Cindy González

Jorge Professeur

• Morin

Gudrun Albrecht

Professeure

Keywords: Algebraic Geometry [math Algebraic-Trigonometric Pythagorean-Hodograph curve, Pythagorean-Hodograph quintic curve, cubic B-spline curve, developable surface, computed tomography Courbe Algébrique Trigonométrique à Hodographe Pythagorien, courbe quintique à Hodographe Pythagorien, courbe B-spline cubique, facette développable, volume odontologique

Algebraic-trigonometric Pythagorean hodograph curves

for solving planar and spatial interpolation problems and their use for visualizing dental information within 3D tomographic volumes

Grâce à ces paramètres de forme, les courbes ATPH se révèlent plus flexibles et plus polyvalentes que leur homologue polynomial, les courbes quintiques à Hodographe Pythagorien (courbes PH) et les courbes polynomiales en général. Ces paramétres sont choisis de maniére optimale grâce à un processus de minimisation de fonctionnelles de forme. Nous réussissons de bien approximer des courbes trigonométriques connues par des courbes ATPH en ajustant les paramètres de forme.

Nous étendons les courbes ATPH planes à l'espace tri-dimensionnel, dans le chapitre List of Tables

Introduction

Free-form curves and surfaces have been widely studied in Computer Aided Geometric Design (CAGD) [START_REF] Farin | Curves and Surfaces for CAGD: A Practical Guide[END_REF][START_REF] Farin | Handbook of CAGD[END_REF]. Among the schemes used in this area to construct curves and surfaces is the concept of interpolation, which consists in designing curves and surfaces that satisfy prescribed geometrical constraints, such as exact points, tangents, curvatures, contact elements (lines/planes), etc. The most commonly used curve types in CAGD for interpolating 2D and 3D point sets are parametric curves such as polynomial splines, Bézier curves and B-spline curves.

Polynomial splines consist of piecewise polynomials that interpolate a set of data points and minimize certain energy functionals. Various continuity conditions can be specified in the data points to impose different degrees of smoothness of the resulting curve. The most popular spline for interpolating a sequence of points is the cubic spline consisting of piecewise cubic polynomials and joining with C 2 continuity in the data points. These curves were developed in the shipbuilding context since they simulate the behavior of wooden splines used for constructing ships, see, e.g., [START_REF] Liming | Practical Analytic Geometry with Applications to Aircraft[END_REF]. Subsequently, they were extensively studied and employed in many applications in CAGD, see, e.g., [START_REF] Liming | Practical Analytic Geometry with Applications to Aircraft[END_REF][START_REF] De Boor | Bicubic spline interpolation[END_REF][START_REF] Schoenberg | Cardinal interpolation and spline functions[END_REF][START_REF] Prenter | Splines and Variational Methods[END_REF][START_REF] Keys | Cubic convolution interpolation for digital image processing[END_REF][START_REF] Sahakian | Spline Functions and Multivariate Interpolations (Mathematics and Its Applications)[END_REF][START_REF] Gasca | Polynomial interpolation in several variables[END_REF].

In the mid 1950's a new class of parametric curves called Bézier curves were introduced by the engineers Pierre Bézier and Paul de Casteljau in the French automobile industry, see [START_REF] Farin | A History of Curves and Surfaces in CAGD[END_REF]. These curves are defined in terms of the Bernstein polynomials and they are controlled by a sequence of points called control points. By moving these points we can modify the shape of the curve. P. de Casteljau designed a recursive method to evaluate Bézier curves called De Casteljau's algorithm. This algorithm has become an important computational tool for the design of Bézier curves. In addition to this feature, Bézier curves possess significant properties such as affine invariance, end point interpolation, convex hull property, variation diminishing property, etc., which make them desirable for modeling smooth curves and surfaces in CAGD and related fields, see, e.g., [START_REF] Gordon | Bernstein-Bézier Methods for the Computer-Aided Design of Free-Form Curves and Surfaces[END_REF][START_REF] Böhm | A survey of curve and surface methods in CAGD[END_REF][START_REF] Farin | The Essentials of CAGD[END_REF].

Furthermore, we can define a Bézier spline curve as a sequence of Bézier segments joined at common end points with a certain continuity depending on the application. These splines have been used to solve planar and spatial interpolation problems see, e.g., [START_REF] Harada | Application of the Bézier curve to data interpolation[END_REF][START_REF] Long | Special Bézier quartics in three dimensional curve design and interpolation[END_REF][START_REF] Sarraga | G 1 interpolation of generally unrestricted cubic Bézier curves[END_REF][START_REF] Lee | Shape preserving and shape control with interpolating Bézier curves[END_REF][START_REF] Forrest | Interactive interpolation and approximation by Bézier polynomials[END_REF]. The Bézier curves were extended to well-known B-spline curves, which are piecewise polynomials defined by a sequence of control points and a knot vector. B-spline curves have similar properties as Bézier curves and B-splines provide more flexibility than the latter. These curves are a powerful tool to represent free-form shapes and they have been used extensively to solve 2D and 3D interpolation problems in CAGD, see, e.g., [START_REF] De Boor | A practical guide to splines[END_REF][START_REF] Farin | Fairing cubic B-spline curves[END_REF][START_REF] Unser | B-spline signal processing. I. Theory[END_REF][START_REF] Unser | B-spline signal processing. II. Efficiency design and applications[END_REF].

Additionally, within the family of parametric curves suitable to solve interpolating problems, are also the well-known Pythagorean Hodograph (PH) curves, which are polynomial parametric curves whose hodograph or derivative satisfies a pythagorean condition.

These curves have been extensively studied, see, e.g., [START_REF] Farouki | Pythagorean hodographs[END_REF][START_REF] Albrecht | Construction of C 2 Pythagorean-hodograph interpolating splines by the homotopy method[END_REF][START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF][START_REF] Pelosi | A control polygon scheme for design of planar C 2 PH quintic spline curves[END_REF][START_REF] Jüttler | Hermite interpolation by Pythagorean hodograph curves of degree seven[END_REF][START_REF] Farouki | Design of rational cam profiles with Pythagorean-hodograph curves[END_REF][START_REF] Šír | Constructing acceleration continuous tool paths using Pythagorean hodograph curves[END_REF][START_REF] Šír | C 2 Hermite interpolation by Pythagorean hodograph space curves[END_REF][START_REF] Farouki | Hermite interpolation by rotationinvariant spatial Pythagorean-hodograph curves[END_REF][START_REF] Farouki | Structural invariance of spatial Pythagorean hodographs[END_REF][START_REF] Farouki | Identification of spatial PH quintic Hermite interpolants with near-optimal shape measures[END_REF]. The PH curves possess significant computational advantages over polynomial curves in general. Among its most important properties are: its arc length can be computed exactly, they have rational offsets and they exhibit smoother curvature profiles than the polynomial curves in the case of interpolation of given points and tangents. A polynomial PH curve is usually obtained by the integration of its hodograph, which is represented by a complex polynomial for a planar PH curve or a quaternion polynomial for a spatial PH curve (see, e.g., [START_REF] Farouki | Pythagorean-Hodograph Curves: Algebra and Geometry Inseparable[END_REF]).

In the case of higher dimensional curves Clifford algebras are involved ( [START_REF] Choi | Clifford algebra, spin representation, and rational parameterization of curves and surfaces[END_REF][START_REF] Sakkalis | Pythagorean-hodograph curves in Euclidean spaces of dimension greater than 3[END_REF]). The lowest degree PH curves with sufficient shape flexibility for interpolation and approximation problems are cubic and quintic curves. Interpolation schemes using planar cubic PH curves can be found in [START_REF] Meek | Hermite interpolation with Tschirnhausen cubic spirals[END_REF][START_REF] Jaklič | Geometric Lagrange interpolation by planar cubic Pythagorean-hodograph curves[END_REF][START_REF] Byrtus | G 1 Hermite interpolation by PH cubics revisited[END_REF][START_REF] Jaklič | On interpolation by planar cubic G 2 Pythagorean-hodograph spline curves[END_REF] and spatial cubics PH curves are considered in [START_REF] Jüttler | Cubic Pythagorean-hodograph spline curves and applications to sweep surface modeling[END_REF][START_REF] Pelosi | Geometric Hermite interpolation by spatial Pythagorean-hodograph cubics[END_REF][START_REF] Kwon | Solvability of G 1 Hermite interpolation by spatial Pythagoreanhodograph cubics and its selection scheme[END_REF][START_REF] Jaklič | An approach to geometric interpolation by Pythagorean-hodograph curves[END_REF]. Lately, cubic PH curves have been extended to cycloidal curves, which were called Pythagorean Hodograph Cycloidal curves (shortly PHC curves), these curves were used to solve the C 1 Hermite interpolation problem, see [START_REF] Kozak | Pythagorean-hodograph cycloidal curves[END_REF]. For planar quintic PH curves, the C 1 Hermite interpolation problem was studied in [START_REF] Farouki | Hermite Interpolation by Pythagorean Hodograph Quintics[END_REF], C 2 splines that interpolate a set of data points were constructed in [START_REF] Albrecht | Construction of C 2 Pythagorean-hodograph interpolating splines by the homotopy method[END_REF][START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF][START_REF] Pelosi | A control polygon scheme for design of planar C 2 PH quintic spline curves[END_REF] and the G 1 Hermite interpolation problem is examined in [START_REF] Choi | Topological criterion for selection of quintic Pythagorean-hodograph hermite interpolants[END_REF]. In the case of spatial PH quintics, the problem of C 1 Hermite interpolation is done in [START_REF] Farouki | Hermite interpolation by rotationinvariant spatial Pythagorean-hodograph curves[END_REF][START_REF] Farouki | Structural invariance of spatial Pythagorean hodographs[END_REF], the construction of C 2 splines that interpolate given point data with prescribed end conditions is discussed in [START_REF] Farouki | Identification of spatial PH quintic Hermite interpolants with near-optimal shape measures[END_REF][START_REF] Šír | Spatial Pythagorean Hodograph Quintics and the Approximation of Pipe Surfaces[END_REF][START_REF] Farouki | Design of C 2 Spatial Pythagorean-Hodograph Quintic Spline Curves by Control Polygons[END_REF] and spatial G 2 continuous curves composed of pairs of PH quintic spiral segments is considered in [START_REF] Meek | G 2 curve design with a pair of Pythagorean hodograph quintic spiral segments[END_REF].

Recently, the PH curves have been extended to trigonometric Pythagorean Hodograph curves defined over a mixed algebraic-trigonometric space by Romani et al, in [START_REF] Romani | Algebraic-Trigonometric Pythagorean-Hodograph curves and their use for Hermite interpolation[END_REF]. This new class of PH curves was called Algebraic-Trigonometric Pythagorean Hodograph (ATPH) curves and they have an analogous property for its hodograph as in the case of PH curves. Besides, these curves enjoy the same computational advantages as their polynomial counterpart. Romani et al. solved the C 1 Hermite interpolation problem and the problem of joining G 2 continuously line segments and circles by ATPH curves of monotone curvature, see [START_REF] Romani | Algebraic-Trigonometric Pythagorean-Hodograph curves and their use for Hermite interpolation[END_REF]. Thanks to the shape parameter of these curves, the ATPH solutions of the aforementioned problems prove to be more flexible and versatile than the PH quintic solutions.

On the other hand, there are many methods that allow to solve the interpolation problems of spatial data sets by surfaces. In the case of interpolating scattered point clouds, parametric and non-parametric methods have been studied, such as, B-splines or NURBS [START_REF] Sarkar | Smooth-surface approximation and reverse engineering[END_REF][START_REF] Lee | Scattered data interpolation with multilevel B-splines[END_REF][START_REF] Kruth | Reverse engineering modelling of free-form surfaces from point clouds subject to boundary conditions[END_REF][START_REF] Caglar | B-spline interpolation compared with finite difference, finite element and finite volume methods which applied to twopoint boundary value problems[END_REF], radial basis functions [START_REF] Carr | Smooth surface reconstruction from noisy range data[END_REF], kernel smoothing techniques [START_REF] Roca-Pardiñas | From laser point clouds to surfaces: Statistical nonparametric methods for three-dimensional reconstruction[END_REF] and triangular meshes [START_REF] Mencl | Interpolation and approximation of surfaces from threedimensional scattered data points[END_REF][START_REF] Várady | Reverse engineering of geometric models-an introduction[END_REF]. In the case of interpolating lines and planes, methods to solve it by ruled surfaces have been proposed in several publications, see, e.g., [START_REF] Chen | Approximation by ruled surfaces[END_REF][START_REF] Hoschek | Interpolation and approximation with ruled surfaces[END_REF][START_REF] Peternell | G 1 Hermite interpolation of ruled surfaces[END_REF]. A surface is called ruled surface if it is generated by a one-parameter family of lines. Ruled surfaces are among the simplest surfaces used in design and modeling. A special case of ruled surfaces are developable surfaces, which can be unfolded into a plane without distortion.

We can find several publications on interpolation and approximation with this kind of surfaces, see, e.g., [START_REF] Aumann | Interpolation with developable bézier patches[END_REF][START_REF] Hoschek | Interpolation and approximation with developable B-spline surfaces[END_REF][START_REF] Pottmann | Approximation algorithms for developable surfaces[END_REF]. Developable surfaces have been used in engineering, design and manufacture, computer graphics and more recently in the medical field.

In the above interpolation contexts, the main objectives of this thesis are to solve practical interpolation problems related to 2D and 3D data with ATPH curves, as well as the interpolation of lines/plane pairs by developable surfaces. The latter surface interpolation aims at visualizing odontological information contained in 3D tomographic volumes. In particular, we solve the following three interpolation problems:

(1) Interpolating a planar sequence of points q 0 , . . . , q N together with prescribed end conditions by a planar C 2 ATPH spline curve.

(2) Constructing a spatial interpolating ATPH curve for any set of spatial C 1 Hermite data, i.e., two points p i , p f with associated first derivatives d i , d f .

(3) Interpolating pairs of consecutive lines/planes of a given sequence {(e k , τ k )} k=1,...,n by developable patches with G 1 continuity along spatial ATPH curves.

Initially, we study the construction of planar C 2 ATPH splines that interpolate a given sequence of points, which is described in Chapter 2. This construction involves the solution of a quadratic system in complex unknowns. In order to determine it, we extend the methods used to solve a similar problem by PH quintic curves (see [START_REF] Albrecht | Construction of C 2 Pythagorean-hodograph interpolating splines by the homotopy method[END_REF][START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF][START_REF] Pelosi | A control polygon scheme for design of planar C 2 PH quintic spline curves[END_REF]) to our case of ATPH solutions. In particular, we first consider a homotopy method that allows to find all possible ATPH spline curves that interpolate a given set of points.

we observe that by adjusting the shape parameters, the planar interpolating C 2 ATPH splines approximate very well, well-known trigonometric curves, such as circles, limaçons, cardiods, deltoids, etc.

As a natural extension of the planar ATPH curves, we study the construction of spatial ATPH curves in Chapter 3. The spatial ATPH curves are defined with respect to four functions in a trigonometric space. We observe that this definition is not very appropriate to attempt solving classical interpolation problems. Therefore, we use a representation in terms of quaternions that facilitates the formulation and analysis of such interpolants, as was done in the case of spatial PH quintic curves in [START_REF] Farouki | Hermite interpolation by rotationinvariant spatial Pythagorean-hodograph curves[END_REF][START_REF] Farouki | Structural invariance of spatial Pythagorean hodographs[END_REF]. In particular, we solve the problem of Hermite interpolation by spatial ATPH curves and we analyze the obtained solutions. The use of the quaternion representation of ATPH curves allows to reduce the aforementioned problem to the solution of a system of three equations in three unknowns. The solution of this system accepts a closed-form and it can be expressed as a three-parameter family, which contains all ATPH interpolants to given spatial Hermite data. In order to fix these three parameters, we consider the minimization of integral shape measures and we compare the resulting ATPH interpolants with their polynomial counterpart, the PH quintic Hermite interpolants. As in the planar case, we observe that calculating the free parameters by minimizing a suitable functional, the ATPH curve interpolating given end points and end tangents from known 3D parametric curves closely agrees with the corresponding parametric curve. Thus, the ATPH Hermite interpolation scheme has immediate application in the smooth interpolation of ordered sequences of points and tangents by piecewise ATPH curves. Therefore, we study the possibility of using these curves for constructing piecewise conical surfaces that allow to visualize information of interest in the odontological area in Chapter 4. For this purpose we revisit various publications that solve analogous problems. A first approach in that direction was presented in the Visible Human Project [START_REF] Hersch | Exploring Anatomic Structures with EPFL's Visible Human web server[END_REF], which is a portal that offers a computational tool for extracting information of organs along planar oblique slices. A more interesting problem is the information extraction along developable curved slices, because of the possibility of unfolding these surfaces into the plane without distortion, which might be useful in surgical planning. Particularly, surface flattening in medical imaging has been extensively studied in recent years, because it can help to detect anatomic abnormalities. For instance, in [START_REF] Haker | Nondistorting flattening maps and the 3D visualization of colon CT images[END_REF][START_REF] Vilanova Bartrolí | Virtual colon flattening[END_REF][START_REF] Hong | Conformal virtual colon flattening[END_REF] several algorithms have been proposed for unfolding surfaces obtained from 3D colon computed tomography (CT), which might be useful for detecting the presence of colon polyps and other pathologies. Angenent et al. proposed in [START_REF] Angenent | On the Laplace-Beltrami operator and brain surface flatenning[END_REF] a technique for flattening the brain surface in order to study its geometry. Another application was proposed in [START_REF] Saroul | Surface extraction and flattening for anatomical visualization[END_REF][START_REF] Hersch | Parallel unfolding visualization of curved surfaces extracted from large three-dimensional volumes[END_REF][START_REF] Paluszny | Between developable surfaces and circular cone splines: curved slices of 3D volumes[END_REF] for extracting information of organs along developable curved slices, which might become a useful tool for the medical professional to make more accurate diagnoses. For instance, the extraction of information about 3D volumes contained in cylinders built on plane curves was proposed by Figuereido and Hersch in [START_REF] Hersch | Parallel unfolding visualization of curved surfaces extracted from large three-dimensional volumes[END_REF]. This has been generalized to developable surfaces by Paluszny in [START_REF] Paluszny | Between developable surfaces and circular cone splines: curved slices of 3D volumes[END_REF]. But the disadvantages of these methods are that the resulting surfaces are not easy to adjust to a prescribed area and the unfolding of these surfaces requires a numerical method. Initially, to overcome these difficulties, we use the method presented by Leopoldseder and Pottmann in [START_REF] Leopoldseder | Approximation of developable surfaces with cone spline surfaces[END_REF] for constructing conical surfaces that model curved slices containing odontological information. The development of these surfaces can be explicitly calculated without numerical integration. However, this method still has some limitations, such as the possible presence of cone vertices within the curved slice or the appearance of the cone segments are not faithful to the anatomical structure of the teeth. Therefore, we study the design of cone segments along spatial C 1 ATPH curves. The cone splines are constructed by joining segments of cones with tangent continuity along given generators. We apply the construction of cone splines from C 1 ATPH curves for modeling a curved slice containing information about all dental pieces of a lower maxilla. This method allows to freely choose the vertex of each cone segment and the conical patches permit to visualize each dental piece of the jaw bone. In this way, the difficulties presented for the above methods are overcome and their advantages are kept. Moreover, the spatial C 1 ATPH interpolants exhibit smoother curvature and torsion profiles compared with those of spatial C 1 PH quintic curves. In this chapter, we are concerned with the construction of a planar C 2 ATPH spline curve that interpolates a sequence of points q 0 , . . . , q N . Different methods to solve a similar problem by PH quintic curves have been discussed in [START_REF] Albrecht | Construction of C 2 Pythagorean-hodograph interpolating splines by the homotopy method[END_REF][START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF][START_REF] Pelosi | A control polygon scheme for design of planar C 2 PH quintic spline curves[END_REF].

In particular, the articles [START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF][START_REF] Pelosi | A control polygon scheme for design of planar C 2 PH quintic spline curves[END_REF] employed the iterative Newton-Raphson method based on initial values estimated by comparing derivatives of the PH quintic spline and the cubic B-spline. We extend this method to the context of ATPH spline curves. The resulting C 2 ATPH spline curve depends on N free shape parameters for the given point sequence, we thus obtain a variety of interpolating splines. For an optimal choice of the values of these parameters we employ minimization techniques involving various fairness functionals.

The remainder of the chapter is organized as follows. In Section 2.2 we review the definition and the properties of ATPH curves. Section 2.3 is dedicated to develop the system of equations required for the construction of the C 2 ATPH interpolating splines which depend on several free shape parameters. In Section 2.4 we first formulate a homotopy method that yields all solutions of this system. We empirically notice that there is a "good" solution not exhibiting undesired loops. We thus use the Newton-Raphson iteration method to find this "good" solution in accordance with the shape of the data and illustrate it for representative examples. We then theoretically and numerically analyze the convergence behaviour of the iteration method by means of the Kantorovich theorem. Various fairness measures to design ATPH curves by automatically adjusting the free shape parameters are presented in Section 2.5. In Section 2.6 we describe how the ATPH curves can be chosen to closely agree with well-known trigonometric curves by adjusting the shape parameters. The investigation presented in this chapter has been published in [START_REF] González | Design of C 2 Algebraic-Trigonometric Pythagorean Hodograph splines with shape parameters[END_REF].

Algebraic Trigonometric Pythagorean Hodograph curves

Romani et al. ( [START_REF] Romani | Algebraic-Trigonometric Pythagorean-Hodograph curves and their use for Hermite interpolation[END_REF]) extended the well-known definition of polynomial PH curves to the algebraic trigonometric case, replacing the spaces of quadratic and quintic polynomials by (algebraic) trigonometric spaces with the normalized B-bases introduced by [START_REF] Sánchez-Reyes | Harmonic rational Bézier curves, p-Bézier curves and trigonometric polynomials[END_REF]. The B-bases have properties analogous to those of the Bernstein basis, yielding parametric curves that are defined by a control polygon as in the case of Bézier curves, see, e.g. [START_REF] Sánchez-Reyes | Harmonic rational Bézier curves, p-Bézier curves and trigonometric polynomials[END_REF][START_REF] Mainar | A general class of Bernstein-like bases[END_REF]. For t ∈ [0, α] and 0 < α < π, we first consider the trigonometric spaces

U 2 = 1,
B 2 0 (t) = cos(α -t) -1 cos(α) -1 , B 2 1 (t) = cos(α) -cos(t) -cos(α -t) + 1 cos(α) -1 , (2.1)
B 2 2 (t) = cos(t) -1 cos(α) -1 ,
1 (t) = 2(cos(α -t) -1)(cos(α) -cos(t) -cos(α -t) + 1) (cos(α) -1) 2 , B 4 2 (t) = 2(cos(α -t) -1)(cos(t) -1) + (cos(α) -cos(t) -cos(α -t) + 1) 2 (cos(α) -1) 2 , 4
B 4 3 (t) = 2(cos(t) -1)(cos(α) -cos(t) -cos(α -t) + 1) (cos(α) -1) 2 , B (2.2)
4 (t) = (cos(t) -1) 2 (cos(α) -1) 2 . 4
For t ∈ [0, α] and 0 < α < 2π we consider the mixed linear-trigonometric function vector space U 5 = 1, t, sin(t), cos(t), sin(2t), cos(2t) . The normalized B-basis of U 5 , denoted by B 5 i i=0,...,5 has the form ( [START_REF] Romani | Algebraic-Trigonometric Pythagorean-Hodograph curves and their use for Hermite interpolation[END_REF]): where we use the notation

B 5 0 (t) = 2 n 0 3(α -t) + sin(α -t)(cos(α -t) -4) , B 5
1 (t) = 4s 1 n 0 n 1 n 0 sin 4 α-t 2 -2s 4 1 3(α -t) + sin(α -t)(cos(α -t) -4) , B 5
c 1 = cos(α/2), s 1 = sin(α/2), c 2 = cos(α), s 2 = sin(α), (2.4)
and the abbreviations

a 0 = 6α + 2s 2 (c 2 -4), a 1 = c 1 (s 2 -3α) + 4s 1 , a 2 = (2 + c 2 )α -3s 2 .
(2.5)

Moreover, as noted in [START_REF] Sánchez-Reyes | Harmonic rational Bézier curves, p-Bézier curves and trigonometric polynomials[END_REF], the normalized B-bases

B 2 i (t) i=0,1,2 , B 4 i (t) i=0,.
..,4 and B 5 i (t) i=0,...,5 tend to the ordinary Bernstein polynomials of degree 2, 4 and 5 respectively, whenever α → 0. In [START_REF] Romani | Algebraic-Trigonometric Pythagorean-Hodograph curves and their use for Hermite interpolation[END_REF], a new class of PH curves has been defined, so-called Algebraic Trigonometric Pythagorean-Hodograph curves. They are defined over the parameter interval t ∈ [0, α] with 0 < α < π, using the fact that if f ∈ U 2 then f 2 ∈ U 4 and f 2 ∈ U 5 . Definition 2.1. Let u(t) and v(t) be nonzero functions from the space U 2 , relatively prime and both non-constant. A parametric curve r(t) = (x(t), y(t)) whose hodograph is of the form

x (t) = u 2 (t) -v 2 (t), y (t) = 2u(t)v(t) (2.6)
is called Algebraic-Trigonometric Pythagorean-Hodograph curve (ATPH), see [START_REF] Romani | Algebraic-Trigonometric Pythagorean-Hodograph curves and their use for Hermite interpolation[END_REF].

In the complex representation [START_REF] Farouki | The conformal map z → z 2 of the hodograph plane[END_REF] an ATPH curve r(t) = x(t) + iy(t) is obtained by squaring the complex function w(t) = u(t) + iv(t) yielding Since w is a complex function in the space U 2 , it can be written as

r (t) = x (t) + iy (t) = u 2 (t) -v 2 (t) + i2u(t)v(t) = w 2 (t).
w(t) = w 0 B 2 0 (t) + w 1 B 2 1 (t) + w 2 B 2 2 (t)
with B 2 k (t) defined in (2.1) and the coefficients w k ∈ C, k = 0, 1, 2. By integrating r (t) = w 2 (t) according to [START_REF] Romani | Algebraic-Trigonometric Pythagorean-Hodograph curves and their use for Hermite interpolation[END_REF], we obtain an ATPH curve in the space U 5 ,

r(t) = 5 k=0 p k B 5 k (t), (2.7)
whose complex control points p k = x k + iy k can be expressed as

p 1 = p 0 + n 0 16s 4 1 w 2 0 , p 2 = p 1 + n 0 -6n 2 8s 4 1 w 0 w 1 , p 3 = p 2 + n 2 4s 4 1 ((1 + c 2 )w 2 1 + w 0 w 2 ), (2.8
)

p 4 = p 3 + n 0 -6n 2 8s 4 1 w 1 w 2 , p 5 = p 4 + n 0 16s 4 1 w 2 2 ,
for an arbitrary constant of integration p 0 and the abbreviations s 1 , c 2 , n 0 , n 1 defined in (2.4) and (2.5).

As in the case of polynomial PH curves the ATPH curves [START_REF] Romani | Algebraic-Trigonometric Pythagorean-Hodograph curves and their use for Hermite interpolation[END_REF] possess important properties that are due to the Pythagorean condition (2.6):

• since the curve's parametric speed is given by σ(t) := (x (t)) 2 + (y (t)) 2 = w 2 (t) , the curve's arc-length can be computed explicitly without numerical quadrature as: where

σ(t)dt =
λ ij = Re( u i u j ) |u j | 2 , i, j ∈ {0, 1, 2} , and u0 = (1 + c2)w1 -w0 -w2 c2 -1 , u1 = s2(w0 -w1) c2 -1 , u2 = c2(w0 -w1) + w2 -w1 c2 -1 .
• For a planar curve r(t), the (signed) curvature has the expression

κ(t) = Im(r (t)r (t)) |r (t)| 3 .
Thus we have

κ(t) = 2 Im( w(t)w (t)) |w(t)| 4 = 2 u(t)v (t) -u (t)v(t) [u 2 (t) + v 2 (t)] 2 .
• The offset curve of an ATPH curve at each (signed) distance d, can be represented exactly by

r d = r(t) + dn(t), t ∈ [0, α],
where the unit normal vector n(t) to r(t), admits a precise rational parametrization given by

n(t) = -iw 2 (

C 2 ATPH spline equations

In this section we will extend the construction of C 2 PH quintic splines from [START_REF] Albrecht | Construction of C 2 Pythagorean-hodograph interpolating splines by the homotopy method[END_REF] to the case of ATPH splines.

Given an arbitrary sequence of parameters {α k } k=1,...,N , with 0 < α k < π, and a sequence of N + 1 complex points {q k } k=0,...,N , we wish to construct a C 2 ATPH spline that interpolates this set of points. This spline will be composed by N segments for k = 1, . . . , N . The segment r k (t) is defined as in (2.7) and its control points p k,i , i = 0, 1, . . . , 5 are expressed as in (2.8).

This construction of C 2 ATPH splines involves the solution of a quadratic system in complex variables. Before deriving these equations, we introduce the following notation and abbreviations for k = 1, . . . , N :

c 1,k = cos(α k /2), s 1,k = sin(α k /2), c 2,k = cos(α k ), s 2,k = sin(α k ), n 0,k = 6α k + 2s 2,k (c 2,k -4), n 1,k = c 1,k (s 2,k -3α k ) + 4s 1,k , n 2,k = (2 + c 2,k )α k -3s 2,k , (2.9)
l 1,k = s 2,k (c 2,k+1 -1), m 1,k = s 2,k (c 2,k-1 -1), x 1,k = l 1,k + m 1,k+1 , y 1,k = l 1,k l 1,k + m 1,k+1
,

where to avoid reference to undefined parameters α 0 , α N +1 , we set α 0 := α 1 and α N +1 := α N . Therefore

y 1,N = l 1,N l 1,N + m 1,N +1 = 1 2 (2.10)
and additionally we consider

y 1,0 = l 1,0 l 1,0 + m 1,1 = 1 2 .
(2.11) By using the above notation and abbreviations, the ATPH segment r k (t), for k = 1, . . . , N and its control points can be written as

r k (t) = 5 i=0 p k,i B 5 i (t),
where

p k,1 = p k,0 + n 0,k 16s 4 1,k w 2 k,0 , p k,2 = p k,1 + n 0,k -6n 2,0 8s 4 1,k w k,0 w k,1 , p k,3 = p k,2 + n 2,k 4s 4 1,k ((1 + c 2,k )w 2 k,1 + w k,0 w k,2 ), p k,4 = p k,3 + n 0,k -6n 2,k 8s 4 1,k w k,1 w k,2 , p k,5 = p k,4 + n 0,k 16s 4 1,k w 2 2,k .
(2.12)

C 2 continuity conditions

We now consider C 2 continuity conditions between two consecutive segments r k and r k+1 at their common point

q k = r k (α k ) = r k+1 (0).
Obviously, we can achieve C 0 continuity taking p k,0 = q k-1 and p k,5 = q k . Writing the hodograph of the segments r k (t) from q k-1 to q k in the form

r k (t) = w k,0 B 2 0 (t) + w k,1 B 2 1 (t) + w k,2 B 2 2 (t) 2 , (2.13)
we can ensure C 1 continuity requiring that w k+1,0 = w k,2 .

By differentiating (2.13) and imposing C 2 continuity at the common point of r k (t) and r k+1 (t), we obtain the condition

w k,2 =y 1,k w k,1 + 1 -y 1,k w k+1,1 , (2.14)
where y 1,k is given in (2.9). Therefore

w k,0 =y 1,k-1 w k-1,1 + 1 -y 1,k-1 w k,1 . (2.15)
By setting b k := w k,1 , the hodograph of the ATPH segment r k (t) can be rewritten in the complex form as

r k (t) = y 1,k-1 b k-1 + (1 -y 1,k-1 )b k B 2 0 (t) + b k B 2 1 (t) + y 1,k b k + (1 -y 1,k )b k+1 B 2 2 (t) 2 . (2.16)
This ensures that segments k and k + 1 exhibit continuity of first and second order derivatives at their common point q k = r k (α k ) = r k+1 (0).

Integrating this hodograph r k (t) subject to the interpolation condition

α k 0 r k (t)dt = ∆q k , with ∆q k = q k -q k-1 , (2.17)
yields the following system of N quadratic equations

f k (b k-1 , b k , b k+1 ) = δ k b 2 k-1 +ψ k b 2 k +θ k b 2 k+1 +φ k b k-1 b k +γ k b k b k+1 +λ k b k-1 b k+1 -k ∆q k = 0, (2.18)
for k = 1, . . . , N , where we introduce the quantities

δ k = n 0,k l 2 1,k-1 x 2 1,k , ψ k = n 0,k m 1,k x 2 1,k m 1,k + 2x 1,k-1 + n 0,k l 1,k x 2 1,k-1 l 1,k + 2x 1,k + 4n 2,k x 1,k x 1,k-1 m 1,k l 1,k + x 1,k x 1,k-1 (c 2,k + 1) -12n 2,k x 1,k x 1,k-1 m 1,k x 1,k + l 1,k x 1,k-1 , θ k = n 0,k m 2 1,k+1 x 2 1,k-1 , φ k = 2n 0,k l 1,k-1 x 2 1,k m 1,k + x 1,k-1 + 4n 2,k l 1,k-1 x 1,k x 1,k-1 l 1,k -3x 1,k , γ k = 2n 0,k m 1,k+1 x 2 1,k-1 l 1,k + x 1,k + 4n 2,k m 1,k+1 x 1,k x 1,k-1 m 1,k -3x 1,k-1 , λ k = 4n 2,k l 1,k-1 m 1,k+1 x 1,k x 1,k-1 , k = 16x 2 1,k x 2 1,k-1 s 4 1,k . (2.19)
We remark that every equation of the system (2. . In order to close the system, as in the case of "ordinary" splines, we need to provide end conditions. The following two possible forms are considered:

(i) complex values d 0 and d N for the end derivatives at the points q 0 and q N are specified, (ii) periodic end conditions are used when a closed C 2 curve is to be created.

In the following subsection we will discuss these conditions in detail.

End conditions

(i) For open curves, the derivatives r 1 (0) = d 0 and r N (α N ) = d N at the end points q 0 and q N are specified. We thus obtain two further equations

f 0 (b 0 , b 1 ) = y 1,0 b 0 + (1 -y 1,0 )b 1 2 -d 0 = 0, f N +1 (b N , b N +1 ) = y 1,N b N + (1 -y 1,N )b N +1 2 -d N = 0. (2.20)
The above equations together with (2.18) for k = 1, . . . , N comprise a system of

N + 2 quadratic equations for N + 2 unknowns b 0 , b 1 , . . . , b N , b N +1 .
(ii) For closed curves, the initial and final interpolating point must be coincident:

q N = q 0 . At this common point q N = q 0 of the first and last span, the continuity conditions r

N (α N ) = r 1 (0), r N (α N ) = r
1 (0), together with the fact that the complex variables b 1 , . . . , b N are considered as a cyclic list, yield the equations

r 1 (t) = ±l 1,N b N + m 1,1 b1 l 1,N + m 1,1 B 2 0 (t) + b1B 2 1 (t) + y 1,1 b1 + (1 -y 1,1 )b2 B 2 2 (t) 2 , r N (t) = y 1,N -1 bN-1 + (1 -y 1,N -1 )bN B 2 0 (t) + bN B 2 1 (t) + ±m 1,1 b1 + l 1,N bN l 1,N + m 1,1 B 2 2 (t) 2 .
The "good" solution may occur for either choice of ± and we will explain how to choose it in Section 2.4.3. By substituting the latter expressions in (2.17) we

obtain f 1 (b N , b 1 , b 2 ) = 0 and f N (b N -1 , b N , b 1 ) = 0, where f 1 (b N , b 1 , b 2 ) = δ 1 b 2 N + ψ 1 b 2 1 + θ 1 b 2 2 ± φ 1 b N b 1 + γ 1 b 1 b 2 ± λ 1 b N b 2 -1 ∆q 1 , (2.21)
f N (b N -1 , b N , b 1 ) = δ N b 2 N -1 + ψ N b 2 N + θ N b 2 1 + φ N b N -1 b N ± γ N b N b 1 ± λ N b N -1 b 1 -N ∆q N .
The above equations, together with (2.18) for k = 2, . . . , N -1, amount to a system of N quadratic equations for the unknowns b 1 , . . . , b N .

Solution of the C 2 ATPH spline equations

Given the sequences of real parameters {α k } k=1,...,N and of interpolating points {q k } k=0,...,N , as in [START_REF] Albrecht | Construction of C 2 Pythagorean-hodograph interpolating splines by the homotopy method[END_REF][START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF][START_REF] Pelosi | A control polygon scheme for design of planar C 2 PH quintic spline curves[END_REF] the total number of distinct C 2 ATPH splines results to be 2 N +1 for the system (2.18), (2.20) and 2 N for the system (2.18), (2.21).

Albrecht and Farouki in [START_REF] Albrecht | Construction of C 2 Pythagorean-hodograph interpolating splines by the homotopy method[END_REF] employed the homotopy method to compute all solutions in the case of PH quintic splines. They observed empirically that all but one of these interpolants exhibit undesirable "looping" behavior. In contexts where a high degree of interactivity is required, and where N is large, the need to examine all solutions in order to select the sole interpolant without self-intersections may seem intolerable. So, Farouki et al. in [START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF][START_REF] Pelosi | A control polygon scheme for design of planar C 2 PH quintic spline curves[END_REF] used the Newton-Raphson iteration as a suitable possibility to directly find this solution. We extend this strategy to the context of ATPH spline curves.

The choice of the appropriate initial values for the Newton-Raphson iteration is important for the convergence to the "good" solution. One possible approach is to compare derivative data from the ATPH spline and the cubic B-spline, both interpolating the same set of points, (as in the context of PH quintic splines, [START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF][START_REF] Pelosi | A control polygon scheme for design of planar C 2 PH quintic spline curves[END_REF]). Therefore, we choose the following moderate constraints on the points q k to be interpolated by the C 2 cubic B-spline

∆q k • ∆q k+1 > 0 and 1 2 ≤ ∆q k+1 |∆q k | ≤ 2 (2.22)
for k = 1, . . . , N -1 for an open curve or k = 1, . . . , N for a closed curve, where q N = q 0 and k is reduced modulo N .

On the other hand, Romani et al. in [START_REF] Romani | Algebraic-Trigonometric Pythagorean-Hodograph curves and their use for Hermite interpolation[END_REF], observed that for increasing values of α k ∈ (0, π), the arc length of the k -th segment increases, which eventually may cause undesirable loops that do not agree with the intuitive shape of the data. To avoid this, we require a "feasible" sequence of parameters α k . Hence the choice of α k is restricted to the interval (0, 2π/3). This choice is supported by the fact that the algebraic-trigonometric space U 5 is an extended Chebyshev space in this interval.

These conditions do not constitute a limitation for most practical applications.

The use of a homotopy method and the Newton-Raphson iterative method to construct the solution ATPH spline, considering the above observations, are described in the next sections.

Homotopy method

In order to examine all solutions to the C 2 ATPH interpolation problem, we have used a homotopy method. This numerical procedure determines all solutions to a given system of polynomial equations. In particular, we need to find the solutions to the system

(2.18), (2.20) of N + 2 quadratic equations in N + 2 complex variables b 0 , . . . , b N +1 .
We have adapted a predictor-corrector path-following program from an acknowledged code [START_REF] Sunyoung | CMPSm: A Continuation Method for Polynomial Systems (MATLAB version) -User's Manual[END_REF] to solve this system. It employs complex arithmetic and linear prediction steps in the homotopy parameter.

We shall use the linear homotopy defined by

h(b, t) = (1 -t)g(b) + tf (b), where b = (b 0 , . . . , b N +1 ). (2.23)
In our case, f is the system (2.18), (2.20) and the system g is obtained from f by omitting some terms from each of the equations

f k (b 0 , . . . , b N +1 ) = 0 for k = 0, . . . , N + 1.
Specifically,

g 0 (b 0 , . . . , b N +1 ) = y 2 0 b 2 0 -d 0 = 0, g k (b 0 , . . . , b N +1 ) = ψ k b 2 k -k ∆q k = 0, for k = 1, . . . , N, g N +1 (b 0 , . . . , b N +1 ) = (1 -y N ) 2 b 2 N +1 -d N = 0.
(2.24)

This initial system g has the closed-form solution, In general, we obtain 2 N +1 different sets of starting values for the variables b 0 , . . . , b N +1 by combining the signs in the above expressions.

b 0 = ± d 0 /|y 0 |, b k = ± k ∆q k /ψ k , for k = 1, . . . , N, b N +1 = ± d N /|(1 -y N )|. ( 2
So, beginning from a known solution of h(b, 0) with t = 0, the homotopy method deforms it continuously to obtain a solution of the objective system h(b, 1) ≡ f (b) = 0 at t = 1, employing a predictor-corrector method.

Empirical results

The homotopy method allows to construct the family of C 2 ATPH splines that interpolate the sequence of points q 0 , . . . , q N . We have observed through numerous examples that there is only one ATPH interpolant that does not exhibit loops and that preserves the shape of the data. To identify the "good" ATPH interpolant, we have employed two measures: the absolute rotation index and the bending energy

R abs = 1 2π |κ|ds E = κ 2 ds (2.26)
where κ and s are the curvature and the arc length, respectively.

Example

We interpolate the points q 0 = 1 + i2, q 1 = 3 + i5, q 2 = 6 + i4, q 3 = 7 + i3, q 4 = 9 + i4 by the family of C 2 ATPH spline curves. There are 32 distinct ATPH solutions that interpolate these five points (see Figure 2.3). Note that the "good" solution to the C 2 interpolation problem, which is traced in red in Figure 2.3, can be identified as the one with the minimum R abs and E (see Table 2.1). Specifically, R abs = 0.72 and E = 4.17.

Once the values b 0 , . . . , b N +1 have been calculated by the homotopy method, we make the assignments (2.15), (2.14) and compute the control points for each ATPH quintic span k = 1, . . . , N using equations (2.12). However, since the homotopy method becomes expensive for big N , we extend in the next section the strategy described in [START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF][START_REF] Pelosi | A control polygon scheme for design of planar C 2 PH quintic spline curves[END_REF] to directly compute the "good" solution in the context of ATPH splines.

Multivariable Newton-Raphson iteration

The Newton-Raphson iteration procedure starts with b (0) and is given by

b (k) = b (k-1) -J -1 (b (k-1) )f (b (k-1) ), (2.27)
for k ≥ 1, J is the Jacobian matrix for the given system f in (2.18) and b

(k) denotes b (k) = (b (k) 0 , . . . , b (k) N +1 ) for open curves or b (k) = (b (k) 1 , . . . , b (k) N )
for closed curves. The Jacobian matrix J for system (2.18) has the non-zero elements:

J k,k-1 = 2δ k b k-1 + φ k b k + λ k b k+1 , J k,k = φ k b k-1 + 2ψ k b k + γ k b k+1 , J k,k+1 = λ k b k-1 + γ k b k + 2θ k b k+1 .
(2.28)

In the case of specified end derivatives, the non-zero elements in the rows k = 0 and

k = N + 1 are J 0,0 = 2l 1,0 l 1,0 b 0 + m 1,1 b 1 /(l 1,0 + m 1,1 ) 2 , J 0,1 = 2m 1,1 l 1,0 b 0 + m 1,1 b 1 /(l 1,0 + m 1,1 ) 2 , J N +1,N = 2l 1,N l 1,N b N + m 1,N +1 b N +1 /(l 1,N + m 1,N +1 ) 2 , J N +1,N +1 = 2m 1,N +1 l 1,N b N + m 1,N +1 b N +1 /(l 1,N + m 1,N +1 ) 2 .
(2.29)

In the case of periodic end conditions the non-zero elements in the rows k = 1 and k = N are given by

J 1,N = 2δ 1 b N ± φ 1 b 1 ± λ 1 b 2 , J N,N -1 = 2δ N b N -1 + φ N b N ± λ N b 1 , J 1,1 = 2ψ 1 b 1 ± φ 1 b N + γ 1 b 2 , J N,N = 2ψ N b N + φ N b N -1 ± γ N b 1 , (2.30)
J 1,2 = 2θ 1 b 2 + γ 1 b 1 ± λ 1 b N , J N,1 = 2θ N b 1 ± γ N b N ± λ N b N -1 .
Regarding the CPU costs, in general, computing the inverse of a matrix with dimensions

N × N requires O(N 3
) operations, see e.g. [START_REF] Conte | Elementary Numerical Analysis: An Algorithmic Approach[END_REF]. The Jacobian matrix for the PH quintic spline is tridiagonal and has constant entries (see [START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF], page 36).

For ATPH splines the Jacobian matrix J for the system (2.18) is also tridiagonal and

depends on the coefficients δ k , φ k , λ k , ψ k , γ k , θ k for k = 1, . . . , N which are given by equa- tion (2.19). The increment vector δb (k-1) = -J -1 (b (k-1) )f (b (k-1)
) can be computed in O(N ) arithmetic operations [START_REF] Conte | Elementary Numerical Analysis: An Algorithmic Approach[END_REF].

The computation cost for ATPH is a bit higher than for the PH splines because of the Considering the quadratic convergence of the Newton-Raphson method [START_REF] Bellman | Methods of Nonlinear Analysis[END_REF] and the tridiagonal nature of J, we can deduce that the method calculates the "good" ATPH spline curve in a very efficient way, if one provides an initial approximation b (0) close to the solution.

As a convergence measure of the Newton-Raphson method, we use

k = N i=1 |f (k) i | 2 (2.31) where f (k) i = f (k) i (b (k) i-1 , b (k) i , b (k) i+1
) is the i-th equation of system (2.18), evaluated at the solution of the k-th Newton-Raphson iteration for k = 1, 2, . . ..

Initial values

To find the solution to the systems (2.18), (2.20) and (2.18), (2.21) with the Newton-Raphson method, which was described in the previous section, we need a starting approximation b (0) . As in [START_REF] Pelosi | A control polygon scheme for design of planar C 2 PH quintic spline curves[END_REF], the initial values are obtained by relating mid-point derivatives of the ATPH spline to mid-point derivatives of the cubic B-spline. For brevity, we restrict our attention to uniformly parametrized cubic B-splines interpolating the same data.

For open curves, we consider the cubic B-spline with knots

t 0 = t 1 = t 2 = t 3 = 0, t k = k -3 for k = 4, . . . , N + 2, and t N +3 = t N +4 = t N +5 = t N +6 = N .
In the case of closed curves, the knots t 0 , . . . , t N +6 for the cubic B-spline are regarded as a cyclical list. Given the knots t 3 < t 4 < . . . < t N +2 < t N +3 , the additional knots t 0 , t 1 , t 2 and t N +4 , t N +5 , t N +6 are determined by the equations:

t N +4 = t N +3 + (t 4 -t 3 ), t N +5 = t N +4 + (t 5 -t 4 ), t N +6 = t N +5 + (t 6 -t 5 ), t 2 = t 3 + (t N +2 -t N +3 ), t 1 = t 2 + (t N +1 -t N +2 ), t 0 = t 1 + (t N -t N +1 ).
By De Boor's algorithm we can calculate the cubic B-spline curve c(t) with nodal points q k = c(t k+3 ), k = 0, . . . , N , as well as their derivatives at the mid-points between the knots, c 1 2 (t k+2 + t k+3 ) for k = 1, . . . , N . Then we compare these derivatives with those from (2.16) evaluated at t = α k 2 and thus obtain the linear system

y 1,k-1 b k-1 + (1 -y 1,k-1 ) + 2c 1,k + y 1,k b k + (1 -y 1,k )b k+1 2(c 1,k + 1) = c 1 2 (t k+2 + t k+3 ) , (2.32)
for k = 1, . . . , N . As in the case of the general system (2.18), we must provide the following end conditions to close this system.

(i) For specified end derivatives the square roots of equations (2.20) provide two more equations:

y 1,0 b 0 + (1 -y 1,0 )b 1 = d 0 and y 1,N b N + (1 -y 1,N )b N +1 = d N . (2.33)
So the solution of the system (2.32), (2.33) for k = 1, . . . , N is the starting appro-

ximation b (0) = (b (0) 0 , b (0) 1 , . . . , b (0) N , b (0) N +1
) for this case.

(ii) The initial and final equations of the system (2.32) in the case of periodic end conditions are given by:

± l 1,N l 1,N +m 1,1 b N + m 1,1 l 1,N +m 1,1 + 2c 1,1 + y 1,1 b 1 + (1 -y 1,1 )b 2 2(c 1,1 + 1) = c 1 2 (t 3 + t 4 ) , (2.34
)

y 1,N -1 b N -1 + (1 -y 1,N -1 ) + 2c 1,N + l 1,N l 1,N +m 1,1 b N ± m 1,1 l 1,N +m 1,1 b 1 2(c 1,N + 1) = c 1 2 (t N + t N +1 ) .
Therefore, the initial value b

(0) = (b (0) 1 , . . . , b (0) N )
is obtained as the solution of the system (2.32), (2.34) for k = 2 . . . , N -1.

Since each complex square root on the right hand side of (2.32),(2.33) or (2.32), (2.34) admits two values, we will use the following procedure to select one as in [START_REF] Pelosi | A control polygon scheme for design of planar C 2 PH quintic spline curves[END_REF]. A particular sign for the square root on the right hand side of the first equation of (2.33), respectively (2.34), is chosen. Then the rest of the values are taken so that the dot product of the right hand side of two consecutive equations is non-negative, by interpreting complex numbers as vectors in R 2 . Finally, the dot product of the two square roots in equations (2.34) determines the sign ±.

Numerical examples

Numerous examples were executed for both open and closed curves and they were implemented in MATLAB. ..,N . We observe that the free parameters {α k } k=1,...,N act as shape parameter for the ATPH interpolants. This can clearly be seen in these figures, where different ATPH splines interpolate the same data.

We remark that according to the choice of the sequence of parameters {α k } k=1,...,N , some ATPH curves have more pleasing distribution of curvature than others. 2.

2).

Curve

N {α k } k=1,...,N

Interpolating points End-derivatives Error

Open curve 3 α 1 = 0.3671π,

q 0 = -2.1 + 1.8i, d 0 = 0.9918 -3.6084i, 0 = 295.3429 α 2 = 0.6667π, q 1 = -0.8333 -0.1i, d 3 = 1.4668 -1.9834i. 1 = 38.3156 α 3 = 0.3766π, q 2 = 0.9833 + 1.4167i, 2 = 1.1791 q 3 = 2.5667 + 0.6i.
3 = 0.013 4 = 0.00000000144 Closed curve 5 α 1 = 0.2228π, q 0 = 6.6670 + 0.6667i, 0 = 0.0982 α 2 = 0.2292π, q 1 = 8.5000 + 4.6670i, 1 = 0.0418 α 3 = 0.2626π, q 2 = 2.1667 + 4.5000i, 2 = 0.0243 α 4 = 0.2578π, q 3 = 0.3333 + 2.5000i, 3 = 0.0084 α 5 = 0.2610π, q 4 = 2.3333 -0.3333i, 4 = 0.00039 q 5 = 6.6670 + 0.6667i. can see that the ATPH interpolants have a slightly more pleasant curvature behaviour than their polynomial counterpart represented by the PH quintic spline curve. From this observation, a method to objectively choose the parameters {α k } k=1,...,N has been developed in section 2.5.

On the other hand, the homotopy method allows to compute the complete family of C 2 ATPH splines that interpolate q 0 , . . . , q N . In the case of open curves there are 2 N +1 distinct ATPH splines and 2 N for closed curves. This method requires to examine all the solutions for selecting a without undesired loop, which becomes very expensive for N > 10. For example, given a sequence of interpolating points q 0 , . . . , q 11 with q 0 = q 11 , we would have to compute 2 12 = 4096 different C 2 ATPH splines and to identify the "good" ATPH spline among all the solutions, as the one with smallest bending energy and smallest absolute rotation index. On the contrary, the Newton-Raphson method computes the "good" ATPH spline in an efficient manner. The convergence behavior of the Newton-Raphson iteration is largely successful. The method works very well even for a big number of interpolating points, for both open and closed curves. Given a sequence of parameters {α k } k=1,...,N , with α k = α for k = 1, . . . N , we have noted that for increasing values of α in (0, 2π/3) the C 2 ATPH spline curve becomes increasingly longer. The reason for this fact lies in the dependency on α of the control points of each ATPH segment curve according to equation (2.12). Numerous experiments suggest that for α k = α the total arc-length of the ATPH spline is monotonically increasing with respect to α. The Figure 2.11 shows an illustrative example.

In the case of a given sequence of parameters {α k } k=1,...,N with α j = α i for j = i, it is more difficult to characterize the behaviour of the C 2 ATPH spline curve for limit values of the shape parameters. This is due to the dependency on α k-1 , α k , α k+1 of the values w 0 , w 2 (see (2.15), (2.14)) of the control points in each segment (2.12).

Kantorovich theorem

We have implemented the Newton-Raphson method to construct C 2 ATPH interpolating splines. In the last section, we presented numerical examples of this construction and we have observed heuristically that this method converges rapidly to the desired solution. Moreover, these results are supported mathematically by the Kantorovich theorem (see [START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF]), which gives a convergence criterion for the Newton-Raphson method.

We will describe the use of this Kantorovich condition in the context of ATPH spline equations.

We employ the infinity norm ( . ∞ ) for vectors and matrices for formulating the convergence conditions. For simplicity, we write . instead of . ∞ . Given the sequence of parameters {α k } k=1,...,N with α k ∈ (0, 2π/3), we suppose:

1. J -1 (b (0) ) ≤ B, B ∈ R + 2. J -1 (b (0) )f (b (0) ) ≤ C, C ∈ R + 3. J(x) -J(y) ≤ K x -y for all x, y ∈ D and K ∈ R + ,
and define In our context the first two conditions of Theorem 2.2 are easy to test. Since the Jacobian J is tridiagonal, we simply calculate the exact inverse J -1 of the Jacobian at the starting approximation b (0) and we set B = J -1 (b (0) ) . In the same way, for condition 2 we explicitly determine the vector

h = BCK r = 1 - √ 1 -2h h C. If the conditions h ≤ 1 2 , S = {b | b -b (0) ≤ r} ⊂ D hold, then the Newton-Raphson iterations b (r+1) = b (r) -J -1 (b (r) )f (b (r) ), r = 1, 2 . . .
J -1 (b (0) )f (b (0) ) and set C = J -1 (b (0) )f (b (0) ) .
Finally, for condition 3 we require a Lipschitz constant K for the Jacobian matrix J, which is obtained according to the following Lemma.

Lemma 2.3. In condition 3 of Theorem 2.2 let D = C N +2 or D = C N . Given a
sequence of parameters {α k } k=1,...,N with α k ∈ (0, 2π/3), for the tridiagonal Jacobian matrix J, we obtain the following Lipschitz constants:

(i) If α k = α for k = 1 . . . , N, we have K = 40072, (ii) If α i = α j for i = j ∈ {1, . . . , N }, we have K = max 1≤k≤N {2, 2 * (| δ k | + | ψ k | + | λ k | + | φ k | + | γ k | + | θ k |)} for open curves or K = max 1≤k≤N {2 * (| δ k | + | ψ k | + | λ k | + | φ k | + | γ k | + | θ k |)} for closed curves.

Proof

(ii) For a given sequence of parameters {α k } k=1,...,N with α i = α j for i = j, x = (x 0 , . . . , x N +1 ), y = (y 0 , . . . , y N +1 ), and the Jacobian matrix defined by (2.28), (2.29) or (2.28), (2.30), we denote the sum of absolute values across row k of J(x) -J(y) by J k , which is

J k =| 2δ k (x k-1 -y k-1 ) + φ k (x k -y k ) + λ k (x k+1 -y k+1 ) | + | φ k (x k-1 -y k-1 ) + 2ψ k (x k -y k ) + γ k (x k+1 -y k+1 ) | + | λ k (x k-1 -y k-1 ) + γ k (x k -y k ) + 2θ k (x k+1 -y k+1 ) | ≤ 2 | δ k | + | φ k | + | λ k | | x k-1 -y k-1 | + (| φ k | +2 | ψ k | + | γ k |) | x k -y k | + | λ k | + | γ k | +2 | θ k | | x k+1 -y k+1 | with 1 ≤ k ≤ N for open curves or 2 ≤ k ≤ N -1 for closed curves
, where we use the triangle inequality.

In the case of open curves, for J 0 and J N +1 , we obtain

J 0 = - 2 l 1,0 + m 1,1 | l 1,0 (x 0 -y 0 ) + m 1,1 (x 1 -y 1 ) | ≤ 2l 1,0 l 1,0 + m 1,1 | x 0 -y 0 | + 2m 1,1 l 1,0 + m 1,1 | x 1 -y 1 | =| x 0 -y 0 | + | x 1 -y 1 |
and

J N +1 = - 2 l 1,N + m 1,N +1 | l 1,N (x N -y N ) + m 1,N +1 (x N +1 -y N +1 ) | ≤ 2l 1,N l 1,N + m 1,N +1 | x N -y N | + 2m 1,N +1 l 1,N + m 1,N +1 | x N +1 -y N +1 | =| x N -y N | + | x N +1 -y N +1 |,
where we have used (2.10) and (2.11).

In the case of closed curves, for J 1 and J N , we obtain

J 1 =| 2ψ 1 (x 1 -y 1 ) + σφ 1 (x N -y N ) + γ 1 (x 2 -y 2 ) | + | 2θ 1 (x 2 -y 2 ) + γ 1 (x 1 -y 1 ) + σλ 1 (x N -y N ) | + | 2δ 1 (x N -y N ) + σφ 1 (x 1 -y 1 ) + σλ 1 (x 2 -y 2 ) | ≤ 2 | ψ 1 | + | σ || φ 1 | + | γ 1 | | x 1 -y 1 | + (2 | θ 1 | + | σ || λ 1 | + | γ 1 |) | x 2 -y 2 | + 2 | δ 1 | + | σ || φ 1 | + | σ || λ 1 | | x N -y N | = 2 | ψ 1 | + | φ 1 | + | γ 1 | | x 1 -y 1 | + (2 | θ 1 | + | λ 1 | + | γ 1 |) | x 2 -y 2 | + 2 | δ 1 | + | φ 1 | + | λ 1 | | x N -y N | .
Analogously

J N =| 2θ N (x 1 -y 1 ) + σγ N (x N -y N ) + σλ N (x N -1 -y N -1 ) | + | 2δ N (x N -1 -y N -1 ) + φ N (x N -y N ) + σλ N (x 1 -y 1 ) | + | 2ψ N (x N -y N ) + σφ N (x N -1 -y N -1 ) + σγ N (x 1 -y 1 ) | ≤ 2 | θ N | + | λ N | + | γ N | | x 1 -y 1 | + (2 | δ N | + | λ N | + | φ N |) | x N -1 -y N -1 | + 2 | ψ N | + | γ N | + | φ N | | x N -y N | .
Thus, noting that

x -y = max 0≤k≤N +1 {| x k -y k |},
we obtain the following inequalities, which proves Lemma 2.

(ii).

For open curves

J(x) -J(y) = max 0≤k≤N +1 {J k } ≤ max 1≤k≤N {2, 2 * (| δ k | + | ψ k | + | λ k | + | φ k | + | γ k | + | θ k |)} x -y ,
and for closed curves

J(x) -J(y) = max 1≤k≤N {J k } ≤ max 1≤k≤N {2 * (| δ k | + | ψ k | + | λ k | + | φ k | + | γ k | + | θ k |)} x -y .
(i) In the case of a given sequence of parameters {α k } k=1,...,N with α k = α for k = 1, . . . , N by considering the fact that α ∈ (0, 2π/3) we have

J k =| 2δ(x k-1 -y k-1 ) + φ(x k -y k ) + λ(x k+1 -y k+1 ) | + =| φ(x k-1 -y k-1 ) + 2ψ(x k -y k ) + γ(x k+1 -y k+1 ) | + =| λ(x k-1 -y k-1 ) + γ(x k -y k ) + 2θ(x k+1 -y k+1 ) | ≤ 2 | δ | + | φ | + | λ | | x k-1 -y k-1 | +(| φ | +2 | ψ | + | γ |) | x k -y k | + | λ | + | γ | +2 | θ | | x k+1 -y k+1 | ≤ 2(| δ | + | γ | + | λ | + | φ | + | ψ | + | θ |) x -y ≤ 32(cos(α) -1) 4 sin 4 (α) 14 | sin(α) | (| cos(α) | +1) + 4α(| cos(α) | +2)(| cos(α) | +1) x -y ≤ 32 * 16 * (28 + 24 * 2π/3) x -y = 40072 x -y .
So we have obtained the inequality of Lemma 2.3 (i).

Empirical results

The global bound K = 40072 works for any α assuming α k = α for k = 1, . . . , N . It is not difficult to calculate Lipschitz constants K α for some specific α in the interval (0, 2π/3). Figure 2.12 and Table 2.5 show these Lipschitz constants K α , which are much smaller than the global bound K = 40072 given in Lemma 2.3 (i). In spite of the choice of each α k was restricted to the interval (0, 2π/3), since the space U 5 is a Chebyshev space in this interval, we calculated Lipschitz constants K α for α ∈ (0, π).

This allows to verify their behavior in the whole interval.

The Newton-Raphson method allows for an initial condition obtained from a cubic B-spline to converge to a "good" solution as considered in [START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF]. In fact, Newton-Raphson method takes the given initial condition to a point for which Lipschitz constants K α guarantees convergence. Tables 2.6, 2.7 illustrate the above assertions for some specific values of α. Note that in these cases Lipschitz constants K α are the same for both open and closed curves, this is because the aforementioned constant only depends on parameters α k = α for k = 1, . . . , N . The Kantorovich condition h < 1/2 is accomplished from the first iteration (h 1 see Table 2.6) for open curves, whereas for closed curves it is fulfilled from the fourth iteration (h 4 see Table 2

.7).

In spite of the fact that in the case of sequence of parameters {α k } k=1,...,N , with α i = α j for i = j the computations are more complex, we may report that the behaviour is similar to the case of α k = α for k = 1, . . . , N . So, to analyse the convergence of the Newton-Raphson method for {α k } k=1,...,N , with α i = α j and i = j, we have selected two different configurations of sequences for both open and closed curves. The first selection comprises 10 sequences of parameters whose values are chosen following a chordal heuristic. This means that for a given α 1 , the choice of the parameter α k with k > 1 depends on the distance between two consecutive interpolating points q k-1 and q k . The second arrangement is composed of 10 arbitrary increasing sequences of {α k } k=1,...,N . The results of the convergence parameter of Theorem 2.2, h i for i ≥ 0, in the case of open curves are listed in Tables 2.8, 2.9 and for closed curves in Tables 2.10, 2.11. We have observed that the Kantorovich condition h < 12 is achieved not only for constant sequences of parameters, but also for any sequences of parameters. Namely, taking chordal sequences of parameters, the Kantorovich condition h < 1 2 is accomplished after the second iteration (h 2 , see Table 2 As proposed in [START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF] we also have observed that the convergence condition h < 1/2 is obtained in less iterations, if we substitute the computed Lipschitz constant by an approximation K * ,

K * = J(b (0) ) -J(b) b (0) -b , (2.35)
where b (0) is the initial point for the Newton-Raphson method and b is the good 1 solution.

Tables 2.13, 2.12 and 2.14, 2.15, summarize the results obtained with K * for chordal and arbitrary sequences of the previously given parameters.

All the experiments performed allow us to conclude that the Newton-Raphson iteration is an efficient practical scheme to construct ATPH splines that interpolate point sequences. 2.5 Designing C 2 ATPH spline curves using fairness functionals

In the last sections, the construction of C 2 ATPH splines by a homotopy method and the Newton-Raphson method was presented. These curves are described by interpolation constraints, end condition constraints and sequences of free parameters for shape adjustment. For different choices of these free parameters {α k } k=1,...,N , we obtain various ATPH splines that interpolate the same data. We have remarked that according to the choice of {α k } k=1,...,N , some ATPH curves are aesthetically more pleasing than others. So, it is useful to develop a method for objectively selecting the sequence {α k } k=1,...,N . In consequence, we will apply nonlinear optimization methods to a pre-established fairness functional to iteratively find the parameters {α k } k=1,...,N that minimize the objective function.

Definition of fairness functionals

We use two traditional measures for fairness of curves, the bending energy (E) (2.26) and the minimum curvature variation (MV), see e.g. [START_REF] Albrecht | Invariante Gütekriterien im Kurvendesign Einige neuere Entwicklungen[END_REF][START_REF] Moreton | Minimum Variation Curves and Surfaces for Computer-Aided Geometric Design[END_REF][START_REF] Roulier | Measures of Fairness for Curves and Surfaces[END_REF]. The MV fairness measure is the integral of the square the derivative of curvature with respect to arc length (κ (s)) 2 ds.

(2.36)

The measure E is interpreted as the measure of the strain energy to bend a thin beam subject to given constraints. The measure (2.36) may be used to minimize the total variation of the curvature.

The value of the functionals E and MV change with a change of scale. So, taking in account this factor, we introduce the following functionals, see [START_REF] Albrecht | Invariante Gütekriterien im Kurvendesign Einige neuere Entwicklungen[END_REF] IM V = ( (ds))

3 (κ (s)) 2 ds, IES = κ 2 (s)ds ds, (2.37)
which are scale independent.

Implementation

The definition of fairness functionals is an essential first step for the design of faired ATPH spline curves. Our approach to design faired ATPH curves relies on the use of standard constrained optimization. We will summarize the generic algorithm as follows:

Algorithm 1. Choose a fairness functional to use as objective function g.

2.

Express g as a function of the parameters {α k } k=1,...,N .

3. Choose a numerical method for the evaluation of the chosen functional g.

4.

Give suitable starting values {α k } k=1,...,N for the algorithm to minimize g.

5.

Choose a multivariate minimization method.

6. The minimization method iteratively performs the following steps:

6a. Calculates the starting approximation as in the subsection 2.4.3.

6b. Uses Newton Raphson's method to construct the ATPH spline, see subsection

2.4.2.
6c. Evaluate the objective function g. Depending on its value the iteration is either stopped or a new set of parameters {α k } k=1,...,N is proposed.

The chosen multivariate minimization method for the examples presented in this section is the interior-point algorithm of the MATLAB "fmincon" constrained-optimization function, where the objective function g is one of the fairness functionals described above.

The numerical method for the evaluation of the fairness functional is the MATLAB function "integral".

Figure 2.17 shows the two ATPH splines that are produced by applying the fairness functional IMV (blue curve) and the bending energy IES (green curve), respectively, which interpolate the same set of points. The curves are visually similar and nearly indistinguishable on some segments, the curvature plots help to identify the shape differences between these curves. For the open ATPH curve, the value of IM V = 116.22 is greater than IES = 25.0506. Similarly, for the closed ATPH curve IM V = 256.6238 is greater than IES = 10.0259. This is because of the cubic factor in the definition of the functional IMV. However, the ATPH curve minimizing the fairness functional IMV has a smoother and lower curvature profile than the one obtained by minimizing the functional IES, indicating a slightly better shape for the first one than the last one.

Reproduction of well-known trigonometric curves

We have considered C 2 ATPH spline curves that interpolate a given set of points sampled from a well-known trigonometric curve such as circle, limaçon, cardioid, deltoid, etc. We have observed that by adjusting the sequences of parameters {α k } k=1,...,N , the interpolating C 2 ATPH spline curve closely agrees with the underlying trigonometric curve. For conciseness, we illustrate this behaviour for the cases of the unit circle and the limaçon, see Figures 2.18 and 2.19. The parametrization of the limaçon is given by the parametric representation x(θ) = (x(θ), y(θ)) = ((0.5 + cos(θ)) cos(θ), (0.5 + cos(θ)) sin(θ)).

We compare these curves with the approximations given by the C 2 PH quintic spline defined with the chordal parametrization given in [START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF] and the C 2 cubic B-spline.

To construct the interpolating C 2 ATPH spline and the C 2 PH quintic spline with the Newton-Raphson method described in section 2.4.2, we need a starting approximation b (0) . These initial values depend on mid-point derivatives of a known curve. In this case, we have used the mid-point derivatives of the respective trigonometric curve.

As a measure of the approximation of the trigonometric curves with C 2 ATPH spline curves and C 2 PH quintic spline curves, we use the approximation errors: maximum square error MSE and root-mean-square error RMSE, defined by

M SE = n k=1 | Ỹk -Y k | 2 n RM SE = √ M SE (2.38)
where Ỹk is the k-th point on the original trigonometric curve (called prediction point) and Y k is the k-th point on the interpolating curve (called observed point).

We consider arcs of planar trigonometric curves defined by parametric equations g(θ) = (x(θ), y(θ)) with θ ∈ [0, 2π], the interpolating C 2 ATPH segment r k (t) with t ∈ [0, α k ] for α k ∈ (0, 2π/3) with k = 1, . . . , N and the C 2 PH quintic spline segment r k,P H (t) defined with the chordal parametrization as in [START_REF] Farouki | Efficient solution of the complex quadratic tridiagonal system for C 2 PH quintic splines[END_REF].

In our case, the observed points are defined by Y k = r k ( • In the case of the unit circle, the sequence of parameters {α k } k=1,...,4 for the ATPH spline curve is given by α k = θ i -θ i-1 , where {θ i } i=0,...,4 = {0, π/4, π/3, π/2, 2π/3} . The four curves, the ATPH spline curve, the PH quintic spline curve, the cubic B-spline and the unit circle, are almost identical on the scale of the plot and their curvature profiles allow to visualize shape differences (see Figure 2.18). The distribution of the curvature (Figure2.18 Right) from the nominal value κ = 1 provides a measure of "smoothness". We can observe that 0.9992 ≤ κ ≤ 1.0019 for the interpolating C 2 ATPH spline whereas for the C 2 PH quintic spline we have 0.8301 ≤ κ ≤ 1.0245 and for the C 2 cubic B-spline we have 0.6104 ≤ κ ≤ 1.0724.

We have used a change of parameter to adjust the curvature plot of the ATPH spline curve to the interval [0, 1]. Quantity (2.31) illustrates the remarkable convergence of the Newton-Raphson method in this example, at the 4-th iteration the value of error (2.31) is already ε 4 =4.3837e-15. In addition the accurate approximation of the unit circle illustrated in Figure 2.18 with the C 2 ATPH spline curve is supported by the errors MSE and RMSE, which are 2.0362e-11 and 4.5124e-6, respectively. On the other hand these quantities for the interpolating C 2 PH quintic spline curve are 1.7036e-6 and 0.0013, respectively.

We have observed through many examples that the C 2 ATPH spline curves can be employed to interpolate large sequences of interpolating points on the unit circle.

In Table 2.16 we present comparative approximation errors: ε r-AT P H , MSE -AT P H and RMSE -AT P H for the ATPH curves and ε r-P H , MSE -P H and RMSE -P H for the PH quintic curves, where r denotes the r-th iteration of the Newton-Raphson method. Both curves interpolate the same data. For these data, we have selected non-uniform and increasingly larger sequences of interpolating points illustrated in Figure 2.20. The corresponding arcs of the unit circle, the interpolating C 2 ATPH spline curves and the interpolating C 2 PH quintic spline curves are in a very good agreement and therefore figures are not included. The approximation errors given by the C 2 ATPH spline curves are smaller than for the C 2 PH quintic spline curves as is clearly apparent in the results. Moreover, we obtain a "good" ATPH spline curve in the first iteration of the Newton-Raphson method and for the PH quintic spline curve more iterations are necessary. This somehow compensates for the higher computational cost for the evaluation of the Jacobian matrix in the case of ATPH spline curves.

• In the case of the limaçon, the interpolating points q i = (x i , y i ) with i = 0, . . . , 4 are given by x i (θ i ) = (0.5 + cos(θ i )) cos(θ i ) and y i (θ i ) = (0.5 + cos(θ i )) sin(θ i ), where {θ i } i=0,...,4 = {π/18, π/4, π/3, π/2, 2π/3} . So we have chosen the sequence of parameters {α k } k=1,...,4 as α k = θ k -θ k-1 . Figure 2.19 shows the C 2 ATPH spline curve, the C 2 PH quintic spline curve, the C 2 cubic B-spline and the respective limaçon, which are barely distinguishable on the scale of the plot. We can observe that the ATPH spline curve and the limaçon have similar curvature profiles (see investigated in [START_REF] Farouki | Hermite Interpolation by Pythagorean Hodograph Quintics[END_REF]. The extension of PH curves to the spatial case has thoroughly been studied in [START_REF] Farouki | Hermite interpolation by rotationinvariant spatial Pythagorean-hodograph curves[END_REF][START_REF] Farouki | Structural invariance of spatial Pythagorean hodographs[END_REF][START_REF] Farouki | Identification of spatial PH quintic Hermite interpolants with near-optimal shape measures[END_REF][START_REF] Farouki | Pythagorean-Hodograph Curves: Algebra and Geometry Inseparable[END_REF][START_REF] Choi | Clifford algebra, spin representation, and rational parameterization of curves and surfaces[END_REF]. Choi et al. [START_REF] Choi | Clifford algebra, spin representation, and rational parameterization of curves and surfaces[END_REF] presented an elegant characterization of spatial PH curves in terms of quaternions. This representation was used by Farouki et al. in [START_REF] Farouki | Hermite interpolation by rotationinvariant spatial Pythagorean-hodograph curves[END_REF][START_REF] Farouki | Pythagorean-Hodograph Curves: Algebra and Geometry Inseparable[END_REF] to solve the C 1 Hermite interpolation problem by spatial PH quintic curves. The resulting solutions have two degrees of freedom. Farouki et al. considered in [START_REF] Farouki | Identification of spatial PH quintic Hermite interpolants with near-optimal shape measures[END_REF] three methods to choose optimally these free angular parameters, which strongly influence the shape of the interpolants. Similarly, we consider the spatial ATPH curves to efficiently solve the C 1 Hermite interpolation problem. As done in [START_REF] Farouki | Hermite interpolation by rotationinvariant spatial Pythagorean-hodograph curves[END_REF][START_REF] Farouki | Structural invariance of spatial Pythagorean hodographs[END_REF] for spatial PH quintic curves, to simplify the construction of the Hermite interpolants a representation in terms of quaternions for spatial ATPH curves is used in this chapter and a general approach exploiting their properties is proposed.

The problem of Hermite interpolation by ATPH curves inherently implies a threeparameter family of solutions. The three parameters that determine spatial ATPH Hermite interpolants are angular variables (θ, β, α). We identify optimal values for them to select a "good" interpolant among the family of solutions.

The remainder of the chapter is organized as follows. In Section 3.2 we recall some results on quaternion algebra. In Section 3.

Preliminaries of quaternion algebra

The theory of complex numbers was extended to three dimensions by W. R. Hamilton [START_REF] Hamilton | On quaternions, or on a new system of imaginaries in algebra[END_REF], calling this new numeral system quaternions. The space of quaternions is denoted by H. Each quaternion number can be uniquely written as a linear combination of the four units 1, i, j, k with coefficients in R. Thus, a quaternion number has the form

A = a + a x i + a y j + a z k = a + a, (3.1)
where a, a x , a y , a z are real numbers and i, j, k denote three imaginary units. The component a of A is called scalar part of the quaternion and a is denominated its vector part.

The 3-dimensional vectors and real numbers are considered as pure vector and pure scalar quaternions, of the form A = a and A = a, respectively. We simply denote such quaternions as a and a.

Quaternion addition or subtraction, and scalar multiplication are defined in the usual manner in R 4 [START_REF] Hamilton | On quaternions, or on a new system of imaginaries in algebra[END_REF], that is

A ± B = (a ± b) + (a x ± b x )i + (a y ± b y )j + (a z ± b z )k = (a ± b) + (a ± b), rA = ra + ra x i + ra y j + ra z k = r(a + a),
for two any quaternion numbers A, B ∈ H and r ∈ R.

The basis elements 1, i, j, k satisfy the following identities:

i 2 = j 2 = k 2 = ijk = -1.
Therefore, all possible products of the basis elements are

ij = -ji = k, jk = -kj = i, ki = -ik = j. (3.2)
Using the above relations, the product of quaternions is given by

AB = (ab -a x b x -a y b y -a z b z ) + (ab x + ba x + a y b z -a z b y )i + (ab y + ba y + a z b x -a x b z )j + (ab z + ba z + a x b y -a y b x )k = ab -a • b + ab + ba + a × b, (3.3)
where the dot and the cross are the usual three dimensional scalar and vector cross products, respectively.

Quaternion multiplication is non-commutative, i.e. AB = BA. However, all other algebraic properties hold, namely, associativity, existence of a neutral element and inverse of elements with respect to multiplication and addition.

For every quaternion A = a + a, there is a conjugate quaternion A * = a -a, and the norm of A is the non-negative quantity |A| defined by

|A| = √ AA * = a 2 + a 2 x + a 2 y + a 2 z = a 2 + |a| 2 .
Besides, for any two quaternions A, B holds |AB| = |BA| and product conjugation satisfies the relation (AB) * = B * A * .

Quaternion equations

The interpolation of first-order Hermite data by spatial ATPH curves, as in the case of PH quintic curves [START_REF] Farouki | Hermite interpolation by rotationinvariant spatial Pythagorean-hodograph curves[END_REF]- [START_REF] Farouki | Structural invariance of spatial Pythagorean hodographs[END_REF], involves equations of the form

AiA * = d, (3.4)
with the unknown quaternion A = a 0 +a x i+a y j+a z k and a given pure vector quaternion

d = d x i + d y j + d z k.
We use the solutions of this system given in [START_REF] Farouki | Hermite interpolation by rotationinvariant spatial Pythagorean-hodograph curves[END_REF]- [START_REF] Farouki | Structural invariance of spatial Pythagorean hodographs[END_REF].

Note that the system of equations (3.4) splits into three equations in four unknowns, namely

a 2 0 + a 2 x -a 2 y -a 2 z = d x , 2(a 0 a z + a x a y ) = d y , 2(a x a z -a 0 a y ) = d z .
So, the solutions of (3.4) have one degree of freedom. Particularly, taking a 0 = 0 and d |d| = (λ x , λ y , λ z ), the real solutions of (3.4) are given by

a x = ± 1 2 (1 + λ x ) | d |, a y = ± λ y 1 + λ x 1 2 (1 + λ x ) | d |, a z = ± λ z 1 + λ x 1 2 (1 + λ x ) | d |.
Therefore,

A = ±(i + λ y 1 + λ x j + λ z 1 + λ x k) 1 2 (1 + λ x ) | d |. (3.5)
On the other hand, consider a quaternion B satisfying the following equation

BiB * = i. (3.6)
Then the quaternion à = AB is also a solution of (3. Note that in the previous equation there is a combination of scalars, vectors and quaternions and the quaternion product is implicit.

Spatial Algebraic Trigonometric Pythagorean Hodograph curves

Romani et al. [START_REF] Romani | Algebraic-Trigonometric Pythagorean-Hodograph curves and their use for Hermite interpolation[END_REF] extended the well-known definition of planar polynomial Pythagorean Hodograph (PH) curves to the algebraic-trigonometric case, calling this new class of planar PH curves, Algebraic Trigonometric Pythagorean Hodograph (ATPH) curves. The complex model for planar ATPH curves was studied in chapter 2, where the hodograph r (t) = (x (t), y (t)) of an ATPH curve r(t) is obtained by squaring a complex function

w(t) = u(t) + iv(t) yielding w 2 (t) = u 2 (t) -v 2 (t) + i2u(t)v(t)
. This hodograph is identified as the real and imaginary parts of w 2 (t). An analogous approach for spatial ATPH curves in terms of quaternions is studied in this section.

Recall that a parametric space curve is called a PH curve (cf. [START_REF] Farouki | Hermite interpolation by rotationinvariant spatial Pythagorean-hodograph curves[END_REF]), if the first derivative is the square of another polynomial. Specifically, r(t) = (x(t), y(t), z(t)) is a space PH curve if

x 2 (t) + y 2 (t) + z 2 (t) = σ 2 (t) (3.10)
for some polynomial σ(t). According to [START_REF] Farouki | Hermite interpolation by rotationinvariant spatial Pythagorean-hodograph curves[END_REF], the condition (3.10) holds if and only if

x (t) = u 2 (t) + v 2 (t) -p 2 (t) -q 2 (t), y (t) = 2u(t)q(t) + 2v(t)p(t), z (t) = 2v(t)q(t) -2u(t)p(t).
for some polynomials u(t), v(t), p(t), q(t).

The spatial ATPH curves are defined similarly, but with respect to four functions u(t), v(t), p(t), q(t) in the trigonometric space U 2 introduced in chapter 2.

Definition 3.1. Let u(t), v(t), p(t), q(t) be non-zero functions in the space U 2 . Then a spatial parametric curve r(t) = (x(t), y(t), z(t)) with t ∈ [0, α] and α ∈ (0, 2π/3) is called spatial Algebraic Trigonometric Pythagorean Hodograph curve or spatial ATPH curve, if its first derivative is of the form

x (t) = u 2 (t) + v 2 (t) -p 2 (t) -q 2 (t), y (t) = 2u(t)q(t) + 2v(t)p(t), z (t) = 2v(t)q(t) -2u(t)p(t).
(3.11)

The representation (3.11) can be reformulated using quaternions [START_REF] Dietz | An algebraic approach to curves and surfaces on the sphere and on other quadrics[END_REF]. Consider the pythagorean hodograph as a pure vector quaternion function in the form

r (t) = [u 2 (t) + v 2 (t) -p 2 (t) -q 2 (t)]i + 2[u(t)q(t) + v(t)p(t)]j + 2[v(t)q(t) -u(t)p(t)]k.
This representation can be expressed in terms of a quaternion function A(t). Namely,

r (t) = A(t)iA * (t), (3.12)
where A(t) = u(t) + v(t)i + p(t)j + q(t)k with u(t), v(t), p(t), q(t) ∈ U 2 and A * (t) denotes the conjugate of A(t), see [START_REF] Choi | Clifford algebra, spin representation, and rational parameterization of curves and surfaces[END_REF].

Since A(t) is a quaternion function in the space U 2 , it can be written as follows,

A(t) = A 0 B 2 0 (t) + A 1 B 2 1 (t) + A 2 B 2 2 (t), (3.13)
with quaternion coefficients A r = u r + v r i + p r j + q r k, r = 0, 1, 2. Note that (3.12) defines r (t) as a pure vector quaternion. By integrating the hodograph defined in (3.12), we obtain a parametric curve in the algebraic-trigonometric space U 5 , which can be expressed in the normalized B-basis of U 5 as

r(t) = 5 k=0 p k B 5 k (t), t ∈ [0, α] (3.14)
with control points p k = x k i + y k j + z k k defined by

p 1 = p 0 + n 0 16s 4 1 A 0 iA * 0 , p 2 = p 1 + n 0 -6n 2 16s 4 1 (A 0 iA * 1 + A 1 iA * 0 ), p 3 = p 2 + n 2 8s 4 1 (A 0 iA * 2 + 2(1 + c 2 )A 1 iA * 1 + A 2 iA * 0 ), p 4 = p 3 + n 0 -6n 2 16s 4 1 (A 1 iA * 2 + A 2 iA * 1 ), p 5 = p 4 + n 0 16s 4 1 A 2 iA * 2 , (3.15)
for an arbitrary constant of integration p 0 and the abbreviations s 1 , c 2 , n 0 , n 1 defined in (2.4) and (2.5).

C 1 Hermite interpolation problem

In this section we are concerned with the Hermite interpolation problem. More precisely, the ATPH curve is used to interpolate prescribed end points p i , p f and tangent vectors at these points. Henceforth, the tangent vectors at p i , p f will be denoted by d i , d f , respectively. We use the quaternion representation for spatial ATPH curves with hodographs of the form (3.12).

ATPH interpolants solving the C 1 Hermite problem

In order to determine conditions under which end points p i , p f and end-derivatives d i , d f can be interpolated by a spatial ATPH curve, we propose a variation of the method proposed in [START_REF] Farouki | Hermite interpolation by rotationinvariant spatial Pythagorean-hodograph curves[END_REF] to solve the Hermite interpolation problem by polynomial PH curves.

To solve this interpolation problem it is convenient to write the data as pure quaternions, namely

p i = p xi i + p yi j + p zi k d i = d xi i + d yi j + d zi k, p f = p xf i + p yf j + p zf k d f = d xf i + d yf j + d zf k.
The interpolation conditions

r (0) = d i , r (α) = d f ,
then yield the system of equations for A 0 and A 2 ,

A 0 iA * 0 = d i A 2 iA * 2 = d f . (3.16)
Writing the vectors d i and d f as

d i |d i | = (λ xi , λ yi , λ zi ) and d f |d f | = (λ xf , λ yf , λ zf
), the equations (3.16) can be solved [START_REF] Farouki | Hermite interpolation by rotationinvariant spatial Pythagorean-hodograph curves[END_REF] as according to (3.8)

A 0 = -sin φ 0 + cos φ 0 i + λ yi cos φ 0 + λ zi sin φ 0 1 + λ xi j + λ zi cos φ 0 -λ yi sin φ 0 1 + λ xi k 1 2 (1 + λ xi ) | d i |, A 2 = -sin φ 2 + cos φ 2 i + λ yf cos φ 2 + λ zf sin φ 2 1 + λ xf j + λ zf cos φ 2 -λ yf sin φ 2 1 + λ xf k 1 2 (1 + λ xf ) | d f |, (3.17)
where D = ρA 0 + ωA 1 + ρA 2 . This equation has the form (3.4), therefore we can deduce a solution from (3.8) as follows

D = -sin φ 1 + cos φ 1 i + λ y cos φ 1 + λ z sin φ 1 1 + λ x j + λ z cos φ 1 -λ y sin φ 1 1 + λ x k 1 2 (1 + λ x ) | d |,
where φ 1 denotes an angular parameter and d |d| = (λ x , λ y , λ z ). Therefore,

A 1 = 1 ω -sin φ 1 + cos φ 1 i + λ y cos φ 1 + λ z sin φ 1 1 + λ x j + λ z cos φ 1 -λ y sin φ 1 1 + λ x k 1 2 (1 + λ x ) | d | - ρ ω (A 0 + A 2 ).
Note that A 1 depends on parameters φ 0 , φ 1 , φ 2 .

The three quaternions coefficients A i , i = 0, 1, 2 can be rewritten in the form (3.9), so

A 0 = | d i |n i exp(φ 0 i), A 1 = - ρ ω (A 0 + A 2 ) + 1 ω | d |n exp(φ 1 i), A 2 = | d f |n f exp(φ 2 i), (3.19)
where

δ i = d i | d i | , δ d = d | d | , δ f = d f | d f | , n i = δ i + i | δ i + i | , n = δ d + i | δ d + i | , n f = δ f + i | δ f + i | and exp(φ r i) = cos φ r + sin φ r i, r = 0, 1, 2.
Note that the hodograph (3.12) depends on differences of the angles φ 0 , φ 1 , φ 2 as follows.

The control points p k , k = 0, . . . , 5, defined by equation (3.15) depend on the product of quaternions A r iA * s for r, s ∈ 0, 1, 2. Each quaternion A r , A s can be expressed as A r = A r (0)(cos φ r + sin φ r i) and A s = A s (0)(cos φ s + sin φ s i) where A r (0), A s (0) is the value when φ r = 0 and φ s = 0, respectively. Therefore

A r iA * s = [A r (0)(cos φ r + sin φ r i)]i[A s (0)(cos φ s + sin φ s i)] * = A r (0)[(cos φ r + sin φ r i)i(cos φ s -sin φ s i)]A * s (0) = A r (0)[(cos φ r sin φ s -cos φ s sin φ r ) + (sin φ r sin φ s + cos φ s cos φ r )i]A * s (0) = A r (0)[sin(φ s -φ r ) + cos(φ s -φ r )i]A * s (0).
This means that A r iA * s depend on the differences of the angular parameters φ 0 , φ 1 , φ 2 . Consequently, without loss of generality, one of the parameters can be taken equal to zero, say φ 1 . Then the three quaternions A 0 , A 1 , A 2 (3.19) can be written in terms of two angles θ = 1 2 (φ 0 + φ 2 ) and β = φ 2 -φ 0 . So, we have that

A 0 = | d i |n i exp((θ -1 2 β)i), A 2 = | d f |n f exp((θ + 1 2 β)i), A 1 = - ρ ω (A 0 + A 2 ) + 1 ω | d |n. (3.20) with (θ, β) ∈ [0, 2π] × [0, 2π].

Empirical results

In this section we present examples of spatial ATPH curves, constructed from given Hermite data p i ,d i and p f ,d f using the method described above. We compare these curves with the polynomial PH quintic solving the same C 1 Hermite problem.

As for the polynomial PH curves, the solution to the ATPH Hermite interpolation problem, depends on the two parameters θ, β. Furthermore in our case, there is one more free parameter α which can be used as shape parameter.

Firstly, for illustrating the method to construct ATPH interpolants, described in section 3.4, we have used the Hermite data given in Table 3. According to the choice of the parameters (θ, β, α) some ATPH curves are more aesthetically pleasing than others. So, to select these parameters such that the resulting ATPH curve satisfies a fairness criterion, we have opted to apply an optimization method to a suitable shape integral for the ATPH Hermite interpolants with respect to these parameters. This method is explained in the next section.

Selection of angular parameters (θ, β, α)

The solution to the ATPH Hermite interpolation problem, specified by (3.19), depends on three free angular parameters (θ, β, α). As one criterion to optimally calculate these parameters, we have considered the minimization of integral shape measures.

We use the concepts of curvature, torsion and arc length to define numerical measures of the fairness of ATPH curves. Recall that the Frenet formulas for spatial parametric curves are defined as (see, e.g., [START_REF] Carmo | Differential Geometry of Curves and Surfaces[END_REF])

t = r | r | n = r × r | r × r | × t b = t × n, (3.21)
where t, n, b are the unit tangent, normal and binormal vectors, respectively, for r. Let κ and τ be the curvature and torsion, respectively, given by

κ = | r × r | | r | 3 τ = (r × r ) • r | r × r | 2 . (3.22)
To calculate the curvature, torsion and their derivatives, it is required to express the first, second, third and fourth order curve derivatives in terms of the quaternion representation of r , so we have that

r = AiA * r = A iA * + AiA * r = A iA * + 2A iA * + AiA * r (4) = A iA * + 3A iA * + 3A iA * + AiA * .
The integral fairness measures that we use for calculating the free parameters (θ, β, α), were studied in [START_REF] Albrecht | Invariante Gütekriterien im Kurvendesign Einige neuere Entwicklungen[END_REF], [START_REF] Roulier | Measures of Fairness for Curves and Surfaces[END_REF], which are

F 1 = κ 2 ds, F 2 = κ 2 ds F 3 = (κ 2 + τ 2 )ds. (3.23)
Note that the functionals F 1 , F 2 , F 3 are invariant under rigid body transformations.

Since the value of these functionals change with a change of scale, we modify them to obtain functionals whose values are independent of scale, by including an arc length term to offset the scaling factor, F1 = ds κ 2 ds,

F2 = ds 3 κ 2 ds F3 = ds 3 (κ 2 + τ 2 )ds, (3.24)
These functionals are scale invariant and they preserve the properties of the functionals F 1 , F 2 , F 3 , i.e., its characteristic shape (geometric invariants of the curve) and parametrization independence. Thus two different parameterizations of the same curve will have the same value of the measure of fairness.

The functionals given in (3.24) are defined in terms of a parametrized function by arc length. To evaluate them, each one must be converted to the parameter t. For this, the differential with respect to s is converted to a differential in t.

Since dt ds = 1 |r (t)| then ds =| r (t) | dt, where | r (t) |= r (t) • r (t) 1/2 .
This yields the following expression for F1 , F2 , F3

F1 = | r | dt κ 2 | r | dt, F2 = | r | dt 3 κ 2 | r | dt, F3 = | r | dt 3 κ 2 + τ 2 1 | r | dt. (3.25)
Finally, we express dκ dt and dτ dt in terms of the derivatives of the ATPH curve r(t). The derivative of κ with respect to t is given by

κ = vdu -udv v 2
where

u = r × r du = r × r v =| r | 3 dv = 3 | r | (r • r )
Taking the derivative of τ with respect to t yields

τ = wdx -xdw w 2 where x = (r × r ) • r dx = (r • r ) • r (4) + (r × r ) • r w = (r × r ) • (r × r ) dw = 2(r × r ) • (r × r ).
For conciseness, the parameter t is omitted in the last equations.

Implementation

As in the planar case, the definition of fairness functionals for designing faired spatial ATPH curves is a fundamental step. This has been considered in the previous section.

Next we apply nonlinear optimization methods to these functionals, in order to objectively select the free parameters of the ATPH curves.

In section 3.4.1, we described the spatial ATPH curves as parametric curves defined by the prescribed end points p i , p f , tangent vectors d i , d f at these points and three free angular parameters (θ, β, α) for shape adjustment. The automatic generation of a faired ATPH curve consists of finding these parameters such that the resulting ATPH curve minimizes one of the fairness functionals described previously. For this, the following algorithm is carried out.

Algorithm

• Choose a fairness functional among F1 , F2 , F3 to be used as objective function g.

• Express the chosen functional as a function of the free parameters (θ, β, α).

• Choose a numerical method for the evaluation of the function g.

• Choose a minimization algorithm and starting values for the algorithm to minimize g. Execute the minimization using the numerical method chosen in the preceding step.

Our approach to implement this algorithm in the MATLAB language, has been to choose as minimization algorithm the constrained-optimization function "fmincon" and as numerical method for evaluating g the function "integral".

Numerical examples

We present examples of C 1 ATPH Hermite interpolants to illustrate the selection criterion of the three angular parameters (θ, β, α) explained above. We compare these curves with their polynomial counterpart, the well-known PH quintic Hermite interpolants. We use the five sets of end derivatives listed in Table 3.3, which were taken from [START_REF] Farouki | Identification of spatial PH quintic Hermite interpolants with near-optimal shape measures[END_REF]. In each case, the end points are p i = (0, 0, 0), We have applied the optimization method to the three functionals F1 , F2 , F3 , to calculate the angular parameters (θ, β, α) for ATPH curves and the parameters (θ, β) for PH

p f = (1,

Approximation of arcs of 3D parametric curves

We consider spatial C 1 ATPH curves that interpolate given end points and associated end-derivatives from known 3D parametric curves. We have observed that calculating the free parameters (θ, β, α) by minimizing an integral shape measure, the interpolating ATPH curve closely approximates the corresponding space parametric curve. To illustrate the performance, we present the process for some kinds of helices. We compare the approximations accomplished through the C 1 ATPH curves with the C 1 PH quintic curves.

We use the approximation errors MSE and RMSE given by the equation (2.38) as measures of the approximation of arcs of 3D parametric curves with ATPH curves and PH quintic curves. Recall that to evaluate the quantities MSE and RMSE, we need to specify the prediction points Ỹk for k = 1, . . . , n (points on the original 3D parametric curve)

and the observed points Y k with k = 1, . . . , n (points on the interpolating 3D curve), see section 2.6.

We consider a discretization {t k } k=1,...,n of n evenly spaced points in the interval [0, α] with α ∈ (0, 2π/3), for the interpolating ATPH curve r(t). Thus the observed points in this case are given by Y k = r(t k ) for k = 1, . . . , n.

Similarly, for the interpolating PH quintic curve r P H (t) with t ∈ [0, 1], the observed points are determined by Y k = r P H (t k ), where {t k } k=1,...,n is an uniform discretization of the interval [0, 1].

Finally, to calculate the prediction points Ỹk , we use the code distance2curve of MATLAB, which allows to determine the nearest point on the original curve to any given point.

For all examples, the free parameters (θ, β, α) for the ATPH curves and the parameters (θ, β) for the PH quintic curves, are calculated with the minimization process of the fairness functional F3 .

• Circular helix: a helix is a smooth 3-dimensional curve, which is parametrized by:

x = r cos t, y = r sin t, z = at,
where r is the radius and a is a constant that give the rise per turn. The curvature and the torsion of a helix are given by

κ = r r 2 + a 2 , τ = a r 2 + a 2 .
In particular, we consider the simplest helix with equations

x = cos t, y = sin t, z = t,
where r = 1 and a = 1. • Hyperbolic helix: a space curve with parametric equations

x = cosh t, y = sinh t, z = t,
is known as a hyperbolic helix. The curvature and torsion of a hyperbolic helix are given by: Chapter 4

κ = 1 2 cosh 2 t , τ = 1 2 cosh 2 t .
Odontological information along conical surfaces

Introduction

Over the past thirty years, the field of medicine has been using different means of medical imaging and these are in constant evolution. In particular, the odontological area is undergoing technological changes. With the discovery of X-rays in 1895, a significant advance in the diagnostic area was obtained; images of internal structures not visible clinically by the human eye were built.

In spite of its great applicability, there is still a limitation with respect to the resulting image and the interpretation of it. Until the 1970's when the diagnostic technique of X-ray computed tomography (CT) was introduced by Godfrey Newbold Hounsfield and Allan MacLeod Cormack, the diagnoses were obtained only in the form of twodimensional radiographic images.

Since the considered structures are three-dimensional (teeth, maxillas, craniums, etc.) this has limited health professionals, including dentists, because the information of width, height and depth of the sagittal, coronal and axial planes is lost in the process of obtaining the CT scan. Taking into consideration that if the radiographic images are flat and based on two of these planes, there is always a third plane of the space sacrificed and not possible to evaluate.

Thanks to the incorporation of tomography and its applicability in dentistry in the 1990's, it is possible to obtain 3D images, based on the three planes of the space (sagittal, coronal and axial) for simultaneously visualizing the structures and performing a multi-planar reconstruction. This has allowed producing images with greater specificity and sensitivity, remarkable properties achieved without super-positioning structures.

Specifically for jaws and teeth, the less explored plane corresponds to the depth or the anterior-posterior relationship of the object of study (tooth) and the adjacent anatomical structures. This relationship is very important to perform an exploration or diagnosis, to propose a treatment plan and of course to avoid complications.

In addition, complete visualization of dental structures and their relationship with adjacent tissues allows exploration beyond the sole purpose of a patient's consultation and highly increases the probability of discovering tomographic indications and specificities that until that moment were not part of the consultation.

The teeth consist of structures of different densities, including enamel, dentin and pulp, which behave differently with the passage of X-rays. Therefore, performing an a priori identification of different dental structures such as anatomical variants, presence of calcifications, dilacerations or root curves, perforations, accessory ducts, pathologies, etc., determines the success or failure of a dental treatment. In order to adequately study the anatomical structure of the teeth, curved slices that follow an area of the 3D volume are highly desirable.

A first approach in that direction was presented in the Visible Human Project [START_REF] Hersch | Exploring Anatomic Structures with EPFL's Visible Human web server[END_REF],

which offers information about human anatomy accessible to the scientific community.

It is a database that consists of 1871 horizontal plane sections from which "photos"of oblique slices -that can be extracted from the data volume-can be reconstructed. This process is mathematically well known: it involves the trilinear interpolation [START_REF] Buss | 3D Computer Graphics: A Mathematical Introduction with OpenGL[END_REF]. In fact, the Visible Human portal offers a computational tool for extracting information of organs along slices, which may be used to reconstruct three-dimensional virtual models of them.

A more interesting problem is the information extraction along curved slices, for instance, along surfaces of interest (in order to determine calcifications or other malformations), or along a jaw bone (with the purpose of making visible the information of several contiguous dental pieces). The deployment of such information is potentially useful in surgical planning. In general, when a surface is displayed in a plane, some of the areas of the slice need to be stretched and this generates a distortion problem: the shape and/or size of the original organ, along a particular section, can differ from the shape and/or size when it is displayed in a plane screen. It is the same problem of deformation that arises when building earth maps.

Particularly, in the medical area, the problem of flattening a surface without stretching was considered by Saroul in his PhD thesis [START_REF] Saroul | Surface extraction and flattening for anatomical visualization[END_REF]. Saroul minimizes the deformation in an area specified by the user at the expense of other areas (which might be less relevant or interesting) where the deformation is not controlled. However, if we consider developable surfaces in this context, the problem of metric distortion can be solved.

Specifically, developable surfaces are a special case of ruled surfaces, which can be unfolded or developed onto a plane without deformation, i.e., without stretching or tearing.

Mathematically speaking, developable surfaces are characterized by the property of possessing the same tangent plane at all points of given ruling. Due to the property to be isometrically mapped into the plane, these surfaces are interesting for visualization purposes. Developable surfaces can be classified into conical surfaces, cylindrical surfaces and tangent surfaces, which have been widely studied in many books [START_REF] Carmo | Differential Geometry of Curves and Surfaces[END_REF][START_REF] Struik | Lectures on Classical Differential Geometry[END_REF][START_REF] Kreyszig | Differential Geometry[END_REF][START_REF] Faux | Computational Geometry for Design and Manufacture[END_REF]. For the applications of this type of surfaces in areas such as engineering and manufacturing the scientific community has been studying them from different perspectives.

A way to describe these surfaces is to use the concept of duality between points and planes in 3D projective space. It means that a developable surface can be seen as the envelope of a one-parametric family of tangent planes see, e.g., [START_REF] Booduluri | Design of developable surfaces using duality between plane and point geometries[END_REF]. This is called the dual representation. In [START_REF] Booduluri | Design of developable surfaces using duality between plane and point geometries[END_REF], a developable surface is designed using control planes with Bézier and B-splines bases. In [START_REF] Pottmann | Developable rational Bézier and B-spline surfaces[END_REF][START_REF] Aumann | Interpolation with developable Bézier patches[END_REF], the authors have combined the dual representation of a developable surface with nonuniform rational B-splines (NURBS) to construct developable NURBS surfaces. In [START_REF] Peternell | Developable surface fitting to point clouds[END_REF], Peternell uses Laguerre geometry to represent developable surfaces as a one-parameter family of planes on the Blaschke model. All these methods define the surface as a NURBS family of planes and consequently as a curve in the dual projective space.

Another approach for dealing with developable surfaces is to use them for interpolating certain constraints. In [START_REF] Aumann | Interpolation with developable Bézier patches[END_REF][START_REF] Lang | Developable (1, n)-Bézier surfaces[END_REF] the authors deal with developable Bézier patches interpolating two 3D boundary curves. Chu and Séquin simplified this method. They showed in [START_REF] Chu | Developable Bézier patches: properties and design[END_REF] that given one 3D boundary curve, five degrees of freedom are available for a second boundary curve of the same degree to completely determine quadratic and cubic developable Bézier patches in 3D. According to this property, Aumann proposes in [START_REF] Aumann | A simple algorithm for designing developable Bézier surfaces[END_REF],

an algorithm to generate developable Bézier patches through a Bézier curve. Fernández-Jambrina in [START_REF] Fernández-Jambrina | B-spline control nets for developable surfaces[END_REF] generalizes that algorithm to B-spline surfaces of arbitrary degree and arbitrary number of pieces. More recently, Bo and Wang present in [START_REF] Bo | Geodesic-controlled developable surfaces for modeling paper bending[END_REF] a method for constructing developable surfaces containing a prescribed 3D curve as geodesic.

On the other hand, in the articles [START_REF] Leopoldseder | Approximation of developable surfaces with cone spline surfaces[END_REF][START_REF] Sun | A technique for constructing developable surfaces[END_REF] the authors construct global developable surfaces by joining several pieces of simple developable surfaces. In particular, Leopoldseder and Pottmann in [START_REF] Leopoldseder | Approximation of developable surfaces with cone spline surfaces[END_REF] based on an article by Fuhs and Stachel [START_REF] Fuhs | Circular pipe connections[END_REF] about circular pipeconnections, approximate general developable surfaces by circular cone segments.

The use of developable surfaces in the medical context has been studied in several articles. Figuereido and Hersch [START_REF] Hersch | Parallel unfolding visualization of curved surfaces extracted from large three-dimensional volumes[END_REF] propose the extraction of information about 3D volumes contained in cylinders built on plane curves. The aforementioned has been generalized

to developable surfaces by Paluszny [START_REF] Paluszny | Between developable surfaces and circular cone splines: curved slices of 3D volumes[END_REF] whose construction, from the computer aided geometric design point of view, has first been studied by Aumann [START_REF] Aumann | A simple algorithm for designing developable Bézier surfaces[END_REF]. Paluszny [START_REF] Paluszny | Between developable surfaces and circular cone splines: curved slices of 3D volumes[END_REF] constructs developable surfaces along a Catmull-Rom-Overhauser curve from [START_REF] Overhauser | Analytic definition of curves and surfaces by parabolic blending[END_REF], cutting the sixteen teeth of the upper human jaw bone. The main disadvantage of these methods is that the resulting surfaces are not easy to adjust to a prescribed curved slice area: a certain amount of manual tweaking is required. It is not clear how to make this process fully automatic in an elegant and numerically stable fashion. Besides, the flat presentation of the developable surface requires a numerical method. For solving these difficulties, we applied in [START_REF] González | Design of C 2 Algebraic-Trigonometric Pythagorean Hodograph splines with shape parameters[END_REF] the method of circular cone splines presented by Leopoldseder and Pottmann in [74]. This construction of conical surfaces is used to model curved slices containing information about some dental pieces of a human jaw bone. The resulting surface crosses the maxilla (upper or lower) and contains the main axis of each dental piece, that is, a straight line that passes through the root and the "center"of the dental crown, which permits to visualize the teeth simultaneously. This method has the advantage that its development can be explicitly calculated without the need to use numerical integration. The conical patches allow for the inspection and measuring of clinically relevant features, which are textured on flattened curved slices.

In this way the difficulty to calculate the development of general developable surfaces generated by Paluszny in [START_REF] Paluszny | Between developable surfaces and circular cone splines: curved slices of 3D volumes[END_REF] and Figuereido and Hersch in [START_REF] Hersch | Parallel unfolding visualization of curved surfaces extracted from large three-dimensional volumes[END_REF] is solved. However, this method of approximation with circular cones has some limitations, such as the possible presence of cone vertices within the curved slice or the appearance of not faithful cone segments to the anatomical structure of the teeth.

To overcome these limitations, we study the design of cone segments along spatial C 1 ATPH curves. These curves allow to solve efficiently the C 1 Hermite problem, as was explained in chapter 3. Besides, the spatial C 1 ATPH interpolants exhibit smoother and more pleasing curvature and torsion profiles compared with those of spatial C 1 PH quintic curves and the cubic B-splines. The cone splines are constructed by joining segments of cones with tangent continuity along given generators (main axis of each tooth) as in the case of circular cone splines. However, for our construction, cone vertices can be excluded of the region of interest for the visualization and the resulting surface has no sharp edges, which is not possible to achieve with the circular cones method described above.

We apply the construction of cone splines from C 1 ATPH curves for modeling a curved slice containing information about all dental pieces of a lower maxilla. These surfaces are compared with those constructed from circular cone segments.

The remainder of this chapter is organized as follows. In Section 4.2, we explain the procedure to deal with tomographic volumes for extracting information of interest. Particularly, within odontological volumes, we define for each dental piece a contact element (e, τ ), which denotes a line and a plane, respectively. These contact elements display clinical information of each tooth. From this, we construct conical surfaces that interpolate sequences of contact elements. These surfaces are related to the choice of a spatial curve. We consider biarc splines and spatial C 1 ATPH splines. The first technique to interpolate consecutive contact elements by circular cones is presented in Section 4.3.

We explain the general algorithm to construct the interpolating cone splines from spatial C 1 ATPH segment curves in Section 4.4, where several computed examples using this construction are compared with the aforementioned method.

Odontological volume

In this section we explain the procedure to deal with odontological volumes, with the purpose of extracting from these the initial information (main axis of the tooth) of each dental piece. This information is necessary for constructing the curved slices along jaw bones within volumes. We start by introducing some terminology about medical imaging.

Digital imaging for the biomedical sciences is commonly referred to as medical imaging. The German physicist W. C. Röntgen produced in 1895 the electromagnetic radiation in the wavelengths, that correspond to X-rays. This discovery allowed medical imaging to become an important tool in medicine, since these can be used in diag-

Object of study

Teeth and jaws are an essential part of the functions of feeding and phonation. Each dental piece is constituted by a soft connective tissue specialized in irrigation and innervation known as the pulp, which is covered by calcified tissues such as dentin, enamel and cement, which are differentiated by their level of mineralization or density. In addition each tooth is attached to the bone through the periodontium or the supporting tissues.

In both jaws (upper and lower), we can find 16 teeth, divided into 4 groups depending on the function that they perform in the masticatory system: incisors, canines, premolars and molars. In the upper jaw, we find anatomical structures such as the nostrils, the floor of the nostrils, the nasopalatine foramen, the intermaxillary suture, the maxillary sinus, the maxillary sinus walls, the maxillary tuberosity, among others. At the mandibular level, the relevant structures are the lingual foramen, the dental foramen, the inferior dental nerve canal, the external oblique line and the basal ridge. Therefore, for the radiographic and tomographic evaluation of the jaws and dental pieces, it is extremely important to know the anatomy of each one of them. Since the radiographies are twodimensional images representing three-dimensional structures, to evaluate its anatomy and variants is difficult for the dentists. So we propose to construct developable curved slices within computed tomography (CT) volumes. For this, we define some elements, which are denominated contact elements and are explained in the next subsection.

Contact element for a dental piece

For each tooth we determine a contact element (e, τ ) (e denotes a line and τ a plane) displaying the information of clinical interest. Given a segmentation of a tooth, the choice of the plane τ depends on the position of the roots that need to be visualized.

The segmentation is manual, it means, for each dental piece we determine the DICOM files that contain information of the tooth structure and we select a set of points on each image. From this dataset, we calculate the best plane that allows to visualize the piece completely. Then we texture it with the corresponding information extracted from the odontological DICOM files using a method of texturization.

The straight line e is a line segment that passes through the root and the "center" of the dental crown1 and is contained in the plane τ .

Curved slices constructed with segments of cones

In this section we construct conical surfaces that interpolate a given sequence of n contact elements {(e k , τ k )} k=1,...,n . These surfaces might be tied up with the choice of a 3D curve. Specifically, we consider biarc splines and spatial C 1 ATPH splines. The first technique deals with making curved slices with circular cone segments. The second one is about the construction of cone segments along a spatial C 1 ATPH curve. In both cases, the resulting surface is a cone spline such that each pair of cones will join with tangent continuity along a common generator. The two techniques implemented allow to construct curved slices with cone segments to interpolate a given sequence of contact elements within a CT volume. Particularly, we apply both constructions for modeling curved slices that allow to visualize relevant information about all the dental pieces of a lower jaw bone. We will explain both methods and their application in the next sections.

Circular cones with prescribed contacts

Leopoldseder and Pottmann presented in [START_REF] Leopoldseder | Approximation of developable surfaces with cone spline surfaces[END_REF] an algorithm for the approximation of a given developable surface by a cone spline surface and it depends on an adecuate choice of generators on this surface. The authors leave the problem of finding the best choice of generators as an open question. Our problem is different, we have "in principle" a fixed sequence of pairs of line/plane as given contact elements. We use the technique from [START_REF] Leopoldseder | Approximation of developable surfaces with cone spline surfaces[END_REF] to build a cone spline with these contact elements. The algorithm works as follows. The sphere Σ is determined from the two consecutive contact elements (e i , τ i ). If m 1 and m 2 are points of the generators e 1 and e 2 , the point m is the intersection of the normal planes

γ i : (x -m i ) • d i = 0, i = 1, 2 (4.1)
with the bisector plane of the two tangent planes

σ : x • (n 1 -n 2 ) -m 1 • n 1 + m 2 • n 1 = 0, (4.2)
where n i and d i are vectors of the orthonormal basis for the plane τ i with i = 1, 2, respectively.

Each one of the cones ∆ 1 , ∆ 2 , will touch the sphere Σ along a circle, c 1 and c 2 respectively, and these circles will be tangentially connected at the point c. The circles c 1 and c 2 allow us to construct a biarc which connects, with tangential continuity in c, a segment of the circle c 1 with a segment of the circle c 2 (see Figure 4.9). This joins generator e 1 with generator e 2 . For a standard rational Bézier representation for the biarc, the generator e i , such that the vector d i is the vector tangent to the spherical biarc at

= (λ 1 + λ 2 ) 2 . (4.3) This is equivalent to (a 2 -a 1 ) 2 -2λ 1 (a 2 -a 1 ) • d 1 -2λ 2 (a 2 -a 1 ) • d 2 + 2λ 1 λ 2 (e 1 • e 2 -1) = 0. (4.4)
Thus, if we choose λ 1 using the equation (4.4) we can calculate λ 2 . For the construction of the biarc, the contact point c is given by:

c = λ 2 b 1 + λ 1 b 2 λ 1 + λ 2 . ( 4.5)
In general, given a set de n + 1 control points b 0 , . . . b n , each one associated with a scalar ω i called weight, a degree-n rational Bézier curve is defined by:

c(t) = ω 0 B n 0 (t)b 0 + . . . + ω n B n n (t)b n ω 0 B n 0 (t) + . . . + ω n B n n (t)b n
, where B n i (t) are well known Bernstein polynomials. So, setting the weights at the end points of the two arcs to 1, we can express the Bézier rational quadratic form of each circular arc as follows:

c 1 (t) = a 1 (1 -t) 2 + ω 11 b 1 2t(1 -t) + ct 2 (1 -t) 2 + ω 11 2t(1 -t) + t 2 (4.6) c 2 (t) = c(1 -t) 2 + ω 12 b 2 2t(1 -t) + a 2 t 2 (1 -t) 2 + ω 12 2t(1 -t) + t 2 , (4.7)
where the weights ω 1i , associated to the internal control points, are given by:

ω 1i = |(b i -a i )(c -a i )| b i -a i c -a i . (4.8)
If λ i is positive, the arc with control points a i , b i , c and weights ω 1i > 0 are used. On the contrary, the complementary arc and a negative weight ω 1i may used.

Finally, once found the biarc as was described above, we can compute the vertices v i of the two cone segments. These are calculated as the intersection of the tangent plane to the sphere Σ at point c with generators e i . The axes of the cone pair are the lines that pass through v i and center m of Σ, respectively.

Numerical examples

Given the two contact elements in Table 4.1, they are interpolated with two circular cone segments, see In order to calculate a cone pair within the one parameter set of solutions, we may choose λ 1 and compute λ 2 using the equation (4.4). For some choices of the parameter λ 1 , it may happen that one or both vertices of the pair of cones lie inside of the region in which it is important to visualize information (see As an example we consider a sequence of 5 contact elements {(e i , τ i )} i=1,...,5 given in These situations constitute a limitation of this method. Sometimes it is not possible to avoid the presence of cone vertices within the approximating curved slice. 4.2.

Dental information along circular cone splines

In spite of the above shortcomings they might be useful to visualize information inside the CT volume. Particularly, we apply the construction of circular cones from biarcs given in the previous section, to modeling a curved slice that contains information about the internal structure of a sequence of neighboring teeth in a jaw bone.

By choosing the planes and lines found for each tooth as contact elements, we construct a curved slice built along segments of circular cones. Specifically, this process involves the following steps:

• Given a pair of contiguous contact elements (e i , τ i ), i = 1, 2 corresponding to two consecutive teeth, we calculate the center m of the sphere Σ (on which we will construct the biarc), as the intersection of three planes whose normal vectors are d 1 , d 2 and n 1 -n 2 , according to equations equations (4.1) and (4.2), where the point m i ∈ e i is the midpoint of the tooth (this is chosen manually) and d i , p i , n i are vectors of the orthonormal basis for the planes τ i with i = 1, 2, respectively.

• We interpolate the given generators e 1 and e 2 by a pair of cones ∆ 1 and ∆ 2 whose axes intersect at the point m. The sphere Σ touches each one of these cones along a circle. We denote them by c 1 and c 2 , respectively.

• The circles c 1 and c 2 generate the spherical biarc that connects the generators e i .

Its endpoints are a i = e i ∩ Σ, and the tangent vectors are p i , with i = 1, 2.

• The inner Bézier points are given by b 1 = a 1 + λ 1 d 1 , c and b 2 = a 2 -λ 2 d 2 , and can be calculated using the equations (4.4) and (4.5), once the parameter λ 1 is chosen freely. which were built for the same pair of contiguous contact elements.

Figure 4.17 illustrates a curved slice constructed with a sequence of segments of circular cones, manually choosing parameter values λ 1 to allow a good overall view of the teeth of the lower jaw bone. Although we can select each segment of cone within one family of circular cones, sometimes occur problems to completely visualize the dental piece under study. For example, in Figure 4.17 note that for some cones its vertex lies inside of the region of interest and others present a sharp edge in the common generator. This is a limitation of the technique, because it is not possible to completely visualize the dental pieces. These issues can be solved applying the technique to construct cone pairs from spatial ATPH curves, which is explained in the next subsection. More precisely, this technique allows to choose the vertices of the cones according to medical visualization needs and avoids the presence of cones with shape edges and s-shaped.

Spatial ATPH curves for designing a G 1 cone pair

Given a sequence of n contact elements {(e k , τ k )} k=1,...,n , we want to construct cone segments that interpolate it with tangent continuity which does not have the problems of circular cone splines. Each cone segment will be defined by a space ATPH curve r(t)

-the directrix -and a fixed point v -the vertex -, lying outside of a region of interest. The steps to construct two cone segments that interpolate two contact elements are given in the following subsection.

General algorithm

Let (e i , τ i ), i = 1, 2 be two consecutive contact elements. The idea is to interpolate Let a i be a point on the line e i and let e i , d i , n i be the vectors constituting the orthonormal basis for each plane τ i , such that n i = e i × d i is the unit normal vector of τ i with i = 1, 2, respectively, see Each ATPH curve can be expressed as where B 5 k (t) is the normalized B-basis of space U 5 (see the equation (2.3)) and p k,i , k = 0, . . . , 5 denote its control points, which are defined by p 1,i = p 0,i + n 0 16s where the abbreviations s 1 , c 2 , n 0 , n 2 are defined in (2.4) and (2.5). The values A 0 , A 1 , A 2 are given by equation (3.20). By construction, we have that p 0,1 = a 1 , p 5,1 = p 0,2 = a c and p 5,2 = a 2 .

(e i , τ i ), i = 1,
r i (t) = 5 k=0 p k,i B 5 k (t), t ∈ [0, α i ], i = 1, 2, (4.9)
Each ATPH curve r i (t) depends on three free parameters (θ i , β i , α i ), i = 1, 2. In order to choose optimally these parameters, we apply a minimization process over one of the following functionals

F1 = | r | dt κ 2 | r | dt, F2 = | r | dt 3 κ 2 | r | dt, F3 = | r | dt 3 κ 2 + τ 2 1 | r | dt, (4.11)
where κ and τ denote the curvature and the torsion of r(t), respectively.

For the sake of clarity, we summarize the steps of the algorithm to construct a G 1 cone pair through C 1 ATPH curves. Each pair of cones will interpolate two given contact elements.

Algorithm

Input: (e i , τ i ), a i , d i , for i = 1, 2. . 3. Compute the three free parameters (θ i , β i , α i ), i = 1, 2 of the ATPH curve r i from (4.9) minimizing a functional given in (4.11).

Output: Control points p k,i , k = 0, . . . , 5, i = 1, 2 of the ATPH curves r 1 , r 2 .

Numerical examples

We consider the data given in Table 4.1. The quantities listed are: the orthonormal base for the two contact elements (e i , τ i ), i = 1, 2. Namely e i , d i , n i for i = 1, 2 and the end points a i lying on the lines e i for i = 1, 2, respectively, which are illustrated in and d c = (0.9386, 0.2636, 0.0931). Finally we compute the three parameters (θ i , β i , α i ), i = 1, 2 for each segment of ATPH curve, minimizing each functional given in (4.11). In Table 4.4 we present the corresponding values of the functionals F1 , F2 , and F3 for the ATPH curves r 1 , r 2 . In this case, the value of F1 is significantly smaller than the values of F2 and F3 . Therefore, we construct the pair of cones through the ATPH curves r 1 , r 2 with the parameters (θ i , β i , α i ), i = 1, 2 given by the minimization of F1 .

In Figure 4.18 we display the points a 1 , a 2 and a c and the contact elements (e and a biarc (see Figure 4.10) are very similar. This is because the spatial ATPH curve closely approximates 3D parametric curves, as circular helices (see Section 3.6.3). On the other hand the data given in Table 4.2 is interpolated with 8 cone segments, see Figure 4.20, which are calculated following the method explained above. This method allows to construct the interpolating cone spline free of singularities inside the region of interest. The resulting surface does not present sharp edges due to the geometric properties of the functional used to calculate the cone segments. In this case, we have used the functional F3 (see equation (4.11)), which minimizes the rate of change of the curvature and the torsion, for each ATPH curve segment that composes the cone spline.

In general the functionals F1 , F2 , F3 (see equation (4.11)) can be used according to the geometric needs of the practical application under consideration.

Even if the circular cones are simpler surfaces than the general cones described in Section 4.4, the construction of the latter ones allows to avoid the presence of singular points in the region of interest. Therefore, for applications in odontological area, using general cones designed from C 1 ATPH curves may be produce more faithful results.

Cone splines from 3D ATPH curves to visualize odontological information

We construct cone slices through 3D ATPH curves, that cross the lower maxilla of a computed tomography (CT) volume, in order to present simultaneously relevant information of all dental pieces. We use the same CT volume used in the section 4.3.3.

For simplicity, let us look at two consecutive contact elements (e i , τ i ), i = 1, 2, for which we compute the corresponding cone pair. The steps of this construction are summarized as follows:

• We consider for each dental piece one contact element. In particular, for the examples in this section we use the sequence of contact elements illustrated in • To construct each cone pair, we choose two points v 1 and v 2 , on the bounding generators e 1 , e 2 , lying outside of the region of interest. These points are the vertices of the cone pair.

• The middle point between v 1 and v 2 will be the contact point a c .

• Let a i be a point on the line e i and e i , d i and n i the vectors forming the orthonormal basis for each plane τ i with i = 1, 2. We calculate d c as the mean vector of d 1 and d 2 .

• We compute C 1 Hermite ATPH curves r 1 and r 2 interpolating the end points a 1 , a c and a c , a 2 with associated end derivatives d 1 , d c and d c , d 2 , respectively.

Each interpolating ATPH curve depends on three parameters (θ, β, α). We minimize a functional F1 , F2 or F3 to determine these free angular parameters. to select the vertex of each cone that composes the developable surface. Additionally, this type of surface has the property of allowing for its planar unfolding, i.e., it may be mapped on the plane without stretching. This process is referred to as isometric Sometimes, the constructed developable surface does not show sufficient information to realize a diagnosis or a treatment plan. For that reason, we use the technique of "offset surfaces". This method permits to construct "parallel" surfaces, i.e, surfaces following a direction into the volume, previously fixed. See, e.g., Figure 4.25.

In addition, the "offset surfaces" can fill up the whole dental volume, which allows to obtain and/or to extend information regarding the dental pieces. In Figure 4.26 can be Chapter 5

Analysis and conclusions

In this thesis, we have solved 2D and 3D interpolation problems with ATPH curves. The solution to the 3D interpolation problem has been used to construct piecewise conical patches, with the purpose of visualizing odontological information contained in 3D CT volumes. The results obtained in each chapter are summarized as follows:

• We presented in chapter 2 a new solution to the problem of interpolating a planar sequence of points by a C 2 continuous spline curve. We obtain the resulting C 2 ATPH splines, which depend on several shape parameters, by solving a nonlinear quadratic system of equations in complex unknowns. Thanks to the shape parameters, we obtain a diversity of ATPH interpolating splines for each given sequence of points. Moreover, these free parameters can be conveniently used to either improve the curvature distribution of the ATPH curve or to adjust the shape according to the application. We compared these splines with the well-known polynomial PH quintic splines and the cubic B-splines. We observe that if the sequences of free parameters are suitably selected, the ATPH interpolants exhibit smoother and more pleasing curvature profiles than those of the PH quintic splines and cubic B-splines. In order to objectively select the sequence of parameters we employed minimization methods of standard pre-established fairness functionals.

• In chapter 3, we have studied spatial ATPH curves, which are characterized in terms of quaternions. After introducing this representation, they are employed to solve the first-order Hermite interpolation problem. This problem is reduced to solving a system of three quadratic equations in three quaternion unknowns.

Three angular variables remain free on solving this system, for which we propose several integral criteria for automatically determining these free parameters. Their performance has been examined through several examples. Furthermore, we have compared the C 1 interpolating ATPH curves with the corresponding PH quintic interpolant. Finally, we successfully approximated known 3D parametric curves with ATPH curves by calculating the three free shape parameters through integral criteria.

• In chapter 4, we introduced two techniques for extracting information along curved slices (which can be flattened without deformation) from dental data volumes.

Under visual inspection both techniques yield results of comparable quality. The main advantage of the approach with cone splines from spatial ATPH curves as compared to circular cone splines is the possibility to avoid the presence of cone vertices within the curved slice, which would be singular points of the cone spline.

In the case of circular cones, once the biarc is found, the vertices v 1 and v 2 of the two cone segments are completely determined, which implies a limitation for this technique.

The main similarities and differences between circular cone splines adapted from [START_REF] Leopoldseder | Approximation of developable surfaces with cone spline surfaces[END_REF] and our new technique of cone splines from ATPH curves are:

-Both cone splines interpolate a sequence of contact elements (plane/line) and each segment can be adjusted independently while preserving the tangent continuity property at common generators. In other words: each method provides local control.

-Both surfaces can be unfolded isometrically onto the plane.

-Cone splines from ATPH curves are faithful to the tooth sequence, because of the vertices can be chosen according to the medical visualization needs. On the contrary, the circular cone splines technique does not have this flexibility because given the biarc, the developable surface is completely determined and it might present sharp edges or s-shapes.

Within the odontological field we consider the extraction of texture information along curved slices of sets of teeth, their internal properties and their placement in the maxillofacial bone structure. Our new technique to construct cone splines from ATPH curves allows for flat display and exact measurements in 2D, because the relationship between the 3D slice and its 2D display is given by an isometry.

As shown in Figure 4.23, the method also allows for the full maxillofacial display in 3D and 2D which might be informative in aspects such as a thorough inspection of dental pieces together with relevant anatomical information such as bone decalcification, the height of the dental floor and the position of the pieces with respect to the palatal vault.

Other possible fields of application of our technique include the construction of curved sections of veins (to study valve function) and arteries (for detecting calcifications).

  Thèse de doctorat Pour obtenir le grade de Docteur de l'Université de VALENCIENNES ET DU HAINAUT-CAMBR ÉSIS Spécialité : MATH ÉMATIQUES APPLIQU ÉES Présentée et soutenue par Cindy GONZ ÁLEZ. Le 25/01/2018, á Valenciennes École doctorale : Sciences Pour l'Ingénieur (SPI) Laboratoire : Laboratoire de Mathématiques et ses Applications de Valenciennes (LAMAV)Algebraic-Trigonometric Pythagorean Hodograph curves for solving planar and spatial interpolation problems and their use for visualizing dental information within 3D tomographic volumes Les courbes Algébriques Trigonométriques à Hodographe Pythagorien pour résoudre des problèmes d'interpolation deux et trois-dimensionnels et leur utilisation pour visualiser les informations dentaires dans des volumes tomographiques 3D JURY Président du jury • ALI-MEHMETI, Felix. Professeur, Université de Valenciennes et du Hainaut-Cambrésis. Résumé Les problèmes d'interpolation ont été largement étudiés dans la Conception Géométrique Assistée par Ordinateur (CGAO). Ces problèmes consistent en la construction de courbes et de surfaces qui passent exactement par un ensemble de données, tel que, des nuages de points, des tangentes, des courbures, des lignes/plans, etc. En général, ces courbes et surfaces sont représentées sous une forme paramétrique. Cette représentation est indépendante du système de coordonnées, elle s'adapte bien aux transformations géométriques et les propriétés géométriques différentielles des courbes et des surfaces sont invariantes sous réparametrization. Dans ce cadre, l'objectif principal de cette thèse est de présenter des méthodes d'interpolation de données 2D et 3D au moyen de courbes Algébriques Trigonométriques à Hodographe Pythagorien (abréviation: courbes ATPH). Celles-ci sont des courbes paramétriques définies dans un espace algébriquetrigonométrique mixte, dont l'hodographe satisfait une condition pythagorienne. Cette représentation permet de calculer analytiquement l'abscisse curviligne aussi bien que la paramétrisation rationnelle-trigonométrique des courbes offset. Ces propriétés sont utilisables pour la conception de modèles géométriques dans de nombreuses applications. En particulier, nous nous intéressons à la modélisation géométrique d'objets odontologiques. À cette fin, nous utilisons les courbes spatiales ATPH pour la construction de surfaces développables dans des volumes 3D odontologiques. Dans l'introduction nous donnons un survol sur comment des problèmes d'interpolation similaires ont été abordés par la communauté scientifique. Ensuite, dans le chapitre 2 nous considérons la construction de courbes planes ATPH avec continuité C 2 qui interpolent une séquence ordonnée de points. Nous employons deux méthodes pour résoudre ce problème. Nous calculons d'abord toutes les solutions avec une méthode d'homotopie. Empiriquement, il a été observé qu'il n'y a qu'une seule "bonne" solution. Par conséquent, la méthode d'itération de Newton-Raphson a été utilisée pour calculer directement cette "bonne" solution. Notez que les courbes C 2 ATPH dépendent de plusieurs paramétres libres, qui permettent d'obtenir une diversité d'interpolants.

3 . 8

 38 ATPH spline curves correspond to minimizing the functionals IES and IMV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.18 Comparison of the C 2 ATPH spline curve, the C 2 PH quintic spline curve, the C 2 cubic B-spline and the unit circle . . . . . . . . . . . . . . . . . . . 2.19 Comparison of the C 2 ATPH spline curve, the C 2 PH quintic spline curve, the C 2 cubic B-spline and the limaçon . . . . . . . . . . . . . . . . . . . . 2.20 Sequences of 10, 15, respectively 20 interpolating points on the unit circle 2.21 Sequences of 10, 15, respectively 20 interpolating points on the limaçon .3.1 The spatial ATPH interpolants to the end points p i , p f and end derivatives d i , d f . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Curvature profiles of the spatial ATPH interpolants . . . . . . . . . . . . 3.3 Torsion profiles of the spatial ATPH interpolants . . . . . . . . . . . . . . 3.4 The spatial ATPH interpolants to the end points p i , p f and end derivatives d i = d f . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 The spatial ATPH curves and PH quintic curve obtained by minimizing F1 3.6 Curvature profiles of the spatial ATPH curves and the corresponding PH quintic curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Torsion profiles of the spatial ATPH curves and the corresponding PH quintic curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Figures ix The spatial ATPH curves and PH quintic curve obtained by minimizing F2 3.9 Curvature profiles of the spatial ATPH curves and the corresponding PH quintic curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Torsion profiles of the spatial ATPH curves and the corresponding PH quintic curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 The spatial ATPH curves and PH quintic curve obtained by minimizing F3 3.12 Curvature profiles of the spatial ATPH curves and the corresponding PH quintic curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Torsion profiles of the spatial ATPH curves and the corresponding PH quintic curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Comparison of the C 1 ATPH curve, the C 1 PH quintic curve and the circular helix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.15 Comparison of the C 1 ATPH curve, the C 1 PH quintic curve and the hyperbolic helix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 MRI slice of the body of a man and CT slice of the head of a woman . . . 4.2 Four axial slices of an odontological volume . . . . . . . . . . . . . . . . . 4.3 Codification of colors of a medical image . . . . . . . . . . . . . . . . . . . 4.4 Contact element (e, τ ) of one tooth dental piece . . . . . . . . . . . . . . . 4.5 3D view of the planes τ i for i = 1, . . . , 16 . . . . . . . . . . . . . . . . . . . 4.6 Two contact elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Circular cones with a common generator and their axes . . . . . . . . . . 4.8 Sphere Σ inscribed in two circular cones . . . . . . . . . . . . . . . . . . . 4.9 Biarc with its control polygon . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Different 3D views of circular cone segments that interpolate two contact elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .4.11 Examples of pair cone segments . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Example of a cone pair where one vertex lies inside of the region of interest 4.13 Example of a cone pair with sharp edges and a cone pair with s-shaped . 4.14 Circular cones interpolating the 5 contact elements . . . . . . . . . . . . . 4.15 Cone spline in 3D composed of four segments of circular cones for one tooth of the lower jaw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Examples of pairs of circular cone segments . . . . . . . . . . . . . . . . . 4.17 Curved slice with segments of cones for some dental pieces of the lower jaw bone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.18 Pair of C 1 ATPH curves r 1 , r 2 that interpolate the end points a 1 , a c , a 2 associated to the end derivatives d 1 , d c , d 2 . . . . . . . . . . . . . . . . . . 4.19 Cone pair defined by the vertices v 1 , v 2 and a pair of C 1 ATPH curves r 1 , r 2 that interpolate the points a 1 , a c , a 2 . . . . . . . . . . . . . . . . . . 4.20 Cone spline designed from C 1 Hermite ATPH curves . . . . . . . . . . . . 4.21 Curved slice for visualizing a tooth of the lower jaw . . . . . . . . . . . . . 4.22 Curved slice with segments of cone splines for all dental pieces of the lower jaw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.23 3D views of a curved slice showing the information related to the lower jaw105 4.24 View of development without distortion of curved slice . . . . . . . . . . . 4.25 Example of "offset" surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . 4.26 Two texturized "offset" spline surfaces . . . . . . . . . . . . . . . . . . . . List of Figures x 4.27 Example of a tooth using the techniques of cone splines from biarcs and 3D ATPH curves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

  ATPH spline and C 2 PH quintic spline interpolating the same sequences of points on the unit circle . . . . . . . . . 2.17 Approximation errors for C 2 ATPH spline and C 2 PH quintic spline interpolating the same sequences of points on the limaçon . . . . . . . . . . 3.1 Hermite data for the examples illustrated in Figures 3.1 and 3.4 . . . . . . 3.2 Various choices of the angular parameters θ, β for the ATPH interpolants illustrated in Figures 3.1 and 3.4. . . . . . . . . . . . . . . . . . . . . . . . 3.3 Derivative data for the five test curves. . . . . . . . . . . . . . . . . . . . . 3.4 Functional values for PH quintic curves. . . . . . . . . . . . . . . . . . . . 3.5 Total arc length corresponding to each functional for the PH quintic curves. 3.6 Functional values for ATPH curves. . . . . . . . . . . . . . . . . . . . . . . 3.7 Total arc length corresponding to each functional for the ATPH curves. . Chapter 1

Finally, in chapter 5

 5 we summarize and analyze the results obtained in this thesis, and give an overview on future work. The Algebraic-Trigonometric Pythagorean Hodograph (ATPH) curves, introduced in [51], are a new class of Pythagorean Hodograph (PH) curves defined over a mixed algebraic-trigonometric space. The well-known polynomial PH curves are characterized by having a closed form representation of their arc length and rational offset curves. They have been used to solve various Hermite interpolation problems, see, e.g. [24-30]. Analogously, the ATPH curves allow a closed form representation of their arc length and have rational algebraic-trigonometric offset curves. Romani et al. [51] solved the C 1 Hermite interpolation problem by ATPH curves and constructed G 2 continuous ATPH transition curves of monotone curvature between basic elements such as line segments and circles.

Figure 2 . 1 :

 21 Figure 2.1: B-basis functions B 2 i (t) i=0,1,2 , B 4 i (t) i=0,...,4 and B 5 i (t) i=0,...,5 for t ∈ [0, α] and α = π/3.

  α -t) -1) 2 (cos(α) -1) 2 , B

Figure 2 . 2 ii=0, 1 , 2 , B 4 i 5 i

 221245 1 shows these B-bases for α = π/3. Note that the parameter α presents a tension-like effect when α is increased in its interval of definition. Moreover, when α tends to π (2π in the case of U 5 ), the normalized B-bases B i=0,...,4 and B i=0,...,5 vanish, except the first and the last functions (see Figure 2.2).

Figure 2 .

 2 Figure 2.2: B-basis functions B 2 i (t) i=0,1,2 , B 4 i (t) i=0,...,4 and B 5 i (t) i=0,...,5 for t ∈ [0, α] and α → π, 2π.

  18) relates three unknowns b k-1 , b k , b k+1 to the known quantity ∆q k . Thus we have a total of N equations for N +2 unknowns b 0 , . . . , b N +1

Figure 2 . 3 :

 23 Figure 2.3: Family of C 2 ATPH splines interpolating five points. The "good" ATPH spline is displayed in red.

  need of evaluation of the entries of the Jacobian matrix of the former. The overhead consists in 690N+50 elementary operations for open and 690N elementary operations for closed ATPH splines, where N is the number of interpolation points. By elementary operations we intend additions/substractions, multiplications/divisions and evaluation of trigonometric functions. In fact, in row k of the Jacobian we have 690 elementary operations for k = 1, ...N , and 25 elementary operations for k = 0 as well as for k = N + 1 in the case of open ATPH splines. This evaluation of the entries of the Jacobian matrix must be computed only once before starting the Newton-Raphson iteration. The higher cost for ATPH splines is offset by the extra flexibility enjoyed by these splines which allow for minimization of shape functionals, such as IMV and IES, considered in section 2.5.

Figure 2 .

 2 4 shows illustrative examples of C 2 ATPH spline curves. Relevant data from these examples are summarized in Table 2.2. Initial values for the Newton-Raphson method are calculated by comparing mid-point derivatives of the ATPH curve with those of the cubic B-spline, (see subsection 2.4.3). Figures 2.5 and 2.6 show the behaviour of different choices of the sequence {α k } k=1,.

Figures 2 .

 2 7 and 2.8 show the comparison of the C 2 cubic B-spline, the C 2 PH quintic spline and the C 2 ATPH spline for open and closed curves, respectively. In these figures, as well as in many others, we use color lines to distinguish among various examples. We

Figure 2 . 4 :

 24 Figure 2.4: Top: Example of an open C 2 ATPH spline curve that interpolates 4 points. Bottom: Example of a closed C 2 ATPH spline curve that interpolates 6 points (see Table2.2).

Table 2 . 2 :

 22 Data for C 2 ATPH splines in Figure 2.4.

Figure 2 . 5 :

 25 Figure 2.5: Open C 2 ATPH splines interpolating the same data, obtained for different sequences of {α k } k=1,2,3 and their corresponding curvature profiles.

Figure 2 . 6 :

 26 Figure 2.6: Variability of a closed C 2 ATPH splines that interpolates q 0 , . . . , q 5 with respect to several choices of the sequences {α k } k=1,...,5 and their corresponding curvature profiles.

Figure 2 . 7 :

 27 Figure 2.7: Comparison of the C 2 ATPH spline, the C 2 PH quintic spline and the C 2 cubic B-spline.

Figure 2 . 8 :

 28 Figure 2.8: Comparison of the C 2 ATPH spline, the C 2 PH quintic spline and the C 2 cubic B-spline.

Figure 2 . 9 :

 29 Figure 2.9: Open C 2 ATPH spline curves for increasing number of input data points.

Figure 2 . 10 :

 210 Figure 2.10: Closed C 2 ATPH spline curves for increasing number of input data points.

Figure 2 . 11 :

 211 Figure 2.11: Examples of ATPH splines with different sequence of parameters {α k } k=1,...,3 , with α k = α and their corresponding total arc-length s T .

Theorem 2 . 2 .

 22 Let f (b) be the map C N +2 → C N +2 for open curves or C N → C N for closed curves, defined by (2.18), (2.20) or (2.18), (2.21), with Jacobian matrix J(b) given by (2.28), (2.29) or (2.28), (2.30), where b denotes b = (b 0 , . . . , b N +1 ) for open curves or b = (b 1 , . . . , b N ) for closed curves. Let b (0) ∈ D ⊆ C N +2 in the case of open curves and respectively, b (0) ∈ D ⊆ C N in the case of closed curves, denote the initial point for the Newton-Raphson method in (2.27).

  always remain inside the ball S with center b (0) and radius r, and converge to a solution b * of f (b) = 0.

  , 2.16 for closed curves. We have used these examples to illustrate the convergence that can be attained in the construction of ATPH splines. Many examples have shown that the convergence condition of Theorem 2.2 is in general satisfied after a small number of iterations for any sequence of parameters.

Figure

  Figure 2.12: Lipschitz constants K α .

Figure 2 . 14 :

 214 Figure 2.14: Open curves corresponding to arbitrary sequences of parameters.

Figure 2 . 15 :

 215 Figure 2.15: Top: closed curves corresponding to chordal sequences of parameters.Left: zoomed curves.

Figure 2 . 16 :

 216 Figure 2.16: Closed curves corresponding to arbitrary sequences of parameters.

Figure 2 . 17 :

 217 Figure 2.17: Left: The C 2 ATPH spline curves correspond to minimizing the functionals IES and IMV. Right: Curvature profiles for both curves.

α k 2 )

 2 for ATPH curves and Y k = r k,P H ( t k-1 +t k 2 ) for PH quintic curves with k = 1 . . . , N . The prediction points Ỹk are computed by using a code called distance2curve of MATLAB, which finds the closest point on a curve to any given point.

Figure 2 . 18 :

 218 Figure 2.18: Left: Comparison of the C 2 ATPH spline curve, the C 2 PH quintic spline curve, the C 2 cubic B-spline and the unit circle. Right: Curvature profiles of the four spline curves.

Figure 2 . 19 :

 219 Figure 2.19: Left: Comparison of the C 2 ATPH spline curve, the C 2 PH quintic spline curve, the C 2 cubic B-spline and the limaçon. Right: curvature profiles of the four curves.

Figure 2 .

 2 Figure 2.19 Right), specifically 1.1151 ≤ κ ≤ 2.3102 for the ATPH spline curve and 1.1156 ≤ κ ≤ 2.3094 for the limaçon. On the other hand the distributions of

Figure 2 . 20 :

 220 Figure 2.20: Sequences of 10, 15, respectively 20 interpolating points on the unit circle.

  We consider in this chapter the construction of smooth spatial, Algebraic Trigonometric Pythagorean Hodograph (ATPH) curves that interpolate given end points and end derivatives. The ATPH curves are parametric curves defined over the mixed algebraic-trigonometric space U 5 (see chapter 2), whose hodograph or derivative satisfies a Pythagorean condition. The solution of the first-order Hermite interpolation problem by planar ATPH curves has been studied by Romani et al. in[START_REF] Romani | Algebraic-Trigonometric Pythagorean-Hodograph curves and their use for Hermite interpolation[END_REF], using the complex representation of them. This characterization facilitates the construction and shape analysis of the ATPH interpolants. The obtained four solutions are analyzed and the best interpolant is identified by means of shape measures. The solution of an analogous problem by the well-known Pythagorean-Hodograph (PH) quintic curves was

C 1

 1 3 we describe the quaternion representation of spatial ATPH curves. In Section 3.4 the spatial ATPH curves are employed to solve the Hermite interpolation problem. Several illustrative examples are presented in Section 3.5. The criteria for selecting the three free parameters of ATPH are then described in Section 3.6.

+ λx k 1 2 ( 1 +. 8 )

 218 4), due to Ãi à * = ABi(AB) * = ABiB * A * = AiA * . By setting B = a + bi + cj + dk, equation (3.6) is equivalent to the system a 2 + b 2 -c 2 -d 2 = 1, 2(ad + bc) = 0, 2(bd -ac) = 0, whose real solutions are of the form (a, b, c, d) = (cos(φ), sin(φ), 0, 0). Thus, we can parametrize all real solutions to (3.6) in terms of a parameter φ as B(φ) = cos(φ) + sin(φ)i. (3.7) Then quaternion à can be written in terms of A and φ as à = -a x sin φ + a y cos φi + (a y cos φ + a z sin φ)j + (a z cos φ -a y sin φ)k and we can deduce a general solution A(φ) to equation (3.4) from the particular solution (3.5) as follows A(φ) = -sin φ + cos φi + λy cos φ + λz sin φ 1 + λx j + λz cos φ -λy sin φ 1 λx) | d |. (3Using the notation exp(φi) = cos φ + sin φi and taking an unit vector δ = d |d| in the direction of d and the unit bisector vector of δ and i, n = δ + i |δ + i| , we can write the above equation as A(φ) = | d |n exp(φi). (3.9)

1 .

 1 The interpolating ATPH curves are displayed in Figures 3.1 and 3.4, together with the corresponding PH quintic curve. Curvature and torsion profiles of the ATPH curves illustrated in Figure 3.1 are displayed in Figures 3.2 and 3.3, respectively. For these examples, we have varied the values of the parameters (θ, β, α) over the domain [0, 2π] × [0, 2π] × [0, 2π/3] (see Tables 3.1,3.2).The selection of the angular parameters (θ, β) has a high influence on the shape of the ATPH interpolants. Moreover, if we fixed these parameters and vary the value of the parameter α, we obtain different ATPH curves interpolating the same data. This can be clearly seen in Figures 3.1 and 3.4. In Figure3.4 we have drawn in blue the ATPH curve that corresponds to the choice of the parameter α = π/10 and in black for the ATPH curve corresponding to α = 2π/3. For the selected values of α given in Table3.1, these curves agree with the shortest and the longest ATPH interpolant. The respective values of the total arc-length are given in the legend (see Figure3.4). Furthermore, based on numerous examples we have observed that, as in the planar case, for increasing values of α ∈ (0, 2π/3) the curves become longer.

  Figure 3.1 p i = (0, 0, 0), p f = (1, 1, 1) α = π/4, π/3, 2π/5, π/2, 3π/5, 2π/3 d i = (-0.8, 0.3, 1.2), d f = (0.5, -1.3, -1.0) Figure 3.4 p i = (5, 0, 1), p f = (-3, -4, 1) α = π/10, π/4, π/3, 2π/5, π/2, 2π/3 d i = (25, -15, 1), d f = (25, -15, 1)

Figures 3 . 4 Figure 3

 343 Figures 3.1(a), 3.4(a) θ = 5π/4 β = 3π/2 Figures 3.1(e), 3.4(c) θ = π/2 β = 0 Figure 3.1(b) θ = 5π/8 β = 7π/4 Figure 3.1(f) θ = 3π/2 β = 0 Figures 3.1(c), 3.4(b) θ = 9π/8 β = 7π/4 Figure 3.1(g) θ = 7π/8 β = 7π/4 Figure 3.1(d) θ = 0 β = 0 Figures 3.1(h), 3.4(d) θ = π/8 β = 3π/4

Figure 3 . 1 :

 31 Figure 3.1: The spatial ATPH interpolants to the end points p i = (0, 0, 0), p f = (1, 1, 1) and end derivatives d i = (-0.8, 0.3, 1.2), d f = (0.5, -1.3, -1.0). The choice of the two parameters is given in Tables 3.1, 3.2. The corresponding PH quintic interpolant is drawn in red.

Figure 3 . 2 :

 32 Figure 3.2: Curvature profiles of the spatial ATPH interpolants illustrated in Figure 3.1. The corresponding curvature profile of the PH quintic interpolant is drawn in red.

Figure 3 . 3 :

 33 Figure 3.3: Torsion profiles of the spatial ATPH interpolants illustrated in Figure 3.1. The corresponding torsion profile of the PH quintic interpolant is drawn in red..

Figure 3 . 4 :

 34 Figure 3.4: The spatial ATPH interpolants to the end points p i = (5, 0, 1), p f = (-3, -4, 1) and end derivatives d i = d f = (25, -15, 1) obtained with the values of the two parameters is given in Tables 3.1 and 3.2. The corresponding PH quintic interpolant is drawn in red. The shortest and longest ATPH interpolant are drawn in blue and black, respectively.

1 , 1 ), excluding case 4 ,

 114 in which we change the point p f to p f = (0.15396, -0.60997, 0.40867). The interpolating C 1 ATPH Hermite curves and C 1 PH quintic Hermite curves are displayed in Figures 3.5-3.11.

case 1 case 2 case 3 case 4 case 5 STable 3 . 7 :

 537 F1 1.8375 3.1418 4.5640 1.0185 1.8714 S F2 1.8720 3.1418 4.5636 1.0216 1.9124 S F3 1.8719 3.1418 4.5647 1.0227 2.2277 Total arc length corresponding to each functional for the ATPH curves.

Figure 3 . 6 :

 36 Figure 3.6: Curvature profiles of the spatial ATPH curves (green line) and the corresponding PH quintic curve (red line) of Figure 3.5.

Figure 3 . 7 :

 37 Figure 3.7: Torsion profiles of the spatial ATPH curves (green line) and the corresponding PH quintic curve (red line) of Figure 3.5.

Figure 3 . 9 :

 39 Figure 3.9: Curvature profiles of the spatial ATPH curves (green line) and the corresponding PH quintic curve (red line) of Figure 3.8.

Figure 3 . 10 :

 310 Figure 3.10: Torsion profiles of the spatial ATPH curves (green line) and the corresponding PH quintic curve (red line) of Figure 3.8.

Figure 3 . 12 :

 312 Figure 3.12: Curvature profiles of the spatial ATPH curves (green line) and the corresponding PH quintic curve (red line) of Figure 3.11.

Figure 3 . 13 :

 313 Figure 3.13: Torsion profiles of the spatial ATPH curves (green line) and the corresponding PH quintic curve (red line) of Figure 3.11.

Figure 3 . 14

 314 Figure 3.14 Left shows the behavior of the C 1 ATPH curve, the C 1 PH quintic curve and the circular helix interpolating the end points p i = (1, 0, 0), p f = (0.7071, 0.7071, 0.7854) and the end derivatives d i = (0, 1, 1), d f = (-0.7071, 0.7071, 1). The three curves are in remarkable agreement and they are

Figure 3 . 14 :

 314 Figure 3.14: Left: Comparison of the C 1 ATPH curve, the C 1 PH quintic curve and the circular helix. Right: Curvature profiles of the three curves. Bottom: The distribution of torsion for each curve.

Figure 3 .

 3 Figure 3.15 illustrates the corresponding arc of the hyperbolic helix, the C 1 ATPH curve and the C 1 PH quintic curve, interpolating the end points p i = (1, 0, 0), p f = (1.1024, 0.4640, 0.4488) and associated end derivatives d i = (0, 1, 1), d f = (0.4640, 1.1024, 1). From their curvature and torsion distribution displayed in Figure 3.15 Right and Bottom, respectively, we can deduce that the ATPH curve achieves a very good approximation of the arc of the hyperbolic helix. This affirmation is supported by the fact that values of the approximation errors MSE and RMSE are low, namely MSE=6.7163e-9 and RMSE=8.1953e-5, while these values in the case of the PH quintic curve are 1.2581e-4 and 0.0112, respectively.

Figure 3 . 15 :

 315 Figure 3.15: Left: Comparison of the C 1 ATPH curve, the C 1 PH quintic curve and the hyperbolic helix. Right: Curvature profiles of the three curves. Bottom: The distribution of torsion for each curve.

  noses, planning treatments, procedures and clinical monitoring. Currently medical images are produced with various techniques, among them are Computed Tomography (CT) and Magnetic Resonance Imaging (MRI).

Figure 4 .

 4 1 shows two slices, of MRI and CT, see https://mri.radiology.uiowa.edu/visible_human_datasets.html for additional examples. Particularly, in this chapter we consider a tomographic volume

Figure 4 . 1 :

 41 Figure 4.1: MRI slice of the body of a man and CT slice of the head of a woman.

  Figure 4.2 shows four slices from this DICOM volume. The CTs or MRIs are composed of images, whichare files with extension ".dcm". The latter is the standard format for manipulation, storage, printing and transmission of medical imaging information. Some programs that allow to read and modify files in format DICOM are OsiriX, MicroDicom, Radiant, Open Dicom Viewer, etc. In our case, we use MATLAB to manipulate these files.

Figure 4 . 2 :

 42 Figure 4.2: Four axial (horizontal) slices of an odontological volume.

Figure 4 .

 4 3 shows an example of the assignment of colors of a tomography. The instruction dicominfo( Image1.dcm ) allows to visualize the additional information about patient and image data, such as position and orientation.

Figure 4 . 3 :

 43 Figure 4.3: Codification of colors of a medical image.

Figure 4 .

 4 4 illustrates the plane chosen for

Figure 4 . 4 :

 44 Figure 4.4: Contact element (e, τ ) of one tooth dental piece, where e, d is the orthonormal basis for τ such that n = e × d is the unit normal vector of τ and e gives the direction of the line e.

Figure 4 .

 4 5 shows all planes chosen along the lower maxilla.

Figure 4 . 5 :

 45 Figure 4.5: 3D view of the planes τ i for i = 1, . . . , 16, which have been textured with the information of each dental piece.

Figure 4 .

 4 Figure 4.6 shows two contact elements (e i , τ i ), i = 1, 2 each one with the corresponding orthonormal basis.

Figure 4 . 6 :

 46 Figure 4.6: Two contact elements (e 1 , τ 1 ), (e 2 , τ 2 ).

4. 3 . 1 . 1

 311 General algorithmGiven two consecutive contact elements (e 1 , τ 1 ) and (e 2 , τ 2 ) we want to find two circular cones ∆ 1 , ∆ 2 , with different vertices v 1 , v 2 , that have a common generator e c and the same tangent plane τ c along this generator. Each ∆ i must also contain the generator e i , and its tangent plane along e i must match up with τ i . The axes of a pair of cones in this position either intersect in a point m or are parallel. We consider the case in which the axes intersect in m (see Figure4.7). It can be shown ([START_REF] Leopoldseder | Approximation of developable surfaces with cone spline surfaces[END_REF]) that the pair of cones

Figure 4 . 7 :

 47 Figure 4.7: Circular cones with a common generator and their axes intersect at a point m.

Figure 4 . 8 :

 48 Figure 4.8: Sphere Σ inscribed in the circular cones ∆ 1 and ∆ 2 .

Figure 4 . 9 :

 49 Figure 4.9: Biarc with its control polygon.

Figure 4 .

 4 10, which are constructed following the algorithm exposed above. e 1 = (-0.1151, 0.2405, 0.9638) e 2 = (-0.1785, 0.1650, 0.9070) n 1 = (0.0427, -0.9682, 0.2466) n 2 = (0.4305, -0.8734, 0.2278) d 1 = (0.9925, 0.0690, 0.1013) d 2 = (0.8847, 0.4583, 0.0849) a 1 = (252.6524, 62.0000, 92.5000) a 2 = (266.5158, 65.1388, 91.5000)Table 4.1: Test data for the cone pair displayed in Figure 4.10.

Figure 4 . 10 :

 410 Figure 4.10: Different 3D views of the circular cone segments that interpolate the two contact elements given in Table 4.1.

Figure 4 .

 4 11 illustrates interpolating circular cone segments for different values of λ 1 , which were built for the same pair of contiguous contact elements. The parameter λ 1 allows for the adjustment of each pair of cone segments according to visualization needs.

Figure 4 . 11 :

 411 Figure 4.11: Pair of segments of cones ∆ 1 and ∆ 2 , which correspond to values λ 1 from left to right to 3.0, 4.5, 6.0, and 7.5, respectively.

Figure 4 . 12 :

 412 Figure 4.12: Example of a cone pair where one vertex lies inside of the region of interest.

Figure 4 . 13 :

 413 Figure 4.13: Left: example of a cone pair with sharp edges. Right: example of a s-shaped cone pair.

Figure 4 .

 4 Figure 4.12). Other times the cones of a pair lie locally on different sides of their common tangent plane, this might lead to either to sharp edges in the cone spline (Figure 4.13 Right) or a s-shaped cone pair (Figure 4.13 Left).

Figure 4 .

 4 [START_REF] Harada | Application of the Bézier curve to data interpolation[END_REF]. Note that the cones ∆ 3 and ∆ 4 present a sharp edge in the common generator. And also the vertex of the cone ∆ 7 lies inside of the region delimited by the line segments e 4 and e 5 .

Figure 4 . 14 :

 414 Figure 4.14: Circular cones interpolating the 5 contact elements given in Table 4.2.

Figure 4 .

 4 Figure 4.15 shows the visualization of one dental piece choosing the parameter λ 1 = 8 for the first pair of cones and λ 1 = 7 for the second.

Figure 4 .Figure 4 . 15 :

 4415 Figure 4.16 illustrates pairs of segments of circular cones for different values of λ 1 ,

Figure 4 . 16 :

 416 Figure 4.16: Pairs of segments of circular cones, which correspond to values λ 1 = 6, 12, 18, 24, respectively, where the blue dash line represents the common generator.

Figure 4 . 17 :

 417 Figure 4.17: 3D view of the construction of a curved slice with segments of cones for some dental pieces of the lower jaw bone.

  2 with two cone segments ∆ 1 , ∆ 2 , which possess the same tangent plane τ c along a common generator e c . The two vertices v i , i = 1, 2 of ∆ 1 , ∆ 2 are chosen on the lines e i , i = 1, 2 such that each one lies outside of the region of interest, see Figure 4.18. The common generator e c of the cone segments ∆ 1 , ∆ 2 is the line joining the points v 1 , v 2 , as it is illustrated in Figure 4.18.

Figure 4 . 6 .

 46 The directrix of each cone will be a C 1 ATPH Hermite curve r i (t) for i = 1, 2, whose end points are a i and a c , and the end tangents at these points are the unit vectors d i and d c for i = 1, 2, respectively, see Figure 4.18. The point a c and the tangent vector d c are calculated as follows. The contact point a c is selected as the midpoint of the vertices v 1 , v 2 and the tangent vector d c at this point is the mean of vectors d i , i = 1, 2, see Figure 4.18. Thus we have the C 1 Hermite data -the end points a 1 , a c , a 2 , with associated end tangents d 1 , d c , d 2 -for calculating the interpolating ATPH curves r 1 , r 2 .

Figure 4 . 18 :

 418 Figure 4.18: Pair of C 1 ATPH curves r 1 , r 2 that interpolate the end points a 1 , a c , a 2 associated to the end derivatives d 1 , d c , d 2 , respectively.

1 . 2 and d c = d 1 +d 2 2

 122 Choose v 1 and v 2 lying outside of the region of interest. 2. Let a c = v 1 +v 2

Figure

4. 18 .

 18 Following the preceding algorithm, we have chosen the vertices v 1 and v 2 , which are v 1 = (256.1051, 54.7844, 63.5861) and v 2 = (254.0189, 76.6916, 159.3996). Then we calculate the contact point a c with associated tangent d c . So, a c = (255.0620,65.7380, 111.4929)

(e 2

 2 , τ 2 ), which are interpolated with the pair of C 1 ATPH curves r 1 , r 2 calculated for the functional F1 illustrated in Figure 4.19. These same contact elements were interpolated with circular cone segments, see Figure 4.10, which were constructed following the algorithm explained in section 4.3.1.1. Note that the resulting cones constructed from a pair of C 1 ATPH curve (see Figure 4.19)

Figure 4 . 19 :

 419 Figure 4.19: Top: Cone pair defined by the vertices v 1 , v 2 and a pair of C 1 ATPH curves r 1 , r 2 that interpolate the points a 1 , a c , a 2 . Bottom: 3D views of this cone pair.

Figure 4 . 20 :

 420 Figure 4.20: Cone spline designed from the C 1 Hermite ATPH curves that interpolate the data given inTable 4.2.

Figure 4

 4 Figure 4.5.

Figure 4 .

 4 Figure 4.21 shows a cone spline constructed as explained above, which allows to visualize a dental piece of the lower maxilla. This spline is composed of four cone segments and contains each straight line that passes through the root and the "center" of the dental crown (black dashed lines in Figure 4.21). The segments of cones are connected G 1 along the common rulings (blue dashed lines in Figure 4.21).

Figure 4 .

 4 Figure 4.22 shows different 3D views of a cone spline in a dental volume that allows

Figure 4 . 21 :

 421 Figure 4.21: 3D view of the construction of a curved slice for visualizing a tooth of the lower jaw.

Figure 4 . 22 :

 422 Figure 4.22: 3D views of the construction of a curved slice with segments of cone splines for all dental pieces of the lower jaw.

Figure 4 .

 4 24 illustrates the corresponding development of the curved slice displayed on Figure 4.23.

Figure 4 . 23 :

 423 Figure 4.23: 3D views of a curved slice showing the information related to the lower jaw.

Figure 4 . 24 :

 424 Figure 4.24: View of development without distortion of the curved slice in Figure 4.23.

Figure 4 . 25 :

 425 Figure 4.25: Example of "offset" surfaces.

Figure 4 . 26 :

 426 Figure 4.26: Two "offset" spline surfaces of the surface illustrated in Figure 4.23.One "offset" in each direction, in order to provide additional information for diagnosis.

Figure 4 . 27 :

 427 Figure 4.27: Example of a tooth using the techniques of cone splines from biarcs (left) and 3D ATPH curves (right).

  Hermite de premier ordre. Les interpolants ATPH obtenues dependent de trois valeurs angulaires libres. Comme dans le cas plan, nous identifions de maniére optimale ces paramètres en minimisant des fonctionnelles de forme. Ce processus est utilisé pour calculer les courbes ATPH interpolantes qui approchent étroitement des courbes paramétriques 3D bien connues. Pour illustrer cette performance, nous présentons le processus pour une sorte de hélices. Dans le chapitre 4, nous utilisons ces splines ATPH spatiales C 1 continues pour guider des facettes développables, qui sont déployées à l'intérieur de volumes tomodensitométriques (CT) odontologiques, afin de visualiser des informations d'intérêt pour le professionnel de santé. En particulier, nous construisons des surfaces coniques par morceaux au long de courbes ATPH pour afficher des informations relatives à la structure anatomique des os de la mâchoire humaine. Basis functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 B-Basis functions for α → π, 2π . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Family of C 2 ATPH splines interpolating five points . . . . . . . . . . . . 2.4 Examples of open and closed ATPH splines . . . . . . . . . . . . . . . . . 2.5 Open ATPH spline curves with different sequences of parameters . . . . . 2.6 Closed ATPH spline curves with different sequences of parameters . . . . 2.7 The comparison of the C 2 ATPH spline, the C 2 PH quintic spline and the C 2 cubic B-spline for open curves . . . . . . . . . . . . . . . . . . . . 2.8 The comparison of the C 2 ATPH spline, the C 2 PH quintic spline and the C 2 cubic B-spline for closed curves . . . . . . . . . . . . . . . . . . . . 2.9 Open C 2 ATPH spline curves for increasing number of input data points . 2.10 Closed C 2 ATPH spline curves for increasing number of input data points 2.11 Examples of ATPH splines with different sequence of parameters α k = α . 2.12 Lipschitz constants K α . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.13 Open curves corresponding to chordal sequences of parameters . . . . . . 2.14 Open curves corresponding to arbitrary sequences of parameters . . . . . 2.15 Closed curves corresponding to chordal sequences of parameters . . . . . . 2.16 Closed curves corresponding to arbitrary sequences of parameters . . . . . 2.17 The C 2

	List of Figures
	2.1 B-

3. Cette caractérisation 3D est donnée en termes de quaternions, car cela permet d'analyser de manière appropriée leurs propriétés et de simplifier les calculs. Ensuite, nous utilisons les courbes spatiales ATPH pour résoudre le problème d'interpolation

  2.1 The absolute rotation index (R abs ) and the bending energy (E) for the family of C 2 ATPH splines in Figure 2.3. . . . . . . . . . . . . . . . . . . 2.2 Data for C 2 ATPH splines in Figure 2.4. . . . . . . . . . . . . . . . . . . . 2.3 Error comparisons for open-curve test cases illustrated in Figure 2.9 . . . 2.4 Error comparisons for closed -curve test cases illustrated in Figure 2.10 . . 2.5 Lipschitz constants K α . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Test data for open ATPH curves . . . . . . . . . . . . . . . . . . . . . . . 2.7 Test data for closed ATPH curves . . . . . . . . . . . . . . . . . . . . . . . 2.8 Test data for open ATPH curves corresponding to chordal α's . . . . . . . 2.9 Test data for open ATPH curves corresponding to arbitrary α's. . . . . . 2.10 Test data for closed ATPH curves chordal α's. . . . . . . . . . . . . . . . 2.11 Test data for closed ATPH curves corresponding to arbitrary α's. . . . . . 2.12 Test data for open ATPH curves for chordal α's . . . . . . . . . . . . . . . 2.13 Test data for open ATPH curves for arbitrary α's . . . . . . . . . . . . . . 2.14 Test data for closed ATPH curves chordal α's. . . . . . . . . . . . . . . . 2.15 Test data for closed ATPH curves arbitrary α's. . . . . . . . . . . . . . . . 2.16 Approximation errors for C 2

Table 2 .

 2

1: The absolute rotation index (R abs ) and the bending energy (E) for the family of C 2 ATPH splines in Figure 2.3.

  In Tables 2.3 and 2.4 we present comparative convergence measures (equation (2.31)) of the Newton-Raphson iteration for ATPH and PH quintic spline curves. The corresponding open and closed C 2 ATPH curves are illustrated in Figures 2.9 and 2.10, respectively.

	Shape	Number of points ε r ATPH curves ε r PH curves
	Semisine	10	ε 6 =2.3668e-16 ε 6 =2.2304e-15
	Mountain	20	ε 6 =2.6414e-16 ε 6 =2.3735e-15
	Spiral	30	ε 5 =5.0446e-15 ε 5 =5.1433e-13
	Quirky	40	ε 6 =3.2643e-16 ε 6 =2.5939e-12

It can be employed to calculate ATPH splines that interpolate large sequences of interpolation points, as is illustrated in Figures 2.9 and 2.10. Moreover, the ATPH curves offer the shape parameters α k , which allow for optimization procedures that can lead to aesthetically more pleasing curves. For example, in Section 2.6, we present the approximation of well-known trigonometric curves with ATPH spline curves, by adjusting the sequence of parameters {α k } k=1,...,N , for k = 1, . . . N .

Table 2 . 3 :

 23 Error comparisons according to(2.31) for open-curve test cases illustrated in Figure2.9.

Table 2 .

 2

4: Error comparisons according to (2.31) for closed -curve test cases illustrated in Figure 2.10.

  .8) for open curves, whereas for closed curves it is satisfied from the fourth iteration (h 4 , see Table 2.10). For arbitrary sequences of parameters the convergence is also achieved after a few iterations, see Table 2.9 for open curves and Table 2.11 for closed curves. The data presented in the aforementioned Tables are by no means "special", they are, in fact, representative of what we have observed in numerous examples. The resulting interpolating ATPH curves are displayed in Figures 2.13, 2.14 for open curves and in Figures 2.15

Table 2 . 5 :

 25 The Lipschitz constants K α corresponding to points *, * and * of Figure 2.12.

			α	K α	α	K α
			π/18	2	5π/9	549.2
			π/9	2	11π/18 1119.4
			π/6	2	2π/3	1717.7
			2π/9	2	13π/18 1965.8
			5π/18	2	7π/9	1620.0
			π/3	2	5π/6	884.2
			7π/18	9.9	8π/9	260.5
			4π/9	53.5 17π/18	21.1
			π/2	201.1	π	2
	α	K α	h 0		h 1		h 2	h 3
	5π/18	2	13.6397	0.4716	6.0255e-04 7.9437e-10
	π/3	2	0.7777	0.0197	1.3149e-05 3.6838e-12
	7π/18	9.9	7.1331	0.5428		0.0030	8.6226e-08
	4π/9	53.5	55.9560	6.9301		0.0993	1.8969e-05
	π/2	201.1	277.6468	46.7705		1.2275	7.6743e-04
	5π/9	549.2	948.3287	199.7283		8.1390	0.0120
	11π/18 1119.4 2.3628e+03	593.1812		34.1615	0.0978
	2π/3	1717.71 4.4283e+03 1.28667e+03	98.6317	0.4862

2.12: Lipschitz constants K α .

Table 2 .

 2

	α	K α	h 0	h 1	h 2	h 3	h 4	h 5
	5π/18	2	155.3650 227.3300 81.0998 4.2599 1.2083 0.3643
	π/3	2	10.0116	17.8138	5.2251 0.7660 0.2189 0.0702
	7π/18	9.9	5.7644	11.4538	3.0132 0.6796 0.2163 0.0736
	4π/9	53.5	5.7564	12.1211	3.0121 0.8746 0.3184 0.1149
	π/2	201.1	5.8799	12.77832 3.0869 1.0745 0.4455 0.1700
	5π/9	549.2	6.1268	13.6497	3.2442 1.3091 0.6142 0.2479
	11π/18 1119.4	6.5297	14.9942	3.5137 1.6184 0.8601 0.3700
	2π/3	1717.71	7.1423	17.2029	3.9430 2.0647 1.2596 0.5085

6:

Test data for open ATPH curves. The quantities listed are: α k = α for k = 1, 2, 3, Lipschitz constants K α and the convergence parameters h 0 at the starting approximation and h i , 1 ≤ i ≤ 3 at the ith iteration.

Table 2 . 7 :

 27 Test data for closed ATPH curves. The quantities listed are: α k = α for k = 1, . . . , 5, Lipschitz constants K α and the convergence parameters h 0 at the starting approximation and h i , 1 ≤ i ≤ 5 at the ith iteration.

	{α k } k=1,2,3	K	h 0	h 1	h 2

{0.3183π, 0.4138π, 0.2546π} 4.187 47.1714 2.0234 0.0139 {0.3254π, 0.4209π, 0.2582π} 5.306 49.3147 2.4345 0.0178 {0.3325π, 0.4279π, 0.2617π} 6.682 51.7001 2.8810 0.0227 {0.3395π, 0.4350π, 0.2653π} 8.364 54.3735 3.3678 0.0285 {0.3466π, 0.4421π, 0.2653π} 10.41 57.3892 3.9018 0.0355 {0.3537π, 0.4492π, 0.2723π} 12.88 60.8117 4.4918 0.0439 {0.3608π, 0.4562π, 0.2759π} 15.84 64.7163 5.1487 0.0539 {0.3678π, 0.4633π, 0.2794π} 19.39 69.1903 5.8858 0.0657 {0.3749π, 0.4704π, 0.2829π} 23.60 74.3346 6.7186 0.0798 {0.3820π, 0.4775π, 0.2865π} 28.57 80.2641 7.6659 0.0965

Table 2 . 8

 28

	{α k } k=1,2,3 K h 0 h 1 h 2 h 3 h 4 h 5	{0.0633π, 0.0926π, 0.0882π} 2 2.3213e+10 2.5265e+8 5.4780e+03 158.4100 11.7970 0.4673	{0.1607π, 0.1301π, 0.1913π} 2 1.5406e+06 1.844e+05 9.4832e+03 19.9027 7.9294e-05 1.6659e-09	{0.1496π, 0.2228π, 0.2728π} 2 3.3926e+05 4.1171e+04 1.7683e+03 2.8234 7.1794e-06 4.7123e-10	{0.3183π, 0.4138π, 0.2706π} 4.8266 25.9829 1.3490 0.0078 2.1490e-07 1.2942e-13 1.3946e-14	{0.3301π, 0.4166π, 0.2706π} 5.5734 28.8287 1.5802 0.0099 3.1082e-07 2.0321e-13 1.9635e-13	{0.3848π, 0.5286π, 0.4749π} 175.0526 155.5343 20.6231 0.3543 1.1605e-04 1.3852e-11 1.6521e-12	{0.4058π, 0.4797π, 0.5433π} 378.8906 489.8410 86.8297 2.4810 0.0018 8.915e-10 5.4512e-12	{0.5139π, 0.4606π, 0.5946π} 687.0410 1.1938e+03 250.6062 9.3172 0.0111 1.5565e-08 9.2800e-12	{0.4673π, 0.5960π, 0.6173π} 1.1631+e03 1.9173e+03 435.1409 20.7548 0.0416 1.5317e-07 2.0242e-11	{0.6667π, 0.4845π, 0.6037π} 1.3720+e03 3.6803e+03 960.3292 57.7008 0.1891 2.1050e-06 2.5398e-11

: Test data for open ATPH curves corresponding to chordal α's.

Figure 2.13: Right: open curves corresponding to chordal sequences of parameters.

Left: zoomed curves.

Table 2 . 9 :

 29 Test data for open ATPH curves corresponding to arbitrary α's.

	{α k } k=1,...,5 K h 0 h 1 h 2 h 3 h 4 h 5	{0.1910π, 0.2228π, 0.2546π, 0.2515π, 0.2578π} 0.0192 76.1675 83.7014 31.3974 5.2977 0.6760 0.1534	{0.1938π, 0.2246π, 0.2557π, 0.2521π, 0.2591π} 0.0202 66.5173 74.5346 27.4078 4.6156 0.5723 0.1309	{0.1966π, 0.2264π, 0.2568π, 0.2527π, 0.2603π} 0.0212 58.2577 66.5427 24.0002 4.0301 0.4866 0.1122	{0.1995π, 0.2281π, 0.2578π, 0.2533π, 0.2616π} 0.0222 51.1685 59.5578 21.0812 3.5263 0.4157 0.0965	{0.2023π, 0.2299π, 0.2589π, 0.2539π, 0.2628π} 0.0233 45.0670 53.4386 18.5737 3.0919 0.3569 0.0835	{0.2051π, 0.2317π, 0.2600π, 0.2545π, 0.2641π} 0.0244 39.8015 48.0653 16.4136 2.7164 0.3081 0.0726	{0.2079π, 0.2334π, 0.2610π, 0.2551π, 0.2653π} 0.0256 35.2456 43.3362 14.5479 2.3911 0.0559 0.0134	{0.2108π, 0.2352π, 0.2621π, 0.2557π, 0.2665π} 0.0268 31.2936 39.1651 12.9321 2.1088 0.2339 0.0559	{0.2136π, 0.2370π, 0.2631π, 0.2563π, 0.2678π} 0.0281 27.8567 35.4781 11.5292 1.8632 0.2058 0.0495	{0.2164π, 0.2387π, 0.2642π, 0.2569π, 0.2690π} 0.0295 24.8606 32.2123 10.1080 1.6491 0.1823 0.0441

Table 2 . 10 :

 210 Test data for closed ATPH curves chordal α's.

	{α k } k=1,...,5 K h 0 h 1 h 2 h 3 h 4 h 5 h 6	{0.1464π, 0.1528π, 0.1910π, 0.1592π, 0.1751π} 7.1410e-05 71.2695 39.8391 43.9982 7.0374 3.2907 0.9337 0.2239	{0.1440π, 0.1607π, 0.1778π, 0.1913π, 0.2069π} 2.3108e-04 148.6130 74.8750 69.8918 10.8287 4.3892 1.1029 0.2492	{0.2164π, 0.2387π, 0.2628π, 0.1614π, 0.2340π} 0.0077 129.5106 107.1060 68.5824 8.7588 2.8788 0.8349 0.2429	{0.2865π, 0.3317π, 0.2762π, 0.2134π, 0.3463π} 0.4881 117.5597 120.8166 49.1497 6.4405 1.5281 0.3982 0.1022	{0.3183π, 0.2292π, 0.4166π, 0.2706π, 0.3820π} 2.3896 176.6208 92.4435 63.9790 8.4589 2.0504 0.7064 0.2493	{0.4138π, 0.3820π, 0.5093π, 0.4775π, 0.5570π} 173.0321 36.1540 48.0362 11.6741 3.4169 1.1357 0.3744 0.1132	{0.4613π, 0.4858π, 0.4797π, 0.5433π, 0.5570π} 390.5769 44.7823 55.6177 15.3843 5.0369 1.8723 0.6788 0.2205	{0.5139π, 0.5411π, 0.4606π, 0.3836π, 0.5946π} 320.8943 27.1107 35.5599 12.5834 4.9199 2.0017 0.7118 0.2228	{0.4673π, 0.4950π, 0.5960π, 0.6173π, 0.5157π} 986.6414 23.1182 56.2846 7.8200 3.4414 1.3212 0.4027 0.1122	{0.6049π, 0.6667π, 0.4845π, 0.5585π, 0.6037π} 1.3379e+03 15.9958 26.6603 8.7193 4.2547 2.2146 0.8672 0.2853

Table 2 . 11 :

 211 Test data for closed ATPH curves corresponding to arbitrary α's.

Table 2 .

 2 12: Test data for open ATPH curves for chordal α's.

	{α k } k=1,2,3		K *	h * 0	h * 1	h * 2
	{0.3183π, 0.4138π, 0.2546π} 1.553 17.4938 0.7504 0.0052
	{0.3254π, 0.4209π, 0.2582π} 2.062 19.1672 0.9462 0.0069
	{0.3325π, 0.4279π, 0.2617π} 2.708 20.9549 1.1677 0.0092
	{0.3395π, 0.4350π, 0.2653π} 3.520 39.4275 3.7657 0.0474
	{0.3466π, 0.4421π, 0.2653π} 4.529 22.8803 1.4172 0.0120
	{0.3537π, 0.4492π, 0.2723π} 5.774 24.9727 1.6979 0.0155
	{0.3608π, 0.4562π, 0.2759π} 7.297 27.2673 2.0141 0.0197
	{0.3678π, 0.4633π, 0.2794π} 9.416 29.8064 2.3714 0.0248
	{0.3749π, 0.4704π, 0.2829π} 11.37 32.6397 2.7765 0.0310
	{0.3820π, 0.4775π, 0.2865π} 14.03 35.8247 3.2380 0.0385
	{α k } k=1,2,3	K *	h * 0	h * 1		h * 2	h * 3
	{0.0633π, 0.0926π, 0.0882π} 2.332e-10 270.6432 29.4568	6.3868	0.1847
	{0.1607π, 0.1301π, 0.1913π} 2.534e-05 19.5222	2.3373		0.1202 2.5221e-04
	{0.1496π, 0.2228π, 0.2728π}	0.0034	570.7164 69.2591	2.9747	0.0047
	{0.3183π, 0.4138π, 0.2706π}	1.9530	10.5136	0.5459		0.0032 8.6955e-08
	{0.3301π, 0.4166π, 0.2706π}	2.3620	12.2178	0.6697		0.0042 1.3173e-07
	{0.3848π, 0.5286π, 0.4749π} 120.2365 106.8301 14.1652	0.2434 7.9713e-05
	{0.4058π, 0.4797π, 0.5433π} 168.0991 217.3235 38.5230	1.1007 8.0024e-04
	{0.5139π, 0.4606π, 0.5946π} 258.6053 449.3682 94.3293	3.5070	0.0042
	{0.4673π, 0.5960π, 0.6173π} 489.6316 807.0998 183.1753 8.7369	0.0175
	{0.6667π, 0.4845π, 0.6037π} 438.3178	1175.8	306.7938 18.4335	0.0604

Table 2 .

 2 13: Test data for open ATPH curves for arbitrary α's.

	{α k } k=1,...,5 K * h * 0 h * 1 h * 2 h * 3 h * 4	{0.1910π, 0.2228π, 0.2546π, 0.2515π, 0.2578π} 0.0059 23.4974 25.8216 9.6860 1.6343 0.2086	{0.1938π, 0.2246π, 0.2557π, 0.2521π, 0.2591π} 0.0062 20.5472 23.0238 8.4663 1.4258 0.1768	{0.1966π, 0.2264π, 0.2568π, 0.2527π, 0.2603π} 0.0065 18.0180 20.5804 7.4228 1.2464 0.1505	{0.1995π, 0.2281π, 0.2578π, 0.2533π, 0.2616π} 0.0069 15.8438 18.4415 6.5276 1.0919 0.1287	{0.2023π, 0.2299π, 0.2589π, 0.2539π, 0.2628π} 0.0072 13.9698 16.5648 5.7574 0.9584 0.1106	{0.2051π, 0.2317π, 0.2600π, 0.2545π, 0.2641π} 0.0076 12.3502 14.9144 5.0931 0.8429 0.0956	{0.2079π, 0.2334π, 0.2610π, 0.2551π, 0.2653π} 0.0079 10.9470 13.4599 4.5185 0.7427 0.0831	{0.2108π, 0.2352π, 0.2621π, 0.2557π, 0.2665π} 0.0083 9.7283 12.1753 4.0202 0.6556 0.0727	{0.2136π, 0.2370π, 0.2631π, 0.2563π, 0.2678π} 0.0087 8.6671 11.0384 3.5871 0.5797 0.0640	{0.2164π, 0.2387π, 0.2642π, 0.2569π, 0.2690π} 0.0092 7.7409 10.0301 3.2096 0.5135 0.0568

Table 2 . 14 :

 214 Test data for closed ATPH curves chordal α's.

	{α k } k=1,...,5 K * h * 0 h * 1 h * 2 h * 3 h * 4 h * 5	{0.1464π, 0.1528π, 0.1910π, 0.1592π, 0.1751π} 2.1094e-05 21.0527 11.7683 12.9968 2.0788 0.9720 0.2758	{0.1440π, 0.1607π, 0.1778π, 0.1913π, 0.2069π} 7.1259e-05 45.8274 23.0890 21.5523 3.3392 1.3535 0.3401	{0.2164π, 0.2387π, 0.2628π, 0.1614π, 0.2340π} 0.0031 52.8686 43.7226 27.9966 3.5755 1.1752 0.3408	{0.2865π, 0.3317π, 0.2762π, 0.2134π, 0.3463π} 0.1972 47.5092 48.8254 19.8628 2.6028 0.6175 0.1609	{0.3183π, 0.2292π, 0.4166π, 0.2706π, 0.3820π} 0.5995 44.3138 23.1939 16.0522 2.1223 0.5144 0.1772	{0.4138π, 0.3820π, 0.5093π, 0.4775π, 0.5570π} 86.5188 18.0776 24.0189 5.8372 1.7085 0.5679 0.1872	{0.4613π, 0.4858π, 0.4797π, 0.5433π, 0.5570π} 124.6401 14.2908 17.7486 4.9094 1.6074 0.5975 0.2166	{0.5139π, 0.5411π, 0.4606π, 0.3836π, 0.5946π} 211.5806 17.8754 26.0837 8.2968 3.2439 1.3198 0.4693	{0.4673π, 0.4950π, 0.5960π, 0.6173π, 0.5157π} 343.8520 8.0569 19.6156 2.7253 1.9993 0.4605 0.1403	{0.6049π, 0.6667π, 0.4845π, 0.5585π, 0.6037π} 656.4706 7.8489 13.0819 4.2785 2.0877 1.0867 0.4255

Table 2 .

 2

15: Test data for closed ATPH curves arbitrary α's.

Table 2 .

 2 [START_REF] Sarraga | G 1 interpolation of generally unrestricted cubic Bézier curves[END_REF]: Comparison of approximation errors for C 2 ATPH spline curves and C 2 PH quintic spline curves interpolating the same sequences of points on the unit circle.

	0347e-

Number of points ε r MSE RMSE 10 ε 1 -AT P H =3.0749e-16 MSE -AT P H =1.0619e-10 RMSE -AT P H =1.0305e-5 ε 1 -P H =44.5434 ε 5 -P H =1.6278e-14 MSE -P H =1.7966e-10 RMSE -P H =1.3404e-5 15 ε 1 -AT P H =3.0696e-16 MSE -AT P H =7.0874e-11 RMSE -AT P H =8.4187e-6 ε 1 -P H =208.0952 ε 6 -P H =3.8881e-14 MSE -P H =7.9633e-11 RMSE -P H =8.9237e-6 20 ε 1 -AT P H =2.8308e-16 MSE -AT P H =5.8221e-11 RMSE -AT P H =7.6302e-6 ε 1 -P H =789.9844 ε 6 -P H =8.0038e-14 MSE -P H =7.1517e-11 RMSE -P H =8.4568e-6 curvature of the PH spline curve and the cubic B-spline are quite different from the distribution of curvature of the limaçon, where for the first and second one we have 0.8541 ≤ κ ≤ 1.7952 and 0.5087 ≤ κ ≤ 2.4232, respectively. In this example, the C 2 ATPH spline approximation errors MSE and RMSE are 1.corresponding values of the approximation errors for the ATPH spline curves and the PH quintic spline curves, both curves interpolating the same sequence of points illustrated in Figure 2.21. All test cases suggest that the ATPH spline curves can be used to reproduce quite accurately arcs of trigonometric curves. Figure 2.21: Sequences of 10, 15, respectively 20 interpolating points on the limaçon. Number of points ε r MSE RMSE 10 ε 5 -AT P H =2.0385e-16 MSE -AT P H =1.9501e-11 RMSE -AT P H =4.4160e-6 ε 5 -P H =5.5070e-4 ε 7 -P H =4.2833e-14 MSE -P H =5.3663e-8 RMSE -P H =2.3165e-4 15 ε 4 -AT P H =1.7580e-16 MSE -AT P H =5.3750e-12 RMSE -AT P H =2.3184e-6 ε 4 -P H =1.3888 ε 7 -P H =7.2197e-14 MSE -P H =2.2188e-8 RMSE -P H =1.4895e-4 20 ε 4 -AT P H =4.4777e-16 MSE -AT P H =2.2622e-13 RMSE -AT P H =4.7562e-7 ε 4 -P H =3.7687 ε 7 -P H =1.0067e-13 MSE -P H =2.6941e-9 RMSE -P H =5.1905e-5

Table 2 .

 2 17: Comparison of approximation errors for C 2 ATPH spline curves and C 2 PH quintic spline curves interpolating the same sequences of points on the limaçon.

Table 3 . 1 :

 31 Hermite data p i ,d i and p f ,d f for the examples illustrated in Figures 3.1 and 3.4.

Table 3 . 2 :

 32 Various choices of the angular parameters θ, β for the ATPH interpolants illustrated in Figures 3.1 and 3.4.

Table 4 .

 4

2: Test data for the cones spline displayed in Figure 4.14.

Table 4 .

 4

[START_REF] Farin | Handbook of CAGD[END_REF]

. The values listed are: the orthonormal base e i , d i , n i for the planes τ i i = 1, . . . , 5 and the end points a i which lie on the lines e i for i = 1, . . . , 5, respectively. This data is interpolated with 8 circular cone segments displayed in

Table 4 . 3

 43

  1 , τ 1 ),

	F1	F2	F3
	r 1 150.9804 4.8725e+06 8.1061e+06
	r 2 150.0417 8.5217e+06 1.5179e+07

Table 4 . 4 :

 44 Values of the functionals F1 , F2 , and F3 , given by the minimization process.

Assuming that b (0) is calculated from an interpolating cubic B-spline without loops (in practice given by the conditions (2.22)). The good solution b is an iteration of b (0) such that its difference with the previous iteration is within a certain tolerance, namely T ol = 1.0 -20 .

Figure 3.5: The spatial ATPH curves obtained by minimizing F1 . The corresponding PH quintic interpolant is drawn in red.

Figure 3.8: The spatial ATPH curves obtained by minimizing F2 . The corresponding PH quintic interpolant is drawn in red.

Figure 3.11: The spatial ATPH curves obtained by minimizing F3 . The corresponding PH quintic interpolant is drawn in red.

If the tooth has more than one root, the straight line must pass through the midpoint of these and the "center" of the dental crown.

where φ 0 , φ 2 are angular parameters. It follows from the interpolation conditions r(0) = p i and r(α) = p f the following equation

then we have that

From the previous equation we would like to find the quaternion coefficient A 1 . After some simplifications, completing squares and making several re-arrangements, we obtain

Substituting (3.16) in the above equation and setting ρ := n 0 -6n 2 and ω := 4n 2 (1+c 2 ), we have respectively. From these Tables, we observe that the minimization method applied to the aforementioned fairness functionals, select "good" spatial ATPH curves among the three-parameter family of solutions.

(1.0,0.0,1.0) (0.0,1.0,1.0) case 2 (-0.8,0.3,1.2) (0.0,1.0,1.0) case 3 (0.4,-1.5,-1.2) (-1.2,-0.6,-1.2) case 4

(-0.8,0.3,1.2) (0.5,-1.3,-1.0) case 5 (10.0,0.0,10.0) (0.0,1.0,1.0) Apart from further exploring the application potential of the presented methods in the medical field, it will be an interesting topic of future research to construct spatial C 2 ATPH splines that interpolate a sequence of points with prescribed end conditions. Note that this leads to a system of quadratic equations in quaternion unknowns, with three degrees of freedom in each subinterval. Suitable strategies for fixing these free parameters based on optimizing shape measures must be addressed.