00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifndef __itkCoxDeBoorBSplineKernelFunction_h
00018 #define __itkCoxDeBoorBSplineKernelFunction_h
00019
00020 #include "itkKernelFunction.h"
00021 #include "itkNumericTraits.h"
00022 #include "vnl/vnl_math.h"
00023 #include "vnl/vnl_real_polynomial.h"
00024 #include "vnl/vnl_matrix.h"
00025
00026 namespace itk
00027 {
00028
00056 template <unsigned int VSplineOrder = 3>
00057 class ITK_EXPORT CoxDeBoorBSplineKernelFunction
00058 : public KernelFunction
00059 {
00060 public:
00062 typedef CoxDeBoorBSplineKernelFunction Self;
00063 typedef KernelFunction Superclass;
00064 typedef SmartPointer<Self> Pointer;
00065 typedef SmartPointer<const Self> ConstPointer;
00066
00068 itkNewMacro(Self);
00069
00071 itkTypeMacro( CoxDeBoorBSplineKernelFunction, KernelFunction );
00072
00073 typedef double RealType;
00074 typedef vnl_vector<RealType> VectorType;
00075 typedef vnl_real_polynomial PolynomialType;
00076 typedef vnl_matrix<RealType> MatrixType;
00077
00079 void SetSplineOrder( unsigned int );
00080 itkGetConstMacro( SplineOrder, unsigned int );
00082
00084 inline RealType Evaluate( const RealType & u ) const
00085 {
00086 RealType absValue = vnl_math_abs( u );
00087 unsigned int which;
00088 if( this->m_SplineOrder % 2 == 0 )
00089 {
00090 which = static_cast<unsigned int>( absValue+0.5 );
00091 }
00092 else
00093 {
00094 which = static_cast<unsigned int>( absValue );
00095 }
00096 if( which < this->m_BSplineShapeFunctions.rows() )
00097 {
00098 return PolynomialType(
00099 this->m_BSplineShapeFunctions.get_row( which ) ).evaluate( absValue );
00100 }
00101 else
00102 {
00103 return NumericTraits<RealType>::Zero;
00104 }
00105 }
00107
00109 inline RealType EvaluateDerivative( const double & u ) const
00110 {
00111 return this->EvaluateNthDerivative( u, 1 );
00112 }
00113
00115 inline RealType EvaluateNthDerivative( const double & u, unsigned int n ) const
00116 {
00117 RealType absValue = vnl_math_abs( u );
00118 unsigned int which;
00119 if( this->m_SplineOrder % 2 == 0 )
00120 {
00121 which = static_cast<unsigned int>( absValue+0.5 );
00122 }
00123 else
00124 {
00125 which = static_cast<unsigned int>( absValue );
00126 }
00127 if( which < this->m_BSplineShapeFunctions.rows() )
00128 {
00129 PolynomialType polynomial(
00130 this->m_BSplineShapeFunctions.get_row( which ) );
00131 for( unsigned int i = 0; i < n; i++ )
00132 {
00133 polynomial = polynomial.derivative();
00134 }
00135 RealType der = polynomial.evaluate( absValue );
00136 if( u < NumericTraits<RealType>::Zero && n % 2 != 0 )
00137 {
00138 return -der;
00139 }
00140 else
00141 {
00142 return der;
00143 }
00144 }
00145 else
00146 {
00147 return NumericTraits<RealType>::Zero;
00148 }
00149 }
00151
00157 MatrixType GetShapeFunctions()
00158 {
00159 return this->m_BSplineShapeFunctions;
00160 }
00161
00166 MatrixType GetShapeFunctionsInZeroToOneInterval();
00167
00168 protected:
00169 CoxDeBoorBSplineKernelFunction();
00170 ~CoxDeBoorBSplineKernelFunction();
00171 void PrintSelf( std::ostream& os, Indent indent ) const;
00172
00173 private:
00174 CoxDeBoorBSplineKernelFunction( const Self& );
00175 void operator=( const Self& );
00176
00181 void GenerateBSplineShapeFunctions( unsigned int );
00182
00189 PolynomialType CoxDeBoor( unsigned short,
00190 VectorType, unsigned int, unsigned int );
00191
00192 MatrixType m_BSplineShapeFunctions;
00193 unsigned int m_SplineOrder;
00194
00195 };
00196
00197 }
00198
00199 #ifndef ITK_MANUAL_INSTANTIATION
00200 #include "itkCoxDeBoorBSplineKernelFunction.txx"
00201 #endif
00202
00203 #endif
00204