An Improved Curvature Circle Algorithm for Orthogonal Projection onto a Planar Algebraic Curve

Abstract

Point orthogonal projection onto planar algebraic curve plays an important role in computer graphics, computer aided design, computer aided geometric design and other fields. For the case where the test point $\mathbf{p}$ is very far from the planar algebraic curve, we propose an improved curvature circle algorithm to find the footpoint. Concretely, the first step is to repeatedly iterate algorithm (the Newton’s steepest gradient descent method) until the iterated point could fall on the planar algebraic curve. Then seek footpoint by using the algorithm (computing footpoint $\mathbf{q}$) where the core technology is the curvature circle method. And the next step is to orthogonally project the footpoint $\mathbf{q}$ onto the planar algebraic curve by using the algorithm (the hybrid tangent vertical foot algorithm). Repeatedly run the algorithm (computing footpoint $\mathbf{q}$) and the algorithm (the hybrid tangent vertical foot algorithm) until the distance between the current footpoint and the previous footpoint is near 0. Furthermore, we propose Second Remedial Algorithm based on Comprehensive Algorithm B. In particular, its robustness is greatly improved than that of Comprehensive Algorithm B and it achieves our expected result. Numerical examples demonstrate that Second Remedial Algorithm could converge accurately and efficiently no matter how far the test point is from the plane algebraic curve and where the initial iteration point is.

1. Introduction

Reconstructing curve/surface is an important work in the field of computer aided geometric design, especially in geometric modeling and processing where it is crucial to fit curve/surface in high accuracy and reduce the error of representation curve/surface. The representation of the four curve types are the explicit-type, implicit-type, parametric-type and subdivision-type. Because implicit representation has unique advantage in the process of computer aided geometric design, it has wide and far-reaching applications. From scattered and unorganized three-dimensional data, Bajaj et al. [1] reconstructed surface and functions on surfaces. They [2,3] have constructed the algebraic B-spline surfaces with least-squares fitting feature using tensor product technique. Schulz et al. [4] constructed an enveloping algebraic surface using gradually approximate algebraization method. Kanatani et al. [5] applied the algebraic curve to construct geometric ellipse fitting using unified strict maximum likelihood estimation method. Mullen et al. [6] reconstructed robust and accurate algebraic surface functions to sign the unsigned from scattered and unorganized three-dimensional data point sets. Upreti et al. [7] used a technique to sign algebraic level sets on NURBS surface and algebraic Boolean level sets on NURBS surfaces. Rouhani et al. [8] applied the algebraic function for polynomial representation system. And L.G. Zagorchev et al. [9] applied the algebraic function for general algebraic surface.

Up to now, there are three main types of methods to solve the problem of point orthogonal projection onto planar algebraic curve: local method, global method and compromise method between these two methods. Here are three typical approaches.

According to the most basic geometric characteristic, orthogonal projection of test point $\mathbf{p}$ onto the planar algebraic curve is actually the point $\mathbf{x}$ on the curve such that cross product of vectors $\overrightarrow{\mathbf{xp}}$ and $\nabla f\left(\mathbf{x}\right)$ is 0.

$$\nabla f\left(\mathbf{x}\right)\times (\mathbf{p}-\mathbf{x})=0.$$

(1)

Equation (1) can be transformed into Newton’s iteration formula (3). Furthermore, Sullivan et al. [10] adopted a hybrid method with Lagrange multiplier and Newton’s iterative method to compute the closest point on the planar algebraic curve for each test point. Some orthogonal projection problems can be transformed into solving system of nonlinear equations. The common characteristic of methods [10,11] is that they converge locally and fast, while methods [10,11] are dependent on the initial points.

The first global method of solving system of nonlinear equations is the Homotopy continuous method [12,13]. They constructed Homotopy continuous formula.

$$H(\mathbf{x},t)=(1-t)P\left(\mathbf{x}\right)+tQ\left(\mathbf{x}\right),t\in \left[0,1\right]$$

(2)

where t is a parameter of continuous transformation from 0 to 1, $P\left(\mathbf{x}\right)=0$ is the original system of nonlinear equations to be solved, $Q\left(\mathbf{x}\right)$ is the objective solution of system of nonlinear equations $P\left(\mathbf{x}\right)=0$. All isolated solutions of system of nonlinear equations $P\left(\mathbf{x}\right)=0$ can be computed by the numerical continuous Homotopy methods [12,13]. So the Homotopy methods [12,13] are global convergence. The Homotopy methods’ robustness is proved by [14] and their high time-consuming property is verified in [15]. Of course, the Homotopy methods [12,13] are ideal in theory, but it is difficult to find or construct the objective system of nonlinear equations $Q\left(\mathbf{x}\right)=0$ in practical engineering applications.

The second global resultant methods convert system of nonlinear equations into the expression of the resultants and then solve the resultants [16,17,18,19]. According to classical elimination theory, system of two nonlinear equations with two variables can be turned into a resultant polynomial with one variable, which is equivalent to the two simultaneous equations. The Sylvester’s resultant and Cayley’s statement of Bézout’s method are the most famous resultant methods [16,17,18,19]. Because the resultant methods [16,17,18,19] can solve all roots if the degree of the planar algebraic curve is less than 4, they are good global methods. However, if the degree of the planar algebraic curve is more than quintic, it becomes harder and harder with increasing degree to solve two-polynomial system with the resultant methods.

The third global method is the adoption of the Bézier clipping technique [20,21,22]. In the first step, solving the nonlinear system of Equation (1) is transformed into solving all roots of Bernstein-Bézier representation with convex hull property. In the second step, if the parts of the domains do not include the solution, we clip the parts of the domains by using convex hull box with Bernstein-Bézier form such that the discarded parts of the region has no solution and all the solutions are in the retained parts of the region. In the third step, the de Casteljau subdivision rule is used to segment the remaining part of the curve obtained by elimination in step 2. Repeat steps 2 and 3 until we can find all the solutions to Equation (1). The advantage of this method is that all solutions of Equation (1) can be found. But this global clipping method has one difficulty: sometimes Equation (1) is difficult or even impossible to convert into Bernstein-Bézier form. For example, specific Equation (1) $\nabla f\left(\mathbf{x}\right)\times (\mathbf{p}-\mathbf{x})=-36p_1{y}^{17}+36x{y}^{17}+6p_2{x}^5-6x^{5}y-4p_{1}x+4p_{2}y+4x^2-4y^2$ is impossible to convert into Bernstein-Bézier form where $f\left(\mathbf{x}\right)=x^6+4xy+2y^{18}-1=0$.

The compromise method is between local and global methods. Consisting of the geometric property with computing the nearest point is proposed by Hartmann [23,24] named as the first compromise method. Repeatedly run the Newton’s steepest gradient descent method (3) until the iterative point falls on the planar algebraic curve, where the initial iterative point is unrestricted.

$${\mathbf{x}}_{n+1}={\mathbf{x}}_n-(f\left({\mathbf{x}}_n\right)/〈\nabla f\left({\mathbf{x}}_n\right),\nabla f\left({\mathbf{x}}_n\right)〉)\nabla f\left({\mathbf{x}}_n\right).$$

(3)

$$\mathbf{q}=\mathbf{p}-(〈\mathbf{p}-{\mathbf{y}}_n,\nabla f\left({\mathbf{y}}_n\right)〉/〈\nabla f\left({\mathbf{y}}_n\right),\nabla f\left({\mathbf{y}}_n\right)〉)\nabla f\left({\mathbf{y}}_n\right).$$

(4)

Running the iterative formula (4) one time, the method [23,24] can obtain the vertical foot point $\mathbf{q}$ where the iterative point ${\mathbf{y}}_n$ of the formula (4) is the final iterative point obtained by formula (3). Continuously iterate the above two steps until the vertical foot point $\mathbf{q}$ is on the planar algebraic curve $f\left(\mathbf{x}\right)$. Unluckily, progressive geometric tangent approximation iteration method with computing vertical foot point $\mathbf{q}$ fails for some planar algebraic curves $f\left(\mathbf{x}\right)$.

The second compromise method is developed by Nicholas [25] who adopted the osculating circle technique to realize orthogonal projection onto the planar algebraic curve. Calculate the corresponding curvature at one point on the planar algebraic curve, and then the radius and center of the curvature circle. The line segment formed by the test point and the center of the curvature circle intersects the curvature circle at footpoint $\mathbf{q}$. Approximately take the footpoint $\mathbf{q}$ as a point on the planar algebraic curve. For the new point on the planar algebraic curve, repeat the above procedure to get a new footpoint and corresponding new approximate point on the planar algebraic curve. Repeat the above behavior until the footpoint $\mathbf{q}$ is the orthogonal projection point ${\mathbf{p}}_{\Gamma }$. Because the planar algebraic curve does not have parametric control like parametric curve, taking the footpoint as an approximate point on the planar algebraic curve will bring about large errors. So it makes the operation of the whole algorithm unstable.

The third compromise method is the circle shrinking technique [26]. Repeatedly run the iterative formula (3) such that the final iterative point ${\mathbf{p}}_c$ falls on the planar algebraic curve as far as possible, where the selection of initial iterative point is arbitrary. The next iterative point on the planar algebraic curve is obtained through a series of combined operations of circle and the planar algebraic curve, where the center and radius of the circle are test point $\mathbf{p}$ and $∥\mathbf{p}-{\mathbf{p}}_c∥$, respectively. A series of combined operations include the two most important steps: Find a point ${\mathbf{p}}^{+}$ on the circle by means of the mean value theorem; Seek the intersection of the line segment $\overline{{\mathbf{p}\mathbf{p}}^{+}}$ and the circle where we call this intersection as the current intersection point ${\mathbf{p}}_c$. Repeatedly run this series of combined operations until the distance between the current point ${\mathbf{p}}_c$ and the previous point ${\mathbf{p}}_c$ is 0. The circle shrinking technique [26] takes a lot of time to seek point ${\mathbf{p}}^{+}$ each time. The algorithm has one difficulty: if the degree of the planar algebraic curve is higher than 5, the intersection point ${\mathbf{p}}_c$ of line segment $\overline{{\mathbf{p}\mathbf{p}}^{+}}$ and the planar algebraic curve cannot be solved directly by formula or the iterative methods to find the intersection ${\mathbf{p}}_c$ will lead to instability.

The four compromise method is a circle double-and-bisect algorithm [27]. The circle doubling algorithm begins with a very small circle where the center is the test point $\mathbf{p}$ and the radius is very small $r_1$. Keep the same center of the circle, take the radius $r_2$ twice of $r_1$ to draw a new circle. If there is no intersection between the new circle and the curve, draw a new circle with radius twice of $r_2$. Continuously repeat the above process until new circle can intersect with the planar algebraic curve and the former circle does not. Naturally, the former circle and the current circle are called interior circle and exterior circle, respectively. Moreover, the bisecting technology implements the rest of the process. Continue to draw a new circle with new radius $r=(r_1+r_2)/2$. If the new current circle whose radius is r intersects with the curve, substitute r for $r_2$, else for $r_1$. Repeatedly run the above progress until the difference between the two radii is approximate zero($\left|r_1-r_2\right|<\epsilon )$. But this method is very difficult to judge whether the exterior circle intersects the planar algebraic curve or not [27].

The fifth compromise method is the integrated hybrid second order algorithm [28]. It includes two sub-algorithms: the hybrid second order algorithm and the initial iterative value estimation algorithm. They mainly exploint three ideas: (1) the tangent orthogonal vertical foot method coupled with calibration method; (2) Newton’s steepest gradient descent iterative method to impel the iteration point to be on the planar implicit curve; (3) Newton’s iterative method to speed up the whole iteration process. Before running the hybrid second order algorithm, the initial iterative value estimation algorithm is used to force the initial iterative value of the formula (17) of the hybrid second order algorithm and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ as close as possible. After a lot of tests, if the distance between the test point $\mathbf{p}$ and the curve is not very far, the advantages of this algorithm are obvious in term of robustness and efficiency. But when the test point is very far from the curve, the integrated hybrid second order algorithm is invalid.

2. The Improved Curvature Circle Algorithm

In Reference [28], when the test point $\mathbf{p}$ is not particularly far away, the integrated hybrid algorithm can have ideal result. But if the test point $\mathbf{p}$ is very far from the curve, the algorithm is invalid where the test point $\mathbf{p}$ can not be robustly and effectively orthogonally projected onto the planar algebraic curve. In order to overcome this difficulty, we propose an improved curvature circle algorithm to ensure robustness and effective convergence with the test point $\mathbf{p}$ being arbitrarily far away. No matter how far the test point $\mathbf{p}$ is from the planar algebraic curve, if the initial iteration point ${\mathbf{x}}_0$ is very close to the orthogonal projection point of the test point $\mathbf{p}$, the preconceived algorithm can converge well. So we attempt to construct an algorithm to find an initial iterative point very close to the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the test point $\mathbf{p}$. The general idea is the following. Repeatedly iterate the formula (3) by utilizing the Newton’s steepest gradient descent method until the iteration point fall on the planar algebraic curve as far as possible, written as ${\mathbf{p}}_c$. This time, the distance between the iteration point ${\mathbf{p}}_c$ and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ is much smaller than that between the original iteration point ${\mathbf{x}}_0$ and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$. The iteration point ${\mathbf{p}}_c$ is closer to the orthogonal projection point ${\mathbf{p}}_{\Gamma }$. In order to further promote the iteration point ${\mathbf{p}}_c$ and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ to be closer, we introduce a key step with curvature circle algorithm. Draw a curvature circle through point ${\mathbf{p}}_c$ on the planar algebraic curve with the radius R determined by the curvature k and the center $\mathbf{m}$ being a normal direction point of point ${\mathbf{p}}_c$ on the planar algebraic curve. Line segment $\overline{\mathbf{m}\mathbf{p}}$ determined by the test point $\mathbf{p}$ and the center $\mathbf{m}$ intersects curvature circle at point $\mathbf{q}$. We take the intersection point $\mathbf{q}$ as the next iteration point for the iteration point ${\mathbf{p}}_c$. Of course, the distance between the intersection point $\mathbf{q}$ and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ is much smaller than the previous one. We use the intersection point $\mathbf{q}$ as the new test point, and run the hybrid algorithm again where the initial iterative point at this moment can be set as $\mathbf{q}-(0.1,0.1)$. Repeatedly iterate until the iteration point falls on the planar algebraic curve $f\left(\mathbf{x}\right)$, written as ${\mathbf{p}}_c$. We repeat the last two key steps in this procedure until the iteration point ${\mathbf{p}}_c$ and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ overlap (See Figure 1).

Let’s elaborate on the general idea. Let $\mathbf{p}$ be a test point on the plane. There is an planar algebraic curve $\Gamma$ on the plane.

$$f(x,y)=0.$$

(5)

The plane algebraic curve (5) can be simply written as

$$f\left(\mathbf{x}\right)=0,$$

(6)

where $\mathbf{x}=(x,y)$. The goal of this paper is to find a point ${\mathbf{p}}_{\Gamma }$ on the planar algebraic curve $f\left(\mathbf{x}\right)$ to satisfy the basic relationship

$$∥\mathbf{p}-{\mathbf{p}}_{\Gamma }∥=\underset{\mathbf{x}\in \Gamma }{min}∥\mathbf{p}-\mathbf{x}∥.$$

(7)

The above problem can be written as

$$\left\{\begin{array}{c}f\left({\mathbf{p}}_{\Gamma }\right)=0,\hfill \\ \nabla f\left({\mathbf{p}}_{\Gamma }\right)\times (\mathbf{p}-{\mathbf{p}}_{\Gamma })=0,\hfill \\ ∥\mathbf{p}-{\mathbf{p}}_{\Gamma }∥=\underset{\mathbf{x}\in \Gamma }{min}∥\mathbf{p}-\mathbf{x}∥,\hfill \end{array}\right.$$

(8)

where ${\displaystyle \nabla f=\left[\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}\right]}$ is Hamiltonian operator and symbol × is cross product. We take s as the arc length parameter of the planar algebraic curve $f\left(\mathbf{x}\right)$ and ${\displaystyle \mathbf{t}=\left[\frac{dx}{ds},\frac{dy}{ds}\right]}$ is unit tangent vector along the planar algebraic curve $f\left(\mathbf{x}\right)$. Take derivative of Equation (6) with respect to arc length parameter s and combine with unit tangent vector condition $∥\mathbf{t}∥=1$, we obtain the following simultaneous system of nonlinear equations,

$$\left\{\begin{array}{c}〈\mathbf{t},\nabla f〉=0,\hfill \\ ∥\mathbf{t}∥=1.\hfill \end{array}\right.$$

(9)

It is easy to get the solution of Equation (9).

$$\mathbf{t}=\left[-\frac{\partial f}{\partial y},\frac{\partial f}{\partial x}\right]/∥\nabla f∥.$$

(10)

Repeatedly iterate Equation (3) called as the Newton’s steepest gradient descent method until until the iterative termination criteria $\left|f\left({\mathbf{x}}_{n+1}\right)\right|<\epsilon$, where the initial iterative point is ${\mathbf{x}}_0=\mathbf{p}-(0.1,0.1)$ and refer to the iterative point ${\mathbf{x}}_{n+1}$ as ${\mathbf{p}}_c$. The first advantage of the Newton’s steepest gradient descent method (3) is to make the iteration point fall on the planar algebraic curve $f\left(\mathbf{x}\right)$ as far as possible. Its second advantage of the Newton’s steepest gradient descent method (3) is that the iteration point fallen on the planar algebraic curve is relatively close to the orthogonal projection point ${\mathbf{p}}_{\Gamma }$, and it brings great convenience to implementation of the subsequent sub-algorithms. The Newton’s steepest gradient descent method (Algorithm 1) can be specifically described as (See Figure 2).

Algorithm 1: The Newton’s steepest gradient descent method.

Input: The test point p and the planar algebraic curve $f\left(\mathbf{x}\right)=0$     Output: The iterative point ${\mathbf{p}}_c$ fallen on planar algebraic curve $f\left(\mathbf{x}\right)=0$     Description:     Step 1:                 ${\mathbf{x}}_{n+1}=\mathbf{p}-(0.1,0.1)$;                 Do {                         ${\mathbf{x}}_n={\mathbf{x}}_{n+1}$ ;                         Update ${\mathbf{x}}_{n+1}$ according to the iterative Equation (3);                 }while ($\left|f\left({\mathbf{x}}_{n+1}\right)\right|>\epsilon \&\&∥{\mathbf{x}}_{n+1}-{\mathbf{x}}_n∥>\epsilon )$;     Step 2:                 ${\mathbf{p}}_c={\mathbf{x}}_{n+1}$ ;                 Return ${\mathbf{p}}_c$;

In this case, if the iterative point ${\mathbf{p}}_c$ fallen on the planar algebraic curve $f\left(\mathbf{x}\right)$ is taken as the initial iterative point of the hybrid algorithm, convergence or divergence may occur where divergence can not improve the algorithm. As for divergence, it can not achieve the purpose of improving the algorithm. From another point of view, the distance between iteration point ${\mathbf{p}}_c$ and orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the test point $\mathbf{p}$ should be closer. It lays a good foundation for the implementation of subsequent sub-algorithms. In order to get the iteration point and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ closer, we adopt curvature circle way to promote the iteration point and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ being closer. Because the iterative point is on the planar algebraic curve, the curvature k at the iterative point ${\mathbf{p}}_c$ fallen on the planar algebraic curve $f\left(\mathbf{x}\right)$ is defined as [29],

$$k=k(x,y)=\frac{\left[-f_y,f_x\right]\mathbf{G}{\left[-f_y,f_x\right]}^T}{{∥\nabla f∥}^3},$$

(11)

where $\mathbf{G}=\left(\begin{array}{cc}{f}_{xx}\hfill & f_{xy}\hfill \\ f_{yx}\hfill & f_{yy}\hfill \end{array}\right)$. The radius R and the center $\mathbf{m}$ of the curvature circle ◯ directed by the curvature k are

$$R=\left|1/k\right|,$$

(12)

and

$$\mathbf{m}={\mathbf{p}}_c+\frac{\overrightarrow{\mathbf{n}}}{k},$$

(13)

where the unit normal vector ${\displaystyle \overrightarrow{\mathbf{n}}}$ is ${\displaystyle \overrightarrow{\mathbf{n}}=\frac{\nabla f}{∥\nabla f∥}}$. The line segment $\overline{\mathbf{mp}}$ determined by the test point $\mathbf{p}$ and the center $\mathbf{m}$ of the curvature circle ◯ intersects the curvature circle ◯ at point $\mathbf{q}$ which is named as footpoint $\mathbf{q}$. From elementary geometric knowledge, the parametric equation of the line segment $\overline{\mathbf{mp}}$ can be expressed as

$$\mathbf{x}=\mathbf{p}+(\mathbf{m}-\mathbf{p})w,$$

(14)

where parametric $0\le w\le 1$ is undetermined. In addition, the equation of the curvature circle ◯ can be written as

$$∥\mathbf{m}-\mathbf{x}∥=R.$$

(15)

By solving Equation (14) and Equation (15) together, the analytic expression of the intersection $\mathbf{q}$ is obtained

$$\mathbf{q}=\mathbf{p}+(\mathbf{m}-\mathbf{p})w,$$

(16)

where the undetermined parameter w is accurately identified as ${\displaystyle w=1-\frac{R}{∥\mathbf{m}-\mathbf{p}∥}}$. The computation of the footpoint $\mathbf{q}$ can be realized through Algorithm 2 (See Figure 3).

Algorithm 2: Computing footpoint $\mathbf{q}$ via the curvature circle ◯ and the line segment $\overline{\mathbf{m}\mathbf{p}}$.

Input: The test point $\mathbf{p}$, the planar algebraic curve $f\left(\mathbf{x}\right)=0$ and the iterative point ${\mathbf{p}}_c$ on the     planar algebraic curve $f\left(\mathbf{x}\right)=0$.     Output: The footpoint $\mathbf{q}$.     Description:     Step 1:                 Compute the curvature k of the iterative point ${\mathbf{p}}_c$ fallen on the planar algebraic curve     $f\left(\mathbf{x}\right)=0$ by the curvature calculation formula (11).     Step 2:                 Calculate the radius R and the center $\mathbf{m}$ of the curvature circle ◯ through the formulas     (12) and (13), respectively.     Step 3:                 Compute the footpoint $\mathbf{q}$ by the formula (16).                 Return $\mathbf{q}$;

Remark 1.

The important formula for computing the curvature k is the formula (11). If the denominator of the curvature k with the formula (11) is 0, the whole iteration process will degenerate. In order to solve this special degeneration, we adopt a small perturbation of the curvature k of the formula (11) in programming implementation of Algorithm 2. Namely, the denominator of the curvature k with the formula (11) could be incremented by a small positive constant ε, the denominator of the curvature k is the denominator of the curvature k +ε, and Algorithm 2 continues to calculate the center and the radius of the curvature circle corresponding to the curvature after disturbance. Of course, in all subsequent formulas or iterative formulas, we also do the same denominators perturbation treatment for the case of the zero denominators of the formulas or the iterative formulas.

Under this circumstance, if the footpoint point $\mathbf{q}$ at this moment is taken as the initial iteration point of the hybrid algorithm, the convergence probability of the hybrid algorithm is much greater than that of using the point ${\mathbf{p}}_c$ in Algorithm 1 as the initial iterative point of the hybrid algorithm. The reason is that the distance $∥\mathbf{q}-{\mathbf{p}}_{\Gamma }∥$ is smaller than the distance $∥{\mathbf{p}}_c-{\mathbf{p}}_{\Gamma }∥$. But divergence may happen in this case. In order to further guarantee the robustness,we orthogonally project the footpoint $\mathbf{q}$ onto the planar algebraic curve $f\left(\mathbf{x}\right)$ by using the hybrid algorithm, instead of directly using the footpoint q as the initial iterative point. At this time we still call the orthogonal projection point of the footpoint $\mathbf{q}$ as the point ${\mathbf{p}}_c$ which is just fallen on the planar algebraic curve $f\left(\mathbf{x}\right)$. Because at this time the footpoint $\mathbf{q}$ is close to the planar algebraic curve $f\left(\mathbf{x}\right)$, the algorithm can ensure complete convergence. The distance between the iterative point ${\mathbf{p}}_c$ and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the test point $\mathbf{p}$ becomes smaller again. The core iterative formula (17) of the hybrid algorithm is as follows (See [28]).

$$\left\{\begin{array}{c}{\mathbf{y}}_n={\mathbf{x}}_n-\left(f\left({\mathbf{x}}_n\right)/〈\nabla f\left({\mathbf{x}}_n\right),\nabla f\left({\mathbf{x}}_n\right)〉\right)\nabla f\left({\mathbf{x}}_n\right),\hfill \\ {\mathbf{z}}_n={\mathbf{y}}_n-\left(\mathbf{F}\left({\mathbf{y}}_n\right)/〈\nabla \mathbf{F}\left({\mathbf{y}}_n\right),\nabla \mathbf{F}\left({\mathbf{y}}_n\right)〉\right)\nabla \mathbf{F}\left({\mathbf{y}}_n\right),\hfill \\ \mathbf{Q}=\mathbf{q}-(〈(\mathbf{q}-{\mathbf{z}}_n),\nabla f{(\mathbf{z}}_n)〉/〈\nabla f{(\mathbf{z}}_n),\nabla f{(\mathbf{z}}_n)〉)\nabla f{(\mathbf{z}}_n),\hfill \\ {\mathbf{u}}_n={\mathbf{z}}_n+sign\left(〈\mathbf{q}-{\mathbf{z}}_n,{\mathbf{t}}_0〉\right){\mathbf{t}}_0\Delta s,\hfill \\ {\mathbf{v}}_n={\mathbf{u}}_n-\left(\mathbf{F}\left({\mathbf{u}}_n\right)/〈\nabla \mathbf{F}\left({\mathbf{u}}_n\right),\nabla \mathbf{F}\left({\mathbf{u}}_n\right)〉\right)\nabla \mathbf{F}\left({\mathbf{u}}_n\right),\hfill \\ {\mathbf{x}}_{n+1}{=\mathbf{v}}_n+\left[-\Delta e,0\right]{\left[\nabla f^T,(\Delta {\mathbf{v}}_n{)}^T\right]}^{-1}(if\left|\left[\nabla f^T,{(\Delta {\mathbf{v}}_n)}^T\right]\right|=0,{\mathbf{x}}_{n+1}={\mathbf{v}}_n),\hfill \end{array}\right.$$

(17)

where ${\mathbf{F}}_0\left(\mathbf{x}\right)=\left[(\mathbf{q}-\mathbf{x})\times \nabla f\left(\mathbf{x}\right)\right]=0$, ${\displaystyle \mathbf{F}\left(\mathbf{x}\right)=\frac{{\mathbf{F}}_0\left(\mathbf{x}\right)}{{\displaystyle \sqrt{〈\nabla f\left(\mathbf{x}\right),\nabla f\left(\mathbf{x}\right)〉}}}}$, ${\displaystyle {\mathbf{t}}_0=\frac{\left[-f_y,f_x\right]}{∥\Delta f∥}}$, $\Delta s=∥\mathbf{Q}-{\mathbf{z}}_n∥$, $f\left({\mathbf{v}}_n\right)=\Delta e$, $\Delta {\mathbf{v}}_n=-\left({\mathbf{F}(\mathbf{u}}_n)/〈\nabla {\mathbf{F}(\mathbf{u}}_n),\nabla {\mathbf{F}(\mathbf{u}}_n)〉\right)\nabla {\mathbf{F}(\mathbf{u}}_n)$.

The iterative formula (17) mainly contains four techniques. The core technology is the tangent foot vertical method with the third step and the fourth step of the iterative formula (17). Draw a tangent line $\mathbf{L}$ from a point on a plane algebraic curve $f\left(\mathbf{x}\right)$. Through the footpoint $\mathbf{q}$ (The footpoint $\mathbf{q}$ at this time is as the test point of iterative formula (17)), make a vertical line of the tangent $\mathbf{L}$ and get its corresponding vertical foot point $\mathbf{Q}$, which is equivalent to the third step in the formula (17). From the fourth step of the iterative formula (17), we get the next iteration point of particular importance for the initial iteration point. When the next iteration point is not very close to the planar algebraic curve $f\left(\mathbf{x}\right)$, we adopt the second important technique with the iteration point correction method, equivalent to the sixth step of the iterative formula (17). The iteration point is to move to the plane algebra curve as close as possible such that the distance between the correction point of the iteration point and the plane algebra curve $f\left(\mathbf{x}\right)$ is as close as possible. These two techniques are pure geometric techniques. When the distance between the test point and the planar algebraic curve is very close, the effect of convergence is obvious. Of course, when the distance between the test point and the planar algebraic curve is relatively long, sometimes there will be non-convergence. In order to improve the robustness of convergence, we add the Newton’s steepest gradient descent method before the first technique with the third step and the fourth step of the iterative formula (17). Its first aim is to bring the initial iteration point closer to the planar algebraic curve $f\left(\mathbf{x}\right)$. Its second aim is to promote the accuracy of subsequent iterations. In order to accelerate the whole iteration process of the iterative formula (17), we once again incorporate the fourth technology of Newton’s iterative method which is closely related to the footpoint $\mathbf{q}$. This technique not only accelerates the convergence rate of the whole iteration process but also improves the iteration robustness. Furthermore, the accuracy of the whole iteration process can be improved by the fourth technique. So we add Newton’s iterative method after the first step with the second technique, and then add it again before the last step with the third technique. Based on the above explanation and illustration, we get the following the hybrid tangent vertical foot algorithm (Algorithm 3).

Algorithm 3: The hybrid tangent vertical foot algorithm (See Figure 4).

Input: The footpoint $\mathbf{q}$ and the planar algebraic curve $f\left(\mathbf{x}\right)=0$.     Output: The point ${\mathbf{p}}_c$ fallen on the planar algebraic curve $f\left(\mathbf{x}\right)=0$.     Description:     Step 1:                 ${\mathbf{x}}_{n+1}=\mathbf{q}-(0.1,0.1)$;                 Do {                         ${\mathbf{x}}_n={\mathbf{x}}_{n+1}$;                         Execute ${\mathbf{x}}_{n+1}$ according to the iterative Equation (17);                 }while (${∥{\mathbf{x}}_{n+1}-{\mathbf{x}}_n∥}^2>\epsilon$&&$\left|f\left({\mathbf{x}}_{n+1}\right)\right|>\epsilon )$     Step 2:                 ${\mathbf{p}}_c={\mathbf{x}}_{n+1}$;                 Return ${\mathbf{p}}_c$;

With the description of the above three algorithms, we propose a comprehensive and complete algorithm (Algorithm 4) closely related to Algorithm 2 (See Figure 4).

Algorithm 4: The first improved curvature circle algorithm (Comprehensive Algorithm A).

Input: Test point $\mathbf{p}$ and the planar algebraic curve $f\left(\mathbf{x}\right)$.     Output: Orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the test point $\mathbf{p}$ which orthogonally projects the test      point $\mathbf{p}$ onto a planar algebraic curve $f\left(\mathbf{x}\right)$.     Description:     Step 1:Starting from the neighbor point of the test point $\mathbf{p}$, calculate the point ${\mathbf{p}}_c$ fallen on the      $f\left(\mathbf{x}\right)$ via Algorithm 1.     Do{           Step 2: Compute the footpoint $\mathbf{q}$ via Algorithm 2.           Step 3: Project footpoint $\mathbf{q}$ onto the planar algebraic curve $f\left(\mathbf{x}\right)$ via Algorithm 3, then get      the new iterative point ${\mathbf{p}}_c$ fallen on the $f\left(\mathbf{x}\right)$.     }while (distance (the current ${\mathbf{p}}_c$, the previous ${\mathbf{p}}_c$)$>ϵ$).     ${\mathbf{p}}_{\Gamma }={\mathbf{p}}_c$;     Return ${\mathbf{p}}_{\Gamma }$;

Through a series of rigorous deductions, Comprehensive Algorithm A is the important algorithm of our paper. No matter how far the test point $\mathbf{p}$ is from the planar algebraic curve $f\left(\mathbf{x}\right)$, test point $\mathbf{p}$ could very robustly orthogonally projects onto the planar algebraic curve $f\left(\mathbf{x}\right)$. This has achieved our desired result. After a lot of testing and observation, when the point on the curve is close to the orthogonal projection point, we find that Comprehensive Algorithm A presents two characteristics: (1) difference between the first distance and the second distance decreases slower and slower, where the first distance and the second distance are the one between the previous iterative point ${\mathbf{p}}_c$ on the planar algebraic curve and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$, and the one between the current iterative point ${\mathbf{p}}_c$ on the planar algebraic curve and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$, respectively; (2) the rate goes even slower at which the absolute value of the inner product gradually approaches zero. These two characteristics are what we don’t want to obtain because they are contrary to the efficiency of computer systems. On the premise of ensuring robustness, we try our best to improve and excavate a certain degree of efficiency for the problem of point orthogonal projection onto planar algebraic curve. We have an ingenious discovery. After each running of Algorithm 3, we run the Newton’s iterative method associated with the original test point $\mathbf{p}$, which can improve the convergence and ensure the orthogonality. Namely, that is to add this step after the last step of the formula (17). Thus the iterative formula (18) is obtained.

$$\left\{\begin{array}{c}{\mathbf{y}}_n={\mathbf{x}}_n-\left(f\left({\mathbf{x}}_n\right)/〈\nabla f\left({\mathbf{x}}_n\right),\nabla f\left({\mathbf{x}}_n\right)〉\right)\nabla f\left({\mathbf{x}}_n\right),\hfill \\ {\mathbf{z}}_n={\mathbf{y}}_n-\left(\mathbf{F}\left({\mathbf{y}}_n\right)/〈\nabla \mathbf{F}\left({\mathbf{y}}_n\right),\nabla \mathbf{F}\left({\mathbf{y}}_n\right)〉\right)\nabla \mathbf{F}\left({\mathbf{y}}_n\right),\hfill \\ \mathbf{Q}=\mathbf{q}-(〈(\mathbf{q}-{\mathbf{z}}_n),\nabla f{(\mathbf{z}}_n)〉/〈\nabla f{(\mathbf{z}}_n),\nabla f{(\mathbf{z}}_n)〉)\nabla f{(\mathbf{z}}_n),\hfill \\ {\mathbf{u}}_n={\mathbf{z}}_n+sign\left(〈\mathbf{q}-{\mathbf{z}}_n,{\mathbf{t}}_0〉\right){\mathbf{t}}_0\Delta s,\hfill \\ {\mathbf{v}}_n={\mathbf{u}}_n-\left(\mathbf{F}\left({\mathbf{u}}_n\right)/〈\nabla \mathbf{F}\left({\mathbf{u}}_n\right),\nabla \mathbf{F}\left({\mathbf{u}}_n\right)〉\right)\nabla \mathbf{F}\left({\mathbf{u}}_n\right),\hfill \\ {\mathbf{w}}_n{=\mathbf{v}}_n+\left[-\Delta e,0\right]{\left[\nabla f^T,(\Delta {\mathbf{v}}_n{)}^T\right]}^{-1}(if\left|\left[\nabla f^T,{(\Delta {\mathbf{v}}_n)}^T\right]\right|=0,then{\mathbf{w}}_n={\mathbf{v}}_n),\hfill \\ {\mathbf{x}}_{n+1}={\mathbf{w}}_n-\left(\mathbf{G}\left({\mathbf{w}}_n\right)/〈\nabla \mathbf{G}\left({\mathbf{w}}_n\right),\nabla \mathbf{G}\left({\mathbf{w}}_n\right)〉\right)\nabla \mathbf{G}\left({\mathbf{w}}_n\right),\hfill \end{array}\right.$$

(18)

where ${\mathbf{F}}_0\left(\mathbf{x}\right)=\left[(\mathbf{q}-\mathbf{x})\times \nabla f\left(\mathbf{x}\right)\right]=0$, ${\displaystyle \mathbf{F}\left(\mathbf{x}\right)=\frac{{\mathbf{F}}_0\left(\mathbf{x}\right)}{{\displaystyle \sqrt{〈\nabla f\left(\mathbf{x}\right),\nabla f\left(\mathbf{x}\right)〉}}}}$, ${\displaystyle {\mathbf{t}}_0=\frac{\left[-f_y,f_x\right]}{∥\Delta f∥}}$, $\Delta s=∥\mathbf{Q}-{\mathbf{z}}_n∥$, $f\left({\mathbf{v}}_n\right)=\Delta e$, $\Delta {\mathbf{v}}_n=-\left({\mathbf{F}(\mathbf{u}}_n)/〈\nabla {\mathbf{F}(\mathbf{u}}_n),\nabla {\mathbf{F}(\mathbf{u}}_n)〉\right)\nabla {\mathbf{F}(\mathbf{u}}_n)$, ${\mathbf{G}}_0\left(\mathbf{x}\right)=\left[(\mathbf{p}-\mathbf{x})\times \nabla f\left(\mathbf{x}\right)\right]=0$, ${\displaystyle \mathbf{G}\left(\mathbf{x}\right)=\frac{{\mathbf{G}}_0\left(\mathbf{x}\right)}{{\displaystyle \sqrt{〈\nabla f\left(\mathbf{x}\right),\nabla f\left(\mathbf{x}\right)〉}}}}$. Because the iterative formula (17) of Algorithm 3 naturally becomes the iterative formula (18), so Algorithm 3 naturally becomes the following Algorithm 5.

Algorithm 5: The hybrid tangent vertical foot algorithm.

Input: The footpoint $\mathbf{q}$ and the planar algebraic curve $f\left(\mathbf{x}\right)=0$.     Output: The point ${\mathbf{p}}_c$ fallen on planar algebraic curve $f\left(\mathbf{x}\right)=0$.     Description:     Step 1:                 ${\mathbf{x}}_{n+1}=\mathbf{q}-(0.1,0.1)$;                 Do {                         ${\mathbf{x}}_n={\mathbf{x}}_{n+1}$;                         Execute ${\mathbf{x}}_{n+1}$ according to the iterative Equation (18);                 }while (${∥{\mathbf{x}}_{n+1}-{\mathbf{x}}_n∥}^2>\epsilon$&&$\left|f\left({\mathbf{x}}_{n+1}\right)\right|>\epsilon )$     Step 2:                 ${\mathbf{p}}_c={\mathbf{x}}_{n+1}$;                 Return ${\mathbf{p}}_c$;

Now let’s replace Algorithm 3 of Comprehensive Algorithms A with Algorithm 5. We get the following Comprehensive Algorithm B (Algorithm 6).

Algorithm 6: The second improved curvature circle algorithm (Comprehensive Algorithm B).

Input: Test point $\mathbf{p}$ and the planar algebraic curve $f\left(\mathbf{x}\right)$.     Output: Orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the test point $\mathbf{p}$ which orthogonally projects the test      point $\mathbf{p}$ onto the planar algebraic curve $f\left(\mathbf{x}\right)$.     Description:     Step 1: Starting from the neighbor point of the test point $\mathbf{p}$, calculate the point ${\mathbf{p}}_c$ fallen on the $f\left(\mathbf{x}\right)$ via Algorithm 1.     Do{               Step 2: Compute the footpoint $\mathbf{q}$ via Algorithm 2.               Step 3: Project footpoint $\mathbf{q}$ onto the planar algebraic curve $f\left(\mathbf{x}\right)$ via Algorithm 5, then get      new point ${\mathbf{p}}_c$ fallen on the $f\left(\mathbf{x}\right)$.     }while(distance(the current ${\mathbf{p}}_c$, the previous ${\mathbf{p}}_c$)$>ϵ$).     ${\mathbf{p}}_{\Gamma }={\mathbf{p}}_c$;     Return ${\mathbf{p}}_{\Gamma }$;

Comprehensive Algorithm A and Comprehensive Algorithm B share common advantage: the robustness of the two algorithms is substantially improved than that of the existing algorithms because our algorithms are not subject to any restrictions on test points and initial iteration points. By comparison, Comprehensive Algorithm B has four advantages over Comprehensive Algorithm A. (1) The last step of the iterative formula (18) in Comprehensive Algorithm B can make corrections continuously; (2) The last step of the iterative formula (18) in Comprehensive Algorithm B accelerates the whole Comprehensive Algorithm B; (3) The last step of the iterative formula (18) in Comprehensive Algorithm B accelerates the inner product of two vectors to 0, where the first vector refers to the vector connecting the test point $\mathbf{p}$ and the iteration point ${\mathbf{z}}_{n+1}$ of Comprehensive Algorithm B and the second vector ${\left[-\frac{\partial f}{\partial y},\frac{\partial f}{\partial x}\right]}_{\left|\mathbf{x}={\mathbf{x}}_{n+1}\right.}$ is the tangent vector derived from the iteration point ${\mathbf{x}}_{n+1}$ on the planar algebraic curve, respectively; (4) Comprehensive Algorithm B overcomes two shortcomings of Comprehensive Algorithm A.

Of course, when the test point is not too far from the plane algebra curve, Comprehensive Algorithm is also convergent for any initial iterative point. However, Comprehensive Algorithm A takes more time than directly using the hybrid second order algorithm. In practical applications such as computer graphics, it’s hard to know if the test point $\mathbf{p}$ is close to or far from a planar algebraic curve. Because the main reason is that the degree and the type of the planar algebraic curve restrict the relative distance between the test point $\mathbf{p}$ and the planar algebraic curve. In order to optimize time efficiency, we take advantage of Comprehensive Algorithm A and the hybrid second order algorithm such that no matter where the test point $\mathbf{p}$ is located, it can be orthogonally projected onto the planar algebraic curve efficiently and robustly. First, the hybrid second order algorithm is iterated. If it does not converge after 100 iterations, it will be changed to Comprehensive Algorithm A to iterate until the iteration point reaches the orthogonal projection point ${\mathbf{p}}_{\Gamma }$. Specific algorithm implementation is the following Comprehensive Integrated Algorithm A (Algorithm 7).

Algorithm 7: The first comprehensive integrated improved curvature circle algorithm (Comprehensive Integrated Algorithm A).

Input:Test point $\mathbf{p}$ and the planar algebraic curve $f\left(\mathbf{x}\right)$.     Output: Orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the test point $\mathbf{p}$.     Description:     Step 1:                 ${\mathbf{x}}_{n+1}=\mathbf{p}-(0.1,0.1)$;                 for($i=0;i<N;i++$) {                         ${\mathbf{x}}_n={\mathbf{x}}_{n+1}$;                         ${\mathbf{x}}_{n+1}$=Hybrid second order algorithm$(f,\mathbf{p},{\mathbf{x}}_n)$;                         if$(∥{\mathbf{x}}_{n+1}-{\mathbf{x}}_n∥<\epsilon )$ break ;                 }     Step 2:                 if($i\ge N\&\&d\ge 1e-15)$ {                         ${\mathbf{x}}_n={\mathbf{x}}_{n+1}$ ;                         ${\mathbf{x}}_{n+1}$=Comprehensive Algorithm A$(f,\mathbf{p},{\mathbf{x}}_n)$;                 }                 ${\mathbf{p}}_{\Gamma }={\mathbf{x}}_{n+1}$;                 Return ${\mathbf{p}}_{\Gamma }$;

Number N is an empirical value of the iterative times where the value N is specified as 5 or 6.

Similar to Comprehensive Algorithm A, by replacing Algorithm 3 with Algorithm 5, the following Comprehensive Integrated Algorithm B (Algorithm 8) can be obtained naturally.

Algorithm 8: The second comprehensive integrated improved curvature circle algorithm (Comprehensive Integrated Algorithm B).

Input: The test point $\mathbf{p}$ and the planar algebraic curve $f\left(\mathbf{x}\right)$.     Output: Orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the test point $\mathbf{p}$.     Description:     Step 1:                 ${\mathbf{x}}_{n+1}=\mathbf{p}-(0.1,0.1)$;                 for($i=0;i<N;i++$) {                         ${\mathbf{x}}_n={\mathbf{x}}_{n+1}$;                         ${\mathbf{x}}_{n+1}$=Hybrid second order algorithm$(f,\mathbf{p},{\mathbf{x}}_n)$;                         if$(∥{\mathbf{x}}_{n+1}-{\mathbf{x}}_n∥<\epsilon )$ break ;                 }     Step 2:                 if($i\ge N\&\&d\ge 1e-15)$ {                         ${\mathbf{x}}_n={\mathbf{x}}_{n+1}$ ;                         ${\mathbf{x}}_{n+1}$=Comprehensive Algorithm B$(f,\mathbf{p},{\mathbf{x}}_n)$;                 }                 ${\mathbf{p}}_{\Gamma }={\mathbf{x}}_{n+1}$;                 Return ${\mathbf{p}}_{\Gamma }$;

Number N is an empirical value of the iterative times where the value N is specified as 5 or 6.

To sum up, we have presented four synthesis algorithms altogether. After analysis and judgment, Comprehensive Algorithm B and Comprehensive Integrated Algorithm B are the most robust and efficient. On the problem of orthogonal projection of point onto planar algebraic curve, if the distance between the test point and the planar algebraic curve is close, we recommend the hybrid second order algorithm, if the distance between the test point and the planar algebraic curve is not close, we recommend Comprehensive Algorithm B. Of course, if the distance between the test point and the planar algebraic curve cannot be known to be very far or close, Comprehensive Integrated Algorithm B is the best choice.

Remark 2.

In sum, Comprehensive Algorithm B has strong superiority over existing algorithms [10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28]. If the distance between the test point and the planar algebraic curve is very far away, the test point can be ideally orthogonally projected onto the planar algebraic curve. But when there are singular points${\displaystyle \frac{\partial f}{\partial x}·\frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}·\frac{\partial f}{\partial y}=0}$in the planar algebraic curve, this case will seriously hinder the correct execution and implementation of Comprehensive Algorithm B. In order to solve the problem in the case of singularities in the planar algebraic curves, we propose a remedy to Comprehensive Algorithm B (Algorithm 9). The specific description is as follows (See Figure 5).

Algorithm 9: The first remedial algorithm of Comprehensive Algorithm B.

Input: Test point $\mathbf{p}$ and the planar algebraic curve $f\left(\mathbf{x}\right)$.     Output: Orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the test point $\mathbf{p}$.     Description:     Step 1.                 Starting from the neighbor point of the test point $\mathbf{p}$, calculate the iterative point ${\mathbf{p}}_c$ fallen     on the planar algebraic curve $f\left(\mathbf{x}\right)$ via Algorithm 1.     Step 2.                 Judge whether to use curvature circle method or tangent method in the next step.     Step 3.                 Find the left endpoint ${\mathbf{L}}_0$ on the other side of $f\left(\mathbf{x}\right)$ relative to the test point $\mathbf{p}$. According     to the result of step 2, if use curvature circle method, then the left endpoint ${\mathbf{L}}_0$ is equal to the     intersection point $\mathbf{q}$ which is computed by the curvature circle method with the formula (16).     If not, then the left endpoint ${\mathbf{L}}_0$ is equal to the vertical foot $\mathbf{Q}$ which is computed by the     tangent method with the third step of the formula (17).     Step 4.                 Calculate the intersection point ${\mathbf{p}}_c$ of the line segment $\overline{{\mathbf{L}}_0\mathbf{p}}$ connecting the current left     endpoint ${\mathbf{L}}_0$ and the test point $\mathbf{p}$ and the planar algebraic curve $f\left(\mathbf{x}\right)$ by the hybrid method of     combining Newton’s iterative method and binary search method. The intersection point ${\mathbf{p}}_c$ is     called as the current iterative point ${\mathbf{p}}_c$;     Step 5.                 Repeat Step 2,Step 3 and Step 4 until the distance between the current iterative point ${\mathbf{p}}_c$     and the previous iterative point ${\mathbf{p}}_c$ is near zero;     Step 6.                 ${\mathbf{p}}_{\Gamma }$= ${\mathbf{p}}_c$;                 Return ${\mathbf{p}}_{\Gamma }$;

Now let’s describe the hybrid method of combining Newton’s iterative method and binary search method in detail. The parameter equation of the line segment $\overline{{\mathbf{L}}_0\mathbf{p}}$ can be expressed as

$$\left\{\begin{array}{c}x=L_1+(p_1-L_1)w,\hfill \\ y=L_2+(p_2-L_2)w,\hfill \end{array}\right.$$

(19)

where ${\mathbf{L}}_0=(L_1,L_2)$, $\mathbf{p}=(p_1,p_2)$, and $0\le w\le 1$ is a parameter of Equation (19). Substitute Equation (19) into Equation (6) of the planar algebraic curve to get a equation on the parameter w,

$$K\left(w\right)=f(x,y)=0,$$

(20)

where the x and y of Equation (20) are completely determined by the x and y of Equation (19). So the most basic Newton’s iterative formula corresponding to Equation (20) is not difficult to write as,

$$w_{n+1}=w_n-\frac{K\left(w_n\right)}{DK\left(w_n\right)},$$

(21)

where $DK\left(w\right)$ is the first derivative of $K\left(w\right)$ about the parameter w. Now we start to iterate the Newton’s iterative formula (21) with the initial iterative value $w_0=0.0$. Based on the actual situation, the intersection of the line segment $\overline{{\mathbf{L}}_0\mathbf{p}}$ and the planar algebraic curve is much closer to the left endpoint ${\mathbf{L}}_0$ and much farther from the original test point $\mathbf{p}$, therefore, the initial interval of the binary search method can be specified as $\left[a,b\right]=\left[0.0,0.5\right]$. The detailed description of the hybrid method of combining Newton’s iterative method and binary search method is as following Algorithm 10.

Algorithm 10: The hybrid method of combining Newton’s iterative method and binary search method.

Input: The planar algebraic curve $f\left(\mathbf{x}\right)$, the original test point $\mathbf{p}=\left(p_1,p_2\right)$, the iterative point ${\mathbf{p}}_c$ via Algorithm 1.     Output: The intersection ${\mathbf{p}}_c$ between the line segment $\overline{{\mathbf{L}}_0\mathbf{p}}$ and the planar algebraic curve $f\left(\mathbf{x}\right)$.     Description:     Step 1:                 The initial interval of the binary search method $[a,b]=[0.0,0.5]$, the initial iterative     value $w=0.0$;       Step 2:                 $w=w-K\left(w\right)/DK\left(w\right)$;                 kmin=min($K\left(a\right)$,$K\left(b\right)$);                 kmax=max($K\left(a\right)$,$K\left(b\right)$);                 if ($K\left(w\right)$ < kmin or $K\left(w\right)$ > kmax)                         $w=\left(a+b\right)/2$;                 $sa$=sign($K\left(a\right)$);                 $sw$=sign($K\left(w\right)$);                 if($sa==sw$)                         $a=w$;                 else                         $b=w$;      Step 3:                 Repeatedly iterate Step 2 until $\left|a-b\right|<\epsilon$;      Step 4:                 ${\mathbf{p}}_c={\mathbf{L}}_0+(\mathbf{p}-{\mathbf{L}}_0)w$;                 Return ${\mathbf{p}}_c$;

The robustness of the first remedial algorithm of Comprehensive Algorithm B is much better than that of Comprehensive Algorithm B while the first remedial algorithm of Comprehensive Algorithm B takes much more time than Comprehensive Algorithm B. The hybrid method of combining Newton’s iterative method and binary search method is a hybrid method which binary search method ensures global convergence and the Newton’s iterative method plays an accelerating role. In order to ensure robustness and improve efficiency, we have fully excavated Comprehensive Algorithm B. We have developed Second Remedial Algorithm (Algorithm 11). The specific description is as follows (See Figure 6).

Algorithm 11: Second Remedial Algorithm.

Input: Test point $\mathbf{p}$ and the planar algebraic curve $f\left(\mathbf{x}\right)$.     Output: Orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the test      point $\mathbf{p}$ which orthogonally projects the test point $\mathbf{p}$ onto the planar algebraic curve $f\left(\mathbf{x}\right)$     Description:     Step 1: Starting from a certain percentage of the test point $\mathbf{p}$, calculate the point ${\mathbf{p}}_c$ fallen on the      $f\left(\mathbf{x}\right)$ via Algorithm 1.     Do{                 Step 2: Compute the footpoint $\mathbf{q}$ via Algorithm 2.                 Step 3: Starting from the footpoint $\mathbf{q}$, compute the iterative point ${\mathbf{p}}_c$ fallen on the $f\left(x\right)$ via      Algorithm 1.     }while(distance(the current ${\mathbf{p}}_c$, the previous ${\mathbf{p}}_c$)$>ϵ$).     Step 4: Compute the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the test point $\mathbf{p}$ via Algorithm 12.     Return ${\mathbf{p}}_{\Gamma }$;

Algorithm 12: The hybrid Newton-type iterative algorithm.

Input: The current iterative point ${\mathbf{p}}_c$ fallen on the planar      algebraic curve $f\left(\mathbf{x}\right)$ and the planar algebraic curve $f\left(\mathbf{x}\right)$.     Output: Orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the test point $\mathbf{p}$ which orthogonally projects the test      point $\mathbf{p}$ onto the planar algebraic curve $f\left(\mathbf{x}\right)$.     Description:     Step 1:                 ${\mathbf{x}}_{n+1}={\mathbf{p}}_c$;                 Do {                         ${\mathbf{x}}_n={\mathbf{x}}_{n+1}$;                         Compute ${\mathbf{x}}_{n+1}$ according to the iterative formula (22);                 }while (${∥{\mathbf{x}}_{n+1}-{\mathbf{x}}_n∥}^2>\epsilon$&&$\left|f\left({\mathbf{x}}_{n+1}\right)\right|>\epsilon )$     Step 2:                 ${\mathbf{p}}_{\Gamma }={\mathbf{x}}_{n+1}$;                 Return ${\mathbf{p}}_{\Gamma }$;

The expression of the iterative formula (22) is as follow,

$$\left\{\begin{array}{c}{\mathbf{y}}_n={\mathbf{x}}_n-\left(f\left({\mathbf{x}}_n\right)/〈\nabla f\left({\mathbf{x}}_n\right),\nabla f\left({\mathbf{x}}_n\right)〉\right)\nabla f\left({\mathbf{x}}_n\right),\hfill \\ {\mathbf{z}}_n={\mathbf{y}}_n-\left(\mathbf{F}\left({\mathbf{y}}_n\right)/〈\nabla \mathbf{F}\left({\mathbf{y}}_n\right),\nabla \mathbf{F}\left({\mathbf{y}}_n\right)〉\right)\nabla \mathbf{F}\left({\mathbf{y}}_n\right),\hfill \\ {\mathbf{x}}_{n+1}={\mathbf{z}}_n+\left[-\Delta e,0\right]{\left[\nabla f^T,(\Delta {\mathbf{z}}_n{)}^T\right]}^{-1}(if\left|\left[\nabla f^T,{(\Delta {\mathbf{z}}_n)}^T\right]\right|=0,then{\mathbf{x}}_{n+1}={\mathbf{z}}_n),\hfill \end{array}\right.$$

(22)

where ${\mathbf{F}}_0\left(\mathbf{x}\right)=\left[(\mathbf{p}-\mathbf{x})\times \nabla f\left(\mathbf{x}\right)\right]=0$, ${\displaystyle \mathbf{F}\left(\mathbf{x}\right)=\frac{{\mathbf{F}}_0\left(\mathbf{x}\right)}{{\displaystyle \sqrt{〈\nabla f\left(\mathbf{x}\right),\nabla f\left(\mathbf{x}\right)〉}}}}$, $f\left({\mathbf{z}}_n\right)=\Delta e$, $\Delta {\mathbf{z}}_n=-\left({\mathbf{F}(\mathbf{y}}_n)/〈\nabla {\mathbf{F}(\mathbf{y}}_n),\nabla {\mathbf{F}(\mathbf{y}}_n)〉\right)\nabla {\mathbf{F}(\mathbf{y}}_n)$.

Remark 3.

In this remark, we present the geometric interpretation of Second Remedial Algorithm. The purpose of the first step is to make the iterative point${\mathbf{p}}_c$fall on the planar algebraic curve as much as possible through Newton’s steepest gradient descent method of Algorithm 1, where the coordinates of the initial iterative point take proportional value of that of the test point$\mathbf{p}$to ensure that Algorithm 1 converges successfully. Otherwise, the distance between the initial iterative point and the planar algebraic curve is very large, which easily leads to the divergence of Algorithm 1. The purpose of Do … While cycle body in Second Remedial Algorithm is to continuously and gradually move the iterative point${\mathbf{p}}_c$to fall on the planar algebraic curve to the orthogonal projection point${\mathbf{p}}_{\Gamma }$. The second step inDo … Whilecycle body in Second Remedial Algorithm has two characteristics. Since the footpoint$\mathbf{q}$is the unique intersection point of the curvature circle and the straight line segment$\overline{\mathbf{mp}}$connecting the centre$\mathbf{m}$of the curvature circle and the test point$\mathbf{q}$, the footpoint$\mathbf{q}$is always closely related to the iterative point${\mathbf{p}}_c$fallen on the planar algebraic curve and the test point$\mathbf{p}$. The first characteristic is that the footpoint$\mathbf{q}$can guarantee the global convergence of the whole algorithm (Second Remedial Algorithm). The second characteristic is that the distance between the footpoint$\mathbf{q}$and the planar algebraic curve is much smaller than the distance between the test point$\mathbf{p}$and the planar algebraic curve. So the third step with Algorithm 1 in Do… Whilecycle body can very robustly iterate the footpoint$\mathbf{q}$onto the planar algebraic curve. The core thought of Do … While cycle body in Second Remedial Algorithm is to keep the iterative point${\mathbf{p}}_c$to fall on the planar algebraic curve and to move towards the orthogonal projection point${\mathbf{p}}_{\Gamma }$such that the distance$∥{\mathbf{p}}_c-{\mathbf{p}}_{\Gamma }∥$between the iterative point${\mathbf{p}}_c$and the orthogonal projection point${\mathbf{p}}_{\Gamma }$becomes smaller and smaller. As the distance$∥{\mathbf{p}}_c-{\mathbf{p}}_{\Gamma }∥$gets smaller and smaller, we have found that there is a defect in Do … While cycle body in Second Remedial Algorithm. The decreasing speed of the distance$∥{\mathbf{p}}_c-{\mathbf{p}}_{\Gamma }∥$is getting slower and slower. Especially the second formula of the formula (8) is very difficult to be satisfied. Namely, it is difficult to orthogonalize the vector$\overrightarrow{\mathbf{p}{\mathbf{p}}_c}$and the vector$\nabla f\left({\mathbf{p}}_c\right)$. In order to overcome the difficulty, we add Algorithm 12 behind Do … While cycle body in Second Remedial Algorithm. Algorithm 12 includes three components: (1) The Newton’s steepest gradient descent method in the first step; (2) The Newton’s iterative method closely associated with the test point$\mathbf{p}$in the second step; (3) Correcting method in the third step. Algorithm 12 has four advantages and important roles: (1) Algorithm 12 plays a role for accelerating the whole algorithm (Second Remedial Algorithm); (2) The first step plays a role for making the iteration point fall on the planar algebraic curve as far as possible; (3) The second step plays a role for accelerating orthogonalization between the vector$\overrightarrow{\mathbf{p}{\mathbf{p}}_c}$and the vector$\nabla f\left({\mathbf{p}}_c\right)$; (4) The third step plays a dual role for the accelerating orthogonalization and the promoting the iterative point to fall on the planar algebraic curve. The numerical tests of Second Remedial Algorithm achieve our expected results. No matter how far the test point$\mathbf{p}$is from the planar algebraic curve$f\left(\mathbf{x}\right)$, Second Remedial Algorithm can converge very robustly and efficiently. Second Remedial Algorithm is the best one in our paper. Of course, the robustness and the efficiency of Second Remedial Algorithm are better than that of the existing algorithms. We are very happy about that.

Remark 4.

In order to further improve the efficiency of the test point$\mathbf{p}$orthogonal projecting onto plane algebraic curve$f\left(\mathbf{x}\right)$, we construct a Comprehensive Integrated Algorithm C which includes two parts: the hybrid second order algorithm in [28] and Second Remedial Algorithm. Firstly run the hybrid second order algorithm in [28]. If the hybrid second order algorithm converges, then it means that Comprehensive Integrated Algorithm C is finished. Otherwise, Second Remedial Algorithm runs until it converges successfully. That is the end of Comprehensive Integrated Algorithm C. The specific description of Comprehensive Integrated Algorithm C is similar to that of Comprehensive Integrated Algorithm B. Here, we are not giving a detailed description of Comprehensive Integrated Algorithm C. When the distance between the test point$\mathbf{p}$and the planar algebraic curve$f\left(\mathbf{x}\right)$is not far, the hybrid second order algorithm in [28] has very high robustness and efficiency. When the distance between the test point$\mathbf{p}$and the planar algebraic curve$f\left(\mathbf{x}\right)$is particularly far, the hybrid second order algorithm does not converge and fails, while Second Remedial Algorithm converges particularly robustly and successfully. To sum up, Comprehensive Integrated Algorithm C absorbs the advantages of two sub-algorithms and overcomes their respective shortcomings such that the robustness and the efficiency of Comprehensive Integrated Algorithm C are maximized.

3. Convergence Analysis

This section proves that several Comprehensive Algorithms do not depend on the initial iteration points.

Theorem 1.

Comprehensive Algorithm A is independent of the initial iterative point.

Proof. 

Firstly, we state the whole operation process of Comprehensive Algorithm A. Comprehensive Algorithm A contains three sub-algorithms (Algorithms 1–3). The role of Algorithm 1 is to repeatedly iterate the iterative formula (3) through Newton’s steepest gradient descent method such that the final iteration point ${\mathbf{x}}_{n+1}$ could fall on the planar algebraic curve where the final iteration point ${\mathbf{x}}_{n+1}$ is denoted as ${\mathbf{p}}_c$. The function of Algorithm 2 is to seek the footpoint $\mathbf{q}$. The curvature circle ◯ determined by the point ${\mathbf{p}}_c$ is obtained from the iterative point ${\mathbf{p}}_c$ on the planar algebraic curve $f\left(\mathbf{x}\right)$ of Algorithm 1, where the curvature k, the radius R and the center $\mathbf{m}$ are determined by formulas (11)–(13), respectively. The intersection of the line segment $\overline{\mathbf{mp}}$ connecting the center $\mathbf{m}$ and the test point $\mathbf{p}$ and the curvature circle ◯ is foot point $\mathbf{q}$. The footpoint $\mathbf{q}$ could be orthogonally projected onto the planar algebraic curve $f\left(\mathbf{x}\right)$ by repeated iteration of Algorithm 3 where at this moment the test point is not the original test point $\mathbf{p}$ but the footpoint point $\mathbf{q}$ solved by Algorithm 2. Repeatedly run Algorithm 2 and Algorithm 3 bound together until the distance between the current footpoint $\mathbf{q}$ and the previous footpoint $\mathbf{q}$ is near zero.

Secondly, the Comprehensive Algorithm A is independent of the initial iterative point. No matter how far the original test point $\mathbf{p}$ is from the planar algebraic curve $f\left(\mathbf{x}\right)$, no matter where the initial iterative point ${\mathbf{x}}_0$ is located, Algorithms 1 can ensure that the final iterative point ${\mathbf{x}}_{n+1}$ or ${\mathbf{p}}_c$ of the initial iterative point can fall on the planar algebraic curve $f\left(\mathbf{x}\right)$. It is obvious that the distance $∥{\mathbf{p}}_c-{\mathbf{p}}_{\Gamma }∥$ between the iteration point ${\mathbf{p}}_c$ and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ is much smaller than the distance $∥\mathbf{p}-{\mathbf{p}}_{\Gamma }∥$ between the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ and the original test point $\mathbf{p}$. From the iterative point ${\mathbf{p}}_c$ fallen on the planar algebraic curve $f\left(\mathbf{x}\right)$, we can calculate the corresponding curvature k and its center $\mathbf{m}$ and radius R. The intersection point $\mathbf{q}$ of the curvature circle ◯ and the line segment $\overline{\mathbf{mp}}$ connecting the original test point $\mathbf{p}$ and the center $\mathbf{m}$ of the curvature circle ◯ is just the footpoint $\mathbf{q}$. That is to say, the footpoint $\mathbf{q}$ is directly generated by the curvature circle ◯ and the line segment $\overline{\mathbf{mp}}$, and the curvature circle ◯ is controlled by the iterative point ${\mathbf{p}}_c$ fallen on the planar algebraic curve $f\left(\mathbf{x}\right)$. So the footpoint $\mathbf{q}$ is directly controlled by the original test point $\mathbf{p}$ and the iterative point ${\mathbf{p}}_c$, while the current footpoint $\mathbf{q}$ is between the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ and the current iterative point ${\mathbf{p}}_c$. It also shows that Algorithm 2 plays a decisive role in the convergence robustness of Comprehensive Algorithm. In addition, we can also know that the distance between the footpoint point $\mathbf{q}$ and the planar algebraic curve $f\left(\mathbf{x}\right)$ is much smaller than the distance between the original test point $\mathbf{p}$ and the planar algebra curve $f\left(\mathbf{x}\right)$. At this point, we keep running Algorithm 3 with the footpoint point $\mathbf{q}$ as the current test point until the current test point can be orthogonally projected onto the plane algebraic curve $f\left(\mathbf{x}\right)$ with guaranteed convergence of Algorithm 3. And now we can call the orthogonal projection point of the footpoint point $\mathbf{q}$ as also the iterative point ${\mathbf{p}}_c$ fallen on the planar algebraic curve $f\left(\mathbf{x}\right)$. The first reason is the distance between the current iterative point ${\mathbf{p}}_c$ and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the original test point point $\mathbf{p}$ is smaller than the one between the previous iterative point ${\mathbf{p}}_c$ and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the original test point $\mathbf{p}$. The second reason is that it establishes a solid foundation for the convergence robustness of the subsequent sub-algorithms implementation. Then according to the requirements of Comprehensive Algorithm A, the second step and third step of Comprehensive Algorithm A are executed once per cycle, the distance $∥{\mathbf{p}}_c-{\mathbf{p}}_{\Gamma }∥$ between the current iterative point ${\mathbf{p}}_c$ on the planar algebraic curve and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the original test point $\mathbf{p}$ of the execution result is smaller than that between the previous iterative point ${\mathbf{p}}_c$ on the planar algebraic curve $f\left(\mathbf{x}\right)$ and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the original test point $\mathbf{p}$. The distance $∥{\mathbf{p}}_c-{\mathbf{p}}_{\Gamma }∥$ between the current iterative point ${\mathbf{p}}_c$ and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the original test point $\mathbf{p}$ is becoming smaller. So repeatedly iterate the second step and the third step of Comprehensive Algorithm A until the distance $∥{\mathbf{p}}_c-{\mathbf{p}}_{\Gamma }∥$ between the current iterative point ${\mathbf{p}}_c$ and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the original test point $\mathbf{p}$ is becoming smaller and smaller. Ultimately, the distance $∥{\mathbf{p}}_c-{\mathbf{p}}_{\Gamma }∥$ between the current iterative point ${\mathbf{p}}_c$ and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the original test point $\mathbf{p}$ is becoming zero. It also demonstrates that Comprehensive Algorithm A is completely convergent. This further proves that Comprehensive Algorithm A can completely converge no matter how far away the original test point $\mathbf{p}$ is from the planar algebraic curve and no matter where the initial iterative point ${\mathbf{x}}_0$ of Comprehensive Algorithm A is on the plane. This means Comprehensive Algorithm A is independent of the initial iterative point. □

Theorem 2.

Comprehensive Algorithm B is independent of the initial iterative point.

Proof. 

In the last step of the iterative formula (18) in Algorithm 5, Newton’s iteration method, which is closely related to the original test point $\mathbf{p}$, is added. In this way, the iterative formula (17) is transformed into the iterative formula (18) in Algorithm 5. Algorithm 5 has several advantages over Algorithm 3. It can speed up the iteration, improve its accuracy and promote the orthogonalization of the tangent vector derived from the iteration point on the planar algebraic curve and the tangent vector connecting the test point and the iterative point. Replace Algorithm 3 of Comprehensive Algorithm A with Algorithm 5 to get Comprehensive Algorithm B. Since Comprehensive Algorithm A is independent of the initial iterative point, so Comprehensive Algorithm B is naturally independent of the initial iterative point. □

Theorem 3.

Comprehensive Integrated Algorithm A is independent of the initial iterative point.

Proof. 

Comprehensive Integrated Algorithm A consists of two parts: the hybrid second order algorithm and Comprehensive Algorithm A. Whether the test point is very far or very close to the planar algebraic curve, the hybrid second order algorithm is executed several times. If this algorithm converges, then it represents that the execution of Comprehensive integrated Algorithm A is over. So Comprehensive Integrated Algorithm A is independent of the initial iterative point. If the hybrid second order algorithm does not converge, then run Comprehensive Algorithm A of the second step of Comprehensive Integrated Algorithm A. Because whether the test point is very far from or very close to the planar algebraic curve, we know from Theorem 1 that Comprehensive Algorithm A is independent of the initial iterative point. To sum up, Comprehensive Integrated Algorithm A is independent of the initial iterative point. □

In a similar way to the proof of Theorem 3, we can state the following result.

Theorem 4.

Comprehensive Integrated Algorithm B is independent of the initial iterative point.

Theorem 5.

The first remedial algorithm of Comprehensive Algorithm B is independent of the initial iterative point.

Proof. 

From the Figure 5, for any initial iterative point, the final iterative point ${\mathbf{p}}_c$ of Algorithm 1 in the first step of the first remedial algorithm of Comprehensive Algorithm B can ensure that it falls on the planar algebraic curve $f\left(\mathbf{x}\right)$. The left endpoint ${\mathbf{L}}_0$ is the only one that can be determined through third step of the first remedial algorithm of Comprehensive Algorithm B. Graphic display shows that the left endpoint ${\mathbf{L}}_0$ and the original test point $\mathbf{p}$ are on both sides of the planar algebraic curve. Namely, there is only one intersection point (also written as ${\mathbf{p}}_c)$ between the line segment $\overline{{\mathbf{L}}_0\mathbf{p}}$ and the planar algebraic curve $f\left(\mathbf{x}\right)$. Because the hybrid method of combining Newton’s iterative method and binary search method is global convergence method, the intersection ${\mathbf{p}}_c$ of the line segment $\overline{{\mathbf{L}}_0\mathbf{p}}$ and the planar algebraic curve $f\left(\mathbf{x}\right)$ can be accurately and uniquely solved by this method. Then repeatedly iterate and run Step 2, Step 3 and Step 4, the distance $∥{\mathbf{p}}_c-{\mathbf{p}}_{\Gamma }∥$ between the current intersection point ${\mathbf{p}}_c$ and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ of the original test point $\mathbf{p}$ continues to shrink to zero. So we have this conclusion that the first remedial algorithm of Comprehensive Algorithm B is independent of the initial iterative point. □

Theorem 6.

Second Remedial Algorithm is independent of the initial iterative point.

Proof. 

In Remark 3, we give a detailed description of the geometric interpretation of Second Remedial Algorithm. In this proof, we only elaborate on the most important geometric significance of Second Remedial Algorithm. The first step of Second Remedial Algorithm is to let the initial iteration point fall on the planar algebraic curve as much as possible through Newton’s steepest gradient descent method of Algorithm 1. Moreover, there is few restriction on the selection of the initial iterative point. The purpose of Do …While cycle body in Second Remedial Algorithm is to continuously and gradually move the iterative point ${\mathbf{p}}_c$ to fall on the planar algebraic curve to the orthogonal projection point ${\mathbf{p}}_{\Gamma }$. The second step in Do…While cycle body in Second Remedial Algorithm has two characteristics. Since the footpoint $\mathbf{q}$ is the unique intersection point of the curvature circle and the straight line segment $\overline{\mathbf{mp}}$ connecting the centre $\mathbf{m}$ of the curvature circle and the test point $\mathbf{p}$, the footpoint $\mathbf{q}$ is always closely related to the iterative point ${\mathbf{p}}_c$ fallen on the planar algebraic curve and the test point $\mathbf{p}$. The first characteristic is that the footpoint $\mathbf{q}$ can guarantee the global convergence of the whole algorithm (Second Remedial Algorithm). The second characteristic is that the distance between the footpoint $\mathbf{q}$ and the planar algebraic curve is much smaller than the distance between the test point $\mathbf{p}$ and the planar algebraic curve. So the third step with Algorithm 1 in Do …While cycle body can very robustly iterate the footpoint $\mathbf{q}$ onto the planar algebraic curve. The core thought of Do …While cycle body in Second Remedial Algorithm is to keep the iterative point ${\mathbf{p}}_c$ fallen on the planar algebraic curve moving towards the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ such that the distance $∥{\mathbf{p}}_c-{\mathbf{p}}_{\Gamma }∥$ between the iterative point ${\mathbf{p}}_c$ and the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ becomes smaller and smaller. If the distance $∥{\mathbf{p}}_c-{\mathbf{p}}_{\Gamma }∥$ gets smaller and smaller, we have found that the decreasing speed of the distance $∥{\mathbf{p}}_c-{\mathbf{p}}_{\Gamma }∥$ is getting slower and slower. Especially the second formula of the formula (8) is very difficult to be satisfied. Algorithm 12 behind the loop body has four advantages and important roles: (1) Algorithm 12 plays a role for accelerating the whole algorithm (Second Remedial Algorithm); (2) The first step plays a role for making the iteration point fall on the planar algebraic curve as far as possible; (3) The second step plays a role for accelerating orthogonalization between the vector $\overrightarrow{{\mathbf{pp}}_c}$ and the vector $\nabla f\left({\mathbf{p}}_c\right)$; (4) The third step plays a dual role for the accelerating orthogonalization and the promoting the iterative point to fall on the planar algebraic curve. No matter how far the test point is from the planar algebraic curve, Second Remedial Algorithm converges very robustly and efficiently. By adding this step, the efficiency and the robustness for Algorithm 12 of Second Remedial Algorithm is further improved. Then the robustness and the efficiency of Second Remedial Algorithm is also further improved. So Second Remedial Algorithm is independent of the initial iterative point. In addition, in a similar way to the proof of Theorem 3, it is not difficult to know that Comprehensive Integrated Algorithm C is also independent of the initial iterative point. □

4. Numerical Comparison Results

We now present some examples to illustrate the efficiency and the comparison of the newly developed method of Comprehensive Algorithm B and Second Remedial Algorithm to show the two algorithms’ high robustness and efficiency for very remote test points. We have three examples to represent closed planar algebraic curve, two sub-closed planar algebraic curves, two branches but not closed planar algebra curves and a single branch not closed the planar algebra curve, respectively. All computations were done using VC++6.0. We used $\epsilon ={10}^{-16}$. The following stopping criteria is used for Comprehensive Algorithm B and Second Remedial Algorithm . In

Table 1. Test results of Comprehensive Algorithm B for Example 1.

p	(1325, 7447)	(779, 325)	(990, 1375)	(0.59, 1377)	(−0.5, 8623)	(−16, 598)	(−3, 231)	(−21, 247)
${\mathbf{p}}_{\Gamma }$	(3.1087, −1.8666)	(3.3647, −1.8184)	(3.3695, −1.8380)	(0.60, 1.140)	(−0.01582, 1.13368)	(−0.3353, 1.13074)	(−0.2268, 1.1328)	(−0.6071, 1.1166)
$\left|f\left({\mathbf{p}}_{\Gamma }\right)\right|$	$3.9·{10}^{-10}$	$7.9·{10}^{-11}$	$1.4·{10}^{-11}$	$2.7·{10}^{-12}$	$2.3·{10}^{-12}$	$4.0·{10}^{-12}$	$1.9·{10}^{-12}$	$2.3·{10}^{-12}$
$\left|〈{\mathbf{V}}_1,{\mathbf{V}}_2〉\right|$	$2.3·{10}^{-10}$	$4.4·{10}^{-11}$	$2.9·{10}^{-09}$	$4.8·{10}^{-12}$	$2.1·{10}^{-13}$	$1.5·{10}^{-12}$	$1.4·{10}^{-11}$	$9.6·{10}^{-12}$
$\mathbf{p}$	(−42, −127)	(−5, −38)	(−9,−579)	(−537, −11)	(78, −123)	(168, −12)	(537, −31)	(91, −221)
${\mathbf{p}}_{\Gamma }$	(0.6633, -0.9941)	(0.4951, −1.0292)	(−0.2183, −1.0445)	(−1.2752, 0.7396)	(3.3519, −1.8685)	(3.3694, −1.8417)	( 3.3694, −1.8414)	(3.3415, −1.8736)
$\left|f\left({\mathbf{p}}_{\Gamma }\right)\right|$	$7.4·{10}^{-12}$	$1.4·{10}^{-12}$	$1.7·{10}^{-12}$	$6.1·{10}^{-12}$	$4.9·{10}^{-12}$	$1.8·{10}^{-13}$	$5.4·{10}^{-13}$	$8.6·{10}^{-12}$
$\left|〈{\mathbf{V}}_1,{\mathbf{V}}_2〉\right|$	$1.7·{10}^{-13}$	$1.2·{10}^{-12}$	$1.1·{10}^{-12}$	$1.7·{10}^{-12}$	$2.3·{10}^{-13}$	$5.7·{10}^{-14}$	$2.8·{10}^{-13}$	$4.5·{10}^{-12}$

,

Table 2. Test results of Comprehensive Algorithm B for Example 2.

p	(565, 945)	(979, 325)	(375, 405)	(1959, 1377)	(−9356, 8623)	(−816, 798)	(−3987, 1231)	(−4821, 647)
${\mathbf{p}}_{\Gamma }$	(0.2978, −0.9108 )	(1.0773, −0.0164 )	(0.6094, 0.5450)	(0.4780, 0.6961 )	(−1.2055, 1.0209 )	(−1.2035, 1.0230 )	(−1.2288, 0.9794 )	(−1.2341, 0.9554 )
$\left|f\left({\mathbf{p}}_{\Gamma }\right)\right|$	$1.1·{10}^{-14}$	$6.0·{10}^{-13}$	$8.0·{10}^{-18}$	$9.2·{10}^{-15}$	$5.8·{10}^{-15}$	$1.1·{10}^{-14}$	$2.5·{10}^{-15}$	$1.1·{10}^{-14}$
$\left|〈{\mathbf{V}}_1,{\mathbf{V}}_2〉\right|$	0	$1.8·{10}^{-12}$	$2.3·{10}^{-13}$	$4.5·{10}^{-12}$	$6.9·{10}^{-10}$	$2.1·{10}^{-11}$	$1.7·{10}^{-10}$	$1.1·{10}^{-10}$
$\mathbf{p}$	(−7942, −275)	(−598, −98)	(−3709, −1979)	(−2937, −1391)	(9708, −323)	(2608, −1912)	(7347, −931)	(5091, −1921)
${\mathbf{p}}_{\Gamma }$	(−1.2367, 0.8930 )	(−1.1847, 0.4851 )	(−1.0035, −0.1601)	(−1.0225, 0.1224 )	(1.1427, −0.9100)	(1.1224, −0.9787)	(1.1414, −0.9273)	(1.1344,−0.9563)
$\left|f\left({\mathbf{p}}_{\Gamma }\right)\right|$	$2.3·{10}^{-15}$	$1.3·{10}^{-11}$	$9.5·{10}^{-15}$	$1.4·{10}^{-13}$	$4.7·{10}^{-15}$	$3.2·{10}^{-15}$	$4.4·{10}^{-15}$	$1.7·{10}^{-17}$
$\left|〈{\mathbf{V}}_1,{\mathbf{V}}_2〉\right|$	$4.0·{10}^{-11}$	$3.6·{10}^{-12}$	$8.7·{10}^{-11}$	$7.2·{10}^{-12}$	$5.8·{10}^{-11}$	$8.7·{10}^{-11}$	$7.2·{10}^{-12}$	$2.1·{10}^{-10}$

and

Table 3. Test results of Comprehensive Algorithm B for Example 3.

p	(1387, 645)	(1879, 395)	(3075, 205)	(1956, 777)	(−9256, 4603)	(−836, 1798)	(−5987, 1031)	(−4181, 1247)
${\mathbf{p}}_{\Gamma }$	(−1.0629, 0.6023)	(0.4232, −0.8650)	(0.4214, −0.8515)	(0.4276, −0.8795)	(−1.1305, 0.9655 )	(−1.0823, 1.0103 )	(−0.4232, 0.8220)	(−1.1369, 0.9490 )
$\left|f\left({\mathbf{p}}_{\Gamma }\right)\right|$	$1.7·{10}^{-13}$	$1.3·{10}^{-16}$	$8.1·{10}^{-17}$	$1.0·{10}^{-16}$	$5.3·{10}^{-17}$	$1.8·{10}^{-15}$	$6.1·{10}^{-17}$	$1.0·{10}^{-17}$
$\left|〈{\mathbf{V}}_1,{\mathbf{V}}_2〉\right|$	$3.6·{10}^{-12}$	$3.2·{10}^{-12}$	$6.3·{10}^{-13}$	$6.3·{10}^{-12}$	$1.7·{10}^{-10}$	$4.3·{10}^{-11}$	$1.8·{10}^{-12}$	$6.5·{10}^{-11}$
$\mathbf{p}$	(−7342, −1275)	(−5098, −918)	(−3217, −2079)	(−2337, −1251)	(9508, −375)	(6608, −712)	(2347, −931)	(1491, −1321)
${\mathbf{p}}_{\Gamma }$	(−0.4226, 0.8617 )	(−0.4227, 0.8623 )	(−1.0274, 0.5370)	(−1.0471, 0.5706)	(1.2367, −0.9267)	(1.2353, −0.9460)	(1.2255, −0.9888)	(1.2068, −1.0194)
$\left|f\left({\mathbf{p}}_{\Gamma }\right)\right|$	$5.7·{10}^{-17}$	$1.3·{10}^{-16}$	$4.4·{10}^{-14}$	$1.1·{10}^{-13}$	$2.2·{10}^{-15}$	$2.8·{10}^{-15}$	$1.8·{10}^{-15}$	$1.1·{10}^{-15}$
$\left|〈{\mathbf{V}}_1,{\mathbf{V}}_2〉\right|$	$1.0·{10}^{-11}$	$5.5·{10}^{-12}$	$7.3·{10}^{-12}$	$9.1·{10}^{-12}$	$2.2·{10}^{-11}$	$3.6·{10}^{-12}$	$5.4·{10}^{-11}$	$2.2·{10}^{-11}$

, the four symbols $\mathbf{p}$, ${\mathbf{p}}_{\Gamma }$, $\left|f\left({\mathbf{p}}_{\Gamma }\right)\right|$ and $\left|〈{\mathbf{V}}_1,{\mathbf{V}}_2〉\right|$ are the original test point, the orthogonal projection point of the original test point, the deviation degree of the orthogonal projection point on the planar algebraic curve and the absolute value of the inner product of two vectors ${\mathbf{V}}_1$ and ${\mathbf{V}}_2$, respectively, where ${\mathbf{V}}_1$ is $\overrightarrow{{\mathbf{pp}}_{\Gamma }}$ and ${\mathbf{V}}_2$ is the tangent vector $\left[-\frac{\partial f}{\partial y},\frac{\partial f}{\partial x}\right]$ of the orthogonal projection point ${\mathbf{p}}_{\Gamma }$ on the planar algebraic curve $f\left(\mathbf{x}\right)$. Thanks to the suggestions by the reviewers, the fourth quadrant result values of the three tables are implemented by Second Remedial Algorithm in Maple 18 environment.

Example 1

(Reference to [28]). Suppose a planar algebraic curve $f(x,y)=x^6+2x^{5}y-2x^3{y}^2+x^4-y^3+2y^8-4=0$ (See Figure 7). In each of the four quadrants, randomly select four distant test points. We calculate the corresponding orthogonal projection point for each test point via computation by Comprehensive Algorithm B and Second Remedial Algorithm. The specific results are shown in Table 1).

Example 2.

Suppose a planar algebraic curve$f(x,y)=x^{10}+6xy+2y^{18}-2=0$(See Figure 8). In each of the four quadrants, randomly select four distant test points. We calculate the corresponding orthogonal projection point for each test point via computation by Comprehensive Algorithm B and Second Remedial Algorithm. The specific results are shown in Table 2.

Example 3.

Suppose a planar algebraic curve$f(x,y)=x^{10}+6xy+2y^{16}+2=0$(See Figure 9). In each of the four quadrants, randomly select four distant test points. We calculate the corresponding orthogonal projection point for each test point via computation by Comprehensive Algorithm B and Second Remedial Algorithm. The specific results are shown in Table 3.

Remark 5.

Besides all test points of the three examples mentioned above are tested by Comprehensive Algorithm B, we have tested them again with the Second Remedial Algorithm. All the test results are consistent with those of Comprehensive Algorithm B and convergent. In addition, in the region$\left[-3000,3000\right]\times \left[-3000,3000\right]$of each example, we randomly select a large number of test points, the probability of non-convergence is particularly low by Second Remedial Algorithm. Further, we use Second Remedial Algorithm other examples with test points in a very large area, and the probability of non-convergence is also very low. Second Remedial Algorithm is verified to be the best one again in our paper. Of course, the robustness and the efficiency of Second Remedial Algorithm is better than that of the existing algorithms.

5. Conclusions and Future Work

In this paper, we have constructed a Comprehensive Algorithm which is an improved curvature circle algorithm for orthogonal projecting onto planar algebraic curve. Based on an integrated hybrid second-order algorithm [28], the Comprehensive Algorithm (the improved curvature circle algorithm) has also incorporated the curvature circle technique and Newton’s gradient steepest descent method such that it can converge robustly and efficiently no matter how far the test point is from the planar algebraic curve and no matter where the initial iterative point is located. Furthermore, we propose Second Remedial Algorithm based on Comprehensive Algorithm B. In particular, its robustness and efficiency is greatly improved than that of Comprehensive Algorithm B and it achieves our expected result. The numerical examples show that the improved curvature circle algorithm is superior to the existing ones. In future work, we try to refine the idea of Comprehensive Algorithm and Second Remedy Algorithm. And the idea is applied to point orthogonal projecting onto spatial algebraic curve and algebraic surface.