KSeExpr  4.0.4.0
Noise.cpp
Go to the documentation of this file.
1 // SPDX-FileCopyrightText: 2011-2019 Disney Enterprises, Inc.
2 // SPDX-License-Identifier: LicenseRef-Apache-2.0
3 // SPDX-FileCopyrightText: 2020 L. E. Segovia <amy@amyspark.me>
4 // SPDX-License-Identifier: GPL-3.0-or-later
5 
6 #include "ExprBuiltins.h"
7 #include <iostream>
8 
9 namespace
10 {
11 #include "NoiseTables.h"
12 }
13 #include "Noise.h"
14 #include "Utils.h"
15 namespace KSeExpr
16 {
18 double s_curve(double t)
19 {
20  return t * t * t * (t * (6 * t - 15) + 10);
21 }
22 
24 template<int d> unsigned char hashReduceChar(int index[d])
25 {
26  uint32_t seed = 0;
27  // blend with seed (constants from Numerical Recipes, attrib. from Knuth)
28  for (int k = 0; k < d; k++) {
29  static const uint32_t M = 1664525, C = 1013904223;
30  seed = seed * M + index[k] + C;
31  }
32  // tempering (from Matsumoto)
33  seed ^= (seed >> 11);
34  seed ^= (seed << 7) & 0x9d2c5680UL;
35  seed ^= (seed << 15) & 0xefc60000UL;
36  seed ^= (seed >> 18);
37  // compute one byte by mixing third and first bytes
38  return (((seed & 0xff0000) >> 4) + (seed & 0xff)) & 0xff;
39 }
40 
42 template<int d> uint32_t hashReduce(uint32_t index[d])
43 {
44  union {
45  uint32_t i;
46  unsigned char c[4];
47  } u1, u2;
48  // blend with seed (constants from Numerical Recipes, attrib. from Knuth)
49  u1.i = 0;
50  for (int k = 0; k < d; k++) {
51  static const uint32_t M = 1664525, C = 1013904223;
52  u1.i = u1.i * M + index[k] + C;
53  }
54  // tempering (from Matsumoto)
55  u1.i ^= (u1.i >> 11);
56  u1.i ^= (u1.i << 7) & 0x9d2c5680U;
57  u1.i ^= (u1.i << 15) & 0xefc60000U;
58  u1.i ^= (u1.i >> 18);
59  // permute bytes (shares perlin noise permutation table)
60  u2.c[3] = p[u1.c[0]];
61  u2.c[2] = p[u1.c[1] + u2.c[3]];
62  u2.c[1] = p[u1.c[2] + u2.c[2]];
63  u2.c[0] = p[u1.c[3] + u2.c[1]];
64 
65  return u2.i;
66 }
67 
69 template<int d, class T, bool periodic> T noiseHelper(const T *X, const int *period = 0)
70 {
71  // find lattice index
72  T weights[2][d]; // lower and upper weights
73  int index[d];
74  for (int k = 0; k < d; k++) {
75  T f = KSeExpr::Utils::floor(X[k]);
76  index[k] = (int)f;
77  if (periodic) {
78  index[k] %= period[k];
79  if (index[k] < 0)
80  index[k] += period[k];
81  }
82  weights[0][k] = X[k] - f;
83  weights[1][k] = weights[0][k] - 1; // dist to cell with index one above
84  }
85  // compute function values propagated from zero from each node
86  const int num = 1 << d;
87  T vals[num];
88  for (int dummy = 0; dummy < num; dummy++) {
89  int latticeIndex[d];
90  int offset[d];
91  for (int k = 0; k < d; k++) {
92  offset[k] = ((dummy & (1 << k)) != 0);
93  latticeIndex[k] = index[k] + offset[k];
94  }
95  // hash to get representative gradient vector
96  int lookup = hashReduceChar<d>(latticeIndex);
97  T val = 0;
98  for (int k = 0; k < d; k++) {
99  double grad = NOISE_TABLES<d>::g[lookup][k];
100  double weight = weights[offset[k]][k];
101  val += grad * weight;
102  }
103  vals[dummy] = val;
104  }
105  // compute linear interpolation coefficients
106  T alphas[d];
107  for (int k = 0; k < d; k++)
108  alphas[k] = s_curve(weights[0][k]);
109  // perform multilinear interpolation (i.e. linear, bilinear, trilinear, quadralinear)
110  for (int newd = d - 1; newd >= 0; newd--) {
111  int newnum = 1 << newd;
112  int k = (d - newd - 1);
113  T alpha = alphas[k];
114  T beta = T(1) - alphas[k];
115  for (int dummy = 0; dummy < newnum; dummy++) {
116  int index = dummy * (1 << (d - newd));
117  int otherIndex = index + (1 << k);
118  // T alpha=s_curve(weights[0][k]);
119  vals[index] = beta * vals[index] + alpha * vals[otherIndex];
120  }
121  }
122  // return reduced version
123  return vals[0];
124 }
125 }
126 
127 namespace KSeExpr
128 {
130 template<int d_in, int d_out, class T> void CellNoise(const T *in, T *out)
131 {
132  uint32_t index[d_in];
133  int dim = 0;
134  for (int k = 0; k < d_in; k++)
135  index[k] = uint32_t(KSeExpr::Utils::floor(in[k]));
136  while (1) {
137  out[dim] = hashReduce<d_in>(index) * (1.0 / 0xffffffffu);
138  if (++dim >= d_out)
139  break;
140  for (int k = 0; k < d_in; k++)
141  index[k] += 1000;
142  }
143 }
144 
146 template<int d_in, int d_out, class T> void Noise(const T *in, T *out)
147 {
148  T P[d_in];
149  for (int i = 0; i < d_in; i++)
150  P[i] = in[i];
151 
152  int i = 0;
153  while (1) {
154  out[i] = noiseHelper<d_in, T, false>(P);
155  if (++i >= d_out)
156  break;
157  for (int k = 0; k < d_out; k++)
158  P[k] += (T)1000;
159  }
160 }
161 
163 template<int d_in, int d_out, class T> void PNoise(const T *in, const int *period, T *out)
164 {
165  T P[d_in];
166  for (int i = 0; i < d_in; i++)
167  P[i] = in[i];
168 
169  int i = 0;
170  while (1) {
171  out[i] = noiseHelper<d_in, T, true>(P, period);
172  if (++i >= d_out)
173  break;
174  for (int k = 0; k < d_out; k++)
175  P[k] += (T)1000;
176  }
177 }
178 
181 template<int d_in, int d_out, bool turbulence, class T> void FBM(const T *in, T *out, int octaves, T lacunarity, T gain)
182 {
183  T P[d_in];
184  for (int i = 0; i < d_in; i++)
185  P[i] = in[i];
186 
187  T scale = 1;
188  for (int k = 0; k < d_out; k++)
189  out[k] = 0;
190  int octave = 0;
191  while (1) {
192  T localResult[d_out];
193  Noise<d_in, d_out>(P, localResult);
194  if (turbulence)
195  for (int k = 0; k < d_out; k++)
196  out[k] += fabs(localResult[k]) * scale;
197  else
198  for (int k = 0; k < d_out; k++)
199  out[k] += localResult[k] * scale;
200  if (++octave >= octaves)
201  break;
202  scale *= gain;
203  for (int k = 0; k < d_in; k++) {
204  P[k] *= lacunarity;
205  P[k] += (T)1234;
206  }
207  }
208 }
209 
210 // Explicit instantiations
211 template void CellNoise<3, 1, double>(const double *, double *);
212 template void CellNoise<3, 3, double>(const double *, double *);
213 template void Noise<1, 1, double>(const double *, double *);
214 template void Noise<2, 1, double>(const double *, double *);
215 template void Noise<3, 1, double>(const double *, double *);
216 template void PNoise<3, 1, double>(const double *, const int *, double *);
217 template void Noise<4, 1, double>(const double *, double *);
218 template void Noise<3, 3, double>(const double *, double *);
219 template void Noise<4, 3, double>(const double *, double *);
220 template void FBM<3, 1, false, double>(const double *, double *, int, double, double);
221 template void FBM<3, 1, true, double>(const double *, double *, int, double, double);
222 template void FBM<3, 3, false, double>(const double *, double *, int, double, double);
223 template void FBM<3, 3, true, double>(const double *, double *, int, double, double);
224 template void FBM<4, 1, false, double>(const double *, double *, int, double, double);
225 template void FBM<4, 3, false, double>(const double *, double *, int, double, double);
226 }
227 
228 #ifdef MAINTEST
229 int main(int argc, char *argv[])
230 {
231  typedef double T;
232  T sum = 0;
233  for (int i = 0; i < 10000000; i++) {
234  T foo[3] = {.3, .3, .3};
235  // for(int k=0;k<3;k++) foo[k]=(double)rand()/double(RAND_MAX)*100.;
236  sum += KSeExpr::noiseHelper<3, T, false>(foo);
237  }
238 }
239 #endif
static constexpr std::array< int, 514 > p
Definition: NoiseTables.h:10
int main(int, char *[])
KSeExpr_DEFAULT double_t floor(double_t val)
Definition: Utils.cpp:168
double s_curve(double t)
This is the Quintic interpolant from Perlin's Improved Noise Paper.
Definition: Noise.cpp:18
template void Noise< 2, 1, double >(const double *, double *)
template void FBM< 4, 1, false, double >(const double *, double *, int, double, double)
template void PNoise< 3, 1, double >(const double *, const int *, double *)
void Noise(const T *in, T *out)
Noise with d_in dimensional domain, d_out dimensional abcissa.
Definition: Noise.cpp:146
template void FBM< 3, 3, false, double >(const double *, double *, int, double, double)
template void FBM< 3, 1, true, double >(const double *, double *, int, double, double)
template void Noise< 1, 1, double >(const double *, double *)
template void FBM< 3, 1, false, double >(const double *, double *, int, double, double)
template void CellNoise< 3, 1, double >(const double *, double *)
void CellNoise(const T *in, T *out)
Computes cellular noise (non-interpolated piecewise constant cell random values)
Definition: Noise.cpp:130
template void CellNoise< 3, 3, double >(const double *, double *)
T noiseHelper(const T *X, const int *period=0)
Noise with d_in dimensional domain, 1 dimensional abcissa.
Definition: Noise.cpp:69
unsigned char hashReduceChar(int index[d])
Does a hash reduce to a character.
Definition: Noise.cpp:24
void PNoise(const T *in, const int *period, T *out)
Periodic Noise with d_in dimensional domain, d_out dimensional abcissa.
Definition: Noise.cpp:163
template void FBM< 4, 3, false, double >(const double *, double *, int, double, double)
double turbulence(int n, const Vec3d *args)
void FBM(const T *in, T *out, int octaves, T lacunarity, T gain)
Fractional Brownian Motion. If turbulence is true then turbulence computed.
Definition: Noise.cpp:181
uint32_t hashReduce(uint32_t index[d])
Does a hash reduce to an integer.
Definition: Noise.cpp:42
template void Noise< 3, 3, double >(const double *, double *)
template void FBM< 3, 3, true, double >(const double *, double *, int, double, double)
template void Noise< 3, 1, double >(const double *, double *)
template void Noise< 4, 1, double >(const double *, double *)
template void Noise< 4, 3, double >(const double *, double *)