Nanospline

Nanospline is a header-only spline library written with modern C++. It is created by Qingnan Zhou as a coding exercise. It supports Bézier, rational Bézier, B-spline and NURBS curves of arbitrary degree in arbitrary dimension. Most of the algorithms are covered by The NURBS Book.

Functionalities

The following functionalities are covered:

Data structure
Creation
Basic information
Evaluation
Derivatives
Curvature
Hodograph
Inverse evaluation
Knot insertion and removal
Split
Degree elevation
Arc length
Inflection
Turning angle
Singularity
Conversion
Intersection

Data structure

Nanospline provide 4 basic data structures for the 4 types of curves: Bézier, rational Bézier, B-spline and NURBS. All of them are templated by 4 parameters:

Scalar: The floating point data type. (e.g. float, double, long double, etc.)
dim: The dimension of the embedding space. (e.g. 2 for 2D curves, and 3 for 3D curves.)
degree: The degree of the curve. (e.g. 2 for quadratic curves, 3 for cubic curves.) The special value -1 is used to indicate dynamic degree.
generic: (Optional) This is a boolean flag indicating whether to treat the curve as a generic curve. Nanospline sometimes provides specialized implementation when the degree of the curve is known. By setting generic to true, we are forcing nanospline to use the default general implementation. If degree is -1 (i.e. dynamic degree), generic should always be true.

Creation

All 4 types of curves can be constructed in the following pattern:

CurveType<Scalar, dim, degree, generic> curve;

where CurveType is one of the following: Bezier, BSpline, RationalBezier and NURBS. Different curve type requires setting different fields:

Fields	Bézier	B-spline	Rational Bézier	NURBS
Control points	Yes	Yes	Yes	Yes
Knots	No	Yes	No	Yes
Weights	No	No	Yes	Yes

All fields are represented using Eigen::Matrix types.

Bézier curve

#include <nanospline/Bezier.h>

// Construct a 2D cubic Bézier curve
nanospline::Bezier<double, 2, 3> curve;

// Setting control points. Assuming `ctrl_pts` is a 4x2 Eigen matrix.
curve.set_control_points(ctrl_pts);

B-spline curve

#include <nanospline/BSpline.h>

// Construct a 2D cubic B-spline curve
nanospline::BSpline<double, 2, 3> curve;

// Setting control points. Assuming `ctrl_pts` is a nx2 Eigen matrix.
curve.set_control_points(ctrl_pts);

// Setting knots.  Assuming `knots` is a mx1 Eigen matrix.
// Where m = n+p+1, and `p` is the degree of the curve.
curve.set_knots(knots);

Rational Bézier curve

#include <nanospline/RationalBezier.h>

// Construct a 2D cubic rational Bézier curve
nanospline::RationalBezier<double, 2, 3> curve;

// Setting control points. Assuming `ctrl_pts` is a 4x2 Eigen matrix.
curve.set_control_points(ctrl_pts);

// Setting weights. Assuming `weights` is a 4x1 Eigen matrix.
curve.set_weights(weights);

// **Important**: RationalBezier requires initialization.
curve.initialize();

NURBS curve

#include <nanospline/NURBS.h>

// Construct a 2D cubic NURBS curve
nanospline::NURBS<double, 2, 3> curve;

// Setting control points. Assuming `ctrl_pts` is a 4x2 Eigen matrix.
curve.set_control_points(ctrl_pts);

// Setting knots.  Assuming `knots` is a mx1 Eigen matrix.
// Where m = n+p+1, and `p` is the degree of the curve.
curve.set_knots(knots);

// Setting weights. Assuming `weights` is a 4x1 Eigen matrix.
curve.set_weights(weights);

// **Important**: NURBS requires initialization.
curve.initialize();

Basic information

Once a curve is initialized, one can query for a number of basic information:

// Polynomial degree.
int degree = curve.get_degree();

// The minimum and maximum parameter value.
auto t_min = curve.get_domain_lower_bound();
auto t_max = curve.get_domain_upper_bound();

// Control points
const auto& ctrl_pts = curve.get_control_points();

// Knots (BSpline and NURBS only).
const auto& knots = curve.get_knots();

// Weights (RationalBezier and NURBS only).
const auto& weights = curve.get_weights();

Evaluation

One can retrieve the point, p, corresponding to a given parameter value, t, by evaluating the curve:

auto p = curve.evaluate(t);

Derivatives

One can compute the first and second derivative vectors, d1 and d2 respectively, at a given parameter value, t, by evaluating the curve derivatives as the following:

auto d1 = curve.evaluate_derivative(t);
auto d2 = curve.evaluate_2nd_derivative(t);

Curvature

Nanospline also provides direct support for computing curvature vector, k, at a given parameter value, t, using the following:

auto k = curve.evaluate_curvature(t);

Hodograph

For Bézier and B-spline, it is well known that their derivative curves are also Bézier or B-spline curves respectively. The first derivative curve is called hodograph, which can be computed with the following:

#include <nanospline/hodograph.h>

auto hodograph = nanospline::compute_hodograph(curve);

Higher order derivative curves can be obtained by computing the hodograph of a hodograph recursively.

Inverse evaluation

Inverse evaluation aims to find a point, parameterized by t, on a given curve that is closest to a query point, q. In the most general case, inverse evaluation can be reduced to finding roots of high degree polynomial. However, closed formula often exist for lower degree curves. Nanospline provide two different methods for inverse evaluation:

// Method 1
auto t = curve.inverse_evaluate(q);

// Method 2
auto t = curve.approximate_inverse_evaluate(q);

// Method 2 complete signature
auto t = curve.approximate_inverse_evaluate(q, t_min, t_max, level);

The first method, inverse_evaluate(), tries to find the exact closest point by solving high degree polynomial or apply closed formula for low degree curves. It is currently under development.

In contrast, the second method, approximate_inverse_evaluate(), uses brute force bisection method to find an approximate closest point. The parameter t_min and t_max specify the search domain, and level is the recursion level. Higher recursion level provides more accurate result.

Knot insertion and removal

For B-spline and NURBS curves, one can insert extra knot, t, with multiplicity, m, using the following:

curve.insert_knot(t, m);

To remove a knot m times:

curve.remove_knot(t, m);

Just to be complete, the full signature of the remove_knot function is

int num_removed = curve.remove_knot(t, m, tolerance);

Where the return value num_removed indicates how many times the knot is removed, and tolerance specifies the max allowed change (L2 distance) in the curve that a valid removal can introduce.

Split

To split a curve into two halves at parameter value t:

#include <nanospline/split.h>

auto halves = nanospline::split(curve, t);

One can also split a B-Spline curve into a sequence of Bézier curves:

#include <nanospline/BSpline.h>

auto r = bspline.convert_to_Bezier();
const auto& bezier_segments = std::get<0>(r);
const auto& parameter_bounds = std::get<1>(r);

Where the ith Bézier curve covers the knot span [parameter_bounds[i], parameter_bounds[i+1]]. It is also possible recombine these Bézier segments to form a B-spline curve:

BSpline<Scalar, dim, degree, generic> curve(
    bezier_segments, parameter_bounds);

// or with uniform knot span

BSpline<Scalar, dim, degree, generic> curve(bezier_segments);

Similarly, NURBS curve can be split into a sequence of rational Bézier curves:

#include <nanospline/BSpline.h>

auto r = nurbs.convert_to_Bezier();
const auto& rational_bezier_segments = std::get<0>(r);
const auto& parameter_bounds = std::get<1>(r);

// Re-combine rational Bézier back into NURBS

NURBS<Scalar, dim, degree, generic> curve(
    rational_bezier_segments, parameter_bounds);

Degree elevation

It is often useful to increase the degree of a curve:

auto curve2 = curve.elevate_degree();
assert(curve2.get_degree() == curve.get_degree()+1);

Arc length

Given a parameter value, t, the following function computes the arc length l at t:

auto l = arc_length(curve, t);

On the other hand, given an arc length l, the following function find the parameter t corresponding to l:

auto t = inverse_arc_length(curve, l);

Inflection

Nanospline also supports computing 2D curve inflection points, i.e. points with zero curvature.

#include <nanospline/inflection.h>

auto inflections = nanospline::compute_inflections(curve);

Where inflections is a vector of parameter values corresponding to inflection points.

Turning angle

Turning angle is the total curvature of a given curve. It represents how much a curve is bending. In Nanospline, turning angle computation is supported for 2D curves:

auto turning_angle = curve.get_turning_angle(t0, t1);

which returns the turning angle (in radians) for the curve segment between t0 and t1. It is often important to determine the locations, tcs, where splitting the curve at tcs will reduce the turning angle by half for each curve piece:

auto turning_angle = curve.get_turning_angle(t0, t1);
std::vector<Scalar> tcs = curve.reduce_turning_angle(t0, t1);

if (!tcs.empty()) {
    auto turning_angle_0 = curve.get_turning_angle(t0, tcs.front());
    assert(std::abs(turning_angle * 0.5 - turning_angle_0) < EPS);
    ...
}

Singularity

Singularity points of a curve are defined as the locations at where the curve has 0 first derivative. Singularity locations between t0 and t1 for 2D curves can be computed in the following way:

std::vector<Scalar> singularites = curve.compute_singularities(t0, t1);

Conversion

It is easy to convert a Bézier curve into a BSpline curve:

#include <nanospline/conversion.h>

auto bspline = nanospline::convert_to_BSpline(bezier);

It is also possible to convert a BSpline curve into a sequence of Bézier curves:

#include <nanospline/conversion.h>

auto beziers = nanospline::convert_to_Bezier(bspline);
for (const auto& curve : bezier) {
    // `curve` is a Bezier<...> curve.
}

To convert a Bézier into rational Bézier curve:

#include <nanospline/conversion.h>

auto rationa_bezier = nanospline::convert_to_RationalBezier(bezier);

To convert a BSpline into NURBS curve:

#include <nanospline/conversion.h>

auto nurbs = nanospline::convert_to_nurbs(bspline);

To convert a NURBS curve into a number of rational Bézier curves:

#include <nanospline/conversion.h>

auto rationa_beziers = nanospline::convert_to_RationalBezier(nurbs);

It is sometimes possible to convert a rational Bézier curve into a plane Bézier curve, and convert a NURBS curves into a BSpline curve. Such conversion is allowed when the all weights are the same. Otherwise, an exception will be raised.

#include <nanospline/conversion.h>

try {
    auto bezier = nanospline::convert_to_Bezier(rational_bezier);
    auto bspline = nanospline::convert_to_BSpline(nurbs);
} catch (const std::exception& e) {
    ...
}

Intersection

Nanospline supports 3 forms of curve-patch intersection computation:

Generic curve-patch intersection

This is the most general form:

#include <nanospline/intersect.h>

Scalar t, u, v;
bool convergdd;

std::tie(t, u, v, converged) = nanospline::intersect(
    curve,           ///< Input curve.
    patch,           ///< Input patch.
    t0,              ///< Initial guess of the intersection on curve.
    u0, v0,          ///< Initial guess of the intersection on patch.
    num_iterations,  ///< Max num iterations.
    tolerance);      ///< Convergence tolerance.

The converged flag will be false if the intersection computation fails. Typical causes are either the input curve does not intersect the input patch or the initial guesses are too off that cause Newton iterations to diverge.

Curve-plane intersection

If the patch is actually a plane, Nanospline offers a faster intersection computation:

#include <nanospline/intersect.h>

std::array<Scalar, 4> plane; // Coefficients of the plane equation.
                             // i.e. [a,b,c,d] -> ax + by + cz + d = 0.
Scalar t;
bool converged;

std::tie(t, converged) = nanospline::intersect(
    curve,          ///< Input curve.
    plane,          ///< Plane eqquation coefficients.
    t0,             ///< Initial guess of the intersection on curve.
    num_iterations, ///< Max num iterations.
    tolerance);     ///< Convergence tolerance.

Embedded curve and plane intersection

Lastly, it is sometimes necessary to compute the intersection of a curve embedded in a patch that is the image of a straight line in parametric space and a plane:

#include <nanospline/intersect.h>

std::array<Scalar, 4> plane; // Coefficients of the plane equation.
                             // i.e. [a,b,c,d] -> ax + by + cz + d = 0.
Scalar u, v;
bool converged;

std::tie(u, v, converged) = nanospline::intersect(
    patch,          ///< The patch that contains the curve.
    pu, pv,         ///< The uv coordinate of the starting point.
    qu, qv,         ///< The uv coordinate of the end point.
    plane,          ///< Plane eqquation coefficients.
    u0, v0,         ///< Initial guess of the intersection on curve.
    num_iterations, ///< Max num iterations.
    tolerance);     ///< Convergence tolerance.

Cheap and reliable Node.js hosting starts at $3/month, and $1/month static HTML hosting

qnzhou / nanospline

Programming Languages