diff --git a/demos/test/simplex.cpp b/demos/test/simplex.cpp index 609a911f..2c5342bd 100644 --- a/demos/test/simplex.cpp +++ b/demos/test/simplex.cpp @@ -17,44 +17,108 @@ using namespace lol; +/* Constants to tweak */ +float const zoom = 0.03f; +int const octaves = 1; + int main(int argc, char **argv) { UNUSED(argc, argv); - float const zoom = 0.03f; - int const octaves = 10; - /* Create an image */ ivec2 const size(1280, 720); Image img(size); array2d &data = img.Lock2D(); - /* Fill image with simplex noise */ + /* Declare plenty of allocators */ simplex_interpolator<2> s2; simplex_interpolator<3> s3; + simplex_interpolator<4> s4; + simplex_interpolator<5> s5; + simplex_interpolator<6> s6; simplex_interpolator<7> s7; + simplex_interpolator<8> s8; + simplex_interpolator<9> s9; + simplex_interpolator<10> s10; + simplex_interpolator<11> s11; + simplex_interpolator<12> s12; + /* Fill image with simplex noise */ for (int j = 0; j < size.y; ++j) for (int i = 0; i < size.x; ++i) { + int cell = j / (size.y / 2) * 3 + i / (size.x / 3); + float x = (float)i, y = (float)j; - float sum = 0.f; - int maxoctave = (j < size.y / 2) ? 1 : (1 << octaves); + float sum = 0.f, coeff = 0.f; + int const max_k = 1 << octaves; - for (int k = 1; k <= maxoctave; k *= 2) + for (int k = 1; k < max_k; k *= 2) { - if (i < size.x / 3) - sum += 0.5f / k * s2.Interp(zoom * k * vec2(x, y)); - else if (i < size.x * 2 / 3) - sum += 0.5f / k * s3.Interp(zoom * k * vec3(x, y, 0.f)); - else - sum += 0.5f / k * s7.Interp(zoom * k * vec7(x, 0.f, 0.f, 0.f, y, 0.f, 0.f)); + float t = 0.f; + + switch (cell) + { + case 0: + t = s2.Interp(zoom * k * vec2(x, y)); + break; + case 1: + t = s3.Interp(zoom * k * vec3(x, 1.f, y)); + break; + case 2: + t = s4.Interp(zoom * k * vec4(x, 1.f, y, 2.f)); + break; + case 3: + t = s6.Interp(zoom * k * vec6(x, 1.f, 2.f, y, 3.f, 4.f)); + break; + case 4: + t = s8.Interp(zoom * k * vec8(x, 1.f, 2.f, 3.f, y, 4.f, + 5.f, 6.f)); + break; + case 5: + t = s12.Interp(zoom * k * vec12(x, 1.f, 2.f, 3.f, 4.f, 5.f, + y, 6.f, 7.f, 8.f, 9.f, 10.f)); + break; + default: + break; + } + + sum += 1.f / k * t; + coeff += 1.f / k; } - float c = saturate(0.5f + sum); + float c = saturate(0.5f + 0.5f * sum / coeff); data[i][j] = vec4(c, c, c, 1.f); + //data[i][j] = Color::HSVToRGB(vec4(c, 1.0f, 0.5f, 1.f)); } +#if 0 + /* Mark simplex vertices */ + vec2 diagonal = normalize(vec2(1.f)); + vec2 dx = mat2::rotate(60.f) * diagonal; + vec2 dy = mat2::rotate(-60.f) * diagonal; + for (int j = -100; j < 100; ++j) + for (int i = -100; i < 100; ++i) + { + auto putpixel = [&](ivec2 p, vec4 c = vec4(1.f, 0.f, 1.f, 1.f)) + { + if (p.x >= 0 && p.x < size.x - 1 && p.y >= 0 && p.y < size.y - 1) + data[p.x][p.y] = c; + }; + + ivec2 coord = ivec2(i / zoom * dx + j / zoom * dy); + + vec2 g = s2.GetGradient((i + 0.1f) * dx + (j + 0.1f) * dy); + for (int t = 0; t < 40; ++t) + putpixel(coord + (ivec2)(g * (t / 2.f)), vec4(0.f, 1.f, 0.f, 1.f)); + + putpixel(coord); + putpixel(coord + ivec2(1, 0)); + putpixel(coord + ivec2(0, 1)); + putpixel(coord + ivec2(1, 1)); + } +#endif + /* Save image */ img.Unlock2D(data); img.Save("simplex.png"); diff --git a/src/lol/math/simplex_interpolator.h b/src/lol/math/simplex_interpolator.h index c13aed3c..84dfc068 100644 --- a/src/lol/math/simplex_interpolator.h +++ b/src/lol/math/simplex_interpolator.h @@ -17,6 +17,27 @@ namespace lol { +/* + * Simplex noise in dimension N + * ---------------------------- + * + * The N-dimensional regular hypercube can be split into N! simplices that + * all have the main diagonal as a shared edge. + * - number of simplices: N! + * - number of vertices per simplex: N+1 + * - number of edges: N(N+1)/2 + * - minimum edge length: 1 (hypercube edges, e.g. [1,0,0,…,0]) + * - maximum edge length: sqrt(N) (hypercube diagonal, i.e. [1,1,1,…,1]) + * + * We skew the simplicial grid along the main diagonal by a factor f = + * sqrt(N+1), which means the diagonal of the initial parallelotope has + * length sqrt(N/(N+1)). The edges of that parallelotope have length + * sqrt(N/(N+1)), too. A formula for the maximum edge length was found + * empirically: + * - minimum edge length: sqrt(N/(N+1)) (parallelotope edges and diagonal) + * - maximum edge length: sqrt(floor((N+1)²/4)/(N+1)) + */ + template class simplex_interpolator { @@ -24,165 +45,218 @@ public: simplex_interpolator(int seed = 0) : m_seed(seed) { - /* Matrix coordinate transformation to skew simplex referential is done - * by inverting the base matrix M which is written as follows: - * - * M = | a b b b … | M^(-1) = | c d d d … | - * | b a b b … | | d c d d … | - * | b b a b … | | d d c d … | - * | … | | … | - * - * where a, b, c, d are computed below ↴ - */ - - float b = (1.f - lol::sqrt(N + 1.f)) / lol::sqrt((float)N * N * N); - float a = b + lol::sqrt((N + 1.f) / N); - - // determinant of matrix M - float det = a * (a + (N - 2) * b) - (b * b) * (N - 1); - - float c = (a + (N - 2) * b) / det; - float d = -b / det; - +#if 0 + // Print some debug information + printf("Simplex Noise of Dimension %d\n", N); + + long long int n = 1; for (int i = 1; i <= N; ++i) n *= i; + printf(" - each hypercube cell has %lld simplices " + "with %d vertices and %d edges\n", n, N + 1, N * (N + 1) / 2); + + vec_t diagonal(1.f); + printf(" - length of hypercube edges is 1, diagonal is %f\n", + length(diagonal)); + printf(" - length of parallelotope edges and diagonal is %f\n", + length(unskew(diagonal))); + + vec_t vertices[N + 1]; + vertices[0] = vec_t(0.f); for (int i = 0; i < N; ++i) { - for (int j = 0; j < N; ++j) + vertices[i + 1] = vertices[i]; + vertices[i + 1][i] += 1.f; + } + float minlength = 10.0f; + float maxlength = 0.0f; + for (int i = 0; i < N + 1; ++i) + for (int j = i + 1; j < N + 1; ++j) { - m_base[i][j] = (i == j ? a : b); - m_base_inverse[i][j] = (i == j ? c : d); + float l = length(unskew(vertices[i] - vertices[j])); + minlength = min(minlength, l); + maxlength = max(maxlength, l); } - } + printf(" - edge lengths between %f and %f\n", minlength, maxlength); + printf(" - predicted: %f and %f\n", sqrt((float)N/(N+1)), + sqrt((N+1)*(N+1)/4/(float)(N+1))); + printf("\n"); +#endif } - // Single interpolation + /* Single interpolation */ inline float Interp(vec_t position) const { - // Computing position in simplex referential - vec_t simplex_position = m_base_inverse * position; - - // Retrieving the closest floor point and decimals + // Retrieve the containing hypercube origin and associated decimals vec_t origin; vec_t pos; + get_origin(skew(position), origin, pos); - this->ExtractFloorDecimal(simplex_position, origin, pos); - - vec_t index_order = this->GetIndexOrder(pos); - - // Resetting decimal point in regular orthonormal coordinates - pos = m_base * pos; + return get_noise(origin, pos); + } - return this->LastInterp(origin, pos, index_order); + /* Only for debug purposes: return the gradient vector */ + inline vec_t GetGradient(vec_t position) const + { + vec_t origin; + vec_t pos; + get_origin(skew(position), origin, pos); + return get_gradient(origin); } protected: - inline float LastInterp(vec_t origin, - vec_t const & pos, - vec_t const & index_order) const + inline float get_noise(vec_t origin, + vec_t const & pos) const { - /* “corner” will traverse the simplex along its edges in - * orthonormal coordinates. */ - vec_t corner(0); - float result = 0; + /* For a given position [0…1]^N inside a regular N-hypercube, find + * the N-simplex which contains that position, and return a path + * along the hypercube edges from (0,0,…,0) to (1,1,…,1) which + * uniquely describes that simplex. */ + vec_t traversal_order; + for (int i = 0; i < N; ++i) + traversal_order[i] = i; + + /* Naïve bubble sort — enough for now */ + for (int i = 0; i < N; ++i) + for (int j = i + 1; j < N; ++j) + if (pos[traversal_order[i]] < pos[traversal_order[j]]) + std::swap(traversal_order[i], traversal_order[j]); + + + /* Get the position in world coordinates, too */ + vec_t world_pos = unskew(pos); + + /* “corner” will traverse the simplex along its edges in world + * coordinates. */ + vec_t world_corner(0.f); + float result = 0.f, coeff = 0.f; for (int i = 0; i < N + 1; ++i) { - vec_t const delta = pos - corner; - vec_t const &gradient = GetGradient(origin); - - /* We use 4/3 because the maximum radius of influence for - * a given corner is sqrt(3/4). FIXME: check whether this - * holds for higher dimensions. */ - float dist = 1.0f - 4.f / 3.f * sqlength(delta); - if (dist > 0) +#if 1 + // FIXME: In “Noise Hardware” (2-17) Perlin uses 0.6 but Gustavson + // makes an exception for dimension 2 and uses 0.5 instead. + // I think m should actually increase with the dimension count. + float const m = N == 2 ? 0.5f : 0.6f; + float d = m - sqlength(world_pos - world_corner); +#else + // DEBUG: this is the linear contribution of each vertex + // in the skewed simplex. Unfortunately it creates artifacts. + float d = ((i == 0) ? 1.f : pos[traversal_order[i - 1]]) + - ((i == N) ? 0.f : pos[traversal_order[i]]); +#endif + + if (d > 0) { - dist *= dist; - result += dist * dist * dot(gradient, delta); + vec_t const &gradient = get_gradient(origin); + + // Perlin uses 8d⁴ whereas Gustavson uses d⁴ and a final + // multiplication factor at the end. + //d = 8.f * d * d * d * d; + d = d * d * d * d; + + //d = (3.f - 2.f * d) * d * d; + //d = ((6 * d - 15) * d + 10) * d * d * d; + + result += d * dot(gradient, world_pos - world_corner); + coeff += d; } if (i < N) { - corner += m_base[index_order[i]]; - origin[index_order[i]] += 1; + vec_t v(0.f); + v[traversal_order[i]] = 1.f; + world_corner += unskew(v); + origin[traversal_order[i]] += 1; } } - // FIXME: another paper uses the value 70 for dimension 2, 32 for - // dimension 3, and 27 for dimension 4; find where this comes from. - return result * 70.f / 16.f; - } - - /* For a given position [0…1]^n inside a square/cube/hypercube etc., - * find the simplex which contains that position, and return the path - * from (0,0,…,0) to (1,1,…,1) that describes that simplex. */ - inline vec_t GetIndexOrder(vec_t const & pos) const - { - vec_t result; - for (int i = 0; i < N; ++i) - result[i] = i; - - /* Naïve bubble sort — enough for now */ - for (int i = 0; i < N; ++i) - for (int j = i + 1; j < N; ++j) - if (pos[result[i]] < pos[result[j]]) - std::swap(result[i], result[j]); - - return result; + // FIXME: Gustavson uses the value 70 for dimension 2, 32 for + // dimension 3, and 27 for dimension 4; find where this comes from + // and maybe find a reasonable formula. + //return 70.f * result / coeff; + float const k = N == 2 ? 70 : N == 3 ? 32 : N == 4 ? 27 : 20; + return k * result; } - inline vec_t GetGradient(vec_t origin) const + inline vec_t get_gradient(vec_t origin) const { /* Quick shuffle table: * strings /dev/urandom | grep . -nm256 | sort -k2 -t: | sed 's|:.*|,|' * Then just replace “256” with “0”. */ static int const shuffle[256] = { - 111, 14, 180, 186, 221, 114, 17, 79, 66, 46, 11, 81, 246, 200, 141, - 172, 85, 244, 112, 92, 34, 106, 218, 205, 236, 7, 121, 115, 109, - 131, 10, 96, 188, 148, 219, 107, 94, 182, 235, 163, 143, 213, 248, - 202, 52, 154, 37, 241, 53, 129, 25, 60, 242, 38, 171, 63, 203, 255, - 193, 6, 42, 209, 28, 176, 210, 159, 54, 144, 3, 71, 89, 116, 12, - 237, 67, 216, 252, 178, 174, 164, 98, 234, 32, 26, 175, 24, 130, - 128, 113, 99, 212, 62, 152, 75, 185, 73, 93, 31, 30, 151, 122, 173, - 139, 91, 136, 162, 194, 208, 56, 101, 68, 69, 211, 44, 97, 55, 83, - 33, 50, 119, 156, 149, 41, 157, 253, 247, 161, 47, 230, 166, 225, - 204, 224, 13, 110, 123, 142, 64, 65, 155, 215, 120, 197, 140, 58, - 77, 214, 126, 195, 179, 220, 232, 125, 147, 8, 39, 187, 27, 217, - 100, 134, 199, 88, 206, 231, 250, 74, 2, 135, 9, 245, 118, 21, 243, - 82, 183, 238, 87, 158, 61, 4, 177, 146, 153, 117, 249, 254, 233, - 90, 222, 207, 48, 15, 18, 20, 239, 133, 0, 165, 138, 127, 169, 72, - 1, 201, 145, 191, 192, 16, 49, 19, 95, 226, 228, 84, 181, 251, 36, - 150, 22, 43, 70, 45, 105, 5, 189, 160, 196, 40, 59, 57, 190, 80, - 104, 167, 78, 124, 103, 240, 184, 170, 137, 29, 23, 223, 108, 102, - 86, 198, 227, 35, 229, 76, 168, 132, 51, + 111, 14, 180, 186, 221, 114, 219, 79, 66, 46, 152, 81, 246, 200, + 141, 172, 85, 244, 112, 92, 34, 106, 218, 205, 236, 7, 121, 115, + 109, 131, 10, 96, 188, 148, 17, 107, 94, 182, 235, 163, 143, 63, + 248, 202, 52, 154, 37, 241, 53, 129, 25, 159, 242, 38, 171, 213, + 6, 203, 255, 193, 42, 209, 28, 176, 210, 60, 54, 144, 3, 71, 89, + 116, 12, 237, 67, 216, 252, 178, 174, 164, 98, 234, 32, 26, 175, + 24, 130, 128, 113, 99, 212, 62, 11, 75, 185, 73, 93, 31, 30, 44, + 122, 173, 139, 91, 136, 162, 194, 41, 56, 101, 68, 69, 211, 151, + 97, 55, 83, 33, 50, 119, 156, 149, 208, 157, 253, 247, 161, 133, + 230, 166, 225, 204, 224, 13, 110, 123, 142, 64, 65, 155, 215, 9, + 197, 140, 58, 77, 214, 126, 195, 179, 220, 232, 125, 147, 8, 39, + 187, 27, 217, 100, 134, 199, 88, 206, 231, 250, 74, 2, 135, 120, + 21, 245, 118, 243, 82, 183, 238, 150, 158, 61, 4, 177, 146, 153, + 117, 249, 254, 233, 90, 222, 207, 48, 15, 18, 20, 16, 47, 0, 51, + 165, 138, 127, 169, 72, 1, 201, 145, 191, 192, 239, 49, 19, 160, + 226, 228, 84, 181, 251, 36, 87, 22, 43, 70, 45, 105, 5, 189, 95, + 40, 196, 59, 57, 190, 80, 104, 167, 78, 124, 103, 240, 184, 170, + 137, 29, 23, 223, 108, 102, 86, 198, 227, 35, 229, 76, 168, 132, }; - /* 16 random vectors; this should be enough for small dimensions */ - auto v = [&]() - { - vec_t ret; - for (int i = 0; i < N; ++i) - ret[i] = rand(-1.f, 1.f); - return normalize(ret); - }; + /* Generate 2^(N+2) random vectors, but at least 2^5 (32) and not + * more than 2^20 (~ 1 million). */ + int const gradient_count = 1 << min(max(N + 2, 5), 20); - static vec_t const gradients[16] = + static auto build_gradients = [&]() { - v(), v(), v(), v(), v(), v(), v(), v(), - v(), v(), v(), v(), v(), v(), v(), v(), + array> ret; + for (int k = 0; k < gradient_count; ++k) + { + vec_t v; + for (int i = 0; i < N; ++i) + v[i] = rand(-1.f, 1.f); + ret << normalize(v); + } + return ret; }; + static array> const gradients = build_gradients(); + int idx = m_seed; for (int i = 0; i < N; ++i) - idx = shuffle[(idx + origin[i]) & 255]; + idx ^= shuffle[(idx + origin[i]) & 255]; + + idx &= (gradient_count - 1); +#if 0 + // DEBUG: only output a few gradients + if (idx > 2) + return vec_t(0); +#endif + return gradients[idx]; + } - return gradients[(idx ^ (idx >> 4)) & 15]; + static inline vec_t skew(vec_t const &v) + { + /* Quoting Perlin in “Hardware Noise” (2-18): + * The “skew factor” f should be set to f = sqrt(N+1), so that + * the point (1,1,...1) is transformed to the point (f,f,...f). */ + float const sum = dot(v, vec_t(1)); + float const f = sqrt(1.f + N); + return v + vec_t(sum * (f - 1) / N); + } + + static inline vec_t unskew(vec_t const &v) + { + float const sum = dot(v, vec_t(1)); + float const f = sqrt(1.f + N); + return v + vec_t(sum * (1 / f - 1) / N); } /* For a given world position, extract grid coordinates (origin) and * the corresponding delta position (pos). */ - inline void ExtractFloorDecimal(vec_t const & simplex_position, - vec_t & origin, - vec_t & pos) const + inline void get_origin(vec_t const & simplex_position, + vec_t & origin, vec_t & pos) const { // Finding floor point index for (int i = 0; i < N; ++i) @@ -192,7 +266,7 @@ protected: pos = simplex_position - (vec_t)origin; } - mat_t m_base, m_base_inverse; + /* A user-provided random seed. Defaults to zero. */ int m_seed; };