00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifndef __itkMultivariateLegendrePolynomial_h
00018 #define __itkMultivariateLegendrePolynomial_h
00019
00020 #include "itkArray.h"
00021
00066 namespace itk {
00067
00068 class MultivariateLegendrePolynomial
00069 {
00070 public:
00071
00072
00073
00074 typedef Array< double > DoubleArrayType ;
00075 typedef Array< unsigned long > ULongArrayType ;
00076 typedef Array< long > LongArrayType ;
00077
00079 typedef DoubleArrayType CoefficientArrayType ;
00080
00083 typedef DoubleArrayType ParametersType ;
00084
00086 typedef ULongArrayType DomainSizeType ;
00087 typedef LongArrayType IndexType ;
00088
00090 MultivariateLegendrePolynomial( unsigned int dimension,
00091 unsigned int degree,
00092 const DomainSizeType & domainSize );
00094 virtual ~MultivariateLegendrePolynomial();
00095
00101 void Initialize(void) ;
00102
00104 unsigned int GetDimension(void) const
00105 { return m_Dimension ; }
00106
00108 unsigned int GetDegree(void) const
00109 { return m_Degree ; }
00110
00116 unsigned int GetNumberOfCoefficients(void) const
00117 { return m_NumberOfCoefficients; }
00118
00120 const DomainSizeType & GetDomainSize( void ) const
00121 { return m_DomainSize; }
00122
00124 class CoefficientVectorSizeMismatch
00125 {
00126 public:
00127 CoefficientVectorSizeMismatch(int given, int required)
00128 {
00129 Required = required ;
00130 Given = given ;
00131 }
00132
00133 int Required;
00134 int Given ;
00135 } ;
00136
00141 void SetCoefficients(const CoefficientArrayType & coef)
00142 throw (CoefficientVectorSizeMismatch) ;
00143
00145 const CoefficientArrayType* GetCoefficients(void) const;
00146
00149 double operator() (IndexType index)
00150 {
00151 if (m_Dimension == 2)
00152 {
00153 if (index[1] != m_PrevY)
00154 {
00155
00156 double norm_y = (*m_NormFactor)[1] *
00157 static_cast<double>( index[1] - 1 );
00158 CalculateXCoef(norm_y, m_CoefficientArray) ;
00159 m_PrevY = index[1] ;
00160 }
00161
00162
00163 double norm_x = (*m_NormFactor)[0] *
00164 static_cast<double>( index[0] - 1 );
00165
00166 return LegendreSum(norm_x, m_Degree, m_CachedXCoef) ;
00167 }
00168 else if (m_Dimension == 3)
00169 {
00170 if (m_PrevZ != index[2])
00171 {
00172
00173 double norm_z = (*m_NormFactor)[2] *
00174 static_cast<double>( index[2] - 1 );
00175 CalculateYCoef(norm_z, m_CoefficientArray) ;
00176 m_PrevZ = index[2] ;
00177 }
00178
00179 if (m_PrevY != index[1])
00180 {
00181
00182 double norm_y = (*m_NormFactor)[1] *
00183 static_cast<double>( index[1] - 1 );
00184 CalculateXCoef(norm_y, m_CachedYCoef) ;
00185 m_PrevY = index[1] ;
00186 }
00187
00188
00189 double norm_x = (*m_NormFactor)[0] *
00190 static_cast<double>( index[0] - 1 );
00191 return LegendreSum(norm_x, m_Degree, m_CachedXCoef) ;
00192 }
00193 return 0 ;
00194 }
00195
00197 unsigned int GetNumberOfCoefficients(void);
00198
00200 unsigned int GetNumberOfCoefficients(unsigned int dimension, unsigned int degree) ;
00201
00204 class SimpleForwardIterator
00205 {
00206 public:
00207 SimpleForwardIterator (MultivariateLegendrePolynomial* polynomial)
00208 {
00209 m_MultivariateLegendrePolynomial = polynomial ;
00210 m_Dimension = m_MultivariateLegendrePolynomial->GetDimension();
00211 m_DomainSize = m_MultivariateLegendrePolynomial->GetDomainSize();
00212 m_Index = IndexType(m_Dimension);
00213 }
00214
00215 void Begin( void )
00216 {
00217 m_Index.Fill( 0 ) ;
00218 m_IsAtEnd = false ;
00219 }
00220
00221 bool IsAtEnd()
00222 { return m_IsAtEnd; }
00223
00224 SimpleForwardIterator& operator++()
00225 {
00226 for (unsigned int dim = 0 ; dim < m_Dimension ; dim++)
00227 {
00228 if (m_Index[dim] < static_cast<int>(m_DomainSize[dim] - 1))
00229 {
00230 m_Index[dim] += 1 ;
00231 return *this ;
00232 }
00233 else
00234 {
00235 if (dim == m_Dimension - 1 )
00236 {
00237 m_IsAtEnd = true ;
00238 break ;
00239 }
00240 else
00241 {
00242 m_Index[dim] = 0 ;
00243 }
00244 }
00245 }
00246 return *this ;
00247 }
00248
00249 double Get()
00250 { return (*m_MultivariateLegendrePolynomial)(m_Index); }
00251
00252 private:
00253 MultivariateLegendrePolynomial* m_MultivariateLegendrePolynomial;
00254 unsigned int m_Dimension;
00255 DomainSizeType m_DomainSize;
00256 IndexType m_Index;
00257 bool m_IsAtEnd;
00258 } ;
00259
00260 protected:
00261 double LegendreSum(const double x, int n, const CoefficientArrayType* coef,
00262 int offset = 0);
00263 void CalculateXCoef(double norm_y, const CoefficientArrayType* coef);
00264 void CalculateYCoef(double norm_z, const CoefficientArrayType* coef);
00265
00266 private:
00267 DomainSizeType m_DomainSize;
00268 unsigned int m_Dimension;
00269 unsigned int m_Degree;
00270 unsigned int m_NumberOfCoefficients;
00271 bool m_MultiplicativeBias;
00272
00273
00274 CoefficientArrayType* m_CoefficientArray;
00275 CoefficientArrayType* m_CachedXCoef;
00276 CoefficientArrayType* m_CachedYCoef;
00277 CoefficientArrayType* m_CachedZCoef;
00278 DoubleArrayType* m_NormFactor;
00279 long m_PrevY ;
00280 long m_PrevZ ;
00281 } ;
00282
00283 }
00284 #endif