basis_spline.hh

#ifndef VTESTBED_GEOMETRY_BASIS_SPLINE_HH
#define VTESTBED_GEOMETRY_BASIS_SPLINE_HH

#include <vtestbed/geometry/types.hh>
#include <vtestbed/linalg/linear_algebra.hh>

namespace vtb {

    namespace geometry {

        template <class T, int N>
        class Basis_spline {

        public:
            using point_type = Vertex<T,N>;
            using element_type = Vertex<T,N+1>;
            using knot_type = T;
            using weight_type = T;
            using point_array = std::vector<point_type>;
            using element_array = std::vector<element_type>;
            using knot_array = std::vector<knot_type>;
            using weight_array = std::vector<weight_type>;
            using coefficient_array = std::vector<T>;
            using scalar_type = T;
            using value_type = element_type;
            using reference = value_type&;
            using const_reference = const value_type&;
            using size_type = size_t;

        private:
            element_array _elements;
            knot_array _knots;
            T _x_min{0}, _x_max{1};

        public:

            Basis_spline() = default;
            ~Basis_spline() = default;
            Basis_spline(const Basis_spline&) = default;
            Basis_spline& operator=(const Basis_spline&) = default;
            Basis_spline(Basis_spline&&) = default;
            Basis_spline& operator=(Basis_spline&&) = default;

            inline explicit
            Basis_spline(element_array elements, knot_array knots):
            _elements(elements), _knots(knots) {}

            inline explicit
            Basis_spline(element_array elements, size_type degree=3):
            Basis_spline(elements, make_knots(elements.size(), degree)) {}

            inline explicit
            Basis_spline(point_array points, weight_array weights, knot_array knots):
            Basis_spline(make_elements(points, weights), knots) {}

            inline explicit
            Basis_spline(point_array points, knot_array knots):
            Basis_spline(make_elements(points), knots) {}

            inline explicit
            Basis_spline(point_array points, weight_array weights, size_type degree=3):
            Basis_spline(make_elements(points, weights),
                    make_knots(points.size(), degree)) {}

            inline explicit
            Basis_spline(point_array points, size_type degree=3):
            Basis_spline(make_elements(points), degree) {}

            inline size_type
            degree() const {
                const auto m = num_knots()-1;
                const auto n = num_points()-1;
                return m-n-1;
            }

            inline size_type
            num_internal_knots() const {
                const auto p = degree();
                const auto m = num_knots()-1;
                return m-p-1-(p+1)+1;
            }

            inline element_type
            operator()(T x) const {
                const auto& t = this->_knots;
                const auto& c = this->_elements;
                const auto p = this->degree();
                size_type k = 0;
                locate(x, k);
                if (k == t.size()) { return c.back(); }
                if (k >= c.size()) { k = c.size()-1; }
                element_array d(p+1);
                for (size_type i=0; i<p+1; ++i) {
                    d[i] = multiply_by_weight(c[i+k-p]);
                }
                for (size_type r=1; r<p+1; ++r) {
                    for (size_type j=p; j>r-1; --j) {
                        T t0 = t[j+k-p];
                        T t1 = t[j+1+k-r];
                        T alpha = (x - t0) / (t1 - t0);
                        d[j] = (T{1} - alpha)*d[j-1] + alpha*d[j];
                    }
                }
                return divide_by_weight(d[p]);
            }

            inline coefficient_array
            coefficients(T x) const {
                const auto& t = this->_knots;
                const auto& c = this->_elements;
                const auto p = this->degree();
                auto npoints = c.size();
                coefficient_array coefs(npoints);
                size_type k = 0;
                locate(x, k);
                if (k == 0) { coefs.front() = 1; return coefs; }
                if (k == t.size()) { coefs.back() = 1; return coefs; } // TODO ???
                coefs[k] = 1;
                for (size_type d=1; d<p+1; ++d) {
                    coefs[k-d] = ((t[k+1]-x) / (t[k+1]-t[k+1-d])) * coefs[k+1-d];
                    for (size_type i=k-d+1; i<k; ++i) {
                        coefs[i] = ((x-t[i]) / (t[i+d]-t[i]))*coefs[i] +
                            ((t[i+d+1]-x) / (t[i+d+1] - t[i+1])) * coefs[i+1];
                    }
                    coefs[k] = ((x-t[k]) / (t[k+d]-t[k]))*coefs[k];
                }
                return coefs;
            }

            inline void
            interpolate() {
                using namespace vtb::linalg;
                auto& c = this->_elements;
                int n = c.size();
                auto delta = T{1}/(n-1);
                for (int k=0; k<N; ++k) {
                    Square_matrix<T> lhs(blitz::shape(n,n));
                    for (int i=0; i<n; ++i) {
                        const auto t = i*delta;
                        auto coefs = coefficients(t);
                        for (int j=0; j<n; ++j) { lhs(i,j) = coefs[j]; }
                    }
                    Vector<T> rhs(n);
                    for (int i=0; i<n; ++i) { rhs(i) = c[i](k); }
                    rhs = solve(lhs, rhs);
                    for (int i=0; i<n; ++i) { c[i](k) = rhs(i); }
                }
            }

            inline size_type size() const { return this->num_points(); }
            inline size_type num_points() const { return this->_elements.size(); }
            inline size_type num_knots() const { return this->_knots.size(); }
            inline const knot_array& knots() const { return this->_knots; }
            inline const element_array& elements() const { return this->_elements; }

            inline void
            range(T min, T max) {
                this->_x_min = min;
                this->_x_max = max;
            }

            inline static knot_array
            make_knots(size_type npoints, size_type degree) {
                const auto n = npoints-1;
                if (degree > n) { degree = n; }
                const auto p = degree;
                const auto m = n + p + 1;
                if (m+1 < p+p) {
                    throw std::invalid_argument("too few points");
                }
                knot_array knots;
                knots.reserve(m+1);
                const auto mreal = m+1-p-p;
                T delta = T{1}/(mreal-1);
                for (size_type i=0; i<p; ++i) { knots.emplace_back(0); }
                for (size_type i=0; i<mreal; ++i) { knots.emplace_back(i*delta); }
                for (size_type i=0; i<p; ++i) { knots.emplace_back(1); }
                return knots;
            }

            inline static knot_array
            make_cyclic_knots(size_type npoints, size_type degree) {
                npoints += degree;
                const auto p = degree;
                const auto n = npoints-1;
                const auto m = n + p + 1;
                if (m+1 < p+p) {
                    throw std::invalid_argument("too few points");
                }
                knot_array knots;
                knots.reserve(m+1);
                const auto mreal = m+1-p-p;
                T delta = T{1}/(mreal-1);
                for (size_type i=0; i<p; ++i) { knots.emplace_back((p-i)*(-delta)); }
                for (size_type i=0; i<mreal; ++i) { knots.emplace_back(i*delta); }
                for (size_type i=0; i<p; ++i) { knots.emplace_back((mreal+i)*delta); }
                return knots;
            }

            inline static element_type
            to_element(const point_type& p, T weight) {
                element_type elem;
                for (int i=0; i<N; ++i) { elem(i) = p(i); }
                elem(N) = weight;
                return elem;
            }

            inline static point_type
            to_point(const element_type& elem) {
                point_type point;
                for (int i=0; i<N; ++i) { point(i) = elem(i); }
                return point;
            }

        private:

            inline void
            locate(T& x, size_type& i) const {
                x = _x_min + (_x_max - _x_min)*x;
                const auto& knots = this->_knots;
                const auto nknots = knots.size();
                for (i=0; i<nknots; ++i) {
                    if (x < knots[i]) {
                        if (i > 0) { --i; }
                        return;
                    }
                }
            }

            inline static element_array
            make_elements(const point_array& points, const weight_array& weights) {
                if (points.size() != weights.size()) {
                    throw std::invalid_argument("bad size");
                }
                const auto npoints = points.size();
                element_array elems;
                if (npoints == 0) { return elems; }
                elems.reserve(npoints);
                for (size_type i=0; i<npoints; ++i) {
                    elems.emplace_back(to_element(points[i], weights[i]));
                }
                return elems;
            }

            inline static element_array
            make_elements(const point_array& points) {
                const auto npoints = points.size();
                element_array elems;
                if (npoints == 0) { return elems; }
                elems.reserve(npoints);
                for (size_type i=0; i<npoints; ++i) {
                    elems.emplace_back(to_element(points[i], T{1}));
                }
                return elems;
            }

            inline static element_type
            multiply_by_weight(element_type elem) {
                if (elem(N) != T{}) {
                    for (int i=0; i<N; ++i) { elem(i) *= elem(N); }
                }
                return elem;
            }

            inline static element_type
            divide_by_weight(element_type elem) {
                if (elem(N) != T{}) {
                    for (int i=0; i<N; ++i) { elem(i) /= elem(N); }
                }
                return elem;
            }

        };

        template <class T, int N>
        class Closed_basis_spline: public Basis_spline<T,N> {

        private:
            using basis_spline = Basis_spline<T,N>;

        public:
            using typename basis_spline::value_type;
            using typename basis_spline::scalar_type;
            using typename basis_spline::size_type;
            using typename basis_spline::element_array;
            using typename basis_spline::knot_array;

        public:
            inline explicit
            Closed_basis_spline(element_array elements, size_type degree=3):
            basis_spline(wrap(elements, degree),
                    basis_spline::make_cyclic_knots(elements.size(), degree)) {}

        private:

            inline element_array
            wrap(element_array c, size_type degree) {
                for (size_type i=0; i<degree; ++i) { c.emplace_back(c[i]); }
                return c;
            }

        };

    }

}

#endif // vim:filetype=cpp