basket.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 "enums.h"
#include "primitive.h"
#include "compute.h"
 
namespace Ponca
{
 
#define BSKNF typename BasketType::NeighborFilter
#define BSKP typename BasketType::DataPoint
 
#ifndef PARSED_WITH_DOXYGEN
    namespace internal
    {
        template <class P, class NF, typename Aggregate, template <class, class, typename> class Ext,
                  template <class, class, typename> class... Exts>
        struct BasketAggregateImpl
        {
            using type = typename BasketAggregateImpl<P, NF, Ext<P, NF, Aggregate>, Exts...>::type;
        };
 
        template <class P, class NF, typename Aggregate, template <class, class, typename> class Ext>
        struct BasketAggregateImpl<P, NF, Aggregate, Ext>
        {
            using type = Ext<P, NF, Aggregate>;
        };
 
        template <class P, class NF, template <class, class, typename> class... Exts>
        struct BasketAggregate : BasketAggregateImpl<P, NF, PrimitiveBase<P, NF>, Exts...>
        {
        };
 
        template <int Type, typename BasketType, template <class, class, int, typename> class Ext,
                  template <class, class, int, typename> class... Exts>
        struct BasketDiffAggregateImpl
        {
            using type = typename BasketDiffAggregateImpl<Type, Ext<BSKP, BSKNF, Type, BasketType>, Exts...>::type;
        };
 
        template <int Type, typename BasketType, template <class, class, int, typename> class Ext>
        struct BasketDiffAggregateImpl<Type, BasketType, Ext>
        {
            using type = Ext<BSKP, BSKNF, Type, BasketType>;
        };
 
        template <typename BasketType, int Type, template <class, class, int, typename> class... Exts>
        struct BasketDiffAggregate : BasketDiffAggregateImpl<Type, BasketType, PrimitiveDer, Exts...>
        {
        };
    } // namespace internal
#endif
 
    template <typename _Derived, typename _Base>
    struct BasketComputeObject : public ComputeObject<_Derived>, public virtual _Base
    {
        using Base    = _Base;    
        using Derived = _Derived; 
        using Scalar  = typename Base::Scalar;
 
    protected:
        using ComputeObject<Derived>::derived;
 
    public:
        using ComputeObject<Derived>::compute; // Makes the default compute(container) accessible when using a CPU
                                               // architecture
 
        template <typename IteratorBegin, typename IteratorEnd>
        PONCA_MULTIARCH inline FIT_RESULT compute(const IteratorBegin& begin, const IteratorEnd& end)
        {
            Base::init();
            FIT_RESULT res = UNDEFINED;
 
            do
            {
                derived().startNewPass();
                for (auto it = begin; it != end; ++it)
                {
                    derived().addNeighbor(*it);
                }
                res = Base::finalize();
            } while (res == NEED_OTHER_PASS);
 
            return res;
        }
 
        template <typename IndexRange, typename PointContainer>
        PONCA_MULTIARCH inline FIT_RESULT computeWithIds(IndexRange ids, const PointContainer& points)
        {
            Base::init();
            FIT_RESULT res = UNDEFINED;
 
            do
            {
                derived().startNewPass();
                for (const auto& i : ids)
                {
                    derived().addNeighbor(points[i]);
                }
                res = Base::finalize();
            } while (res == NEED_OTHER_PASS);
 
            return res;
        }
 
    protected:
        template <typename Func>
        PONCA_MULTIARCH FIT_RESULT computeMLSImpl(Func&& computeFunc, const int mlsIter, const Scalar epsilon)
        {
            FIT_RESULT res = UNDEFINED;
            auto lastPos   = Base::getNeighborFilter().evalPos();
 
            for (int mm = 0; mm < mlsIter; ++mm)
            {
                Base::m_nFilter.changeNeighborhoodFrame(lastPos);
                res = computeFunc();
 
                if (Base::isStable())
                {
                    auto newPos = Base::project(lastPos);
                    if (newPos.isApprox(lastPos, epsilon))
                        return res;
                    lastPos = newPos;
                }
                else
                {
                    return res;
                }
            }
            return res;
        }
 
    public:
        template <typename PointContainer>
        PONCA_MULTIARCH FIT_RESULT computeMLS(const PointContainer& points, const int mlsIter = 5,
                                              const Scalar epsilon = Eigen::NumTraits<Scalar>::dummy_precision())
        {
            return computeMLSImpl([&]() { return compute(points); }, mlsIter, epsilon);
        }
 
        template <typename IndexRange, typename PointContainer>
        PONCA_MULTIARCH FIT_RESULT computeWithIdsMLS(const IndexRange& ids, const PointContainer& points,
                                                     const int mlsIter    = 5,
                                                     const Scalar epsilon = Eigen::NumTraits<Scalar>::dummy_precision())
        {
            return computeMLSImpl([&]() { return computeWithIds(ids, points); }, mlsIter, epsilon);
        }
    };
 
#define WRITE_COMPUTE_FUNCTIONS                            \
    using BasketComputeObject<Self, Base>::compute;        \
    using BasketComputeObject<Self, Base>::computeWithIds; \
    using BasketComputeObject<Self, Base>::computeMLS;     \
    using BasketComputeObject<Self, Base>::computeWithIdsMLS;
 
    template <typename BasketType, int Type, template <class, class, int, typename> class Ext0,
              template <class, class, int, typename> class... Exts>
    class BasketDiff
        : public BasketComputeObject<BasketDiff<BasketType, Type, Ext0, Exts...>,
                                     typename internal::BasketDiffAggregate<BasketType, Type, Ext0, Exts...>::type>
    {
    private:
        using Self = BasketDiff;
 
    public:
        using Base = typename internal::BasketDiffAggregate<BasketType, Type, Ext0, Exts...>::type;
        using NeighborFilter = BSKNF;
        using DataPoint = BSKP;
        using Scalar = typename BasketType::Scalar;
        WRITE_COMPUTE_FUNCTIONS
 
        PONCA_MULTIARCH inline bool addNeighbor(const DataPoint& _nei)
        {
            // compute weight
            auto neiFilterOutput = Base::getNeighborFilter()(_nei);
            typename Base::ScalarArray dw;
 
            if (neiFilterOutput.first > Scalar(0.))
            {
                Base::addLocalNeighbor(neiFilterOutput.first, neiFilterOutput.second, _nei, dw);
                return true;
            }
            return false;
        }
    };
 
    template <class P, class NF, template <class, class, typename> class Ext0,
              template <class, class, typename> class... Exts>
    class Basket : public BasketComputeObject<Basket<P, NF, Ext0, Exts...>,
                                              typename internal::BasketAggregate<P, NF, Ext0, Exts...>::type>
    {
    private:
        using Self = Basket;
 
    public:
        using Base = typename internal::BasketAggregate<P, NF, Ext0, Exts...>::type;
        using Scalar     = typename P::Scalar;
        using VectorType = typename P::VectorType;
        using DataPoint = P;
        using NeighborFilter = NF;
 
        WRITE_COMPUTE_FUNCTIONS
 
        PONCA_MULTIARCH inline bool addNeighbor(const DataPoint& _nei)
        {
            // compute weight
            auto neiFilterOutput = Base::getNeighborFilter()(_nei);
 
            if (neiFilterOutput.first > Scalar(0.))
            {
                Base::addLocalNeighbor(neiFilterOutput.first, neiFilterOutput.second, _nei);
                return true;
            }
            return false;
        }
 
        PONCA_MULTIARCH [[nodiscard]] inline VectorType projectDescent(const VectorType& _q, int nbIter = 16) const
        {
            PONCA_MULTIARCH_STD_MATH(min)
 
            // turn to centered basis
            const VectorType lq = Base::getNeighborFilter().convertToLocalBasis(_q);
 
            VectorType grad;
            VectorType dir  = Base::primitiveGradientLocal(lq);
            Scalar ilg      = Scalar(1.) / dir.norm();
            dir             = dir * ilg;
            Scalar ad       = Base::potentialLocal(lq);
            Scalar delta    = -ad * min(ilg, Scalar(1.));
            VectorType proj = lq + dir * delta;
 
            for (int i = 0; i < nbIter; ++i)
            {
                grad  = Base::primitiveGradientLocal(proj);
                ilg   = Scalar(1.) / grad.norm();
                delta = -Base::potentialLocal(proj) * min(ilg, Scalar(1.));
                proj += dir * delta;
            }
            return Base::getNeighborFilter().convertToGlobalBasis(proj);
        }
    }; // class Basket
 
#undef WRITE_COMPUTE_FUNCTIONS
} // namespace Ponca