ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkLinearInterpolateImageFunction.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 #ifndef __itkLinearInterpolateImageFunction_h
00019 #define __itkLinearInterpolateImageFunction_h
00020 
00021 #include "itkInterpolateImageFunction.h"
00022 
00023 namespace itk
00024 {
00047 template< class TInputImage, class TCoordRep = double >
00048 class ITK_EXPORT LinearInterpolateImageFunction:
00049   public InterpolateImageFunction< TInputImage, TCoordRep >
00050 {
00051 public:
00053   typedef LinearInterpolateImageFunction                     Self;
00054   typedef InterpolateImageFunction< TInputImage, TCoordRep > Superclass;
00055   typedef SmartPointer< Self >                               Pointer;
00056   typedef SmartPointer< const Self >                         ConstPointer;
00057 
00059   itkTypeMacro(LinearInterpolateImageFunction, InterpolateImageFunction);
00060 
00062   itkNewMacro(Self);
00063 
00065   typedef typename Superclass::OutputType OutputType;
00066 
00068   typedef typename Superclass::InputImageType InputImageType;
00069 
00071   typedef typename Superclass::InputPixelType InputPixelType;
00072 
00074   typedef typename Superclass::RealType RealType;
00075 
00077   itkStaticConstMacro(ImageDimension, unsigned int, Superclass::ImageDimension);
00078 
00080   typedef typename Superclass::IndexType      IndexType;
00081 
00083   typedef typename Superclass::ContinuousIndexType ContinuousIndexType;
00084 
00093   virtual inline OutputType EvaluateAtContinuousIndex(const
00094                                                       ContinuousIndexType &
00095                                                       index) const
00096   {
00097     return this->EvaluateOptimized(Dispatch< ImageDimension >(), index);
00098   }
00099 
00100 protected:
00101   LinearInterpolateImageFunction();
00102   ~LinearInterpolateImageFunction();
00103   void PrintSelf(std::ostream & os, Indent indent) const;
00104 
00105 private:
00106   LinearInterpolateImageFunction(const Self &); //purposely not implemented
00107   void operator=(const Self &);                 //purposely not implemented
00108 
00110   static const unsigned long m_Neighbors;
00111 
00112   struct DispatchBase {};
00113   template< unsigned int >
00114   struct Dispatch: public DispatchBase {};
00115 
00116   inline OutputType EvaluateOptimized(const Dispatch< 0 > &,
00117                                       const ContinuousIndexType & ) const
00118   {
00119     return 0;
00120   }
00121 
00122   inline OutputType EvaluateOptimized(const Dispatch< 1 > &,
00123                                       const ContinuousIndexType & index) const
00124   {
00125     IndexType basei;
00126 
00127     basei[0] = Math::Floor< IndexValueType >(index[0]);
00128     if ( basei[0] < this->m_StartIndex[0] )
00129       {
00130       basei[0] = this->m_StartIndex[0];
00131       }
00132 
00133     const double distance = index[0] - static_cast< double >( basei[0] );
00134 
00135     const RealType val0 = this->GetInputImage()->GetPixel(basei);
00136 
00137     if ( distance <= 0. )
00138       {
00139       return ( static_cast< OutputType >( val0 ) );
00140       }
00141 
00142     ++basei[0];
00143     if ( basei[0] > this->m_EndIndex[0] )
00144       {
00145       return ( static_cast< OutputType >( val0 ) );
00146       }
00147     const RealType val1 = this->GetInputImage()->GetPixel(basei);
00148 
00149     return ( static_cast< OutputType >( val0 + ( val1 - val0 ) * distance ) );
00150   }
00151 
00152   inline OutputType EvaluateOptimized(const Dispatch< 2 > &,
00153                                       const ContinuousIndexType & index) const
00154   {
00155     IndexType basei;
00156 
00157     basei[0] = Math::Floor< IndexValueType >(index[0]);
00158     if ( basei[0] < this->m_StartIndex[0] )
00159       {
00160       basei[0] = this->m_StartIndex[0];
00161       }
00162     const double distance0 = index[0] - static_cast< double >( basei[0] );
00163 
00164     basei[1] = Math::Floor< IndexValueType >(index[1]);
00165     if ( basei[1] < this->m_StartIndex[1] )
00166       {
00167       basei[1] = this->m_StartIndex[1];
00168       }
00169     const double distance1 = index[1] - static_cast< double >( basei[1] );
00170 
00171     const RealType val00 = this->GetInputImage()->GetPixel(basei);
00172     if ( distance0 <= 0. && distance1 <= 0. )
00173       {
00174       return ( static_cast< OutputType >( val00 ) );
00175       }
00176     else if ( distance1 <= 0. ) // if they have the same "y"
00177       {
00178       ++basei[0];  // then interpolate across "x"
00179       if ( basei[0] > this->m_EndIndex[0] )
00180         {
00181         return ( static_cast< OutputType >( val00 ) );
00182         }
00183       const RealType val10 = this->GetInputImage()->GetPixel(basei);
00184       return ( static_cast< OutputType >( val00 + ( val10 - val00 ) * distance0 ) );
00185       }
00186     else if ( distance0 <= 0. ) // if they have the same "x"
00187       {
00188       ++basei[1];  // then interpolate across "y"
00189       if ( basei[1] > this->m_EndIndex[1] )
00190         {
00191         return ( static_cast< OutputType >( val00 ) );
00192         }
00193       const RealType val01 = this->GetInputImage()->GetPixel(basei);
00194       return ( static_cast< OutputType >( val00 + ( val01 - val00 ) * distance1 ) );
00195       }
00196     // fall-through case:
00197     // interpolate across "xy"
00198     ++basei[0];
00199     if ( basei[0] > this->m_EndIndex[0] ) // interpolate across "y"
00200       {
00201       --basei[0];
00202       ++basei[1];
00203       if ( basei[1] > this->m_EndIndex[1] )
00204         {
00205         return ( static_cast< OutputType >( val00 ) );
00206         }
00207       const RealType val01 = this->GetInputImage()->GetPixel(basei);
00208       return ( static_cast< OutputType >( val00 + ( val01 - val00 ) * distance1 ) );
00209       }
00210     const RealType val10 = this->GetInputImage()->GetPixel(basei);
00211 
00212     const RealType valx0 = val00 + ( val10 - val00 ) * distance0;
00213 
00214     ++basei[1];
00215     if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "x"
00216       {
00217       return ( static_cast< OutputType >( valx0 ) );
00218       }
00219     const RealType val11 = this->GetInputImage()->GetPixel(basei);
00220     --basei[0];
00221     const RealType val01 = this->GetInputImage()->GetPixel(basei);
00222 
00223     const RealType valx1 = val01 + ( val11 - val01 ) * distance0;
00224 
00225     return ( static_cast< OutputType >( valx0 + ( valx1 - valx0 ) * distance1 ) );
00226   }
00227 
00228   inline OutputType EvaluateOptimized(const Dispatch< 3 > &,
00229                                       const ContinuousIndexType & index) const
00230   {
00231     IndexType basei;
00232 
00233     basei[0] = Math::Floor< IndexValueType >(index[0]);
00234     if ( basei[0] < this->m_StartIndex[0] )
00235       {
00236       basei[0] = this->m_StartIndex[0];
00237       }
00238     const double distance0 = index[0] - static_cast< double >( basei[0] );
00239 
00240     basei[1] = Math::Floor< IndexValueType >(index[1]);
00241     if ( basei[1] < this->m_StartIndex[1] )
00242       {
00243       basei[1] = this->m_StartIndex[1];
00244       }
00245     const double distance1 = index[1] - static_cast< double >( basei[1] );
00246 
00247     basei[2] = Math::Floor< IndexValueType >(index[2]);
00248     if ( basei[2] < this->m_StartIndex[2] )
00249       {
00250       basei[2] = this->m_StartIndex[2];
00251       }
00252     const double distance2 = index[2] - static_cast< double >( basei[2] );
00253 
00254     if ( distance0 <= 0. && distance1 <= 0. && distance2 <= 0. )
00255       {
00256       return ( static_cast< OutputType >( this->GetInputImage()->GetPixel(basei) ) );
00257       }
00258 
00259     const RealType val000 = this->GetInputImage()->GetPixel(basei);
00260 
00261     if ( distance2 <= 0. )
00262       {
00263       if ( distance1 <= 0. ) // interpolate across "x"
00264         {
00265         ++basei[0];
00266         if ( basei[0] > this->m_EndIndex[0] )
00267           {
00268           return ( static_cast< OutputType >( val000 ) );
00269           }
00270         const RealType val100 = this->GetInputImage()->GetPixel(basei);
00271 
00272         return static_cast< OutputType >( val000 + ( val100 - val000 ) * distance0 );
00273         }
00274       else if ( distance0 <= 0. ) // interpolate across "y"
00275         {
00276         ++basei[1];
00277         if ( basei[1] > this->m_EndIndex[1] )
00278           {
00279           return ( static_cast< OutputType >( val000 ) );
00280           }
00281         const RealType val010 = this->GetInputImage()->GetPixel(basei);
00282 
00283         return static_cast< OutputType >( val000 + ( val010 - val000 ) * distance1 );
00284         }
00285       else  // interpolate across "xy"
00286         {
00287         ++basei[0];
00288         if ( basei[0] > this->m_EndIndex[0] ) // interpolate across "y"
00289           {
00290           --basei[0];
00291           ++basei[1];
00292           if ( basei[1] > this->m_EndIndex[1] )
00293             {
00294             return ( static_cast< OutputType >( val000 ) );
00295             }
00296           const RealType val010 = this->GetInputImage()->GetPixel(basei);
00297 
00298           return static_cast< OutputType >( val000 + ( val010 - val000 ) * distance1 );
00299           }
00300         const RealType val100 = this->GetInputImage()->GetPixel(basei);
00301 
00302         const RealType valx00 = val000 + ( val100 - val000 ) * distance0;
00303 
00304         ++basei[1];
00305         if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "x"
00306           {
00307           return ( static_cast< OutputType >( valx00 ) );
00308           }
00309         const RealType val110 = this->GetInputImage()->GetPixel(basei);
00310 
00311         --basei[0];
00312         const RealType val010 = this->GetInputImage()->GetPixel(basei);
00313 
00314         const RealType valx10 = val010 + ( val110 - val010 ) * distance0;
00315 
00316         return static_cast< OutputType >( valx00 + ( valx10 - valx00 ) * distance1 );
00317         }
00318       }
00319     else
00320       {
00321       if ( distance1 <= 0. )
00322         {
00323         if ( distance0 <= 0. ) // interpolate across "z"
00324           {
00325           ++basei[2];
00326           if ( basei[2] > this->m_EndIndex[2] )
00327             {
00328             return ( static_cast< OutputType >( val000 ) );
00329             }
00330           const RealType val001 = this->GetInputImage()->GetPixel(basei);
00331 
00332           return static_cast< OutputType >( val000 + ( val001 - val000 ) * distance2 );
00333           }
00334         else // interpolate across "xz"
00335           {
00336           ++basei[0];
00337           if ( basei[0] > this->m_EndIndex[0] ) // interpolate across "z"
00338             {
00339             --basei[0];
00340             ++basei[2];
00341             if ( basei[2] > this->m_EndIndex[2] )
00342               {
00343               return ( static_cast< OutputType >( val000 ) );
00344               }
00345             const RealType val001 = this->GetInputImage()->GetPixel(basei);
00346 
00347             return static_cast< OutputType >( val000 + ( val001 - val000 ) * distance2 );
00348             }
00349           const RealType val100 = this->GetInputImage()->GetPixel(basei);
00350 
00351           const RealType valx00 = val000 + ( val100 - val000 ) * distance0;
00352 
00353           ++basei[2];
00354           if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "x"
00355             {
00356             return ( static_cast< OutputType >( valx00 ) );
00357             }
00358           const RealType val101 = this->GetInputImage()->GetPixel(basei);
00359 
00360           --basei[0];
00361           const RealType val001 = this->GetInputImage()->GetPixel(basei);
00362 
00363           const RealType valx01 = val001 + ( val101 - val001 ) * distance0;
00364 
00365           return static_cast< OutputType >( valx00 + ( valx01 - valx00 ) * distance2 );
00366           }
00367         }
00368       else if ( distance0 <= 0. ) // interpolate across "yz"
00369         {
00370         ++basei[1];
00371         if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "z"
00372           {
00373           --basei[1];
00374           ++basei[2];
00375           if ( basei[2] > this->m_EndIndex[2] )
00376             {
00377             return ( static_cast< OutputType >( val000 ) );
00378             }
00379           const RealType val001 = this->GetInputImage()->GetPixel(basei);
00380 
00381           return static_cast< OutputType >( val000 + ( val001 - val000 ) * distance2 );
00382           }
00383         const RealType val010 = this->GetInputImage()->GetPixel(basei);
00384 
00385         const RealType val0x0 = val000 + ( val010 - val000 ) * distance1;
00386 
00387         ++basei[2];
00388         if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "y"
00389           {
00390           return ( static_cast< OutputType >( val0x0 ) );
00391           }
00392         const RealType val011 = this->GetInputImage()->GetPixel(basei);
00393 
00394         --basei[1];
00395         const RealType val001 = this->GetInputImage()->GetPixel(basei);
00396 
00397         const RealType val0x1 = val001 + ( val011 - val001 ) * distance1;
00398 
00399         return static_cast< OutputType >( val0x0 + ( val0x1 - val0x0 ) * distance2 );
00400         }
00401       else // interpolate across "xyz"
00402         {
00403         ++basei[0];
00404         if ( basei[0] > this->m_EndIndex[0] ) // interpolate across "yz"
00405           {
00406           --basei[0];
00407           ++basei[1];
00408           if ( basei[1] > this->m_EndIndex[1] )  // interpolate across "z"
00409             {
00410             --basei[1];
00411             ++basei[2];
00412             if ( basei[2] > this->m_EndIndex[2] )
00413               {
00414               return ( static_cast< OutputType >( val000 ) );
00415               }
00416             const RealType val001 = this->GetInputImage()->GetPixel(basei);
00417 
00418             return static_cast< OutputType >( val000 + ( val001 - val000 ) * distance2 );
00419             }
00420           const RealType val010 = this->GetInputImage()->GetPixel(basei);
00421 
00422           const RealType val0x0 = val000 + ( val010 - val000 ) * distance1;
00423 
00424           ++basei[2];
00425           if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "y"
00426             {
00427             return ( static_cast< OutputType >( val0x0 ) );
00428             }
00429           const RealType val011 = this->GetInputImage()->GetPixel(basei);
00430 
00431           --basei[1];
00432           const RealType val001 = this->GetInputImage()->GetPixel(basei);
00433 
00434           const RealType val0x1 = val001 + ( val011 - val001 ) * distance1;
00435 
00436           return static_cast< OutputType >( val0x0 + ( val0x1 - val0x0 ) * distance2 );
00437           }
00438         const RealType val100 = this->GetInputImage()->GetPixel(basei);
00439 
00440         const RealType valx00 = val000 + ( val100 - val000 ) * distance0;
00441 
00442         ++basei[1];
00443         if ( basei[1] > this->m_EndIndex[1] ) // interpolate across "xz"
00444           {
00445           --basei[1];
00446           ++basei[2];
00447           if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "x"
00448             {
00449             return ( static_cast< OutputType >( valx00 ) );
00450             }
00451           const RealType val101 = this->GetInputImage()->GetPixel(basei);
00452 
00453           --basei[0];
00454           const RealType val001 = this->GetInputImage()->GetPixel(basei);
00455 
00456           const RealType valx01 = val001 + ( val101 - val001 ) * distance0;
00457 
00458           return static_cast< OutputType >( valx00 + ( valx01 - valx00 ) * distance2 );
00459           }
00460         const RealType val110 = this->GetInputImage()->GetPixel(basei);
00461 
00462         --basei[0];
00463         const RealType val010 = this->GetInputImage()->GetPixel(basei);
00464 
00465         const RealType valx10 = val010 + ( val110 - val010 ) * distance0;
00466 
00467         const RealType valxx0 = valx00 + ( valx10 - valx00 ) * distance1;
00468 
00469         ++basei[2];
00470         if ( basei[2] > this->m_EndIndex[2] ) // interpolate across "xy"
00471           {
00472           return ( static_cast< OutputType >( valxx0 ) );
00473           }
00474         const RealType val011 = this->GetInputImage()->GetPixel(basei);
00475 
00476         ++basei[0];
00477         const RealType val111 = this->GetInputImage()->GetPixel(basei);
00478 
00479         --basei[1];
00480         const RealType val101 = this->GetInputImage()->GetPixel(basei);
00481 
00482         --basei[0];
00483         const RealType val001 = this->GetInputImage()->GetPixel(basei);
00484 
00485         const RealType valx01 = val001 + ( val101 - val001 ) * distance0;
00486         const RealType valx11 = val011 + ( val111 - val011 ) * distance0;
00487         const RealType valxx1 = valx01 + ( valx11 - valx01 ) * distance1;
00488 
00489         return ( static_cast< OutputType >( valxx0 + ( valxx1 - valxx0 ) * distance2 ) );
00490         }
00491       }
00492   }
00493 
00494   inline OutputType EvaluateOptimized(const DispatchBase &,
00495                                       const ContinuousIndexType & index) const
00496   {
00497     return this->EvaluateUnoptimized(index);
00498   }
00499 
00500   virtual inline OutputType EvaluateUnoptimized(
00501     const ContinuousIndexType & index) const;
00502 };
00503 } // end namespace itk
00504 
00505 // Define instantiation macro for this template.
00506 #define ITK_TEMPLATE_LinearInterpolateImageFunction(_, EXPORT, TypeX, TypeY)     \
00507   namespace itk                                                                  \
00508   {                                                                              \
00509   _( 2 ( class EXPORT LinearInterpolateImageFunction< ITK_TEMPLATE_2 TypeX > ) ) \
00510   namespace Templates                                                            \
00511   {                                                                              \
00512   typedef LinearInterpolateImageFunction< ITK_TEMPLATE_2 TypeX >                 \
00513   LinearInterpolateImageFunction##TypeY;                                       \
00514   }                                                                              \
00515   }
00516 
00517 #if ITK_TEMPLATE_EXPLICIT
00518 #include "Templates/itkLinearInterpolateImageFunction+-.h"
00519 #endif
00520 
00521 #if ITK_TEMPLATE_TXX
00522 #include "itkLinearInterpolateImageFunction.hxx"
00523 #endif
00524 
00525 #endif
00526