NURBS Modeling and Curve Interpolation Optimization of 3D Graphics

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.


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

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][7][8][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 as v 0 , v 1 , and v 2 . The subscript is defined as the number of the angle.
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 (v 2 ) and the number of ∠ABD (v 0 ), 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. 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 are Q 1 ðx 1 ; y 1 ; z 1 Þ, Q 2 ðx 2 ; y 2 ; z 2 Þ, and Q 3 ðx 3 ; y 3 ; z 3 Þ. The section plane intersects the triangular patch. The intersection point of the plane and the edge Q 1 Q 2 is Q k ðx k ; y k ; z k Þ. The point satisfies Eq. (2): The coordinates of intersection point Q k can be obtained by solving Eq. (3).
where h is 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.
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.

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 Q k , k = 0,1,2…,n. It is necessary to create a p-th NURBS curve to interpolate these points. For the sake of simplicity, in the weight vector, let all w i ¼ 1. If a parameter value u k is specified for each data point Q k , and a reasonable node vector U ¼ u 0 ; u 1 ; …; u m f g is designed, then Eq. (4) can be obtained.
It can be seen that Eq. (4) is a system of linear equations with n + 1 control points P t as the independent variable and matrix of (n + 1) × (n + 1) as the coefficient matrix.  where R i,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 value u k and the node vector U = {u 0 , u 1 ,…u m }. It is assumed that all of the parameters satisfy u ∈ [0,1]. Generally, the chord length parameter method can be used to select the parameter value.
Similar to NURBS curves, it is necessary to calculate the reasonable parameters ðu k ; v l Þ, node vectors U and V at first. Using u k as an example, the algorithm is as follows: Step 1: Calculate the mean value u 0 l ; u 1 l ; …; u n l in V direction.
The algorithm of solving v l is similar to that of u k . Once u k and v l are obtained, the node vectors U and V can be obtained. For the control points: where It can be found that Eq. (8) is the interpolation equation for curve interpolation of data point Q k,j (k = 0,1, Similarly, when i is fixed and l changes, Eq. (9) is an interpolation equation for data point K i,l (l = 0,1,…,m). P i,j on 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 of U, solve the curve interpolation equation according to the parameter u k and the data point Q k,j (k = 0,1,…,n), then get K i,l .
Step 2: In the direction of V, solve the curve interpolation equation according to the parameters v l and the data point K i,l (l = 0,1,…,m), then get P i,j .

Method of Finding the Position Parameter Based on the Optimized Adams Algorithm
NURBS interpolation is the process of calculating the tool position P j+1 = {x(u j+1 ), y(u j+1 ), z(u j+1 )} after an interpolation cycle based on the federate v(t), the interpolation cycleT , and the current tool position P j = {x(u j ), y(u j ), z(u j )}. In the NURBS parameter curve C(u), the position parameter u j+1, which is after an interpolation period, needs to be computed based on the current position parameter u j in the parameter domain and the feed speed v(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 calculate u j+1 quickly.
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].
Once the function value y(x n ) of the node x n is 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 of y(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).
where y n ¼ yðx n Þ, y n 0 ¼ f n ¼ f ðx n ; y n Þ. r is the order of the Adams formula.
The common three-step fourth-order Adams explicit formula is shown in Eq. (13).
where h is the step value and the truncation error is O(t 5 ).
Applying the idea of an ordinary differential equation to the NURBS curve model, the following settings can be made: Make interpolation cycle T ¼ h, parameter node u j ¼ y n , then: According to the above settings, and referring to the most commonly used three-step four-order Adams formula, the method of calculating u j+1 is 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 are O(t 5 ). 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: þ 37 vðt jÀ2 Þ Cðu jÀ2 Þ À Cðu jÀ3 Þ u jÀ2 À u jÀ3 À 9 vðt jÀ3 Þ Cðu jÀ3 Þ À Cðu jÀ4 Þ u jÀ3 À u jÀ4 It can be seen from Eq. (17) that the value u j+1 of the position parameter can be calculated by the current position parameter and the four position parameters u j , u j-1 , u j-2 , u j-3 and u j-4 before the current position. At the same time, the previous multi-step position parameters u j , u j-1 , u j-2 …… are needed to calculate the current position parameter u j+1 by 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, let y n = y(x n ), y n+1 = y(x n+1 ), there is y nþ1 ¼ y n þ hy n 0 ¼ y n þ hf ðx n ; y n Þ; n ¼ 0; 1; … Eq. (18) is the Euler formula. f(x n , y n ) is the slope at point x n . 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 [x n , x n+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 parameters u j+1 of 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 parameters u 1 , u 2 , u 3 and u 4 are 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 speed v(t) can refer to the relevant discussion and calculation [15].

Hardening Implementation of Parallel Computing
It can be seen from the previous description that, for a NURBS curve with n parameters, 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][17][18].
The hardware structure of implementing the Runge-Kutta and Adams formulas by parallel computing is shown in Figs. 5 and 6. In Figs. 5 and 6, module cu is used to calculate C(u), module vt is used to calculate v(t), module ui_out is used to calculate the next position parameter u i+1 , module cuk is used to calculate C(u j-1 + (T/2)K i ), and module Kout is used to calculate K i . 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 Δt and the calculation time of the module cu is T, the parallel structure calculation of the Adams and Runge-Kutta formulas require T + Δt and 4T + 5Δt, respectively. Considering T >>Δ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. 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.

Simulation and Experiment
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: Degree p = 2. The control points are P 0 (1,2), P 1 (1.5,1), P 2 (3,3), P 3 (4,3.5), P 4 (5,3), P 5 (6.5,1), and P 6 (7,2). Node vector U = {0,0,0,1/5,2/5,3/5,4/5,1,1,1}, weight W = {1,1,1,1,1,1,1}. The interpolation path is shown in Fig. 8. The processing parameters are as follows: interpolation period T = 1 ms, maximum feed speed v m = 4 mm/s, and initial speed v s = 0 mm/s. Maximum bow height error ε = 1 μm. Reference [15] introduced a model to calculate v(t). The relationship between the position parameter u and the processing time is shown in Fig. 9. The relationship between the step value Δu of u and 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 7: NURBS modeling of a 3D image 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 cu T = 205 μs. The time required to calculate the first 12 position parameters u j+1 by serial calculation and parallel calculation is shown in Fig. 11. It can be seen that the time to calculate the position parameters u j+1 by 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.

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.  Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.