Integrated Hybrid Second Order Algorithm for Orthogonal Projection onto a Planar Implicit Curve

The computation of the minimum distance between a point and a planar implicit curve is a very important problem in geometric modeling and graphics. An integrated hybrid second order algorithm to facilitate the computation is presented. The proofs indicate that the convergence of the algorithm is independent of the initial value and demonstrate that its convergence order is up to two. Some numerical examples further confirm that the algorithm is more robust and efficient than the existing methods.


Introduction
Due to its great properties, the implicit curve has many applications.As a result, how to render implicit curves and surfaces is an important topic in computer graphics [1], which usually adopts four techniques: (1) representation conversion; (2) curve tracking; (3) space subdivision; and (4) symbolic computation.Using approximate distance tests to replace the Euclidean distance test, a practical rendering algorithm is proposed to rasterize algebraic curves in [2].Employing the idea that field functions can be combined both on their values and gradients, a set of binary composition operators is developed to tackle four major problems in constructive modeling in [3].As a powerful tool for implicit shape modeling, a new type of bivariate spline function is applied in [4], and it can be created from any given set of 2D polygons that divides the 2D plane into any required degree of smoothness.Furthermore, the spline basis functions created by the proposed procedure are piecewise polynomials and explicit in an analytical form.
Aside from rendering of computer graphics, implicit curves also play an important role in other aspects of computer graphics.To facilitate applications, it is important to compute the intersection of parametric and algebraic curves.Elimination theory and matrix determinant expression of the resultant in the intersection equations are used in [5].Some researchers try to transform the problem of intersection into that of computing the eigenvalues and eigenvectors of a numeric matrix.Similar to elimination theory and matrix determinant expression, combining the marching methods with the algebraic formulation generates an efficient algorithm to compute the intersection of algebraic and NURBSsurfaces in [6].For the cases with a degenerate intersection of two quadric surfaces, which are frequently applied in geometric and solid modeling, a simple method is proposed to determine the conic types without actually computing the intersection and to enumerate all possible conic types in [7].M.Aizenshtein et al. [8] present a solver to robustly solve well-constrained n × n transcendental systems, which applies to curve-curve, curve-surface intersections, ray-trap and geometric constraint problems.
To improve implicit modeling, many techniques have been developed to compute the distance between a point and an implicit curve or surface.In order to compute the bounded Hausdorff distance between two real space algebraic curves, a theoretical result can reduce the bound of the Hausdorff distance of algebraic curves from the spatial to the planar case in [9].Ron [10] discusses and analyzes formulas to calculate the curvature of implicit planar curves, the curvature and torsion of implicit space curves and the mean and Gaussian curvature for implicit surfaces, as well as curvature formulas to higher dimensions.Using parametric approximation of an implicit curve or surface, Thomas et al. [11] introduce a relatively small number of low-degree curve segments or surface patches to approximate an implicit curve or surface accurately and further constructs monoid curves and surfaces after eliminating the undesirable singularities and the undesirable branches normally associated with implicit representation.Slightly different from ref. [11], Eva et al. [12] use support function representation to identify and approximate monotonous segments of algebraic curves.Anderson et al. [13] present an efficient and robust algorithm to compute the foot points for planar implicit curves.
Contribution: An integrated hybrid second order algorithm is presented for orthogonal projection onto planar implicit curves.For any test point p, any planar implicit curve with or without singular points and any order of the planar implicit curve, any distance between the test point and the planar implicit curve, the algorithm could be convergent.It consists of two parts: the hybrid second order algorithm and the initial iterative value estimation algorithm.
The hybrid second order algorithm fuses the three basic ideas: (1) the tangent line orthogonal iteration method with one correction; (2) the steepest descent method to force the iteration point to fall on the planar implicit curve as much as it can; (3) Newton-Raphson's iterative method to accelerate iteration.
Therefore, the hybrid second order algorithm is composed of six steps.The first step uses the steepest descent method of Newton's iterative method to force the iterative value of the initial value to lie on the planar implicit curve, which is not associated with the test point p.In the second step, Newton's iterative method employs the relationship determined by the test point p to accelerate the iteration process.The third step finds the orthogonal projection point q on the tangent line, which goes through the initial iterative point, of a test point p.The fourth step gets the linear orthogonal increment value.The same relationship in the second step is used once more to accelerate the iteration process in the fifth step.The final step gives some correction to the result of the iterative value in the fourth and fifth step.
One problem for the hybrid second order algorithm is that it appears divergent if the test point p lies particularly far away from the planar implicit curve.Since it has been found that when the initial iterative point is close to the orthogonal projection point p Γ , no matter how far away the test point p is from the planar implicit curve, it will be convergent, an algorithm, named the initial iterative value estimation algorithm, is proposed to drive the initial iterative value toward the orthogonal projection point p Γ as much as possible.Accordingly, the second order algorithm with the initial iterative value estimation algorithm is named as the integrated hybrid second order algorithm.
The rest of this paper is organized as follows.Section 2 presents related work for orthogonal projection onto the planar implicit curve.Section 3 presents the integrated hybrid second order algorithm for orthogonal projection onto the planar implicit curve.In Section 4, convergent analysis for the integrated hybrid second order algorithm is described.The experimental results including the evaluation of performance data are given in Section 5. Finally, Sections 6 and 7 conclude the paper.

Related Work
The existing methods can be divided into three categories: local methods, global methods and compromise methods between local and global methods.

Local Methods
The first one is Newton's iterative method.Let point x = (x, y) be on a plane, and let Γ : f (x) = 0 be a smooth planar implicit curve; its specific form can be represented as: Let p = (p 1 , p 2 ) be a point in the vicinity of curve Γ.The orthogonal projection point p Γ satisfies the relationships: where ∧ is the difference-product ( [14]).The nonlinear system (2) can be solved using Newton's iterative method [15]: where ) is the Jacobian matrix of partial derivatives of L(x) with respect to x. Sullivan et al. [16] used Lagrange multipliers and Newton's algorithm to compute the closest point on the curve for each point.
where λ is the Lagrange multiplier.It will converge after repeated iteration of Equation ( 4) for increment δx, δy, δλ with the initial iterative point and λ = 0.However, Newton's iterative method or Newton-type's iterative method is of local convergence, i.e., it sometimes failed to converge even with a reasonably good initial guess.On the other hand, once a good initial guess lies in the vicinity of the solution, two advantages emerge: its fast convergence speed and high convergence accuracy.In this paper, these two advantages to improve the accuracy and effectiveness of convergence for the integrated hybrid second order algorithm will be employed.The second one is the homotopy method [17,18].In order to solve the target system of nonlinear Equation (2), they start with nonlinear equations g(x) = 0, where g(x) = g 1 (x, y) = 0, g 2 (x, y) = 0, and give a homotopy formula: where t is a continuous parameter and ranges from 0-1, deforming the starting system of nonlinear equations g(x) into the target system of nonlinear equations L(x).The numerical continuous homotopy methods can compute all isolated solutions of a polynomial system and are globally convergent and exhaustive solvers, i.e., their robustness is summarized by [19], and their high computational cost is confirmed in [20].

Global Methods
Firstly, a global method is a resultants' one.For the algebraic curve with low degree no more than quartic, the resultant methods are a good choice.With classical elimination theory, it will yield a resultant polynomial from two polynomial equations with two unknown variables where the roots of the resultant polynomial in one variable correspond to the solution of the two simultaneous equations [21][22][23][24].Assume two polynomial equations f (x, y) = 0 and g(x, y) = 0 with respective degree m and n.Let p(p ≤ m) and q(q ≤ n) are two integers, respectively.To facilitate the resultant method calculation, y is a constant in this case.It indicates that: where α is the zero of the x-coordinate and (x, y) is a common solution of f = 0 and g = 0, δ(x, α) = f (x,y)g(α,y)− f (α,y)g(x,y) x−α , l = max(m, n) − 1 and B is the Bézout matrix of f (x, y) and g(x, y) (with elements consisting of polynomials in y) that has nothing to do with variables α and x.Therefore, the determinant det(B(y)) = 0 of the Bézout matrix is a polynomial in y, and there are at most mn intersections of f = 0 and g = 0. Therefore, the roots of this polynomial give the y-coordinates of mn possible closest points on the curve.The best known results, such as the Sylvester's resultant and Cayley's statement of Bézout's method [21,22,24], obviously indicate that, if the algebraic curve has a high degree more than quintic, it is very difficult to use the resultant method to solve a two-polynomial system.
Secondly, the global method uses the Bézier clipping technique.Using the convex hull property of Bernstein-Bézier representations, footpoints can be found by solving the nonlinear system of Equation ( 2) [25][26][27].Transformation of (2) into Bernstein-Bézier form eliminates parts of the domains outside the convex hull box not including a solution.Elimination rules are repeatedly applied using the de Casteljau subdivision algorithm.Once it meets a certain accuracy requirement, the algorithm will end.The Bézier clipping method can find all solutions of the system (2), especially the singular points on the implicit curve.Certainly, the Bézier clipping method is of global convergence, but with relatively expensive computation due to many subdivision steps.
Based on and more efficient than [26], a hybrid parallel algorithm in [28] is proposed to solve systems with multivariate constraints by exploiting both the CPU and the GPU multi-core architectures.In addition, their GPU-based subdivision method utilizes the inherent parallelism in the multivariate polynomial subdivision.Their hybrid parallel algorithm can be geometrically applied and improves the performance greatly compared to the state-of-the-art subdivision-based CPU.Two blending schemes presented in [29] efficiently eliminate domains without any root and therefore greatly cut down the number of subdivisions.Through a simple linear blend of functions of the given polynomial system, this seek function would satisfy two conditions: no-root contributing and exhausting all control points of its Bernstein-Bézier representation of the same sign.It continually keeps generating these functions so as to eliminate the no-root domain during the subdivision process.
Van Sosin et al. [30] decompose and efficiently solve a wide variety of complex piecewise polynomial constraint systems with both zero constraints and inequality constraints with zero-dimensional or univariate solution spaces.The algorithm contains two parts: a subdivision-based polynomial solver and a decomposition algorithm, which can deal with large complex systems.It confirms that its performance is more effective than the existing ones.

Compromise Methods between Local and Global Methods
Firstly, the compromise method lies between the local and global method and uses successive tangent approximation techniques.The geometrically iterative method for orthogonal projection onto the implicit curve presented by Hartmann [31,32] consists of two steps for the whole process of iteration, while the local approximation of the curve by its tangent or tangent parabola is the key step.
where f (x) = f (x, y).The first step in [31,32] repeatedly iterates the iterative Formula (7) in the steepest way such that the iterative point is as close as possible to the curve f (x) = 0 with arbitrary initial iterative point.This is called the first step of [31,32].Then, the second step of [31,32] will get the vertical point q by the iterative Formula (8).
Repeatedly run these two steps until the vertical point q falls on the curve f (x).Unfortunately, the successive tangent approximation method fails for the planar implicit curves.
Secondly, the compromise method between the local and global method uses the successive circular approximation technique.Similar to [31,32], Nicholas [33] uses the osculating circle to develop another geometric iteration method.For a planar implicit curve f (x) = 0, the curvature formula of a point on an implicit curve could be defined as K (see [10]), and the radius r of curvature could be expressed as: The osculating circle determined by the point x k has its center at x r where the radius r is from Formula (9).For the current point x k from the orthogonal projection point p Γ on the curve to the test point p, then the next iterative point x k+1 will be the intersection point determined by the line segment px r and the current osculating circle.Repeatedly iterate until the distance between the former iterative point x k and the latter iterative point x k+1 is almost zero.The geometric iteration method [33] may fail in three cases in which there is no intersection, or the intersection is difficult to solve when the algebraic curve is of high degree more than quintic, or the new estimated iterative point lies very far from the planar implicit curve.The third geometric iteration method uses the curvature information to orthogonally project point onto the parametric osculating circle and osculating sphere [34].Although the method in [34] handles point orthogonal projection onto the parametric curve and surface, the basic idea is the same as that in [33] to deal with the problem for the planar implicit curve.The convergence analysis for the method in [34] is provided in [35].The third geometric iteration method [34] is more robust than the existing methods, but it is time consuming.
Thirdly, the compromise method between the local and global method also uses the circle shrinking technique [14].It repeatedly iterates Equation (7) in the steepest way such that the iterative point is as close as possible to the curve f (x) with the arbitrary initial iterative point.This time, the iterative point falling on the curve is called the point p c , and then, two points p and p c define a circle with center p. Compute the point p + by calculating the (local) maximum of curve f (x) along the circle with center p, starting from p c .The intersection between the line segment pp + and the planar implicit curve f (x) will be the next iterative point p' c .Replace the point p c with the point p' c , repeat the above process until the distance between these two points approaches zero.Hu et al. [36] proposed a circle double-and-bisect algorithm to reliably evaluate the accurate geometric distance between a point and an implicit curve.The circle doubling algorithm begins with a very small circle centered at the test point p.Extend the circle by doubling its radius if the circle does not intersect with the implicit curve f (x).The extending process continues until the circle intersects with the implicit circle where the former circle does not hit the curve, but the current one does.At the same time, the former radius and the current radius are called interior radius r 1 and exterior radius r 2 , respectively.The bisecting process continues yielding new radius r = r 1 +r 2 2 .If a circle with its radius r intersects with the curve, replace r with r 2 , else with r 1 .Repeatedly iterate the above procedure until |r 1 − r 2 | < ε.Similar to the circle shrinking method, Chen et al. [37,38] made some contribution to the orthogonal projection onto the parametric curve and surface.Given a test point p = (p 1 , p 2 ) and a planar implicit curve f (x) = 0, the footpoint p Γ has to be a solution of a 2 × 2 well-constrained system: where this formula is from [38].The efficient algebraic solvers can solve this system, and one just needs to take the minimum over all possible footpoints.It uses a circular/spherical clipping technique to eliminate the curve parts/surface patches outside a circle/sphere with the test point as its center point, where the objective squared distance function for judging whether a curve/surface is outside a circle/sphere is the key technique.The radius of the elimination circle/sphere gets smaller and smaller during the subdivision process.Once the radius of the circle/sphere can no longer become smaller, the iteration will end.With the advantage of high robustness, the algorithm still faces two difficulties: it is time consuming and had difficulty calculating the intersection between the circle and the planar implicit curve with a relatively high degree (more than five).

Integrated Hybrid Second Order Algorithm
Let Γ : f (x) = f (x, y) = 0 be a smooth planar implicit curve, and let p = (p 1 , p 2 ) be a point in the vicinity of curve Γ (test point).Assume that s is the arc length parameter for the planar implicit curve Γ : f (x) = f (x, y) = 0. t = [dx/ds, dy/ds] is the tangent vector along the implicit curve Γ : f (x) = 0.The orthogonal projection point p Γ to satisfy this relationship: where ∧ is the difference-product ( [14]).

Orthogonal Tangent Vector Method
The derivative of the planar implicit curve f (x) with respect to parameter s is, where ∇ = ∂ ∂x , ∂ ∂y is the Hamiltonian operator and the symbol is the inner product.Its geometric meaning is that the tangent vector t is orthogonal to the corresponding gradient ∇ f .The combination of the tangent vector t and Formula (12) will generate: From (13), it is not difficult to know that the unit tangent vector of t is: The following first order iterative algorithm determines the foot point of p on Γ.
where t 0 = − f y , f x ∇ f , ∆s = q − x n .q is the corresponding orthogonal projection point of test point p at the tangent line determined by the initial iterative point x n (see Figure 1).Formula ( 14) can be expressed as, q = p − ( (p where x n is the initial iterative point.Many numerical tests illustrate that iterative Formula (15) depends on the initial iterative point, namely it is very difficult for the iterative value y n to fall on the planar implicit curve.
Graphic demonstration for the hybrid second order algorithm.

Steepest Descent Method
To move the iterative value y n to fall on the planar implicit curve f (x) as much as possible, a method of preprocessing is introduced.Before the implementation of the iterative Formula ( 15), the steepest descent method will be adopted, namely a basic Newton's iterative formula is added such that the iterative value y n falls on the planar implicit curve f (x) as much as possible.
where t 0 = − f y , f x ∇ f , ∆s = q − y n .

Linear Calibrating Method
Although more robust than the iterative Formula (15) to a certain extent, the iterative Formula (16) will often change convergence if the test point p or the initial iterative point x 0 takes different values.Especially for large ∆s, iterative point z n will deviate from the planar implicit curve greatly, namely In this case, a correction for the deviation of iterative point z n is proposed as follows.
is used for correction.That is to say, z n = z n + δz n , z n and z n are the iteration values before and after correction, respectively, and | f (z n )| < ε.The correction aims to make the deviation of the iteration value z n from the planar implicit curve as small as possible.Let δz n be perpendicular to increment value ∆z n = sign(∆s)t 0 ∆s and orthogonal to the planar implicit curve such that δz n , ∆z n = 0 and ∇ f , δz n = −∆e, where ∇ f and ∆e take the value at z n .Then, it is easy . The corresponding iterative formula for correction will be, where Obviously the stability and efficiency of the iterative Formula ( 17) improve greatly, compared with the previous iterative Formulas ( 15) and ( 16).

Newton's Accelerated Method
Many tests for the iterative Formula ( 17) conducted indicate that it is sometimes not convergent when the test point lies far from the planar implicit curve.Newton's accelerated method is then adopted to correct the problem.For the classic Newton second order iterative method, its iterative expression is: where ∇F 0 (x), ∇F 0 (x) is inner product of the gradient of the function F 0 (x) with itself.The function where the symbol [ ] denotes the determinant of a matrix (p − x) × ∇f(x).In order to improve the stability and rate of convergence, based on the iterative Formula ( 17), the hybrid second order algorithm is proposed to orthogonally project onto the planar implicit curve f (x).Between Step 1 and Step 2 of the iterative Formula ( 17) and between Step 3 and Step 4 of the same formula, the iterative Formula ( 18) is inserted twice.After this, the stability, the rapidity, the efficiency and the numerical iterative accuracy of the iterative algorithm (17) all improve.Then, the iterative formula becomes, where Iterative termination for the iterative Formula (20) satisfies: x n+1 − x n < ε.The robustness and the stability of the iterative Formula (20) improves, compared with the previous iteration formulas.
That is to say, even for test point p being far away from the planar implicit curve, the iterative Formula ( 20) is still convergent.After normalization of the second equation and the fifth equation in the iterative Formula (20), it becomes, where . The iterative Formula ( 21) can be implemented in six steps.The first step computes the point x n on the planar implicit curve using the basic Newton's iterative formula, which is not associated with test point p for any initial iterative point.The second step uses Newton's iterative method to accelerate the whole iteration process and get the new iterative point z n , which is associated with test point p.The third step gets the orthogonal projection point q (footpoint) at the tangent line to f (x).The fourth equation in iterative Formula (21) yields the new iterative point u n .The third step and the fourth step compute the linear orthogonal increment, which is the core component (including linear calibrating method of sixth step) of the iterative Formula ( 21).The fifth step accelerates the previous steps again and yields the iterative point, which is associated with test point p.The sixth step corrects the iterative result for the previous three steps.Therefore, the whole six steps ensure the robustness of the whole iteration process.The above procedure is repeated until the iterative point coincides with the orthogonal projection point p Γ (see Figure 1 and the detailed explanation of Remark 3).

Remark 1.
In the actual implementation of the iterative Formula ( 21) of the hybrid second order algorithm (Algorithm 1), three techniques are used to optimize the process.On the right-hand side of Step 1, Step 2, Step 3 and step 5, the part in parentheses is calculated firstly and then the part outside the parentheses to prevent overflow of the intermediate calculation process.Error handling for the second term is added in the right-hand side of Step 4 in the iterative Formula (21).Namely, if p − z n , t 0 = 0, sign( p − z n , t 0 ) = 1.For the second term of the right-hand side in Step 6 of the iterative Formula ( 21), if the determinant of ∇ f T , (∆v n ) T is zero, equals zero, then substitute the sixth step with x n+1 = v n to avoid the overflow problem.
According to the analyses above, the hybrid second order algorithm is presented as follows.
Algorithm 1: Hybrid second order algorithm.
Input: Initial iterative value x 0 , test point p and planar implicit curve f (x) = 0. Output: The orthogonal projection point p Γ .Description: Step 1: Remark 2. Many tests demonstrate that if the test point p is not far away from the planar implicit curve, Algorithm 1 will converge for any initial iterative point x 0 .For instance, assume a planar implicit curve f (x, y) = x 6 + 2x 5 y − 2x 3 y 2 + x 4 − y 3 + 2y 8 − 4 and four different test points (13, 7), (3, −4), (−2, 2), (−7, −3); Algorithm 1 converges efficiently for the given initial iterative value.See Table 1 for details, where p is the test point, x 0 is the initial iterative point, iterations is the number of iterations, | f (p Γ )| is the absolute function value with the orthogonal projection point p Γ and Error_2 =| However, when the test point p is far away from the planar implicit curve, no matter whether the initial iterative point p is close to the planar implicit curve, Algorithm 1 sometimes produces oscillation such that subsequent iterations could not ensure convergence.For example, for the same planar implicit curve with test point p = (17, −11) and initial iterative point x 0 = (2, −2), it constantly produces oscillation such that subsequent iterations could not ensure convergence after 838 iterations (see Table 2).

Initial Iterative Value Estimation Algorithm
Through Remark 2, when the test point p is not far away from the planar implicit curve, with the initial iterative point x 0 in any position, Algorithm 1 could ensure convergence.However, when the test point p is far away from the planar implicit curve, even if the initial iterative point x 0 is close to the planar implicit curve, Algorithm 1 sometimes produces oscillation such that subsequent iterations could not ensure convergence.This is essentially a problem for any Newton-based method.
Consider high nonlinearity, which cannot be captured just by f (x, y) = 0, i.e, the surface [x, y, f (x, y)] is very oscillatory in the neighborhood of the z = 0 plane, but it does intersect the z = 0 plane in a single closed branch.Under this case, for far-away test point p from the planar implicit curve, any Newton-based method sometimes will produce oscillation to cause non-convergence.We will give a counter example in Remark 2. To solve the problem of non-convergence, some method is proposed to put the initial iterative point x 0 close to the orthogonal projection point p Γ .Therefore, the task changes to construct an algorithm such that the initial iterative value x 0 of the iterative Formula ( 21) and the orthogonal projection point p Γ are as close as possible.The algorithm can be summarized as follows.Input an initial iterative point x 0 , and repeatedly iterate with the basic Newton's iterative formula y = x − ( f (x)/ ∇ f (x), ∇ f (x) ) ∇ f (x) such that the iterative point lies on the planar implicit curve f (x) (see Figure 2c).After that, iterate once through the formula where the blue point denotes the initial iterative value (see Figure 2c).After the first round iteration in Figure 2, then replace the initial iterative value with the iterated value q, and do the second round iteration (see Figure 3).After the second round iteration, replace the initial iterative value with the iterated value q, and do the third round iteration (see Figure 4).The detailed algorithm is the following.
Firstly, the notations for Figures 2-4 are clarified.Black and green points represent test point p and orthogonal projection point p Γ , respectively.The blue point denotes Step 2 in Algorithm 2, whether it is on the planar implicit curve f (x) or not.The footpoint q(red point) denotes q in Step 3 of Algorithm 2, and the brown curve describes the planar implicit curve f (x) = 0.
Step 2 in Algorithm 2 uses basic Newton's iterative method.That is to say, it repeatedly iterates using the steepest descent method in Section 3.2 until the blue point Step 2 lies on the planar implicit curve f (x).At the same time, through Step 3 in Algorithm 2, it yields footpoint q (see Figure 2).The integer n of the iteration round counts one after the first round iteration of Algorithm 2. When the blue point is on the planar implicit curve f (x), at this time, replace the initial iterative value with the iterated value q, and do the second round iteration; the integer n of the iteration round counts two (see Figure 3).Replace the initial iterative value with the iterated value q again after the second round iteration, and do the third round iteration; the integer n of the iteration round counts three (see Figure 4).When n = 3 in Step 4, then exit Algorithm 2. At this time, the current footpoint q from Algorithm 2 will be the initial iterative value for Algorithm 1.

Algorithm 2: Initial iterative value estimation algorithm.
Input: Initial iterative value x 0 , test point p and planar implicit curve f (x) = 0. Output: The footpoint point q.Description: Step 1: n = 0; x n+1 = x 0 ; Step 2: Thirdly, the reason for choosing n = 3 in Algorithm 2 is explained.Many cases are tested for planar implicit curves with no singular point.As long as n = 2, the output value from Algorithm 2 could be used as the initial iterative value of Algorithm 1 to get convergence.However, if the planar implicit curve has singular points or big fluctuation and oscillation appear, n = 3 can guarantee the convergence.In a future study, a more optimized and efficient algorithm needs to be developed to automatically specify the integer n.

Integrated Hybrid Second Order Algorithm
Algorithm 2 can optimize the initial iterative value for Algorithm 1.Then, Algorithm 1 can project the test point p onto planar implicit curve f (x).The integrated hybrid second order algorithm (Algorithm 3) is presented to take advantage of Algorithms 1 and 2, which are denoted as Algorithm 1 ((x 0 , p, f (x)) and Algorithm 2 (q, p, f (x)) for convenience, respectively.Algorithm 3 can be described as follows (see Figure 5).Firstly, the notations for Figure 5 are clarified, which describes the entire iterative process in Algorithm 3. The black point is test point p; the green point is orthogonal projection point p Γ ; the blue point is the left-hand side value of the equality of the first step of the iterative Formula (21) in Algorithm 1; footpoint q (red point) is the left-hand side value of the equality of the third step of the iterative Formula (21) in Algorithm 1; and the brown curve represents the planar implicit curve f (x).
Algorithm 3: Integrated hybrid second order algorithm.
Input: Initial iterative value x 0 , test point p and planar implicit curve f (x) = 0. Output: The orthogonal projection point p Γ .Description: Step 1: q = Algorithm 2 (x 0 , p, f (x)); Step 2: p Γ = Algorithm 1 (q, p, f (x)); Step 3: return p Γ ; Secondly, Algorithm 3 is interpreted.The output from Algorithm 2 is taken as the initial iterative value for Algorithm 1 (see footpoint q or the red point in Figure 4c).Algorithm 1 repeatedly iterates until it satisfies the termination criteria ( x n+1 − x n < ε) (see Figure 5).The six subgraphs in Figure 5 represent successive steps in the entire iterative process of Algorithm 1.In the end, three points of green, blue and red merge into orthogonal projection point p Γ (see Figure 5f).Remark 3. Algorithm 3 with two sub-algorithms is interpreted geometrically, where Algorithms 1 and 2 are graphically demonstrated by Figures 6 and 7, respectively.In Figures 6a and 7a, several closed loops represent the orthogonal projection of the contour lines on the surface z = f (x, y) onto the horizontal plane x − y, respectively.In Figure 7b,e, several closed loops also represent orthogonal projection of the contour lines on the surface z = F(x, y) onto the horizontal plane x − y, respectively.In Figure 6a, the vector starting with point x 0 is gradient −∇ f (x 0 ), and the length of the vector is f (x 0 )/ ∇ f (x 0 ), ∇ f (x 0 ) .For arbitrary initial iterative point x 0 , the iterative formula Step 2 of Algorithm 2) from the steepest descent method repeatedly iterates until the iterative point x n lies on the planar implicit curve f (x).
In Figure 6b, the footpoint q, i.e., the intersection of tangent line (from the point x n on the planar implicit curve f (x)) and perpendicular line (from test point p) is acquired by Step 3 in Algorithm 2. After the first round iteration of Algorithm 2, replace the initial iterative point x 0 with the footpoint q, and then, do the second round and the third round iteration.The three rounds of iteration constitute Algorithm 2 and part of Algorithm 3.
In each sub-figure of Figure 7, points p and p Γ are the test point and the corresponding orthogonal projective point, respectively.In Figure 7a, the vector starting with point x n is gradient −∇ f (x n ), and the length of the vector is f (x n )/ ∇ f (x n ), ∇ f (x n ) .For the initial iterative point x n from Algorithm 2, the iterative formula y of Algorithm 1) from the steepest descent method iterates once.In Figure 7b, the vector starting with point y n is gradient −∇F(y n ), and the length of the vector is F(y n )/ ∇F(y n ), ∇F(y n ) .

For the initial iterative point y
, the iterative Step 2 of Algorithm 1) from the steepest descent method iterates once.In Figure 7c, the footpoint q, i.e., the intersection of tangent line (from the point z' n on the planar implicit curve f (x)) and the perpendicular line (from test point p), is acquired by Step 3 in Algorithm 1.In the actual iterative process, point z' n is approximately equivalent to the point z n .In Figure 7d, point u n comes form the fourth step of Algorithm 1, which aims to obtain a linear orthogonal increment.In Figure 7e, the vector starting with point u n is gradient −∇F(u n ), and the length of the vector is Step 5 of Algorithm 1) from the steepest descent method iterates once more.In Figure 7f, the iterative point x n+1 from the sixth step in Algorithm 1 gives a correction for the iterative point v n from the fifth step in Algorithm 1. Repeatedly iterate the above six steps until the iteration exit criteria are met.In the end, three points of footpoint q, the iterative point x n+1 and the orthogonal projecting point p Γ merge into orthogonal projection point p Γ .These six steps constitute Algorithm 1 and part of Algorithm 3.   21); (c) The third step of the iterative Formula ( 21); (d) The fourth step of the iterative Formula ( 21); (e) The fifth step of the iterative Formula ( 21); (f) The sixth step of the iterative Formula (21).

Convergence Analysis
In this section, the convergence analysis for the integrated hybrid second order algorithm is presented.Proofs indicate the convergence order of the algorithm is up to two, and Algorithm 3 is independent of the initial value.
Theorem 1.Given an implicit function f (x) that can be parameterized, the convergence order of the iterative Formula ( 21) is up to two.
Proof.Without loss of generality, assume that the parametric representation of the planar implicit curve Γ : The first part will derive that the order of convergence of the first step for the iterative Formula ( 21) is up to two.It is not difficult to know the iteration equation in the corresponding Newton's second order parameterized iterative method, i.e., the first step for the iterative Formula (21): Taylor expansion around α generates: where e n = t n − α and Thus, it is easy to have: From ( 22)-( 24), the error iteration can be expressed as, where The second part will prove that the order of convergence of the second step for the iterative Formula ( 21) is two.It is easy to get the corresponding parameterized iterative equation for Newton's second-order iterative method, essentially the second step for the iterative Formula (21), where: Using Taylor expansion around α, it is easy to get: where e n = t n − α and b i = (1/i!)(F(α)), i = 0, 1, 2. Thus, it is easy to get: According to Formula ( 26)-( 29), after Taylor expansion and simplifying, the error relationship can be expressed as follows, where Because the fifth step is completely equal to the second step of the iterative Formula (21) and outputs from Newton's iterative method are closely related with test point p, the order of convergence for the fifth step of the iterative Formula ( 21) is also two.The third part will derive that the order of convergence of the third step and fourth step for iterative Formula ( 21) is one.According to the first order method for orthogonal projection onto the parametric curve [32,39,40], the footpoint q = (q 1 , q 2 ) of the parameterized iterative equation of the third step of the iterative Formula ( 21) can be expressed in the following way, q = c(t n ) + ∆tc (t n ). ( From the iterative Equation ( 31) and combining with the fourth step of the iterative Formula ( 21), it is easy to have: where x, y denotes the scalar product of vectors x, y ∈ R 2 .Let t n + ∆t → t n , and repeat the procedure (32) until ∆t is less than a given tolerance ε.Because parameter α is the orthogonal projection point of test point p = (p 1 , p 2 ) onto the parametric curve Because the footpoint q is the intersection of the tangent line of the parametric curve c(t) at t = t n and the perpendicular line − → pq determined by the test point p, the equation of the tangent line of the parametric curve c(t) at t = t n is: At the same time, the vector of the line segment connected by the test point p and the point c(t n ) is: ( The vector (35) and the tangent vector c (t n ) = ( f 1 (t n ), f 2 (t n )) of the tangent line (34) are mutually orthogonal, so the parameter value s 0 of the tangent line (34) is: Substituting ( 36) into (34) and simplifying, it is not difficult to get the footpoint q = (q 1 , q 2 ), Substituting ( 37) into (32) and simplifying, it is easy to obtain, From (33) and combined with (38), using Taylor expansion by the symbolic computation software Maple 18, it is easy to get: Simplifying (30), it is easy to obtain: where the symbol C 2 denotes the coefficient in the first order error e n of the right-hand side of Formula (40).The result shows that the third step and the fourth step of the iterative Formula ( 21) comprise the first order convergence.According to the iterative Formula (21) and combined with three error iteration relationships ( 25), ( 30) and ( 40), the convergent order of each iterative formula is not more than two.Then, the iterative error relationship of the iterative Formula ( 21) can be expressed as follows: To sum up, the convergence order of the iterative Formula ( 21) is up to two.
Theorem 2. The convergence of the hybrid second order algorithm (Algorithm 1) is a compromise method between the local and global method.
Proof.The third step and fourth step of the iterative Formula ( 21) of Algorithm 1 are equivalent to the foot point algorithm for implicit curves in [32].The work in [14] has explained that the convergence of the foot point algorithm for the implicit curve proposed in [14] is a compromise method between the local and global method.Then, the convergence of Algorithm 1 is also a compromise method between the local and global method.Namely, if a test point is close to the foot point of the planar implicit curve, the convergence of Algorithm 1 is independent of the initial iterative value, and if not, the convergence of Algorithm 1 is dependent on the initial iterative value.The sixth step in Algorithm 1 promotes the robustness.However, the third step, the fourth step and the sixth step in Algorithm 1 still constitute a compromise method between the local and global ones.Certainly, the first step (steepest descent method) of Algorithm 1 can make the iterative point fall on the planar implicit curve and improves its robustness.The second step and the fifth step constitute the classical Newton's iterative method to accelerate convergence and improve robustness in some way.The steepest descent method of the first step and Newton's iterative method of the second step and the fifth step in Algorithm 1 are more robust and efficient, but they can change the fact that Algorithm 1 is the compromise method between the local and global ones.To sum up, Algorithm 1 is the compromise method between the local and global ones.
Theorem 3. The convergence of the integrated hybrid second order algorithm (Algorithm 3) is independent of the initial iterative value.
Proof.The integrated hybrid second order algorithm (Algorithm 3) is composed of two parts sub-algorithms (Algorithm 1 and Algorithm 2).From Theorem 2, Algorithm 1 is a compromise method between the local and global method.Of course, whether the test point p is very far away or not far away from the planar implicit curve f (x), if the initial iterative value lies close to the orthogonal projection point p Γ , Algorithm 1 could be convergent.In any case, Algorithm 2 can change the initial iterative value of Algorithm 1 sufficiently close to the orthogonal projection point p Γ to ensure the convergence of Algorithm 1.In this way, Algorithm 3 can converge for any initial iterative value.Therefore, the convergence of the integrated hybrid second order algorithm (Algorithm 3) is independent of the initial value.
In Examples 2-6, if the degree of every planar implicit curve is more than five, it is difficult to get the intersection between the line segment determined by test point p and p + and the planar implicit curve by using the circle shrinking algorithm in [14].The running time comparison for Algorithm in [14] was not done, and it was not done for the circle double-and-bisect algorithm in [36] due to the same reason.The running time comparison test by using the circle double-and-bisect algorithm in [36] has not been done because it is difficult to solve the intersection between the circle and the planar implicit curve by using the circle double-and-bisect algorithm.In addition, many methods (Newton's method, the geometrically-motivated method [31,32], the osculating circle algorithm [33], the Bézer clipping method [25][26][27], etc.) cannot guarantee complete convergence for Examples 2-5.The running time comparison test for those methods in [25][26][27][31][32][33] has not been done yet.From Table 2 in [36], the circle shrinking algorithm in [14] is faster than the existing methods, while Algorithm 3 is faster than the circle shrinking algorithm in [14] in our Example 1.Then, Algorithm 3 is faster than the existing methods.Furthermore, Algorithm 3 is more robust and efficient than the existing methods.Besides, it is not difficult to find that if test point p is close to the planar implicit curve and initial iterative point x 0 is close to the test point p, for a lower degree of and fewer terms in the planar implicit curve and lower precision of the iteration, Algorithm 3 will use less total average running time.Otherwise, Algorithm 3 will use more time.
To find as many solutions as possible, a larger value of m is taken.

Remark 6.
In Example 1, for the test points (−0.1,1.0),(0.2,1.0), (0.1,0.1), (0.45,0.5), by using Algorithm 3, the corresponding orthogonal projection points p Γ are (−0.47144354751227009,0.70879213227958752), (−0.42011639143389254, 0.63408011508207950), (−0.33334322619432892, 0.099785192603767206), (−0.34352305539212918,0.401230229163152532), respectively (see Figure 14 and Table 16).In addition to the six test examples, many other examples have also been tested.According to these results, if test point p is close to the planar implicit curve f (x), for different initial iterative values x 0 , which are also close to the corresponding orthogonal projection point p Γ , it can converge to the corresponding orthogonal projection point p Γ by using Algorithm 3, namely the test point p and its corresponding orthogonal projection point p Γ satisfy the inequality relationships: Thus, it illustrates that the convergence of Algorithm 3 is independent of the initial value and Algorithm 3 is efficient.In sum, the algorithm can meet the top two of the ten challenges proposed by Professor Les A. Piegl [41] in terms of robustness and efficiency.Remark 7. From the authors' six test examples, Algorithm 3 is robust and efficient.If test point p is very far away from the planar implicit curve and the degree of the planar implicit curve is very high, Algorithm 3 also converges.However, inequality relationships (42) could not be satisfied simultaneously.In addition, if the planar implicit curve contains singular points, Algorithm 3 only works for test point p in a suitable position.Namely, for any initial iterative point x 0 , test point p can be orthogonally projected onto the planar implicit curve, but with a larger distance p − p Γ than the minimum distance p − p s between the test point and the orthogonal projection point, where p s is the singular point.For example, for the test point (1.0,0.01),(0.6,0.1), (0.5,−0.15), (0.8,−0.1),Algorithm 3 gives the corresponding orthogonal projection points p Γ as (0.66370473801453017, 0.092784537693334545), (0.66704812931370775, 0.097528910436113817), (0.663704738014530, 0.13435089298485379), (0.66418591136724639, −0.090702201378858334), respectively.However, the actual corresponding orthogonal projection point of four test points is (0.66666666666666667, 0.0) (see Figure 14 and Table 16).Iterations 10 6.69 × 10 −1 15 6.69 × 10 −1 18 6.69 × 10 −1 22 6.69 × 10 −1 32 6.69 × 10 −1 14 6.69 × 10 −1

Conclusions
This paper investigates the problem related to a point projection onto a planar implicit curve.The integrated hybrid second order algorithm is proposed, which is composed of two sub-algorithms (hybrid second order algorithm and initial iterative value estimation algorithm).For any test point p, any planar implicit curve containing singular points, whether test point p is close to or very far away from the planar implicit curve, the integrated hybrid second order algorithm could be convergent.It is proven that the convergence of Algorithm 3 is independent of the initial value.Convergence analysis of the integrated hybrid second order algorithm demonstrates that the convergence order is second order.Numerical examples illustrate that the algorithm is robust and efficient.

Future Work
For any initial iterative point and test point in any position of the plane, for any planar implicit curve (including containing singular points, the degree of the planar implicit curve being arbitrarily high), the future work is to construct a brand new algorithm to meet three requirements: (1) it does converge, and the orthogonal projection point does simultaneously satisfy three relationships of Formula (11); (2) it is very effective at tackling singularity; (3) it takes less time than the current Algorithm 3. Of course, it will be very challenging to find this kind of algorithm in the future.
Another potential topic for future research is to develop a more efficient method to compute the minimum distance between a point and a spatial implicit curve or a spatial implicit surface.

Figure 2 .Figure 3 .
Figure 2. The entire graphical demonstration of the first round iteration in Algorithm 2. (a) Initial status; (b) Intermediate status; (c) Final status.

Figure 4 .
Figure 4.The entire graphical demonstration of the third round iteration in Algorithm 2. (a) Initial status; (b) Intermediate status; (c) Final status.

Figure 5 .
Figure 5.The entire graphical demonstration for the whole iterative process of Algorithm 3. (a) Initial status; (b) First intermediate status; (c) Second intermediate status; (d) Third intermediate status; (e) Fourth intermediate status; (f) Final status.

Figure 7 .
Figure 7.The entire graphical demonstration of Algorithm 1.(a) The first step of the iterative Formula (21); (b) The second step of the iterative Formula (21); (c) The third step of the iterative Formula (21); (d) The fourth step of the iterative Formula (21); (e) The fifth step of the iterative Formula (21); (f) The sixth step of the iterative Formula (21).

Table 1 .
Convergence of the hybrid second order algorithm for four given test points.

Table 2 .
Oscillation of the hybrid second order algorithm for the planar implicit curve with a far-away test point.

Table 3 .
Running time for different initial iterative values by Algorithm 3 in Example 1.

Table 5 .
Running times for different initial iterative values by Algorithm 3 in Example 1.

Table 6 .
Running times for different initial iterative values by Algorithm 3 in Example 1.
, the average running times of Algorithm 3 for eight different initial iterative values are 4,487,449, 4,202,203, 4,555,396, 4,533,326, 4,304,781, 4,163,107, 4,268,792 and 4,378,470 nanoseconds, respectively.In the end, the overall average running time is 4,361,691 nanoseconds (see Figure

Table 7 .
Running times for different initial iterative values by Algorithm 3 in Example 2.

Table 9 .
Running times for different initial iterative values by Algorithm 3 in Example 3.

Table 10 .
The error analysis of the iteration process of Algorithm 3 in Example 3.

Table 11 .
Running times for different initial iterative values by Algorithm 3 in Example 4.

Table 13 .
Running times for different initial iterative values by Algorithm 3 in Example 5.

Table 14 .
The error analysis of the iteration process of Algorithm 3 in Example 5.

Table 15 .
The error analysis of the iteration process of Algorithm 3 in Example 6.

Table 21 .
The error ratios for each iteration in Example 5 of Algorithm 3.

Table 22 .
The error ratios for each iteration in Example 6 of Algorithm 3.