Hybrid Second-Order Iterative Algorithm for Orthogonal Projection onto a Parametric Surface

To compute the minimum distance between a point and a parametric surface, three well-known first-order algorithms have been proposed by Hartmann (1999), Hoschek, et al. (1993) and Hu, et al. (2000) (hereafter, the First-Order method). In this paper, we prove the method’s first-order convergence and its independence of the initial value. We also give some numerical examples to illustrate its faster convergence than the existing methods. For some special cases where the First-Order method does not converge, we combine it with Newton’s second-order iterative method to present the hybrid second-order algorithm. Our method essentially exploits hybrid iteration, thus it performs very well with a second-order convergence, it is faster than the existing methods and it is independent of the initial value. Some numerical examples confirm our conclusion.


Introduction
In this paper, we discuss how to compute the minimum distance between a point and a parametric surface, and to return the nearest point (footpoint) on the surface as well as its corresponding parameter, which is also called the point projection problem (or the point inversion problem) of a parametric surface.It is a very interesting problem in geometric modeling, computer graphics and computer vision [1].Both projection and inversion are essential for interactively selecting curves and surfaces [1,2], for the curve fitting problem [1,2], for reconstructing surfaces [3][4][5] and for projecting a space curve onto a surface for surface curve design [6,7].It is also a key issue in the ICP (iterative closest point) algorithm for shape registration [8].Mortenson (1985) [9] turns the projection problem into finding the root of a polynomial, then finds the root by using the Newton-Raphson method.Zhou et al. (1993) [10] present an algorithm for computation of the stationary points of the squared distance functions between two point sets.The problem is reformulated in terms of solution of n polynomial equations with n variables expressed in the tensor product Bernstein basis.Johnson and Cohen (2005) [11] present a robust search for distance extrema from a point to a curve or a surface.The robustness comes from using geometric operations with tangent cones rather than numerical methods to find all local extrema.Limaien and Trochu (1995) [12] compute the orthogonal projection of a point onto parametric curves and surfaces by constructing an auxiliary function and finding its zeros.Polak et al. (2003) [13] present a new feedback precision-adjustment rule with a smoothing technique and standard unconstrained minimization algorithms in the solution of finite minimax problems.Patrikalakis and Maekawa (2001) [14] reduce the distance function problem to solving systems of nonlinear polynomial equations.Based on Ma et al. [1], Selimovic (2006) [15] presents improved algorithms for the projection of points on NURBS curves and surfaces.Cohen et al. (1980) [16] provide classical subdivision algorithms which have been widely applied in computer-aided geometric design, computer graphics, and numerical analysis.Based on the subdividing concept [16], Piegl and Tiller (1995) [17] present an algorithm for point projection on NURBS surfaces by subdividing a NURBS surface into quadrilaterals, projecting the test point onto the closest quadrilateral, and then recover the parameter from the closest quadrilateral.Liu et al. (2009) [18] propose a local surface approximation technique-torus patch approximation-and prove that the approximation torus patch and the original surface are second-order osculating.By using the tangent line and the rectifying plane, Li et al. (2013) [19] present an algorithm for finding the intersection between the two spatial parametric curves.Hu et al. (2005) [8] use curvature iteration information for solving the projection problem.Scholars (Ku-Jin Kim (2003) [20], Li et al. (2004) [21], Chen et al. (2010) [22], Chen et al. (2009) [23], Bharath Ram Sundar et al. (2014) [24]) have fully analyzed and discussed the intersection curve between two surfaces, the minimum distance between two curves, the minimum distance between curve and surface and the minimum distance between two surfaces.By clipping circle technology, Chen et al. (2008) [25] provide a method for computing the minimum distance between a point and a NURBS curve.Then, based on clipping sphere strategy, Chen et al. (2009) [26] propose a method for computing the minimum distance between a point and a clamped B-spline surface.Being analogous to [25,26], based on an efficient culling technique that eliminates redundant curves and surfaces which obviously contain no projection from the given point, Young-Taek Oh et al. (2012) [27] present an efficient algorithm for projecting a given point to its closest point on a family of freeform curves and surfaces.Song et al. (2011) [7] propose an algorithm for calculating the orthogonal projection of parametric curves onto B-spline surfaces.It uses a second-order tracing method to construct a polyline to approximate the pre-image curve of the orthogonal projection curve in the parametric domain of the base surface.Regarding the projection problem, Kwanghee Ko et al. (2014) [28] give a detailed review on literatures before 2014.To sum up, those algorithms employ various techniques such as turning the problem into finding the root of the system of nonlinear equations, geometric methods, subdivision methods and circular clipping algorithms.It is well known that there are three classical first-order algorithms for computing the minimum distance between a point and a parametric surface [29][30][31].However, they did not prove convergence for the First-Order method.In this paper, we contribute in two aspects.Firstly, we prove the method's first-order convergence and its independence of the initial value.We also give some numerical examples to illustrate its faster convergence than the existing ones.Secondly, for some special cases where the First-Order method does not converge, we combine it with Newton's second-order iterative method to present the hybrid second-order algorithm.Our method essentially exploits hybrid iteration, thus it performs very well with a second-order convergence, it is faster than the existing methods, and it is independent of the initial value.Some numerical examples confirm our conclusion.
The rest of this paper is organized as follows.Section 2 presents a convergence analysis for the First-Order method for orthogonal projection onto a parametric surface.In Section 3, some numerical examples illustrate that it converges faster than the existing methods.In Section 4, for some special cases where the First-Order method is not convergent, an improved hybrid second-order algorithm is presented.Convergence analysis and experimental results for the hybrid second-order algorithm are also presented in this section.Finally, Section 5 concludes the paper.

Convergence Analysis for the First-Order Method
In this section, we prove that the method defined by ( 5) is the first-order convergent and its convergence is independent of the initial value.
Assume a regular parametric surface s(u, v) = ( f 1 (u, v), f 2 (u, v), f 3 (u, v)), i.e., p = (p 1 , p 2 , p 3 ) is a test point.The first-order geometric iteration [29][30][31] to compute the footpoint q of the test point p is the following.Projecting the test point p onto a tangent plane of the parametric surface s(u, v) at (u, v) = (u n , v n ) yields a point q determined by s(u n , v n ) and the partial derivatives s u (u n , v n ), s v (u n , v n ).(see the following Formulas (2) and ( 3)).The footpoint can be approximatively expressed in the following way where Multiplying with s u (u n , v n ) and s v (u n , v n ), respectively, we obtain where x, y denotes the scalar product of vectors x, y ∈ R 3 , and x denotes the norm of a vector x.
The corresponding solution of a regular system of linear Equation ( 4) is where We update u n , v n by adding ∆u, ∆v, and repeat the above procedure until |∆u| < ε and |∆v| < ε where ε is a given tolerance.This is the first-order geometric iteration method in [29][30][31] (See Figure 1).Furthermore, convergence of the First-Order method is independent of the initial value.Theorem 1.The method defined by (5) is the first-order convergent.Convergence of the iterative Formula (5) is independent of the initial value.
Proof.We firstly derive the expression of footpoint q.Assume that parameter (α, β) takes the value so that the test point p = (p 1 , p 2 , p 3 ) orthogonally projects onto the parametric surface s(u, v) = ( f 1 (u, v), f 2 (u, v), f 3 (u, v)).It is not difficult to show that there is a relational expression where h = ( The relationship can also be expressed in the following way, Because the footpoint q is the intersection of the tangent plane of the parametric surface s(u, v) at (u, v) = (u n , v n ) and the perpendicular line determined by test point p, the equation of the tangent plane of the parametric surface s At the same time, the vector of the line segment connected by the test point p and the point Because the vector (9) and the tangent vectors are orthogonal to each other, respectively, we naturally obtain a system of nonlinear equations of the parameters u, v, So the corresponding solutions of parameters u, v of the Formula (10) are where Substituting ( 11) into ( 8), and simplifying, we have So the footpoint q = (q 1 , q 2 , q 3 ) is the Formula (12).Substituting ( 12) into (5), and simplifying, we obtain the relationship, where Using Taylor's expansion, we obtain where Thus, we have and From ( 7) and ( 14), then we have By ( 14)-( 16), and using Taylor's expansion, Formula (13) can be transformed into the following form, where From (17), we can obtain uA = 0 and vA = 0.So Formula ( 18) can be transformed into the following form where e 1(n+1 . Formula (20) can be further simplified into the following form, where The result of Formula ( 21) implies that the iterative Formula ( 5) is the first-order convergent.In the following, we continue to illustrate that convergence of the iterative Formula ( 5) is independent of the initial value(See Figure 2).
Graphic demonstration for convergence analysis.
Our proof method is analogous to those methods in references [32,33].We project all points, curves and the surface of Figure 1 onto the y − z plane, and this yields Figure 2. When the iterative Formula (5) starts to iterate, according to the graphical demonstration, the corresponding parametric value of the footpoint q is (u n+1 , v n+1 ).The middle point of point (u n+1 , v n+1 ) and point ).From the graphic demonstration, there are two inequality , where e n 2 = (e 1n − α) 2 + (e 2n − β) 2 .To sum up, we can verify that convergence of the iterative Formula ( 5) is independent of the initial value.

Numerical Examples
Example 1.There is a general parametric surface s(u, v) = ( Using the First-Order method, the corresponding orthogonal projection parametric value is (α, β) =(1.9142974355881810, −0.80714871779409050), the initial iterative values (u 0 , v 0 ) are (1, −2), (−2, 2), (−2, −2), (0, 0), (1, 1), and(−2, −1), respectively.Each initial iterative value repeatedly iterates 10 times, respectively, yielding 10 different iteration times in the time unit of nanoseconds.In Table 1, the mean running time of the first-order iterative method is 134,760, 141,798, 41,033, 140,051, 137,059, and 42,399 nanoseconds for six different initial iterative values, respectively.In the end, the overall average running time in Table 1 is 106183.33nanoseconds (≈0.10618 ms), while the overall average running time of Tables 1 and  2 in [18] is 0.3565 ms.So the First-Order method is faster than the algorithm in [18].(See Figure 3).Example 2. There is a quasi-B-spline parametric surface s(u, v) = ( with a test point p = (p 1 , p 2 , p 3 ) = (15.0,20.0, 25.0).Using the First-Order method, the corresponding orthogonal projection parametric value is (α, β) =(1.0199334308624865, 1.3569112459785527), the initial iterative values (u 0 , v 0 ) are (1, 1), (2, 2), (2, 0), (1, 2), (2, 1), and(0, 2.0), respectively.Each initial iterative value repeatedly iterates 10 times, respectively, yielding 10 different iteration times in nanoseconds.In Table 2, the mean running time of the First-Order method is 500,831, 440,815, 480,969, 445,755, 426,737, and 488,092 nanoseconds for six different initial iterative values, respectively.In the end, the overall average running time is 463,866.50nanoseconds (≈0.463865 ms), while the overall average running time of the Example 2 for the algorithm in [18] is 1.705624977 ms under the same initial iteration condition.So the First-Order method is faster than the algorithm in [18].(See Figure 4).In Tables 1 and 2, the convergence tolerance required such that (u The approximate zero (α, β) found up to the 17th decimal place is displayed.All computations were done under g++ in Fedora linux 8 environment using our personal computer with T2080 1.73 GHz CPU and 2.5 GB memory.The overall average running time of 0.336222 ms in Examples 1 and 2 for their algorithm [18] means the First-Order method is faster than the algorithm in [18].In [18], authors point out that their algorithm is faster than that in [8], which means the First-Order method is faster than the one in [8].At the same time, the overall average running time of 61.81167 ms in three Examples for their algorithm [26] is obtained, then the First-Order method is faster than the algorithm in [26].However, in [26], the authors indicate that their algorithm is faster than those in [1,15], which means the First-Order method is faster than the one in [1,15].To sum up, the First-Order method converges faster than the existing methods in [1,8,15,18,26]

Counterexamples
In Sections 2 and 3, convergence of the First-Order method is independent of the initial value and some numerical examples illustrate that it converges faster than the existing methods.To show some special cases where the First-Order method is not convergent, we create five counterexamples as follows.

The Improved Algorithm
Since the First-Order method is not convergent for some special cases, we present the improved algorithm which will converge for any parametric surface, test point and initial iterative value.For simplicity, we briefly write (u, v) T as t, namely, t = (u, v) T , t n = (u n , v n ) T .It is well known that the classic iterative method to solve a system of nonlinear equations is Newton's method, whose iterative expression is where the system F(t) = 0 is expressed as follows: Then, the more specific expression of Newton's iterative Formula ( 22) is the following, where This method is quadratically convergent in a neighborhood of the solution where the Jacobian matrix is non-singular.Its convergence rate is greater than that of the First-Order method.However, a limitation of this method is that its convergence depends on the choice of the initial value.Only if the convergence condition for Newton's second-order iterative method is satisfied, is the method effective.In order to improve the robustness of convergence, based on the First-Order convergence method, we present Algorithm 1 (the hybrid second-order algorithm) for orthogonal projection onto the parametric surface.The hybrid second-order algorithm essentially combines the respective advantages of these two methods: if the iterative parametric value from the First-Order geometric iteration method satisfies the convergence condition for Newton's second-order iterative method, we then turn to the method to accelerate the convergence process.If not, we continue the First-Order geometric iteration method until its iterative parametric value satisfies the convergence condition for Newton's second-order iterative method, and we then turn to it as above.Then comes the end of the whole algorithm.The procedure that we proposed guarantees that the whole iterative process is convergent, and its convergence is independent of the choice of the initial value.In addition, the procedure can increase its rapidity.The hybrid second-order iterative algorithm for computing the minimum distance between a point and a parametric surface can be realized as follows.
Algorithm 1: The hybrid second-order algorithm.
Input: Input the initial iterative parametric value t 0 , parametric surface s(u, v) and test point p.
Output: Output the final iterative parametric value.
Step 1. Input the initial iterative parametric value t 0 .
Step 2. Use the iterative Formula (5), compute the parametric incremental value ∆t, and update t 0 + ∆t to t 0 , namely, t 0 = t 0 + ∆t.Step 3. Judge whether the norm of difference between the former t 0 and the latter t 0 is near 0( ∆t < ε).If so, end this algorithm.
Step 4. Substitute new t 0 into { Using Newton's second-order iterative Formula ( 22), compute until the norm of difference between the former t 0 and the latter t 0 is near 0( ∆t < ε), then end this algorithm.

}
Remark 1.We give a geometric interpretation of Newton's iterative method with two variables in Figure 5, where abscissa axis t represents the two-dimensional coordinate (u, v), the vertical coordinate z represents function value F 1 (t) or F 2 (t).Curve F 1 (t) and F 2 (t) actually denote two surfaces in three-dimensional space, which are determined by the first and the second Formula in Equation (23), respectively.γ is the intersection point of two lines, the first of which is the intersection line of F 1 (t) and plane t, and the second of which is the intersection line of F 2 (t) and plane t.In another way, γ is the solution to Equation (23).Through point t 0 , draw a line perpendicular to the plane t, and then the perpendicular line intersects the surface F 1 (t) and F 2 (t), respectively.The intersection points are designated as the first and the second intersection point, respectively.Through the first and the second intersection point, we make two tangent planes, respectively.These two tangent planes intersect with plane t and generate two intersection lines, which are designated as the first and the second intersection lines.Two intersection lines intersect to obtain the point t 1 .Given the intersection point t 1 , we repeat the procedure above to obtain the new first and the second intersection point, and then the new first and the second intersection lines.Finally, we can obtain the new intersection point t 2 .Repeat the above steps until the iterative value converges to γ=(α, β).Of course, it must satisfy the condition for the fixed point theorem of Newton's iterative method before Newton's iterative method starts to work.
Graphic demonstration for the hybrid second-order algorithm.
Remark 2. For some special cases, in Section 4, where the First-Order method is not convergent, our hybrid second-order algorithm could converge.In addition, we find that the method is convergent for any initial iterative value, any test point and parametric surface in many tested examples.

Convergence Analysis for the Improved Algorithm
Definition 1.In Reference [34] Theorem 2. In Reference [34] (Fixed Point Theorem), let D = {(x 1 , x 2 , ..., Then, these sequence {x k } ∞ k=0 defined by an arbitrarily selected x 0 in D and generated by This converges to the unique fixed point p ∈ D and satisfies the iterative relationship In the following, we directly present the fixed point theorem for Newton's iterative method. for some collection of constants a 1 , a 2 and b 1 , b 2 .Suppose G (see Equation ( 24)) is a continuous function from D ⊂ R 2 into R < L 2 , whenever t ∈ D, for each j = 1, 2 and each component function G i .Then, these sequence {t k } ∞ k=0 defined by an arbitrarily selected t 0 in D and generated by This converges to the unique fixed point p ∈ D and satisfies the iterative relationship Theorem 4. The convergence order of the hybrid second-order algorithm for orthogonal projection onto the parametric surface is 2. Convergence of the hybrid second-order algorithm is independent of the initial value.
Proof.Let γ = (α, β) T be a root of system F(t) = 0, t = (u, v) T , where the system F(t) = 0 is expressed by Formula (23).Define e n = t n − γ and use Taylor's expansion around γ, we have Since γ is a simple zero of the Formula ( 23), so F(γ) = 0.Then, we have where Using (21), and ( 29)-( 31), we obtain The result means that the hybrid algorithm is second-order convergent.According to the procedure of the hybrid second-order algorithm, if the iterative parametric value satisfies the convergence condition for Newton's second-order iterative method, we then turn to Newton's second-order iterative method.If not, we continue the First-Order geometric iteration method until its iterative parametric value satisfies the convergence condition of Newton's second-order iterative method.Then comes the end of the whole algorithm.When the hybrid second-order iterative algorithm implements the First-Order iterative method, its independence of the initial value can be guaranteed by Theorem 1.When the hybrid second-order iterative algorithm implements Newton's second-order iterative method, its independence of the initial value can be guaranteed by Theorem 3. To sum up, convergence of the hybrid second-order iterative algorithm is independent of the initial value throughout the whole algorithm operating process.

Numerical Experiments
Example 3. We replicate Example 1 using the hybrid second-order algorithm in Table 3 and compare it with results in Table 1.In Table 3, the mean running time of the hybrid second-order algorithm is 303,210.7,319,045.5, 92,325.2,315,116.9,308,383.9, and 95,399.5 nanoseconds for six different initial iterative values, respectively.In the end, the overall average running time in Table 3 is 238,913.62nanoseconds (≈0.2389 ms), while the overall average running time of Tables 1 and 2 in [18] is 0.3565 ms.So our hybrid second-order algorithm is faster than the algorithm in [18].On the other hand, from Table 3 = 0.1180339887.In the end, the overall average running time is 1,043,700 nanoseconds(≈1.0437ms), while the overall average running time of the Example 2 for the algorithm in [18] is 1.705624977 ms under the same initial iteration condition.So the hybrid second-order algorithm is faster than the algorithm proposed in [18].On the other hand, from Table 4, In Tables 3 and 4, the convergence tolerance, the decimal place of approximate zero and the computation environment are set as those in Tables 1 and 2. The overall average running time of 0.336222 ms in Examples 3 and 4 for the algorithm in [18] implies that our hybrid second-order algorithm is faster than the algorithm in [18] .In [18], the authors point out that their algorithm is faster than the algorithm in [8], which means the hybrid second-order algorithm is faster than the one in [8].At the same time, the overall average running time of 61.81167 ms in three Examples for their algorithm implies that our hybrid second-order algorithm is faster than the algorithm in [26].However, in [26], the authors indicate that their algorithm is faster than the ones in [1,15], which means our hybrid second-order algorithm is faster than the ones in [1,15].To sum up, with the exception of the First-Order method, our hybrid second-order algorithm converges faster than the existing methods in [1,8,15,18,26].In Tables 3 and 4, the unit of time is nanoseconds.Remark 3. The hybrid second-order algorithm presented can be applied for one projection point case for orthogonal projection of a test point onto a parametric surface s(u, v).For the multiple orthogonal projection points situation, the basic idea of our approach is as follows: (2) Randomly select an initial iterative parametric value in each sub-region.
(3) For each initial iterative parametric value, holding other initial iterative parametric values fixed, use the hybrid second-order iterative algorithm and iterate until it converges.Let us assume that the converged iterative parametric values are (α i , β j ), i, j = 0, 1, 2, ..., M -1, respectively.(4) Compute the local minimum distances d ij , i, j = 0, 1, 2, ..., M -1, where d ij = p − s(α i , β j ) .(u 0 , v 0 )  ( This indicates that these two inequality relationships satisfy Formula ( 6) or (7).Thus, it illustrates that convergence of the hybrid second-order algorithm is independent of the initial value.Furthermore, the hybrid second-order algorithm is robust and efficient as shown by the previous two of ten challenges proposed by [35].

Conclusions
This paper investigates the problem related to a point projection onto a parametric surface.To compute the minimum distance between a point and a parametric surface, three well-known first-order algorithms have been proposed by [29][30][31] (hereafter, the First-Order method).In this paper, we prove the method's first-order convergence and its independence of the initial value.We also give some numerical examples to illustrate its faster convergence than the existing ones.For some special cases where the First-Order method does not converge, we combine it with Newton's second-order iterative method to present the hybrid second-order algorithm.Our method essentially exploits hybrid iteration, thus it performs very well with a second-order convergence, it is faster than the existing methods and it is independent of the initial value.Some numerical examples confirm our conclusion.An area for future research is to develop a method for computing the minimum distance between a point and a higher-dimensional parametric surface.

Figure 1 .
Figure 1.Graphic demonstration of the First-Order method.
2, ..., n} for some collection of constants a 1 , a 2 , ..., a n and b 1 , b 2 , ..., b n .Suppose G is a continuous function from D ⊂ R n into R n with the property that G(x) ∈ D whenever x ∈ D.Then, G has a fixed point in D.Moreover, suppose that all the component functions of G have continuous partial derivatives and a constant 0 < L < 1 exists with ∂g i (x) ∂x j ≤ L n , whenever x ∈ D, for each j = 1, 2, ..., n and each component function g i .

Remark 4 .p
In addition to two test examples, we have also tested many other examples.According to these test results, for different initial iterative values, it can converge to the corresponding orthogonal projection point by using the hybrid second-order algorithm.Namely, if the initial iterative value is (u0 , v 0 ) ∈ [a, b] × [c, d],which belongs to the parametric region of a parametric surface s(u, v), the corresponding orthogonal projection parametric value for the test point p = (p 1 , p 2 , p 3 ) is (α, β), and then the test point p and its corresponding orthogonal projection parametric value (α, β) satisfy two inequality relationships − s(α, β), ∂s(α, β) ∂u < 1E − 16, p − s(α, β), ∂s(α, β) ∂v < 1E − 16.

Table 1 .
Running time (in nanoseconds) for different initial iterative values by the First-Order method.

Table 2 .
Running time for (in nanoseconds) for different initial iterative values by the First-Order method.

Table 3 .
Running time (in nanoseconds) for different initial iterative values by the hybrid second-order algorithm.Example 4. We replicate Example 2 using the hybrid second-order algorithm in Table4and compare it with results in

Table 2 .
In Table 4, the mean running time of the hybrid second-order algorithm is 1,126,871, 991,834.4,1,082,181, 1,002,949, 960,158.4,and 1,098,208 nanoseconds for six different initial iterative values, respectively.

Table 4 .
Running time (in nanoseconds) for different initial iterative values using the hybrid second-order algorithm.