SincInterpolateImageFunction.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_SincInterpolateImageFunction_HXX
21 #define MIRTK_SincInterpolateImageFunction_HXX
22 
23 #include "mirtk/SincInterpolateImageFunction.h"
24 #include "mirtk/InterpolateImageFunction.hxx"
25 
26 #include "mirtk/Math.h"
27 #include "mirtk/Sinc.h"
28 #include "mirtk/VoxelCast.h"
29 
30 
31 namespace mirtk {
32 
33 
34 // =============================================================================
35 // Construction/Destruction
36 // =============================================================================
37 
38 // -----------------------------------------------------------------------------
39 template <class TImage>
40 GenericSincInterpolateImageFunction<TImage>
41 ::GenericSincInterpolateImageFunction()
42 :
43  _Epsilon(1e-6)
44 {
45 }
46 
47 // -----------------------------------------------------------------------------
48 template <class TImage>
49 void GenericSincInterpolateImageFunction<TImage>::Initialize(bool coeff)
50 {
51  // Initialize base class
52  Superclass::Initialize(coeff);
53 
54  // Initialize lookup table of Sinc function values
55  Kernel::Initialize();
56 
57  // Domain on which the Sinc interpolation is defined
58  switch (this->NumberOfDimensions()) {
59  case 4:
60  this->_t1 = Kernel::Radius;
61  this->_t2 = this->Input()->T() - Kernel::Radius - 1;
62  case 3:
63  this->_z1 = Kernel::Radius;
64  this->_z2 = this->Input()->Z() - Kernel::Radius - 1;
65  default:
66  this->_y1 = Kernel::Radius;
67  this->_y2 = this->Input()->Y() - Kernel::Radius - 1;
68  this->_x1 = Kernel::Radius;
69  this->_x2 = this->Input()->X() - Kernel::Radius - 1;
70  }
71 }
72 
73 // =============================================================================
74 // Domain checks
75 // =============================================================================
76 
77 // -----------------------------------------------------------------------------
78 template <class TImage>
79 void GenericSincInterpolateImageFunction<TImage>
80 ::BoundingInterval(double x, int &i, int &I) const
81 {
82  const int c = iround(x);
83  i = c - Kernel::Radius;
84  I = c + Kernel::Radius;
85 }
86 
87 // =============================================================================
88 // Evaluation
89 // =============================================================================
90 
91 // -----------------------------------------------------------------------------
92 template <class TImage>
93 inline typename GenericSincInterpolateImageFunction<TImage>::VoxelType
94 GenericSincInterpolateImageFunction<TImage>
95 ::Get2D(double x, double y, double z, double t) const
96 {
97  int i = iround(x);
98  int j = iround(y);
99  int k = iround(z);
100  int l = iround(t);
101 
102  if (k < 0 || k >= this->Input()->Z() ||
103  l < 0 || l >= this->Input()->T()) {
104  return voxel_cast<VoxelType>(this->DefaultValue());
105  }
106 
107  // Return value of nearest neighbor if distance is negligible
108  if (fequal(x, i, _Epsilon) &&
109  fequal(y, j, _Epsilon)) {
110  if (0 <= i && i < this->Input()->X() &&
111  0 <= j && j < this->Input()->Y()) {
112  return this->Input()->Get(i, j, k, l);
113  } else {
114  return voxel_cast<VoxelType>(this->DefaultValue());
115  }
116  }
117 
118  // Truncated Sinc using Hanning window, H(dx/R)*Sinc(dx), R=6 where
119  // Sinc(dx) = sin(pi*dx)/(pi*dx), H(dx/R) = 0.5*(1+cos(pi*dx/R))
120  const int i1 = i - Kernel::Radius;
121  const int j1 = j - Kernel::Radius;
122  const int i2 = i + Kernel::Radius;
123  const int j2 = j + Kernel::Radius;
124 
125  RealType val = voxel_cast<RealType>(0);
126  Real nrm(0), wy, w;
127 
128  for (j = j1; j <= j2; ++j) {
129  if (0 <= j && j < this->Input()->Y()) {
130  wy = Kernel::Lookup(Real(y - j));
131  for (i = i1; i <= i2; ++i) {
132  if (0 <= i && i < this->Input()->X()) {
133  w = Kernel::Lookup(Real(x - i)) * wy;
134  val += w * voxel_cast<RealType>(this->Input()->Get(i, j, k, l));
135  nrm += w;
136  }
137  }
138  }
139  }
140 
141  if (nrm > 1e-3) val /= nrm;
142  else val = voxel_cast<RealType>(this->DefaultValue());
143 
144  return voxel_cast<VoxelType>(val);
145 }
146 
147 // -----------------------------------------------------------------------------
148 template <class TImage>
149 inline typename GenericSincInterpolateImageFunction<TImage>::VoxelType
150 GenericSincInterpolateImageFunction<TImage>
151 ::GetWithPadding2D(double x, double y, double z, double t) const
152 {
153  int i = iround(x);
154  int j = iround(y);
155  int k = iround(z);
156  int l = iround(t);
157 
158  if (k < 0 || k >= this->Input()->Z() ||
159  l < 0 || l >= this->Input()->T()) {
160  return voxel_cast<VoxelType>(this->DefaultValue());
161  }
162 
163  // Return value of nearest neighbor if distance is negligible
164  if (fequal(x, i, _Epsilon) &&
165  fequal(y, j, _Epsilon)) {
166  if (this->Input()->IsInsideForeground(i, j, k, l)) {
167  return this->Input()->Get(i, j, k, l);
168  } else {
169  return voxel_cast<VoxelType>(this->DefaultValue());
170  }
171  }
172 
173  // Truncated Sinc using Hanning window, H(dx/R)*Sinc(dx), R=6 where
174  // Sinc(dx) = sin(pi*dx)/(pi*dx), H(dx/R) = 0.5*(1+cos(pi*dx/R))
175  const int i1 = i - Kernel::Radius;
176  const int j1 = j - Kernel::Radius;
177  const int i2 = i + Kernel::Radius;
178  const int j2 = j + Kernel::Radius;
179 
180  RealType val = voxel_cast<RealType>(0);
181  Real fgw(0), bgw(0), wy, w;
182 
183  for (j = j1; j <= j2; ++j) {
184  wy = Kernel::Lookup(Real(y - j));
185  for (i = i1; i <= i2; ++i) {
186  w = Kernel::Lookup(Real(x - i)) * wy;
187  if (this->Input()->IsInsideForeground(i, j, k, l)) {
188  val += w * voxel_cast<RealType>(this->Input()->Get(i, j, k, l));
189  fgw += w;
190  } else {
191  bgw += w;
192  }
193  }
194  }
195 
196  if (bgw > fgw || AreEqual(bgw, fgw, 1e-3)) {
197  return voxel_cast<VoxelType>(this->DefaultValue());
198  }
199  return voxel_cast<VoxelType>(val / fgw);
200 }
201 
202 // -----------------------------------------------------------------------------
203 template <class TImage> template <class TOtherImage>
204 inline typename TOtherImage::VoxelType
205 GenericSincInterpolateImageFunction<TImage>
206 ::Get2D(const TOtherImage *input, double x, double y, double z, double t) const
207 {
208  typedef typename TOtherImage::VoxelType VoxelType;
209  typedef typename TOtherImage::RealType RealType;
210 
211  int i = iround(x);
212  int j = iround(y);
213  int k = iround(z);
214  int l = iround(t);
215 
216  // Return value of nearest neighbor if distance is negligible
217  if (fequal(x, i, _Epsilon) &&
218  fequal(y, j, _Epsilon)) {
219  return input->Get(i, j, k, l);
220  }
221 
222  // Truncated Sinc using Hanning window, H(dx/R)*Sinc(dx), R=6 where
223  // Sinc(dx) = sin(pi*dx)/(pi*dx), H(dx/R) = 0.5*(1+cos(pi*dx/R))
224  const int i1 = i - Kernel::Radius;
225  const int j1 = j - Kernel::Radius;
226  const int i2 = i + Kernel::Radius;
227  const int j2 = j + Kernel::Radius;
228 
229  RealType val = voxel_cast<RealType>(0);
230  Real nrm(0), wy, w;
231 
232  for (j = j1; j <= j2; ++j) {
233  wy = Kernel::Lookup(Real(y - j));
234  for (i = i1; i <= i2; ++i) {
235  w = Kernel::Lookup(Real(x - i)) * wy;
236  val += w * voxel_cast<RealType>(input->Get(i, j, k, l));
237  nrm += w;
238  }
239  }
240 
241  if (nrm > 1e-3) val /= nrm;
242  else val = voxel_cast<RealType>(this->DefaultValue());
243 
244  return voxel_cast<VoxelType>(val);
245 }
246 
247 // -----------------------------------------------------------------------------
248 template <class TImage> template <class TOtherImage>
249 inline typename TOtherImage::VoxelType
250 GenericSincInterpolateImageFunction<TImage>
251 ::GetWithPadding2D(const TOtherImage *input, double x, double y, double z, double t) const
252 {
253  typedef typename TOtherImage::VoxelType VoxelType;
254  typedef typename TOtherImage::RealType RealType;
255 
256  int i = iround(x);
257  int j = iround(y);
258  int k = iround(z);
259  int l = iround(t);
260 
261  // Return value of nearest neighbor if distance is negligible
262  if (fequal(x, i, _Epsilon) &&
263  fequal(y, j, _Epsilon)) {
264  if (input->IsForeground(i, j, k, l)) {
265  return input->Get(i, j, k, l);
266  } else {
267  return voxel_cast<VoxelType>(this->DefaultValue());
268  }
269  }
270 
271  // Truncated Sinc using Hanning window, H(dx/R)*Sinc(dx), R=6 where
272  // Sinc(dx) = sin(pi*dx)/(pi*dx), H(dx/R) = 0.5*(1+cos(pi*dx/R))
273  const int i1 = i - Kernel::Radius;
274  const int j1 = j - Kernel::Radius;
275  const int i2 = i + Kernel::Radius;
276  const int j2 = j + Kernel::Radius;
277 
278  RealType val = voxel_cast<RealType>(0);
279  Real fgw(0), bgw(0), wy, w;
280 
281  for (j = j1; j <= j2; ++j) {
282  wy = Kernel::Lookup(Real(y - j));
283  for (i = i1; i <= i2; ++i) {
284  w = Kernel::Lookup(Real(x - i)) * wy;
285  if (input->IsForeground(i, j, k, l)) {
286  val += w * voxel_cast<RealType>(input->Get(i, j, k, l));
287  fgw += w;
288  } else {
289  bgw += w;
290  }
291  }
292  }
293 
294  if (bgw > fgw || AreEqual(bgw, fgw, 1e-3)) {
295  return voxel_cast<VoxelType>(this->DefaultValue());
296  }
297  return voxel_cast<VoxelType>(val / fgw);
298 }
299 
300 // -----------------------------------------------------------------------------
301 template <class TImage>
302 inline typename GenericSincInterpolateImageFunction<TImage>::VoxelType
303 GenericSincInterpolateImageFunction<TImage>
304 ::Get3D(double x, double y, double z, double t) const
305 {
306  int i = iround(x);
307  int j = iround(y);
308  int k = iround(z);
309  int l = iround(t);
310 
311  if (l < 0 || l >= this->Input()->T()) {
312  return voxel_cast<VoxelType>(this->DefaultValue());
313  }
314 
315  // Return value of nearest neighbor if distance is negligible
316  if (fequal(x, i, _Epsilon) &&
317  fequal(y, j, _Epsilon) &&
318  fequal(z, k, _Epsilon)) {
319  if (0 <= i && i < this->Input()->X() &&
320  0 <= j && j < this->Input()->Y() &&
321  0 <= k && k < this->Input()->Z()) {
322  return this->Input()->Get(i, j, k, l);
323  } else {
324  return voxel_cast<VoxelType>(this->DefaultValue());
325  }
326  }
327 
328  // Truncated Sinc using Hanning window, H(dx/R)*Sinc(dx), R=6 where
329  // Sinc(dx) = sin(pi*dx)/(pi*dx), H(dx/R) = 0.5*(1+cos(pi*dx/R))
330  const int i1 = i - Kernel::Radius;
331  const int j1 = j - Kernel::Radius;
332  const int k1 = k - Kernel::Radius;
333  const int i2 = i + Kernel::Radius;
334  const int j2 = j + Kernel::Radius;
335  const int k2 = k + Kernel::Radius;
336 
337  RealType val = voxel_cast<RealType>(0);
338  Real nrm(0), wz, wyz, w;
339 
340  for (k = k1; k <= k2; ++k) {
341  if (0 <= k && k < this->Input()->Z()) {
342  wz = Kernel::Lookup(Real(z - k));
343  for (j = j1; j <= j2; ++j) {
344  if (0 <= j && j < this->Input()->Y()) {
345  wyz = Kernel::Lookup(Real(y - j)) * wz;
346  for (i = i1; i <= i2; ++i) {
347  if (0 <= i && i < this->Input()->X()) {
348  w = Kernel::Lookup(Real(x - i)) * wyz;
349  val += w * voxel_cast<RealType>(this->Input()->Get(i, j, k, l));
350  nrm += w;
351  }
352  }
353  }
354  }
355  }
356  }
357 
358  if (nrm > 1e-3) val /= nrm;
359  else val = voxel_cast<RealType>(this->DefaultValue());
360 
361  return voxel_cast<VoxelType>(val);
362 }
363 
364 // -----------------------------------------------------------------------------
365 template <class TImage>
366 inline typename GenericSincInterpolateImageFunction<TImage>::VoxelType
367 GenericSincInterpolateImageFunction<TImage>
368 ::GetWithPadding3D(double x, double y, double z, double t) const
369 {
370  int i = iround(x);
371  int j = iround(y);
372  int k = iround(z);
373  int l = iround(t);
374 
375  if (l < 0 || l >= this->Input()->T()) {
376  return voxel_cast<VoxelType>(this->DefaultValue());
377  }
378 
379  // Return value of nearest neighbor if distance is negligible
380  if (fequal(x, i, _Epsilon) &&
381  fequal(y, j, _Epsilon) &&
382  fequal(z, k, _Epsilon)) {
383  if (this->Input()->IsInsideForeground(i, j, k, l)) {
384  return this->Input()->Get(i, j, k, l);
385  } else {
386  return voxel_cast<VoxelType>(this->DefaultValue());
387  }
388  }
389 
390  // Truncated Sinc using Hanning window, H(dx/R)*Sinc(dx), R=6 where
391  // Sinc(dx) = sin(pi*dx)/(pi*dx), H(dx/R) = 0.5*(1+cos(pi*dx/R))
392  const int i1 = i - Kernel::Radius;
393  const int j1 = j - Kernel::Radius;
394  const int k1 = k - Kernel::Radius;
395  const int i2 = i + Kernel::Radius;
396  const int j2 = j + Kernel::Radius;
397  const int k2 = k + Kernel::Radius;
398 
399  RealType val = voxel_cast<RealType>(0);
400  Real fgw(0), bgw(0), wz, wyz, w;
401 
402  for (k = k1; k <= k2; ++k) {
403  wz = Kernel::Lookup(Real(z - k));
404  for (j = j1; j <= j2; ++j) {
405  wyz = Kernel::Lookup(Real(y - j)) * wz;
406  for (i = i1; i <= i2; ++i) {
407  w = Kernel::Lookup(Real(x - i)) * wyz;
408  if (this->Input()->IsInsideForeground(i, j, k, l)) {
409  val += w * voxel_cast<RealType>(this->Input()->Get(i, j, k, l));
410  fgw += w;
411  } else {
412  bgw += w;
413  }
414  }
415  }
416  }
417 
418  if (bgw > fgw || AreEqual(bgw, fgw, 1e-3)) {
419  return voxel_cast<VoxelType>(this->DefaultValue());
420  }
421  return voxel_cast<VoxelType>(val / fgw);
422 }
423 
424 // -----------------------------------------------------------------------------
425 template <class TImage> template <class TOtherImage>
426 inline typename TOtherImage::VoxelType
427 GenericSincInterpolateImageFunction<TImage>
428 ::Get3D(const TOtherImage *input, double x, double y, double z, double t) const
429 {
430  typedef typename TOtherImage::VoxelType VoxelType;
431  typedef typename TOtherImage::RealType RealType;
432 
433  int i = iround(x);
434  int j = iround(y);
435  int k = iround(z);
436  int l = iround(t);
437 
438  // Return value of nearest neighbor if distance is negligible
439  if (fequal(x, i, _Epsilon) &&
440  fequal(y, j, _Epsilon) &&
441  fequal(z, k, _Epsilon)) {
442  return input->Get(i, j, k, l);
443  }
444 
445  // Truncated Sinc using Hanning window, H(dx/R)*Sinc(dx), R=6 where
446  // Sinc(dx) = sin(pi*dx)/(pi*dx), H(dx/R) = 0.5*(1+cos(pi*dx/R))
447  const int i1 = i - Kernel::Radius;
448  const int j1 = j - Kernel::Radius;
449  const int k1 = k - Kernel::Radius;
450  const int i2 = i + Kernel::Radius;
451  const int j2 = j + Kernel::Radius;
452  const int k2 = k + Kernel::Radius;
453 
454  RealType val = voxel_cast<RealType>(0);
455  Real nrm(0), wz, wyz, w;
456 
457  for (k = k1; k <= k2; ++k) {
458  wz = Kernel::Lookup(Real(z - k));
459  for (j = j1; j <= j2; ++j) {
460  wyz = Kernel::Lookup(Real(y - j)) * wz;
461  for (i = i1; i <= i2; ++i) {
462  w = Kernel::Lookup(Real(x - i)) * wyz;
463  val += w * voxel_cast<RealType>(input->Get(i, j, k, l));
464  nrm += w;
465  }
466  }
467  }
468 
469  if (nrm > 1e-3) val /= nrm;
470  else val = voxel_cast<RealType>(this->DefaultValue());
471 
472  return voxel_cast<VoxelType>(val);
473 }
474 
475 // -----------------------------------------------------------------------------
476 template <class TImage> template <class TOtherImage>
477 inline typename TOtherImage::VoxelType
478 GenericSincInterpolateImageFunction<TImage>
479 ::GetWithPadding3D(const TOtherImage *input, double x, double y, double z, double t) const
480 {
481  typedef typename TOtherImage::VoxelType VoxelType;
482  typedef typename TOtherImage::RealType RealType;
483 
484  int i = iround(x);
485  int j = iround(y);
486  int k = iround(z);
487  int l = iround(t);
488 
489  // Return value of nearest neighbor if distance is negligible
490  if (fequal(x, i, _Epsilon) &&
491  fequal(y, j, _Epsilon) &&
492  fequal(z, k, _Epsilon)) {
493  if (input->IsForeground(i, j, k, l)) {
494  return input->Get(i, j, k, l);
495  } else {
496  return voxel_cast<VoxelType>(this->DefaultValue());
497  }
498  }
499 
500  // Truncated Sinc using Hanning window, H(dx/R)*Sinc(dx), R=6 where
501  // Sinc(dx) = sin(pi*dx)/(pi*dx), H(dx/R) = 0.5*(1+cos(pi*dx/R))
502  const int i1 = i - Kernel::Radius;
503  const int j1 = j - Kernel::Radius;
504  const int k1 = k - Kernel::Radius;
505  const int i2 = i + Kernel::Radius;
506  const int j2 = j + Kernel::Radius;
507  const int k2 = k + Kernel::Radius;
508 
509  RealType val = voxel_cast<RealType>(0);
510  Real fgw(0), bgw(0), wz, wyz, w;
511 
512  for (k = k1; k <= k2; ++k) {
513  wz = Kernel::Lookup(Real(z - k));
514  for (j = j1; j <= j2; ++j) {
515  wyz = Kernel::Lookup(Real(y - j)) * wz;
516  for (i = i1; i <= i2; ++i) {
517  w = Kernel::Lookup(Real(x - i)) * wyz;
518  if (input->IsForeground(i, j, k, l)) {
519  val += w * voxel_cast<RealType>(input->Get(i, j, k, l));
520  fgw += w;
521  } else {
522  bgw += w;
523  }
524  }
525  }
526  }
527 
528  if (bgw > fgw || AreEqual(bgw, fgw, 1e-3)) {
529  return voxel_cast<VoxelType>(this->DefaultValue());
530  }
531  return voxel_cast<VoxelType>(val / fgw);
532 }
533 
534 // -----------------------------------------------------------------------------
535 template <class TImage>
536 inline typename GenericSincInterpolateImageFunction<TImage>::VoxelType
537 GenericSincInterpolateImageFunction<TImage>
538 ::Get4D(double x, double y, double z, double t) const
539 {
540  int i = iround(x);
541  int j = iround(y);
542  int k = iround(z);
543  int l = iround(t);
544 
545  // Return value of nearest neighbor if distance is negligible
546  if (fequal(x, i, _Epsilon) &&
547  fequal(y, j, _Epsilon) &&
548  fequal(z, k, _Epsilon) &&
549  fequal(t, l, _Epsilon)) {
550  if (this->Input()->IsInside(i, j, k, l)) {
551  return this->Input()->Get(i, j, k, l);
552  } else {
553  return voxel_cast<VoxelType>(this->DefaultValue());
554  }
555  }
556 
557  // Truncated Sinc using Hanning window, H(dx/R)*Sinc(dx), R=6 where
558  // Sinc(dx) = sin(pi*dx)/(pi*dx), H(dx/R) = 0.5*(1+cos(pi*dx/R))
559  const int i1 = i - Kernel::Radius;
560  const int j1 = j - Kernel::Radius;
561  const int k1 = k - Kernel::Radius;
562  const int l1 = l - Kernel::Radius;
563  const int i2 = i + Kernel::Radius;
564  const int j2 = j + Kernel::Radius;
565  const int k2 = k + Kernel::Radius;
566  const int l2 = l + Kernel::Radius;
567 
568  RealType val = voxel_cast<RealType>(0);
569  Real nrm(0), wt, wzt, wyzt, w;
570 
571  for (l = l1; l <= l2; ++l) {
572  if (0 <= l && l < this->Input()->T()) {
573  wt = Kernel::Lookup(Real(t - l));
574  for (k = k1; k <= k2; ++k) {
575  if (0 <= k && k < this->Input()->Z()) {
576  wzt = Kernel::Lookup(Real(z - k)) * wt;
577  for (j = j1; j <= j2; ++j) {
578  if (0 <= j && j < this->Input()->Y()) {
579  wyzt = Kernel::Lookup(Real(y - j)) * wzt;
580  for (i = i1; i <= i2; ++i) {
581  if (0 <= i && i < this->Input()->X()) {
582  w = Kernel::Lookup(Real(x - i)) * wyzt;
583  val += w * voxel_cast<RealType>(this->Input()->Get(i, j, k, l));
584  nrm += w;
585  }
586  }
587  }
588  }
589  }
590  }
591  }
592  }
593 
594  if (nrm > 1e-3) val /= nrm;
595  else val = voxel_cast<RealType>(this->DefaultValue());
596 
597  return voxel_cast<VoxelType>(val);
598 }
599 
600 // -----------------------------------------------------------------------------
601 template <class TImage>
602 inline typename GenericSincInterpolateImageFunction<TImage>::VoxelType
603 GenericSincInterpolateImageFunction<TImage>
604 ::GetWithPadding4D(double x, double y, double z, double t) const
605 {
606  int i = iround(x);
607  int j = iround(y);
608  int k = iround(z);
609  int l = iround(t);
610 
611  // Return value of nearest neighbor if distance is negligible
612  if (fequal(x, i, _Epsilon) &&
613  fequal(y, j, _Epsilon) &&
614  fequal(z, k, _Epsilon) &&
615  fequal(t, l, _Epsilon)) {
616  if (this->Input()->IsInsideForeground(i, j, k, l)) {
617  return this->Input()->Get(i, j, k, l);
618  } else {
619  return voxel_cast<VoxelType>(this->DefaultValue());
620  }
621  }
622 
623  // Truncated Sinc using Hanning window, H(dx/R)*Sinc(dx), R=6 where
624  // Sinc(dx) = sin(pi*dx)/(pi*dx), H(dx/R) = 0.5*(1+cos(pi*dx/R))
625  const int i1 = i - Kernel::Radius;
626  const int j1 = j - Kernel::Radius;
627  const int k1 = k - Kernel::Radius;
628  const int l1 = l - Kernel::Radius;
629  const int i2 = i + Kernel::Radius;
630  const int j2 = j + Kernel::Radius;
631  const int k2 = k + Kernel::Radius;
632  const int l2 = l + Kernel::Radius;
633 
634  RealType val = voxel_cast<RealType>(0);
635  Real fgw(0), bgw(0), wt, wzt, wyzt, w;
636 
637  for (l = l1; l <= l2; ++l) {
638  wt = Kernel::Lookup(Real(t - l));
639  for (k = k1; k <= k2; ++k) {
640  wzt = Kernel::Lookup(Real(z - k)) * wt;
641  for (j = j1; j <= j2; ++j) {
642  wyzt = Kernel::Lookup(Real(y - j)) * wzt;
643  for (i = i1; i <= i2; ++i) {
644  w = Kernel::Lookup(Real(x - i)) * wyzt;
645  if (this->Input()->IsInsideForeground(i, j, k, l)) {
646  val += w * voxel_cast<RealType>(this->Input()->Get(i, j, k, l));
647  fgw += w;
648  } else {
649  bgw += w;
650  }
651  }
652  }
653  }
654  }
655 
656  if (bgw > fgw || AreEqual(bgw, fgw, 1e-3)) {
657  return voxel_cast<VoxelType>(this->DefaultValue());
658  }
659  return voxel_cast<VoxelType>(val / fgw);
660 }
661 
662 // -----------------------------------------------------------------------------
663 template <class TImage> template <class TOtherImage>
664 inline typename TOtherImage::VoxelType
665 GenericSincInterpolateImageFunction<TImage>
666 ::Get4D(const TOtherImage *input, double x, double y, double z, double t) const
667 {
668  typedef typename TOtherImage::VoxelType VoxelType;
669  typedef typename TOtherImage::RealType RealType;
670 
671  int i = iround(x);
672  int j = iround(y);
673  int k = iround(z);
674  int l = iround(t);
675 
676  // Return value of nearest neighbor if distance is negligible
677  if (fequal(x, i, _Epsilon) &&
678  fequal(y, j, _Epsilon) &&
679  fequal(z, k, _Epsilon) &&
680  fequal(t, l, _Epsilon)) {
681  return input->Get(i, j, k, l);
682  }
683 
684  // Truncated Sinc using Hanning window, H(dx/R)*Sinc(dx), R=6 where
685  // Sinc(dx) = sin(pi*dx)/(pi*dx), H(dx/R) = 0.5*(1+cos(pi*dx/R))
686  const int i1 = i - Kernel::Radius;
687  const int j1 = j - Kernel::Radius;
688  const int k1 = k - Kernel::Radius;
689  const int l1 = l - Kernel::Radius;
690  const int i2 = i + Kernel::Radius;
691  const int j2 = j + Kernel::Radius;
692  const int k2 = k + Kernel::Radius;
693  const int l2 = l + Kernel::Radius;
694 
695  RealType val = voxel_cast<RealType>(0);
696  Real nrm(0), wt, wzt, wyzt, w;
697 
698  for (l = l1; l <= l2; ++l) {
699  wt = Kernel::Lookup(Real(t - l));
700  for (k = k1; k <= k2; ++k) {
701  wzt = Kernel::Lookup(Real(z - k)) * wt;
702  for (j = j1; j <= j2; ++j) {
703  wyzt = Kernel::Lookup(Real(y - j)) * wzt;
704  for (i = i1; i <= i2; ++i) {
705  w = Kernel::Lookup(Real(x - i)) * wyzt;
706  val += w * voxel_cast<RealType>(input->Get(i, j, k, l));
707  nrm += w;
708  }
709  }
710  }
711  }
712 
713  if (nrm > 1e-3) val /= nrm;
714  else val = voxel_cast<RealType>(this->DefaultValue());
715 
716  return voxel_cast<VoxelType>(val);
717 }
718 
719 // -----------------------------------------------------------------------------
720 template <class TImage> template <class TOtherImage>
721 inline typename TOtherImage::VoxelType
722 GenericSincInterpolateImageFunction<TImage>
723 ::GetWithPadding4D(const TOtherImage *input, double x, double y, double z, double t) const
724 {
725  typedef typename TOtherImage::VoxelType VoxelType;
726  typedef typename TOtherImage::RealType RealType;
727 
728  int i = iround(x);
729  int j = iround(y);
730  int k = iround(z);
731  int l = iround(t);
732 
733  // Return value of nearest neighbor if distance is negligible
734  if (fequal(x, i, _Epsilon) &&
735  fequal(y, j, _Epsilon) &&
736  fequal(z, k, _Epsilon) &&
737  fequal(t, l, _Epsilon)) {
738  if (input->IsForeground(i, j, k, l)) {
739  return input->Get(i, j, k, l);
740  } else {
741  return voxel_cast<VoxelType>(this->DefaultValue());
742  }
743  }
744 
745  // Truncated Sinc using Hanning window, H(dx/R)*Sinc(dx), R=6 where
746  // Sinc(dx) = sin(pi*dx)/(pi*dx), H(dx/R) = 0.5*(1+cos(pi*dx/R))
747  const int i1 = i - Kernel::Radius;
748  const int j1 = j - Kernel::Radius;
749  const int k1 = k - Kernel::Radius;
750  const int l1 = l - Kernel::Radius;
751  const int i2 = i + Kernel::Radius;
752  const int j2 = j + Kernel::Radius;
753  const int k2 = k + Kernel::Radius;
754  const int l2 = l + Kernel::Radius;
755 
756  RealType val = voxel_cast<RealType>(0);
757  Real fgw(0), bgw(0), wt, wzt, wyzt, w;
758 
759  for (l = l1; l <= l2; ++l) {
760  wt = Kernel::Lookup(Real(t - l));
761  for (k = k1; k <= k2; ++k) {
762  wzt = Kernel::Lookup(Real(z - k)) * wt;
763  for (j = j1; j <= j2; ++j) {
764  wyzt = Kernel::Lookup(Real(y - j)) * wzt;
765  for (i = i1; i <= i2; ++i) {
766  w = Kernel::Lookup(Real(x - i)) * wyzt;
767  if (input->IsForeground(i, j, k, l)) {
768  val += w * voxel_cast<RealType>(input->Get(i, j, k, l));
769  fgw += w;
770  } else {
771  bgw += w;
772  }
773  }
774  }
775  }
776  }
777 
778  if (bgw > fgw || AreEqual(bgw, fgw, 1e-3)) {
779  return voxel_cast<VoxelType>(this->DefaultValue());
780  }
781  return voxel_cast<VoxelType>(val / fgw);
782 }
783 
784 // -----------------------------------------------------------------------------
785 template <class TImage>
786 inline typename GenericSincInterpolateImageFunction<TImage>::VoxelType
787 GenericSincInterpolateImageFunction<TImage>
788 ::Get(double x, double y, double z, double t) const
789 {
790  switch (this->NumberOfDimensions()) {
791  case 3: return Get3D(x, y, z, t);
792  case 2: return Get2D(x, y, z, t);
793  default: return Get4D(x, y, z, t);
794  }
795 }
796 
797 // -----------------------------------------------------------------------------
798 template <class TImage>
799 inline typename GenericSincInterpolateImageFunction<TImage>::VoxelType
800 GenericSincInterpolateImageFunction<TImage>
801 ::GetWithPadding(double x, double y, double z, double t) const
802 {
803  switch (this->NumberOfDimensions()) {
804  case 3: return GetWithPadding3D(x, y, z, t);
805  case 2: return GetWithPadding2D(x, y, z, t);
806  default: return GetWithPadding4D(x, y, z, t);
807  }
808 }
809 
810 // -----------------------------------------------------------------------------
811 template <class TImage> template <class TOtherImage>
812 inline typename TOtherImage::VoxelType
813 GenericSincInterpolateImageFunction<TImage>
814 ::Get(const TOtherImage *input, double x, double y, double z, double t) const
815 {
816  switch (this->NumberOfDimensions()) {
817  case 3: return Get3D(input, x, y, z, t);
818  case 2: return Get2D(input, x, y, z, t);
819  default: return Get4D(input, x, y, z, t);
820  }
821 }
822 
823 // -----------------------------------------------------------------------------
824 template <class TImage> template <class TOtherImage>
825 inline typename TOtherImage::VoxelType
826 GenericSincInterpolateImageFunction<TImage>
827 ::GetWithPadding(const TOtherImage *input, double x, double y, double z, double t) const
828 {
829  switch (this->NumberOfDimensions()) {
830  case 3: return GetWithPadding3D(input, x, y, z, t);
831  case 2: return GetWithPadding2D(input, x, y, z, t);
832  default: return GetWithPadding4D(input, x, y, z, t);
833  }
834 }
835 
836 // -----------------------------------------------------------------------------
837 template <class TImage>
838 inline typename GenericSincInterpolateImageFunction<TImage>::VoxelType
839 GenericSincInterpolateImageFunction<TImage>
840 ::GetInside(double x, double y, double z, double t) const
841 {
842  return Get(this->Input(), x, y, z, t);
843 }
844 
845 // -----------------------------------------------------------------------------
846 template <class TImage>
847 inline typename GenericSincInterpolateImageFunction<TImage>::VoxelType
848 GenericSincInterpolateImageFunction<TImage>
849 ::GetOutside(double x, double y, double z, double t) const
850 {
851  if (this->Extrapolator()) {
852  return Get(this->Extrapolator(), x, y, z, t);
853  } else {
854  return Get(x, y, z, t);
855  }
856 }
857 
858 // -----------------------------------------------------------------------------
859 template <class TImage>
860 inline typename GenericSincInterpolateImageFunction<TImage>::VoxelType
861 GenericSincInterpolateImageFunction<TImage>
862 ::GetWithPaddingInside(double x, double y, double z, double t) const
863 {
864  return GetWithPadding(this->Input(), x, y, z, t);
865 }
866 
867 // -----------------------------------------------------------------------------
868 template <class TImage>
869 inline typename GenericSincInterpolateImageFunction<TImage>::VoxelType
870 GenericSincInterpolateImageFunction<TImage>
871 ::GetWithPaddingOutside(double x, double y, double z, double t) const
872 {
873  if (this->Extrapolator()) {
874  return GetWithPadding(this->Extrapolator(), x, y, z, t);
875  } else {
876  return GetWithPadding(x, y, z, t);
877  }
878 }
879 
880 
881 } // namespace mirtk
882 
883 #endif // MIRTK_SincInterpolateImageFunction_HXX
mirtk::Get
string Get(const ParameterList &params, string name)
Get parameter value from parameters list.
Definition: Object.h:202
mirtk::GenericSincInterpolateImageFunction::GetWithPadding2D
VoxelType GetWithPadding2D(double, double, double=0, double=0) const
Definition: SincInterpolateImageFunction.hxx:151
mirtk::GenericSincInterpolateImageFunction::Get4D
VoxelType Get4D(double, double, double=0, double=0) const
Definition: SincInterpolateImageFunction.hxx:538
mirtk::GenericSincInterpolateImageFunction::GetWithPadding4D
VoxelType GetWithPadding4D(double, double, double=0, double=0) const
Definition: SincInterpolateImageFunction.hxx:604
mirtk::GenericSincInterpolateImageFunction::GetWithPaddingInside
virtual VoxelType GetWithPaddingInside(double, double, double=0, double=0) const
Definition: SincInterpolateImageFunction.hxx:862
mirtk::fequal
MIRTKCU_API bool fequal(double a, double b, double tol=1e-12)
Definition: Math.h:138
mirtk
Definition: IOConfig.h:41
mirtk::GenericSincInterpolateImageFunction::Get2D
VoxelType Get2D(double, double, double=0, double=0) const
Definition: SincInterpolateImageFunction.hxx:95
mirtk::GenericSincInterpolateImageFunction::GetOutside
virtual VoxelType GetOutside(double, double, double=0, double=0) const
Evaluate generic image at an arbitrary location (in pixels)
Definition: SincInterpolateImageFunction.hxx:849
mirtk::GenericSincInterpolateImageFunction::GetInside
virtual VoxelType GetInside(double, double, double=0, double=0) const
Definition: SincInterpolateImageFunction.hxx:840
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::GenericSincInterpolateImageFunction::GetWithPadding
virtual VoxelType GetWithPadding(double, double, double=0, double=0) const
Definition: SincInterpolateImageFunction.hxx:801
mirtk::GenericSincInterpolateImageFunction::Get
VoxelType Get(double, double, double=0, double=0) const
Definition: SincInterpolateImageFunction.hxx:788
mirtk::iround
MIRTKCU_API int iround(T x)
Round floating-point value and cast to int.
Definition: Math.h:170
mirtk::GenericSincInterpolateImageFunction::GetWithPaddingOutside
virtual VoxelType GetWithPaddingOutside(double, double, double=0, double=0) const
Definition: SincInterpolateImageFunction.hxx:871
mirtk::GenericSincInterpolateImageFunction::BoundingInterval
virtual void BoundingInterval(double, int &, int &) const
Definition: SincInterpolateImageFunction.hxx:80
mirtk::GenericSincInterpolateImageFunction::GetWithPadding3D
VoxelType GetWithPadding3D(double, double, double=0, double=0) const
Definition: SincInterpolateImageFunction.hxx:368
mirtk::GenericSincInterpolateImageFunction::Get3D
VoxelType Get3D(double, double, double=0, double=0) const
Definition: SincInterpolateImageFunction.hxx:304
© 2013-2018 Imperial College London, Andreas Schuh.