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