Convergence Analysis on a Second Order Algorithm for Orthogonal Projection onto Curves

Regarding the point projection and inversion problem, a classical algorithm for orthogonal projection onto curves and surfaces has been presented by Hu and Wallner (2005). The objective of this paper is to give a convergence analysis of the projection algorithm. On the point projection problem, we give a formal proof that it is second order convergent and independent of the initial value to project a point onto a planar parameter curve. Meantime, for the point inversion problem, we then give a formal proof that it is third order convergent and independent of the initial value.

To find the orthogonal projection of a point onto surfaces or parametric curves, Limaien and Trochu [1] constructed an auxiliary function and obtain all its zeros.Based on the geometric method [4], Li et al. [2] use the geometric iterative method with second order approximation properties.The common features of [2] and [4] are that they both use the curvature information, and they are independent of the initial value, respectively.If the algorithm in Hu et al. [4] is applied in the spatial parametric curve, the order of convergence will be one rather than two.In addition, the line segment between the center of circle, test point and osculating circle will not intersect in some cases.
Pottmann et al. [3] adopt the idea of the ICP algorithm and provide a new method to geometrically align a point cloud to a surface and to the related registration problems.An alternative notion has been presented in [3], which depends on both instantaneous kinematics and the geometry of the squared distance function corresponding to a surface.Evidenced by local convergence analysis, as well as by experiment results, the newly proposed ICP method is faster to converge than the classic one.Piegl et al. [5] present a method to project a point onto NURBS surfaces with three stpes: decompose a NURBS surface into quadrilaterals, project the test point onto the nearest quadrilateral, and finally get the parameter from the nearest quadrilateral.Mortenson [21] converted the projection problem into the problem of obtaining a polynomial's root, and then adopted the Newton-Raphson approach to get the root.
Chen et al. [18] used clipping circle technology and proposed an approach to find the minimum distance between a point and a NURBS curve.In the same vein, Chen et al. [19] used clipping sphere technology and provided an approach to find the minimum distance between a clamped B-Spline surface and a point.Similar to [18,19], Young-Taek Oh et al. [22] used an efficient culling technique to reduce redundant curves and surfaces and proposed an efficient method to project a point to its nearest point on a set of freeform surfaces and curves.Chen et al. [20] proposed a new quadratic clipping approach to calculate a root of a polynomial f (t) of degree n within an interval.Distinct from the classic method in R 1 space, it approximates (t, f (t)) and gets three quadratic curves in R 2 space, which results in a higher order of approximation.Among various methods to find the roots to a polynomial equation in computer graphics and computer aided geometric design, clipping methods with the Bernstein-Bézier form demonstrate great numerical stability.To search the roots of a polynomial f (t) in an interval, Chen et al. [11] proposed a rational cubic clipping method, of which the approximation order and the corresponding convergence rate to search a single root are both seven.Based on their previous method in [11], Chen et al. [12] revise a rational cubic clipping approach to get two bounding cubic within O(n) time, which could get a faster convergence rate of five instead of four in a previous attempt.For Bézier curves, Chen et al. [23] improved the algebraic pruning method to prevent invalid roots after conversion of the point projection problem into a root finding problem, which is applicable for NURBS curves because, in NURBS curves, is not hard to be changed into Bézier form.
Song et al. [24] presented an algorithm to solve the problem of orthogonal projection of parametric curves onto B-spline surfaces.They use a tracing method and then create a polyline approximating the pre-image curve of the projection curve.Bharath Ram Sundar et al. [25] have constructed a measure of distance computation between surfaces and curves.To find point inversion for a parametric surface, Wang et al. [26] proposed a new method to give a simple procedure to get the parameters determined by a given point on a surface.They essentially formulate ordinary differential equation systems, which are determined by the intersection curve segment or by the projection curve segment.Li et al. [27] combined a first order algorithm with Newton's second order algorithm to propose the hybrid second order approach, which goes faster than the current methods and does not depend on the initial value.The adoption of ADMM in [28] contributes to refine the estimate in each iteration because it can incorporate information about the direction of estimates gotten in previous steps.In sum, those methods use lots of techniques, for instance, converting the projection problem into a root searching problem for a nonlinear equation system, subdivision approaches, and geometric approaches as well as circular clipping approaches.
One widely used, fast and robust method is introduced by the paper [4] and consists of two steps: calculate parameter values by projecting points onto curvature circles in the first step and then compute parameter increments using Taylor's expansion of the surface or curve in the second step.The paper [4] claims to be quadratic in its convergence rate but does not prove that.In this paper, regarding the point projection problem, we formally prove that projecting a point onto a planar parameter curve is of second order convergence and independent of the initial value.On the point inversion problem, we formally prove that it is of third order convergence and independent of the initial value.

Convergence Analysis
In this part, we conduct a convergence analysis of iteration, presented in formula Equation (5) in paper [4].We parameterize the curvature circle c such that it shares the same Taylor's polynomial with the curve c.Then this equality holds: In R 2 , we deal with the Equation (1) using the method in the paper [4] and get det(q On the point projection problem, to verify that the algorithm defined by Equation (2) (hereafter, the iterative method (2)) is quadratically convergent, we firstly derive the expression of the footpoint q.The corresponding parameter is denoted as α when we orthogonally project a test point p onto a parameter curve c(t), where p = (p 1 , p 2 ) and c(t) = ( f 1 (t), f 2 (t)).By the definition of the minimum distance between a curve and a point, we will get the relation as follows where ).The unit normal vector, the curvature, the relative curvature, and the radius of curvature circle of curve c(t) can be expressed as follows, respectively, In this paper, we prove the case where Additionally, the proof method also holds for the case only with an opposite sign to the previous case.So the formula Equation ( 7) can be written as The center of curvature circle of the planar curve c(t) is Now, inserting Equations ( 4) and ( 6) into Equation ( 9), we get and The line connecting the test point p and the center of curvature circle m can be parameterized as where u is a parameter.Also, the equation of the curvature circle is As the footpoint q = (q 1 , q 2 ) is the intersection of the line Equation ( 12) and the curvature circle Equation ( 13), substituting Equation ( 12) into Equation ( 13), we get Solving Equation ( 14), we have Since the footpoint q lies between the center of curvature circle m and the test point p, the parameter u must belong to (0, 1).By Equation ( 15), the parameter u will be From Equations ( 3), ( 10), ( 12) and ( 16), we get horizontal coordinate q 1 of the intersection q From Equations ( 11), ( 12) and ( 16), we obtain vertical coordinate q 2 of the intersection q From Equation ( 17), we have From Equation (18), we get Till now, we provide the detailed illustration of the iterative method (2).Next, we will use a method similar to those in the literature [29][30][31] to prove the theorem.The standard convergence analysis method in numerical analysis in [29][30][31] are used, namely, e n+1 = C 0 e k n + o(e k+1 n ), where e n+1 = t n+1 − α, e n = t n − α, k is a positive integer, C 0 is constant.
) be a real parametric curve function in two dimensions, suppose q(t) − c(t) has first and second derivatives within an interval I.When q(t) − c(t) has only a single root α ∈ I, and t n stays close sufficiently to α, then the iterative method (2) will be of second order convergence for the point projection problem.
Furthermore, we have and For the sake of simplicity, we use the Maple 18 package to calculate the Taylor's expansion in the following way.Combine the derivation of Equations ( 4)- (20) with the orthogonal projection condition for Equation (3) and Equation ( 21)- (26), and then simplify these equations, we acquire parameterized increment ∆t as follows: where Substituting Equations ( 28) into (27), we can get the following error equation similar to that in [32], This implies that the algorithm defined by the iterative method ( 2) is quadratically convergent for the point projection problem.The proof is completed.
Then, we will illustrate that the algorithm expressed by the iterative method (2) is third order convergent for the point inversion problem.Theorem 2. Let q(t) − c(t) = (q 1 (t) − f 1 (t), q 2 (t) − f 2 (t)) be a real parametric curve function in two dimensions, suppose q(t) − c(t) has first and second derivatives within an interval.When q(t) − c(t) has only one single root α ∈ I, and t n stays close sufficiently to α, so the iterative method (2) will be of cubic convergence for point the inversion problem.
Proof.Given the point inversion problem is a special case of the point projection problem, the proof of this theorem will use the results in the proof of Theorem 1.For the same part, we will not repeat.Since the test point p is on the parametric curve c(t), from Equation (3), we can know that (p 1 , p 2 ) = (a 0 , b 0 ).From Equation ( 29) and (p 1 , p 2 ) = (a 0 , b 0 ), we acquire the following error equation in [32], This implies that the algorithm expressed by the iterative method ( 2) is third order convergent for the point inversion problem.The proof is completed.Theorem 3. The iterative method (2) is independent of the initial value for point projection and inversion problem.
Proof.Firstly, we presents the interpretation for Figure 1.Assume there exists a horizontal axis t and planar parametric curve c(t).There are two points on the planar parametric curve c(t): the first point c(t n ) with t n on the horizontal axis, and the second point c(α) is the corresponding orthogonal projection point of the test point p onto the planar parametric curve c(t), where α is the corresponding parameter value on the horizontal axis.While, according to the iterative method (2), the line segment determined by point p and point c(α) and tangent line of planar parametric curve c(t) at t = α are mutually perpendicular.The tangent line of the planar parametric curve c(t) through point c(t n ) will determine a footpoint q, which is obtained by orthogonally projecting test point p onto the tangent line of the curve c(t) at t = t n .The corresponding parametric value t n+1 of footpoint q on the horizontal axis is obviously the next iterative value used by the iterative method (2).The corresponding parametric value of the middle point of the line segment determined by the point c(t n ) and the footpoint q is M.
Secondly, we present the proof of Theorem 3. When the iterative method (2) only has one solution, it is independent of the initial value for point projection and inversion problem.The proof method is similar to the ones in references [33][34][35].It is not difficult to show that the corresponding parameter of the first dimensional for the planar parametric curve on the two-dimensional horizontal plane is t.When the iterative method (2) begins to iterate, as the iterative parameter value satisfies the inequality relationship t n < α, by the graphical demonstration, the parameter corresponding to the footpoint q will be t n+1 .The middle point of point (t n+1 , 0) and point (t n , 0) is (M, 0), which implies the relationship M = t n + t n+1 2 .Additionally, because of 0 < ∆t = t n+1 − t n , there is an inequality relationship t n < M < α.This result sufficiently yields the inequality  2) is independent of the initial iterative parametric value.Alternatively, we will present a descriptive proof that the iterative method ( 2) is independent of the initial iterative parametric value.This proof method is also similar to the one in the literature [33][34][35].Suppose that the starting point lies on the left of α, then footpoint q will lie on the right to the starting point, and hence ∆t will be positive.By the iterative method (2), the iterative sequence is expressed as follows, where If t n < α, the sequence t n is strictly monotonically increasing.If t n > α, by at most three iterations, the sequence t n converge to α.Notice that this iterative sequence performs like an attenuated pendulum and converges to α.Furthermore, when the starting point lies to the right of α, the convergence holds in the same way.To summarize, the footpoint q depends on the test point p all the time, rather than on the starting iteration point.So the iterative method (2) is independent of the initial value.The proof is completed.
Remark 1.The iterative method (2) is independent of the initial value, but it only requires that c(t) is C 2 , while the methods in the literature [33][34][35] need to satisfy the sufficient conditions for global convergence.In addition, the iterative method (2) fits in with robustness and stabilization in Professor Les A. Piegl's view [36].
Remark 2. If k(t) = 0 of the Equation ( 5), then t = t 0 .Hence, we calculate the distance between the point c(t 0 ) on the curve c(t) and the test point p, i.e., d 1 = p − c(t 0 ) .At the same time, in order to deal with the possible failure of the projection or inversion, we use the parameter perturbation method and increment the parameter t 0 by a small positive constant ε, t 0 = t 0 + ε, so that the iteration can continue.During the iteration process, when the derivative of the curve c(t) at t = t 0 becomes zero, the same perturbation approach is applied.Eventually, d =min{d 1 , d 2 , d 3 , ...} becomes the minimum distance between the curve c(t) and the test point p.

Numerical Experiments
This section provides numerical evidence to compare the behavior of the Newton iterative algorithm and the iterative method (2) discussed above.
Example 1.We consider the curve c(t) = (t, sin(t)) with test point p = (2, 2) and the corresponding parameter of projection point α = 1.7838126561069.Table 1 shows the results of Newton iteration and the iterative method (2).The results verify that the iterative method (2) has better convergence and robustness, while the Newton algorithm is dependent on the initial value but is unstable (See Figure 2).  1 to represent that the Newton method does not converge, and the same notation applies in Table 2 to Table 6 .Example 2. We use the curve c(t) = (t, cos(t)) with test point p = (2, 5) and the corresponding parameter of projection point α = 0.402360707683495.Table 2 shows the results of Newton iteration and the iterative method (2).The results also verify that the iterative method (2) has better convergence and robustness, while Newton algorithm is dependent on the initial value but is unstable (See Figure 3).Example 3. We consider the curve c(t) = (t 2 + sin(t), sin(sin(t)) + cos(t)) with test point is p = (−1, 1) and the corresponding parameter of projection point α = −0.26523161027243385.Table 3 shows the number of iterations for the Newton iteration and the iterative method (2).The results show the iterative method (2) has better convergence and robustness, at the same time, the Newton's method sometimes converges slowly and depends on the initial value but is unstable (See Figure 4).Example 4. We use the curve c(t) = (t, sin(t) + cos(t)), t ∈ [3,6] .Test point p = (−2, −6) and the corresponding parameter of projection point is 3.1213051310788399.Table 4 shows the number of iterations for Newton iteration the iterative method (2).Results show the iterative method (2) has better convergence and robustness, at the same time, Newton's method is sometimes dependent on the initial value but is unstable (See Figure 5).5 shows number of iterations for Newton iteration and the iterative method (2).Results show the the iterative method (2) has better convergence and robustness, at the same time, Newton's method is sometimes dependent on the initial value but is unstable (See Figure 6).6 shows the number of iterations for Newton iteration and the iterative method (2).Results show the iterative method (2) has better convergence and robustness, at the same time, Newton's method is sometimes dependent on the initial value but is unstable (See Figure 7).Remark 3. The iterative method (2) orthogonally projects a test point onto a parametric curve c(t).For the case with multiple orthogonal points, our approach is changed in the following way: (1) Divide a parameter interval [a, b] parametric curve c(t) into M subintervals with equal length.
(2) Randomly select an initial iterative parametric value in each subinterval.
(3) Using the iterative method (2) and using each initial iterative parametric value, iterate, respectively.Suppose that the iterative parametric values are α 1 ,α 2 ,...,α M , respectively.If we attempt to search all the solutions as fast as possible, divide a parameter interval [a, b] into M subintervals with equal length such that M is sufficiently large.
According to Figures 5 and 7, either blue point is the only orthogonal projection point of the test point, but not the point with minimum distance.For t ∈ [−6, 6] in Figure 5, four parameter values are −5.8827783434357889,−2.3086073340017088, 1.0869042777921642, 3.1213051310788399, respectively.According to Remark 3, it is easy to find that the orthogonal projection point determined by the parameter value of −2.3086073340017088 is the one with minimum distance, but the orthogonal projection points determined by the other parameter values are not.For t ∈ [−10, 3] in Figure 7, eight parameter values are −8.1893433622682510,−5.9663756035363462, −4.8785206470437198, −3.0152895587637027, −1.5956295339830861, −5.8855753296468201 × 10 −2 , 1.6840779491863052, 2.8953236261013390, respectively.In the same way, only the orthogonal projection point determined by the parameter value of −5.8855753296468201 × 10 −2 is the one with minimum distance.

Conclusions
In this paper, regarding the point projection problem, we formally prove that projection of a point onto a planar parameter curve is quadratically convergent and independent of the initial value.On the point inversion problem, we formally prove that the inversion problem is cubically convergent and independent of the initial value.In the future, we would like to present a formal proof of convergence on projecting a point onto a spatial parametric curve and surface in literature [4].
namely, there is an iterative error expression |e n+1 | < |e n |, where e n = t n − α.If t n > α, then the proving method is similar to the case with t n < α.Thus, for the planar parametric curve c(t), there is an iterative error expression |e n+1 | < |e n | in the two-dimensional horizontal plane, namely, the iterative method (

Table 3 .
Comparison of the robustness and effectiveness for the Newton's method and the iterative method (2).

Table 5 .
Comparison of the robustness and effectiveness for the Newton's method and the iterative method (2).

Table 6 .
Comparison of the robustness and effectiveness for the Newton's method and the iterative method (2).