itkMultivariateLegendrePolynomial.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  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
 *
 *         https://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();
 
  MultivariateLegendrePolynomial(const MultivariateLegendrePolynomial &) = default;
  MultivariateLegendrePolynomial(MultivariateLegendrePolynomial &&) = default;
  MultivariateLegendrePolynomial &
  operator=(const MultivariateLegendrePolynomial &) = default;
  MultivariateLegendrePolynomial &
  operator=(MultivariateLegendrePolynomial &&) = default;

  [[nodiscard]] unsigned int
  GetDimension() const
  {
    return m_Dimension;
  }

  [[nodiscard]] unsigned int
  GetDegree() const
  {
    return m_Degree;
  }

  [[nodiscard]] unsigned int
  GetNumberOfCoefficients() const
  {
    return m_NumberOfCoefficients;
  }

  [[nodiscard]] const DomainSizeType &
  GetDomainSize() const
  {
    return m_DomainSize;
  }

  class CoefficientVectorSizeMismatch
  {
  public:
    CoefficientVectorSizeMismatch(const size_t given, const size_t required)
      : m_Required{ required }
      , m_Given{ given }
    {}
 
    size_t m_Required;
    size_t m_Given;
  };

  void
  SetCoefficients(const CoefficientArrayType & coefficients);
 
  void
  SetCoefficients(const ParametersType & coefficients);

  [[nodiscard]] const CoefficientArrayType &
  GetCoefficients() const;

  double
  Evaluate(IndexType & index)
  {
    if (m_Dimension == 2)
    {
      if (index[1] != m_PrevY)
      {
        // normalized y [-1, 1]
        const 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]
      const double norm_x = m_NormFactor[0] * static_cast<double>(index[0] - 1);
 
      return LegendreSum(norm_x, m_Degree, m_CachedXCoef);
    }
    if (m_Dimension == 3)
    {
      if (index[2] != m_PrevZ)
      {
        // normalized z [-1, 1]
        const 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]
        const 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]
      const 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 dimension, unsigned int degree);

  class SimpleForwardIterator
  {
  public:
    SimpleForwardIterator(MultivariateLegendrePolynomial * polynomial)
      : m_MultivariateLegendrePolynomial(polynomial)
      , m_Dimension(m_MultivariateLegendrePolynomial->GetDimension())
      , m_DomainSize(m_MultivariateLegendrePolynomial->GetDomainSize())
 
    {
      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;
      }
    }
 
    [[nodiscard]] 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;
        }
 
        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{ false };
  }; // end of class Iterator
 
  void
  Print(std::ostream & os) const;
 
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, const MultivariateLegendrePolynomial & poly);
} // end of namespace itk
#endif