cubic_spline.hh

#ifndef VTESTBED_GEOMETRY_CUBIC_SPLINE_HH
#define VTESTBED_GEOMETRY_CUBIC_SPLINE_HH

#include <vector>

#include <vtestbed/base/blitz.hh>
#include <vtestbed/geometry/curve_segments.hh>
#include <vtestbed/geometry/math.hh>

namespace vtb {

    namespace geometry {

        template <class V, class T>
        inline V
        difference(const V& a, const V& b, T delta) {
            V result = a-b;
            if (delta > 0) {
                result /= delta;
            }
            return result;
        }

        template <class V>
        inline V
        difference(const V& a, const V& b) {
            auto delta = distance(a,b);
            V result = a-b;
            if (delta > 0) {
                result /= delta;
            }
            return result;
        }

        template <class T, int N=3>
        class Cubic_spline {

        public:
            using value_type = Vertex<T,N>;
            using reference = value_type&;
            using const_reference = const value_type&;
            using vertex_array = std::vector<value_type>;
            using size_type = size_t;
            using segments_type = Curve_segments<T,N>;

        private:
            vertex_array _points;
            segments_type _segments;
            vertex_array _derivatives;

        public:

            inline explicit
            Cubic_spline(const vertex_array& values):
            _points(values), _segments(values) {
                const auto n = values.size();
                if (n <= 1) { return; }
                const auto& p = this->_points;
                auto& dt = this->_segments;
                auto& m = this->_derivatives;
                m.reserve(n);
                m.emplace_back(difference(p[1],p[0],dt[0]));
                for (size_t i=1; i<n-1; ++i) {
                    m.emplace_back(T{0.5}*(
                        difference(p[i+1],p[i],dt[i]) +
                        difference(p[i],p[i-1],dt[i-1])
                    ));
                }
                m.emplace_back(difference(p[n-1],p[n-2],dt[n-2]));
                dt.normalise();
            }

            inline value_type
            operator()(T t) const {
                const auto& segments = this->_segments;
                const auto& points = this->_points;
                const auto& derivatives = this->_derivatives;
                if (points.size() == 1) {
                    return points.front();
                }
                size_type i0 = 0;
                segments.locate(t, i0);
                if (i0 == segments.size()) {
                    return points.back();
                }
                size_type i1 = i0+1;
                const auto& p0 = points[i0];
                const auto& p1 = points[i1];
                const auto& m0 = derivatives[i0];
                const auto& m1 = derivatives[i1];
                T t1 = (T{1} - t);
                t1 *= t1;
                T t2 = t*t;
                T h00 = (T{1} + T{2}*t)*t1;
                T h10 = t*t1;
                T h01 = t2*(T{3} - T{2}*t);
                T h11 = t2*(t - T{1});
                return h00*p0 + h10*m0 + h01*p1 + h11*m1;
            }

            inline const_reference
            front() const noexcept {
                return this->_points.front();
            }

            inline const_reference
            back() const noexcept {
                return this->_points.back();
            }

            inline size_type
            size() const {
                return this->_points.size();
            }

            inline const vertex_array&
            derivatives() const {
                return this->_derivatives;
            }

            inline const segments_type&
            segments() const {
                return this->_segments;
            }

            inline segments_type&
            segments() {
                return this->_segments;
            }

            inline void
            range(T min, T max) {
                this->_segments.range(min, max);
            }

            inline T
            length() const {
                return this->_segments.length();
            }

            inline void
            smooth_left(value_type p_left) {
                auto& m = this->_derivatives;
                const auto& p = this->_points;
                m[0] = T{0.5}*(difference(p[1],p[0]) + difference(p[0], p_left));
            }

            inline void
            smooth_right(value_type p_right) {
                auto& m = this->_derivatives;
                const auto& p = this->_points;
                const auto n = p.size();
                m[n-1] = T{0.5}*(difference(p[n-1],p[n-2]) + difference(p[n-1], p_right));
            }

        };

    }

}

#endif // vim:filetype=cpp