Total Least Squares Fitting Of Bezier And B-Spline Curves To . - CORE

View metadata, citation and similar papers at core.ac.ukbrought to you byCOREprovided by Calhoun, Institutional Archive of the Naval Postgraduate SchoolCalhoun: The NPS Institutional ArchiveFaculty and Researcher PublicationsFaculty and Researcher Publications2001Total Least Squares Fitting of Bezierand B-Spline Curves to Ordered Data.Computer Aided Geometric DesignBorges, C.F.C.F. Borges and T.A. Pastva, Total Least Squares Fitting of Bezier and B-Spline Curves toOrdered Data. Computer Aided Geometric Designhttp://hdl.handle.net/10945/38066

Computer Aided Geometric Design 19 (2002) 275–289www.elsevier.com/locate/comaidTotal least squares fitting of Bézier andB-spline curves to ordered dataCarlos F. Borges , Tim PastvaNaval Postgraduate School, Code NA/BC, 93943 Monterey, CA, USAReceived 1 October 2000; received in revised form 1 November 2001AbstractWe begin by considering the problem of fitting a single Bézier curve segment to a set of ordered data so thatthe error is minimized in the total least squares sense. We develop an algorithm for applying the Gauss–Newtonmethod to this problem with a direct method for evaluating the Jacobian based on implicitly differentiating apseudo-inverse. We then demonstrate the simple extension of this algorithm to B-spline curves. We present someexperimental results for both cases. 2002 Elsevier Science B.V. All rights reserved.Keywords: Data fitting; Total least-squares; B-splines1. IntroductionA degree n Bézier curve segment is a parametric polynomial curve whose shape is determined by anordered set of control points (xj , yj ) for j 0, 1, 2, . . . , n. In particular, for 0 t 1, the componentsof the curve are given byx(t) n Bjn (t)xjj 0andy(t) n Bjn (t)yj ,j 0where Bjn (t) is the j th Bernstein polynomial of degree n, that is n jnt (1 t)n j .Bj (t) j* Corresponding author.E-mail address: borges@nps.navy.mil (C.F. Borges).0167-8396/02/ – see front matter 2002 Elsevier Science B.V. All rights reserved.PII: S 0 1 6 7 - 8 3 9 6 ( 0 2 ) 0 0 0 8 8 - 2

276C.F. Borges, T. Pastva / Computer Aided Geometric Design 19 (2002) 275–289The problem we are interested in solving can be stated as follows. Given an ordered set of data pointsmn{(ui , vi )}mi 1 and a non-negative integer n m, find nodes {ti }i 1 and control points {(xj , yj )}j 1 thatminimize 2 2 mnn nnui .(1)Bj (ti )xj vi Bj (ti )yjj 1i 1j 1In short, we are looking for the best total least-squares fit of a Bézier curve segment to a set of ordereddata points in the plane.We will find it much easier to proceed if we recast the problem in a more convenient linear algebrasetting. To this end we introduce the following definition:Definition 1. Given a vector t Rm and a non-negative integer n we define the Bernstein matrix of degreen, which we denote by Bn (t), to be the real m (n 1) matrix whose i, j element is Bjn (ti ). In particular nB0 (t1 ) B1n (t1 ) . . . Bnn (t1 ) B0n (t2 ) B1n (t2 ) . . . Bnn (t2 ) .Bn (t) . . B0n (tm )B1n (tm ). . . Bnn (tm )Be warned that for ease of notation we will generally omit the arguments of matrix functions when noconfusion is possible.This definition allows us to rewrite the objective function in expression (1) as follows 2 ux Bn , (2) v Bn y 2where x, y, u, and v are vectors whose elements are the xi , yi , ui , and vi respectively.It is of particular interest that the variables in this problem are of two distinct types—those that appearlinearly (the control points) and those that appear non-linearly (the nodes). Note that for a fixed t theoptimal associated control points are found by solving a linear least squares problem. Formally we maywrite uxBn,(3) Bn vywhere Bn is a generalized inverse (see (Rao and Mitra, 1971)) of Bn . Specifically, Bn is any (n 1) mmatrix function satisfying the following equationsBn Bn Bn Bn ,(4)(Bn Bn )T(5) Bn Bn .We note that the Moore–Penrose generalized inverse, Bn , is one such matrix although there can beothers.11 This is a formal derivation and it is important to note that one should not form generalized inverses as a method forcomputing solutions to linear least-squares problems.

C.F. Borges, T. Pastva / Computer Aided Geometric Design 19 (2002) 275–289277This structure leads to a useful approach to solving the problem. In particular, one can substitute (3)into the objective function (2) and formally uncouple the variables. It is then possible to solve for theoptimal t independently of the control points by minimizing the variable projection functional 2 uu Bn Bn . (6)r(t) Bn Bnv 2vIt is important to note that r(t) rT r, where the residual vector, r, is also a function of t and is givenby r(t) PB n uPB n v ,(7)where PBn Bn Bn is the orthogonal projector onto the range of B and PB n I PBn is the orthogonalprojector onto the null-space of BnT . Note that both of these are also matrix functions although we havesuppressed their arguments.1.1. Numerical evaluation of r(t)From a computational standpoint, we may reliably evaluate the variable projection functional byexploiting the QR factorization (see (Golub and Van Loan, 1996)) of the Bernstein matrix. In particular,there exists a permutation matrix Π such that R1,1 R1,2,Bn Π [Q1 Q2 ]00where [Q1 Q2 ] is an m m orthogonal matrix and R1,1 is r r, upper-triangular, and invertible wherer rank(Bn ). Note that R1,2 vanishes whenever Bn is full rank. It is easily verified that 1 T R1,1 Q1 (8)Bn Π0satisfies Eqs. (4) and (5), from which it follows directly thatPBn Q1 QT1 ,PB n Q2 QT2 ,whence 2 2 r(t) QT2 u 2 QT2 v 2 .2. Minimizing the variable projection functionalIn this section we will concern ourselves with minimizing the variable projection functional whichamounts to solving a non-linear least-squares problem. Although a variety of methods exist for suchproblems we will explore the Gauss–Newton approach.The Gauss–Newton approach involves using an affine model at the current point tk as followsmk (t) r(tk ) J (tk )(t tk ),

278C.F. Borges, T. Pastva / Computer Aided Geometric Design 19 (2002) 275–289where the i, j element of the Jacobian matrix isJ (t)i,j ri (t). tjMinimizing mk (t) yieldstk 1 tk J (tk )r(tk ).For a variety of reasons it is prudent to use some form of stepsize control which leads to the followingalgorithmtk 1 tk αk J (tk )r(tk ),where J (tk )r(tk ) gives a descent direction and αk is usually chosen to guarantee that the step does infact decrease the residual.The critical step for us is computing the Jacobian since this involves differentiating an expressioninvolving a generalized inverse of the Bernstein matrix. In particular, the j th column of the Jacobianmatrix is given by P Bn r(t) tj u . PB n tjv tjNotice that the Jacobian has an upper and a lower half, and furthermore, that to evaluate them we needto be able to evaluate PB nx, tjwhere x can be either u or v. It is essential therefore that we are able to differentiate a generalized inverseand we will find the following theorem will quite useful in this regard (see also (Hanson and Lawson,1969) or (Golub and Pereyra, 1973)).Theorem 1. Let A be an m n matrix function and let A be an n m matrix function such thatAA A A and (AA )T AA . Then T PA A A PAA PAA.(9) tk tk tkProof. First we note that by the product rule A PA A PA A PA tk tk tkand since PA A A we have A A PA A PA, tk tk tkwhence A A A PAA PA PA . tk tk tk tk

C.F. Borges, T. Pastva / Computer Aided Geometric Design 19 (2002) 275–289279Multiplying both sides on the right by A and using the fact that PA AA yields A PAPA PA A . tk tkTransposing both sides and invoking the symmetry of PA we also have T PA A PAA.PA tk tkNow, since PA is idempotent PA2 PA PA PA PA PA tk tk tk tk A A T PA A PA A. tk tkFinally, we note that PA I PA and hence PA PA tk tkfrom which the result follows. Applying this to our case yields 1 T 1 T T PB nR1,1 Q1R1,1 Q1 Bn B T Q2 QT2ΠΠ T n Q2 QT2 . tj tj tj00It is clear that Bn tjis all zeros with the exception of its j th row which is nB0 (tj ) B1n (tj ) . . . Bnn (tj ) tjand this structure may be exploited to reduce the total number of operations required to evaluate theJacobian. In particular, we can simultaneously compute all of the columns in the following way. LetBn (t) be the derivative of the Bernstein matrix which can be easily computed with Bn (t) n 0 Bn 1 (t) Bn 1 (t) 0 , where 0 is a single column of zeros. Then the Jacobian is given by Q2 QT2 diag(P u) P T diag(Q2 QT2 u),(10)Q2 QT2 diag(P v) P T diag(Q2 QT2 v)where 1 T R1,1Q1.P Bn Π0 Careful implementation of this formula is important to maximize economies.

280C.F. Borges, T. Pastva / Computer Aided Geometric Design 19 (2002) 275–2893. Implementation detailsAs we are fitting a single Bézier segment the problem may be simplified, without loss of generality, byfixing the nodes associated with the first and last data points to occur at the ends of the curve segment, inparticular we set t1 0 and tm 1 respectively. This has the effect of preventing the curve segment fromgrowing excessively outside the region of the data set. Although it does not force the endpoints of thecurve segment to coincide with the first and last points of the data set, it does restrict them to be tangentsto circles centered at these points.Indeed, experience shows that failure to enforce this restriction generally leads to serious problemsas the curve segment tends to expand without bound. With this restriction it is useful to reduce thedimension of the search vector to m 2 by deleting the first and last elements (t1 and tm ). The Jacobianmatrix becomes (2m) (m 2) and may be computed by replacing the first and last rows of Bn withzeros, then evaluating formula (10), and finally deleting the first and last columns from the resultingmatrix.In our implementation we have used the stepsize control parameter αk to enforce the condition thatthe nodal values satisfy 0 ti 1 for all i. If taking a full step would violate this constraint, then we setαk to a value that will preserve it. Having done this, we then use a primitive line search before taking thestep. In particular, we check to see if the current step will in fact decrease the residual. If not, then werepeatedly halve αk until it does.As with any iterative method, a good initial guess can be most valuable. The simplest is the chordlength parameterization. In particular, we let di [ui vi ]T and make our initial guess for the nodal values ij 1 di di 1 2ti mj 1 di di 1 2for i 1, 2, . . . , m with t0 0.This works reasonably well but we have noted better results using an affine invariant norm (see(Nielson, 1987)) in place of the 2-norm in the formula given above. In particular, we use the normdefined byx2V xT V x,where V is the inverse of the data covariance matrix, which is given by 1 m1 Tdi di.V m i 1Another excellent method is the affine invariant angle metric which is developed in the very elegantpaper of Foley and Nielson (1989). It is somewhat more involved in its computation but worked verywell in our experiments.Finally, we use a traditional test for convergence and stop whenrcurrent rlast tol,rlastwhere rcurrent is the squared residual at the current step, rlast is the squared residual at the last step, and tolis an arbitrary tolerance. In order to push the algorithm to the extreme we set tol 10 12 for the examplespresented in the next section.