BezierSpline.hh

#ifndef GRAIL_BEZIER_SPLINE_H
#define GRAIL_BEZIER_SPLINE_H
 
#include "Curve.hh"
#include "../../libs/vector2/vector2.h"
 
#include <cassert>
 
namespace grail
{
    namespace curves
    {
        template <typename ContextType>
        class BezierSpline final : public Curve<ContextType>
        {
        public:
            BezierSpline(std::shared_ptr<utility::Evaluator<ContextType>> childEvaluator,
                         const std::vector<Vector2>& points,
                         const std::vector<Vector2>& tangents,
                         float epsilon = 0.001f)
                : Curve<ContextType>{childEvaluator}, points(points), tangents{ tangents }, epsilon{ epsilon }
            {
                Initialize();
            }
 
            BezierSpline(std::shared_ptr<utility::Evaluator<ContextType>> childEvaluator,
                         std::vector<Vector2>&& points,
                         std::vector<Vector2>&& tangents,
                         float epsilon = 0.001f)
                : Curve<ContextType>{childEvaluator}, points(std::move(points)), tangents(std::move(tangents)), epsilon{ epsilon }
            {
                Initialize();
            }
 
            virtual float Sample(float argument) const override final
            {
                if (argument <= points.front().x)
                    return points.front().y;
                else if (argument > points.back().x)
                    return points.back().y;
 
                for (std::size_t i = 1; i < points.size(); ++i)
                {
                    if (argument == points[i].x)
                        return points[i].y;
 
                    if (argument > points[i - 1].x && argument < points[i].x)
                    {
                        float begin = 0;
                        float end = 1;
                        Vector2 result;
 
                        do
                        {
                            float middle = (begin + end) / 2.0f;
                            result = BezierInterpolation(i - 1, i, middle);
 
                            if (result.x < argument)
                                begin = middle;
                            else
                                end = middle;
                        }
                        while (end - begin > epsilon);
 
                        return result.y;
                    }
                }
                return 0;
            }
 
            std::vector<Vector2>& GetPoints() { return points; }
            const std::vector<Vector2>& GetPoints() const { return points; }
 
            std::vector<Vector2>& GetTangents() { return tangents; }
            const std::vector<Vector2>& GetTangents() const { return tangents; }
 
            EvaluatorType GetEvaluatorType() const override final { return EvaluatorType::CURVE_BEZIER; }
 
        private:
            void Initialize()
            {
                assert(points.size() > 1);
                assert(points.size() == tangents.size());
                Vector2 leftVelocity = tangents.front();
                leftGradient = leftVelocity.x != 0 ? leftVelocity.y / leftVelocity.x : std::numeric_limits<float>::quiet_NaN();
                leftIntercept = leftVelocity.y - leftGradient * leftVelocity.x;
 
                Vector2 rightVelocity = tangents.back();
                rightGradient = rightVelocity.x != 0 ? rightVelocity.y / rightVelocity.x : std::numeric_limits<float>::quiet_NaN();
                rightIntercept = rightVelocity.y - rightGradient * rightVelocity.x;
            }
 
            Vector2 BezierInterpolation(const std::size_t beginIndex, const std::size_t endIndex, float t) const
            {
                static const float BINOMIAL_COEFFICIENTS[4] = { 1, 3, 3, 1 };
                Vector2 p0 = points[beginIndex];
                Vector2 p3 = points[endIndex];
                Vector2 p1 = p0 + tangents[beginIndex];
                Vector2 p2 = p3 - tangents[endIndex];
                Vector2 interpPoints[4]{ p0, p1, p2, p3 };
 
                float tComponent = 1; //to start with t^0
                float oneMinusT = 1 - t;
                float oneMinusTComponent = oneMinusT * oneMinusT * oneMinusT; //to start with (1 - t)^3
                Vector2 result{ 0, 0 };
                for (int pointIndex = 0; pointIndex < 4; ++pointIndex)
                {
                    result += BINOMIAL_COEFFICIENTS[pointIndex] *
                        tComponent *
                        oneMinusTComponent *
                        interpPoints[pointIndex]; //p(i) * t^(3-i) * (1-t)^(i)
                    tComponent *= t; //increment exponent
                    oneMinusTComponent /= (1 - t); //decrement exponent
                }
 
                return result;
            }
 
            std::vector<Vector2> points{};
            std::vector<Vector2> tangents{};
            float epsilon = 0.0f;
 
            float leftGradient = 0.0f;
            float leftIntercept = 0.0f;
            float rightGradient = 0.0f;
            float rightIntercept = 0.0f;
        };
    }
}
 
#endif // GRAIL_BEZIER_SPLINE_H