1 module atelier.core.simplex; 2 3 import std.range : iota; 4 import std.math; 5 6 private { 7 struct Grad { 8 double x, y, z; 9 } 10 11 static Grad[] _grads = [ 12 Grad(1, 1, 0), Grad(-1, 1, 0), Grad(1, -1, 0), Grad(-1, -1, 0), 13 Grad(1, 0, 1), Grad(-1, 0, 1), Grad(1, 0, -1), Grad(-1, 0, -1), 14 Grad(0, 1, 1), Grad(0, -1, 1), Grad(0, 1, -1), Grad(0, -1, -1) 15 ]; 16 17 immutable ubyte[] _p = [ 18 151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7, 19 225, 140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23, 190, 6, 20 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 21 11, 32, 57, 177, 33, 88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 22 68, 175, 74, 165, 71, 134, 139, 48, 27, 166, 77, 146, 158, 231, 83, 23 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 24 244, 102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 25 208, 89, 18, 169, 200, 196, 135, 130, 116, 188, 159, 86, 164, 100, 109, 26 198, 173, 186, 3, 64, 52, 217, 226, 250, 124, 123, 5, 202, 38, 147, 27 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 28 189, 28, 42, 223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 29 221, 153, 101, 155, 167, 43, 172, 9, 129, 22, 39, 253, 19, 98, 108, 30 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, 228, 251, 31 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235, 32 249, 14, 239, 107, 49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 33 176, 115, 121, 50, 45, 127, 4, 150, 254, 138, 236, 205, 93, 222, 114, 34 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180 35 ]; 36 37 enum _permSize = 512; 38 static short[_permSize] _perm, _permMod12; 39 40 static double F2 = 0.5 * (sqrt(3.0) - 1.0); 41 static double G2 = (3.0 - sqrt(3.0)) / 6.0; 42 static double F3 = 1.0 / 3.0; 43 static double G3 = 1.0 / 6.0; 44 static double F4 = (sqrt(5.0) - 1.0) / 4.0; 45 static double G4 = (5.0 - sqrt(5.0)) / 20.0; 46 } 47 48 static this() { 49 static foreach (i; iota(0, _permSize)) { 50 _perm[i] = _p[i & 255]; 51 _permMod12[i] = cast(short)(_perm[i] % 12); 52 } 53 } 54 55 private int fastfloor(double x) { 56 const int xi = cast(int) x; 57 return x < xi ? xi - 1 : xi; 58 } 59 60 private double dot(Grad g, double x, double y) { 61 return g.x * x + g.y * y; 62 } 63 64 private double dot(Grad g, double x, double y, double z) { 65 return g.x * x + g.y * y + g.z * z; 66 } 67 68 /// 1D simplex noise 69 double noise(double x) { 70 return noise(x, 0.0); 71 } 72 73 /// 2D simplex noise 74 double noise(double xin, double yin) { 75 double n0, n1, n2; // Noise contributions from the three corners 76 // Skew the input space to determine which simplex cell we're in 77 const double s = (xin + yin) * F2; // Hairy factor for 2D 78 const int i = fastfloor(xin + s); 79 const int j = fastfloor(yin + s); 80 const double t = (i + j) * G2; 81 const double X0 = i - t; // Unskew the cell origin back to (x,y) space 82 const double Y0 = j - t; 83 double x0 = xin - X0; // The x,y distances from the cell origin 84 double y0 = yin - Y0; 85 // For the 2D case, the simplex shape is an equilateral triangle. 86 // Determine which simplex we are in. 87 int i1, j1; // Offsets for second (middle) corner of simplex in (i,j) coords 88 if (x0 > y0) { 89 i1 = 1; 90 j1 = 0; 91 } // lower triangle, XY order: (0,0)->(1,0)->(1,1) 92 else { 93 i1 = 0; 94 j1 = 1; 95 } // upper triangle, YX order: (0,0)->(0,1)->(1,1) 96 // A step of (1,0) in (i,j) means a step of (1-c,-c) in (x,y), and 97 // a step of (0,1) in (i,j) means a step of (-c,1-c) in (x,y), where 98 // c = (3-sqrt(3))/6 99 double x1 = x0 - i1 + G2; // Offsets for middle corner in (x,y) unskewed coords 100 double y1 = y0 - j1 + G2; 101 double x2 = x0 - 1.0 + 2.0 * G2; // Offsets for last corner in (x,y) unskewed coords 102 double y2 = y0 - 1.0 + 2.0 * G2; 103 // Work out the hashed gradient indices of the three simplex corners 104 const int ii = i & 255; 105 const int jj = j & 255; 106 int gi0 = _permMod12[ii + _perm[jj]]; 107 int gi1 = _permMod12[ii + i1 + _perm[jj + j1]]; 108 int gi2 = _permMod12[ii + 1 + _perm[jj + 1]]; 109 // Calculate the contribution from the three corners 110 double t0 = 0.5 - x0 * x0 - y0 * y0; 111 if (t0 < 0) 112 n0 = 0.0; 113 else { 114 t0 *= t0; 115 n0 = t0 * t0 * dot(_grads[gi0], x0, y0); // (x,y) of _grads used for 2D gradient 116 } 117 double t1 = 0.5 - x1 * x1 - y1 * y1; 118 if (t1 < 0) 119 n1 = 0.0; 120 else { 121 t1 *= t1; 122 n1 = t1 * t1 * dot(_grads[gi1], x1, y1); 123 } 124 double t2 = 0.5 - x2 * x2 - y2 * y2; 125 if (t2 < 0) 126 n2 = 0.0; 127 else { 128 t2 *= t2; 129 n2 = t2 * t2 * dot(_grads[gi2], x2, y2); 130 } 131 // Add contributions from each corner to get the final noise value. 132 // The result is scaled to return values in the interval [-1,1]. 133 return 70.0 * (n0 + n1 + n2); 134 } 135 136 /// 3D simplex noise 137 double noise(double xin, double yin, double zin) { 138 double n0, n1, n2, n3; // Noise contributions from the four corners 139 // Skew the input space to determine which simplex cell we're in 140 const double s = (xin + yin + zin) * F3; // Very nice and simple skew factor for 3D 141 const int i = fastfloor(xin + s); 142 const int j = fastfloor(yin + s); 143 const int k = fastfloor(zin + s); 144 const double t = (i + j + k) * G3; 145 const double X0 = i - t; // Unskew the cell origin back to (x,y,z) space 146 const double Y0 = j - t; 147 const double Z0 = k - t; 148 double x0 = xin - X0; // The x,y,z distances from the cell origin 149 double y0 = yin - Y0; 150 double z0 = zin - Z0; 151 // For the 3D case, the simplex shape is a slightly irregular tetrahedron. 152 // Determine which simplex we are in. 153 int i1, j1, k1; // Offsets for second corner of simplex in (i,j,k) coords 154 int i2, j2, k2; // Offsets for third corner of simplex in (i,j,k) coords 155 if (x0 >= y0) { 156 if (y0 >= z0) { 157 i1 = 1; 158 j1 = 0; 159 k1 = 0; 160 i2 = 1; 161 j2 = 1; 162 k2 = 0; 163 } // X Y Z order 164 else if (x0 >= z0) { 165 i1 = 1; 166 j1 = 0; 167 k1 = 0; 168 i2 = 1; 169 j2 = 0; 170 k2 = 1; 171 } // X Z Y order 172 else { 173 i1 = 0; 174 j1 = 0; 175 k1 = 1; 176 i2 = 1; 177 j2 = 0; 178 k2 = 1; 179 } // Z X Y order 180 } 181 else { // x0<y0 182 if (y0 < z0) { 183 i1 = 0; 184 j1 = 0; 185 k1 = 1; 186 i2 = 0; 187 j2 = 1; 188 k2 = 1; 189 } // Z Y X order 190 else if (x0 < z0) { 191 i1 = 0; 192 j1 = 1; 193 k1 = 0; 194 i2 = 0; 195 j2 = 1; 196 k2 = 1; 197 } // Y Z X order 198 else { 199 i1 = 0; 200 j1 = 1; 201 k1 = 0; 202 i2 = 1; 203 j2 = 1; 204 k2 = 0; 205 } // Y X Z order 206 } 207 // A step of (1,0,0) in (i,j,k) means a step of (1-c,-c,-c) in (x,y,z), 208 // a step of (0,1,0) in (i,j,k) means a step of (-c,1-c,-c) in (x,y,z), and 209 // a step of (0,0,1) in (i,j,k) means a step of (-c,-c,1-c) in (x,y,z), where 210 // c = 1/6. 211 double x1 = x0 - i1 + G3; // Offsets for second corner in (x,y,z) coords 212 double y1 = y0 - j1 + G3; 213 double z1 = z0 - k1 + G3; 214 double x2 = x0 - i2 + 2.0 * G3; // Offsets for third corner in (x,y,z) coords 215 double y2 = y0 - j2 + 2.0 * G3; 216 double z2 = z0 - k2 + 2.0 * G3; 217 double x3 = x0 - 1.0 + 3.0 * G3; // Offsets for last corner in (x,y,z) coords 218 double y3 = y0 - 1.0 + 3.0 * G3; 219 double z3 = z0 - 1.0 + 3.0 * G3; 220 // Work out the hashed gradient indices of the four simplex corners 221 const int ii = i & 255; 222 const int jj = j & 255; 223 const int kk = k & 255; 224 int gi0 = _permMod12[ii + _perm[jj + _perm[kk]]]; 225 int gi1 = _permMod12[ii + i1 + _perm[jj + j1 + _perm[kk + k1]]]; 226 int gi2 = _permMod12[ii + i2 + _perm[jj + j2 + _perm[kk + k2]]]; 227 int gi3 = _permMod12[ii + 1 + _perm[jj + 1 + _perm[kk + 1]]]; 228 // Calculate the contribution from the four corners 229 double t0 = 0.6 - x0 * x0 - y0 * y0 - z0 * z0; 230 if (t0 < 0) 231 n0 = 0.0; 232 else { 233 t0 *= t0; 234 n0 = t0 * t0 * dot(_grads[gi0], x0, y0, z0); 235 } 236 double t1 = 0.6 - x1 * x1 - y1 * y1 - z1 * z1; 237 if (t1 < 0) 238 n1 = 0.0; 239 else { 240 t1 *= t1; 241 n1 = t1 * t1 * dot(_grads[gi1], x1, y1, z1); 242 } 243 double t2 = 0.6 - x2 * x2 - y2 * y2 - z2 * z2; 244 if (t2 < 0) 245 n2 = 0.0; 246 else { 247 t2 *= t2; 248 n2 = t2 * t2 * dot(_grads[gi2], x2, y2, z2); 249 } 250 double t3 = 0.6 - x3 * x3 - y3 * y3 - z3 * z3; 251 if (t3 < 0) 252 n3 = 0.0; 253 else { 254 t3 *= t3; 255 n3 = t3 * t3 * dot(_grads[gi3], x3, y3, z3); 256 } 257 // Add contributions from each corner to get the final noise value. 258 // The result is scaled to stay just inside [-1,1] 259 return 32.0 * (n0 + n1 + n2 + n3); 260 }