weingarten.hpp

#include <Eigen/Eigenvalues>
#include <Eigen/Core>
 
namespace Ponca
{
    template <class DataPoint, class _NFilter, typename T>
    typename FundamentalFormWeingartenEstimator<DataPoint, _NFilter, T>::Matrix2 FundamentalFormWeingartenEstimator<
        DataPoint, _NFilter, T>::firstFundamentalForm() const
    {
        Matrix2 first;
        firstFundamentalForm(first);
        return first;
    }
 
    template <class DataPoint, class _NFilter, typename T>
    template <typename Matrix2Derived>
    void FundamentalFormWeingartenEstimator<DataPoint, _NFilter, T>::firstFundamentalForm(Matrix2Derived& first) const
    {
        Base::firstFundamentalFormComponents(first(0, 0), first(1, 0), first(1, 1));
        first(0, 1) = first(1, 0); // diagonal
    }
 
    template <class DataPoint, class _NFilter, typename T>
    typename FundamentalFormWeingartenEstimator<DataPoint, _NFilter, T>::Matrix2 FundamentalFormWeingartenEstimator<
        DataPoint, _NFilter, T>::secondFundamentalForm() const
    {
        Matrix2 second;
        secondFundamentalForm(second);
        return second;
    }
 
    template <class DataPoint, class _NFilter, typename T>
    template <typename Matrix2Derived>
    void FundamentalFormWeingartenEstimator<DataPoint, _NFilter, T>::secondFundamentalForm(Matrix2Derived& second) const
    {
        Base::secondFundamentalFormComponents(second(0, 0), second(1, 0), second(1, 1));
        second(0, 1) = second(1, 0); // diagonal
    }
 
    template <class DataPoint, class _NFilter, typename T>
    typename FundamentalFormWeingartenEstimator<DataPoint, _NFilter, T>::Matrix2 FundamentalFormWeingartenEstimator<
        DataPoint, _NFilter, T>::weingartenMap() const
    {
        Matrix2 w;
        weingartenMap(w);
        return w;
    }
 
    template <class DataPoint, class _NFilter, typename T>
    template <typename Matrix2Derived>
    void FundamentalFormWeingartenEstimator<DataPoint, _NFilter, T>::weingartenMap(Matrix2Derived& w) const
    {
        w = firstFundamentalForm().inverse() * secondFundamentalForm();
    }
 
    template <class DataPoint, class _NFilter, typename T>
    typename FundamentalFormWeingartenEstimator<DataPoint, _NFilter, T>::Scalar FundamentalFormWeingartenEstimator<
        DataPoint, _NFilter, T>::kMean() const
    {
        Scalar E, F, G, L, M, N;
        Base::firstFundamentalFormComponents(E, F, G);
        Base::secondFundamentalFormComponents(L, M, N);
        return (G * L - Scalar(2) * F * M + E * N) / (Scalar(2) * (E * G - F * F));
    }
 
    template <class DataPoint, class _NFilter, typename T>
    typename FundamentalFormWeingartenEstimator<DataPoint, _NFilter, T>::Scalar FundamentalFormWeingartenEstimator<
        DataPoint, _NFilter, T>::GaussianCurvature() const
    {
        Scalar E, F, G, L, M, N;
        Base::firstFundamentalFormComponents(E, F, G);
        Base::secondFundamentalFormComponents(L, M, N);
        return (L * N - M * M) / (E * G - F * F);
    }
 
    template <class DataPoint, class _NFilter, int DiffType, typename T>
    typename NormalDerivativeWeingartenEstimator<DataPoint, _NFilter, DiffType, T>::Matrix2
    NormalDerivativeWeingartenEstimator<DataPoint, _NFilter, DiffType, T>::weingartenMap() const
    {
        Matrix2 w;
        weingartenMap(w);
        return w;
    }
 
    template <class DataPoint, class _NFilter, int DiffType, typename T>
    FIT_RESULT NormalDerivativeWeingartenEstimator<DataPoint, _NFilter, DiffType, T>::finalize()
    {
 
        PONCA_MULTIARCH_STD_MATH(abs);
        PONCA_MULTIARCH_STD_MATH(sqrt);
 
        using Index = typename VectorType::Index;
 
        Base::finalize();
        // Test if base finalize end on a viable case (stable / unstable)
        if (this->isReady())
        {
            Index i0 = Index(-1), i1 = Index(-1), i2 = Index(-1);
 
            MatrixType dN = Base::dNormal().template middleCols<DataPoint::Dim>(Base::isScaleDer() ? 1 : 0);
 
            VectorType n = Base::primitiveGradient();
            n.array().abs().minCoeff(&i0); // i0: dimension where n extends the least
            i1 = (i0 + 1) % 3;
            i2 = (i0 + 2) % 3;
 
            m_tangentBasis.col(0) = n;
 
            m_tangentBasis.col(1)[i0] = 0;
            m_tangentBasis.col(1)[i1] = n[i2];
            m_tangentBasis.col(1)[i2] = -n[i1];
 
            m_tangentBasis.col(1).normalize();
            m_tangentBasis.col(2) = m_tangentBasis.col(1).cross(n);
        }
        return Base::m_eCurrentState;
    }
 
    template <class DataPoint, class _NFilter, int DiffType, typename T>
    template <typename Matrix2Derived>
    void NormalDerivativeWeingartenEstimator<DataPoint, _NFilter, DiffType, T>::weingartenMap(Matrix2Derived& W) const
    {
        PONCA_MULTIARCH_STD_MATH(abs);
        PONCA_MULTIARCH_STD_MATH(sqrt);
 
        using Index    = typename VectorType::Index;
        using Matrix32 = Eigen::Matrix<Scalar, 3, 2>;
 
        // Get the object space Weingarten map dN
        MatrixType dN = Base::dNormal().template middleCols<DataPoint::Dim>(Base::isScaleDer() ? 1 : 0);
 
        // Compute tangent-space change of basis function
        auto B = m_tangentBasis.template rightCols<2>();
 
        // Compute the 2x2 matrix representing the shape operator by transforming dN to the basis B.
        // Recall that dN is a bilinear form, it thus transforms as follows:
        W = B.transpose() * dN * B;
 
        // Recall that at this stage, the shape operator represented by W describes the normal curvature K_n(u) in the
        // direction u \in R^2 as follows:
        //   K_n(u) = u^T W u
        // The principal curvatures are fully defined by the values and the directions of the extrema of K_n.
        //
        // If the normal field N(x) comes from the gradient of a scalar field, then N(x) is curl-free, and dN and W are
        // symmetric matrices. In this case, the extrema of the previous quadratic form are directly obtained by the
        // eigenvalue decomposition of W. However, if N(x) is only an approximation of the normal field of a surface,
        // then N(x) is not necessarily curl-free, and in this case W is not symmetric. In this case, we first have to
        // find an equivalent symmetric matrix W' such that:
        //   K_n(u) = u^T W' u,
        // for any u \in R^2.
        // It is easy to see that such a W' is simply obtained as:
        //   W' = (W + W^T)/2
        W(0, 1) = W(1, 0) = (W(0, 1) + W(1, 0)) / Scalar(2);
    }
 
    template <class DataPoint, class _NFilter, int DiffType, typename T>
    typename NormalDerivativeWeingartenEstimator<DataPoint, _NFilter, DiffType, T>::VectorType
    NormalDerivativeWeingartenEstimator<DataPoint, _NFilter, DiffType, T>::worldToTangentPlane(
        const VectorType& _q, bool _isPositionVector) const
    {
        return m_tangentBasis.normalized().transpose() *
               Base::getNeighborFilter().convertToLocalBasis(_q, _isPositionVector);
    }
 
    template <class DataPoint, class _NFilter, int DiffType, typename T>
    typename NormalDerivativeWeingartenEstimator<DataPoint, _NFilter, DiffType, T>::VectorType
    NormalDerivativeWeingartenEstimator<DataPoint, _NFilter, DiffType, T>::tangentPlaneToWorld(
        const VectorType& _lq, bool _isPositionVector) const
    {
        return Base::getNeighborFilter().convertToGlobalBasis(m_tangentBasis.normalized().transpose().inverse() * _lq,
                                                              _isPositionVector);
    }
 
    namespace internal
    {
        template <class DataPoint, class _NFilter, typename T>
        FIT_RESULT WeingartenCurvatureEstimatorBase<DataPoint, _NFilter, T>::finalize()
        {
 
            if (Base::finalize() != STABLE)
                return Base::m_eCurrentState;
 
            Matrix2 w;
            Base::weingartenMap(w);
 
            // w is self adjoint by construction
            Eigen::SelfAdjointEigenSolver<Matrix2> solver;
            solver.computeDirect(w);
 
            Scalar kmin = solver.eigenvalues().x();
            Scalar kmax = solver.eigenvalues().y();
            VectorType vmin, vmax; // maxi directions
 
            vmin(0)                       = Scalar(0); // set height
            vmax(0)                       = Scalar(0); // set height
            vmin.template bottomRows<2>() = solver.eigenvectors().col(0);
            vmax.template bottomRows<2>() = solver.eigenvectors().col(1);
 
            vmin = Base::tangentPlaneToWorld(vmin, false);
            vmax = Base::tangentPlaneToWorld(vmax, false);
 
            Base::setCurvatureValues(kmin, kmax, vmin, vmax);
 
            return Base::m_eCurrentState;
        }
    } // namespace internal
} // namespace Ponca