ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkMultivariateLegendrePolynomial.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 __itkMultivariateLegendrePolynomial_h
00019 #define __itkMultivariateLegendrePolynomial_h
00020 
00021 #include "itkIntTypes.h"
00022 #include "itkIndent.h"
00023 #include <vector>
00024 #include "itkArray.h"
00025 
00026 namespace itk
00027 {
00074 class ITK_EXPORT MultivariateLegendrePolynomial
00075 {
00076 public:
00077   typedef MultivariateLegendrePolynomial Self;
00078 
00079   typedef std::vector< double >        DoubleArrayType;
00080   typedef std::vector< unsigned long > ULongArrayType;
00081   typedef std::vector< long >          LongArrayType;
00082 
00084   typedef DoubleArrayType CoefficientArrayType;
00085 
00088   typedef Array< double > ParametersType;
00089 
00091   typedef ULongArrayType DomainSizeType;
00092   typedef LongArrayType  IndexType;
00093 
00095   MultivariateLegendrePolynomial(unsigned int dimension,
00096                                  unsigned int degree,
00097                                  const DomainSizeType & domainSize);
00098 
00100   virtual ~MultivariateLegendrePolynomial();
00101 
00103   unsigned int GetDimension(void) const
00104   { return m_Dimension; }
00105 
00107   unsigned int GetDegree(void) const
00108   { return m_Degree; }
00109 
00116   unsigned int GetNumberOfCoefficients(void) const
00117   { return m_NumberOfCoefficients; }
00118 
00120   const DomainSizeType & GetDomainSize(void) const
00121   { return m_DomainSize; }
00122 
00127   class CoefficientVectorSizeMismatch
00128   {
00129 public:
00130     CoefficientVectorSizeMismatch(int given, int required)
00131     {
00132       m_Required = required;
00133       m_Given = given;
00134     }
00135 
00136     int m_Required;
00137     int m_Given;
00138   };
00139 
00144   void SetCoefficients(const CoefficientArrayType & coef)
00145   throw ( CoefficientVectorSizeMismatch );
00146 
00147   void SetCoefficients(const ParametersType & coef)
00148   throw ( CoefficientVectorSizeMismatch );
00149 
00151   const CoefficientArrayType & GetCoefficients(void) const;
00152 
00155   double Evaluate(IndexType & index)
00156   {
00157     if ( m_Dimension == 2 )
00158       {
00159       if ( index[1] != m_PrevY )
00160         {
00161         // normalized y [-1, 1]
00162         double norm_y =  m_NormFactor[1]
00163                         * static_cast< double >( index[1] - 1 );
00164         this->CalculateXCoef(norm_y, m_CoefficientArray);
00165         m_PrevY = index[1];
00166         }
00168 
00169       // normalized x [-1, 1]
00170       double norm_x =  m_NormFactor[0]
00171                       * static_cast< double >( index[0] - 1 );
00172 
00173       return LegendreSum(norm_x, m_Degree, m_CachedXCoef);
00174       }
00175     else if ( m_Dimension == 3 )
00176       {
00177       if ( index[2] != m_PrevZ )
00178         {
00179         // normalized z [-1, 1]
00180         double norm_z =  m_NormFactor[2]
00181                         * static_cast< double >( index[2] - 1 );
00182         this->CalculateYCoef(norm_z, m_CoefficientArray);
00183         m_PrevZ = index[2];
00184         }
00185 
00186       if ( index[1] != m_PrevY )
00187         {
00188         // normalized y [-1, 1]
00189         double norm_y =  m_NormFactor[1]
00190                         * static_cast< double >( index[1] - 1 );
00191         this->CalculateXCoef(norm_y, m_CachedYCoef);
00192         m_PrevY = index[1];
00193         }
00194 
00195       // normalized x [-1, 1]
00196       double norm_x =  m_NormFactor[0]
00197                       * static_cast< double >( index[0] - 1 );
00198       return this->LegendreSum(norm_x, m_Degree, m_CachedXCoef);
00199       }
00200     return 0;
00201   }
00202 
00204   unsigned int GetNumberOfCoefficients();
00205 
00207   unsigned int GetNumberOfCoefficients(unsigned int dimension, unsigned int degree);
00208 
00215   class SimpleForwardIterator
00216   {
00217 public:
00218     SimpleForwardIterator (MultivariateLegendrePolynomial *polynomial)
00219     {
00220       m_MultivariateLegendrePolynomial = polynomial;
00221       m_Dimension   = m_MultivariateLegendrePolynomial->GetDimension();
00222       m_DomainSize  = m_MultivariateLegendrePolynomial->GetDomainSize();
00223       m_Index.resize(m_Dimension);
00224       std::fill(m_Index.begin(), m_Index.end(), 0);
00225     }
00226 
00227     void Begin(void)
00228     {
00229       m_IsAtEnd = false;
00230       for ( unsigned int dim = 0; dim < m_Dimension; dim++ )
00231         {
00232         m_Index[dim] = 0;
00233         }
00234     }
00235 
00236     bool IsAtEnd()
00237     { return m_IsAtEnd; }
00238 
00239     SimpleForwardIterator & operator++()
00240     {
00241       for ( unsigned int dim = 0; dim < m_Dimension; dim++ )
00242         {
00243         if ( m_Index[dim] < static_cast< int >( m_DomainSize[dim] - 1 ) )
00244           {
00245           m_Index[dim] += 1;
00246           return *this;
00247           }
00248         else
00249           {
00250           if ( dim == m_Dimension - 1 )
00251             {
00252             m_IsAtEnd = true;
00253             break;
00254             }
00255           else
00256             {
00257             m_Index[dim] = 0;
00258             }
00259           }
00260         }
00261       return *this;
00262     }
00263 
00264     double Get()
00265     { return m_MultivariateLegendrePolynomial->Evaluate(m_Index); }
00266 private:
00267     MultivariateLegendrePolynomial *m_MultivariateLegendrePolynomial;
00268     unsigned int                    m_Dimension;
00269     DomainSizeType                  m_DomainSize;
00270     IndexType                       m_Index;
00271     bool                            m_IsAtEnd;
00272   };   // end of class Iterator
00273 
00274   void Print(std::ostream & os);
00275 
00276 protected:
00277   void PrintSelf(std::ostream & os, Indent indent) const;
00278 
00279   double LegendreSum(const double x, int n,
00280                      const CoefficientArrayType & coef,
00281                      int offset = 0);
00282 
00283   void CalculateXCoef(double norm_y, const CoefficientArrayType & coef);
00284 
00285   void CalculateYCoef(double norm_z, const CoefficientArrayType & coef);
00286 
00287 private:
00288   DomainSizeType m_DomainSize;
00289   unsigned int   m_Dimension;
00290   unsigned int   m_Degree;
00291   unsigned int   m_NumberOfCoefficients;
00292   bool           m_MultiplicativeBias;
00293 
00294   CoefficientArrayType m_CoefficientArray;
00295   CoefficientArrayType m_CachedXCoef;
00296   CoefficientArrayType m_CachedYCoef;
00297   CoefficientArrayType m_CachedZCoef;
00298 
00299   DoubleArrayType m_NormFactor;
00300   IndexValueType  m_PrevY;
00301   IndexValueType  m_PrevZ;
00302 }; // end of class
00303 
00304 ITK_EXPORT std::ostream & operator<<(std::ostream & os,
00305                           MultivariateLegendrePolynomial & poly);
00306 } // end of namespace itk
00307 #endif
00308