1 /* 2 * A speed-improved simplex noise algorithm for 2D, 3D and 4D in Java. 3 * 4 * Based on example code by Stefan Gustavson (stegu@itn.liu.se). 5 * Optimisations by Peter Eastman (peastman@drizzle.stanford.edu). 6 * Better rank ordering method for 4D by Stefan Gustavson in 2012. 7 * 8 * This could be speeded up even further, but it's useful as it is. 9 * 10 * Version 2012-03-09 11 * 12 * This code was placed in the public domain by its original author, 13 * Stefan Gustavson. You may use it as you see fit, but 14 * attribution is appreciated. 15 * 16 * http://webstaff.itn.liu.se/~stegu/simplexnoise/SimplexNoise.java 17 */ 18 19 module voxelman.math.simplex; 20 21 import std.math; 22 23 // Simplex noise in 2D, 3D and 4D 24 struct SimplexNoise { 25 // 2D simplex noise 26 static double noise(double xin, double yin) pure @nogc nothrow { 27 double n0, n1, n2; // Noise contributions from the three corners 28 // Skew the input space to determine which simplex cell we're in 29 double s = (xin+yin)*F2; // Hairy factor for 2D 30 int i = fastfloor(xin+s); 31 int j = fastfloor(yin+s); 32 double t = (i+j)*G2; 33 double X0 = i-t; // Unskew the cell origin back to (x,y) space 34 double Y0 = j-t; 35 double x0 = xin-X0; // The x,y distances from the cell origin 36 double y0 = yin-Y0; 37 // For the 2D case, the simplex shape is an equilateral triangle. 38 // Determine which simplex we are in. 39 int i1, j1; // Offsets for second (middle) corner of simplex in (i,j) coords 40 if(x0>y0) {i1=1; j1=0;} // lower triangle, XY order: (0,0)->(1,0)->(1,1) 41 else {i1=0; j1=1;} // upper triangle, YX order: (0,0)->(0,1)->(1,1) 42 // A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and 43 // a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where 44 // c = (3-sqrt(3))/6 45 double x1 = x0 - i1 + G2; // Offsets for middle corner in (x,y) unskewed coords 46 double y1 = y0 - j1 + G2; 47 double x2 = x0 - 1.0 + 2.0 * G2; // Offsets for last corner in (x,y) unskewed coords 48 double y2 = y0 - 1.0 + 2.0 * G2; 49 // Work out the hashed gradient indices of the three simplex corners 50 int ii = i & 255; 51 int jj = j & 255; 52 int gi0 = permMod12[ii+perm[jj]]; 53 int gi1 = permMod12[ii+i1+perm[jj+j1]]; 54 int gi2 = permMod12[ii+1+perm[jj+1]]; 55 // Calculate the contribution from the three corners 56 double t0 = 0.5 - x0*x0-y0*y0; 57 if(t0<0) n0 = 0.0; 58 else { 59 t0 *= t0; 60 n0 = t0 * t0 * dot(grad3[gi0], x0, y0); // (x,y) of grad3 used for 2D gradient 61 } 62 double t1 = 0.5 - x1*x1-y1*y1; 63 if(t1<0) n1 = 0.0; 64 else { 65 t1 *= t1; 66 n1 = t1 * t1 * dot(grad3[gi1], x1, y1); 67 } 68 double t2 = 0.5 - x2*x2-y2*y2; 69 if(t2<0) n2 = 0.0; 70 else { 71 t2 *= t2; 72 n2 = t2 * t2 * dot(grad3[gi2], x2, y2); 73 } 74 // Add contributions from each corner to get the final noise value. 75 // The result is scaled to return values in the interval [-1,1]. 76 return 70.0 * (n0 + n1 + n2); 77 } 78 79 // 3D simplex noise 80 static double noise(double xin, double yin, double zin) pure @nogc nothrow { 81 double n0, n1, n2, n3; // Noise contributions from the four corners 82 // Skew the input space to determine which simplex cell we're in 83 double s = (xin+yin+zin)*F3; // Very nice and simple skew factor for 3D 84 int i = fastfloor(xin+s); 85 int j = fastfloor(yin+s); 86 int k = fastfloor(zin+s); 87 double t = (i+j+k)*G3; 88 double X0 = i-t; // Unskew the cell origin back to (x,y,z) space 89 double Y0 = j-t; 90 double Z0 = k-t; 91 double x0 = xin-X0; // The x,y,z distances from the cell origin 92 double y0 = yin-Y0; 93 double z0 = zin-Z0; 94 // For the 3D case, the simplex shape is a slightly irregular tetrahedron. 95 // Determine which simplex we are in. 96 int i1, j1, k1; // Offsets for second corner of simplex in (i,j,k) coords 97 int i2, j2, k2; // Offsets for third corner of simplex in (i,j,k) coords 98 if(x0>=y0) { 99 if(y0>=z0) 100 { i1=1; j1=0; k1=0; i2=1; j2=1; k2=0; } // X Y Z order 101 else if(x0>=z0) { i1=1; j1=0; k1=0; i2=1; j2=0; k2=1; } // X Z Y order 102 else { i1=0; j1=0; k1=1; i2=1; j2=0; k2=1; } // Z X Y order 103 } 104 else { // x0<y0 105 if(y0<z0) { i1=0; j1=0; k1=1; i2=0; j2=1; k2=1; } // Z Y X order 106 else if(x0<z0) { i1=0; j1=1; k1=0; i2=0; j2=1; k2=1; } // Y Z X order 107 else { i1=0; j1=1; k1=0; i2=1; j2=1; k2=0; } // Y X Z order 108 } 109 // A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in (x,y,z), 110 // a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in (x,y,z), and 111 // a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in (x,y,z), where 112 // c = 1/6. 113 double x1 = x0 - i1 + G3; // Offsets for second corner in (x,y,z) coords 114 double y1 = y0 - j1 + G3; 115 double z1 = z0 - k1 + G3; 116 double x2 = x0 - i2 + 2.0*G3; // Offsets for third corner in (x,y,z) coords 117 double y2 = y0 - j2 + 2.0*G3; 118 double z2 = z0 - k2 + 2.0*G3; 119 double x3 = x0 - 1.0 + 3.0*G3; // Offsets for last corner in (x,y,z) coords 120 double y3 = y0 - 1.0 + 3.0*G3; 121 double z3 = z0 - 1.0 + 3.0*G3; 122 // Work out the hashed gradient indices of the four simplex corners 123 int ii = i & 255; 124 int jj = j & 255; 125 int kk = k & 255; 126 int gi0 = permMod12[ii+perm[jj+perm[kk]]]; 127 int gi1 = permMod12[ii+i1+perm[jj+j1+perm[kk+k1]]]; 128 int gi2 = permMod12[ii+i2+perm[jj+j2+perm[kk+k2]]]; 129 int gi3 = permMod12[ii+1+perm[jj+1+perm[kk+1]]]; 130 // Calculate the contribution from the four corners 131 double t0 = 0.6 - x0*x0 - y0*y0 - z0*z0; 132 if(t0<0) n0 = 0.0; 133 else { 134 t0 *= t0; 135 n0 = t0 * t0 * dot(grad3[gi0], x0, y0, z0); 136 } 137 double t1 = 0.6 - x1*x1 - y1*y1 - z1*z1; 138 if(t1<0) n1 = 0.0; 139 else { 140 t1 *= t1; 141 n1 = t1 * t1 * dot(grad3[gi1], x1, y1, z1); 142 } 143 double t2 = 0.6 - x2*x2 - y2*y2 - z2*z2; 144 if(t2<0) n2 = 0.0; 145 else { 146 t2 *= t2; 147 n2 = t2 * t2 * dot(grad3[gi2], x2, y2, z2); 148 } 149 double t3 = 0.6 - x3*x3 - y3*y3 - z3*z3; 150 if(t3<0) n3 = 0.0; 151 else { 152 t3 *= t3; 153 n3 = t3 * t3 * dot(grad3[gi3], x3, y3, z3); 154 } 155 // Add contributions from each corner to get the final noise value. 156 // The result is scaled to stay just inside [-1,1] 157 return 32.0*(n0 + n1 + n2 + n3); 158 } 159 160 // 4D simplex noise, better simplex rank ordering method 2012-03-09 161 static double noise(double x, double y, double z, double w) pure @nogc nothrow { 162 163 double n0, n1, n2, n3, n4; // Noise contributions from the five corners 164 // Skew the (x,y,z,w) space to determine which cell of 24 simplices we're in 165 double s = (x + y + z + w) * F4; // Factor for 4D skewing 166 int i = fastfloor(x + s); 167 int j = fastfloor(y + s); 168 int k = fastfloor(z + s); 169 int l = fastfloor(w + s); 170 double t = (i + j + k + l) * G4; // Factor for 4D unskewing 171 double X0 = i - t; // Unskew the cell origin back to (x,y,z,w) space 172 double Y0 = j - t; 173 double Z0 = k - t; 174 double W0 = l - t; 175 double x0 = x - X0; // The x,y,z,w distances from the cell origin 176 double y0 = y - Y0; 177 double z0 = z - Z0; 178 double w0 = w - W0; 179 // For the 4D case, the simplex is a 4D shape I won't even try to describe. 180 // To find out which of the 24 possible simplices we're in, we need to 181 // determine the magnitude ordering of x0, y0, z0 and w0. 182 // Six pair-wise comparisons are performed between each possible pair 183 // of the four coordinates, and the results are used to rank the numbers. 184 int rankx = 0; 185 int ranky = 0; 186 int rankz = 0; 187 int rankw = 0; 188 if(x0 > y0) rankx++; else ranky++; 189 if(x0 > z0) rankx++; else rankz++; 190 if(x0 > w0) rankx++; else rankw++; 191 if(y0 > z0) ranky++; else rankz++; 192 if(y0 > w0) ranky++; else rankw++; 193 if(z0 > w0) rankz++; else rankw++; 194 int i1, j1, k1, l1; // The integer offsets for the second simplex corner 195 int i2, j2, k2, l2; // The integer offsets for the third simplex corner 196 int i3, j3, k3, l3; // The integer offsets for the fourth simplex corner 197 // [rankx, ranky, rankz, rankw] is a 4-vector with the numbers 0, 1, 2 and 3 198 // in some order. We use a thresholding to set the coordinates in turn. 199 // Rank 3 denotes the largest coordinate. 200 i1 = rankx >= 3 ? 1 : 0; 201 j1 = ranky >= 3 ? 1 : 0; 202 k1 = rankz >= 3 ? 1 : 0; 203 l1 = rankw >= 3 ? 1 : 0; 204 // Rank 2 denotes the second largest coordinate. 205 i2 = rankx >= 2 ? 1 : 0; 206 j2 = ranky >= 2 ? 1 : 0; 207 k2 = rankz >= 2 ? 1 : 0; 208 l2 = rankw >= 2 ? 1 : 0; 209 // Rank 1 denotes the second smallest coordinate. 210 i3 = rankx >= 1 ? 1 : 0; 211 j3 = ranky >= 1 ? 1 : 0; 212 k3 = rankz >= 1 ? 1 : 0; 213 l3 = rankw >= 1 ? 1 : 0; 214 // The fifth corner has all coordinate offsets = 1, so no need to compute that. 215 double x1 = x0 - i1 + G4; // Offsets for second corner in (x,y,z,w) coords 216 double y1 = y0 - j1 + G4; 217 double z1 = z0 - k1 + G4; 218 double w1 = w0 - l1 + G4; 219 double x2 = x0 - i2 + 2.0*G4; // Offsets for third corner in (x,y,z,w) coords 220 double y2 = y0 - j2 + 2.0*G4; 221 double z2 = z0 - k2 + 2.0*G4; 222 double w2 = w0 - l2 + 2.0*G4; 223 double x3 = x0 - i3 + 3.0*G4; // Offsets for fourth corner in (x,y,z,w) coords 224 double y3 = y0 - j3 + 3.0*G4; 225 double z3 = z0 - k3 + 3.0*G4; 226 double w3 = w0 - l3 + 3.0*G4; 227 double x4 = x0 - 1.0 + 4.0*G4; // Offsets for last corner in (x,y,z,w) coords 228 double y4 = y0 - 1.0 + 4.0*G4; 229 double z4 = z0 - 1.0 + 4.0*G4; 230 double w4 = w0 - 1.0 + 4.0*G4; 231 // Work out the hashed gradient indices of the five simplex corners 232 int ii = i & 255; 233 int jj = j & 255; 234 int kk = k & 255; 235 int ll = l & 255; 236 int gi0 = perm[ii+perm[jj+perm[kk+perm[ll]]]] % 32; 237 int gi1 = perm[ii+i1+perm[jj+j1+perm[kk+k1+perm[ll+l1]]]] % 32; 238 int gi2 = perm[ii+i2+perm[jj+j2+perm[kk+k2+perm[ll+l2]]]] % 32; 239 int gi3 = perm[ii+i3+perm[jj+j3+perm[kk+k3+perm[ll+l3]]]] % 32; 240 int gi4 = perm[ii+1+perm[jj+1+perm[kk+1+perm[ll+1]]]] % 32; 241 // Calculate the contribution from the five corners 242 double t0 = 0.6 - x0*x0 - y0*y0 - z0*z0 - w0*w0; 243 if(t0<0) n0 = 0.0; 244 else { 245 t0 *= t0; 246 n0 = t0 * t0 * dot(grad4[gi0], x0, y0, z0, w0); 247 } 248 double t1 = 0.6 - x1*x1 - y1*y1 - z1*z1 - w1*w1; 249 if(t1<0) n1 = 0.0; 250 else { 251 t1 *= t1; 252 n1 = t1 * t1 * dot(grad4[gi1], x1, y1, z1, w1); 253 } 254 double t2 = 0.6 - x2*x2 - y2*y2 - z2*z2 - w2*w2; 255 if(t2<0) n2 = 0.0; 256 else { 257 t2 *= t2; 258 n2 = t2 * t2 * dot(grad4[gi2], x2, y2, z2, w2); 259 } 260 double t3 = 0.6 - x3*x3 - y3*y3 - z3*z3 - w3*w3; 261 if(t3<0) n3 = 0.0; 262 else { 263 t3 *= t3; 264 n3 = t3 * t3 * dot(grad4[gi3], x3, y3, z3, w3); 265 } 266 double t4 = 0.6 - x4*x4 - y4*y4 - z4*z4 - w4*w4; 267 if(t4<0) n4 = 0.0; 268 else { 269 t4 *= t4; 270 n4 = t4 * t4 * dot(grad4[gi4], x4, y4, z4, w4); 271 } 272 // Sum up and scale the result to cover the range [-1,1] 273 return 27.0 * (n0 + n1 + n2 + n3 + n4); 274 } 275 } 276 277 private: 278 279 // Skewing and unskewing factors for 2, 3, and 4 dimensions 280 enum double F2 = 0.5 * (sqrt(3.0) - 1.0); 281 enum double G2 = (3.0 - sqrt(3.0)) / 6.0; 282 enum double F3 = 1.0 / 3.0; 283 enum double G3 = 1.0 / 6.0; 284 enum double F4 = (sqrt(5.0) - 1.0) / 4.0; 285 enum double G4 = (5.0 - sqrt(5.0)) / 20.0; 286 287 pure const @nogc nothrow: 288 pragma(inline, true): 289 290 // This method is a *lot* faster than using (int)floor(x) 291 int fastfloor(double x) { 292 int xi = cast(int)x; 293 return x<xi ? xi-1 : xi; 294 } 295 296 double dot(const byte[] g, double x, const double y) { 297 return g[0]*x + g[1]*y; 298 } 299 300 double dot(const byte[] g, double x, double y, double z) { 301 return g[0]*x + g[1]*y + g[2]*z; 302 } 303 304 double dot(const byte[] g, double x, double y, double z, double w) { 305 return g[0]*x + g[1]*y + g[2]*z + g[3]*w; 306 } 307 308 immutable byte[3][12] grad3 = [ 309 [1, 1, 0], [-1, 1, 0], [1,-1, 0], [-1,-1, 0], 310 [1, 0, 1], [-1, 0, 1], [1, 0,-1], [-1, 0,-1], 311 [0, 1, 1], [ 0,-1, 1], [0, 1,-1], [ 0,-1,-1]]; 312 immutable byte[4][32] grad4 = [ 313 [ 0, 1, 1, 1], [ 0, 1, 1,-1], [ 0, 1,-1, 1], [ 0, 1,-1,-1], 314 [ 0,-1, 1, 1], [ 0,-1, 1,-1], [ 0,-1,-1, 1], [ 0,-1,-1,-1], 315 [ 1, 0, 1, 1], [ 1, 0, 1,-1], [ 1, 0,-1, 1], [ 1, 0,-1,-1], 316 [-1, 0, 1, 1], [-1, 0, 1,-1], [-1, 0,-1, 1], [-1, 0,-1,-1], 317 [ 1, 1, 0, 1], [ 1, 1, 0,-1], [ 1,-1, 0, 1], [ 1,-1, 0,-1], 318 [-1, 1, 0, 1], [-1, 1, 0,-1], [-1,-1, 0, 1], [-1,-1, 0,-1], 319 [ 1, 1, 1, 0], [ 1, 1,-1, 0], [ 1,-1, 1, 0], [ 1,-1,-1, 0], 320 [-1, 1, 1, 0], [-1, 1,-1, 0], [-1,-1, 1, 0], [-1,-1,-1, 0]]; 321 322 immutable ushort[512] perm = [ 323 151,160,137,91, 90,15,131,13,201, 95,96, 53,194,233,7,225,140, 36,103,30,69, 324 142, 8, 99,37,240,21,10,23, 190,6,148,247,120,234, 75, 0, 26,197,62, 94,252, 325 219,203,117, 35,11,32, 57,177,33, 88,237,149, 56,87,174, 20,125,136,171,168, 326 68,175, 74,165,71,134,139, 48,27,166,77,146,158,231, 83,111,229,122, 60,211, 327 133,230,220,105, 92,41,55,46,245, 40,244,102,143, 54,65,25, 63,161,1,216,80, 328 73,209, 76,132,187,208, 89, 18,169,200,196,135,130,116, 188,159, 86,164,100, 329 109,198,173,186, 3,64, 52,217,226,250,124,123, 5,202, 38,147,118,126,255,82, 330 85,212, 207,206,59,227,47, 16,58,17, 182,189,28,42,223,183,170, 213,119,248, 331 152, 2,44,154,163, 70,221,153,101,155,167, 43,172, 9,129, 22, 39,253, 19,98, 332 108,110, 79,113, 224,232, 178,185,112, 104,218,246, 97,228, 251, 34,242,193, 333 238,210,144, 12,191,179,162,241, 81, 51,145,235,249, 14,239,107, 49,192,214, 334 31,181,199, 106,157,184, 84,204,176,115,121, 50, 45,127, 4,150, 254,138,236, 335 205, 93, 222,114, 67, 29, 24, 72,243, 141,128, 195, 78, 66,215, 61, 156,180, 336 337 151,160,137,91, 90,15,131,13,201, 95,96, 53,194,233,7,225,140, 36,103,30,69, 338 142, 8, 99,37,240,21,10,23, 190,6,148,247,120,234, 75, 0, 26,197,62, 94,252, 339 219,203,117, 35,11,32, 57,177,33, 88,237,149, 56,87,174, 20,125,136,171,168, 340 68,175, 74,165,71,134,139, 48,27,166,77,146,158,231, 83,111,229,122, 60,211, 341 133,230,220,105, 92,41,55,46,245, 40,244,102,143, 54,65,25, 63,161,1,216,80, 342 73,209, 76,132,187,208, 89, 18,169,200,196,135,130,116, 188,159, 86,164,100, 343 109,198,173,186, 3,64, 52,217,226,250,124,123, 5,202, 38,147,118,126,255,82, 344 85,212, 207,206,59,227,47, 16,58,17, 182,189,28,42,223,183,170, 213,119,248, 345 152, 2,44,154,163, 70,221,153,101,155,167, 43,172, 9,129, 22, 39,253, 19,98, 346 108,110, 79,113, 224,232, 178,185,112, 104,218,246, 97,228, 251, 34,242,193, 347 238,210,144, 12,191,179,162,241, 81, 51,145,235,249, 14,239,107, 49,192,214, 348 31,181,199, 106,157,184, 84,204,176,115,121, 50, 45,127, 4,150, 254,138,236, 349 205, 93, 222,114, 67, 29, 24, 72,243, 141,128, 195, 78, 66,215, 61, 156,180]; 350 351 // permMod12[i] = (short)(perm[i] % 12); 352 immutable ushort[512] permMod12 = [ 353 7,4,5,7,6,3,11,1,9,11,0,5,2,5,7,9,8,0,7,6,9,10,8,3,1,0,9,10,11,10,6,4,7,0,6, 354 3,0,2,5,2,10,0,3,11, 9,11,11,8,9,9,9,4,9,5,8,3,6,8,5,4,3,0,8,7,2,9,11,2,7,0, 355 3,10,5,2,2,3,11,3,1,2,0,7,1,2,4,9,8,5,7,10,5,4,4,6,11,6,5,1,3,5,1,0,8,1,5,4, 356 0,7,4,5,6,1,8,4,3,10,8,8,3,2,8,4,1,6,5,6,3,4,4, 1,10,10,4,3,5,10,2,3,10,6,3, 357 10,1,8,3,2,11,11,11,4,10,5,2,9,4,6,7,3,2,9,11,8,8,2,8,10,7,10,5,9,5,11,11,7, 358 4,9,9,10,3,1,7,2,0,2,7,5,8, 4,10,5,4,8,2,6,1,0,11,10,2,1,10,6,0,0,11,11,6,1, 359 9,3,1,7,9,2,11,11,1,0,10,7,1,7,10,1,4,0,0,8,7,1,2,9,7,4,6,2,6,8,1,9,6,6,7,5, 360 0,0,3,9,8,3,6,6,11,1,0,0,7,4,5,7,6,3,11,1,9,11,0,5,2,5,7,9,8,0,7,6,9,10,8,3, 361 1,0,9,10,11,10,6,4,7,0,6,3,0,2,5,2,10,0,3,11,9,11,11,8,9, 9,9,4,9,5,8,3,6,8, 362 5,4,3,0,8,7,2,9,11,2,7,0,3,10,5,2,2,3,11,3,1,2,0,7,1,2,4,9,8,5,7,10,5,4,4,6, 363 11,6,5,1,3, 5,1,0,8,1,5,4,0, 7,4,5,6,1,8,4,3,10,8,8,3,2,8,4,1,6,5,6,3,4,4,1, 364 10,10,4,3,5,10,2,3,10,6,3,10,1,8,3,2,11,11,11,4,10,5,2,9,4,6,7,3,2,9,11,8,8, 365 2,8,10,7,10,5,9,5,11,11,7,4,9,9,10,3,1,7,2, 0,2,7,5,8,4,10,5,4,8,2,6,1,0,11, 366 10,2,1,10,6,0,0,11,11,6,1,9,3,1,7,9,2,11,11,1,0,10,7,1,7,10,1,4,0,0,8,7,1,2, 367 9,7,4,6,2,6,8,1,9,6,6,7,5,0,0,3,9,8,3,6,6,11,1,0,0];