LinearInterpolateImageFunction.hxx
1 /*
2  * Medical Image Registration ToolKit (MIRTK)
3  *
4  * Copyright 2013-2017 Imperial College London
5  * Copyright 2013-2017 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_LinearInterpolateImageFunction_HXX
21 #define MIRTK_LinearInterpolateImageFunction_HXX
22 
23 #include "mirtk/LinearInterpolateImageFunction.h"
24 #include "mirtk/InterpolateImageFunction.hxx"
25 
26 #include "mirtk/Math.h"
27 #include "mirtk/VoxelCast.h"
28 
29 
30 namespace mirtk {
31 
32 
33 // =============================================================================
34 // Construction/Destruction
35 // =============================================================================
36 
37 // -----------------------------------------------------------------------------
38 template <class TImage>
39 GenericLinearInterpolateImageFunction<TImage>
40 ::GenericLinearInterpolateImageFunction()
41 {
42 }
43 
44 // -----------------------------------------------------------------------------
45 template <class TImage>
46 GenericLinearInterpolateImageFunction<TImage>
47 ::~GenericLinearInterpolateImageFunction()
48 {
49 }
50 
51 // -----------------------------------------------------------------------------
52 template <class TImage>
53 void GenericLinearInterpolateImageFunction<TImage>::Initialize(bool coeff)
54 {
55  /// Initialize base class
56  Superclass::Initialize(coeff);
57 
58  // Domain on which linear interpolation is defined: [0, x)
59  switch (this->NumberOfDimensions()) {
60  case 4:
61  this->_t1 = 0.;
62  this->_t2 = fdec(this->Input()->T() - 1);
63  case 3:
64  this->_z1 = 0.;
65  this->_z2 = fdec(this->Input()->Z() - 1);
66  default:
67  this->_y1 = 0.;
68  this->_y2 = fdec(this->Input()->Y() - 1);
69  this->_x1 = 0.;
70  this->_x2 = fdec(this->Input()->X() - 1);
71  }
72 
73  // Calculate offsets for fast pixel access
74  const int xoffset[2] = {0, 1};
75  const int yoffset[2] = {0, this->Input()->X()};
76  const int zoffset[2] = {0, this->Input()->X() * this->Input()->Y()};
77  const int toffset[2] = {0, this->Input()->X() * this->Input()->Y() * this->Input()->Z()};
78 
79  int idx = 0;
80  for (int d = 0; d <= 1; ++d)
81  for (int c = 0; c <= 1; ++c)
82  for (int b = 0; b <= 1; ++b)
83  for (int a = 0; a <= 1; ++a, ++idx) {
84  _Offset[idx] = xoffset[a] + yoffset[b] + zoffset[c] + toffset[d];
85  }
86 }
87 
88 // =============================================================================
89 // Interpolation weights
90 // =============================================================================
91 
92 // -----------------------------------------------------------------------------
93 template <class TImage>
94 inline int GenericLinearInterpolateImageFunction<TImage>
95 ::ComputeWeights(double x, Real w[2])
96 {
97  const int i = ifloor(x);
98  w[1] = static_cast<Real>(x) - static_cast<Real>(i);
99  w[0] = static_cast<Real>(1) - w[1];
100  return i;
101 }
102 
103 // =============================================================================
104 // Domain checks
105 // =============================================================================
106 
107 // -----------------------------------------------------------------------------
108 template <class TImage>
109 void GenericLinearInterpolateImageFunction<TImage>
110 ::BoundingInterval(double x, int &i, int &I) const
111 {
112  i = ifloor(x), I = i + 1;
113 }
114 
115 // =============================================================================
116 // Evaluation
117 // =============================================================================
118 
119 // -----------------------------------------------------------------------------
120 template <class TImage>
121 inline typename GenericLinearInterpolateImageFunction<TImage>::VoxelType
122 GenericLinearInterpolateImageFunction<TImage>
123 ::Get2D(double x, double y, double z, double t) const
124 {
125  const int k = iround(z);
126  const int l = iround(t);
127  if (k < 0 || k >= this->Input()->Z() ||
128  l < 0 || l >= this->Input()->T()) {
129  return voxel_cast<VoxelType>(this->DefaultValue());
130  }
131 
132  Real wx[2], wy[2];
133  const int i = ComputeWeights(x, wx);
134  const int j = ComputeWeights(y, wy);
135 
136  RealType val = voxel_cast<RealType>(0);
137  Real nrm(0), w;
138 
139  int ia, jb;
140  for (int b = 0; b <= 1; ++b) {
141  jb = j + b;
142  if (0 <= jb && jb < this->Input()->Y()) {
143  for (int a = 0; a <= 1; ++a) {
144  ia = i + a;
145  if (0 <= ia && ia < this->Input()->X()) {
146  w = wx[a] * wy[b];
147  val += w * voxel_cast<RealType>(this->Input()->Get(ia, jb, k, l));
148  nrm += w;
149  }
150  }
151  }
152  }
153 
154  if (nrm > Real(1e-3)) val /= nrm;
155  else val = voxel_cast<RealType>(this->DefaultValue());
156 
157  return voxel_cast<VoxelType>(val);
158 }
159 
160 // -----------------------------------------------------------------------------
161 template <class TImage>
162 inline typename GenericLinearInterpolateImageFunction<TImage>::VoxelType
163 GenericLinearInterpolateImageFunction<TImage>
164 ::GetWithPadding2D(double x, double y, double z, double t) const
165 {
166  const int k = iround(z);
167  const int l = iround(t);
168  if (k < 0 || k >= this->Input()->Z() ||
169  l < 0 || l >= this->Input()->T()) {
170  return voxel_cast<VoxelType>(this->DefaultValue());
171  }
172 
173  Real wx[2], wy[2];
174  const int i = ComputeWeights(x, wx);
175  const int j = ComputeWeights(y, wy);
176 
177  RealType val = voxel_cast<RealType>(0);
178  Real fgw(0), w;
179 
180  int ia, jb;
181  for (int b = 0; b <= 1; ++b) {
182  jb = j + b;
183  for (int a = 0; a <= 1; ++a) {
184  ia = i + a;
185  w = wx[a] * wy[b];
186  if (this->Input()->IsInsideForeground(ia, jb, k, l)) {
187  val += w * voxel_cast<RealType>(this->Input()->Get(ia, jb, k, l));
188  fgw += w;
189  }
190  }
191  }
192 
193  if (ifloor(fgw * 1.e3) > 500) val /= fgw;
194  else val = voxel_cast<RealType>(this->DefaultValue());
195 
196  return voxel_cast<VoxelType>(val);
197 }
198 
199 // -----------------------------------------------------------------------------
200 template <class TImage> template <class TOtherImage>
201 inline typename TOtherImage::VoxelType
202 GenericLinearInterpolateImageFunction<TImage>
203 ::Get2D(const TOtherImage *input, double x, double y, double z, double t) const
204 {
205  typedef typename TOtherImage::VoxelType VoxelType;
206  typedef typename TOtherImage::RealType RealType;
207 
208  Real wx[2], wy[2];
209  const int i = ComputeWeights(x, wx), I = i + 1;
210  const int j = ComputeWeights(y, wy), J = j + 1;
211  const int k = iround(z);
212  const int l = iround(t);
213 
214  const auto v00 = voxel_cast<RealType>(input->Get(i, j, k, l));
215  const auto v01 = voxel_cast<RealType>(input->Get(I, j, k, l));
216  const auto v10 = voxel_cast<RealType>(input->Get(i, J, k, l));
217  const auto v11 = voxel_cast<RealType>(input->Get(I, J, k, l));
218 
219  return voxel_cast<VoxelType>(wy[0] * (wx[0] * v00 + wx[1] * v01) +
220  wy[1] * (wx[0] * v10 + wx[1] * v11));
221 }
222 
223 // -----------------------------------------------------------------------------
224 template <class TImage> template <class TOtherImage>
225 inline typename TOtherImage::VoxelType
226 GenericLinearInterpolateImageFunction<TImage>
227 ::GetWithPadding2D(const TOtherImage *input, double x, double y, double z, double t) const
228 {
229  typedef typename TOtherImage::VoxelType VoxelType;
230  typedef typename TOtherImage::RealType RealType;
231 
232  const int k = iround(z);
233  const int l = iround(t);
234  if (k < 0 || k >= input->Z() ||
235  l < 0 || l >= input->T()) {
236  return voxel_cast<VoxelType>(this->DefaultValue());
237  }
238 
239  Real wx[2], wy[2];
240  const int i = ComputeWeights(x, wx);
241  const int j = ComputeWeights(y, wy);
242 
243  RealType val = voxel_cast<RealType>(0);
244  Real fgw(0), w;
245 
246  int ia, jb;
247  for (int b = 0; b <= 1; ++b) {
248  jb = j + b;
249  for (int a = 0; a <= 1; ++a) {
250  ia = i + a;
251  w = wx[a] * wy[b];
252  if (input->IsForeground(ia, jb, k, l)) {
253  val += w * voxel_cast<RealType>(input->Get(ia, jb, k, l));
254  fgw += w;
255  }
256  }
257  }
258 
259  if (ifloor(fgw * 1.e3) > 500) val /= fgw;
260  else val = voxel_cast<RealType>(this->DefaultValue());
261 
262  return voxel_cast<VoxelType>(val);
263 }
264 
265 // -----------------------------------------------------------------------------
266 template <class TImage>
267 inline typename GenericLinearInterpolateImageFunction<TImage>::VoxelType
268 GenericLinearInterpolateImageFunction<TImage>
269 ::Get3D(double x, double y, double z, double t) const
270 {
271  const int l = iround(t);
272  if (l < 0 || l >= this->Input()->T()) {
273  return voxel_cast<VoxelType>(this->DefaultValue());
274  }
275 
276  Real wx[2], wy[2], wz[2];
277  const int i = ComputeWeights(x, wx);
278  const int j = ComputeWeights(y, wy);
279  const int k = ComputeWeights(z, wz);
280 
281  RealType val = voxel_cast<RealType>(0);
282  Real nrm(0), w;
283 
284  int ia, jb, kc;
285  for (int c = 0; c <= 1; ++c) {
286  kc = k + c;
287  if (0 <= kc && kc < this->Input()->Z()) {
288  for (int b = 0; b <= 1; ++b) {
289  jb = j + b;
290  if (0 <= jb && jb < this->Input()->Y()) {
291  for (int a = 0; a <= 1; ++a) {
292  ia = i + a;
293  if (0 <= ia && ia < this->Input()->X()) {
294  w = wx[a] * wy[b] * wz[c];
295  val += w * voxel_cast<RealType>(this->Input()->Get(ia, jb, kc, l));
296  nrm += w;
297  }
298  }
299  }
300  }
301  }
302  }
303 
304  if (nrm > Real(1e-3)) val /= nrm;
305  else val = voxel_cast<RealType>(this->DefaultValue());
306 
307  return voxel_cast<VoxelType>(val);
308 }
309 
310 // -----------------------------------------------------------------------------
311 template <class TImage>
312 inline typename GenericLinearInterpolateImageFunction<TImage>::VoxelType
313 GenericLinearInterpolateImageFunction<TImage>
314 ::GetWithPadding3D(double x, double y, double z, double t) const
315 {
316  const int l = iround(t);
317  if (l < 0 || l >= this->Input()->T()) {
318  return voxel_cast<VoxelType>(this->DefaultValue());
319  }
320 
321  Real wx[2], wy[2], wz[2];
322  const int i = ComputeWeights(x, wx);
323  const int j = ComputeWeights(y, wy);
324  const int k = ComputeWeights(z, wz);
325 
326  RealType val = voxel_cast<RealType>(0);
327  Real fgw(0), w;
328 
329  int ia, jb, kc;
330  for (int c = 0; c <= 1; ++c) {
331  kc = k + c;
332  for (int b = 0; b <= 1; ++b) {
333  jb = j + b;
334  for (int a = 0; a <= 1; ++a) {
335  ia = i + a;
336  w = wx[a] * wy[b] * wz[c];
337  if (this->Input()->IsInsideForeground(ia, jb, kc, l)) {
338  val += w * voxel_cast<RealType>(this->Input()->Get(ia, jb, kc, l));
339  fgw += w;
340  }
341  }
342  }
343  }
344 
345  if (ifloor(fgw * 1.e3) > 500) val /= fgw;
346  else val = voxel_cast<RealType>(this->DefaultValue());
347 
348  return voxel_cast<VoxelType>(val);
349 }
350 
351 // -----------------------------------------------------------------------------
352 template <class TImage> template <class TOtherImage>
353 inline typename TOtherImage::VoxelType
354 GenericLinearInterpolateImageFunction<TImage>
355 ::Get3D(const TOtherImage *input, double x, double y, double z, double t) const
356 {
357  typedef typename TOtherImage::VoxelType VoxelType;
358  typedef typename TOtherImage::RealType RealType;
359 
360  Real wx[2], wy[2], wz[2];
361  const int i = ComputeWeights(x, wx), I = i + 1;
362  const int j = ComputeWeights(y, wy), J = j + 1;
363  const int k = ComputeWeights(z, wz), K = k + 1;
364  const int l = iround(t);
365 
366  const auto v000 = voxel_cast<RealType>(input->Get(i, j, k, l));
367  const auto v001 = voxel_cast<RealType>(input->Get(I, j, k, l));
368  const auto v010 = voxel_cast<RealType>(input->Get(i, J, k, l));
369  const auto v011 = voxel_cast<RealType>(input->Get(I, J, k, l));
370  const auto v100 = voxel_cast<RealType>(input->Get(i, j, K, l));
371  const auto v101 = voxel_cast<RealType>(input->Get(I, j, K, l));
372  const auto v110 = voxel_cast<RealType>(input->Get(i, J, K, l));
373  const auto v111 = voxel_cast<RealType>(input->Get(I, J, K, l));
374 
375  return voxel_cast<VoxelType>(wz[0] * (wy[0] * (wx[0] * v000 + wx[1] * v001) +
376  wy[1] * (wx[0] * v010 + wx[1] * v011)) +
377  wz[1] * (wy[0] * (wx[0] * v100 + wx[1] * v101) +
378  wy[1] * (wx[0] * v110 + wx[1] * v111)));
379 }
380 
381 // -----------------------------------------------------------------------------
382 template <class TImage> template <class TOtherImage>
383 inline typename TOtherImage::VoxelType
384 GenericLinearInterpolateImageFunction<TImage>
385 ::GetWithPadding3D(const TOtherImage *input, double x, double y, double z, double t) const
386 {
387  typedef typename TOtherImage::VoxelType VoxelType;
388  typedef typename TOtherImage::RealType RealType;
389 
390  const int l = iround(t);
391  if (l < 0 || l >= input->T()) {
392  return voxel_cast<VoxelType>(this->DefaultValue());
393  }
394 
395  Real wx[2], wy[2], wz[2];
396  const int i = ComputeWeights(x, wx);
397  const int j = ComputeWeights(y, wy);
398  const int k = ComputeWeights(z, wz);
399 
400  RealType val = voxel_cast<RealType>(0);
401  Real fgw(0), wyz, w;
402 
403  int ia, jb, kc;
404  for (int c = 0; c <= 1; ++c) {
405  kc = k + c;
406  for (int b = 0; b <= 1; ++b) {
407  jb = j + b;
408  wyz = wy[b] * wz[c];
409  for (int a = 0; a <= 1; ++a) {
410  ia = i + a;
411  w = wx[a] * wyz;
412  if (input->IsForeground(ia, jb, kc, l)) {
413  val += w * static_cast<RealType>(input->Get(ia, jb, kc, l));
414  fgw += w;
415  }
416  }
417  }
418  }
419 
420  if (ifloor(fgw * 1.e3) > 500) val /= fgw;
421  else val = voxel_cast<RealType>(this->DefaultValue());
422 
423  return voxel_cast<VoxelType>(val);
424 }
425 
426 // -----------------------------------------------------------------------------
427 template <class TImage>
428 inline typename GenericLinearInterpolateImageFunction<TImage>::VoxelType
429 GenericLinearInterpolateImageFunction<TImage>
430 ::Get4D(double x, double y, double z, double t) const
431 {
432  Real wx[2], wy[2], wz[2], wt[2];
433  const int i = ComputeWeights(x, wx);
434  const int j = ComputeWeights(y, wy);
435  const int k = ComputeWeights(z, wz);
436  const int l = ComputeWeights(t, wt);
437 
438  RealType val = voxel_cast<RealType>(0);
439  Real nrm(0), w;
440 
441  int ia, jb, kc, ld;
442  for (int d = 0; d <= 1; ++d) {
443  ld = l + d;
444  if (0 <= ld && ld < this->Input()->T()) {
445  for (int c = 0; c <= 1; ++c) {
446  kc = k + c;
447  if (0 <= kc && kc < this->Input()->Z()) {
448  for (int b = 0; b <= 1; ++b) {
449  jb = j + b;
450  if (0 <= jb && jb < this->Input()->Y()) {
451  for (int a = 0; a <= 1; ++a) {
452  ia = i + a;
453  if (0 <= ia && ia < this->Input()->X()) {
454  w = wx[a] * wy[b] * wz[c] * wt[d];
455  val += w * voxel_cast<RealType>(this->Input()->Get(ia, jb, kc, ld));
456  nrm += w;
457  }
458  }
459  }
460  }
461  }
462  }
463  }
464  }
465 
466  if (nrm > Real(1e-3)) val /= nrm;
467  else val = voxel_cast<RealType>(this->DefaultValue());
468 
469  return voxel_cast<VoxelType>(val);
470 }
471 
472 // -----------------------------------------------------------------------------
473 template <class TImage>
474 inline typename GenericLinearInterpolateImageFunction<TImage>::VoxelType
475 GenericLinearInterpolateImageFunction<TImage>
476 ::GetWithPadding4D(double x, double y, double z, double t) const
477 {
478  Real wx[2], wy[2], wz[2], wt[2];
479  const int i = ComputeWeights(x, wx);
480  const int j = ComputeWeights(y, wy);
481  const int k = ComputeWeights(z, wz);
482  const int l = ComputeWeights(t, wt);
483 
484  RealType val = voxel_cast<RealType>(0);
485  Real fgw(0), w;
486 
487  int ia, jb, kc, ld;
488  for (int d = 0; d <= 1; ++d) {
489  ld = l + d;
490  for (int c = 0; c <= 1; ++c) {
491  kc = k + c;
492  for (int b = 0; b <= 1; ++b) {
493  jb = j + b;
494  for (int a = 0; a <= 1; ++a) {
495  ia = i + a;
496  w = wx[a] * wy[b] * wz[c] * wt[d];
497  if (this->Input()->IsInsideForeground(ia, jb, kc, ld)) {
498  val += w * voxel_cast<RealType>(this->Input()->Get(ia, jb, kc, ld));
499  fgw += w;
500  }
501  }
502  }
503  }
504  }
505 
506  if (ifloor(fgw * 1.e3) > 500) val /= fgw;
507  else val = voxel_cast<RealType>(this->DefaultValue());
508 
509  return voxel_cast<VoxelType>(val);
510 }
511 
512 // -----------------------------------------------------------------------------
513 template <class TImage> template <class TOtherImage>
514 inline typename TOtherImage::VoxelType
515 GenericLinearInterpolateImageFunction<TImage>
516 ::Get4D(const TOtherImage *input, double x, double y, double z, double t) const
517 {
518  typedef typename TOtherImage::VoxelType VoxelType;
519  typedef typename TOtherImage::RealType RealType;
520 
521  Real wx[2], wy[2], wz[2], wt[2];
522  const int i = ComputeWeights(x, wx), I = i + 1;
523  const int j = ComputeWeights(y, wy), J = j + 1;
524  const int k = ComputeWeights(z, wz), K = k + 1;
525  const int l = ComputeWeights(t, wt), L = l + 1;
526 
527  const auto v0000 = voxel_cast<RealType>(input->Get(i, j, k, l));
528  const auto v0001 = voxel_cast<RealType>(input->Get(I, j, k, l));
529  const auto v0010 = voxel_cast<RealType>(input->Get(i, J, k, l));
530  const auto v0011 = voxel_cast<RealType>(input->Get(I, J, k, l));
531  const auto v0100 = voxel_cast<RealType>(input->Get(i, j, K, l));
532  const auto v0101 = voxel_cast<RealType>(input->Get(I, j, K, l));
533  const auto v0110 = voxel_cast<RealType>(input->Get(i, J, K, l));
534  const auto v0111 = voxel_cast<RealType>(input->Get(I, J, K, l));
535  const auto v1000 = voxel_cast<RealType>(input->Get(i, j, k, L));
536  const auto v1001 = voxel_cast<RealType>(input->Get(I, j, k, L));
537  const auto v1010 = voxel_cast<RealType>(input->Get(i, J, k, L));
538  const auto v1011 = voxel_cast<RealType>(input->Get(I, J, k, L));
539  const auto v1100 = voxel_cast<RealType>(input->Get(i, j, K, L));
540  const auto v1101 = voxel_cast<RealType>(input->Get(I, j, K, L));
541  const auto v1110 = voxel_cast<RealType>(input->Get(i, J, K, L));
542  const auto v1111 = voxel_cast<RealType>(input->Get(I, J, K, L));
543 
544  return voxel_cast<VoxelType>(wt[0] * (wz[0] * (wy[0] * (wx[0] * v0000 + wx[1] * v0001) +
545  wy[1] * (wx[0] * v0010 + wx[1] * v0011)) +
546  wz[1] * (wy[0] * (wx[0] * v0100 + wx[1] * v0101) +
547  wy[1] * (wx[0] * v0110 + wx[1] * v0111))) +
548  wt[1] * (wz[0] * (wy[0] * (wx[0] * v1000 + wx[1] * v1001) +
549  wy[1] * (wx[0] * v1010 + wx[1] * v1011)) +
550  wz[1] * (wy[0] * (wx[0] * v1100 + wx[1] * v1101) +
551  wy[1] * (wx[0] * v1110 + wx[1] * v1111))));
552 }
553 
554 // -----------------------------------------------------------------------------
555 template <class TImage> template <class TOtherImage>
556 inline typename TOtherImage::VoxelType
557 GenericLinearInterpolateImageFunction<TImage>
558 ::GetWithPadding4D(const TOtherImage *input, double x, double y, double z, double t) const
559 {
560  typedef typename TOtherImage::VoxelType VoxelType;
561  typedef typename TOtherImage::RealType RealType;
562 
563  Real wx[2], wy[2], wz[2], wt[2];
564  const int i = ComputeWeights(x, wx);
565  const int j = ComputeWeights(y, wy);
566  const int k = ComputeWeights(z, wz);
567  const int l = ComputeWeights(t, wt);
568 
569  RealType val = voxel_cast<RealType>(0);
570  Real fgw(0), wzt, wyzt, w;
571 
572  int ia, jb, kc, ld;
573  for (int d = 0; d <= 1; ++d) {
574  ld = l + d;
575  for (int c = 0; c <= 1; ++c) {
576  kc = k + c;
577  wzt = wz[c] * wt[d];
578  for (int b = 0; b <= 1; ++b) {
579  jb = j + b;
580  wyzt = wy[b] * wzt;
581  for (int a = 0; a <= 1; ++a) {
582  ia = i + a;
583  w = wx[a] * wyzt;
584  if (input->IsForeground(ia, jb, kc, ld)) {
585  val += w * voxel_cast<RealType>(input->Get(ia, jb, kc, ld));
586  fgw += w;
587  }
588  }
589  }
590  }
591  }
592 
593  if (ifloor(fgw * 1.e3) > 500) val /= fgw;
594  else val = voxel_cast<RealType>(this->DefaultValue());
595 
596  return voxel_cast<VoxelType>(val);
597 }
598 
599 // -----------------------------------------------------------------------------
600 template <class TImage>
601 inline typename GenericLinearInterpolateImageFunction<TImage>::VoxelType
602 GenericLinearInterpolateImageFunction<TImage>
603 ::Get(double x, double y, double z, double t) const
604 {
605  switch (this->NumberOfDimensions()) {
606  case 3: return Get3D(x, y, z, t);
607  case 2: return Get2D(x, y, z, t);
608  default: return Get4D(x, y, z, t);
609  }
610 }
611 
612 // -----------------------------------------------------------------------------
613 template <class TImage>
614 inline typename GenericLinearInterpolateImageFunction<TImage>::VoxelType
615 GenericLinearInterpolateImageFunction<TImage>
616 ::GetWithPadding(double x, double y, double z, double t) const
617 {
618  switch (this->NumberOfDimensions()) {
619  case 3: return GetWithPadding3D(x, y, z, t);
620  case 2: return GetWithPadding2D(x, y, z, t);
621  default: return GetWithPadding4D(x, y, z, t);
622  }
623 }
624 
625 // -----------------------------------------------------------------------------
626 template <class TImage> template <class TOtherImage>
627 inline typename TOtherImage::VoxelType
628 GenericLinearInterpolateImageFunction<TImage>
629 ::Get(const TOtherImage *input, double x, double y, double z, double t) const
630 {
631  switch (this->NumberOfDimensions()) {
632  case 3: return Get3D(input, x, y, z, t);
633  case 2: return Get2D(input, x, y, z, t);
634  default: return Get4D(input, x, y, z, t);
635  }
636 }
637 
638 // -----------------------------------------------------------------------------
639 template <class TImage> template <class TOtherImage>
640 inline typename TOtherImage::VoxelType
641 GenericLinearInterpolateImageFunction<TImage>
642 ::GetWithPadding(const TOtherImage *input, double x, double y, double z, double t) const
643 {
644  switch (this->NumberOfDimensions()) {
645  case 3: return GetWithPadding3D(input, x, y, z, t);
646  case 2: return GetWithPadding2D(input, x, y, z, t);
647  default: return GetWithPadding4D(input, x, y, z, t);
648  }
649 }
650 
651 // -----------------------------------------------------------------------------
652 template <class TImage>
653 inline typename GenericLinearInterpolateImageFunction<TImage>::VoxelType
654 GenericLinearInterpolateImageFunction<TImage>
655 ::GetInside2D(double x, double y, double z, double t) const
656 {
657  Real wx[2], wy[2];
658  const int i = ComputeWeights(x, wx);
659  const int j = ComputeWeights(y, wy);
660  const int k = iround(z);
661  const int l = iround(t);
662 
663  auto img = reinterpret_cast<const VoxelType *>(this->Input()->GetDataPointer(i, j, k, l));
664 
665  const auto v00 = voxel_cast<RealType>(img[_Offset[0]]);
666  const auto v01 = voxel_cast<RealType>(img[_Offset[1]]);
667  const auto v10 = voxel_cast<RealType>(img[_Offset[2]]);
668  const auto v11 = voxel_cast<RealType>(img[_Offset[3]]);
669 
670  return voxel_cast<VoxelType>(wy[0] * (wx[0] * v00 + wx[1] * v01) +
671  wy[1] * (wx[0] * v10 + wx[1] * v11));
672 }
673 
674 // -----------------------------------------------------------------------------
675 template <class TImage>
676 inline typename GenericLinearInterpolateImageFunction<TImage>::VoxelType
677 GenericLinearInterpolateImageFunction<TImage>
678 ::GetInside3D(double x, double y, double z, double t) const
679 {
680  Real wx[2], wy[2], wz[2];
681  const int i = ComputeWeights(x, wx);
682  const int j = ComputeWeights(y, wy);
683  const int k = ComputeWeights(z, wz);
684  const int l = iround(t);
685 
686  auto img = reinterpret_cast<const VoxelType *>(this->Input()->GetDataPointer(i, j, k, l));
687 
688  const auto v000 = voxel_cast<RealType>(img[_Offset[0]]);
689  const auto v001 = voxel_cast<RealType>(img[_Offset[1]]);
690  const auto v010 = voxel_cast<RealType>(img[_Offset[2]]);
691  const auto v011 = voxel_cast<RealType>(img[_Offset[3]]);
692  const auto v100 = voxel_cast<RealType>(img[_Offset[4]]);
693  const auto v101 = voxel_cast<RealType>(img[_Offset[5]]);
694  const auto v110 = voxel_cast<RealType>(img[_Offset[6]]);
695  const auto v111 = voxel_cast<RealType>(img[_Offset[7]]);
696 
697  return voxel_cast<VoxelType>(wz[0] * (wy[0] * (wx[0] * v000 + wx[1] * v001) +
698  wy[1] * (wx[0] * v010 + wx[1] * v011)) +
699  wz[1] * (wy[0] * (wx[0] * v100 + wx[1] * v101) +
700  wy[1] * (wx[0] * v110 + wx[1] * v111)));
701 }
702 
703 // -----------------------------------------------------------------------------
704 template <class TImage>
705 inline typename GenericLinearInterpolateImageFunction<TImage>::VoxelType
706 GenericLinearInterpolateImageFunction<TImage>
707 ::GetInside4D(double x, double y, double z, double t) const
708 {
709  Real wx[2], wy[2], wz[2], wt[2];
710  const int i = ComputeWeights(x, wx);
711  const int j = ComputeWeights(y, wy);
712  const int k = ComputeWeights(z, wz);
713  const int l = ComputeWeights(t, wt);
714 
715  auto img = reinterpret_cast<const VoxelType *>(this->Input()->GetDataPointer(i, j, k, l));
716 
717  const auto v0000 = voxel_cast<RealType>(img[_Offset[ 0]]);
718  const auto v0001 = voxel_cast<RealType>(img[_Offset[ 1]]);
719  const auto v0010 = voxel_cast<RealType>(img[_Offset[ 2]]);
720  const auto v0011 = voxel_cast<RealType>(img[_Offset[ 3]]);
721  const auto v0100 = voxel_cast<RealType>(img[_Offset[ 4]]);
722  const auto v0101 = voxel_cast<RealType>(img[_Offset[ 5]]);
723  const auto v0110 = voxel_cast<RealType>(img[_Offset[ 6]]);
724  const auto v0111 = voxel_cast<RealType>(img[_Offset[ 7]]);
725  const auto v1000 = voxel_cast<RealType>(img[_Offset[ 8]]);
726  const auto v1001 = voxel_cast<RealType>(img[_Offset[ 9]]);
727  const auto v1010 = voxel_cast<RealType>(img[_Offset[10]]);
728  const auto v1011 = voxel_cast<RealType>(img[_Offset[11]]);
729  const auto v1100 = voxel_cast<RealType>(img[_Offset[12]]);
730  const auto v1101 = voxel_cast<RealType>(img[_Offset[13]]);
731  const auto v1110 = voxel_cast<RealType>(img[_Offset[14]]);
732  const auto v1111 = voxel_cast<RealType>(img[_Offset[15]]);
733 
734  return voxel_cast<VoxelType>(wt[0] * (wz[0] * (wy[0] * (wx[0] * v0000 + wx[1] * v0001) +
735  wy[1] * (wx[0] * v0010 + wx[1] * v0011)) +
736  wz[1] * (wy[0] * (wx[0] * v0100 + wx[1] * v0101) +
737  wy[1] * (wx[0] * v0110 + wx[1] * v0111))) +
738  wt[1] * (wz[0] * (wy[0] * (wx[0] * v1000 + wx[1] * v1001) +
739  wy[1] * (wx[0] * v1010 + wx[1] * v1011)) +
740  wz[1] * (wy[0] * (wx[0] * v1100 + wx[1] * v1101) +
741  wy[1] * (wx[0] * v1110 + wx[1] * v1111))));
742 }
743 
744 // -----------------------------------------------------------------------------
745 template <class TImage>
746 inline typename GenericLinearInterpolateImageFunction<TImage>::VoxelType
747 GenericLinearInterpolateImageFunction<TImage>
748 ::GetInside(double x, double y, double z, double t) const
749 {
750  // Use faster pixel access than Get(Input(), x, y, z, t), which requires
751  // a 4D array lookup with three indirections instead of a direct 1D lookup
752  switch (this->NumberOfDimensions()) {
753  case 3: return GetInside3D(x, y, z, t);
754  case 2: return GetInside2D(x, y, z, t);
755  default: return GetInside4D(x, y, z, t);
756  }
757 }
758 
759 // -----------------------------------------------------------------------------
760 template <>
761 inline typename GenericLinearInterpolateImageFunction<BaseImage>::VoxelType
762 GenericLinearInterpolateImageFunction<BaseImage>
763 ::GetInside(double x, double y, double z, double t) const
764 {
765  return Get(this->Input(), x, y, z, t);
766 }
767 
768 // -----------------------------------------------------------------------------
769 template <class TImage>
770 inline typename GenericLinearInterpolateImageFunction<TImage>::VoxelType
771 GenericLinearInterpolateImageFunction<TImage>
772 ::GetOutside(double x, double y, double z, double t) const
773 {
774  if (this->Extrapolator()) {
775  return Get(this->Extrapolator(), x, y, z, t);
776  } else {
777  return Get(x, y, z, t);
778  }
779 }
780 
781 // -----------------------------------------------------------------------------
782 template <class TImage>
783 inline typename GenericLinearInterpolateImageFunction<TImage>::VoxelType
784 GenericLinearInterpolateImageFunction<TImage>
785 ::GetWithPaddingInside(double x, double y, double z, double t) const
786 {
787  return GetWithPadding(this->Input(), x, y, z, t);
788 }
789 
790 // -----------------------------------------------------------------------------
791 template <class TImage>
792 inline typename GenericLinearInterpolateImageFunction<TImage>::VoxelType
793 GenericLinearInterpolateImageFunction<TImage>
794 ::GetWithPaddingOutside(double x, double y, double z, double t) const
795 {
796  if (this->Extrapolator()) {
797  return GetWithPadding(this->Extrapolator(), x, y, z, t);
798  } else {
799  return GetWithPadding(x, y, z, t);
800  }
801 }
802 
803 // =============================================================================
804 // Evaluation of derivatives
805 // =============================================================================
806 
807 // -----------------------------------------------------------------------------
808 template <class TImage>
809 void GenericLinearInterpolateImageFunction<TImage>
810 ::EvaluateJacobianInside(Matrix &jac, double x, double y, double z, double t) const
811 {
812  Jacobian(jac, this->Input(), x, y, z, t);
813 }
814 
815 // -----------------------------------------------------------------------------
816 template <class TImage>
817 void GenericLinearInterpolateImageFunction<TImage>
818 ::EvaluateJacobianOutside(Matrix &jac, double x, double y, double z, double t) const
819 {
820  if (this->Extrapolator()) {
821  Jacobian(jac, this->Extrapolator(), x, y, z, t);
822  } else {
823  Jacobian(jac, x, y, z, t);
824  }
825 }
826 
827 // -----------------------------------------------------------------------------
828 template <class TImage>
829 void GenericLinearInterpolateImageFunction<TImage>
830 ::EvaluateJacobianWithPaddingInside(Matrix &jac, double x, double y, double z, double t) const
831 {
832  JacobianWithPadding(jac, x, y, z, t);
833 }
834 
835 // -----------------------------------------------------------------------------
836 template <class TImage>
837 void GenericLinearInterpolateImageFunction<TImage>
838 ::EvaluateJacobianWithPaddingOutside(Matrix &jac, double x, double y, double z, double t) const
839 {
840  JacobianWithPadding(jac, x, y, z, t);
841 }
842 
843 // -----------------------------------------------------------------------------
844 // without use of extrapolation
845 // -----------------------------------------------------------------------------
846 
847 // -----------------------------------------------------------------------------
848 template <class TImage>
849 void GenericLinearInterpolateImageFunction<TImage>
850 ::Jacobian2D(Matrix &jac, double x, double y, double z, double t) const
851 {
852  const TImage * const input = this->Input();
853 
854  Real wx[2], wy[2], wd[2] = {Real(-1), Real(1)};
855  const int i = ComputeWeights(x, wx);
856  const int j = ComputeWeights(y, wy);
857  const int k = iround(z);
858 
859  if (IsNaN(t) && IsZero(input->TSize()) && input->N() == 1) {
860 
861  jac.Initialize(input->T(), 2);
862 
863  if (i < 0 || i >= input->X() - 1 ||
864  j < 0 || j >= input->Y() - 1 ||
865  k < 0 || k >= input->Z()) {
866  return; // zero
867  }
868 
869  Real dv, w[2];
870  int ia, jb;
871  for (int b = 0; b <= 1; ++b) {
872  jb = j + b;
873  for (int a = 0; a <= 1; ++a) {
874  ia = i + a;
875  w[0] = wd[a] * wy[b];
876  w[1] = wx[a] * wd[b];
877  for (int l = 0; l < input->T(); ++l) {
878  dv = voxel_cast<Real>(input->Get(ia, jb, k, l));
879  jac(l, 0) += w[0] * dv;
880  jac(l, 1) += w[1] * dv;
881  }
882  }
883  }
884 
885  } else {
886 
887  jac.Initialize(input->N(), 2);
888 
889  const int l = (IsNaN(t) ? 0 : iround(t));
890  if (i < 0 || i >= input->X() - 1 ||
891  j < 0 || j >= input->Y() - 1 ||
892  k < 0 || k >= input->Z() ||
893  l < 0 || l >= input->T()) {
894  return; // zero
895  }
896 
897  RealType dv; // input value at discrete point
898  RealType dx = voxel_cast<RealType>(0); // derivative(s) w.r.t. x
899  RealType dy = voxel_cast<RealType>(0); // derivative(s) w.r.t. y
900 
901  int ia, jb;
902  for (int b = 0; b <= 1; ++b) {
903  jb = j + b;
904  for (int a = 0; a <= 1; ++a) {
905  ia = i + a;
906  dv = voxel_cast<RealType>(input->Get(ia, jb, k, l));
907  dx += (wd[a] * wy[b]) * dv;
908  dy += (wx[a] * wd[b]) * dv;
909  }
910  }
911 
912  for (int r = 0; r < input->N(); ++r) {
913  jac(r, 0) = get(dx, r);
914  jac(r, 1) = get(dy, r);
915  }
916 
917  }
918 }
919 
920 // -----------------------------------------------------------------------------
921 template <class TImage>
922 void GenericLinearInterpolateImageFunction<TImage>
923 ::Jacobian3D(Matrix &jac, double x, double y, double z, double t) const
924 {
925  const TImage * const input = this->Input();
926 
927  Real wx[2], wy[2], wz[2], wd[2] = {Real(-1), Real(1)};
928  const int i = ComputeWeights(x, wx);
929  const int j = ComputeWeights(y, wy);
930  const int k = ComputeWeights(z, wz);
931 
932  if (IsNaN(t) && IsZero(input->TSize()) && input->N() == 1) {
933 
934  jac.Initialize(input->T(), 3);
935 
936  if (i < 0 || i >= input->X() - 1 ||
937  j < 0 || j >= input->Y() - 1 ||
938  k < 0 || k >= input->Z() - 1) {
939  return; // zero
940  }
941 
942  Real w[3], dv;
943  int ia, jb, kc;
944  for (int c = 0; c <= 1; ++c) {
945  kc = k + c;
946  for (int b = 0; b <= 1; ++b) {
947  jb = j + b;
948  for (int a = 0; a <= 1; ++a) {
949  ia = i + a;
950  w[0] = wd[a] * wy[b] * wz[c];
951  w[1] = wx[a] * wd[b] * wz[c];
952  w[2] = wx[a] * wy[b] * wd[c];
953  for (int l = 0; l < input->T(); ++l) {
954  dv = voxel_cast<Real>(input->Get(ia, jb, kc, l));
955  jac(l, 0) += w[0] * dv;
956  jac(l, 1) += w[1] * dv;
957  jac(l, 2) += w[2] * dv;
958  }
959  }
960  }
961  }
962 
963  } else {
964 
965  jac.Initialize(input->N(), 3);
966 
967  const int l = iround(t);
968  if (i < 0 || i >= input->X() - 1 ||
969  j < 0 || j >= input->Y() - 1 ||
970  k < 0 || k >= input->Z() - 1 ||
971  l < 0 || l >= input->T()) {
972  return; // zero
973  }
974 
975  RealType dv; // input value at discrete point
976  RealType dx = voxel_cast<RealType>(0); // derivative(s) w.r.t. x
977  RealType dy = voxel_cast<RealType>(0); // derivative(s) w.r.t. y
978  RealType dz = voxel_cast<RealType>(0); // derivative(s) w.r.t. z
979 
980  int ia, jb, kc;
981  for (int c = 0; c <= 1; ++c) {
982  kc = k + c;
983  for (int b = 0; b <= 1; ++b) {
984  jb = j + b;
985  for (int a = 0; a <= 1; ++a) {
986  ia = i + a;
987  dv = voxel_cast<RealType>(input->Get(ia, jb, kc, l));
988  dx += (wd[a] * wy[b] * wz[c]) * dv;
989  dy += (wx[a] * wd[b] * wz[c]) * dv;
990  dz += (wx[a] * wy[b] * wd[c]) * dv;
991  }
992  }
993  }
994 
995  for (int r = 0; r < input->N(); ++r) {
996  jac(r, 0) = get(dx, r);
997  jac(r, 1) = get(dy, r);
998  jac(r, 2) = get(dz, r);
999  }
1000 
1001  }
1002 }
1003 
1004 // -----------------------------------------------------------------------------
1005 template <class TImage>
1006 void GenericLinearInterpolateImageFunction<TImage>
1007 ::Jacobian4D(Matrix &jac, double x, double y, double z, double t) const
1008 {
1009  const TImage * const input = this->Input();
1010 
1011  jac.Initialize(input->N(), 4);
1012 
1013  Real wx[2], wy[2], wz[2], wt[2], wd[2] = {Real(-1), Real(1)};
1014  const int i = ComputeWeights(x, wx);
1015  const int j = ComputeWeights(y, wy);
1016  const int k = ComputeWeights(z, wz);
1017  const int l = ComputeWeights(IsNaN(t) ? 0. : t, wt);
1018 
1019  if (i < 0 || i >= input->X() - 1 ||
1020  j < 0 || j >= input->Y() - 1 ||
1021  k < 0 || k >= input->Z() - 1 ||
1022  l < 0 || l >= input->T() - 1) {
1023  return; // zero
1024  }
1025 
1026  RealType dv; // input value at discrete point
1027  RealType dx = voxel_cast<RealType>(0); // derivative(s) w.r.t. x
1028  RealType dy = voxel_cast<RealType>(0); // derivative(s) w.r.t. y
1029  RealType dz = voxel_cast<RealType>(0); // derivative(s) w.r.t. z
1030  RealType dt = voxel_cast<RealType>(0); // derivative(s) w.r.t. t
1031 
1032  int ia, jb, kc, ld;
1033  for (int d = 0; d <= 1; ++d) {
1034  ld = l + d;
1035  for (int c = 0; c <= 1; ++c) {
1036  kc = k + c;
1037  for (int b = 0; b <= 1; ++b) {
1038  jb = j + b;
1039  for (int a = 0; a <= 1; ++a) {
1040  ia = i + a;
1041  dv = voxel_cast<RealType>(input->Get(ia, jb, kc, ld));
1042  dx += (wd[a] * wy[b] * wz[c] * wt[d]) * dv;
1043  dy += (wx[a] * wd[b] * wz[c] * wt[d]) * dv;
1044  dz += (wx[a] * wy[b] * wd[c] * wt[d]) * dv;
1045  dt += (wx[a] * wy[b] * wz[c] * wd[d]) * dv;
1046  }
1047  }
1048  }
1049  }
1050 
1051  for (int r = 0; r < input->N(); ++r) {
1052  jac(r, 0) = get(dx, r);
1053  jac(r, 1) = get(dy, r);
1054  jac(r, 2) = get(dz, r);
1055  jac(r, 3) = get(dt, r);
1056  }
1057 }
1058 
1059 // -----------------------------------------------------------------------------
1060 template <class TImage>
1061 void GenericLinearInterpolateImageFunction<TImage>
1062 ::Jacobian(Matrix &jac, double x, double y, double z, double t) const
1063 {
1064  switch (this->NumberOfDimensions()) {
1065  case 3: return Jacobian3D(jac, x, y, z, t);
1066  case 2: return Jacobian2D(jac, x, y, z, t);
1067  default: return Jacobian4D(jac, x, y, z, t);
1068  }
1069 }
1070 
1071 
1072 // -----------------------------------------------------------------------------
1073 // without use of extrapolation, inside image foreground
1074 // -----------------------------------------------------------------------------
1075 
1076 // -----------------------------------------------------------------------------
1077 template <class TImage>
1078 void GenericLinearInterpolateImageFunction<TImage>
1079 ::JacobianWithPadding2D(Matrix &jac, double x, double y, double z, double t) const
1080 {
1081  const TImage * const input = this->Input();
1082 
1083  Real wx[2], wy[2], wd[2] = {Real(-1), Real(1)};
1084  const int i = ComputeWeights(x, wx);
1085  const int j = ComputeWeights(y, wy);
1086  const int k = iround(z);
1087 
1088  if (IsNaN(t) && IsZero(input->TSize()) && input->N() == 1) {
1089 
1090  jac.Initialize(input->T(), 2);
1091 
1092  if (i < 0 || i >= input->X() - 1 ||
1093  j < 0 || j >= input->Y() - 1 ||
1094  k < 0 || k >= input->Z()) {
1095  return; // zero
1096  }
1097 
1098  Real w[2], dv;
1099  int ia, jb;
1100  for (int b = 0; b <= 1; ++b) {
1101  jb = j + b;
1102  for (int a = 0; a <= 1; ++a) {
1103  ia = i + a;
1104  if (!input->IsForeground(ia, jb, k)) {
1105  jac.Initialize(input->T(), 2);
1106  return; // zero
1107  }
1108  w[0] = wd[a] * wy[b];
1109  w[1] = wx[a] * wd[b];
1110  for (int l = 0; l < input->T(); ++l) {
1111  dv = voxel_cast<Real>(input->Get(ia, jb, k, l));
1112  jac(l, 0) += w[0] * dv;
1113  jac(l, 1) += w[1] * dv;
1114  }
1115  }
1116  }
1117 
1118  } else {
1119 
1120  jac.Initialize(input->N(), 2);
1121 
1122  const int l = (IsNaN(t) ? 0 : iround(t));
1123  if (i < 0 || i >= input->X() - 1 ||
1124  j < 0 || j >= input->Y() - 1 ||
1125  k < 0 || k >= input->Z() ||
1126  l < 0 || l >= input->T()) {
1127  return; // zero
1128  }
1129 
1130  RealType dv; // input value at discrete point
1131  RealType dx = voxel_cast<RealType>(0); // derivative(s) w.r.t. x
1132  RealType dy = voxel_cast<RealType>(0); // derivative(s) w.r.t. y
1133 
1134  int ia, jb;
1135  for (int b = 0; b <= 1; ++b) {
1136  jb = j + b;
1137  for (int a = 0; a <= 1; ++a) {
1138  ia = i + a;
1139  if (!input->IsForeground(ia, jb, k, l)) {
1140  return; // zero
1141  }
1142  dv = voxel_cast<RealType>(input->Get(ia, jb, k, l));
1143  dx += (wd[a] * wy[b]) * dv;
1144  dy += (wx[a] * wd[b]) * dv;
1145  }
1146  }
1147 
1148  for (int r = 0; r < input->N(); ++r) {
1149  jac(r, 0) = get(dx, r);
1150  jac(r, 1) = get(dy, r);
1151  }
1152 
1153  }
1154 }
1155 
1156 // -----------------------------------------------------------------------------
1157 template <class TImage>
1158 void GenericLinearInterpolateImageFunction<TImage>
1159 ::JacobianWithPadding3D(Matrix &jac, double x, double y, double z, double t) const
1160 {
1161  const TImage * const input = this->Input();
1162 
1163  Real wx[2], wy[2], wz[2], wd[2] = {Real(-1), Real(1)};
1164  const int i = ComputeWeights(x, wx);
1165  const int j = ComputeWeights(y, wy);
1166  const int k = ComputeWeights(z, wz);
1167 
1168  if (IsNaN(t) && IsZero(input->TSize()) && input->N() == 1) {
1169 
1170  jac.Initialize(input->T(), 3);
1171 
1172  if (i < 0 || i >= input->X() - 1 ||
1173  j < 0 || j >= input->Y() - 1 ||
1174  k < 0 || k >= input->Z() - 1) {
1175  return; // zero
1176  }
1177 
1178  Real w[3], dv;
1179  int ia, jb, kc;
1180  for (int c = 0; c <= 1; ++c) {
1181  kc = k + c;
1182  for (int b = 0; b <= 1; ++b) {
1183  jb = j + b;
1184  for (int a = 0; a <= 1; ++a) {
1185  ia = i + a;
1186  if (!input->IsForeground(ia, jb, kc)) {
1187  jac.Initialize(input->T(), 3);
1188  return; // zero
1189  }
1190  w[0] = wd[a] * wy[b] * wz[c];
1191  w[1] = wx[a] * wd[b] * wz[c];
1192  w[2] = wx[a] * wy[b] * wd[c];
1193  for (int l = 0; l < input->T(); ++l) {
1194  dv = voxel_cast<Real>(input->Get(ia, jb, kc, l));
1195  jac(l, 0) += w[0] * dv;
1196  jac(l, 1) += w[1] * dv;
1197  jac(l, 2) += w[2] * dv;
1198  }
1199  }
1200  }
1201  }
1202 
1203  } else {
1204 
1205  jac.Initialize(input->N(), 3);
1206 
1207  const int l = (IsNaN(t) ? 0 : iround(t));
1208  if (i < 0 || i >= input->X() - 1 ||
1209  j < 0 || j >= input->Y() - 1 ||
1210  k < 0 || k >= input->Z() - 1 ||
1211  l < 0 || l >= input->T()) {
1212  return; // zero
1213  }
1214 
1215  RealType dv; // input value at discrete point
1216  RealType dx = voxel_cast<RealType>(0); // derivative(s) w.r.t. x
1217  RealType dy = voxel_cast<RealType>(0); // derivative(s) w.r.t. y
1218  RealType dz = voxel_cast<RealType>(0); // derivative(s) w.r.t. z
1219 
1220  int ia, jb, kc;
1221  for (int c = 0; c <= 1; ++c) {
1222  kc = k + c;
1223  for (int b = 0; b <= 1; ++b) {
1224  jb = j + b;
1225  for (int a = 0; a <= 1; ++a) {
1226  ia = i + a;
1227  if (!input->IsForeground(ia, jb, kc, l)) {
1228  return; // zero
1229  }
1230  dv = voxel_cast<RealType>(input->Get(ia, jb, kc, l));
1231  dx += (wd[a] * wy[b] * wz[c]) * dv;
1232  dy += (wx[a] * wd[b] * wz[c]) * dv;
1233  dz += (wx[a] * wy[b] * wd[c]) * dv;
1234  }
1235  }
1236  }
1237 
1238  for (int r = 0; r < input->N(); ++r) {
1239  jac(r, 0) = get(dx, r);
1240  jac(r, 1) = get(dy, r);
1241  jac(r, 2) = get(dz, r);
1242  }
1243 
1244  }
1245 }
1246 
1247 // -----------------------------------------------------------------------------
1248 template <class TImage>
1249 void GenericLinearInterpolateImageFunction<TImage>
1250 ::JacobianWithPadding4D(Matrix &jac, double x, double y, double z, double t) const
1251 {
1252  const TImage * const input = this->Input();
1253 
1254  jac.Initialize(input->N(), 4);
1255 
1256  Real wx[2], wy[2], wz[2], wt[2], wd[2] = {Real(-1), Real(1)};
1257  const int i = ComputeWeights(x, wx);
1258  const int j = ComputeWeights(y, wy);
1259  const int k = ComputeWeights(z, wz);
1260  const int l = ComputeWeights(IsNaN(t) ? 0. : t, wt);
1261 
1262  if (i < 0 || i >= input->X() - 1 ||
1263  j < 0 || j >= input->Y() - 1 ||
1264  k < 0 || k >= input->Z() - 1 ||
1265  l < 0 || l >= input->T() - 1) {
1266  return; // zero
1267  }
1268 
1269  RealType dv; // input value at discrete point
1270  RealType dx = voxel_cast<RealType>(0); // derivative(s) w.r.t. x
1271  RealType dy = voxel_cast<RealType>(0); // derivative(s) w.r.t. y
1272  RealType dz = voxel_cast<RealType>(0); // derivative(s) w.r.t. z
1273  RealType dt = voxel_cast<RealType>(0); // derivative(s) w.r.t. t
1274 
1275  int ia, jb, kc, ld;
1276  for (int d = 0; d <= 1; ++d) {
1277  ld = l + d;
1278  for (int c = 0; c <= 1; ++c) {
1279  kc = k + c;
1280  for (int b = 0; b <= 1; ++b) {
1281  jb = j + b;
1282  for (int a = 0; a <= 1; ++a) {
1283  ia = i + a;
1284  if (!input->IsForeground(ia, jb, kc, ld)) {
1285  return; // zero
1286  }
1287  dv = voxel_cast<RealType>(input->Get(ia, jb, kc, ld));
1288  dx += (wd[a] * wy[b] * wz[c] * wt[d]) * dv;
1289  dy += (wx[a] * wd[b] * wz[c] * wt[d]) * dv;
1290  dz += (wx[a] * wy[b] * wd[c] * wt[d]) * dv;
1291  dt += (wx[a] * wy[b] * wz[c] * wd[d]) * dv;
1292  }
1293  }
1294  }
1295  }
1296 
1297  for (int r = 0; r < input->N(); ++r) {
1298  jac(r, 0) = get(dx, r);
1299  jac(r, 1) = get(dy, r);
1300  jac(r, 2) = get(dz, r);
1301  jac(r, 3) = get(dt, r);
1302  }
1303 }
1304 
1305 // -----------------------------------------------------------------------------
1306 template <class TImage>
1307 void GenericLinearInterpolateImageFunction<TImage>
1308 ::JacobianWithPadding(Matrix &jac, double x, double y, double z, double t) const
1309 {
1310  switch (this->NumberOfDimensions()) {
1311  case 3: return JacobianWithPadding3D(jac, x, y, z, t);
1312  case 2: return JacobianWithPadding2D(jac, x, y, z, t);
1313  default: return JacobianWithPadding4D(jac, x, y, z, t);
1314  }
1315 }
1316 
1317 
1318 // -----------------------------------------------------------------------------
1319 // inside image domain or with extrapolation function
1320 // -----------------------------------------------------------------------------
1321 
1322 // -----------------------------------------------------------------------------
1323 template <class TImage> template <class TOtherImage>
1324 void GenericLinearInterpolateImageFunction<TImage>
1325 ::Jacobian2D(Matrix &jac, const TOtherImage *input, double x, double y, double z, double t) const
1326 {
1327  Real wx[2], wy[2], wd[2] = {Real(-1), Real(1)};
1328  const int i = ComputeWeights(x, wx);
1329  const int j = ComputeWeights(y, wy);
1330  const int k = iround(z);
1331 
1332  if (IsNaN(t) && IsZero(input->TSize()) && input->N() == 1) {
1333 
1334  jac.Initialize(input->T(), 2);
1335 
1336  if (k < 0 || k >= input->Z()) {
1337  return; // zero
1338  }
1339 
1340  Real w[2], dv;
1341  int ia, jb;
1342  for (int b = 0; b <= 1; ++b) {
1343  jb = j + b;
1344  for (int a = 0; a <= 1; ++a) {
1345  ia = i + a;
1346  w[0] = wd[a] * wy[b];
1347  w[1] = wx[a] * wd[b];
1348  for (int l = 0; l < input->T(); ++l) {
1349  dv = voxel_cast<Real>(input->Get(ia, jb, k, l));
1350  jac(l, 0) += w[0] * dv;
1351  jac(l, 1) += w[1] * dv;
1352  }
1353  }
1354  }
1355 
1356  } else {
1357 
1358  jac.Initialize(input->N(), 2);
1359 
1360  const int l = (IsNaN(t) ? 0 : iround(t));
1361  if (k < 0 || k >= input->Z() ||
1362  l < 0 || l >= input->T()) {
1363  return; // zero
1364  }
1365 
1366  RealType dv; // input value at discrete point
1367  RealType dx = voxel_cast<RealType>(0); // derivative(s) w.r.t. x
1368  RealType dy = voxel_cast<RealType>(0); // derivative(s) w.r.t. y
1369 
1370  int ia, jb;
1371  for (int b = 0; b <= 1; ++b) {
1372  jb = j + b;
1373  for (int a = 0; a <= 1; ++a) {
1374  ia = i + a;
1375  dv = voxel_cast<RealType>(input->Get(ia, jb, k, l));
1376  dx += (wd[a] * wy[b]) * dv;
1377  dy += (wx[a] * wd[b]) * dv;
1378  }
1379  }
1380 
1381  for (int r = 0; r < input->N(); ++r) {
1382  jac(r, 0) = get(dx, r);
1383  jac(r, 1) = get(dy, r);
1384  }
1385 
1386  }
1387 }
1388 
1389 // -----------------------------------------------------------------------------
1390 template <class TImage> template <class TOtherImage>
1391 void GenericLinearInterpolateImageFunction<TImage>
1392 ::Jacobian3D(Matrix &jac, const TOtherImage *input,
1393  double x, double y, double z, double t) const
1394 {
1395  Real wx[2], wy[2], wz[2], wd[2] = {Real(-1), Real(1)};
1396  const int i = ComputeWeights(x, wx);
1397  const int j = ComputeWeights(y, wy);
1398  const int k = ComputeWeights(z, wz);
1399 
1400  if (IsNaN(t) && IsZero(input->TSize()) && input->N() == 1) {
1401 
1402  jac.Initialize(input->T(), 3);
1403 
1404  Real w[3], dv;
1405  int ia, jb, kc;
1406  for (int c = 0; c <= 1; ++c) {
1407  kc = k + c;
1408  for (int b = 0; b <= 1; ++b) {
1409  jb = j + b;
1410  for (int a = 0; a <= 1; ++a) {
1411  ia = i + a;
1412  w[0] = wd[a] * wy[b] * wz[c];
1413  w[1] = wx[a] * wd[b] * wz[c];
1414  w[2] = wx[a] * wy[b] * wd[c];
1415  for (int l = 0; l < input->T(); ++l) {
1416  dv = voxel_cast<Real>(input->Get(ia, jb, kc, l));
1417  jac(l, 0) += w[0] * dv;
1418  jac(l, 1) += w[1] * dv;
1419  jac(l, 2) += w[2] * dv;
1420  }
1421  }
1422  }
1423  }
1424 
1425  } else {
1426 
1427  jac.Initialize(input->N(), 3);
1428 
1429  const int l = (IsNaN(t) ? 0 : iround(t));
1430  if (l < 0 || l >= input->T()) {
1431  return; // zero
1432  }
1433 
1434  RealType dv; // input value at discrete point
1435  RealType dx = voxel_cast<RealType>(0); // derivative(s) w.r.t. x
1436  RealType dy = voxel_cast<RealType>(0); // derivative(s) w.r.t. y
1437  RealType dz = voxel_cast<RealType>(0); // derivative(s) w.r.t. z
1438 
1439  int ia, jb, kc;
1440  for (int c = 0; c <= 1; ++c) {
1441  kc = k + c;
1442  for (int b = 0; b <= 1; ++b) {
1443  jb = j + b;
1444  for (int a = 0; a <= 1; ++a) {
1445  ia = i + a;
1446  dv = voxel_cast<RealType>(input->Get(ia, jb, kc, l));
1447  dx += (wd[a] * wy[b] * wz[c]) * dv;
1448  dy += (wx[a] * wd[b] * wz[c]) * dv;
1449  dz += (wx[a] * wy[b] * wd[c]) * dv;
1450  }
1451  }
1452  }
1453 
1454  for (int r = 0; r < input->N(); ++r) {
1455  jac(r, 0) = get(dx, r);
1456  jac(r, 1) = get(dy, r);
1457  jac(r, 2) = get(dz, r);
1458  }
1459 
1460  }
1461 }
1462 
1463 // -----------------------------------------------------------------------------
1464 template <class TImage> template <class TOtherImage>
1465 void GenericLinearInterpolateImageFunction<TImage>
1466 ::Jacobian4D(Matrix &jac, const TOtherImage *input,
1467  double x, double y, double z, double t) const
1468 {
1469  jac.Initialize(input->N(), 4);
1470 
1471  Real wx[2], wy[2], wz[2], wt[2], wd[2] = {Real(-1), Real(1)};
1472  const int i = ComputeWeights(x, wx);
1473  const int j = ComputeWeights(y, wy);
1474  const int k = ComputeWeights(z, wz);
1475  const int l = ComputeWeights(IsNaN(t) ? 0. : t, wt);
1476 
1477  RealType dv; // input value at discrete point
1478  RealType dx = voxel_cast<RealType>(0); // derivative(s) w.r.t. x
1479  RealType dy = voxel_cast<RealType>(0); // derivative(s) w.r.t. y
1480  RealType dz = voxel_cast<RealType>(0); // derivative(s) w.r.t. z
1481  RealType dt = voxel_cast<RealType>(0); // derivative(s) w.r.t. t
1482 
1483  int ia, jb, kc, ld;
1484  for (int d = 0; d <= 1; ++d) {
1485  ld = l + d;
1486  for (int c = 0; c <= 1; ++c) {
1487  kc = k + c;
1488  for (int b = 0; b <= 1; ++b) {
1489  jb = j + b;
1490  for (int a = 0; a <= 1; ++a) {
1491  ia = i + a;
1492  dv = voxel_cast<RealType>(input->Get(ia, jb, kc, ld));
1493  dx += (wd[a] * wy[b] * wz[c] * wt[d]) * dv;
1494  dy += (wx[a] * wd[b] * wz[c] * wt[d]) * dv;
1495  dz += (wx[a] * wy[b] * wd[c] * wt[d]) * dv;
1496  dt += (wx[a] * wy[b] * wz[c] * wd[d]) * dv;
1497  }
1498  }
1499  }
1500  }
1501 
1502  for (int r = 0; r < input->N(); ++r) {
1503  jac(r, 0) = get(dx, r);
1504  jac(r, 1) = get(dy, r);
1505  jac(r, 2) = get(dz, r);
1506  jac(r, 3) = get(dt, r);
1507  }
1508 }
1509 
1510 // -----------------------------------------------------------------------------
1511 template <class TImage> template <class TOtherImage>
1512 void GenericLinearInterpolateImageFunction<TImage>
1513 ::Jacobian(Matrix &jac, const TOtherImage *input,
1514  double x, double y, double z, double t) const
1515 {
1516  switch (this->NumberOfDimensions()) {
1517  case 3: return Jacobian3D(jac, input, x, y, z, t);
1518  case 2: return Jacobian2D(jac, input, x, y, z, t);
1519  default: return Jacobian4D(jac, input, x, y, z, t);
1520  }
1521 }
1522 
1523 
1524 } // namespace mirtk
1525 
1526 #endif // MIRTK_LinearInterpolateImageFunction_HXX
mirtk::GenericLinearInterpolateImageFunction::GetWithPadding4D
VoxelType GetWithPadding4D(double, double, double=0, double=0) const
Definition: LinearInterpolateImageFunction.hxx:476
mirtk::GenericLinearInterpolateImageFunction::Get
VoxelType Get(double, double, double=0, double=0) const
Definition: LinearInterpolateImageFunction.hxx:603
mirtk::GenericLinearInterpolateImageFunction::JacobianWithPadding3D
void JacobianWithPadding3D(Matrix &, double, double, double=0., double=NaN) const
Definition: LinearInterpolateImageFunction.hxx:1159
mirtk::GenericLinearInterpolateImageFunction::Jacobian4D
void Jacobian4D(Matrix &, double, double, double=0., double=NaN) const
Definition: LinearInterpolateImageFunction.hxx:1007
mirtk::IsZero
MIRTKCU_API bool IsZero(double a, double tol=1e-12)
Determine equality of a floating point number with zero.
Definition: Math.h:131
mirtk::Get
string Get(const ParameterList &params, string name)
Get parameter value from parameters list.
Definition: Object.h:202
mirtk::GenericLinearInterpolateImageFunction::EvaluateJacobianWithPaddingInside
virtual void EvaluateJacobianWithPaddingInside(Matrix &, double, double, double=0, double=NaN) const
Definition: LinearInterpolateImageFunction.hxx:830
mirtk::GenericLinearInterpolateImageFunction::Get3D
VoxelType Get3D(double, double, double=0, double=0) const
Definition: LinearInterpolateImageFunction.hxx:269
mirtk::GenericLinearInterpolateImageFunction::GetInside
virtual VoxelType GetInside(double, double, double=0, double=0) const
Definition: LinearInterpolateImageFunction.hxx:748
mirtk::IsNaN
MIRTKCU_API bool IsNaN(double x)
Check if floating point value is not a number (NaN)
Definition: Math.h:91
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::GenericLinearInterpolateImageFunction::EvaluateJacobianInside
virtual void EvaluateJacobianInside(Matrix &, double, double, double=0, double=NaN) const
Definition: LinearInterpolateImageFunction.hxx:810
mirtk::fdec
MIRTKCU_API double fdec(double f)
Definition: Math.h:188
mirtk::Matrix::Initialize
void Initialize(int, int=-1, double *=NULL)
Initialize matrix with number of rows and columns.
mirtk::GenericLinearInterpolateImageFunction::EvaluateJacobianOutside
virtual void EvaluateJacobianOutside(Matrix &, double, double, double=0, double=NaN) const
Definition: LinearInterpolateImageFunction.hxx:818
mirtk
Definition: IOConfig.h:41
mirtk::GenericLinearInterpolateImageFunction::GetInside3D
virtual VoxelType GetInside3D(double, double, double=0, double=0) const
Evaluate generic 3D image without handling boundary conditions.
Definition: LinearInterpolateImageFunction.hxx:678
mirtk::GenericLinearInterpolateImageFunction::GetWithPaddingOutside
virtual VoxelType GetWithPaddingOutside(double, double, double=0, double=0) const
Definition: LinearInterpolateImageFunction.hxx:794
mirtk::GenericLinearInterpolateImageFunction::BoundingInterval
virtual void BoundingInterval(double, int &, int &) const
Definition: LinearInterpolateImageFunction.hxx:110
mirtk::GenericLinearInterpolateImageFunction::Jacobian
void Jacobian(Matrix &, double, double, double=0., double=NaN) const
Get 1st order derivatives of given image at arbitrary location (in pixels)
Definition: LinearInterpolateImageFunction.hxx:1062
mirtk::GenericLinearInterpolateImageFunction::Get4D
VoxelType Get4D(double, double, double=0, double=0) const
Definition: LinearInterpolateImageFunction.hxx:430
mirtk::GenericLinearInterpolateImageFunction::GetInside2D
virtual VoxelType GetInside2D(double, double, double=0, double=0) const
Evaluate generic 2D image without handling boundary conditions.
Definition: LinearInterpolateImageFunction.hxx:655
mirtk::GenericLinearInterpolateImageFunction::Jacobian3D
void Jacobian3D(Matrix &, double, double, double=0., double=NaN) const
Definition: LinearInterpolateImageFunction.hxx:923
mirtk::GenericLinearInterpolateImageFunction::Get2D
VoxelType Get2D(double, double, double=0, double=0) const
Definition: LinearInterpolateImageFunction.hxx:123
mirtk::GenericLinearInterpolateImageFunction::GetWithPadding3D
VoxelType GetWithPadding3D(double, double, double=0, double=0) const
Definition: LinearInterpolateImageFunction.hxx:314
mirtk::GenericLinearInterpolateImageFunction::GetWithPaddingInside
virtual VoxelType GetWithPaddingInside(double, double, double=0, double=0) const
Definition: LinearInterpolateImageFunction.hxx:785
mirtk::GenericLinearInterpolateImageFunction::JacobianWithPadding4D
void JacobianWithPadding4D(Matrix &, double, double, double=0., double=NaN) const
Definition: LinearInterpolateImageFunction.hxx:1250
mirtk::GenericLinearInterpolateImageFunction::Jacobian2D
void Jacobian2D(Matrix &, double, double, double=0., double=NaN) const
Definition: LinearInterpolateImageFunction.hxx:850
mirtk::GenericLinearInterpolateImageFunction::GetWithPadding
virtual VoxelType GetWithPadding(double, double, double=0, double=0) const
Definition: LinearInterpolateImageFunction.hxx:616
mirtk::iround
MIRTKCU_API int iround(T x)
Round floating-point value and cast to int.
Definition: Math.h:170
mirtk::GenericLinearInterpolateImageFunction::EvaluateJacobianWithPaddingOutside
virtual void EvaluateJacobianWithPaddingOutside(Matrix &, double, double, double=0, double=NaN) const
Definition: LinearInterpolateImageFunction.hxx:838
mirtk::GenericLinearInterpolateImageFunction::GetInside4D
virtual VoxelType GetInside4D(double, double, double=0, double=0) const
Evaluate generic 4D image without handling boundary conditions.
Definition: LinearInterpolateImageFunction.hxx:707
mirtk::GenericLinearInterpolateImageFunction::JacobianWithPadding
void JacobianWithPadding(Matrix &, double, double, double=0, double=NaN) const
Get 1st order derivatives of given image at arbitrary location (in pixels)
Definition: LinearInterpolateImageFunction.hxx:1308
mirtk::GenericLinearInterpolateImageFunction::GetOutside
virtual VoxelType GetOutside(double, double, double=0, double=0) const
Evaluate generic image at an arbitrary location (in pixels)
Definition: LinearInterpolateImageFunction.hxx:772
mirtk::GenericLinearInterpolateImageFunction::GetWithPadding2D
VoxelType GetWithPadding2D(double, double, double=0, double=0) const
Definition: LinearInterpolateImageFunction.hxx:164
mirtk::GenericLinearInterpolateImageFunction::JacobianWithPadding2D
void JacobianWithPadding2D(Matrix &, double, double, double=0., double=NaN) const
Definition: LinearInterpolateImageFunction.hxx:1079
© 2013-2018 Imperial College London, Andreas Schuh.