1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkImageBSplineInternals.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15
16 // This code has been modified from the original C code from Thevenaz.
17 // The functions have been converted into static C++ class methods,
18 // and two new boundary conditions have been added: clamped and repeat.
19
20 /*****************************************************************************
21 * Date: January 5, 2009
22 *----------------------------------------------------------------------------
23 * This C program is based on the following paper:
24 * P. Thevenaz, T. Blu, M. Unser, "Interpolation Revisited,"
25 * IEEE Transactions on Medical Imaging,
26 * vol. 19, no. 7, pp. 739-758, July 2000.
27 *----------------------------------------------------------------------------
28 * Philippe Thevenaz
29 * EPFL/STI/IMT/LIB/BM.4.137
30 * Station 17
31 * CH-1015 Lausanne VD
32 *----------------------------------------------------------------------------
33 * phone (CET): +41(21)693.51.61
34 * fax: +41(21)693.37.01
35 * RFC-822: philippe.thevenaz@epfl.ch
36 * X-400: /C=ch/A=400net/P=switch/O=epfl/S=thevenaz/G=philippe/
37 * URL: http://bigwww.epfl.ch/
38 ****************************************************************************/
39
40 #include "vtkImageBSplineInternals.h"
41 #include "vtkAbstractImageInterpolator.h"
42 #include <cstddef>
43 #include <cmath>
44
45 /*--------------------------------------------------------------------------*/
46 void vtkImageBSplineInternals::
ConvertToInterpolationCoefficients(double c[],long DataLength,long Border,double z[],long NbPoles,double Tolerance)47 ConvertToInterpolationCoefficients
48 (
49 double c[], /* input samples --> output coefficients */
50 long DataLength, /* number of samples or coefficients */
51 long Border, /* border mode */
52 double z[], /* poles */
53 long NbPoles, /* number of poles */
54 double Tolerance /* admissible relative error */
55 )
56
57 { /* begin ConvertToInterpolationCoefficients */
58
59 double Lambda = 1.0;
60 long n, k;
61
62 /* special case required by mirror boundaries */
63 if (DataLength == 1L) {
64 return;
65 }
66 /* compute the overall gain */
67 for (k = 0L; k < NbPoles; k++) {
68 Lambda = Lambda * (1.0 - z[k]) * (1.0 - 1.0 / z[k]);
69 }
70 /* apply the gain */
71 for (n = 0L; n < DataLength; n++) {
72 c[n] *= Lambda;
73 }
74 /* loop over all poles */
75 for (k = 0L; k < NbPoles; k++) {
76 /* causal initialization */
77 c[0] = InitialCausalCoefficient(c, DataLength, Border, z[k], Tolerance);
78 /* causal recursion */
79 for (n = 1L; n < DataLength; n++) {
80 c[n] += z[k] * c[n - 1L];
81 }
82 /* anticausal initialization */
83 c[DataLength - 1L] = InitialAntiCausalCoefficient(c, DataLength, Border, z[k], Tolerance);
84 /* anticausal recursion */
85 for (n = DataLength - 2L; 0 <= n; n--) {
86 c[n] = z[k] * (c[n + 1L] - c[n]);
87 }
88 }
89 } /* end ConvertToInterpolationCoefficients */
90
91 /*--------------------------------------------------------------------------*/
92 double vtkImageBSplineInternals::
InitialCausalCoefficient(double c[],long DataLength,long Border,double z,double Tolerance)93 InitialCausalCoefficient
94 (
95 double c[], /* coefficients */
96 long DataLength, /* number of coefficients */
97 long Border, /* border mode */
98 double z, /* actual pole */
99 double Tolerance /* admissible relative error */
100 )
101
102 { /* begin InitialCausalCoefficient */
103
104 double Sum, zn, z2n, iz;
105 long n, Horizon;
106
107 switch (Border) {
108 case VTK_IMAGE_BORDER_CLAMP:
109 /* this initialization corresponds to repeating edge pixels */
110 Horizon = DataLength;
111 if (Tolerance > 0.0 && DataLength > 16) {
112 Horizon = (long)(ceil(log(Tolerance) / log(fabs(z))));
113 }
114 if (Horizon < DataLength) {
115 /* accelerated loop */
116 zn = z;
117 Sum = c[0];
118 for (n = 0; n < Horizon; n++) {
119 Sum += zn * c[n];
120 zn *= z;
121 }
122 return(Sum);
123 }
124 else {
125 /* full loop */
126 zn = z;
127 iz = 1.0 / z;
128 z2n = pow(z, (double)DataLength);
129 Sum = z * c[0] + z2n * z2n * c[0];
130 z2n *= z2n * iz;
131 for (n = 1L; n <= DataLength - 1L; n++) {
132 zn *= z;
133 Sum += (zn + z2n) * c[n];
134 z2n *= iz;
135 }
136 return(c[0] + Sum / (1.0 - zn * zn));
137 }
138 /* break; (unnecessary because of return) */
139
140 case VTK_IMAGE_BORDER_MIRROR:
141 /* this initialization corresponds to mirror boundaries */
142 Horizon = DataLength;
143 if (Tolerance > 0.0 && DataLength > 16) {
144 Horizon = (long)(ceil(log(Tolerance) / log(fabs(z))));
145 }
146 if (Horizon < DataLength) {
147 /* accelerated loop */
148 zn = z;
149 Sum = c[0];
150 for (n = 1L; n < Horizon; n++) {
151 Sum += zn * c[n];
152 zn *= z;
153 }
154 return(Sum);
155 }
156 else {
157 /* full loop */
158 zn = z;
159 iz = 1.0 / z;
160 z2n = pow(z, (double)(DataLength - 1L));
161 Sum = c[0] + z2n * c[DataLength - 1L];
162 z2n *= z2n * iz;
163 for (n = 1L; n <= DataLength - 2L; n++) {
164 Sum += (zn + z2n) * c[n];
165 zn *= z;
166 z2n *= iz;
167 }
168 return(Sum / (1.0 - zn * zn));
169 }
170 /* break; (unnecessary because of return) */
171
172 case VTK_IMAGE_BORDER_REPEAT:
173 /* this initialization corresponds to periodic boundaries */
174 Horizon = DataLength;
175 if (Tolerance > 0.0 && DataLength > 16) {
176 Horizon = (long)(ceil(log(Tolerance) / log(fabs(z))));
177 }
178 if (Horizon < DataLength) {
179 /* accelerated loop */
180 zn = z;
181 Sum = c[0];
182 for (n = 1L; n < Horizon; n++) {
183 Sum += zn * c[DataLength-n];
184 zn *= z;
185 }
186 return(Sum);
187 }
188 else {
189 /* full loop */
190 zn = z;
191 Sum = c[0];
192 for (n = 1L; n < DataLength; n++) {
193 Sum += zn * c[DataLength-n];
194 zn *= z;
195 }
196 return(Sum / (1.0 - zn));
197 }
198 /* break; (unnecessary because of return) */
199 }
200
201 return(0.0);
202 } /* end InitialCausalCoefficient */
203
204 /*--------------------------------------------------------------------------*/
205 double vtkImageBSplineInternals::
InitialAntiCausalCoefficient(double c[],long DataLength,long Border,double z,double Tolerance)206 InitialAntiCausalCoefficient
207 (
208 double c[], /* coefficients */
209 long DataLength, /* number of samples or coefficients */
210 long Border, /* border mode */
211 double z, /* actual pole */
212 double Tolerance /* admissible relative error */
213 )
214
215 { /* begin InitialAntiCausalCoefficient */
216 double Sum;
217 double zn;
218 long n, Horizon;
219
220 switch (Border) {
221 case VTK_IMAGE_BORDER_CLAMP:
222 /* this initialization corresponds to repeating edge pixels */
223 return((z / (z - 1.0)) * c[DataLength - 1L]);
224 case VTK_IMAGE_BORDER_MIRROR:
225 /* this initialization corresponds to mirror boundaries */
226 return((z / (z * z - 1.0)) * (z * c[DataLength - 2L] + c[DataLength - 1L]));
227 case VTK_IMAGE_BORDER_REPEAT:
228 /* this initialization corresponds to periodic boundaries */
229 Horizon = DataLength;
230 if (Tolerance > 0.0 && DataLength > 16) {
231 Horizon = (long)(ceil(log(Tolerance) / log(fabs(z))));
232 }
233 if (Horizon < DataLength) {
234 /* accelerated loop */
235 zn = z;
236 Sum = c[0];
237 for (n = 1L; n < Horizon; n++) {
238 Sum += zn * c[n];
239 zn *= z;
240 }
241 return(-Sum * z * z - z * c[DataLength - 1L]);
242 }
243 else {
244 /* full loop */
245 zn = z;
246 Sum = c[0];
247 for (n = 1L; n < DataLength; n++) {
248 Sum += zn * c[n];
249 zn *= z;
250 }
251 return(Sum * z * z / (zn - 1.0) - z * c[DataLength - 1L]);
252 }
253 }
254 return(0.0);
255 } /* end InitialAntiCausalCoefficient */
256
257 /*--------------------------------------------------------------------------*/
258 int vtkImageBSplineInternals::
GetPoleValues(double Pole[4],long & NbPoles,long SplineDegree)259 GetPoleValues
260 (
261 double Pole[4], /* poles for b-spline */
262 long &NbPoles, /* number of poles (will be four, at maximum) */
263 long SplineDegree /* degree of the spline model */
264 )
265
266 { /* begin GetPoleValues */
267
268 /* recover the poles from a lookup table */
269 switch (SplineDegree) {
270 case 0L:
271 NbPoles = 0L;
272 break;
273 case 1L:
274 NbPoles = 0L;
275 break;
276 case 2L:
277 NbPoles = 1L;
278 Pole[0] = sqrt(8.0) - 3.0;
279 break;
280 case 3L:
281 NbPoles = 1L;
282 Pole[0] = sqrt(3.0) - 2.0;
283 break;
284 case 4L:
285 NbPoles = 2L;
286 Pole[0] = sqrt(664.0 - sqrt(438976.0)) + sqrt(304.0) - 19.0;
287 Pole[1] = sqrt(664.0 + sqrt(438976.0)) - sqrt(304.0) - 19.0;
288 break;
289 case 5L:
290 NbPoles = 2L;
291 Pole[0] = sqrt(135.0 / 2.0 - sqrt(17745.0 / 4.0)) + sqrt(105.0 / 4.0)
292 - 13.0 / 2.0;
293 Pole[1] = sqrt(135.0 / 2.0 + sqrt(17745.0 / 4.0)) - sqrt(105.0 / 4.0)
294 - 13.0 / 2.0;
295 break;
296 case 6L:
297 NbPoles = 3L;
298 Pole[0] = -0.48829458930304475513011803888378906211227916123938;
299 Pole[1] = -0.081679271076237512597937765737059080653379610398148;
300 Pole[2] = -0.0014141518083258177510872439765585925278641690553467;
301 break;
302 case 7L:
303 NbPoles = 3L;
304 Pole[0] = -0.53528043079643816554240378168164607183392315234269;
305 Pole[1] = -0.12255461519232669051527226435935734360548654942730;
306 Pole[2] = -0.0091486948096082769285930216516478534156925639545994;
307 break;
308 case 8L:
309 NbPoles = 4L;
310 Pole[0] = -0.57468690924876543053013930412874542429066157804125;
311 Pole[1] = -0.16303526929728093524055189686073705223476814550830;
312 Pole[2] = -0.023632294694844850023403919296361320612665920854629;
313 Pole[3] = -0.00015382131064169091173935253018402160762964054070043;
314 break;
315 case 9L:
316 NbPoles = 4L;
317 Pole[0] = -0.60799738916862577900772082395428976943963471853991;
318 Pole[1] = -0.20175052019315323879606468505597043468089886575747;
319 Pole[2] = -0.043222608540481752133321142979429688265852380231497;
320 Pole[3] = -0.0021213069031808184203048965578486234220548560988624;
321 break;
322 default:
323 NbPoles = 0L;
324 return(1);
325 }
326
327 return(0);
328 } /* end GetPoleValues */
329
330 /*--------------------------------------------------------------------------*/
331 namespace {
332
333 template<class T>
vtkImageBSplineGetInterpolationWeights(T xWeight[10],double w,long SplineDegree)334 int vtkImageBSplineGetInterpolationWeights
335 (
336 T xWeight[10], /* weights, size is SplineDegree + 1 */
337 double w, /* offset, value between 0 and 1 */
338 long SplineDegree /* degree of spline, between 2 and 9 */
339 )
340 { /* begin GetInterpolationWeights */
341 double w2, w4, t, t0, t1;
342
343 /* compute the interpolation weights */
344 switch (SplineDegree) {
345 case 0L:
346 xWeight[0] = 1.0;
347 break;
348 case 1L:
349 xWeight[0] = 1.0 - w;
350 xWeight[1] = w;
351 break;
352 case 2L:
353 xWeight[1] = 3.0 / 4.0 - w * w;
354 xWeight[2] = (1.0 / 2.0) * (w - xWeight[1] + 1.0);
355 xWeight[0] = 1.0 - xWeight[1] - xWeight[2];
356 break;
357 case 3L:
358 xWeight[3] = (1.0 / 6.0) * w * w * w;
359 xWeight[0] = (1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - xWeight[3];
360 xWeight[2] = w + xWeight[0] - 2.0 * xWeight[3];
361 xWeight[1] = 1.0 - xWeight[0] - xWeight[2] - xWeight[3];
362 break;
363 case 4L:
364 w2 = w * w;
365 t = (1.0 / 6.0) * w2;
366 xWeight[0] = 1.0 / 2.0 - w;
367 xWeight[0] *= xWeight[0];
368 xWeight[0] *= (1.0 / 24.0) * xWeight[0];
369 t0 = w * (t - 11.0 / 24.0);
370 t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
371 xWeight[1] = t1 + t0;
372 xWeight[3] = t1 - t0;
373 xWeight[4] = xWeight[0] + t0 + (1.0 / 2.0) * w;
374 xWeight[2] = 1.0 - xWeight[0] - xWeight[1] - xWeight[3] - xWeight[4];
375 break;
376 case 5L:
377 w2 = w * w;
378 xWeight[5] = (1.0 / 120.0) * w * w2 * w2;
379 w2 -= w;
380 w4 = w2 * w2;
381 w -= 1.0 / 2.0;
382 t = w2 * (w2 - 3.0);
383 xWeight[0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - xWeight[5];
384 t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
385 t1 = (-1.0 / 12.0) * w * (t + 4.0);
386 xWeight[2] = t0 + t1;
387 xWeight[3] = t0 - t1;
388 t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
389 t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
390 xWeight[1] = t0 + t1;
391 xWeight[4] = t0 - t1;
392 break;
393 case 6L:
394 xWeight[0] = 1.0 / 2.0 - w;
395 xWeight[0] *= xWeight[0] * xWeight[0];
396 xWeight[0] *= xWeight[0] / 720.0;
397 xWeight[1] = (361.0 / 192.0 - w * (59.0 / 8.0 + w
398 * (-185.0 / 16.0 + w * (25.0 / 3.0 + w * (-5.0 / 2.0 + w)
399 * (1.0 / 2.0 + w))))) / 120.0;
400 xWeight[2] = (10543.0 / 960.0 + w * (-289.0 / 16.0 + w
401 * (79.0 / 16.0 + w * (43.0 / 6.0 + w * (-17.0 / 4.0 + w
402 * (-1.0 + w)))))) / 48.0;
403 w2 = w * w;
404 xWeight[3] = (5887.0 / 320.0 - w2 * (231.0 / 16.0 - w2
405 * (21.0 / 4.0 - w2))) / 36.0;
406 xWeight[4] = (10543.0 / 960.0 + w * (289.0 / 16.0 + w
407 * (79.0 / 16.0 + w * (-43.0 / 6.0 + w * (-17.0 / 4.0 + w
408 * (1.0 + w)))))) / 48.0;
409 xWeight[6] = 1.0 / 2.0 + w;
410 xWeight[6] *= xWeight[6] * xWeight[6];
411 xWeight[6] *= xWeight[6] / 720.0;
412 xWeight[5] = 1.0 - xWeight[0] - xWeight[1] - xWeight[2] - xWeight[3]
413 - xWeight[4] - xWeight[6];
414 break;
415 case 7L:
416 xWeight[0] = 1.0 - w;
417 xWeight[0] *= xWeight[0];
418 xWeight[0] *= xWeight[0] * xWeight[0];
419 xWeight[0] *= (1.0 - w) / 5040.0;
420 w2 = w * w;
421 xWeight[1] = (120.0 / 7.0 + w * (-56.0 + w * (72.0 + w
422 * (-40.0 + w2 * (12.0 + w * (-6.0 + w)))))) / 720.0;
423 xWeight[2] = (397.0 / 7.0 - w * (245.0 / 3.0 + w * (-15.0 + w
424 * (-95.0 / 3.0 + w * (15.0 + w * (5.0 + w
425 * (-5.0 + w))))))) / 240.0;
426 xWeight[3] = (2416.0 / 35.0 + w2 * (-48.0 + w2 * (16.0 + w2
427 * (-4.0 + w)))) / 144.0;
428 xWeight[4] = (1191.0 / 35.0 - w * (-49.0 + w * (-9.0 + w
429 * (19.0 + w * (-3.0 + w) * (-3.0 + w2))))) / 144.0;
430 xWeight[5] = (40.0 / 7.0 + w * (56.0 / 3.0 + w * (24.0 + w
431 * (40.0 / 3.0 + w2 * (-4.0 + w * (-2.0 + w)))))) / 240.0;
432 xWeight[7] = w2;
433 xWeight[7] *= xWeight[7] * xWeight[7];
434 xWeight[7] *= w / 5040.0;
435 xWeight[6] = 1.0 - xWeight[0] - xWeight[1] - xWeight[2] - xWeight[3]
436 - xWeight[4] - xWeight[5] - xWeight[7];
437 break;
438 case 8L:
439 xWeight[0] = 1.0 / 2.0 - w;
440 xWeight[0] *= xWeight[0];
441 xWeight[0] *= xWeight[0];
442 xWeight[0] *= xWeight[0] / 40320.0;
443 w2 = w * w;
444 xWeight[1] = (39.0 / 16.0 - w * (6.0 + w * (-9.0 / 2.0 + w2)))
445 * (21.0 / 16.0 + w * (-15.0 / 4.0 + w * (9.0 / 2.0 + w
446 * (-3.0 + w)))) / 5040.0;
447 xWeight[2] = (82903.0 / 1792.0 + w * (-4177.0 / 32.0 + w
448 * (2275.0 / 16.0 + w * (-487.0 / 8.0 + w * (-85.0 / 8.0 + w
449 * (41.0 / 2.0 + w * (-5.0 + w * (-2.0 + w)))))))) / 1440.0;
450 xWeight[3] = (310661.0 / 1792.0 - w * (14219.0 / 64.0 + w
451 * (-199.0 / 8.0 + w * (-1327.0 / 16.0 + w * (245.0 / 8.0 + w
452 * (53.0 / 4.0 + w * (-8.0 + w * (-1.0 + w)))))))) / 720.0;
453 xWeight[4] = (2337507.0 / 8960.0 + w2 * (-2601.0 / 16.0 + w2
454 * (387.0 / 8.0 + w2 * (-9.0 + w2)))) / 576.0;
455 xWeight[5] = (310661.0 / 1792.0 - w * (-14219.0 / 64.0 + w
456 * (-199.0 / 8.0 + w * (1327.0 / 16.0 + w * (245.0 / 8.0 + w
457 * (-53.0 / 4.0 + w * (-8.0 + w * (1.0 + w)))))))) / 720.0;
458 xWeight[7] = (39.0 / 16.0 - w * (-6.0 + w * (-9.0 / 2.0 + w2)))
459 * (21.0 / 16.0 + w * (15.0 / 4.0 + w * (9.0 / 2.0 + w
460 * (3.0 + w)))) / 5040.0;
461 xWeight[8] = 1.0 / 2.0 + w;
462 xWeight[8] *= xWeight[8];
463 xWeight[8] *= xWeight[8];
464 xWeight[8] *= xWeight[8] / 40320.0;
465 xWeight[6] = 1.0 - xWeight[0] - xWeight[1] - xWeight[2] - xWeight[3]
466 - xWeight[4] - xWeight[5] - xWeight[7] - xWeight[8];
467 break;
468 case 9L:
469 xWeight[0] = 1.0 - w;
470 xWeight[0] *= xWeight[0];
471 xWeight[0] *= xWeight[0];
472 xWeight[0] *= xWeight[0] * (1.0 - w) / 362880.0;
473 xWeight[1] = (502.0 / 9.0 + w * (-246.0 + w * (472.0 + w
474 * (-504.0 + w * (308.0 + w * (-84.0 + w * (-56.0 / 3.0 + w
475 * (24.0 + w * (-8.0 + w))))))))) / 40320.0;
476 xWeight[2] = (3652.0 / 9.0 - w * (2023.0 / 2.0 + w * (-952.0 + w
477 * (938.0 / 3.0 + w * (112.0 + w * (-119.0 + w * (56.0 / 3.0 + w
478 * (14.0 + w * (-7.0 + w))))))))) / 10080.0;
479 xWeight[3] = (44117.0 / 42.0 + w * (-2427.0 / 2.0 + w * (66.0 + w
480 * (434.0 + w * (-129.0 + w * (-69.0 + w * (34.0 + w * (6.0 + w
481 * (-6.0 + w))))))))) / 4320.0;
482 w2 = w * w;
483 xWeight[4] = (78095.0 / 63.0 - w2 * (700.0 + w2 * (-190.0 + w2
484 * (100.0 / 3.0 + w2 * (-5.0 + w))))) / 2880.0;
485 xWeight[5] = (44117.0 / 63.0 + w * (809.0 + w * (44.0 + w
486 * (-868.0 / 3.0 + w * (-86.0 + w * (46.0 + w * (68.0 / 3.0 + w
487 * (-4.0 + w * (-4.0 + w))))))))) / 2880.0;
488 xWeight[6] = (3652.0 / 21.0 - w * (-867.0 / 2.0 + w * (-408.0 + w
489 * (-134.0 + w * (48.0 + w * (51.0 + w * (-4.0 + w) * (-1.0 + w)
490 * (2.0 + w))))))) / 4320.0;
491 xWeight[7] = (251.0 / 18.0 + w * (123.0 / 2.0 + w * (118.0 + w
492 * (126.0 + w * (77.0 + w * (21.0 + w * (-14.0 / 3.0 + w
493 * (-6.0 + w * (-2.0 + w))))))))) / 10080.0;
494 xWeight[9] = w2 * w2;
495 xWeight[9] *= xWeight[9] * w / 362880.0;
496 xWeight[8] = 1.0 - xWeight[0] - xWeight[1] - xWeight[2] - xWeight[3]
497 - xWeight[4] - xWeight[5] - xWeight[6] - xWeight[7] - xWeight[9];
498 break;
499 default:
500 return(0);
501 }
502 return(1);
503 } /* end GetInterpolationWeights */
504
505 /*--------------------------------------------------------------------------*/
506 template<class T>
vtkImageBSplineInterpolatedValue(const T * Bcoeff,T * v,long Width,long Height,long Slices,long Depth,double x,double y,double z,long SplineDegree,long Border)507 int vtkImageBSplineInterpolatedValue
508 (
509 const T *Bcoeff, /* input B-spline array of coefficients */
510 T *v, /* resulting pixel value */
511 long Width, /* width of the image */
512 long Height, /* height of the image */
513 long Slices, /* number of slices of the image */
514 long Depth, /* number of samples per pixel */
515 double x, /* x coordinate where to interpolate */
516 double y, /* y coordinate where to interpolate */
517 double z, /* y coordinate where to interpolate */
518 long SplineDegree, /* degree of the spline model */
519 long Border /* what to do at the border */
520 )
521
522 { /* begin InterpolatedValue */
523
524 const T *p1, *p2, *p3;
525 T xWeight[10], yWeight[10], zWeight[10];
526 double interpolated;
527 double w, u;
528 double s, t, r;
529 long xIndex[10], yIndex[10], zIndex[10];
530 ptrdiff_t xIncrement, yIncrement, zIncrement;
531 long Width2 = 2L * Width - 2L;
532 long Height2 = 2L * Height - 2L;
533 long Slices2 = 2L * Slices - 2L;
534 long CentralIndex = SplineDegree / 2L;
535 long i, j, k, l, c;
536 long imax, jmax, kmax;
537
538 if (SplineDegree < 0 || SplineDegree > 9)
539 {
540 return 0;
541 }
542
543 /* check for 1D and 2D images */
544 imax = SplineDegree;
545 jmax = SplineDegree;
546 kmax = SplineDegree;
547 if (Width == 1)
548 {
549 imax = 0;
550 }
551 if (Height == 1)
552 {
553 jmax = 0;
554 }
555 if (Slices == 1)
556 {
557 kmax = 0;
558 }
559
560 /* compute the interpolation indexes and fractional offsets*/
561 if (SplineDegree & 1L) {
562 i = (long)(floor(x));
563 j = (long)(floor(y));
564 k = (long)(floor(z));
565 }
566 else {
567 i = (long)(floor(x + 0.5));
568 j = (long)(floor(y + 0.5));
569 k = (long)(floor(z + 0.5));
570 }
571
572 s = x - i;
573 t = y - j;
574 r = z - k;
575
576 i -= CentralIndex;
577 j -= CentralIndex;
578 k -= CentralIndex;
579
580 for (l = 0L; l <= SplineDegree; l++) {
581 xIndex[l] = i++;
582 yIndex[l] = j++;
583 zIndex[l] = k++;
584 }
585
586 /* get the interpolation weights */
587 xWeight[0] = 1.0;
588 yWeight[0] = 1.0;
589 zWeight[0] = 1.0;
590
591 if (Width > 1) {
592 vtkImageBSplineInternals::
593 GetInterpolationWeights(xWeight, s, SplineDegree);
594 }
595 if (Height > 1) {
596 vtkImageBSplineInternals::
597 GetInterpolationWeights(yWeight, t, SplineDegree);
598 }
599 if (Slices > 1) {
600 vtkImageBSplineInternals::
601 GetInterpolationWeights(zWeight, r, SplineDegree);
602 }
603
604 switch (Border)
605 {
606 case VTK_IMAGE_BORDER_CLAMP:
607 /* apply the constant boundary conditions */
608 for (l = 0L; l <= SplineDegree; l++) {
609 if (xIndex[l] < 0) {
610 xIndex[l] = 0;
611 }
612 else if (xIndex[l] >= Width) {
613 xIndex[l] = Width - 1;
614 }
615 }
616 for (l = 0L; l <= SplineDegree; l++) {
617 if (yIndex[l] < 0) {
618 yIndex[l] = 0;
619 }
620 else if (yIndex[l] >= Height) {
621 yIndex[l] = Height - 1;
622 }
623 }
624 for (l = 0L; l <= SplineDegree; l++) {
625 if (zIndex[l] < 0) {
626 zIndex[l] = 0;
627 }
628 else if (zIndex[l] >= Slices) {
629 zIndex[l] = Slices - 1;
630 }
631 }
632 break;
633
634 case VTK_IMAGE_BORDER_MIRROR:
635 /* apply the mirror boundary conditions */
636 for (l = 0L; l <= SplineDegree; l++) {
637 xIndex[l] = (Width == 1L) ? (0L) : ((xIndex[l] < 0L) ?
638 (-xIndex[l] - Width2 * ((-xIndex[l]) / Width2))
639 : (xIndex[l] - Width2 * (xIndex[l] / Width2)));
640 if (Width <= xIndex[l]) {
641 xIndex[l] = Width2 - xIndex[l];
642 }
643 yIndex[l] = (Height == 1L) ? (0L) : ((yIndex[l] < 0L) ?
644 (-yIndex[l] - Height2 * ((-yIndex[l]) / Height2))
645 : (yIndex[l] - Height2 * (yIndex[l] / Height2)));
646 if (Height <= yIndex[l]) {
647 yIndex[l] = Height2 - yIndex[l];
648 }
649 zIndex[l] = (Slices == 1L) ? (0L) : ((zIndex[l] < 0L) ?
650 (-zIndex[l] - Slices2 * ((-zIndex[l]) / Slices2))
651 : (zIndex[l] - Slices2 * (zIndex[l] / Slices2)));
652 if (Slices <= zIndex[l]) {
653 zIndex[l] = Slices2 - zIndex[l];
654 }
655 }
656 break;
657
658 case VTK_IMAGE_BORDER_REPEAT:
659 /* apply the repeat boundary conditions */
660 for (l = 0L; l <= SplineDegree; l++) {
661 if ((xIndex[l] %= Width) < 0) {
662 xIndex[l] += Width;
663 }
664 if ((yIndex[l] %= Height) < 0) {
665 yIndex[l] += Height;
666 }
667 if ((zIndex[l] %= Slices) < 0) {
668 zIndex[l] += Slices;
669 }
670 }
671 break;
672 }
673
674 /* precompute the increments in each direction */
675 xIncrement = Depth;
676 yIncrement = xIncrement * Width;
677 zIncrement = yIncrement * Height;
678
679 /* perform interpolation */
680 for (c = 0L; c < Depth; c++) {
681 p1 = Bcoeff + c;
682 interpolated = 0.0;
683 for (k = 0L; k <= kmax; k++) {
684 p2 = p1 + (zIndex[k] * zIncrement);
685 u = 0.0;
686 for (j = 0L; j <= jmax; j++) {
687 p3 = p2 + (yIndex[j] * yIncrement);
688 w = 0.0;
689 for (i = 0L; i <= imax; i++) {
690 w += xWeight[i] * p3[xIndex[i] * xIncrement];
691 }
692 u += yWeight[j] * w;
693 }
694 interpolated += zWeight[k] * u;
695 }
696 v[c] = interpolated;
697 }
698
699 return(1);
700 } /* end InterpolateValue */
701
702 }
703
704 /*--------------------------------------------------------------------------*/
GetInterpolationWeights(float weights[10],double w,long degree)705 int vtkImageBSplineInternals::GetInterpolationWeights(
706 float weights[10], double w, long degree)
707 {
708 return vtkImageBSplineGetInterpolationWeights(weights, w, degree);
709 }
710
GetInterpolationWeights(double weights[10],double w,long degree)711 int vtkImageBSplineInternals::GetInterpolationWeights(
712 double weights[10], double w, long degree)
713 {
714 return vtkImageBSplineGetInterpolationWeights(weights, w, degree);
715 }
716
717 /*--------------------------------------------------------------------------*/
InterpolatedValue(const double * coeffs,double * value,long width,long height,long slices,long depth,double x,double y,double z,long degree,long border)718 int vtkImageBSplineInternals::InterpolatedValue(
719 const double *coeffs, double *value,
720 long width, long height, long slices, long depth,
721 double x, double y, double z, long degree, long border)
722 {
723 return vtkImageBSplineInterpolatedValue(
724 coeffs, value, width, height, slices, depth, x, y, z, degree, border);
725 }
726
InterpolatedValue(const float * coeffs,float * value,long width,long height,long slices,long depth,double x,double y,double z,long degree,long border)727 int vtkImageBSplineInternals::InterpolatedValue(
728 const float *coeffs, float *value,
729 long width, long height, long slices, long depth,
730 double x, double y, double z, long degree, long border)
731 {
732 return vtkImageBSplineInterpolatedValue(
733 coeffs, value, width, height, slices, depth, x, y, z, degree, border);
734 }
735