svn merge -r 23207:23528 https://svn.blender.org/svnroot/bf-blender/trunk/blender
[blender.git] / source / blender / blenlib / intern / voxel.c
1 /**
2  *
3  * ***** BEGIN GPL LICENSE BLOCK *****
4  *
5  * This program is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU General Public License
7  * as published by the Free Software Foundation; either version 2
8  * of the License, or (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program; if not, write to the Free Software Foundation,
17  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
18  *
19  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
20  * All rights reserved.
21  *
22  * The Original Code is: all of this file.
23  *
24  * Contributor(s): Matt Ebb, Raul Fernandez Hernandez (Farsthary).
25  *
26  * ***** END GPL LICENSE BLOCK *****
27  */
28 #include <math.h>
29
30 #include "BLI_voxel.h"
31
32 #include "BKE_utildefines.h"
33
34
35 #if defined( _MSC_VER ) && !defined( __cplusplus )
36 # define inline __inline
37 #endif // defined( _MSC_VER ) && !defined( __cplusplus )
38
39 static inline float D(float *data,  int *res, int x, int y, int z)
40 {
41         CLAMP(x, 0, res[0]-1);
42         CLAMP(y, 0, res[1]-1);
43         CLAMP(z, 0, res[2]-1);
44         return data[ V_I(x, y, z, res) ];
45 }
46
47 /* *** nearest neighbour *** */
48 /* input coordinates must be in bounding box 0.0 - 1.0 */
49 float voxel_sample_nearest(float *data, int *res, float *co)
50 {
51         int xi, yi, zi;
52         
53         xi = co[0] * res[0];
54         yi = co[1] * res[1];
55         zi = co[2] * res[2];
56         
57         return D(data, res, xi, yi, zi);
58 }
59
60 // returns highest integer <= x as integer (slightly faster than floor())
61 inline int FLOORI(float x)
62 {
63         const int r = (int)x;
64         return ((x >= 0.f) || (float)r == x) ? r : (r - 1);
65 }
66
67 // clamp function, cannot use the CLAMPIS macro, it sometimes returns unwanted results apparently related to gcc optimization flag -fstrict-overflow which is enabled at -O2
68 // this causes the test (x + 2) < 0 with int x == 2147483647 to return false (x being an integer, x + 2 should wrap around to -2147483647 so the test < 0 should return true, which it doesn't)
69 inline int _clamp(int a, int b, int c)
70 {
71         return (a < b) ? b : ((a > c) ? c : a);
72 }
73
74 float voxel_sample_trilinear(float *data, int *res, float *co)
75 {
76         if (data) {
77         
78                 const float xf = co[0] * res[0] - 0.5f;
79                 const float yf = co[1] * res[1] - 0.5f;
80                 const float zf = co[2] * res[2] - 0.5f;
81                 
82                 const int x = FLOORI(xf), y = FLOORI(yf), z = FLOORI(zf);
83         
84                 const int xc[2] = {_clamp(x, 0, res[0] - 1), _clamp(x + 1, 0, res[0] - 1)};
85                 const int yc[2] = {res[0] * _clamp(y, 0, res[1] - 1), res[0] * _clamp(y + 1, 0, res[1] - 1)};
86                 const int zc[2] = {res[0] * res[1] * _clamp(z, 0, res[2] - 1), res[0] * res[1] * _clamp(z + 1, 0, res[2] - 1)};
87         
88                 const float dx = xf - (float)x;
89                 const float dy = yf - (float)y;
90                 const float dz = zf - (float)z;
91                 
92                 const float u[2] = {1.f - dx, dx};
93                 const float v[2] = {1.f - dy, dy};
94                 const float w[2] = {1.f - dz, dz};
95         
96                 return w[0] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[0]] + u[1] * data[xc[1] + yc[0] + zc[0]] )
97                                 + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[0]] + u[1] * data[xc[1] + yc[1] + zc[0]] ) )
98                      + w[1] * (   v[0] * ( u[0] * data[xc[0] + yc[0] + zc[1]] + u[1] * data[xc[1] + yc[0] + zc[1]] )
99                                 + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[1]] + u[1] * data[xc[1] + yc[1] + zc[1]] ) );
100         
101         }
102         return 0.f;
103 }
104         
105
106 float voxel_sample_triquadratic(float *data, int *res, float *co)
107 {
108         if (data) {
109
110                 const float xf = co[0] * res[0], yf = co[1] * res[1], zf = co[2] * res[2];
111                 const int x = FLOORI(xf), y = FLOORI(yf), z = FLOORI(zf);
112
113                 const int xc[3] = {_clamp(x - 1, 0, res[0] - 1), _clamp(x, 0, res[0] - 1), _clamp(x + 1, 0, res[0] - 1)};
114                 const int yc[3] = {res[0] * _clamp(y - 1, 0, res[1] - 1), res[0] * _clamp(y, 0, res[1] - 1), res[0] * _clamp(y + 1, 0, res[1] - 1)};
115                 const int zc[3] = {res[0] * res[1] * _clamp(z - 1, 0, res[2] - 1), res[0] * res[1] * _clamp(z, 0, res[2] - 1), res[0] * res[1] * _clamp(z + 1, 0, res[2] - 1)};
116
117                 const float dx = xf - (float)x, dy = yf - (float)y, dz = zf - (float)z;
118                 const float u[3] = {dx*(0.5f*dx - 1.f) + 0.5f, dx*(1.f - dx) + 0.5f, 0.5f*dx*dx};
119                 const float v[3] = {dy*(0.5f*dy - 1.f) + 0.5f, dy*(1.f - dy) + 0.5f, 0.5f*dy*dy};
120                 const float w[3] = {dz*(0.5f*dz - 1.f) + 0.5f, dz*(1.f - dz) + 0.5f, 0.5f*dz*dz};
121
122                 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]] )
123                                 + 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]] )
124                                 + 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]] ) )
125                      + 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]] )
126                                 + 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]] )
127                                 + 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]] ) )
128                      + 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]] )
129                                 + 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]] )
130                                 + 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]] ) );
131
132 }
133         return 0.f;
134 }
135
136 float voxel_sample_tricubic(float *data, int *res, float *co, int bspline)
137 {
138         if (data) {
139
140                 const float xf = co[0] * res[0] - 0.5f, yf = co[1] * res[1] - 0.5f, zf = co[2] * res[2] - 0.5f;
141                 const int x = FLOORI(xf), y = FLOORI(yf), z = FLOORI(zf);
142
143                 const int xc[4] = {_clamp(x - 1, 0, res[0] - 1), _clamp(x, 0, res[0] - 1), _clamp(x + 1, 0, res[0] - 1), _clamp(x + 2, 0, res[0] - 1)};
144                 const int yc[4] = {res[0] * _clamp(y - 1, 0, res[1] - 1), res[0] * _clamp(y, 0, res[1] - 1), res[0] * _clamp(y + 1, 0, res[1] - 1), res[0] * _clamp(y + 2, 0, res[1] - 1)};
145                 const int zc[4] = {res[0] * res[1] * _clamp(z - 1, 0, res[2] - 1), res[0] * res[1] * _clamp(z, 0, res[2] - 1), res[0] * res[1] * _clamp(z + 1, 0, res[2] - 1), res[0] * res[1] * _clamp(z + 2, 0, res[2] - 1)};
146
147                 const float dx = xf - (float)x, dy = yf - (float)y, dz = zf - (float)z;
148
149                 float u[4], v[4], w[4];
150                 if (bspline) {  // B-Spline
151                         u[0] = (((-1.f/6.f)*dx + 0.5f)*dx - 0.5f)*dx + (1.f/6.f);
152                         u[1] =  ((     0.5f*dx - 1.f )*dx       )*dx + (2.f/3.f);
153                         u[2] =  ((    -0.5f*dx + 0.5f)*dx + 0.5f)*dx + (1.f/6.f);
154                         u[3] =   ( 1.f/6.f)*dx*dx*dx;
155                         v[0] = (((-1.f/6.f)*dy + 0.5f)*dy - 0.5f)*dy + (1.f/6.f);
156                         v[1] =  ((     0.5f*dy - 1.f )*dy       )*dy + (2.f/3.f);
157                         v[2] =  ((    -0.5f*dy + 0.5f)*dy + 0.5f)*dy + (1.f/6.f);
158                         v[3] =  ( 1.f/6.f)*dy*dy*dy;
159                         w[0] = (((-1.f/6.f)*dz + 0.5f)*dz - 0.5f)*dz + (1.f/6.f);
160                         w[1] =  ((     0.5f*dz - 1.f )*dz       )*dz + (2.f/3.f);
161                         w[2] =  ((    -0.5f*dz + 0.5f)*dz + 0.5f)*dz + (1.f/6.f);
162                         w[3] =   ( 1.f/6.f)*dz*dz*dz;
163                 }
164                 else {  // Catmull-Rom
165                         u[0] = ((-0.5f*dx + 1.0f)*dx - 0.5f)*dx;
166                         u[1] = (( 1.5f*dx - 2.5f)*dx       )*dx + 1.0f;
167                         u[2] = ((-1.5f*dx + 2.0f)*dx + 0.5f)*dx;
168                         u[3] = (( 0.5f*dx - 0.5f)*dx       )*dx;
169                         v[0] = ((-0.5f*dy + 1.0f)*dy - 0.5f)*dy;
170                         v[1] = (( 1.5f*dy - 2.5f)*dy       )*dy + 1.0f;
171                         v[2] = ((-1.5f*dy + 2.0f)*dy + 0.5f)*dy;
172                         v[3] = (( 0.5f*dy - 0.5f)*dy       )*dy;
173                         w[0] = ((-0.5f*dz + 1.0f)*dz - 0.5f)*dz;
174                         w[1] = (( 1.5f*dz - 2.5f)*dz       )*dz + 1.0f;
175                         w[2] = ((-1.5f*dz + 2.0f)*dz + 0.5f)*dz;
176                         w[3] = (( 0.5f*dz - 0.5f)*dz       )*dz;
177                 }
178
179                 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]] )
180                                 + 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]] )
181                                 + 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]] )
182                                 + 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]] ) )
183                      + 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]] )
184                                 + 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]] )
185                                 + 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]] )
186                                 + 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]] ) )
187                      + 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]] )
188                                 + 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]] )
189                                 + 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]] )
190                                 + 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]] ) )
191                      + 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]] )
192                                 + 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]] )
193                                 + 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]] )
194                                 + 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]] ) );
195
196         }
197         return 0.f;
198 }