project.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 "../concepts.h"
 
#include <algorithm> // std::min
 
namespace Ponca
{
    struct DirectProjectionOperator
    {
    public:
        template <ProvidesProjectionOperator ComputeObject>
        PONCA_MULTIARCH typename ComputeObject::VectorType project(const ComputeObject& _f,
                                                                   const typename ComputeObject::VectorType& _pos) const
        {
            return _f.project(_pos);
        }
 
        template <typename ComputeObject>
        PONCA_MULTIARCH typename ComputeObject::VectorType operator()(
            const ComputeObject& _f, const typename ComputeObject::VectorType& _pos) const
        {
            return project(_f, _pos);
        }
    };
 
    struct GradientDescentProjectionOperator
    {
    public:
        GradientDescentProjectionOperator(unsigned int _nbIter = 16) : nbIter(_nbIter) {}
 
        template <ProvidesImplicitPrimitive ComputeObject>
        PONCA_MULTIARCH typename ComputeObject::VectorType project(const ComputeObject& _f,
                                                                   const typename ComputeObject::VectorType& _pos) const
        {
            PONCA_MULTIARCH_STD_MATH(min)
            using VectorType = typename ComputeObject::VectorType;
            using Scalar     = typename ComputeObject::Scalar;
 
            VectorType grad;
            VectorType dir  = _f.primitiveGradient(_pos);
            Scalar ilg      = Scalar(1.) / dir.norm();
            dir             = dir * ilg;
            Scalar ad       = _f.potential(_pos);
            Scalar delta    = -ad * min(ilg, Scalar(1.));
            VectorType proj = _pos + dir * delta;
 
            for (unsigned int i = 0; i < nbIter; ++i)
            {
                grad  = _f.primitiveGradient(proj);
                ilg   = Scalar(1.) / grad.norm();
                delta = -_f.potential(proj) * min(ilg, Scalar(1.));
                proj += dir * delta;
            }
            return proj;
        }
 
        template <typename ComputeObject>
        PONCA_MULTIARCH typename ComputeObject::VectorType operator()(
            const ComputeObject& f, const typename ComputeObject::VectorType& pos) const
        {
            return project(f, pos);
        }
 
    public:
        unsigned int nbIter;
    };
}; // namespace Ponca