algebraicSphere.h

/*
 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*/
 
#pragma once
 
#include "../../defines.h"
#include "../concepts.h"
#include "../../Filters/concepts.h" // HasLocalFrame
 
#include PONCA_MULTIARCH_INCLUDE_STD(cmath)
#include PONCA_MULTIARCH_INCLUDE_CU_STD(limits)
 
#include <Eigen/Core>
 
#define ALGEBRAIC_SPHERE_REQUIREMENTS ProvidesBasketUnitBase<T>&& HasLocalFrame<_NFilter>
 
namespace Ponca
{

        by the evaluation position, the public methods involving a query (such as
        project, potential, gradient) have to be defined in global
        coordinates (e.g. you don't need to convert your query in the current locale
        frame).
 
        This primitive satisfies: ProvidesAlgebraicSphere
 
        \todo Deal with planar case (_uq == 0) and what about _ul == 0 ?
    */
 
    template <class DataPoint, class _NFilter, typename T>
        requires ALGEBRAIC_SPHERE_REQUIREMENTS
    class AlgebraicSphere : public T
    {
        PONCA_FITTING_DECLARE_DEFAULT_TYPES
 
    protected:
        bool m_isNormalized{false};
 
        // results
    public:
        Scalar m_uc{0},                      
            m_uq{0};                         
        VectorType m_ul{VectorType::Zero()}; 
 
    public:
        PONCA_EXPLICIT_BASKET_CAST_OPERATORS(AlgebraicSphere, algebraicSphere)
        PONCA_EXPLICIT_BASKET_CAST_OPERATORS(AlgebraicSphere, implicitPrimitive)
        PONCA_EXPLICIT_BASKET_CAST_OPERATORS(AlgebraicSphere, projectionOperator)
        PONCA_EXPLICIT_BASKET_CAST_OPERATORS(AlgebraicSphere, meanCurvature)

        PONCA_MULTIARCH inline void init()
        {
            Base::init();
 
            m_uc = Scalar(0);
            m_ul = VectorType::Zero();
            m_uq = Scalar(0);
 
            m_isNormalized = false;
        }

        PONCA_MULTIARCH [[nodiscard]] inline bool isValid() const
        {
            return !(m_ul.isApprox(VectorType::Zero()) && m_uc == Scalar(0) && m_uq == Scalar(0));
        }

        PONCA_MULTIARCH inline bool operator==(const AlgebraicSphere<DataPoint, NeighborFilter, T>& other) const
        {
            PONCA_MULTIARCH_STD_MATH(pow);
            const Scalar epsilon        = Eigen::NumTraits<Scalar>::dummy_precision();
            const Scalar squaredEpsilon = epsilon * epsilon;
            return pow(m_uc - other.m_uc, Scalar(2)) < squaredEpsilon &&
                   pow(m_uq - other.m_uq, Scalar(2)) < squaredEpsilon && m_ul.isApprox(other.m_ul);
        }

        template <typename Other>
        PONCA_MULTIARCH [[nodiscard]] inline bool isApprox(
            const Other& other, const Scalar& epsilon = Eigen::NumTraits<Scalar>::dummy_precision()) const
        {
            PONCA_MULTIARCH_STD_MATH(pow);
            const Scalar squaredEpsilon = epsilon * epsilon;
            return pow(m_uc - other.m_uc, Scalar(2)) < squaredEpsilon &&
                   pow(m_uq - other.m_uq, Scalar(2)) < squaredEpsilon && m_ul.isApprox(other.m_ul);
        }

        PONCA_MULTIARCH inline bool operator!=(const AlgebraicSphere<DataPoint, NeighborFilter, T>& other) const
        {
            return !((*this) == other);
        }

        PONCA_MULTIARCH inline void changeBasis(const VectorType& newbasis)
        {
            auto& frame = Base::getNeighborFrame();
 
            VectorType diff = frame.center() - newbasis;
            frame.changeNeighborhoodFrame(newbasis);
            Base::init();
            m_uc = m_uc - m_ul.dot(diff) + m_uq * diff.dot(diff);
            m_ul = m_ul - Scalar(2.) * m_uq * diff;
            // m_uq is not changed
            m_isNormalized = false;
            applyPrattNorm();
        }

        PONCA_MULTIARCH [[nodiscard]] inline Scalar prattNorm() const
        {
            PONCA_MULTIARCH_STD_MATH(sqrt);
            return sqrt(prattNorm2());
        }

        PONCA_MULTIARCH [[nodiscard]] inline Scalar prattNorm2() const
        {
            return m_ul.squaredNorm() - Scalar(4.) * m_uc * m_uq;
        }

        PONCA_MULTIARCH inline bool applyPrattNorm()
        {
            if (!m_isNormalized)
            {
                Scalar pn = prattNorm();
                m_uc /= pn;
                m_ul *= Scalar(1.) / pn;
                m_uq /= pn;
 
                m_isNormalized = true;
            }
            return true;
        }

        PONCA_MULTIARCH [[nodiscard]] inline Scalar radius() const
        {
            PONCA_MULTIARCH_CU_STD_FUNC(numeric_limits);
            if (isPlane())
                return numeric_limits<Scalar>::infinity(); // non-sense value
 
            PONCA_MULTIARCH_STD_MATH(sqrt);
            Scalar b = Scalar(1.) / m_uq;
            return Scalar(sqrt(((Scalar(-0.5) * b) * m_ul).squaredNorm() - m_uc * b));
        }
 
        PONCA_MULTIARCH [[nodiscard]] inline Scalar kMean() const { return isPlane() ? Scalar(0) : Scalar(2) * m_uq; }

        PONCA_MULTIARCH [[nodiscard]] inline VectorType center() const
        {
            PONCA_MULTIARCH_CU_STD_FUNC(numeric_limits);
            if (isPlane())
                return VectorType::Constant(numeric_limits<Scalar>::infinity()); // non-sense value
 
            Scalar b = Scalar(1.) / m_uq;
            return Base::getNeighborFrame().convertToGlobalBasis((Scalar(-0.5) * b) * m_ul);
        }

        PONCA_MULTIARCH [[nodiscard]] inline bool isNormalized() const { return m_isNormalized; }

        PONCA_MULTIARCH [[nodiscard]] inline Scalar potential(const VectorType& _q) const
        {
            // Turn to centered basis
            const VectorType lq = Base::getNeighborFrame().convertToLocalBasis(_q);
            return potentialLocal(lq);
        }

        PONCA_MULTIARCH [[nodiscard]] inline Scalar potential() const { return m_uc; }

        PONCA_MULTIARCH [[nodiscard]] inline VectorType project(const VectorType& _q) const;

        PONCA_MULTIARCH [[nodiscard]] inline VectorType primitiveGradient(const VectorType& _q) const
        {
            // Turn to centered basis
            const VectorType lq = Base::getNeighborFrame().convertToLocalBasis(_q);
            return primitiveGradientLocal(lq);
        }

        PONCA_MULTIARCH [[nodiscard]] inline const VectorType& primitiveGradient() const { return m_ul; }

        PONCA_MULTIARCH [[nodiscard]] inline bool isPlane() const
        {
            PONCA_MULTIARCH_STD_MATH(abs);
            bool bPlanar = Eigen::internal::isApprox(m_uq, Scalar(0));
            bool bReady  = Base::isReady();
            return bReady && bPlanar;
        }
 
    protected:
        PONCA_MULTIARCH [[nodiscard]] inline Scalar potentialLocal(const VectorType& _lq) const;

        PONCA_MULTIARCH [[nodiscard]] inline VectorType primitiveGradientLocal(const VectorType& _lq) const;
    }; // class AlgebraicSphere
 
#include "algebraicSphere.hpp"
} // namespace Ponca