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>
GenericSincInterpolateImageFunction()41 ::GenericSincInterpolateImageFunction()
42 :
43 _Epsilon(1e-6)
44 {
45 }
46
47 // -----------------------------------------------------------------------------
48 template <class TImage>
Initialize(bool coeff)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>
BoundingInterval(double x,int & i,int & I) const80 ::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>
Get2D(double x,double y,double z,double t) const95 ::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>
GetWithPadding2D(double x,double y,double z,double t) const151 ::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>
Get2D(const TOtherImage * input,double x,double y,double z,double t) const206 ::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>
GetWithPadding2D(const TOtherImage * input,double x,double y,double z,double t) const251 ::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>
Get3D(double x,double y,double z,double t) const304 ::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>
GetWithPadding3D(double x,double y,double z,double t) const368 ::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>
Get3D(const TOtherImage * input,double x,double y,double z,double t) const428 ::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>
GetWithPadding3D(const TOtherImage * input,double x,double y,double z,double t) const479 ::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>
Get4D(double x,double y,double z,double t) const538 ::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>
GetWithPadding4D(double x,double y,double z,double t) const604 ::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>
Get4D(const TOtherImage * input,double x,double y,double z,double t) const666 ::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>
GetWithPadding4D(const TOtherImage * input,double x,double y,double z,double t) const723 ::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>
Get(double x,double y,double z,double t) const788 ::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>
GetWithPadding(double x,double y,double z,double t) const801 ::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>
Get(const TOtherImage * input,double x,double y,double z,double t) const814 ::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>
GetWithPadding(const TOtherImage * input,double x,double y,double z,double t) const827 ::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>
GetInside(double x,double y,double z,double t) const840 ::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>
GetOutside(double x,double y,double z,double t) const849 ::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>
GetWithPaddingInside(double x,double y,double z,double t) const862 ::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>
GetWithPaddingOutside(double x,double y,double z,double t) const871 ::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
884