SeExpr
Noise.cpp
Go to the documentation of this file.
1 /*
2  Copyright Disney Enterprises, Inc. All rights reserved.
3 
4  Licensed under the Apache License, Version 2.0 (the "License");
5  you may not use this file except in compliance with the License
6  and the following modification to it: Section 6 Trademarks.
7  deleted and replaced with:
8 
9  6. Trademarks. This License does not grant permission to use the
10  trade names, trademarks, service marks, or product names of the
11  Licensor and its affiliates, except as required for reproducing
12  the content of the NOTICE file.
13 
14  You may obtain a copy of the License at
15  http://www.apache.org/licenses/LICENSE-2.0
16 */
17 
18 #include <iostream>
19 #ifdef __SSE4_1__
20 #include <smmintrin.h>
21 #endif
22 #include "ExprBuiltins.h"
23 
24 namespace {
25 #include "NoiseTables.h"
26 }
27 #include "Noise.h"
28 namespace SeExpr2 {
29 
30 #ifdef __SSE4_1__
31 inline double floorSSE(double val)
32 {
33  return _mm_cvtsd_f64(_mm_floor_sd(_mm_set_sd(0.0), _mm_set_sd(val)));
34 }
35 
36 inline double roundSSE(double val)
37 {
38  return _mm_cvtsd_f64(_mm_round_sd(_mm_set_sd(0.0), _mm_set_sd(val), _MM_FROUND_TO_NEAREST_INT));
39 }
40 #else
41 #define floorSSE floor
42 #define roundSSE round
43 #endif
44 
46 double s_curve(double t) { return t * t * t * (t * (6 * t - 15) + 10); }
47 
49 template <int d>
50 unsigned char hashReduceChar(int index[d]) {
51  uint32_t seed = 0;
52  // blend with seed (constants from Numerical Recipes, attrib. from Knuth)
53  for (int k = 0; k < d; k++) {
54  static const uint32_t M = 1664525, C = 1013904223;
55  seed = seed * M + index[k] + C;
56  }
57  // tempering (from Matsumoto)
58  seed ^= (seed >> 11);
59  seed ^= (seed << 7) & 0x9d2c5680UL;
60  seed ^= (seed << 15) & 0xefc60000UL;
61  seed ^= (seed >> 18);
62  // compute one byte by mixing third and first bytes
63  return (((seed & 0xff0000) >> 4) + (seed & 0xff)) & 0xff;
64 }
65 
67 template <int d>
68 uint32_t hashReduce(uint32_t index[d]) {
69  union {
70  uint32_t i;
71  unsigned char c[4];
72  } u1, u2;
73  // blend with seed (constants from Numerical Recipes, attrib. from Knuth)
74  u1.i = 0;
75  for (int k = 0; k < d; k++) {
76  static const uint32_t M = 1664525, C = 1013904223;
77  u1.i = u1.i * M + index[k] + C;
78  }
79  // tempering (from Matsumoto)
80  u1.i ^= (u1.i >> 11);
81  u1.i ^= (u1.i << 7) & 0x9d2c5680U;
82  u1.i ^= (u1.i << 15) & 0xefc60000U;
83  u1.i ^= (u1.i >> 18);
84  // permute bytes (shares perlin noise permutation table)
85  u2.c[3] = p[u1.c[0]];
86  u2.c[2] = p[u1.c[1] + u2.c[3]];
87  u2.c[1] = p[u1.c[2] + u2.c[2]];
88  u2.c[0] = p[u1.c[3] + u2.c[1]];
89 
90  return u2.i;
91 }
92 
94 template <int d, class T, bool periodic>
95 T noiseHelper(const T* X, const int* period = 0) {
96  // find lattice index
97  T weights[2][d]; // lower and upper weights
98  int index[d];
99  for (int k = 0; k < d; k++) {
100  T f = floorSSE(X[k]);
101  index[k] = (int)f;
102  if (periodic) {
103  index[k] %= period[k];
104  if (index[k] < 0) index[k] += period[k];
105  }
106  weights[0][k] = X[k] - f;
107  weights[1][k] = weights[0][k] - 1; // dist to cell with index one above
108  }
109  // compute function values propagated from zero from each node
110  const int num = 1 << d;
111  T vals[num];
112  for (int dummy = 0; dummy < num; dummy++) {
113  int latticeIndex[d];
114  int offset[d];
115  for (int k = 0; k < d; k++) {
116  offset[k] = ((dummy & (1 << k)) != 0);
117  latticeIndex[k] = index[k] + offset[k];
118  }
119  // hash to get representative gradient vector
120  int lookup = hashReduceChar<d>(latticeIndex);
121  T val = 0;
122  for (int k = 0; k < d; k++) {
123  double grad = NOISE_TABLES<d>::g[lookup][k];
124  double weight = weights[offset[k]][k];
125  val += grad * weight;
126  }
127  vals[dummy] = val;
128  }
129  // compute linear interpolation coefficients
130  T alphas[d];
131  for (int k = 0; k < d; k++) alphas[k] = s_curve(weights[0][k]);
132  // perform multilinear interpolation (i.e. linear, bilinear, trilinear, quadralinear)
133  for (int newd = d - 1; newd >= 0; newd--) {
134  int newnum = 1 << newd;
135  int k = (d - newd - 1);
136  T alpha = alphas[k];
137  T beta = T(1) - alphas[k];
138  for (int dummy = 0; dummy < newnum; dummy++) {
139  int index = dummy * (1 << (d - newd));
140  int otherIndex = index + (1 << k);
141  // T alpha=s_curve(weights[0][k]);
142  vals[index] = beta * vals[index] + alpha * vals[otherIndex];
143  }
144  }
145  // return reduced version
146  return vals[0];
147 }
148 }
149 
150 namespace SeExpr2 {
151 
153 template <int d_in, int d_out, class T>
154 void CellNoise(const T* in, T* out) {
155  uint32_t index[d_in];
156  int dim = 0;
157  for (int k = 0; k < d_in; k++) index[k] = uint32_t(floorSSE(in[k]));
158  while (1) {
159  out[dim] = hashReduce<d_in>(index) * (1.0 / 0xffffffffu);
160  if (++dim >= d_out) break;
161  for (int k = 0; k < d_in; k++) index[k] += 1000;
162  }
163 }
164 
166 template <int d_in, int d_out, class T>
167 void Noise(const T* in, T* out) {
168  T P[d_in];
169  for (int i = 0; i < d_in; i++) P[i] = in[i];
170 
171  int i = 0;
172  while (1) {
173  out[i] = noiseHelper<d_in, T, false>(P);
174  if (++i >= d_out) break;
175  for (int k = 0; k < d_out; k++) P[k] += (T)1000;
176  }
177 }
178 
180 template <int d_in, int d_out, class T>
181 void PNoise(const T* in, const int* period, T* out) {
182  T P[d_in];
183  for (int i = 0; i < d_in; i++) P[i] = in[i];
184 
185  int i = 0;
186  while (1) {
187  out[i] = noiseHelper<d_in, T, true>(P, period);
188  if (++i >= d_out) break;
189  for (int k = 0; k < d_out; k++) P[k] += (T)1000;
190  }
191 }
192 
195 template <int d_in, int d_out, bool turbulence, class T>
196 void FBM(const T* in, T* out, int octaves, T lacunarity, T gain) {
197  T P[d_in];
198  for (int i = 0; i < d_in; i++) P[i] = in[i];
199 
200  T scale = 1;
201  for (int k = 0; k < d_out; k++) out[k] = 0;
202  int octave = 0;
203  while (1) {
204  T localResult[d_out];
205  Noise<d_in, d_out>(P, localResult);
206  if (turbulence)
207  for (int k = 0; k < d_out; k++) out[k] += fabs(localResult[k]) * scale;
208  else
209  for (int k = 0; k < d_out; k++) out[k] += localResult[k] * scale;
210  if (++octave >= octaves) break;
211  scale *= gain;
212  for (int k = 0; k < d_in; k++) {
213  P[k] *= lacunarity;
214  P[k] += (T)1234;
215  }
216  }
217 }
218 
219 // Explicit instantiations
220 template void CellNoise<3, 1, double>(const double*, double*);
221 template void CellNoise<3, 3, double>(const double*, double*);
222 template void Noise<1, 1, double>(const double*, double*);
223 template void Noise<2, 1, double>(const double*, double*);
224 template void Noise<3, 1, double>(const double*, double*);
225 template void PNoise<3, 1, double>(const double*, const int*, double*);
226 template void Noise<4, 1, double>(const double*, double*);
227 template void Noise<3, 3, double>(const double*, double*);
228 template void Noise<4, 3, double>(const double*, double*);
229 template void FBM<3, 1, false, double>(const double*, double*, int, double, double);
230 template void FBM<3, 1, true, double>(const double*, double*, int, double, double);
231 template void FBM<3, 3, false, double>(const double*, double*, int, double, double);
232 template void FBM<3, 3, true, double>(const double*, double*, int, double, double);
233 template void FBM<4, 1, false, double>(const double*, double*, int, double, double);
234 template void FBM<4, 3, false, double>(const double*, double*, int, double, double);
235 }
236 
237 #ifdef MAINTEST
238 int main(int argc, char* argv[]) {
239  typedef double T;
240  T sum = 0;
241  for (int i = 0; i < 10000000; i++) {
242  T foo[3] = {.3, .3, .3};
243  // for(int k=0;k<3;k++) foo[k]=(double)rand()/double(RAND_MAX)*100.;
244  sum += SeExpr2::noiseHelper<3, T, false>(foo);
245  }
246 }
247 #endif
index
The result is computed int int< br >< div style="margin-left: 40px;"> Picks values randomly between loRange and hiRange based on supplied index(which is automatically hashed). &nbsp
f
with numParticles numAttributes A variable block contains variable names and types but doesn t care what the values are< pre > void f(const std::string &s, MyParticleData *p, int outputDim=3)
Definition: varblocks.txt:35
SeExpr2::FBM< 4, 3, false, double >
template void FBM< 4, 3, false, double >(const double *, double *, int, double, double)
SeExpr2::noiseHelper
T noiseHelper(const T *X, const int *period=0)
Noise with d_in dimensional domain, 1 dimensional abcissa.
Definition: Noise.cpp:95
roundSSE
#define roundSSE
Definition: Noise.cpp:42
SeExpr2::CellNoise
void CellNoise(const T *in, T *out)
Computes cellular noise (non-interpolated piecewise constant cell random values)
Definition: Noise.cpp:154
SeExpr2::FBM< 3, 1, false, double >
template void FBM< 3, 1, false, double >(const double *, double *, int, double, double)
SeExpr2::PNoise< 3, 1, double >
template void PNoise< 3, 1, double >(const double *, const int *, double *)
SeExpr2::PNoise
void PNoise(const T *in, const int *period, T *out)
Periodic Noise with d_in dimensional domain, d_out dimensional abcissa.
Definition: Noise.cpp:181
NoiseTables.h
NOISE_TABLES
Definition: NoiseTables.h:46
SeExpr2
Definition: Context.h:22
SeExpr2::FBM< 3, 1, true, double >
template void FBM< 3, 1, true, double >(const double *, double *, int, double, double)
SeExpr2::CellNoise< 3, 1, double >
template void CellNoise< 3, 1, double >(const double *, double *)
ExprBuiltins.h
SeExpr2::Noise< 3, 3, double >
template void Noise< 3, 3, double >(const double *, double *)
SeExpr2::Noise
void Noise(const T *in, T *out)
Noise with d_in dimensional domain, d_out dimensional abcissa.
Definition: Noise.cpp:167
floorSSE
#define floorSSE
Definition: Noise.cpp:41
SeExpr2::FBM< 3, 3, false, double >
template void FBM< 3, 3, false, double >(const double *, double *, int, double, double)
SeExpr2::FBM< 3, 3, true, double >
template void FBM< 3, 3, true, double >(const double *, double *, int, double, double)
SeExpr2::hashReduce
uint32_t hashReduce(uint32_t index[d])
Does a hash reduce to an integer.
Definition: Noise.cpp:68
SeExpr2::FBM
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:196
SeExpr2::Noise< 1, 1, double >
template void Noise< 1, 1, double >(const double *, double *)
SeExpr2::Noise< 4, 1, double >
template void Noise< 4, 1, double >(const double *, double *)
p
static const int p[514]
Definition: NoiseTables.h:20
SeExpr2::Noise< 4, 3, double >
template void Noise< 4, 3, double >(const double *, double *)
SeExpr2::FBM< 4, 1, false, double >
template void FBM< 4, 1, false, double >(const double *, double *, int, double, double)
SeExpr2::Noise< 3, 1, double >
template void Noise< 3, 1, double >(const double *, double *)
main
int main(int argc, char *argv[])
Definition: EditMain.cpp:24
Noise.h
SeExpr2::turbulence
double turbulence(int n, const Vec3d *args)
Definition: ExprBuiltins.cpp:551
SeExpr2::s_curve
double s_curve(double t)
This is the Quintic interpolant from Perlin's Improved Noise Paper.
Definition: Noise.cpp:46
SeExpr2::Noise< 2, 1, double >
template void Noise< 2, 1, double >(const double *, double *)
SeExpr2::CellNoise< 3, 3, double >
template void CellNoise< 3, 3, double >(const double *, double *)
SeExpr2::hashReduceChar
unsigned char hashReduceChar(int index[d])
Does a hash reduce to a character.
Definition: Noise.cpp:50