itkQuadEdgeMeshParamMatrixCoefficients.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 itkQuadEdgeMeshParamMatrixCoefficients_h
#define itkQuadEdgeMeshParamMatrixCoefficients_h

#include "itkQuadEdgeMesh.h"
#include "itkTriangleHelper.h"
#include "itkMath.h"

namespace itk
{
template< typename TInputMesh >
class MatrixCoefficients
{
public:
  using InputMeshType = TInputMesh;
  using InputCoordRepType = typename InputMeshType::CoordRepType;
  using InputQEType = typename InputMeshType::QEType;

  MatrixCoefficients()= default;
  virtual ~MatrixCoefficients() = default;

  virtual InputCoordRepType operator()
    (const InputMeshType *iMesh, InputQEType *iEdge) const = 0;
};

template< typename TInputMesh >
class OnesMatrixCoefficients:public MatrixCoefficients< TInputMesh >
{
public:
  using Superclass = MatrixCoefficients< TInputMesh >;

  using InputMeshType = TInputMesh;
  using InputCoordRepType = typename InputMeshType::CoordRepType;
  using InputQEType = typename InputMeshType::QEType;

  OnesMatrixCoefficients() = default;

  InputCoordRepType operator()( const InputMeshType *itkNotUsed(iMesh),
                                InputQEType *itkNotUsed(iEdge) ) const override
  {
    return 1.0;
  }
};

template< typename TInputMesh >
class InverseEuclideanDistanceMatrixCoefficients:
  public MatrixCoefficients< TInputMesh >
{
public:
  using Superclass = MatrixCoefficients< TInputMesh >;

  using InputMeshType = TInputMesh;
  using InputCoordRepType = typename InputMeshType::CoordRepType;
  using InputPointType = typename InputMeshType::PointType;
  using InputPointIdentifier = typename InputMeshType::PointIdentifier;
  using InputQEType = typename InputMeshType::QEType;
  using InputVectorType = typename InputMeshType::VectorType;

  InverseEuclideanDistanceMatrixCoefficients() = default;

  InputCoordRepType operator()(const InputMeshType *iMesh, InputQEType *iEdge) const override
  {
    InputPointIdentifier id1 = iEdge->GetOrigin();
    InputPointIdentifier id2 = iEdge->GetDestination();

    InputPointType pt1 = iMesh->GetPoint(id1);
    InputPointType pt2 = iMesh->GetPoint(id2);

    InputCoordRepType oValue = 1.0 / pt1.EuclideanDistanceTo(pt2);

    return oValue;
  }
};

template< typename TInputMesh >
class ConformalMatrixCoefficients:public MatrixCoefficients< TInputMesh >
{
public:
  using Superclass = MatrixCoefficients< TInputMesh >;

  using InputMeshType = TInputMesh;
  using InputCoordRepType = typename InputMeshType::CoordRepType;
  using InputPointType = typename InputMeshType::PointType;
  using InputPointIdentifier = typename InputMeshType::PointIdentifier;
  using InputQEType = typename InputMeshType::QEType;

  ConformalMatrixCoefficients() = default;

  InputCoordRepType operator()(const InputMeshType *iMesh, InputQEType *iEdge) const override
  {
    InputPointIdentifier id1 = iEdge->GetOrigin();
    InputPointIdentifier id2 = iEdge->GetDestination();
    InputPointType       pt1 = iMesh->GetPoint(id1);
    InputPointType       pt2 = iMesh->GetPoint(id2);

    InputCoordRepType oValue(0.0);

    if ( iEdge->IsLeftSet() )
      {
      InputPointIdentifier idA = iEdge->GetLnext()->GetDestination();
      InputPointType       ptA = iMesh->GetPoint(idA);
      oValue += TriangleHelper< InputPointType >::Cotangent(pt1, ptA, pt2);
      }
    if ( iEdge->IsRightSet() )
      {
      InputPointIdentifier idB = iEdge->GetRnext()->GetOrigin();
      InputPointType       ptB = iMesh->GetPoint(idB);
      oValue += TriangleHelper< InputPointType >::Cotangent(pt1, ptB, pt2);
      }

    return std::max( NumericTraits< InputCoordRepType >::ZeroValue(), oValue);
  }
};

template< typename TInputMesh >
class AuthalicMatrixCoefficients:public MatrixCoefficients< TInputMesh >
{
public:
  using Superclass = MatrixCoefficients< TInputMesh >;

  using InputMeshType = TInputMesh;
  using InputCoordRepType = typename InputMeshType::CoordRepType;
  using InputPointType = typename InputMeshType::PointType;
  using InputPointIdentifier = typename InputMeshType::PointIdentifier;
  using InputQEType = typename InputMeshType::QEType;

  AuthalicMatrixCoefficients() = default;

  InputCoordRepType operator()(const InputMeshType *iMesh, InputQEType *iEdge) const override
  {
    InputPointIdentifier id1 = iEdge->GetOrigin();
    InputPointType       pt1 = iMesh->GetPoint(id1);

    InputPointIdentifier id2 = iEdge->GetDestination();
    InputPointType       pt2 = iMesh->GetPoint(id2);

    InputCoordRepType oValue = NumericTraits< InputCoordRepType >::ZeroValue();

    if ( iEdge->IsLeftSet() )
      {
      InputPointIdentifier idA = iEdge->GetLnext()->GetDestination();
      InputPointType       ptA = iMesh->GetPoint(idA);
      oValue +=
        TriangleHelper< InputPointType >::Cotangent(pt1, pt2, ptA);
      }

    if ( iEdge->IsRightSet() )
      {
      InputPointIdentifier idB = iEdge->GetRnext()->GetOrigin();
      InputPointType       ptB = iMesh->GetPoint(idB);
      oValue += TriangleHelper< InputPointType >::Cotangent(pt1, pt2, ptB);
      }

    return oValue / pt1.SquaredEuclideanDistanceTo(pt2);
  }
};

template< typename TInputMesh >
class IntrinsicMatrixCoefficients:public MatrixCoefficients< TInputMesh >
{
public:
  using Superclass = MatrixCoefficients< TInputMesh >;

  using InputMeshType = TInputMesh;
  using InputCoordRepType = typename InputMeshType::CoordRepType;
  using InputQEType = typename InputMeshType::QEType;

  InputCoordRepType m_Lambda;

  IntrinsicMatrixCoefficients(const InputCoordRepType & iLambda):
    m_Lambda(iLambda)
  {}

  InputCoordRepType operator()(const InputMeshType *iMesh,
                               InputQEType *iEdge) const
  {
    AuthalicMatrixCoefficients< TInputMesh >  authalic;
    ConformalMatrixCoefficients< TInputMesh > conformal;

    InputCoordRepType oValue = m_Lambda * conformal(iMesh, iEdge)
                               + ( 1.0 - m_Lambda ) * authalic(iMesh, iEdge);

    return oValue;
  }
};

template< typename TInputMesh >
class HarmonicMatrixCoefficients:public MatrixCoefficients< TInputMesh >
{
public:
  using Superclass = MatrixCoefficients< TInputMesh >;

  using InputMeshType = TInputMesh;
  using InputCoordRepType = typename InputMeshType::CoordRepType;
  using InputPointType = typename InputMeshType::PointType;
  using InputVectorType = typename InputPointType::VectorType;
  using InputPointIdentifier = typename InputMeshType::PointIdentifier;
  using InputQEType = typename InputMeshType::QEType;

  static constexpr unsigned int PointDimension = InputPointType::PointDimension;

  HarmonicMatrixCoefficients() = default;

  InputCoordRepType operator()(const InputMeshType *iMesh, InputQEType *iEdge) const override
  {
    InputPointIdentifier id1 = iEdge->GetOrigin();
    InputPointIdentifier id2 = iEdge->GetDestination();

    InputPointIdentifier idA = iEdge->GetLnext()->GetDestination();
    InputPointIdentifier idB = iEdge->GetRnext()->GetOrigin();

    InputPointType pt1 = iMesh->GetPoint(id1);
    InputPointType pt2 = iMesh->GetPoint(id2);
    InputPointType ptA = iMesh->GetPoint(idA);
    InputPointType ptB = iMesh->GetPoint(idB);

    InputVectorType v1A = ptA - pt1;
    InputVectorType v1B = ptB - pt1;
    InputVectorType v12 = pt2 - pt1;

    InputCoordRepType L1A = v1A * v1A;
    InputCoordRepType L1B = v1B * v1B;
    InputCoordRepType L12 = v12 * v12;

    InputCoordRepType L2A = pt2.SquaredEuclideanDistanceTo(ptA);
    InputCoordRepType L2B = pt2.SquaredEuclideanDistanceTo(ptB);

    CrossHelper< InputVectorType > cross;

    InputCoordRepType AreaA = 0.5 * ( cross(v1A, v12).GetNorm() );
    InputCoordRepType AreaB = 0.5 * ( cross(v1B, v12).GetNorm() );

    InputCoordRepType
      oValue = ( L1A + L2A - L12 ) / AreaA + ( L1B + L2B - L12 ) / AreaB;

    return oValue;
  }
};
}
#endif