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