BLI_math_rotation: properly name the quaternion power function.
[blender.git] / source / blender / blenlib / intern / voxel.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program 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 this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): Matt Ebb, Raul Fernandez Hernandez (Farsthary).
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27
28 /** \file blender/blenlib/intern/voxel.c
29  *  \ingroup bli
30  */
31
32 #include "BLI_utildefines.h"
33 #include "BLI_voxel.h"
34
35 #include "BLI_strict_flags.h"
36
37
38 BLI_INLINE float D(float *data, const int res[3], int x, int y, int z)
39 {
40         CLAMP(x, 0, res[0] - 1);
41         CLAMP(y, 0, res[1] - 1);
42         CLAMP(z, 0, res[2] - 1);
43         return data[BLI_VOXEL_INDEX(x, y, z, res)];
44 }
45
46 /* *** nearest neighbor *** */
47 /* input coordinates must be in bounding box 0.0 - 1.0 */
48 float BLI_voxel_sample_nearest(float *data, const int res[3], const float co[3])
49 {
50         int xi, yi, zi;
51
52         xi = (int)(co[0] * (float)res[0]);
53         yi = (int)(co[1] * (float)res[1]);
54         zi = (int)(co[2] * (float)res[2]);
55
56         return D(data, res, xi, yi, zi);
57 }
58
59 /* returns highest integer <= x as integer (slightly faster than floor()) */
60 BLI_INLINE int FLOORI(float x)
61 {
62         const int r = (int)x;
63         return ((x >= 0.f) || (float)r == x) ? r : (r - 1);
64 }
65
66 /* clamp function, cannot use the CLAMPIS macro, it sometimes returns unwanted results apparently related to
67  * gcc optimization flag -fstrict-overflow which is enabled at -O2
68  *
69  * this causes the test (x + 2) < 0 with int x == 2147483647 to return false (x being an integer,
70  * x + 2 should wrap around to -2147483647 so the test < 0 should return true, which it doesn't) */
71 BLI_INLINE int64_t _clamp(int a, int b, int c)
72 {
73         return (a < b) ? b : ((a > c) ? c : a);
74 }
75
76 float BLI_voxel_sample_trilinear(float *data, const int res[3], const float co[3])
77 {
78         if (data) {
79
80                 const float xf = co[0] * (float)res[0] - 0.5f;
81                 const float yf = co[1] * (float)res[1] - 0.5f;
82                 const float zf = co[2] * (float)res[2] - 0.5f;
83
84                 const int x = FLOORI(xf), y = FLOORI(yf), z = FLOORI(zf);
85
86                 const int64_t xc[2] = {
87                     _clamp(x,     0, res[0] - 1),
88                     _clamp(x + 1, 0, res[0] - 1),
89                 };
90                 const int64_t yc[2] = {
91                     _clamp(y,     0, res[1] - 1) * res[0],
92                     _clamp(y + 1, 0, res[1] - 1) * res[0],
93                 };
94                 const int64_t zc[2] = {
95                     _clamp(z,     0, res[2] - 1) * res[0] * res[1],
96                     _clamp(z + 1, 0, res[2] - 1) * res[0] * res[1],
97                 };
98
99                 const float dx = xf - (float)x;
100                 const float dy = yf - (float)y;
101                 const float dz = zf - (float)z;
102
103                 const float u[2] = {1.f - dx, dx};
104                 const float v[2] = {1.f - dy, dy};
105                 const float w[2] = {1.f - dz, dz};
106
107                 return w[0] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[0]] + u[1] * data[xc[1] + yc[0] + zc[0]] )
108                                 + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[0]] + u[1] * data[xc[1] + yc[1] + zc[0]] ) )
109                      + w[1] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[1]] + u[1] * data[xc[1] + yc[0] + zc[1]] )
110                                 + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[1]] + u[1] * data[xc[1] + yc[1] + zc[1]] ) );
111
112         }
113         return 0.f;
114 }
115
116 float BLI_voxel_sample_triquadratic(float *data, const int res[3], const float co[3])
117 {
118         if (data) {
119
120                 const float xf = co[0] * (float)res[0];
121                 const float yf = co[1] * (float)res[1];
122                 const float zf = co[2] * (float)res[2];
123                 const int x = FLOORI(xf), y = FLOORI(yf), z = FLOORI(zf);
124
125                 const int64_t xc[3] = {
126                     _clamp(x - 1, 0, res[0] - 1),
127                     _clamp(x,     0, res[0] - 1),
128                     _clamp(x + 1, 0, res[0] - 1),
129                 };
130                 const int64_t yc[3] = {
131                     _clamp(y - 1, 0, res[1] - 1) * res[0],
132                     _clamp(y,     0, res[1] - 1) * res[0],
133                     _clamp(y + 1, 0, res[1] - 1) * res[0],
134                 };
135                 const int64_t zc[3] = {
136                     _clamp(z - 1, 0, res[2] - 1) * res[0] * res[1],
137                     _clamp(z,     0, res[2] - 1) * res[0] * res[1],
138                     _clamp(z + 1, 0, res[2] - 1) * res[0] * res[1],
139                 };
140
141                 const float dx = xf - (float)x, dy = yf - (float)y, dz = zf - (float)z;
142                 const float u[3] = {dx * (0.5f * dx - 1.f) + 0.5f, dx * (1.0f - dx) + 0.5f, 0.5f * dx * dx};
143                 const float v[3] = {dy * (0.5f * dy - 1.f) + 0.5f, dy * (1.0f - dy) + 0.5f, 0.5f * dy * dy};
144                 const float w[3] = {dz * (0.5f * dz - 1.f) + 0.5f, dz * (1.0f - dz) + 0.5f, 0.5f * dz * dz};
145
146                 return w[0] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[0]] + u[1] * data[xc[1] + yc[0] + zc[0]] + u[2] * data[xc[2] + yc[0] + zc[0]] )
147                                 + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[0]] + u[1] * data[xc[1] + yc[1] + zc[0]] + u[2] * data[xc[2] + yc[1] + zc[0]] )
148                                 + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[0]] + u[1] * data[xc[1] + yc[2] + zc[0]] + u[2] * data[xc[2] + yc[2] + zc[0]] ) )
149                      + w[1] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[1]] + u[1] * data[xc[1] + yc[0] + zc[1]] + u[2] * data[xc[2] + yc[0] + zc[1]] )
150                                 + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[1]] + u[1] * data[xc[1] + yc[1] + zc[1]] + u[2] * data[xc[2] + yc[1] + zc[1]] )
151                                 + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[1]] + u[1] * data[xc[1] + yc[2] + zc[1]] + u[2] * data[xc[2] + yc[2] + zc[1]] ) )
152                      + w[2] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[2]] + u[1] * data[xc[1] + yc[0] + zc[2]] + u[2] * data[xc[2] + yc[0] + zc[2]] )
153                                 + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[2]] + u[1] * data[xc[1] + yc[1] + zc[2]] + u[2] * data[xc[2] + yc[1] + zc[2]] )
154                                 + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[2]] + u[1] * data[xc[1] + yc[2] + zc[2]] + u[2] * data[xc[2] + yc[2] + zc[2]] ) );
155
156         }
157         return 0.f;
158 }
159
160 float BLI_voxel_sample_tricubic(float *data, const int res[3], const float co[3], int bspline)
161 {
162         if (data) {
163
164                 const float xf = co[0] * (float)res[0] - 0.5f;
165                 const float yf = co[1] * (float)res[1] - 0.5f;
166                 const float zf = co[2] * (float)res[2] - 0.5f;
167                 const int x = FLOORI(xf), y = FLOORI(yf), z = FLOORI(zf);
168
169                 const int64_t xc[4] = {
170                     _clamp(x - 1, 0, res[0] - 1),
171                     _clamp(x,     0, res[0] - 1),
172                     _clamp(x + 1, 0, res[0] - 1),
173                     _clamp(x + 2, 0, res[0] - 1),
174                 };
175                 const int64_t yc[4] = {
176                     _clamp(y - 1, 0, res[1] - 1) * res[0],
177                     _clamp(y,     0, res[1] - 1) * res[0],
178                     _clamp(y + 1, 0, res[1] - 1) * res[0],
179                     _clamp(y + 2, 0, res[1] - 1) * res[0],
180                 };
181                 const int64_t zc[4] = {
182                     _clamp(z - 1, 0, res[2] - 1) * res[0] * res[1],
183                     _clamp(z,     0, res[2] - 1) * res[0] * res[1],
184                     _clamp(z + 1, 0, res[2] - 1) * res[0] * res[1],
185                     _clamp(z + 2, 0, res[2] - 1) * res[0] * res[1],
186                 };
187                 const float dx = xf - (float)x, dy = yf - (float)y, dz = zf - (float)z;
188
189                 float u[4], v[4], w[4];
190                 if (bspline) {  // B-Spline
191                         u[0] = (((-1.f/6.f)*dx + 0.5f)*dx - 0.5f)*dx + (1.f/6.f);
192                         u[1] =  ((     0.5f*dx - 1.f )*dx       )*dx + (2.f/3.f);
193                         u[2] =  ((    -0.5f*dx + 0.5f)*dx + 0.5f)*dx + (1.f/6.f);
194                         u[3] =   ( 1.f/6.f)*dx*dx*dx;
195                         v[0] = (((-1.f/6.f)*dy + 0.5f)*dy - 0.5f)*dy + (1.f/6.f);
196                         v[1] =  ((     0.5f*dy - 1.f )*dy       )*dy + (2.f/3.f);
197                         v[2] =  ((    -0.5f*dy + 0.5f)*dy + 0.5f)*dy + (1.f/6.f);
198                         v[3] =  ( 1.f/6.f)*dy*dy*dy;
199                         w[0] = (((-1.f/6.f)*dz + 0.5f)*dz - 0.5f)*dz + (1.f/6.f);
200                         w[1] =  ((     0.5f*dz - 1.f )*dz       )*dz + (2.f/3.f);
201                         w[2] =  ((    -0.5f*dz + 0.5f)*dz + 0.5f)*dz + (1.f/6.f);
202                         w[3] =   ( 1.f/6.f)*dz*dz*dz;
203                 }
204                 else {  // Catmull-Rom
205                         u[0] = ((-0.5f*dx + 1.0f)*dx - 0.5f)*dx;
206                         u[1] = (( 1.5f*dx - 2.5f)*dx       )*dx + 1.0f;
207                         u[2] = ((-1.5f*dx + 2.0f)*dx + 0.5f)*dx;
208                         u[3] = (( 0.5f*dx - 0.5f)*dx       )*dx;
209                         v[0] = ((-0.5f*dy + 1.0f)*dy - 0.5f)*dy;
210                         v[1] = (( 1.5f*dy - 2.5f)*dy       )*dy + 1.0f;
211                         v[2] = ((-1.5f*dy + 2.0f)*dy + 0.5f)*dy;
212                         v[3] = (( 0.5f*dy - 0.5f)*dy       )*dy;
213                         w[0] = ((-0.5f*dz + 1.0f)*dz - 0.5f)*dz;
214                         w[1] = (( 1.5f*dz - 2.5f)*dz       )*dz + 1.0f;
215                         w[2] = ((-1.5f*dz + 2.0f)*dz + 0.5f)*dz;
216                         w[3] = (( 0.5f*dz - 0.5f)*dz       )*dz;
217                 }
218
219                 return w[0] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[0]] + u[1] * data[xc[1] + yc[0] + zc[0]] + u[2] * data[xc[2] + yc[0] + zc[0]] + u[3] * data[xc[3] + yc[0] + zc[0]] )
220                                 + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[0]] + u[1] * data[xc[1] + yc[1] + zc[0]] + u[2] * data[xc[2] + yc[1] + zc[0]] + u[3] * data[xc[3] + yc[1] + zc[0]] )
221                                 + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[0]] + u[1] * data[xc[1] + yc[2] + zc[0]] + u[2] * data[xc[2] + yc[2] + zc[0]] + u[3] * data[xc[3] + yc[2] + zc[0]] )
222                                 + v[3] * ( u[0] * data[xc[0] + yc[3] + zc[0]] + u[1] * data[xc[1] + yc[3] + zc[0]] + u[2] * data[xc[2] + yc[3] + zc[0]] + u[3] * data[xc[3] + yc[3] + zc[0]] ) )
223                      + w[1] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[1]] + u[1] * data[xc[1] + yc[0] + zc[1]] + u[2] * data[xc[2] + yc[0] + zc[1]] + u[3] * data[xc[3] + yc[0] + zc[1]] )
224                                 + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[1]] + u[1] * data[xc[1] + yc[1] + zc[1]] + u[2] * data[xc[2] + yc[1] + zc[1]] + u[3] * data[xc[3] + yc[1] + zc[1]] )
225                                 + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[1]] + u[1] * data[xc[1] + yc[2] + zc[1]] + u[2] * data[xc[2] + yc[2] + zc[1]] + u[3] * data[xc[3] + yc[2] + zc[1]] )
226                                 + v[3] * ( u[0] * data[xc[0] + yc[3] + zc[1]] + u[1] * data[xc[1] + yc[3] + zc[1]] + u[2] * data[xc[2] + yc[3] + zc[1]] + u[3] * data[xc[3] + yc[3] + zc[1]] ) )
227                      + w[2] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[2]] + u[1] * data[xc[1] + yc[0] + zc[2]] + u[2] * data[xc[2] + yc[0] + zc[2]] + u[3] * data[xc[3] + yc[0] + zc[2]] )
228                                 + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[2]] + u[1] * data[xc[1] + yc[1] + zc[2]] + u[2] * data[xc[2] + yc[1] + zc[2]] + u[3] * data[xc[3] + yc[1] + zc[2]] )
229                                 + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[2]] + u[1] * data[xc[1] + yc[2] + zc[2]] + u[2] * data[xc[2] + yc[2] + zc[2]] + u[3] * data[xc[3] + yc[2] + zc[2]] )
230                                 + v[3] * ( u[0] * data[xc[0] + yc[3] + zc[2]] + u[1] * data[xc[1] + yc[3] + zc[2]] + u[2] * data[xc[2] + yc[3] + zc[2]] + u[3] * data[xc[3] + yc[3] + zc[2]] ) )
231                      + w[3] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[3]] + u[1] * data[xc[1] + yc[0] + zc[3]] + u[2] * data[xc[2] + yc[0] + zc[3]] + u[3] * data[xc[3] + yc[0] + zc[3]] )
232                                 + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[3]] + u[1] * data[xc[1] + yc[1] + zc[3]] + u[2] * data[xc[2] + yc[1] + zc[3]] + u[3] * data[xc[3] + yc[1] + zc[3]] )
233                                 + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[3]] + u[1] * data[xc[1] + yc[2] + zc[3]] + u[2] * data[xc[2] + yc[2] + zc[3]] + u[3] * data[xc[3] + yc[2] + zc[3]] )
234                                 + v[3] * ( u[0] * data[xc[0] + yc[3] + zc[3]] + u[1] * data[xc[1] + yc[3] + zc[3]] + u[2] * data[xc[2] + yc[3] + zc[3]] + u[3] * data[xc[3] + yc[3] + zc[3]] ) );
235
236         }
237         return 0.f;
238 }