BSpline.h
1 /*
2  * Medical Image Registration ToolKit (MIRTK)
3  *
4  * Copyright 2013-2017 Imperial College London
5  * Copyright 2013-2018 Andreas Schuh
6  *
7  * Licensed under the Apache License, Version 2.0 (the "License");
8  * you may not use this file except in compliance with the License.
9  * You may obtain a copy of the License at
10  *
11  * http://www.apache.org/licenses/LICENSE-2.0
12  *
13  * Unless required by applicable law or agreed to in writing, software
14  * distributed under the License is distributed on an "AS IS" BASIS,
15  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16  * See the License for the specific language governing permissions and
17  * limitations under the License.
18  */
19 
20 #ifndef MIRTK_BSpline_H
21 #define MIRTK_BSpline_H
22 
23 #include "mirtk/Config.h"
24 #include "mirtk/Math.h"
25 #include "mirtk/Stream.h"
26 
27 #include "mirtk/NumericsExport.h"
28 
29 
30 namespace mirtk {
31 
32 
33 /**
34  * Cubic B-spline, its basis functions, and their derivatives.
35  *
36  * This static class provides methods to evaluate the cubic B-spline function,
37  * each of its four basis functions, and the corresponding first and second
38  * derivatives in the interval [0, 1]. It is in particular intended for use
39  * by classes implementing the irtkFreeFormTransformation interface using
40  * a cubic B-spline basis to represent the displacement or velocity field,
41  * respectively.
42  *
43  * The cubic B-spline pieces are indexed here in reverse order!
44  * This is especially important when using the derivatives, because otherwise
45  * the sign may be incorrect when the intervals are exchanged.
46  * Spline function B(t) consists of pieces:
47  * - B3(t + 2) for t in [-2, -1)
48  * - B2(t + 1) for t in [-1, 0)
49  * - B1(t - 0) for t in [ 0, 1)
50  * - B0(t - 1) for t in [ 1, 2]
51  *
52  * \note Though not explicitly enforced by protecting these members,
53  * do <b>not</b> modify the lookup tables of precomputed function
54  * values as it will affect all code which makes use of them!
55  * Moreover, make sure to call the static Initialize() method at least
56  * once in your program before accessing any of these precomputed values.
57  *
58  * Usage:
59  * \code
60  * BSplineFunction kernel;
61  * // evaluate function explicitly
62  * double exact = kernel.B1(0.2);
63  * // or use lookup table; note that the Initialize() method must only be
64  * // called once during the life-time of the program
65  * BSplineFunction::Initialize();
66  * double approx = kernel.LookupTable[kernel.VariableToIndex(0.2)][1];
67  * \endcode
68  */
69 
70 template <class TReal>
71 class BSpline
72 {
73 public:
74 
75  // ---------------------------------------------------------------------------
76  // Generel 1D interpolation kernel interface
77 
78  /// Floating point precision type used
79  typedef TReal RealType;
80 
81  /// Lookup table of B-spline basis function values for t=0
82  /// (i.e., kernel centered at lattice coordinate)
83  static const TReal LatticeWeights[4];
84 
85  /// Lookup table of B-spline basis function 1st derivative values for t=0
86  /// (i.e., kernel centered at lattice coordinate)
87  static const TReal LatticeWeights_I[4];
88 
89  /// Lookup table of B-spline basis function 2nd derivative values for t=0
90  /// (i.e., kernel centered at lattice coordinate)
91  static const TReal LatticeWeights_II[4];
92 
93  /// Returns the value of the B-spline function
94  MIRTKCU_API static TReal Weight(TReal);
95 
96  /// Returns the values of the B-spline function for
97  /// - [0]: t in ]-2, -1[
98  /// - [1]: t in [-1, 0]
99  /// - [2]: t in ]0, 1]
100  /// - [3]: t in ]1, 2[
101  MIRTKCU_API static void Weights(TReal, TReal[4]);
102 
103  /// Size of lookup tables of pre-computed B-spline values
104  static const unsigned int LookupTableSize = 1000;
105 
106  /// Initialize lookup tables
107  static void Initialize(bool = false);
108 
109  /// Returns the lookup table index for a given value in [0,1]
110  static int VariableToIndex(TReal);
111 
112  /// Returns the lookup table indices for any t value
113  ///
114  /// \param[in] t Cubic B-spline function parameter.
115  /// \param[out] i First lookup table index in [0, LookupTableSize - 1].
116  /// \param[out] j Second lookup table index in [0, 3].
117  ///
118  /// \returns Whether t is within local support of cubic B-spline function.
119  /// When the return value is false, lookup table indices corresponding
120  /// to a boundary of the local support region is returned.
121  static bool VariableToIndex(TReal, int &i, int &j);
122 
123  /// Lookup table of B-spline function values
124  MIRTK_Numerics_EXPORT static TReal WeightLookupTable[LookupTableSize];
125 
126  /// Lookup table of B-spline basis function values
127  MIRTK_Numerics_EXPORT static TReal LookupTable[LookupTableSize][4];
128 
129  /// Lookup table of B-spline basis function 1st derivative values
130  MIRTK_Numerics_EXPORT static TReal LookupTable_I[LookupTableSize][4];
131 
132  /// Lookup table of B-spline basis function 2nd derivative values
133  MIRTK_Numerics_EXPORT static TReal LookupTable_II[LookupTableSize][4];
134 
135  // ---------------------------------------------------------------------------
136  // B-spline basis functions
137 
138  /// Returns the value of the B-spline function
139  MIRTKCU_API static TReal B(TReal);
140 
141  /// Returns the value of the i-th B-spline basis function
142  MIRTKCU_API static TReal B(int, TReal);
143 
144  /// Returns the value of the first B-spline basis function
145  MIRTKCU_API static TReal B0(TReal);
146 
147  /// Returns the value of the second B-spline basis function
148  MIRTKCU_API static TReal B1(TReal);
149 
150  /// Returns the value of the third B-spline basis function
151  MIRTKCU_API static TReal B2(TReal);
152 
153  /// Returns the value of the fourth B-spline basis function
154  MIRTKCU_API static TReal B3(TReal);
155 
156  /// Returns the 1st derivative value of the B-spline function
157  MIRTKCU_API static TReal B_I(TReal);
158 
159  /// Returns the 1st derivative value of the i-th B-spline basis function
160  MIRTKCU_API static TReal B_I(int, TReal);
161 
162  /// Returns the 1st derivative value of the first B-spline basis function
163  MIRTKCU_API static TReal B0_I(TReal);
164 
165  /// Returns the 1st derivative value of the second B-spline basis function
166  MIRTKCU_API static TReal B1_I(TReal);
167 
168  /// Returns the 1st derivative value of the third B-spline basis function
169  MIRTKCU_API static TReal B2_I(TReal);
170 
171  /// Returns the 1st derivative value of the fourth B-spline basis function
172  MIRTKCU_API static TReal B3_I(TReal);
173 
174  /// Returns the 2nd derivative value of the B-spline function
175  MIRTKCU_API static TReal B_II(TReal);
176 
177  /// Returns the 2nd derivative value of the i-th B-spline basis function
178  MIRTKCU_API static TReal B_II(int, TReal);
179 
180  /// Returns the 2nd derivative value of the first B-spline basis function
181  MIRTKCU_API static TReal B0_II(TReal);
182 
183  /// Returns the 2nd derivative value of the second B-spline basis function
184  MIRTKCU_API static TReal B1_II(TReal);
185 
186  /// Returns the 2nd derivative value of the third B-spline basis function
187  MIRTKCU_API static TReal B2_II(TReal);
188 
189  /// Returns the 2nd derivative value of the fourth B-spline basis function
190  MIRTKCU_API static TReal B3_II(TReal);
191 
192  /// Returns the 3rd derivative value of the B-spline function
193  MIRTKCU_API static TReal B_III(TReal);
194 
195  /// Returns the 3rd derivative value of the i-th B-spline basis function
196  MIRTKCU_API static TReal B_III(int, TReal);
197 
198  /// Returns the 3rd derivative value of the first B-spline basis function
199  MIRTKCU_API static TReal B0_III(TReal);
200 
201  /// Returns the 3rd derivative value of the second B-spline basis function
202  MIRTKCU_API static TReal B1_III(TReal);
203 
204  /// Returns the 3rd derivative value of the third B-spline basis function
205  MIRTKCU_API static TReal B2_III(TReal);
206 
207  /// Returns the 3rd derivative value of the fourth B-spline basis function
208  MIRTKCU_API static TReal B3_III(TReal);
209 
210  /// Returns the n-th derivative value of the B-spline function
211  MIRTKCU_API static TReal B_nI(int, TReal);
212 
213  /// Returns the n-th derivative value of the i-th B-spline basis function
214  MIRTKCU_API static TReal B_nI(int, int, TReal);
215 
216  /// Returns the n-th derivative value of the first B-spline basis function
217  MIRTKCU_API static TReal B0_nI(int, TReal);
218 
219  /// Returns the n-th derivative value of the second B-spline basis function
220  MIRTKCU_API static TReal B1_nI(int, TReal);
221 
222  /// Returns the n-th derivative value of the third B-spline basis function
223  MIRTKCU_API static TReal B2_nI(int, TReal);
224 
225  /// Returns the n-th derivative value of the fourth B-spline basis function
226  MIRTKCU_API static TReal B3_nI(int, TReal);
227 
228 protected:
229 
230  /// Flag which indicates whether the lookup tables are initialized
231  MIRTK_Numerics_EXPORT static bool _initialized;
232 };
233 
234 // -----------------------------------------------------------------------------
235 // Fix Clang warning regarding missing definition of static template members
236 #ifdef __clang__
237 
238 // Lookup table of B-spline function values
239 template <class TReal> TReal BSpline<TReal>::WeightLookupTable[BSpline<TReal>::LookupTableSize];
240 
241 // Lookup table of B-spline basis function values
242 template <class TReal> TReal BSpline<TReal>::LookupTable [BSpline<TReal>::LookupTableSize][4];
243 template <class TReal> TReal BSpline<TReal>::LookupTable_I [BSpline<TReal>::LookupTableSize][4];
244 template <class TReal> TReal BSpline<TReal>::LookupTable_II[BSpline<TReal>::LookupTableSize][4];
245 
246 // Wether lookup tables of B-spline kernel were initialized
247 template <class TReal> bool BSpline<TReal>::_initialized = false;
248 
249 #endif // defined(__clang__)
250 
251 // -----------------------------------------------------------------------------
252 // Weights often used by B-spline transformations for evaluation of derivatives
253 // at lattice points only where B-spline parameter t is zero
254 template <class TReal>
255 const TReal BSpline<TReal>::LatticeWeights[4] = {
256  TReal(1.0/6.0), TReal(2.0/3.0), TReal(1.0/6.0), TReal(0.0)
257 };
258 
259 template <class TReal>
260 const TReal BSpline<TReal>::LatticeWeights_I[4] = {
261  TReal(-0.5), TReal(0.0), TReal(0.5), TReal(0.0)
262 };
263 
264 template <class TReal>
265 const TReal BSpline<TReal>::LatticeWeights_II[4] = {
266  TReal(1.0), TReal(-2.0), TReal(1.0), TReal(0.0)
267 };
268 
269 // -----------------------------------------------------------------------------
270 typedef BSpline<double> BSplineFunction;
271 
272 ////////////////////////////////////////////////////////////////////////////////
273 // Inline definitions
274 ////////////////////////////////////////////////////////////////////////////////
275 
276 // -----------------------------------------------------------------------------
277 template <class TReal>
278 inline void BSpline<TReal>::Weights(TReal t, TReal value[4])
279 {
280  TReal c = static_cast<TReal>(1.0 - t);
281  TReal b = t * t;
282  TReal a = b * t;
283  value[0] = static_cast<TReal>(c * c * c / 6.0);
284  value[1] = static_cast<TReal>(( 3.0 * a - 6.0 * b + 4.0) / 6.0);
285  value[2] = static_cast<TReal>((-3.0 * a + 3.0 * b + 3.0 * t + 1.0) / 6.0);
286  value[3] = static_cast<TReal>(a / 6.0);
287 }
288 
289 // -----------------------------------------------------------------------------
290 template <class TReal>
291 inline TReal BSpline<TReal>::Weight(TReal t)
292 {
293  double value = 0.0;
294  t = fabs(t);
295  if (t < 2.0) {
296  if (t < 1.0) {
297  value = 2.0/3.0 + (0.5*t-1.0)*t*t; // B1
298  } else {
299  t -= 2.0;
300  value = -t*t*t/6.0; // B3
301  }
302  }
303  return static_cast<TReal>(value);
304 }
305 
306 // -----------------------------------------------------------------------------
307 template <class TReal>
308 inline TReal BSpline<TReal>::B(TReal t)
309 {
310  return Weight(t);
311 }
312 
313 // -----------------------------------------------------------------------------
314 template <class TReal>
315 inline TReal BSpline<TReal>::B(int i, TReal t)
316 {
317  switch (i) {
318  case 0: return B0(t);
319  case 1: return B1(t);
320  case 2: return B2(t);
321  case 3: return B3(t);
322  default: return TReal(0);
323  }
324 }
325 
326 // -----------------------------------------------------------------------------
327 template <class TReal>
328 inline TReal BSpline<TReal>::B0(TReal t)
329 {
330  return static_cast<TReal>((1.0-t)*(1.0-t)*(1.0-t)/6.0);
331 }
332 
333 // -----------------------------------------------------------------------------
334 template <class TReal>
335 inline TReal BSpline<TReal>::B1(TReal t)
336 {
337  return static_cast<TReal>((3.0*t*t*t - 6.0*t*t + 4.0)/6.0);
338 }
339 
340 // -----------------------------------------------------------------------------
341 template <class TReal>
342 inline TReal BSpline<TReal>::B2(TReal t)
343 {
344  return static_cast<TReal>((-3.0*t*t*t + 3.0*t*t + 3.0*t + 1.0)/6.0);
345 }
346 
347 // -----------------------------------------------------------------------------
348 template <class TReal>
349 inline TReal BSpline<TReal>::B3(TReal t)
350 {
351  return static_cast<TReal>((t*t*t)/6.0);
352 }
353 
354 // -----------------------------------------------------------------------------
355 template <class TReal>
356 inline TReal BSpline<TReal>::B_I(TReal t)
357 {
358  double y = abs(t);
359  double value = 0.0;
360  if (y < 2.0) {
361  if (y < 1.0) {
362  value = (1.5 * y - 2.0) * t;
363  } else {
364  y -= 2.0;
365  value = -0.5 * y * y;
366  if (t < 0.0) value = -value;
367  }
368  }
369  return static_cast<TReal>(value);
370 }
371 
372 // -----------------------------------------------------------------------------
373 template <class TReal>
374 inline TReal BSpline<TReal>::B_I(int i, TReal t)
375 {
376  switch (i) {
377  case 0: return B0_I(t);
378  case 1: return B1_I(t);
379  case 2: return B2_I(t);
380  case 3: return B3_I(t);
381  default: return TReal(0);
382  }
383 }
384 
385 // -----------------------------------------------------------------------------
386 template <class TReal>
387 inline TReal BSpline<TReal>::B0_I(TReal t)
388 {
389  return static_cast<TReal>(-(1.0-t)*(1.0-t)/2.0);
390 }
391 
392 // -----------------------------------------------------------------------------
393 template <class TReal>
394 inline TReal BSpline<TReal>::B1_I(TReal t)
395 {
396  return static_cast<TReal>((9.0*t*t - 12.0*t)/6.0);
397 }
398 
399 // -----------------------------------------------------------------------------
400 template <class TReal>
401 inline TReal BSpline<TReal>::B2_I(TReal t)
402 {
403  return static_cast<TReal>((-9.0*t*t + 6.0*t + 3.0)/6.0);
404 }
405 
406 // -----------------------------------------------------------------------------
407 template <class TReal>
408 inline TReal BSpline<TReal>::B3_I(TReal t)
409 {
410  return static_cast<TReal>((t*t)/2.0);
411 }
412 
413 // -----------------------------------------------------------------------------
414 template <class TReal>
415 inline TReal BSpline<TReal>::B_II(TReal t)
416 {
417  double value = 0.0;
418  t = fabs(t);
419  if (t < 2.0) {
420  if (t < 1.0) {
421  value = 3.0 * t - 2.0;
422  } else {
423  value = -t + 2.0;
424  }
425  }
426  return static_cast<TReal>(value);
427 }
428 
429 // -----------------------------------------------------------------------------
430 template <class TReal>
431 inline TReal BSpline<TReal>::B_II(int i, TReal t)
432 {
433  switch (i) {
434  case 0: return B0_II(t);
435  case 1: return B1_II(t);
436  case 2: return B2_II(t);
437  case 3: return B3_II(t);
438  default: return TReal(0);
439  }
440 }
441 
442 // -----------------------------------------------------------------------------
443 template <class TReal>
444 inline TReal BSpline<TReal>::B0_II(TReal t)
445 {
446  return static_cast<TReal>(1.0 - t);
447 }
448 
449 // -----------------------------------------------------------------------------
450 template <class TReal>
451 inline TReal BSpline<TReal>::B1_II(TReal t)
452 {
453  return static_cast<TReal>(3.0*t - 2.0);
454 }
455 
456 // -----------------------------------------------------------------------------
457 template <class TReal>
458 inline TReal BSpline<TReal>::B2_II(TReal t)
459 {
460  return static_cast<TReal>(-3.0*t + 1.0);
461 }
462 
463 // -----------------------------------------------------------------------------
464 template <class TReal>
465 inline TReal BSpline<TReal>::B3_II(TReal t)
466 {
467  return t;
468 }
469 
470 // -----------------------------------------------------------------------------
471 template <class TReal>
472 inline TReal BSpline<TReal>::B_III(TReal t)
473 {
474  TReal y = fabs(t), value;
475  if (y < 2.0) {
476  if (y < 1.0) {
477  value = copysign(TReal(3.0), t);
478  } else {
479  value = copysign(TReal(1.0), -t);
480  }
481  } else {
482  value = TReal(0);
483  }
484  return value;
485 }
486 
487 // -----------------------------------------------------------------------------
488 template <class TReal>
489 inline TReal BSpline<TReal>::B_III(int i, TReal t)
490 {
491  switch (i) {
492  case 0: return B0_III(t);
493  case 1: return B1_III(t);
494  case 2: return B2_III(t);
495  case 3: return B3_III(t);
496  default: return TReal(0);
497  }
498 }
499 
500 // -----------------------------------------------------------------------------
501 template <class TReal>
502 inline TReal BSpline<TReal>::B0_III(TReal)
503 {
504  return TReal(-1.0);
505 }
506 
507 // -----------------------------------------------------------------------------
508 template <class TReal>
509 inline TReal BSpline<TReal>::B1_III(TReal)
510 {
511  return TReal(3.0);
512 }
513 
514 // -----------------------------------------------------------------------------
515 template <class TReal>
516 inline TReal BSpline<TReal>::B2_III(TReal)
517 {
518  return TReal(-3.0);
519 }
520 
521 // -----------------------------------------------------------------------------
522 template <class TReal>
523 inline TReal BSpline<TReal>::B3_III(TReal)
524 {
525  return TReal(1.0);
526 }
527 
528 // -----------------------------------------------------------------------------
529 template <class TReal>
530 inline TReal BSpline<TReal>::B_nI(int n, TReal t)
531 {
532  switch(n) {
533  case 0: return B (t);
534  case 1: return B_I (t);
535  case 2: return B_II (t);
536  case 3: return B_III(t);
537  default: return TReal(0);
538  }
539 }
540 
541 // -----------------------------------------------------------------------------
542 template <class TReal>
543 inline TReal BSpline<TReal>::B_nI(int i, int n, TReal t)
544 {
545  switch(n) {
546  case 0: return B (i, t);
547  case 1: return B_I (i, t);
548  case 2: return B_II (i, t);
549  case 3: return B_III(i, t);
550  default: return TReal(0);
551  }
552 }
553 
554 // -----------------------------------------------------------------------------
555 template <class TReal>
556 inline TReal BSpline<TReal>::B0_nI(int n, TReal t)
557 {
558  switch(n) {
559  case 0: return B0 (t);
560  case 1: return B0_I (t);
561  case 2: return B0_II (t);
562  case 3: return B0_III(t);
563  default: return TReal(0);
564  }
565 }
566 
567 // -----------------------------------------------------------------------------
568 template <class TReal>
569 inline TReal BSpline<TReal>::B1_nI(int n, TReal t)
570 {
571  switch(n) {
572  case 0: return B1 (t);
573  case 1: return B1_I (t);
574  case 2: return B1_II (t);
575  case 3: return B1_III(t);
576  default: return TReal(0);
577  }
578 }
579 
580 // -----------------------------------------------------------------------------
581 template <class TReal>
582 inline TReal BSpline<TReal>::B2_nI(int n, TReal t)
583 {
584  switch(n) {
585  case 0: return B2 (t);
586  case 1: return B2_I (t);
587  case 2: return B2_II (t);
588  case 3: return B2_III(t);
589  default: return TReal(0);
590  }
591 }
592 
593 // -----------------------------------------------------------------------------
594 template <class TReal>
595 inline TReal BSpline<TReal>::B3_nI(int n, TReal t)
596 {
597  switch(n) {
598  case 0: return B3 (t);
599  case 1: return B3_I (t);
600  case 2: return B3_II (t);
601  case 3: return B3_III(t);
602  default: return TReal(0);
603  }
604 }
605 
606 // -----------------------------------------------------------------------------
607 template <class TReal>
608 inline void BSpline<TReal>::Initialize(bool force)
609 {
610  if (!_initialized || force) {
611  for (unsigned int i = 0; i < LookupTableSize; i++) {
612  WeightLookupTable[i] = Weight(static_cast<TReal>(i) / static_cast<TReal>(LookupTableSize-1));
613  for (int l = 0; l < 4; l++) {
614  LookupTable [i][l] = B (l, static_cast<TReal>(i) / static_cast<TReal>(LookupTableSize-1));
615  LookupTable_I [i][l] = B_I (l, static_cast<TReal>(i) / static_cast<TReal>(LookupTableSize-1));
616  LookupTable_II[i][l] = B_II(l, static_cast<TReal>(i) / static_cast<TReal>(LookupTableSize-1));
617  }
618  }
619  _initialized = true;
620  }
621 }
622 
623 // -----------------------------------------------------------------------------
624 template <class TReal>
625 inline int BSpline<TReal>::VariableToIndex(TReal t)
626 {
627  if (t < .0) return 0;
628  else if (t > 1.0) return LookupTableSize-1;
629  else return iround(t * static_cast<TReal>(LookupTableSize-1));
630 }
631 
632 // -----------------------------------------------------------------------------
633 template <class TReal>
634 inline bool BSpline<TReal>::VariableToIndex(TReal t, int &i, int &j)
635 {
636  if (t < 0.) {
637  if (t <= -2.) {
638  j = 3;
639  i = 0;
640  return false;
641  }
642  if (t < -1.) {
643  j = 3;
644  i = iround((t + 2.) * static_cast<TReal>(LookupTableSize - 1));
645  } else {
646  j = 2;
647  i = iround((t + 1.) * static_cast<TReal>(LookupTableSize - 1));
648  }
649  } else {
650  if (t >= 2.) {
651  j = 0;
652  i = LookupTableSize - 1;
653  return false;
654  }
655  if (t < 1.) {
656  j = 1;
657  i = iround(t * static_cast<TReal>(LookupTableSize - 1));
658  } else {
659  j = 0;
660  i = iround((t - 1.) * static_cast<TReal>(LookupTableSize - 1));
661  }
662  }
663  return true;
664 }
665 
666 ////////////////////////////////////////////////////////////////////////////////
667 // B-spline interpolation weights
668 ////////////////////////////////////////////////////////////////////////////////
669 
670 // -----------------------------------------------------------------------------
671 /// Compute indices and weights required for 2D B-spline interpolation
672 template <class TReal>
673 inline void ComputeBSplineIndicesAndWeights(double x, double y, int spline_degree,
674  int xIndex [6], int yIndex [6],
675  TReal xWeight[6], TReal yWeight[6])
676 {
677  // Compute the interpolation indices
678  int i, j;
679  if (spline_degree & 1) {
680  i = static_cast<int>(floor(x )) - spline_degree / 2;
681  j = static_cast<int>(floor(y )) - spline_degree / 2;
682  } else {
683  i = static_cast<int>(floor(x + 0.5)) - spline_degree / 2;
684  j = static_cast<int>(floor(y + 0.5)) - spline_degree / 2;
685  }
686  for (int m = 0; m <= spline_degree; m++) {
687  xIndex[m] = i++;
688  yIndex[m] = j++;
689  }
690 
691  // Compute the interpolation weights
692  double w, w2, w4, t, t0, t1;
693 
694  switch (spline_degree) {
695  case 2:
696  /* x */
697  w = x - static_cast<double>(xIndex[1]);
698  xWeight[1] = TReal(3.0 / 4.0 - w * w);
699  xWeight[2] = TReal((1.0 / 2.0) * (w - xWeight[1] + 1.0));
700  xWeight[0] = TReal(1.0 - xWeight[1] - xWeight[2]);
701  /* y */
702  w = y - static_cast<double>(yIndex[1]);
703  yWeight[1] = TReal(3.0 / 4.0 - w * w);
704  yWeight[2] = TReal((1.0 / 2.0) * (w - yWeight[1] + 1.0));
705  yWeight[0] = TReal(1.0 - yWeight[1] - yWeight[2]);
706  break;
707  case 3:
708  /* x */
709  w = x - static_cast<double>(xIndex[1]);
710  xWeight[3] = TReal((1.0 / 6.0) * w * w * w);
711  xWeight[0] = TReal((1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - xWeight[3]);
712  xWeight[2] = TReal(w + xWeight[0] - 2.0 * xWeight[3]);
713  xWeight[1] = TReal(1.0 - xWeight[0] - xWeight[2] - xWeight[3]);
714  /* y */
715  w = y - static_cast<double>(yIndex[1]);
716  yWeight[3] = TReal((1.0 / 6.0) * w * w * w);
717  yWeight[0] = TReal((1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - yWeight[3]);
718  yWeight[2] = TReal(w + yWeight[0] - 2.0 * yWeight[3]);
719  yWeight[1] = TReal(1.0 - yWeight[0] - yWeight[2] - yWeight[3]);
720  break;
721  case 4:
722  /* x */
723  w = x - static_cast<double>(xIndex[2]);
724  w2 = w * w;
725  t = (1.0 / 6.0) * w2;
726  xWeight[0] = TReal(1.0 / 2.0 - w);
727  xWeight[0] *= xWeight[0];
728  xWeight[0] *= TReal(1.0 / 24.0) * xWeight[0];
729  t0 = w * (t - 11.0 / 24.0);
730  t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
731  xWeight[1] = TReal(t1 + t0);
732  xWeight[3] = TReal(t1 - t0);
733  xWeight[4] = TReal(xWeight[0] + t0 + (1.0 / 2.0) * w);
734  xWeight[2] = TReal(1.0 - xWeight[0] - xWeight[1] - xWeight[3] - xWeight[4]);
735  /* y */
736  w = y - static_cast<double>(yIndex[2]);
737  w2 = w * w;
738  t = (1.0 / 6.0) * w2;
739  yWeight[0] = TReal(1.0 / 2.0 - w);
740  yWeight[0] *= yWeight[0];
741  yWeight[0] *= TReal(1.0 / 24.0) * yWeight[0];
742  t0 = w * (t - 11.0 / 24.0);
743  t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
744  yWeight[1] = TReal(t1 + t0);
745  yWeight[3] = TReal(t1 - t0);
746  yWeight[4] = TReal(yWeight[0] + t0 + (1.0 / 2.0) * w);
747  yWeight[2] = TReal(1.0 - yWeight[0] - yWeight[1] - yWeight[3] - yWeight[4]);
748  break;
749  case 5:
750  /* x */
751  w = x - static_cast<double>(xIndex[2]);
752  w2 = w * w;
753  xWeight[5] = TReal((1.0 / 120.0) * w * w2 * w2);
754  w2 -= w;
755  w4 = w2 * w2;
756  w -= 1.0 / 2.0;
757  t = w2 * (w2 - 3.0);
758  xWeight[0] = TReal((1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - xWeight[5]);
759  t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
760  t1 = (-1.0 / 12.0) * w * (t + 4.0);
761  xWeight[2] = TReal(t0 + t1);
762  xWeight[3] = TReal(t0 - t1);
763  t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
764  t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
765  xWeight[1] = TReal(t0 + t1);
766  xWeight[4] = TReal(t0 - t1);
767  /* y */
768  w = y - static_cast<double>(yIndex[2]);
769  w2 = w * w;
770  yWeight[5] = TReal((1.0 / 120.0) * w * w2 * w2);
771  w2 -= w;
772  w4 = w2 * w2;
773  w -= 1.0 / 2.0;
774  t = w2 * (w2 - 3.0);
775  yWeight[0] = TReal((1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - yWeight[5]);
776  t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
777  t1 = (-1.0 / 12.0) * w * (t + 4.0);
778  yWeight[2] = TReal(t0 + t1);
779  yWeight[3] = TReal(t0 - t1);
780  t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
781  t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
782  yWeight[1] = TReal(t0 + t1);
783  yWeight[4] = TReal(t0 - t1);
784  break;
785  default:
786  cerr << "ComputeBSplineIndicesAndWeights: Unsupported B-spline degree: " << spline_degree << endl;
787  exit(1);
788  }
789 }
790 
791 // -----------------------------------------------------------------------------
792 /// Compute indices and weights required for 3D B-spline interpolation
793 template <class TReal>
794 inline void ComputeBSplineIndicesAndWeights(double x, double y, double z, int spline_degree,
795  int xIndex [6], int yIndex [6], int zIndex [6],
796  TReal xWeight[6], TReal yWeight[6], TReal zWeight[6])
797 {
798  // Compute the interpolation indices
799  int i, j, k, m;
800 
801  if (spline_degree & 1) {
802  i = static_cast<int>(floor(x)) - spline_degree / 2;
803  j = static_cast<int>(floor(y)) - spline_degree / 2;
804  k = static_cast<int>(floor(z)) - spline_degree / 2;
805  for (m = 0; m <= spline_degree; m++) {
806  xIndex[m] = i++;
807  yIndex[m] = j++;
808  zIndex[m] = k++;
809  }
810  } else {
811  i = static_cast<int>(floor(x + 0.5)) - spline_degree / 2;
812  j = static_cast<int>(floor(y + 0.5)) - spline_degree / 2;
813  k = static_cast<int>(floor(z + 0.5)) - spline_degree / 2;
814  for (m = 0; m <= spline_degree; m++) {
815  xIndex[m] = i++;
816  yIndex[m] = j++;
817  zIndex[m] = k++;
818  }
819  }
820 
821  // Compute the interpolation weights
822  double w, w2, w4, t, t0, t1;
823 
824  switch (spline_degree) {
825  case 2:
826  /* x */
827  w = x - static_cast<double>(xIndex[1]);
828  xWeight[1] = TReal(3.0 / 4.0 - w * w);
829  xWeight[2] = TReal((1.0 / 2.0) * (w - xWeight[1] + 1.0));
830  xWeight[0] = TReal(1.0 - xWeight[1] - xWeight[2]);
831  /* y */
832  w = y - static_cast<double>(yIndex[1]);
833  yWeight[1] = TReal(3.0 / 4.0 - w * w);
834  yWeight[2] = TReal((1.0 / 2.0) * (w - yWeight[1] + 1.0));
835  yWeight[0] = TReal(1.0 - yWeight[1] - yWeight[2]);
836  /* z */
837  w = z - static_cast<double>(zIndex[1]);
838  zWeight[1] = TReal(3.0 / 4.0 - w * w);
839  zWeight[2] = TReal((1.0 / 2.0) * (w - zWeight[1] + 1.0));
840  zWeight[0] = TReal(1.0 - zWeight[1] - zWeight[2]);
841  break;
842  case 3:
843  /* x */
844  w = x - static_cast<double>(xIndex[1]);
845  xWeight[3] = TReal((1.0 / 6.0) * w * w * w);
846  xWeight[0] = TReal((1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - xWeight[3]);
847  xWeight[2] = TReal(w + xWeight[0] - 2.0 * xWeight[3]);
848  xWeight[1] = TReal(1.0 - xWeight[0] - xWeight[2] - xWeight[3]);
849  /* y */
850  w = y - static_cast<double>(yIndex[1]);
851  yWeight[3] = TReal((1.0 / 6.0) * w * w * w);
852  yWeight[0] = TReal((1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - yWeight[3]);
853  yWeight[2] = TReal(w + yWeight[0] - 2.0 * yWeight[3]);
854  yWeight[1] = TReal(1.0 - yWeight[0] - yWeight[2] - yWeight[3]);
855  /* z */
856  w = z - static_cast<double>(zIndex[1]);
857  zWeight[3] = TReal((1.0 / 6.0) * w * w * w);
858  zWeight[0] = TReal((1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - zWeight[3]);
859  zWeight[2] = TReal(w + zWeight[0] - 2.0 * zWeight[3]);
860  zWeight[1] = TReal(1.0 - zWeight[0] - zWeight[2] - zWeight[3]);
861  break;
862  case 4:
863  /* x */
864  w = x - static_cast<double>(xIndex[2]);
865  w2 = w * w;
866  t = (1.0 / 6.0) * w2;
867  xWeight[0] = TReal(1.0 / 2.0 - w);
868  xWeight[0] *= xWeight[0];
869  xWeight[0] *= TReal(1.0 / 24.0) * xWeight[0];
870  t0 = w * (t - 11.0 / 24.0);
871  t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
872  xWeight[1] = TReal(t1 + t0);
873  xWeight[3] = TReal(t1 - t0);
874  xWeight[4] = TReal(xWeight[0] + t0 + (1.0 / 2.0) * w);
875  xWeight[2] = TReal(1.0 - xWeight[0] - xWeight[1] - xWeight[3] - xWeight[4]);
876  /* y */
877  w = y - static_cast<double>(yIndex[2]);
878  w2 = w * w;
879  t = (1.0 / 6.0) * w2;
880  yWeight[0] = TReal(1.0 / 2.0 - w);
881  yWeight[0] *= yWeight[0];
882  yWeight[0] *= TReal(1.0 / 24.0) * yWeight[0];
883  t0 = w * (t - 11.0 / 24.0);
884  t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
885  yWeight[1] = TReal(t1 + t0);
886  yWeight[3] = TReal(t1 - t0);
887  yWeight[4] = TReal(yWeight[0] + t0 + (1.0 / 2.0) * w);
888  yWeight[2] = TReal(1.0 - yWeight[0] - yWeight[1] - yWeight[3] - yWeight[4]);
889  /* z */
890  w = z - static_cast<double>(zIndex[2]);
891  w2 = w * w;
892  t = (1.0 / 6.0) * w2;
893  zWeight[0] = TReal(1.0 / 2.0 - w);
894  zWeight[0] *= zWeight[0];
895  zWeight[0] *= TReal(1.0 / 24.0) * zWeight[0];
896  t0 = w * (t - 11.0 / 24.0);
897  t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
898  zWeight[1] = TReal(t1 + t0);
899  zWeight[3] = TReal(t1 - t0);
900  zWeight[4] = TReal(zWeight[0] + t0 + (1.0 / 2.0) * w);
901  zWeight[2] = TReal(1.0 - zWeight[0] - zWeight[1] - zWeight[3] - zWeight[4]);
902  break;
903  case 5:
904  /* x */
905  w = x - static_cast<double>(xIndex[2]);
906  w2 = w * w;
907  xWeight[5] = TReal((1.0 / 120.0) * w * w2 * w2);
908  w2 -= w;
909  w4 = w2 * w2;
910  w -= 1.0 / 2.0;
911  t = w2 * (w2 - 3.0);
912  xWeight[0] = TReal((1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - xWeight[5]);
913  t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
914  t1 = (-1.0 / 12.0) * w * (t + 4.0);
915  xWeight[2] = TReal(t0 + t1);
916  xWeight[3] = TReal(t0 - t1);
917  t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
918  t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
919  xWeight[1] = TReal(t0 + t1);
920  xWeight[4] = TReal(t0 - t1);
921  /* y */
922  w = y - static_cast<double>(yIndex[2]);
923  w2 = w * w;
924  yWeight[5] = TReal((1.0 / 120.0) * w * w2 * w2);
925  w2 -= w;
926  w4 = w2 * w2;
927  w -= 1.0 / 2.0;
928  t = w2 * (w2 - 3.0);
929  yWeight[0] = TReal((1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - yWeight[5]);
930  t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
931  t1 = (-1.0 / 12.0) * w * (t + 4.0);
932  yWeight[2] = TReal(t0 + t1);
933  yWeight[3] = TReal(t0 - t1);
934  t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
935  t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
936  yWeight[1] = TReal(t0 + t1);
937  yWeight[4] = TReal(t0 - t1);
938  /* z */
939  w = z - static_cast<double>(zIndex[2]);
940  w2 = w * w;
941  zWeight[5] = TReal((1.0 / 120.0) * w * w2 * w2);
942  w2 -= w;
943  w4 = w2 * w2;
944  w -= 1.0 / 2.0;
945  t = w2 * (w2 - 3.0);
946  zWeight[0] = TReal((1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - zWeight[5]);
947  t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
948  t1 = (-1.0 / 12.0) * w * (t + 4.0);
949  zWeight[2] = TReal(t0 + t1);
950  zWeight[3] = TReal(t0 - t1);
951  t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
952  t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
953  zWeight[1] = TReal(t0 + t1);
954  zWeight[4] = TReal(t0 - t1);
955  break;
956  default:
957  cerr << "ComputeBSplineIndicesAndWeights: Unsupported B-spline degree: " << spline_degree << endl;
958  exit(1);
959  }
960 }
961 
962 // -----------------------------------------------------------------------------
963 /// Compute indices and weights required for 4D B-spline interpolation
964 template <class TReal>
965 inline void ComputeBSplineIndicesAndWeights(double x, double y, double z, double t, int spline_degree,
966  int xIndex [6], int yIndex [6], int zIndex [6], int tIndex[6],
967  TReal xWeight[6], TReal yWeight[6], TReal zWeight[6], TReal tWeight[6])
968 {
969  // Compute the interpolation indices
970  int i, j, k, l, m;
971 
972  if (spline_degree & 1) {
973  i = ifloor(x) - spline_degree / 2;
974  j = ifloor(y) - spline_degree / 2;
975  k = ifloor(z) - spline_degree / 2;
976  l = ifloor(t) - spline_degree / 2;
977  for (m = 0; m <= spline_degree; m++) {
978  xIndex[m] = i++;
979  yIndex[m] = j++;
980  zIndex[m] = k++;
981  tIndex[m] = l++;
982  }
983  } else {
984  i = ifloor(x + .5) - spline_degree / 2;
985  j = ifloor(y + .5) - spline_degree / 2;
986  k = ifloor(z + .5) - spline_degree / 2;
987  l = ifloor(t + .5) - spline_degree / 2;
988  for (m = 0; m <= spline_degree; m++) {
989  xIndex[m] = i++;
990  yIndex[m] = j++;
991  zIndex[m] = k++;
992  tIndex[m] = l++;
993  }
994  }
995 
996  // Compute the interpolation weights
997  double w, w2, w4, t0, t1;
998 
999  switch (spline_degree) {
1000  case 2:
1001  /* x */
1002  w = x - static_cast<double>(xIndex[1]);
1003  xWeight[1] = TReal(3.0 / 4.0 - w * w);
1004  xWeight[2] = TReal((1.0 / 2.0) * (w - xWeight[1] + 1.0));
1005  xWeight[0] = TReal(1.0 - xWeight[1] - xWeight[2]);
1006  /* y */
1007  w = y - static_cast<double>(yIndex[1]);
1008  yWeight[1] = TReal(3.0 / 4.0 - w * w);
1009  yWeight[2] = TReal((1.0 / 2.0) * (w - yWeight[1] + 1.0));
1010  yWeight[0] = TReal(1.0 - yWeight[1] - yWeight[2]);
1011  /* z */
1012  w = z - static_cast<double>(zIndex[1]);
1013  zWeight[1] = TReal(3.0 / 4.0 - w * w);
1014  zWeight[2] = TReal((1.0 / 2.0) * (w - zWeight[1] + 1.0));
1015  zWeight[0] = TReal(1.0 - zWeight[1] - zWeight[2]);
1016  /* t */
1017  w = t - static_cast<double>(tIndex[1]);
1018  tWeight[1] = TReal(3.0 / 4.0 - w * w);
1019  tWeight[2] = TReal((1.0 / 2.0) * (w - tWeight[1] + 1.0));
1020  tWeight[0] = TReal(1.0 - tWeight[1] - tWeight[2]);
1021  break;
1022  case 3:
1023  /* x */
1024  w = x - static_cast<double>(xIndex[1]);
1025  xWeight[3] = TReal((1.0 / 6.0) * w * w * w);
1026  xWeight[0] = TReal((1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - xWeight[3]);
1027  xWeight[2] = TReal(w + xWeight[0] - 2.0 * xWeight[3]);
1028  xWeight[1] = TReal(1.0 - xWeight[0] - xWeight[2] - xWeight[3]);
1029  /* y */
1030  w = y - static_cast<double>(yIndex[1]);
1031  yWeight[3] = TReal((1.0 / 6.0) * w * w * w);
1032  yWeight[0] = TReal((1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - yWeight[3]);
1033  yWeight[2] = TReal(w + yWeight[0] - 2.0 * yWeight[3]);
1034  yWeight[1] = TReal(1.0 - yWeight[0] - yWeight[2] - yWeight[3]);
1035  /* z */
1036  w = z - static_cast<double>(zIndex[1]);
1037  zWeight[3] = TReal((1.0 / 6.0) * w * w * w);
1038  zWeight[0] = TReal((1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - zWeight[3]);
1039  zWeight[2] = TReal(w + zWeight[0] - 2.0 * zWeight[3]);
1040  zWeight[1] = TReal(1.0 - zWeight[0] - zWeight[2] - zWeight[3]);
1041  /* t */
1042  w = t - static_cast<double>(tIndex[1]);
1043  tWeight[3] = TReal((1.0 / 6.0) * w * w * w);
1044  tWeight[0] = TReal((1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - tWeight[3]);
1045  tWeight[2] = TReal(w + tWeight[0] - 2.0 * tWeight[3]);
1046  tWeight[1] = TReal(1.0 - tWeight[0] - tWeight[2] - tWeight[3]);
1047  break;
1048  case 4:
1049  /* x */
1050  w = x - static_cast<double>(xIndex[2]);
1051  w2 = w * w;
1052  t = (1.0 / 6.0) * w2;
1053  xWeight[0] = TReal(1.0 / 2.0 - w);
1054  xWeight[0] *= xWeight[0];
1055  xWeight[0] *= TReal(1.0 / 24.0) * xWeight[0];
1056  t0 = w * (t - 11.0 / 24.0);
1057  t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
1058  xWeight[1] = TReal(t1 + t0);
1059  xWeight[3] = TReal(t1 - t0);
1060  xWeight[4] = TReal(xWeight[0] + t0 + (1.0 / 2.0) * w);
1061  xWeight[2] = TReal(1.0 - xWeight[0] - xWeight[1] - xWeight[3] - xWeight[4]);
1062  /* y */
1063  w = y - static_cast<double>(yIndex[2]);
1064  w2 = w * w;
1065  t = (1.0 / 6.0) * w2;
1066  yWeight[0] = TReal(1.0 / 2.0 - w);
1067  yWeight[0] *= yWeight[0];
1068  yWeight[0] *= TReal(1.0 / 24.0) * yWeight[0];
1069  t0 = w * (t - 11.0 / 24.0);
1070  t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
1071  yWeight[1] = TReal(t1 + t0);
1072  yWeight[3] = TReal(t1 - t0);
1073  yWeight[4] = TReal(yWeight[0] + t0 + (1.0 / 2.0) * w);
1074  yWeight[2] = TReal(1.0 - yWeight[0] - yWeight[1] - yWeight[3] - yWeight[4]);
1075  /* z */
1076  w = z - static_cast<double>(zIndex[2]);
1077  w2 = w * w;
1078  t = (1.0 / 6.0) * w2;
1079  zWeight[0] = TReal(1.0 / 2.0 - w);
1080  zWeight[0] *= zWeight[0];
1081  zWeight[0] *= TReal(1.0 / 24.0) * zWeight[0];
1082  t0 = w * (t - 11.0 / 24.0);
1083  t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
1084  zWeight[1] = TReal(t1 + t0);
1085  zWeight[3] = TReal(t1 - t0);
1086  zWeight[4] = TReal(zWeight[0] + t0 + (1.0 / 2.0) * w);
1087  zWeight[2] = TReal(1.0 - zWeight[0] - zWeight[1] - zWeight[3] - zWeight[4]);
1088  /* t */
1089  w = t - static_cast<double>(tIndex[2]);
1090  w2 = w * w;
1091  t = (1.0 / 6.0) * w2;
1092  tWeight[0] = TReal(1.0 / 2.0 - w);
1093  tWeight[0] *= tWeight[0];
1094  tWeight[0] *= TReal(1.0 / 24.0) * tWeight[0];
1095  t0 = w * (t - 11.0 / 24.0);
1096  t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
1097  tWeight[1] = TReal(t1 + t0);
1098  tWeight[3] = TReal(t1 - t0);
1099  tWeight[4] = TReal(tWeight[0] + t0 + (1.0 / 2.0) * w);
1100  tWeight[2] = TReal(1.0 - tWeight[0] - tWeight[1] - tWeight[3] - tWeight[4]);
1101  break;
1102  case 5:
1103  /* x */
1104  w = x - static_cast<double>(xIndex[2]);
1105  w2 = w * w;
1106  xWeight[5] = TReal((1.0 / 120.0) * w * w2 * w2);
1107  w2 -= w;
1108  w4 = w2 * w2;
1109  w -= 1.0 / 2.0;
1110  t = w2 * (w2 - 3.0);
1111  xWeight[0] = TReal((1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - xWeight[5]);
1112  t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
1113  t1 = (-1.0 / 12.0) * w * (t + 4.0);
1114  xWeight[2] = TReal(t0 + t1);
1115  xWeight[3] = TReal(t0 - t1);
1116  t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
1117  t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
1118  xWeight[1] = TReal(t0 + t1);
1119  xWeight[4] = TReal(t0 - t1);
1120  /* y */
1121  w = y - static_cast<double>(yIndex[2]);
1122  w2 = w * w;
1123  yWeight[5] = TReal((1.0 / 120.0) * w * w2 * w2);
1124  w2 -= w;
1125  w4 = w2 * w2;
1126  w -= 1.0 / 2.0;
1127  t = w2 * (w2 - 3.0);
1128  yWeight[0] = TReal((1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - yWeight[5]);
1129  t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
1130  t1 = (-1.0 / 12.0) * w * (t + 4.0);
1131  yWeight[2] = TReal(t0 + t1);
1132  yWeight[3] = TReal(t0 - t1);
1133  t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
1134  t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
1135  yWeight[1] = TReal(t0 + t1);
1136  yWeight[4] = TReal(t0 - t1);
1137  /* z */
1138  w = z - static_cast<double>(zIndex[2]);
1139  w2 = w * w;
1140  zWeight[5] = TReal((1.0 / 120.0) * w * w2 * w2);
1141  w2 -= w;
1142  w4 = w2 * w2;
1143  w -= 1.0 / 2.0;
1144  t = w2 * (w2 - 3.0);
1145  zWeight[0] = TReal((1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - zWeight[5]);
1146  t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
1147  t1 = (-1.0 / 12.0) * w * (t + 4.0);
1148  zWeight[2] = TReal(t0 + t1);
1149  zWeight[3] = TReal(t0 - t1);
1150  t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
1151  t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
1152  zWeight[1] = TReal(t0 + t1);
1153  zWeight[4] = TReal(t0 - t1);
1154  /* t */
1155  w = t - static_cast<double>(tIndex[2]);
1156  w2 = w * w;
1157  tWeight[5] = TReal((1.0 / 120.0) * w * w2 * w2);
1158  w2 -= w;
1159  w4 = w2 * w2;
1160  w -= 1.0 / 2.0;
1161  t = w2 * (w2 - 3.0);
1162  tWeight[0] = TReal((1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - tWeight[5]);
1163  t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
1164  t1 = (-1.0 / 12.0) * w * (t + 4.0);
1165  tWeight[2] = TReal(t0 + t1);
1166  tWeight[3] = TReal(t0 - t1);
1167  t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
1168  t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
1169  tWeight[1] = TReal(t0 + t1);
1170  tWeight[4] = TReal(t0 - t1);
1171  break;
1172  default:
1173  cerr << "ComputeBSplineIndicesAndWeights: Unsupported B-spline degree: " << spline_degree << endl;
1174  exit(1);
1175  }
1176 }
1177 
1178 
1179 } // namespace mirtk
1180 
1181 #endif // MIRTK_BSpline_H
mirtk::BSpline::VariableToIndex
static int VariableToIndex(TReal)
Returns the lookup table index for a given value in [0,1].
Definition: BSpline.h:625
mirtk::BSpline::Weight
static MIRTKCU_API TReal Weight(TReal)
Returns the value of the B-spline function.
Definition: BSpline.h:291
mirtk::BSpline::B_I
static MIRTKCU_API TReal B_I(TReal)
Returns the 1st derivative value of the B-spline function.
Definition: BSpline.h:356
mirtk::BSpline::B
static MIRTKCU_API TReal B(TReal)
Returns the value of the B-spline function.
Definition: BSpline.h:308
mirtk::BSpline::B3
static MIRTKCU_API TReal B3(TReal)
Returns the value of the fourth B-spline basis function.
Definition: BSpline.h:349
mirtk::BSpline::LookupTable_II
static MIRTK_Numerics_EXPORT TReal LookupTable_II[LookupTableSize][4]
Lookup table of B-spline basis function 2nd derivative values.
Definition: BSpline.h:133
mirtk::BSpline::B2
static MIRTKCU_API TReal B2(TReal)
Returns the value of the third B-spline basis function.
Definition: BSpline.h:342
mirtk::BSpline::B2_II
static MIRTKCU_API TReal B2_II(TReal)
Returns the 2nd derivative value of the third B-spline basis function.
Definition: BSpline.h:458
mirtk::BSpline::B1_III
static MIRTKCU_API TReal B1_III(TReal)
Returns the 3rd derivative value of the second B-spline basis function.
Definition: BSpline.h:509
mirtk::BSpline::WeightLookupTable
static MIRTK_Numerics_EXPORT TReal WeightLookupTable[LookupTableSize]
Lookup table of B-spline function values.
Definition: BSpline.h:124
mirtk::BSpline::RealType
TReal RealType
Floating point precision type used.
Definition: BSpline.h:79
mirtk::BSpline::B2_III
static MIRTKCU_API TReal B2_III(TReal)
Returns the 3rd derivative value of the third B-spline basis function.
Definition: BSpline.h:516
mirtk::BSpline::B0_nI
static MIRTKCU_API TReal B0_nI(int, TReal)
Returns the n-th derivative value of the first B-spline basis function.
Definition: BSpline.h:556
mirtk::BSpline::B1
static MIRTKCU_API TReal B1(TReal)
Returns the value of the second B-spline basis function.
Definition: BSpline.h:335
mirtk::ifloor
MIRTKCU_API int ifloor(T x)
Round floating-point value to next smaller integer and cast to int.
Definition: Math.h:154
mirtk::BSpline::B3_nI
static MIRTKCU_API TReal B3_nI(int, TReal)
Returns the n-th derivative value of the fourth B-spline basis function.
Definition: BSpline.h:595
mirtk::BSpline::B0_I
static MIRTKCU_API TReal B0_I(TReal)
Returns the 1st derivative value of the first B-spline basis function.
Definition: BSpline.h:387
mirtk::BSpline::B_II
static MIRTKCU_API TReal B_II(TReal)
Returns the 2nd derivative value of the B-spline function.
Definition: BSpline.h:415
mirtk::BSpline::LookupTableSize
static const unsigned int LookupTableSize
Size of lookup tables of pre-computed B-spline values.
Definition: BSpline.h:104
mirtk::BSpline::LookupTable
static MIRTK_Numerics_EXPORT TReal LookupTable[LookupTableSize][4]
Lookup table of B-spline basis function values.
Definition: BSpline.h:127
mirtk
Definition: IOConfig.h:41
mirtk::BSpline::B_III
static MIRTKCU_API TReal B_III(TReal)
Returns the 3rd derivative value of the B-spline function.
Definition: BSpline.h:472
mirtk::BSpline::B2_nI
static MIRTKCU_API TReal B2_nI(int, TReal)
Returns the n-th derivative value of the third B-spline basis function.
Definition: BSpline.h:582
mirtk::BSpline::LatticeWeights_I
static const TReal LatticeWeights_I[4]
Definition: BSpline.h:87
mirtk::BSpline::LatticeWeights_II
static const TReal LatticeWeights_II[4]
Definition: BSpline.h:91
mirtk::BSpline::B3_II
static MIRTKCU_API TReal B3_II(TReal)
Returns the 2nd derivative value of the fourth B-spline basis function.
Definition: BSpline.h:465
mirtk::ComputeBSplineIndicesAndWeights
void ComputeBSplineIndicesAndWeights(double x, double y, int spline_degree, int xIndex [6], int yIndex [6], TReal xWeight[6], TReal yWeight[6])
Compute indices and weights required for 2D B-spline interpolation.
Definition: BSpline.h:673
mirtk::BSpline::B1_II
static MIRTKCU_API TReal B1_II(TReal)
Returns the 2nd derivative value of the second B-spline basis function.
Definition: BSpline.h:451
mirtk::BSpline::Weights
static MIRTKCU_API void Weights(TReal, TReal[4])
Definition: BSpline.h:278
mirtk::BSpline::B_nI
static MIRTKCU_API TReal B_nI(int, TReal)
Returns the n-th derivative value of the B-spline function.
Definition: BSpline.h:530
mirtk::BSpline::Initialize
static void Initialize(bool=false)
Initialize lookup tables.
Definition: BSpline.h:608
mirtk::BSpline::B1_I
static MIRTKCU_API TReal B1_I(TReal)
Returns the 1st derivative value of the second B-spline basis function.
Definition: BSpline.h:394
mirtk::BSpline::_initialized
static MIRTK_Numerics_EXPORT bool _initialized
Flag which indicates whether the lookup tables are initialized.
Definition: BSpline.h:231
mirtk::BSpline::B0_II
static MIRTKCU_API TReal B0_II(TReal)
Returns the 2nd derivative value of the first B-spline basis function.
Definition: BSpline.h:444
mirtk::iround
MIRTKCU_API int iround(T x)
Round floating-point value and cast to int.
Definition: Math.h:170
mirtk::BSpline::B2_I
static MIRTKCU_API TReal B2_I(TReal)
Returns the 1st derivative value of the third B-spline basis function.
Definition: BSpline.h:401
mirtk::BSpline::B0_III
static MIRTKCU_API TReal B0_III(TReal)
Returns the 3rd derivative value of the first B-spline basis function.
Definition: BSpline.h:502
mirtk::BSpline::B0
static MIRTKCU_API TReal B0(TReal)
Returns the value of the first B-spline basis function.
Definition: BSpline.h:328
mirtk::BSpline::B1_nI
static MIRTKCU_API TReal B1_nI(int, TReal)
Returns the n-th derivative value of the second B-spline basis function.
Definition: BSpline.h:569
mirtk::BSpline::B3_I
static MIRTKCU_API TReal B3_I(TReal)
Returns the 1st derivative value of the fourth B-spline basis function.
Definition: BSpline.h:408
mirtk::BSpline::LatticeWeights
static const TReal LatticeWeights[4]
Definition: BSpline.h:83
mirtk::BSpline::B3_III
static MIRTKCU_API TReal B3_III(TReal)
Returns the 3rd derivative value of the fourth B-spline basis function.
Definition: BSpline.h:523
mirtk::BSpline::LookupTable_I
static MIRTK_Numerics_EXPORT TReal LookupTable_I[LookupTableSize][4]
Lookup table of B-spline basis function 1st derivative values.
Definition: BSpline.h:130
© 2013-2018 Imperial College London, Andreas Schuh.