itkMultivariateLegendrePolynomial.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef itkMultivariateLegendrePolynomial_h
#define itkMultivariateLegendrePolynomial_h

#include "itkIntTypes.h"
#include "itkIndent.h"
#include <vector>
#include "itkArray.h"
#include "ITKPolynomialsExport.h"

namespace itk
{
class ITKPolynomials_EXPORT MultivariateLegendrePolynomial
{
public:
  using Self = MultivariateLegendrePolynomial;

  using DoubleArrayType = std::vector< double >;
  using ULongArrayType = std::vector< unsigned long >;
  using LongArrayType = std::vector< long >;

  using CoefficientArrayType = DoubleArrayType;

  using ParametersType = Array< double >;

  using DomainSizeType = ULongArrayType;
  using IndexType = LongArrayType;

  MultivariateLegendrePolynomial(unsigned int dimension,
                                 unsigned int degree,
                                 const DomainSizeType & domainSize);

  virtual ~MultivariateLegendrePolynomial();

  unsigned int GetDimension() const
  { return m_Dimension; }

  unsigned int GetDegree() const
  { return m_Degree; }

  unsigned int GetNumberOfCoefficients() const
  { return m_NumberOfCoefficients; }

  const DomainSizeType & GetDomainSize() const
  { return m_DomainSize; }

  class CoefficientVectorSizeMismatch
  {
public:
    CoefficientVectorSizeMismatch(int given, int required)
    {
      m_Required = required;
      m_Given = given;
    }

    int m_Required;
    int m_Given;
  };

  void SetCoefficients(const CoefficientArrayType & coef);

  void SetCoefficients(const ParametersType & coef);

  const CoefficientArrayType & GetCoefficients() const;

  double Evaluate(IndexType & index)
  {
    if ( m_Dimension == 2 )
      {
      if ( index[1] != m_PrevY )
        {
        // normalized y [-1, 1]
        double norm_y =  m_NormFactor[1]
                        * static_cast< double >( index[1] - 1 );
        this->CalculateXCoef(norm_y, m_CoefficientArray);
        m_PrevY = index[1];
        }

      // normalized x [-1, 1]
      double norm_x =  m_NormFactor[0]
                      * static_cast< double >( index[0] - 1 );

      return LegendreSum(norm_x, m_Degree, m_CachedXCoef);
      }
    else if ( m_Dimension == 3 )
      {
      if ( index[2] != m_PrevZ )
        {
        // normalized z [-1, 1]
        double norm_z =  m_NormFactor[2]
                        * static_cast< double >( index[2] - 1 );
        this->CalculateYCoef(norm_z, m_CoefficientArray);
        m_PrevZ = index[2];
        }

      if ( index[1] != m_PrevY )
        {
        // normalized y [-1, 1]
        double norm_y =  m_NormFactor[1]
                        * static_cast< double >( index[1] - 1 );
        this->CalculateXCoef(norm_y, m_CachedYCoef);
        m_PrevY = index[1];
        }

      // normalized x [-1, 1]
      double norm_x =  m_NormFactor[0]
                      * static_cast< double >( index[0] - 1 );
      return this->LegendreSum(norm_x, m_Degree, m_CachedXCoef);
      }
    return 0;
  }

  unsigned int GetNumberOfCoefficients();

  unsigned int GetNumberOfCoefficients(unsigned int dimension, unsigned int degree);

  class SimpleForwardIterator
  {
public:
    SimpleForwardIterator (MultivariateLegendrePolynomial *polynomial) :
      m_MultivariateLegendrePolynomial(polynomial),
      m_Dimension  (m_MultivariateLegendrePolynomial->GetDimension()),
      m_DomainSize (m_MultivariateLegendrePolynomial->GetDomainSize()),
      m_IsAtEnd(false)
    {
      m_Index.resize(m_Dimension);
      std::fill(m_Index.begin(), m_Index.end(), 0);
    }

    void Begin()
    {
      m_IsAtEnd = false;
      for ( unsigned int dim = 0; dim < m_Dimension; dim++ )
        {
        m_Index[dim] = 0;
        }
    }

    bool IsAtEnd() const
    { return m_IsAtEnd; }

    SimpleForwardIterator & operator++()
    {
      for ( unsigned int dim = 0; dim < m_Dimension; dim++ )
        {
        if ( m_Index[dim] < static_cast< int >( m_DomainSize[dim] - 1 ) )
          {
          m_Index[dim] += 1;
          return *this;
          }
        else
          {
          if ( dim == m_Dimension - 1 )
            {
            m_IsAtEnd = true;
            break;
            }
          else
            {
            m_Index[dim] = 0;
            }
          }
        }
      return *this;
    }

    double Get()
    { return m_MultivariateLegendrePolynomial->Evaluate(m_Index); }

private:
    MultivariateLegendrePolynomial *m_MultivariateLegendrePolynomial;
    unsigned int                    m_Dimension;
    DomainSizeType                  m_DomainSize;
    IndexType                       m_Index;
    bool                            m_IsAtEnd;
  };   // end of class Iterator

  void Print(std::ostream & os);

protected:
  void PrintSelf(std::ostream & os, Indent indent) const;

  double LegendreSum(const double x, int n,
                     const CoefficientArrayType & coef,
                     int offset = 0);

  void CalculateXCoef(double norm_y, const CoefficientArrayType & coef);

  void CalculateYCoef(double norm_z, const CoefficientArrayType & coef);

private:
  DomainSizeType m_DomainSize;
  unsigned int   m_Dimension;
  unsigned int   m_Degree;
  unsigned int   m_NumberOfCoefficients;

  CoefficientArrayType m_CoefficientArray;
  CoefficientArrayType m_CachedXCoef;
  CoefficientArrayType m_CachedYCoef;
  CoefficientArrayType m_CachedZCoef;

  DoubleArrayType m_NormFactor;
  IndexValueType  m_PrevY;
  IndexValueType  m_PrevZ;
}; // end of class

 ITKPolynomials_EXPORT std::ostream & operator<<(std::ostream & os,
                          MultivariateLegendrePolynomial & poly);
} // end of namespace itk
#endif