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 PONCA_MULTIARCH_INCLUDE_STD(cmath)
#include PONCA_MULTIARCH_INCLUDE_STD(limits)
 
#include <Eigen/Core>
 
namespace Ponca
{
 
template < class DataPoint, class _WFunctor, typename T >
class AlgebraicSphere : public T
{
    PONCA_FITTING_DECLARE_DEFAULT_TYPES
 
protected:
    enum
    {
        check = Base::PROVIDES_PRIMITIVE_BASE,  
        PROVIDES_ALGEBRAIC_SPHERE               
    };
 
protected:
    bool m_isNormalized;
 
// results
public:
    Scalar m_uc {0},       
           m_uq {0};       
    VectorType m_ul {VectorType::Zero()};   
public:
    PONCA_EXPLICIT_CAST_OPERATORS(AlgebraicSphere,algebraicSphere)
 
    
    PONCA_MULTIARCH inline void init()
    {
        Base::init();
 
        m_uc = Scalar(0);
        m_ul = VectorType::Zero();
        m_uq = Scalar(0);
 
        m_isNormalized = false;
    }
 
    PONCA_MULTIARCH 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, WFunctor, 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 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, WFunctor, T>& other) const{
        return ! ((*this) == other);
    }
 
    PONCA_MULTIARCH inline void changeBasis(const VectorType& newbasis)
    {
        VectorType diff = Base::m_w.basisCenter() - newbasis;
        Base::setWeightFunc({newbasis, Base::m_w.evalScale()});
        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 inline Scalar prattNorm() const
    {
        PONCA_MULTIARCH_STD_MATH(sqrt);
        return sqrt(prattNorm2());
    }
 
    PONCA_MULTIARCH 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 inline Scalar radius() const
    {
        if(isPlane())
        {
            //return infinity (non-sense value)
#ifdef __CUDACC__
            Scalar inf = 0.;
            return Scalar(1.)/inf;
#else
            return std::numeric_limits<Scalar>::infinity();
#endif
        }
 
        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 inline VectorType center() const
    {
        if(isPlane())
        {
            //return infinity (non-sense value)
#ifdef __CUDACC__
            return VectorType::Constant(Scalar(1.)/Scalar(0));
#else
            return VectorType::Constant(std::numeric_limits<Scalar>::infinity());
#endif
        }
 
        Scalar b = Scalar(1.)/m_uq;
        return Base::m_w.convertToGlobalBasis((Scalar(-0.5)*b)*m_ul);
    }
 
    PONCA_MULTIARCH inline bool isNormalized() const { return m_isNormalized; }
 
    PONCA_MULTIARCH inline Scalar potential (const VectorType& _q) const;
 
    PONCA_MULTIARCH inline Scalar potential() const { return m_uc; }
 
    PONCA_MULTIARCH inline VectorType project (const VectorType& _q) const;
 
    PONCA_MULTIARCH inline VectorType projectDescent (const VectorType& _q, int nbIter = 16) const;
 
    PONCA_MULTIARCH inline VectorType primitiveGradient (const VectorType& _q) const;
 
    PONCA_MULTIARCH inline const VectorType& primitiveGradient () const { return m_ul; }
 
    PONCA_MULTIARCH 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;
    }
 
}; //class AlgebraicSphere
 
#include "algebraicSphere.hpp"
}