spline.h

Go to the documentation of this file.

 
 
#pragma once
 
#include <scimesh/types.h>
#include <glm/geometric.hpp>
 
#include <algorithm>
#include <cmath>
#include <vector>
 
namespace scimesh {
 
namespace detail {
 
constexpr size_t kMinPathPoints = 2u;
 
constexpr double kMinKnotInterval = 1e-6;
 
constexpr double kArcLengthEpsilon = 1e-9;
 
inline Vec3 hermite_segment(const Vec3 &p0, const Vec3 &p1, const Vec3 &m0,
                            const Vec3 &m1, double h, double s) {
    const double s2 = s * s;
    const double s3 = s2 * s;
 
    const double w00 = 2.0 * s3 - 3.0 * s2 + 1.0;
    const double w10 = s3 - 2.0 * s2 + s;
    const double w01 = -2.0 * s3 + 3.0 * s2;
    const double w11 = s3 - s2;
 
    const Vec3 ends = p0 * static_cast<float>(w00) + p1 * static_cast<float>(w01);
    const Vec3 tangents = m0 * static_cast<float>(w10 * h) +
                          m1 * static_cast<float>(w11 * h);
    return ends + tangents;
}
 
inline size_t segment_count(size_t num_points, bool closed) {
    if (num_points < kMinPathPoints) {
        return 0u;
    }
    return closed ? num_points : num_points - 1u;
}
 
inline size_t next_index(size_t i, size_t num_points, bool closed) {
    const size_t n = i + 1u;
    if (n < num_points) {
        return n;
    }
    return closed ? 0u : num_points - 1u;
}
 
inline int sample_count(int samples_per_segment) {
    return std::max(1, samples_per_segment);
}
 
} // namespace detail
 
// ---------------------------------------------------------------------------
//  Path queries
// ---------------------------------------------------------------------------
 
inline float path_length(const std::vector<Vec3> &path, bool closed = false) {
    const size_t num_segments = detail::segment_count(path.size(), closed);
    double total = 0.0;
    for (size_t i = 0; i < num_segments; ++i) {
        total += glm::length(path[detail::next_index(i, path.size(), closed)] -
                             path[i]);
    }
    return static_cast<float>(total);
}
 
inline std::vector<Vec3> remove_duplicate_points(const std::vector<Vec3> &path,
                                                 float epsilon = 1e-6f) {
    std::vector<Vec3> cleaned;
    cleaned.reserve(path.size());
    for (const Vec3 &p : path) {
        if (cleaned.empty() || glm::length(p - cleaned.back()) > epsilon) {
            cleaned.push_back(p);
        }
    }
    return cleaned;
}
 
inline std::vector<Vec3> path_tangents(const std::vector<Vec3> &path,
                                       bool closed = false) {
    std::vector<Vec3> tangents;
    const size_t n = path.size();
    if (static_cast<int>(n) < static_cast<int>(detail::kMinPathPoints)) {
        return tangents;
    }
    tangents.resize(n);
 
    for (size_t i = 0; i < n; ++i) {
        Vec3 dir;
        if (closed) {
            const Vec3 prev = path[(i + n - 1u) % n];
            const Vec3 next = path[(i + 1u) % n];
            dir = next - prev;
        } else if (i == 0u) {
            dir = path[1] - path[0];
        } else if (i + 1u == n) {
            dir = path[n - 1u] - path[n - 2u];
        } else {
            dir = path[i + 1u] - path[i - 1u];
        }
 
        const float len = glm::length(dir);
        if (len > 1e-8f) {
            tangents[i] = dir / len;
        } else if (i > 0u) {
            tangents[i] = tangents[i - 1u];
        } else {
            tangents[i] = Vec3(0.0f, 0.0f, 1.0f);
        }
    }
    return tangents;
}
 
inline std::vector<float> path_curvature(const std::vector<Vec3> &path,
                                         bool closed = false) {
    const size_t n = path.size();
    if (n < 3u) {
        return std::vector<float>();
    }
 
    std::vector<float> curvature(n, 0.0f);
    for (size_t i = 0; i < n; ++i) {
        if (!closed && (i == 0u || i + 1u == n)) {
            continue;   // filled in below from the neighbouring interior point.
        }
        const Vec3 prev = path[closed ? (i + n - 1u) % n : i - 1u];
        const Vec3 next = path[closed ? (i + 1u) % n : i + 1u];
 
        const Vec3 v = next - prev;                       // ~ x' / (2h)
        const Vec3 a = next - 2.0f * path[i] + prev;      // ~ x'' / h^2
        const float v_len = glm::length(v);
        if (v_len < 1e-8f) {
            continue;
        }
        const float cross_len = glm::length(glm::cross(v, a));
        curvature[i] = (4.0f * cross_len) /
                       (v_len * v_len * v_len);
    }
 
    if (!closed && n >= 3u) {
        curvature[0] = curvature[1];
        curvature[n - 1u] = curvature[n - 2u];
    }
    return curvature;
}
 
// ---------------------------------------------------------------------------
//  Resampling
// ---------------------------------------------------------------------------
 
inline std::vector<Vec3> resample_by_arclength(const std::vector<Vec3> &path,
                                               float step,
                                               bool closed = false) {
    const size_t n = path.size();
    if (step <= 0.0f || n < 2u) {
        return path;
    }
 
    const double total = path_length(path, closed);
    if (total <= 0.0) {
        return std::vector<Vec3>(1u, path[0]);
    }
 
    // A closed loop is resampled at the step that covers it in a whole number
    // of steps, so that the seam is a regular step like every other one.
    double step_used = static_cast<double>(step);
    if (closed) {
        const double steps = std::max(
            3.0, std::round(total / static_cast<double>(step)));
        step_used = total / steps;
    }
 
    std::vector<Vec3> out;
    out.reserve(static_cast<size_t>(std::ceil(total / step_used)) + 2u);
    out.push_back(path[0]);
 
    const double kStep = step_used;
    const double eps = detail::kArcLengthEpsilon * total;
    const size_t num_segments = detail::segment_count(n, closed);
 
    double travelled = 0.0;   // arc length consumed before the current segment
    double next_at = kStep;   // arc length of the next sample
    for (size_t i = 0; i < num_segments; ++i) {
        const Vec3 a = path[i];
        const Vec3 b = path[detail::next_index(i, n, closed)];
        const double seg_len = glm::length(b - a);
 
        if (seg_len > 0.0) {
            while (next_at <= travelled + seg_len && next_at < total - eps) {
                const double u = (next_at - travelled) / seg_len;
                out.push_back(a + static_cast<float>(u) * (b - a));
                next_at += kStep;
            }
        }
        travelled += seg_len;
    }
 
    // The end of an open path is part of the shape, even when it does not fall
    // on a sample boundary; a closed path closes on its starting point.
    out.push_back(closed ? path[0] : path[n - 1u]);
    return out;
}
 
// ---------------------------------------------------------------------------
//  Curve sampling
// ---------------------------------------------------------------------------
 
inline std::vector<Vec3> hermite_path(const std::vector<Vec3> &points,
                                      const std::vector<Vec3> &tangents,
                                      int samples_per_segment = 8,
                                      bool closed = false) {
    std::vector<Vec3> path;
    const std::vector<Vec3> pts = remove_duplicate_points(points);
    if (pts.size() != points.size() || tangents.size() != points.size()) {
        return path;   // caller passed degenerate points or mismatched sizes.
    }
    if (pts.size() < detail::kMinPathPoints ||
        (closed && pts.size() < 3u)) {
        return path;
    }
 
    const int samples = detail::sample_count(samples_per_segment);
    const size_t num_segments = detail::segment_count(pts.size(), closed);
    path.reserve(num_segments * static_cast<size_t>(samples) + 1u);
 
    // Each segment spans one unit of the curve parameter, so h == 1 and the
    // caller's tangents are used as-is.
    for (size_t i = 0; i < num_segments; ++i) {
        const size_t j = detail::next_index(i, pts.size(), closed);
        for (int k = 0; k < samples; ++k) {
            const double s = static_cast<double>(k) / static_cast<double>(samples);
            path.push_back(detail::hermite_segment(pts[i], pts[j], tangents[i],
                                                   tangents[j], 1.0, s));
        }
    }
    if (closed) {
        path.push_back(pts[0]);   // the loop closes on the first point.
    } else {
        path.push_back(pts.back());
    }
    return path;
}
 
inline std::vector<Vec3> catmull_rom_path(const std::vector<Vec3> &points,
                                          int samples_per_segment = 8,
                                          bool closed = false,
                                          float alpha = 0.5f) {
    std::vector<Vec3> path;
    const std::vector<Vec3> pts = remove_duplicate_points(points);
    const size_t n = pts.size();
    if (n < detail::kMinPathPoints || (closed && n < 3u)) {
        return path;
    }
 
    // Knot intervals: the parameter distance between two neighbouring points.
    // For alpha = 0 every interval is 1 (the uniform case), otherwise it grows
    // with the chord length to the power of alpha.
    std::vector<double> interval(n, 1.0);
    for (size_t i = 0; i < n; ++i) {
        const size_t j = detail::next_index(i, n, closed);
        if (!closed && j == i) {
            break;   // open curve: the last point has no outgoing segment.
        }
        const double chord = glm::length(pts[j] - pts[i]);
        interval[i] = std::max(
            (alpha == 0.0f) ? 1.0 : std::pow(chord, static_cast<double>(alpha)),
            detail::kMinKnotInterval);
    }
 
    // Cumulative knot values along the curve.
    std::vector<double> knot(n, 0.0);
    for (size_t i = 1; i < n; ++i) {
        knot[i] = knot[i - 1u] + interval[i - 1u];
    }
 
    // Tangent at each control point: the central difference of its neighbours,
    // divided by their knot distance (which is what makes the formula work for
    // the non-uniform parameterizations).
    std::vector<Vec3> tangents(n);
    for (size_t i = 0; i < n; ++i) {
        Vec3 prev, next;
        double t_prev, t_next;
        if (closed) {
            prev = pts[(i + n - 1u) % n];
            next = pts[(i + 1u) % n];
            t_prev = knot[i] - interval[(i + n - 1u) % n];
            t_next = knot[i] + interval[i];
        } else if (i == 0u) {
            // Reflected phantom point: keeps the curve starting at pts[0] with
            // the tangent direction of the first segment.
            prev = 2.0f * pts[0] - pts[1];
            next = pts[1];
            t_prev = knot[0] - interval[0];
            t_next = knot[1];
        } else if (i + 1u == n) {
            next = 2.0f * pts[n - 1u] - pts[n - 2u];
            prev = pts[n - 2u];
            t_prev = knot[n - 2u];
            t_next = knot[n - 1u] + interval[n - 2u];
        } else {
            prev = pts[i - 1u];
            next = pts[i + 1u];
            t_prev = knot[i - 1u];
            t_next = knot[i + 1u];
        }
 
        const double span = t_next - t_prev;
        if (span > detail::kMinKnotInterval) {
            tangents[i] = (next - prev) / static_cast<float>(span);
        }
    }
 
    const int samples = detail::sample_count(samples_per_segment);
    const size_t num_segments = detail::segment_count(n, closed);
    path.reserve(num_segments * static_cast<size_t>(samples) + 1u);
 
    for (size_t i = 0; i < num_segments; ++i) {
        const size_t j = detail::next_index(i, n, closed);
        const double h = (j == 0u) ? (interval[n - 1u]) : (knot[j] - knot[i]);
        for (int k = 0; k < samples; ++k) {
            const double s = static_cast<double>(k) / static_cast<double>(samples);
            path.push_back(detail::hermite_segment(pts[i], pts[j], tangents[i],
                                                   tangents[j], h, s));
        }
    }
    path.push_back(closed ? pts[0] : pts[n - 1u]);
    return path;
}
 
inline std::vector<Vec3> bspline_path(const std::vector<Vec3> &points,
                                      int samples_per_segment = 8,
                                      bool closed = false) {
    std::vector<Vec3> path;
    const std::vector<Vec3> pts = remove_duplicate_points(points);
    const size_t n = pts.size();
    if (n < 4u) {
        return path;   // a cubic B-spline needs four control points per piece.
    }
 
    const int samples = detail::sample_count(samples_per_segment);
 
    // Uniform cubic B-spline basis, as a function of the local parameter s in
    // [0, 1] within one knot span.
    const auto basis = [](double s, double &b0, double &b1, double &b2, double &b3) {
        const double s2 = s * s;
        const double s3 = s2 * s;
        b0 = (1.0 - 3.0 * s + 3.0 * s2 - s3) / 6.0;
        b1 = (4.0 - 6.0 * s2 + 3.0 * s3) / 6.0;
        b2 = (1.0 + 3.0 * s + 3.0 * s2 - 3.0 * s3) / 6.0;
        b3 = s3 / 6.0;
    };
 
    // Segment `i` is controlled by the points i-3 .. i (wrapped for a closed
    // curve): a uniform cubic B-spline is defined on the knot span [3, n] of
    // the uniform knot vector, which is n - 3 segments for an open curve and n
    // for a closed one.  For the open curve, segment i therefore uses the
    // points i .. i+3 directly.
    const size_t num_segments = closed ? n : (n - 3u);
    path.reserve(num_segments * static_cast<size_t>(samples) + 1u);
 
    for (size_t i = 0; i < num_segments; ++i) {
        const Vec3 &p0 = pts[closed ? (i + n - 3u) % n : i];
        const Vec3 &p1 = pts[closed ? (i + n - 2u) % n : i + 1u];
        const Vec3 &p2 = pts[closed ? (i + n - 1u) % n : i + 2u];
        const Vec3 &p3 = pts[closed ? i : i + 3u];
 
        // The last segment of an open curve also emits its endpoint (s = 1),
        // so that the path ends where the curve does.
        const int emit = (!closed && i + 1u == num_segments) ? samples + 1 : samples;
        for (int k = 0; k < emit; ++k) {
            const double s = static_cast<double>(k) / static_cast<double>(samples);
            double b0, b1, b2, b3;
            basis(s, b0, b1, b2, b3);
            path.push_back(static_cast<float>(b0) * p0 +
                           static_cast<float>(b1) * p1 +
                           static_cast<float>(b2) * p2 +
                           static_cast<float>(b3) * p3);
        }
    }
 
    // A periodic B-spline returns to its starting point after n knot spans, so
    // the loop closes on the first sample.
    if (closed) {
        path.push_back(path[0]);
    }
    return path;
}
 
inline std::vector<Vec3> bezier_path(const std::vector<Vec3> &control_points,
                                     int samples = 64) {
    std::vector<Vec3> path;
    const size_t n = control_points.size();
    if (n < 2u) {
        return path;
    }
 
    const int count = std::max(2, samples);
    path.reserve(static_cast<size_t>(count));
 
    std::vector<Vec3> work(n);
    for (int k = 0; k < count; ++k) {
        // De Casteljau: repeated linear interpolation between neighbouring
        // points converges to the curve point for the current parameter.
        const float u = static_cast<float>(k) / static_cast<float>(count - 1);
        work = control_points;
        for (size_t level = n - 1u; level > 0u; --level) {
            for (size_t i = 0; i < level; ++i) {
                work[i] = (1.0f - u) * work[i] + u * work[i + 1u];
            }
        }
        path.push_back(work[0]);
    }
    return path;
}
 
} // namespace scimesh