A Computational Method with Maple for Finding the Maximum Curvature of a Bézier-Spline Curve

Abstract

In this paper, we propose two Maple procedures and some related utilities to determine the maximum curvature of a cubic Bézier-spline curve that interpolates an ordered set of points in ${\mathbb{R}}^2$ or ${\mathbb{R}}^3$. The procedures are designed from closed-form formulas for such open and closed curves.

1. Introduction

We recall here the definition of curvature of a smooth curve. This is a fundamental concept in differential geometry that has been studied deeply in applied mathematics, engineering, and computer graphics. The problem of finding the curvature extremum has been investigated by many authors. We can point out some of their results on establishing necessary and sufficient conditions for the regularity of offset curves or tubular surfaces, designing various types of aesthetic curves from constrained conditions, representing and modifying curves to adapt principles of interpolation and animation, etc. In the following, we give a short description of the method for solving this problem and its limitations.

If $\mathbf{r}:[a,b]⟶{\mathbb{R}}^3$ is the position vector of a smooth curve $\mathcal{L}$, then the point $\mathbf{r}\left(t\right)$ of $\mathcal{L}$ at $t\in [a,b]$ is written as a 3-tuple $\left(x\right(t),y(t),z(t\left)\right)$, or a column vector $\langle x\left(t\right),y\left(t\right),z\left(t\right)\rangle ={(x\left(t\right),y\left(t\right),z\left(t\right))}^T$ when used in a matrix expression; hence $\mathbf{r}{\left(t\right)}^T$ is the row matrix $\left[\begin{array}{ccc}x\left(t\right)& y\left(t\right)& z\left(t\right)\end{array}\right]$. The image of $\mathbf{r}$ or the map $\mathbf{r}$ itself is called a parametrization of $\mathcal{L}$. Along $\mathcal{L}$, we consider the coordinates of the Frenet frame $\{\widehat{\mathbf{T}},\widehat{\mathbf{N}},\widehat{\mathbf{B}}\}$ whose origin is located at points of $\mathcal{L}$. From differential geometry and calculus (see [1,2]), we know that

$$\widehat{\mathbf{T}}=\frac{\mathbf{v}}{\parallel \mathbf{v}\parallel },\widehat{\mathbf{B}}=\frac{\mathbf{v}\times \mathbf{a}}{\parallel \mathbf{v}\times \mathbf{a}\parallel },\widehat{\mathbf{N}}=\widehat{\mathbf{B}}\times \widehat{\mathbf{T}},$$

where $\mathbf{v}\left(t\right)={\mathbf{r}}^{\prime }\left(t\right)={(x^{\prime }\left(t\right),y^{\prime }\left(t\right),z^{\prime }\left(t\right))}^T$ and $\mathbf{a}\left(t\right)={\mathbf{r}}^{″}\left(t\right)={(x^{″}\left(t\right),y^{″}\left(t\right),z^{″}\left(t\right))}^T$ are usually called the velocity vector and the acceleration vector of $\mathcal{L}$ at t (or at $\mathbf{r}\left(t\right)$).

Let $\mathbf{r}=\mathbf{r}\left(t\right)={(x\left(t\right),y\left(t\right),z\left(t\right))}^T$ be a parametrization of a smooth space curve $\mathcal{L}$ with $\mathbf{a}={\mathbf{r}}^{″}\in C\left([a,b]\right)$. Then, the curvature of $\mathcal{L}$ at t is given by

$$\kappa =\frac{\parallel \mathbf{v}\times \mathbf{a}\parallel }{v^3},$$

(1)

where $\mathbf{v}=\mathbf{v}\left(t\right)={\mathbf{r}}^{\prime }\left(t\right)$, $\mathbf{a}=\mathbf{a}\left(t\right)={\mathbf{r}}^{″}\left(t\right)$, and $v=\parallel \mathbf{v}\left(t\right)\parallel$. Under this assumption, we usually find the maximum value of $\kappa \left(t\right)$ on $[a,b]$ by choosing the largest value of $\kappa$ at the points where ${\kappa }^{\prime }\left(t\right)=0$. However, how to solve the equation ${\kappa }^{\prime }\left(t\right)=0$ exactly or approximately? In general, we cannot do it. Moreover, apart from this, the expression for ${\kappa }^{\prime }\left(t\right)$ is complicated! Let us take $K={\kappa }^2$ and write it out in the components of $\mathbf{v}\times \mathbf{a}$. This gives

$$K=\frac{{(y^{\prime }{z}^{″}-z^{\prime }{y}^{″})}^2+{(z^{\prime }{x}^{″}-x^{\prime }{z}^{″})}^2+{(x^{\prime }{y}^{″}-y^{\prime }{x}^{″})}^2}{{[{\left(x^{\prime }\right)}^2+{\left(y^{\prime }\right)}^2+{\left(z^{\prime }\right)}^2]}^3}.$$

On the other hand, if $\mathbf{r}=\mathbf{r}\left(s\right)$ is the parametrization of $\mathcal{L}$ by the arc length s, then from the relation ${\widehat{\mathbf{T}}}^{\prime }\left(s\right)=\kappa \left(s\right)\widehat{\mathbf{N}}\left(s\right)$ we obtain

$${\widehat{\mathbf{T}}}^{″}\left(s\right)={\kappa }^{\prime }\left(s\right)\widehat{\mathbf{N}}\left(s\right)+\kappa \left(s\right){\widehat{\mathbf{N}}}^{\prime }\left(s\right).$$

This gives a simple expression for ${\kappa }^{\prime }\left(s\right)$:

$${\widehat{\mathbf{T}}}^{″}\left(s\right)·\widehat{\mathbf{N}}\left(s\right)={\kappa }^{\prime }\left(s\right).$$

However, we again encounter another hard problem: how to convert a parametrization of a smooth space curve $\mathcal{L}$ by a general parameter t into the parametrization by the arc length s. In general, this is impossible.

To overcome those obstacles, many researchers restricted their attention to the class of Bézier curves and their variants. Recently, the papers related to maximizing or minimizing the curvature of these curves provided a lot of theoretical results and useful algorithms, as well as practical tools. We highlight some representative papers with a brief note. Ref. [3] presents a unique design on a piecewise quadratic Bézier curve that interpolates its local maximum curvature points that are also its control points. Ref. [4], adopted from [3], proposes new methods to modify local curvature at the interpolation points by taking basis functions of higher degree. Ref. [5] provides conditions for the curvature of a quadratic rational Bézier curve to be monotone or to have a local minimum and maximum. Ref. [6] establishes conditions for Bézier plane curves generated by a matrix to have monotone curvature. Ref. [7] also establishes conditions for Bézier curves to have monotone curvature, based on control points of the position vector of the curve and its derivatives. Ref. [8] treats typical Bézier plane curves with one curvature extremum that can be easily calculated, which can help to divide the curve into two typical curves with monotone curvature.

Traditionally, the papers noted above and many others paid much attention to control points and polygons. Actually, these objects have direct effects on the shape of the curves, so they have been modified in order to obtain a curve with properties needed in design applications. However, we have some changes in mind when relating this widespread trend to our result in [9]. The formula in [9] (Theorem 3.1) can be seen as a way of approximation by interpolation with Bézier-spline curves. Therefore, we prefer to place emphasis on interpolated points. These points can be chosen in a way to design curves in ${\mathbb{R}}^2$ or ${\mathbb{R}}^3$ with desired shapes or can be taken from special partitions of the parameter interval of a smooth curve with given parametrization to approximate its curvature extremum.

Now, we go back to our main purpose: making a Maple procedure to compute the maximum value ${\kappa }_{max}$ of the curvature. We restrict our attention to Bézier-spline curves. This objective is based on the power of Maple on symbolic computation and on solving polynomial equations of high degree and on the explicit piecewise cubic parametrization of these special curves.

The present paper is organized as follows. In Section 2, we construct the Maple procedures to represent Bézier-spline space curves for both open and closed curves. In Section 3, we propose a pseudo-algorithm for computing ${\kappa }_{max}$, then we provide the full code of the procedures corresponding to the algorithm. In Section 4, we discuss some modifications to obtain procedures to represent Bézier-spline plane curves and to compute their maximum curvature. In Section 5, we give some concluding remarks.

2. Bézier-Spline Space Curves with Maple Parametrization

In this section, we consider a Bézier-spline curve, which is obtained from a closed-form solution to the inverse problem, which interpolates an ordered set of points ${\mathbf{s}}_0$, ${\mathbf{s}}_1$, …, ${\mathbf{s}}_n$, given in [9]. Such a plane curve can be obviously extended to a space curve $\mathcal{C}$ given by a piecewise cubic function $\mathbf{f}$ of $C^2\left([0,n]\right)$. According to the construction of such curves in [9], we present here a more convenient way to derive their parametrization.

First, $\mathbf{f}\left(t\right)$ is composed of the cubic functions ${\mathbf{f}}_k$ given by

$${\mathbf{f}}_k\left(t\right)={(1-t)}^3{\mathbf{s}}_{k-1}+3t{(1-t)}^2{\mathbf{p}}_{k-1}+3t^2(1-t){\mathbf{q}}_k+t^3{\mathbf{s}}_k,$$

(2)

where $k=1,\dots ,n$ and $t\in [0,1]$, and ${\mathbf{f}}_k$ is the parametrization of the cubic Bézier curve ${\mathcal{C}}_k$ with the control points ${\mathbf{s}}_{k-1}$, ${\mathbf{p}}_{k-1}$, ${\mathbf{q}}_k$, and ${\mathbf{s}}_k$. These points satisfy the known relations

$${\mathbf{p}}_{k-1}=\frac{2}{3}{\mathbf{b}}_{k-1}+\frac{1}{3}{\mathbf{b}}_k,{\mathbf{q}}_k=\frac{1}{3}{\mathbf{b}}_{k-1}+\frac{2}{3}{\mathbf{b}}_k,$$

(3)

and

$${\mathbf{s}}_k=\frac{{\mathbf{q}}_k+{\mathbf{p}}_k}{2}.$$

(4)

On the other hand, from [9] (Theorem 3.1), the points ${\mathbf{b}}_k,{\mathbf{s}}_k$, $k=0,\dots ,n$, are now in ${\mathbb{R}}^3$ with ${\mathbf{b}}_0={\mathbf{s}}_0$ and ${\mathbf{b}}_n={\mathbf{s}}_n$, and for $k=1,\dots ,n-1$, we have

$${\mathbf{b}}_k=\frac{{\beta }_{n-1-k}}{{\beta }_{n-1}}\left[{(-1)}^k{\mathbf{s}}_0+6\sum _{j=1}^{k-1}{(-1)}^{k-j}{\beta }_{j-1}{\mathbf{s}}_j\right]+\frac{{\beta }_{k-1}}{{\beta }_{n-1}}\left[6\sum _{j=k}^{n-1}{(-1)}^{j-k}{\beta }_{n-1-j}{\mathbf{s}}_j+{(-1)}^{n-k}{\mathbf{s}}_n\right],$$

(5)

where ${\beta }_k$ ($k=0,\dots ,n-1$) is evaluated by the formula

$${\beta }_k=2^k\sum _{m=0}^{\lfloor k/2\rfloor }\left(\begin{array}{c}k+1\\ k-2m\end{array}\right){(3/4)}^m.$$

(6)

Finally, we can give a simple process to obtain a so-called relaxed, uniform B-spline space curve $\mathcal{C}$ that interpolates an ordered set of points ${\mathbf{s}}_0$, ${\mathbf{s}}_1$, …, ${\mathbf{s}}_n$ (see [10]). The parametrization $\mathbf{f}$ of $\mathcal{C}$ is a piecewise cubic function on $[0,n]$ whose components ${\mathbf{f}}_k$, $k=1,\dots ,n$, derived from (2) and (3), can be now given by

$$\begin{array}{cc}\hfill {\mathbf{f}}_k\left(t\right)& ={(k-t)}^3{\mathbf{s}}_{k-1}+(t-k+1){(k-t)}^2(2{\mathbf{b}}_{k-1}+{\mathbf{b}}_k)\hfill \\ \hfill & +{(t-k+1)}^2(k-t)({\mathbf{b}}_{k-1}+2{\mathbf{b}}_k)+{(t-k+1)}^3{\mathbf{s}}_k,\hfill \end{array}$$

(7)

where $t\in [k-1,k]$, and the ${\mathbf{b}}_k$ are obtained from (5) for $k=1,\dots ,n-1$, and ${\mathbf{b}}_0={\mathbf{s}}_0$ and ${\mathbf{b}}_n={\mathbf{s}}_n$. Since ${\mathbf{f}}_1^{″}\left(0\right)={\mathbf{f}}_n^{″}\left(n\right)=\mathbf{0}$, the curvatures of $\mathcal{C}$ at $t=0$ and $t=n$ are both zero and we call $\mathcal{C}$ a relaxed Bézier-spline space curve.

We are interested in the implementation of the above parametrization by a Maple procedure. We list here some Maple commands that will appear in our procedures. They are all very important and frequently used in graphic and computation programming: args, nargs, op, nops, ERROR, RETURN, convert, evalf, diff, expand, floor, for, fsolve, map, max, min, piecewise, plot, plot3d, proc, seq, solve, unapply, and while; in addition, LinearAlgebra and plots are the great packages containing many procedures for specific purposes. A declaration to create a function, e.g., f, such as f:=x->F(x) or f:=unapply(F(x),x) (sub-procedures in F(x) are evaluated first), where F(x) is an expression or a list of expressions in x, is a very useful and convenient tool. In addition, the conditional structures if-then-else and if-then-elif-else are indispensable in branch programming, whereas the type “list” is a flexible ordered arrangement of operands (or components, elements) inside the square brackets [, ]. See [11,12] and Maple help pages in each session to know more details about meaning, syntax, and usage of these commands, structures, and types. The implementation of some specific task by calling a procedure name together with appropriate arguments is usually said to be a calling sequence. As a convention, we choose type list for elements of ${\mathbb{R}}^2$ or ${\mathbb{R}}^3$.

Now, let us make a procedure to compute $\mathbf{f}\left(t\right)$ on $[0,n]$ from (5)–(7), with ${\mathbf{b}}_0={\mathbf{s}}_0$ and ${\mathbf{b}}_n={\mathbf{s}}_n$. It takes $S=\{{\mathbf{s}}_0,{\mathbf{s}}_1,\dots ,{\mathbf{s}}_n\}$ as its input and gives $\mathbf{f}\left(t\right)$ as its output in the form of $\left[F\right(t),G(t),H(t\left)\right]$ such that $F,G,H$ are the piecewise functions in t on $[0,n]$. This procedure is called BScurve3d and its full code is given in the following.

BScurve3d

BScurve3d:=proc(Lst::list(list(realcons)))

local n,S,b,B,f,F,G,H,j,k;

n:=nops(Lst)-1:

for k from 0 to n do

S[k]:=Lst[k+1]:

end do:

for j from 0 to n-1 do

B[j]:=2^j*add(binomial(j+1,j-2*m)*(3/4)^m,m=0..floor(j/2)):

end do:

for k from 1 to n-1 do

b[k]:=(B[n-1-k]/B[n-1])*((-1)^k*S[0]+6*add((-1)^(k-j)*B[j-1]*S[j],j=1..k-1))

+(B[k-1]/B[n-1])*((-1)^(n-k)*S[n]+6*add((-1)^(j-k)*B[n-1-j]*S[j],j=k..n-1)):

end do:

b[0]:=S[0]:

b[n]:=S[n]:

for k from 1 to n do

f[k]:=unapply(expand((k-t)^3*S[k-1]+(t-k+1)*(k-t)^2*(2*b[k-1]+b[k])

+(t-k+1)^2*(k-t)*(b[k-1]+2*b[k])+(t-k+1)^3*S[k]),t):

end do:

F:=t->piecewise(-1<=t and t<1,f[1](t)[1],

seq([(k-1)<=t and t<k,f[k](t)[1]][],k=2..n-1),(n-1)<=t and t<n+1,f[n](t)[1]):

G:=t->piecewise(-1<=t and t<1,f[1](t)[2],

seq([(k-1)<=t and t<k,f[k](t)[2]][],k=2..n-1),(n-1)<=t and t<n+1,f[n](t)[2]):

H:=t->piecewise(-1<=t and t<1,f[1](t)[3],

seq([(k-1)<=t and t<k,f[k](t)[3]][],k=2..n-1),(n-1)<=t and t<n+1,f[n](t)[3]):

RETURN(unapply([F(t),G(t),H(t)],t));

end proc:

To declare a finite sequence of indexed expressions, we can use the operator [ ] to extract the contents of a list. For example, the command [$1,2,3$][ ] results in $1,2,3$, and the declaration of the sequence

$$1\le t\mathrm{and}t<2,f_2\left(t\right),2\le t\mathrm{and}t<3,f_3\left(t\right),\dots ,n-2\le t \mathrm{and}t<n-1,f_{n-1}\left(t\right)$$

can be written as: seq([k-1<=t and t<k,f[k](t)][],k=2..n-1). We have used this declaration in BScurve3d.

If we have an ordered set of $(n+1)$ distinct points in ${\mathbb{R}}^3$ that are declared in Maple as a list of lists

$$S:=\left[\right[a_0,b_0,c_0],[a_1,b_1,c_1],\dots ,[a_n,b_n,c_n\left]\right],$$

then we can obtain the position vector of the Bézier−spline curve $\mathcal{C}$ that interpolates these points by calling f:=BScurve3d(S).The plot of $\mathcal{C}$ can be made by the procedure spacecurve in the plots package with the command

$$\mathtt{plots}\left[\mathtt{spacecurve}\right]\left(\mathtt{f}\right(\mathtt{t}),\mathtt{t}=\mathtt{0}\dots \mathtt{n},\mathtt{opts});$$

The ‘opts’ means plotting options. To implement these steps, for instance, we set

$$L:=\left[\right[-2,5,1],[0,-1,0],[-3,1,2],[1,2,3],[-5,-2,1],[2,-4,-5],[0,1,7],[-6,3,2\left]\right]$$

and $\mathbf{f}:=\mathtt{BScurve}\mathtt{3}\mathtt{d}\left(L\right)$. Then, the plot of the Bézier-spline curve $\mathcal{C}$ that interpolates L is given by the calling sequence

$$\mathtt{plots}\left[\mathtt{spacecurve}\right]\left(f\right(t),t = 0..7,\mathtt{opts});$$

Figure 1 shows the curve $\mathcal{C}$ and the points of L. We take one more example of interpolating an ordered set of points on a given space curve to see how well a Bézier-spline curve fits this curve. Let ${\mathcal{L}}_1$ be a curve whose parametrization is $\mathbf{r}\left(t\right)={(tcost,0.8tsint,tcost-sint)}^T$, $t\in [-2.4,5]$, with the declaration

$$\mathbf{r}:=\mathtt{t}->[\mathtt{t}*\mathtt{cos}(\mathtt{t}),(\mathtt{0}.\mathtt{8})*\mathtt{t}*\mathtt{sin}(\mathtt{t}),\mathtt{t}*\mathtt{cos}(\mathtt{t})-\mathtt{sin}(\mathtt{t}\left)\right]:$$

Consider a partition of $[-2.4,5]$ by the points $t_0:=-2.4$, $t_1:=-1.2$, $t_2:=0.0$, $t_3:=1.4$, $t_4:=2.7$, $t_5:=3.8$, $t_6:=5.0$, and set

$$L_1:=[r\left(t0\right),r\left(t1\right),\cdots ,r\left(t6\right)]:$$

The curve ${\mathcal{L}}_1$ and its approximation by a Bézier-spline curve ${\mathcal{C}}_1$ are also given in Figure 1.

We will make some changes to obtain a so-called closed, uniform B-spline space curve $\mathcal{C}$ that interpolates an ordered set of points ${\mathbf{s}}_0$, ${\mathbf{s}}_1$, …, ${\mathbf{s}}_n$, with ${\mathbf{s}}_n={\mathbf{s}}_0$ (see [10]). We call such a curve a closed Bézier-spline space curve. In this case, we choose appropriate settings to have again the relations (3) and (4) at the common point. The first setting should be ${\mathbf{b}}_n={\mathbf{b}}_0$. Then, $\mathcal{C}$ is still composed of the cubic Bézier curves ${\mathcal{C}}_k$, $k=1,\dots ,n$, as above. Specifically, at the interpolated point ${\mathbf{s}}_n={\mathbf{s}}_0$, (4) becomes

$${\mathbf{s}}_0=\frac{{\mathbf{p}}_0+{\mathbf{q}}_n}{2},$$

(8)

where we have from (3) that

$${\mathbf{p}}_0=\frac{2}{3}{\mathbf{b}}_0+\frac{1}{3}{\mathbf{b}}_1,{\mathbf{q}}_n=\frac{1}{3}{\mathbf{b}}_{n-1}+\frac{2}{3}{\mathbf{b}}_n=\frac{1}{3}{\mathbf{b}}_{n-1}+\frac{2}{3}{\mathbf{b}}_0.$$

(9)

Now, from (8) and (9), we get the last setting

$${\mathbf{b}}_0=\frac{1}{4}(6{\mathbf{s}}_0-{\mathbf{b}}_{n-1}-{\mathbf{b}}_1).$$

(10)

It is easy to have another Maple procedure, say BScurve3dC, for representing a closed Bézier-spline space curve $\mathcal{C}$ from the dataset $\{{\mathbf{s}}_0,{\mathbf{s}}_1,\dots ,{\mathbf{s}}_n\}$ with ${\mathbf{s}}_n={\mathbf{s}}_0$, and the settings (10) and ${\mathbf{b}}_n={\mathbf{b}}_0$. From the above discussion, the parametrization $\mathbf{f}$ of $\mathcal{C}$ has the components ${\mathbf{f}}_k$ given by (7), and we can check that

$${\mathbf{f}}_1^{\prime }\left(0\right)={\mathbf{f}}_n^{\prime }\left(n\right)=-3{\mathbf{s}}_0+2{\mathbf{b}}_0+{\mathbf{b}}_1,{\mathbf{f}}_1^{″}\left(0\right)={\mathbf{f}}_n^{″}\left(n\right)=6({\mathbf{s}}_0-{\mathbf{b}}_0).$$

Therefore, $\mathbf{f}$ is in $C^2\left([0,n]\right)$ again and ${\mathcal{C}}_1,{\mathcal{C}}_n$ have the same curvature at their common point.

Let us take some examples on using BScurve3dC. As the steps to display closed Bézier-spline space curves are the same as for BScurve3d, we just give the graphical results of these examples. Note that the initial and terminal points of the input list for BScurve3dC have to be the same. Let L be the list

$$L:=\left[\right[0,3,-1],[1,-2,2],[-2,0,-5],[2,4,3],[-1,2,4],[6,-1,0],[0,3,-1\left]\right].$$

The display of the closed Bézier-spline curve ${\mathcal{C}}_L$ that interpolates L is given in Figure 2. Let $\mathbf{r}$ be the parametrization of a closed space curve $\mathcal{L}$ called a trefoil knot (see [1] (p. 897)) with

$$\mathbf{r}\left(t\right):=\left[\right(1+0.3cos\left(3t\right))cos(2t),(1+0.3cos\left(3t\right))sin(2t),0.35sin(3t\left)\right](t\in [0,2\pi \left]\right).$$

The endpoints of $\mathcal{L}$ coincide with the point $[1.3,0,0]$. Let T be a set of points on $\mathcal{L}$ such that

$$T:=\left[\right[1.3,0,0],\mathbf{r}(0.5),\mathbf{r}(1.2),\mathbf{r}(1.8),\mathbf{r}(2.5),\mathbf{r}(3.2),\mathbf{r}(3.9),\mathbf{r}(4.5),\mathbf{r}(5.1),\mathbf{r}(5.8),[1.3,0,0\left]\right].$$

In Figure 2, we also give the display of $\mathcal{L}$ together with its approximation ${\mathcal{C}}_T$ that interpolates T.

3. Computation of the Maximum Curvature of a Bézier-Spline Curve

Let $\mathcal{C}$ be a Bézier-spline curve that interpolates an ordered set T of points ${\mathbf{s}}_0$, ${\mathbf{s}}_1$, …, ${\mathbf{s}}_n$ in ${\mathbb{R}}^3$. Then $\mathbf{r}\left(t\right)$, the position vector of $\mathcal{C}$, is a piecewise cubic function of $C^2\left([0,n]\right)$, given by its components ${\mathbf{f}}_k\left(t\right)$ in (7).

Avoiding the square root function, we have from (1):

$$K:={\kappa }^2=\frac{{\parallel \mathbf{v}\times \mathbf{a}\parallel }^2}{v^6}=\frac{(\mathbf{v}\times \mathbf{a})·(\mathbf{v}\times \mathbf{a})}{{(\mathbf{v}·\mathbf{v})}^3}.$$

Then, we have that   

$$\begin{array}{cc}\hfill K^{\prime }& =\frac{{[(\mathbf{v}\times \mathbf{a})·(\mathbf{v}\times \mathbf{a})]}^{\prime }{(\mathbf{v}·\mathbf{v})}^3-[(\mathbf{v}\times \mathbf{a})·(\mathbf{v}\times \mathbf{a})]3{(\mathbf{v}·\mathbf{v})}^{2}2(\mathbf{v}·{\mathbf{v}}^{\prime })}{{(\mathbf{v}·\mathbf{v})}^6}\hfill \\ & =\frac{2[(\mathbf{v}\times \mathbf{a})·(\mathbf{v}\times {\mathbf{a}}^{\prime })](\mathbf{v}·\mathbf{v})-6[(\mathbf{v}\times \mathbf{a})·(\mathbf{v}\times \mathbf{a})](\mathbf{v}·\mathbf{a})}{{(\mathbf{v}·\mathbf{v})}^4}.\hfill \end{array}$$

As $\kappa$ attains its maximum value on $[0,n]$ only at solutions of the equation $K^{\prime }=0$ or

$$(\mathbf{v}·\mathbf{v})[(\mathbf{v}\times \mathbf{a})·(\mathbf{v}\times {\mathbf{a}}^{\prime })]-3[(\mathbf{v}\times \mathbf{a})·(\mathbf{v}\times \mathbf{a})](\mathbf{v}·\mathbf{a})=0$$

(11)

in the intervals $(i-1,i)$ and at their endpoints $i-1,i$, with $i=1,\dots ,n$, we can find ${\kappa }_{\max }$ on $[0,n]$ by the procedure MaxCurvature3d. However, at first, we present the procedure in the form of a pseudo-code algorithm. It would be easy to translate statements in such algorithms into Maple codes or other programming languages. Moreover, our discussion on how to use appropriate commands for a specific purpose will give a clear description of our procedures. In addition, we sometimes use built-in Maple procedures in those algorithms with their most simple form for convenience and simplicity.

Note that the left-hand side of (11) is a polynomial of degree at most 7. Letting Q be such a polynomial (in one variable, say, t), we will use a powerful tool of Maple to find numerically all the zeros of Q in a given interval. That is the procedure fsolve and it has been called in Algorithm 1 by the command: fsolve(Q,t,$\alpha$…$\beta$). The output of this calling is a sequence of all real zeros of Q in $[\alpha ,\beta ]$. Moreover, the expressions of dot and cross products are given by the great package of Maple: LinearAlgebra. We also select points in $[0,n]$ at which $\kappa$ attains its maximum value. Thus, the output of Algorithm 1 consists of ${\kappa }_{max}$ and the set ${\kappa }^{-1}\left({\kappa }_{max}\right)$ in $[0,n]$. Now, for relaxed Bézier-spline curves, we give the full code of MaxCurvature3d at the end of this section.

Algorithm 1: Finding the maximum curvature of a Bézier-spline curve

Input:  a set T of $(n+1)$ points in ${\mathbb{R}}^3$; Output:  The maximum curvature ${\kappa }_{max}$ of the Bézier-spline curve interpolating T; 1: $S_p:=\left\{\right\}$; 2: for$i=1$ndo 3:    $\mathbf{v}:={\mathbf{f}}_i^{\prime }\left(t\right)$, $\mathbf{a}:={\mathbf{f}}_i^{″}\left(t\right)$     // ${\mathbf{f}}_i$ from BScurve3d; 4:    $Q:=\mathrm{the}\mathrm{left}-\mathrm{hand}\mathrm{side}\mathrm{of}\left(\right)$; 5:    $P:=\left\{\mathtt{fsolve}\right($ Q,t,i-1 .. i$\left)\right\}$ // Solving (11) for $t\in [i-1,i]$; 6:    $S_p:=S_p\cup P$; 7:    $S:=\{F_i\left(t\right):t\in P\}\cup \{F_i(i-1),F_i\left(i\right)\}$ // $F_i:=t\mapsto \kappa \left(t\right)=\parallel \mathbf{v}\times \mathbf{a}\parallel /v^3$, $v:=\parallel \mathbf{v}\parallel$; 8:    $m_i:=maxS$; 9: end for 10: ${\kappa }_{max}:=max\{m_1,\dots ,m_n\}$; 11: $M:=S_p\cup \{0,1,\dots ,n\}$; $T_p:=F^{-1}\left({\kappa }_{max}\right)\subset M$ // $F:=t\mapsto \kappa \left(t\right)$ on $[0,n]$, $F=F_i$ on $[i-1,i]$; 12: return ${\kappa }_{max}$ and $T_p$;

MaxCurvature3d

MaxCurvature3d:=proc(L::list(list(realcons)))

local a,ad,A,AD,m,n,S,b,B,f,H,E,F,G,R,k,i,j,v,V,M,Kmax,N,P,Q,Sp,Tp,Tpoint;

n:=nops(L)-1:

for k from 0 to n do

S[k]:=L[k+1]:

end do:

for j from 0 to n-1 do

B[j]:=2^j*add(binomial(j+1,j-2*m)*(3/4)^m,m=0..floor(j/2)):

end do:

for k from 1 to n-1 do

b[k]:=(B[n-1-k]/B[n-1])*((-1)^k*S[0]+6*add((-1)^(k-j)*B[j-1]*S[j],j=1..k-1))

+(B[k-1]/B[n-1])*((-1)^(n-k)*S[n]+6*add((-1)^(j-k)*B[n-1-j]*S[j],j=k..n-1)):

end do:

b[0]:=S[0]:

b[n]:=S[n]:

for k from 1 to n do

f[k]:=unapply(expand((k-t)^3*S[k-1]+(t-k+1)*(k-t)^2*(2*b[k-1]+b[k])

+(t-k+1)^2*(k-t)*(b[k-1]+2*b[k])+(t-k+1)^3*S[k]),t):

end do:

Sp:={}:

for i from 1 to n do

v:=unapply(diff(f[i](t),t),t):

a:=unapply(diff(f[i](t),t$2),t):

ad:=unapply(diff(f[i](t),t$3),t):

V:=convert(v(t),Vector):

A:=convert(a(t),Vector):

AD:=convert(ad(t),Vector):

M:=expand(LinearAlgebra[DotProduct](V,V,conjugate=false)):

N:=expand(LinearAlgebra[DotProduct](V,A,conjugate=false)):

G:=map(expand,LinearAlgebra[CrossProduct](V,A)):

H:=LinearAlgebra[CrossProduct](V,AD):

E:=expand(LinearAlgebra[DotProduct](G,H,conjugate=false)):

R:=expand(LinearAlgebra[DotProduct](G,G,conjugate=false)):

Q:=expand(M*E-3*R*N):

F[i]:=unapply(sqrt(abs(R))/abs(M)^(3/2),t):

P:={fsolve(Q,t,i-1..i)}:

m[i]:=max(seq(F[i](P[j]),j=1..nops(P)),F[i](i-1),F[i](i)):

Sp:=Sp union P:

end do:

Kmax:=max(seq(m[j],j=1..n)):

Tpoint:={}:

for i from 1 to n do

if (abs(F[i](i-1)-Kmax)=0) or (abs(F[i](i-1)-Kmax)=0.) then

Tpoint:= Tpoint union {i-1}:

elif (abs(F[i](i)-Kmax)=0) or (abs(F[i](i)-Kmax)=0.) then

Tpoint:= Tpoint union {i}:

end if:

end do:

if nops(Sp)=0 then

RETURN(Kmax,Tpoint);

end if:

for j from 1 to nops(Sp) do

if (floor(Sp[j])=n) then

if (abs(F[n](Sp[j])-Kmax)=0 or abs(F[n](Sp[j])-Kmax)=0.) then

Tpoint:= Tpoint union {n}:

end if:

elif (abs(F[floor(Sp[j])+1](Sp[j])-Kmax)=0 or

abs(F[floor(Sp[j])+1](Sp[j])-Kmax)=0.) then

Tpoint:= Tpoint union {Sp[j]}:

end if:

end do:

RETURN(Kmax,Tpoint);

end proc:

Here we give an explanation of how to determine the set ${\kappa }^{-1}\left({\kappa }_{max}\right)=\{t\in [0,n]:\kappa \left(t\right)=$${\kappa }_{max}\}$. We first check whether i ($\in \{0,1,\dots ,n\}$) belongs to this set, then we check the same for all solutions of the equation $Q=0$. In addition, when we need the value $\mathbf{f}\left(t\right)$ for a point $t\in [0,n]$, we should take an appropriate component ${\mathbf{f}}_k$ of $\mathbf{f}$; if $t=n$, then we take $k=n$, else we take $k=\lfloor t\rfloor +1$, since $\lfloor t\rfloor \le t<\lfloor t\rfloor +1$.

We use MaxCurvature3d to determine ${\kappa }_{max}$ in the examples whose graphical results are given in Figure 1. From the lists L and $L_1$ in these examples, the calling sequences MaxCurvature3d(L) and MaxCurvature3d($L_1$) result in

$$11.87724489,\left\{2.913967861\right\}\mathrm{and}3.135008280,\left\{2\right\},$$

respectively. If the parametrization of ${\mathcal{C}}_1$ is $\mathbf{h}$, then

$$\begin{array}{cc}\hfill \mathbf{f}\left(2.913967861\right)& =[1.06704678,2.13642213,2.94912425],\hfill \\ \hfill \mathbf{h}\left(2\right)& =[-0.8322936730,1.454875883,-1.741591100]\hfill \end{array}$$

are the maximum curvature points on the curves $\mathcal{C}$ and ${\mathcal{C}}_1$, respectively.

The new version of MaxCurvature3d for closed Bézier-spline curves, say, MaxCurvature3dC, can be derived easily from Algorithm 1 with the modification “${\mathbf{f}}_i$ from BScurve3dC”. Then, we use MaxCurvature3dC to find ${\kappa }_{max}$ of ${\mathcal{C}}_L$ displayed in Figure 2 and we obtain the result

$$1.652746390,\left\{1.979287163\right\}.$$

Let $\mathbf{g}$ be the parametrization of ${\mathcal{C}}_L$. The point

$$\mathbf{g}\left(1.979287163\right)=[-2.02583824,-0.094336298,-4.99371046]$$

is given in Figure 3 and this maximum point of curvature of ${\mathcal{C}}_L$ is very close to the interpolated point $[-2,0,-5]$.

To avoid a comparison error between fractions and decimal numbers, we may use the decimal point ‘.’ for at least one component of the points in the list argument of MaxCurvature3d. Note that the result of fsolve only contains decimal numbers, so it will give us ‘$2.$’, for instance, if it contains the integer ‘2’. To cover this case, we add a condition such as ‘$A-B=0.$’ in the definition of MaxCurvature3d, and it should be sometimes ‘$|A-B|<{10}^{-m}$’ with some positive integer m when we need to obtain an expected result.

4. Remarks on the Two-Dimensional Case

For a relaxed Bézier-spline plane curve $\mathcal{L}$ that interpolates a dataset M of points ${\mathbf{s}}_0$, ${\mathbf{s}}_1$, …, ${\mathbf{s}}_n$ in ${\mathbb{R}}^2$, we have already a procedure to obtain its position vector $\mathbf{f}\left(t\right)$. That is just removing the lines H:=t->⋯ and modifying the lines RETURN(t->[F(t),G(t),H(t)]) to RETURN(t->[F(t),G(t)]) in BScurve3d, and the remaining part is for BScurve2d, the Maple parametrization of $\mathcal{L}$. Maple provides the plot procedure to display plane curves with their parametrization $\left[x\right(t),y(t\left)\right]$, $t\in [a,b]$, by the declaration

$$\mathtt{plot}\left(\right[\left(t\right),y\left(t\right),t = a\dots b],\mathrm{opts});$$

Accordingly, we first set $\mathbf{f}:=\mathtt{BScurve}\mathtt{2}\mathtt{d}$(M), then we call

$$\mathtt{plot}\left(\right[\mathtt{op}\left(f\right(t\left)\right),t = 0\dots n],\mathrm{opts});\mathrm{or}\mathtt{plot}\left(\right[f\left(t\right)\left[\right],t = 0\dots n],\mathrm{opts});$$

to display $\mathcal{L}$ . Similarly, we get BScurve2dC, the new version of BScurve2d for closed Bézier−spline plane curves, from BScurve3dC.

MaxCurvature3d can be modified to use only sets of points in ${\mathbb{R}}^2$ and we call its new version MaxCurvature2d. Let $\mathcal{L}$ be a smooth plane curve with parametrization $\mathbf{r}$. From the formula $\kappa \left(s\right)={\widehat{\mathbf{T}}}^{\prime }\left(s\right)·\widehat{\mathbf{N}}\left(s\right)$ with arc length parameter s, we can write

$$\kappa =\frac{{\widehat{\mathbf{T}}}^{\prime }}{v}·\widehat{\mathbf{N}}=\frac{1}{v}\left[{\left(\frac{1}{v}\right)}^{\prime }\mathbf{v}+\frac{1}{v}\mathbf{a}\right]·\widehat{\mathbf{N}}=\frac{1}{v^2}\mathbf{a}·\left(J\widehat{\mathbf{T}}\right)=\frac{\mathbf{a}·\mathbf{u}}{v^3}$$

for a general parameter t, where

$$\mathbf{v}={\mathbf{r}}^{\prime },v=\parallel \mathbf{v}\parallel ,\mathbf{u}=J\mathbf{v},\mathbf{a}={\mathbf{v}}^{\prime },J=\left[\begin{array}{cc}0& \hfill -1\\ 1& \hfill 0\end{array}\right].$$

We also consider $K:={\kappa }^2={(\mathbf{a}·\mathbf{u})}^2/{(\mathbf{v}·\mathbf{v})}^3$ and derive the following equation from $K^{\prime }=0$:

$$(\mathbf{v}·\mathbf{v})(\mathbf{a}·\mathbf{u})({\mathbf{a}}^{\prime }·\mathbf{u})-3{(\mathbf{a}·\mathbf{u})}^2(\mathbf{v}·\mathbf{a})=0.$$

(12)

This equation has the same role as (11), so we can make a new version of MaxCurvature3d, say, MaxCurvature2d, following the steps in Algorithm 1 with some modification: there is no cross product in this version. The expression of the local variable Q in the definition of MaxCurvature2d is given by the left-hand side of (12). Note that the function ${\mathbf{F}}_i$ in the pseudo-code algorithm of MaxCurvature2d is now

$${\mathbf{F}}_i:=\frac{|\mathbf{a}·\mathbf{u}|}{v^3}.$$

In the definition of MaxCurvature2d, the matrix J is declared at right above Sp:={} by

$$\mathtt{J}:=\mathtt{Matrix}(\mathtt{2},\mathtt{2},[\mathtt{0},-\mathtt{1},\mathtt{1},\mathtt{0}\left]\right):$$

and, then, we set U:=J.V inside the for loop at right below.

We give two more examples on getting the maximum curvature of a relaxed Bézier-spline plane curve and its maximum curvature points. Let M be the list

$$M:=\left[\right[0.5,1],[0,2],[0.8,3],[1.7,2.25],[2.5,1.8],[3.5,2.2],[3.2,2.8\left]\right]$$

and $\mathcal{C}$ be the Bézier-spline curve that interpolates M. Then, we obtain the parametrization of $\mathcal{C}$ by setting

$$\mathbf{f}:=\mathtt{BScurve}\mathtt{2}\mathtt{d}\left(M\right).$$

On the other hand, the calling sequence MaxCurvature2d(M) gives us the result

$$5.011177204,\left\{5.168088495\right\}.$$

Thus, the maximum curvature of $\mathcal{C}$ is ${\kappa }_{max}=5.011177204$ and this value is attained at the point $\mathbf{f}\left(5.168088495\right)=[3.53273441,2.298980995]$.

Let $\mathcal{L}$ be a curve with the parametrization $\mathbf{r}:t\mapsto [3sint,tcos(3t\left)\right]$, $t\in [-1,2]$, and let

$$t_0:=-1,t_1:=-0.6,t_2:=-0.2,t_3:=0.2,t_4:=0.6,t_5:=0.9,t_6:=1.3,t_7:=1.7,t_8:=2$$

be a partition of $[-1,2]$. We set $N:=[\mathbf{r}\left(t_0\right),\mathbf{r}\left(t_1\right),\mathbf{r}\left(t_2\right),\dots ,\mathbf{r}\left(t_8\right)]$. Letting ${\mathcal{C}}_1$ be the Bézier-spline curve that interpolates N, we get its parametrization

$$\mathbf{h}:=\mathtt{BScurve}\mathtt{2}\mathtt{d}\left(N\right).$$

Then, we derive

$$7.637332459,\left\{5.720743580\right\}$$

from the calling sequence MaxCurvature2d(N). It follows that ${\mathcal{C}}_1$ attains the maximum curvature ${\kappa }_{max}=7.637332459$ at the point $\mathbf{h}\left(5.720743580\right)=[2.76963535,-1.07851249]$.

The results from the two examples above are given in Figure 4.

The next examples are dealing with closed Bézier-spline curves. The datasets A and B in the first two examples are chosen to be symmetric to an axis, namely

$$\begin{array}{cc}\hfill A& :=\left[\right[3,1],[5,3],[4,4],[3,3],[2,4],[1,3],[3,1\left]\right];\hfill \\ \hfill B& :=\left[\right[1,4],[0.6,2],[2,0.4],[3.4,1],[2.6,2.8],[2.2,2.4],[4,1.6],[4.6,3],[3,4.4],[1,4\left]\right].\hfill \end{array}$$

The curves ${\mathcal{C}}_A$ and ${\mathcal{C}}_B$ that interpolate A and B, respectively, are given in Figure 5. The shapes of ${\mathcal{C}}_A$ and ${\mathcal{C}}_B$ can already be seen if we first display the interpolated points on the plane.

In the last two examples, we compute the maximum curvature of two closed Bézier-spline curves and show the maximum curvature points on these curves. The results are given in Figure 6. The interpolated points (on the right) in Figure 6 are chosen on the ellipse $x^2/9+y^2/4=1$ (depicted in cyan).

We sometimes want to compute the curvature of a Bézier-spline curve at a $p\in [0,n]$, so we should have a tool to do that. Algorithm 2 can be used to make such a tool.

Algorithm 2: Finding the curvature of a Bézier-spline curve at its given point

Input:  a set T of $(n+1)$ points in ${\mathbb{R}}^2$, a point $p\in [0,n]$; Output:  The curvature $\kappa$ at p of the Bézier-spline curve interpolating T; 1: if$p=n$then 2:    $i:=n$; 3: else 4:    $i:=\lfloor p\rfloor +1$; 5: end if 6: $\mathbf{v}:={\mathbf{f}}_i^{\prime }\left(t\right)$, $\mathbf{a}:={\mathbf{f}}_i^{″}\left(t\right)$     // ${\mathbf{f}}_i$ from BScurve2d or BScurve2dC; 7: $\mathbf{u}:=J\mathbf{v}$     // $J=\left[\begin{array}{cc}0& \hfill -1\\ 1& 0\end{array}\right]$; 8: $F:=|\mathbf{a}·\mathbf{u}|/v^3$    // $v:=\parallel \mathbf{v}\parallel$, $F:t\mapsto \kappa \left(t\right)$; 9; return $F\left(p\right)$;

It is easy to derive the Maple procedure from Algorithm 2 that we call Curvature2d or Curvature2dC, depending on whether ${\mathbf{f}}_i$ is from BScurve2d or from BScurve2dC. Similarly, the extension of this procedure for the three-dimensional case can be obtained with the curvature function from Algorithm 1.

5. Conclusions

Our procedures may also be convenient tools for non-Maple users to estimate the maximum curvature of parametric smooth curves and to approximate a curve from its interpolated points by a Bézier-spline curve, as similarly done in [13,14,15,16]. They could be used for doing calculations to establish conditions for non-singularity of tubular surfaces and offset curves associated to a Bézier-spline curve, according to the issues profoundly presented in [17,18]. Moreover, our parametrization derived from the closed-form solution to a linear system for the interpolation problem would be better than those derived from the existing direct or iterative methods to solve the system, as noted in [19,20].