Cleanup: remove redundant doxygen \file argument
[blender.git] / intern / smoke / intern / FFT_NOISE.h
1 /** \file \ingroup smoke
2  */
3 //////////////////////////////////////////////////////////////////////
4 // This file is part of Wavelet Turbulence.
5 // 
6 // Wavelet Turbulence is free software: you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published by
8 // the Free Software Foundation, either version 3 of the License, or
9 // (at your option) any later version.
10 // 
11 // Wavelet Turbulence is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 // GNU General Public License for more details.
15 // 
16 // You should have received a copy of the GNU General Public License
17 // along with Wavelet Turbulence.  If not, see <http://www.gnu.org/licenses/>.
18 // 
19 // Copyright 2008 Theodore Kim and Nils Thuerey
20 // 
21 /////////////////////////////////////////////////////////////////////////
22 //
23
24 #ifndef FFT_NOISE_H_
25 #define FFT_NOISE_H_
26
27 #ifdef WITH_FFTW3
28 #include <iostream>
29 #include <fftw3.h>
30 #include <MERSENNETWISTER.h>
31
32 #include "WAVELET_NOISE.h"
33
34 #ifndef M_PI
35 #define M_PI 3.14159265
36 #endif
37
38 /////////////////////////////////////////////////////////////////////////
39 // shift spectrum to the format that FFTW expects
40 /////////////////////////////////////////////////////////////////////////
41 static void shift3D(float*& field, int xRes, int yRes, int zRes)
42 {
43   int xHalf = xRes / 2;
44   int yHalf = yRes / 2;
45   int zHalf = zRes / 2;
46  // int slabSize = xRes * yRes;
47   for (int z = 0; z < zHalf; z++)
48     for (int y = 0; y < yHalf; y++)
49       for (int x = 0; x < xHalf; x++)
50       {
51         int index = x + y * xRes + z * xRes * yRes;
52         float temp;
53         int xSwap = xHalf;
54         int ySwap = yHalf * xRes;
55         int zSwap = zHalf * xRes * yRes;
56         
57         // [0,0,0] to [1,1,1]
58         temp = field[index];
59         field[index] = field[index + xSwap + ySwap + zSwap];
60         field[index + xSwap + ySwap + zSwap] = temp;
61
62         // [1,0,0] to [0,1,1]
63         temp = field[index + xSwap];
64         field[index + xSwap] = field[index + ySwap + zSwap];
65         field[index + ySwap + zSwap] = temp;
66
67         // [0,1,0] to [1,0,1]
68         temp = field[index + ySwap];
69         field[index + ySwap] = field[index + xSwap + zSwap];
70         field[index + xSwap + zSwap] = temp;
71         
72         // [0,0,1] to [1,1,0]
73         temp = field[index + zSwap];
74         field[index + zSwap] = field[index + xSwap + ySwap];
75         field[index + xSwap + ySwap] = temp;
76       }
77 }
78
79 static void generatTile_FFT(float* const noiseTileData, std::string filename)
80 {
81         if (loadTile(noiseTileData, filename)) return;
82         
83         int res = NOISE_TILE_SIZE;
84         int xRes = res;
85         int yRes = res;
86         int zRes = res;
87         int totalCells = xRes * yRes * zRes;
88         
89         // create and shift the filter
90         float* filter = new float[totalCells];
91         for (int z = 0; z < zRes; z++)
92                 for (int y = 0; y < yRes; y++)
93                         for (int x = 0; x < xRes; x++)
94                         {
95                                 int index = x + y * xRes + z * xRes * yRes;
96                                 float diff[] = {(float)abs(x - xRes / 2),
97                                                 (float)abs(y - yRes / 2),
98                                                 (float)abs(z - zRes / 2)};
99                                 float radius = sqrtf(diff[0] * diff[0] + 
100                                 diff[1] * diff[1] + 
101                                 diff[2] * diff[2]) / (xRes / 2);
102                                 radius *= M_PI;
103                                 float H = cos((M_PI / 2.0f) * log(4.0f * radius / M_PI) / log(2.0f));
104                                 H = H * H;
105                                 float filtered = H;
106                                 
107                                 // clamp everything outside the wanted band
108                                 if (radius >= M_PI / 2.0f)
109                                         filtered = 0.0f;
110                                 
111                                 // make sure to capture all low frequencies
112                                 if (radius <= M_PI / 4.0f)
113                                         filtered = 1.0f;
114                                 
115                                 filter[index] = filtered;
116                         }
117         shift3D(filter, xRes, yRes, zRes);
118         
119         // create the noise
120         float* noise = new float[totalCells];
121         int index = 0;
122         MTRand twister;
123         for (int z = 0; z < zRes; z++)
124         for (int y = 0; y < yRes; y++)
125           for (int x = 0; x < xRes; x++, index++)
126                 noise[index] = twister.randNorm();
127         
128         // create padded field
129         fftw_complex* forward = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * totalCells);
130         
131         // init padded field
132         index = 0;
133         for (int z = 0; z < zRes; z++)
134         for (int y = 0; y < yRes; y++)
135           for (int x = 0; x < xRes; x++, index++)
136           {
137                 forward[index][0] = noise[index];
138                 forward[index][1] = 0.0f;
139           }
140         
141         // forward FFT 
142         fftw_complex* backward = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * totalCells);
143         fftw_plan forwardPlan = fftw_plan_dft_3d(xRes, yRes, zRes, forward, backward, FFTW_FORWARD, FFTW_ESTIMATE);  
144         fftw_execute(forwardPlan);
145         fftw_destroy_plan(forwardPlan);
146         
147         // apply filter
148         index = 0;
149         for (int z = 0; z < zRes; z++)
150                 for (int y = 0; y < yRes; y++)
151                         for (int x = 0; x < xRes; x++, index++)
152                         {
153                                 backward[index][0] *= filter[index];
154                                 backward[index][1] *= filter[index];
155                         }
156         
157         // backward FFT
158         fftw_plan backwardPlan = fftw_plan_dft_3d(xRes, yRes, zRes, backward, forward, FFTW_BACKWARD, FFTW_ESTIMATE);  
159         fftw_execute(backwardPlan);
160         fftw_destroy_plan(backwardPlan);
161         
162         // subtract out the low frequency components
163         index = 0;
164         for (int z = 0; z < zRes; z++)
165                 for (int y = 0; y < yRes; y++)
166                         for (int x = 0; x < xRes; x++, index++)
167                                 noise[index] -= forward[index][0] / totalCells;
168
169         // fill noiseTileData
170         memcpy(noiseTileData, noise, sizeof(float) * totalCells);
171         // save out the noise tile
172         saveTile(noise, filename);
173         
174         fftw_free(forward);
175         fftw_free(backward);
176         delete[] filter;
177         delete[] noise;
178 }
179
180 #endif
181
182 #endif /* FFT_NOISE_H_ */