# A Modified Complete Spline Interpolation and Exponential Parameterization

Ryszard Kozera, Lyle Noakes, Magdalena Wilkolazka

► To cite this version:

Ryszard Kozera, Lyle Noakes, Magdalena Wilkolazka. A Modified Complete Spline Interpolation and Exponential Parameterization. 14th Computer Information Systems and Industrial Management (CISIM), Sep 2015, Warsaw, Poland. pp.98-110, 10.1007/978-3-319-24369-6_8. hal-01444508

HAL Id: hal-01444508

https://inria.hal.science/hal-01444508

Submitted on 24 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Modified Complete Spline Interpolation and Exponential Parameterization

Ryszard Kozera¹,³, Lyle Noakes², and Magdalena Wilkolazka³

¹ Warsaw University of Life Sciences - SGGW

Faculty of Applied Informatics and Mathematics

Nowoursynowska str. 159, 02-776 Warsaw, Poland

ryszard.kozera@gmail.com ryszard.kozera@sggw.pl

² The University of Western Australia

School of Mathematics and Statistics

35 Stirling Highway, Crawley W.A. 6009, Perth, Australia

lyle.noakes@maths.uwa.edu.au

³ John Paul II Catholic University of Lublin

Faculty of Mathematics, IT and Landscape Architecture

Konstantynów str. 1H, 20-708 Lublin, Poland,

magda.wilkolazka@gmail.com

**Abstract.** In this paper a modified complete spline interpolation based on reduced data is examined in the context of trajectory approximation. Reduced data constitute an ordered collection of interpolation points in arbitrary euclidean space, stripped from the corresponding interpolation knots. The exponential parameterization (controlled by $\lambda \in [0, 1]$) compensates the above loss of information and provides specific scheme to approximate the distribution of the missing knots. This approach is commonly used in computer graphics or computer vision in curve modeling and image segmentation or in biometrics for feature extraction. The numerical verification of asymptotic orders $\alpha(\lambda)$ in trajectory estimation by modified complete spline interpolation is performed here for regular curves sampled more-or-less uniformly with the missing knots parameterized according to exponential parameterization. Our approach is equally applicable to either sparse or dense data. The numerical experiments confirm a slow linear convergence orders $\alpha(\lambda) = 1$ holding for all $\lambda \in [0, 1)$ and a quartic one $\alpha(1) = 4$ once modified complete spline is used. The paper closes with an example of medical image segmentation.

**Keywords:** Spline interpolation, curve approximation and modeling, reduced data, biometrics and feature extraction, computer graphics and vision, medical image processing.

## 1 Problem formulation

Let $\gamma : [0, T] \to E^n$ be a smooth regular parametric curve, i.e. the curve with $\dot{\gamma}(t) \neq 0$ over $t \in [0, T]$ (here $T < \infty$). Reduced data represent a sequence of $m+1$ interpolation points $Q_m = \{q_i\}_{i=0}^m$ in arbitrary euclidean space $E^n$ satisfying
$q_i = \gamma(t_i)$ and $q_{i+1} \neq q_i$. The corresponding interpolation knots $\{t_i\}_{i=0}^m$ fulfilling $t_0 < \dots < t_i < \dots < t_m$ are assumed here to be unknown. Any data fitting scheme $\hat{\gamma}$ based on reduced data $Q_m$ is called non-parametric interpolation. In order to construct $\hat{\gamma}$ explicitly, first the knot estimates $\{\hat{t}_i\}_{i=0}^m \approx \{t_i\}_{i=0}^m$ need to be somehow guessed (here one naturally sets $\hat{\gamma}(\hat{t}_i) = q_i$). Upon selecting a specific interpolation scheme $\hat{\gamma}: [0, \hat{T}] \to E^n$ and the substitutes $\{\hat{t}_i\}_{i=0}^m$ of the missing knots $\{t_i\}_{i=0}^m$, the analysis yielding the order $\alpha$ in $\gamma$ approximation by $\hat{\gamma}$ needs to be carried out (for $m \to \infty$). The appropriate choice of $\{\hat{t}_i\}_{i=0}^m$ should ensure convergence of the interpolant $\hat{\gamma}$ to the unknown curve $\gamma$ with possibly fast order $\alpha$.

We recall now a necessary background information (see e.g. [1]). In fact, reduced data $Q_m = \{\gamma(t_i)\}_{i=0}^m$ are formed from the set of admissible samplings:

**Definition 1.** The interpolation knots $\{t_i\}_{i=0}^m$ are called admissible if they satisfy:

$$ \lim_{m \to \infty} \delta_m \to 0^+, \quad \text{where} \quad \delta_m = \max_{1 \le i \le m} \{t_i - t_{i-1} : i = 1, 2, \dots, m\}. \quad (1) $$

In this paper a special subfamily of admissible samplings i.e. the so-called more-or-less uniform samplings is considered (see also [3]):

**Definition 2.** The sampling $\{t_i\}_{i=0}^m$ is more-or-less uniform if for some constants $0 < K_l \le K_u$ and sufficiently large $m$ the following holds:

$$ \frac{K_l}{m} \le t_i - t_{i-1} \le \frac{K_u}{m}, \quad (2) $$

for all $i = 1, 2, \dots, m$. Alternatively, condition (2) can be replaced by the equivalent inequality $\beta\delta_m \le t_{i+1} - t_i \le \delta_m$ satisfied for some $0 < \beta \le 1$ and sufficiently large $m$.

Recall now the next definition (see e.g. [1] or [2]):

**Definition 3.** Consider a family $\{f_{\delta_m}, \delta_m > 0\}$ of functions $f_{\delta_m}: [0, T] \to E$. We say that $f_{\delta_m}$ is of order $O(\delta_m^\alpha)$ (denoted as $f_{\delta_m} = O(\delta_m^\alpha)$), if there is a constant $K > 0$ such that, for some $\bar{\delta} > 0$ the inequality $|f_{\delta_m}(t)| < K\delta_m^\alpha$ holds for all $\delta_m \in (0, \bar{\delta})$, uniformly over $[0, T]$. In case of vector-valued functions $F_{\delta_m}: [0, T] \to E^n$ by $F_{\delta_m} = O(\delta_m^\alpha)$ it is understood that $\|F_{\delta_m}\| = O(\delta_m^\alpha)$ (here $\|\cdot\|$ denotes a standard euclidean norm).

To examine the asymptotics in trajectory estimation (i.e. the coefficient $\alpha$ from Def. 3) in case of classical parametric interpolation $\hat{\gamma}: [0, T] \to E^n$, where both $q_i = \tilde{\gamma}(t_i)$ and $\{t_i\}_{i=0}^m$ are given, one sets $F_{\delta_m} = \gamma - \tilde{\gamma}$. On the other hand, for non-parametric interpolation a slight adjustment in the last expression for $F_{\delta_m}$ is required (see [1]). Indeed, the latter stems from the fact that both domains $[0, T]$ of $\gamma$ and $[0, \hat{T}]$ of $\hat{\gamma}$ do not generically coincide (here $T = t_m$ and $\hat{T} = \hat{t}_m$). Consequently, for the non-parametric interpolant $\hat{\gamma}$ (given $\{\hat{t}_i\}_{i=0}^m$ are
somehow guessed) a reparameterization $\psi : [0, T] \to [0, \hat{T}]$ is needed so that the asymptotics

$$ (\hat{\gamma} \circ \psi)(t) - \gamma(t) = O(\delta_m^\alpha) \qquad (3) $$

can be examined over a common domain $[0, T]$. Noticeably, the reparametrization issue (given fixed $\hat{\gamma}$ and $\{\hat{t}_i\}_{i=0}^m$) is essential to both theoretical and numerical examination determining an intrinsic asymptotics built in (3). In addition, preferably $\psi$ should be a genuine reparameterization (i.e. $\dot{\psi} > 0$) e.g. if length of $\gamma$ is to be estimated by the length of $\hat{\gamma}$. Evidently, on its own the construction of the interpolant $\hat{\gamma}$ does rely on explicit formula standing for $\psi$. Independently from $\psi$, the derivation of any non-parametric interpolant $\hat{\gamma}$ stipulates the appropriate choice of the estimates $\{\hat{t}_i\}_{i=0}^m$ mimicking the missing knots $\{t_i\}_{i=0}^m$. In doing so, recall now a definition of exponential parameterization (e.g. in [4]):

$$ \hat{t}_0^\lambda = 0 \quad \text{and} \quad \hat{t}_j^\lambda = \hat{t}_{j-1}^\lambda + \|q_j - q_{j-1}\|^\lambda, \qquad (4) $$

where $j = 1, 2, \dots, m$ and $\lambda \in [0, 1]$. If $\lambda = 0$ a blind guess yielding uniform knots $\hat{t}_i^0 = i$ follows. On the other hand, the case of $\lambda = 1$ results in a cumulative chord parameterization $\hat{t}_i^1 = \hat{t}_{j-1}^1 + \|q_j - q_{j-1}\|$ (see [4] or [5]). From now on we suppress the superscript notation with $\lambda$ in (4), unless needed otherwise. The term exponential parameterization stands for the determination of a discrete set of knots $\{\hat{t}_i\}_{i=0}^m \approx \{t_i\}_{i=0}^m$, whereas a similar term i.e. a reparameterization represents a piecewise-smooth mapping $\psi : [0, T] \to [0, \hat{T}]$.

Previous result [6] proved that for $\lambda = 1$ and for an arbitrary admissible sampling (1) a Lagrange piecewise-quadratic(-cubic) $\hat{\gamma}_r$ ($r = 2, 3$ - see e.g. [2]) interpolation combined with (4) yields $\alpha_r(1) = r + 1$. Hence for cumulative chords the interpolant $\hat{\gamma}_r$ ($r = 2, 3$) renders either cubic or quartic convergence orders in trajectory estimation (see (3)). Interestingly, opposite to the parametric interpolation $\hat{\gamma}_r$, the convergence orders in question do not necessarily increase for $r > 3$ and $\lambda = 1$ - see [1] or [7]. In addition, a recent result by [8] (see also [9]) proves that $\hat{\gamma}_2$ combined with (4) and more-or-less uniformly sampled (2) reduced data $Q_m$ yields $\alpha(1) = 3$ and $\alpha(\lambda) = 1$ for $\lambda \in [0, 1)$. The latter demonstrates an unexpected left-hand side discontinuity of $\alpha(\lambda)$ at $\lambda = 1$. Interestingly, such trend continues once (4) is combined with piecewise-cubics $\hat{\gamma}_3$. Indeed the following holds (see [10]):

**Theorem 1.** Suppose $\gamma$ is a regular $C^4([0, T])$ curve in $E^n$ sampled more-or-less uniformly (2). Assume that $\{\hat{t}_i^\lambda\}_{i=0}^m$ are computed from $Q_m$ according to (4). Then there exists a piecewise-cubic $C^\infty$ mapping $\psi : [0, T] \to [0, \hat{T}]$, such that over $[0, T]$, we have for either $\lambda \in [0, 1)$:

$$ \hat{\gamma}_3 \circ \psi - \gamma = O(\delta_m) \qquad (5) $$

or for $\lambda = 1$ (and (1)):

$$ \hat{\gamma}_3 \circ \psi - \gamma = O(\delta_m^4). \qquad (6) $$

Undesirably, the interpolants $\hat{\gamma}_r$ ($r = 2, 3$), are generically non-smooth at junction points, where both neighboring local quadratics (cubics) are glued together over two consecutive segments $[t_i, t_{i+r}]$ and $[t_{i+r}, t_{i+2r}]$ (with $r = 2, 3$). In
order to alleviate such deficiency, a modified $C^1$ Hermite interpolation $\hat{\gamma}_3^H$ based on $Q_m$, cumulative chords and general admissible samplings (1) is introduced and examined in [11] or [12]. Here the unknown derivatives at all interpolation points $\{q_i\}_{i=0}^m$ are approximated with high accuracy via special procedure (see [11]). This permits to obtain quartic order $\alpha(1) = 4$ in trajectory estimation once $\hat{\gamma}_3^H$ and (1) are coupled together. Analogously to Th. 1, the latter extends to all remaining $\lambda \in [0, 1)$ (for samplings (2)) resulting in $\alpha(\lambda) = 1$ (see [13]). Recurrent left-hand side discontinuity in convergence order $\alpha(\lambda)$ at $\lambda = 1$ is here manifested again.

For certain applications (e.g. approximation of curvature of $\gamma$, image segmentation or other feature extraction in biometrics) the interpolant $\hat{\gamma}$ should be at least continuously twice differentiable. Such constraint is not generically fulfilled by so-far discussed interpolants at any junction point. The remedy guaranteeing $C^2$ smoothness is met upon applying various hybrids of $C^2$ cubic spline interpolants $\hat{\gamma}_3^S$ (see [2]) based on $Q_m$ and (4). One of them (called complete cubic spline $\hat{\gamma}_3^C$) relies on the provision of initial and terminal velocities $\gamma'(t_0 = 0) = v_0$ and $\gamma'(t_m = T) = v_m$ usually not accompanying reduced data $Q_m$. This special case is discussed in [14] (also limited exclusively to $\lambda = 1$), where quartic order $\alpha(1) = 4$ for trajectory estimation by $\hat{\gamma}_3^C$ is established.

In this paper we extend the latter (at least with the aid of numerical tests) to twofold more general situation. Similarly to $\hat{\gamma}_3^H$, we estimate first both missing velocities $\gamma'(t_0) \approx v_0^a$ and $\gamma'(t_m) \approx v_m^a$. Next a modified complete spline interpolant $\hat{\gamma}_3^C$ based on $Q_m$, $v_0^a$, $v_m^a$ and (4) is introduced for all $\lambda \in [0, 1]$ - see Section 2. The conjectured asymptotics reads as:

**Theorem 2.** Let $\gamma$ be a regular $C^4([0, T])$ curve in $E^n$ sampled more-or-less-uniformly (4). Approximate $(\gamma'(t_0), \gamma'(T))$ with $v_0^a = \hat{\gamma}_3'(0)$ and $v_m^a = \hat{\gamma}_3'(\bar{T})$, where $\hat{\gamma}_3$ defines a piecewise cubic based on $Q_m$ and (4) with $\lambda \in [0, 1]$. Assume also that $\hat{\gamma}_3^C : [0, \bar{T}] \to E^n$ define a modified complete spline constructed on $Q_m$, estimated velocities $(v_0^a, v_m^a)$ and exponential parameterization (4). Then there is a piecewise-$C^\infty$ mapping $\psi : [0, T] \to [0, \bar{T}]$ such that over $[0, T]$ we either have for all $\lambda \in [0, 1)$:

$$ \hat{\gamma}_3^C \circ \psi - \gamma = O(\delta_m) \qquad (7) $$

or for $\lambda = 1$:

$$ \hat{\gamma}_3^C \circ \psi - \gamma = O(\delta_m^4). \qquad (8) $$

In Section 3 the asymptotics from Th. 2 is numerically verified as sharp and specific application of modified $C^2$ complete spline $\hat{\gamma}_3^C$ is given. In addition, we compare our interpolant $\hat{\gamma}_3^C$ against $\hat{\gamma}_3^H$. Finally, our paper concludes with hints for possible extension of this work. Extra literature references concerning related work and spin-off applications are also provided.

## 2 Modified Complete Spline on Reduced Data

A modified complete spline interpolant $\hat{\gamma}_3^C$ based on reduced data $Q_m$ (see also [2]) and exponential parameterization (4) is introduced below. This scheme applicable to both dense and sparse $Q_m$ falls into the following steps:
1. Calculate the estimates {$\hat{t}_i$}$_{i=0}^m$ of the missing knots {$t_i$}$_{i=0}^m$ according to the exponential parameterization (4) (with $\lambda \in [0, 1]$).

2. The so-called general $C^2$ piecewise-cubic spline $\hat{\gamma}_{3,i}^S$ interpolant (a sum-track of cubics {$\hat{\gamma}_{3,i}^S$}$_{i=0}^{m-1}$ - see [2]) fulfills the following constraints over each segment $[\hat{t}_i, \hat{t}_{i+1}]$:

$$ 
\begin{aligned}
\hat{\gamma}_{3,i}^{S}(\hat{t}_i) &= q_i, & \hat{\gamma}_{3,i}^{S}(\hat{t}_{i+1}) &= q_{i+1}, \\
\hat{\gamma}_{3,i}^{S'}(\hat{t}_i) &= \mathbf{v}_i, & \hat{\gamma}_{3,i}^{S'}(\hat{t}_{i+1}) &= \mathbf{v}_{i+1},
\end{aligned}
\qquad (9) $$

where $\mathbf{v}_0, \dots, \mathbf{v}_m$ represent the unknown slopes $\mathbf{v}_i \in \mathbb{R}^n$. The internal velocities {$\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_{m-1}$} can be uniquely computed from $C^2$ constraints imposed on $\hat{\gamma}_3^S$ at junction points {$q_1, \dots, q_{m-1}$} i.e. by enforcing:

$$ \hat{\gamma}_{3,i-1}^{S''}(\hat{t}_i) = \hat{\gamma}_{3,i}^{S''}(\hat{t}_i), \qquad (10) $$

provided both $\mathbf{v}_0$ and $\mathbf{v}_m$ are somehow computed (or a priori given). The computational method to determine all slopes {$\mathbf{v}_i$}$_{i=0}^m$ (including initial and terminal ones) is discussed in next.

3. Assuming temporarily the provision of all velocities {$\mathbf{v}_i$}$_{i=0}^m$, each cubic $\hat{\gamma}_{3,i}^S$ over $\hat{t} \in [\hat{t}_i, \hat{t}_{i+1}]$ reads as:

$$ \hat{\gamma}_{3,i}^{S}(\hat{t}) = c_{1,i} + c_{2,i}(\hat{t} - \hat{t}_i) + c_{3,i}(\hat{t} - \hat{t}_i)^2 + c_{4,i}(\hat{t} - \hat{t}_i)^3, \qquad (11) $$

where its respective coefficients (with $\Delta\hat{t}_i = \hat{t}_{i+1} - \hat{t}_i$) are equal to:

$$ c_{1,i} = q_i, \quad c_{2,i} = \mathbf{v}_i, $$

$$ c_{3,i} = \frac{\frac{q_{i+1}-q_i}{\Delta\hat{t}_i} - \mathbf{v}_i}{\Delta\hat{t}_i} - c_{4,i}\Delta\hat{t}_i, \quad c_{4,i} = \frac{\mathbf{v}_i + \mathbf{v}_{i+1} - 2\frac{q_{i+1}-q_i}{\Delta\hat{t}_i}}{(\Delta\hat{t}_i)^2}. \qquad (12) $$

If additionally $\mathbf{v}_i = \gamma'(t_i)$ are given then formulas (11) and (12) yield a well-known $C^1$ Hermite spline. However, the required velocities {$\mathbf{v}_0, \mathbf{v}_1, \dots, \mathbf{v}_m$} are not usually supplemented to $Q_m$. A scheme for computing the corresponding missing internal velocities {$\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_{m-1}$} is recalled next (see [2]). Following the latter a method of estimating {$\mathbf{v}_0, \mathbf{v}_m$} is given. It is inspired by the approach adopted in [11].

4. Formulas (11) and (12) render $\hat{\gamma}_{3,i}^{S''}(\hat{t}_i) = 2c_{3,i}$ and $\hat{\gamma}_{3,i-1}^{S''}(\hat{t}_i) = 2c_{3,i-1} + 6c_{4,i-1}(\hat{t}_i - \hat{t}_{i-1})$ which combined with (10) leads to the linear system:

$$ \mathbf{v}_{i-1}\Delta\hat{t}_i + 2\mathbf{v}_i(\Delta\hat{t}_{i-1} + \Delta\hat{t}_i) + \mathbf{v}_{i+1}\Delta\hat{t}_{i-1} = b_i, \qquad (13) $$

where

$$ b_i = 3 \left( \Delta \hat{t}_i \frac{q_i - q_{i-1}}{\Delta \hat{t}_{i-1}} + \Delta \hat{t}_{i-1} \frac{q_{i+1} - q_i}{\Delta \hat{t}_i} \right). $$

Assuming that the end-slopes $\mathbf{v}_0$ and $\mathbf{v}_m$ are somehow given the system (13) solves uniquely in {$\mathbf{v}_i$}$_{i=1}^{m-1}$. The latter yields a $C^2$ spline $\hat{\gamma}_3^S$ (which fits
reduced data $Q_m$) defined as a track-sum of ${\hat{\gamma}_{3,i}^S}_{i=0}^{m-1}$ introduced in (11).
If extra conditions hold, i.e. $\gamma'(t_0) = v_0$ and $\gamma'(T) = v_m$ then $\hat{\gamma}_3^S$ is called a
complete cubic spline (denoted here as $\hat{\gamma}_3^C$).

5. Since $Q_m$ are usually deprived from both initial and terminal velocities $\{\gamma'(t_0) = v_0, \gamma'(T) = v_m\}$ a good estimate $\{v_0^a, v_m^a\}$ is therefore required. Of course, any choice of $\{v_0^a, v_m^a\}$ renders a unique explicit formula for $\hat{\gamma}_3^C$. This however is insufficient for our consideration. Indeed, still a proper estimate of these two velocities is needed so that (7) and (8) follow. In doing so, we invoke Lagrange cubic $\hat{\gamma}_{3,0}^L : [0, \hat{t}_4^\lambda] \to E^n$ (and $\hat{\gamma}_{3,m-3}^L : [\hat{t}_{m-3}^\lambda, \hat{T}] \to E^n$), satisfying $\hat{\gamma}_{3,0}^L(\hat{t}_i^\lambda) = q_i$ (and $\hat{\gamma}_{3,m-3}^L(\hat{t}_{m-3+i}^\lambda) = q_{m-3+i}$), with $i = 0, 1, 2, 3$ - here the same $\lambda \in [0, 1]$ is applied in the derivation of $\hat{\gamma}_{3,0}^L, \hat{\gamma}_{3,m-3}^L$ and $\hat{\gamma}_3^C$. Set now for $v_0^a = \hat{\gamma}_{3,0}^{L'}(0)$ and for $v_m^a = \hat{\gamma}_{3,m-3}^{L'}(\hat{T})$, respectively.

This completes a description of a modified $C^2$ complete spline $\hat{\gamma}_3^C$ based on
reduced data $Q_m$ and exponential parameterization (4).

However, to verify the asymptotics from (7) and (8) (either numerically or theoretically) a candidate for a reparameterization $\psi : [0, T] \to [0, \hat{T}]$ is still required, as justified in Section 1. In doing so, consider a $C^2$ complete spline $\psi = \psi_3^C : [0, T] \to [0, \hat{T}]$ satisfying the knots' interpolation constraints $\psi_3^C(t_i) = \hat{t}_i$, where $\{\hat{t}_i\}_{i=0}^m$ are defined according to (4). In addition, the initial and terminal velocities of $s_0 = \psi_3^{C'}(0)$ and $s_m = \psi_3^{C'}(T)$ are set similarly to the construction from above. More specifically, define two Lagrange cubics $\psi_{3,0} : [0, t_{i+3}] \to [0, \hat{t}_{i+3}^\lambda]$ and $\psi_{3,m-3} : [t_{m-3+i}, T] \to [\hat{t}_{m-3+i}^\lambda, \hat{T}]$ satisfying interpolation conditions $\psi_{3,0}(t_i) = \hat{t}_i^\lambda$ and $\psi_{3,m-3}(t_{m-3+i}) = \hat{t}_{m-3+i}^\lambda$ (with $i = 0, 1, 2, 3$ and the same $\lambda \in [0, 1]$ as for the construction of $\hat{\gamma}_3^C$), respectively. One sets here for $s_0 = \psi_3^{C'}(0) = \psi'_{3,0}(0)$ and for $s_m = \psi_3^{C'}(T) = \psi'_{3,m-3}(T)$.

We pass now to the experimental section of this paper which tests the asymptotics from Th. 2. As already indicated, a sole derivation of a modified $C^2$ complete spline $\hat{\gamma}_3^C$ relies exclusively on reduced data $Q_m$ (either dense or sparse) and (4). On the other hand, any numerical verification or theoretical proof of the asymptotics $\alpha(\lambda)$ involved (e.g. from Th. 2), requires an extra introduction of reparameterization $\psi$ (proposed here as $\psi_3^C$) as well as an admittance of sufficiently densely more-or-less uniformly sampled points $Q_m$. The latter enables to assess a desired asymptotics controlling the decrease in difference $\hat{\gamma}_3 \circ \psi_3^C - \gamma$, uniformly over $[0, T]$ (once $m \to \infty$).

3 Experiments

In this section, a numerical verification of the asymptotics $\alpha(\lambda)$ (and its sharpness) claimed in Th. 2 is conducted. Recall that, given fixed $\lambda \in [0, 1]$, by sharpness we understand the existence of at least one curve $\gamma \in C^4(0, T]$) and one special family of more-or-less uniform sampling (2) such that the asymptotics in differences $\hat{\gamma}_3^C \circ \psi_3^C - \gamma$ (over $[0, T]$) is not faster than predicted $\alpha(\lambda)$. A positive verification of (7) and (8) would point out again to a bizarre phenomenon. Namely, the existence of the left-hand side discontinuity in $\alpha(\lambda)$ at $\lambda = 1$.
All tests for this paper are carried out in Mathematica 8.0⁴ (see also [15]) and resort to two types of skew-symmetric more-or-less uniform samplings (2). The first one selected (for $t_i \in [0, 1]$) reads as:

$$ t_i = \begin{cases} \frac{i}{m} + \frac{1}{2m}, & \text{for } i = 4k + 1; \\ \frac{i}{m} - \frac{1}{2m}, & \text{for } i = 4k + 3; \\ \frac{i}{m}, & \text{for } i \text{ even;} \end{cases} \qquad (14) $$

with $K_l = (1/2)$ and $K_u = 1$ as introduced in (2). The second one is defined according to:

$$ t_i = \frac{i}{m} + \frac{(-1)^{i+1}}{3m}, \qquad (15) $$

with constants $K_l = (1/2)$ and $K_u = (5/3)$ from (2). For a given $m$, the error $E_m$, between $\gamma$ and reparameterized modified complete spline $\hat{\gamma}_3^C \circ \psi_3^C$ reads as:

$$ E_m = \max_{t \in [0,1]} \|(\hat{\gamma}_3^C \circ \psi_3^C)(t) - \gamma(t)\|. \qquad (16) $$

The latter is computed over each sub-interval $[t_i, t_{i+1}]$ (for $i = 0, \dots, m-1$) by using Mathematica function - FindMaximum and then upon taking the maximal values from all segments' optima. In order to approximate $\alpha(\lambda)$ we calculate first $E_m$ for $m_{min} \le m \le m_{max}$, where $m_{min}$ and $m_{max}$ are sufficiently large fixed constants. Then a linear regression yielding a function $y(x) = \bar{\alpha}(\lambda)x + b$ is applied to $\{(\log(m), -\log(E_m))\}_{m_{min}}^{m_{max}}$. Mathematica built-in function Linear-ModelFit extracts a coefficient $\bar{\alpha}(\lambda) \approx \alpha(\lambda)$. A full justification of this procedure to approximate $\alpha(\lambda)$ by $\bar{\alpha}(\lambda)$ is given in [1]. Note also that since both (7) and (8) have asymptotic character the constants $m_{min} < m_{max}$ should be taken as sufficiently large. On other hand, a potential negative impact of machine rounding-off errors stipulates these two constants not to exceed big values. In practice, the appropriate choices for $m_{min} < m_{max}$ are adjusted each time during the experimental phase. The tests conducted here employ three types of $C^\infty$ regular curves: a spiral $\gamma_{sp}$ and a cubic $\gamma_c$ both in $E^2$ as well as a helix $\gamma_h$ in $E^3$. They are sampled more-or-less uniformly according to either (14) or (15). For comparison reasons we also test here the asymptotic orders $\alpha_H(\lambda)$ in trajectory estimation for modified $C^1$ Hermite interpolant $\hat{\gamma}_3^H$ examined in [11] and [12] (here $\alpha_H(1) = 4$ and $\alpha_H(\lambda) = 1$ for $\lambda \in [0, 1)$). However, since the interpolant $\hat{\gamma}_3^H \in C^1$ (over $Q_m$) it does not permit to approximate the curvature of $\gamma$ at interpolation points. However, the latter can be accomplished with the aid of $\hat{\gamma}_3^C$ due it is higher order of smoothness (i.e. $\hat{\gamma}_3^C \in C^2$ over $\hat{t} \in [0, \hat{T}]$).

**Example 1.** Consider a regular planar spiral $\gamma_{sp} : [0, 1] \to E^2$,

$$ \gamma_{sp}(t) = ((0.2 + t) \cos(\pi(1 - t)), (0.2 + t) \sin(\pi(1 - t))). \qquad (17) $$

Figure 1 (or Figure 2) contains the plots of $\gamma_{sp}$ (or of $\hat{\gamma}_3^C$) with $\lambda = 0$ sampled (here $m = 15$) according to either (14) or (15).
<figure><img src="image_2.png" /><figcaption>Fig. 1. A spiral $\gamma_{sp}$ from (17) sampled along (dotted): a) (14) or b) (15), for $m = 15$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. A spiral $\gamma_{sp}$ from (17) for: a) (14) b) (15) fitted by $\hat{\gamma}_3^C$ (here $m = 15$ and $\lambda = 0$).</figcaption></figure>

The respective linear regression based estimates $\bar{\alpha}(\lambda) \approx \alpha(\lambda)$ (for various $\lambda \in [0, 1]$) are computed here for $m_{min} = 60 \le m \le m_{max} = 120$. The numerical results contained in Table 1 confirm the sharpness of (7) and (8) for $\lambda \in \{0.0, 0.1, 0.3, 0.5, 0.7\}$ and yield marginally faster (though still consistent with asymptotics from Th. (2)) $\alpha(\lambda)$ for $\lambda \approx 1$. For comparison reasons, Table 1 contains also the corresponding numerical results established for estimating $\gamma$ with modified $C^1$ Hermite interpolant $\hat{\gamma}_3^H$ based on the same reduced data $Q_m$ and exponential parameterization.

<table><caption>Table 1. Computed $\bar{\alpha}(\lambda) \approx \alpha(\lambda)$ in (7) & (8) for $\gamma_{sp}$ from (17) and various $\lambda \in [0, 1]$.</caption><thead><tr><th>λ</th><th>0.0</th><th>0.1</th><th>0.3</th><th>0.5</th><th>0.7</th><th>0.9</th><th>1.0</th></tr></thead><tbody><tr><td>α(λ) for (14)</td><td>1.0067</td><td>1.0085</td><td>1.0134</td><td>1.0218</td><td>1.0409</td><td>1.1463</td><td>4.2537</td></tr><tr><td>α(λ) for (15)</td><td>1.0121</td><td>1.0128</td><td>1.0160</td><td>1.0248</td><td>1.0506</td><td>1.2099</td><td>3.9912</td></tr><tr><td>α(λ) in Th. 2</td><td>1.0</td><td>1.0</td><td>1.0</td><td>1.0</td><td>1.0</td><td>1.0</td><td>4.0</td></tr><tr><td>α<sub>H</sub>(λ) for (14)</td><td>1.0070</td><td>1.0084</td><td>1.0129</td><td>1.0205</td><td>1.0371</td><td>1.1282</td><td>3.9192</td></tr><tr><td>α<sub>H</sub>(λ) for (15)</td><td>1.0009</td><td>1.0023</td><td>1.0113</td><td>1.0484</td><td>1.0499</td><td>4.8304</td><td>4.0584</td></tr></tbody></table>

We pass now to the example with a helix having a trajectory in $E^3$.

*Example 2.* Let $\gamma_h : [0, 1] \to E^3$ be defined as

$$ \gamma_h(t) = (1.5 \cos(2\pi t), \sin(2\pi t), 2\pi t/4). \quad (18) $$

<footnote n="4">⁴ This research was supported in part by computing resources of ACC Cyfronet AGH.</footnote>
Figure 3 (or Figure 4) illustrates the trajectories of $\gamma_h$ (or of $\hat{\gamma}_3^C$) for $\lambda = 0.3$ sampled according to either (14) or (15), with $m = 15$. As previously, a linear regression estimating $\bar{\alpha}(\lambda) \approx \alpha(\lambda)$ from Th. 2 is used here, for $m$ ranging over $60 \le m \le 120$ with various $\lambda \in [0, 1]$. The coefficients $\bar{\alpha}(\lambda)$ (see Table 2)

<figure><img src="image_4.png" /><figcaption>Fig. 3. A helix from (18) sampled along (dotted): a) (14) b) (15), for $m = 15$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. A helix $\gamma_h$ from (18) for: a) (14) b) (15) fitted by $\hat{\gamma}_3^C$ (here $m = 15$ and $\lambda = 0.3$).</figcaption></figure>

computed numerically all sharply coincide with those specified in (7) and (8). Again, for comparison reasons, Table 2 presents the corresponding numerical results derived for estimating $\gamma$ with modified $C^1$ Hermite interpolant $\hat{\gamma}_3^H$ based on the same reduced data $Q_m$ and exponential parameterization.

Finally, a planar cubic $\gamma_c$ is tested.

*Example 3.* Let $\gamma_c : [0, 1] \to E^2$ be defined as follows:

$$ \gamma_c(t) = (\pi t, (\pi t + 1)^3(\pi + 1)^{-3}). \qquad (19) $$

Figure 5 (or Figure 6) contains the plots of $\gamma_c$ (or of $\hat{\gamma}_3^C$) sampled along either (14) or (15), with $\lambda = 1$ and $m = 15$. In order to compute $\bar{\alpha}(\lambda) \approx \alpha(\lambda)$ estimating the asymptotics from Th. 2 again a linear regression is used (as explained at the beginning of this section) for $60 \le m \le 120$ and varying $\lambda \in [0, 1]$. Table
<table><caption>**Table 2.** Computed $\bar{\alpha}(\lambda) \approx \alpha(\lambda)$ in (7) & (8) for $\gamma_h$ from (18) and various $\lambda \in [0, 1]$.</caption><thead><tr><th>λ</th><th>0.0</th><th>0.1</th><th>0.3</th><th>0.5</th><th>0.7</th><th>0.9</th><th>1.0</th></tr></thead><tbody><tr><td>ᾱ(λ) for (14)</td><td>1.0000</td><td>1.0000</td><td>1.0002</td><td>1.0006</td><td>1.0009</td><td>1.0065</td><td>3.9949</td></tr><tr><td>ᾱ(λ) for (15)</td><td>1.0000</td><td>1.0000</td><td>1.0001</td><td>1.0005</td><td>1.0015</td><td>1.0127</td><td>3.9992</td></tr><tr><td>α(λ) in Th. 2</td><td>1.0</td><td>1.0</td><td>1.0</td><td>1.0</td><td>1.0</td><td>1.0</td><td>4.0</td></tr><tr><td>ᾱ<sub>H</sub>(λ) for (14)</td><td>1.0000</td><td>1.0000</td><td>1.0002</td><td>1.0003</td><td>1.0008</td><td>1.0049</td><td>3.9833</td></tr><tr><td>ᾱ<sub>H</sub>(λ) for (15)</td><td>1.0002</td><td>1.0001</td><td>1.0001</td><td>1.0026</td><td>1.0043</td><td>1.0317</td><td>3.9888</td></tr></tbody></table>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** A cubic planar curve (19) sampled along (dotted): a) (14) b) (15), for $m = 15$.</figcaption></figure>

3 enlists numerically computed estimates $\bar{\alpha}(\lambda) \approx \alpha(\lambda)$ for various $\lambda \in [0, 1]$ and samplings (14) and (15). Evidently these numerical results re-emphasize the sharpness of the asymptotics determined by (7) and (8), with marginally faster case of $\alpha(1)$. Similarly to the previous examples, Table 3 contains also the corresponding numerical results obtained for estimating $\gamma$ with modified $C^1$ Hermite interpolant $\hat{\gamma}_3^H$ based on the same reduced data $Q_m$ and exponential parameterization.

<table><caption>**Table 3.** Computed $\bar{\alpha}(\lambda) \approx \alpha(\lambda)$ in (7) & (8) for $\gamma_c$ from (19) and various $\lambda \in [0, 1]$.</caption><thead><tr><th>λ</th><th>0.0</th><th>0.1</th><th>0.3</th><th>0.5</th><th>0.7</th><th>0.9</th><th>1.0</th></tr></thead><tbody><tr><td>ᾱ(λ) for (14)</td><td>1.0001</td><td>1.0001</td><td>1.0001</td><td>1.0002</td><td>1.0003</td><td>1.0011</td><td>4.1612</td></tr><tr><td>ᾱ(λ) for (15)</td><td>1.0001</td><td>1.0001</td><td>1.0002</td><td>1.0002</td><td>1.0003</td><td>1.0017</td><td>4.1196</td></tr><tr><td>α(λ) in Th. 2</td><td>1.0</td><td>1.0</td><td>1.0</td><td>1.0</td><td>1.0</td><td>1.0</td><td>4.0</td></tr><tr><td>ᾱ<sub>H</sub>(λ) for (14)</td><td>1.0001</td><td>1.0001</td><td>1.0001</td><td>1.0002</td><td>1.0003</td><td>1.0010</td><td>4.2868</td></tr><tr><td>ᾱ<sub>H</sub>(λ) for (15)</td><td>0.9999</td><td>1.0000</td><td>1.0001</td><td>1.0002</td><td>0.9998</td><td>0.9991</td><td>4.3044</td></tr></tbody></table>

The examples presented herein demonstrate the sharpness of (7) and (8) resulting in a left-hand side discontinuity of $\alpha(\lambda)$ at $\lambda = 1$ which is consistent with Th. 2. We close this section with an application of $\hat{\gamma}_3^C$ to medical image processing.

*Example 4.* A medical image of a kidney is shown in Figure 7. A segmentation of an image of any human organ from its image background (e.g. from a digital image) permits to focus on vital geometrical or other properties (like $\gamma$ perimeter, section internal area, average curvature) of the examined organ. This ultimately
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** A cubic from $\gamma_c$ (19) for: a) (14) b) (15) fitted by $\hat{\gamma}_3^C$ (here $m = 15$ and $\lambda = 1$).</figcaption></figure>

can be exploited in medical diagnosis and further treatment. Indeed, a physician
can mark $m+1$ selected consecutive points on the kidney’s boundary (represent-
ing the trajectory of the unknown curve $\gamma$). Such input points, positioned along
trajectory of $\gamma$, form the set of available interpolation points $Q_m$. Naturally, the
corresponding knots $\{t_i\}_{i=0}^m$ parameterizing $Q_m$ are here defaulted. A modified
complete spline $\hat{\gamma}_3^C$ based on (4) and $Q_m$ can be applied now. The relevant points’
coordinates are determined here by using *Get Coordinate Tool in Mathematica*.
Figure 7 contains of a plot a modified complete spline $\hat{\gamma}_3^C$ based on 67 marked
points (here as $q_0 = q_{67}$ we have 67 different points) with either $\lambda = 0$ or $\lambda = 1$
set in (4) - see Figure 7 a) or b), respectively. Note that the boundary of the
kidney forms a loop which re-translates e.g. into $q_0 = \gamma(0) = \gamma(T) = q_m$. Con-
sequently the interpolant $\hat{\gamma}_3^C$ is generically not smooth at a single point $q_0 = q_m$
unless $v_0 = v_m$, for which $C^1$ class follows. This weakness can be removed e.g.
by taking the average of $v_0$ and $v_m$ at both overlapping “ends” of the curve $\gamma$.
Finally, for comparison reasons, Fig. 8 a) or b) presents the trajectory of the cor-
responding modified Hermite interpolant $\hat{\gamma}_3^H$ constructed on the same reduced
data $Q_m$ and exponential parameterization.

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** The shape of a kidney determined by $\hat{\gamma}_3^C$ with a) $\lambda = 0$ b) $\lambda = 1$, for $m = 67$.</figcaption></figure>

# 4 Conclusion

The tests in Section 3 confirm *the sharpness* of the asymptotics from Th. 2
to approximate $\gamma$ via modified complete spline $\hat{\gamma}_3^C$ based on reduced data $Q_m$,
<figure><img src="image_10.png" /><figcaption>Fig. 8. The shape of a kidney determined by $\hat{\gamma}_3^H$ with a) $\lambda = 0$ b) $\lambda = 1$, for $m = 67$.</figcaption></figure>

more-or-less uniform samplings (2) and exponential parameterization (4). A possible extension of this work includes e.g. an analytical proof of Th. 2 (including investigation of asymptotic constants) or determination of sufficient conditions imposed on samplings $\{t_i\}_{i=0}^m$ to render $\psi_3^C$ as a genuine piecewise-$C^\infty$ reparameterization of $[0, T]$ into $[0, \hat{T}]$. The investigation of the asymptotics in curvature estimation by $\hat{\gamma}_3^C$ is also an open problem. The case with $\lambda = 1$ offers the fastest quartic asymptotics in trajectory estimation for $\hat{\gamma}_3^C$ and (4). However, one can also focus on enforcing specific geometrical properties or constraints by selecting the best $\hat{\gamma}_3^C$ (depending on $\lambda \in [0, 1]$ from (4)) to optimize newly adopted criteria (criterion). This paper shows that if the speed of $\gamma$ approximation is not the main issue then a decisive factor in choosing optimal $\hat{\gamma}_3^C$ should stem from such extra requirement(s) as almost all $\hat{\gamma}_3^C$ have an identical $\alpha(\lambda)$ for $\gamma$ approximation. Related work on $\varepsilon$-uniform samplings combined with (4) can be found in [16]. More specific applications on interpolating (or approximating) reduced data are provided e.g. in [4], [17], [18], [19] or [20]. Splines can also be used in trajectory planning [21], finding algebraic and implicit curves [22] and [23] or in bifurcating surfaces [24]. To supplement (4), there are also other parameterizations applied predominantly on sparse data (applicable also on dense $Q_m$) - see e.g. the so-called *blending parameterization* [25], *monotonicity or convexity preserving ones* [4] or [26].

## References

1. Kozera, R.: Curve modeling via interpolation based on multidimensional reduced data. Studia Informatica 25(4B-61), 1-140 (2004).

2. de Boor, C.: A Practical Guide to Spline. Springer-Verlag, New York Heidelberg, Berlin, (1985).

3. Noakes, L., Kozera, R.: More-or-less uniform samplings and length of curves. Quarterly of Applied Mathematics 61(3), 475-484 (2003).

4. Kvasov, B.I.: Methods of Shape-Preserving Spline Approximation. World Scientific, Singapore, (2000).

5. Lee, E.T.Y.: Choosing nodes in parametric curve interpolation. Computer-Aided Design 21(6), 363-370 (1989).

6. Noakes, L., Kozera, R.: Cumulative chords piecewise-quadratics and piecewise-cubics. In Geometric Properties of Incomplete Data, Eds. R. Klette, R. Kozera,
L. Noakes and J. Weickert, Computational Imaging and Vision Vol. 31, 59–75, Klu-
ver Academic Publishers, The Netherlands (2006)

7. Kozera, R.: Asymptotics for length and trajectory from cumulative chord piecewise-quartics. Fundamenta Informaticae 61(3-4), 267–283 (2004).

8. Kozera, R., Noakes, L.: Piecewise-quadratics and exponential parameterization for reduced data. Applied Mathematics and Computation 221, 620–638 (2013).

9. Kozera, R., Noakes, L., Szmielew, P.: Trajectory estimation for exponential pa-
rameterization and different samplings. Lecture Notes in Computer Science 8104,
430–441, Springer-Verlag (2013).

10. Kozera, R., Noakes, L., Wilkołazka, M.: Piecewise-cubics and exponential param-
eterization for reduced data. Submitted.

11. Kozera, R., Noakes, L.: $C^1$ interpolation with cumulative chord cubics. Funda-
menta Informaticae 61(3-4), 285–301 (2004).

12. Kozera, R., Noakes, L. Szmielew, P.: Quartic orders and sharpness in trajectory estimation for smooth cumulative chord cubics. Lecture Notes in Computer Science 8671, 9–16, Springer (2014).

13. Kozera, R., Noakes, L., Wilkołazka, M.: $C^1$ interpolation with cumulative chord cubics and exponential parameterization. Submitted.

14. Floater, M.S.: Chordal cubic spline interpolation is fourth order accurate. IMA Journal of Numerical Analysis 26, 25–33 (2006).

15. Wolfram Mathematica 9, Documentation Center, URL link: reference.wolfram.com/mathematica/guide/Mathematica.html.

16. Noakes, L., Kozera, R., Klette, R.: Length estimation for curves with $\varepsilon$-uniform samplings. In Digital Image Geometry, Eds. G. Bertrand, A. Imiya, and R. Klette, Lecture Notes in Computer Science 2243, 339–351, Springer-Verlag, Berlin Heidelberg (2001).

17. Janik, M., Kozera, R., Koziol, P.: Reduced data for curve modeling - applications in graphics, computer vision and physics. Advances in Science and Technology 7(18), 28–35 (2013)

18. Piegl, L., Tiller, W.: The NURBS Book. Springer-Verlag, Berlin Heidelberg, New York (1997).

19. Yan-Bin, J.: Interpolation by Spline Functions. Lecture Notes in Computer Science 477/577, Fall, Iowa State University (2014).

20. Jupp, D.L.P.: Approximation to data by splines with free knots. Journal on Numerical Analysis, 15(2), 328–343 (1978).

21. Budzko, D.A., Prokopenya, A. N.: Symbolic-numerical methods for searching equilibriums states in a restricted four-body problem. Programming and Computer Software 39(2), 74–80 (2013).

22. Kozera, R.: Uniqueness in shape from shading revisited. Journal of Mathematical Imaging and Vision 7(2), 123–138 (1997).

23. Taubin, G.: Estimation of planar and non-planar space curves defined by implicit equations with applications to edge and range segmentation. IEEE Transactions in Pattern Analysis and Machine Intelligence 13(11), 1115-1138 (1991).

24. Circularly symmetrical eikonal equations and non-uniqueness in computer vision. Journal of Mathematical Analysis and Applications. 165, 192–215 (1992).

25. Kocić, Lj.M., Simoncelli, A.C., Della Vecchia B.: Blending parameterization of polynomial and spline interpolants. Facta Universitatis (NIŠ), Series Mathematics and Informatics 5, 95–107 (1990).

26. Mørken, K., Scherer, K.: A general framework for high-accuracy parametric interpolation. Mathematics of Computation 66(217), 237–260 (1997).