ITK  4.6.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 >
50  public InterpolateImageFunction< TInputImage, TCoordRep >
51 {
52 public:
58 
61 
63  itkNewMacro(Self);
64 
66  typedef typename Superclass::OutputType OutputType;
67 
70 
73 
75  typedef typename Superclass::RealType RealType;
76 
78  itkStaticConstMacro(ImageDimension, unsigned int, Superclass::ImageDimension);
79 
81  typedef typename Superclass::IndexType IndexType;
82 
85  typedef typename ContinuousIndexType::ValueType InternalComputationType;
86 
97  index) const
98  {
99  return this->EvaluateOptimized(Dispatch< ImageDimension >(), index);
100  }
101 
102 protected:
105  void PrintSelf(std::ostream & os, Indent indent) const;
106 
107 private:
108  LinearInterpolateImageFunction(const Self &); //purposely not implemented
109  void operator=(const Self &); //purposely not implemented
110 
112  static const unsigned long m_Neighbors;
113 
114  struct DispatchBase {};
115  template< unsigned int >
116  struct Dispatch: public DispatchBase {};
117 
119  const ContinuousIndexType & ) const
120  {
121  return 0;
122  }
123 
125  const ContinuousIndexType & index) const
126  {
127  IndexType basei;
128  basei[0] = Math::Floor< IndexValueType >(index[0]);
129  if ( basei[0] < this->m_StartIndex[0] )
130  {
131  basei[0] = this->m_StartIndex[0];
132  }
133 
134  const InternalComputationType & distance = index[0] - static_cast< InternalComputationType >( basei[0] );
135 
136  const TInputImage * const inputImagePtr = this->GetInputImage();
137  const RealType & val0 = inputImagePtr->GetPixel(basei);
138  if ( distance <= 0. )
139  {
140  return ( static_cast< OutputType >( val0 ) );
141  }
142 
143  ++basei[0];
144  if ( basei[0] > this->m_EndIndex[0] )
145  {
146  return ( static_cast< OutputType >( val0 ) );
147  }
148  const RealType & val1 = inputImagePtr->GetPixel(basei);
149 
150  return ( static_cast< OutputType >( val0 + ( val1 - val0 ) * distance ) );
151  }
152 
154  const ContinuousIndexType & index) const
155  {
156  IndexType basei;
157 
158  basei[0] = Math::Floor< IndexValueType >(index[0]);
159  if ( basei[0] < this->m_StartIndex[0] )
160  {
161  basei[0] = this->m_StartIndex[0];
162  }
163  const InternalComputationType & distance0 = index[0] - static_cast< InternalComputationType >( basei[0] );
164 
165  basei[1] = Math::Floor< IndexValueType >(index[1]);
166  if ( basei[1] < this->m_StartIndex[1] )
167  {
168  basei[1] = this->m_StartIndex[1];
169  }
170  const InternalComputationType & distance1 = index[1] - static_cast< InternalComputationType >( basei[1] );
171 
172  const TInputImage * const inputImagePtr = this->GetInputImage();
173  const RealType & val00 = inputImagePtr->GetPixel(basei);
174  if ( distance0 <= 0. && distance1 <= 0. )
175  {
176  return ( static_cast< OutputType >( val00 ) );
177  }
178  else if ( distance1 <= 0. ) // if they have the same "y"
179  {
180  ++basei[0]; // then interpolate across "x"
181  if ( basei[0] > this->m_EndIndex[0] )
182  {
183  return ( static_cast< OutputType >( val00 ) );
184  }
185  const RealType & val10 = inputImagePtr->GetPixel(basei);
186  return ( static_cast< OutputType >( val00 + ( val10 - val00 ) * distance0 ) );
187  }
188  else if ( distance0 <= 0. ) // if they have the same "x"
189  {
190  ++basei[1]; // then interpolate across "y"
191  if ( basei[1] > this->m_EndIndex[1] )
192  {
193  return ( static_cast< OutputType >( val00 ) );
194  }
195  const RealType & val01 = inputImagePtr->GetPixel(basei);
196  return ( static_cast< OutputType >( val00 + ( val01 - val00 ) * distance1 ) );
197  }
198  // fall-through case:
199  // interpolate across "xy"
200  ++basei[0];
201  if ( basei[0] > this->m_EndIndex[0] ) // interpolate across "y"
202  {
203  --basei[0];
204  ++basei[1];
205  if ( basei[1] > this->m_EndIndex[1] )
206  {
207  return ( static_cast< OutputType >( val00 ) );
208  }
209  const RealType & val01 = inputImagePtr->GetPixel(basei);
210  return ( static_cast< OutputType >( val00 + ( val01 - val00 ) * distance1 ) );
211  }
212  const RealType & val10 = inputImagePtr->GetPixel(basei);
213 
214  const RealType & valx0 = val00 + ( val10 - val00 ) * distance0;
215 
216  ++basei[1];
217  if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "x"
218  {
219  return ( static_cast< OutputType >( valx0 ) );
220  }
221  const RealType & val11 = inputImagePtr->GetPixel(basei);
222  --basei[0];
223  const RealType & val01 = inputImagePtr->GetPixel(basei);
224 
225  const RealType & valx1 = val01 + ( val11 - val01 ) * distance0;
226 
227  return ( static_cast< OutputType >( valx0 + ( valx1 - valx0 ) * distance1 ) );
228  }
229 
231  const ContinuousIndexType & index) const
232  {
233  IndexType basei;
234  basei[0] = Math::Floor< IndexValueType >(index[0]);
235  if ( basei[0] < this->m_StartIndex[0] )
236  {
237  basei[0] = this->m_StartIndex[0];
238  }
239  const InternalComputationType & distance0 = index[0] - static_cast< InternalComputationType >( basei[0] );
240 
241  basei[1] = Math::Floor< IndexValueType >(index[1]);
242  if ( basei[1] < this->m_StartIndex[1] )
243  {
244  basei[1] = this->m_StartIndex[1];
245  }
246  const InternalComputationType & distance1 = index[1] - static_cast< InternalComputationType >( basei[1] );
247 
248  basei[2] = Math::Floor< IndexValueType >(index[2]);
249  if ( basei[2] < this->m_StartIndex[2] )
250  {
251  basei[2] = this->m_StartIndex[2];
252  }
253  const InternalComputationType & distance2 = index[2] - static_cast< InternalComputationType >( basei[2] );
254 
255  const TInputImage * const inputImagePtr = this->GetInputImage();
256  const RealType & val000 = inputImagePtr->GetPixel(basei);
257  if ( distance0 <= 0. && distance1 <= 0. && distance2 <= 0. )
258  {
259  return ( static_cast< OutputType >( val000 ) );
260  }
261 
262  if ( distance2 <= 0. )
263  {
264  if ( distance1 <= 0. ) // interpolate across "x"
265  {
266  ++basei[0];
267  if ( basei[0] > this->m_EndIndex[0] )
268  {
269  return ( static_cast< OutputType >( val000 ) );
270  }
271  const RealType & val100 = inputImagePtr->GetPixel(basei);
272 
273  return static_cast< OutputType >( val000 + ( val100 - val000 ) * distance0 );
274  }
275  else if ( distance0 <= 0. ) // interpolate across "y"
276  {
277  ++basei[1];
278  if ( basei[1] > this->m_EndIndex[1] )
279  {
280  return ( static_cast< OutputType >( val000 ) );
281  }
282  const RealType & val010 = inputImagePtr->GetPixel(basei);
283 
284  return static_cast< OutputType >( val000 + ( val010 - val000 ) * distance1 );
285  }
286  else // interpolate across "xy"
287  {
288  ++basei[0];
289  if ( basei[0] > this->m_EndIndex[0] ) // interpolate across "y"
290  {
291  --basei[0];
292  ++basei[1];
293  if ( basei[1] > this->m_EndIndex[1] )
294  {
295  return ( static_cast< OutputType >( val000 ) );
296  }
297  const RealType & val010 = inputImagePtr->GetPixel(basei);
298  return static_cast< OutputType >( val000 + ( val010 - val000 ) * distance1 );
299  }
300  const RealType & val100 = inputImagePtr->GetPixel(basei);
301  const RealType & valx00 = val000 + ( val100 - val000 ) * distance0;
302 
303  ++basei[1];
304  if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "x"
305  {
306  return ( static_cast< OutputType >( valx00 ) );
307  }
308  const RealType & val110 = inputImagePtr->GetPixel(basei);
309 
310  --basei[0];
311  const RealType & val010 = inputImagePtr->GetPixel(basei);
312  const RealType & valx10 = val010 + ( val110 - val010 ) * distance0;
313 
314  return static_cast< OutputType >( valx00 + ( valx10 - valx00 ) * distance1 );
315  }
316  }
317  else
318  {
319  if ( distance1 <= 0. )
320  {
321  if ( distance0 <= 0. ) // interpolate across "z"
322  {
323  ++basei[2];
324  if ( basei[2] > this->m_EndIndex[2] )
325  {
326  return ( static_cast< OutputType >( val000 ) );
327  }
328  const RealType & val001 = inputImagePtr->GetPixel(basei);
329 
330  return static_cast< OutputType >( val000 + ( val001 - val000 ) * distance2 );
331  }
332  else // interpolate across "xz"
333  {
334  ++basei[0];
335  if ( basei[0] > this->m_EndIndex[0] ) // interpolate across "z"
336  {
337  --basei[0];
338  ++basei[2];
339  if ( basei[2] > this->m_EndIndex[2] )
340  {
341  return ( static_cast< OutputType >( val000 ) );
342  }
343  const RealType & val001 = inputImagePtr->GetPixel(basei);
344 
345  return static_cast< OutputType >( val000 + ( val001 - val000 ) * distance2 );
346  }
347  const RealType & val100 = inputImagePtr->GetPixel(basei);
348 
349  const RealType & valx00 = val000 + ( val100 - val000 ) * distance0;
350 
351  ++basei[2];
352  if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "x"
353  {
354  return ( static_cast< OutputType >( valx00 ) );
355  }
356  const RealType & val101 = inputImagePtr->GetPixel(basei);
357 
358  --basei[0];
359  const RealType & val001 = inputImagePtr->GetPixel(basei);
360 
361  const RealType & valx01 = val001 + ( val101 - val001 ) * distance0;
362 
363  return static_cast< OutputType >( valx00 + ( valx01 - valx00 ) * distance2 );
364  }
365  }
366  else if ( distance0 <= 0. ) // interpolate across "yz"
367  {
368  ++basei[1];
369  if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "z"
370  {
371  --basei[1];
372  ++basei[2];
373  if ( basei[2] > this->m_EndIndex[2] )
374  {
375  return ( static_cast< OutputType >( val000 ) );
376  }
377  const RealType & val001 = inputImagePtr->GetPixel(basei);
378 
379  return static_cast< OutputType >( val000 + ( val001 - val000 ) * distance2 );
380  }
381  const RealType & val010 = inputImagePtr->GetPixel(basei);
382 
383  const RealType & val0x0 = val000 + ( val010 - val000 ) * distance1;
384 
385  ++basei[2];
386  if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "y"
387  {
388  return ( static_cast< OutputType >( val0x0 ) );
389  }
390  const RealType & val011 = inputImagePtr->GetPixel(basei);
391 
392  --basei[1];
393  const RealType & val001 = inputImagePtr->GetPixel(basei);
394 
395  const RealType & val0x1 = val001 + ( val011 - val001 ) * distance1;
396 
397  return static_cast< OutputType >( val0x0 + ( val0x1 - val0x0 ) * distance2 );
398  }
399  else // interpolate across "xyz"
400  {
401  ++basei[0];
402  if ( basei[0] > this->m_EndIndex[0] ) // interpolate across "yz"
403  {
404  --basei[0];
405  ++basei[1];
406  if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "z"
407  {
408  --basei[1];
409  ++basei[2];
410  if ( basei[2] > this->m_EndIndex[2] )
411  {
412  return ( static_cast< OutputType >( val000 ) );
413  }
414  const RealType & val001 = inputImagePtr->GetPixel(basei);
415 
416  return static_cast< OutputType >( val000 + ( val001 - val000 ) * distance2 );
417  }
418  const RealType & val010 = inputImagePtr->GetPixel(basei);
419  const RealType & val0x0 = val000 + ( val010 - val000 ) * distance1;
420 
421  ++basei[2];
422  if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "y"
423  {
424  return ( static_cast< OutputType >( val0x0 ) );
425  }
426  const RealType & val011 = inputImagePtr->GetPixel(basei);
427 
428  --basei[1];
429  const RealType & val001 = inputImagePtr->GetPixel(basei);
430 
431  const RealType & val0x1 = val001 + ( val011 - val001 ) * distance1;
432 
433  return static_cast< OutputType >( val0x0 + ( val0x1 - val0x0 ) * distance2 );
434  }
435  const RealType & val100 = inputImagePtr->GetPixel(basei);
436 
437  const RealType & valx00 = val000 + ( val100 - val000 ) * distance0;
438 
439  ++basei[1];
440  if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "xz"
441  {
442  --basei[1];
443  ++basei[2];
444  if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "x"
445  {
446  return ( static_cast< OutputType >( valx00 ) );
447  }
448  const RealType & val101 = inputImagePtr->GetPixel(basei);
449 
450  --basei[0];
451  const RealType & val001 = inputImagePtr->GetPixel(basei);
452 
453  const RealType & valx01 = val001 + ( val101 - val001 ) * distance0;
454 
455  return static_cast< OutputType >( valx00 + ( valx01 - valx00 ) * distance2 );
456  }
457  const RealType & val110 = inputImagePtr->GetPixel(basei);
458 
459  --basei[0];
460  const RealType & val010 = inputImagePtr->GetPixel(basei);
461 
462  const RealType & valx10 = val010 + ( val110 - val010 ) * distance0;
463 
464  const RealType & valxx0 = valx00 + ( valx10 - valx00 ) * distance1;
465 
466  ++basei[2];
467  if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "xy"
468  {
469  return ( static_cast< OutputType >( valxx0 ) );
470  }
471  const RealType & val011 = inputImagePtr->GetPixel(basei);
472 
473  ++basei[0];
474  const RealType & val111 = inputImagePtr->GetPixel(basei);
475 
476  --basei[1];
477  const RealType & val101 = inputImagePtr->GetPixel(basei);
478 
479  --basei[0];
480  const RealType & val001 = inputImagePtr->GetPixel(basei);
481 
482  const RealType & valx01 = val001 + ( val101 - val001 ) * distance0;
483  const RealType & valx11 = val011 + ( val111 - val011 ) * distance0;
484  const RealType & valxx1 = valx01 + ( valx11 - valx01 ) * distance1;
485 
486  return ( static_cast< OutputType >( valxx0 + ( valxx1 - valxx0 ) * distance2 ) );
487  }
488  }
489  }
490 
492  const ContinuousIndexType & index) const
493  {
494  return this->EvaluateUnoptimized(index);
495  }
496 
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
void MakeZeroInitializer(const TInputImage *const inputImagePtr, VariableLengthVector< RealTypeScalarRealType > &tempZeros) const
A method to generically set all components to zero.
virtual OutputType EvaluateUnoptimized(const ContinuousIndexType &index) const
ContinuousIndexType::ValueType InternalComputationType
void PrintSelf(std::ostream &os, Indent indent) const
virtual OutputType EvaluateAtContinuousIndex(const ContinuousIndexType &index) const
static const unsigned int ImageDimension
Represents an array whose length can be defined at run-time.
Superclass::InputImageType InputImageType
Linearly interpolate an image at specified positions.
Superclass::ContinuousIndexType ContinuousIndexType
Base class for all image interpolaters.
void SetSize(unsigned int sz, bool destroyExistingData=true)
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
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