lolengine
/
lol
miroir de https://github.com/lolengine/lol
1
0
Bifurcation 0
Code Tickets 0 Versions 0 Wiki Activité
Parcourir la source

math: tweak simplex noise scale according to dimension.

undefined
Sam Hocevar il y a 10 ans
Parent
7faf5d912a
révision
b0b5bcc6fa
3 fichiers modifiés avec 159 ajouts et 385 suppressions
Vue séparée
Option de Diff
Voir les Statistiques Télécharger le Fichier Patch Télécharger le Fichier des Différences
  1. +33
    -14
      demos/test/simplex.cpp
  2. +125
    -27
      src/lol/math/simplex_interpolator.h
  3. +1
    -344
      src/t/math/simplex_interpolator.cpp

+ 33
- 14
demos/test/simplex.cpp Voir le fichier

@@ -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<vec4> &data = img.Lock2D<PixelFormat::RGBA_F32>();

@@ -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));



+ 125
- 27
src/lol/math/simplex_interpolator.h Voir le fichier

@@ -50,8 +50,8 @@ public:
#endif
}

/* Single interpolation */
inline float Interp(vec_t<float, N> position) const
/* Evaluate noise at a given point */
inline float eval(vec_t<float, N> position) const
{
// Retrieve the containing hypercube origin and associated decimals
vec_t<int, N> origin;
@@ -61,8 +61,9 @@ public:
return get_noise(origin, pos);
}

/* Only for debug purposes: return the gradient vector */
inline vec_t<float, N> GetGradient(vec_t<float, N> position) const
/* Only for debug purposes: return the gradient vector of the given
* point’s simplex origin. */
inline vec_t<float, N> gradient(vec_t<float, N> position) const
{
vec_t<int, N> origin;
vec_t<float, N> 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<float, N> 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<float, N> get_gradient(vec_t<int, N> 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<float, N, N> 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<float, N> 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<vec_t<float, N>> deltas;
for (int i = 0; i < N; ++i)
{
vec_t<float, N> v(0.f);
v[i] = 1.f;
deltas << v << -v;
}
for (int run = 0; run < 1000; ++run)
{
/* Pick a random vector */
vec_t<float, N> 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");
}



+ 1
- 344
src/t/math/simplex_interpolator.cpp Voir le fichier

@@ -19,354 +19,11 @@
namespace lol
{

template<int N>
class test_interpolator : public simplex_interpolator<N>
{
public:

test_interpolator() :
simplex_interpolator<N>()
{
}

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<int, N> GetIndexOrder(vec_t<float, N> const & decimal_point)
{
return simplex_interpolator<N>::GetIndexOrder(decimal_point);
}
};

lolunit_declare_fixture(SimplexInterpolatorTest)
lolunit_declare_fixture(SimplexevalolatorTest)
{
void SetUp() {}

void TearDown() {}

template<int N>
void check_base_matrix()
{
test_interpolator<N> 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<int N>
void check_index_ordering()
{
static int gen = 12345678;
test_interpolator<N> s;

vec_t<float, N> vec_ref;
for (int i = 0 ; i < N ; ++i)
vec_ref[i] = (gen = gen * gen + gen);

vec_t<int, N> 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<float, 2> > gradients(vec_t<ptrdiff_t, 2>({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<float, 2> 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<float, 2>({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
};

}

Écrire Aperçu
Chargement…
Annuler
Enregistrer