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