# An Approach to Computing Multipoint Inversion and Multiray Surface Intersection on Parametric Surface

Authors: Pei Jingyu; Wang Xiaoping; Zhang Leen
Journal: Mathematical Problems in Engineering
 (2019)

---

## 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.


---

## Body

## 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 O(log⁡N) 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 On 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 surfaceS(u,v),u,v∈D⊂R2 (a mapping from D→R3) with at least 2-order derivative existing. For any u,v∈D there are Su≠0 and Sv≠0. The problem is as follows: P is a point on the surface S, with its Cartesian coordinates known. See Figure 1. We want to compute its parametric coordinates u,v∈D. Firstly take a point P0(u0,v0) arbitrarily on the given surface S. Construct a line segment L(t) with P and P0 as follows:(1)Lt=tP+1-tP0,0≤t≤1Calculate a normal vector N0 of the known point P0 on the surface S. Then define a vector V as V=[N0×(P-P0)]×(P-P0), which represents the project direction of line segment L(t). Let the projected curve of L(t) onto the surface S be P(t), which can be characterized as follows(2)Pt-Lt×V=0.Computing the derivative on both sides of Equation (2) with respect to t, it follows that(3)Su×VSv×Vutvt=Lt×V.Taking the dot product of both sides of (3) with the vector Su and Sv separately yields(4)Su·Su×VSu·Sv×VSv·Su×VSv·Sv×Vutvt=Su·P-P0×VSv·P-P0×V.Thus, if the mixed product V,Su,Sv is not equal to 0, by solving the matrix equation (4), we have ut and vt as follows:(5)utvt=01Sv·Su×V1Su·Sv×V0Su·P-P0×VSv·P-P0×V=Sv·V×P-P0V·Su×SvSu·P-P0×VV·Su×Sv.(6)P-Sun,vn<ϵ1.Equation (5) is the govern equation of multipoint inversion. To solve the 1st order matrix differential equation (5), we will obtain initial value for the Newton-Raphson method. Inequality (6) is the convergence condition.Figure 1
Geometry diagram of projection method.


### 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 Figure2. Similarly to point inversion, we reduce the problem to the projection of straight line segment onto the related surface.Figure 2
Geometry diagram of multiray surface tracing.The ray line can be expressed as(7)Rt=O+kt,t≥0Firstly, choose an initial point P0u0,v0 on the base surface S, calculate the foot point Pj projects from P0 to Rt as follows:(8)Pj=O+P0-O·k·kk·kConstruct a line segment Lt as follows:(9)Lt=1-tP0+tPj,0≤t≤1Then we have the following characteristic equation:(10)Lt-O×Pt-O=0Taking the derivative of (10) with respect to t, we get(11)Su×Lt-Odudt+Sv×Lt-Odvdt=L′t×Pt-OAssume(12)Su×Sv·Lt-O≠0along the projective curve segment. Finally we have(13)dudt=Su,v-O×Pj-P0·SvSu×Sv·Lt-Odvdt=Pj-P0×Su,v-O·SuSu×Sv·Lt-OAdding initial value conditions(14)u0=u0,v0=v0By solving (13), we can obtain the parameters of point P.


### 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:(15)dxdt=ft,x,xt0=x0where t is the independent variable, x is a vector of dependent variables to be found and f(t,x) is a function of t and x. The mathematical problem is specified when f(t,x) on the right-hand side of (15) is set and the initial conditions x(t0)=x0 at time t0 are given.The disadvantage of using ODE45 is slow convergence because of the variable time step. When the tracing pointPn is very close to the P, 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 4th order Runge-Kutta method to trace along the base surface. Equations (5) and (13) can be converted as follows:(16)ut=Ft,uwhere ut=uttvtt, u=utvt.The explicit4th-order fixed-step Runge-Kutta calculation formulas are as follows:(17)un+1=un+∆t6K1+2K2+2K3+K4K1=Ftn,unK2=Ftn+∆t2,un+∆t2K1K3=Ftn+∆t2,un+∆t2K2K4=Ftn+∆t,un+∆tK3where ∆t is fixed step satisfying(18)ti+1-ti=∆t,i=0,1,2,…,n,(19)t0=0,tn=1The end computing conditions for multipoint inversion are inequality (6) and 0≤t≤1.For multiray surface intersection, when the projection distanced from surface point P to the ray Rt is small enough, the point P is considered as exact solution. It can be described as follows:(20)d=OP·kk≤ϵ1So, the convergence condition for multiray surface intersection algorithm is inequality (20) and 0≤t≤1. The two conditions can guarantee the algorithm not go into an infinite loop.The precision of the4th-order Runge-Kutta is c∆t5 in every computation step, where c is a positive arbitrary constant. After tracing along the base surface S, the accumulated error e will be very great. Paper [30] analysed systematically the accuracy of Runge-Kutta method. For a known surface S, the truncation error between exact solution and resulting solution E in each calculation step has an upper bound as follows:(21)E≤73720ML4h5where L and M are positive constants independent of t,u, and h.


### 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 P0, a linear combination of tangent vectors is obtained. Connect the point P with the centre of normal curvature circle Q. The intersection point of the line segment PQ and the normal curvature circle is X. Project the point X along the normal vector at P0 onto the base surface. The foot of perpendicular from X to the base surface is approximated by a known function f(Δt). With Δt 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 P is too far from P0 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 Pn will be found, which is not close enough to point P. The point Pn 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:(22)ru,v=Su,v-Pand two scalar equations(23)fu,v=ru,v·Suu,v(24)gu,v=ru,v·Svu,v.Equations (23) and (24) can be written as a matrix form(25)δi=∆u∆v=ui+1-uivi+1-vi(26)Ji=fufvgugv=Su2+r·SuuSu·Sv+r·SuvSu·Sv+r·SvuSv2+r·Svv(27)κi=-fgwhere all items in the matrix Ji and κi are calculated at ui,vi. At the ith iteration we need to solve a 2×2 system of linear equations in the unknown δi given by(28)Jiδi=κiConvergence criteria are given by(29)P-Sui,vi≤ϵ2However, for multiray surface intersection, the computation process can be described as(30)Gu,v,t=Su,v-Rt=0in which Su,v is a parametric surface and Rt is a ray in 3D space. Given a point P 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 c, 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 c. 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 x=u,v,tT, the amount of update for x is obtained by solving equation as follows:(31)Δx=-J-1GxThe Jacobian matrix J is given by(32)J=SuSvkIn which Su, Sv, and k are column vectors.Consider a pointP 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 t as follows:(33)t=P-O·Nk·Nin which N is the unit surface normal vector of P.Convergence criteria are given by(34)Sui,vi-O×kk≤ϵ2The tolerance ϵ2 is much smaller than the tolerance ϵ1. If the inequality (29) or (34) is satisfied, then the calculation is considered successful.


## 2.1. Multipoint Inversion

Let us consider a parametric surfaceS(u,v),u,v∈D⊂R2 (a mapping from D→R3) with at least 2-order derivative existing. For any u,v∈D there are Su≠0 and Sv≠0. The problem is as follows: P is a point on the surface S, with its Cartesian coordinates known. See Figure 1. We want to compute its parametric coordinates u,v∈D. Firstly take a point P0(u0,v0) arbitrarily on the given surface S. Construct a line segment L(t) with P and P0 as follows:(1)Lt=tP+1-tP0,0≤t≤1Calculate a normal vector N0 of the known point P0 on the surface S. Then define a vector V as V=[N0×(P-P0)]×(P-P0), which represents the project direction of line segment L(t). Let the projected curve of L(t) onto the surface S be P(t), which can be characterized as follows(2)Pt-Lt×V=0.Computing the derivative on both sides of Equation (2) with respect to t, it follows that(3)Su×VSv×Vutvt=Lt×V.Taking the dot product of both sides of (3) with the vector Su and Sv separately yields(4)Su·Su×VSu·Sv×VSv·Su×VSv·Sv×Vutvt=Su·P-P0×VSv·P-P0×V.Thus, if the mixed product V,Su,Sv is not equal to 0, by solving the matrix equation (4), we have ut and vt as follows:(5)utvt=01Sv·Su×V1Su·Sv×V0Su·P-P0×VSv·P-P0×V=Sv·V×P-P0V·Su×SvSu·P-P0×VV·Su×Sv.(6)P-Sun,vn<ϵ1.Equation (5) is the govern equation of multipoint inversion. To solve the 1st order matrix differential equation (5), we will obtain initial value for the Newton-Raphson method. Inequality (6) is the convergence condition.Figure 1
Geometry diagram of projection method.


## 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 Figure2. Similarly to point inversion, we reduce the problem to the projection of straight line segment onto the related surface.Figure 2
Geometry diagram of multiray surface tracing.The ray line can be expressed as(7)Rt=O+kt,t≥0Firstly, choose an initial point P0u0,v0 on the base surface S, calculate the foot point Pj projects from P0 to Rt as follows:(8)Pj=O+P0-O·k·kk·kConstruct a line segment Lt as follows:(9)Lt=1-tP0+tPj,0≤t≤1Then we have the following characteristic equation:(10)Lt-O×Pt-O=0Taking the derivative of (10) with respect to t, we get(11)Su×Lt-Odudt+Sv×Lt-Odvdt=L′t×Pt-OAssume(12)Su×Sv·Lt-O≠0along the projective curve segment. Finally we have(13)dudt=Su,v-O×Pj-P0·SvSu×Sv·Lt-Odvdt=Pj-P0×Su,v-O·SuSu×Sv·Lt-OAdding initial value conditions(14)u0=u0,v0=v0By solving (13), we can obtain the parameters of point P.


## 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:(15)dxdt=ft,x,xt0=x0where t is the independent variable, x is a vector of dependent variables to be found and f(t,x) is a function of t and x. The mathematical problem is specified when f(t,x) on the right-hand side of (15) is set and the initial conditions x(t0)=x0 at time t0 are given.The disadvantage of using ODE45 is slow convergence because of the variable time step. When the tracing pointPn is very close to the P, 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 4th order Runge-Kutta method to trace along the base surface. Equations (5) and (13) can be converted as follows:(16)ut=Ft,uwhere ut=uttvtt, u=utvt.The explicit4th-order fixed-step Runge-Kutta calculation formulas are as follows:(17)un+1=un+∆t6K1+2K2+2K3+K4K1=Ftn,unK2=Ftn+∆t2,un+∆t2K1K3=Ftn+∆t2,un+∆t2K2K4=Ftn+∆t,un+∆tK3where ∆t is fixed step satisfying(18)ti+1-ti=∆t,i=0,1,2,…,n,(19)t0=0,tn=1The end computing conditions for multipoint inversion are inequality (6) and 0≤t≤1.For multiray surface intersection, when the projection distanced from surface point P to the ray Rt is small enough, the point P is considered as exact solution. It can be described as follows:(20)d=OP·kk≤ϵ1So, the convergence condition for multiray surface intersection algorithm is inequality (20) and 0≤t≤1. The two conditions can guarantee the algorithm not go into an infinite loop.The precision of the4th-order Runge-Kutta is c∆t5 in every computation step, where c is a positive arbitrary constant. After tracing along the base surface S, the accumulated error e will be very great. Paper [30] analysed systematically the accuracy of Runge-Kutta method. For a known surface S, the truncation error between exact solution and resulting solution E in each calculation step has an upper bound as follows:(21)E≤73720ML4h5where L and M are positive constants independent of t,u, and h.


## 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 P0, a linear combination of tangent vectors is obtained. Connect the point P with the centre of normal curvature circle Q. The intersection point of the line segment PQ and the normal curvature circle is X. Project the point X along the normal vector at P0 onto the base surface. The foot of perpendicular from X to the base surface is approximated by a known function f(Δt). With Δt 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 P is too far from P0 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 Pn will be found, which is not close enough to point P. The point Pn 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:(22)ru,v=Su,v-Pand two scalar equations(23)fu,v=ru,v·Suu,v(24)gu,v=ru,v·Svu,v.Equations (23) and (24) can be written as a matrix form(25)δi=∆u∆v=ui+1-uivi+1-vi(26)Ji=fufvgugv=Su2+r·SuuSu·Sv+r·SuvSu·Sv+r·SvuSv2+r·Svv(27)κi=-fgwhere all items in the matrix Ji and κi are calculated at ui,vi. At the ith iteration we need to solve a 2×2 system of linear equations in the unknown δi given by(28)Jiδi=κiConvergence criteria are given by(29)P-Sui,vi≤ϵ2However, for multiray surface intersection, the computation process can be described as(30)Gu,v,t=Su,v-Rt=0in which Su,v is a parametric surface and Rt is a ray in 3D space. Given a point P 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 c, 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 c. 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 x=u,v,tT, the amount of update for x is obtained by solving equation as follows:(31)Δx=-J-1GxThe Jacobian matrix J is given by(32)J=SuSvkIn which Su, Sv, and k are column vectors.Consider a pointP 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 t as follows:(33)t=P-O·Nk·Nin which N is the unit surface normal vector of P.Convergence criteria are given by(34)Sui,vi-O×kk≤ϵ2The tolerance ϵ2 is much smaller than the tolerance ϵ1. 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. Figure3 shows the whole calculation flowchart for the new algorithm.Figure 3
The calculation flowchart.The main pseudocode of new method with Newton-Raphson iteration are showed in Algorithm1. When solving multipoints inversion, the default O coincides with P0.Algorithm 1:The Main algorithm.
Require:  Su,v⋁P⋁P0⋁u0 1: procedure Main 2:if  O-P0≠0  then 3:k←(O-P0)/O-P0 4:end if 5:t←0 6:h←h0 7:u←ODE(h,u0,P,P0,O,k) 8:Pn←S(u[1],u[2]) 9:t←t+h 10:if  O-P0=0  then 11:while  Pn-P>ϵ1&t≤1  do 12:u0←u 13:u←ODE(h,u0,P,P0,O,k) 14:Pn←S(u[1],u[2]) 15:t←t+h 16:end while 17:else 18:while  OP⋅k/k>ϵ1&t≤1  do 19:u0←u 20:u←ODE(h,u0,P,P0,O,k) 21:Pn←S(u[1],u[2]) 22:t←t+h 23:end while 24:end if 25:if  O-P0=0  then 26:while  Pn-P>ϵ2  do 27:u0←u 28:u←NewtonIteration(u0,P.O,k) 29:Pn←S(u[1],u[2]) 30:end while 31:else 32:while  |OP⋅k|/k>ϵ2  do 33:u0←u 34:u←NewtonIteration(u0,P.O,k) 35:←S(u[1],u[2]) 36:end while 37:end if 38: end procedureThe pseudocodes of the ODE solver are showed in Algorithm2. The pseudocode of Newton-Raphson iterations are showed in Algorithm 3.Algorithm 2:The ODE Solver.
Require:  Su←∂S/∂u⋁Sv←∂S/∂v 1: function ODE 2:for  i=0→3  do 3:if  O-P0=0  then 4:dui←Equ(5)(P,P0,ut,h) 5:else 6:dui←Equ(13)(P,P0,ut,h,k) 7:end if 8:end for 9:K1←du0 10:u1←du0+h×K1/2 11:K2←du1 12:u2←du1+h×K2/2 13:K3←du2 14:u3←du2+h×K3 15:K4←du3 16:u3←du3 17:y←u+h/6(K1+2K2+2K3+K4)  return  y 18: end functionAlgorithm 3:The Newton Iteration.
Require:  Su←∂S/∂u⋁Sv←(∂S/∂v)Su←∂S/∂u⋁Sv←∂S/∂v⋁Suu←∂2S/∂u2⋁Svv←∂2S/∂v2⋁Suv← ∂ 2 S / ∂ u ∂ v ⋁ S v u ← ∂ 2 S / ∂ v ∂ u 1: function NewtonIteration 2:if  O-P0=0  then 3:J←Equ(26)(P,u0) 4:delta←Equ(25)(P,u0) 5:u←u+delta 6:else 7:dui←Equ(32)(k,u0) 8:δx←Equ(31)(k,u0) 9:u[1]←u[1]+δx[1] 10:u[2]←u[2]+δx[2] 11:end ifreturn  u 12: end functionIn 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 F:D⊂Rn→Rn is G-differentiable in an open neighbourhood S0⊂D of point Fx∗=0, x∗∈D. Assuming F′ is continuous at x∗and F′x∗ is not singular. Then x∗ is an attraction point of Newton iteration(35)J:xk+1=xk-F′xk-1Fxk,k=0,1,…,In addition, if there is a constant α<+∞ and p∈0,1 satisfying(36)F′x-F′x∗≤αx-x∗p,∀x∈S0then ORJ,x∗≥OQJ,x∗≥1+p. At last, if F is continuously differentiable on S0, and its second-order F-derivative exists which satisfies(37)F′′x∗hh≠0,∀h∈Rn,h≠0then ORJ,x∗=OQJ,x∗=2.The theorem shows that there is always an attraction domainS0. If the initial approximation is in S0, the sequence xk generated by Newton-Raphson iteration always lies in S0 and converges to x∗. 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, xk+1 only depends on F and xk. 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ϵ2. 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 Pn1, then Pn1 will be the initial value of the second line segment. The line segment of PPn1 is described as L2(t) as follows:(38)L2t=tP+1-tPn1,0≤t≤1Afterwards, trace the base surface along L2 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[0.5,0.5] as initial value. Proceed as follows:(1) Divide NURBS surface into(n-k+1)×(m-k+1) pieces of Bezier surfaces via inserting nodes algorithm.(2) For each piece of Bezier surface, discrete the Bezier surface ifP 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:(39)z=x24-y29Convert (39) to parametric surface as the following characteristic equation:(40)Su,v=u,v,u24-v29The calculation initial point P0 is 0,0,0. 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.Figure 4
(a) Points on hyperbolic paraboloid. (b) Single point calculation step on hyperbolic paraboloid for new method.
(a) (b)By comparing with Wang’s method, the results of calculating hyperbolic paraboloid are showed in Tables1–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.Table 1
Two different methods in solving multipoint inversion on hyperbolic paraboloid.
New Method Wang’s Method Number of points 2000 2000 CPU Time(s) 29.322629806545390 273.7031542564848 Error 10 - 5 10 - 5 Failure 0 0 Maximum iteration 200 200 Parametric range [ - 300,300 ] [ - 300,300 ]Table 2
Experiment of new method for hyperbolic paraboloid.
New Method CPU Time(s) Failure Experiment 1 (2000 points) 31.163845291790800 0 Experiment 2 (2000 points) 30.898390837119763 0 Experiment 3 (2000 points) 30.926868387619347 0Table 3
Experiment of Wang’s Method for hyperbolic paraboloid.
Wang’s Method CPU Time(s) Failure Experiment 1 (2000 points) 255.5892196051562 0 Experiment 2 (2000 points) 259.5492855696442 0 Experiment 3 (2000 points) 265.7260784417069 0Tables2 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 asSu,v=∑i=0m∑j=0nbi,jBi,muBj,nv. The control vertices bi,j is in Table 4.Table 4
Control Points of Bezier Surface.
j = 0 1 2 3 i=0 0,0 , 0 0,1.116,0.067 0,2.232,0.134 0.433,3.132,0.076 1 - 0.616,0.433,1.217 - 0.183,1.333,1.158 - 0.183,2.449,1.225 1.116,2.915,0.917 2 - 0.366,0.433,2.183 - 0.366,1.549,2.25 - 0.366,2.665,2.317 0.933,3.132,2.009 3 0.75,0 , 2.900 0.75,1.116,2.967 0.75,2.232,3.033 2.049,2.699,2.725 4 1.866 , - 0.433,3.616 1.433,0.900,3.808 1.433,2.016,3.875 2.732,2.483,3.567The initial pointP0 is 0,0,0 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.Figure 5
(a) Points on Bezier surface. (b) Single point calculation step on Bezier surface for new method.
(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.Table5 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.Table 5
Three different methods in solving multipoint inversion on Bezier surface.
New Method Wang’s Method Hu’s second order method Number of points 2000 2000 2000 CPU Time(s) 6.926983581000674 69.154373822805740 714.2510553485881 Error 10 - 5 10 - 5 10 - 5 Failure 0 0 112 Maximum iteration 200 200 200 Point range [ 0,1 ] [ 0,1 ] [ 0,1 ]By randomized numerical experiment in three different methods, Tables6–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.Table 6
Experiment of new method for Bezier surface.
New method CPU Time(s) Failure Experiment 1 (2000 points) 6.962026113531086 0 Experiment 2 (2000 points) 6.546980250216993 0 Experiment 3 (2000 points) 6.702481947743534 0Table 7
Experiment of Wang’s Method for Bezier surface.
Wang’s Method CPU Time(s) Failure Experiment 1 (2000 points) 68.42192473628629 0 Experiment 2 (2000 points) 66.11651475984002 0 Experiment 3 (2000 points) 66.90059448796946 0Table 8
Experiment of Hu’s second order method for Bezier surface.
Hu’s second order method CPU Time(s) Failure Experiment 1 (2000 points) 691.0428227734579 132 Experiment 2 (2000 points) 667.2941246635530 109 Experiment 3 (2000 points) 665.4855344774451 123Three comparative results of new method, Wang’s method through multisegment tracing and Hu’s second order method with subdivision, are displayed on Tables6–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.Figure 6
Rays reflecting from freeform surface.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 Figure7. Project the object onto a Bezier surface screen by a point light source. The surface screen is a 3×4-order Bezier surface. The control points are showed in Table 9.Table 9
Control points of Bezier surface screen.
j = 0 1 2 3 i=0 0,0 , 0 0,20,10 0,40,20 0,60,20 1 20,0 , 20 20,20,20 20,40,30 20,60,10 2 40,0 , 4 40,20,30 40,40,40 40,60,20 3 60,0 , 0 60,20,10 60,40,20 60,60,0 4 80,0 , 0 80,20,0 80,40,10 80,60,0Figure 7
Mesh surface for human hand.The calculation process for single point and the relative position of the point light source, object and surface screen are showed in Figure8. The object locates in 40,50,-35. The light source position is 39.59,96.44,-19.78. 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.Table 10
Algorithm parameters and results for the case.
Number of points CPU Time(s) Error Failure Maximum Iteration Step of ODE Solver 66848 246.0136530198419 10 - 5 0 200 0.1Figure 8
Single point calculation process.Figure 9
Image of the hand mesh object.All the algorithm examples are tested in a MATLAB R2016a on Windows 7 with an Intel 4-core 3.7GHz CPU and 8G RAM.


## 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 Pn1, then Pn1 will be the initial value of the second line segment. The line segment of PPn1 is described as L2(t) as follows:(38)L2t=tP+1-tPn1,0≤t≤1Afterwards, trace the base surface along L2 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[0.5,0.5] as initial value. Proceed as follows:(1) Divide NURBS surface into(n-k+1)×(m-k+1) pieces of Bezier surfaces via inserting nodes algorithm.(2) For each piece of Bezier surface, discrete the Bezier surface ifP 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:(39)z=x24-y29Convert (39) to parametric surface as the following characteristic equation:(40)Su,v=u,v,u24-v29The calculation initial point P0 is 0,0,0. 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.Figure 4
(a) Points on hyperbolic paraboloid. (b) Single point calculation step on hyperbolic paraboloid for new method.
(a) (b)By comparing with Wang’s method, the results of calculating hyperbolic paraboloid are showed in Tables1–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.Table 1
Two different methods in solving multipoint inversion on hyperbolic paraboloid.
New Method Wang’s Method Number of points 2000 2000 CPU Time(s) 29.322629806545390 273.7031542564848 Error 10 - 5 10 - 5 Failure 0 0 Maximum iteration 200 200 Parametric range [ - 300,300 ] [ - 300,300 ]Table 2
Experiment of new method for hyperbolic paraboloid.
New Method CPU Time(s) Failure Experiment 1 (2000 points) 31.163845291790800 0 Experiment 2 (2000 points) 30.898390837119763 0 Experiment 3 (2000 points) 30.926868387619347 0Table 3
Experiment of Wang’s Method for hyperbolic paraboloid.
Wang’s Method CPU Time(s) Failure Experiment 1 (2000 points) 255.5892196051562 0 Experiment 2 (2000 points) 259.5492855696442 0 Experiment 3 (2000 points) 265.7260784417069 0Tables2 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 asSu,v=∑i=0m∑j=0nbi,jBi,muBj,nv. The control vertices bi,j is in Table 4.Table 4
Control Points of Bezier Surface.
j = 0 1 2 3 i=0 0,0 , 0 0,1.116,0.067 0,2.232,0.134 0.433,3.132,0.076 1 - 0.616,0.433,1.217 - 0.183,1.333,1.158 - 0.183,2.449,1.225 1.116,2.915,0.917 2 - 0.366,0.433,2.183 - 0.366,1.549,2.25 - 0.366,2.665,2.317 0.933,3.132,2.009 3 0.75,0 , 2.900 0.75,1.116,2.967 0.75,2.232,3.033 2.049,2.699,2.725 4 1.866 , - 0.433,3.616 1.433,0.900,3.808 1.433,2.016,3.875 2.732,2.483,3.567The initial pointP0 is 0,0,0 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.Figure 5
(a) Points on Bezier surface. (b) Single point calculation step on Bezier surface for new method.
(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.Table5 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.Table 5
Three different methods in solving multipoint inversion on Bezier surface.
New Method Wang’s Method Hu’s second order method Number of points 2000 2000 2000 CPU Time(s) 6.926983581000674 69.154373822805740 714.2510553485881 Error 10 - 5 10 - 5 10 - 5 Failure 0 0 112 Maximum iteration 200 200 200 Point range [ 0,1 ] [ 0,1 ] [ 0,1 ]By randomized numerical experiment in three different methods, Tables6–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.Table 6
Experiment of new method for Bezier surface.
New method CPU Time(s) Failure Experiment 1 (2000 points) 6.962026113531086 0 Experiment 2 (2000 points) 6.546980250216993 0 Experiment 3 (2000 points) 6.702481947743534 0Table 7
Experiment of Wang’s Method for Bezier surface.
Wang’s Method CPU Time(s) Failure Experiment 1 (2000 points) 68.42192473628629 0 Experiment 2 (2000 points) 66.11651475984002 0 Experiment 3 (2000 points) 66.90059448796946 0Table 8
Experiment of Hu’s second order method for Bezier surface.
Hu’s second order method CPU Time(s) Failure Experiment 1 (2000 points) 691.0428227734579 132 Experiment 2 (2000 points) 667.2941246635530 109 Experiment 3 (2000 points) 665.4855344774451 123Three comparative results of new method, Wang’s method through multisegment tracing and Hu’s second order method with subdivision, are displayed on Tables6–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.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:(39)z=x24-y29Convert (39) to parametric surface as the following characteristic equation:(40)Su,v=u,v,u24-v29The calculation initial point P0 is 0,0,0. 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.Figure 4
(a) Points on hyperbolic paraboloid. (b) Single point calculation step on hyperbolic paraboloid for new method.
(a) (b)By comparing with Wang’s method, the results of calculating hyperbolic paraboloid are showed in Tables1–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.Table 1
Two different methods in solving multipoint inversion on hyperbolic paraboloid.
New Method Wang’s Method Number of points 2000 2000 CPU Time(s) 29.322629806545390 273.7031542564848 Error 10 - 5 10 - 5 Failure 0 0 Maximum iteration 200 200 Parametric range [ - 300,300 ] [ - 300,300 ]Table 2
Experiment of new method for hyperbolic paraboloid.
New Method CPU Time(s) Failure Experiment 1 (2000 points) 31.163845291790800 0 Experiment 2 (2000 points) 30.898390837119763 0 Experiment 3 (2000 points) 30.926868387619347 0Table 3
Experiment of Wang’s Method for hyperbolic paraboloid.
Wang’s Method CPU Time(s) Failure Experiment 1 (2000 points) 255.5892196051562 0 Experiment 2 (2000 points) 259.5492855696442 0 Experiment 3 (2000 points) 265.7260784417069 0Tables2 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 asSu,v=∑i=0m∑j=0nbi,jBi,muBj,nv. The control vertices bi,j is in Table 4.Table 4
Control Points of Bezier Surface.
j = 0 1 2 3 i=0 0,0 , 0 0,1.116,0.067 0,2.232,0.134 0.433,3.132,0.076 1 - 0.616,0.433,1.217 - 0.183,1.333,1.158 - 0.183,2.449,1.225 1.116,2.915,0.917 2 - 0.366,0.433,2.183 - 0.366,1.549,2.25 - 0.366,2.665,2.317 0.933,3.132,2.009 3 0.75,0 , 2.900 0.75,1.116,2.967 0.75,2.232,3.033 2.049,2.699,2.725 4 1.866 , - 0.433,3.616 1.433,0.900,3.808 1.433,2.016,3.875 2.732,2.483,3.567The initial pointP0 is 0,0,0 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.Figure 5
(a) Points on Bezier surface. (b) Single point calculation step on Bezier surface for new method.
(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.Table5 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.Table 5
Three different methods in solving multipoint inversion on Bezier surface.
New Method Wang’s Method Hu’s second order method Number of points 2000 2000 2000 CPU Time(s) 6.926983581000674 69.154373822805740 714.2510553485881 Error 10 - 5 10 - 5 10 - 5 Failure 0 0 112 Maximum iteration 200 200 200 Point range [ 0,1 ] [ 0,1 ] [ 0,1 ]By randomized numerical experiment in three different methods, Tables6–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.Table 6
Experiment of new method for Bezier surface.
New method CPU Time(s) Failure Experiment 1 (2000 points) 6.962026113531086 0 Experiment 2 (2000 points) 6.546980250216993 0 Experiment 3 (2000 points) 6.702481947743534 0Table 7
Experiment of Wang’s Method for Bezier surface.
Wang’s Method CPU Time(s) Failure Experiment 1 (2000 points) 68.42192473628629 0 Experiment 2 (2000 points) 66.11651475984002 0 Experiment 3 (2000 points) 66.90059448796946 0Table 8
Experiment of Hu’s second order method for Bezier surface.
Hu’s second order method CPU Time(s) Failure Experiment 1 (2000 points) 691.0428227734579 132 Experiment 2 (2000 points) 667.2941246635530 109 Experiment 3 (2000 points) 665.4855344774451 123Three comparative results of new method, Wang’s method through multisegment tracing and Hu’s second order method with subdivision, are displayed on Tables6–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.Figure 6
Rays reflecting from freeform surface.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 Figure7. Project the object onto a Bezier surface screen by a point light source. The surface screen is a 3×4-order Bezier surface. The control points are showed in Table 9.Table 9
Control points of Bezier surface screen.
j = 0 1 2 3 i=0 0,0 , 0 0,20,10 0,40,20 0,60,20 1 20,0 , 20 20,20,20 20,40,30 20,60,10 2 40,0 , 4 40,20,30 40,40,40 40,60,20 3 60,0 , 0 60,20,10 60,40,20 60,60,0 4 80,0 , 0 80,20,0 80,40,10 80,60,0Figure 7
Mesh surface for human hand.The calculation process for single point and the relative position of the point light source, object and surface screen are showed in Figure8. The object locates in 40,50,-35. The light source position is 39.59,96.44,-19.78. 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.Table 10
Algorithm parameters and results for the case.
Number of points CPU Time(s) Error Failure Maximum Iteration Step of ODE Solver 66848 246.0136530198419 10 - 5 0 200 0.1Figure 8
Single point calculation process.Figure 9
Image of the hand mesh object.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 Table8 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.



---

*Source: 3790762-2019-06-11.xml*