The London Perl and Raku Workshop takes place on 26th Oct 2024. If your company depends on Perl, please consider sponsoring and/or attending.

NAME

Math::BSpline::Basis - B-spline basis functions

VERSION

version 0.002

SYNOPSIS

    use Math::BSpline::Basis;

    my $bspline = Math::BSpline::Basis->new(
        degree      => 3,
        knot_vector => [0, 0, 0, 0, 1, 1, 1, 1],
    );
    my $P = [-1, 2, 5, -3];

    my $u = 0.6;
    my $s = $bspline->find_knot_span($u);
    my $N = $bspline->evaluate_basis_functions($s, $u);

    my $value = 0;
    my $p     = $bspline->degree;
    for (my $i=0;$i<=$p;$i++) {
        $value += $N->[$i] * $P->[$s-$p+$i];
    }

    my $D = $bspline->evaluate_basis_derivatives($s, $u, 3);

DESCRIPTION

Introduction

A spline S of degree p is a piecewise polynomial function from a real interval [a, b] to the real numbers. This means that given a strictly increasing sequence of breakpoints a = u_0 < u_1 < ... < u_m = b, S is a polynomial function of degree less or equal to p on each interval [u_i, u_{i+1}].

At the breakpoints u_i, S can be discontinuous or continuously differentiable to a certain degree. If we specify minimum continuity conditions for each breakpoint then the set of splines of degree p satisfying these conditions forms a real vector space.

Given the degree, sequence of breakpoints, and continuity conditions, one can define a certain set of splines with convenient properties (local support, partition of unity etc) that form a basis of the above vector space. These are called basis splines or B-splines.

It turns out that the breakpoints and continuity conditions can be neatly specified in one go by defining a non-decreasing sequence of knots (also called knot vector) a = u_0 <= u_1 <= ... <= u_m = b. Note the difference to the breakpoints above: while the breakpoints are distinct, the knots only need to be non-decreasing. The breakpoints correspond to the set of distinct knots and the continuity conditions are encoded in the multiplicity of breakpoints. If a breakpoint has multiplicity k then the B-splines are p-k times continuously differentiable at this breakpoint.

B-Spline Support in Math::BSpline::Basis

As discussed above, a degree p and a knot vector U determine a vector space of splines and a basis of B-splines for this vector space. This module allows you to evaluate the non-zero basis functions and derivatives at a given position u in [a, b].

An important property of B-splines is that on a given knot span [u_i, u_{i+1}[ only a small number of B-splines is non-zero. This property is called local support. More specifically, the B-spline function N_{i,p} is equal to 0 outside the interval [u_i, u_{i+p+1}]. This implies that on a given knot span [u_i, u_{i+1}[, at most the p+1 B-splines N_{i-p,p},...,N_{p,p} can be non-zero.

Hence, the first step in evaluating B-splines at a given position is to determine the knot span the position falls into. After that, it is possible to evaluate specifically the non-zero basis functions.

    my $u = 0.6;
    my $s = $bspline->find_knot_span($u);
    my $N = $bspline->evaluate_basis_functions($s, $u);

Special attention has to be given to the first and the last breakpoint. A knot vector is called clamped (or non-periodic or open) if the first and the last breakpoint each have multiplicity p+1. This corresponds to no continuity restrictions at the boundaries of the interval [a, b]. This module only supports clamped knot vectors.

For a clamped knot vector of length m+1 (u_0,...,u_m), the dimension of the vector space of splines is n+1 with n = m - p - 1, hence we have basis splines N_{0,p},...,N_{n,p}. In order to form a B-spline function as a linear combination of the basis splines, we need n+1 so-called control points P_0,...,P_n. (The control points can also be vectors, which results in a B-spline curve.) In order to evaluate the resulting B-spline function (or curve) at a given position u we need to evaluate the non-vanishing basis functions as described above and then form the linear combination making sure that we get the indexing right:

    my $u = 0.6;
    my $s = $bspline->find_knot_span($u);
    my $N = $bspline->evaluate_basis_functions($s, $u);

    my $value = 0;
    for (my $i=0;$i<=$p;$i++) {
        $value += $N->[$i] * $P->[$s-$p+$i];
    }

Have a look at Math::BSpline::Curve if you would like to work with B-spline curves (or even functions).

CONSTRUCTORS

new

    $bspline = Math::BSpline::Basis->new(
        degree      => 3,
        knot_vector => [0, 0, 0, 0, 1, 1, 1, 1],
    );

Parameters:

degree (mandatory)

The degree of the B-splines.

knot_vector

The knot vector as array reference. A full clamped knot vector must be a sequence of non-decreasing numbers with p+1 copies of the same number at the beginning and p+1 copies of the same number at the end. In order to achieve a valid knot vector, some automatic trimming is applied on a copy of the knot vector while the original value remains unchanged. However, no fool-proof validation is performed, specifically, the knot vector is not validated against undefined or non-numeric values.

sorting

In order to protect against decreasing knot values (e.g. due to rounding errors in upstream calculations) the knot vector is sorted.

clamping

Mostly to relieve the user of this task, the knot vector is extended on both ends, such that the first p+1 values are identical and the last p+1 values are identical.

breakpoint multiplicity

The multiplicity of the first and last breakpoint is limited to p+1, the multiplicity of each internal breakpoint is limited to p. Increasing breakpoint multiplicity any further would lead to basis functions that are discontinous at the respective breakpoint. This is currently not supported.

If not specified at all, the knot vector defaults to [0,...,0, 1,...,1] with p+1 copies each. This results in a Bezier spline of degree p.

ATTRIBUTES

degree

The degree of the B-splines. Must be set at construction and cannot be modified afterwards.

knot_vector

The knot vector of the B-splines. Can only be set at construction. Defaults to [0,...,0, 1,...,1] resulting in a Bezier spline.

METHODS

find_knot_span

  $s = $bspline->find_knot_span($u)

If the knot vector is (u_0,...u_m) then this method returns the index s such that u_s <= u < u_{s+1}. A special case is u = u_m. In this case, the method returns s such that u_s < u <= u_{s+1}, i.e. u is the upper boundary of the knot span.

Outside of [u_p, u_{m-p}], all basis functions vanish. This is supported by evaluate_basis_functions and evaluate_basis_derivatives. This method returns p for u < u_p and m-p-1 for u > u_{m-p}. If you use the return value exclusively for feeding it into evaluate_basis_functions or evaluate_basis_derivatives then you do not need to care about this.

evaluate_basis_functions

    $P = [-1, 2, 5, -3];
    $u = 0.6;
    $s = $bspline->find_knot_span($u);
    $N = $bspline->evaluate_basis_functions($s, $u);

    $value = 0;
    $p     = $bspline->degree;
    for (my $i=0;$i<=$p;$i++) {
        $value += $N->[$i] * $P->[$s-$p+$i];
    }

Expects a knot span index s as returned by find_knot_span and a parameter value u and returns the values of the p+1 (potentially) non-vanishing B-splines N_{s-p,p},...,N_{p,p} as an array reference.

evaluate_basis_derivatives

    $P = [-1, 2, 5, -3];
    $u = 0.6;
    $d = 3;
    $s = $bspline->find_knot_span($u);
    $D = $bspline->evaluate_basis_derivatives($s, $u, $d);

    $value = [];
    $p     = $bspline->degree;
    for (my $k=0;$k<=$d;$k++) {
        $value->[$k] = 0;

        if ($k <= $p) {
            for (my $i=0;$i<=$p;$i++) {
                $value->[$k] += $D->[$k]->[$i] * $P->[$s-$p+$i];
            }
        }
    }

ACKNOWLEDGEMENTS

This implementation is based on the theory and algorithms presented in the NURBS book [1] and (to a much lesser degree) on the theory presented in [2].

[1] Piegl, L., Tiller, W.: The NURBS book, 2nd Edition. Springer, 1997.
[2] de Boor, C.: A Practical Guide to Splines, Revised Edition. Springer 2001.

AUTHOR

Lutz Gehlen <perl@lutzgehlen.de>

COPYRIGHT AND LICENSE

This software is copyright (c) 2020 by Lutz Gehlen.

This is free software; you can redistribute it and/or modify it under the same terms as the Perl 5 programming language system itself.