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 "itkIndent.h"
00021 #include <vector>
00022 #include "itkArray.h"
00023
00024 namespace itk {
00025
00070 class MultivariateLegendrePolynomial
00071 {
00072 public:
00073 typedef MultivariateLegendrePolynomial Self ;
00074
00075 typedef std::vector< double > DoubleArrayType ;
00076 typedef std::vector< unsigned long > ULongArrayType ;
00077 typedef std::vector< long > LongArrayType ;
00078
00080 typedef DoubleArrayType CoefficientArrayType ;
00081
00084 typedef Array< double > ParametersType ;
00085
00087 typedef ULongArrayType DomainSizeType ;
00088 typedef LongArrayType IndexType ;
00089
00091 MultivariateLegendrePolynomial( unsigned int dimension,
00092 unsigned int degree,
00093 const DomainSizeType & domainSize );
00094
00096 virtual ~MultivariateLegendrePolynomial();
00097
00099 unsigned int GetDimension(void) const
00100 { return m_Dimension ; }
00101
00103 unsigned int GetDegree(void) const
00104 { return m_Degree ; }
00105
00111 unsigned int GetNumberOfCoefficients(void) const
00112 { return m_NumberOfCoefficients; }
00113
00115 const DomainSizeType & GetDomainSize( void ) const
00116 { return m_DomainSize; }
00117
00119 class CoefficientVectorSizeMismatch
00120 {
00121 public:
00122 CoefficientVectorSizeMismatch(int given, int required)
00123 {
00124 Required = required ;
00125 Given = given ;
00126 }
00127
00128 int Required;
00129 int Given ;
00130 } ;
00131
00136 void SetCoefficients(const CoefficientArrayType& coef)
00137 throw (CoefficientVectorSizeMismatch) ;
00138
00139 void SetCoefficients(const ParametersType& coef)
00140 throw (CoefficientVectorSizeMismatch) ;
00141
00143 const CoefficientArrayType& GetCoefficients(void) const;
00144
00147 double Evaluate(IndexType& index)
00148 {
00149 if (m_Dimension == 2)
00150 {
00151 if (index[1] != m_PrevY)
00152 {
00153
00154 double norm_y = m_NormFactor[1] *
00155 static_cast<double>( index[1] - 1 );
00156 this->CalculateXCoef(norm_y, m_CoefficientArray) ;
00157 m_PrevY = index[1] ;
00158 }
00160
00161
00162 double norm_x = m_NormFactor[0] *
00163 static_cast<double>( index[0] - 1 );
00164
00165 return LegendreSum(norm_x, m_Degree, m_CachedXCoef) ;
00166 }
00167 else if (m_Dimension == 3)
00168 {
00169 if (index[2] != m_PrevZ )
00170 {
00171
00172 double norm_z = m_NormFactor[2] *
00173 static_cast<double>( index[2] - 1 );
00174 this->CalculateYCoef(norm_z, m_CoefficientArray) ;
00175 m_PrevZ = index[2] ;
00176 }
00177
00178 if (index[1] != m_PrevY)
00179 {
00180
00181 double norm_y = m_NormFactor[1] *
00182 static_cast<double>( index[1] - 1 );
00183 this->CalculateXCoef(norm_y, m_CachedYCoef) ;
00184 m_PrevY = index[1] ;
00185 }
00186
00187
00188 double norm_x = m_NormFactor[0] *
00189 static_cast<double>( index[0] - 1 );
00190 return this->LegendreSum(norm_x, m_Degree, m_CachedXCoef) ;
00191 }
00192 return 0 ;
00193 }
00194
00196 unsigned int GetNumberOfCoefficients();
00197
00199 unsigned int GetNumberOfCoefficients(unsigned int dimension, unsigned int degree) ;
00200
00203 class SimpleForwardIterator
00204 {
00205 public:
00206 SimpleForwardIterator (MultivariateLegendrePolynomial* polynomial)
00207 {
00208 m_MultivariateLegendrePolynomial = polynomial ;
00209 m_Dimension = m_MultivariateLegendrePolynomial->GetDimension() ;
00210 m_DomainSize = m_MultivariateLegendrePolynomial->GetDomainSize() ;
00211 m_Index.resize(m_Dimension) ;
00212 std::fill(m_Index.begin(), m_Index.end(), 0) ;
00213 }
00215
00216 void Begin( void )
00217 {
00218 m_IsAtEnd = false ;
00219 for (unsigned int dim = 0 ; dim < m_Dimension ; dim++)
00220 {
00221 m_Index[dim] = 0 ;
00222 }
00223 }
00224
00225 bool IsAtEnd()
00226 { return m_IsAtEnd; }
00227
00228 SimpleForwardIterator& operator++()
00229 {
00230 for (unsigned int dim = 0 ; dim < m_Dimension ; dim++)
00231 {
00232 if (m_Index[dim] < static_cast<int>(m_DomainSize[dim] - 1))
00233 {
00234 m_Index[dim] += 1 ;
00235 return *this ;
00236 }
00237 else
00238 {
00239 if (dim == m_Dimension - 1 )
00240 {
00241 m_IsAtEnd = true ;
00242 break ;
00243 }
00244 else
00245 {
00246 m_Index[dim] = 0 ;
00247 }
00248 }
00249 }
00250 return *this ;
00251 }
00252
00253 double Get()
00254 { return m_MultivariateLegendrePolynomial->Evaluate(m_Index); }
00255
00256 private:
00257 MultivariateLegendrePolynomial* m_MultivariateLegendrePolynomial;
00258 unsigned int m_Dimension;
00259 DomainSizeType m_DomainSize;
00260 IndexType m_Index;
00261 bool m_IsAtEnd;
00262 } ;
00263
00264 void Print(std::ostream& os) ;
00265
00266 protected:
00267 void PrintSelf(std::ostream& os, Indent indent) const;
00268 double LegendreSum(const double x, int n, const CoefficientArrayType& coef,
00269 int offset = 0);
00270 void CalculateXCoef(double norm_y, const CoefficientArrayType& coef);
00271 void CalculateYCoef(double norm_z, const CoefficientArrayType& coef);
00272
00273 private:
00274 DomainSizeType m_DomainSize;
00275 unsigned int m_Dimension;
00276 unsigned int m_Degree;
00277 unsigned int m_NumberOfCoefficients;
00278 bool m_MultiplicativeBias;
00279
00280 CoefficientArrayType m_CoefficientArray;
00281 CoefficientArrayType m_CachedXCoef;
00282 CoefficientArrayType m_CachedYCoef;
00283 CoefficientArrayType m_CachedZCoef;
00284 DoubleArrayType m_NormFactor;
00285 long m_PrevY ;
00286 long m_PrevZ ;
00287 } ;
00288
00289 std::ostream& operator<< (std::ostream& os,
00290 MultivariateLegendrePolynomial& poly) ;
00291 }
00292 #endif
00293