SourceAccessorCubicInterpolate.h

Go to the documentation of this file.
00001 // $Id$
00002 
00003 #ifndef __SOURCEACCESSORCUBICINTERPOLATE_H__
00004 #define __SOURCEACCESSORCUBICINTERPOLATE_H__
00005 
00006 #include "SourceGenerics.h"
00007 
00008 #include "Promotion.h"
00009 #include <math.h>
00010 
00011 bool DoSpecialDebug();
00012 
00013 //============================================================================
00014 // forward declarations of interpolation functions
00015 //============================================================================
00016 template <class DATA>
00017 inline PROMOTION(Real,DATA) CubicInterpolate (
00018     const Real& x,
00019     const DATA& data0, const DATA& data1, const DATA& data2, const DATA& data3
00020 );
00021 
00022 //============================================================================
00023 // base class
00024 //============================================================================
00025 template<class DATA, int DIMENSIONALITY, class PRECISION, class SOURCE>
00026 class SourceAccessorCubicInterpolateOf: public SourceAccessorOf<DATA,
00027     DIMENSIONALITY, PRECISION, SOURCE>
00028 {
00029 };
00030 
00031 //============================================================================
00032 // 2D
00033 //============================================================================
00034 template<class DATA, class PRECISION, class SOURCE>
00035 class SourceAccessorCubicInterpolateOf<DATA, 2, PRECISION, SOURCE> : public SourceAccessorOf<
00036     DATA, 2, PRECISION, SOURCE>
00037 {
00038 public:
00039 SOURCE_ACTUALS_2D  ;
00040 
00041   virtual String Describe() const
00042     {
00043     return (
00044         _LINE + "SourceAccessorCubicInterpolateOf (2D)" + LINE_ +
00045         this->DescribeCommon()
00046     );
00047     }
00048 
00049 private:
00050   inline PROMOTION(DATA,PRECISION) InterpolateOnY(
00051       const PRECISION &rX, const PRECISION &rY
00052   ) const
00053     {
00054     // note center is at pxl11, not pxl00
00055     DATA data0; ::Get(*(this->m_psource), data0, (PRECISION)((int)(rX-1)),(PRECISION)((int)(rY)));
00056     DATA data1; ::Get(*(this->m_psource), data1, (PRECISION)((int)(rX+0)),(PRECISION)((int)(rY)));
00057     DATA data2; ::Get(*(this->m_psource), data2, (PRECISION)((int)(rX+1)),(PRECISION)((int)(rY)));
00058     DATA data3; ::Get(*(this->m_psource), data3, (PRECISION)((int)(rX+2)),(PRECISION)((int)(rY)));
00059 
00060     return CubicInterpolate(rX-int(rX), data0,data1,data2,data3);
00061     }
00062 
00063 public:
00064 
00065   inline void Get(DATA& dataOut, const PRECISION &rX, const PRECISION &rY) const
00066     {
00067     PROMOTION(DATA,PRECISION) data0 = InterpolateOnY((PRECISION)((int)(rX)),(PRECISION)((int)(rY-1)));
00068     PROMOTION(DATA,PRECISION) data1 = InterpolateOnY((PRECISION)((int)(rX)),(PRECISION)((int)(rY)) );
00069     PROMOTION(DATA,PRECISION) data2 = InterpolateOnY((PRECISION)((int)(rX)),(PRECISION)((int)(rY+1)));
00070     PROMOTION(DATA,PRECISION) data3 = InterpolateOnY((PRECISION)((int)(rX)),(PRECISION)((int)(rY+2)));
00071 
00072     dataOut=CubicInterpolate(rY-int(rY), data0,data1,data2,data3);
00073     }
00074 
00075   void Set(const PRECISION &rX, const PRECISION &rY, const DATA& data)
00076     {
00077     ERROR("NOT IMPLEMENTED");
00078     }
00079 
00080   };
00081 
00082 //============================================================================
00083 // 3d
00084 //============================================================================
00085 template<class DATA, class PRECISION, class SOURCE>
00086 class SourceAccessorCubicInterpolateOf<DATA, 3, PRECISION, SOURCE> : public SourceAccessorOf<
00087     DATA, 3, PRECISION, SOURCE>
00088 {
00089 public:
00090 SOURCE_ACTUALS_3D  ;
00091 
00092   virtual String Describe() const
00093     {
00094     return (
00095         _LINE + "SourceAccessorCubicInterpolateOf (3D)" + LINE_ +
00096         this->DescribeCommon()
00097     );
00098     }
00099 
00100 private:
00101   inline PROMOTION(DATA,Real) InterpolateOnYZ(
00102       const PRECISION &rX, const PRECISION &rY, const PRECISION &rZ
00103   ) const
00104     {
00105     // note center is at pxl11, not pxl00
00106     DATA data0; ::Get( *(this->m_psource), data0, (PRECISION)((int)(rX-1)), (PRECISION)((int)(rY)), (PRECISION)((int)(rZ)));
00107     DATA data1; ::Get( *(this->m_psource), data1, (PRECISION)((int)(rX+0)), (PRECISION)((int)(rY)), (PRECISION)((int)(rZ)));
00108     DATA data2; ::Get( *(this->m_psource), data2, (PRECISION)((int)(rX+1)), (PRECISION)((int)(rY)), (PRECISION)((int)(rZ)));
00109     DATA data3; ::Get( *(this->m_psource), data3, (PRECISION)((int)(rX+2)), (PRECISION)((int)(rY)), (PRECISION)((int)(rZ)));
00110 
00111     return CubicInterpolate(rX-int(rX), data0,data1,data2,data3);
00112     }
00113 
00114   inline PROMOTION(DATA,Real) InterpolateOnZ(
00115       const PRECISION &rX, const PRECISION &rY, const PRECISION &rZ
00116   ) const
00117     {
00118     PROMOTION(DATA,Real) data0 = InterpolateOnYZ(rX,rY-1,rZ);
00119     PROMOTION(DATA,Real) data1 = InterpolateOnYZ(rX,rY ,rZ);
00120     PROMOTION(DATA,Real) data2 = InterpolateOnYZ(rX,rY+1,rZ);
00121     PROMOTION(DATA,Real) data3 = InterpolateOnYZ(rX,rY+2,rZ);
00122 
00123     return CubicInterpolate(rY-int(rY), data0,data1,data2,data3);
00124     }
00125 
00126 public:
00127   inline void Get(DATA& dataOut, const PRECISION &rX, const PRECISION &rY, const PRECISION &rZ) const
00128     {
00129     PROMOTION(DATA,Real) data0 = InterpolateOnZ(rX,rY,rZ-1);
00130     PROMOTION(DATA,Real) data1 = InterpolateOnZ(rX,rY,rZ );
00131     PROMOTION(DATA,Real) data2 = InterpolateOnZ(rX,rY,rZ+1);
00132     PROMOTION(DATA,Real) data3 = InterpolateOnZ(rX,rY,rZ+2);
00133 
00134     dataOut= CubicInterpolate(rZ-int(rZ), data0, data1, data2, data3);
00135     }
00136 
00137   void Set(const PRECISION &rX, const PRECISION &rY, const PRECISION &rZ, const DATA& data)
00138     {
00139     ERROR("NOT IMPLEMENTED");
00140     }
00141   };
00142 
00143 // using formulae from
00144 //  http://en.wikipedia.org/wiki/Cubic_interpolation#Interpolation_on_the_unit_interval_without_exact_derivatives
00145 // computes cubic interpolation of 4 points at point x,y
00146 // note: variable indicies off by +1, so centered at data1 == f(x=0)
00147 template <class DATA>
00148 inline PROMOTION(Real,DATA) CubicInterpolate (
00149     const Real& x,
00150     const DATA& data0, const DATA& data1, const DATA& data2, const DATA& data3
00151 ) __attribute__((always_inline));
00152 
00153 template <class DATA>
00154 inline PROMOTION(Real,DATA) CubicInterpolate (
00155     const Real& x,
00156     const DATA& data0, const DATA& data1, const DATA& data2, const DATA& data3
00157 )
00158   {
00159   const PROMOTION(Real,DATA) x2 = x * x;
00160 
00161   return Real(.5) * (
00162       + (x*((2-x)*x-1)) * data0
00163       + (x2*(3*x-5)+2) * data1
00164       + (x*((4-3*x)*x+1)) * data2
00165       + ((x-1)*x2) * data3);
00166 
00167   }
00168 
00169 MAKE_INTERPOLANT_ASSISTANT( CubicInterpolation, SourceAccessorCubicInterpolateOf );
00170 
00171 #endif // __SOURCEACCESSORCUBICINTERPOLATE_H__

Generated on 6 Apr 2011 for Slicer3 by  doxygen 1.6.1