493ed1642fc123c389af20204d8dcb86c30b34dc
[blender.git] / intern / smoke / intern / WAVELET_NOISE.h
1 //////////////////////////////////////////////////////////////////////
2 // This file is part of Wavelet Turbulence.
3 // 
4 // Wavelet Turbulence is free software: you can redistribute it and/or modify
5 // it under the terms of the GNU General Public License as published by
6 // the Free Software Foundation, either version 3 of the License, or
7 // (at your option) any later version.
8 // 
9 // Wavelet Turbulence is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12 // GNU General Public License for more details.
13 // 
14 // You should have received a copy of the GNU General Public License
15 // along with Wavelet Turbulence.  If not, see <http://www.gnu.org/licenses/>.
16 // 
17 // Copyright 2008 Theodore Kim and Nils Thuerey
18 // 
19 //////////////////////////////////////////////////////////////////////////////////////////
20 // Wavelet noise functions
21 //
22 // This code is based on the C code provided in the appendices of:
23 //
24 // @article{1073264,
25 //  author = {Robert L. Cook and Tony DeRose},
26 //  title = {Wavelet noise},
27 //  journal = {ACM Trans. Graph.},
28 //  volume = {24},
29 //  number = {3},
30 //  year = {2005},
31 //  issn = {0730-0301},
32 //  pages = {803--811},
33 //  doi = {http://doi.acm.org/10.1145/1073204.1073264},
34 //  publisher = {ACM},
35 //  address = {New York, NY, USA},
36 //  }
37 //
38 //////////////////////////////////////////////////////////////////////////////////////////
39
40 #ifndef WAVELET_NOISE_H
41 #define WAVELET_NOISE_H
42
43 #include <MERSENNETWISTER.h>
44
45 #define NOISE_TILE_SIZE 128
46 static const int noiseTileSize = NOISE_TILE_SIZE;
47
48 // warning - noiseTileSize has to be 128^3!
49 #define modFast128(x) ((x) & 127)
50 #define modFast64(x)  ((x) & 63)
51 #define DOWNCOEFFS 0.000334f,-0.001528f, 0.000410f, 0.003545f,-0.000938f,-0.008233f, 0.002172f, 0.019120f, \
52                   -0.005040f,-0.044412f, 0.011655f, 0.103311f,-0.025936f,-0.243780f, 0.033979f, 0.655340f, \
53                    0.655340f, 0.033979f,-0.243780f,-0.025936f, 0.103311f, 0.011655f,-0.044412f,-0.005040f, \
54                    0.019120f, 0.002172f,-0.008233f,-0.000938f, 0.003546f, 0.000410f,-0.001528f, 0.000334f
55
56 //////////////////////////////////////////////////////////////////////////////////////////
57 // Wavelet downsampling -- periodic boundary conditions
58 //////////////////////////////////////////////////////////////////////////////////////////
59 static void downsampleX(float *from, float *to, int n){
60   // if these values are not local incorrect results are generated
61   float downCoeffs[32] = { DOWNCOEFFS };
62         const float *a = &downCoeffs[16];
63         for (int i = 0; i < n / 2; i++) {
64                 to[i] = 0;
65                 for (int k = 2 * i - 16; k <= 2 * i + 16; k++)
66                         to[i] += a[k - 2 * i] * from[modFast128(k)];
67         }
68 }
69 static void downsampleY(float *from, float *to, int n){
70   // if these values are not local incorrect results are generated
71   float downCoeffs[32] = { DOWNCOEFFS };
72         const float *a = &downCoeffs[16];
73         for (int i = 0; i < n / 2; i++) {
74                 to[i * n] = 0;
75                 for (int k = 2 * i - 16; k <= 2 * i + 16; k++)
76                         to[i * n] += a[k - 2 * i] * from[modFast128(k) * n];
77         }
78 }
79 static void downsampleZ(float *from, float *to, int n){
80   // if these values are not local incorrect results are generated
81   float downCoeffs[32] = { DOWNCOEFFS };
82         const float *a = &downCoeffs[16];
83         for (int i = 0; i < n / 2; i++) {
84                 to[i * n * n] = 0;
85                 for (int k = 2 * i - 16; k <= 2 * i + 16; k++)
86                         to[i * n * n] += a[k - 2 * i] * from[modFast128(k) * n * n];
87         }
88 }
89
90 //////////////////////////////////////////////////////////////////////////////////////////
91 // Wavelet downsampling -- Neumann boundary conditions
92 //////////////////////////////////////////////////////////////////////////////////////////
93 static void downsampleNeumann(const float *from, float *to, int n, int stride)
94 {
95   // if these values are not local incorrect results are generated
96   float downCoeffs[32] = { DOWNCOEFFS };
97   static const float *const aCoCenter= &downCoeffs[16];
98         for (int i = 0; i < n / 2; i++) {
99                 to[i * stride] = 0;
100                 for (int k = 2 * i - 16; k < 2 * i + 16; k++) { 
101                         // handle boundary
102                         float fromval; 
103                         if (k < 0) {
104                                 fromval = from[0];
105                         } else if(k > n - 1) {
106                                 fromval = from[(n - 1) * stride];
107                         } else {
108                                 fromval = from[k * stride]; 
109                         } 
110                         to[i * stride] += aCoCenter[k - 2 * i] * fromval; 
111                 }
112         }
113 }
114 static void downsampleXNeumann(float* to, const float* from, int sx,int sy, int sz) {
115         for (int iy = 0; iy < sy; iy++) 
116                 for (int iz = 0; iz < sz; iz++) {
117                         const int i = iy * sx + iz*sx*sy;
118                         downsampleNeumann(&from[i], &to[i], sx, 1);
119                 }
120 }
121 static void downsampleYNeumann(float* to, const float* from, int sx,int sy, int sz) {
122         for (int ix = 0; ix < sx; ix++) 
123                 for (int iz = 0; iz < sz; iz++) {
124                         const int i = ix + iz*sx*sy;
125                         downsampleNeumann(&from[i], &to[i], sy, sx);
126     }
127 }
128 static void downsampleZNeumann(float* to, const float* from, int sx,int sy, int sz) {
129         for (int ix = 0; ix < sx; ix++) 
130                 for (int iy = 0; iy < sy; iy++) {
131                         const int i = ix + iy*sx;
132                         downsampleNeumann(&from[i], &to[i], sz, sx*sy);
133     }
134 }
135
136 //////////////////////////////////////////////////////////////////////////////////////////
137 // Wavelet upsampling - periodic boundary conditions
138 //////////////////////////////////////////////////////////////////////////////////////////
139 static float _upCoeffs[4] = {0.25f, 0.75f, 0.75f, 0.25f};
140 static void upsampleX(float *from, float *to, int n) {
141         const float *p = &_upCoeffs[2];
142
143         for (int i = 0; i < n; i++) {
144                 to[i] = 0;
145                 for (int k = i / 2; k <= i / 2 + 1; k++)
146                         to[i] += p[i - 2 * k] * from[modFast64(k)];
147         }
148 }
149 static void upsampleY(float *from, float *to, int n) {
150         const float *p = &_upCoeffs[2];
151
152         for (int i = 0; i < n; i++) {
153                 to[i * n] = 0;
154                 for (int k = i / 2; k <= i / 2 + 1; k++)
155                         to[i * n] += p[i - 2 * k] * from[modFast64(k) * n];
156         }
157 }
158 static void upsampleZ(float *from, float *to, int n) {
159         const float *p = &_upCoeffs[2];
160
161         for (int i = 0; i < n; i++) {
162                 to[i * n * n] = 0;
163                 for (int k = i / 2; k <= i / 2 + 1; k++)
164                         to[i * n * n] += p[i - 2 * k] * from[modFast64(k) * n * n];
165         }
166 }
167
168 //////////////////////////////////////////////////////////////////////////////////////////
169 // Wavelet upsampling - Neumann boundary conditions
170 //////////////////////////////////////////////////////////////////////////////////////////
171 static void upsampleNeumann(const float *from, float *to, int n, int stride) {
172   static const float *const pCoCenter = &_upCoeffs[2];
173         for (int i = 0; i < n; i++) {
174                 to[i * stride] = 0;
175                 for (int k = i / 2; k <= i / 2 + 1; k++) {
176                         float fromval;
177                         if(k>n/2) {
178                                 fromval = from[(n/2) * stride];
179                         } else {
180                                 fromval = from[k * stride]; 
181                         }  
182                         to[i * stride] += pCoCenter[i - 2 * k] * fromval; 
183                 }
184         }
185 }
186 static void upsampleXNeumann(float* to, const float* from, int sx, int sy, int sz) {
187         for (int iy = 0; iy < sy; iy++) 
188                 for (int iz = 0; iz < sz; iz++) {
189                         const int i = iy * sx + iz*sx*sy;
190                         upsampleNeumann(&from[i], &to[i], sx, 1);
191                 }
192 }
193 static void upsampleYNeumann(float* to, const float* from, int sx, int sy, int sz) {
194         for (int ix = 0; ix < sx; ix++) 
195                 for (int iz = 0; iz < sz; iz++) {
196                         const int i = ix + iz*sx*sy;
197                         upsampleNeumann(&from[i], &to[i], sy, sx);
198                 }
199 }
200 static void upsampleZNeumann(float* to, const float* from, int sx, int sy, int sz) {
201         for (int ix = 0; ix < sx; ix++) 
202                 for (int iy = 0; iy < sy; iy++) {
203                         const int i = ix + iy*sx;
204                         upsampleNeumann(&from[i], &to[i], sz, sx*sy);
205                 }
206 }
207
208
209 //////////////////////////////////////////////////////////////////////////////////////////
210 // load in an existing noise tile
211 //////////////////////////////////////////////////////////////////////////////////////////
212 static bool loadTile(float* const noiseTileData, std::string filename)
213 {
214         FILE* file;
215         file = fopen(filename.c_str(), "rb");
216
217         if (file == NULL) {
218                 printf("loadTile: No noise tile '%s' found.\n", filename.c_str());
219                 return false;
220         }
221
222         // dimensions
223         int gridSize = noiseTileSize * noiseTileSize * noiseTileSize;
224
225         // noiseTileData memory is managed by caller
226         size_t bread = fread((void*)noiseTileData, sizeof(float), gridSize, file);
227         fclose(file);
228         printf("Noise tile file '%s' loaded.\n", filename.c_str());
229
230         if (bread != gridSize) {
231                 printf("loadTile: Noise tile '%s' is wrong size %d.\n", filename.c_str(), bread);
232                 return false;
233         } 
234         return true;
235 }
236
237 //////////////////////////////////////////////////////////////////////////////////////////
238 // write out an existing noise tile
239 //////////////////////////////////////////////////////////////////////////////////////////
240 static void saveTile(float* const noiseTileData, std::string filename)
241 {
242         FILE* file;
243         file = fopen(filename.c_str(), "wb");
244
245         if (file == NULL) {
246                 printf("saveTile: Noise tile '%s' could not be saved.\n", filename.c_str());
247                 return;
248         } 
249
250         fwrite((void*)noiseTileData, sizeof(float), noiseTileSize * noiseTileSize * noiseTileSize, file);
251         fclose(file);
252
253         printf("saveTile: Noise tile file '%s' saved.\n", filename.c_str());
254 }
255
256 //////////////////////////////////////////////////////////////////////////////////////////
257 // create a new noise tile if necessary
258 //////////////////////////////////////////////////////////////////////////////////////////
259 static void generateTile_WAVELET(float* const noiseTileData, std::string filename) {
260         // if a tile already exists, just use that
261         if (loadTile(noiseTileData, filename)) return;
262
263         const int n = noiseTileSize;
264         const int n3 = n*n*n;
265         std::cout <<"Generating new 3d noise tile size="<<n<<"^3 \n";
266         MTRand twister;
267
268         float *temp13 = new float[n3];
269         float *temp23 = new float[n3];
270         float *noise3 = new float[n3];
271
272         // initialize
273         for (int i = 0; i < n3; i++) {
274                 temp13[i] = temp23[i] = noise3[i] = 0.;
275         }
276
277         // Step 1. Fill the tile with random numbers in the range -1 to 1.
278         for (int i = 0; i < n3; i++) 
279                 noise3[i] = twister.randNorm();
280
281         // Steps 2 and 3. Downsample and upsample the tile
282         for (int iy = 0; iy < n; iy++) 
283                 for (int iz = 0; iz < n; iz++) {
284                         const int i = iy * n + iz*n*n;
285                         downsampleX(&noise3[i], &temp13[i], n);
286                         upsampleX  (&temp13[i], &temp23[i], n);
287                 }
288         for (int ix = 0; ix < n; ix++) 
289                 for (int iz = 0; iz < n; iz++) {
290                         const int i = ix + iz*n*n;
291                         downsampleY(&temp23[i], &temp13[i], n);
292                         upsampleY  (&temp13[i], &temp23[i], n);
293                 }
294         for (int ix = 0; ix < n; ix++) 
295                 for (int iy = 0; iy < n; iy++) {
296                         const int i = ix + iy*n;
297                         downsampleZ(&temp23[i], &temp13[i], n);
298                         upsampleZ  (&temp13[i], &temp23[i], n);
299                 }
300
301         // Step 4. Subtract out the coarse-scale contribution
302         for (int i = 0; i < n3; i++) 
303                 noise3[i] -= temp23[i];
304
305         // Avoid even/odd variance difference by adding odd-offset version of noise to itself.
306         int offset = n / 2;
307         if (offset % 2 == 0) offset++;
308
309         int icnt=0;
310         for (int ix = 0; ix < n; ix++)
311                 for (int iy = 0; iy < n; iy++)
312                         for (int iz = 0; iz < n; iz++) { 
313                                 temp13[icnt] = noise3[modFast128(ix+offset) + modFast128(iy+offset)*n + modFast128(iz+offset)*n*n];
314                                 icnt++;
315                         }
316
317         for (int i = 0; i < n3; i++) 
318                 noise3[i] += temp13[i];
319
320         for (int i = 0; i < n3; i++) 
321                 noiseTileData[i] = noise3[i];
322
323         saveTile(noise3, filename); 
324         delete[] temp13;
325         delete[] temp23;
326         delete[] noise3;
327         std::cout <<"Generating new 3d noise done\n";
328 }
329
330 //////////////////////////////////////////////////////////////////////////////////////////
331 // x derivative of noise
332 //////////////////////////////////////////////////////////////////////////////////////////
333 static inline float WNoiseDx(Vec3 p, float* data) { 
334   int c[3], mid[3], n = noiseTileSize;
335   float w[3][3], t, result = 0;
336   
337   mid[0] = (int)ceil(p[0] - 0.5); 
338   t = mid[0] - (p[0] - 0.5);
339         w[0][0] = -t;
340         w[0][2] = (1.f - t);
341         w[0][1] = 2.0f * t - 1.0f;
342   
343   mid[1] = (int)ceil(p[1] - 0.5); 
344   t = mid[1] - (p[1] - 0.5);
345   w[1][0] = t * t / 2; 
346   w[1][2] = (1 - t) * (1 - t) / 2;
347   w[1][1] = 1 - w[1][0] - w[1][2];
348
349   mid[2] = (int)ceil(p[2] - 0.5); 
350   t = mid[2] - (p[2] - 0.5);
351   w[2][0] = t * t / 2; 
352   w[2][2] = (1 - t) * (1 - t)/2; 
353   w[2][1] = 1 - w[2][0] - w[2][2];
354  
355   // to optimize, explicitly unroll this loop
356   for (int z = -1; z <=1; z++)
357     for (int y = -1; y <=1; y++)
358       for (int x = -1; x <=1; x++)
359       {
360         float weight = 1.0f;
361         c[0] = modFast128(mid[0] + x);
362         weight *= w[0][x+1];
363         c[1] = modFast128(mid[1] + y);
364         weight *= w[1][y+1];
365         c[2] = modFast128(mid[2] + z);
366         weight *= w[2][z+1];
367         result += weight * data[c[2]*n*n+c[1]*n+c[0]];
368       }
369  return result;
370 }
371
372 //////////////////////////////////////////////////////////////////////////////////////////
373 // y derivative of noise
374 //////////////////////////////////////////////////////////////////////////////////////////
375 static inline float WNoiseDy(Vec3 p, float* data) { 
376   int c[3], mid[3], n=noiseTileSize; 
377   float w[3][3], t, result =0;
378   
379   mid[0] = (int)ceil(p[0] - 0.5); 
380   t = mid[0]-(p[0] - 0.5);
381   w[0][0] = t * t / 2; 
382   w[0][2] = (1 - t) * (1 - t) / 2;
383   w[0][1] = 1 - w[0][0] - w[0][2];
384   
385   mid[1] = (int)ceil(p[1] - 0.5); 
386   t = mid[1]-(p[1] - 0.5);
387         w[1][0] = -t;
388         w[1][2] = (1.f - t);
389         w[1][1] = 2.0f * t - 1.0f;
390
391   mid[2] = (int)ceil(p[2] - 0.5); 
392   t = mid[2] - (p[2] - 0.5);
393   w[2][0] = t * t / 2; 
394   w[2][2] = (1 - t) * (1 - t)/2; 
395   w[2][1] = 1 - w[2][0] - w[2][2];
396   
397   // to optimize, explicitly unroll this loop
398   for (int z = -1; z <=1; z++)
399     for (int y = -1; y <=1; y++)
400       for (int x = -1; x <=1; x++)
401       {
402         float weight = 1.0f;
403         c[0] = modFast128(mid[0] + x);
404         weight *= w[0][x+1];
405         c[1] = modFast128(mid[1] + y);
406         weight *= w[1][y+1];
407         c[2] = modFast128(mid[2] + z);
408         weight *= w[2][z+1];
409         result += weight * data[c[2]*n*n+c[1]*n+c[0]];
410       }
411
412   return result;
413 }
414
415 //////////////////////////////////////////////////////////////////////////////////////////
416 // z derivative of noise
417 //////////////////////////////////////////////////////////////////////////////////////////
418 static inline float WNoiseDz(Vec3 p, float* data) { 
419   int c[3], mid[3], n=noiseTileSize; 
420   float w[3][3], t, result =0;
421
422   mid[0] = (int)ceil(p[0] - 0.5); 
423   t = mid[0]-(p[0] - 0.5);
424   w[0][0] = t * t / 2; 
425   w[0][2] = (1 - t) * (1 - t) / 2;
426   w[0][1] = 1 - w[0][0] - w[0][2];
427   
428   mid[1] = (int)ceil(p[1] - 0.5); 
429   t = mid[1]-(p[1] - 0.5);
430   w[1][0] = t * t / 2; 
431   w[1][2] = (1 - t) * (1 - t) / 2;
432   w[1][1] = 1 - w[1][0] - w[1][2];
433
434   mid[2] = (int)ceil(p[2] - 0.5); 
435   t = mid[2] - (p[2] - 0.5);
436         w[2][0] = -t;
437         w[2][2] = (1.f - t);
438         w[2][1] = 2.0f * t - 1.0f;
439
440   // to optimize, explicitly unroll this loop
441   for (int z = -1; z <=1; z++)
442     for (int y = -1; y <=1; y++)
443       for (int x = -1; x <=1; x++)
444       {
445         float weight = 1.0f;
446         c[0] = modFast128(mid[0] + x);
447         weight *= w[0][x+1];
448         c[1] = modFast128(mid[1] + y);
449         weight *= w[1][y+1];
450         c[2] = modFast128(mid[2] + z);
451         weight *= w[2][z+1];
452         result += weight * data[c[2]*n*n+c[1]*n+c[0]];
453       }
454   return result;
455 }
456
457 #endif
458