itkQuadEdgeMeshDecimationQuadricElementHelper.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 itkQuadEdgeMeshDecimationQuadricElementHelper_h
#define itkQuadEdgeMeshDecimationQuadricElementHelper_h
 
#include "itkPoint.h"
#include "itkMathSVD.h"
#include "vnl/vnl_vector_fixed.h"
#include "vnl/vnl_matrix.h"
#include <algorithm>
 
#include "itkTriangleHelper.h"
 
namespace itk
{
template <typename TPoint>
class QuadEdgeMeshDecimationQuadricElementHelper
{
public:
  using Self = QuadEdgeMeshDecimationQuadricElementHelper;
  using PointType = TPoint;
  using CoordType = typename PointType::CoordinateType;
 
  static constexpr unsigned int PointDimension = PointType::PointDimension;
  static constexpr unsigned int NumberOfCoefficients = PointDimension * (PointDimension + 1 / 2 + PointDimension + 1);
 
  using VectorType = typename PointType::VectorType;
  using VNLMatrixType = vnl_matrix<CoordType>;
  using VNLVectorType = vnl_vector_fixed<CoordType, Self::PointDimension>;
  using CoefficientVectorType = vnl_vector_fixed<CoordType, Self::NumberOfCoefficients>;
  using TriangleType = TriangleHelper<PointType>;
 
  // *****************************************************************
  QuadEdgeMeshDecimationQuadricElementHelper()
    : m_Coefficients(CoordType{})
    , m_A(PointDimension, PointDimension, CoordType{})
    , m_B(CoordType{})
    , m_Rank(PointDimension)
    , m_SVDAbsoluteThreshold(static_cast<CoordType>(1e-6))
    , m_SVDRelativeThreshold(static_cast<CoordType>(1e-3))
  {}
 
  QuadEdgeMeshDecimationQuadricElementHelper(const CoefficientVectorType & iCoefficients)
    : m_Coefficients(iCoefficients)
    , m_A(PointDimension, PointDimension, CoordType{})
    , m_B(CoordType{})
    , m_Rank(PointDimension)
    , m_SVDAbsoluteThreshold(static_cast<CoordType>(1e-3))
    , m_SVDRelativeThreshold(static_cast<CoordType>(1e-3))
  {
    this->ComputeAMatrixAndBVector();
  }
 
  ~QuadEdgeMeshDecimationQuadricElementHelper() = default;
 
  [[nodiscard]] CoefficientVectorType
  GetCoefficients() const
  {
    return this->m_Coefficients;
  }
 
  VNLMatrixType
  GetAMatrix()
  {
    this->ComputeAMatrixAndBVector();
    return m_A;
  }
 
  VNLVectorType
  GetBVector()
  {
    ComputeAMatrixAndBVector();
    return m_B;
  }
 
  [[nodiscard]] unsigned int
  GetRank() const
  {
    return m_Rank;
  }

  [[nodiscard]] inline CoordType
  ComputeError(const PointType & iP) const
  {
    //     ComputeAMatrixAndBVector();
    const auto      svd = itk::Math::SVD(m_A, /*canonicalizeSigns=*/false);
    const CoordType wmax = svd.W[0];
    // Fold the absolute and relative singular-value tolerances into one rcond.
    const CoordType rcond =
      (wmax > CoordType{}) ? std::max(m_SVDAbsoluteThreshold / wmax, m_SVDRelativeThreshold) : m_SVDRelativeThreshold;
    const CoordType oError = inner_product(iP.GetVnlVector(), svd.Recompose(rcond) * iP.GetVnlVector());
 
    return this->m_Coefficients.back() - oError;
    /*
    CoordType oError( 0. );
 
    std::vector< CoordType > pt( PointDimension + 1, 1. );
 
    unsigned int dim1( 0 ), dim2, k( 0 );
 
    while( dim1 < PointDimension )
      {
      pt[dim1] = iP[dim1];
      ++dim1;
      }
 
    for( dim1 = 0; dim1 < PointDimension + 1; ++dim1 )
      {
      oError += this->m_Coefficients[k++] * pt[dim1] * pt[dim1];
 
      for( dim2 = dim1 + 1; dim2 < PointDimension + 1; ++dim2 )
        {
        oError += 2. * this->m_Coefficients[k++] * pt[dim1] * pt[dim2];
        }
      }
    oError += this->m_Coefficients[k++];
 
    return oError;*/
  }

  inline CoordType
  ComputeErrorAtOptimalLocation(const PointType & iP)
  {
    PointType optimal_location = ComputeOptimalLocation(iP);
 
    return ComputeError(optimal_location);
  }
 
  PointType
  ComputeOptimalLocation(const PointType & iP)
  {
    ComputeAMatrixAndBVector();
 
    const auto      svd = itk::Math::SVD(m_A, /*canonicalizeSigns=*/false);
    const CoordType wmax = svd.W[0];
    const CoordType rcond =
      (wmax > CoordType{}) ? std::max(m_SVDAbsoluteThreshold / wmax, m_SVDRelativeThreshold) : m_SVDRelativeThreshold;
 
    m_Rank = svd.Rank(rcond);
 
    const auto y = (m_B.as_vector() - m_A * iP.GetVnlVector());
 
    VNLVectorType displacement = svd.Solve(y, rcond);
 
    PointType oP;
    for (unsigned int dim = 0; dim < PointDimension; ++dim)
    {
      oP[dim] = iP[dim] + displacement[dim];
    }
    return oP;
  }

  PointType
  ComputeOptimalLocation(const unsigned int)
  {}
 
  void
  AddTriangle(const PointType & iP1,
              const PointType & iP2,
              const PointType & iP3,
              const CoordType & iWeight = static_cast<CoordType>(1.))
  {
    const VectorType N = TriangleType::ComputeNormal(iP1, iP2, iP3);
 
    AddPoint(iP1, N, iWeight);
  }
 
  void
  AddPoint(const PointType & iP, const VectorType & iN, const CoordType & iWeight = static_cast<CoordType>(1.))
  {
    unsigned int    k(0);
    const CoordType d = -iN * iP.GetVectorFromOrigin();
    for (unsigned int dim1 = 0; dim1 < PointDimension; ++dim1)
    {
      for (unsigned int dim2 = dim1; dim2 < PointDimension; ++dim2)
      {
        this->m_Coefficients[k++] += iWeight * iN[dim1] * iN[dim2];
      }
      this->m_Coefficients[k++] += iWeight * iN[dim1] * d;
    }
 
    this->m_Coefficients[k++] += iWeight * d * d;
  }
 
  // ***********************************************************************
  // operators
  Self &
  operator=(const Self & iRight)
  {
    if (this != &iRight)
    {
      this->m_Coefficients = iRight.m_Coefficients;
    }
    return *this;
  }
 
  Self
  operator+(const Self & iRight) const
  {
    return Self(this->m_Coefficients + iRight.m_Coefficients);
  }
 
  Self &
  operator+=(const Self & iRight)
  {
    this->m_Coefficients += iRight.m_Coefficients;
    return *this;
  }
 
  Self
  operator-(const Self & iRight) const
  {
    return Self(this->m_Coefficients - iRight.m_Coefficients);
  }
 
  Self &
  operator-=(const Self & iRight)
  {
    this->m_Coefficients -= iRight.m_Coefficients;
    return *this;
  }
 
  Self
  operator*(const CoordType & iV) const
  {
    Self oElement(this->m_Coefficients * iV);
 
    return oElement;
  }
 
  Self &
  operator*=(const CoordType & iV)
  {
    this->m_Coefficients *= iV;
    return *this;
  }
 
protected:
  CoefficientVectorType m_Coefficients;
  VNLMatrixType         m_A;
  VNLVectorType         m_B;
  unsigned int          m_Rank;
  CoordType             m_SVDAbsoluteThreshold;
  CoordType             m_SVDRelativeThreshold;
  // bool                        m_MatrixFilled;
 
  void
  ComputeAMatrixAndBVector()
  {
    unsigned int k(0);
    for (unsigned int dim1 = 0; dim1 < PointDimension; ++dim1)
    {
      for (unsigned int dim2 = dim1; dim2 < PointDimension; ++dim2)
      {
        m_A[dim1][dim2] = m_A[dim2][dim1] = m_Coefficients[k++];
      }
      m_B[dim1] = -m_Coefficients[k++];
    }
    // m_MatrixFilled = true;
  }
};
} // namespace itk
#endif