ITK  5.4.0
Insight Toolkit
itkLinearInterpolateImageFunction.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright NumFOCUS
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * https://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef itkLinearInterpolateImageFunction_h
19 #define itkLinearInterpolateImageFunction_h
20 
23 #include <algorithm> // For max.
24 
25 namespace itk
26 {
50 template <typename TInputImage, typename TCoordRep = double>
51 class ITK_TEMPLATE_EXPORT LinearInterpolateImageFunction : public InterpolateImageFunction<TInputImage, TCoordRep>
52 {
53 public:
54  ITK_DISALLOW_COPY_AND_MOVE(LinearInterpolateImageFunction);
55 
61 
63  itkOverrideGetNameOfClassMacro(LinearInterpolateImageFunction);
64 
66  itkNewMacro(Self);
67 
69  using typename Superclass::OutputType;
70 
72  using typename Superclass::InputImageType;
73 
75  using typename Superclass::InputPixelType;
76 
78  using typename Superclass::RealType;
79 
81  static constexpr unsigned int ImageDimension = Superclass::ImageDimension;
82 
84  using typename Superclass::IndexType;
85 
87  using typename Superclass::SizeType;
88 
90  using typename Superclass::ContinuousIndexType;
92 
101  OutputType
102  EvaluateAtContinuousIndex(const ContinuousIndexType & index) const override
103  {
104  return this->EvaluateOptimized(Dispatch<ImageDimension>(), index);
105  }
106 
107  SizeType
108  GetRadius() const override
109  {
110  return SizeType::Filled(1);
111  }
112 
113 protected:
114  LinearInterpolateImageFunction() = default;
115  ~LinearInterpolateImageFunction() override = default;
116  void
117  PrintSelf(std::ostream & os, Indent indent) const override;
118 
119 private:
121  {};
122  template <unsigned int>
123  struct Dispatch : public DispatchBase
124  {};
125 
126  inline OutputType
128  {
129  return 0;
130  }
131 
132  inline OutputType
133  EvaluateOptimized(const Dispatch<1> &, const ContinuousIndexType & index) const
134  {
135  IndexType basei;
136  basei[0] = std::max(Math::Floor<IndexValueType>(index[0]), this->m_StartIndex[0]);
137 
138  const InternalComputationType distance = index[0] - static_cast<InternalComputationType>(basei[0]);
139 
140  const TInputImage * const inputImagePtr = this->GetInputImage();
141  const RealType & val0 = inputImagePtr->GetPixel(basei);
142  if (distance <= 0.)
143  {
144  return (static_cast<OutputType>(val0));
145  }
146 
147  ++basei[0];
148  if (basei[0] > this->m_EndIndex[0])
149  {
150  return (static_cast<OutputType>(val0));
151  }
152  const RealType & val1 = inputImagePtr->GetPixel(basei);
153 
154  return (static_cast<OutputType>(val0 + (val1 - val0) * distance));
155  }
156 
157  inline OutputType
158  EvaluateOptimized(const Dispatch<2> &, const ContinuousIndexType & index) const
159  {
160  IndexType basei;
161 
162  basei[0] = std::max(Math::Floor<IndexValueType>(index[0]), this->m_StartIndex[0]);
163  const InternalComputationType distance0 = index[0] - static_cast<InternalComputationType>(basei[0]);
164 
165  basei[1] = std::max(Math::Floor<IndexValueType>(index[1]), this->m_StartIndex[1]);
166  const InternalComputationType distance1 = index[1] - static_cast<InternalComputationType>(basei[1]);
167 
168  const TInputImage * const inputImagePtr = this->GetInputImage();
169  const RealType & val00 = inputImagePtr->GetPixel(basei);
170  if (distance0 <= 0. && distance1 <= 0.)
171  {
172  return (static_cast<OutputType>(val00));
173  }
174  else if (distance1 <= 0.) // if they have the same "y"
175  {
176  ++basei[0]; // then interpolate across "x"
177  if (basei[0] > this->m_EndIndex[0])
178  {
179  return (static_cast<OutputType>(val00));
180  }
181  const RealType & val10 = inputImagePtr->GetPixel(basei);
182  return (static_cast<OutputType>(val00 + (val10 - val00) * distance0));
183  }
184  else if (distance0 <= 0.) // if they have the same "x"
185  {
186  ++basei[1]; // then interpolate across "y"
187  if (basei[1] > this->m_EndIndex[1])
188  {
189  return (static_cast<OutputType>(val00));
190  }
191  const RealType & val01 = inputImagePtr->GetPixel(basei);
192  return (static_cast<OutputType>(val00 + (val01 - val00) * distance1));
193  }
194  // fall-through case:
195  // interpolate across "xy"
196  ++basei[0];
197  if (basei[0] > this->m_EndIndex[0]) // interpolate across "y"
198  {
199  --basei[0];
200  ++basei[1];
201  if (basei[1] > this->m_EndIndex[1])
202  {
203  return (static_cast<OutputType>(val00));
204  }
205  const RealType & val01 = inputImagePtr->GetPixel(basei);
206  return (static_cast<OutputType>(val00 + (val01 - val00) * distance1));
207  }
208  const RealType & val10 = inputImagePtr->GetPixel(basei);
209 
210  const RealType & valx0 = val00 + (val10 - val00) * distance0;
211 
212  ++basei[1];
213  if (basei[1] > this->m_EndIndex[1]) // interpolate across "x"
214  {
215  return (static_cast<OutputType>(valx0));
216  }
217  const RealType & val11 = inputImagePtr->GetPixel(basei);
218  --basei[0];
219  const RealType & val01 = inputImagePtr->GetPixel(basei);
220 
221  const RealType & valx1 = val01 + (val11 - val01) * distance0;
222 
223  return (static_cast<OutputType>(valx0 + (valx1 - valx0) * distance1));
224  }
225 
226  inline OutputType
227  EvaluateOptimized(const Dispatch<3> &, const ContinuousIndexType & index) const
228  {
229  IndexType basei;
230  basei[0] = std::max(Math::Floor<IndexValueType>(index[0]), this->m_StartIndex[0]);
231  const InternalComputationType distance0 = index[0] - static_cast<InternalComputationType>(basei[0]);
232 
233  basei[1] = std::max(Math::Floor<IndexValueType>(index[1]), this->m_StartIndex[1]);
234  const InternalComputationType distance1 = index[1] - static_cast<InternalComputationType>(basei[1]);
235 
236  basei[2] = std::max(Math::Floor<IndexValueType>(index[2]), this->m_StartIndex[2]);
237  const InternalComputationType distance2 = index[2] - static_cast<InternalComputationType>(basei[2]);
238 
239  const TInputImage * const inputImagePtr = this->GetInputImage();
240  const RealType & val000 = inputImagePtr->GetPixel(basei);
241  if (distance0 <= 0. && distance1 <= 0. && distance2 <= 0.)
242  {
243  return (static_cast<OutputType>(val000));
244  }
245 
246  if (distance2 <= 0.)
247  {
248  if (distance1 <= 0.) // interpolate across "x"
249  {
250  ++basei[0];
251  if (basei[0] > this->m_EndIndex[0])
252  {
253  return (static_cast<OutputType>(val000));
254  }
255  const RealType & val100 = inputImagePtr->GetPixel(basei);
256 
257  return static_cast<OutputType>(val000 + (val100 - val000) * distance0);
258  }
259  else if (distance0 <= 0.) // interpolate across "y"
260  {
261  ++basei[1];
262  if (basei[1] > this->m_EndIndex[1])
263  {
264  return (static_cast<OutputType>(val000));
265  }
266  const RealType & val010 = inputImagePtr->GetPixel(basei);
267 
268  return static_cast<OutputType>(val000 + (val010 - val000) * distance1);
269  }
270  else // interpolate across "xy"
271  {
272  ++basei[0];
273  if (basei[0] > this->m_EndIndex[0]) // interpolate across "y"
274  {
275  --basei[0];
276  ++basei[1];
277  if (basei[1] > this->m_EndIndex[1])
278  {
279  return (static_cast<OutputType>(val000));
280  }
281  const RealType & val010 = inputImagePtr->GetPixel(basei);
282  return static_cast<OutputType>(val000 + (val010 - val000) * distance1);
283  }
284  const RealType & val100 = inputImagePtr->GetPixel(basei);
285  const RealType & valx00 = val000 + (val100 - val000) * distance0;
286 
287  ++basei[1];
288  if (basei[1] > this->m_EndIndex[1]) // interpolate across "x"
289  {
290  return (static_cast<OutputType>(valx00));
291  }
292  const RealType & val110 = inputImagePtr->GetPixel(basei);
293 
294  --basei[0];
295  const RealType & val010 = inputImagePtr->GetPixel(basei);
296  const RealType & valx10 = val010 + (val110 - val010) * distance0;
297 
298  return static_cast<OutputType>(valx00 + (valx10 - valx00) * distance1);
299  }
300  }
301  else
302  {
303  if (distance1 <= 0.)
304  {
305  if (distance0 <= 0.) // interpolate across "z"
306  {
307  ++basei[2];
308  if (basei[2] > this->m_EndIndex[2])
309  {
310  return (static_cast<OutputType>(val000));
311  }
312  const RealType & val001 = inputImagePtr->GetPixel(basei);
313 
314  return static_cast<OutputType>(val000 + (val001 - val000) * distance2);
315  }
316  else // interpolate across "xz"
317  {
318  ++basei[0];
319  if (basei[0] > this->m_EndIndex[0]) // interpolate across "z"
320  {
321  --basei[0];
322  ++basei[2];
323  if (basei[2] > this->m_EndIndex[2])
324  {
325  return (static_cast<OutputType>(val000));
326  }
327  const RealType & val001 = inputImagePtr->GetPixel(basei);
328 
329  return static_cast<OutputType>(val000 + (val001 - val000) * distance2);
330  }
331  const RealType & val100 = inputImagePtr->GetPixel(basei);
332 
333  const RealType & valx00 = val000 + (val100 - val000) * distance0;
334 
335  ++basei[2];
336  if (basei[2] > this->m_EndIndex[2]) // interpolate across "x"
337  {
338  return (static_cast<OutputType>(valx00));
339  }
340  const RealType & val101 = inputImagePtr->GetPixel(basei);
341 
342  --basei[0];
343  const RealType & val001 = inputImagePtr->GetPixel(basei);
344 
345  const RealType & valx01 = val001 + (val101 - val001) * distance0;
346 
347  return static_cast<OutputType>(valx00 + (valx01 - valx00) * distance2);
348  }
349  }
350  else if (distance0 <= 0.) // interpolate across "yz"
351  {
352  ++basei[1];
353  if (basei[1] > this->m_EndIndex[1]) // interpolate across "z"
354  {
355  --basei[1];
356  ++basei[2];
357  if (basei[2] > this->m_EndIndex[2])
358  {
359  return (static_cast<OutputType>(val000));
360  }
361  const RealType & val001 = inputImagePtr->GetPixel(basei);
362 
363  return static_cast<OutputType>(val000 + (val001 - val000) * distance2);
364  }
365  const RealType & val010 = inputImagePtr->GetPixel(basei);
366 
367  const RealType & val0x0 = val000 + (val010 - val000) * distance1;
368 
369  ++basei[2];
370  if (basei[2] > this->m_EndIndex[2]) // interpolate across "y"
371  {
372  return (static_cast<OutputType>(val0x0));
373  }
374  const RealType & val011 = inputImagePtr->GetPixel(basei);
375 
376  --basei[1];
377  const RealType & val001 = inputImagePtr->GetPixel(basei);
378 
379  const RealType & val0x1 = val001 + (val011 - val001) * distance1;
380 
381  return static_cast<OutputType>(val0x0 + (val0x1 - val0x0) * distance2);
382  }
383  else // interpolate across "xyz"
384  {
385  ++basei[0];
386  if (basei[0] > this->m_EndIndex[0]) // interpolate across "yz"
387  {
388  --basei[0];
389  ++basei[1];
390  if (basei[1] > this->m_EndIndex[1]) // interpolate across "z"
391  {
392  --basei[1];
393  ++basei[2];
394  if (basei[2] > this->m_EndIndex[2])
395  {
396  return (static_cast<OutputType>(val000));
397  }
398  const RealType & val001 = inputImagePtr->GetPixel(basei);
399 
400  return static_cast<OutputType>(val000 + (val001 - val000) * distance2);
401  }
402  const RealType & val010 = inputImagePtr->GetPixel(basei);
403  const RealType & val0x0 = val000 + (val010 - val000) * distance1;
404 
405  ++basei[2];
406  if (basei[2] > this->m_EndIndex[2]) // interpolate across "y"
407  {
408  return (static_cast<OutputType>(val0x0));
409  }
410  const RealType & val011 = inputImagePtr->GetPixel(basei);
411 
412  --basei[1];
413  const RealType & val001 = inputImagePtr->GetPixel(basei);
414 
415  const RealType & val0x1 = val001 + (val011 - val001) * distance1;
416 
417  return static_cast<OutputType>(val0x0 + (val0x1 - val0x0) * distance2);
418  }
419  const RealType & val100 = inputImagePtr->GetPixel(basei);
420 
421  const RealType & valx00 = val000 + (val100 - val000) * distance0;
422 
423  ++basei[1];
424  if (basei[1] > this->m_EndIndex[1]) // interpolate across "xz"
425  {
426  --basei[1];
427  ++basei[2];
428  if (basei[2] > this->m_EndIndex[2]) // interpolate across "x"
429  {
430  return (static_cast<OutputType>(valx00));
431  }
432  const RealType & val101 = inputImagePtr->GetPixel(basei);
433 
434  --basei[0];
435  const RealType & val001 = inputImagePtr->GetPixel(basei);
436 
437  const RealType & valx01 = val001 + (val101 - val001) * distance0;
438 
439  return static_cast<OutputType>(valx00 + (valx01 - valx00) * distance2);
440  }
441  const RealType & val110 = inputImagePtr->GetPixel(basei);
442 
443  --basei[0];
444  const RealType & val010 = inputImagePtr->GetPixel(basei);
445 
446  const RealType & valx10 = val010 + (val110 - val010) * distance0;
447 
448  const RealType & valxx0 = valx00 + (valx10 - valx00) * distance1;
449 
450  ++basei[2];
451  if (basei[2] > this->m_EndIndex[2]) // interpolate across "xy"
452  {
453  return (static_cast<OutputType>(valxx0));
454  }
455  const RealType & val011 = inputImagePtr->GetPixel(basei);
456 
457  ++basei[0];
458  const RealType & val111 = inputImagePtr->GetPixel(basei);
459 
460  --basei[1];
461  const RealType & val101 = inputImagePtr->GetPixel(basei);
462 
463  --basei[0];
464  const RealType & val001 = inputImagePtr->GetPixel(basei);
465 
466  const RealType & valx01 = val001 + (val101 - val001) * distance0;
467  const RealType & valx11 = val011 + (val111 - val011) * distance0;
468  const RealType & valxx1 = valx01 + (valx11 - valx01) * distance1;
469 
470  return (static_cast<OutputType>(valxx0 + (valxx1 - valxx0) * distance2));
471  }
472  }
473  }
474 
475  inline OutputType
476  EvaluateOptimized(const DispatchBase &, const ContinuousIndexType & index) const
477  {
478  return this->EvaluateUnoptimized(index);
479  }
480 
482  virtual inline OutputType
483  EvaluateUnoptimized(const ContinuousIndexType & index) const;
484 
487  template <typename RealTypeScalarRealType>
488  void
489  MakeZeroInitializer(const TInputImage * const inputImagePtr,
491  {
492  // Variable length vector version to get the size of the pixel correct.
493  constexpr typename TInputImage::IndexType idx = { { 0 } };
494  const typename TInputImage::PixelType & tempPixel = inputImagePtr->GetPixel(idx);
495  const unsigned int sizeOfVarLengthVector = tempPixel.GetSize();
496  tempZeros.SetSize(sizeOfVarLengthVector);
497  tempZeros.Fill(RealTypeScalarRealType{});
498  }
499 
500  template <typename RealTypeScalarRealType>
501  void
502  MakeZeroInitializer(const TInputImage * const itkNotUsed(inputImagePtr), RealTypeScalarRealType & tempZeros) const
503  {
504  // All other cases
505  tempZeros = RealTypeScalarRealType{};
506  }
507 };
508 } // end namespace itk
509 
510 #ifndef ITK_MANUAL_INSTANTIATION
511 # include "itkLinearInterpolateImageFunction.hxx"
512 #endif
513 
514 #endif
itk::VariableLengthVector::Fill
void Fill(TValue const &v)
itk::LinearInterpolateImageFunction::EvaluateOptimized
OutputType EvaluateOptimized(const DispatchBase &, const ContinuousIndexType &index) const
Definition: itkLinearInterpolateImageFunction.h:476
itk::ImageFunction< TInputImage, NumericTraits< TInputImage::PixelType >::RealType, TCoordRep >::IndexType
typename InputImageType::IndexType IndexType
Definition: itkImageFunction.h:90
itkVariableLengthVector.h
itk::LinearInterpolateImageFunction::EvaluateOptimized
OutputType EvaluateOptimized(const Dispatch< 3 > &, const ContinuousIndexType &index) const
Definition: itkLinearInterpolateImageFunction.h:227
itk::VariableLengthVector::SetSize
void SetSize(unsigned int sz, TReallocatePolicy reallocatePolicy, TKeepValuesPolicy keepValues)
itk::GTest::TypedefsAndConstructors::Dimension2::SizeType
ImageBaseType::SizeType SizeType
Definition: itkGTestTypedefsAndConstructors.h:49
itk::LinearInterpolateImageFunction::EvaluateOptimized
OutputType EvaluateOptimized(const Dispatch< 2 > &, const ContinuousIndexType &index) const
Definition: itkLinearInterpolateImageFunction.h:158
itk::SmartPointer< Self >
itk::Indent
Control indentation during Print() invocation.
Definition: itkIndent.h:49
itk::FunctionBase< Point< TCoordRep, TInputImage::ImageDimension >, NumericTraits< TInputImage::PixelType >::RealType >::OutputType
NumericTraits< TInputImage::PixelType >::RealType OutputType
Definition: itkFunctionBase.h:62
itk::LinearInterpolateImageFunction::MakeZeroInitializer
void MakeZeroInitializer(const TInputImage *const, RealTypeScalarRealType &tempZeros) const
Definition: itkLinearInterpolateImageFunction.h:502
itk::GTest::TypedefsAndConstructors::Dimension2::IndexType
ImageBaseType::IndexType IndexType
Definition: itkGTestTypedefsAndConstructors.h:50
itk::LinearInterpolateImageFunction
Linearly interpolate an image at specified positions.
Definition: itkLinearInterpolateImageFunction.h:51
itk::LightObject
Light weight base class for most itk classes.
Definition: itkLightObject.h:55
itk::LinearInterpolateImageFunction::InternalComputationType
typename ContinuousIndexType::ValueType InternalComputationType
Definition: itkLinearInterpolateImageFunction.h:91
itk::LinearInterpolateImageFunction::MakeZeroInitializer
void MakeZeroInitializer(const TInputImage *const inputImagePtr, VariableLengthVector< RealTypeScalarRealType > &tempZeros) const
A method to generically set all components to zero.
Definition: itkLinearInterpolateImageFunction.h:489
itk::LinearInterpolateImageFunction::DispatchBase
Definition: itkLinearInterpolateImageFunction.h:120
itk::VariableLengthVector
Represents an array whose length can be defined at run-time.
Definition: itkConstantBoundaryCondition.h:28
itk::LinearInterpolateImageFunction::Dispatch
Definition: itkLinearInterpolateImageFunction.h:123
itk
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....
Definition: itkAnnulusOperator.h:24
itk::ContinuousIndex< TCoordRep, Self::ImageDimension >
itk::ContinuousIndex< TCoordRep, Self::ImageDimension >::ValueType
TCoordRep ValueType
Definition: itkContinuousIndex.h:58
itk::LinearInterpolateImageFunction::GetRadius
SizeType GetRadius() const override
Definition: itkLinearInterpolateImageFunction.h:108
itkInterpolateImageFunction.h
itk::LinearInterpolateImageFunction::EvaluateAtContinuousIndex
OutputType EvaluateAtContinuousIndex(const ContinuousIndexType &index) const override
Definition: itkLinearInterpolateImageFunction.h:102
itk::LinearInterpolateImageFunction::EvaluateOptimized
OutputType EvaluateOptimized(const Dispatch< 0 > &, const ContinuousIndexType &) const
Definition: itkLinearInterpolateImageFunction.h:127
itk::InterpolateImageFunction
Base class for all image interpolators.
Definition: itkInterpolateImageFunction.h:45
itk::LinearInterpolateImageFunction::EvaluateOptimized
OutputType EvaluateOptimized(const Dispatch< 1 > &, const ContinuousIndexType &index) const
Definition: itkLinearInterpolateImageFunction.h:133
itk::InterpolateImageFunction::RealType
typename NumericTraits< typename TInputImage::PixelType >::RealType RealType
Definition: itkInterpolateImageFunction.h:85