1 
2 
3 // DO NOT EDIT !
4 // This file is generated using the MantaFlow preprocessor (prep generate).
5 
6 /******************************************************************************
7  *
8  * MantaFlow fluid solver framework
9  * Copyright 2011 Tobias Pfaff, Nils Thuerey
10  *
11  * This program is free software, distributed under the terms of the
12  * Apache License, Version 2.0
13  * http://www.apache.org/licenses/LICENSE-2.0
14  *
15  * Wavelet noise field
16  *
17  ******************************************************************************/
18 
19 #ifndef _NOISEFIELD_H_
20 #define _NOISEFIELD_H_
21 
22 #include "vectorbase.h"
23 #include "manta.h"
24 #include "grid.h"
25 #include <atomic>
26 
27 namespace Manta {
28 
29 #define NOISE_TILE_SIZE 128
30 
31 // wrapper for a parametrized field of wavelet noise
32 
33 class WaveletNoiseField : public PbClass {
34  public:
35   WaveletNoiseField(FluidSolver *parent, int fixedSeed = -1, int loadFromFile = false);
_W_0(PyObject * _self,PyObject * _linargs,PyObject * _kwds)36   static int _W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
37   {
38     PbClass *obj = Pb::objFromPy(_self);
39     if (obj)
40       delete obj;
41     try {
42       PbArgs _args(_linargs, _kwds);
43       bool noTiming = _args.getOpt<bool>("notiming", -1, 0);
44       pbPreparePlugin(0, "WaveletNoiseField::WaveletNoiseField", !noTiming);
45       {
46         ArgLocker _lock;
47         FluidSolver *parent = _args.getPtr<FluidSolver>("parent", 0, &_lock);
48         int fixedSeed = _args.getOpt<int>("fixedSeed", 1, -1, &_lock);
49         int loadFromFile = _args.getOpt<int>("loadFromFile", 2, false, &_lock);
50         obj = new WaveletNoiseField(parent, fixedSeed, loadFromFile);
51         obj->registerObject(_self, &_args);
52         _args.check();
53       }
54       pbFinalizePlugin(obj->getParent(), "WaveletNoiseField::WaveletNoiseField", !noTiming);
55       return 0;
56     }
57     catch (std::exception &e) {
58       pbSetError("WaveletNoiseField::WaveletNoiseField", e.what());
59       return -1;
60     }
61   }
62 
~WaveletNoiseField()63   ~WaveletNoiseField()
64   {
65     if (mNoiseTile && !mNoiseReferenceCount) {
66       delete mNoiseTile;
67       mNoiseTile = NULL;
68     }
69   };
70 
71   //! evaluate noise
72   inline Real evaluate(Vec3 pos, int tile = 0) const;
73   //! evaluate noise as a vector
74   inline Vec3 evaluateVec(Vec3 pos, int tile = 0) const;
75   //! evaluate curl noise
76   inline Vec3 evaluateCurl(Vec3 pos) const;
77 
78   //! direct data access
data()79   Real *data()
80   {
81     return mNoiseTile;
82   }
83 
84   //! compute wavelet decomposition of an input grid (stores residual coefficients)
85   static void computeCoefficients(Grid<Real> &input, Grid<Real> &tempIn1, Grid<Real> &tempIn2);
86 
87   // helper
88   std::string toString();
89 
90   // texcoord position and scale
91   Vec3 mPosOffset;
_GET_mPosOffset(PyObject * self,void * cl)92   static PyObject *_GET_mPosOffset(PyObject *self, void *cl)
93   {
94     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
95     return toPy(pbo->mPosOffset);
96   }
_SET_mPosOffset(PyObject * self,PyObject * val,void * cl)97   static int _SET_mPosOffset(PyObject *self, PyObject *val, void *cl)
98   {
99     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
100     pbo->mPosOffset = fromPy<Vec3>(val);
101     return 0;
102   }
103 
104   Vec3 mPosScale;
_GET_mPosScale(PyObject * self,void * cl)105   static PyObject *_GET_mPosScale(PyObject *self, void *cl)
106   {
107     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
108     return toPy(pbo->mPosScale);
109   }
_SET_mPosScale(PyObject * self,PyObject * val,void * cl)110   static int _SET_mPosScale(PyObject *self, PyObject *val, void *cl)
111   {
112     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
113     pbo->mPosScale = fromPy<Vec3>(val);
114     return 0;
115   }
116 
117   // value offset & scale
118   Real mValOffset;
_GET_mValOffset(PyObject * self,void * cl)119   static PyObject *_GET_mValOffset(PyObject *self, void *cl)
120   {
121     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
122     return toPy(pbo->mValOffset);
123   }
_SET_mValOffset(PyObject * self,PyObject * val,void * cl)124   static int _SET_mValOffset(PyObject *self, PyObject *val, void *cl)
125   {
126     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
127     pbo->mValOffset = fromPy<Real>(val);
128     return 0;
129   }
130 
131   Real mValScale;
_GET_mValScale(PyObject * self,void * cl)132   static PyObject *_GET_mValScale(PyObject *self, void *cl)
133   {
134     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
135     return toPy(pbo->mValScale);
136   }
_SET_mValScale(PyObject * self,PyObject * val,void * cl)137   static int _SET_mValScale(PyObject *self, PyObject *val, void *cl)
138   {
139     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
140     pbo->mValScale = fromPy<Real>(val);
141     return 0;
142   }
143 
144   // clamp? (default 0-1)
145   bool mClamp;
_GET_mClamp(PyObject * self,void * cl)146   static PyObject *_GET_mClamp(PyObject *self, void *cl)
147   {
148     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
149     return toPy(pbo->mClamp);
150   }
_SET_mClamp(PyObject * self,PyObject * val,void * cl)151   static int _SET_mClamp(PyObject *self, PyObject *val, void *cl)
152   {
153     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
154     pbo->mClamp = fromPy<bool>(val);
155     return 0;
156   }
157 
158   Real mClampNeg;
_GET_mClampNeg(PyObject * self,void * cl)159   static PyObject *_GET_mClampNeg(PyObject *self, void *cl)
160   {
161     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
162     return toPy(pbo->mClampNeg);
163   }
_SET_mClampNeg(PyObject * self,PyObject * val,void * cl)164   static int _SET_mClampNeg(PyObject *self, PyObject *val, void *cl)
165   {
166     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
167     pbo->mClampNeg = fromPy<Real>(val);
168     return 0;
169   }
170 
171   Real mClampPos;
_GET_mClampPos(PyObject * self,void * cl)172   static PyObject *_GET_mClampPos(PyObject *self, void *cl)
173   {
174     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
175     return toPy(pbo->mClampPos);
176   }
_SET_mClampPos(PyObject * self,PyObject * val,void * cl)177   static int _SET_mClampPos(PyObject *self, PyObject *val, void *cl)
178   {
179     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
180     pbo->mClampPos = fromPy<Real>(val);
181     return 0;
182   }
183 
184   // animated over time
185   Real mTimeAnim;
_GET_mTimeAnim(PyObject * self,void * cl)186   static PyObject *_GET_mTimeAnim(PyObject *self, void *cl)
187   {
188     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
189     return toPy(pbo->mTimeAnim);
190   }
_SET_mTimeAnim(PyObject * self,PyObject * val,void * cl)191   static int _SET_mTimeAnim(PyObject *self, PyObject *val, void *cl)
192   {
193     WaveletNoiseField *pbo = dynamic_cast<WaveletNoiseField *>(Pb::objFromPy(self));
194     pbo->mTimeAnim = fromPy<Real>(val);
195     return 0;
196   }
197 
198  protected:
199   // noise evaluation functions
200   static inline Real WNoiseDx(const Vec3 &p, Real *data);
201   static inline Vec3 WNoiseVec(const Vec3 &p, Real *data);
202   static inline Real WNoise(const Vec3 &p, Real *data);
203 
204   // helpers for tile generation , for periodic 128 grids only
205   static void downsample(Real *from, Real *to, int n, int stride);
206   static void upsample(Real *from, Real *to, int n, int stride);
207 
208   // for grids with arbitrary sizes, and neumann boundary conditions
209   static void downsampleNeumann(const Real *from, Real *to, int n, int stride);
210   static void upsampleNeumann(const Real *from, Real *to, int n, int stride);
211 
modSlow(int x,int n)212   static inline int modSlow(int x, int n)
213   {
214     int m = x % n;
215     return (m < 0) ? m + n : m;
216   }
217 // warning - noiseTileSize has to be 128^3!
218 #define modFast128(x) ((x)&127)
219 
getTime()220   inline Real getTime() const
221   {
222     return mParent->getTime() * mParent->getDx() * mTimeAnim;
223   }
224 
225   // pre-compute tile data for wavelet noise
226   void generateTile(int loadFromFile);
227 
228   // animation over time
229   // grid size normalization (inverse size)
230   Real mGsInvX, mGsInvY, mGsInvZ;
231   // random offset into tile to simulate different random seeds
232   Vec3 mSeedOffset;
233 
234   static Real *mNoiseTile;
235   // global random seed storage
236   static int randomSeed;
237   // global reference count for noise tile
238   static std::atomic<int> mNoiseReferenceCount;
239  public:
240   PbArgs _args;
241 }
242 #define _C_WaveletNoiseField
243 ;
244 
245 // **************************************************************************
246 // Implementation
247 
248 #define ADD_WEIGHTED(x, y, z) \
249   weight = 1.0f; \
250   xC = modFast128(midX + (x)); \
251   weight *= w[0][(x) + 1]; \
252   yC = modFast128(midY + (y)); \
253   weight *= w[1][(y) + 1]; \
254   zC = modFast128(midZ + (z)); \
255   weight *= w[2][(z) + 1]; \
256   result += weight * data[(zC * NOISE_TILE_SIZE + yC) * NOISE_TILE_SIZE + xC];
257 
258 //////////////////////////////////////////////////////////////////////////////////////////
259 // derivatives of 3D noise - unrolled for performance
260 //////////////////////////////////////////////////////////////////////////////////////////
WNoiseDx(const Vec3 & p,Real * data)261 inline Real WaveletNoiseField::WNoiseDx(const Vec3 &p, Real *data)
262 {
263   Real w[3][3], t, result = 0;
264 
265   // Evaluate quadratic B-spline basis functions
266   int midX = (int)ceil(p[0] - 0.5f);
267   t = midX - (p[0] - 0.5f);
268   w[0][0] = -t;
269   w[0][2] = (1.f - t);
270   w[0][1] = 2.0f * t - 1.0f;
271 
272   int midY = (int)ceil(p[1] - 0.5f);
273   t = midY - (p[1] - 0.5f);
274   w[1][0] = t * t * 0.5f;
275   w[1][2] = (1.f - t) * (1.f - t) * 0.5f;
276   w[1][1] = 1.f - w[1][0] - w[1][2];
277 
278   int midZ = (int)ceil(p[2] - 0.5f);
279   t = midZ - (p[2] - 0.5f);
280   w[2][0] = t * t * 0.5f;
281   w[2][2] = (1.f - t) * (1.f - t) * 0.5f;
282   w[2][1] = 1.f - w[2][0] - w[2][2];
283 
284   // Evaluate noise by weighting noise coefficients by basis function values
285   int xC, yC, zC;
286   Real weight = 1;
287 
288   ADD_WEIGHTED(-1, -1, -1);
289   ADD_WEIGHTED(0, -1, -1);
290   ADD_WEIGHTED(1, -1, -1);
291   ADD_WEIGHTED(-1, 0, -1);
292   ADD_WEIGHTED(0, 0, -1);
293   ADD_WEIGHTED(1, 0, -1);
294   ADD_WEIGHTED(-1, 1, -1);
295   ADD_WEIGHTED(0, 1, -1);
296   ADD_WEIGHTED(1, 1, -1);
297 
298   ADD_WEIGHTED(-1, -1, 0);
299   ADD_WEIGHTED(0, -1, 0);
300   ADD_WEIGHTED(1, -1, 0);
301   ADD_WEIGHTED(-1, 0, 0);
302   ADD_WEIGHTED(0, 0, 0);
303   ADD_WEIGHTED(1, 0, 0);
304   ADD_WEIGHTED(-1, 1, 0);
305   ADD_WEIGHTED(0, 1, 0);
306   ADD_WEIGHTED(1, 1, 0);
307 
308   ADD_WEIGHTED(-1, -1, 1);
309   ADD_WEIGHTED(0, -1, 1);
310   ADD_WEIGHTED(1, -1, 1);
311   ADD_WEIGHTED(-1, 0, 1);
312   ADD_WEIGHTED(0, 0, 1);
313   ADD_WEIGHTED(1, 0, 1);
314   ADD_WEIGHTED(-1, 1, 1);
315   ADD_WEIGHTED(0, 1, 1);
316   ADD_WEIGHTED(1, 1, 1);
317 
318   return result;
319 }
320 
WNoise(const Vec3 & p,Real * data)321 inline Real WaveletNoiseField::WNoise(const Vec3 &p, Real *data)
322 {
323   Real w[3][3], t, result = 0;
324 
325   // Evaluate quadratic B-spline basis functions
326   int midX = (int)ceilf(p[0] - 0.5f);
327   t = midX - (p[0] - 0.5f);
328   w[0][0] = t * t * 0.5f;
329   w[0][2] = (1.f - t) * (1.f - t) * 0.5f;
330   w[0][1] = 1.f - w[0][0] - w[0][2];
331 
332   int midY = (int)ceilf(p[1] - 0.5f);
333   t = midY - (p[1] - 0.5f);
334   w[1][0] = t * t * 0.5f;
335   w[1][2] = (1.f - t) * (1.f - t) * 0.5f;
336   w[1][1] = 1.f - w[1][0] - w[1][2];
337 
338   int midZ = (int)ceilf(p[2] - 0.5f);
339   t = midZ - (p[2] - 0.5f);
340   w[2][0] = t * t * 0.5f;
341   w[2][2] = (1.f - t) * (1.f - t) * 0.5f;
342   w[2][1] = 1.f - w[2][0] - w[2][2];
343 
344   // Evaluate noise by weighting noise coefficients by basis function values
345   int xC, yC, zC;
346   Real weight = 1;
347 
348   ADD_WEIGHTED(-1, -1, -1);
349   ADD_WEIGHTED(0, -1, -1);
350   ADD_WEIGHTED(1, -1, -1);
351   ADD_WEIGHTED(-1, 0, -1);
352   ADD_WEIGHTED(0, 0, -1);
353   ADD_WEIGHTED(1, 0, -1);
354   ADD_WEIGHTED(-1, 1, -1);
355   ADD_WEIGHTED(0, 1, -1);
356   ADD_WEIGHTED(1, 1, -1);
357 
358   ADD_WEIGHTED(-1, -1, 0);
359   ADD_WEIGHTED(0, -1, 0);
360   ADD_WEIGHTED(1, -1, 0);
361   ADD_WEIGHTED(-1, 0, 0);
362   ADD_WEIGHTED(0, 0, 0);
363   ADD_WEIGHTED(1, 0, 0);
364   ADD_WEIGHTED(-1, 1, 0);
365   ADD_WEIGHTED(0, 1, 0);
366   ADD_WEIGHTED(1, 1, 0);
367 
368   ADD_WEIGHTED(-1, -1, 1);
369   ADD_WEIGHTED(0, -1, 1);
370   ADD_WEIGHTED(1, -1, 1);
371   ADD_WEIGHTED(-1, 0, 1);
372   ADD_WEIGHTED(0, 0, 1);
373   ADD_WEIGHTED(1, 0, 1);
374   ADD_WEIGHTED(-1, 1, 1);
375   ADD_WEIGHTED(0, 1, 1);
376   ADD_WEIGHTED(1, 1, 1);
377 
378   return result;
379 }
380 
381 #define ADD_WEIGHTEDX(x, y, z) \
382   weight = dw[0][(x) + 1] * w[1][(y) + 1] * w[2][(z) + 1]; \
383   result += weight * neighbors[x + 1][y + 1][z + 1];
384 
385 #define ADD_WEIGHTEDY(x, y, z) \
386   weight = w[0][(x) + 1] * dw[1][(y) + 1] * w[2][(z) + 1]; \
387   result += weight * neighbors[x + 1][y + 1][z + 1];
388 
389 #define ADD_WEIGHTEDZ(x, y, z) \
390   weight = w[0][(x) + 1] * w[1][(y) + 1] * dw[2][(z) + 1]; \
391   result += weight * neighbors[x + 1][y + 1][z + 1];
392 
393 //////////////////////////////////////////////////////////////////////////////////////////
394 // compute all derivatives in at once
395 //////////////////////////////////////////////////////////////////////////////////////////
WNoiseVec(const Vec3 & p,Real * data)396 inline Vec3 WaveletNoiseField::WNoiseVec(const Vec3 &p, Real *data)
397 {
398   Vec3 final(0.);
399   Real w[3][3];
400   Real dw[3][3];
401   Real result = 0;
402   int xC, yC, zC;
403   Real weight;
404 
405   int midX = (int)ceil(p[0] - 0.5f);
406   int midY = (int)ceil(p[1] - 0.5f);
407   int midZ = (int)ceil(p[2] - 0.5f);
408 
409   Real t0 = midX - (p[0] - 0.5f);
410   Real t1 = midY - (p[1] - 0.5f);
411   Real t2 = midZ - (p[2] - 0.5f);
412 
413   // precache all the neighbors for fast access
414   Real neighbors[3][3][3];
415   for (int z = -1; z <= 1; z++)
416     for (int y = -1; y <= 1; y++)
417       for (int x = -1; x <= 1; x++) {
418         xC = modFast128(midX + (x));
419         yC = modFast128(midY + (y));
420         zC = modFast128(midZ + (z));
421         neighbors[x + 1][y + 1][z + 1] =
422             data[zC * NOISE_TILE_SIZE * NOISE_TILE_SIZE + yC * NOISE_TILE_SIZE + xC];
423       }
424 
425   ///////////////////////////////////////////////////////////////////////////////////////
426   // evaluate splines
427   ///////////////////////////////////////////////////////////////////////////////////////
428   dw[0][0] = -t0;
429   dw[0][2] = (1.f - t0);
430   dw[0][1] = 2.0f * t0 - 1.0f;
431 
432   dw[1][0] = -t1;
433   dw[1][2] = (1.0f - t1);
434   dw[1][1] = 2.0f * t1 - 1.0f;
435 
436   dw[2][0] = -t2;
437   dw[2][2] = (1.0f - t2);
438   dw[2][1] = 2.0f * t2 - 1.0f;
439 
440   w[0][0] = t0 * t0 * 0.5f;
441   w[0][2] = (1.f - t0) * (1.f - t0) * 0.5f;
442   w[0][1] = 1.f - w[0][0] - w[0][2];
443 
444   w[1][0] = t1 * t1 * 0.5f;
445   w[1][2] = (1.f - t1) * (1.f - t1) * 0.5f;
446   w[1][1] = 1.f - w[1][0] - w[1][2];
447 
448   w[2][0] = t2 * t2 * 0.5f;
449   w[2][2] = (1.f - t2) * (1.f - t2) * 0.5f;
450   w[2][1] = 1.f - w[2][0] - w[2][2];
451 
452   ///////////////////////////////////////////////////////////////////////////////////////
453   // x derivative
454   ///////////////////////////////////////////////////////////////////////////////////////
455   result = 0.0f;
456   ADD_WEIGHTEDX(-1, -1, -1);
457   ADD_WEIGHTEDX(0, -1, -1);
458   ADD_WEIGHTEDX(1, -1, -1);
459   ADD_WEIGHTEDX(-1, 0, -1);
460   ADD_WEIGHTEDX(0, 0, -1);
461   ADD_WEIGHTEDX(1, 0, -1);
462   ADD_WEIGHTEDX(-1, 1, -1);
463   ADD_WEIGHTEDX(0, 1, -1);
464   ADD_WEIGHTEDX(1, 1, -1);
465 
466   ADD_WEIGHTEDX(-1, -1, 0);
467   ADD_WEIGHTEDX(0, -1, 0);
468   ADD_WEIGHTEDX(1, -1, 0);
469   ADD_WEIGHTEDX(-1, 0, 0);
470   ADD_WEIGHTEDX(0, 0, 0);
471   ADD_WEIGHTEDX(1, 0, 0);
472   ADD_WEIGHTEDX(-1, 1, 0);
473   ADD_WEIGHTEDX(0, 1, 0);
474   ADD_WEIGHTEDX(1, 1, 0);
475 
476   ADD_WEIGHTEDX(-1, -1, 1);
477   ADD_WEIGHTEDX(0, -1, 1);
478   ADD_WEIGHTEDX(1, -1, 1);
479   ADD_WEIGHTEDX(-1, 0, 1);
480   ADD_WEIGHTEDX(0, 0, 1);
481   ADD_WEIGHTEDX(1, 0, 1);
482   ADD_WEIGHTEDX(-1, 1, 1);
483   ADD_WEIGHTEDX(0, 1, 1);
484   ADD_WEIGHTEDX(1, 1, 1);
485   final[0] = result;
486 
487   ///////////////////////////////////////////////////////////////////////////////////////
488   // y derivative
489   ///////////////////////////////////////////////////////////////////////////////////////
490   result = 0.0f;
491   ADD_WEIGHTEDY(-1, -1, -1);
492   ADD_WEIGHTEDY(0, -1, -1);
493   ADD_WEIGHTEDY(1, -1, -1);
494   ADD_WEIGHTEDY(-1, 0, -1);
495   ADD_WEIGHTEDY(0, 0, -1);
496   ADD_WEIGHTEDY(1, 0, -1);
497   ADD_WEIGHTEDY(-1, 1, -1);
498   ADD_WEIGHTEDY(0, 1, -1);
499   ADD_WEIGHTEDY(1, 1, -1);
500 
501   ADD_WEIGHTEDY(-1, -1, 0);
502   ADD_WEIGHTEDY(0, -1, 0);
503   ADD_WEIGHTEDY(1, -1, 0);
504   ADD_WEIGHTEDY(-1, 0, 0);
505   ADD_WEIGHTEDY(0, 0, 0);
506   ADD_WEIGHTEDY(1, 0, 0);
507   ADD_WEIGHTEDY(-1, 1, 0);
508   ADD_WEIGHTEDY(0, 1, 0);
509   ADD_WEIGHTEDY(1, 1, 0);
510 
511   ADD_WEIGHTEDY(-1, -1, 1);
512   ADD_WEIGHTEDY(0, -1, 1);
513   ADD_WEIGHTEDY(1, -1, 1);
514   ADD_WEIGHTEDY(-1, 0, 1);
515   ADD_WEIGHTEDY(0, 0, 1);
516   ADD_WEIGHTEDY(1, 0, 1);
517   ADD_WEIGHTEDY(-1, 1, 1);
518   ADD_WEIGHTEDY(0, 1, 1);
519   ADD_WEIGHTEDY(1, 1, 1);
520   final[1] = result;
521 
522   ///////////////////////////////////////////////////////////////////////////////////////
523   // z derivative
524   ///////////////////////////////////////////////////////////////////////////////////////
525   result = 0.0f;
526   ADD_WEIGHTEDZ(-1, -1, -1);
527   ADD_WEIGHTEDZ(0, -1, -1);
528   ADD_WEIGHTEDZ(1, -1, -1);
529   ADD_WEIGHTEDZ(-1, 0, -1);
530   ADD_WEIGHTEDZ(0, 0, -1);
531   ADD_WEIGHTEDZ(1, 0, -1);
532   ADD_WEIGHTEDZ(-1, 1, -1);
533   ADD_WEIGHTEDZ(0, 1, -1);
534   ADD_WEIGHTEDZ(1, 1, -1);
535 
536   ADD_WEIGHTEDZ(-1, -1, 0);
537   ADD_WEIGHTEDZ(0, -1, 0);
538   ADD_WEIGHTEDZ(1, -1, 0);
539   ADD_WEIGHTEDZ(-1, 0, 0);
540   ADD_WEIGHTEDZ(0, 0, 0);
541   ADD_WEIGHTEDZ(1, 0, 0);
542   ADD_WEIGHTEDZ(-1, 1, 0);
543   ADD_WEIGHTEDZ(0, 1, 0);
544   ADD_WEIGHTEDZ(1, 1, 0);
545 
546   ADD_WEIGHTEDZ(-1, -1, 1);
547   ADD_WEIGHTEDZ(0, -1, 1);
548   ADD_WEIGHTEDZ(1, -1, 1);
549   ADD_WEIGHTEDZ(-1, 0, 1);
550   ADD_WEIGHTEDZ(0, 0, 1);
551   ADD_WEIGHTEDZ(1, 0, 1);
552   ADD_WEIGHTEDZ(-1, 1, 1);
553   ADD_WEIGHTEDZ(0, 1, 1);
554   ADD_WEIGHTEDZ(1, 1, 1);
555   final[2] = result;
556 
557   // debMsg("FINAL","at "<<p<<" = "<<final); // DEBUG
558   return final;
559 }
560 #undef ADD_WEIGHTEDX
561 #undef ADD_WEIGHTEDY
562 #undef ADD_WEIGHTEDZ
563 
evaluate(Vec3 pos,int tile)564 inline Real WaveletNoiseField::evaluate(Vec3 pos, int tile) const
565 {
566   pos[0] *= mGsInvX;
567   pos[1] *= mGsInvY;
568   pos[2] *= mGsInvZ;
569   pos += mSeedOffset;
570 
571   // time anim
572   pos += Vec3(getTime());
573 
574   pos[0] *= mPosScale[0];
575   pos[1] *= mPosScale[1];
576   pos[2] *= mPosScale[2];
577   pos += mPosOffset;
578 
579   const int n3 = square(NOISE_TILE_SIZE) * NOISE_TILE_SIZE;
580   Real v = WNoise(pos, &mNoiseTile[tile * n3]);
581 
582   v += mValOffset;
583   v *= mValScale;
584   if (mClamp) {
585     if (v < mClampNeg)
586       v = mClampNeg;
587     if (v > mClampPos)
588       v = mClampPos;
589   }
590   return v;
591 }
592 
evaluateVec(Vec3 pos,int tile)593 inline Vec3 WaveletNoiseField::evaluateVec(Vec3 pos, int tile) const
594 {
595   pos[0] *= mGsInvX;
596   pos[1] *= mGsInvY;
597   pos[2] *= mGsInvZ;
598   pos += mSeedOffset;
599 
600   // time anim
601   pos += Vec3(getTime());
602 
603   pos[0] *= mPosScale[0];
604   pos[1] *= mPosScale[1];
605   pos[2] *= mPosScale[2];
606   pos += mPosOffset;
607 
608   const int n3 = square(NOISE_TILE_SIZE) * NOISE_TILE_SIZE;
609   Vec3 v = WNoiseVec(pos, &mNoiseTile[tile * n3]);
610 
611   v += Vec3(mValOffset);
612   v *= mValScale;
613 
614   if (mClamp) {
615     for (int i = 0; i < 3; i++) {
616       if (v[i] < mClampNeg)
617         v[i] = mClampNeg;
618       if (v[i] > mClampPos)
619         v[i] = mClampPos;
620     }
621   }
622   return v;
623 }
624 
evaluateCurl(Vec3 pos)625 inline Vec3 WaveletNoiseField::evaluateCurl(Vec3 pos) const
626 {
627   // gradients of w0-w2
628   Vec3 d0 = evaluateVec(pos, 0), d1 = evaluateVec(pos, 1), d2 = evaluateVec(pos, 2);
629 
630   return Vec3(d0.y - d1.z, d2.z - d0.x, d1.x - d2.y);
631 }
632 
633 }  // namespace Manta
634 
635 #endif
636