Abstract
This article presents a method for multipoint inversion and multiray surface intersection problem on the parametric surface. By combining tracing along the surface and classical Newton iteration, it can solve point inversion and ray-surface intersection issues concerning a large number of points or rays in a stable and high-speed way. What is more, the computation result can approximate to exact solutions with arbitrary precision because of the self-correction of Newton-Raphson iteration. The main ideas are adopting a scheme tracing along the surface to obtain a good initial point, which is close to the desired point with any prescribed precision, and conducting Newton iteration process with the point as a start point to compute desired parameters. The new method enhances greatly iterative convergence rate compared with traditional Newton’s iteration related methods. In addition, it has a better performance than traditional methods, especially in dealing with multipoint inversion and multiray surface intersection problems. The result shows that the new method is superior to them in both speed and stability and can be widely applied to industrial and research field related to CAD and CG.
1. Introduction
Point inversion and ray-surface intersection are significant problems, which are to compute the parametric coordinates value of a point by its Cartesian coordinates or ray tangent vectors. Many computational geometry problems, computer vision problems, and computer graphics applications are closely related to surface manipulation. For a parametric surface, the most ordinary manipulation might be to find the parameters of points on the surface in a fast and stable way. Several efforts have been made towards point inversion issue of polynomial parametric surface [1–5]. For approximate methods, Hoschek et al. [6] and Hartmann et al. [7] proposed the first-order geometric iteration method. Hu and Wallner [8] gave out a second-order iteration by computing a normal curvature circle and a line segment intersection point. Liu et al. [9] enhanced the stability and speeded up the convergence by changing the normal curvature circle into a torus patch. Some classic methods are provided by Piegl and Tiller [10] by decomposing a NURBS surface into quadrilaterals, projecting the test point onto the closest quadrilateral, and then recovering the parameter from the closest quadrilateral. Ma and Hewitt [11] put forward an algorithm to obtain a good initial value for the Newton-Raphson method. Xu [5] developed an applicable method by subdividing the parametric domain and checking the signs of all coefficients of a Bernstein polynomial. For NURBS curve and surface, Johnson and Cohen [12] concentrated on the polygonal model for finding the minimum distance between two geometric objects. Chin [13] and Edelsbruner [14] found an algorithm with complexity for finding the minimum distance between two convex polygons. For ray-surface intersection, Taezoon et al. [15] provided a second order geometric iteration which provides better approximation to the local shape of the surface. Wang et al. [16] provided a method combining Bezier clipping and Newton’s method. Wang et al. [17] presented an algorithm for enhancing the performance of both numerical and subdivision methods. Saini [18] gave an application for inverse geometric reconstruction of conics in 3D space by using a ray-surface intersection. For parametric model, Mortensen [19] gave a numerical approach for this problem through finding a vector from the point on the curve that is perpendicular to the tangent vector. Song et al. [20] constructed a biarc that locally approximates a segment on the original curve starting from the current projective point. Ko [21] reviewed methods for computing orthogonal projection of points onto curves and surfaces, which are given on implicit or parametric form or as points cloud. Through parameterizing the progenitor curve, a novel curve projection scheme is proposed by Xu [22] for computing the orthogonal projection. Limaiem [23] presented another approach to find the minimum distance to convex parametric curves and surfaces. Lin [24] provided the approach for finding the minimum distance for concave surfaces. For further reading, papers [25, 26] presented several algorithms for computing the minimum distance between a point and a surface. With the Newton-Raphson iteration initial value Piegl and Tiller [10] provided another method to solve the point projection problem for NURBS surface.
Wang et al. [27] developed several methods for point inversion and ray surface intersection issues by constructing and projecting a line segment and tracing the parameters along the projected curve or intersection curve from an arbitrarily prescribed initial point. To some extent, the approach avoids fully and truly the sensitivity to the choice of the initial value in those Newton-Raphson iteration related methods. Young-Taek Oh et al. [28] used an efficient culling technique to reduce redundant curves and surfaces and proposed an efficient method to project a point to its nearest point on a set of freedom surfaces and curves. Chen et al. [29] revised a rational cubic clipping approach to get two bounding cubes within time, which could get a faster convergence rate.
So far, these existing methods are good at solving point inversion and ray-surface intersection problems only of a single point or several points. However, many engineering issues such as composite material fibre placement and nonplane optical imaging often need to deal with point inversion or ray-surface intersection issues of a large number points simultaneously. Also it needs high computational efficiency (high precision and high success rate). There are two demerits in using traditional Newton related methods to handle the multipoint inversion and multiray surface intersection problems. One is totally time-consumption in finding initial value usually through subdivision for every point. Another is the sensitivity of Newton’s method to an initial guess. As a matter of fact, when sensitivity occurs the iteration process fails to converge. Thus user must restart the whole computational process. So the two disadvantages lead to huge computational costs. In practice, even 5% failure rates of iteration will result in the failure of a path planning process for robotic fibre placement. Therefore we need to develop new method towards the multipoint inversion problems.
Enlightened by the method [27], which can trace solution with any prescribed initial point and Newton iteration’s fast computation, we try to resolve the new problem by combining the two existing methods and improving them. Our main contributions are as follows: overcome low efficiency of traditional methods when dealing with multipoint inversion issue; reduce and even eliminate the sensitivity of traditional Newton’s methods to an initial guess; propose a new algorithm good at multipoint inversion and multiray surface intersection problems.
2. Fundamentals of Algorithm
Here algorithm foundations are introduced and deduced, which are the indispensable parts of the method. Actually the new algorithm of multipoints inversion and multiray surface intersection can be divided into two steps. First, we improve the geometric method to obtain a good initial value. The second important step is to enhance the accuracy of solving the differential equations through the classic Newton-Raphson iteration.
2.1. Multipoint Inversion
Let us consider a parametric surface (a mapping from ) with at least 2-order derivative existing. For any there are and . The problem is as follows: is a point on the surface , with its Cartesian coordinates known. See Figure 1. We want to compute its parametric coordinates . Firstly take a point arbitrarily on the given surface . Construct a line segment with and as follows:Calculate a normal vector of the known point on the surface . Then define a vector as , which represents the project direction of line segment . Let the projected curve of onto the surface be , which can be characterized as followsComputing the derivative on both sides of Equation (2) with respect to , it follows thatTaking the dot product of both sides of (3) with the vector and separately yieldsThus, if the mixed product is not equal to 0, by solving the matrix equation (4), we have and as follows:Equation (5) is the govern equation of multipoint inversion. To solve the order matrix differential equation (5), we will obtain initial value for the Newton-Raphson method. Inequality (6) is the convergence condition.
2.2. Multiray Surface Tracing
The computing ideas can also apply to the multiray surface tracing, which provides the basic algorithm for ray-surface intersection. Traditional methods for ray-surface intersection problem can be roughly classified into Newton-iteration-based ones and subdivision-based ones. Here we give a new method to solve this problem. See Figure 2. Similarly to point inversion, we reduce the problem to the projection of straight line segment onto the related surface.
The ray line can be expressed asFirstly, choose an initial point on the base surface , calculate the foot point projects from to as follows:Construct a line segment as follows:Then we have the following characteristic equation:Taking the derivative of (10) with respect to , we getAssumealong the projective curve segment. Finally we haveAdding initial value conditionsBy solving (13), we can obtain the parameters of point .
2.3. Improvement on the ODE Solver
ODE45 is a function implementing a Runge-Kutta method with a variable time step for efficient computation. ODE45 is designed to handle the following general problem:where is the independent variable, is a vector of dependent variables to be found and is a function of and . The mathematical problem is specified when on the right-hand side of (15) is set and the initial conditions at time are given.
The disadvantage of using ODE45 is slow convergence because of the variable time step. When the tracing point is very close to the , ODE45 has to change the step size very frequently. It takes plenty of time to satisfy convergence conditions. In order to overcome the disadvantages above, new method will use fixed step order Runge-Kutta method to trace along the base surface. Equations (5) and (13) can be converted as follows:where , .
The explicit -order fixed-step Runge-Kutta calculation formulas are as follows:where is fixed step satisfyingThe end computing conditions for multipoint inversion are inequality (6) and .
For multiray surface intersection, when the projection distance from surface point to the ray is small enough, the point is considered as exact solution. It can be described as follows:So, the convergence condition for multiray surface intersection algorithm is inequality (20) and . The two conditions can guarantee the algorithm not go into an infinite loop.
The precision of the -order Runge-Kutta is in every computation step, where is a positive arbitrary constant. After tracing along the base surface , the accumulated error will be very great. Paper [30] analysed systematically the accuracy of Runge-Kutta method. For a known surface , the truncation error between exact solution and resulting solution in each calculation step has an upper bound as follows:where and are positive constants independent of , and .
2.4. Newton-Raphson Iteration
The second-order algorithm of Hu et al. [31] is a good representation of Newton-Raphson method. By calculating the normal curvature circle of base surface at , a linear combination of tangent vectors is obtained. Connect the point with the centre of normal curvature circle . The intersection point of the line segment and the normal curvature circle is . Project the point along the normal vector at onto the base surface. The foot of perpendicular from to the base surface is approximated by a known function . With computed, a new set of parameters can be obtained for next iteration. However, the problem of the method is the stiffness equations which consist of two tangent vectors and a normal vector. According to the numerical experiment, if the point is too far from or the curvature of base surface is too small, the coefficients of tangent vectors will be much greater than the coefficients of normal vector. The phenomenon will cause the linear combination equation to have no solution, especially when solving multipoints inversion.
The new method combining the Runge-Kutta algorithm and the Newton-Raphson algorithm can avoid the problem in the second order of Hu et al. [8]. Meanwhile, it can also save lots of time, especially in calculating multipoints inversion. The total CPU time will be economized violently, with the calculation points increasing. For the calculation process, after tracing along the given surface, a result point will be found, which is not close enough to point . The point is taken as the initial value of the Newton-Raphson iteration. Usually, Newton-Raphson iteration needs a good initial value, so that the computation can converge. According to equation (19), selecting a suitable step size of ODE solver can assure the initial value good enough to the Newton-Raphson iteration. Piegl [32] shows a calculation process to solve point inversion by Newton-Raphson iteration. Without a good initial value, the method often causes numerical oscillation. There computation process is listed as follows:and two scalar equationsEquations (23) and (24) can be written as a matrix formwhere all items in the matrix and are calculated at . At the th iteration we need to solve a system of linear equations in the unknown given byConvergence criteria are given byHowever, for multiray surface intersection, the computation process can be described asin which is a parametric surface and is a ray in 3D space. Given a point and a ray, a plane containing them can be created as illustrated in Figure 2. Then the intersection of the generated plane and the surface products a curve , which is indicated as a dotted curve in Figure 2. Since the intersection curve is on the plane, it becomes a planar curve and subsequently the intersection point between the ray and the surface should lie on the intersection curve . So, the same method developed for the ray and 2D planar curve intersection case is directly applicable for this problem.
The basic approach to solve (30) is identical to the case of finding intersection between a ray and a planar curve. Given , the amount of update for is obtained by solving equation as follows:The Jacobian matrix is given byIn which , , and are column vectors.
Consider a point on the surface. The intersection between the ray and the tangent plane is computed by plugging the ray equation into the tangent plane equations, yielding the value of as follows:in which is the unit surface normal vector of .
Convergence criteria are given byThe tolerance is much smaller than the tolerance . If the inequality (29) or (34) is satisfied, then the calculation is considered successful.
3. Computational Process of the Proposed New Method
A single point calculation is taken as an example in the computational process to explain the whole process which combines tracing along the base surface and the Newton-Raphson iteration. Figure 3 shows the whole calculation flowchart for the new algorithm.
The main pseudocode of new method with Newton-Raphson iteration are showed in Algorithm 1. When solving multipoints inversion, the default coincides with .
| Require: | |
| 1: procedure MAIN | |
| 2:if then | |
| 3: | |
| 4:end if | |
| 5: | |
| 6: | |
| 7: | |
| 8: | |
| 9: | |
| 10:if then | |
| 11:while do | |
| 12: | |
| 13: | |
| 14: | |
| 15: | |
| 16:end while | |
| 17:else | |
| 18:while do | |
| 19: | |
| 20: | |
| 21: | |
| 22: | |
| 23:end while | |
| 24:end if | |
| 25:if then | |
| 26:while do | |
| 27: | |
| 28: | |
| 29: | |
| 30:end while | |
| 31:else | |
| 32:while do | |
| 33: | |
| 34: | |
| 35: | |
| 36:end while | |
| 37:end if | |
| 38: end procedure |
The pseudocodes of the ODE solver are showed in Algorithm 2. The pseudocode of Newton-Raphson iterations are showed in Algorithm 3.
In practical, the algorithm transforms the initial value problem of Newton iteration into the design problem of ODE solver error. By reasonably designing the error of the ODE solver, the Newton iteration solver can obtain an excellent initial value.
4. Convergence Rate Analysis
In this chapter we will analyse the convergence rate. There is a theorem about convergence of Newton iteration. The proof of the theorem can be found in section 10.2.2 of book [33].
Theorem. Assuming is -differentiable in an open neighbourhood of point , . Assuming is continuous at and is not singular. Then is an attraction point of Newton iterationIn addition, if there is a constant and satisfyingthen . At last, if is continuously differentiable on , and its second-order -derivative exists which satisfiesthen .
The theorem shows that there is always an attraction domain . If the initial approximation is in , the sequence generated by Newton-Raphson iteration always lies in and converges to . Roughly speaking, iterating once with Newton-Raphson iteration method can double the number of significant digits. That means Newton-Raphson iteration converges quickly. It also should be noted that the Newton-Raphson iteration method is self-correcting. That is to say, only depends on and . So, the rounding error generated by the previous iteration will not be passed down step by step. That is the reason why the combination of Newton-Raphson iteration and Runge-Kutta method can obtain precise result.
In this paper, the Runge-Kutta method and the Newton iteration we provided are both first-order convergent. However, Wang’s method needs to build multiline tracing in order to make the error satisfy the convergence condition. That will make the step size of ODE far less than the error . Therefore, Wang’s method has a precision limitation when using the fix-step Runge-Kutta method.
Moreover, there are two reasons for the proposed method being faster than Hu’s second order method. One is Hu’s second-order method needs to subdivide the base surface when calculating each point. That will cost a lot of time. The other one is related to its success rate. According to the experiment, Hu’s second-order often brings a failure rate around 6%. That will slow down the algorithm.
For single point calculation, Hu’s second order method has a higher convergence order than others. But when calculating large number points on base surface, Hu’s second-order method needs to subdivide the base surface in order to obtain a good initial point. Strictly speaking, when calculating a failure point, Hu’s second-order method needs to resubdivide the base surface for a better initial point. That is time wasting.
5. Implementation Examples
5.1. Multipoint Inversion Cases
Here we will use practical examples to demonstrate why the proposed method is better than the latest method [27] and traditional Newton iteration-related methods (we use method [8] as a representative). Generally speaking, the method [27] avoids the sensitivity to initial value of traditional Newton iteration-related methods. The process of calculation of Wang et al. [27] method should be divided into a lot of tracing segments. For example, if the output of the first line tracing along the base surface is , then will be the initial value of the second line segment. The line segment of is described as as follows:Afterwards, trace the base surface along with a small iteration step and repeat the tracing loop, until inequality (29) is satisfied. The repeated tracing along a curve on surface is time-killing.
For the trajectory of composite fibre path planning, a large number of points often need to be processed continuously on a given surface. That is a typical multipoint inversion problem. So successful calculation at one stroke is very important. We cannot afford even a few iteration failures, which, as we know, may not be avoided in most Newton-iteration-related methods. For example, in order to improve success rate of Hu’s second-order method, when the base surface is NURBS, a subdivision strategy to optimize the initial value for this method will be adopted. The main process of segmentation algorithm is to divide the NURBS surface into a lot of subsurfaces. In each subsurface select the parametric point as initial value. Proceed as follows:
(1) Divide NURBS surface into pieces of Bezier surfaces via inserting nodes algorithm.
(2) For each piece of Bezier surface, discrete the Bezier surface if is in its convex hull. Usually, the surface needs to be divided three times.
5.1.1. Hyperbolic Paraboloid Surface
Hyperbolic paraboloid surfaces are widely used in construction, manufacturing, aerospace and nuclear industries. The hyperbolic paraboloid surface can be described in Cartesian coordinate equation as follows:Convert (39) to parametric surface as the following characteristic equation:The calculation initial point is . The 2000 points need to be computed on base surface. See Figure 4(a). Figure 4(b) directly gives the result of single point calculation step for new method. It takes 106 iterations to complete the calculation of this point.
(a)
(b)
By comparing with Wang’s method, the results of calculating hyperbolic paraboloid are showed in Tables 1–3. Table 1 shows the time cost and failure number of points in calculating same 2000 points. From the table we can see that new method costs less time than Wang’s method.
Tables 2 and 3 give three times experiment in calculating random 2000 points in parametric domain for three different methods. In each experiment the new method costs less time than other methods and the calculation time has very small fluctuations. Wang’s method through multisegment tracing can also have a high success rate.
According to the results, both Wang’s method and new method have no failure points. However, the new method is an order of magnitude faster than Wang’s method.
5.1.2. Bezier Surface
Another important surface in modern industry is the Bezier surface, it is widely used in surface modelling. A Bezier surface can be expressed as . The control vertices is in Table 4.
The initial point is and the 2000 points will be computed on Bezier surface. See Figure 5(a). Figure 5(b) directly gives the result of single point calculation step for new method. It takes 14 iterations to complete the calculation of this point.
(a)
(b)
As explained in the beginning of this section, Hu et.al [8] elaborates an approach to segment the NURBS surface into Bezier surfaces. A Bezier surface also needs to be segmented into 64 small pieces to optimize initial point for the Newton-Raphson iteration. For this calculation example, the deCasteljau midpoint segmentation algorithm will be applied to solve surface segmentation. According to Hu et al. [31], the detachment times are generally selected as 3. However, the algorithm will cause segmentation error. In order to make sure that every point is in the convex hull of divided surface, the surface convex hull scope should be expanded.
Table 5 shows that the new method is superior to Wang’s method through multisegment tracing and Hu’s second-order method with subdivision. Another important phenomenon is that, for Hu’s method, the calculation success rate has a great improvement after surface subdivision. But, the subdivision algorithm will cost much more time in every point calculating process.
By randomized numerical experiment in three different methods, Tables 6–8 show that the presented method (new method) costs the least time among the three methods. It also has a great successful rate. Meanwhile, Table 7 shows that Wang’s method through multisegment tracing also has no failure in three numerical experiments. But the CPU time of Wang’s method through multisegment tracing is approximately 10 times as much as that of new method. At last, Hu’s second order method needs to conduct multilevel subdivision calculation. This is the reason why the time cost by Hu’s method is approximately 100 times than that of new method takes. Moreover the failure rate of Hu’s method with subdivision is about 6.07%. See Table 8.
Three comparative results of new method, Wang’s method through multisegment tracing and Hu’s second order method with subdivision, are displayed on Tables 6–8. With the total number of points growing, the CPU time cost using the three algorithms shows different properties. Here we choose a Bezier surface and a hyperbolic paraboloid surface as base surfaces to test the CPU time cost.
5.2. Multiray Surface Intersection
For multiray surface intersection, one of the industrial applications in optical engineering is freeform optical system. However, freeform optics that achieve excellent optical performance have broad application prospects in various areas, such as green energy, manufacturing, illumination, and biomedical [34–38]. The typical feature of such an imaging system is that the object plane or the phase plane is freeform surface. See Figure 6.
The aberration parameters such as NURBS surface are obtained by controlling the surface shape [39], because it is formed from piecewise splines. According to the geometric properties of the surface and the law of reflection the incident rays and reflected rays can be determined. Therefore, the core algorithm for freeform optical system is the calculation of feature image point which is closely related to multiray surface intersection.
The current case takes a human hand mesh surface as the object which has 66848 points and 133692 faces. See Figure 7. Project the object onto a Bezier surface screen by a point light source. The surface screen is a -order Bezier surface. The control points are showed in Table 9.
The calculation process for single point and the relative position of the point light source, object and surface screen are showed in Figure 8. The object locates in . The light source position is . The surface image of the object is shown in Figure 9. The parameters and results of multiray surface intersection algorithm for the case are showed in Table 10.
All the algorithm examples are tested in a MATLAB R2016a on Windows 7 with an Intel 4-core 3.7GHz CPU and 8G RAM.
6. Discussion
The results show that the algorithm combining the Newton-Raphson iteration and the Runge-Kutta method can solve multipoint inversion faster than Wang’s and Hu’s method. The new method also has a strong stability because the distance between initial point and point to be calculated has an upper bound for a known surface. However, the algorithm through multisegment tracing also shows a great success rate on calculating multipoint inversion. The advantage of this method is no sensitivity to initial value. The method can satisfy tolerance conditions through optimizing the step size of Runge-Kutta algorithm in each time tracing line segment. The second order of Hu’s method with subdivision can only be used for calculating point inversion on NURBS surface. The reason for high failure rate is 3 times detachment. If the failure points are recomputed with more detachment times, the failure points may be calculated successfully. But, it will cost more time than Table 8 showed. For multiray surface intersection, the result of present application supports that the algorithm can accurately calculate the imaging under multiray irradiation. In detail, the following advantages and precautions can be presented with the comparison of three algorithms:
(1) The point inversion and ray-surface intersection are very basic problems in computer graphics. There are a lot of algorithms to solve this problem. But, in industry application, engineers often have to face multipoint inversion and multiray surface intersection on base surface. This paper provides a new method to fill this gap.
(2) There is no need to subdivide surface for better initial value. Theoretically speaking, the algorithm absorbs the merit of the classic Newton-Raphson iteration and tracing along the base surface. The cumulative error of Runge-Kutta algorithm determines the initial value of the classic Newton-Raphson iteration. The suitable step size of Runge-Kutta algorithm is the essential key to guarantee the Newton-Raphson iteration convergence.
(3) The algorithm has a large convergence scope. It can solve multipoint inversion and multiray surface intersection not only on the NURBS surface but also on most small curvature open surfaces. In most cases, it can obtain a high success rate.
(4) Similarly to Wang’s method, the method needs to be given a curve from initial point to the point to be calculated. If the base surface has a complex topology, we can insert some middle points to be traced on the surface until the initial point is good enough for the Newton-Raphson iteration. For example, if the base surface is a tour patch surface, some intermediate points should be traced when calculating points with no line segment projection on base surface.
(5) The algorithm can also combine tracing along the base surface and Hu’s second-order method. The author thinks that the method tracing along the base surface can provide good initial value instead of surface subdivision in contrast to Hu’s second-order method.
(6) If the point lies on the boundary of base surface or coincident with the endpoints, the algorithms can also have an excellent calculation success rate.
7. Conclusions
In this paper, we have presented a novel method to solve multipoint inversion and multiray surface intersection problem. This method achieves better results in both accuracy and effectiveness, especially on small curvature open parametric surface. Different from Wang’s method, this method combines the tracing along base surface and the classic Newton-Raphson iteration. Better than Hu’s method [32], this method need not select initial value in computing every single point inversion. The new approach can save a lot of calculation time and can be used in many industrial applications.
On the other hand, for a NURBS surface, this method significantly increases the success rate of the algorithm compared with Hu’s method with three times subdivision. Another important advantage is that the time required to complete calculation is approximately proportional to the number of points on Bezier surface.
For algorithm stabilization, new method has a strong stability on both Bezier surface and hyperbolic parabolic surface. In calculating multipoints inversion problem, Wang’s method through multisegmentation tracing also has a good stability, but it costs much more time. Yet, second order of Hu’s method with subdivision is the most unstable one on Bezier surface among the three methods.
Finally, the algorithm can also provide excellent support for industrial application. Compared with traditional algorithms, our algorithm can acquire high precision phase in designing freeform optical system.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The work in this paper was supported by National Natural Science Foundation of China under grant No. 51575266.