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, and3
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 settinggeneric
totrue
, we are forcing nanospline to use the default general implementation. Ifdegree
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 i
th 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.