Spline Curve Distance Minimization: Robust and performant computation of distance queries for spline curves on B-spline form

Minimization

Robust and performant computation of distance queries for spline curves on B-spline form

Torgeir Baraldsnes

Master’s Thesis Autumn 2016

Torgeir Baraldsnes 1st August 2016

This thesis investigates the problem of efficiently determining the shortest distance from a point to a spline curve, and the location of the nearest point on the spline. Splines are versatile and popular approximation functions used in a plethora of applications, including geometric modeling software, physical simulations and differential equation solvers, computer graphics, image analysis, statistical tools, and more. Distance queries are an important part of the basic operations on which many other operations are dependent. It is therefore essential that these operations are robust and highly performant.

A robust and efficient algorithm for solving this problem is presented, which combines established methods with recent techniques developed for spline curves and functions. The algorithm is conceptually simple with clearly separated constituent parts. It can therefore be easily adapted to more specialized applications. The approach can also be extended to related, more general, problems.

I would like to thank my supervisors Martin Reimers and Fritz Albregtsen for their guidance, help and enthusiasm. I would also like to thank my family and friends for their invaluable support during all the work leading up to this thesis.

1 Introduction 1

1.1 Problem statement . . . 2

1.2 Notation and definitions . . . 5

1.2.1 Optimization . . . 5

1.2.2 Splines and B-splines . . . 6

1.2.3 Operations on splines . . . 10

1.3 Formal description of the problem . . . 13

2 Preliminary Theory 15 2.1 Tools and techniques . . . 15

2.1.1 Newton’s method . . . 15

2.1.2 Preprocessing . . . 16

2.1.3 Culling and heuristic methods. . . 17

2.2 Existing approaches . . . 17

2.2.1 Tangent cone methods . . . 18

2.2.2 Curvature based second order iteration . . . 18

2.2.3 Hyperplane clipping. . . 18

2.2.4 Circle clipping. . . 19

2.2.5 Separating axes and k-DOPs. . . 19

2.3 Products of splines . . . 20

2.3.1 Relevance to distance computation with splines . . . 20

2.3.2 Construction of the direct B-spline product . . . 20

2.3.3 Disadvantages of direct products . . . 21

2.3.4 Solution: Products of Bernstein polynomials . . . 21

2.4 Converting to Bernstein Bézier form . . . 24

2.4.1 Advantages of the BB-form . . . 24

2.4.2 Various methods of conversion . . . 24

2.5 Root finding for splines and polynomials . . . 26

2.5.1 Application to spline optimization . . . 26

2.5.2 Algorithm outline . . . 26

2.5.4 Curve generalization. . . 28

2.5.5 Challenges with storage and data structures . . . 28

3 Algorithms 31 3.1 High level description of the algorithms . . . 31

3.1.1 Basic algorithm . . . 31

3.1.2 Adaptive algorithm . . . 32

3.2 A more detailed look . . . 32

3.2.1 The main algorithms . . . 32

3.2.2 Sub-routines. . . 36

3.2.3 Adaptive algorithm example . . . 40

4 Numerical tests and analysis 45 4.1 About the tests . . . 45

4.2 Results and plots . . . 45

4.3 Analysis . . . 48

4.4 Profiling . . . 49

4.5 Analysis summary . . . 50

5 Conclusion 53 5.1 Summary . . . 53

5.2 Further Improvements . . . 55

5.3 Related problems . . . 56

A Source Code 61

1.1 Problem illustration . . . 1

1.2 Application Example . . . 3

1.3 B-spline functions plot . . . 7

1.4 The control polygon of a spline curve . . . 9

1.5 Knot insertion . . . 12

2.1 Bézier product summation . . . 23

2.2 Finding the roots of a spline function . . . 26

2.3 The convex hull of a spline curve. . . 29

3.1 Algorithm outline . . . 35

3.2 Adaptive algorithm example: The spline curve . . . 41

3.3 Adaptive algorithm example: The BB-segments . . . 42

3.4 Adaptive algorithm example: The squared distance function . . . 43

4.1 Timings of algorithms by degrees . . . 46

4.2 Timings of algorithms by number of segments . . . 47

4.3 Profiling results . . . 51

A.1 De Casteljau’s algorithm . . . 62

A.2 De Boor’s Algorithm . . . 62

A.3 Böhm’s algorithm . . . 63

A.4 Auxiliary code . . . 64

A.5 Bézier Product algorithms . . . 65

INTRODUCTION

This first chapter will present the prerequisite mathematical framework and problem statement, setting the stage for the rest of the thesis and serving as reference for later chapters. The first section provides some background and intuition for the problem in an informal fashion, while the formal treatment is deferred until after the necessary notation and definitions have been introduced.

P f

Figure 1.1: The smallest distance from a point to a curve

1.1 Problem statement

The goal of this thesis is to construct an algorithm for the efficient computation of the smallest distance between a point and a spline curve. This problem is a practical and useful one in its own right, as well as a stepping stone towards more general optimization problems for splines. The investigation into this problem begins with an introduction of the central objects and terminology.

Curves, points and distances In order to speak meaningfully about distance it is necessary to consider a pair of objects between which the distance is measured. In this case, the objects are a curve along with a point residing in the same space. Furthermore, it will be assumed that the space under consideration is Euclidean and endowed with a Cartesian coordinate system. These assumptions admit the use of the familiar notions of Euclidean distance and norm∥ · ∥as induced by the standard inner product.

Let f be a curve andPbe an arbitrary point. Then, the problem is to compute the value of distance(f,P) =inf

x ∥^f(_x)−^P∥ ^(1.1)

The pointPwill be referred to as thequery point. The nature of these objects will be more properly defined insection 1.2.

Curves can be defined in various ways depending on the application and what degree of generality is preferred. The type of curves considered here will be analogous to the physical trajectory of a particle: An arbitrary path through space such that for any given point in time, x, there is a single point in space, f(_x), associated with it. Such a curve is called a parametric curve. For a parametric curve, computing the smallest distance to a point is closely related to the problems of finding the point, or points, on the curve where the smallest distance is attained and the associated argument values. Because distance is relatively easily recovered from the argument of the minimum, by evaluation, the emphasis will be on finding the argument and essentially getting the distance estimate as a byproduct. Determining the argument of points on the curve is commonly referred to aspoint inversion. The argument where the smallest distance is attained may not be unique. The test point may be equidistant from separate parts of the curve, or even lie on a self intersection.

Any argument value corresponding to aglobal minimumwill be considered a valid solution for the purposes of this discussion.

When designing algorithms for practical applications, it is usually beneficial to restrict attention to a class of widely used curves. One such class is that ofspline curves, particularly inB-splineform.

Splines, what are they? Splines are mathematical functions used for modeling curves and surfaces, as well as other types of approximation and interpolation problems where control over smoothness and local changes is required. They have their name from the tool used by drafters in the shipbuilding-, aeroplane- and automotive-industry as a flexible ruler to draw smooth, streamlined curves through a set of points. The mathematical counterparts have somewhat analogous properties which has found them a wide variety of uses across several disciplines.

Why splines? Splines are ubiquitous in geometric modeling software tools and graphical design applications. They can be found in music creation software, text rendering, computer graphics

applications, and all manner of scientific and engineering tools. Some of the reason for their widespread use can be attributed to the fact that they are general and versatile function approximators which lend themselves nicely to both intuitive editing by hand and automatic fitting to data.

Motivation With splines being so pervasive in so many applications follows the need for high quality algorithms for transforming and querying them. Locating the closest point on a spline is an example of such a frequently occurring operation. It has natural connections to other fundamental geometrical operations such as intersection and separation queries.

Figure 1.2: Pixel center distance to a shape defined by a curved outline. Note that while the figure shows the underlying shape ideal, the overlapping pixels can only take on a single intensity value, which may lead to aliasing in simple sampling schemes. Distance queries can be used to estimate coverage, leading to a more accurate rendition of the shape.

A concrete example Distance queries for splines can be used in vector graphics applications, where each shape is represented by an outline of spline curves. To be viewed on a screen, the shapes need to be sampled and mapped to pixels on the screen with color and intensity values determined in part by how much the shape covers the pixel. If the pixel naively sampled in the center of each pixel and the output is given by binary value based on whether the sample is inside the the shape or not, the output will suffer from jagged edges caused byaliasing. A sampling scheme which mitigates this problem can be devised by estimating pixel coverage in terms the distance of the pixel center to the shape. This information is then used to give pixels which are close to the boundary of the shape

intermediate values based on the estimated coverage. The distance query requires a method to find the closest point on a planar spline curve to the point representing the pixel center. The scheme is not perfect, however, as it will overestimate the coverage thin features that are small compared to the pixel yet still stretch close to the sample point. Still, it should serve as a good enough estimate, or serve as a starting point for more elaborate methods. Alternatively, a high resolution distance map can be precomputed and downsampled. An example of such an approach can be found in [12].

Anticipated challenges The principal challenge regarding minimization of spline curves is to deal with the non-convexity of the curve in an efficient manner. The most common methods for root- finding used in such problems require an initial estimate of the minimum to satisfy certain conditions so that convergence can be ensured. Finding such an estimate and ensuring that the conditions are met is not altogether straight forward. There have been developed methods that work by finding estimates and checking for sufficient conditions to begin iteration. Unfortunately, they are often plagued by the need for laborious subdivision and preprocessing steps. In fact, much of the progress made on the problem has been in the form of techniques which reduce the number of subdivision steps required. This is discussed further inchapter 2.

Before formulating and investigating the formal problem statement it is necessary to review some basic definitions and terminology from mathematical optimization and spline theory, and to introduce the notation employed throughout the rest of the text.

1.2 Notation and definitions

This section covers the basic theory of splines and B-splines and establishes the prerequisite definitions and notation.

1.2.1 Optimization

The problem of finding the smallest distance is an instance of anoptimization problem. It is therefore useful to briefly review some optimization terminology. The notation used mostly follows the treatments in [3] and [9].

Definition 1.2.1(Optimization terminology).

Let f be a function with domainXand rangeY⊆R. A pointxis said to be Astationary pointif f is differentiable atxand f^′(_x) =0

Acritical pointifxis either a stationary point or f is not differentiable atx Aglobal minimum pointif f(_x)≤ f(_y)for ally∈ X

Aglobal maximum pointif f(x)≥ f(y)for ally ∈X

Alocal minimum pointif there exists anϵ>0such that f(_x)≤ f(_y)for allywithin a distanceϵof x, that is|y−x|<_ϵ.

Alocal maximum pointif there exists anϵ> 0such that f(_x)≥ f(_y)for allywithin a distanceϵof x.

Ifxis either a minimum point or maximum point of f then the function value, then f(x) is called either aminimum value or maximum value, or simply a minimum or maximum, respectively. Anextremumis a function value which is either a minimum or maximum of the function.

Convexity

The notion of convexity is central in the field of mathematical optimization and is the primary focus of the important sub-field calledconvex optimization. It is used to guarantee global optimality and as a condition for efficiency for many important methods. Informally, convexity refers to the property of being without ’holes’ or ’dents’ such that a convex set is one where any pair of its points can be connected by a straight path entirely contained in the set. This intuition is formalized in the following definition:

Definition 1.2.2(Convexity).

A setCof elements from a vector spaceVis said to beconvexif and only if for each pair of elementsu,v∈C, each pointw(_t) = (1−t)u+_tvis also inCfor anyt∈[0, 1]⊂R

A vector w(_t)inDefinition 1.2.2is called aconvex combination of the vectorsuandv. This, and related terms are summarized here.

Definition 1.2.3(Convexity cont.).

Aconvex combinationof vectors{^vi}ⁿ_i=1is a linear combination

∑

λiv_i

whereλi ∈ Ralso satisfies∑ⁿi=1λi =1and allλi ≥^0.

Theconvex hullof a set is the set of all possible convex combinations of its elements.

Aconvex functionis a function whose epigraph is a convex set, or equivalently a function f : X → R such f((1−t)_x1+_tx2) ≤ (1−t)_f(_x0) +_{t f}(_x1) for every pair x0,x1 ∈X.

1.2.2 Splines and B-splines

Knots and degree The parametric continuity (number of continuous derivatives)C^kof the curves will be dictated by two parameters which are used in the definition of B-splines.

Definition 1.2.4(Knot vectors and polynomial degree).

• The polynomial degree d ∈ N is the highest exponent occurring in a univariate polynomial. This notion is extended to splines to denote the degree of all basis functions used to represent the spline.

• Aknot vector tis an non-decreasing sequence of real numbers(_ti)ⁿ_i=⁺₁^d+1calledknots.

The knot vector partitions an interval of the real line into contiguous, half-open, sub- intervals[_ti,ti+1). The intervals are said to beemptyif and only ifti =_ti+1for some interval. A knottiinthasmultiplicitymifti occursmtimes int. This text assumes knot multiplicitydat the endpoints. This has the effect of clamping the spline in each end. This is assumption is convenient and simplifies the discussion without causing loss of generality.

Basis functions The algorithms presented in this thesis will all deal with splines represented in B- spline form. B-splines are a family of functions used as a basis for representing splines for which the control points have particularly nice properties. Splines are typically viewed as linear combinations of B-splines, which are piecewise polynomials that form a basis for a spline space with a given domain partition and degree. B-splines are usually defined by a recurrence relation, as follows:

Definition 1.2.5(B-Spline).

A B-Spline is a function B_i,d,t : R → R, with index i ∈ {1, . . . ,n} where n ∈ N, polynomialdegree d ∈ Nand a knot vector t = (_ti)ⁿ_i=⁺₁^d+1. The B-splines satisfies the recurrence relation

Bi,d,t(_x) = ^x−ti

t_i+d−tiBi,d−1,t(_x) + ^ti+1−x

t_i+d+1−ti+1Bi+1,d−1,t(_x) with base case

B_i,0,t(x) =

{1, ifti ≤x <_ti+1

0, otherwise .

Figure 1.3: B-splines of degree 1, 2 and 3 from top to bottom. Each bump (B-spline) is a basis function, and together they span a spline space of a given degree and knot vector. The figure shows the effect of clamping the splines at the interval endpoints

Curves The space of polynomial spline curves can now be defined in terms of B-splines. All spline curves f considered here will be elements of this space. The following is adapted from Definition 2.13 and 2.14 in [17]

Definition 1.2.6(Spline curves).

Lett =⁽_tj^)n+d+1

j=1 be a nondecreasing sequence of real numbers, and letq≥²be an integer.

The space of all spline curves inR^qof degreedand with knotstis defined as S^q_d,t=

{ n j

∑

c_jBj,d,t

c_j ∈R^q for 1≤ j<_n }

Thencontrol points,(c_i)ⁿ_i=1, of a spline are the coefficients in the linear combination of basis functions used to represent the spline. Scalar control points may alternatively be referred to as control coefficientsor simply coefficients to distinguish them from higher-dimensional points.

The control points of a given spline depend upon which basis is chosen, and will therefore be prefixed with the name of the basis wherever necessary to avoid confusion.

An element ofS^q_d,t is called a spline vector function or aparametric spline curveof degreed with knot vectortand control points(

c_j)_n

j=1, or justspline curvefor short.

Using the control points to approximate the curve by a piecewise affine (degree one) spline is a central idea in spline theory. It essentially represents an un-smoothed view of the curve that interpolates the control points. Each B-splineBi,d,t in a spline has a corresponding coefficient c_i and set ofinterior knots ti+1, . . . ,t_i+dassociated with it. The interior knots are the only knots on which the B-spline is non-zero. The average of the interior knots are thei’th knot average oft. This is information is used to parametrize the approximate spline. The construction is formalized in the following definition adapted from those given in [22] and [24].

Definition 1.2.7(Control polygon of spline curves).

Lettanddbe a knot vector and the degree of a spline curve f inS^q_d,tas inDefinition 1.2.6.

Thecontrol polygonof f is the piecewise affine curve obtained by connecting neighbouring control points (c_j)ⁿ_j=1by straight line segments. Each control point has an associatedknot average(_t^¯_i)ⁿ_i=1, wheret¯i = ^ti+1+···_d^+ti+d. The control polygon of f ontwill be denoted

Γf,t(_x) =_Γj,f,t(_x), for tj ≤x <_tj+d (1.2) where

Γj,f,t=

(¯tj+1−x t¯j+1−t^¯j

c_j+ ^x−t^¯j

t¯j+1−t^¯j

c_j+1

)

= (1−λ)c_j+_λcj+1 andλ= _t ^x−t¯j

j+d+1−t_j+1. For simplicity, the notationΓandΓj, will be used where the parameters are clear from the context.

Figure 1.4: A spline curve and its control polygon consisting of connected line segments Bézier curves and the Bernstein polynomial basis Bézier curves are special cases of spline-curves which consist of only one polynomial segment. This means that they are just regular polynomial curves, but on a special basis, the Bernstein polynomial basis, which possesses many properties equal or similar to those of B-splines. Many of these properties hold only over a an associated interval, analogous to how B-spline representations are associated to their knot vector. For instance, the non-negativity and convex hull properties only hold over the specified interval. The interval can be arbitrary, but any such curve can be obtained by translation of a curve over the domain[0, 1].

Definition 1.2.8(Bernstein Bézier form of a polynomial).

Letdbe a non-negative integer and let(c_i)^d_i=0bed+1coefficients inR^m And let thei’th Bernstein basis polynomialb_i,d:R→Rof degreedbe defined by

b_i,d(x) = (d

i )

xⁱ(1−x)^d−i. Then a Bézier curve f is the the linear combination

f =

∑

c_ibi,d

restricted to the interval[0, 1].

Bézier curves can be subdivided and evaluated using thede Casteljau algorithm(seeListing A.1).

Because they are a special case, Bézier curves can be viewed as a spline with the a knot vector of [a^d+1,b^d+1]where the superscript denotes multiplicity. The canonical knot vector isa = 0and b=1. Any spline algorithms for B-splines therefore work on Bézier curves. Because the degree can be inferred from the number of control points, it is possible to drop the knot vector and only store the interval boundsaandb. The knot vector can then be retrieved from this data. The coefficients will be said to be onBernstein-Bézier-form, or BB-form for short. For the knot vector[_a^d+1,b^d+1], they are equal to the B-spline coefficients of the curve.

1.2.3 Operations on splines

Evaluation While splines can be evaluated by computing the value of the B-splines through the recurrence relation inDefinition 1.2.5, it is more common to employ a method known asde Boor’s algorithm[10]. The algorithm works by performing repeated convex combinations of the control points. The coefficients of the convex combinations are determined by the function argument and the knots. De Boor’s algorithm is a generalization of de Casteljau’s algorithm which is used to evaluate Bézier curves in the same fashion.

Algorithm 1.2.1(de Boor’s Algorithm).

Letx be a point in the domain of a spline f such that x ∈ [_tµ,t_µ+1)for some µ, where d ≤ µ ≤ n. Then the only B-splines on the knot vector which are non-zero at x are (B_i)^µ_i=µ−d. By lettingc⁰_i = c_i, the value f(x) = c^d_µcan be computed from the recurrence relation

c^k_i = (1−λi,k)c_i^k₋⁻₁¹+_λi,kc^k_i⁻¹ forµ−d+_k≤ i≤µ, and whereλi,k = _t ^x−ti

i+d−k+1−ti.

Knot Insertion Another important operation is the process of computing the coefficients of a given spline on a refined knot vector, known asknot insertion. The new coefficients that emerge

when inserting new knots will cause the control polygon to more closely approximate the spline near the parameter value of the new knots, seeFigure 1.5. This property can be used as an alternate way of plotting the spline by inserting enough knots such that the control polygon can be plotted directly.

The refinement property of knot insertion can be leveraged to perform iterative spline queries, such as the root-finding method developed in [22], which many queries such as intersection and optimisation can be transformed into, this will be discussed further inchapter 2.

There are two main methods for inserting knots: Boehm’s method [1] and the Oslo Algorithm [8]. Boehm’s method is relatively simple and computes the result of inserting a single knot into the knot vector. The Oslo Algorithm is somewhat more involved, but can be used to insert an arbitrary number of knots.

Algorithm 1.2.2(Böhm’s Method).

Letτ = (_τ)ⁿ_j=⁺₁^d+1 be a given knot vector and lett = (_t)ⁿ_j=⁺₁^d+2be a knot vector obtained by inserting a knotzinτin the interval[_τµ,τ_µ+1). If

f =

∑

c_jBj,d,τ =

n+1 i

∑

a_iBi,d,t (1.3)

then(a_i)ⁿ_i=⁺₁¹can be expressed in terms of(c_j)ⁿ_j=1through the formulas

a_i =







c_i, if1≤i≤µ−d,

z−τi

τi+d−τic_i+_τ^τi+d−z

i+d−τic_i−1, ifµ−d+1≤i≤µ,

c_i−1, ifµ+1≤i≤n+1.

(1.4)

Derivatives of splines Differentiation is another central operation in optimization, so it is useful to review some basic results for computing derivatives of splines. As with regular polynomials the derivatives of splines can be written as splines of lower degree. For non-differentiable points it is natural to consider one-sided derivatives, subderivatives, or jump discontinuities. These only occur at the breakpoints of the spline. It is possible for a curve to have cusp singularities, but these correspond to zeroes of the derivative with sign changes in the vector components of the curve, and not necessarily discontinuities in the derivative. From Theorem 3.16 in [17]

Figure 1.5: The effect of inserting a knot on a cubic spline. The spline is represented on a refined knot vector, in this case a single additional knot, which causes the control polygon to contain an additional point and to be brought closer to the spline. The dotted lines represents the new edges of the control polygon along which the corners of the old control polygon is clipped.

Theorem 1.2.1(Derivative of B-splines).

The derivative of thej’th B-spline of degreedontis given by DBj,d,t(_x) =_d

(Bj,d−1,t(_x)

tj+d−tj − ^Bj+1,d−1,t(_x) tj+d+1−tj+1

)

ford≥¹and for any real numberx. The derivative ofBj,d,tcan also be expressed as

DBj,d,t(_x) = ^d d−¹

( x−tj

t_j+d−tjDBj,d−1,t(_x)− ^tj+d+1−x t_j+d+1−tj+1

DBj+1,d−1,t(_x) )

ford≥2

Similarly, using the notation from expression (2) in [22] the derivative of the spline f ∈S^q_d,tcan be written written

f^′ =

∑

∆cjBj,d−1,t

where

∆c_j = ^cj−c_j−1

¯tj−t^¯j−1

= ^cj−c_j−1 t_j+d−tj

andc_j = 0for anyj<0,j≥ n, and∆cj =0if the interval[_tj,tj+d)is empty, and consequently the corresponding B-spline is identically zero.

1.3 Formal description of the problem

Now that the necessary basic definitions have been established, it is time to recast the problem in a more formal presentation.

Revisiting the problem Consider a spline curve f : R →R^mof degreedwhose representation on B-spline form has n control points (c_i)ⁿ_i=1, each in R^m, on the d+1-extended knot vector t =⁽_tj^)n+d+1

j=1 with eachtj ∈R. Then the curve can be written on the form f =

∑

c_iBi,d,t (1.5)

Then, for a query pointP∈R^mthe distance minimization problem minimize

x∈[t1,t_n+d+1) ∥f(x)−P∥ (1.6) is equivalent to solving the squared distance minimization problem

minimize

x∈[t1,t_n+d+1) ∥^f(x)∥^{2 (1.7)} for a translated spline

f =

∑

(c_i−^P)_Bi,d,t

A solution to this problem is on the form (_x^∗,∥^f(_x^∗)∥), for a parameter value x^∗ for which

∥^f(_x^∗)∥ ≤ ∥^f(_x)∥for all pointsxin the domain of f. Thenx^∗ is aglobal minimum. It is also worth repeating that solutions are not necessarily unique in general. In such cases the convention adopted here will be to return a single solution pair, although one could equally well choose to return all such solutions or specify solution intervals where applicable.

A common approach to solving the problem is by searching stationary points, which is anyx satisfying∥^f(_x)∥^′ = 0and separately checking the points in the domain where the spline is non- differentiable, which together constitute the critical points.

By expanding the equation for the stationary points of the squared distance function, we get

Dx∥f(x)∥²=0 (1.8)

=⇒ Dx(f(_x)· ^f(_x)) =_Dx ( n

i

∑

c_iBi,d,t(_x)·

∑

i=j

c_jBj,d,t(_x) )

(1.9)

= _Dx ( n

∑ ∑

)

(1.10)

In going from (1.6) toEquation 1.8, some information is lost resulting in solutions admitting local and global maxima as well as local minima. This can be remedied by requiring the second derivative to be non-negative. Unfortunately, this limits the applicability of the procedure to convex partitions of the objective function, which requires a method to determine such partitions. It would be more ideal to employ a method that does not waste time by computing redundant stationary points.

The product of B-splines, as appearing in (1.10) can be represented by higher degree B-splines on a suitable knot vector. The resulting spline is one-dimensional due to taking the scalar product of the coefficients. It can then be optimized just like a regular one-dimensional spline function.

Using Leibniz rule from Equation 1.8 gives the following simple form f(_x)·Dxf(_x) = 0.

The interpretation of this is that the stationary points are those for which the position and tangent vectors are orthogonal. A procedure implementation based on this form might seem appealing:

Computing the derivative before the product for the purpose of reducing the degree of one of the factors. However, this somewhat breaks the symmetry of the square product, taking away a possible source of efficiency. Thesavings per segment, in the number of computed product terms, in the symmetrical case is ^(d+₂^1)d versus justdfor the early differentiation. The early differentiation also forces the use of root finding instead of a direct minimization procedure. This would work, but requires additional detection and processing of the non-differentiable point of the spline as well as distinguishing between maxima and minima by a second derivative test.

It is important to remember that it is the spline curve f which is given in this problem, represented by a tuple((c)ⁿ_i=1,d,(_tj)ⁿ_j=⁺₁^d+1), andnot g(_x) =∥^f(_x)∥which is the actual objective function. If the functiongwas given on spline form, for instance, the problem could be solved by applying a suitable existing optimization method for one-dimensional spline functions. The current challenge lies in the combined effort of obtaining selected parts of, and minimizing,gefficiently by somehow finding a better approach than computingg(_x)via evaluation of f(_x).

PRELIMINARY THEORY

This chapter contains a review of earlier literature connected to spline curve optimization problems and as well as a discussion on each of the central building-blocks used in both existing algorithms and the algorithms of this thesis.

2.1 Tools and techniques

There are a number of frequently occurring established techniques used in problems of this type which can largely be divided into two main categories: Preprocessing and optimization. The optimization methods come in a variety of forms but most of the relevant ones are Newton’s method and other methods similar to it. These methods, which are based on fixed point iteration, usually share a great deal of strengths and weaknesses.

2.1.1 Newton’s method

In its standard formulation, Newton’s method is a simple and efficient method for finding roots of functions. It is simple, at least conceptually, to adapt it to locate stationary points of functions, for the purpose of optimization, by replacing the objective function with its derivative.

The standard Newton update step of an estimatex_k isx_k+1 = _xk− _f^f′⁽(^xx⁾). To locate a root of an objective function on the form(

∥f(_x)∥^2)′the update step specializes to x_k+1= _xk− ^f′(_x)· ^f(_x)

f^′′(x)·f(x) + f^′(x)²

which requires three spline evaluations each step. Alternatively, the product g = _f · f can be computed, and differentiated so that the update term becomes−_g^g′⁽(^xx⁾). This form requires only two evaluations, as is usual, but each evaluation might be more computationally expensive, and obtaining the product itself may incur a significant cost. This will be discussed further insection 2.3.

Convergence concerns Newton’s method is a numerical method designed to iteratively produce better approximations to the position of a root of the function from an initial estimate. As such, questions with regards to stability and convergence characteristics are primary concerns. In short

satisfied, there exists a neighborhood around the root for which all points converges to the root.

This is what is often meant by requiring the initial estimate to be ’close’ to the root. In many cases though it is more useful to find a sufficient condition for the class of functions in question than to check all the necessary ones.

Whenever Newton’s method converges, the convergence rate is quadratic except when the derivative of the function is zero in the neighbourhood of the root. This means, roughly speaking, that the error is squared in each step, producing twice as many correct significant digits in each step as it converges towards the root. Consequently, the method typically only needs a very low number of iterations to reach full precision for both 32-bit and 64-bit floating point numbers. Because the evaluations required for Newtons method are fairly inexpensive, it is likely more worthwhile to focus efforts on the methods used to obtain a starting-point, which is considered next.

2.1.2 Preprocessing

In order to ensure convergence of Newton’s and related methods, it is often necessary to perform some manner of preprocessing to locate a good initial estimate. For spline functions, this usually involves some local or global refinement procedure and a sufficient condition on the control points such that global optimality is ensured. This usually depends properties such as convexity or monotonicity in some form.

Sampling One of the most simple preprocessing methods is to sample the objective function densely to obtain approximate local information. To minimize the function, one may simply choose to start Newtons method at the closest sample point. This is not guaranteed to return the global minimum, however, and is not even guaranteed to converge at all. For spline curves, it should be possible to find a sufficient sampling density to get correct convergence in most cases. However, the number of samples required to get decent results in such a scheme is usually large enough to render it inefficient compared to more sophisticated methods. One advantageous aspect to sampling beyond simplicity is the fact that it requires very little in the way of memory and consequently memory allocation. This in contrast to the related process of subdivision, which uses essentially the same algorithm, but gives better local information by storing the refined control points near the sample point.

Subdivision A very common form of preprocessing for spline curves, especially Bézier-splines, is subdivision. When a curve is subdivided, the previous curve-representation is transformed into an equivalent form where the curve has been split into a pair of independent pieces (sub-curves). For splines on B-spline form, this is usually achieved using the de Boor algorithm. For splines on BB- form, the de Casteljau algorithm is used. Performing a subdivision step is a divide and conquer approach to locating a starting point for Newtons method. The process can be recursively applied to each of the sub-curves until one or more fulfils a sufficient condition to start Newtons method.

Sub-curves that do not contain a minimum can be removed from further consideration.

2.1.3 Culling and heuristic methods

The process of removing portions of the search space which need not be considered further goes by several names in different disciplines, but the idea is the same: Reducing the search space can cut down on superfluous work down the line. Alternatively, one can see the reduction of the search space as a complementary way of increasing the precision of the solution: By removing extraneous parts the remaining space shrinks, causing the bounds on the solution becomesharper.

Criteria When designing a culling method it is important to choose the culling criteria carefully, not only to guarantee correctness, but to ensure that testing the criteria is significantly less computationally expensive than just processing the discarded elements.

Culling strategies are usually geared towards average case complexity or ’typical use-case’

performance. They rely on the fact that worst case scenarios, where the method may perform worse with culling than without, are rare in practice. Rather than performing a rigorous analysis on the effect of the culling criteria in the average case, they are usually selected based on observed usage patterns and application specific knowledge. For spline curves, it is natural to wish to capitalize on the convex hull property to bound curve segments which may then be tested for culling against existing bounds.

Heuristics Whenever new bounds are acquired, and there are multiple candidates which might contain the solution, a decision must be made regarding which candidates are considered in the next step. The choice is made on the basis of aheuristic function, which ranks the available candidates. The choice will often have a significant impact on the performance of the algorithm on different types of problem instances. More elaborate heuristics specializes the algorithm towards certain applications.

There are a wide range of useful heuristics to consider for the spline curve minimization. For Computer Aided Geometric Design and Manufacturing the curves are mostly of low polynomial degree and dimension and often exhibit local structure which may be factored in when choosing heuristics. In data-science, statistics and image analysis applications, on the other hand, generally deals with very high dimensional data which may render techniques designed for CAD/CAM applications less efficient. Likewise one can imagine cases where the splines are of much higher polynomial degree than usually found in CAD/CAM applications, for instance where degree- elevating operations like products, and especially function composition, are frequently used [11].

2.2 Existing approaches

Numerous papers have been published on this problem for different classes of curves, including the different types of curves in the spline family, for example NURBS [6], Bezier curves [7], and special cases like planar, fixed degree (usually cubic) [30] or cardinal spline curves. Ko and Sakkali provide a fairly recent and concise survey of developments on orthogonal projection onto curves and surfaces in CAD [15]. Papers on minimization of arbitrary degree non-uniform B-spline curves such as those considered here are not as common, but they are a special case of NURBS curves they may therefore benefit from any method used to optimize those.

A number of advances have been made over the last decade, some of which are discussed next.

2.2.1 Tangent cone methods

In the paper ’Distance Extrema for Spline Models Using Tangent Cones’ [14], Johnson and Cohen present a method for locating extrema from a point to a spline curve or surface. Crucially, the method is designed to locateall distance extrema of a spline model. The method is therefore centered around first obtaining all solutions toEquation 1.8. To that end, tangent cone hierarchies are employed to cull away the parts of the spline which can not contain an orthogonal tangent. The idea is to bound the tangents of an interval of the spline by a cone and test against this cone rather than the set of tangents on the interval. The interval is in turn is bounded by a sphere around the local convex hull.

The parts of the spline which satisfy the orthogonality condition are kept and further subdivided, while the rest are culled. Finally, when the angular precision of the tangent cone method reaches a user specified threshold, Newton’s method can be applied to rapidly increase the precision. The method was also extended to dynamic queries.

The method turned out to be robust and especially performant for dynamic queries. Some potential drawbacks seems to be the need for recursive subdivision of Bézier patches, which is expensive. The method, as described, also does not distinguish between the global minimum and other extrema though the modification is likely to be straightforward.

2.2.2 Curvature based second order iteration

Rather than directly treating distance minimization, the paper ’A second order algorithm for orthogonal projection onto curves and surfaces’ [13] by Hu and Wallner tackles the problem of computing orthogonal projections by the means of an iterative scheme based on using curvature of the osculating circle to approximate the curve in each step. The circle is assumed to be parametrized such that the second order Taylor expansion coincides with that of the curve. By sampling the curve and its first two derivatives at the current parameter estimate, the curvature of the circle is computed, via the area spanned by the first and second derivative vectors, and used to obtain a new parameter estimate. A similar approach is used for the surface case, where the curvature circle is aligned with the vector from the estimate point to the projection of the query point onto the tangent space of the surface.

Like Newton’s method, which it resembles, the curvature based approach requires an initial estimate to begin iteration. The algorithm is reported to be fairly insensitive to the choice of estimate.

However, the authors do not mention any requirements to guarantee convergence, though it is reasonable to expect conditions similar to those of Newton’s method. By the nature of the method, it does not by itself find more than one of the possible projections onto the curve, and which of the projections is obtained depends on the initial parameter estimate.

Advantages of the approach include second order convergence, insensitivity to the initial estimate, and parametrization independent iteration. Although it does not directly compute global minima, it can be used to speed up the final iteration of subdivision based methods after a sufficient condition for convergence to a global minima is met.

2.2.3 Hyperplane clipping

Selimovic details a hyperplane clipping method in the paper ’Improved algorithms for the projection of points on NURBS curves and surfaces’ [27]. The hyperplane elimination criterion works by

determining whether all control points of a NURBS curve lies of the far side of a hyperplane separating the segment and query point. Using the vector from the query point to the closest endpoint, each control point of the segment needs to be bounded away from the query point by a hyperplane orthogonal to the vector and containing the endpoint. If the condition is met, the sub- curve can be safely culled. As noted in [23], the condition is equivalent to checking whether the query point lies in the exterior Voronoi cell of the convex hull mapping to the endpoint, without explicitly computing the convex hull as in [18]. Because the method relies on the endpoint interpolation property the NURBS curve must be fully subdivided at each step, after which the elimination criterion is then applied. The paper also makes use of tangent cones in a somewhat similar fashion to [14]

2.2.4 Circle clipping

Chen et al describe an approach based on the squared distance function in the paper ’Computing the minimum distance between a point and a NURBS curve’[6]. The resulting technique is referred to as a circular clipping method. The curve is converted to Bernstein Bézier-form and the squared distance function is computed. An upper bound on the distance to the curve is maintained by the algorithm, and is conceptualized as the squared radius of a clipping circle. Each segment of the squared distance function with all its control points further than the clipping radius can be safely eliminated. Then, the curve is recursively subdivided and clipped until Newton’s method can be applied to each remaining segment. The method demonstrated superior performance compared to earlier methods in the provided test cases. However, only fairly short Bézier curves were tested, so it is not immediately clear how well the method scales in practice.

2.2.5 Separating axes and k-DOPs

The paper ’Efficient point-projection to freeform curves and surfaces’[23], by Oh et al, combines the approaches of [6] and [27]. It is noted that computing the squared distance function can be expensive, which motivates the use of a preceding clipping procedure. The idea is to construct what is essentially a discrete approximation of a clipping sphere, before the product is computed, to eliminate segments. A predetermined set ofkdirections are chosen for the purpose of performing a separating axis test. The directions can be regarded as the normal vectors of hyperplanes enclosing a convex polytope. This construction is frequently used in computer graphics, and often called a k-Discrete Oriented Polytope. In this case, the k-DOP approximates a clipping sphere of radius equal to the smallest available upper bound on the distance, which can be obtained from interpolating endpoints or taking a small number of samples. Remaining sub-curves are then clipped approximately against the hyperplanes. Further segments are eliminated using the criterion of [27], and the radius can be updated. The squared distance function is computed for the remaining parts of the curve, which are subsequently subdivided until a Newton-type method can be applied.

A common trait in most of the recent approaches have been to leverage usefulness of the squared distance function in clipping and minimization by constructing it explicitly on spline form. The next section looks into the construction of this function for spline curves.

2.3 Products of splines

2.3.1 Relevance to distance computation with splines

In order to work directly with the squared distance function, instead of via the spline curve, the symbolic product needs to have been computed beforehand. This can be advantageous as the control points of the product spline more directly conveys information about the squared distance function.

The idea is to leverage the additional information granted by the product spline control polygon to capitalize on specialized optimization methods for one-dimensional spline functions. Unfortunately, the direct spline product is known to be rather tricky to implement efficiently. A brief investigation into the theory of spline products is followed by a possible remedy to the efficiency problem.

2.3.2 Construction of the direct B-spline product

The theory for direct multiplication of B-splines was first developed by Mørken in [21]. The results were initially presented via the theory of dual functionals for B-splines, using Marsden’s Identity [20] and discrete B-splines [19]. The results were later formulated using blossoming by Ueda [29].

A more recent paper by Chen et al [5] gives a more efficient algorithm based on the blossoming approach.

The degree Each polynomial segment of the product spline is the product of segments of the factor splines which overlap in the domain, meaning that the process can be thought of as regular polynomial multiplication. The difficulty lies in efficiently obtaining the B-spline representation of the resulting function. The result of the multiplication is a segment with increased degree. Therefore the degree of the B-splines used to represent the product must also be elevated. Regardless of representation, the product of two factor polynomials of degreed1andd2isd1+_d2, so the product B-splines must be of this degree as well.

Forming a new knot vector The order of continuity of a spline at a knot is equal to the degree of the spline minus the multiplicity of the knot. It can be shown, by using Leibniz’ formula, that the order of continuity of a product of two functions f₁ ∈ C^a and f₂ ∈ C^b is min(_a,b). Using this the multiplicity of the knot in the product can be determined. By considering a product of two splines f = f₁f2, where f, f₁and f2are of degreeD,d1andd2respectively, and a knot occurs with multiplicitym1inf1andm2inf₂. Then, the unknown multiplicitymof that knot infcan be solved for by requiring thatD−m=min(_d1−m1,d2−m2), givingm=_D−^min(_d1−m1,d2−m2). This also holds for knots with multiplicity zero by generalizing toc= min(_c1,c2), wherec,c1,c2

is the order of continuity for the splines in a given point. For squaring the B-spline the multiplicity of the knots are simply increased by the degree. Such symmetries are a potential source of increased efficiency for this special case, as the number of knots added is the same everywhere.