6800c3b84b02345f6cc90ba5ba76391da468dcef
[blender.git] / intern / smoke / intern / INTERPOLATE.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 #ifndef INTERPOLATE_H
21 #define INTERPOLATE_H
22
23 #include <iostream>
24 #include <VEC3.h>
25
26 namespace INTERPOLATE {
27
28 //////////////////////////////////////////////////////////////////////
29 // linear interpolators
30 //////////////////////////////////////////////////////////////////////
31 static inline float lerp(float t, float a, float b) {
32         return ( a + t * (b - a) );
33 }
34
35 static inline float lerp(float* field, float x, float y, int res) {
36         // clamp backtrace to grid boundaries
37         if (x < 0.5f) x = 0.5f;
38         if (x > res - 1.5f) x = res - 1.5f;
39         if (y < 0.5f) y = 0.5f;
40         if (y > res - 1.5f) y = res - 1.5f;
41
42         const int x0 = (int)x;
43         const int y0 = (int)y;
44         x -= x0;
45         y -= y0;
46         float d00, d10, d01, d11;
47
48         // lerp the velocities
49         d00 = field[x0 + y0 * res];
50         d10 = field[(x0 + 1) + y0 * res];
51         d01 = field[x0 + (y0 + 1) * res];
52         d11 = field[(x0 + 1) + (y0 + 1) * res];
53         return lerp(y, lerp(x, d00, d10),
54                         lerp(x, d01, d11));
55 }
56
57 //////////////////////////////////////////////////////////////////////////////////////////
58 // 3d linear interpolation
59 ////////////////////////////////////////////////////////////////////////////////////////// 
60 static inline float lerp3d(float* field, float x, float y, float z,  int xres, int yres, int zres) {
61         // clamp pos to grid boundaries
62         if (x < 0.5) x = 0.5;
63         if (x > xres - 1.5) x = xres - 1.5;
64         if (y < 0.5) y = 0.5;
65         if (y > yres - 1.5) y = yres - 1.5;
66         if (z < 0.5) z = 0.5;
67         if (z > zres - 1.5) z = zres - 1.5;
68
69         // locate neighbors to interpolate
70         const int x0 = (int)x;
71         const int x1 = x0 + 1;
72         const int y0 = (int)y;
73         const int y1 = y0 + 1;
74         const int z0 = (int)z;
75         const int z1 = z0 + 1;
76
77         // get interpolation weights
78         const float s1 = x - (float)x0;
79         const float s0 = 1.0f - s1;
80         const float t1 = y - (float)y0;
81         const float t0 = 1.0f - t1;
82         const float u1 = z - (float)z0;
83         const float u0 = 1.0f - u1;
84
85         const int slabSize = xres*yres;
86         const int i000 = x0 + y0 * xres + z0 * slabSize;
87         const int i010 = x0 + y1 * xres + z0 * slabSize;
88         const int i100 = x1 + y0 * xres + z0 * slabSize;
89         const int i110 = x1 + y1 * xres + z0 * slabSize;
90         const int i001 = x0 + y0 * xres + z1 * slabSize;
91         const int i011 = x0 + y1 * xres + z1 * slabSize;
92         const int i101 = x1 + y0 * xres + z1 * slabSize;
93         const int i111 = x1 + y1 * xres + z1 * slabSize;
94
95         // interpolate (indices could be computed once)
96         return ( u0 * (s0 * (t0 * field[i000] +
97                 t1 * field[i010]) +
98                 s1 * (t0 * field[i100] +
99                 t1 * field[i110])) +
100                 u1 * (s0 * (t0 * field[i001] +
101                 t1 * field[i011]) +
102                 s1 * (t0 * field[i101] +
103                 t1 * field[i111])) );
104 }
105
106 //////////////////////////////////////////////////////////////////////////////////////////
107 // convert field entries of type T to floats, then interpolate
108 //////////////////////////////////////////////////////////////////////////////////////////
109 template <class T> 
110 static inline float lerp3dToFloat(T* field1,
111                 float x, float y, float z,  int xres, int yres, int zres) {
112         // clamp pos to grid boundaries
113         if (x < 0.5) x = 0.5;
114         if (x > xres - 1.5) x = xres - 1.5;
115         if (y < 0.5) y = 0.5;
116         if (y > yres - 1.5) y = yres - 1.5;
117         if (z < 0.5) z = 0.5;
118         if (z > zres - 1.5) z = zres - 1.5;
119
120         // locate neighbors to interpolate
121         const int x0 = (int)x;
122         const int x1 = x0 + 1;
123         const int y0 = (int)y;
124         const int y1 = y0 + 1;
125         const int z0 = (int)z;
126         const int z1 = z0 + 1;
127
128         // get interpolation weights
129         const float s1 = x - (float)x0;
130         const float s0 = 1.0f - s1;
131         const float t1 = y - (float)y0;
132         const float t0 = 1.0f - t1;
133         const float u1 = z - (float)z0;
134         const float u0 = 1.0f - u1;
135
136         const int slabSize = xres*yres;
137         const int i000 = x0 + y0 * xres + z0 * slabSize;
138         const int i010 = x0 + y1 * xres + z0 * slabSize;
139         const int i100 = x1 + y0 * xres + z0 * slabSize;
140         const int i110 = x1 + y1 * xres + z0 * slabSize;
141         const int i001 = x0 + y0 * xres + z1 * slabSize;
142         const int i011 = x0 + y1 * xres + z1 * slabSize;
143         const int i101 = x1 + y0 * xres + z1 * slabSize;
144         const int i111 = x1 + y1 * xres + z1 * slabSize;
145
146         // interpolate (indices could be computed once)
147         return (float)(
148                         ( u0 * (s0 * (t0 * (float)field1[i000] +
149                                 t1 * (float)field1[i010]) +
150                                 s1 * (t0 * (float)field1[i100] +
151                                 t1 * (float)field1[i110])) +
152                                 u1 * (s0 * (t0 * (float)field1[i001] +
153                                 t1 * (float)field1[i011]) +
154                                 s1 * (t0 * (float)field1[i101] +
155                                 t1 * (float)field1[i111])) ) );
156 }
157
158 //////////////////////////////////////////////////////////////////////////////////////////
159 // interpolate a vector from 3 fields
160 //////////////////////////////////////////////////////////////////////////////////////////
161 static inline Vec3 lerp3dVec(float* field1, float* field2, float* field3, 
162                 float x, float y, float z,  int xres, int yres, int zres) {
163         // clamp pos to grid boundaries
164         if (x < 0.5) x = 0.5;
165         if (x > xres - 1.5) x = xres - 1.5;
166         if (y < 0.5) y = 0.5;
167         if (y > yres - 1.5) y = yres - 1.5;
168         if (z < 0.5) z = 0.5;
169         if (z > zres - 1.5) z = zres - 1.5;
170
171         // locate neighbors to interpolate
172         const int x0 = (int)x;
173         const int x1 = x0 + 1;
174         const int y0 = (int)y;
175         const int y1 = y0 + 1;
176         const int z0 = (int)z;
177         const int z1 = z0 + 1;
178
179         // get interpolation weights
180         const float s1 = x - (float)x0;
181         const float s0 = 1.0f - s1;
182         const float t1 = y - (float)y0;
183         const float t0 = 1.0f - t1;
184         const float u1 = z - (float)z0;
185         const float u0 = 1.0f - u1;
186
187         const int slabSize = xres*yres;
188         const int i000 = x0 + y0 * xres + z0 * slabSize;
189         const int i010 = x0 + y1 * xres + z0 * slabSize;
190         const int i100 = x1 + y0 * xres + z0 * slabSize;
191         const int i110 = x1 + y1 * xres + z0 * slabSize;
192         const int i001 = x0 + y0 * xres + z1 * slabSize;
193         const int i011 = x0 + y1 * xres + z1 * slabSize;
194         const int i101 = x1 + y0 * xres + z1 * slabSize;
195         const int i111 = x1 + y1 * xres + z1 * slabSize;
196
197         // interpolate (indices could be computed once)
198         return Vec3(
199                         ( u0 * (s0 * (t0 * field1[i000] +
200                                 t1 * field1[i010]) +
201                                 s1 * (t0 * field1[i100] +
202                                 t1 * field1[i110])) +
203                                 u1 * (s0 * (t0 * field1[i001] +
204                                 t1 * field1[i011]) +
205                                 s1 * (t0 * field1[i101] +
206                                 t1 * field1[i111])) ) , 
207                         ( u0 * (s0 * (t0 * field2[i000] +
208                                 t1 * field2[i010]) +
209                                 s1 * (t0 * field2[i100] +
210                                 t1 * field2[i110])) +
211                                 u1 * (s0 * (t0 * field2[i001] +
212                                 t1 * field2[i011]) +
213                                 s1 * (t0 * field2[i101] +
214                                 t1 * field2[i111])) ) , 
215                         ( u0 * (s0 * (t0 * field3[i000] +
216                                 t1 * field3[i010]) +
217                                 s1 * (t0 * field3[i100] +
218                                 t1 * field3[i110])) +
219                                 u1 * (s0 * (t0 * field3[i001] +
220                                 t1 * field3[i011]) +
221                                 s1 * (t0 * field3[i101] +
222                                 t1 * field3[i111])) ) 
223                         );
224 }
225
226 };
227 #endif