FastCubicBSplineInterpolateImageFunction.hxx
1 /*
2  * Medical Image Registration ToolKit (MIRTK)
3  *
4  * Copyright 2013-2015 Imperial College London
5  * Copyright 2013-2015 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_FastCubicBSplineInterpolateImageFunction_HXX
21 #define MIRTK_FastCubicBSplineInterpolateImageFunction_HXX
22 
23 #include "mirtk/FastCubicBSplineInterpolateImageFunction.h"
24 
25 #include "mirtk/Math.h"
26 #include "mirtk/BSpline.h"
27 #include "mirtk/Matrix.h"
28 #include "mirtk/Vector3D.h" // get
29 #include "mirtk/Vector4D.h" // get
30 
31 #include "mirtk/InterpolateImageFunction.hxx"
32 #include "mirtk/ImageToInterpolationCoefficients.h"
33 
34 
35 namespace mirtk {
36 
37 
38 // =============================================================================
39 // Construction/Destruction
40 // =============================================================================
41 
42 // -----------------------------------------------------------------------------
43 template <class TImage>
44 GenericFastCubicBSplineInterpolateImageFunction<TImage>
45 ::GenericFastCubicBSplineInterpolateImageFunction()
46 :
47  _InfiniteCoefficient(NULL)
48 {
49  // Default extrapolation mode is to apply the mirror boundary condition
50  // which is also assumed when converting an input image to spline coefficients
51  this->Extrapolator(ExtrapolatorType::New(Extrapolation_Mirror), true);
52 }
53 
54 // -----------------------------------------------------------------------------
55 template <class TImage>
56 GenericFastCubicBSplineInterpolateImageFunction<TImage>
57 ::~GenericFastCubicBSplineInterpolateImageFunction()
58 {
59  delete _InfiniteCoefficient;
60 }
61 
62 // -----------------------------------------------------------------------------
63 template <class TImage>
64 void GenericFastCubicBSplineInterpolateImageFunction<TImage>
65 ::Initialize(bool coeff)
66 {
67  // Initialize base class
68  Superclass::Initialize(coeff);
69 
70  // Initialize kernel lookup table
71  Kernel::Initialize();
72 
73  // Domain on which cubic B-spline interpolation is defined
74  const double margin = 2.0;
75  switch (this->NumberOfDimensions()) {
76  case 4:
77  this->_t1 = fdec(margin);
78  this->_t2 = this->Input()->T() - margin - 1;
79  case 3:
80  this->_z1 = fdec(margin);
81  this->_z2 = this->Input()->Z() - margin - 1;
82  default:
83  this->_y1 = fdec(margin);
84  this->_y2 = this->Input()->Y() - margin - 1;
85  this->_x1 = fdec(margin);
86  this->_x2 = this->Input()->X() - margin - 1;
87  }
88 
89  // Initialize coefficient image
90  RealType *data = nullptr;
91  _UseInputCoefficients = coeff;
92  if (_UseInputCoefficients && this->Input()->GetDataType() == voxel_info<RealType>::type()) {
93  data = reinterpret_cast<RealType *>(const_cast<void *>(this->Input()->GetDataPointer()));
94  }
95  _Coefficient.Initialize(this->Input()->Attributes(), data);
96  this->Update();
97 
98  // Initialize infinite coefficient image (i.e., extrapolator)
99  if (!_InfiniteCoefficient || _InfiniteCoefficient->ExtrapolationMode() != this->ExtrapolationMode()) {
100  delete _InfiniteCoefficient;
101  _InfiniteCoefficient = CoefficientExtrapolator::New(this->ExtrapolationMode(), &_Coefficient);
102  }
103  if (_InfiniteCoefficient) {
104  _InfiniteCoefficient->Input(&_Coefficient);
105  _InfiniteCoefficient->Initialize();
106  }
107 
108  // Compute strides for fast iteration over coefficient image
109  //_s1 = 1;
110  _s2 = this->Input()->X() - 4;
111  _s3 = (this->Input()->Y() - 4) * this->Input()->X();
112  _s4 = (this->Input()->Z() - 4) * this->Input()->Y() * this->Input()->X();
113 }
114 
115 // -----------------------------------------------------------------------------
116 template <class TImage>
117 void GenericFastCubicBSplineInterpolateImageFunction<TImage>::Update()
118 {
119  if (_Coefficient.GetDataPointer() != this->Input()->GetDataPointer()) {
120  _Coefficient = *(this->Input());
121  if (!_UseInputCoefficients) {
122  Real poles[2];
123  int npoles;
124  SplinePoles(3, poles, npoles);
125  FillBackgroundBeforeConversionToSplineCoefficients(_Coefficient);
126  switch (this->NumberOfDimensions()) {
127  case 4: ConvertToInterpolationCoefficientsT(_Coefficient, poles, npoles);
128  case 3: ConvertToInterpolationCoefficientsZ(_Coefficient, poles, npoles);
129  default: ConvertToInterpolationCoefficientsY(_Coefficient, poles, npoles);
130  ConvertToInterpolationCoefficientsX(_Coefficient, poles, npoles);
131  }
132  }
133  }
134 }
135 
136 // =============================================================================
137 // Domain checks
138 // =============================================================================
139 
140 // -----------------------------------------------------------------------------
141 template <class TImage>
142 void GenericFastCubicBSplineInterpolateImageFunction<TImage>
143 ::BoundingInterval(double x, int &i, int &I) const
144 {
145  i = ifloor(x), I = i + 3;
146 }
147 
148 // =============================================================================
149 // Evaluation
150 // =============================================================================
151 
152 // -----------------------------------------------------------------------------
153 template <class TImage>
154 inline typename GenericFastCubicBSplineInterpolateImageFunction<TImage>::VoxelType
155 GenericFastCubicBSplineInterpolateImageFunction<TImage>
156 ::Get2D(double x, double y, double z, double t) const
157 {
158  int i = ifloor(x);
159  int j = ifloor(y);
160  int k = iround(z);
161  int l = iround(t);
162 
163  if (k < 0 || k >= _Coefficient.Z() ||
164  l < 0 || l >= _Coefficient.T()) {
165  return voxel_cast<VoxelType>(this->DefaultValue());
166  }
167 
168  const int A = Kernel::VariableToIndex(Real(x - i));
169  const int B = Kernel::VariableToIndex(Real(y - j));
170 
171  --i, --j;
172 
173  RealType val = voxel_cast<RealType>(0);
174  Real w;
175 
176  int ia, jb;
177  for (int b = 0; b <= 3; ++b) {
178  jb = j + b;
179  DefaultExtrapolator::Apply(jb, _Coefficient.Y() - 1);
180  for (int a = 0; a <= 3; ++a) {
181  ia = i + a;
182  DefaultExtrapolator::Apply(ia, _Coefficient.X() - 1);
183  w = Kernel::LookupTable[A][a] * Kernel::LookupTable[B][b];
184  val += w * _Coefficient(ia, jb, k, l);
185  }
186  }
187 
188  return voxel_cast<VoxelType>(val);
189 }
190 
191 // -----------------------------------------------------------------------------
192 template <class TImage>
193 inline typename GenericFastCubicBSplineInterpolateImageFunction<TImage>::VoxelType
194 GenericFastCubicBSplineInterpolateImageFunction<TImage>
195 ::GetWithPadding2D(double x, double y, double z, double t) const
196 {
197  int i = ifloor(x);
198  int j = ifloor(y);
199  int k = iround(z);
200  int l = iround(t);
201 
202  if (k < 0 || k >= _Coefficient.Z() ||
203  l < 0 || l >= _Coefficient.T()) {
204  return voxel_cast<VoxelType>(this->DefaultValue());
205  }
206 
207  const int A = Kernel::VariableToIndex(Real(x - i));
208  const int B = Kernel::VariableToIndex(Real(y - j));
209 
210  --i, --j;
211 
212  RealType val = voxel_cast<RealType>(0);
213  Real fgw(0), bgw(0), w;
214 
215  int ia, jb;
216  for (int b = 0; b <= 3; ++b) {
217  jb = j + b;
218  DefaultExtrapolator::Apply(jb, _Coefficient.Y() - 1);
219  for (int a = 0; a <= 3; ++a) {
220  ia = i + a;
221  DefaultExtrapolator::Apply(ia, _Coefficient.X() - 1);
222  w = Kernel::LookupTable[A][a] * Kernel::LookupTable[B][b];
223  val += w * _Coefficient(ia, jb, k, l);
224  if (this->Input()->IsInsideForeground(i + a, j + b, k, l)) {
225  fgw += w;
226  } else {
227  bgw += w;
228  }
229  }
230  }
231 
232  if (bgw > fgw || AreEqual(bgw, fgw, 1e-3)) {
233  val = voxel_cast<RealType>(this->DefaultValue());
234  }
235  return voxel_cast<VoxelType>(val);
236 }
237 
238 // -----------------------------------------------------------------------------
239 template <class TImage> template <class TCoefficient>
240 inline typename TCoefficient::VoxelType
241 GenericFastCubicBSplineInterpolateImageFunction<TImage>
242 ::Get2D(const TCoefficient *coeff, double x, double y, double z, double t) const
243 {
244  typedef typename TCoefficient::VoxelType VoxelType;
245  typedef typename TCoefficient::RealType RealType;
246 
247  int i = ifloor(x);
248  int j = ifloor(y);
249  int k = iround(z);
250  int l = iround(t);
251 
252  const int A = Kernel::VariableToIndex(Real(x - i));
253  const int B = Kernel::VariableToIndex(Real(y - j));
254 
255  --i, --j;
256 
257  RealType val = voxel_cast<RealType>(0);
258 
259  int jb;
260  for (int b = 0; b <= 3; ++b) {
261  jb = j + b;
262  val += Kernel::LookupTable[A][0] * Kernel::LookupTable[B][b] * coeff->Get(i, jb, k, l);
263  val += Kernel::LookupTable[A][1] * Kernel::LookupTable[B][b] * coeff->Get(i+1, jb, k, l);
264  val += Kernel::LookupTable[A][2] * Kernel::LookupTable[B][b] * coeff->Get(i+2, jb, k, l);
265  val += Kernel::LookupTable[A][3] * Kernel::LookupTable[B][b] * coeff->Get(i+3, jb, k, l);
266  }
267 
268  return voxel_cast<VoxelType>(val);
269 }
270 
271 // -----------------------------------------------------------------------------
272 template <class TImage> template <class TOtherImage, class TCoefficient>
273 inline typename TCoefficient::VoxelType
274 GenericFastCubicBSplineInterpolateImageFunction<TImage>
275 ::GetWithPadding2D(const TOtherImage *input, const TCoefficient *coeff,
276  double x, double y, double z, double t) const
277 {
278  typedef typename TCoefficient::VoxelType VoxelType;
279  typedef typename TCoefficient::RealType RealType;
280 
281  int i = ifloor(x);
282  int j = ifloor(y);
283  int k = iround(z);
284  int l = iround(t);
285 
286  const int A = Kernel::VariableToIndex(Real(x - i));
287  const int B = Kernel::VariableToIndex(Real(y - j));
288 
289  --i, --j;
290 
291  RealType val = voxel_cast<RealType>(0);
292  Real fgw(0), bgw(0), w;
293 
294  int ia, jb;
295  for (int b = 0; b <= 3; ++b) {
296  jb = j + b;
297  for (int a = 0; a <= 3; ++a) {
298  ia = i + a;
299  w = Kernel::LookupTable[A][a] * Kernel::LookupTable[B][b];
300  val += w * voxel_cast<RealType>(coeff->Get(ia, jb, k, l));
301  if (input->IsForeground(ia, jb, k, l)) {
302  fgw += w;
303  } else {
304  bgw += w;
305  }
306  }
307  }
308 
309  if (bgw > fgw || AreEqual(bgw, fgw, 1e-3)) {
310  val = voxel_cast<RealType>(this->DefaultValue());
311  }
312  return voxel_cast<VoxelType>(val);
313 }
314 
315 // -----------------------------------------------------------------------------
316 template <class TImage>
317 inline typename GenericFastCubicBSplineInterpolateImageFunction<TImage>::VoxelType
318 GenericFastCubicBSplineInterpolateImageFunction<TImage>
319 ::Get3D(double x, double y, double z, double t) const
320 {
321  int i = ifloor(x);
322  int j = ifloor(y);
323  int k = ifloor(z);
324  int l = iround(t);
325 
326  if (l < 0 || l >= _Coefficient.T()) {
327  return voxel_cast<VoxelType>(this->DefaultValue());
328  }
329 
330  const int A = Kernel::VariableToIndex(Real(x - i));
331  const int B = Kernel::VariableToIndex(Real(y - j));
332  const int C = Kernel::VariableToIndex(Real(z - k));
333 
334  --i, --j, --k;
335 
336  RealType val = voxel_cast<RealType>(0);
337  Real wyz, w;
338 
339  int ia, jb, kc;
340  for (int c = 0; c <= 3; ++c) {
341  kc = k + c;
342  DefaultExtrapolator::Apply(kc, _Coefficient.Z() - 1);
343  for (int b = 0; b <= 3; ++b) {
344  jb = j + b;
345  DefaultExtrapolator::Apply(jb, _Coefficient.Y() - 1);
346  wyz = Kernel::LookupTable[B][b] * Kernel::LookupTable[C][c];
347  for (int a = 0; a <= 3; ++a) {
348  ia = i + a;
349  DefaultExtrapolator::Apply(ia, _Coefficient.X() - 1);
350  w = Kernel::LookupTable[A][a] * wyz;
351  val += w * _Coefficient(ia, jb, kc, l);
352  }
353  }
354  }
355 
356  return voxel_cast<VoxelType>(val);
357 }
358 
359 // -----------------------------------------------------------------------------
360 template <class TImage>
361 inline typename GenericFastCubicBSplineInterpolateImageFunction<TImage>::VoxelType
362 GenericFastCubicBSplineInterpolateImageFunction<TImage>
363 ::GetWithPadding3D(double x, double y, double z, double t) const
364 {
365  int i = ifloor(x);
366  int j = ifloor(y);
367  int k = ifloor(z);
368  int l = iround(t);
369 
370  if (l < 0 || l >= _Coefficient.T()) {
371  return voxel_cast<VoxelType>(this->DefaultValue());
372  }
373 
374  const int A = Kernel::VariableToIndex(Real(x - i));
375  const int B = Kernel::VariableToIndex(Real(y - j));
376  const int C = Kernel::VariableToIndex(Real(z - k));
377 
378  --i, --j, --k;
379 
380  RealType val = voxel_cast<RealType>(0);
381  Real fgw(0), bgw(0), wyz, w;
382 
383  int ia, jb, kc;
384  for (int c = 0; c <= 3; ++c) {
385  kc = k + c;
386  DefaultExtrapolator::Apply(kc, _Coefficient.Z() - 1);
387  for (int b = 0; b <= 3; ++b) {
388  jb = j + b;
389  DefaultExtrapolator::Apply(jb, _Coefficient.Y() - 1);
390  wyz = Kernel::LookupTable[B][b] * Kernel::LookupTable[C][c];
391  for (int a = 0; a <= 3; ++a) {
392  ia = i + a;
393  DefaultExtrapolator::Apply(ia, _Coefficient.X() - 1);
394  w = Kernel::LookupTable[A][a] * wyz;
395  val += w * _Coefficient(ia, jb, kc, l);
396  if (this->Input()->IsInsideForeground(i + a, j + b, k + c, l)) {
397  fgw += w;
398  } else {
399  bgw += w;
400  }
401  }
402  }
403  }
404 
405  if (bgw > fgw || AreEqual(bgw, fgw, 1e-3)) {
406  val = voxel_cast<RealType>(this->DefaultValue());
407  }
408  return voxel_cast<VoxelType>(val);
409 }
410 
411 // -----------------------------------------------------------------------------
412 template <class TImage> template <class TCoefficient>
413 inline typename TCoefficient::VoxelType
414 GenericFastCubicBSplineInterpolateImageFunction<TImage>
415 ::Get3D(const TCoefficient *coeff, double x, double y, double z, double t) const
416 {
417  typedef typename TCoefficient::VoxelType VoxelType;
418  typedef typename TCoefficient::RealType RealType;
419 
420  int i = ifloor(x);
421  int j = ifloor(y);
422  int k = ifloor(z);
423  int l = iround(t);
424 
425  const int A = Kernel::VariableToIndex(Real(x - i));
426  const int B = Kernel::VariableToIndex(Real(y - j));
427  const int C = Kernel::VariableToIndex(Real(z - k));
428 
429  --i, --j, --k;
430 
431  RealType val = voxel_cast<RealType>(0);
432 
433  int jb, kc;
434  Real wyz;
435  for (int c = 0; c <= 3; ++c) {
436  kc = k + c;
437  for (int b = 0; b <= 3; ++b) {
438  jb = j + b;
439  wyz = Kernel::LookupTable[B][b] * Kernel::LookupTable[C][c];
440  val += Kernel::LookupTable[A][0] * wyz * voxel_cast<RealType>(coeff->Get(i, jb, kc, l));
441  val += Kernel::LookupTable[A][1] * wyz * voxel_cast<RealType>(coeff->Get(i+1, jb, kc, l));
442  val += Kernel::LookupTable[A][2] * wyz * voxel_cast<RealType>(coeff->Get(i+2, jb, kc, l));
443  val += Kernel::LookupTable[A][3] * wyz * voxel_cast<RealType>(coeff->Get(i+3, jb, kc, l));
444  }
445  }
446 
447  return voxel_cast<VoxelType>(val);
448 }
449 
450 // -----------------------------------------------------------------------------
451 template <class TImage> template <class TOtherImage, class TCoefficient>
452 inline typename TCoefficient::VoxelType
453 GenericFastCubicBSplineInterpolateImageFunction<TImage>
454 ::GetWithPadding3D(const TOtherImage *input, const TCoefficient *coeff,
455  double x, double y, double z, double t) const
456 {
457  typedef typename TCoefficient::VoxelType VoxelType;
458  typedef typename TCoefficient::RealType RealType;
459 
460  int i = ifloor(x);
461  int j = ifloor(y);
462  int k = ifloor(z);
463  int l = iround(t);
464 
465  const int A = Kernel::VariableToIndex(Real(x - i));
466  const int B = Kernel::VariableToIndex(Real(y - j));
467  const int C = Kernel::VariableToIndex(Real(z - k));
468 
469  --i, --j, --k;
470 
471  RealType val = voxel_cast<RealType>(0);
472  Real fgw(0), bgw(0), wyz, w;
473 
474  int ia, jb, kc;
475  for (int c = 0; c <= 3; ++c) {
476  kc = k + c;
477  DefaultExtrapolator::Apply(kc, _Coefficient.Z() - 1);
478  for (int b = 0; b <= 3; ++b) {
479  jb = j + b;
480  DefaultExtrapolator::Apply(jb, _Coefficient.Y() - 1);
481  wyz = Kernel::LookupTable[B][b] * Kernel::LookupTable[C][c];
482  for (int a = 0; a <= 3; ++a) {
483  ia = i + a;
484  DefaultExtrapolator::Apply(ia, _Coefficient.X() - 1);
485  w = Kernel::LookupTable[A][a] * wyz;
486  val += w * voxel_cast<RealType>(coeff->Get(ia, jb, kc, l));
487  if (input->IsForeground(i + a, j + b, k + c, l)) {
488  fgw += w;
489  } else {
490  bgw += w;
491  }
492  }
493  }
494  }
495 
496  if (bgw > fgw || AreEqual(bgw, fgw, 1e-3)) {
497  val = voxel_cast<RealType>(this->DefaultValue());
498  }
499  return voxel_cast<VoxelType>(val);
500 }
501 
502 // -----------------------------------------------------------------------------
503 template <class TImage>
504 inline typename GenericFastCubicBSplineInterpolateImageFunction<TImage>::VoxelType
505 GenericFastCubicBSplineInterpolateImageFunction<TImage>
506 ::Get4D(double x, double y, double z, double t) const
507 {
508  int i = ifloor(x);
509  int j = ifloor(y);
510  int k = ifloor(z);
511  int l = ifloor(t);
512 
513  const int A = Kernel::VariableToIndex(Real(x - i));
514  const int B = Kernel::VariableToIndex(Real(y - j));
515  const int C = Kernel::VariableToIndex(Real(z - k));
516  const int D = Kernel::VariableToIndex(Real(t - l));
517 
518  --i, --j, --k, --l;
519 
520  RealType val = voxel_cast<RealType>(0);
521  Real wzt, wyzt, w;
522 
523  int ia, jb, kc, ld;
524  for (int d = 0; d <= 3; ++d) {
525  ld = l + d;
526  DefaultExtrapolator::Apply(ld, _Coefficient.T() - 1);
527  for (int c = 0; c <= 3; ++c) {
528  kc = k + c;
529  DefaultExtrapolator::Apply(kc, _Coefficient.Z() - 1);
530  wzt = Kernel::LookupTable[C][c] * Kernel::LookupTable[D][d];
531  for (int b = 0; b <= 3; ++b) {
532  jb = j + b;
533  DefaultExtrapolator::Apply(jb, _Coefficient.Y() - 1);
534  wyzt = Kernel::LookupTable[B][b] * wzt;
535  for (int a = 0; a <= 3; ++a) {
536  ia = i + a;
537  DefaultExtrapolator::Apply(ia, _Coefficient.X() - 1);
538  w = Kernel::LookupTable[A][a] * wyzt;
539  val += w * _Coefficient(ia, jb, kc, ld);
540  }
541  }
542  }
543  }
544 
545  return voxel_cast<VoxelType>(val);
546 }
547 
548 // -----------------------------------------------------------------------------
549 template <class TImage>
550 inline typename GenericFastCubicBSplineInterpolateImageFunction<TImage>::VoxelType
551 GenericFastCubicBSplineInterpolateImageFunction<TImage>
552 ::GetWithPadding4D(double x, double y, double z, double t) const
553 {
554  int i = ifloor(x);
555  int j = ifloor(y);
556  int k = ifloor(z);
557  int l = ifloor(t);
558 
559  const int A = Kernel::VariableToIndex(Real(x - i));
560  const int B = Kernel::VariableToIndex(Real(y - j));
561  const int C = Kernel::VariableToIndex(Real(z - k));
562  const int D = Kernel::VariableToIndex(Real(t - l));
563 
564  --i, --j, --k, --l;
565 
566  RealType val = voxel_cast<RealType>(0);
567  Real fgw(0), bgw(0), wzt, wyzt, w;
568 
569  int ia, jb, kc, ld;
570  for (int d = 0; d <= 3; ++d) {
571  ld = l + d;
572  DefaultExtrapolator::Apply(ld, _Coefficient.T() - 1);
573  for (int c = 0; c <= 3; ++c) {
574  kc = k + c;
575  DefaultExtrapolator::Apply(kc, _Coefficient.Z() - 1);
576  wzt = Kernel::LookupTable[C][c] * Kernel::LookupTable[D][d];
577  for (int b = 0; b <= 3; ++b) {
578  jb = j + b;
579  DefaultExtrapolator::Apply(jb, _Coefficient.Y() - 1);
580  wyzt = Kernel::LookupTable[B][b] * wzt;
581  for (int a = 0; a <= 3; ++a) {
582  ia = i + a;
583  DefaultExtrapolator::Apply(ia, _Coefficient.X() - 1);
584  w = Kernel::LookupTable[A][a] * wyzt;
585  val += w * _Coefficient(ia, jb, kc, ld);
586  if (this->Input()->IsInsideForeground(i + a, j + b, k + c, l + d)) {
587  fgw += w;
588  } else {
589  bgw += w;
590  }
591  }
592  }
593  }
594  }
595 
596  if (bgw > fgw || AreEqual(bgw, fgw, 1e-3)) {
597  val = voxel_cast<RealType>(this->DefaultValue());
598  }
599  return voxel_cast<VoxelType>(val);
600 }
601 
602 // -----------------------------------------------------------------------------
603 template <class TImage> template <class TCoefficient>
604 inline typename TCoefficient::VoxelType
605 GenericFastCubicBSplineInterpolateImageFunction<TImage>
606 ::Get4D(const TCoefficient *coeff, double x, double y, double z, double t) const
607 {
608  typedef typename TCoefficient::VoxelType VoxelType;
609  typedef typename TCoefficient::RealType RealType;
610 
611  int i = ifloor(x);
612  int j = ifloor(y);
613  int k = ifloor(z);
614  int l = ifloor(t);
615 
616  const int A = Kernel::VariableToIndex(Real(x - i));
617  const int B = Kernel::VariableToIndex(Real(y - j));
618  const int C = Kernel::VariableToIndex(Real(z - k));
619  const int D = Kernel::VariableToIndex(Real(t - l));
620 
621  --i, --j, --k, --l;
622 
623  RealType val = voxel_cast<RealType>(0);
624 
625  int jb, kc, ld;
626  Real wzt, wyzt;
627  for (int d = 0; d <= 3; ++d) {
628  ld = l + d;
629  for (int c = 0; c <= 3; ++c) {
630  kc = k + c;
631  wzt = Kernel::LookupTable[C][c] * Kernel::LookupTable[D][d];
632  for (int b = 0; b <= 3; ++b) {
633  jb = j + b;
634  wyzt = Kernel::LookupTable[B][b] * wzt;
635  val += Kernel::LookupTable[A][0] * wyzt * voxel_cast<RealType>(coeff->Get(i, jb, kc, ld));
636  val += Kernel::LookupTable[A][1] * wyzt * voxel_cast<RealType>(coeff->Get(i+1, jb, kc, ld));
637  val += Kernel::LookupTable[A][2] * wyzt * voxel_cast<RealType>(coeff->Get(i+2, jb, kc, ld));
638  val += Kernel::LookupTable[A][3] * wyzt * voxel_cast<RealType>(coeff->Get(i+3, jb, kc, ld));
639  }
640  }
641  }
642 
643  return voxel_cast<VoxelType>(val);
644 }
645 
646 // -----------------------------------------------------------------------------
647 template <class TImage> template <class TOtherImage, class TCoefficient>
648 inline typename TCoefficient::VoxelType
649 GenericFastCubicBSplineInterpolateImageFunction<TImage>
650 ::GetWithPadding4D(const TOtherImage *input, const TCoefficient *coeff,
651  double x, double y, double z, double t) const
652 {
653  typedef typename TCoefficient::VoxelType VoxelType;
654  typedef typename TCoefficient::RealType RealType;
655 
656  int i = ifloor(x);
657  int j = ifloor(y);
658  int k = ifloor(z);
659  int l = ifloor(t);
660 
661  const int A = Kernel::VariableToIndex(Real(x - i));
662  const int B = Kernel::VariableToIndex(Real(y - j));
663  const int C = Kernel::VariableToIndex(Real(z - k));
664  const int D = Kernel::VariableToIndex(Real(t - l));
665 
666  --i, --j, --k, --l;
667 
668  RealType val = voxel_cast<RealType>(0);
669  Real fgw(0), bgw(0), wzt, wyzt, w;
670 
671  int ia, jb, kc, ld;
672  for (int d = 0; d <= 3; ++d) {
673  ld = l + d;
674  for (int c = 0; c <= 3; ++c) {
675  kc = k + c;
676  wzt = Kernel::LookupTable[C][c] * Kernel::LookupTable[D][d];
677  for (int b = 0; b <= 3; ++b) {
678  jb = j + b;
679  wyzt = Kernel::LookupTable[B][b] * wzt;
680  for (int a = 0; a <= 3; ++a) {
681  ia = i + a;
682  w = Kernel::LookupTable[A][a] * wyzt;
683  val += w * voxel_cast<RealType>(coeff->Get(ia, jb, kc, ld));
684  if (input->IsForeground(ia, jb, kc, ld)) {
685  fgw += w;
686  } else {
687  bgw += w;
688  }
689  }
690  }
691  }
692  }
693 
694  if (bgw > fgw || AreEqual(bgw, fgw, 1e-3)) {
695  val = voxel_cast<RealType>(this->DefaultValue());
696  }
697  return voxel_cast<VoxelType>(val);
698 }
699 
700 // -----------------------------------------------------------------------------
701 template <class TImage>
702 inline typename GenericFastCubicBSplineInterpolateImageFunction<TImage>::VoxelType
703 GenericFastCubicBSplineInterpolateImageFunction<TImage>
704 ::Get(double x, double y, double z, double t) const
705 {
706  switch (this->NumberOfDimensions()) {
707  case 3: return Get3D(x, y, z, t);
708  case 2: return Get2D(x, y, z, t);
709  default: return Get4D(x, y, z, t);
710  }
711 }
712 
713 // -----------------------------------------------------------------------------
714 template <class TImage>
715 inline typename GenericFastCubicBSplineInterpolateImageFunction<TImage>::VoxelType
716 GenericFastCubicBSplineInterpolateImageFunction<TImage>
717 ::GetWithPadding(double x, double y, double z, double t) const
718 {
719  switch (this->NumberOfDimensions()) {
720  case 3: return GetWithPadding3D(x, y, z, t);
721  case 2: return GetWithPadding2D(x, y, z, t);
722  default: return GetWithPadding4D(x, y, z, t);
723  }
724 }
725 
726 // -----------------------------------------------------------------------------
727 template <class TImage> template <class TOtherImage>
728 inline typename TOtherImage::VoxelType
729 GenericFastCubicBSplineInterpolateImageFunction<TImage>
730 ::Get(const TOtherImage *coeff, double x, double y, double z, double t) const
731 {
732  switch (this->NumberOfDimensions()) {
733  case 3: return Get3D(coeff, x, y, z, t);
734  case 2: return Get2D(coeff, x, y, z, t);
735  default: return Get4D(coeff, x, y, z, t);
736  }
737 }
738 
739 // -----------------------------------------------------------------------------
740 template <class TImage> template <class TOtherImage, class TCoefficient>
741 inline typename TCoefficient::VoxelType
742 GenericFastCubicBSplineInterpolateImageFunction<TImage>
743 ::GetWithPadding(const TOtherImage *input, const TCoefficient *coeff,
744  double x, double y, double z, double t) const
745 {
746  switch (this->NumberOfDimensions()) {
747  case 3: return GetWithPadding3D(input, coeff, x, y, z, t);
748  case 2: return GetWithPadding2D(input, coeff, x, y, z, t);
749  default: return GetWithPadding4D(input, coeff, x, y, z, t);
750  }
751 }
752 
753 // -----------------------------------------------------------------------------
754 template <class TImage>
755 inline typename GenericFastCubicBSplineInterpolateImageFunction<TImage>::VoxelType
756 GenericFastCubicBSplineInterpolateImageFunction<TImage>
757 ::GetInside2D(double x, double y, double z, double t) const
758 {
759  int i = ifloor(x);
760  int j = ifloor(y);
761  int k = iround(z);
762  int l = iround(t);
763 
764  const int A = Kernel::VariableToIndex(Real(x - i));
765  const int B = Kernel::VariableToIndex(Real(y - j));
766 
767  --i, --j;
768 
769  RealType val = voxel_cast<RealType>(0);
770  const RealType *coeff = _Coefficient.Data(i, j, k, l);
771 
772  for (int b = 0; b <= 3; ++b, coeff += _s2) {
773  val += Kernel::LookupTable[A][0] * Kernel::LookupTable[B][b] * (*coeff), ++coeff;
774  val += Kernel::LookupTable[A][1] * Kernel::LookupTable[B][b] * (*coeff), ++coeff;
775  val += Kernel::LookupTable[A][2] * Kernel::LookupTable[B][b] * (*coeff), ++coeff;
776  val += Kernel::LookupTable[A][3] * Kernel::LookupTable[B][b] * (*coeff), ++coeff;
777  }
778 
779  return voxel_cast<VoxelType>(val);
780 }
781 
782 // -----------------------------------------------------------------------------
783 template <class TImage>
784 inline typename GenericFastCubicBSplineInterpolateImageFunction<TImage>::VoxelType
785 GenericFastCubicBSplineInterpolateImageFunction<TImage>
786 ::GetInside3D(double x, double y, double z, double t) const
787 {
788  int i = ifloor(x);
789  int j = ifloor(y);
790  int k = ifloor(z);
791  int l = iround(t);
792 
793  const int A = Kernel::VariableToIndex(Real(x - i));
794  const int B = Kernel::VariableToIndex(Real(y - j));
795  const int C = Kernel::VariableToIndex(Real(z - k));
796 
797  --i, --j, --k;
798 
799  RealType val = voxel_cast<RealType>(0);
800  const RealType *coeff = _Coefficient.Data(i, j, k, l);
801 
802  Real wyz;
803  for (int c = 0; c <= 3; ++c, coeff += _s3) {
804  for (int b = 0; b <= 3; ++b, coeff += _s2) {
805  wyz = Kernel::LookupTable[B][b] * Kernel::LookupTable[C][c];
806  val += Kernel::LookupTable[A][0] * wyz * (*coeff), ++coeff;
807  val += Kernel::LookupTable[A][1] * wyz * (*coeff), ++coeff;
808  val += Kernel::LookupTable[A][2] * wyz * (*coeff), ++coeff;
809  val += Kernel::LookupTable[A][3] * wyz * (*coeff), ++coeff;
810  }
811  }
812 
813  return voxel_cast<VoxelType>(val);
814 }
815 
816 // -----------------------------------------------------------------------------
817 template <class TImage>
818 inline typename GenericFastCubicBSplineInterpolateImageFunction<TImage>::VoxelType
819 GenericFastCubicBSplineInterpolateImageFunction<TImage>
820 ::GetInside4D(double x, double y, double z, double t) const
821 {
822  int i = ifloor(x);
823  int j = ifloor(y);
824  int k = ifloor(z);
825  int l = ifloor(t);
826 
827  const int A = Kernel::VariableToIndex(Real(x - i));
828  const int B = Kernel::VariableToIndex(Real(y - j));
829  const int C = Kernel::VariableToIndex(Real(z - k));
830  const int D = Kernel::VariableToIndex(Real(t - l));
831 
832  --i, --j, --k, --l;
833 
834  RealType val = voxel_cast<RealType>(0);
835  const RealType *coeff = _Coefficient.Data(i, j, k, l);
836 
837  Real wzt, wyzt;
838  for (int d = 0; d <= 3; ++d, coeff += _s4) {
839  for (int c = 0; c <= 3; ++c, coeff += _s3) {
840  wzt = Kernel::LookupTable[C][c] * Kernel::LookupTable[D][d];
841  for (int b = 0; b <= 3; ++b, coeff += _s2) {
842  wyzt = Kernel::LookupTable[B][b] * wzt;
843  val += Kernel::LookupTable[A][0] * wyzt * (*coeff), ++coeff;
844  val += Kernel::LookupTable[A][1] * wyzt * (*coeff), ++coeff;
845  val += Kernel::LookupTable[A][2] * wyzt * (*coeff), ++coeff;
846  val += Kernel::LookupTable[A][3] * wyzt * (*coeff), ++coeff;
847  }
848  }
849  }
850 
851  return voxel_cast<VoxelType>(val);
852 }
853 
854 // -----------------------------------------------------------------------------
855 template <class TImage>
856 inline typename GenericFastCubicBSplineInterpolateImageFunction<TImage>::VoxelType
857 GenericFastCubicBSplineInterpolateImageFunction<TImage>
858 ::GetInside(double x, double y, double z, double t) const
859 {
860  // Use faster coefficient iteration than Get(Coefficient(), x, y, z, t)
861  switch (this->NumberOfDimensions()) {
862  case 3: return GetInside3D(x, y, z, t);
863  case 2: return GetInside2D(x, y, z, t);
864  default: return GetInside4D(x, y, z, t);
865  }
866 }
867 
868 // -----------------------------------------------------------------------------
869 template <class TImage>
870 inline typename GenericFastCubicBSplineInterpolateImageFunction<TImage>::VoxelType
871 GenericFastCubicBSplineInterpolateImageFunction<TImage>
872 ::GetOutside(double x, double y, double z, double t) const
873 {
874  if (_InfiniteCoefficient) {
875  return voxel_cast<VoxelType>(Get(_InfiniteCoefficient, x, y, z, t));
876  } else {
877  return Get(x, y, z, t);
878  }
879 }
880 
881 // -----------------------------------------------------------------------------
882 template <class TImage>
883 inline typename GenericFastCubicBSplineInterpolateImageFunction<TImage>::VoxelType
884 GenericFastCubicBSplineInterpolateImageFunction<TImage>
885 ::GetWithPaddingInside(double x, double y, double z, double t) const
886 {
887  switch (this->NumberOfDimensions()) {
888  case 3: return voxel_cast<VoxelType>(GetWithPadding3D(this->Input(), &_Coefficient, x, y, z, t));
889  case 2: return voxel_cast<VoxelType>(GetWithPadding2D(this->Input(), &_Coefficient, x, y, z, t));
890  default: return voxel_cast<VoxelType>(GetWithPadding4D(this->Input(), &_Coefficient, x, y, z, t));
891  }
892 }
893 
894 // -----------------------------------------------------------------------------
895 template <class TImage>
896 inline typename GenericFastCubicBSplineInterpolateImageFunction<TImage>::VoxelType
897 GenericFastCubicBSplineInterpolateImageFunction<TImage>
898 ::GetWithPaddingOutside(double x, double y, double z, double t) const
899 {
900  if (this->Extrapolator() && _InfiniteCoefficient) {
901  return voxel_cast<VoxelType>(GetWithPadding(this->Extrapolator(), _InfiniteCoefficient, x, y, z, t));
902  } else {
903  return GetWithPadding(x, y, z, t);
904  }
905 }
906 
907 // =============================================================================
908 // Evaluation of derivatives
909 // =============================================================================
910 
911 // -----------------------------------------------------------------------------
912 template <class TImage>
913 void GenericFastCubicBSplineInterpolateImageFunction<TImage>
914 ::Jacobian2D(Matrix &jac, double x, double y, double z, double t) const
915 {
916  jac.Initialize(_Coefficient.N(), 2);
917 
918  int i = ifloor(x);
919  int j = ifloor(y);
920  int k = iround(z);
921  int l = iround(t);
922 
923  if (k < 0 || k >= _Coefficient.Z() ||
924  l < 0 || l >= _Coefficient.T()) {
925  return;
926  }
927 
928  const int A = Kernel::VariableToIndex(Real(x - i));
929  const int B = Kernel::VariableToIndex(Real(y - j));
930 
931  --i, --j;
932 
933  RealType dx = voxel_cast<RealType>(0);
934  RealType dy = voxel_cast<RealType>(0);
935  Real wx[2], wy[2];
936 
937  int ia, jb;
938  for (int b = 0; b <= 3; ++b) {
939  jb = j + b;
940  DefaultExtrapolator::Apply(jb, _Coefficient.Y() - 1);
941  wy[0] = Kernel::LookupTable [B][b];
942  wy[1] = Kernel::LookupTable_I[B][b];
943  for (int a = 0; a <= 3; ++a) {
944  ia = i + a;
945  DefaultExtrapolator::Apply(ia, _Coefficient.X() - 1);
946  wx[0] = Kernel::LookupTable [A][a];
947  wx[1] = Kernel::LookupTable_I[A][a];
948  dx += (wx[1] * wy[0]) * _Coefficient(ia, jb, k, l);
949  dy += (wx[0] * wy[1]) * _Coefficient(ia, jb, k, l);
950  }
951  }
952 
953  for (int r = 0; r < _Coefficient.N(); ++r) {
954  jac(r, 0) = get(dx, r);
955  jac(r, 1) = get(dy, r);
956  }
957 }
958 
959 // -----------------------------------------------------------------------------
960 template <class TImage> template <class TCoefficient>
961 void GenericFastCubicBSplineInterpolateImageFunction<TImage>
962 ::Jacobian2D(Matrix &jac, const TCoefficient *coeff,
963  double x, double y, double z, double t) const
964 {
965  jac.Initialize(coeff->N(), 2);
966 
967  int i = ifloor(x);
968  int j = ifloor(y);
969  int k = iround(z);
970  int l = iround(t);
971 
972  const int A = Kernel::VariableToIndex(Real(x - i));
973  const int B = Kernel::VariableToIndex(Real(y - j));
974 
975  --i, --j;
976 
977  RealType dx = voxel_cast<RealType>(0);
978  RealType dy = voxel_cast<RealType>(0);
979  Real wx[2], wy[2];
980 
981  int ia, jb;
982  for (int b = 0; b <= 3; ++b) {
983  jb = j + b;
984  wy[0] = Kernel::LookupTable [B][b];
985  wy[1] = Kernel::LookupTable_I[B][b];
986  for (int a = 0; a <= 3; ++a) {
987  ia = i + a;
988  wx[0] = Kernel::LookupTable [A][a];
989  wx[1] = Kernel::LookupTable_I[A][a];
990  dx += (wx[1] * wy[0]) * voxel_cast<RealType>(coeff->Get(ia, jb, k, l));
991  dy += (wx[0] * wy[1]) * voxel_cast<RealType>(coeff->Get(ia, jb, k, l));
992  }
993  }
994 
995  for (int r = 0; r < coeff->N(); ++r) {
996  jac(r, 0) = get(dx, r);
997  jac(r, 1) = get(dy, r);
998  }
999 }
1000 
1001 // -----------------------------------------------------------------------------
1002 template <class TImage>
1003 void GenericFastCubicBSplineInterpolateImageFunction<TImage>
1004 ::Jacobian3D(Matrix &jac, double x, double y, double z, double t) const
1005 {
1006  jac.Initialize(_Coefficient.N(), 3);
1007 
1008  int i = ifloor(x);
1009  int j = ifloor(y);
1010  int k = ifloor(z);
1011  int l = iround(t);
1012 
1013  if (l < 0 || l >= _Coefficient.T()) {
1014  return;
1015  }
1016 
1017  const int A = Kernel::VariableToIndex(Real(x - i));
1018  const int B = Kernel::VariableToIndex(Real(y - j));
1019  const int C = Kernel::VariableToIndex(Real(z - k));
1020 
1021  --i, --j, --k;
1022 
1023  RealType dx = voxel_cast<RealType>(0);
1024  RealType dy = voxel_cast<RealType>(0);
1025  RealType dz = voxel_cast<RealType>(0);
1026  Real wx[2], wy[2], wz[2];
1027 
1028  int ia, jb, kc;
1029  for (int c = 0; c <= 3; ++c) {
1030  kc = k + c;
1031  DefaultExtrapolator::Apply(kc, _Coefficient.Z() - 1);
1032  wz[0] = Kernel::LookupTable [C][c];
1033  wz[1] = Kernel::LookupTable_I[C][c];
1034  for (int b = 0; b <= 3; ++b) {
1035  jb = j + b;
1036  DefaultExtrapolator::Apply(jb, _Coefficient.Y() - 1);
1037  wy[0] = Kernel::LookupTable [B][b];
1038  wy[1] = Kernel::LookupTable_I[B][b];
1039  for (int a = 0; a <= 3; ++a) {
1040  ia = i + a;
1041  DefaultExtrapolator::Apply(ia, _Coefficient.X() - 1);
1042  wx[0] = Kernel::LookupTable [A][a];
1043  wx[1] = Kernel::LookupTable_I[A][a];
1044  dx += (wx[1] * wy[0] * wz[0]) * _Coefficient(ia, jb, kc, l);
1045  dy += (wx[0] * wy[1] * wz[0]) * _Coefficient(ia, jb, kc, l);
1046  dz += (wx[0] * wy[0] * wz[1]) * _Coefficient(ia, jb, kc, l);
1047  }
1048  }
1049  }
1050 
1051  for (int r = 0; r < _Coefficient.N(); ++r) {
1052  jac(r, 0) = get(dx, r);
1053  jac(r, 1) = get(dy, r);
1054  jac(r, 2) = get(dz, r);
1055  }
1056 }
1057 
1058 // -----------------------------------------------------------------------------
1059 template <class TImage> template <class TOtherImage>
1060 void GenericFastCubicBSplineInterpolateImageFunction<TImage>
1061 ::Jacobian3D(Matrix &jac, const TOtherImage *coeff,
1062  double x, double y, double z, double t) const
1063 {
1064  jac.Initialize(coeff->N(), 3);
1065 
1066  int i = ifloor(x);
1067  int j = ifloor(y);
1068  int k = ifloor(z);
1069  int l = iround(t);
1070 
1071  const int A = Kernel::VariableToIndex(Real(x - i));
1072  const int B = Kernel::VariableToIndex(Real(y - j));
1073  const int C = Kernel::VariableToIndex(Real(z - k));
1074 
1075  --i, --j, --k;
1076 
1077  RealType dx = voxel_cast<RealType>(0);
1078  RealType dy = voxel_cast<RealType>(0);
1079  RealType dz = voxel_cast<RealType>(0);
1080  Real wx[2], wy[2], wz[2];
1081 
1082  int ia, jb, kc;
1083  for (int c = 0; c <= 3; ++c) {
1084  kc = k + c;
1085  wz[0] = Kernel::LookupTable [C][c];
1086  wz[1] = Kernel::LookupTable_I[C][c];
1087  for (int b = 0; b <= 3; ++b) {
1088  jb = j + b;
1089  wy[0] = Kernel::LookupTable [B][b];
1090  wy[1] = Kernel::LookupTable_I[B][b];
1091  for (int a = 0; a <= 3; ++a) {
1092  ia = i + a;
1093  wx[0] = Kernel::LookupTable [A][a];
1094  wx[1] = Kernel::LookupTable_I[A][a];
1095  dx += (wx[1] * wy[0] * wz[0]) * voxel_cast<RealType>(coeff->Get(ia, jb, kc, l));
1096  dy += (wx[0] * wy[1] * wz[0]) * voxel_cast<RealType>(coeff->Get(ia, jb, kc, l));
1097  dz += (wx[0] * wy[0] * wz[1]) * voxel_cast<RealType>(coeff->Get(ia, jb, kc, l));
1098  }
1099  }
1100  }
1101 
1102  for (int r = 0; r < coeff->N(); ++r) {
1103  jac(r, 0) = get(dx, r);
1104  jac(r, 1) = get(dy, r);
1105  jac(r, 2) = get(dz, r);
1106  }
1107 }
1108 
1109 // -----------------------------------------------------------------------------
1110 template <class TImage>
1111 void GenericFastCubicBSplineInterpolateImageFunction<TImage>
1112 ::Jacobian4D(Matrix &jac, double x, double y, double z, double t) const
1113 {
1114  jac.Initialize(_Coefficient.N(), 4);
1115 
1116  int i = ifloor(x);
1117  int j = ifloor(y);
1118  int k = ifloor(z);
1119  int l = ifloor(t);
1120 
1121  const int A = Kernel::VariableToIndex(Real(x - i));
1122  const int B = Kernel::VariableToIndex(Real(y - j));
1123  const int C = Kernel::VariableToIndex(Real(z - k));
1124  const int D = Kernel::VariableToIndex(Real(t - l));
1125 
1126  --i, --j, --k, --l;
1127 
1128  RealType dx = voxel_cast<RealType>(0);
1129  RealType dy = voxel_cast<RealType>(0);
1130  RealType dz = voxel_cast<RealType>(0);
1131  RealType dt = voxel_cast<RealType>(0);
1132  Real wx[2], wy[2], wz[2], wt[2];
1133 
1134  int ia, jb, kc, ld;
1135  for (int d = 0; d <= 3; ++d) {
1136  ld = l + d;
1137  DefaultExtrapolator::Apply(ld, _Coefficient.T() - 1);
1138  wt[0] = Kernel::LookupTable [D][d];
1139  wt[1] = Kernel::LookupTable_I[D][d];
1140  for (int c = 0; c <= 3; ++c) {
1141  kc = k + c;
1142  DefaultExtrapolator::Apply(kc, _Coefficient.Z() - 1);
1143  wz[0] = Kernel::LookupTable [C][c];
1144  wz[1] = Kernel::LookupTable_I[C][c];
1145  for (int b = 0; b <= 3; ++b) {
1146  jb = j + b;
1147  DefaultExtrapolator::Apply(jb, _Coefficient.Y() - 1);
1148  wy[0] = Kernel::LookupTable [B][b];
1149  wy[1] = Kernel::LookupTable_I[B][b];
1150  for (int a = 0; a <= 3; ++a) {
1151  ia = i + a;
1152  DefaultExtrapolator::Apply(ia, _Coefficient.X() - 1);
1153  wx[0] = Kernel::LookupTable [A][a];
1154  wx[1] = Kernel::LookupTable_I[A][a];
1155  dx += (wx[1] * wy[0] * wz[0] * wt[0]) * _Coefficient(ia, jb, kc, ld);
1156  dy += (wx[0] * wy[1] * wz[0] * wt[0]) * _Coefficient(ia, jb, kc, ld);
1157  dz += (wx[0] * wy[0] * wz[1] * wt[0]) * _Coefficient(ia, jb, kc, ld);
1158  dt += (wx[0] * wy[0] * wz[0] * wt[1]) * _Coefficient(ia, jb, kc, ld);
1159  }
1160  }
1161  }
1162  }
1163 
1164  for (int r = 0; r < _Coefficient.N(); ++r) {
1165  jac(r, 0) = get(dx, r);
1166  jac(r, 1) = get(dy, r);
1167  jac(r, 2) = get(dz, r);
1168  jac(r, 3) = get(dt, r);
1169  }
1170 }
1171 
1172 // -----------------------------------------------------------------------------
1173 template <class TImage> template <class TOtherImage>
1174 void GenericFastCubicBSplineInterpolateImageFunction<TImage>
1175 ::Jacobian4D(Matrix &jac, const TOtherImage *coeff,
1176  double x, double y, double z, double t) const
1177 {
1178  jac.Initialize(coeff->N(), 4);
1179 
1180  int i = ifloor(x);
1181  int j = ifloor(y);
1182  int k = ifloor(z);
1183  int l = ifloor(t);
1184 
1185  const int A = Kernel::VariableToIndex(Real(x - i));
1186  const int B = Kernel::VariableToIndex(Real(y - j));
1187  const int C = Kernel::VariableToIndex(Real(z - k));
1188  const int D = Kernel::VariableToIndex(Real(t - l));
1189 
1190  --i, --j, --k, --l;
1191 
1192  RealType dx = voxel_cast<RealType>(0);
1193  RealType dy = voxel_cast<RealType>(0);
1194  RealType dz = voxel_cast<RealType>(0);
1195  RealType dt = voxel_cast<RealType>(0);
1196  Real wx[2], wy[2], wz[2], wt[2];
1197 
1198  int ia, jb, kc, ld;
1199  for (int d = 0; d <= 3; ++d) {
1200  ld = l + d;
1201  wt[0] = Kernel::LookupTable [D][d];
1202  wt[1] = Kernel::LookupTable_I[D][d];
1203  for (int c = 0; c <= 3; ++c) {
1204  kc = k + c;
1205  wz[0] = Kernel::LookupTable [C][c];
1206  wz[1] = Kernel::LookupTable_I[C][c];
1207  for (int b = 0; b <= 3; ++b) {
1208  jb = j + b;
1209  wy[0] = Kernel::LookupTable [B][b];
1210  wy[1] = Kernel::LookupTable_I[B][b];
1211  for (int a = 0; a <= 3; ++a) {
1212  ia = i + a;
1213  wx[0] = Kernel::LookupTable [A][a];
1214  wx[1] = Kernel::LookupTable_I[A][a];
1215  dx += (wx[1] * wy[0] * wz[0] * wt[0]) * voxel_cast<RealType>(coeff->Get(ia, jb, kc, ld));
1216  dy += (wx[0] * wy[1] * wz[0] * wt[0]) * voxel_cast<RealType>(coeff->Get(ia, jb, kc, ld));
1217  dz += (wx[0] * wy[0] * wz[1] * wt[0]) * voxel_cast<RealType>(coeff->Get(ia, jb, kc, ld));
1218  dt += (wx[0] * wy[0] * wz[0] * wt[1]) * voxel_cast<RealType>(coeff->Get(ia, jb, kc, ld));
1219  }
1220  }
1221  }
1222  }
1223 
1224  for (int r = 0; r < _Coefficient.N(); ++r) {
1225  jac(r, 0) = get(dx, r);
1226  jac(r, 1) = get(dy, r);
1227  jac(r, 2) = get(dz, r);
1228  jac(r, 3) = get(dt, r);
1229  }
1230 }
1231 
1232 // -----------------------------------------------------------------------------
1233 template <class TImage>
1234 void GenericFastCubicBSplineInterpolateImageFunction<TImage>
1235 ::Jacobian(Matrix &jac, double x, double y, double z, double t) const
1236 {
1237  switch (this->NumberOfDimensions()) {
1238  case 3: return Jacobian3D(jac, x, y, z, t);
1239  case 2: return Jacobian2D(jac, x, y, z, t);
1240  default: return Jacobian4D(jac, x, y, z, t);
1241  }
1242 }
1243 
1244 // -----------------------------------------------------------------------------
1245 template <class TImage> template <class TOtherImage>
1246 void GenericFastCubicBSplineInterpolateImageFunction<TImage>
1247 ::Jacobian(Matrix &jac, const TOtherImage *coeff,
1248  double x, double y, double z, double t) const
1249 {
1250  switch (this->NumberOfDimensions()) {
1251  case 3: return Jacobian3D(jac, coeff, x, y, z, t);
1252  case 2: return Jacobian2D(jac, coeff, x, y, z, t);
1253  default: return Jacobian4D(jac, coeff, x, y, z, t);
1254  }
1255 }
1256 
1257 // -----------------------------------------------------------------------------
1258 template <class TImage>
1259 void GenericFastCubicBSplineInterpolateImageFunction<TImage>
1260 ::EvaluateJacobianInside(Matrix &jac, double x, double y, double z, double t) const
1261 {
1262  Jacobian(jac, &_Coefficient, x, y, z, t);
1263 }
1264 
1265 // -----------------------------------------------------------------------------
1266 template <class TImage>
1267 void GenericFastCubicBSplineInterpolateImageFunction<TImage>
1268 ::EvaluateJacobianOutside(Matrix &jac, double x, double y, double z, double t) const
1269 {
1270  Jacobian(jac, _InfiniteCoefficient, x, y, z, t);
1271 }
1272 
1273 
1274 } // namespace mirtk
1275 
1276 #endif // MIRTK_FastCubicBSplineInterpolateImageFunction_HXX
mirtk::GenericFastCubicBSplineInterpolateImageFunction::EvaluateJacobianOutside
virtual void EvaluateJacobianOutside(Matrix &, double, double, double=0, double=0) const
Get 1st order derivatives of given image at arbitrary location (in pixels)
Definition: FastCubicBSplineInterpolateImageFunction.hxx:1268
mirtk::GenericFastCubicBSplineInterpolateImageFunction::Jacobian4D
void Jacobian4D(Matrix &, double, double, double=0, double=0) const
Definition: FastCubicBSplineInterpolateImageFunction.hxx:1112
mirtk::GenericFastCubicBSplineInterpolateImageFunction::GetWithPadding2D
VoxelType GetWithPadding2D(double, double, double=0, double=0) const
Definition: FastCubicBSplineInterpolateImageFunction.hxx:195
mirtk::SplinePoles
void SplinePoles(int degree, TReal pole[2], int &npoles)
Recover spline poles from a lookup table.
Definition: ImageToInterpolationCoefficients.h:379
mirtk::GenericFastCubicBSplineInterpolateImageFunction::Jacobian3D
void Jacobian3D(Matrix &, double, double, double=0, double=0) const
Definition: FastCubicBSplineInterpolateImageFunction.hxx:1004
mirtk::GenericFastCubicBSplineInterpolateImageFunction::GetInside4D
virtual VoxelType GetInside4D(double, double, double=0, double=0) const
Evaluate generic 4D image without handling boundary conditions.
Definition: FastCubicBSplineInterpolateImageFunction.hxx:820
mirtk::Get
string Get(const ParameterList &params, string name)
Get parameter value from parameters list.
Definition: Object.h:202
mirtk::GenericFastCubicBSplineInterpolateImageFunction::GetWithPadding4D
VoxelType GetWithPadding4D(double, double, double=0, double=0) const
Definition: FastCubicBSplineInterpolateImageFunction.hxx:552
mirtk::GenericFastCubicBSplineInterpolateImageFunction::EvaluateJacobianInside
virtual void EvaluateJacobianInside(Matrix &, double, double, double=0, double=0) const
Get 1st order derivatives of given image at arbitrary location (in pixels)
Definition: FastCubicBSplineInterpolateImageFunction.hxx:1260
mirtk::GenericFastCubicBSplineInterpolateImageFunction::GetWithPaddingInside
virtual VoxelType GetWithPaddingInside(double, double, double=0, double=0) const
Definition: FastCubicBSplineInterpolateImageFunction.hxx:885
mirtk::FillBackgroundBeforeConversionToSplineCoefficients
void FillBackgroundBeforeConversionToSplineCoefficients(GenericImage< TData > &image)
Fill background by front propagation of foreground.
Definition: ImageToInterpolationCoefficients.h:458
mirtk::GenericFastCubicBSplineInterpolateImageFunction::GetInside3D
virtual VoxelType GetInside3D(double, double, double=0, double=0) const
Evaluate generic 3D image without handling boundary conditions.
Definition: FastCubicBSplineInterpolateImageFunction.hxx:786
mirtk::GenericFastCubicBSplineInterpolateImageFunction::Jacobian
void Jacobian(Matrix &, double, double, double=0, double=0) const
Get 1st order derivatives of given image at arbitrary location (in pixels)
Definition: FastCubicBSplineInterpolateImageFunction.hxx:1235
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::fdec
MIRTKCU_API double fdec(double f)
Definition: Math.h:188
mirtk::GenericFastCubicBSplineInterpolateImageFunction::GetWithPadding
virtual VoxelType GetWithPadding(double, double, double=0, double=0) const
Definition: FastCubicBSplineInterpolateImageFunction.hxx:717
mirtk::GenericFastCubicBSplineInterpolateImageFunction::GetWithPaddingOutside
virtual VoxelType GetWithPaddingOutside(double, double, double=0, double=0) const
Definition: FastCubicBSplineInterpolateImageFunction.hxx:898
mirtk::Matrix::Initialize
void Initialize(int, int=-1, double *=NULL)
Initialize matrix with number of rows and columns.
mirtk
Definition: IOConfig.h:41
mirtk::GenericFastCubicBSplineInterpolateImageFunction::GetOutside
virtual VoxelType GetOutside(double, double, double=0, double=0) const
Evaluate generic image at an arbitrary location (in pixels)
Definition: FastCubicBSplineInterpolateImageFunction.hxx:872
mirtk::ExtrapolationMode
ExtrapolationMode
Image extrapolation modes.
Definition: ExtrapolationMode.h:31
mirtk::GenericFastCubicBSplineInterpolateImageFunction::GetWithPadding3D
VoxelType GetWithPadding3D(double, double, double=0, double=0) const
Definition: FastCubicBSplineInterpolateImageFunction.hxx:363
mirtk::AreEqual
MIRTKCU_API bool AreEqual(double a, double b, double tol=1e-12)
Determine equality of two floating point numbers.
Definition: Math.h:115
mirtk::GenericFastCubicBSplineInterpolateImageFunction::Jacobian2D
void Jacobian2D(Matrix &, double, double, double=0, double=0) const
Definition: FastCubicBSplineInterpolateImageFunction.hxx:914
mirtk::GenericFastCubicBSplineInterpolateImageFunction::GetInside
virtual VoxelType GetInside(double, double, double=0, double=0) const
Definition: FastCubicBSplineInterpolateImageFunction.hxx:858
mirtk::iround
MIRTKCU_API int iround(T x)
Round floating-point value and cast to int.
Definition: Math.h:170
mirtk::GenericFastCubicBSplineInterpolateImageFunction::GetInside2D
virtual VoxelType GetInside2D(double, double, double=0, double=0) const
Evaluate generic 2D image without handling boundary conditions.
Definition: FastCubicBSplineInterpolateImageFunction.hxx:757
© 2013-2018 Imperial College London, Andreas Schuh.