ITK  4.13.0
Insight Segmentation and Registration Toolkit
itkLinearInterpolateImageFunction.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
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  * http://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 
24 namespace itk
25 {
48 template< typename TInputImage, typename TCoordRep = double >
49 class ITK_TEMPLATE_EXPORT LinearInterpolateImageFunction:
50  public InterpolateImageFunction< TInputImage, TCoordRep >
51 {
52 public:
58 
61 
63  itkNewMacro(Self);
64 
66  typedef typename Superclass::OutputType OutputType;
67 
69  typedef typename Superclass::InputImageType InputImageType;
70 
72  typedef typename Superclass::InputPixelType InputPixelType;
73 
75  typedef typename Superclass::RealType RealType;
76 
78  itkStaticConstMacro(ImageDimension, unsigned int, Superclass::ImageDimension);
79 
81  typedef typename Superclass::IndexType IndexType;
82 
84  typedef typename Superclass::ContinuousIndexType ContinuousIndexType;
85  typedef typename ContinuousIndexType::ValueType InternalComputationType;
86 
97  index) const ITK_OVERRIDE
98  {
99  return this->EvaluateOptimized(Dispatch< ImageDimension >(), index);
100  }
101 
102 protected:
104  ~LinearInterpolateImageFunction() ITK_OVERRIDE;
105  void PrintSelf(std::ostream & os, Indent indent) const ITK_OVERRIDE;
106 
107 private:
108  ITK_DISALLOW_COPY_AND_ASSIGN(LinearInterpolateImageFunction);
109 
111  static const unsigned long m_Neighbors;
112 
113  struct DispatchBase {};
114  template< unsigned int >
115  struct Dispatch: public DispatchBase {};
116 
118  const ContinuousIndexType & ) const
119  {
120  return 0;
121  }
122 
124  const ContinuousIndexType & index) const
125  {
126  IndexType basei;
127  basei[0] = Math::Floor< IndexValueType >(index[0]);
128  if ( basei[0] < this->m_StartIndex[0] )
129  {
130  basei[0] = this->m_StartIndex[0];
131  }
132 
133  const InternalComputationType & distance = index[0] - static_cast< InternalComputationType >( basei[0] );
134 
135  const TInputImage * const inputImagePtr = this->GetInputImage();
136  const RealType & val0 = inputImagePtr->GetPixel(basei);
137  if ( distance <= 0. )
138  {
139  return ( static_cast< OutputType >( val0 ) );
140  }
141 
142  ++basei[0];
143  if ( basei[0] > this->m_EndIndex[0] )
144  {
145  return ( static_cast< OutputType >( val0 ) );
146  }
147  const RealType & val1 = inputImagePtr->GetPixel(basei);
148 
149  return ( static_cast< OutputType >( val0 + ( val1 - val0 ) * distance ) );
150  }
151 
153  const ContinuousIndexType & index) const
154  {
155  IndexType basei;
156 
157  basei[0] = Math::Floor< IndexValueType >(index[0]);
158  if ( basei[0] < this->m_StartIndex[0] )
159  {
160  basei[0] = this->m_StartIndex[0];
161  }
162  const InternalComputationType & distance0 = index[0] - static_cast< InternalComputationType >( basei[0] );
163 
164  basei[1] = Math::Floor< IndexValueType >(index[1]);
165  if ( basei[1] < this->m_StartIndex[1] )
166  {
167  basei[1] = this->m_StartIndex[1];
168  }
169  const InternalComputationType & distance1 = index[1] - static_cast< InternalComputationType >( basei[1] );
170 
171  const TInputImage * const inputImagePtr = this->GetInputImage();
172  const RealType & val00 = inputImagePtr->GetPixel(basei);
173  if ( distance0 <= 0. && distance1 <= 0. )
174  {
175  return ( static_cast< OutputType >( val00 ) );
176  }
177  else if ( distance1 <= 0. ) // if they have the same "y"
178  {
179  ++basei[0]; // then interpolate across "x"
180  if ( basei[0] > this->m_EndIndex[0] )
181  {
182  return ( static_cast< OutputType >( val00 ) );
183  }
184  const RealType & val10 = inputImagePtr->GetPixel(basei);
185  return ( static_cast< OutputType >( val00 + ( val10 - val00 ) * distance0 ) );
186  }
187  else if ( distance0 <= 0. ) // if they have the same "x"
188  {
189  ++basei[1]; // then interpolate across "y"
190  if ( basei[1] > this->m_EndIndex[1] )
191  {
192  return ( static_cast< OutputType >( val00 ) );
193  }
194  const RealType & val01 = inputImagePtr->GetPixel(basei);
195  return ( static_cast< OutputType >( val00 + ( val01 - val00 ) * distance1 ) );
196  }
197  // fall-through case:
198  // interpolate across "xy"
199  ++basei[0];
200  if ( basei[0] > this->m_EndIndex[0] ) // interpolate across "y"
201  {
202  --basei[0];
203  ++basei[1];
204  if ( basei[1] > this->m_EndIndex[1] )
205  {
206  return ( static_cast< OutputType >( val00 ) );
207  }
208  const RealType & val01 = inputImagePtr->GetPixel(basei);
209  return ( static_cast< OutputType >( val00 + ( val01 - val00 ) * distance1 ) );
210  }
211  const RealType & val10 = inputImagePtr->GetPixel(basei);
212 
213  const RealType & valx0 = val00 + ( val10 - val00 ) * distance0;
214 
215  ++basei[1];
216  if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "x"
217  {
218  return ( static_cast< OutputType >( valx0 ) );
219  }
220  const RealType & val11 = inputImagePtr->GetPixel(basei);
221  --basei[0];
222  const RealType & val01 = inputImagePtr->GetPixel(basei);
223 
224  const RealType & valx1 = val01 + ( val11 - val01 ) * distance0;
225 
226  return ( static_cast< OutputType >( valx0 + ( valx1 - valx0 ) * distance1 ) );
227  }
228 
230  const ContinuousIndexType & index) const
231  {
232  IndexType basei;
233  basei[0] = Math::Floor< IndexValueType >(index[0]);
234  if ( basei[0] < this->m_StartIndex[0] )
235  {
236  basei[0] = this->m_StartIndex[0];
237  }
238  const InternalComputationType & distance0 = index[0] - static_cast< InternalComputationType >( basei[0] );
239 
240  basei[1] = Math::Floor< IndexValueType >(index[1]);
241  if ( basei[1] < this->m_StartIndex[1] )
242  {
243  basei[1] = this->m_StartIndex[1];
244  }
245  const InternalComputationType & distance1 = index[1] - static_cast< InternalComputationType >( basei[1] );
246 
247  basei[2] = Math::Floor< IndexValueType >(index[2]);
248  if ( basei[2] < this->m_StartIndex[2] )
249  {
250  basei[2] = this->m_StartIndex[2];
251  }
252  const InternalComputationType & distance2 = index[2] - static_cast< InternalComputationType >( basei[2] );
253 
254  const TInputImage * const inputImagePtr = this->GetInputImage();
255  const RealType & val000 = inputImagePtr->GetPixel(basei);
256  if ( distance0 <= 0. && distance1 <= 0. && distance2 <= 0. )
257  {
258  return ( static_cast< OutputType >( val000 ) );
259  }
260 
261  if ( distance2 <= 0. )
262  {
263  if ( distance1 <= 0. ) // interpolate across "x"
264  {
265  ++basei[0];
266  if ( basei[0] > this->m_EndIndex[0] )
267  {
268  return ( static_cast< OutputType >( val000 ) );
269  }
270  const RealType & val100 = inputImagePtr->GetPixel(basei);
271 
272  return static_cast< OutputType >( val000 + ( val100 - val000 ) * distance0 );
273  }
274  else if ( distance0 <= 0. ) // interpolate across "y"
275  {
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 
283  return static_cast< OutputType >( val000 + ( val010 - val000 ) * distance1 );
284  }
285  else // interpolate across "xy"
286  {
287  ++basei[0];
288  if ( basei[0] > this->m_EndIndex[0] ) // interpolate across "y"
289  {
290  --basei[0];
291  ++basei[1];
292  if ( basei[1] > this->m_EndIndex[1] )
293  {
294  return ( static_cast< OutputType >( val000 ) );
295  }
296  const RealType & val010 = inputImagePtr->GetPixel(basei);
297  return static_cast< OutputType >( val000 + ( val010 - val000 ) * distance1 );
298  }
299  const RealType & val100 = inputImagePtr->GetPixel(basei);
300  const RealType & valx00 = val000 + ( val100 - val000 ) * distance0;
301 
302  ++basei[1];
303  if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "x"
304  {
305  return ( static_cast< OutputType >( valx00 ) );
306  }
307  const RealType & val110 = inputImagePtr->GetPixel(basei);
308 
309  --basei[0];
310  const RealType & val010 = inputImagePtr->GetPixel(basei);
311  const RealType & valx10 = val010 + ( val110 - val010 ) * distance0;
312 
313  return static_cast< OutputType >( valx00 + ( valx10 - valx00 ) * distance1 );
314  }
315  }
316  else
317  {
318  if ( distance1 <= 0. )
319  {
320  if ( distance0 <= 0. ) // interpolate across "z"
321  {
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  else // interpolate across "xz"
332  {
333  ++basei[0];
334  if ( basei[0] > this->m_EndIndex[0] ) // interpolate across "z"
335  {
336  --basei[0];
337  ++basei[2];
338  if ( basei[2] > this->m_EndIndex[2] )
339  {
340  return ( static_cast< OutputType >( val000 ) );
341  }
342  const RealType & val001 = inputImagePtr->GetPixel(basei);
343 
344  return static_cast< OutputType >( val000 + ( val001 - val000 ) * distance2 );
345  }
346  const RealType & val100 = inputImagePtr->GetPixel(basei);
347 
348  const RealType & valx00 = val000 + ( val100 - val000 ) * distance0;
349 
350  ++basei[2];
351  if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "x"
352  {
353  return ( static_cast< OutputType >( valx00 ) );
354  }
355  const RealType & val101 = inputImagePtr->GetPixel(basei);
356 
357  --basei[0];
358  const RealType & val001 = inputImagePtr->GetPixel(basei);
359 
360  const RealType & valx01 = val001 + ( val101 - val001 ) * distance0;
361 
362  return static_cast< OutputType >( valx00 + ( valx01 - valx00 ) * distance2 );
363  }
364  }
365  else if ( distance0 <= 0. ) // interpolate across "yz"
366  {
367  ++basei[1];
368  if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "z"
369  {
370  --basei[1];
371  ++basei[2];
372  if ( basei[2] > this->m_EndIndex[2] )
373  {
374  return ( static_cast< OutputType >( val000 ) );
375  }
376  const RealType & val001 = inputImagePtr->GetPixel(basei);
377 
378  return static_cast< OutputType >( val000 + ( val001 - val000 ) * distance2 );
379  }
380  const RealType & val010 = inputImagePtr->GetPixel(basei);
381 
382  const RealType & val0x0 = val000 + ( val010 - val000 ) * distance1;
383 
384  ++basei[2];
385  if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "y"
386  {
387  return ( static_cast< OutputType >( val0x0 ) );
388  }
389  const RealType & val011 = inputImagePtr->GetPixel(basei);
390 
391  --basei[1];
392  const RealType & val001 = inputImagePtr->GetPixel(basei);
393 
394  const RealType & val0x1 = val001 + ( val011 - val001 ) * distance1;
395 
396  return static_cast< OutputType >( val0x0 + ( val0x1 - val0x0 ) * distance2 );
397  }
398  else // interpolate across "xyz"
399  {
400  ++basei[0];
401  if ( basei[0] > this->m_EndIndex[0] ) // interpolate across "yz"
402  {
403  --basei[0];
404  ++basei[1];
405  if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "z"
406  {
407  --basei[1];
408  ++basei[2];
409  if ( basei[2] > this->m_EndIndex[2] )
410  {
411  return ( static_cast< OutputType >( val000 ) );
412  }
413  const RealType & val001 = inputImagePtr->GetPixel(basei);
414 
415  return static_cast< OutputType >( val000 + ( val001 - val000 ) * distance2 );
416  }
417  const RealType & val010 = inputImagePtr->GetPixel(basei);
418  const RealType & val0x0 = val000 + ( val010 - val000 ) * distance1;
419 
420  ++basei[2];
421  if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "y"
422  {
423  return ( static_cast< OutputType >( val0x0 ) );
424  }
425  const RealType & val011 = inputImagePtr->GetPixel(basei);
426 
427  --basei[1];
428  const RealType & val001 = inputImagePtr->GetPixel(basei);
429 
430  const RealType & val0x1 = val001 + ( val011 - val001 ) * distance1;
431 
432  return static_cast< OutputType >( val0x0 + ( val0x1 - val0x0 ) * distance2 );
433  }
434  const RealType & val100 = inputImagePtr->GetPixel(basei);
435 
436  const RealType & valx00 = val000 + ( val100 - val000 ) * distance0;
437 
438  ++basei[1];
439  if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "xz"
440  {
441  --basei[1];
442  ++basei[2];
443  if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "x"
444  {
445  return ( static_cast< OutputType >( valx00 ) );
446  }
447  const RealType & val101 = inputImagePtr->GetPixel(basei);
448 
449  --basei[0];
450  const RealType & val001 = inputImagePtr->GetPixel(basei);
451 
452  const RealType & valx01 = val001 + ( val101 - val001 ) * distance0;
453 
454  return static_cast< OutputType >( valx00 + ( valx01 - valx00 ) * distance2 );
455  }
456  const RealType & val110 = inputImagePtr->GetPixel(basei);
457 
458  --basei[0];
459  const RealType & val010 = inputImagePtr->GetPixel(basei);
460 
461  const RealType & valx10 = val010 + ( val110 - val010 ) * distance0;
462 
463  const RealType & valxx0 = valx00 + ( valx10 - valx00 ) * distance1;
464 
465  ++basei[2];
466  if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "xy"
467  {
468  return ( static_cast< OutputType >( valxx0 ) );
469  }
470  const RealType & val011 = inputImagePtr->GetPixel(basei);
471 
472  ++basei[0];
473  const RealType & val111 = inputImagePtr->GetPixel(basei);
474 
475  --basei[1];
476  const RealType & val101 = inputImagePtr->GetPixel(basei);
477 
478  --basei[0];
479  const RealType & val001 = inputImagePtr->GetPixel(basei);
480 
481  const RealType & valx01 = val001 + ( val101 - val001 ) * distance0;
482  const RealType & valx11 = val011 + ( val111 - val011 ) * distance0;
483  const RealType & valxx1 = valx01 + ( valx11 - valx01 ) * distance1;
484 
485  return ( static_cast< OutputType >( valxx0 + ( valxx1 - valxx0 ) * distance2 ) );
486  }
487  }
488  }
489 
491  const ContinuousIndexType & index) const
492  {
493  return this->EvaluateUnoptimized(index);
494  }
495 
497  virtual inline OutputType EvaluateUnoptimized(
498  const ContinuousIndexType & index) const;
499 
502  template<typename RealTypeScalarRealType>
503  void
504  MakeZeroInitializer(const TInputImage * const inputImagePtr,
506  {
507  // Variable length vector version to get the size of the pixel correct.
508  typename TInputImage::IndexType idx;
509  idx.Fill(0);
510  const typename TInputImage::PixelType & tempPixel = inputImagePtr->GetPixel(idx);
511  const unsigned int sizeOfVarLengthVector = tempPixel.GetSize();
512  tempZeros.SetSize(sizeOfVarLengthVector);
514  }
515 
516  template<typename RealTypeScalarRealType>
517  void
518  MakeZeroInitializer(const TInputImage * const itkNotUsed( inputImagePtr ),
519  RealTypeScalarRealType & tempZeros) const
520  {
521  // All other cases
523  }
524 
525 };
526 } // end namespace itk
527 
528 #ifndef ITK_MANUAL_INSTANTIATION
529 #include "itkLinearInterpolateImageFunction.hxx"
530 #endif
531 
532 #endif
OutputType EvaluateOptimized(const DispatchBase &, const ContinuousIndexType &index) const
Light weight base class for most itk classes.
void MakeZeroInitializer(const TInputImage *const , RealTypeScalarRealType &tempZeros) const
Superclass::ContinuousIndexType ContinuousIndexType
virtual OutputType EvaluateAtContinuousIndex(const ContinuousIndexType &index) const override
void MakeZeroInitializer(const TInputImage *const inputImagePtr, VariableLengthVector< RealTypeScalarRealType > &tempZeros) const
A method to generically set all components to zero.
ContinuousIndexType::ValueType InternalComputationType
Represents an array whose length can be defined at run-time.
void Fill(IndexValueType value)
Definition: itkIndex.h:293
Linearly interpolate an image at specified positions.
Superclass::ContinuousIndexType ContinuousIndexType
Base class for all image interpolaters.
Control indentation during Print() invocation.
Definition: itkIndent.h:49
Define additional traits for native types such as int or float.
OutputType EvaluateOptimized(const Dispatch< 1 > &, const ContinuousIndexType &index) const
void SetSize(unsigned int sz, TReallocatePolicy reallocatePolicy, TKeepValuesPolicy keepValues)
OutputType EvaluateOptimized(const Dispatch< 0 > &, const ContinuousIndexType &) const
void Fill(TValue const &v)
OutputType EvaluateOptimized(const Dispatch< 3 > &, const ContinuousIndexType &index) const
InterpolateImageFunction< TInputImage, TCoordRep > Superclass
OutputType EvaluateOptimized(const Dispatch< 2 > &, const ContinuousIndexType &index) const
NumericTraits< typename TInputImage::PixelType >::RealType RealType