001package squidpony.squidmath;
002
003/**
004 * This is Ken Perlin's third revision of his noise function. It is sometimes
005 * referred to as "Simplex Noise". Results are bound by (-1, 1) inclusive.
006 *
007 *
008 * It is significantly faster than his earlier versions. This particular version
009 * was originally from Stefan Gustavson. This is much preferred to the earlier
010 * versions of Perlin Noise due to the reasons noted in the articles:
011 * <ul>
012 * <li>http://www.itn.liu.se/~stegu/simplexnoise/simplexnoise.pdf</li>
013 * <li>http://webstaff.itn.liu.se/~stegu/TNM022-2005/perlinnoiselinks/ch02.pdf</li>
014 * </ul>
015 *
016 */
017public class PerlinNoise {
018
019    private static final int grad3[][] = {{1, 1, 0}, {-1, 1, 0}, {1, -1, 0},
020    {-1, -1, 0}, {1, 0, 1}, {-1, 0, 1},
021    {1, 0, -1}, {-1, 0, -1}, {0, 1, 1},
022    {0, -1, 1}, {0, 1, -1}, {0, -1, -1}};
023    private static final int grad4[][] = {{0, 1, 1, 1}, {0, 1, 1, -1},
024    {0, 1, -1, 1}, {0, 1, -1, -1},
025    {0, -1, 1, 1}, {0, -1, 1, -1},
026    {0, -1, -1, 1}, {0, -1, -1, -1},
027    {1, 0, 1, 1}, {1, 0, 1, -1},
028    {1, 0, -1, 1}, {1, 0, -1, -1},
029    {-1, 0, 1, 1}, {-1, 0, 1, -1},
030    {-1, 0, -1, 1}, {-1, 0, -1, -1},
031    {1, 1, 0, 1}, {1, 1, 0, -1},
032    {1, -1, 0, 1}, {1, -1, 0, -1},
033    {-1, 1, 0, 1}, {-1, 1, 0, -1},
034    {-1, -1, 0, 1}, {-1, -1, 0, -1},
035    {1, 1, 1, 0}, {1, 1, -1, 0},
036    {1, -1, 1, 0}, {1, -1, -1, 0},
037    {-1, 1, 1, 0}, {-1, 1, -1, 0},
038    {-1, -1, 1, 0}, {-1, -1, -1, 0}};
039    private static final int p[] = {151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96,
040        53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142,
041        8, 99, 37, 240, 21, 10, 23, 190, 6, 148, 247,
042        120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203,
043        117, 35, 11, 32, 57, 177, 33, 88, 237, 149, 56,
044        87, 174, 20, 125, 136, 171, 168, 68, 175, 74,
045        165, 71, 134, 139, 48, 27, 166, 77, 146, 158,
046        231, 83, 111, 229, 122, 60, 211, 133, 230, 220,
047        105, 92, 41, 55, 46, 245, 40, 244, 102, 143, 54,
048        65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132,
049        187, 208, 89, 18, 169, 200, 196, 135, 130, 116,
050        188, 159, 86, 164, 100, 109, 198, 173, 186, 3,
051        64, 52, 217, 226, 250, 124, 123, 5, 202, 38, 147,
052        118, 126, 255, 82, 85, 212, 207, 206, 59, 227,
053        47, 16, 58, 17, 182, 189, 28, 42, 223, 183, 170,
054        213, 119, 248, 152, 2, 44, 154, 163, 70, 221,
055        153, 101, 155, 167, 43, 172, 9, 129, 22, 39, 253,
056        19, 98, 108, 110, 79, 113, 224, 232, 178, 185,
057        112, 104, 218, 246, 97, 228, 251, 34, 242, 193,
058        238, 210, 144, 12, 191, 179, 162, 241, 81, 51,
059        145, 235, 249, 14, 239, 107, 49, 192, 214, 31,
060        181, 199, 106, 157, 184, 84, 204, 176, 115, 121,
061        50, 45, 127, 4, 150, 254, 138, 236, 205, 93, 222,
062        114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66,
063        215, 61, 156, 180};
064    private static final double F2 = 0.5 * (Math.sqrt(3.0) - 1.0);
065    private static final double G2 = (3.0 - Math.sqrt(3.0)) / 6.0;
066    private static final double F3 = 1.0 / 3.0;
067    private static final double G3 = 1.0 / 6.0;
068    private static final double F4 = (Math.sqrt(5.0) - 1.0) / 4.0;
069    private static final double G4 = (5.0 - Math.sqrt(5.0)) / 20.0;
070    // To remove the need for index wrapping, double the permutation table
071    // length
072    private static final int perm[] = new int[512];
073
074    static {
075        for (int i = 0; i < 512; i++) {
076            perm[i] = p[i & 255];
077        }
078    }
079    private PerlinNoise()
080    {
081
082    }
083    // A lookup table to traverse the simplex around a given point in 4D.
084    // Details can be found where this table is used, in the 4D noise method.
085    private static final int simplex[][]
086            = {{0, 1, 2, 3}, {0, 1, 3, 2}, {0, 0, 0, 0}, {0, 2, 3, 1},
087            {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 2, 3, 0},
088            {0, 2, 1, 3}, {0, 0, 0, 0}, {0, 3, 1, 2}, {0, 3, 2, 1},
089            {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {1, 3, 2, 0},
090            {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
091            {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
092            {1, 2, 0, 3}, {0, 0, 0, 0}, {1, 3, 0, 2}, {0, 0, 0, 0},
093            {0, 0, 0, 0}, {0, 0, 0, 0}, {2, 3, 0, 1}, {2, 3, 1, 0},
094            {1, 0, 2, 3}, {1, 0, 3, 2}, {0, 0, 0, 0}, {0, 0, 0, 0},
095            {0, 0, 0, 0}, {2, 0, 3, 1}, {0, 0, 0, 0}, {2, 1, 3, 0},
096            {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
097            {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
098            {2, 0, 1, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
099            {3, 0, 1, 2}, {3, 0, 2, 1}, {0, 0, 0, 0}, {3, 1, 2, 0},
100            {2, 1, 0, 3}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0},
101            {3, 1, 0, 2}, {0, 0, 0, 0}, {3, 2, 0, 1}, {3, 2, 1, 0}};
102
103    private static double dot(int g[], double x, double y) {
104        return g[0] * x + g[1] * y;
105    }
106
107    private static double dot(int g[], double x, double y, double z) {
108        return g[0] * x + g[1] * y + g[2] * z;
109    }
110
111    private static double dot(int g[], double x, double y, double z, double w) {
112        return g[0] * x + g[1] * y + g[2] * z + g[3] * w;
113    }
114
115    /**
116     * 2D simplex noise.
117     *
118     * @param xin x input
119     * @param yin y input
120     * @return noise from -1.0 to 1.0, inclusive
121     */
122    public static double noise(double xin, double yin) {
123        double noise0, noise1, noise2; // from the three corners
124        // Skew the input space to determine which simplex cell we're in
125        double skew = (xin + yin) * F2; // Hairy factor for 2D
126        int i = (int) Math.floor(xin + skew);
127        int j = (int) Math.floor(yin + skew);
128        double t = (i + j) * G2;
129        double X0 = i - t; // Unskew the cell origin back to (x,y) space
130        double Y0 = j - t;
131        double x0 = xin - X0; // The x,y distances from the cell origin
132        double y0 = yin - Y0;
133        // For the 2D case, the simplex shape is an equilateral triangle.
134        // Determine which simplex we are in.
135        int i1, j1; // Offsets for second (middle) corner of simplex in (i,j)
136        // coords
137        if (x0 > y0) {
138            i1 = 1;
139            j1 = 0;
140        } // lower triangle, XY order: (0,0)->(1,0)->(1,1)
141        else {
142            i1 = 0;
143            j1 = 1;
144        } // upper triangle, YX order: (0,0)->(0,1)->(1,1)
145        // A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and
146        // a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y),
147        // where
148        // c = (3-sqrt(3))/6
149        double x1 = x0 - i1 + G2; // Offsets for middle corner in (x,y)
150        // unskewed coords
151        double y1 = y0 - j1 + G2;
152        double x2 = x0 - 1.0 + 2.0 * G2; // Offsets for last corner in (x,y)
153        // unskewed coords
154        double y2 = y0 - 1.0 + 2.0 * G2;
155        // Work out the hashed gradient indices of the three simplex corners
156        int ii = i & 255;
157        int jj = j & 255;
158        int gi0 = perm[ii + perm[jj]] % 12;
159        int gi1 = perm[ii + i1 + perm[jj + j1]] % 12;
160        int gi2 = perm[ii + 1 + perm[jj + 1]] % 12;
161        // Calculate the contribution from the three corners
162        double t0 = 0.5 - x0 * x0 - y0 * y0;
163        if (t0 < 0) {
164            noise0 = 0.0;
165        } else {
166            t0 *= t0;
167            noise0 = t0 * t0 * dot(grad3[gi0], x0, y0); // (x,y) of grad3 used
168            // for 2D gradient
169        }
170        double t1 = 0.5 - x1 * x1 - y1 * y1;
171        if (t1 < 0) {
172            noise1 = 0.0;
173        } else {
174            t1 *= t1;
175            noise1 = t1 * t1 * dot(grad3[gi1], x1, y1);
176        }
177        double t2 = 0.5 - x2 * x2 - y2 * y2;
178        if (t2 < 0) {
179            noise2 = 0.0;
180        } else {
181            t2 *= t2;
182            noise2 = t2 * t2 * dot(grad3[gi2], x2, y2);
183        }
184        // Add contributions from each corner to get the final noise value.
185        // The result is scaled to return values in the interval [-1,1].
186        return 70.0 * (noise0 + noise1 + noise2);
187    }
188
189    /**
190     * 3D simplex noise.
191     *
192     * @param xin X input
193     * @param yin Y input
194     * @param zin Z input
195     * @return noise from -1.0 to 1.0, inclusive
196     */
197    public static double noise(double xin, double yin, double zin) {
198        double n0, n1, n2, n3; // Noise contributions from the four corners
199        // Skew the input space to determine which simplex cell we're in
200        double s = (xin + yin + zin) * F3; // Very nice and simple skew
201        // factor for 3D
202        int i = (int) Math.floor(xin + s);
203        int j = (int) Math.floor(yin + s);
204        int k = (int) Math.floor(zin + s);
205        double t = (i + j + k) * G3;
206        double X0 = i - t; // Unskew the cell origin back to (x,y,z) space
207        double Y0 = j - t;
208        double Z0 = k - t;
209        double x0 = xin - X0; // The x,y,z distances from the cell origin
210        double y0 = yin - Y0;
211        double z0 = zin - Z0;
212        // For the 3D case, the simplex shape is a slightly irregular
213        // tetrahedron.
214        // Determine which simplex we are in.
215        int i1, j1, k1; // Offsets for second corner of simplex in (i,j,k)
216        // coords
217        int i2, j2, k2; // Offsets for third corner of simplex in (i,j,k)
218        // coords
219        if (x0 >= y0) {
220            if (y0 >= z0) {
221                i1 = 1;
222                j1 = 0;
223                k1 = 0;
224                i2 = 1;
225                j2 = 1;
226                k2 = 0;
227            } // X Y Z order
228            else if (x0 >= z0) {
229                i1 = 1;
230                j1 = 0;
231                k1 = 0;
232                i2 = 1;
233                j2 = 0;
234                k2 = 1;
235            } // X Z Y order
236            else {
237                i1 = 0;
238                j1 = 0;
239                k1 = 1;
240                i2 = 1;
241                j2 = 0;
242                k2 = 1;
243            } // Z X Y order
244        } else { // x0<y0
245            if (y0 < z0) {
246                i1 = 0;
247                j1 = 0;
248                k1 = 1;
249                i2 = 0;
250                j2 = 1;
251                k2 = 1;
252            } // Z Y X order
253            else if (x0 < z0) {
254                i1 = 0;
255                j1 = 1;
256                k1 = 0;
257                i2 = 0;
258                j2 = 1;
259                k2 = 1;
260            } // Y Z X order
261            else {
262                i1 = 0;
263                j1 = 1;
264                k1 = 0;
265                i2 = 1;
266                j2 = 1;
267                k2 = 0;
268            } // Y X Z order
269        }
270        // A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in
271        // (x,y,z),
272        // a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in
273        // (x,y,z), and
274        // a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in
275        // (x,y,z), where
276        // c = 1/6.
277        double x1 = x0 - i1 + G3; // Offsets for second corner in (x,y,z)
278        // coords
279        double y1 = y0 - j1 + G3;
280        double z1 = z0 - k1 + G3;
281        double x2 = x0 - i2 + 2.0 * G3; // Offsets for third corner in
282        // (x,y,z) coords
283        double y2 = y0 - j2 + 2.0 * G3;
284        double z2 = z0 - k2 + 2.0 * G3;
285        double x3 = x0 - 1.0 + 3.0 * G3; // Offsets for last corner in
286        // (x,y,z) coords
287        double y3 = y0 - 1.0 + 3.0 * G3;
288        double z3 = z0 - 1.0 + 3.0 * G3;
289        // Work out the hashed gradient indices of the four simplex corners
290        int ii = i & 255;
291        int jj = j & 255;
292        int kk = k & 255;
293        int gi0 = perm[ii + perm[jj + perm[kk]]] % 12;
294        int gi1 = perm[ii + i1 + perm[jj + j1 + perm[kk + k1]]] % 12;
295        int gi2 = perm[ii + i2 + perm[jj + j2 + perm[kk + k2]]] % 12;
296        int gi3 = perm[ii + 1 + perm[jj + 1 + perm[kk + 1]]] % 12;
297        // Calculate the contribution from the four corners
298        double t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0;
299        if (t0 < 0) {
300            n0 = 0.0;
301        } else {
302            t0 *= t0;
303            n0 = t0 * t0 * dot(grad3[gi0], x0, y0, z0);
304        }
305        double t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1;
306        if (t1 < 0) {
307            n1 = 0.0;
308        } else {
309            t1 *= t1;
310            n1 = t1 * t1 * dot(grad3[gi1], x1, y1, z1);
311        }
312        double t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2;
313        if (t2 < 0) {
314            n2 = 0.0;
315        } else {
316            t2 *= t2;
317            n2 = t2 * t2 * dot(grad3[gi2], x2, y2, z2);
318        }
319        double t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3;
320        if (t3 < 0) {
321            n3 = 0.0;
322        } else {
323            t3 *= t3;
324            n3 = t3 * t3 * dot(grad3[gi3], x3, y3, z3);
325        }
326        // Add contributions from each corner to get the final noise value.
327        // The result is scaled to stay just inside [-1,1]
328        return 32.0 * (n0 + n1 + n2 + n3);
329    }
330
331    /**
332     * 4D simplex noise.
333     *
334     * @param x X position
335     * @param y Y position
336     * @param z Z position
337     * @param w Fourth-dimensional position. It is I, the Fourth-Dimensional Ziltoid the Omniscient!
338     * @return noise from -1.0 to 1.0, inclusive
339     */
340    public static double noise(double x, double y, double z, double w) {
341        // The skewing and unskewing factors are hairy again for the 4D case
342        double n0, n1, n2, n3, n4; // Noise contributions from the five
343        // corners
344        // Skew the (x,y,z,w) space to determine which cell of 24 simplices
345        // we're in
346        double s = (x + y + z + w) * F4; // Factor for 4D skewing
347        int i = (int) Math.floor(x + s);
348        int j = (int) Math.floor(y + s);
349        int k = (int) Math.floor(z + s);
350        int l = (int) Math.floor(w + s);
351        double t = (i + j + k + l) * G4; // Factor for 4D unskewing
352        double X0 = i - t; // Unskew the cell origin back to (x,y,z,w) space
353        double Y0 = j - t;
354        double Z0 = k - t;
355        double W0 = l - t;
356        double x0 = x - X0; // The x,y,z,w distances from the cell origin
357        double y0 = y - Y0;
358        double z0 = z - Z0;
359        double w0 = w - W0;
360        // For the 4D case, the simplex is a 4D shape I won't even try to
361        // describe.
362        // To find out which of the 24 possible simplices we're in, we need
363        // to
364        // determine the magnitude ordering of x0, y0, z0 and w0.
365        // The method below is a good way of finding the ordering of x,y,z,w
366        // and
367        // then find the correct traversal order for the simplex we’re in.
368        // First, six pair-wise comparisons are performed between each
369        // possible pair
370        // of the four coordinates, and the results are used to add up binary
371        // bits
372        // for an integer index.
373        int c1 = (x0 > y0) ? 32 : 0;
374        int c2 = (x0 > z0) ? 16 : 0;
375        int c3 = (y0 > z0) ? 8 : 0;
376        int c4 = (x0 > w0) ? 4 : 0;
377        int c5 = (y0 > w0) ? 2 : 0;
378        int c6 = (z0 > w0) ? 1 : 0;
379        int c = c1 + c2 + c3 + c4 + c5 + c6;
380        int i1, j1, k1, l1; // The integer offsets for the second simplex
381        // corner
382        int i2, j2, k2, l2; // The integer offsets for the third simplex
383        // corner
384        int i3, j3, k3, l3; // The integer offsets for the fourth simplex
385        // corner
386        // simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some
387        // order.
388        // Many values of c will never occur, since e.g. x>y>z>w makes x<z,
389        // y<w and x<w
390        // impossible. Only the 24 indices which have non-zero entries make
391        // any sense.
392        // We use a thresholding to set the coordinates in turn from the
393        // largest magnitude.
394        // The number 3 in the "simplex" array is at the position of the
395        // largest coordinate.
396        i1 = simplex[c][0] >= 3 ? 1 : 0;
397        j1 = simplex[c][1] >= 3 ? 1 : 0;
398        k1 = simplex[c][2] >= 3 ? 1 : 0;
399        l1 = simplex[c][3] >= 3 ? 1 : 0;
400        // The number 2 in the "simplex" array is at the second largest
401        // coordinate.
402        i2 = simplex[c][0] >= 2 ? 1 : 0;
403        j2 = simplex[c][1] >= 2 ? 1 : 0;
404        k2 = simplex[c][2] >= 2 ? 1 : 0;
405        l2 = simplex[c][3] >= 2 ? 1 : 0;
406        // The number 1 in the "simplex" array is at the second smallest
407        // coordinate.
408        i3 = simplex[c][0] >= 1 ? 1 : 0;
409        j3 = simplex[c][1] >= 1 ? 1 : 0;
410        k3 = simplex[c][2] >= 1 ? 1 : 0;
411        l3 = simplex[c][3] >= 1 ? 1 : 0;
412        // The fifth corner has all coordinate offsets = 1, so no need to
413        // look that up.
414        double x1 = x0 - i1 + G4; // Offsets for second corner in (x,y,z,w)
415        // coords
416        double y1 = y0 - j1 + G4;
417        double z1 = z0 - k1 + G4;
418        double w1 = w0 - l1 + G4;
419        double x2 = x0 - i2 + 2.0 * G4; // Offsets for third corner in
420        // (x,y,z,w) coords
421        double y2 = y0 - j2 + 2.0 * G4;
422        double z2 = z0 - k2 + 2.0 * G4;
423        double w2 = w0 - l2 + 2.0 * G4;
424        double x3 = x0 - i3 + 3.0 * G4; // Offsets for fourth corner in
425        // (x,y,z,w) coords
426        double y3 = y0 - j3 + 3.0 * G4;
427        double z3 = z0 - k3 + 3.0 * G4;
428        double w3 = w0 - l3 + 3.0 * G4;
429        double x4 = x0 - 1.0 + 4.0 * G4; // Offsets for last corner in
430        // (x,y,z,w) coords
431        double y4 = y0 - 1.0 + 4.0 * G4;
432        double z4 = z0 - 1.0 + 4.0 * G4;
433        double w4 = w0 - 1.0 + 4.0 * G4;
434        // Work out the hashed gradient indices of the five simplex corners
435        int ii = i & 255;
436        int jj = j & 255;
437        int kk = k & 255;
438        int ll = l & 255;
439        int gi0 = perm[ii + perm[jj + perm[kk + perm[ll]]]] % 32;
440        int gi1
441                = perm[ii + i1 + perm[jj + j1 + perm[kk + k1 + perm[ll + l1]]]] % 32;
442        int gi2
443                = perm[ii + i2 + perm[jj + j2 + perm[kk + k2 + perm[ll + l2]]]] % 32;
444        int gi3
445                = perm[ii + i3 + perm[jj + j3 + perm[kk + k3 + perm[ll + l3]]]] % 32;
446        int gi4
447                = perm[ii + 1 + perm[jj + 1 + perm[kk + 1 + perm[ll + 1]]]] % 32;
448        // Calculate the contribution from the five corners
449        double t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0 - w0 * w0;
450        if (t0 < 0) {
451            n0 = 0.0;
452        } else {
453            t0 *= t0;
454            n0 = t0 * t0 * dot(grad4[gi0], x0, y0, z0, w0);
455        }
456        double t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1 - w1 * w1;
457        if (t1 < 0) {
458            n1 = 0.0;
459        } else {
460            t1 *= t1;
461            n1 = t1 * t1 * dot(grad4[gi1], x1, y1, z1, w1);
462        }
463        double t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2 - w2 * w2;
464        if (t2 < 0) {
465            n2 = 0.0;
466        } else {
467            t2 *= t2;
468            n2 = t2 * t2 * dot(grad4[gi2], x2, y2, z2, w2);
469        }
470        double t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3 - w3 * w3;
471        if (t3 < 0) {
472            n3 = 0.0;
473        } else {
474            t3 *= t3;
475            n3 = t3 * t3 * dot(grad4[gi3], x3, y3, z3, w3);
476        }
477        double t4 = 0.6 - x4 * x4 - y4 * y4 - z4 * z4 - w4 * w4;
478        if (t4 < 0) {
479            n4 = 0.0;
480        } else {
481            t4 *= t4;
482            n4 = t4 * t4 * dot(grad4[gi4], x4, y4, z4, w4);
483        }
484        // Sum up and scale the result to cover the range [-1,1]
485        return 27.0 * (n0 + n1 + n2 + n3 + n4);
486    }
487}