BezierSpline.hh

// Copyright QED Software 2023.
 
#ifndef GRAIL_BEZIER_SPLINE_H
#define GRAIL_BEZIER_SPLINE_H
 
#include "Curve.hh"
#include "../../libs/vector2/vector2.h"
 
#include <cassert>
 
namespace grail
{
namespace evaluator
{
 template <typename ContextType>
 class BezierSpline final : public Curve<ContextType>
 {
 public:
 BezierSpline(std::shared_ptr<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<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; }
 
 virtual data::EvaluatorType GetEvaluatorType() const override final { return data::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