Fix #27855: crash on enabling high resolution smoke.
[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
266         // check for invalid nan tile data that could be generated. bug is now
267         // fixed, but invalid files may still hang around
268         if (isnan(noiseTileData[0])) {
269                 printf("loadTile: Noise tile '%s' contains nan values.\n", filename.c_str());
270                 return false;
271         }
272
273         return true;
274 }
275
276 //////////////////////////////////////////////////////////////////////////////////////////
277 // write out an existing noise tile
278 //////////////////////////////////////////////////////////////////////////////////////////
279 static void saveTile(float* const noiseTileData, std::string filename)
280 {
281         FILE* file;
282         file = fopen(filename.c_str(), "wb");
283         int endiantest=1;
284         char longsize;
285
286         if (file == NULL) {
287                 printf("saveTile: Noise tile '%s' could not be saved.\n", filename.c_str());
288                 return;
289         } 
290
291         //Write file header
292         fwrite(tilefile_headerstring, strlen(tilefile_headerstring), 1, file);
293         fwrite(tilefile_revision, strlen(tilefile_revision), 1, file);
294         //Endianness
295         if (*((unsigned char*)&endiantest) == 1)
296                 fwrite("L", 1, 1, file); //Little endian
297         else
298                 fwrite("B",1,1,file); //Big endian
299         //32/64bit
300         longsize = (char)sizeof(long)+'0';
301         fwrite(&longsize, 1, 1, file);
302         
303         
304         fwrite((void*)noiseTileData, sizeof(float), noiseTileSize * noiseTileSize * noiseTileSize, file);
305         fclose(file);
306
307         printf("saveTile: Noise tile file '%s' saved.\n", filename.c_str());
308 }
309
310 //////////////////////////////////////////////////////////////////////////////////////////
311 // create a new noise tile if necessary
312 //////////////////////////////////////////////////////////////////////////////////////////
313 static void generateTile_WAVELET(float* const noiseTileData, std::string filename) {
314         // if a tile already exists, just use that
315         if (loadTile(noiseTileData, filename)) return;
316
317         const int n = noiseTileSize;
318         const int n3 = n*n*n;
319         std::cout <<"Generating new 3d noise tile size="<<n<<"^3 \n";
320         MTRand twister;
321
322         float *temp13 = new float[n3];
323         float *temp23 = new float[n3];
324         float *noise3 = new float[n3];
325
326         // initialize
327         for (int i = 0; i < n3; i++) {
328                 temp13[i] = temp23[i] = noise3[i] = 0.;
329         }
330
331         // Step 1. Fill the tile with random numbers in the range -1 to 1.
332         for (int i = 0; i < n3; i++) 
333                 noise3[i] = twister.randNorm();
334
335         // Steps 2 and 3. Downsample and upsample the tile
336         for (int iy = 0; iy < n; iy++) 
337                 for (int iz = 0; iz < n; iz++) {
338                         const int i = iy * n + iz*n*n;
339                         downsampleX(&noise3[i], &temp13[i], n);
340                         upsampleX  (&temp13[i], &temp23[i], n);
341                 }
342         for (int ix = 0; ix < n; ix++) 
343                 for (int iz = 0; iz < n; iz++) {
344                         const int i = ix + iz*n*n;
345                         downsampleY(&temp23[i], &temp13[i], n);
346                         upsampleY  (&temp13[i], &temp23[i], n);
347                 }
348         for (int ix = 0; ix < n; ix++) 
349                 for (int iy = 0; iy < n; iy++) {
350                         const int i = ix + iy*n;
351                         downsampleZ(&temp23[i], &temp13[i], n);
352                         upsampleZ  (&temp13[i], &temp23[i], n);
353                 }
354
355         // Step 4. Subtract out the coarse-scale contribution
356         for (int i = 0; i < n3; i++) 
357                 noise3[i] -= temp23[i];
358
359         // Avoid even/odd variance difference by adding odd-offset version of noise to itself.
360         int offset = n / 2;
361         if (offset % 2 == 0) offset++;
362
363         int icnt=0;
364         for (int ix = 0; ix < n; ix++)
365                 for (int iy = 0; iy < n; iy++)
366                         for (int iz = 0; iz < n; iz++) { 
367                                 temp13[icnt] = noise3[modFast128(ix+offset) + modFast128(iy+offset)*n + modFast128(iz+offset)*n*n];
368                                 icnt++;
369                         }
370
371         for (int i = 0; i < n3; i++) 
372                 noise3[i] += temp13[i];
373
374         for (int i = 0; i < n3; i++) 
375                 noiseTileData[i] = noise3[i];
376
377         saveTile(noise3, filename); 
378         delete[] temp13;
379         delete[] temp23;
380         delete[] noise3;
381         std::cout <<"Generating new 3d noise done\n";
382 }
383
384 //////////////////////////////////////////////////////////////////////////////////////////
385 // x derivative of noise
386 //////////////////////////////////////////////////////////////////////////////////////////
387 static inline float WNoiseDx(Vec3 p, float* data) { 
388   int c[3], mid[3], n = noiseTileSize;
389   float w[3][3], t, result = 0;
390   
391   mid[0] = (int)ceil(p[0] - 0.5); 
392   t = mid[0] - (p[0] - 0.5);
393         w[0][0] = -t;
394         w[0][2] = (1.f - t);
395         w[0][1] = 2.0f * t - 1.0f;
396   
397   mid[1] = (int)ceil(p[1] - 0.5); 
398   t = mid[1] - (p[1] - 0.5);
399   w[1][0] = t * t / 2; 
400   w[1][2] = (1 - t) * (1 - t) / 2;
401   w[1][1] = 1 - w[1][0] - w[1][2];
402
403   mid[2] = (int)ceil(p[2] - 0.5); 
404   t = mid[2] - (p[2] - 0.5);
405   w[2][0] = t * t / 2; 
406   w[2][2] = (1 - t) * (1 - t)/2; 
407   w[2][1] = 1 - w[2][0] - w[2][2];
408  
409   // to optimize, explicitly unroll this loop
410   for (int z = -1; z <=1; z++)
411     for (int y = -1; y <=1; y++)
412       for (int x = -1; x <=1; x++)
413       {
414         float weight = 1.0f;
415         c[0] = modFast128(mid[0] + x);
416         weight *= w[0][x+1];
417         c[1] = modFast128(mid[1] + y);
418         weight *= w[1][y+1];
419         c[2] = modFast128(mid[2] + z);
420         weight *= w[2][z+1];
421         result += weight * data[c[2]*n*n+c[1]*n+c[0]];
422       }
423  return result;
424 }
425
426 //////////////////////////////////////////////////////////////////////////////////////////
427 // y derivative of noise
428 //////////////////////////////////////////////////////////////////////////////////////////
429 static inline float WNoiseDy(Vec3 p, float* data) { 
430   int c[3], mid[3], n=noiseTileSize; 
431   float w[3][3], t, result =0;
432   
433   mid[0] = (int)ceil(p[0] - 0.5); 
434   t = mid[0]-(p[0] - 0.5);
435   w[0][0] = t * t / 2; 
436   w[0][2] = (1 - t) * (1 - t) / 2;
437   w[0][1] = 1 - w[0][0] - w[0][2];
438   
439   mid[1] = (int)ceil(p[1] - 0.5); 
440   t = mid[1]-(p[1] - 0.5);
441         w[1][0] = -t;
442         w[1][2] = (1.f - t);
443         w[1][1] = 2.0f * t - 1.0f;
444
445   mid[2] = (int)ceil(p[2] - 0.5); 
446   t = mid[2] - (p[2] - 0.5);
447   w[2][0] = t * t / 2; 
448   w[2][2] = (1 - t) * (1 - t)/2; 
449   w[2][1] = 1 - w[2][0] - w[2][2];
450   
451   // to optimize, explicitly unroll this loop
452   for (int z = -1; z <=1; z++)
453     for (int y = -1; y <=1; y++)
454       for (int x = -1; x <=1; x++)
455       {
456         float weight = 1.0f;
457         c[0] = modFast128(mid[0] + x);
458         weight *= w[0][x+1];
459         c[1] = modFast128(mid[1] + y);
460         weight *= w[1][y+1];
461         c[2] = modFast128(mid[2] + z);
462         weight *= w[2][z+1];
463         result += weight * data[c[2]*n*n+c[1]*n+c[0]];
464       }
465
466   return result;
467 }
468
469 //////////////////////////////////////////////////////////////////////////////////////////
470 // z derivative of noise
471 //////////////////////////////////////////////////////////////////////////////////////////
472 static inline float WNoiseDz(Vec3 p, float* data) { 
473   int c[3], mid[3], n=noiseTileSize; 
474   float w[3][3], t, result =0;
475
476   mid[0] = (int)ceil(p[0] - 0.5); 
477   t = mid[0]-(p[0] - 0.5);
478   w[0][0] = t * t / 2; 
479   w[0][2] = (1 - t) * (1 - t) / 2;
480   w[0][1] = 1 - w[0][0] - w[0][2];
481   
482   mid[1] = (int)ceil(p[1] - 0.5); 
483   t = mid[1]-(p[1] - 0.5);
484   w[1][0] = t * t / 2; 
485   w[1][2] = (1 - t) * (1 - t) / 2;
486   w[1][1] = 1 - w[1][0] - w[1][2];
487
488   mid[2] = (int)ceil(p[2] - 0.5); 
489   t = mid[2] - (p[2] - 0.5);
490         w[2][0] = -t;
491         w[2][2] = (1.f - t);
492         w[2][1] = 2.0f * t - 1.0f;
493
494   // to optimize, explicitly unroll this loop
495   for (int z = -1; z <=1; z++)
496     for (int y = -1; y <=1; y++)
497       for (int x = -1; x <=1; x++)
498       {
499         float weight = 1.0f;
500         c[0] = modFast128(mid[0] + x);
501         weight *= w[0][x+1];
502         c[1] = modFast128(mid[1] + y);
503         weight *= w[1][y+1];
504         c[2] = modFast128(mid[2] + z);
505         weight *= w[2][z+1];
506         result += weight * data[c[2]*n*n+c[1]*n+c[0]];
507       }
508   return result;
509 }
510
511 #endif
512