NURBS Modeling and Curve Interpolation Optimization of 3D Graphics

Hao Zhu,Mulan Wang,Kun Liu and Weiye Xu

1School of Information and Communication Engineering, Nanjing Institute of Technology, Nanjing, 211167,China

2Jiangsu Key Laboratory of Advanced Numerical Control Technology, Nanjing Institute of Technology, Nanjing, 211167,China

3Department of informatics, University of Leicester, Leicester, LE1 7RH, UK

Abstract: In order to solve the problem of complicated Non-Uniform Rational B-Splines (NURBS) modeling and improve the real-time performance of the high-order derivative of the curve interpolation process, the method of NURBS modeling based on the slicing and layering of triangular mesh is introduced.The research and design of NURBS curve interpolation are carried out from the two aspects of software algorithm and hardware structure.Based on the analysis of the characteristics of traditional computing methods with Taylor series expansion, the Adams formula and the Runge-Kutta formula are used in the NURBS curve interpolation process,and the process is then optimized according to the characteristics of NURBS interpolation.This can ensure accuracy, and avoid the calculation of higher-order derivatives.Furthermore,the hardware modules for the Adams and Runge-Kutta formulas are designed by using the parallel hardware construction technology of Field Programmable Gate Array (FPGA)chips.The parallel computing process using FPGA is compared with the traditional serial computing process using CPUs.Simulation and experimental results show that this scheme can improve the computational speed of the system and that the algorithm is feasible.

Keywords: NURBS modeling;interpolation;Adams;Runge-Kutta;FPGA

1 Introduction

The Non-Uniform Rational B-Splines (NURBS) model was proposed by Piegl et al.[1].This model provides a unified mathematical description method for analytic curves, analytic surfaces, free curves, and free surfaces.Due to its universality and simplicity of calculation, the NURBS model has become the only mathematical description method used in computer-aided design (CAD) to define product geometry.Having the NURBS curve interpolation function has become a key index for determining the performance of computer numerical control (CNC) systems [2].There are many ways to describe the NURBS curve.One of the most commonly used is the parameter description [3].

2 NURBS Surface Generation Based on Triangle Mesh Model

2.1 Slicing and Layering of Triangle Mesh Model

3D reconstruction technology based on triangular patch is a classic reverse engineering technology.Reference [4] applied artificial neural network technology to the reconstruction of triangular patches.Reference [5] proposed a method of generating triangular patch model based on shape from silhouette technology.The research of slicing and layering technology based on the triangle mesh model has enjoyed many successes [6-9].There are three categories of layered algorithms:The layered algorithm based on the location of triangle meshes, the layered algorithm based on the topological relationship of triangle meshes, and the layered algorithm based on out-of-core technology.In this paper, because the subsequent modeling requires an ordered point set, the layered algorithm based on the topological relationship of triangle meshes is used to obtain it.Reference [10] proposed a recursive slicing method based on the directed weighted graph.Reference [11] modified this method to improve the layering speed.In this improved method, a directed weighted sorting for the adjacent triangular patches was established.Using the triangular pyramid ABCD shown in Fig.1 as an example, its four faces are numbered 0, 1, 2, and 3.Face 0 is the under-side triangular patch BCD.Face 1 is the left-side triangular patch ABD.Face 2 is the right-side triangular patch ACD.Face 3 is the front-side triangular patch ABC.The three angles of each triangular patch are defined asv0,v1, andv2.The subscript is defined as the number of the angle.

Figure 1:Weighted sorting of triangle meshes

The information for each face includes its number and the three weights of the three edges.The weight of edge is defined as follows.For an edge of the current face, the weight is the sum of two angles’ numbers.These two angles satisfy the two following conditions:They are in the adjacent face (the face that shares the edge with the current face); and they are adjacent to the current face.As shown in Fig.1, edge AB is the edge shared by face 3 (the current face) and face 1 (the adjacent face).The weight of the edge AB is the sum of the number of ∠BAD (v2) and the number of ∠ABD (v0), i.e., 2 + 0 = 2.The rules for adding angle numbers are shown in Eq.(1).

The directed weighted vector graph is shown in Fig.2.The circles represent four triangular patches.The arrows point from the current face to the adjacent face.The value of the arrow indicates the weight.

Figure 2:Directed weighted vector graph of triangular patches

According to the above definition,every edge of each triangle patch can be recorded.The information of the triangular patch can be stored in the linked list structure.Each triangular patch corresponds to a chain node.Every node member includes the triangular patch number,the normal vector,the division mark,and the vertex coordinate.In addition,each edge is a child node of the triangular patch node.Each child node member includes the number of the triangle patch that shares this edge,the weight of this edge,and the address of the next edge.According to this structure,each edge of each triangular patch can be found in order when the triangular mesh model is layered.If the current triangular patch has been layered,it can be marked at the same time.

The core process of triangle mesh model layering is to use a series of planes perpendicular to the z axis of the rectangular coordinate to intersect with the triangular patch and record the intersection point set.As shown in Fig.3, the three vertices of the triangular patch areQ1(x1,y1,z1),Q2(x2,y2,z2), andQ3(x3,y3,z3).The section plane intersects the triangular patch.The intersection point of the plane and the edgeQ1Q2isQk(xk,yk,zk).The point satisfies Eq.(2):

Figure 3:Schematic diagram of intersection of section plane and triangular patch

The coordinates of intersection pointQkcan be obtained by solving Eq.(3).

wherehis the height of the current section plane.

As shown in Fig.4, when the section plane intersects with each triangle patch, two intersections are obtained.For two adjacent triangle patches, the intersection of the section plane and the common edge produces the repeating node.In the link list node structure, the segmentation mark can be used to record the layered state of the triangular patch by the current section plane.The repeating node can be eliminated through the above state.

Figure 4:Repeating node of adjacent triangular patches

The layering algorithm is described as follows:

Step 1:Prescreen the triangle patch according to the Z coordinate of each triangle, and filter out the triangular patches that do not intersect with the current section plane.

Step 2:Weighted sort the remaining triangular patches to determine the order.

Step 3:Traverse each triangle in order of weight,and store the calculated intersection points in the chain list.

The resulting linked list data are the ordered intersection data.

2.2 NURBS Curve Reconstruction

As described in Section 2.1, a series of ordered intersection points can be obtained.Reference [1]presented a method of interpolating non-rational B-spline curves and surfaces based on the ordered data points.The algorithm of interpolating NURBS curve and surface can be designed based on this idea.

Set the interpolation object as a series of ordered data points Qk,k=0,1,2…,n.It is necessary to create ap-th NURBS curve to interpolate these points.For the sake of simplicity,in the weight vector,let allwi=1.If a parameter valueis specified for each data point Qk,and a reasonable node vectorU={u0,u1,…,um}is designed,then Eq.(4)can be obtained.

where

It can be seen that Eq.(4)is a system of linear equations withn+1 control points Ptas the independent variable and matrix of(n+ 1) ×(n+ 1) as the coefficient matrix.

whereRi,p(u) can be calculated by the Cox-de Boor recurrence formula.Therefore, the key to solving this system of equations is how to select the appropriate parameter valueand the node vectorU= {u0,u1,…um}.It is assumed that all of the parameters satisfyu∈[0,1].Generally, the chord length parameter method can be used to select the parameter value.

2.3 NURBS Surface Interpolation

Suppose there are (n+ 1) × (m+ 1) data points Qk,l,k= 0,1,…n,l= 0,1,…,m.Surface interpolation creates a non-rational(p, q)degree B-spline surface to interpolate these points.

Step 1:Calculate the mean valueinVdirection.

The algorithm of solvingis similar to that of.Onceandare obtained,the node vectorsUandVcan be obtained.For the control points:

where

It can be found that Eq.(8)is the interpolation equation for curve interpolation of data point Qk,j(k=0,1,…,n)in theUdirection.Ki,lis the isoparametric curve control point of surface S(u,v)whenv=.Similarly,wheniis fixed andlchanges,Eq.(9)is an interpolation equation for data point Ki,l(l=0,1,…,m).Pi,jon the right side of the equation is the control point to be solved.Therefore,the algorithm to solve the control point is as follows:

Step 1:In the direction ofU,solve the curve interpolation equation according to the parameterand the data point Qk,j(k= 0,1,…,n), then get Ki,l.

Step 2:In the direction ofV,solve the curve interpolation equation according to the parametersand the data point Ki,l(l= 0,1,…,m), then get Pi,j.

No sooner had she reached her room than the Princess exclaimed, Now let us see what these fine plums can add to my beauty, and throwing off her hood52, she picked up a couple and ate them

3 NURBS Curve Interpolation Optimization Method

3.1 Method of Finding the Position Parameter Based on the Optimized Adams Algorithm

NURBS interpolation is the process of calculating the tool position Pj+1={x(uj+1),y(uj+1),z(uj+1)}after an interpolation cycle based on the federatev(t), the interpolation cycleT, and the current tool position Pj= {x(uj),y(uj),z(uj)}.In the NURBS parameter curve C(u), the position parameteruj+1,which is after an interpolation period, needs to be computed based on the current position parameterujin the parameter domain and the feed speedv(t) in the track domain.This process is commonly implemented by using the Taylor series expansion algorithm [12].

The disadvantage of the Taylor series expansion method is that if we want more accurate results, we need higher-order Taylor series expansion.That is to say, we need to compute the higher-order derivatives of the NURBS curves, which is a complex process.However, reducing the Taylor series expansion order will reduce the computational accuracy.In view of this problem, we can design an optimization algorithm based on the Adams formula and calculateuj+1quickly.

The Adams formula is a numerical integration method.Its basic purpose is to solve general ordinary differential equations such as Eq.(10)with a linear multi-step method[13].

In this equation,the step size ish, letxn=x0+nh(n =1, 2, …).

Once the function valuey(xn)of the nodexnis obtained according to a certain recursion principle,Eq.(10) can be solved.

In the Adams algorithm,take the slope values of the first few steps ofy(x).The weighted average of these values is taken as the slope value to estimate the current function value.When the selected point is the node before the current point,it is the explicit Adams format.When the selected points include the node before the current point and after the current point, it is the implicit Adams format.The explicit and implicit Adams formats are shown as Eqs.(11)and (12).

whereyn=y(xn),yn′=fn=f(xn,yn).ris the order of the Adams formula.

The common three-step fourth-order Adams explicit formula is shown in Eq.(13).

wherehis the step value and the truncation error isO(t5).

Applying the idea of an ordinary differential equation to the NURBS curve model,the following settings can be made:

Make interpolation cycleT=h,parameter nodeuj=yn,then:

According to the above settings,and referring to the most commonly used three-step four-order Adams formula, the method of calculatinguj+1is shown in Eq.(15).

It can be seen from Eq.(15) that for the same order Adams algorithm and Taylor series expansion algorithm, their truncation error is the same.Using the fourth-order algorithms as an example, both of them areO(t5).By comparison, the Adams explicit formula avoids the high-order derivative calculation encountered by Taylor series expansion.At the same time, it can greatly improve the operation speed.However, the first derivative still exists in the Adams formula.Therefore, we continue to optimize the algorithm by replacing the derivative with difference quotient of C(u).

Take Eqs.(14) and (16)into Eq.(15),we can obtain:

It can be seen from Eq.(17)that the valueuj+1of the position parameter can be calculated by the current position parameter and the four position parametersuj,uj-1,uj-2,uj-3anduj-4before the current position.At the same time,the previous multi-step position parametersuj,uj-1,uj-2……are needed to calculate the current position parameteruj+1by the Adams formula.The higher the order is, the more that pre-order position parameters are needed.Therefore, on the premise that the truncation error is small enough, the system cannot start itself only by the Adams algorithm.To solve this problem, we need to use the linear single step algorithm to obtain the initial n position parameters.Of the many single-step algorithms, the Runge-Kutta algorithm is a calculation method with high accuracy [14].In this paper, the algorithm is used to achieve the self-starting of the system.

The Runge-Kutta algorithm evolved from the Euler algorithm.For the first-order ordinary differential equation shown in Eq.(10), we can see that

Therefore, letyn=y(xn),yn+1=y(xn+1), there isyn+1=yn+hyn′=yn+hf(xn,yn),n=0,1,…

Eq.(18)is the Euler formula.f(xn,yn)is the slope at pointxn.Because the function value of the latter step is predicted only by the single point slope,the error of the Euler formula is large.The idea of the Runge-Kutta algorithm is to take several points in the interval[xn,xn+1]and take the weighted average value of its slope as the function value after the average slope prediction.This gives it higher accuracy.

The commonly used fourth-order Runge-Kutta formula is shown in Eq.(19).

By using the definitions of Eq.(14)and substituting the derivative with difference quotient,the fourthorder Runge-Kutta formula can be applied to the NURBS curve model to obtain Eq.(20).

According to the principle of the Runge-Kutta formula, the parametersuj+1of NURBS curve can be calculated by using Eq.(20) recursively.However, by comparing Eqs.(20) and (17), it can be seen that the Runge-Kutta method needs to use the recursive operation, which requires a large amount of calculation.Therefore, according to the Runge-Kutta method, the first four position parametersu1,u2,u3andu4are calculated by Eq.(20).According to the Adams method, with a relatively small amount of calculation, the rest of the position parameters are calculated by Eq.(17).Of these, the calculation method of the feed speedv(t) can refer to the relevant discussion and calculation [15].

3.2 Hardening Implementation of Parallel Computing

It can be seen from the previous description that, for a NURBS curve withnparameters, the Runge-Kutta formula of 4 times and the Adams formula of n-4 times need to be calculated.If the conventional CPU such as a microcontroller unit (MCU) or advanced reduced instruction set machine (ARM) is used for calculation, using the fourth-order Runge-Kutta and Adams formula as an example, the next parameter needs to be obtained by 8 times of serial calculation of C(u).This is undoubtedly a bottleneck for real-time CNC machining processes.For this type of operation, using a large number of logic units integrated in the FPGA to build parallel computing logic can greatly improve the system operation speed.At the same time, after the FPGA program is hardened, the calculation process is completed by hardware digital logic.Compared with the process of executing a software program with the CPU,the speed will be greatly improved[16-18].

The hardware structure of implementing the Runge-Kutta and Adams formulas by parallel computing is shown in Figs.5 and 6.

Figure 5:Parallel computing module structure of the Runge-Kutta formula

Figure 6:Parallel computing module structure of the Adams formula

In Figs.5 and 6,module cu is used to calculate C(u),module vt is used to calculatev(t),module ui_out is used to calculate the next position parameterui+1, module cuk is used to calculate C(uj-1+ (T/2)Ki), and module Kout is used to calculateKi.The calculation of module cu requires multiple clock cycles, which will take up most of the calculation time, while the Kout and ui_out modules can be calculated in one clock cycle by using the parallel multiplier and divider in the FPGA.If the clock period is set as Δtand the calculation time of the module cu isT, the parallel structure calculation of the Adams and Runge-Kutta formulas requireT+ Δtand 4T+ 5Δt, respectively.ConsideringT＞＞Δt, for NURBS curves with n parameters, the times for parallel and serial calculation are shown in Tab.1.It can be seen that the parallel computing will greatly shorten the computing time.

Table 1:Time comparison between serial computing and parallel computing

4 Simulation and Experiment

Fig.7 shows,from left to right,the point cloud model,the triangular patch model,the NURBS layered curve model,and the surface reconstruction model of the target.

Figure 7:NURBS modeling of a 3D image

To verify the correctness of the interpolation algorithm,this paper simulates the quadratic NURBS curve in MATLAB.The parameters of NURBS curve are as follows:Degreep=2.The control points are P0(1,2),P1(1.5,1),P2(3,3),P3(4,3.5),P4(5,3),P5(6.5,1),and P6(7,2).Node vectorU={0,0,0,1/5,2/5,3/5,4/5,1,1,1},weightW= {1,1,1,1,1,1,1}.The interpolation path is shown in Fig.8.The processing parameters are as follows:interpolation periodT= 1ms, maximum feed speedvm= 4mm/s, and initial speedvs= 0mm/s.Maximum bow height error ε = 1 μm.Reference [15] introduced a model to calculatev(t).The relationship between the position parameteruand the processing time is shown in Fig.9.The relationship between the step value Δuofuand the processing time is shown in Fig.10.As the curvature radius of the machining path changes, the step value of the position parameter also changes, so the position parameter does not have a linear relationship with the machining time.

Figure 8:NURBS path to be machined

Figure 9:Position parameter u

Next,the time consumption of serial computing and parallel computing was studied.The FPGA model was ep2c8q208,the system clock was 50 MHz,and then the clock cycle Δt=0.02 μs,with the computing time of cuT=205 μs.The time required to calculate the first 12 position parametersuj+1by serial calculation and parallel calculation is shown in Fig.11.It can be seen that the time to calculate the position parametersuj+1by serial calculation was more than the actual interpolation period of the machine tool,and therefore,it cannot meet the real-time requirements of machining.However,the time for the parallel calculation was far less than the actual interpolation period of the machine tool.

Figure 10:Position parameter step value

Figure 11:Position parameter calculation time

5 Conclusion

This paper introduces a method of NURBS surface reconstruction based on triangle mesh layering.Based on an analysis of existing algorithms, the Adams and Runge-Kutta algorithms were introduced into the process of NURBS interpolation, thereby avoiding the complex high-order calculation of the NURBS curve and reducing the amount of calculation from the perspective of the software algorithm.At the same time, the hardware parallel computing module was designed based on the FPGA parallel hardware architecture.This module replaces the common software computing mode based on the CPU and improves the computing speed.The analysis results show that the algorithm proposed in this paper can correctly obtain the position parameters of NURBS curve.Therefore, when combined with the hardware structure of parallel computing,this method can improve the real-time performance of the system.

Funding Statement:This work was supported by Natural Science Foundation of Universities in Jiangsu Province (Grant No.16KJA460001), Creation Foundation Project of Nanjing Institute of Technology(Grant No.:ZKJ201611).

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.