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