ITK  4.3.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 
21 #include "itkInterpolateImageFunction.h"
22 
23 namespace itk
24 {
47 template< class TInputImage, class TCoordRep = double >
48 class ITK_EXPORT LinearInterpolateImageFunction:
49  public InterpolateImageFunction< TInputImage, TCoordRep >
50 {
51 public:
53  typedef LinearInterpolateImageFunction Self;
54  typedef InterpolateImageFunction< TInputImage, TCoordRep > Superclass;
55  typedef SmartPointer< Self > Pointer;
56  typedef SmartPointer< const Self > ConstPointer;
57 
59  itkTypeMacro(LinearInterpolateImageFunction, InterpolateImageFunction);
60 
62  itkNewMacro(Self);
63 
65  typedef typename Superclass::OutputType OutputType;
66 
68  typedef typename Superclass::InputImageType InputImageType;
69 
71  typedef typename Superclass::InputPixelType InputPixelType;
72 
74  typedef typename Superclass::RealType RealType;
75 
77  itkStaticConstMacro(ImageDimension, unsigned int, Superclass::ImageDimension);
78 
80  typedef typename Superclass::IndexType IndexType;
81 
83  typedef typename Superclass::ContinuousIndexType ContinuousIndexType;
84 
93  virtual inline OutputType EvaluateAtContinuousIndex(const
94  ContinuousIndexType &
95  index) const
96  {
97  return this->EvaluateOptimized(Dispatch< ImageDimension >(), index);
98  }
99 
100 protected:
101  LinearInterpolateImageFunction();
102  ~LinearInterpolateImageFunction();
103  void PrintSelf(std::ostream & os, Indent indent) const;
104 
105 private:
106  LinearInterpolateImageFunction(const Self &); //purposely not implemented
107  void operator=(const Self &); //purposely not implemented
108 
110  static const unsigned long m_Neighbors;
111 
112  struct DispatchBase {};
113  template< unsigned int >
114  struct Dispatch: public DispatchBase {};
115 
116  inline OutputType EvaluateOptimized(const Dispatch< 0 > &,
117  const ContinuousIndexType & ) const
118  {
119  return 0;
120  }
121 
122  inline OutputType EvaluateOptimized(const Dispatch< 1 > &,
123  const ContinuousIndexType & index) const
124  {
125  IndexType basei;
126 
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 double distance = index[0] - static_cast< double >( basei[0] );
134 
135  const RealType val0 = this->GetInputImage()->GetPixel(basei);
136 
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 = this->GetInputImage()->GetPixel(basei);
148 
149  return ( static_cast< OutputType >( val0 + ( val1 - val0 ) * distance ) );
150  }
151 
152  inline OutputType EvaluateOptimized(const Dispatch< 2 > &,
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 double distance0 = index[0] - static_cast< double >( 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 double distance1 = index[1] - static_cast< double >( basei[1] );
170 
171  const RealType val00 = this->GetInputImage()->GetPixel(basei);
172  if ( distance0 <= 0. && distance1 <= 0. )
173  {
174  return ( static_cast< OutputType >( val00 ) );
175  }
176  else if ( distance1 <= 0. ) // if they have the same "y"
177  {
178  ++basei[0]; // then interpolate across "x"
179  if ( basei[0] > this->m_EndIndex[0] )
180  {
181  return ( static_cast< OutputType >( val00 ) );
182  }
183  const RealType val10 = this->GetInputImage()->GetPixel(basei);
184  return ( static_cast< OutputType >( val00 + ( val10 - val00 ) * distance0 ) );
185  }
186  else if ( distance0 <= 0. ) // if they have the same "x"
187  {
188  ++basei[1]; // then interpolate across "y"
189  if ( basei[1] > this->m_EndIndex[1] )
190  {
191  return ( static_cast< OutputType >( val00 ) );
192  }
193  const RealType val01 = this->GetInputImage()->GetPixel(basei);
194  return ( static_cast< OutputType >( val00 + ( val01 - val00 ) * distance1 ) );
195  }
196  // fall-through case:
197  // interpolate across "xy"
198  ++basei[0];
199  if ( basei[0] > this->m_EndIndex[0] ) // interpolate across "y"
200  {
201  --basei[0];
202  ++basei[1];
203  if ( basei[1] > this->m_EndIndex[1] )
204  {
205  return ( static_cast< OutputType >( val00 ) );
206  }
207  const RealType val01 = this->GetInputImage()->GetPixel(basei);
208  return ( static_cast< OutputType >( val00 + ( val01 - val00 ) * distance1 ) );
209  }
210  const RealType val10 = this->GetInputImage()->GetPixel(basei);
211 
212  const RealType valx0 = val00 + ( val10 - val00 ) * distance0;
213 
214  ++basei[1];
215  if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "x"
216  {
217  return ( static_cast< OutputType >( valx0 ) );
218  }
219  const RealType val11 = this->GetInputImage()->GetPixel(basei);
220  --basei[0];
221  const RealType val01 = this->GetInputImage()->GetPixel(basei);
222 
223  const RealType valx1 = val01 + ( val11 - val01 ) * distance0;
224 
225  return ( static_cast< OutputType >( valx0 + ( valx1 - valx0 ) * distance1 ) );
226  }
227 
228  inline OutputType EvaluateOptimized(const Dispatch< 3 > &,
229  const ContinuousIndexType & index) const
230  {
231  IndexType basei;
232 
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 double distance0 = index[0] - static_cast< double >( 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 double distance1 = index[1] - static_cast< double >( 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 double distance2 = index[2] - static_cast< double >( basei[2] );
253 
254  if ( distance0 <= 0. && distance1 <= 0. && distance2 <= 0. )
255  {
256  return ( static_cast< OutputType >( this->GetInputImage()->GetPixel(basei) ) );
257  }
258 
259  const RealType val000 = this->GetInputImage()->GetPixel(basei);
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 = this->GetInputImage()->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 = this->GetInputImage()->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 = this->GetInputImage()->GetPixel(basei);
297 
298  return static_cast< OutputType >( val000 + ( val010 - val000 ) * distance1 );
299  }
300  const RealType val100 = this->GetInputImage()->GetPixel(basei);
301 
302  const RealType valx00 = val000 + ( val100 - val000 ) * distance0;
303 
304  ++basei[1];
305  if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "x"
306  {
307  return ( static_cast< OutputType >( valx00 ) );
308  }
309  const RealType val110 = this->GetInputImage()->GetPixel(basei);
310 
311  --basei[0];
312  const RealType val010 = this->GetInputImage()->GetPixel(basei);
313 
314  const RealType valx10 = val010 + ( val110 - val010 ) * distance0;
315 
316  return static_cast< OutputType >( valx00 + ( valx10 - valx00 ) * distance1 );
317  }
318  }
319  else
320  {
321  if ( distance1 <= 0. )
322  {
323  if ( distance0 <= 0. ) // interpolate across "z"
324  {
325  ++basei[2];
326  if ( basei[2] > this->m_EndIndex[2] )
327  {
328  return ( static_cast< OutputType >( val000 ) );
329  }
330  const RealType val001 = this->GetInputImage()->GetPixel(basei);
331 
332  return static_cast< OutputType >( val000 + ( val001 - val000 ) * distance2 );
333  }
334  else // interpolate across "xz"
335  {
336  ++basei[0];
337  if ( basei[0] > this->m_EndIndex[0] ) // interpolate across "z"
338  {
339  --basei[0];
340  ++basei[2];
341  if ( basei[2] > this->m_EndIndex[2] )
342  {
343  return ( static_cast< OutputType >( val000 ) );
344  }
345  const RealType val001 = this->GetInputImage()->GetPixel(basei);
346 
347  return static_cast< OutputType >( val000 + ( val001 - val000 ) * distance2 );
348  }
349  const RealType val100 = this->GetInputImage()->GetPixel(basei);
350 
351  const RealType valx00 = val000 + ( val100 - val000 ) * distance0;
352 
353  ++basei[2];
354  if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "x"
355  {
356  return ( static_cast< OutputType >( valx00 ) );
357  }
358  const RealType val101 = this->GetInputImage()->GetPixel(basei);
359 
360  --basei[0];
361  const RealType val001 = this->GetInputImage()->GetPixel(basei);
362 
363  const RealType valx01 = val001 + ( val101 - val001 ) * distance0;
364 
365  return static_cast< OutputType >( valx00 + ( valx01 - valx00 ) * distance2 );
366  }
367  }
368  else if ( distance0 <= 0. ) // interpolate across "yz"
369  {
370  ++basei[1];
371  if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "z"
372  {
373  --basei[1];
374  ++basei[2];
375  if ( basei[2] > this->m_EndIndex[2] )
376  {
377  return ( static_cast< OutputType >( val000 ) );
378  }
379  const RealType val001 = this->GetInputImage()->GetPixel(basei);
380 
381  return static_cast< OutputType >( val000 + ( val001 - val000 ) * distance2 );
382  }
383  const RealType val010 = this->GetInputImage()->GetPixel(basei);
384 
385  const RealType val0x0 = val000 + ( val010 - val000 ) * distance1;
386 
387  ++basei[2];
388  if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "y"
389  {
390  return ( static_cast< OutputType >( val0x0 ) );
391  }
392  const RealType val011 = this->GetInputImage()->GetPixel(basei);
393 
394  --basei[1];
395  const RealType val001 = this->GetInputImage()->GetPixel(basei);
396 
397  const RealType val0x1 = val001 + ( val011 - val001 ) * distance1;
398 
399  return static_cast< OutputType >( val0x0 + ( val0x1 - val0x0 ) * distance2 );
400  }
401  else // interpolate across "xyz"
402  {
403  ++basei[0];
404  if ( basei[0] > this->m_EndIndex[0] ) // interpolate across "yz"
405  {
406  --basei[0];
407  ++basei[1];
408  if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "z"
409  {
410  --basei[1];
411  ++basei[2];
412  if ( basei[2] > this->m_EndIndex[2] )
413  {
414  return ( static_cast< OutputType >( val000 ) );
415  }
416  const RealType val001 = this->GetInputImage()->GetPixel(basei);
417 
418  return static_cast< OutputType >( val000 + ( val001 - val000 ) * distance2 );
419  }
420  const RealType val010 = this->GetInputImage()->GetPixel(basei);
421 
422  const RealType val0x0 = val000 + ( val010 - val000 ) * distance1;
423 
424  ++basei[2];
425  if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "y"
426  {
427  return ( static_cast< OutputType >( val0x0 ) );
428  }
429  const RealType val011 = this->GetInputImage()->GetPixel(basei);
430 
431  --basei[1];
432  const RealType val001 = this->GetInputImage()->GetPixel(basei);
433 
434  const RealType val0x1 = val001 + ( val011 - val001 ) * distance1;
435 
436  return static_cast< OutputType >( val0x0 + ( val0x1 - val0x0 ) * distance2 );
437  }
438  const RealType val100 = this->GetInputImage()->GetPixel(basei);
439 
440  const RealType valx00 = val000 + ( val100 - val000 ) * distance0;
441 
442  ++basei[1];
443  if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "xz"
444  {
445  --basei[1];
446  ++basei[2];
447  if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "x"
448  {
449  return ( static_cast< OutputType >( valx00 ) );
450  }
451  const RealType val101 = this->GetInputImage()->GetPixel(basei);
452 
453  --basei[0];
454  const RealType val001 = this->GetInputImage()->GetPixel(basei);
455 
456  const RealType valx01 = val001 + ( val101 - val001 ) * distance0;
457 
458  return static_cast< OutputType >( valx00 + ( valx01 - valx00 ) * distance2 );
459  }
460  const RealType val110 = this->GetInputImage()->GetPixel(basei);
461 
462  --basei[0];
463  const RealType val010 = this->GetInputImage()->GetPixel(basei);
464 
465  const RealType valx10 = val010 + ( val110 - val010 ) * distance0;
466 
467  const RealType valxx0 = valx00 + ( valx10 - valx00 ) * distance1;
468 
469  ++basei[2];
470  if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "xy"
471  {
472  return ( static_cast< OutputType >( valxx0 ) );
473  }
474  const RealType val011 = this->GetInputImage()->GetPixel(basei);
475 
476  ++basei[0];
477  const RealType val111 = this->GetInputImage()->GetPixel(basei);
478 
479  --basei[1];
480  const RealType val101 = this->GetInputImage()->GetPixel(basei);
481 
482  --basei[0];
483  const RealType val001 = this->GetInputImage()->GetPixel(basei);
484 
485  const RealType valx01 = val001 + ( val101 - val001 ) * distance0;
486  const RealType valx11 = val011 + ( val111 - val011 ) * distance0;
487  const RealType valxx1 = valx01 + ( valx11 - valx01 ) * distance1;
488 
489  return ( static_cast< OutputType >( valxx0 + ( valxx1 - valxx0 ) * distance2 ) );
490  }
491  }
492  }
493 
494  inline OutputType EvaluateOptimized(const DispatchBase &,
495  const ContinuousIndexType & index) const
496  {
497  return this->EvaluateUnoptimized(index);
498  }
499 
500  virtual inline OutputType EvaluateUnoptimized(
501  const ContinuousIndexType & index) const;
502 };
503 } // end namespace itk
504 
505 #ifndef ITK_MANUAL_INSTANTIATION
506 #include "itkLinearInterpolateImageFunction.hxx"
507 #endif
508 
509 #endif
510 
Tarballs of the nightly generated Doxygen documentation are available for the html, xml, and tag file.
Generated on Sun Dec 9 2012 01:17:56 for ITK by   doxygen 1.8.2