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"
 
#include PONCA_MULTIARCH_INCLUDE_STD(iterator)
 
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...>
    {
    };
}
#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>
        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::getNeighborFilter().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>
        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>
        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