Generation of spherical non-uniform rational basis spline curves and its application in five-axis machining

Abstract

A method of generating spherical non-uniform rational basis spline curves based on De Boor’s algorithm is presented in this article. The spherical curve preserves many good properties from non-uniform rational basis spline curves in Euclidean space, such as local modification property, convex hull property, rotation invariant property, knot insertion property, and so on. Construction of closed spherical non-uniform rational basis spline curve will be discussed too. Furthermore, a progressive iterative approximation scheme is proposed to approximate the given points on unit sphere using spherical non-uniform rational basis spline curves. The convergence and curve continuity will be discussed. Since the analytical expression of the spherical non-uniform rational basis spline curve is messy, numerical differentiation is used to test its continuity at interior knots. Finally, several examples are presented to verify the effectiveness of the proposed method. The proposed method is applied on tool path generation of five-axis machining to avoid interference by adjusting fewer directions and obtain smooth tool path simultaneously.

Introduction

Generating smooth spherical curves on unit sphere can be applied in many situations, such as designing spherical curve, generating smooth rigid motion in computer animation, and tool path generation of multi-axis machining.^1–4 Since B-spline has many good properties, it has already been widely used in curves and surface designs.^5–9 Therefore, much effort has been made to extend B-spline curve in Euclidean space to the construction of spherical B-spline curves.

In 1979, Parker and Denham proposed a method to interpolate in the ambient space E ^{m+ 1} and normalize to produce a curve in S^m . This method may cause extra geometrical distortion. ¹⁰ In 1985, in order to solve the problem of rigid motion interpolation in computer animation, Shoemake ¹ proposed the spherical linear interpolation for two points by using spherical distance instead of Euclidean distance. Furthermore, the classical De Casteljau’s algorithm was generalized to 3-sphere S ³. The interpolation spherical curve has C ¹ continuity at the curve joint. In 1994, Pobegailo ² proposed a method for designing spherical curves by two weighted spatial rotations. The generated curve has C ¹ continuity, local control, and invariance orthogonal transformations of coordinate system. In 1992, Nielson and colleagues^11,12 proposed a spherical B-spline method to construct spherical curve. The inner points of the cubic segment are computed at 1/3 and 2/3 of the geodesic distance between the control points on sphere. In his article, Nielson ¹¹ pointed out that this is the smoothest method he observed. But Kim et al.^13,14 thought that the constructed spherical B-spline curve is not C ² continuous, and furthermore, they proposed a general construction method for unit quaternion curves. The spline curve is reformulated in a cumulative basis form, and the corresponding quaternion curve is generated by converting each vector addition into the quaternion multiplication. This method is easy to evaluate high-order derivatives. But if the order of control points is reversed, the shape of spherical curve will change. ¹⁵ Therefore, it cannot be used to construct closed spherical curve. Also in 1995, Kim and Nam ¹⁶ presented a method interpolating a given sequence of solid orientation using circular blending quaternion curves. In 2001, Wang and Tang ¹⁷ tried to construct a spherical non-uniform rational basis spline (NURBS) curve, but it cannot guarantee that the generated curve lies on the unit sphere. In 1997, Park and Ravani ¹⁸ interpolated multiple points in SO(3) by regarding SO(3) as a Lie group equipped with a natural Riemannian metric. Russ and Fillmore ¹⁵ proposed a method for computing weighted averages on unit sphere based on the least square method. The existence and uniqueness properties of the weighted averages are proved. Two iterative algorithms are presented to generate smooth average quaternion. In 2005, Séquin et al. ¹⁹ proposed a class of circle spline that produced fair-looking G ²-continuous curves without any cusps, which can be used on sphere for spherical curve design. Based on Shoemake’s work, Popiel and Noakes ²⁰ proposed a method to piece spherical Bézier curves together into C ² splines. Since the velocity and acceleration are considered, this method has good characteristics of controlling their variation. It can also be used to interpolate sparse points on the unit sphere.

In this article, a simple and practical method will be proposed to generate smooth spherical NURBS curves based on De Boor’s algorithm. The structure of this article is organized as follows. In the ‘Definition of spherical NURBS curve based on De Boor’s algorithm’ section, the definition of spherical NURBS curves based on De Boor’s algorithm is proposed. In the ‘PIA of points on a sphere using spherical NURBS curve’ section, a progressive iterative approximation (PIA) scheme is proposed to approximate the given points on unit sphere. In the ‘Case study and discussion’ section, several examples are shown to discuss the continuity of spherical NURBS curve and verify the effectiveness of the proposed method. The ‘Conclusion’ section includes conclusion and discussion.

Definition of spherical NURBS curve based on De Boor’s algorithm

It is supposed that the point on the unit sphere $\mathbf{S}$ can be expressed as a unit vector $\mathbf{P}=(x,y,z)$ $(\left|\mathbf{P}\right|=1)$ . Spherical linear interpolation of two points ${\mathbf{P}}_1$ and ${\mathbf{P}}_2$ on $\mathbf{S}$ can be expressed as follows ¹

$$\mathbf{P}(\alpha )=\frac{\sin [(1-\alpha )\theta ]}{\sin \theta }{\mathbf{P}}_1+\frac{\sin (\alpha \theta )}{\sin \theta }{\mathbf{P}}_2,\alpha \in [0,1]$$

(1)

where $\theta ={\cos }^{-1}({\mathbf{P}}_1•{\mathbf{P}}_2)$ . According to De Casteljau’s algorithm for calculating the Bézier curve, equation (1) can be used to generate spherical Bézier curve, which likes a corner cutting process. De Boor’s algorithm is the generalization of De Casteljau’s algorithm for NURBS. The basis of this method is knot insertion algorithm, which is one of the most important of NURBS algorithms. ⁷ Similar to Shoemake’s method, ¹ De Boor’s algorithm can be extended to the generation of spherical NURBS curves as follows.

Given a set of control points $({\mathbf{P}}_i;w_i)(i=0\dots n)$ on unit sphere $\mathbf{S}$ . $w_i>0$ denotes weights corresponding to points ${\mathbf{P}}_i$ . The total number of control points is $n+1$ . $p$ denotes degree, $s$ denotes multiplicity of inserted knot $t$ , and $U=\{u_0,\dots ,u_p,u_{p+1},\dots ,u_{m-p-1},u_{m-p},\dots ,u_m\}$ denotes knot vector. We can insert a new knot $t$ into the knot vector $U$ . ⁷ Suppose $t$ belongs to the interval $[u_k,u_{k+1})$ , after inserting $t$ with $r$ times, the new control point can be expressed as

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\mathbf{P}}_{i,r}=\frac{\sin \left[(1-{\beta }_{i,r}){\theta }_{i,r-1}\right]}{\sin {\theta }_{i,r-1}}\hfill \\ \multicolumn{1}{c}{}& {\mathbf{P}}_{i-1,r-1}+\frac{\sin \left({\beta }_{i,r}{\theta }_{i,r-1}\right)}{\sin {\theta }_{i,r-1}}{\mathbf{P}}_{i,r-1}\hfill \end{array}$$

(2)

where

$$w_{i,r}=w_{i-1,r-1}(1-{\alpha }_{i,r})+w_{i,r-1}{\alpha }_{i,r}$$

(3)

$$1-{\beta }_{i,r}=\frac{w_{i-1,r-1}(1-{\alpha }_{i,r})}{w_{i,r}}$$

(4)

$${\beta }_{i,r}=\frac{w_{i,r-1}{\alpha }_{i,r}}{w_{i,r}}$$

(5)

$${\theta }_{i,r-1}={\cos }^{-1}({\mathbf{P}}_{i-1,r-1}•{\mathbf{P}}_{i,r-1})$$

(6)

$${\alpha }_{i,r}=\{\begin{array}{cc}\hfill 1\hfill & \hfill i\le k-p+r-1\hfill \\ \multicolumn{1}{c}{\frac{t-u_i}{u_{i+p-r+1}-u_i}}& \hfill k-p+r\le i\le k-s\hfill \\ \multicolumn{1}{c}{0}& \hfill i\ge k-s+1\hfill \end{array}$$

(7)

It should be noted that we still use linear interpolation to get weight $w_{i,r}$ in equation (3), which is the same as NURBS curve in Euclidean space. Comparing equation (2) with equation (1), we know that $\left|{\mathbf{P}}_{i,r}\right|=1$ , which indicates that the new interpolation point ${\mathbf{P}}_{i,r}$ must lie on the unit sphere $\mathbf{S}$ . If $s=0$ and $r=p$ , equations (2) and (3) generate two triangular tables of control points and weights, as shown in Figure 1. It will be helpful to understand the calculation process.

Figure 1.

Two triangular tables of control points and weights: (a) control points and (b) weights.

In addition, suppose ${\mathbf{P}}_i$ denotes a unit quaternion and $\mathbf{S}$ just denotes four-dimensional (4D) unit hypersphere, equations (2)–(7) still work. Considering equation (5), equation (2) also can be rewritten using the multiplication of quaternion as follows

$${\mathbf{P}}_{i,r}={\mathbf{P}}_{i-1,r-1}{\left({\mathbf{P}}_{i-1,r-1}^{-1}{\mathbf{P}}_{i,r-1}\right)}^{{\beta }_{i,r}}$$

(8)

According to the knot insertion algorithm of NURBS, ⁷ after inserting $t$ with $r=p-s$ times, the new point just denotes the point on the NURBS curve with parameter $t$ . Therefore, ${\mathbf{P}}_{i,p-s}$ can be looked as the point on the generated spherical NURBS curve. Suppose $t\in [0,1]$ , the whole spherical NURBS curve can be calculated.

According to the definition, we know that the constructed spherical curve has many good properties, such as local modification property and convex hull property. In addition, when the control points rotate around the centre of the unit sphere, the shape of spherical curve will not change, which can be proved as follows: apply rotation matrix $\mathbf{M}$ to equation (2), we have

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathbf{MP}}_{i,r}=}& \frac{\sin \left[(1-{\beta }_{i,r}){\theta }_{i,r-1}\right]}{\sin {\theta }_{i,r-1}}{\mathbf{MP}}_{i-1,r-1}\hfill \\ \multicolumn{1}{c}{}& +\frac{\sin \left({\beta }_{i,r}{\theta }_{i,r-1}\right)}{\sin {\theta }_{i,r-1}}{\mathbf{MP}}_{i,r-1}\hfill \end{array}$$

(9)

Equation (9) shows that the rotation of control points results in the same rotation of the point on the spherical NURBS curve. Therefore, the shape of curve will not change. It is called as rotation invariant property. Since the generation of spherical NURBS curves is based on knot insertion algorithm in Euclidean space, new control points can be added to increase the flexibility of shape control. But it should be noted that the shape of curve will be changed, which is different from NURBS curve in Euclidean space.

PIA of points on a sphere using spherical NURBS curve

According to equation (2), we know that the spherical NURBS curve is not a linear function with respect to control points for $p\ge 2$ . We cannot interpolate the given points on unit sphere by solving a group of linear equations like interpolation of NURBS curve in Euclidean space. Recently, PIA^21–23 was studied to approximate the given points in Euclidean space, which can be used to interpolate curves and surfaces approximately. It can bring more flexibility for shape control in data fitting. Similarly, a method called geometric algorithm was proposed by Maekawa et al., ²⁴ which can be used to interpolate a triangular or a quadrilateral mesh of arbitrary topological type formed by n vertices. The difference between geometric algorithm and PIA is the registration method of points. PIA uses the assigned parameter of to-be-interpolate points at each iteration. Geometric algorithm gets the parameter of to-be-interpolate points by computing the point-surface distance. Strictly speaking, the two methods are not an interpolation but an approximation since the generated curve does not precisely pass through the given points. But they do not need solving linear system of equations. It is easy to calculate. Following this basic idea, a simple iterative method is developed to approximate the points on sphere using spherical NURBS curves in this section.

The computing scheme is shown in Table 1, which can be described in detail as follows:

• Input: degree $p$ of spherical NURBS curve, points on unit sphere ${\mathbf{Q}}_i$ , weights $w_i$ , tolerance $\epsilon >0$ , maximum iteration number $\mathrm{MINum}$ , and initial error $E_i^{M=0}$ ( $M$ denotes the number of iterations);

• Output: spherical NURBS curve.

1. There are several methods to determine knot vector $U$ . We can choose uniform or non-uniform knot vector. In this article, non-uniform knot vector is calculated from ${\mathbf{Q}}_i$ . First, calculate the parameter $t_i$ of ${\mathbf{Q}}_i$ corresponding to the initial spherical NURBS curve. By connecting ${\mathbf{Q}}_i$ in turn using geodesic arc, a geodesic arc polygon can be obtained. Suppose the total length of the polygon is 1 and $D_i$ denotes the length of geodesic between ${\mathbf{Q}}_i$ and ${\mathbf{Q}}_{i-1}$ , parameter $t_i$ of ${\mathbf{Q}}_i$ can be expressed as

   $$t_i=\{\begin{array}{c}\hfill \frac{\sum _{k=1}^i{D}_k}{\sum _{j=1}^n{D}_i}\begin{array}{cc}\hfill \hfill & \hfill (i=1\dots n)\hfill \end{array}\hfill \\ \multicolumn{1}{c}{0\begin{array}{ccc}\hfill \hfill & \hfill \hfill & \hfill i=0\hfill \end{array}}\end{array}$$

   (10)

   where

   $$D_i=\arccos \frac{{\mathbf{Q}}_i•{\mathbf{Q}}_{i-1}}{\left|{\mathbf{Q}}_i\right|•\left|{\mathbf{Q}}_{i-1}\right|}=\arccos ({\mathbf{Q}}_i•{\mathbf{Q}}_{i-1})(i=1\dots n)$$

   (11)

   Equation (11) just means that $D_i$ is described using the angle between ${\mathbf{Q}}_i$ and ${\mathbf{Q}}_{i-1}$ .

   By using parameter $t_i$ , knot vector $U$ can be obtained ⁷

   $$\{\begin{array}{c}\hfill u_0=\dots =u_p=0\hfill \\ \multicolumn{1}{c}{u_{j+p}=\frac{1}{p}\sum _{i=j}^{j+p-1}{t}_i\begin{array}{cc}\hfill \hfill & \hfill j=1,\dots ,n-p\hfill \end{array}}\\ \multicolumn{1}{c}{u_{m-p}=\dots =u_m=1}\end{array}$$

   (12)

   Suppose ${\mathbf{P}}_i^M={\mathbf{Q}}_i(M=0)$ are the control points, we can calculate the initial spherical NURBS curve ${\mathbf{R}}^M(t)$ using control points ${\mathbf{P}}_i^M$ and knot vector $U$ according to the ‘Definition of spherical NURBS curve based on De Boor’s algorithm’ section. For the convenience of presentation, geodesic arcs on the unit sphere are all expressed as line segments in Figure 2.

2. In order to find the error $E_i^M$ between points ${\mathbf{Q}}_i$ and curve ${\mathbf{R}}^M(t)$ , we need to find the closest point ${\mathbf{R}}^M(t_i^{min})$ to ${\mathbf{Q}}_i$ from ${\mathbf{R}}^M(t)$ within $t\in [u_i,u_{i+p+1})$ . Since the distance can be measured using the length of geodesic arc, the error can be described as the length of geodesic arc between ${\mathbf{R}}^M(t_i^{min})$ and ${\mathbf{Q}}_i$ as follows

   $$E_i^M=\arccos \frac{{\mathbf{Q}}_i•{\mathbf{R}}^M(t_i^{min})}{\left|{\mathbf{Q}}_i\right|•|{\mathbf{R}}^M(t_i^{min})|}=\arccos ({\mathbf{Q}}_i•{\mathbf{R}}^M(t_i^{min}))$$

   (13)

3. We hope to adjust control points ${\mathbf{P}}_i^M$ to ${\mathbf{P}}_i^{M+1}$ and make the new spherical NURBS curve ${\mathbf{R}}^{M+1}(t)$ closer to ${\mathbf{Q}}_i$ than ${\mathbf{R}}^M(t)$ . The new control point ${\mathbf{P}}_i^{M+1}$ can be written as

   $$\begin{array}{cc}\multicolumn{1}{c}{{\mathbf{P}}_i^{M+1}=}& ({\mathbf{Q}}_i\times {\mathbf{R}}^M(t_i^{min}))\times {\mathbf{P}}_i^M\sin \left(E_i^M\right)\hfill \\ \multicolumn{1}{c}{}& +{\mathbf{P}}_i^M\cos \left(E_i^M\right)\hfill \end{array}$$

   (14)

   Obviously, ${\mathbf{P}}_i^{M+1}$ still lies on the unit sphere. Furthermore, since the shape of spherical NURBS curve changes after each iteration calculation, the adjusted direction of ${\mathbf{P}}_i^M$ on the unit sphere also changes.

4. Let the new spherical NURBS curve be looked as the initial spherical NURBS and repeat steps 2 and 3 until $\max (E_i^M)<\epsilon$ , we can obtain the approximation spherical NURBS curve. If $M>\mathrm{MINum}$ , the approximation spherical NURBS curve cannot be found. In this situation, the iteration calculation process may be divergent. We will discuss it using several examples in the ‘Case study and discussion’ section.

Table 1.

Approximation of the points using spherical NURBS curve.

Input:	Degree P	Points on unit sphere Q i maximum iteration number $\mathrm{MINum}$
	Tolerance ε
	Initial error $E_i^{M=0}$
Output:	Spherical NURBS curve
Initialization:	Calculate the initial spherical NURBS curve ${\mathbf{R}}^M(t)$ (M = 0) using Q i as initial control points
If $max(E_i^M)>\epsilon$ and iteration number $M<\mathrm{MINum}$
	Find the nearest point $R^M(t_i^{min})$ of Q i on the spherical NURBS curve
	Calculate errors $E_i^M$ between $R(t_i^{min})$ and Q i
	Adjust control points Q i
	M = M + 1
elseIf $max(E_i^M)<\epsilon$
	Output control points and knot vector of spherical NURBS curve
elseIf $M<\mathrm{MINum}$
Print
(spherical NURBS curve cannot be found)
End

Figure 2.

Approximation of the given points ${\mathbf{Q}}_i$ on unit sphere using spherical NURBS curve: (a) calculation of the initial spherical NURBS curve, (b) calculation of the errors between ${\mathbf{Q}}_i$ and the initial spherical NURBS curve, and (c) moving ${\mathbf{P}}_i^0$ to ${\mathbf{P}}_i^1$ to make the initial spherical NURBS curve approximate ${\mathbf{Q}}_i$ .

As mentioned in step 2, we need to find the closest point ${\mathbf{R}}^M(t_i^{min})$ corresponding to ${\mathbf{Q}}_i$ for defining the error $E_i^M$ . Therefore, we need to know how the closest point behaves. Suppose the distance between ${\mathbf{Q}}_i$ and a point on a sphere curve $\mathbf{R}(t)$ can be written as

$$\theta (t)=\arccos (\mathbf{R}(t),•{\mathbf{Q}}_i)\theta \in [0,\pi )$$

(15)

Suppose $\mathbf{R}(t)$ has at least $C^1$ continuity, differentiating equation (15) with respect to $t$ , we can obtain

$$\frac{d\theta (t)}{\mathrm{dt}}=-\frac{1}{\sqrt{1-{(\mathbf{R}(t)•{\mathbf{Q}}_i)}^2}}\frac{d\mathbf{R}(t)}{\mathrm{dt}}•{\mathbf{Q}}_i$$

(16)

According to extremum condition and equation (16), the closest point meets the following necessary condition

$$\frac{d\theta (t)}{\mathrm{dt}}=0\Rightarrow \frac{d\mathbf{R}(t)}{\mathrm{dt}}•{\mathbf{Q}}_i=0$$

(17)

As shown in Figure 3, equation (17) means that ${\mathbf{Q}}_i$ should be perpendicular to the tangent ${\mathbf{T}}_c$ at the closest point. Obviously, for farthest point, equation (17) also exists. This property is helpful to define the errors in equation (13).

Figure 3.

Extremum condition of distance between point and spherical curve on the unit sphere.

Case study and discussion

In this section, the proposed method is implemented using MATLAB language. Several cases are shown for different purposes.

Case 1

There are some methods to construct spherical curves as mentioned in the ‘Introduction’ section. Usually, we always hope that the spherical NURBS curve inherit all good properties of NURBS curve in Euclidean space. But it is very difficult, even impossible. Compared with NURBS curve in Euclidean space, we can find that the proposed spherical NURBS curve inherits some good properties. In this case, we will show several typical properties with several examples as follows:

1. Local modification property: moving ${\mathbf{P}}_i$ only changes the spherical curve in the interval $[u_i,u_{i+p+1})$ . As shown in Figure 4, $U=\{0,0,0,0,0.25,0.37,1,1,1,1\}$ , $w_i=1$ , moving ${\mathbf{P}}_1$ to ${\mathbf{P}\mathbf{\text{'}}}_1$ only changes the curve in interval $[0,0.25)$ ;

2. Similar to NURBS curve in Euclidean space, ⁶ closed spherical NURBS curve can be constructed as follows: design a uniform knot vector $U=\{u_0=0,u_1=1/m,u_2=2/m,\dots ,u_{m-1}=(m-1)/m,u_m=1\}$ . The curve is defined within $[u_p,u_{n-p}]$ . Then, let ${\mathbf{P}}_0={\mathbf{P}}_{n-p+1}$ , ${\mathbf{P}}_1={\mathbf{P}}_{n-p+2},\dots ,$ ${\mathbf{P}}_{p-2}={\mathbf{P}}_{n-1}$ , and ${\mathbf{P}}_{p-1}={\mathbf{P}}_n$ . Figure 5 shows two examples of closed spherical NURBS curve.

Figure 4.

Local modification property of the spherical NURBS curve.

Figure 5.

Construction of closed spherical NURBS curves.

Case 2

In this case, the continuity of constructed spherical NURBS curve will be discussed. As mentioned in the ‘Definition of spherical NURBS curve based on De Boor’s algorithm’ section, equation (2) just describes a spherical linear interpolation process. For knot multiplicity $s=0$ , the point on the spherical NURBS curve can be obtained after inserting knot $p$ times. Obviously, the spherical NURBS curve is infinitely differentiable between interior knots. But we do not know its continuity at interior knots. Theoretically, we can get the analytical expression of the spherical NURBS curve and prove its continuity at the knot. Unfortunately, its expression is very messy due to the complicated expression of ${\theta }_{i,r-1}$ in equation (6). In this situation, the points are sampled from spherical NURBS curve with equal parameters. Then, the numerical differentiation function in MATLAB will be used to test its parameter continuity for some examples.

For example 1, control points and weights are shown in Table 2. For $p=4$ , the knot vector $U=\{0,0,0,0,0,0.623,1,1,1,1,1\}$ is calculated according to equation (12). The constructed spherical NURBS curve is shown in Figure 6. A total of 1000 points are sampled from spherical NURBS curve. The first- to fourth-order approximate derivatives are shown in Figure 7. We can find that the fourth-order approximate derivatives are not continuous at knot 0.6225 obviously.

Table 2.

Control points and weights for example 1.

	x	y	z	$w_i$
${\mathbf{Q}}_0$	0	1	0	1
${\mathbf{Q}}_1$	0.7387	0.3764	0.5592	1
${\mathbf{Q}}_2$	0	0	1	1
${\mathbf{Q}}_3$	0	0.7071	0.7071	1
${\mathbf{Q}}_4$	−0.4111	0.7417	0.5299	1
${\mathbf{Q}}_5$	−0.6076	0.3368	0.7193	1

Figure 6.

Spherical NURBS curve for example 1.

Figure 7.

First- to fourth-order approximate derivatives of spherical NURBS curve for example 1: (a) first-order approximate derivatives, (b) second-order approximate derivatives, (c) third-order approximate derivatives, and (d) fourth-order approximate derivatives (number of sampling points: 1000, blue curve: x component, green curve: y component, and red curve: z component).

For example 2, control points and weights are shown in Table 3. For $p=3$ , the calculated knot vector is $U=\{0,0,0,0,0.533,0.723,1,1,1,1\}$ . Two curves with different weights are shown in Figure 8. First- and second-order approximate derivatives are still smooth at two interior knots $u_4=0.533$ and $u_5=0.723$ . Figure 9 shows the discontinuity of third-order approximate derivatives of the spherical NURBS curve.

Table 3.

Control points and weights for example 2.

	x	y	z	$w_i$
${\mathbf{Q}}_0$	0.9992	−0.0405	0	1	1
${\mathbf{Q}}_1$	−0.9952	0.0497	0.0837	1	10
${\mathbf{Q}}_2$	0.8885	−0.3059	−0.3420	1	10
${\mathbf{Q}}_3$	0.9803	−0.1204	−0.1564	1	10
${\mathbf{Q}}_4$	−0.5050	0.5609	0.6561	1	10
${\mathbf{Q}}_5$	−0.9659	0.2414	−0.0926	1	1

Figure 8.

Spherical NURBS curve for example 2: (a) $w_i=1$ and (b) $w_0=w_5=1$ , $w_i=10(i=1-4)$ .

Figure 9.

Third-order approximate derivatives are not continuous at two interior knots for example 2 (blue curve: x component, green curve: y component, and red curve: z component).

It seems that this spherical NURBS curve has $p-1$ continuity at interior knots. But there is still another possibility: the discontinuity at interior knot may be very small that it is not detected by numerical procedures. Furthermore, the number of sampling points is changed to 20,000. The first- and second-order approximate derivatives are still smooth, but the third-order approximate derivatives have a small break at knot 0.623, as shown in Figure 10. It shows that the spherical NURBS curve may be not $C^{p-1}$ continuous. But the break is very small and difficult to detect. In order to study this property further, three additional tests are conducted. All spherical NURBS curves will be sampled into 20,000 points to calculate approximate derivatives. For the convenience of the study, the relative error of the continuity at interior knots R_E_C is defined to describe the degree of continuity.

$$\mathrm{R}\_\mathrm{E}\_\mathrm{C}=\frac{\mathrm{h}}{\mathrm{H}}\times 100\%$$

Figure 10.

Third-order approximate derivatives of spherical NURBS curve for example 1 at knot 0.623 (number of sampling points: 20,000; blue curve: x component, green curve: y component, and red curve: z component).

As shown in Figure 11, $h_i$ denotes the difference of the derivatives between two adjacent curves at interior knots. H denotes the variation range of the derivatives for the whole spherical curve. Let $\mathrm{D}\_\mathrm{N}$ denote the derivative order, we can calculate $\mathrm{R}\_\mathrm{E}\_\mathrm{C}$ for the three examples. Results are shown in Appendix 1, from which we know that (1) all examples seem to be $C^1$ continuous; (2) when $\mathrm{D}\_\mathrm{N}\le p-1$ , $\mathrm{R}\_\mathrm{E}\_\mathrm{C}$ is very small and usually less than 1%; (3) when $\mathrm{D}\_\mathrm{N}=p$ , $\mathrm{R}\_\mathrm{E}\_\mathrm{C}$ increases obviously; and (4) since the number of sampling points is huge and the difference between adjacent sampling points is very tiny, when calculating the fourth-order derivatives, numerical calculation is unstable for two of the three examples. These tests show that the spherical NURBS curve may have $C^1$ continuity and have no $C^{p-1}$ continuity for $p\ge 3$ .

Figure 11.

Definition of relative error of the continuity at interior knots.

Case 3

In this case, the approximation of points on sphere using spherical NURBS curve and its convergence will be discussed. Two examples are shown in Figures 12 and 13. Viviani’s curve is a sphere–cylinder intersection. In this example, nine points are sampled from Viviani’s curve on the unit sphere, which is written as

$$\{\begin{array}{c}\hfill \mathrm{x}=0.5(1+\mathrm{cost})\hfill \\ \multicolumn{1}{c}{\mathrm{y}=0.5\mathrm{sint}}\\ \multicolumn{1}{c}{\mathrm{z}=\sin (0.5\mathrm{t})}\end{array}\mathrm{t}\in [-2\pi ,2\pi )$$

Figure 12.

Approximation of nine sampled points from Viviani’s curve: (a) initial spherical NURBS curve defined using sampled point from Viviani’s curve and (b) $M=10$ , $max(E_i^{10})=0.00087$ .

Figure 13.

Approximation of 30 sampled points from a sphere spiral: (a) the initial spherical NURBS curve and (b) $p=3$ , $M=6$ , $max(E_i^6)=0.00080$ .

Let sampling points be the initial control points; for $p=4$ , the initial spherical NURBS curve is generated in Figure 12(a). When the number of iterations $M=10$ , the maximum error $E_i^{10}=0.00087$ . The approximation of spherical NURBS curve is shown in Figure 12(b). Similarly, 30 points are sampled from a sphere spiral, which are written as

$$\{\begin{array}{c}\hfill x=\cos (-\frac{\pi }{2}+0.05t)\cos (t)\hfill \\ \multicolumn{1}{c}{y=\cos (-\frac{\pi }{2}+0.05t)\sin (t)}\\ \multicolumn{1}{c}{z=\sin (-\frac{\pi }{2}+0.05t)}\end{array}\mathrm{t}\in [0,29]$$

For $p=3$ , these points are approximated, as shown in Figure 13. When the number of iterations $M=6$ , the maximum error $E_i^6=0.00080$ . From the two examples, we know that the proposed method can approximate the given points on the unit sphere effectively.

However, it is still needed to note that the two cases have different convergence rates. Since the approximation process is iterative, it is necessary to discuss its convergence further. As shown in Figure 14, 30 sampled points from sphere spiral are approximated with different degree $p(p=5,6,10)$ . We can find that the convergence rate will decrease with the increase in $p$ . As shown in Figure 15, degree $p=2-5$ , points ${\mathbf{Q}}_i$ , and weights $w_i$ are shown in Table 2. When $M=1-20$ and $p=5$ , the calculation process is divergent. Form Figure 15(b), we can see that the spherical NURBS curve is deformed much. The two examples show that the convergence is affected by distribution of the points greatly. Since this process is complex, we cannot prove strict convergence conditions mathematically at present. But through many tests, we can summarize some of the practical convergence conditions: (1) the angle ${\theta }_{i,r-1}$ is within $[0,\pi /2)$ ; (2) the angle ${\delta }_i$ between geodesic arcs at ${\mathbf{Q}}_i$ is within $[\pi /2,\pi ]$ , as shown in Figure 6. Furthermore, if possible, smaller $p$ should be used to reduce the calculation time and the possibility of divergence. Usually, $p\le 4$ is enough in most situations.

Figure 14.

Analysis of convergence for approximating sampling points from sphere spiral: (a) $p$  = 10 M = 20 and (b) relationship between the number of iterations $M$ and $max(E_i^M)$ for different ps.

Figure 15.

Approximation of the points in Table 2 using different ps: (a) $p=3$ , $M=20$ ; (b) $p=5$ , $M=20$ ; and (c) relationship between the number of iterations $M$ and $max(E_i^M)$ .

Case 4

Smoothness of tool path is very important for dynamic performance of the machine tool. ²⁵ To overcome the common weakness of dramatic change of tool orientations, a new method has been developed for smooth tool motions by searching the alternative tool orientation in the C-space with the smallest change of tool orientations. ²⁶ But this method cannot smooth the tool path globally. Langeron et al. ²⁷ proposed a new format of tool path polynomial interpolation in five-axis machining using two B-spline curves, which can be used to avoid five-axis singularities ²⁸ and keep the tool path smooth at the same time. The difference between this method and our method is that the tool axis orientation is described as a spherical NURBS curve in this article, which makes it easy to see the variation of tool axis orientation on a sphere. In this case, the proposed method is used to generate smooth tool path without interference in five-axis machining. As shown in Figure 16(a), six typical tool positions and tool axis directions are defined. In computer-aided manufacturing (CAM) system, they are usually defined using man–machine interaction and free of interference. Furthermore, we need to interpolate the six tool axes and get the whole tool path. Since the six directions can be described as six points on a unit sphere, the proposed method is used to approximate these points and get the directions of the whole tool path. The trajectory of the tip of the tool can be interpolated using NURBS curve in Euclidean space. It is obvious that interference happens between the directions of the tool axes and the obstacle. To avoid interference, tool axis A is chosen and adjusted parallel to the observing direction. A new smooth tool path is obtained, as shown in Figure 16(b). Obviously, the interference is eliminated. In addition, the adjustment of tool axis A only affects the initial tool path locally. Thus, this method can be used to modify five-axis tool path locally.

Figure 16.

Avoiding interference by adjusting fewer directions and obtaining smooth tool path interference in five-axis machining: (a) tool path with interference and (b) tool path without interference.

Conclusion

A novel and simple method of generating spherical NURBS curve and its application have been studied in this article. The definition of spherical NURBS curves is based on De Boor’s algorithm. We can get the following conclusions:

1. The generated spherical NURBS curve preserves some good properties from NURBS curve in Euclidean space, such as local modification property, convex hull property, rotation invariant property, and knot insertion property. But we need to know that knot insertion will change the shape of the curve. It is easy to construct closed spherical NURBS curve using the proposed method. This method will be very useful to design spherical curve interactively.

2. For the sequence of the given points on the unit sphere, a PIA scheme has been proposed to approximate them using smooth spherical NURBS curves. Results show that the position relationship of points will affect the convergence. A practical but not strict convergence condition for approximation is presented.

3. Since the analytical expression of spherical NURBS curves is messy, numerical differentiation is used to analyse its continuity at interior knots. Results show that spherical NURBS curves may have $C^1$ continuity for $p\ge 2$ and do not have $C^{p-1}$ continuity at interior knots like NURBS curve in Euclidean space. But when the derivative order is less than $p$ , the difference of derivatives between two adjacent curves at the interior knots is very small. When derivative order is equal to $p$ , the difference of derivatives between two adjacent curves at the interior knots becomes obvious.

4. The proposed method is applied on tool path generation of five-axis machining to avoid interference by adjusting fewer directions and obtain smooth tool path simultaneously.