curve.h

/********************************************************************
 * Copyright (C) 2015 Liangliang Nan <liangliang.nan@gmail.com>
 * https://3d.bk.tudelft.nl/liangliang/
 *
 * This file is part of Easy3D. If it is useful in your research/work,
 * I would be grateful if you show your appreciation by citing it:
 * ------------------------------------------------------------------
 *      Liangliang Nan.
 *      Easy3D: a lightweight, easy-to-use, and efficient C++ library
 *      for processing and rendering 3D data.
 *      Journal of Open Source Software, 6(64), 3255, 2021.
 * ------------------------------------------------------------------
 *
 * Easy3D is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License Version 3
 * as published by the Free Software Foundation.
 *
 * Easy3D is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program. If not, see <http://www.gnu.org/licenses/>.
 ********************************************************************/
 
#ifndef EASY3D_CORE_CURVE_H
#define EASY3D_CORE_CURVE_H
 
 
#include <vector>
 
 
namespace easy3d {

    namespace curve {

        template <template <size_t, class> class Point, size_t N, typename T>
        void quadratic(const Point<N, T>& A, const Point<N, T>& B, const Point<N, T>& C,
                       std::vector<Point<N, T>>& curve, unsigned int bezier_steps = 4,
                       bool include_end = false)
        {
            using Point_t = Point<N, T>;
            for (unsigned int i = 0; i <= bezier_steps; i++) {   // i <= bezier_steps: to fully sample the curve
                const T t = static_cast<T>(i) / bezier_steps;
                const Point_t U = (1.0 - t) * A + t * B;
                const Point_t V = (1.0 - t) * B + t * C;
                curve.push_back((1.0 - t) * U + t * V);
            }
 
            if (include_end)
                curve.push_back(C);
        }

        template <template <size_t, class> class Point, size_t N, typename T>
        void cubic(const Point<N, T> &A, const Point<N, T> &B, const Point<N, T> &C, const Point<N, T> &D,
                   std::vector<Point<N, T>> &curve,
                   unsigned int bezier_steps = 4,
                   bool include_end = false)
        {
            using Point_t = Point<N, T>;
            for (unsigned int i = 0; i <= bezier_steps; i++) {  // i <= bezier_steps: to fully sample the curve
                const T t = static_cast<T>(i) / bezier_steps;
                const Point_t U = (1.0 - t) * A + t * B;
                const Point_t V = (1.0 - t) * B + t * C;
                const Point_t W = (1.0 - t) * C + t * D;
                const Point_t P = (1.0 - t) * U + t * V;
                const Point_t Q = (1.0 - t) * V + t * W;
                curve.push_back((1.0 - t) * P + t * Q);
            }
            if (include_end)
                curve.push_back(D);
        }
 
    }
 
    // =============================================================================
 

    template<template <size_t, class> class Point, size_t N, typename T>
    class Curve {
    public:
        using Point_t = Point<N, T>;

        virtual std::vector<Point_t> generate(const std::vector<Point_t> &control_points, size_t num_samples) const = 0;
    };

    template<template <size_t, class> class Point, size_t N, typename T>
    class Bezier : public Curve<Point, N, T> {
    public:
        using Point_t = Point<N, T>;

        Bezier() = default;

        std::vector<Point_t> generate(const std::vector<Point_t> &control_points, size_t num_samples) const override {
            std::vector<Point_t> curve;
            curve.reserve(num_samples + 1); // Reserve space to avoid reallocations
 
            for (size_t i = 0; i <= num_samples; ++i) {
                T t = static_cast<T>(i) / num_samples;
                curve.push_back(interpolate(control_points, t));
            }
            return curve;
        }
 
    private:
        Point_t interpolate(const std::vector<Point_t> &control_points, T t) const {
            std::vector<Point_t> points = control_points; // Copy control points
            size_t n = points.size();
            for (size_t k = 1; k < n; ++k) {
                for (size_t i = 0; i < n - k; ++i) {
                    points[i] = points[i] * (1 - t) + points[i + 1] * t;
                }
            }
            return points[0];
        }
    };

    template<template <size_t, class> class Point, size_t N, typename T>
    class BSpline : public Curve<Point, N, T> {
    public:
        using Point_t = Point<N, T>;

        BSpline() = default;

        std::vector<Point_t> generate(const std::vector<Point_t> &control_points, size_t num_samples) const override {
            std::vector<Point_t> curve;
            size_t n = control_points.size();
            if (n < 4)
                return curve;
 
            curve.reserve(num_samples + 1); // Reserve space to avoid reallocations
 
            // Sample the curve
            for (size_t i = 0; i <= num_samples; ++i) {
                T t = static_cast<T>(i) / num_samples * (n - 3);
                size_t k = static_cast<size_t>(t);
                t -= k;
 
                // Ensure k is within the valid range
                if (k >= n - 3) {
                    k = n - 4;
                    t = 1.0;
                }
 
                // Evaluate using a proper B-spline basis
                curve.push_back(interpolate(control_points, t, k));
            }
            return curve;
        }
 
    private:
        Point_t interpolate(const std::vector<Point_t> &control_points, T t, size_t k) const {
            const T one_sixth = static_cast<T>(1) / 6;
            Point_t P0 = control_points[k];
            Point_t P1 = control_points[k + 1];
            Point_t P2 = control_points[k + 2];
            Point_t P3 = control_points[k + 3];
 
            T t2 = t * t;
            T t3 = t2 * t;
 
            return (P0 * (-t3 + 3 * t2 - 3 * t + 1) * one_sixth) +
                   (P1 * (3 * t3 - 6 * t2 + 4) * one_sixth) +
                   (P2 * (-3 * t3 + 3 * t2 + 3 * t + 1) * one_sixth) +
                   (P3 * (t3) * one_sixth);
        }
    };

    template<template <size_t, class> class Point, size_t N, typename T>
    class CatmullRom : public Curve<Point, N, T> {
    public:
        using Point_t = Point<N, T>;

        CatmullRom() = default;

        std::vector<Point_t> generate(const std::vector<Point_t> &control_points, size_t num_samples) const override {
            std::vector<Point_t> curve;
            curve.reserve(num_samples + 1); // Reserve space to avoid reallocations
 
            for (size_t i = 0; i <= num_samples; ++i) {
                T t = static_cast<T>(i) / num_samples;
                curve.push_back(interpolate(control_points, t));
            }
            return curve;
        }
 
    private:
        Point_t interpolate(const std::vector<Point_t> &control_points, T t) const {
            size_t n = control_points.size() - 3;
            size_t i = std::min(static_cast<size_t>(t * n), n - 1);
            T u = t * n - i;
            Point_t P0 = control_points[i], P1 = control_points[i + 1], P2 = control_points[i + 2], P3 = control_points[
                i + 3];
            return ((P0 * (-1) + P1 * 3 - P2 * 3 + P3) * (u * u * u) / 2 +
                    (P0 * 2 - P1 * 5 + P2 * 4 - P3) * (u * u) / 2 +
                    (P0 * (-1) + P2) * u / 2 +
                    P1);
        }
    };
 
}   // namespace easy3d
 
 
#endif  // EASY3D_CORE_CURVE_H