diff --git a/demos/test/simplex.cpp b/demos/test/simplex.cpp index 2c5342bd..f9bcae61 100644 --- a/demos/test/simplex.cpp +++ b/demos/test/simplex.cpp @@ -18,15 +18,17 @@ using namespace lol; /* Constants to tweak */ -float const zoom = 0.03f; +ivec2 const size(1280 * 1, 720 * 1); +float const zoom = 0.03f / 1; int const octaves = 1; int main(int argc, char **argv) { UNUSED(argc, argv); + srand(time(nullptr)); + /* Create an image */ - ivec2 const size(1280, 720); Image img(size); array2d &data = img.Lock2D(); @@ -52,6 +54,7 @@ int main(int argc, char **argv) float x = (float)i, y = (float)j; float sum = 0.f, coeff = 0.f; int const max_k = 1 << octaves; + bool b_centre = false, b_sphere1 = false, b_sphere2 = false; for (int k = 1; k < max_k; k *= 2) { @@ -60,36 +63,52 @@ int main(int argc, char **argv) switch (cell) { case 0: - t = s2.Interp(zoom * k * vec2(x, y)); + t = s2.eval(zoom * k * vec2(x, y)); break; case 1: - t = s3.Interp(zoom * k * vec3(x, 1.f, y)); + t = s3.eval(zoom * k * vec3(x, 1.f, y)); break; case 2: - t = s4.Interp(zoom * k * vec4(x, 1.f, y, 2.f)); + t = s4.eval(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)); + t = s6.eval(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)); + t = s8.eval(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)); + //t = s12.eval(zoom * k * vec12(x / 2, -x / 2, y / 2, -y / 2, + // -x / 2, x / 2, -y / 2, y / 2, + // 7.f, 8.f, 9.f, 10.f)); + t = s12.eval(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; } + if (t == -2.f) b_centre = true; + if (t == -3.f) b_sphere1 = true; + if (t == -4.f) b_sphere2 = true; + sum += 1.f / k * t; coeff += 1.f / k; } - 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 (b_centre) + data[i][j] = vec4(1.f, 0.f, 1.f, 1.f); + else if (b_sphere1) + data[i][j] = vec4(0.f, 1.f, 0.f, 1.f); + else if (b_sphere2) + data[i][j] = vec4(0.f, 0.f, 1.f, 1.f); + else + { + 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 @@ -108,7 +127,7 @@ int main(int argc, char **argv) ivec2 coord = ivec2(i / zoom * dx + j / zoom * dy); - vec2 g = s2.GetGradient((i + 0.1f) * dx + (j + 0.1f) * dy); + vec2 g = s2.gradient((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)); diff --git a/src/lol/math/simplex_interpolator.h b/src/lol/math/simplex_interpolator.h index 1ba7bbfc..666c016a 100644 --- a/src/lol/math/simplex_interpolator.h +++ b/src/lol/math/simplex_interpolator.h @@ -50,8 +50,8 @@ public: #endif } - /* Single interpolation */ - inline float Interp(vec_t position) const + /* Evaluate noise at a given point */ + inline float eval(vec_t position) const { // Retrieve the containing hypercube origin and associated decimals vec_t origin; @@ -61,8 +61,9 @@ public: return get_noise(origin, pos); } - /* Only for debug purposes: return the gradient vector */ - inline vec_t GetGradient(vec_t position) const + /* Only for debug purposes: return the gradient vector of the given + * point’s simplex origin. */ + inline vec_t gradient(vec_t position) const { vec_t origin; vec_t pos; @@ -83,7 +84,8 @@ protected: for (int i = 0; i < N; ++i) traversal_order[i] = i; - /* Naïve bubble sort — enough for now */ + /* Naïve bubble sort — enough for now since the general complexity + * of our algorithm is O(N²). Sorting in O(N²) will not hurt more! */ for (int i = 0; i < N; ++i) for (int j = i + 1; j < N; ++j) if (pos[traversal_order[i]] < pos[traversal_order[j]]) @@ -96,19 +98,26 @@ protected: /* “corner” will traverse the simplex along its edges in world * coordinates. */ vec_t world_corner(0.f); - float result = 0.f, sum = 0.f; + float result = 0.f, sum = 0.f, special = 0.f; + UNUSED(sum, special); for (int i = 0; i < N + 1; ++i) { #if 1 - // In “Noise Hardware” (2-17) Perlin uses 0.6 but Gustavson uses - // 0.5 instead, saying “else the noise is not continuous at - // simplex boundaries”. - // And indeed, the distance between any given simplex vertex and - // the opposite hyperplane is 1/sqrt(2), so the contribution of - // that vertex should never be > 0 for points further than - // 1/sqrt(2). Hence 0.5 - d². - float d = 0.5f - sqlength(world_pos - world_corner); + /* In “Noise Hardware” (2-17) Perlin uses 0.6 - d². + * + * In an errata to “Simplex noise demystified”, Gustavson uses + * 0.5 - d² instead, saying “else the noise is not continuous at + * simplex boundaries”. + * And indeed, the distance between any given simplex vertex and + * the opposite hyperplane is 1/sqrt(2), so the contribution of + * that vertex should never be > 0 for points further than + * 1/sqrt(2). Hence 0.5 - d². + * + * I prefer to use 1 - 2d² and compensate for the d⁴ below by + * dividing the final result by 2⁴, because manipulating values + * between 0 and 1 is more convenient. */ + float d = 1.0f - 2.f * sqlength(world_pos - world_corner); #else // DEBUG: this is the linear contribution of each vertex // in the skewed simplex. Unfortunately it creates artifacts. @@ -116,6 +125,16 @@ protected: - ((i == N) ? 0.f : pos[traversal_order[i]]); #endif +#if 0 + // DEBUG: identify simplex features: + // -4.f: centre (-2.f), + // -3.f: r=0.38 sphere of influence (contribution = 1/4) + // -2.f: r=0.52 sphere of influence (contribution = 1/24) + if (d > 0.99f) special = min(special, -4.f); + if (d > 0.7f && d < 0.72f) special = min(special, -3.f); + if (d > 0.44f && d < 0.46f) special = min(special, -2.f); +#endif + if (d > 0) { // Perlin uses 8d⁴ whereas Gustavson uses d⁴ and a final @@ -140,16 +159,30 @@ protected: } } - // FIXME: Gustavson uses the value 70 for dimension 2, 32 for - // dimension 3, and 27 for dimension 4, and uses non-unit gradients - // of length sqrt(2), sqrt(2) and sqrt(3). Find out where this comes - // from and maybe find a more generic formula. - float const k = N == 2 ? (70.f * sqrt(2.f)) // approx. 99 - : N == 3 ? (70.f * sqrt(2.f)) // approx. 45 - : N == 4 ? (70.f * sqrt(3.f)) // approx. 47 - : 50.f; - //return k * result / sum; - return k * result; +#if 0 + if (special < 0.f) + return special; +#endif + + return get_scale() * result; + } + + static inline float get_scale() + { + /* FIXME: Gustavson uses the value 70 for dimension 2, 32 for + * dimension 3, and 27 for dimension 4, and uses non-unit gradients + * of length sqrt(2), sqrt(2) and sqrt(3), saying “The result is + * scaled to stay just inside [-1,1]” which honestly is just not + * true. + * Experiments show that the scaling factor is remarkably close + * to 6.7958 for all high dimensions (measured up to 12). */ + return N == 2 ? 6.2003f + : N == 3 ? 6.7297f + : N == 4 ? 6.7861f + : N == 5 ? 6.7950f + : N == 6 ? 6.7958f + : N == 7 ? 6.7958f + : /* 7+ */ 6.7958f; } inline vec_t get_gradient(vec_t origin) const @@ -272,12 +305,13 @@ private: for (int i = 0; i < N + 1; ++i) vertices[i] = unskew(vertices[i]); + /* Output information for each vertex */ for (int i = 0; i < N + 1; ++i) { printf(" - vertex %d\n", i); -#if 0 /* Coordinates for debugging purposes */ +#if 0 printf(" · ["); for (int k = 0; k < N; ++k) printf(" %f", vertices[i][k]); @@ -285,6 +319,7 @@ private: #endif /* Analyze edge lengths from that vertex */ +#if 0 float minlength = 1.0f; float maxlength = 0.0f; for (int j = 0; j < N + 1; ++j) @@ -298,12 +333,13 @@ private: } printf(" · edge lengths between %f and %f\n", minlength, maxlength); +#endif -#if 0 /* Experimental calculation of the distance to the opposite * hyperplane, by picking random points. Works reasonably * well up to dimension 6. After that, we’d need something * better such as gradient walk. */ +#if 0 float mindist = 1.0f; for (int run = 0; run < 10000000; ++run) { @@ -322,7 +358,6 @@ private: printf(" · approx. dist. to opposite hyperplane: %f\n", mindist); #endif -#if 0 /* Find a normal vector to the opposite hyperplane. First, pick * any point p(i0) on the hyperplane. We just need i0 != i. Then, * build a matrix where each row is p(i0)p(j) for all j != i0. @@ -330,6 +365,7 @@ private: * full of zeroes except at position i. So we build a vector * full of zeroes except at position i, and multiply it by the * matrix inverse. */ +#if 0 int i0 = (i == 0) ? 1 : 0; mat_t m; for (int j = 0; j < N; ++j) @@ -349,6 +385,68 @@ private: printf(" · distance to opposite hyperplane: %f\n", dist); #endif } + + /* Compute some statistics about the noise. TODO: histogram */ +#if 0 + vec_t input(0.f); + for (int run = 0; run < 1000000; ++run) + { + float t = eval(input); + + input[run % N] = rand(1000.f); + } +#endif + + /* Try to find max noise value by climbing gradient */ + float minval = 0.f, maxval = 0.f; + array> deltas; + for (int i = 0; i < N; ++i) + { + vec_t v(0.f); + v[i] = 1.f; + deltas << v << -v; + } + for (int run = 0; run < 1000; ++run) + { + /* Pick a random vector */ + vec_t v; + for (int i = 0; i < N; ++i) + v[i] = rand(-100.f, 100.f); + float t = eval(v); + float e = 0.1f; + /* Climb up gradient in all dimensions */ + while (e > 1e-6f) + { + int best_delta = -1; + float best_t2 = t; + for (int i = 0; i < deltas.Count(); ++i) + { + float t2 = eval(v + e * deltas[i]); + if (abs(t2) > abs(best_t2)) + { + best_delta = i; + best_t2 = t2; + } + } + if (best_delta == -1) + e *= 0.5f; + else + { + v += e * deltas[best_delta]; + t = best_t2; + } + } + minval = min(t, minval); + maxval = max(t, maxval); + } + printf(" - noise value min/max: %f %f\n", minval, maxval); + float newscale = 1.f / max(-minval, maxval); + if (newscale < 1.f) + printf(" - could replace scale %f with %f\n", + get_scale(), newscale * get_scale()); + else + printf(" - scale looks OK\n"); + printf("\n"); } diff --git a/src/t/math/simplex_interpolator.cpp b/src/t/math/simplex_interpolator.cpp index 50d3a629..e828ebc6 100644 --- a/src/t/math/simplex_interpolator.cpp +++ b/src/t/math/simplex_interpolator.cpp @@ -19,354 +19,11 @@ namespace lol { -template -class test_interpolator : public simplex_interpolator -{ -public: - - test_interpolator() : - simplex_interpolator() - { - } - - void DumpMatrix() - { - std::cout << std::endl; - for (int i = 0; i < N; ++i) - { - for (int j = 0; j < N; ++j) - { - std::cout << this->m_base[i][j] << ", "; - } - std::cout << ";"; - } - - std::cout << std::endl; - - for (int i = 0; i < N; ++i) - { - for (int j = 0; j < N; ++j) - { - std::cout << this->m_base_inverse[i][j] << ", "; - } - std::cout << ";"; - } - - std::cout << std::endl; - } - - void DumpMatrix(float result[N][N]) - { - for (int i = 0; i < N; ++i) - { - for (int j = 0; j < N; ++j) - { - result[i][j] = this->m_base[i][j]; - } - } - } - - void DumpCheckInverse(float result[N][N]) - { - for (int i = 0; i < N; ++i) - for (int j = 0; j < N; ++j) - result[i][j] = 0; - - for (int i = 0; i < N; ++i) - { - for (int j = 0; j < N; ++j) - { - for (int k = 0; k < N; ++k) - { - result[i][j] += this->m_base[i][k] * this->m_base_inverse[k][j]; - } - } - } - } - - vec_t GetIndexOrder(vec_t const & decimal_point) - { - return simplex_interpolator::GetIndexOrder(decimal_point); - } -}; - -lolunit_declare_fixture(SimplexInterpolatorTest) +lolunit_declare_fixture(SimplexevalolatorTest) { void SetUp() {} void TearDown() {} - - template - void check_base_matrix() - { - test_interpolator s; - - float result[N][N]; - s.DumpCheckInverse(result); - - // Check base matrix inverse - for (int i = 0; i < N; ++i) - { - for (int j = 0; j < N; ++j) - { - if (i == j) - lolunit_assert_doubles_equal(1, result[i][j], 1e-6); - else - lolunit_assert_doubles_equal(0, result[i][j], 1e-6); - } - } - - s.DumpMatrix(result); - - // Check base vectors’ norms - for (int i = 0; i < N; ++i) - { - float norm = 0; - - for (int j = 0; j < N; ++j) - { - norm += result[i][j] * result[i][j]; - } - - lolunit_assert_doubles_equal(1, norm, 1e-6); - } - - // Check result of sum of base vectors => must have norm = 1 - float vec_result[N]; - for (int i = 0; i < N; ++i) - { - vec_result[i] = 0; - for (int j = 0; j < N; ++j) - { - vec_result[i] += result[i][j]; - } - } - - float norm = 0; - for (int i = 0 ; i < N ; ++i) - { - norm += vec_result[i] * vec_result[i]; - } - lolunit_assert_doubles_equal(1, norm, 1e-6); - } - - lolunit_declare_test(CoordinateMatrixTest) - { - check_base_matrix<1>(); - check_base_matrix<2>(); - check_base_matrix<3>(); - check_base_matrix<4>(); - check_base_matrix<5>(); - check_base_matrix<6>(); - check_base_matrix<7>(); - check_base_matrix<8>(); - check_base_matrix<9>(); - check_base_matrix<10>(); - } - - template - void check_index_ordering() - { - static int gen = 12345678; - test_interpolator s; - - vec_t vec_ref; - for (int i = 0 ; i < N ; ++i) - vec_ref[i] = (gen = gen * gen + gen); - - vec_t result = s.GetIndexOrder(vec_ref); - - for (int i = 1 ; i < N ; ++i) - lolunit_assert(vec_ref[result[i]] < vec_ref[result[i-1]]); - } - - lolunit_declare_test(IndexOrderTest) - { - for (int i = 1 ; i < 10 ; ++i) - check_index_ordering<1>(); - for (int i = 1 ; i < 10 ; ++i) - check_index_ordering<2>(); - for (int i = 1 ; i < 10 ; ++i) - check_index_ordering<3>(); - for (int i = 1 ; i < 10 ; ++i) - check_index_ordering<4>(); - for (int i = 1 ; i < 10 ; ++i) - check_index_ordering<5>(); - for (int i = 1 ; i < 10 ; ++i) - check_index_ordering<6>(); - for (int i = 1 ; i < 10 ; ++i) - check_index_ordering<7>(); - for (int i = 1 ; i < 10 ; ++i) - check_index_ordering<8>(); - for (int i = 1 ; i < 10 ; ++i) - check_index_ordering<9>(); - for (int i = 1 ; i < 10 ; ++i) - check_index_ordering<10>(); - } - - lolunit_declare_test(CheckSampleOrder2) - { - static int gen = 12345678; - - arraynd<2, vec_t > gradients(vec_t({10, 10})); - for (int i = 0 ; i < 10 ; ++i) - { - for (int j = 0 ; j < 10 ; ++j) - { - float x1 = (gen = gen * gen + gen) % 255; - float x2 = (gen = gen * gen + gen) % 255; - - vec_t v = {x1, x2}; - - gradients[i][j] = v; - } - } - - simplex_interpolator<2> s; - s.SetGradients(gradients); - - std::cout << std::endl; - for (int i = -64 ; i < 64 ; ++i) - { - for (int j = -64 ; j < 64 ; ++j) - { - std::cout << s.Interp(vec_t({i * 0.1f, j * 0.1f})) << ", "; - } - std::cout << std::endl; - } - std::cout << std::endl; - } - -#if 0 - lolunit_declare_test(FloatGridPoints2D1x1) - { - simplex_interpolator<2> s({{1.f}}); - float val; - - val = s.Interp(GridPoint2D(-1, -1)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(0, -1)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(1, -1)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(-1, 0)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(0, 0)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(1, 0)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(-1, 1)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(0, 1)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(1, 1)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - } - - lolunit_declare_test(VectorGridPoints2D1x1) - { - simplex_interpolator<2, vec3> s({{vec3(1.f, 2.f, 3.f)}}); - - vec3 val = s.Interp(GridPoint2D(-1, 0)); - lolunit_assert_doubles_equal(1.f, val.x, 1e-5f); - lolunit_assert_doubles_equal(2.f, val.y, 1e-5f); - lolunit_assert_doubles_equal(3.f, val.z, 1e-5f); - } - - lolunit_declare_test(GridPoints2D2x2) - { - simplex_interpolator<2> s({{1.f, 1.f}, {1.f, 2.f}}); - float val; - - val = s.Interp(GridPoint2D(0, 0)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(1, 0)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(0, 1)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(1, 1)); - lolunit_assert_doubles_equal(2.f, val, 1e-5f); - } - - lolunit_declare_test(MoreGridPoints2D2x2) - { - simplex_interpolator<2> s({{1.f, 1.f}, {1.f, 2.f}}); - float val; - - val = s.Interp(GridPoint2D(-1, -1)); - lolunit_assert_doubles_equal(2.f, val, 1e-5f); - val = s.Interp(GridPoint2D(0, -1)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(1, -1)); - lolunit_assert_doubles_equal(2.f, val, 1e-5f); - val = s.Interp(GridPoint2D(2, -1)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(2, 0)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(2, 1)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(2, 2)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(1, 2)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(0, 2)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(-1, 2)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(-1, 1)); - lolunit_assert_doubles_equal(2.f, val, 1e-5f); - val = s.Interp(GridPoint2D(-1, 0)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - } - - lolunit_declare_test(GridPoints2D3x3) - { - simplex_interpolator<2> s({{1, 1, 2}, {1, 2, 2}, {2, 2, 2}}); - float val; - - val = s.Interp(GridPoint2D(0, 0)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(1, 0)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(0, 1)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - val = s.Interp(GridPoint2D(1, 1)); - lolunit_assert_doubles_equal(2.f, val, 1e-5f); - } - - lolunit_declare_test(OtherPoints2D3x3) - { - simplex_interpolator<2> s({{1, 1, 2}, {1, 2, 2}, {2, 2, 2}}); - float val; - - /* 1/1 triangle edge */ - val = s.Interp(vec2(1.f, lol::sin(F_PI / 3))); - lolunit_assert_doubles_equal(1.5f, val, 1e-5f); - - /* 1/1/1 triangle interior */ - val = s.Interp(vec2(0.5f, 0.5f)); - lolunit_assert_doubles_equal(1.f, val, 1e-5f); - - /* 1/1/2 triangle interior */ - val = s.Interp(vec2(0.f, 0.5f)); - lolunit_assert(val < 2.f); - lolunit_assert(val > 1.f); - - /* Another 1/1/2 triangle interior */ - val = s.Interp(vec2(1.f, 0.5f)); - lolunit_assert(val < 2.f); - lolunit_assert(val > 1.f); - } - -private: - vec2 GridPoint2D(int i, int j) - { - static vec2 vi(1.f, 0.f); - static vec2 vj(lol::cos(F_PI / 3), lol::sin(F_PI / 3)); - - return (float)i * vi + (float)j * vj; - } -#endif }; }