diff --git a/demos/test/simplex.cpp b/demos/test/simplex.cpp index e71b2cbb..6b5a2777 100644 --- a/demos/test/simplex.cpp +++ b/demos/test/simplex.cpp @@ -21,51 +21,37 @@ int main(int argc, char **argv) { UNUSED(argc, argv); - int const period = 10; - float const zoom = 0.02f; - - /* Create some gradients */ - array2d gradients(vec_t({period, period})); - for (int i = 0 ; i < period ; ++i) - for (int j = 0 ; j < period ; ++j) - gradients[i][j] = normalize(vec2(rand(-1.f, 1.f), rand(-1.f, 1.f))); - simplex_interpolator<2> s; - s.SetGradients(gradients); + float const zoom = 0.03f; + int const octaves = 10; /* Create an image */ - ivec2 const size(800, 600); + ivec2 const size(1280, 720); Image img(size); array2d &data = img.Lock2D(); /* Fill image with simplex noise */ + simplex_interpolator<2> s2; + simplex_interpolator<3> s3; + simplex_interpolator<4> s4; + for (int j = 0; j < size.y; ++j) for (int i = 0; i < size.x; ++i) { - float p = 0.5f * s.Interp(zoom * vec2(i, j)); -#if 1 - for (int k = 2; k < 128; k *= 2) - p += 0.5f / k * s.Interp(zoom * k * vec2(i, j)); - //p += 0.5f / k * lol::abs(s.Interp(zoom * k * vec2(i, j))); -#endif - data[i][j] = vec4(saturate(0.5f + p), 0.f, 0.f, 1.f); - } + float sum = 0.f; + int maxoctave = (j < size.y / 2) ? 1 : (1 << octaves); - /* 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) - { - ivec2 coord = ivec2(i / zoom * dx + j / zoom * dy); - if (coord.x >= 0 && coord.x < size.x - 1 - && coord.y >= 0 && coord.y < size.y - 1) + for (int k = 1; k <= maxoctave; k *= 2) { - data[coord.x][coord.y] = vec4(1.f, 0.f, 1.f, 1.f); - data[coord.x + 1][coord.y] = vec4(1.f, 0.f, 1.f, 1.f); - data[coord.x][coord.y + 1] = vec4(1.f, 0.f, 1.f, 1.f); - data[coord.x + 1][coord.y + 1] = vec4(1.f, 0.f, 1.f, 1.f); + if (i < size.x / 3) + sum += 0.5f / k * s2.Interp(zoom * k * vec2(i, j)); + else if (i < size.x * 2 / 3) + sum += 0.5f / k * s3.Interp(zoom * k * vec3(i, j, 0.0f)); + else + sum += 0.5f / k * s4.Interp(zoom * k * vec4(i, j, 0.0f, 0.0f)); } + + float c = saturate(0.5f + sum); + data[i][j] = vec4(c, c, c, 1.f); } /* Save image */ diff --git a/src/lol/math/simplex_interpolator.h b/src/lol/math/simplex_interpolator.h index 7160641f..c13aed3c 100644 --- a/src/lol/math/simplex_interpolator.h +++ b/src/lol/math/simplex_interpolator.h @@ -12,6 +12,8 @@ #pragma once +#include + namespace lol { @@ -19,21 +21,37 @@ template class simplex_interpolator { public: - - simplex_interpolator() : - m_gradients() + simplex_interpolator(int seed = 0) + : m_seed(seed) { - this->InitBase(); - } + /* 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 ↴ + */ - inline void SetGradients(arraynd > const & gradients) - { - this->m_gradients = gradients; - } + float b = (1.f - lol::sqrt(N + 1.f)) / lol::sqrt((float)N * N * N); + float a = b + lol::sqrt((N + 1.f) / N); - inline arraynd > const & GetGradients() const - { - return this->m_gradients; + // 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; + + for (int i = 0; i < N; ++i) + { + for (int j = 0; j < N; ++j) + { + m_base[i][j] = (i == j ? a : b); + m_base_inverse[i][j] = (i == j ? c : d); + } + } } // Single interpolation @@ -57,20 +75,19 @@ public: } protected: - inline float LastInterp(vec_t origin, vec_t const & pos, vec_t const & index_order) const { - // “corner” will traverse the simplex along its edges in - // orthonormal coordinates + /* “corner” will traverse the simplex along its edges in + * orthonormal coordinates. */ vec_t corner(0); float result = 0; for (int i = 0; i < N + 1; ++i) { vec_t const delta = pos - corner; - vec_t const &gradient = this->m_gradients[origin]; + 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 @@ -85,8 +102,7 @@ protected: if (i < N) { corner += m_base[index_order[i]]; - origin[index_order[i]] = this->Mod(origin[index_order[i]] + 1, - index_order[i]); + origin[index_order[i]] += 1; } } @@ -104,6 +120,7 @@ protected: 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]]) @@ -112,11 +129,53 @@ protected: return result; } - inline int Mod(int value, int index) const + inline vec_t GetGradient(vec_t origin) const { - int const dim = this->m_gradients.GetSize()[index]; - int const ret = value % dim; - return ret >= 0 ? ret : ret + dim; + /* 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, + }; + + /* 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); + }; + + static vec_t const gradients[16] = + { + v(), v(), v(), v(), v(), v(), v(), v(), + v(), v(), v(), v(), v(), v(), v(), v(), + }; + + int idx = m_seed; + for (int i = 0; i < N; ++i) + idx = shuffle[(idx + origin[i]) & 255]; + + return gradients[(idx ^ (idx >> 4)) & 15]; } /* For a given world position, extract grid coordinates (origin) and @@ -131,46 +190,11 @@ protected: // Extracting decimal part from simplex sample pos = simplex_position - (vec_t)origin; - - // Never exceed the size of the gradients and loop on it - for (int i = 0; i < N; ++i) - origin[i] = this->Mod(origin[i], i); - } - - inline void InitBase() - { - /* Matrix coordinate transformation to skew simplex referential is done - by inversing 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; - - for (int i = 0; i < N; ++i) - { - for (int j = 0; j < N; ++j) - { - m_base[i][j] = (i == j ? a : b); - m_base_inverse[i][j] = (i == j ? c : d); - } - } } mat_t m_base, m_base_inverse; - arraynd > m_gradients; + int m_seed; }; } +