Equiareal parameterizations of NURBS surfaces

Graphical Models 76 (2014) 43–55

Contents lists available at ScienceDirect

Graphical Models journal homepage: www.elsevier.com/locate/gmod

Equiareal parameterizations of NURBS surfaces q Yi-Jun Yang a,⇑, Wei Zeng b, Jian-Feng Chen c a b c

School of Computer Science and Technology, Shandong University, Jinan, China School of Computing & Information Sciences, Florida International University, USA School of Software, Shandong University, Jinan, China

a r t i c l e

i n f o

Article history: Received 30 May 2013 Received in revised form 20 October 2013 Accepted 24 October 2013 Available online 14 November 2013 Keywords: NURBS surfaces Möbius transformations Equiareality

a b s t r a c t The equiareality of NURBS surfaces greatly affects the results of visualization and tessellation applications, especially when dealing with extruding and intruding shapes. To improve the equiareality of given NURBS surfaces, an optimization algorithm using the Möbius transformations is presented in this paper. The optimal Möbius transformation is obtained by computing the intersection of two planar algebraic curves, whose coefficients are computed explicitly for Bézier and B-spline surfaces, while numerically for NURBS surfaces. Examples are given to show the performance of our algorithm for visualization and tessellation applications. Ó 2013 Elsevier Inc. All rights reserved.

1. Introduction NURBS surfaces play an increasingly important role in contemporary Computer Aided Design (CAD). The results of most surface algorithms such as surface registration, surface visualization, surface tessellation, and so on [24–26,30,31], highly depend on the surface parameterization. A NURBS surface has an intrinsic rational polynomial mapping (see Fig. 1) from the 3D surface to the 2D parameter domain (a unit square). By surface reparameterizations [12,24–26], the surface may have infinitely many different parameterizations. Depending on where and how it will used, one may need to find a suitable or optimal parameterization out of the infinitely many, or to convert the given parameterization into another (more) suitable one [12,24–26]. In many applications, such as the quantitative analysis and visualization of brain, joint bone and mechanical surfaces with rich extruding and intruding fea-

q This paper has been recommended for acceptance by Alla Sheffer and Peter Lindstrom. ⇑ Corresponding author. Address: 1500 Shunhua Road, Gaoxin District, Shandong University, Jinan, Shandong 250100, China. E-mail address: [email protected] (Y.-J. Yang).

1524-0703/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.gmod.2013.10.007

tures, it is highly desirable that the parameterization can preserve area elements of the 3D surfaces, so that many area-related patterns, e.g., intestinal density, activation extent, thickness, etc., can be better represented in the parameter domain [31]. Equiareality is a natural and direct 2D surface extension of arc-length parameterization for curve cases [2,7,8,11,15,18–23,29]. Moreover, an equiareal parameterization will lead to more robust and stable computations for derivative based algorithms such as surface intersection, curvature computation, and so on [1,14]. 1.1. Related works We focus on review on NURBS curve and surface parameterization methods. For triangle meshes and subdivision surfaces, we concentrate on analyzing the difference between their parameterization and freeform surface parameterization, a detailed review of mesh and subdivision surface parameterization techniques being beyond the scope of this article. For the details of mesh parameterization, the reader is recommended to see the survey paper by Floater and Hormann [8] and the references therein. For the details of subdivision surface parameterization, see the papers by He et al. [9] for the latest progress.

44

Y.-J. Yang et al. / Graphical Models 76 (2014) 43–55

Fig. 1. The intrinsic mapping of a NURBS surface.

In the past twenty years, how to achieve optimal parameterization of Bézier curves has been studied extensively in the literature such as [2,7,11,15,18–23,29]. Farouki [7] identified arc-length parameterization as the optimal parameterization of Bézier curves. By minimizing an integral which measures the deviation from arc-length parameterization, the optimal representation is obtained by solving a quadratic equation. Jüttler [11] presented a simplified approach to Farouki’s result by using a back substitution in the integral. Costantini et al. [2] obtained closer approximations to the arc-length parameterization by applying composite reparameterizations to Bézier curves. Farin [6] showed that for a circular arc represented by its rational quadratic representation, its parameterization is a chord length one. Motivated by Farin’s work, Lu [13] identified a family of curves that can be parameterized by chord length. The parameterization of triangle meshes has been studied extensively in the last decade and still remains as a hot topic until now [8,16,31]. The main purpose of the research on the parameterization of triangle meshes is to construct a suitable, bijective mapping between the triangle mesh embedded in 3D and a simple 2D domain, referred to as the parameter space or parameter domain. To minimize the parameterization distortion in either angles or areas, many different algorithms have been proposed in the literature [8,16,31]. Also He et al. presented algorithms for parameterizing subdivision surfaces in [9]. As the NURBS surface already has a rational polynomial mapping (see Fig. 1) from the 2D parameter domain (a unit square) to the 3D surface, its parameterization has some specific properties different from the parameterization of triangle meshes and subdivision surfaces. First the NURBS surface has an intrinsic mapping already and we do not need to construct an initial surface mapping from the 3D surface to the 2D parameter domain, which is the case for triangle meshes and subdivision surfaces. Second the parameterization of NURBS surface is a continuous rational polynomial mapping while those of triangle meshes and subdivision surfaces are discrete, usually defined by the correspondence between their vertices and the correspondence of points inside the triangles/quads is obtained from the vertices correspondence by interpolation techniques. If we convert the NURBS surface into a triangle mesh and apply the mesh parameterization method, some parameterization results can be obtained subsequently. However, there is one main drawback for this kind of methods. The resultant surface representation is not NURBS anymore,

which is problematic for subsequent CAD algorithms designed for freeform surfaces. Though we can reconstruct the NURBS surface by traditional least-square fitting methods from the triangle parameterization, neither the surface shape nor the triangle parameterization are preserved precisely during the fitting procedure, which is not allowed for high accurate CAD applications. While most of the successes have been reported about the parameterization of NURBS curve, triangle mesh and subdivision surface parameterization [2,7–11,15,16,18– 23,29,31], the NURBS surface parameterization has not met with similar achievements. The results of rendering and tessellation applications for NURBS surfaces largely depend on the parameterization quality. He et al. [10] gave a rational bicubic reparameterization method to improve the parameterization of the approximate Gregory patches such that the new parameterization conforms better to that of the given subdivision surface. Both the explicit representation of the reparameterized surface and the equiareality of the final surface are not considered therein. Yang et al. [24] presented an algorithm to optimize the uniformity of iso-parameter curves for Bézier surfaces based on Möbius transformations [5,24], which can change only the distribution of iso-parameter curves but not the shape of them. However, the above method [24] minimize the uniformity only on sampled iso-parameter curves, not on the whole surface. Furthermore, the equiareality energy describes the uniformity of the surface stretch better than the formulated energies utilized in [24]. To the author’s knowledge, there is no equiareality related work about NURBS surfaces until now. The aim of this paper is to optimize the equiareality of given NURBS surfaces for CAD applications such as surface visualization, surface tessellation, surface intersection, curvature computation and so on. Most of all, we want to keep the surface geometry (the surface shape) unchanged and let the resultant surface represented as a NURBS surface, which is preferable and convenient for the algorithms designed for the above mentioned Computer Aided Design applications. The method presented in this paper utilizes the Möbius transformation, which satisfies both the surface shape and representation requirements and serves as a startup research to optimize the equiareality of given NURBS surfaces. Due to the ability of Möbius transformations to change surface parameterizations, there is still a lot of room to improve, which is left as a future work. 1.2. Algorithm overview The input to our algorithm is a NURBS surface and our goal is to find an optimal Möbius transformation to optimize its surface equiareality. A NURBS surface in three dimensional space is defined by

Pnu 1 Pnv 1

p q i¼0 j¼0 N i ðuÞN j ðv Þxi;j Pi;j Xðu; v Þ ¼ Pnu 1 Pnv 1 p ; q i¼0 j¼0 N i ðuÞN j ðv Þxi;j

u; v 2 ½0; 1; ð1Þ

where Pi;j are the control points, xi;j are the weights, and N pi ðuÞ, N qj ðv Þ are the pth-degree and qth-degree B-spline basis functions defined on the knot vectors

45

Y.-J. Yang et al. / Graphical Models 76 (2014) 43–55

8 9 > > < = U ¼ 0; . . . ; 0; upþ1 ; . . . ; unu 1 ; 1; . . . ; 1 ; |fflfflfflffl{zfflfflfflffl}> > :|fflfflfflffl{zfflfflfflffl} ;

where E and G give the square length of the two partial derivatives and F measures the orthogonality of the two partial derivatives. Then, we have

and

ds ¼ ðdu dv Þ I

pþ1

pþ1

2

8 9 > > < = V ¼ 0; . . . ; 0; v qþ1 ; . . . ; v nv 1 ; 1; . . . ; 1 ; |fflfflfflffl{zfflfflfflffl} |fflfflfflffl{zfflfflfflffl} > > : ; qþ1

qþ1

1. Compute the coefficients of two planar curves which correspond to two planar equations obtained by setting the partial derivatives of formulated energy to be zero (see Section 2.4). 2. Obtain the intersection of the two planar curves, whose coordinates correspond to the optimal Möbius transformation parameters. 3. Get the optimal equiareal NURBS surface by reparameterizing given NURBS surface using the obtained Möbius transformation. The paper is organized as follows. Section 2 describes how to minimize the equiareality of given NURBS surfaces by Möbius transformations. In Section 3, several examples are given to show the performance of our algorithm and Section 4 concludes the paper. 2. Area preserving NURBS surfaces

The parameterization in Eq. (1) is characterized by its first fundamental form [4]

where Xu ¼ @X and Xv ¼ @X are the two partial derivative @u @v vectors of the surface X. The first fundamental form describes the metric on the surface X. Let

G ¼ Xv  Xv ;

and rewrite the coefficients in a symmetric matrix

 I¼

E

F

F

G

 ;

the discriminant of the quadratic form. It can be easily verified that

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi EG  F 2 ¼ kXu  Xv k;

u; v 2 ½0; 1;

always holds across the whole surface. Let A denote the area of the given surface,

A¼

Z

1

Z

1

Z

0

¼

Z

0

1

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k EG  F 2 kdudv

0 1

0

kXu  Xv kdudv :

ð2Þ

As the area of the parameter domain (a unit square) is always equal to 1, it is impossible to obtain an area-preserving parameterization for a surface with larger or smaller surface areas (i.e., A – 1). Thus we relax the equiareality to uniform area stretch across the freeform surfaces, i.e., we call surface X area preserving if

gðu; v Þ  A2 ;

u; v 2 ½0; 1;

ð3Þ

holds. Condition (3) can also be rewritten as follows

kXu  Xv k  A;

u; v 2 ½0; 1:

Without explicit notation, in the remainder of this paper, equiareality always means that the surface is uniformly stretched across the whole freeform surface. In the next subsection, we describe the Möbius transformations, based on which the equiareality of given freeform surfaces is optimized.

The Möbius transformations of freeform surfaces described in Eq. (1) are as follows.

ða  1Þs ; 2as  s  a

a 2 ð0; 1Þ;

ð4Þ

ðb  1Þt ; 2bt  t  b

b 2 ð0; 1Þ:

ð5Þ

and

v ¼ f2 ðtÞ ¼

2

ds ¼ Xu  Xu ðduÞ þ 2Xu  Xv dudv þ Xv  Xv ðdv Þ ;

F ¼ Xu  Xv ;

g ¼ det I ¼ EG  F 2 ;

u ¼ f1 ðsÞ ¼

2.1. Differential geometry of NURBS surfaces

E ¼ Xu  Xu ;

dv

 :

2.2. Möbius transformations of NURBS surfaces

In this section, we first study the differential geometry of NURBS surfaces and then try to optimize the equiareality of given NURBS surfaces using Möbius transformations.

2

du

Under the assumption of regularity, the matrix has a strictly positive determinant

respectively. In the optimization procedure, a nonlinear energy measuring the equiareality deviations for NURBS surfaces is formulated. By setting the two partial derivatives of the formulated energy to be zero, the optimal Möbius transformation can be obtained by computing the intersection of two corresponding planar algebraic curves. The coefficients of the two planar curves are explicitly computed for Bézier surface and B-spline surfaces. In addition, the algorithm is generalized to NURBS surfaces by computing the curve coefficients using a numerical quadrature method. The main algorithm flow is described as follows.

2



For the two parameters, the two transformations map the unit intervals onto themselves, respectively. In particular, they satisfy



f1 ð0Þ ¼ 0; f1 ð1Þ ¼ 1 and f 1 ð0:5Þ ¼ 1  a f2 ð0Þ ¼ 0; f2 ð1Þ ¼ 1 and f 2 ð0:5Þ ¼ 1  b

:

The inverse transformations of Eqs. (4) and (5) evaluate to

s ¼ f3 ðuÞ ¼

au ; 1  a þ 2au  u

46

Y.-J. Yang et al. / Graphical Models 76 (2014) 43–55

and

becomes as small as possible. The integral measures the deviation of the current surface parameterization Xðs; tÞ from its equiareal parameterizations. From Eqs. (6) and (7), we have

bv t ¼ f4 ðv Þ ¼ : 1  b þ 2bv  v Applying the transformations (4) and (5) to surface (1) results in the reparameterized surface [12]

Pnu 1 Pnv 1

p q ~ i¼0 j¼0 N i ðsÞN j ðtÞ i;j Pi;j Pnu 1 Pnv 1 p q ~ i;j N ðsÞN ðtÞ i¼0 j¼0 i j

Xðs; tÞ ¼

x x

;

s; t 2 ½0; 1;

Qn

0

Z

ðkXs  Xt k  AÞ2 dsdt

0

Z

1

Z

1

kXs  Xt k2 dsdt  2A

Z

0

0

1

Z

1

kXs

0

 Xt kdsdt þ A2 Z 1Z 1 ¼ kXs  Xt k2 dsdt  2A2 þ A2

v

l¼1 ð1  b þ 2b jþl

 v jþl Þ

0

¼

Z

0

1

Z

0

and the two new knot vectors

8 > < aupþ1 e ¼ 0; . . . ; 0; U ;...; |fflfflfflffl{zfflfflfflffl} > 1  a þ 2 aupþ1  upþ1 : pþ1 9 > = aunu 1 ; 1; . . . ; 1 ; 1  a þ 2aunu 1  unu 1 |fflfflfflffl{zfflfflfflffl}> ;

1

kXs  Xt k2 dsdt  A2

ð8Þ

0

As an immediate consequence from Eq. (8), this is equivalent to finding the minimum of the integral function

Jða; bÞ ¼

Z

1

0

Z

1

kXs  Xt k2 dsdt:

ð9Þ

0

Let

pþ1

ci;j ¼

and 8

e¼ V

1

0

xi;j

k¼1 ð1  a þ 2auiþk  uiþk Þ

1

¼

with the same control points, the new weights

~ i;j ¼ Qm x

Z

Z

1

Z

0

> < bv qþ1 0; . . . ; 0; ;...; |fflfflfflffl{zfflfflfflffl} > 1  b þ 2b v qþ1  v qþ1 : qþ1 9 > = bv nv 1 ; 1; . . . ; 1 : 1  b þ 2bv nv 1  v nv 1 |fflfflfflffl{zfflfflfflffl}> ;

1

0

kXu  Xv k2 B2i ðuÞB2j ðv Þdudv ;

i ¼ 0; 1; 2;

j ¼ 0; 1; 2;

ð10Þ

where B2i ðuÞ and B2j ðv Þ are quadratic Bernstein basis of variables u and v, respectively. The integral term in Eq. (10) is always non-negative for points on the surface. Thus the coefficients ci;j are positive for any freeform surface. By the chain rule, we obtain from Eq. (9)

qþ1

The Möbius transformations only change the magnitude of the partial derivative vectors, not the direction of them. Thus the Möbius transformations keep the angle between the u=v partial derivative vectors as well as the orthogonality of iso-parameter curves. After describing the Möbius transformations, we will present a method to optimize the equiareality energy of given freeform surfaces based on the Möbius transformations in the next subsection.

Jða; bÞ ¼

Z

1

0

¼

Z

¼

1

Z

1

Z

Z

1

0

1

0

2.3. Optimization method of equiareality energy

1

0

0

¼

1

kXs  Xt k2 dsdt

0

0

Z

Z

Z

0

1

kXu  Xv k2 

 2  2 @u @v  dsdt @s @t

kXu  Xv k2 

@u @ v  dudv @s @t

kXu  Xv k2 

ð1  a þ 2ua  uÞ2 að1  aÞ

ð11Þ

ð1  b þ 2v b  v Þ2 dudv bð1  bÞ 2 2 X 2 X ci;j Bi ðaÞB2j ðbÞ    ¼ 2 2 i¼0 j¼0 að1  aÞbð1  bÞ j i 

Given a NURBS surface, here we try to optimize its equiareality using Möbius transformations. The Möbius transformations only change the surface parameterization, but not the shape of the surface. Thus the surface area is preserved during the Möbius transformation, which leads to the following identity equation

A

Z

1

0

Z

1

kXs  Xt kdsdt;

ð6Þ

0

1

0

0¼

Z

1

0¼

@Jða; bÞ ; @a

ð12Þ

@Jða; bÞ : @b

ð13Þ

and

which is the area of the surface. Here we want to choose the Möbius transformations (4) and (5) such that the value of

Z

The solution satisfies the following two equations

2

ðkXs  Xt k  AÞ dsdt

ð7Þ

0¼

From Eqs. (11) and (12), we have

0

@Jða; bÞ B22 ðaÞðc2;0 B20 ðbÞ þ c2;1 B21 ðbÞ=2 þ c2;2 B22 ðbÞÞ  B20 ðaÞðc0;0 B20 ðbÞ þ c0;1 B21 ðbÞ=2 þ c0;2 B22 ðbÞÞ ¼ : @a a2 ð1  aÞ2 bð1  bÞ

ð14Þ

47

Y.-J. Yang et al. / Graphical Models 76 (2014) 43–55

To satisfy Condition (12), the numerator of Eq. (15) must be set to zero, which leads to

B22 ðaÞ B20 ðaÞ

¼

c0;0 B20 ðbÞ þ c0;1 B21 ðbÞ=2 þ c0;2 B22 ðbÞ c2;0 B20 ðbÞ þ c2;1 B21 ðbÞ=2 þ c2;2 B22 ðbÞ

;

ð15Þ

where B2i ðaÞ and B2i ðbÞ are the Bernstein polynomial basis of a and b, respectively. Similarly, from Eq. (13), we can obtain

B22 ðbÞ B20 ðbÞ

c0;0 B20 ðaÞ þ c1;0 B21 ðaÞ=2 þ c2;0 B22 ðaÞ

¼

c0;2 B20 ðaÞ þ c1;2 B21 ðaÞ=2 þ c2;2 B22 ðaÞ

:

ð16Þ

The curves expressed by Eqs. (15) and (16) are planar algebraic curves. We first examine the planar curve Y 1 expressed by Eq. (15). Substituting b0 ¼ 0 and b1 ¼ 1 into Eq. (15), we have

0 < a0 ¼

1þ

1 qffiffiffiffiffiffi and a1 ¼ c 2;0

c0;0

1þ

1 qffiffiffiffiffiffi < 1: c2;2 c0;2

ð17Þ

The right term of Eq. (15) can be expressed as a degree-two Bézier form as follows.

R¼ ¼

c0;0 B20 ðbÞ þ c0;1 B21 ðbÞ=2 þ c0;2 B22 ðbÞ c2;0 B20 ðbÞ þ c2;1 B21 ðbÞ=2 þ c2;2 B22 ðbÞ P0 x0 B20 ðbÞ þ P1 x1 B21 ðbÞ þ P2 x2 B22 ðbÞ

x0 B20 ðbÞ þ x1 B21 ðbÞ þ x2 B22 ðbÞ

2.4.1. Coefficients computation for Bézier surfaces A Bézier surface can be represented in the following form

8 c P ¼ 0;0 > > > 0 c2;0 > > c > > P1 ¼ c0;1 > > 2;1 > < c P2 ¼ c0;2 : 2;2 > > x0 ¼ c2;0 > > > > > x1 ¼ c2;1 =2 > > > : x2 ¼ c2;2

Xðu; v Þ ¼

u 2 ½0; 1 and

v

2 ½0; 1;

ð20Þ

n where Pi;j are the control points, Bm i ðuÞ and Bj ðv Þ are the Bernstein polynomials. Its u=v partial derivatives are as follows.

m 1X n X

1 pffiffiffiffiffiffiffiffi < 1; 1 þ 1=R

ð19Þ

The shape of the curve Y 1 expressed by Eq. (15) is illustrated in Fig. 2. Similar statements hold for the curve Y 2 , which passes points ð0; b2 Þ and ð1; b3 Þ where 0;2

c0;0

i¼0 j¼0

ð21Þ and

minðP0 ; P1 ; P2 Þ 6 R 6 maxðP0 ; P1 ; P2 Þ:

1 qffiffiffiffiffiffi and b3 ¼ c

Bm1 ðuÞBnj ðv ÞQ i;j ; u 2 ½0;1 and v 2 ½0; 1; i

Xu ðu; v Þ ¼ m

ð18Þ

where

1þ

m X n X n Bm i ðuÞBj ðv ÞPi;j ; i¼0 j¼0

For each 0 6 b 6 1, there exists a unique a inside ð0; 1Þ

0 < b2 ¼

The intersection points can be obtained robustly using the algorithm presented by Sederberg [17]. The main limitation of this paper is that there is no guarantee that there is only one intersection point for the two planar curves, i.e., energy (11) is convex. However, in practice, we observe that the algorithm always has a unique solution (only one intersection point), as discussed in more detail in Section 4. 2.4. Computation of coefficients

;

where

0
Fig. 2. The illustration of the two planar curves Y 1 ; Y 2 and their intersection.

1þ

1 qffiffiffiffiffiffi < 1: c2;2 c2;0

Xv ðu; v Þ ¼ n

m X n1 X n1 Bm ðv ÞRi;j ; u 2 ½0;1 and v 2 ½0;1; i ðuÞBj i¼0 j¼0

ð22Þ

where



Q i;j ¼ Piþ1;j  Pi;j Ri;j ¼ Pi;jþ1  Pi;j

:

From Eqs. (21) and (22), the first fundamental form terms are computed as follows.

8 E ¼ Xu  Xu ¼ > > Pm1 Pn Pm1 Pn ð2m2ikÞ!ðiþkÞ!ð2njlÞ!ðjþlÞ! 2m2 > 2n > m!m!n!n! > > i¼0 j¼0 k¼0 l¼0 ðm1iÞ!i!ðm1kÞ!k!ðnjÞ!j!ðnlÞ!l! Biþk ðuÞBjþl ðv ÞQ i;j  Q k;l ð2m2Þ!ð2nÞ! > > > < G ¼ Xv  Xv ¼ Pm Pn1 Pm Pn1 ð2mikÞ!ðiþkÞ!ð2n2jlÞ!ðjþlÞ! 2m m!m!n!n! B ðuÞB2n2 > jþl ðv ÞR i;j  R k;l > > ð2n2Þ!ð2mÞ! i¼0 j¼0 k¼0 l¼0 ðn1jÞ!j!ðn1lÞ!l!ðmiÞ!i!ðmkÞ!k! iþk > > > ¼ > > F ¼ Xu  Xv P > m1 Pn Pm Pn1 ð2m1ikÞ!ðiþkÞ!ð2n1jlÞ!ðjþlÞ! 2m1 : m!m!n!n! 2n1 i¼0 j¼0 k¼0 l¼0 ðm1iÞ!i!ðmkÞ!k!ðnjÞ!j!ðn1lÞ!l! Biþk ðuÞBjþl ðv ÞQ i;j  R k;l ð2m1Þ!ð2n1Þ!

ð23Þ

48

Y.-J. Yang et al. / Graphical Models 76 (2014) 43–55

Then, the coefficients ci;j in Eq. (10) can be obtained as follows

Z

1

Z

ci;j ¼

X ncr;i;j ; r¼1

1

kXu  Xv k2 B2i ðuÞB2j ðv Þdudv    Z 1Z 1 2 m!m!n!n!m!m!n!n! Xm1 Xn Xm1 Xn Xm Xn1 Xm Xn1 2 ¼ ðEG  F 2 ÞB2i ðuÞB2j ðv Þdudv ¼ i1 ¼0 j1 ¼0 k1 ¼0 l1 ¼0 i2 ¼0 j2 ¼0 k2 ¼0 l2 ¼0 j ð4m þ 1Þ!ð4n þ 1Þ! i 0 0 ð4m  i1  i2  k1  k2  iÞ!ði1 þ i2 þ k1 þ k2 þ iÞ! ð4n  j1  j2  l1  l2  jÞ!ðj1 þ j2 þ l1 þ l2 þ jÞ! ðm  1  i1 Þ!i1 !ðm  1  k1 Þ!k1 !ðn  j1 Þ!j1 !ðn  l1 Þ!l1 ! ðn  1  j2 Þ!j2 !ðn  1  l2 Þ!l2 !ðm  i2 Þ!i2 !ðm  k2 Þ!k2 !

ci;j ¼

0

0

ðQ i1 ;j  Q k1 ;l1 ÞðRi2 ;j2  Rk2 ;l2 Þ  1  2 2 m!m!n!n!m!m!n!n! Xm1 Xn Xm Xn1 Xm1 Xn Xm Xn1 i1 ¼0 j1 ¼0 k1 ¼0 l1 ¼0 i2 ¼0 j2 ¼0 k2 ¼0 l2 ¼0 i j ð4m þ 1Þ!ð4n þ 1Þ! ð4m  i1  i2  k1  k2  iÞ!ði1 þ i2 þ k1 þ k2 þ iÞ! ð4n  j1  j2  l1  l2  jÞ!ðj1 þ j2 þ l1 þ l2 þ jÞ! ðm  1  i1 Þ!i1 !ðm  k1 Þ!k1 !ðn  j1 Þ!j1 !ðn  1  l1 Þ!l1 ! ðm  1  i2 Þ!i2 !ðm  k2 Þ!k2 !ðn  j2 Þ!j2 !ðn  1  l2 Þ!l2 ! ðQ i1 ;j1  Rk1 ;l1 ÞðQ i2 ;j2  Rk2 ;l2 Þ    2 m!m!n!n!m!m!n!n! Xm1 Xm1 Xn1 Xn1 Xm Xm Xn Xn 2 ¼ i1 ¼0 i2 ¼0 j1 ¼0 j2 ¼0 k1 ¼0 k2 ¼0 l1 ¼0 l2 ¼0 j ð4m þ 1Þ!ð4n þ 1Þ! i ð4m  i1  k1  i2  k2  iÞ!ði1 þ k1 þ i2 þ k2 þ iÞ! ð4n  l2  j2  l1  j1  jÞ!ðl2 þ j2 þ l1 þ j1 þ jÞ! ðm  1  i1 Þ!i1 !ðm  1  i2 Þ!i2 !ðn  l2 Þ!l2 !ðn  l1 Þ!l1 ! ðn  1  j2 Þ!j2 !ðn  1  j1 Þ!j1 !ðm  k1 Þ!k1 !ðm  k2 Þ!k2 ! h i ðQ i1 ;l1  Q i2 ;l2 ÞðRk1 ;j1  Rk2 ;j2 Þ  ðQ i1 ;l1  Rk1 ;j1 ÞðQ i2 ;l2  Rk2 ;j2 Þ :

2.4.2. Coefficients computation for B-spline surfaces We can extend the above coefficient computation to the B-spline surfaces as follows. First subdivide the B-spline surface into n Bézier surfaces by knot insertion (see Fig. 3). Denote the located parameter domain of the r Bézier surface as ður0 ; v r0 Þ ! ður1 ; v r1 Þ. Then the coefficients ci;j can be computed as follows

where

Z

cr;i;j ¼

ur 1

Z vr 1 v r0

ur 0

kXu  Xv k2 B2i ðuÞB2j ðv Þdudv :

ð24Þ

The coefficient in Eq. (24) can be computed on the corresponding Bézier surface by the following parameter transformations.



u ¼ ð1  1Þur0 þ 1ur1

v ¼ ð1  sÞv r

0

þ sv r 1

ð25Þ

:

Substituting Eqs. (25) into Eq. (24), we have cr;i;j ¼

Z

ur 1 ur 0

R1 R1 ¼

0

0

Z vr 1 v r0

kX1  Xs k2 ð

@1 2 @s 2 2 Þ ð Þ B ðuÞB2j ðv Þdudv @u @ v i

kX1  Xs k2 B2i ðð1  1Þur0 þ 1ur1 ÞB2j ðð1  sÞv r0 þ sv r1 Þd1ds : ður1  ur0 Þðv r1  v r0 Þ

The explicit representation of the above coefficient can be further derived similarly as the computation for Bézier surfaces (see Appendix A). After we obtain the above coefficients, the left optimization procedure is the same as the Bézier surface. 2.4.3. Numerical quadrature method for NURBS surfaces For rational surfaces, the exact integration via partialfraction decomposition is too expensive. Here we use 2D compound trapezoidal [3] to evaluate the coefficients in Eq. (10) as follows. Sample 2l points uniformly in the u interval and 2k points uniformly in the v interval

 Fig. 3. B-spline surface subdivision: (a) the original B-spline surface (venus); and (b) the resultant Bézier surfaces.

0 ¼ u0 < u1 < u2 <    < u2l ¼ 1 0 ¼ v 0 < v 1 < v 2 <    < v 2k ¼ 1

:

Y.-J. Yang et al. / Graphical Models 76 (2014) 43–55

(a)

(b)

(c)

(d)

(e)

(f)

49

Fig. 4. Area preserving optimization for a Bézier surface: (a) texture mapping of the original surface parameterization; (b) texture mapping of the resultant surface parameterization with a ¼ 0:47 and b ¼ 0:29; (c) iso-parameter curves of the original parameterization; (d) iso-parameter curves of the resultant parameterization; (e) the color-coding image of the original surface; and (f) the color-coding image of the resultant surface.

(a)

(b)

(c)

(d)

Fig. 5. The intersection point of two planar algebraic curves derived from (a) the Bézier surface in Fig. 4, (b) the face surface in Fig. 6, (c) the joint bone surface in Fig. 7, and (d) the vase surface in Fig. 8.

50

Y.-J. Yang et al. / Graphical Models 76 (2014) 43–55

These subdivisions determine the step sizes h ¼ 1=ð2lÞ and p ¼ 1=ð2kÞ. Let

Hðu; v Þ ¼ kXu  Xv k2 B2i ðuÞB2j ðv Þ: For the coefficient ci;j , the numerical approximation has the form

ci;j ¼ T þ Eerror ; where

8 P P2k T ¼ h14g1 2l > i1 ¼0 j1 ¼0 Ai1 Aj1 Hðui1 ; v j1 Þ > > > 4 Hðu 4 ^ ;^ >  ; @ v Þ vÞ 4 4 @ Hðu > > > Eerror ¼  h @u4 þp @v 4 > > 180 < h1 ¼ 2l1 ; > > 1 > g 1 ¼ 2k > > > > > ui1 ¼ i1 h1 > > : v j1 ¼ j1 g 1 and

Ai1 ¼



1 i1 ¼ 0; 2l ; Aj1 ¼ 2 i1 ¼ 1; 2; . . . ; 2l  1



1 j1 ¼ 0; 2k ; 2 j1 ¼ 1; 2; . . . ; 2k  1

^; v ; v ^ Þ and ðu  Þ in the parameter domain. We can for some ðu increase the sampling density l and k to decrease the error term Eerror . Without explicit explanation, we use l; k ¼ 10 in the remainder of this paper such that the error term has the following form

Eerror ¼ 

 ;v Þ @ 4 Hðu @u4

4

dðu; v Þ ¼ jkXu  Xv k  Aj: The parameterization Xðu; v Þ is area-preserving at a point ðu; v Þ if and only if dðu; v Þ ¼ 0. In Fig. 4, we can see that the original Bézier surface is far from being equiareal. By using the derived expression in Section 2.4.1, we explicitly compute the coefficients in Eq. (10). Substituting the obtained coefficients into Eqs. (15) and (16), we obtain the intersection point in ð0; 0Þ ! ð1; 1Þ of the two planar curves (see Fig. 5(a)) ð0:47; 0:29Þ. From the texture mapping, iso-parameter curves and color-coding area stretch in Fig. 4, we can see

(a)

(b)

(c)

(d)

(e)

(f)

u ;v Þ þ @ Hð @v 4 : 7 2:88  10 ^^

For a given NURBS surface, the numerator of Eerror is bounded and Eerror is seven order of magnitude smaller than that bound with our sampling density l ¼ k ¼ 10. The algorithm flow of the double numerical integration for the coefficient computation is as follows. 1. 2. 3. 4. 5.

The non-uniform local area stretching induced by the parameterization can be characterized by the surface/ parameter area ratio associated with each point, denoted as ki . If the parameterization is conformal, ki is also called the conformal factor. In the color-coding images, we measure for a parameterization the area stretch kXu  Xv k, and use colors to illustrate the absolute value of the difference between this area stretch and the surface area. In other words, the color-coding images illustrate the following function across the surface

Set the initial value l0 ¼ k0 ¼ 10. Compute T 0 ¼ T in Equation (26). Double the sampling intensity lr ¼ 2lr1 and kr ¼ 2kr1 . Compute the new approximation T r in Equation (26). If jT r  T r1 j < , exit; otherwise perform r ¼ r þ 1 and goto Step 3.

After we obtain the coefficient ci;j for the NURBS surface, the optimization algorithm is the same as the Bézier and B-spline surfaces. 3. Experimental results Our algorithm is implemented on a PC with Intel Duo CPU 3.06 GHZ, 2G Memory and Microsoft Visual Studio 2008. To show the performance of our algorithm, several examples are given as follows. Fig. 4 illustrates the optimization result of a Bézier surface of degree 3  3. To show the performance, each surface is illustrated using three methods: texture mapping, iso-parameter curves and color-coding area distortion. The iso-parameter curve networks can also be seen as quad tessellations of the given freeform surfaces.

Fig. 6. Area preserving optimization for a B-spline face surface: (a) texture mapping of the original surface parameterization; (b) texture mapping of the resultant surface parameterization with a ¼ 0:59 and b ¼ 0:36; (c) iso-parameter curves of the original parameterization; (d) iso-parameter curves of the resultant parameterization; (e) the colorcoding image of the original surface; and (f) the color-coding image of the resultant surface.

Y.-J. Yang et al. / Graphical Models 76 (2014) 43–55

that equiareality of the resultant surface has been improved. Another example is given in Fig. 6, where a B-spline face surface with 17 by 17 control points is optimized using our method. The coefficients are computed using the method presented in Section 2.4.2. The two planar algebraic curves as well as the intersection point are illustrated in Fig. 5(b). From Fig. 6, we can see that our method can improve the equiareality of the given face B-spline surface.

51

In the third example, a NURBS joint bone surface is optimized using our method (see Fig. 7). The coefficients are computed numerically using the method in Section 2.4.3. Two planar algebraic curves as well as their intersection are given in Fig. 5(c). Another NURBS vase surface is also given in Fig. 8. The planar curve intersection are given in Fig. 5(d). Also, to show the performance of our algorithm to optimize the surface generated by the NURBS modeler, an example of a NURBS car top surface is given in Fig. 9.

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 7. Area preserving optimization for a NURBS bone surface: (a) texture mapping of the original surface parameterization; (b) texture mapping of the resultant surface parameterization with a ¼ 0:32 and b ¼ 0:37; (c) iso-parameter curves of the original parameterization; (d) iso-parameter curves of the resultant parameterization; (e) the color-coding image of the original surface; and (f) the color-coding image of the resultant surface.

52

Y.-J. Yang et al. / Graphical Models 76 (2014) 43–55

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 8. Area preserving optimization for a NURBS vase surface: (a) texture mapping of the original surface parameterization; (b) texture mapping of the resultant surface parameterization with a ¼ 0:24 and b ¼ 0:77; (c) iso-parameter curves of the original parameterization; (d) iso-parameter curves of the resultant parameterization; (e) the color-coding image of the original surface; and (f) the color-coding image of the resultant surface.

4. Conclusions and future works In order to improve equiareality of given freeform surfaces, an optimization algorithm is presented in this paper

based on Möbius transformations. By setting the partial derivatives of formulated equiareality energy to be zero, the optimization problem is reduced to solving the intersection of two planar algebraic curves, whose coefficients

53

Y.-J. Yang et al. / Graphical Models 76 (2014) 43–55

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 9. Area preserving optimization for a NURBS car top surface: (a) texture mapping of the original surface parameterization; (b) texture mapping of the resultant surface parameterization with a ¼ 0:70 and b ¼ 0:60; (c) iso-parameter curves of the original parameterization; (d) iso-parameter curves of the resultant parameterization; (e) the color-coding image of the original surface; and (f) the color-coding image of the resultant surface.

are computed explicitly for Bézier and B-spline surfaces, while numerically for NURBS surfaces. Four examples are given to show the performance of our algorithm. Our method can improve the equiareality of given freeform surfaces. Sometimes, the parameterization can only be improved to a limited degree. The reason is as follows. The Möbius transformations only change the distribution of iso-parameter curves, not the shape of them. With Möbius transformations, the degrees of freedom that we gain may still be not enough. In order to achieve better results, we may consider higher order composite reparameterization in future. From Fig. 5, we can see that there always exists a unique intersection point for the examples we tested. The convexity demonstration of the energy in Eq. (11) is left as a future work. Besides the surface rendering, another potential area that may benefit from this work is the sampling of freeform surfaces for the inspecting planning on CMM [27,28]. To sample the freeform surfaces, a major issue is to control the distribution of the sampling points. Based on the deviation measurement, optimization algorithms using iterative methods are presented in [27,28], which involve intensive computations. By generating more area-preserving parameterizations using Möbius transformations, the optimal sampling of the surface can be reduced to uniformly sample the parameter domain of the optimal parameterization. Unlike the methods in [27,28], the computation using Möbius transformations is deterministic (the energy is convex).

and Technology Support Program (2013BAH39F00), Shandong Province Outstanding Young Scientist Research Award Foundation (BS2012DX014), Independent Innovation Foundation of Shandong University, IIFSDU (2012TB012) and Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry (SRF for ROCS, SEM). Appendix A. The coefficient computation for B-spline surfaces Let

8 f ¼ 1  ur0 > > > > > g ¼ 1  v r0 > > > < i ¼ 1  ur1 : > > j ¼ 1  v r1 > > > > > l ¼ ður0 i þ ur1 fÞ > : m ¼ ðv r0 j þ v r1 gÞ Using the above notations, the Bernstein basis of can be written as follows.

1 and s

8 2 B0 ðð1  1Þur0 þ 1ur1 Þ ¼ B20 ð1Þf2 þ B21 ð1Þfi þ B22 ð1Þi2 > > > > > > > B21 ðð1  1Þur0 þ 1ur1 Þ ¼ 2B20 ð1Þur0 f þ 2B21 ð1Þl þ 2B22 ð1Þur1 i > > > > > B2 ðð1  1Þu þ 1u Þ ¼ B2 ð1Þu2 þ B2 ð1Þu u þ B2 ð1Þu2 < 2

r0

r1

0

r0

1

r0

r1

2

r1

Acknowledgments

> > B20 ðð1  sÞv r0 þ sv r1 Þ ¼ B20 ðsÞg2 þ B21 ðsÞgj þ B22 ðsÞj2 > > > > > 2 > B1 ðð1  sÞv r0 þ sv r1 Þ ¼ 2B20 ðsÞv r0 g þ 2B21 ðsÞm þ 2B22 ðsÞv r1 j > > > > : 2 B2 ðð1  sÞv r0 þ sv r1 Þ ¼ B20 ðsÞv 2r0 þ B21 ðsÞv r0 v r1 þ B22 ðsÞv 2r1

This work was supported by the China National Natural Science Foundation (61202146, 61272243, 61070093, 61332015 and U1035004), the China National Science

From the above equation, B2i ðð1  1Þur0 þ 1ur1 ÞB2j ðð1  sÞ v r0 þ sv r1 Þ; i ¼ 0; 1; 2; j ¼ 0; 1; 2 can be expressed as a matrix form

:

54

0

Y.-J. Yang et al. / Graphical Models 76 (2014) 43–55

B20 ðð1  1Þur0 þ 1ur1 ÞB20 ðð1  sÞv r0 þ sv r1 Þ

1

0

B20 ð1ÞB20 ðsÞ

1

C C B 2 B 2 B B0 ðð1  1Þur0 þ 1ur1 ÞB21 ðð1  sÞv r0 þ sv r1 Þ C B B0 ð1ÞB21 ðsÞ C C C B B C C B 2 B 2 2 2 B B0 ðð1  1Þur0 þ 1ur1 ÞB2 ðð1  sÞv r0 þ sv r1 Þ C B B0 ð1ÞB2 ðsÞ C C C B B B B2 ðð1  1Þu þ 1u ÞB2 ðð1  sÞv þ sv Þ C B B2 ð1ÞB2 ðsÞ C C B 1 B 1 r0 r1 r0 r1 C 0 0 C C B 2 B B B ðð1  1Þur þ 1ur ÞB2 ðð1  sÞv r þ sv r Þ C ¼ AT B B2 ð1ÞB2 ðsÞ C i;j B 1 C B 1 1 1 0 1 0 1 C C C B 2 B B B1 ðð1  1Þur0 þ 1ur1 ÞB22 ðð1  sÞv r0 þ sv r1 Þ C B B21 ð1ÞB22 ðsÞ C C C B B C C B 2 B 2 B B2 ðð1  1Þur0 þ 1ur1 ÞB20 ðð1  sÞv r0 þ sv r1 Þ C B B2 ð1ÞB20 ðsÞ C C C B B C B 2 B B2 ð1ÞB2 ðsÞ C 2 A @ B2 ðð1  1Þur0 þ 1ur1 ÞB1 ðð1  sÞv r0 þ sv r1 Þ A @ 2 1 B22 ð1ÞB22 ðsÞ

B22 ðð1  1Þur0 þ 1ur1 ÞB22 ðð1  sÞv r0 þ sv r1 Þ

where column vectors Ai;j ¼ ½A0;0 A0;1 A0;2 A1;0 A1;1 A1;2 A2;0 A2;1 A2;2  are defined as follows

0

f2 g2

2f2 v r0 g

f2 v 2r0

B 2 B f gj 2f2 m f2 v r 0 v r 1 B B 2 2 2 Bf j 2f v r1 j f2 v 2r1 B B B fig2 2fiv r0 g fiv 2r0 B B figj 2fim fiv r0 v r1 B B B fij2 2fiv r1 j fiv 2r1 B B 2 2 2i2 v r0 g i2 v 2r0 Bi g B B i2 gj 2 2 2i m i v r0 v r 1 @

i2 j2 2i2 v r1 j

i2 v 2r1

2ur0 fg2

4ur0 fv r0 g

2ur0 fv 2r0

u2r0 g2

2u2r0 v r0 g

2ur0 fgj

4ur0 fm

2ur0 fv r0 v r1

u2r0 gj

2u2r0 m

2ur0 fj2

4ur0 fv r1 j

2ur0 fv 2r1

u2r0 j2

2l v

R1 R1 0

0

2lg

4l v r 0 g 4lm

2lv r0 v r1

ur0 ur1 gj

2ur0 ur1 m

2lj2

4l v r 1 j

2lv 2r1

ur0 ur1 j2

2ur0 ur1 v r1 j

2ur1 ig2

4ur1 iv r0 g

2ur1 iv 2r0

u2r1 g2

2u2r1 v r0 g

2ur1 igj

4ur1 im

2ur1 iv r0 v r1

2ur1 ij2

4ur1 iv r1 j

2ur1 iv 2r1

2 r0

kX1  Xs k2 B2i ðð1  1Þur0 þ 1ur1 ÞB2j ðð1  sÞv r0 þ sv r1 Þd1ds ður1  ur0 Þðv r1  v r0 Þ

1 B20 ð1ÞB20 ðsÞ C B 2 B B0 ð1ÞB21 ðsÞ C C B C B 2 B B0 ð1ÞB22 ðsÞ C C B B B2 ð1ÞB2 ðsÞ C C B 1 0 C B R1 R1 2 T B 2 kX1  Xs k Ai;j B B1 ð1ÞB21 ðsÞ C 0 0 Cd1ds C B 2 B B1 ð1ÞB22 ðsÞ C C B C B 2 B B2 ð1ÞB20 ðsÞ C C B B B2 ð1ÞB2 ðsÞ C A @ 2 1 0

¼

2u2r0 v r1 j

2lgj

2

Finally we obtain the coefficient for the B-spline surface as follows. cr;i;j ¼

[3] G. Dahlquist, A. Björck, Numerical Methods in Scientific Computing, SIAM press, 2008. [4] M.P. do Carmo, Differential Geometry of Curves and Surfaces, Prentice Hall, 1976. [5] G. Farin, NURBS from Projective Geometry to Practical Use, second ed., A K Peters, 1999. [6] G. Farin, Rational quadratic circles are parametrized by chord length, Computer Aided Geometric Design 23 (9) (2006) 722–724. [7] R.T. Farouki, Optimal parameterizations, Computer Aided Geometric Design 14 (2) (1997) 153–168. [8] M.S. Floater, K. Hormann, Surface parameterization: a tutorial and survey, in: N.A. Dodgson, M.S. Floater, M.A. Sabin (Eds.), Advances in Multiresolution for Geometric Modelling, Mathematics and Visualization, Springer, Berlin, Heidelberg, 2005, pp. 157–186. [9] L. He, S. Schaefer, K. Hormann, Parameterizing subdivision surfaces, ACM Transactions on Graphics (TOG) 29 (4) (2010). [10] L. He, C. Loop, S. Schaefer, Improving the parameterization of approximate subdivision surfaces 31(7) (2012) 2127-2134.

B22 ð1ÞB22 ðsÞ ður1  ur0 Þðv r1  v r0 Þ

1 cr;i;j;0;0 C Bc B r;i;j;0;1 C C B B cr;i;j;0;2 C C B C Bc B r;i;j;1;0 C C B T B Ai;j B cr;i;j;1;1 C C C B B cr;i;j;1;2 C C B B cr;i;j;2;0 C C B C B @ cr;i;j;2;1 A cr;i;j;2;2 ; ¼ ður1  ur0 Þðv r1  v r0 Þ 0

where cr;i;j;m;n ; m ¼ 0; 1; 2; n ¼ 0; 1; 2 are the coefficients in Eq. (10) for the rth sub Bézier surface, which can be computed explicitly by the method in Section 2.4.1. References [1] G. Cheung, R. Lau, F. Li, Incremental rendering of deforming NURBS surfaces, in: VRST 03: Proceedings of the ACM Symposium on Virtual Reality Software and Technology, ACM Press, New York, NY, USA, 2003, pp. 48–55. [2] P. Costantini, R.T. Farouki, C. Manni, A. Sestini, Computation of optimal composite re-parameterizations, Computer Aided Geometric Design 18 (9) (2001) 875–897.

2

ur0 ur1 g

u2r1

gj 2 2 ur1 j

2ur0 ur1 v r0 g

2u2r1 m 2u2r1 v r1 j

u2r0 v 2r0

1

C u2r0 v r0 v r1 C C C C u2r0 v 2r1 C C 2 ur0 ur1 v r0 C C ur0 ur1 v r0 v r1 C C C ur0 ur1 v 2r1 C C C u2r1 v 2r0 C C 2 ur1 v r0 v r1 C A u2r1 v 2r1

[11] B. Jüttler, A vegetarian approach to optimal parameterizations, Computer Aided Geometric Design 14 (9) (1997) 887–890. [12] E.T.Y. Lee, M.L. Lucian, Möbius reparametrizations of rational Bsplines, Computer Aided Geometric Design 8 (3) (1991) 213–215. [13] W. Lu, Curves with chord length parameterization, Computer Aided Geometric Design 26 (3) (2009) 342–350. [14] W.M.M. Ng, S.T. Tan, Incremental tessellation of trimmed parametric surfaces, Computer-Aided Design 32 (4) (2000) 279–294. [15] B.H. Ong, An extraction of almost arc-length parameterization from parametric curves, Annals of Numerical Mathematics 3 (1996) 305– 316. [16] N. Ray, W.C. Li, B. Levy, A. Sheffer, P. Alliez, Periodic global parameterization, ACM Transactions on Graphics 25 (4) (2006) 1460–1485. [17] T.W. Sederberg, Algorithm for algebraic curve intersection, Computer-Aided Design 21 (9) (1989) 547–554. [18] M. Shpitalni, Y. Koren, C.C. Lo, Realtime curve interpolators, Computer-Aided Design 26 (11) (1994) 832–838. [19] F.C. Wang, D.C.H. Yang, Nearly arc-length parameterized quintic spline interpolation for precision machining, Computer-Aided Design 25 (5) (1993) 281–288. [20] F.C. Wang, P.K. Wright, B.A. Barsky, D.C.H. Yang, Approximately arclength parameterized C 3 quintic interpolatory splines, ASME Journal of Mechanical Design 121 (1999) 430–439. [21] U. Wever, Optimal parameterization for cubic spline, ComputerAided Design 23 (9) (1991) 641–644. [22] D.C.H. Yang, T. Kong, Parametric interpolator versus linear interpolator for precision CNC machining, Computer-Aided Design 26 (3) (1994) 225–234. [23] D.C.H. Yang, F.C. Wang, A quintic spline interpolator for motion command generation of computer-controlled machines, ASME Journal of Mechanical Design 116 (1994) 226–231. [24] Y.J. Yang, J.H. Yong, H. Zhang, J.C. Paul, J.G. Sun, Optimal parameterizations of Bézier surfaces, in: 2nd International Symposium on Visual Computing(ISVC06), 2006, pp. 672–681. [25] Y.J. Yang, W. Zeng, H. Zhang, J.C. Paul, J.H. Yong, Projection of curves on B-spline surfaces using quadratic reparameterization, Journal of Graphical Models 72 (5) (2010) 47–59.

Y.-J. Yang et al. / Graphical Models 76 (2014) 43–55 [26] Y.J. Yang, W. Zeng, C.L. Yang, X.X. Meng, J.H. Yong, B.L. Deng, G1 continuous approximate curves on NURBS surfaces, Computer-Aided Design 44 (9) (2012) 824–834. [27] M.R. Yu, Y.J. Zhang, Y.L. Li, D. Zhang, Adaptive sampling method for inspection planning on CMM for free-form surfaces, The International Journal of Advanced Manufacturing Technology 67 (2013) 1967–1975. [28] S.M. Obeidat, S. Raman, An intelligent sampling method for inspecting free-form surfaces, The International Journal of Advanced Manufacturing Technology 40 (11-12) (2009) 1125–1136.

55

[29] S.-S. Yeh, P.-L. Hsu, The speed-controlled interpolator for machining parametric curves, Computer-Aided Design 31 (5) (1999) 349–357. [30] G. Zou, J. Hua, O. Muzik, Non-rigid surace registration using spherical thin-plate splines, in: proceeding of the 10th International Conference on Medical Image Computing and Computer Assisted Intervention, 2007, pp. 367–374. [31] G. Zou, J.X. Hu, X.F. Gu, J. Hua, Authalic parameterization of general surfaces using Lie advection, IEEE Transactions on Visualization and Computer Graphics 17 (12) (2011) 2005–2014.