Dynamic evaluation of exponential polynomial curves and surfaces via basis transformation

TL;DR

This paper introduces a robust and efficient basis transformation algorithm for the dynamic evaluation of exponential polynomial curves and surfaces, improving accuracy and reducing computational costs compared to traditional differential system methods.

Contribution

It presents an explicit basis transformation approach that enables robust, accurate, and computationally efficient dynamic evaluation of exponential polynomial curves and surfaces with changing parameters.

Findings

The new method provides robust and accurate evaluations for any parameter steps.

It reduces computational time by using a constant matrix for polynomial curves.

The approach outperforms conventional differential system solutions in efficiency.

Abstract

It is shown in "SIAM J. Sci. Comput. 39 (2017):B424-B441" that free-form curves used in computer aided geometric design can usually be represented as the solutions of linear differential systems and points and derivatives on the curves can be evaluated dynamically by solving the differential systems numerically. In this paper we present an even more robust and efficient algorithm for dynamic evaluation of exponential polynomial curves and surfaces. Based on properties that spaces spanned by general exponential polynomials are translation invariant and polynomial spaces are invariant with respect to a linear transform of the parameter, the transformation matrices between bases with or without translated or linearly transformed parameters are explicitly computed. Points on curves or surfaces with equal or changing parameter steps can then be evaluated dynamically from a start point using a pre-computed matrix. Like former dynamic evaluation algorithms, the newly proposed approach needs only arithmetic operations for evaluating exponential polynomial curves and surfaces. Unlike conventional numerical methods that solve a linear differential system, the new method can give robust and accurate evaluation results for any chosen parameter steps. Basis transformation technique also enables dynamic evaluation of polynomial curves with changing parameter steps using a constant matrix, which reduces time costs significantly than computing each point individually by classical algorithms.

Figures

Tables

Table 1. Table 1: Basis vectors and their elementary decompositions.

Basis vector	Elementary decomposition
${(1,t,\cos t,\sin t)}^T$ [32]	$U_1\left(t\right)\bigsqcup V\left(t\right)$
${(1,\cos t,\sin t,\dots ,\cos nt,\sin nt)}^T$ [26]	$U_0\left(t\right)\bigsqcup V\left(t\right)\bigsqcup \dots \bigsqcup V\left(nt\right)$
${(1,\cosh t,\sinh t,\dots ,\cosh nt,\sinh nt)}^T$ [28]	$U_0\left(t\right)\bigsqcup W\left(t\right)\bigsqcup \dots \bigsqcup W\left(nt\right)$
${(1,t,\cos t,\sin t,t\cos t,t\sin t)}^T$ [17]	$U_1\left(t\right)\bigsqcup U_1\left(t\right)\otimes V\left(t\right)$
${(1,t,\dots ,t^{n-2},\cos t,\sin t)}^T$ [5]	$U_{n-2}\left(t\right)\bigsqcup V\left(t\right)$
${(1,t,\dots ,t^{n-2},\cosh t,\sinh t)}^T$ [14]	$U_{n-2}\left(t\right)\bigsqcup W\left(t\right)$
${(1,\cosh t,\sinh t,\cos t,\sin t)}^T$ [4]	$U_0\left(t\right)\bigsqcup W\left(t\right)\bigsqcup V\left(t\right)$
${(1,t,\dots ,t^{n-5},\cosh t,\sinh t,\cos t,\sin t)}^T$ [30]	$U_{n-5}\left(t\right)\bigsqcup W\left(t\right)\bigsqcup V\left(t\right)$
${(1,\cos t,\sin t,t\cos t,t\sin t,\dots ,t^n\cos t,t^n\sin t)}^T$ [29]	$U_0\left(t\right)\bigsqcup U_n\left(t\right)\otimes V\left(t\right)$

Table 2. Table 2: Maximum deviations for evaluating the intrinsically defined curve.

#points	basis transformation	Taylor’s method
10	4.261E-14	2.916E+1
20	5.153E-14	3.026E-1
100	1.196E-13	1.971E-5
200	2.160E-13	3.087E-7
1000	6.407E-13	1.967E-11
2000	1.467E-12	3.884E-13
10000	4.606E-12	1.125E-12
20000	3.954E-12	1.306E-11

Equations

$$L=\left(\frac{d}{dt}-{\lambda }_0\right)\cdots \left(\frac{d}{dt}-{\lambda }_n\right),$$

$${\Phi }_a(t)\bigsqcup {\Phi }_b(t)={\left(a_0(t),\dots ,a_n(t),b_0(t),\dots ,b_m(t)\right)}^T.$$

$${\Phi }_a(t)\otimes {\Phi }_b(t)={\bigsqcup }_{i=0}^n{a}_i(t){\Phi }_b(t).$$

$$\frac{d}{dt}({\Phi }_a(t)\bigsqcup {\Phi }_b(t))=\mathrm{diag}(A,B)({\Phi }_a(t)\bigsqcup {\Phi }_b(t)),$$

$$\begin{array}[]{lcl}\frac{d}{dt}(\Phi_{a}(t)\otimes\Phi_{b}(t))&=&\Phi_{a}^{\prime}(t)\otimes\Phi_{b}(t)+\Phi_{a}(t)\otimes\Phi_{b}^{\prime}(t)\\ &=&(A\otimes I_{m+1})(\Phi_{a}(t)\otimes\Phi_{b}(t))+(I_{n+1}\otimes B)(\Phi_{a}(t)\otimes\Phi_{b}(t))\\ &=&(A\otimes I_{m+1}+I_{n+1}\otimes B)(\Phi_{a}(t)\otimes\Phi_{b}(t)).\end{array}$$

$$\left\{\begin{array}[]{lcl}\Phi^{\prime}(t)&=&A\Phi(t),\ \ \ t\in\mathbb{R}.\\ \Phi(t_{0})&=&\Phi_{0},\end{array}\right.$$

$$M_{U_{n}}^{h}=\left(\begin{array}[]{ccccc}1&0&0&\cdots&0\\ h&1&0&\cdots&0\\ h^{2}&2h&1&\cdots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ h^{n}&\binom{n}{1}h^{n-1}&\binom{n}{2}h^{n-2}&\cdots&1\end{array}\right).$$

$$M_{V}^{h}=\left(\begin{array}[]{cc}\cos h&-\sin h\\ \sin h&\cos h\end{array}\right)$$

$$M_{W}^{h}=\left(\begin{array}[]{cc}\cosh h&\sinh h\\ \sinh h&\cosh h\end{array}\right),$$

$$\begin{array}[]{lll}\sin h&=&h-\frac{h^{3}}{3!}+\frac{h^{5}}{5!}-\frac{h^{7}}{7!}+\frac{h^{9}}{9!}-\frac{h^{11}}{11!}+\frac{h^{13}}{13!}-\frac{h^{15}}{15!}\\ \cos h&=&1-\frac{h^{2}}{2!}+\frac{h^{4}}{4!}-\frac{h^{6}}{6!}+\frac{h^{8}}{8!}-\frac{h^{10}}{10!}+\frac{h^{12}}{12!}-\frac{h^{14}}{14!}-\frac{h^{16}}{16!}\\ \sinh h&=&h+\frac{h^{3}}{3!}+\frac{h^{5}}{5!}+\frac{h^{7}}{7!}+\frac{h^{9}}{9!}+\frac{h^{11}}{11!}+\frac{h^{13}}{13!}+\frac{h^{15}}{15!}\\ \cosh h&=&1+\frac{h^{2}}{2!}+\frac{h^{4}}{4!}+\frac{h^{6}}{6!}+\frac{h^{8}}{8!}+\frac{h^{10}}{10!}+\frac{h^{12}}{12!}+\frac{h^{14}}{14!}+\frac{h^{16}}{16!}\end{array}$$

$${\Phi }_1(t+h)\bigsqcup {\Phi }_2(t+h)=\mathrm{diag}(M_1,M_2)({\Phi }_1(t)\bigsqcup {\Phi }_2(t)).$$

$${\Phi }_1(t+h)\otimes {\Phi }_2(t+h)=(M_1\otimes M_2)({\Phi }_1(t)\otimes {\Phi }_2(t)),$$

$$\begin{array}[]{lcl}\Phi(u+h_{1},v+h_{2})&=&\Phi_{1}(u+h_{1})\otimes\Phi_{2}(v+h_{2})\\ &=&(M_{1}\otimes M_{2})(\Phi_{1}(u)\otimes\Phi_{2}(v))\\ &=&(M_{1}\otimes M_{2})\Phi(u,v).\end{array}$$

$$C_{A}=\left(\begin{array}[]{ccccc}1&0&0&\ldots&0\\ a_{1}&a_{0}&0&\ldots&0\\ a_{1}^{2}&2a_{0}a_{1}&a_{0}^{2}&\ldots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ a_{1}^{n}&\binom{n}{1}a_{0}a_{1}^{n-1}&\binom{n}{2}a_{0}^{2}a_{1}^{n-2}&\ldots&a_{0}^{n}\end{array}\right).$$

$${\Phi }_B((1-t)a+tb)=C_B{\Phi }_B(t),$$

$$C_{B}=\left(\begin{array}[]{cccc}c_{00}&c_{01}&\ldots&c_{0n}\\ c_{10}&c_{11}&\ldots&c_{1n}\\ \vdots&\vdots&\ddots&\vdots\\ c_{n0}&c_{n1}&\ldots&c_{nn}\end{array}\right)$$

$$\begin{array}[]{lcl}\Phi_{B}((1-t)a+tb)&=&\{[1-(1-t)a-tb]I+[(1-t)a+tb]E\}^{n}\mathbf{e}_{0}\\ &=&\{[(1-a)I+aE](1-t)+[(1-b)I+bE]t\}^{n}\mathbf{e}_{0}\\ &=&\sum_{l=0}^{n}[(1-b)I+bE]^{l}[(1-a)I+aE]^{n-l}\mathbf{e}_{0}B_{l,n}(t)\\ &=&\sum_{l=0}^{n}\mathbf{q}_{l}B_{l,n}(t),\end{array}$$

$$\begin{array}[]{lcl}\mathbf{q}_{l}&=&[(1-b)I+bE]^{l}[(1-a)I+aE]^{n-l}\mathbf{e}_{0}\\ &=&\sum_{i=0}^{l}E^{i}B_{i,l}(b)\sum_{j=0}^{n-l}E^{j}B_{j,n-l}(a)\mathbf{e}_{0}\\ &=&\sum_{k=0}^{n}\sum^{i+j=k}_{0\leq i\leq l,0\leq j\leq n-l}B_{i,l}(b)B_{j,n-l}(a)E^{k}\mathbf{e}_{0}\\ &=&\sum_{k=0}^{n}c_{kl}\mathbf{e}_{k}.\end{array}$$

$$\begin{array}[]{lcl}\Phi_{B}((1-t)a+tb)&=&(\begin{array}[]{cccc}\mathbf{q}_{0}&\mathbf{q}_{1}&\ldots&\mathbf{q}_{n}\end{array})\left(\begin{array}[]{c}B_{0,n}(t)\\ B_{1,n}(t)\\ \vdots\\ B_{n,n}(t)\end{array}\right)\\ &=&\left(\begin{array}[]{cccc}c_{00}&c_{01}&\ldots&c_{0n}\\ c_{10}&c_{11}&\ldots&c_{1n}\\ \vdots&\vdots&\ddots&\vdots\\ c_{n0}&c_{n1}&\ldots&c_{nn}\end{array}\right)\left(\begin{array}[]{c}B_{0,n}(t)\\ B_{1,n}(t)\\ \vdots\\ B_{n,n}(t)\end{array}\right)\\ &=&C_{B}\Phi_{B}(t).\end{array}$$

$$\left\{\begin{array}[]{lcl}X(t_{0})&=&X_{orig},\\ X(t_{i})&=&M_{h}X(t_{i-1}),\ \ \ i=1,2,\ldots\end{array}\right.$$

$$X(t_i)=M_{B}X(t_{i-1}),i=1,2,\dots$$

$$\begin{array}[]{lcl}t_{i}&=&a+(b-a)t_{i-1}\\ &=&a+a(b-a)+\cdots+a(b-a)^{i-1}+t_{0}(b-a)^{i}\\ &=&\left\{\begin{array}[]{ll}t_{0}+ia&\text{if}\ (b-a)=1\\ a\frac{1-(b-a)^{i}}{1-(b-a)}+t_{0}(b-a)^{i}&\text{otherwise}.\end{array}\right.\end{array}$$

$$X(u+h_1,v)=M_u^{h_1}X(u,v)$$

$$X(u,v+h_2)=M_v^{h_2}X(u,v).$$

$$\begin{array}[]{lcl}\frac{\partial X(u,v)}{\partial u}&=&A_{u}X(u,v),\\ \frac{\partial X(u,v)}{\partial v}&=&A_{v}X(u,v),\end{array}$$

$$\begin{array}[]{lcl}Q^{\prime}(t)&=&\frac{dX(u(t),v(t))}{dt}\\ &=&u^{\prime}(t)\frac{\partial X(u,v)}{\partial u}+v^{\prime}(t)\frac{\partial X(u,v)}{\partial v}\\ &=&M_{X}[u^{\prime}(t)A_{1}+v^{\prime}(t)A_{2}]\Phi(u(t),v(t)).\end{array}$$

$$\left\{\begin{array}[]{lcl}Q^{\prime}(t)&=&A_{\gamma}Q(t),\ \ \ t\in[\alpha,\beta],\\ Q(\alpha)&=&X(u(\alpha),v(\alpha)),\end{array}\right.$$

$$\begin{array}[]{lcl}Q(t_{i}+h)&=&X(u_{i}+\delta h,v_{i}+\eta h)\\ &=&M_{X}\Phi(u_{i}+\delta h,v_{i}+\eta h)\\ &=&M_{X}C_{u,v}^{\delta h,\eta h}\Phi(u_{i},v_{i})\\ &=&M_{X}C_{u,v}^{\delta h,\eta h}M_{X}^{-1}Q(t_{i})\\ &=&M_{u,v}^{\delta h,\eta h}Q(t_{i}),\end{array}$$

$$\displaystyle\textbf{r}(\theta)=\Bigg{(}\begin{array}[]{c}x(\theta)\\ y(\theta)\end{array}\Bigg{)}=\Bigg{(}\begin{array}[]{c}\int_{0}^{\theta}\rho(t)\cos tdt\\ \int_{0}^{\theta}\rho(t)\sin tdt\end{array}\Bigg{)}$$

$$\begin{array}[]{lcl}\textbf{r}(\theta)&=&\Big{(}\begin{array}[]{c}-1.494\\ 0.52\end{array}\Big{)}+\Big{(}\begin{array}[]{c}1.494\\ -0.52\end{array}\Big{)}\cos\theta+\Big{(}\begin{array}[]{c}0.52\\ 1.494\end{array}\Big{)}\sin\theta\\ &&+\Big{(}\begin{array}[]{c}-0.12\\ -1.494\end{array}\Big{)}\theta\cos\theta+\Big{(}\begin{array}[]{c}1.494\\ -0.12\end{array}\Big{)}\theta\sin\theta+\Big{(}\begin{array}[]{c}0.003\\ 0.06\end{array}\Big{)}\theta^{2}\cos\theta\\ &&+\Big{(}\begin{array}[]{c}-0.06\\ 0.003\end{array}\Big{)}\theta^{2}\sin\theta+\Big{(}\begin{array}[]{c}0\\ -0.001\end{array}\Big{)}\theta^{3}\cos\theta+\Big{(}\begin{array}[]{c}0.001\\ 0\end{array}\Big{)}\theta^{3}\sin\theta.\end{array}$$

Full text

{AMS}

65D17, 65D18, 65D25, 65L05

1 Introduction

The exponential polynomials that lie in the null spaces of constant coefficient linear differential operators have nice properties and they have often been used for the construction of curves and surfaces in the fields of CAGD (computer aided geometric design) [1, 9, 12]. The frequently used exponential polynomials are polynomials, trigonometric functions, hyperbolic functions or their mixtures. Besides polynomial curves and surfaces, typical curves and surfaces such as ellipses, cycloids, involutes, helices, etc. can be represented by exponential polynomials exactly [17, 21, 24, 32]. By choosing a proper parameter interval, normalized B-bases that are useful for optimal shape design can be obtained from the exponential polynomials [5, 18, 25, 27]. Algebraic trigonometric polynomials can be used to define curves intrinsically or design curves with Pythagorean hodographs [29, 23, 22].

Many algorithms have been given in the literature to evaluate polynomial curves and surfaces. The de Casteljau algorithm or the rational de Casteljau algorithm can be employed to robustly evaluate single points on Bézier or rational Bézier curves [9, 6]. The Horner algorithm, the VS algorithm, et al. can be used to evaluate polynomial or Bézier curves with even lower complexity [7, 3]. If for rendering or machining purposes, sequences of points have to be evaluated in an efficient way [8, 10]. Particularly, the forward differencing approaches have been successfully employed for fast rendering of Bézier or NURBS (Non-Uniform Rational B-spline) curves and surfaces [15, 16].

Differently from polynomial curves and surfaces that can be evaluated by arithmetic operations, points on curves and surfaces which are constructed by transcendental functions or mixtures of polynomials and transcendental functions have to be evaluated by inquiring pre-computed special function tables or loading special mathematical libraries. Though this seems feasible for many modern computing machines [13, 20], evaluating general exponential polynomial curves and surfaces by only arithmetic operations without any pre-computed tables or special math library can have its own advantages. Particularly, the speed and efficiency of evaluation play important roles in the fields of CNC machining and interactive rendering.

Recently, we have shown that free-form curves defined in various spaces in CAGD are the solutions of linear differential systems and points and derivatives on the curves can be obtained by solving the linear differential systems numerically [31]. Particularly, when the parameter step is fixed, points on a free-form curve can be evaluated dynamically by multiplying a pre-computed constant matrix with prior points. Iso-parameter curves on a surface can also be dynamically evaluated by establishing a linear differential system for each iso-parameter curve. The method is simple and universal and points on polynomial as well as transcendental curves and surfaces can be evaluated with only arithmetic operations. However, the evaluation accuracy varies much when the differential systems have been solved by different numerical methods or the parameter step has been chosen different values. If the constant matrix for dynamic evaluation is given by the exponential of the coefficient matrix of a linear differential system, careful attention should be paid for robust computation of the matrix [19].

Instead of solving linear differential systems numerically, in this paper we derive the constant matrices for dynamic evaluation of curves and surfaces using basis transformation. This is based on the fact that spaces spanned by the exponential polynomial basis are invariant with respect to the translation of the parameter while curves and surfaces used for CAGD are usually constructed by exponential polynomials. By using identities of exponential polynomials, we compute explicitly the transformation matrices between bases with or without the translation of the parameter. Combined with control points, a constant matrix for evaluating points on an exponential polynomial curve with equal parameter steps is derived. It is also noticed that a polynomial space of degree no more than a given number is even invariant with respect to any linear transformation of the parameter. The matrix for polynomial basis transformation can then be used to evaluate points on polynomial curves with changing parameter steps. Based on basis transformation, a family of surface curves with various iso-parameters can be evaluated dynamically using a single matrix and surface curves with skew parametrization can also be evaluated dynamically with a pre-computed constant matrix.

The rest of the paper is structured as follows. In Section 2 we present explicit formulae for computing basis transformation for exponential polynomials. In Section 3 robust algorithm for evaluating points on general exponential polynomial curves with a fixed parameter step, dynamic algorithm for evaluating polynomial curves with changing parameter steps, and dynamic algorithms for evaluating iso-parameter curves on surfaces or surface curves with skew parametrization will be given. Examples and comparisons with some known methods for curve and surface evaluation are given in Section 4. Section 5 concludes the paper with a brief summary of our work.

2 Basis transformation for spaces composed of exponential polynomials

As parametric curves and surfaces are usually defined by basis functions together with coefficients or control points, a parametric curve or surface can then be evaluated efficiently by exploring distinguished properties of the basis. This section presents explicit transformation formulae for exponential polynomial basis which will be used for robust and efficient curve or surface evaluation in next section.

2.1 Spaces spanned by exponential polynomials

Suppose that a linear differential operator $L$ with constant coefficients is given by

[TABLE]

where $\lambda_{i}\in\mathbb{C}$, $i=0,1,\ldots,n$, and $\Lambda=\{\lambda_{0},\lambda_{1},\ldots,\lambda_{n}\}$ is closed under conjugation. A function $f(t)$ that satisfies $Lf(t)=0$ is referred an exponential polynomial. Let $\Omega$ be the null space of the linear differential operator. From the knowledge of differential equation [2] we know that $\Omega=\text{span}\{\phi_{0}(t),\phi_{1}(t),\ldots,\phi_{n}(t)\}$, where $\phi_{i}(t)$, $i=0,1,\ldots,n$, are the basis functions of the space. Based on the definition of exponential polynomials we have the following proposition.

Proposition 2.1**.**

*Suppose that $L$ is a constant coefficient linear differential operator and $\Omega$ is the null space of the operator. Let $h$ be an arbitrary given real number. If function $f(t)$ satisfies $Lf(t)=0$, it yields that $f^{\prime}(t)\in\Omega$ and $f(t+h)\in\Omega$. *

Proposition 2.1 states that the null space of a linear differential operator is closed with respect to a differentiation and the space is also invariant with respect to any translation of the parameter.

Before deriving formulae for basis transformation, we present the definition of union or product of two sets of bases. Assume $\Phi_{a}(t)=\left(a_{0}(t),a_{1}(t),\ldots,a_{n}(t)\right)^{T}$ and $\Phi_{b}(t)=\left(b_{0}(t),b_{1}(t),\ldots,b_{m}(t)\right)^{T}$, where the capital ’T’ means the transpose of a vector or matrix. The (ordered) union of $\Phi_{a}(t)$ and $\Phi_{b}(t)$ is given by

[TABLE]

Let $a(t)\Phi_{b}(t)=\left(a(t)b_{0}(t),\ldots,a(t)b_{m}(t)\right)^{T}$. The product of $\Phi_{a}(t)$ and $\Phi_{b}(t)$ is obtained as

[TABLE]

Just like the precedence of ’$\times$’ over ’$+$’, we assume the operation ’$\otimes$’ has precedence over ’$\sqcup$’. Thus, $\Phi_{1}(t)\sqcup\Phi_{2}(t)\otimes\Phi_{3}(t)$ has the same meaning as $\Phi_{1}(t)\sqcup(\Phi_{2}(t)\otimes\Phi_{3}(t))$.

Proposition 2.2**.**

*Suppose spaces spanned by basis $\Phi_{a}(t)$ or $\Phi_{b}(t)$ are closed with respect to a differentiation. Then, spaces spanned by the union $\Phi_{a}(t)\sqcup\Phi_{b}(t)$ or by the product $\Phi_{a}(t)\otimes\Phi_{b}(t)$ are also closed with respect to a differentiation. *

Proof 2.3**.**

Because spaces spanned by $\Phi_{a}(t)$ or $\Phi_{b}(t)$ are closed with respect to a differentiation, there exist matrix $A$ of order $n+1$ and matrix $B$ of order $m+1$ such that $\Phi_{a}^{\prime}(t)=A\Phi_{a}(t)$ and $\Phi_{b}^{\prime}(t)=B\Phi_{b}(t)$. Then, the derivative of $\Phi_{a}(t)\sqcup\Phi_{b}(t)$ can be computed as

[TABLE]

where \operatorname{diag}(A,B)=\left(\begin{array}[]{cc}A&0\\ 0&B\end{array}\right). Let $I_{n+1}$ and $I_{m+1}$ be the identity matrices of order $n+1$ or order $m+1$, respectively. The derivative of $\Phi_{a}(t)\otimes\Phi_{b}(t)$ is computed by

[TABLE]

*Therefore, the spaces spanned by $\Phi_{a}(t)\sqcup\Phi_{b}(t)$ or by $\Phi_{a}(t)\otimes\Phi_{b}(t)$ are also closed under differentiation. This completes the proof. *

In the following we assume that the basis functions are real. Particularly, we assume that the bases are obtained by unions or products of a few elementary basis vectors. Let $U_{n}(t)=(1,t,t^{2},\ldots,t^{n})^{T}$, $V(t)=(\cos(t),\sin(t))^{T}$ and $W(t)=(\cosh(t),\sinh(t))^{T}$. The basis functions for free-form curves and surfaces in CAGD can usually be obtained by recursive compositions or tensor products of the elementary bases. Several popular basis vectors for construction of free-form curves and surfaces in CAGD and their elementary decomposition can be found in Table 1.

2.2 Transformation of general exponential polynomial basis with parameter translation

In this subsection we show that the space spanned by general exponential polynomial basis is invariant with respect to a translation of the parameter and a simple and robust method for computing the transformation matrix between different bases will be presented.

Proposition 2.4**.**

*Suppose $\Omega=\text{span}\{\phi_{0}(t),\phi_{1}(t),\ldots,\phi_{n}(t)\}$ is closed with respect to a differentiation. Let $h$ be an arbitrary real number. It yields that $\Omega=\text{span}\{\phi_{0}(t+h),\phi_{1}(t+h),\ldots,\phi_{n}(t+h)\}$. Meanwhile, there exists a matrix $C_{h}$ such that $\Phi(t+h)=C_{h}\Phi(t)$. *

Proof 2.5**.**

Because $\phi_{i}^{\prime}(t)\in\Omega$, $i=0,1,\ldots,n$, there exists a matrix $A$ such that $\Phi(t)$ satisfies a linear differential system

[TABLE]

*From Equation (3), $\Phi(t)$ can be represented as $\Phi(t)=e^{A(t-t_{0})}\Phi_{0}$. Therefore, we have $\Phi(t+h)=e^{Ah}\Phi(t)$. Let $C_{h}=e^{Ah}$. Since $\text{det}(C_{h})\neq 0$ and $\Phi(t)=C_{h}^{-1}\Phi(t+h)$, we have $\Omega=\text{span}\{\phi_{0}(t+h),\phi_{1}(t+h),\ldots,\phi_{n}(t+h)\}$. This proves the proposition. *

Proposition 2.6**.**

*Suppose $\Omega=\textrm{span}\{\phi_{0}(u,v),\phi_{1}(u,v),\ldots,\phi_{L}(u,v)\}$ and $\frac{\partial\phi_{l}(u,v)}{\partial u}\in\Omega$, $\frac{\partial\phi_{l}(u,v)}{\partial v}\in\Omega$, $l=0,1,\ldots,L$. The space $\Omega$ is invariant when parameter $u$ or parameter $v$ or both of the two parameters have been translated. *

Proof 2.7**.**

We first prove that the space $\Omega$ is invariant when the parameter $u$ or $v$ has been translated. Let $\Phi(u,v)=(\phi_{0}(u,v),\phi_{1}(u,v),\ldots,\phi_{L}(u,v))^{T}$. Because $\frac{\partial\phi_{l}(u,v)}{\partial u}\in\Omega$, $l=0,1,\ldots,L$, there exists matrix $A_{1}$ such that the basis vector $\Phi(u,v)$ satisfy $\frac{\partial\Phi(u,v)}{\partial u}=A_{1}\Phi(u,v)$. Therefore, we have $\Phi(u,v)=e^{A_{1}(u-u_{0})}\Phi(u_{0},v)$ for any selected real number $u_{0}$. From this expression of $\Phi(u,v)$, we have $\Phi(u+h_{1},v)=e^{A_{1}h_{1}}\Phi(u,v)$. By the same reason, we have $\frac{\partial\Phi(u,v)}{\partial v}=A_{2}\Phi(u,v)$ and $\Phi(u,v+h_{2})=e^{A_{2}h_{2}}\Phi(u,v)$. Denote $C_{u}^{h_{1}}=e^{A_{1}h_{1}}$ and $C_{v}^{h_{2}}=e^{A_{2}h_{2}}$. Because matrices $C_{u}^{h_{1}}$ and $C_{v}^{h_{2}}$ are non-singular, it implies that both $\Phi(u+h_{1},v)$ and $\Phi(u,v+h_{2})$ are basis vectors of the space $\Omega$.

*Now we prove that the space $\Omega$ is invariant when both parameters $u$ and $v$ have been translated. In fact, $\Phi(u+h_{1},v+h_{2})=C_{u}^{h_{1}}\Phi(u,v+h_{2})=C_{u}^{h_{1}}C_{v}^{h_{2}}\Phi(u,v)$. Let $C_{u,v}^{h_{1},h_{2}}=C_{u}^{h_{1}}C_{v}^{h_{2}}$. Since $\det(C_{u,v}^{h_{1},h_{2}})=\det(C_{u}^{h_{1}})\det(C_{v}^{h_{2}})\neq 0$, we conclude that $\Phi(u+h_{1},v+h_{2})$ is also the basis vector of space $\Omega$. *

Though the transformation matrices between bases with or without translation of the parameters are defined by exponentials of constant matrices, accurate and efficient evaluation of exponentials of matrices is not a trivial task [19]. We propose to compute the transformation matrices for exponential polynomial bases with translated parameters directly. Particularly, we derive the transformation matrices for polynomials, trigonometric functions or hyperbolic functions using the identities of the functions first and then compute the transformation matrices for even more general basis through their elementary decompositions.

Assume $U_{n}(t)$, $V(t)$ and $W(t)$ are the basis vectors as given in Subsection 2.1. With simple computation, the transformation for basis vector $U_{n}(t)$ can be obtained as $U_{n}(t+h)=M_{U_{n}}^{h}U_{n}(t)$, where

[TABLE]

Similarly, the transformation matrices for $V(t)$ and $W(t)$ are

[TABLE]

and

[TABLE]

respectively. It yields that $V(t+h)=M_{V}^{h}V(t)$ and $W(t+h)=M_{W}^{h}W(t)$.

For a small parameter step $h$, the values of $\cos h$, $\sin h$, $\cosh h$ and $\sinh h$ can be computed efficiently with only arithmetic operations by Taylor expansion:

[TABLE]

The above expressions can be evaluated by Horner algorithm in practice. From the definition of $M_{V}^{h}$ and $M_{W}^{h}$ we know that $M_{V}^{h}=(M_{V}^{h/K})^{K}$ and $M_{W}^{h}=(M_{W}^{h/K})^{K}$. If the parameter step $h$ is larger than a threshold, for example $0.1$, we can choose a proper integer $K$ and compute the elements of matrix $M_{V}^{h/K}$ or $M_{W}^{h/K}$ first, and then compute the matrix $M_{V}^{h}$ or $M_{W}^{h}$ by matrix multiplication. Similarly, transformation matrices with other parameter steps can be computed by $M^{2h}=(M^{h})^{2}$, $M^{3h}=(M^{h})^{3}$, etc.

Assume the transformation matrices for $\Phi_{1}(t)$ or $\Phi_{2}(t)$ are $M_{1}$ and $M_{2}$, respectively. The transformation matrix for $\Phi_{1}(t)\sqcup\Phi_{2}(t)$ can be computed as

[TABLE]

The transformation matrix for $\Phi_{1}(t)\otimes\Phi_{2}(t)$ is obtained as follows

[TABLE]

where $M_{1}\otimes M_{2}$ is also known as the Kronecker product of two matrices.

Similar to the product of two univariate bases, the transformation matrix for the tensor product basis of a surface can be computed easily. Suppose that $\Phi(u,v)=\Phi_{1}(u)\otimes\Phi_{2}(v)$, $\Phi_{1}(u+h_{1})=M_{1}\Phi_{1}(u)$ and $\Phi_{2}(v+h_{2})=M_{2}\Phi_{2}(v)$, the transformation matrix for $\Phi(u,v)$ is computed by

[TABLE]

By Equations (4), (5) and (6), the transformation matrices for even more complicated basis will be computed accurately and robustly.

2.3 Transformation of polynomial basis with linear transformation of the parameter

In this subsection we show that space of polynomials of degree no more than a given number is invariant with respect to a linear transformation of the parameter. Particularly, the transformation matrix for Bernstein basis will be given.

Proposition 2.8**.**

*Let $\mathbb{P}_{n}(t)=\text{span}\{1,t,t^{2},\ldots,t^{n}\}$. The space $\mathbb{P}_{n}(t)$ is invariant with respect to a linear transformation of the parameter. *

Proof 2.9**.**

Let $U_{n}(t)=(1,t,t^{2},\ldots,t^{n})^{T}$. Assume $a_{0}\neq 0$ and $a_{1}$ are real numbers. Replacing $t$ within $U_{n}(t)$ by $a_{0}t+a_{1}$, we have $U_{n}(a_{0}t+a_{1})=C_{A}U_{n}(t)$, where

[TABLE]

*Because $\det(C_{A})=a_{0}^{n(n+1)/2}\neq 0$, we have $U_{n}(t)=C_{A}^{-1}U_{n}(a_{0}t+a_{1})$. Therefore, $U_{n}(a_{0}t+a_{1})$ is another set of basis of the space $\mathbb{P}_{n}(t)$. The proposition is proven. *

Besides power basis, another popular basis used for polynomial curve and surface modeling is Bernstein basis. We derive transformation matrix for Bernstein basis under a linear transformation of the parameter. The transformation matrix will be used for dynamic evaluation of Bézier curves and surfaces with constant or changing parameter steps.

Proposition 2.10**.**

Assume $\Phi_{B}(t)=(B_{0,n}(t),B_{1,n}(t),\ldots,B_{n,n}(t))^{T}$, where $B_{i,n}(t)=\frac{n!}{i!(n-i)!}t^{i}(1-t)^{n-i}$, $i=0,1,\ldots,n$, are Bernstein basis functions. Assume $a\neq b$ are real numbers. Then the basis vector satisfies

[TABLE]

where

[TABLE]

*and $c_{kl}=\sum\limits^{i+j=k}_{0\leq i\leq l,0\leq j\leq n-l}B_{i,l}(b)B_{j,n-l}(a)$, $0\leq k,l\leq n$. *

Proof 2.11**.**

Let $\mathbf{e}_{0}=(1,0,0,\ldots,0)^{T}$, $\mathbf{e}_{1}=(0,1,0,\ldots,0)^{T}$, $\ldots$, $\mathbf{e}_{n}=(0,0,0,\ldots,1)^{T}$. The basis vector $\Phi_{B}(t)$ can be represented as $\Phi_{B}(t)=\sum_{i=0}^{n}\mathbf{e}_{i}B_{i,n}(t)$. Assume $I$ and $E$ are identity or shift operators which satisfy $I\mathbf{e}_{i}=\mathbf{e}_{i}$ and $E\mathbf{e}_{i}=\mathbf{e}_{i+1}$. The basis vector can be reformulated as $\Phi_{B}(t)=[(1-t)I+tE]^{n}\mathbf{e}_{0}$. Then we have

[TABLE]

where

[TABLE]

As $\mathbf{q}_{l}=\sum_{k=0}^{n}c_{kl}\mathbf{e}_{k}=(c_{0l},c_{1l},\ldots,c_{nl})^{T}$, $l=0,1,\ldots,n$, representing $\Phi_{B}((1-t)a+tb)$ in matrix form, we have

[TABLE]

*This proves the proposition. *

3 Dynamic evaluation of exponential polynomial curves and surfaces

In this section we show that curves and surfaces defined by basis and control points can be evaluated dynamically by applying the basis transformation recursively. If a linear differential system is available, the derivatives of the curve or surface at the evaluated points can be obtained simultaneously.

3.1 Dynamic evaluation of exponential polynomial curves

Suppose a free-form curve is defined by $P(t)=\sum_{i=0}^{n}P_{i}\phi_{i}(t)$, where $P_{i}\in\mathbb{R}^{d}$, $i=0,1,\ldots,n$, are the control points. Let $\Phi(t)=(\phi_{0}(t),\phi_{1}(t),\ldots,\phi_{n}(t))^{T}$. The curve can be represented as $P(t)=M_{P}\Phi(t)$, where $M_{P}=(P_{0},P_{1},\ldots,P_{n})$ and $\Phi(t)$ is also referred as the normal curve in $\mathbb{R}^{n+1}$ [11]. If $d<n+1$, we first lift all control points in space $\mathbb{R}^{n+1}$ by adding additional coordinates as that presented in [31]. Assume the lifted curve is $X(t)=\sum_{i=0}^{n}X_{i}\phi_{i}(t)$, where $X_{i}\in\mathbb{R}^{n+1}$, $i=0,1,\ldots,n$. When a point on $X(t)$ has been evaluated, the point on $P(t)$ will be obtained just by choosing the first few coordinates.

Suppose that the space $\Omega=\text{span}\{\phi_{0}(t),\phi_{1}(t),\ldots,\phi_{n}(t)\}$ is closed under a differentiation and the matrix $M_{X}=(X_{0},X_{1},\ldots,X_{n})$ is nonsingular. The curve $X(t)$ is formulated as $X(t)=M_{X}\Phi(t)$. From Proposition 2.4 we know that the basis vector satisfies $\Phi(t+h)=C_{h}\Phi(t)$. Therefore, all points on curve $X(t)$ satisfy $X(t+h)=M_{X}C_{h}\Phi(t)$. Substituting $\Phi(t)=M_{X}^{-1}X(t)$, we have $X(t+h)=M_{h}X(t)$, where $M_{h}=M_{X}C_{h}M_{X}^{-1}$. Suppose an initial point at parameter $t=t_{0}$ is known, points with a constant time step will be computed dynamically as follows

[TABLE]

where $t_{i}=t_{0}+hi$. From Equation (3) we know that $\Phi^{\prime}(t)=A\Phi(t)$. Then we have $X^{\prime}(t)=A_{c}X(t)$, where $A_{c}=M_{X}AM_{X}^{-1}$. The derivatives at the evaluated points are obtained as $X^{\prime}(t_{i})=A_{c}X(t_{i})$, $X^{\prime\prime}(t_{i})=A_{c}X^{\prime}(t_{i})=A_{c}^{2}X(t_{i})$, etc.

As discussed in Subsection 2.3 the polynomial spaces are invariant under a linear transformation of the parameter. Points with changing parameter steps on a polynomial curve can be evaluated by using a fixed iteration matrix. Assume $\Phi_{B}(t)$ be the basis vector as given in Proposition 2.10. By applying a linear transformation that maps interval $[0,1]$ to $[a,b]$, the basis vector becomes $\Phi_{B}((1-t)a+tb)=C_{B}\Phi_{B}(t)$. A polynomial curve $X(t)=M_{X}\Phi_{B}(t)$ satisfies $X((1-t)a+tb)=M_{B}X(t)$, where $M_{B}=M_{X}C_{B}M_{X}^{-1}$. Starting from an initial point $X(t_{0})$, a sequence of points on curve $X(t)$ will be computed by

[TABLE]

where

[TABLE]

A linear differential system that represents the Bézier curves has been given in [31]. By the differential system, the derivatives at any evaluated point on a Bézier curve can be obtained just by multiplying the coefficient matrix with the point.

3.2 Dynamic evaluation of exponential polynomial surfaces

Similar to free-form curves, a bivariate surface can also be reformulated as the solution to a linear differential system when the space spanned by the basis functions is closed with respect to the partial differentials. Iso-parameter curves and surface curves with skew parametrization can be computed dynamically via basis transformation.

$\bullet$ Dynamic evaluation of iso-curves of free-form surfaces. Suppose a surface is given by $X(u,v)=\sum_{l=0}^{L}X_{l}\phi_{l}(u,v)$, where $X_{l}\in\mathbb{R}^{L+1}$, $l=0,1,\ldots,L$, are the control points and $\phi_{l}(u,v)$, $l=0,1,\ldots,L$, are the basis or blending functions. Assume the matrix $M_{X}=(X_{0},X_{1},\ldots,X_{L})$ is nonsingular and the space $\Omega=\text{span}\{\phi_{0}(u,v),\phi_{1}(u,v),\ldots,\phi_{L}(u,v)\}$ is closed with respect to the partial differentiations $\frac{\partial}{\partial u}$ and $\frac{\partial}{\partial v}$. Let $\Phi(u,v)=(\phi_{0}(u,v),\phi_{1}(u,v),\ldots,\phi_{L}(u,v))^{T}$. The surface is reformulated as $X(u,v)=M_{X}\Phi(u,v)$. From Proposition 2.6 we know that $\Phi(u+h_{1},v)=C_{u}^{h_{1}}\Phi(u,v)$ and $\Phi(u,v+h_{2})=C_{v}^{h_{2}}\Phi(u,v)$. Let $M_{u}^{h_{1}}=M_{X}C_{u}^{h_{1}}M_{X}^{-1}$ and $M_{v}^{h_{2}}=M_{X}C_{v}^{h_{2}}M_{X}^{-1}$. The iso-parameter curves or points on all $u$-curves or all $v$-curves can be dynamically evaluated by

[TABLE]

or

[TABLE]

From the proof of Proposition 2.6 we also know that $\frac{\partial\Phi(u,v)}{\partial u}=A_{1}\Phi(u,v)$ and $\frac{\partial\Phi(u,v)}{\partial v}=A_{2}\Phi(u,v)$. Therefore, the derivatives of $X(u,v)$ can be computed by

[TABLE]

where $A_{u}=M_{X}A_{1}M_{X}^{-1}$ and $A_{v}=M_{X}A_{2}M_{X}^{-1}$. Besides the first order derivatives, higher order derivatives can also be computed by multiplying these two matrices repeatedly. For example, $\frac{\partial^{2}X(u,v)}{\partial u^{2}}=A_{u}^{2}X(u,v)$, $\frac{\partial^{2}X(u,v)}{\partial u\partial v}=A_{v}A_{u}X(u,v)$.

$\bullet$ Dynamic evaluation of surface curves with skew parametrization. In addition to the iso-parameter curves, curves with skew parametrization on a surface can satisfy differential equations and can be evaluated dynamically too. Assume $\gamma(t)=(u(t),v(t))$ is a tangent smooth curve in the parameter domain. A surface curve is obtained as $Q(t)=X(u(t),v(t))$. The derivative of the surface curve is

[TABLE]

Substituting $\Phi(u(t),v(t))=M_{X}^{-1}X(u(t),v(t))$ into above equation, we have

[TABLE]

where $A_{\gamma}=M_{X}[u^{\prime}(t)A_{1}+v^{\prime}(t)A_{2}]M_{X}^{-1}$. In particular, if $u(t)$ and $v(t)$ are linear functions, i.e., $u(t)=u_{0}+\delta(t-\alpha)$ and $v(t)=v_{0}+\eta(t-\alpha)$, it yields that $u^{\prime}(t)=\delta$, $v^{\prime}(t)=\eta$, and $A_{\gamma}$ is a constant matrix.

Assume $u(t)$ and $v(t)$ are linear functions, we evaluate a sequence of points on curve $Q(t)=X(u(t),v(t))$ by basis transformation. Suppose that the point $Q(t_{i})=X(u_{i},v_{i})$ is known, we compute $Q(t_{i}+h)$ as follows

[TABLE]

where $C_{u,v}^{\delta h,\eta h}=C_{u}^{\delta h}C_{v}^{\eta h}$ is the basis transformation matrix as defined in Proposition 2.6. When the matrices $C_{u,v}^{\delta h,\eta h}$ and $M_{u,v}^{\delta h,\eta h}=M_{X}C_{u,v}^{\delta h,\eta h}M_{X}^{-1}$ have been obtained, points and derivatives of the curve $Q(t)$ will be evaluated dynamically by Equations (13) and (12).

3.3 Evaluation of curves and surfaces with combined parameter steps

For efficiency of rendering and machining, points on a curve or surface may be evaluated with non-constant or adaptive parameter steps. For polynomial curves and surfaces, Equation (9) can be employed to evaluate points with changing parameter steps. For general exponential curves and surfaces one may compute transformation matrices with different parameter steps first and then compute points on curves or surfaces using combinations of these transformation matrices.

From Proposition 2.4 we know that the transformation of a basis vector satisfies $\Phi(t+h_{1}+h_{2})=C_{h_{1}+h_{2}}\Phi(t)=C_{h_{1}}C_{h_{2}}\Phi(t)=C_{h_{2}}C_{h_{1}}\Phi(t)$. It is also verified that the transformation matrices given by Equation (8) satisfies $M_{h_{1}+h_{2}}=M_{h_{1}}M_{h_{2}}=M_{h_{2}}M_{h_{1}}$. Then, if one or more points should be added within a curve segment, we can begin with any point on the curve and compute additional points using transformation matrices with smaller parameter steps.

A surface $X(u,v)=M_{X}\Phi(u,v)$ can be evaluated along $u$-curves, $v$-curves, or curves with skew parametrization using different transformation matrices. From Proposition 2.6 we know that the transformation matrices for $\Phi(u,v)$ satisfy $C_{u,v}^{h_{1},h_{2}}=C_{u}^{h_{1}}C_{v}^{h_{2}}=C_{v}^{h_{2}}C_{u}^{h_{1}}$. Then the transformation matrices for surface points satisfy $M_{u,v}^{h_{1},h_{2}}=M_{u}^{h_{1}}M_{v}^{h_{2}}=M_{v}^{h_{2}}M_{u}^{h_{1}}$. This implies that $X(u+h_{1},v+h_{2})$ can be computed by a transformation from $X(u,v)$ directly or through intermediate points like $X(u+h_{1},v)$ or $X(u,v+h_{2})$. Using combinations of matrices $M_{u}^{h_{1}}$, $M_{v}^{h_{2}}$ and $M_{u,v}^{h_{1},h_{2}}$, one can compute curves with even more complicated parameter steps on a surface.

4 Examples and comparisons

The proposed algorithm for curve and surface evaluation was implemented using C++ on a laptop with Intel(R) Core(TM) i7-4910MQ CPU@2.90GHz 2.89GHz and 8G RAM. All numbers are represented in double precision. Comparisons between the proposed method and some known algorithms will be given.

Example 1. In the first example we evaluate an intrinsically defined planar curve

[TABLE]

where $\rho(t)$ represents the curvature radius of a planar curve and $\theta$ is the angle between the tangent direction of the curve and the positive direction of $x$-axis. Just as in [29], we choose $\rho(t)=0.001t^{3}-0.06t^{2}+1.5t+0.4$. The Cartesian coordinates of the curve are given by

[TABLE]

Just like that presented in [31], we lift the curve from $\mathbb{R}^{2}$ to $\mathbb{R}^{9}$. Assume the lifted coefficient vectors are $X_{i}\in\mathbb{R}^{9}$, $i=0,1,\ldots,8$. The lifted curve is obtained as $X(\theta)=M_{X}\Phi(\theta)$, where $M_{X}=(X_{0},X_{1},\ldots,X_{8})$ and $\Phi(\theta)$ is the corresponding basis vector. From Section 2 we know that $\Phi(\theta)=U_{0}(\theta)\sqcup U_{3}(\theta)\otimes V(\theta)$, where $U_{n}(\theta)$ and $V(\theta)$ are elementary bases as defined in Subsection 2.1. Based on Equations (4) and (5), the transformation of the basis vector with a translated parameter step is obtained as $\Phi(\theta+h)=C_{h}\Phi(\theta)$, where $C_{h}=\operatorname{diag}(M_{U_{0}}^{h},M_{U_{3}}^{h}\otimes M_{V}^{h})$. Furthermore, we compute $M_{h}=M_{X}C_{h}M_{X}^{-1}$. By this matrix, points with a fixed parameter step $h$ are computed consequently according to Equation (8). When a point $X(\theta)$ has been computed, the point $\textbf{r}(\theta)$ is obtained by choosing the first two coordinates.

In our experiments, we compute points on the curve that is defined on the parameter interval $[0,8\pi]$. Particularly, the start point is obtained as $X(0)=X_{0}+X_{1}$. The parameter step is chosen as $h=\frac{8\pi}{m}$ when $m$ points are to be computed on the curve. As $\rho(\theta)$ increases when $\theta$ changes from 0 to $8\pi$, the deviations from the evaluated points to the exact ones may increase too. To measure the accuracy of the proposed evaluation method, we compute the Euclidean distance from the last evaluated point to point $\mathbf{r}(8\pi)$ which is computed by loading the math library. From [31] we know that the lifted curve $X(\theta)$ satisfies a linear differential system and the points on the curve can be evaluated by employing Taylor’s method or implicit mid-point scheme to solve the differential system. Since the implicit mid-point scheme has only quadratic precision for solving a linear differential system, we compare the results by the proposed method only with Taylor’s method. The maximum deviations for evaluating points by basis transformation or by Taylor’s method with various choices of number $m$ are given in Table 2. From the table we can see that the evaluation accuracy by Taylor’s method (expansion order $s=6$) depends heavily on the parameter steps while the presented basis transformation approach can always give high accuracy results for various choices of point numbers.

Example 2. The second example is about dynamic evaluation of a Bézier curve with fixed or changing parameter steps.

Figure 1 illustrates a planar Bézier curve of degree 8. To evaluate the curve by the proposed dynamic algorithm, we lift the curve from $\mathbb{R}^{2}$ to $\mathbb{R}^{9}$. Assume the lifted Bézier curve is $X(t)=\sum_{i=0}^{8}X_{i}B_{i,8}(t)$. It can then be represented as $X(t)=M_{X}\Phi_{B}(t)$, where $M_{X}$ and $\Phi_{B}(t)$ are the coefficient matrix or the basis vector, respectively. For any two distinctive real numbers $a$ and $b$, we compute a basis transformation matrix $C_{B}$ by Equation (7) and then a curve transformation matrix $M_{B}$ by Equation (9). By matrix $M_{B}$ we compute points in $\mathbb{R}^{9}$ and obtain points in the plane.

We first compute points by choosing $X(t_{0})=X_{0}$, $a=0.01$ and $b=1.01$. The obtained points are plotted in Figure 1(a). Because $b-a=1$, a sequence of points with a fixed forward parameter step have been obtained. By choosing $a=0.01$ and $b=1.005$, a sequence of points with decreasing parameter steps have been obtained by applying Equation (9) recursively. See Figure 1(b) for the computed points starting from $X_{0}$. Similarly, points with increasing parameter steps can be obtained when we choose $a=0.01$ and $b=1.015$. See Figure 1(c) for the evaluated points. When we choose $a=-0.005$, $b=0.99$ and $X(t_{0})=X_{8}$, a sequence of points with decreasing parameter steps can be computed starting from the right end point of the curve. See Figure 1(d).

From the above results we can see that different basis transformations can lead to different sampling speeds or directions on the curve. One can then tune the sampling speed or sampling direction adaptively by choosing various transformation matrices. Because a linear differential system can be constructed from each given Bézier curve [31], derivatives at the sampled points can be evaluated directly using the differential system. As analyzed in [31], dynamic evaluation of Bézier curves needs much less time than evaluating points individually using classical de Casteljau algorithm, even though both of the two algorithms have $O(n^{2})$ time complexity.

Example 3. In the third example we show how to evaluate piecewise smooth curves on a tensor product Bézier surface by the proposed algorithm.

Assume a Bézier surface of degree $5\times 7$ is given by

[TABLE]

where $P_{ij}\in\mathbb{R}^{3}$ are the control points. To evaluate points on the surface we first reformulate the surface in matrix form

[TABLE]

Let $M_{L}=(P_{0,0},\ldots,P_{0,7},\ldots,P_{5,0},\ldots,P_{5,7})$, $\Phi_{B_{n}}(t)=(B_{0,n}(t),\ldots,B_{n,n}(t))^{T}$ and $\Phi(u,v)=\Phi_{B_{5}}(u)\otimes\Phi_{B_{7}}(v)$. The Bézier surface can be further represented as $S(u,v)=M_{L}\Phi(u,v)$. As $M_{L}$ is a $3\times 48$ matrix, using the technique presented in [31], we lift the Bézier surface from $\mathbb{R}^{3}$ to $\mathbb{R}^{48}$. Assume $M_{X}$ is the lifted matrix of order 48 and $X_{ij}$ are the lifted control points, the lifted surface becomes $X(u,v)=M_{X}\Phi(u,v)$. When $X(u,v)$ has been evaluated, the point $S(u,v)$ is obtained by choosing the first three coordinates of $X(u,v)$.

To derive the matrices for dynamic evaluation of curves on the Bézier surface, we first compute the transformation matrices for bases $\Phi_{B_{5}}(u)$ and $\Phi_{B_{7}}(v)$. Assume the translated parameter step is $h$, we choose $a=h$ and $b=1+h$. From Equation (7) we have the transformation matrices $C_{B_{5}}$ or $C_{B_{7}}$ for the basis $\Phi_{B_{5}}(u)$ or $\Phi_{B_{7}}(v)$. Then, the transformation matrices for the basis $\Phi(u,v)$ with translated parameter $u$ or $v$ are obtained as $C_{u}^{h}=C_{B_{5}}\otimes I_{8}$ or $C_{v}^{h}=I_{6}\otimes C_{B_{7}}$, respectively. Now, the matrices for dynamic evaluation of points on $u$-curves or $v$-curves on the lifted surface $X(u,v)$ with a fixed parameter step $h$ are obtained as $M_{u}^{h}=M_{X}C_{u}^{h}M_{X}^{-1}$ and $M_{v}^{h}=M_{X}C_{v}^{h}M_{X}^{-1}$. If we replace the parameter step $h$ by $-h$, we have matrices $M_{u}^{-h}$ and $M_{v}^{-h}$ for dynamic evaluation of $u$-curves or $v$-curves in the opposite directions.

In our experiments, we choose $h=\frac{1}{80}$ and $X_{0,0}$ as the initial point for dynamic evaluation of a piecewise smooth curve that is consisting of 33 pieces of full or partial iso-parameter curves. Particularly, we evaluate points on $u$-curves with parameter step $h$, $v$-curves with parameter step $h$, $u$-curves with parameter step $-h$ and $v$-curves with parameter step $-h$, alternately. Assume the curve segments are numbered as $j=0,1,\ldots,32$. The point number for each curve segment is chosen as $m_{j}=80-5\times[(j-1)/2]$, where $[(j-1)/2]$ means the integer part of a real number. When points on a specified curve segment have been evaluated, the obtained last point is chosen as the start point for dynamic evaluation of next curve segment. See Figure 2(a) for the evaluated points and Figure 2(b) for the obtained piecewise curve. We note that the lastly evaluated point by the proposed technique is corresponding to the center of the surface. Assume the distance between two corner control points $P_{0,0}$ and $P_{5,7}$ is 1. The absolute error between the last point obtained by the proposed algorithm and $S(0.5,0.5)$ computed by conventional de Casteljau algorithm is $7.931\times 10^{-13}$.

Example 4. In the fourth example we evaluate a family of iso-parameter curves on a helicoidal patch. Suppose the surface patch is given by

[TABLE]

Let $\Phi(u,v)=(1,v,\cos v,\sin v,u\cos v,u\sin v)^{T}$. It is easily verified that the space spanned by the basis $\Phi(u,v)$ is closed with respect to partial differentiations $\frac{\partial}{\partial u}$ and $\frac{\partial}{\partial v}$. To evaluate the surface by the dynamic algorithm, we lift the surface from $\mathbb{R}^{3}$ to $\mathbb{R}^{6}$ by adding three more coordinates to the coefficients. The lifted surface represented in matrix form is

[TABLE]

Denote the coefficient matrix of $X(u,v)$ as $M_{X}$. It yields that $X(u,v)=M_{X}\Phi(u,v)$. When $X(u,v)$ has been evaluated, the point $(x(u,v),y(u,v),z(u,v))$ is obtained by choosing the first three coordinates.

Because the basis vector $\Phi(u,v)$ can be decomposed as $\Phi(u,v)=U_{1}(v)\sqcup U_{1}(u)\otimes V(v)$, the transformation matrix for the basis vector $\Phi(u,v)$ with respect to the translation of parameter $v$ is obtained as $C_{v}^{h}=\operatorname{diag}(M_{U_{1}}^{h},I_{2}\otimes M_{V}^{h})$. Because $\det(M_{X})=-4\neq 0$, we compute a transformation matrix as $M_{v}^{h}=M_{X}C_{v}^{h}M_{X}^{-1}$. The points on any surface curve with a fixed parameter $u$ are then computed by

[TABLE]

According to Equation (20), points on a family of $v$-curves on the surface are obtained iteratively starting from a set of points on the boundary line. Figure 3 illustrates the evaluated results after 10, 120 or 200 steps of evaluation, where $u=0,0.2,0.4,\ldots,2$ and the parameter step is chosen as $h=\frac{4\pi}{200}$.

Example 5. Lastly, we evaluate curves with skew parametrization on a Dupin-Cyclide. Let $a=6$, $b=4\sqrt{2}$, $c=2$ and $\mu=3$. The Cartesian coordinates of a Dupin-Cyclide are given by [24]

[TABLE]

where

[TABLE]

To evaluate the Cartesian coordinates of the surface, we should compute the homogeneous coordinates first. Let

[TABLE]

and $\Phi(u,v)=(1,\sin v,\cos v,\sin u,\sin u\sin v,\sin u\cos v,\cos u,\cos u\sin v,\cos u\cos v)^{T}$. The homogeneous coordinates of the surface are represented as $X_{H}(u,v)=M_{H}\Phi(u,v)$. To evaluate the homogeneous coordinates dynamically, we lift $X_{H}(u,v)$ from $\mathbb{R}^{4}$ to $\mathbb{R}^{9}$. Assume $M_{H}=(H_{1},H_{2})$, where $H_{1}$ and $H_{2}$ are the $4\times 5$, $4\times 4$ sub-matrices, respectively. Let

[TABLE]

where $I_{5}$ is the identity matrix of order 5. The lifted homogeneous surface is obtained as $X(u,v)=M_{X}\Phi(u,v)$.

Let $\Theta(t)=(1,\sin t,\cos t)^{T}$. It yields that $\Phi(u,v)=\Theta(u)\otimes\Theta(v)$. It is easily verified that spaces spanned by $\Theta(t)$ or $\Phi(u,v)$ are closed with respect to a differentiation or translation of the parameters. We have $\Theta(t+h)=C_{h}\Theta(t)$ and $\Phi(u+h_{1},v+h_{2})=C_{u,v}^{h_{1},h_{2}}\Phi(u,v)$, where

[TABLE]

and $C_{u,v}^{h_{1},h_{2}}=C_{h_{1}}\otimes C_{h_{2}}$. As the inverse of matrix $M_{X}$ is

[TABLE]

we have $X(u+h_{1},v+h_{2})=M_{u,v}^{h_{1},h_{2}}X(u,v)$, where $M_{u,v}^{h_{1},h_{2}}=M_{X}C_{u,v}^{h_{1},h_{2}}M_{X}^{-1}$. Starting from any point $X(u_{0},v_{0})$, a sequence of points on the surface will be computed by

[TABLE]

In our experiments, we choose $u_{0}=0$, $v_{0}=\pi$ and the start point is obtained as the outer intersection point between the surface and the $x$-axis (the red arrow line) in the positive direction; see Figure 4(a). To evaluate points on piecewise surface curves with skew parametrization, we choose the parameter steps as $h_{1}=\frac{0.18\pi}{m}$ and $h_{2}=\frac{\pi}{m}$, where $m$ is the number of points that will be computed on each piece of surface curve. By choosing $m=100$, two iteration matrices $M_{u,v}^{h_{1},h_{2}}$ and $M_{u,v}^{h_{1},-h_{2}}$ are computed first. Given $X(u_{0},v_{0})$, the points on the first piece of surface curve are computed by Equation (22) using matrix $M_{u,v}^{h_{1},h_{2}}$. Starting from the end of first piece of curve we compute points on the second piece using matrix $M_{u,v}^{h_{1},-h_{2}}$. We continue this process by using $M_{u,v}^{h_{1},h_{2}}$ and $M_{u,v}^{h_{1},-h_{2}}$ for evaluating odd number or even number pieces of curves alternately. Figure 4(a) illustrates the result with 22 pieces of surface curves while Figure 4(b) illustrates the surface curve with 100 pieces. Due to periodicity, the last point on the 100th piece is theoretically the same as the start point of the first one. Practically, the distance between these two points is $6.526\times 10^{-13}$ even after 10000 times of matrix-point multiplication in the presence of truncation errors of irrational numbers.

Numerical stability of dynamic evaluation. Since all dynamically evaluated points on a curve or surface are computed by one or a few constant iteration matrices and a start point, the accuracy of the constant iteration matrices play a key role for the accuracy of the evaluated points. The results computed by the Taylor method in Example 1 show that an inexact iteration matrix may cause deviations in the following evaluated points. All examples employing the basis transformation technique demonstrate that the new method can be used to compute the iteration matrices and the points on curves or surfaces accurately enough.

Even the iteration matrix is accurate, the noise at the initial point can propagate to the dynamically evaluated points. Fortunately, the propagated errors are bounded and controlled when points on a curve segment or a surface patch are dynamically evaluated. Suppose $X(t)=M_{X}\Phi(t)$ is an exponential polynomial curve as defined in Section 3.1 and a sequence of points on the curve are computed by Equation (8). If the start point has been changed as $\tilde{X}(t_{0})=X(t_{0})+X_{\varepsilon}$, the dynamically evaluated points become

[TABLE]

where $M_{h}=M_{X}C_{h}M_{X}^{-1}$ and $C_{h}=e^{Ah}$ are as defined in Equation (8) or Proposition 2.4. The error magnitude for the $i$th point is estimated as

[TABLE]

Since we compute points on a curve segment or a surface patch in practice, the parameter $ih$ lies on a limited interval. Therefore, the norm of the matrix $M_{X}e^{Aih}M_{X}^{-1}$ and the noise magnitudes of the evaluated points are bounded. We have recomputed points for above examples using start points with added noise. It is found that the deviation magnitudes for the dynamically evaluated points are around the same or a few times larger than the magnitudes of added noise.

5 Conclusions

This paper has presented a robust and efficient algorithm for dynamic evaluation of free-form curves and surfaces constructed by general exponential polynomials. By explicit computation of transformation matrices between exponential polynomial bases with or without translation of the parameter, points on curves or surfaces with equal parameter steps can be evaluated dynamically with only arithmetic operations. The proposed technique suffers no shortcomings of classical numerical algorithms for solving linear differential systems any more and it can be used for accurate and stable evaluation of general exponential polynomial curves with any parameter steps. Besides evaluating points with fixed parameter steps or families of iso-parameter curves on surfaces, the basis transformation technique can also be used for evaluating polynomial curves with changing parameter steps or dynamic evaluation of skew-parameterized curves on surfaces in a simple and efficient way.

Acknowledgment

We owe thanks to referees for their invaluable comments and suggestions which helped to improve the presentation of the paper greatly.