Smoke: decoupling of wavelet #2, new noise strength option on gui, fftw3 option in...
[blender.git] / intern / smoke / intern / FFT_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 //
21
22 #ifndef FFT_NOISE_H_
23 #define FFT_NOISE_H_
24
25 #if FFTW3==1
26 #include <iostream>
27 #include <fftw3.h>
28 #include <MERSENNETWISTER.h>
29
30 #include "WAVELET_NOISE.h"
31
32 #ifndef M_PI
33 #define M_PI 3.14159265
34 #endif
35
36 /////////////////////////////////////////////////////////////////////////
37 // shift spectrum to the format that FFTW expects
38 /////////////////////////////////////////////////////////////////////////
39 static void shift3D(float*& field, int xRes, int yRes, int zRes)
40 {
41   int xHalf = xRes / 2;
42   int yHalf = yRes / 2;
43   int zHalf = zRes / 2;
44   int slabSize = xRes * yRes;
45   for (int z = 0; z < zHalf; z++)
46     for (int y = 0; y < yHalf; y++)
47       for (int x = 0; x < xHalf; x++)
48       {
49         int index = x + y * xRes + z * xRes * yRes;
50         float temp;
51         int xSwap = xHalf;
52         int ySwap = yHalf * xRes;
53         int zSwap = zHalf * xRes * yRes;
54         
55         // [0,0,0] to [1,1,1]
56         temp = field[index];
57         field[index] = field[index + xSwap + ySwap + zSwap];
58         field[index + xSwap + ySwap + zSwap] = temp;
59
60         // [1,0,0] to [0,1,1]
61         temp = field[index + xSwap];
62         field[index + xSwap] = field[index + ySwap + zSwap];
63         field[index + ySwap + zSwap] = temp;
64
65         // [0,1,0] to [1,0,1]
66         temp = field[index + ySwap];
67         field[index + ySwap] = field[index + xSwap + zSwap];
68         field[index + xSwap + zSwap] = temp;
69         
70         // [0,0,1] to [1,1,0]
71         temp = field[index + zSwap];
72         field[index + zSwap] = field[index + xSwap + ySwap];
73         field[index + xSwap + ySwap] = temp;
74       }
75 }
76
77 static void generatTile_FFT(float* const noiseTileData, std::string filename)
78 {
79         if (loadTile(noiseTileData, filename)) return;
80         
81         int res = NOISE_TILE_SIZE;
82         int xRes = res;
83         int yRes = res;
84         int zRes = res;
85         int totalCells = xRes * yRes * zRes;
86         
87         // create and shift the filter
88         float* filter = new float[totalCells];
89         for (int z = 0; z < zRes; z++)
90                 for (int y = 0; y < yRes; y++)
91                         for (int x = 0; x < xRes; x++)
92                         {
93                                 int index = x + y * xRes + z * xRes * yRes;
94                                 float diff[] = {abs(x - xRes/2), 
95                                 abs(y - yRes/2), 
96                                 abs(z - zRes/2)};
97                                 float radius = sqrtf(diff[0] * diff[0] + 
98                                 diff[1] * diff[1] + 
99                                 diff[2] * diff[2]) / (xRes / 2);
100                                 radius *= M_PI;
101                                 float H = cos((M_PI / 2.0f) * log(4.0f * radius / M_PI) / log(2.0f));
102                                 H = H * H;
103                                 float filtered = H;
104                                 
105                                 // clamp everything outside the wanted band
106                                 if (radius >= M_PI / 2.0f)
107                                         filtered = 0.0f;
108                                 
109                                 // make sure to capture all low frequencies
110                                 if (radius <= M_PI / 4.0f)
111                                         filtered = 1.0f;
112                                 
113                                 filter[index] = filtered;
114                         }
115         shift3D(filter, xRes, yRes, zRes);
116         
117         // create the noise
118         float* noise = new float[totalCells];
119         int index = 0;
120         MTRand twister;
121         for (int z = 0; z < zRes; z++)
122         for (int y = 0; y < yRes; y++)
123           for (int x = 0; x < xRes; x++, index++)
124                 noise[index] = twister.randNorm();
125         
126         // create padded field
127         fftw_complex* forward = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * totalCells);
128         
129         // init padded field
130         index = 0;
131         for (int z = 0; z < zRes; z++)
132         for (int y = 0; y < yRes; y++)
133           for (int x = 0; x < xRes; x++, index++)
134           {
135                 forward[index][0] = noise[index];
136                 forward[index][1] = 0.0f;
137           }
138         
139         // forward FFT 
140         fftw_complex* backward = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * totalCells);
141         fftw_plan forwardPlan = fftw_plan_dft_3d(xRes, yRes, zRes, forward, backward, FFTW_FORWARD, FFTW_ESTIMATE);  
142         fftw_execute(forwardPlan);
143         fftw_destroy_plan(forwardPlan);
144         
145         // apply filter
146         index = 0;
147         for (int z = 0; z < zRes; z++)
148                 for (int y = 0; y < yRes; y++)
149                         for (int x = 0; x < xRes; x++, index++)
150                         {
151                                 backward[index][0] *= filter[index];
152                                 backward[index][1] *= filter[index];
153                         }
154         
155         // backward FFT
156         fftw_plan backwardPlan = fftw_plan_dft_3d(xRes, yRes, zRes, backward, forward, FFTW_BACKWARD, FFTW_ESTIMATE);  
157         fftw_execute(backwardPlan);
158         fftw_destroy_plan(backwardPlan);
159         
160         // subtract out the low frequency components
161         index = 0;
162         for (int z = 0; z < zRes; z++)
163                 for (int y = 0; y < yRes; y++)
164                         for (int x = 0; x < xRes; x++, index++)
165                                 noise[index] -= forward[index][0] / totalCells;
166
167         // save out the noise tile
168         saveTile(noise, filename);
169         
170         fftw_free(forward);
171         fftw_free(backward);
172         delete[] filter;
173         delete[] noise;
174 }
175
176 #endif
177
178 #endif /* FFT_NOISE_H_ */