Point Orthogonal Projection onto a Spatial Algebraic Curve

: Point orthogonal projection onto a spatial algebraic curve plays an important role in computer graphics, computer-aided geometric design, etc. We propose an algorithm for point orthogonal projection onto a spatial algebraic curve based on Newton’s steepest gradient descent method and geometric correction method. The purpose of Algorithm 1 in the ﬁrst step of Algorithm 4 is to let the initial iteration point fall on the spatial algebraic curve completely and successfully. On the basis of ensuring that the iteration point fallen on the spatial algebraic curve, the purpose of the intermediate for loop body including Step 2 and Step 3 is to let the iteration point gradually approach the orthogonal projection point (the closest point) such that the distance between them is very small. robust and efﬁcient which it achieves the expected and ideal result.


Introduction
Orthogonal projection is widely used and plays an important role in geometric modeling, computer graphics and computer aided geometric design. Pegna et al. [1] first proposed the concept of orthogonal projection, and discussed the calculation projecting problem of spatial parametric curve orthogonal projection onto a parametric surface and algebraic surface. The so-called orthogonal projection is to find a point on the curve such that the line segment connected to this point and the given point is perpendicular to the tangent line of the curve at this point. There is a close relationship between the orthogonal projection problem and distance projection problem [1]; therefore, the study of the distance projection problem is helpful for the study of the orthogonal projection problem to some extent.
Many scholars have done a lot of work in these two aspects. Hartmann [2] proposed a first-order tangent line method to calculate point orthogonal projection onto parametric curve and surface. For a few cases of non-convergence and as supplement and perfection of the first-order tangent line method [2], Liang et al. [3] and Li et al. [4] presented hybrid second-order method for orthogonal projection onto parametric curve and surface, respectively. Hu et al. [5] proposed a second-order geometric iterative algorithm with curvature information to approximate the orthogonal projection point of the given point on parametric curves and surfaces. Based on their work [5], Li et al. [6] provided an improved method for point orthogonal projection onto general parametric surface such that the efficiency in [6] is greatly improved compared with the existing methods. Ma et al. [7] studied the problem for point orthogonal projection onto NURBS curves and surfaces for which they adopted four-step technique: subdividing curve or surface into curve segments or surface patches, finding out the relationship between the control polygon of curve segment or surface patch and the test point, candidate curve segment or surface patch and approximate projection points being found out by comparison and the final projective point being obtained by comparing the distance between the test point and these approximate projective points. Since the minimum distance between two geometries occurred between a pair of special points, they studied the minimum distance between various specific geometry by using the property and obtained some satisfactory results [8][9][10].
Since the algebraic curve and surface are difficult to control accurately, it is difficult to find the orthogonal projection point on the algebraic curve and surface. At the same time, geometric modeling, computer graphics, computer aided geometric design and many other problems can come down to point orthogonal projection onto curve and surface problems. Therefore, it is of great significance to study point orthogonal projection onto algebraic curves and surfaces.
Up to the present day, the study of point orthogonal projection onto algebraic curves and surfaces is mainly divided into three types: orthogonal projection onto planar algebraic curve, spatial algebraic curve and algebraic surface. Since algebraic curves and algebraic surfaces have no parametric control, it is not easy to solve such problems directly. However, many scholars turned such problems into root-finding problems. There are three types for point orthogonal projection onto planar algebraic curve: local method, global method and compromise method between these two methods. According to the most basic geometric characteristics, the problem of point orthogonal projection onto planar algebraic curve can be understood as solving the equation which the cross product of gradient ∇ f (X) of the planar algebraic curve f (X) at point X and the vector − → PX is zero. The equation can be specifically expressed as the following, ∇ f (X) × (P − X) = 0 (1) In fact, the geometric meaning is that the inner product of the tangent vector of the planar algebraic curve f (X) at point X and the vector − → PX is zero. Equation (1) can be transformed into the corresponding iterative formula of Newton's steepest gradient descent method.
William et al. [11] used composed method including Lagrange multiplier and Newton's steepest gradient descent method to compute the orthogonal projection point on the planar algebraic curve for each test point. The characteristic of the method [11] is that it is fast but local and depends on the initial iterative point.
In addition to the Newton's local iterative method to solve the equation for point orthogonal projection onto the planar algebraic curve, the first global method converting into solving system of nonlinear equations is the Homotopy continuous method [12,13]. They used the most classical Homotopy transformation method, where t is a continuous transformational parameter from 0 to 1, P(X) = 0 and Q(X) = 0 are original system of nonlinear equations to be solved and objective solution system of nonlinear equations, respectively. All isolated solutions of original system of nonlinear equations P(X) = 0 can be obtained by the numerical continuous Homotopy method [12,13], where all the isolated solutions of P(X) = 0 are exactly the same as the target solutions of system of nonlinear equations Q(X) = 0. The robustness of the numerical continuous Homotopy methods [12,13] with global convergence is expounded and verified by [14] and their high time consumption property is summarized in [15]. The major difficulty of Homotopy method is how to find a very satisfactory objective solutions of system of nonlinear equations Q(X) = 0. Turning the projection problem into a resultant problem and then solving the resultants [16][17][18][19] is called the second global resultant method. Using the classical elimination theory, a system of two nonlinear equations with two variables can be transformed into a resultant polynomial with one variable being equivalent to the two original equations. The Sylvester's resultant and Cayley's statement of Bézout's method are the most classical resultant methods [16][17][18][19]. If the degree of the planar algebraic curve is no more than 5, all roots of the nonlinear polynomial equation formed by the resultant methods can be solved. However, if the degree of the planar algebraic curve is more than 5, it is very difficult to directly solve all roots of two-polynomial system with the resultant method.
The third global method is Bézier clipping technique [20][21][22]. In the first step, the nonlinear system of Equation (1) is turned into Bézier form with convex hull property. The rest of the treatment is exactly the same as that of Ma et al. [7] so we do not give a detailed description here. The superiority of global Bézier clipping method is that all roots can be solved, or all orthogonal projection points can be obtained this way. There are two disadvantages of the global clipping method: the first disadvantage is that it needs a lot of time to subdivide, seek, judge and compare; the second disadvantage is that Equation (1) is difficult to directly convert into Bernstein-Bézier form for a small number of planar algebraic curves. Of course, an indirect solution to the difficulty which the idea is come from reviewer's suggestion is that one can use a multivariate solver for this problem. Equation (7) gives a 3×3 algebraic system and its solution gives the stationary (candidate) points.
The first compromise method that is between local and global methods is proposed by Hartmann [2,23] which is composed of the geometric tangent property for computing orthogonal projection point. Continue to execute the iterative formula (2) until the iterative point falls on the planar algebra curve. Taking the iterative point fallen on the planar algebraic curve as being the initial iterative point, compute the foot point once by the iterative formula (4), Repeat the above two steps until the foot point Q and the orthogonal projection point completely overlap. Unfortunately, for a small part of the planar algebraic curve, it is invalid to approach the orthogonal projective point with progressive tangent foot point. Based on Hartmann's work in [2,23], Li et al. [24] proposed an improved method for calculating footpoint which the effect is good if the test point is not far away, but the effect is not ideal if the test point is particularly far away. In addition, on the important roles and applications of the algebraic curves and algebraic surfaces, some introduction and explanation are presented in these papers [2,9,10,23,24]. Here we will not explain and elaborate upon them in detail. The second compromise method is proposed by Redding [25] who used the osculating circle technique to accomplish point orthogonal projection onto the planar algebraic curve. The osculating circle technique mainly includes three steps: (1) Compute curvature at point on the planar algebraic curve and the corresponding radius and center of the curvature circle determined by curvature. (2) Get the footpoint Q intersected with the line segment formed by the test point and the center of the curvature circle and curvature circle. (3) Take the footpoint Q approximately as the point fallen on the planar algebraic curve. Repeatedly execute the above three steps until the foorpoint Q and the orthogonal projection point P Γ completely coincide. Since there is a certain deviation in the third step of the osculating circle technique for footpoint Q approximating to the orthogonal projection point on the planar algebraic curve and the planar algebraic curve has no parametric control like parametric curve, the robustness of the osculating circle technique cannot be achieved very well. Greatly inspired by these works [2,23,25], Wu et al. [26] proposed an improved curvature circle algorithm for orthogonal projection onto a planar algebraic curve. Not only in robustness but also in efficiency, the Second Remedial Algorithm in [26] is very satisfactory and ideal.
The third compromise method is the circle shrinking method [27]. Repeatedly iterate the iterative formula (2) such that the iterative point can fall on the planar algebraic curve as much as possible. Draw a circle where the center and radius are the test point P and P − P I , respectively. Identify a point P + on the circle by means of the mean value theorem and find the intersection named the current iterative point P I which is the intersection between the line segment PP + and the circle. Repeatedly run the above behavior until the distance between the current point P I and the previous point P I is almost zero. The circle shrinking technique [27] consumes more time to find point P + each time and it is difficult to directly compute the intersection P I intersected of line segment PP + and the planar algebraic curve if the degree of the planar algebraic curve is more than 5.
The fourth compromise method with circle double-and-bisect algorithm is proposed by Hu et al. [28]. Starting with a very small circle with center being test point P and radius r 1 being arbitrarily small, draw a new circle with the same center being test point P and radius r 2 = 2r 1 (After that, the center of all circles is the test point P). If the new circle does not intersect the planar algebraic curve, draw a new circle with radius twice of r 2 . Keep repeating the above action, until the latest circle can intersect with the planar algebraic curve. The previous circle and the latest circle are tagged as the interior circle and the exterior circle, respectively. Furthermore, implement the rest of the process by adopting the bisecting technology. Continue to draw a new circle with new radius r = (r 1 + r 2 )/2. If the current circle with radius r intersects with the planar algebraic curve, substitute r for r 2 , else for r 1 . Repeatedly execute the above progress until the difference between the current radius and the previous radius is almost 0. However, with this method it is not easy to judge whether the exterior circle intersects the planar algebraic curve or not [28]. Moreover, the circle double-and-bisect algorithm [28] consumes more time to seek the intersection between the exterior circle and the planar algebraic curve.
On the problem of point orthogonal projection onto curve and surface, there exists a classical method which is transforming the problem of point orthogonal projection onto curve and surface into a specific solver problem. Bartoň M. [29] proposed two blending schemes solvers for solving the problem of point orthogonal projection onto curve and surface. As a system of nonlinear equations, throwing away no-root domain can be realized by a simple linear combination of functions and then determine all control points for its Bernstein-Bézier basis having the same sign, which must be in accord with the seeking function. It can continuously obtain these types of functions to get rid of the no-root domain through the continuous subdivision process. The efficiency of the geometric-based algorithm [29] is higher than that of the existing subdivision mode based GPU. Therefore, two blending schemes in [29] can efficiently reduce the number of subdivisions by using continuous eliminating no-root domains. Result from the consequence in [29], van Sosin and Elber [30] constructed a variety of complex piecewise polynomial systems having zero or inequality constraints in zero-dimensional or one-dimensional solution spaces. Finding out all the root solutions is the advantage of these methods [29,30]; however, more computation consume-time and the need for many subdivision steps are their disadvantages.

Newton's Steepest Gradient Descent Method
Let us elaborate on the general idea. Let P be a test point. There is a spatial algebraic curve C which can be written as, where X is equal to (x, y, z), and f 1 (X) = 0 and f 2 (X) = 0 are two algebraic surfaces, respectively. Their intersection set constitutes a spatial algebra curve C. The aim of this paper is to seek a spatial point P Γ on the spatial algebraic curve C to satisfy the basic relationship (See Figure 1), On the other hand, the vector − → PX determined by the test point P and point X on the spatial algebraic curve C is orthogonal to the tangent vector T, where T = ∇ f 1 (X) × ∇ f 2 (X). The corresponding equation formed by the inner product of the vector − → PX and vector T can be expressed as, The above two questions can be summed up as, is Hamiltonian operator, symbol × is cross product and symbol , is inner product. We take s as the arc length parameter of the spatial algebraic curve C and Let be tangent vector along the spatial algebraic curve C. There are many ways to construct the Newton's steepest descent methods of the spatial algebraic curve with Equation (5). After many tests and comparisons, we find that the following model of Newton's steepest descent method with Equation (9) is the best choice for the Newton's steepest descent methods of the spatial algebraic curve.
Repeatedly run Equation (9), called Newton's steepest gradient descent method, until the iterative termination criteria are satisfied | f 1 (X n+1 )| < ε and | f 2 (X n+1 )| < ε, where the initial iterative point is X 0 = P − (0.1, 0.1, 0.1) and the iterative point X n+1 fallen on the spatial algebraic curve C is called P I . The Newton's steepest gradient descent method (Algorithm 1) can be specifically described as (See

Algorithm 1:
The Newton's steepest gradient descent method. Input: The test point P and the spatial algebraic curve C. Output: The iterative point P I fallen on the spatial algebraic curve C. Description: Step 1: Update X n+1 according to the iterative Equation (9); Step 2:

Remark 1.
Starting from the initial iterative point X 0 , the iterative point of the first equation of the Newton's steepest gradient descent method with Equation (9) causes the iterative point to be closer to the first algebraic surface f 1 (X). Then through the iterating of the second equation with Equation (9), the iterative point of the second equation with Equation (9) causes the iterative point to be closer to the second algebraic surface f 2 (X). In this way, Equation (9) is iterated repeatedly such that the distance between the finally iterative point and every algebraic surface is almost zero. Namely, through alternating iteration of two sub-equations with Equation (9), the distance between the finally iterative point and the spatial algebraic curve determined by the intersection set of two algebraic surfaces is almost zero. The Newton's steepest gradient descent method (9) has two advantages. The first advantage of the Newton's steepest gradient descent method (9) is to promote the iterative point to fall on the spatial algebraic curve C as much as possible. The second advantage of the Newton's steepest gradient descent method (9) is that the iteration point fallen on the spatial algebraic curve is closer to the orthogonal projection point P Γ , it is very convenient for implementation of the subsequent sub-algorithms.

Computing Foot-Point Q
Although the iterative point P I falls on the spatial algebraic curve, there is still a certain small distance between the iterative point P I and the orthogonal projection point P Γ . There is still a small certain gap between this and our ideal and goal. Our idea is to find a way to gradually reduce the distance between the iteration point P I and the orthogonal projection point P Γ to zero. In order to promote the iteration point P I and the orthogonal projection point P Γ closer, we use the tangent line perpendicular foot method to promote the iteration point P I and the orthogonal projection point P Γ being closer. The basic idea of this method can be realized as the following form. Draw a tangent line L of the spatial algebraic curve C at the iterative point P I fallen the the spatial algebraic curve C. This tangent line L can be written as, where t is a tangent line parameter to be determined, and T is the tangent vector of the spatial algebraic curve C at the iterative point P I . Through the test point P, let us construct a plane such that the tangent line L is perpendicular to the plane. From primary geometric intuitiveness, it is not difficult to find that the intersection Q 0 named the foot point intersected with the tangent line L and the plane can be expressed as where t 0 = P − P I , T / T, T . In this way, the position of the foot point Q 0 is between the iteration point P I and the orthogonal projection P Γ . At this time, taking foot point Q 0 as the initial iterative point of Algorithm 1, we obtain the new iterative point P I fallen on the spatial algebraic curve C which the new iterative point P I is called the current iterative point P I . According to our observations and the basic geometry characteristic, we can also choose the foot point closer to the spatial algebraic curve.
Let us choose a point Q on the line segment By solving Equation (12), foot-point Q can be specifically expressed as The computation of the foot-point Q can be achieved via Algorithm 2 (See Figure 3).

Algorithm 2:
Computing foot-point Q. Input: The test point P, the spatial algebraic curve C and the iterative point P I on the spatial algebraic curve C. Output: The foot-point Q. Description: Step 1: Compute the foot-point Q by the formula (13).

Accelerating and Correcting Method
By Algorithm 1, the initial iteration point X 0 can fall on the spatial algebraic curve C and a final iterative point P I is generated. Then the foot-point Q is obtained by Algorithm 2 derived from the final iterative point P I . So the foot-point Q is taken as the initial iterative point of Algorithm 1, and a new current final iterative point P I fallen on the spatial algebraic curve C is generated. Run Algorithm 1 and Algorithm 2 alternately and repeatedly, the current iterative point P I is getting closer to the orthogonal projection point P Γ . However, the rate of the iterative point P I approaching the orthogonal projection point P Γ becomes slower and slower. And orthogonalization is hard to achieve. In other words, the first three sub-formulas of Formula (8) are difficult to satisfy which this process is not up to the ideal.
In order to overcome these two disadvantages, we add accelerating and correcting method after Algorithm 2 such that the first three sub-formulas of formula (8) can be satisfied as soon as possible.
Since T = ∇ f 1 × ∇ f 2 is the tangent vector of the spatial algebraic curve C at point X, so the relationship between the tangent vector T = ∇ f 1 × ∇ f 2 and the spatial algebraic curve C can be written as, Furthermore, let c = d 2 x dt 2 , d 2 y dt 2 , d 2 z dt 2 be the curvature vector of the spatial algebraic curve C at the point X. Taking the first-order derivative of Equation (14) with respect to arc length parameter s, get the corresponding equation about curvature vector, where ∂x∂y , etc. Solving Equation (15), the specific curvature vector c can be expressed as where ∂x , x s = dx ds , etc. For the case of the spatial algebraic curve, the iterative point P I may deviate from the spatial algebraic curve C. That is, the iterative point P I is not on two algebraic surfaces, f 1 (P I ) = err1 = 0, f 2 (P I ) = err2 = 0. The basic idea of error correcting can be described as follows. When the iterative point P I deviates from two algebraic surfaces, use correction vector δV = [δx, δy, δz] to correct iteration point P I such that | f 1 (P I | < Le and | f 2 (P I )| < Le (Le is a particularly small value of the allowable error range). The corrected iteration point P I is guaranteed to be on the spatial algebraic curve and the deviations between the corrected iterative point P I and two algebraic surfaces are almost zero. ∆V is on the plane formed by the tangent vector T and curvature vector c derived from the iteration point P I before correction and is perpendicular to the vector T × c. Let δV be on the plane spanned by the vector ∆V and T × c, and δV points to two algebraic surfaces. From the above description, we can get the equation about the correction vector δV, From Equation (17), it is not difficult to get the expression of the correction vector δV, where err1 = f 1 (P I ), err2 = f 2 (P I ), at this moment, the iteration point P I is the iteration point before correction, ∆V is the difference vector of the corrected iteration point P I subtracting the iteration point P I before correction. After a lot of testing, we have obtained a preliminary result. When the correction vector δV is added to the iteration point P I after iterating, it can promote the robustness and high efficiency of convergence from iteration point P I to orthogonal projection point P Γ . However, the rate which the iteration point P I converges to the orthogonal projection point P Γ is still very slow, at the same time, the rate of the orthogonalization of the tangent vector of the spatial algebraic curve at the iteration point P I and the vector connected to the test point P and the iteration point P I is also very slow. There is still a big gap between this situation and our expected goal. In order to accelerate the convergence of the iteration point to the orthogonal projection point P Γ and speed up the orthogonalization, we add three iterative formulas related to Newton's steepest gradient descent method in front of the correction vector δV. Therefore, the corresponding iterative formula for accelerating and orthogonalizing is as follows, where err1 = f 1 (Z n ), err2 = f 2 (Z n ), ∇ f 1 , ∇ f 2 , T, c assigned by Z n , ∆V = U n − Z n .

Remark 2.
We present a complete geometric interpretation for Equation (19). The first two formulas in Equation (19) that are not related to the test point P are Newton's steepest gradient descent method. The ultimate goal is to make the initial iteration point fall on the spatial algebraic curve C. Please refer to Remark 1 for the specific explanation of geometric meaning. The third formula of Equation (19) is the Newton's steepest gradient descent method related to the test point P. This iterative formula is the iterative formula of Newton's steepest gradient descent method formed by the direct deformation of the formula (7), which the formula (7) is the inner product of the vector − → PX formed by the test point P and the point X on the spatial algebra curve C with the tangent vector T at that point X is zero. Therefore, the third iteration formula of Equation (19) is mainly used to accelerate the orthogonalization of the vector − → PX and the tangent vector T. The last formula of Equation (19) is to use the geometric correction vector δV to promote the iteration point fall on the spatial algebraic curve C as much as possible once again. Therefore, Equation (19) plays a dual roles of accelerating the iteration point to fall on the spatial algebraic curve C and orthogonalization of the vector T and the vector − − → PP Γ . Then accelerating and correcting the method with Equation (19) makes the first three formulas of the formula (8) achieves very satisfactory result. 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 accelerating and correcting method with Equation (19).
Based on the above illustration, we obtain the following Algorithm 3 (See Figure 4).

Algorithm 3:
Accelerating and geometric correcting method for orthogonal projection onto a spatial algebraic curve.
Input: The iterative point P I fallen on the spatial algebraic curve C and the spatial algebraic curve C. Output: The final orthogonal projection point P Γ fallen on the spatial algebraic curve C. Description: Step 1: Step 2:  Based on the above description, we give the complete algorithm for point orthogonal projection onto a spatial algebraic curve(Algorithm 4) (See Figure 5).

Algorithm 4:
A complete algorithm for point orthogonal projection onto a spatial algebraic curve.
Input: Test point P and the spatial algebraic curve C. Output: The final orthogonal projection point P Γ of the test point P. Description: Step 1: Starting from the neighbor point of the test point P, calculate the iterative point P I fallen on the spatial algebraic C via Algorithm 1. for(i = 0;i < 5;i + +) { Step 2: Compute the foot-point Q via Algorithm 2.
Step 3: Project foot-point Q onto the spatial algebraic curve C via Algorithm 1, then get new point P I fallen on the spatial algebraic curve C. } Step 4: The iterative point P I fallen on the spatial algebraic curve C being as the initial iterative point, calculate the finally orthogonal projection point P Γ fallen on the spatial algebraic C via Algorithm 3.

Remark 3.
In this remark, we present a geometric interpretation for Algorithm 4. The purpose of the first step is to promote the initial iterative point completely iterate to the spatial algebraic curve. In Step 2, the foot-point Q is derived from the iterative point P I fallen on the spatial algebraic curve of the first step. Now, the position of foot-point Q is between the iteration point P I and the orthogonal projection point P. In other words, the foot-point Q is closer to the orthogonal projection point than the test point P. Furthermore, by using Algorithm 1 of the third step, the foot-point Q can be completely iterated to the spatial algebra curve, and the point at which this iteration point falls on the spatial algebraic curve is called the current iteration point and the last iteration point to fall on the spatial algebra curve is named the previous iteration point. It is obvious that the current iteration point is closer to the orthogonal projection point than the previous iterative point. After joint running Step 2 and Step 3 five times in this way, the distance between the latest current iterative point and the orthogonal projection point is very small. If we keep running in this way, the rate of the iteration point getting closer to the orthogonal projection point becomes slower and slower. Moreover, orthogonalization of the tangent vector T and the vector − − → PP Γ is also difficult to realize. However, Algorithm 3 in the fourth step can make up for these two deficiencies. At this time, running Algorithm 3 plays double important roles for accelerating convergence and orthogonalization. Thus, Algorithm 4 is a robust and efficient algorithm which it achieves a very satisfactory result for point orthogonal projection onto the spatial algebraic curves. In addition, the five cycles of the intermediate for loop body are an empirical value. When the distance between the test point and the spatial algebra curve is not very far, the intermediate for loop body only needs 5 cycles. However, when the distance between the test point and the spatial algebraic curve is very large, the number of cycles of the for loop body needs to be increased. The purpose is to promote the iterative point after the loop of for loop body closer to the orthogonal projection point. Therefore, Algorithm 3 of the fourth step can conver the orthogonal projection point P Γ successfully and completely.
We have done a lot of tests for Algorithm 4. Whether the test point P is close to or far from the spatial algebra curve C, the test point can be orthogonally projected to the spatial algebraic curve in a very robust and efficient way. We have another surprise and wonderful discovery. When the test point is not far from the spatial algebraic curve, Algorithm 4 can be expressed in simplified version form, where the for loop body in the middle can be omitted and the fourth formula with geometric correction method of Equation (19) can also be deleted. So Equation (19) naturally becomes Equation (20).
The simplification version of Algorithm 4 can be described as the following Algorithm 5, where Algorithm 5 maintains the same robustness as Algorithm 4, and the time running efficiency of Algorithm 5 is much better than that of Algorithm 4.

Algorithm 5:
A simplified version algorithm for point orthogonal projection onto a spatial algebraic curve.
Input: Test point P and the spatial algebraic curve C. Output: Orthogonal projection point P Γ of the test point P which orthogonally projects the test point P onto the spatial algebraic curve C. Description: Step1: Calculate the iterative point P I fallen on the the spatial algebraic curve C via Algorithm 1.

Remark 4.
In the actual algorithms implementation process, we take two optimized techniques. Firstly, if the test point P is very far away from the spatial algebraic curve C, the initial iterative X 0 of Algorithm 1 is a very small proportion number of the test point P which ensures that the initial iterative point X 0 converges to the spatial algebraic curve C very robustly. However, the initial iterative point X 0 of Algorithm 1 in Step 3 of Algorithm 4 remains unchanged(X 0 =Q ). Secondly, due to the spatial algebraic curve being composed of the intersection of two algebraic surfaces, the singularity of any algebraic surface may lead to the singularity of spatial algebraic curve. That is to say, the singularity of a spatial algebraic curve will cause the denominators of all formulas and iterative formulas to be zero. In order to avoid the degradation of all iterative formulas with zero denominators, we add a small perturbation number ε to the denominator of each iteration formula such that every iterative formula will not degenerate. Thus, every algorithm can run and iterate normally. Furthermore, no matter how far away the test point is from the spatial algebraic curve and whether there are singular points or not in the spatial algebra curve, Algorithm 4 can run very robustly and efficiently which it meets our expectations and ideals.

Finding out Initial Iterative Candidate Points
This subsection comes from the reviewer's previous suggestion. In the field of practical computer graphics, computer-aided geometric design and other applications, it is necessary to solve not only a single orthogonal projection point but also all orthogonal projection points as much as possible. Our preliminary idea is to outline the spatial algebraic curve which the specific expression is presented by the formula (5). We try to find out several 3D bounding boxes within the specified region of the spatial algebraic curve such that every curve segment of the spatial algebraic curve is contained in one 3D bounding box. Randomly select a point in each 3D bounding box as the initial iteration point of Algorithm 4 or Algorithm 5. In this way, we can yield as many orthogonal projection points as possible by using Algorithm 4 or Algorithm 5. Then, through computing the distance between the test point P and each orthogonal projection point and by comparing with each distance one by one, the orthogonal projection point corresponding to the shortest distance is what we are looking for. However, it is very difficult to express the spatial algebraic curve directly. On the other hand, the geometric essence of point orthogonal projection onto the spatial algebraic curve is to find a point X on the spatial algebraic curve such that the vector − → PX and the tangent vector T(T = ∇ f 1 (X) × ∇ f 2 (X)) on the spatial algebraic curve at point X are perpendicular to each other. Namely, the vector − → PX is orthogonal to the tangent vector T. The corresponding equation formed by the inner product of the vector − → PX and the vector T is expressed as the form of Equation (7). The geometric meaning of Equation (7) is algebraic surface. So, outlining the spatial algebraic curve with Equation (5) is transforming into outlining the algebraic surface with Equation (7). We are gong to find several 3D bounding boxes within the specified region of the algebraic surface with Equation (7) such that every surface patch of the algebraic surface is contained in one 3D bounding box. Randomly select a point in each 3D bounding box as the initial iteration point of Algorithm 4 or Algorithm 5. In this way, we can yield as many orthogonal projection points as possible by using Algorithm 4 or Algorithm 5. Then through computing the distance between the test point P and each orthogonal projection point and by comparing with each distance one by one, the orthogonal projection point (the closest point) corresponding to the shortest distance is what we are looking for.Therefore, the orthogonal projection point (the closest point) satisfied the shortest distance completely conforms to every formula of Equation (8).
Let us suppose that the region Ω of the spatial algebraic curve with Equation (5) is . We may decide that the region of the algebraic surface is the same as that of the spatial algebraic curve. So the region of the algebraic surface is also Ω. We adopt adaptive affine arithmetic method [31][32][33] to identify a series of 3D bounding boxes such that every algebraic surface patch is contained in every 3D bounding box. We orthogonally project the algebraic surface F(X) = 0 onto y − z plane, the x − z plane and the x − y plane at the point (x 0 , y 0 , z 0 ), respectively, and we obtain the planar algebraic curves F(x 0 , y, z) = 0, F(x, y 0 , z) = 0 and F(x, y, z0) = 0. respectively. For the sake of convenience, the planar algebraic curves F(x 0 , y, z) = 0, F(x, y 0 , z) = 0 and F(x, y, z0) = 0 can be written as F 1 (y, z) = 0, F 2 (x, z) = 0, F 3 (x, y) = 0, respectively. We set up a crucial judgment function for the planar algebraic curve F 1 (y, z) = 0, Similarly, we also set up two other crucial judgment functions for the planar algebraic curve F 2 (x, z) = 0 and the planar algebraic curve F 3 (x, y) = 0, and where F 1y = ∂F 1 (y,z) The each substitution value of variable X for nine partial derivative functions F 1y , F 1z , F 2x , F 2z , F 3x , F 3y , F 1yz , F 2xz , F 3xy is (x 0 , y 0 , z 0 ). Afterwards, we will call this point (x 0 , y 0 , z 0 ) as the center point of the 3D bounding box. Combining with Equation (21), Equation (22) and Equation (23), we get the crucial judgment function Ψ with Equation (24), We specify a threshold for the crucial judgment function Ψ with Equation (24). The adaptive approach method which we have adopted for solving a series of 3D bounding boxes of the algebraic surface has absorbed the idea of the papers [31][32][33], which is for solving a series of 2D bounding boxes of the planar algebraic curve. The detailed explanation of Equation (24) is exactly similar to the explanation in [34]. On affine arithmetic, we mainly absorbed the ideas of this paper. As for the more complex topological structure of the algebraic surface, we have to absorb the idea of the paper [35] to solve a series of 3D bounding boxes of the algebraic surface. The specific description for solving a series of 3D bounding boxes of the algebraic surface can be expressed as Algorithm 6. Step 1: Subdivide this 3D bounding box into 8 3D sub-bounding boxes by dividing 2 on each axis.
Step 2: Evaluate the value Ψ of each 3D sub-bounding box according to Equation (24).

Convergence Analysis
This section proves that Algorithm 4 does not depend on the initial iteration point.

Theorem 1. Algorithm 4 is independent of the initial iterative point.
Proof. According to Remark 4, if the test point P is very far away from the spatial algebraic curve C, the initial iterative X 0 of Algorithm 1 is a very small proportion number of the test point P. That is to say, as long as the distance between the initial iteration point X 0 and the spatial algebraic curve is relatively near, the initial iteration point can converge to the spatial algebraic curve by using Algorithm 1, or the initial iteration point can successfully and completely iterate and fall on the spatial algebraic curve. Based on the iterative result of Algorithm 1, the general idea of the intermediate for loop body is to ensure that the iteration point falls on the spatial algebraic curve under the premise such that the iterative point P I can gradually move to a very close position of the orthogonal projection point. According to Remark 2, Algorithm 3 with Equation (19) plays a dual roles of accelerating the iteration point to fall on the spatial algebraic curve C and orthogonalization of the vector T and the vector −→ PP Γ . Then accelerating and correcting method with Equation (19) makes the first three formulas of the formula (8) achieve a very satisfactory result. 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 accelerating and correcting the method with Equation (19). Therefore, for the initial iteration point X 0 at any position, by using Algorithm 4, it can be finally iterated to spatial algebraic curve such that the finally iterative point fallen on the spatial algebraic curve is the orthogonal projection point. Namely, the final iteration point X n+1 that falls on the spatial algebraic curve satisfies the first three terms of the Formula (8). Specially, the test point P and the final iteration point X n+1 satisfy the relationship −−−→ PX n+1 , T = 0 or −→ PP Γ , T = 0. It is possible that the orthogonal projection points obtained by Algorithm 4 may be different for different initial iteration points. Once the initial iteration point is determined, the orthogonal projection point is also uniquely determined accordingly. It can converge to the orthogonal projection point robustly and efficiently. It is demonstrated that Algorithm 4 is independent of the initial iterative point.
In a similar way to the proof of Theorem 1, we can state the following result. Theorem 2. Algorithm 5 is independent of the initial iterative point.

Numerical Results
In this section, we present two examples to explain that Algorithm 4 is robust and efficient. All calculation were implemented by using Maple 18 environment with ε = 10 −20 . In Tables 1 and 2, the five symbols P, P Γ , | f 1 (P Γ )|, | f 2 (P Γ )|, |F(P Γ )| are the test point, the orthogonal projection point of the test point, the absolute value of the deviation degree of the orthogonal projection point on the first algebraic surface, the absolute value of the deviation degree of the orthogonal projection point on the second algebraic surface and the absolute value of the expression of Equation (7), respectively.
In this area [−300, 300]×[−300, 300]×[−300, 300]. We randomly select a large number of test points, the probability of convergence can achieve 100% by Algorithm 4 which is verified to be very high robustness and efficiency for remote test points. In addition, in each of the eight quadrants, we randomly select two different test points. We calculate the corresponding orthogonal projection point, the absolute value of the deviation degree of the orthogonal projection point on the first algebraic surface, the absolute value of the deviation degree of the orthogonal projection point on the second algebraic surface and the absolute value of the expression of Equation (7) for each test point via Algorithm 4. The specific results are shown in Table 1.
Example 2. Assume that there is a spatial algebraic curve f 1 (x, y, z) = x 2 9 + y 2 16 + z 2 9 + 4xy − 1 f 2 (x, y, z) = x 2 + y 2 + x − 4z (See Figure 7). In this area [−300, 300]×[-300, 300]×[-300, 300]. We randomly select a large number of test points, the probability of convergence can achieve 100% by Algorithm 4 which is verified to be very high robustness and efficiency for remote test points. In addition, in each of the eight quadrants, randomly select two different test points. We calculate the corresponding orthogonal projection point, the absolute value of the deviation degree of the orthogonal projection point on the first algebraic surface, the absolute value of the deviation degree of the orthogonal projection point on the second algebraic surface and the absolute value of the expression of Equation (7) for each test point via Algorithm 4. The specific results are shown in Table 2.

Remark 5.
If the test point is farther away from the spatial algebraic curve or if the spatial algebraic curve is smooth or if there are fewer singularities in the spatial algebraic curve or if the curvature of the spatial algebraic curve is locally non-high curvature, the algorithm provided by us has certain advantages. However, if the spatial algebraic curve has one of the following cases, any case will hinder and affect the normal operation of Algorithm 4. The robustness and the efficiency of Algorithm 4 are greatly reduced.
(1) The test point is particularly farther away from the spatial algebraic curve.
(2) The spatial algebraic curve is not smooth.
(3) There are many singularities in the spatial algebraic curve.
(4) The curvature of the spatial algebraic curve is locally high.
(5) There are many symmetries in spatial algebraic curves. (6) The spatial algebraic curve are extremely or particularly oscillatory. For example, there is a planar curve f (x) = cos 4 ( 1 x ) + sin 19 (x) + sin 8 ( 1 x 2 ), x ∈ [−4, 4] (See Figure 8). From the graphic display results, this curve is particularly oscillatory. If the initial iteration point is slightly further away from the curve, it is difficult to implement by using Algorithm 1 because the inherent essence of the iterative formula (9) for Algorithm 1 is exactly the same as the classical Newton's iterative method. In view of these non-convergent cases, we are going to take three approaches. The first approach is transforming the test point orthogonal projection onto the spatial algebraic curve problem into the root-finding problem [29,30,36,37]. The second approach is to explore a strategy such that the initial iteration point should be as close as possible to the spatial algebraic curve which the goal is to minimize the sensitivity of the initial iteration point. Absorb the very robust methods in the each part of the algorithm that we already have. Then optimize them to the maximum. The third approach is to explore a one or more brand new algorithms such that every algorithm could solve all or part of the above difficult situations. Of course, it is necessary to study the singularity, non-smoothness, local high curvature and oscillation of spatial algebraic curve. Then find out how they relate to each other. Finally, some ideal algorithms that can meet the requirements of robustness and high efficiency are explored based on their intrinsic relevance.

Conclusions and Future Work
In this paper, we have developed an algorithm for point orthogonal projection onto a spatial algebraic curve. Algorithm 4 mainly includes three core technologies: Newton's steepest gradient descent method, computing foot-point method and geometric correction method. Algorithm 1 impels the iterative point to fall on the spatial algebraic curve. The middle for loop body is to let the iteration point move to the position closed to the orthogonal projection point. The algorithm 3 of the fourth step plays an important double accelerating and orthogonalization role. Numerical examples show that our algorithm is robust and efficient and it achieves our expected result. In future work, we will try to study the singularity, non-smoothness, local high curvature and oscillation of spatial algebraic curve. Then we will find out how they relate to each other. Finally, some ideal algorithms that can meet the requirements of robustness and high efficiency are explored based on their intrinsic relevance. We also try to extract and purify the idea of Algorithm 4. Moreover, the idea is applied to point orthogonal projecting onto algebraic surface.