1 Introduction

B-spline free-form
deformation of polygonal
object as trimmed Bézier
surfaces
Jieqing Feng1 , Tomoyuki Nishita2 ,Xiaogang Jin1 , Qunsheng Peng1
1 State Key Laboratory of CAD&CG, Zhejiang

University, Hangzhou, 310027, P.R. China
E-mail: {jqfeng, jin,peng}@cad.zju.edu.cn
2 Department of Complexity Science andEngineering, University of Tokyo, Hongo 7-3-1,
Bunkyo, Tokyo 113-0033, Japan
E-mail: nis@is.s.u-tokyo.ac.jp

Published online: November 6, 2002c Springer-Verlag 2002
Free-form deformation is a powerful shapemodification tool. How to approximate orcompute the real deformation of a polygonal object is still problematic. In this paper,a new solution is proposed for this problem. First, a special initial B-spline volume
is defined whose Jacobian is an identity matrix. The accurate deformation is as trimmed
tensor-product Bézier surfaces. The description of the trimmed surfaces is consistent
with that in the industrial standard, STEP.
The degrees of the Bézier surfaces are lower
than the theoretical results. Compared withprevious algorithms, the proposed algorithm
has the advantages of both storage and runtime.
Key words: Free-form deformation –B-spline – Polygons – Bézier surface –Trimmed surface

The Visual Computer (2002) 18:493–510
Digital Object Identifier (DOI) 10.1007/s00371-002-0171-1

Free-form deformation, abbreviated as FFD, wasfirst proposed by Sederberg and Parry [23]. It is aneffective tool for geometric shape modification, and
it is independent of the representation and topologyof object. The FFD can be concisely described as follows. Firstly, the object is frozen into an intermediate
flexible space, such that when the shape of the intermediate space is changed, the object will follow itschange. The 3D tensor-product Bézier volume wasadopted as an appropriate flexible space by Sederberg and Parry. Its initial shape is an orthogonal parallelepiped. Later the FFD was extended to include
the B-spline FFD [13], extended FFD [4], rational
FFD [16], direct manipulation FFD [15], NURBS
FFD [17], continuous FFD [1], arbitrary topologicallattice FFD [19], and Dirichlet FFD [20]. It has been
applied to geometric modelling [6, 14] and computer
animation [3, 5, 16]. Bechmann produced a survey
about the various space deformation techniques [2].
The problem of sampling for the polygonal object
deformation is unsolved as yet. Theoretically, the deformation should act upon the whole polygonal object. However, this can not be achieved, since thepolygonal object can be only represented discretely.
Thus, the deformation can only act on the sampling
points of the polygonal object. If the number of sampling points is insufficient, aliasing phenomena will
occur in the deformation result, and this will degrade
the final deformation effect. Heretofore, there were
three types of solutions to the problem. They are
the uniform subdivision approach, the error-test approach and the accurate FFD approach.

1.1 Uniform subdivision approach
The polygonal object is uniformly subdivided on the
midpoint of each edge recursively. This subdivisionoccurs regardless of the geometric shape of the object, and the stop criterion is decided by the user. Allof the meshes in the object will be treated equally,
and therefore it is a brute-force solution; but it is thesolution that is adopted for most popular 3D computer animation software. The drawback of the solution is obvious: many tiny meshes must be generated
if a perfect deformation result is required, and thisincreases the rendering and storage burden.

1.2 Error-test approach
Parry first proposed an adaptive subdivision algorithm for solving the sampling problem [22]. The

494

J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

two triangles that share a common “long side” arehalved. The long side for a triangle is defined as follows. Initially, an arbitrary edge of the triangle islabelled, and then, after subdivision, the other two
undivided edges are labelled as the long sides forthe two newly generated triangles. Parry adopted two
criteria to decide the subdivision of a triangle: thesize of the triangle as projected on a screen and thetriangle’s “curvature”. Griessmair and Purgathofer
proposed another solution for the polygonal object
deformation [13]. In their method, the middle pointof an edge after deformation is checked to determine whether it is far away from its true position.
If the disparity is greater than a predefined threshold, the edge will be halved. Inspired by the existing
work on domain discretisation in finite element analysis, Nimscheck proposed an advancing-front meshing algorithm for FFD [21]. Each edge of the initial
polygon is recursively subdivided based on either
its post-deformation length or the disparity between
its endpoint tangents. Then a new vertex is inserted
for triangulation, which approximates the ideal deformation result. Besides subdivision, Gain merged
the triangles that are approximately coplanar [12].
This method can decrease the number of meshes inthe deformed object. All of the four methods described above follow a similar strategy: sampling,
testing and subdivision. According to the sampling
theory, some singular cases cannot be avoided [12].
These are inappropriate for processing by the above
methods. For example, for a wave-shaped curve, if
the normals to both its ends are the same, the above
methods may detect it falsely.

lated into triangular meshes or converted into threetensor-product Bézier surfaces in advance, which isan additional cost. Another weakness is the degree ofthe triangular Bézier surface. It is equal to the degreeof the B-spline volume. When a triangular Bézier
surface is converted into three tensor-product Bézier
surfaces, the degree of the resulting surfaces is twice
the degree of the triangular Bézier surface. This is
the maximal degree according to the results in [7, 8].
For example, when the degree of the B-spline volume in the accurate FFD is 2 × 2 × 2, the degree ofthe triangular Bézier surface is 6. The degree of thecorresponding tensor-product Bézier surface will be
6 × 6. This high degree means high computational,
storage and display costs.
In this paper, a new accurate FFD method is proposed. First, a new definition of the initial B-spline
volume is given. It is an end-interpolation B-spline
volume, whose Jacobian matrix is an identity matrix. Then, a robust subdivision algorithm for the
polygonal object is proposed, which is free of degenerate polygons. The deformation result is trimmed
tensor-product Bézier surfaces whose description
is consistent with that in the industrial standard,
STEP [24]. Furthermore the degree of the tensorproduct Bézier surface is always lower than thetheoretical results [8]. Experimental results and theoretical analysis show that the new accurate FFD has
advantages both in computational cost and storage,
and can make use of graphics hardware to render thetrimmed surfaces.

2 Overview of the proposed algorithm

1.3 Accurate FFD approach
The accurate FFD methods proposed by the present
authors can avoid the problems of the error-test
approach [10, 11]. Traditionally, the FFD acts ona point, whereas accurate FFD acts on each triangle. The deformation result is a triangular Bézier
surface according to the Bernstein polynomial composition [8]. There are two weaknesses in accurate
FFD methods. Firstly, a triangular Bézier surface is
not directly supported by current modelling systems
and rendering pipelines. In general, these support
tensor-product Bézier surfaces. Sometimes the rendering procedure can be accelerated by the hardware.
To make full use of the graphics hardware ability andmake it compatible with the current modelling system, the triangular Bézier surface needs to be tessel-

This research is inspired by the Bernstein polynomial composition [8]. The composition between
a planar quadrilateral and a Bézier volume is a tensorproduct Bézier surface. With this result, the sampling
problem in the first FFD method can be solved [23].
Since the Bézier volume can only produce global
deformation, the B-spline volume is preferable inpractical use. As we know, the B-spline volume
can be converted into piecewise continuous Bézier
volumes [9]. Thus, the Bernstein polynomial composition can be employed here for the B-spline FFDof the polygonal object. With proper subdivision of
the polygonal object in the parametric space of the
B-spline volume, the B-spline FFD can be converted
as Bézier FFDs. According to the results in [8], thedegree of the resultant tensor-product Bézier surface

J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

Inputing
Polygonal
ObjectInitial
B-Spline
VolumePolygonal
ObjectMarking
Co-Planar

Subdivision

Polygons

Determination ofTrimmed Surface
DegreeSampling
andInterpolation

495

OutputingDeformed
ObjectInteractively
Editing
B-Spline
VolumeFig. 1. The flow chart of the proposed algorithm

is equal to the tensor product of the degree of theBézier volume. According to our analysis, this degree will never be reached in the work in this paper;that is, the degrees of the resultant tensor-productBézier surfaces are always lower than this. Unlike
the previous algorithms in [8, 18], we proposed a fast
polynomial interpolation algorithm to compute thecontrol points of the tensor-product Bézier surface.
The flow chart of the proposed algorithm is shownin Fig. 1.
The initial shape of the B-spline volume can be either arbitrary or regular (orthogonal parallelepiped).
An initial B-spline volume with an arbitrary shape isconvenient for interactive shape modification. However, it is discommodious to define or automaticallygenerate such an initial B-spline volume with an arbitrary shape. In practical use, an orthogonal parallelepiped is preferable, whose extension is identicalwith the object bounding box. Many computer animation systems adopt such an approach. In this paper
we consider the regular case.
First, a new initial B-spline volume is introduced inthis paper. The proposed initial B-spline volume has
two desirable properties. To start with, the point coordinates in the object space are identical with its
local coordinates corresponding to the B-spline volume, and the Jacobian of the B-spline is an identity
matrix. After the initial B-spline volume is defined,
the polygonal object is subdivided in the parametric space according to the knot vectors distribution.
Each sub-polygon in the object after subdivision lies
within a “knot box” (see definition in Sect. 3.2) of the
B-spline volume. In fact, the B-spline volume limited in a knot box is a Bézier volume. The next stepis to classify the coplanar sub-polygons in each knot
box. At the same time, the degree of the resultant
tensor-product Bézier surface is determined. For the
coplanar sub-polygons in a knot box, it is not necessary to compute each sub-polygon deformation. We

will only compute the deformation of the planar rectangular bounding box of the coplanar sub-polygons.
The deformation of the coplanar sub-polygons ina knot box is the trimmed tensor-product Bézier surfaces. The trimmed curves are the edges of the coplanar sub-polygons. The control points of the tensorproduct Bézier surface are evaluated through polynomial interpolation, where the interpolation matrix isconstant.

3 B-spline FFD polygonal object astrimmed Bézier surfaces
The proposed algorithm is described in detail in thissection. The section is organized according to theflow chart in Fig. 1.

3.1 A new definition of initial B-spline
volume3.1.1 The initial B-spline volume and local
coordinates
From a mathematical point of view, FFD can be
described as two mappings: 3 → 3 . The first mapping is from the object space to the parametric space,and it freezes the object in the B-spline volume. Let
the mapping be noted as L : (x, y, z) → (u, v, w).
The mapping L is computed once after the initial B-spline volume is defined. During the deformation, (u, v, w) is fixed for the point (x, y, z).
The second mapping is the B-spline volume itself
R : (u, v, w) → (x̃, ỹ, z̃). The relationship between(x, y, z) and (x̃, ỹ, z̃) is linked by the triple (u, v, w),
namely, the local coordinates of point (x, y, z).
For the given initial B-spline volume, there are two
ways to compute the local coordinates (u, v, w). One
way is called forward mapping, which is a linear

496

J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

a

b

c

d

Fig. 2a–d. The degree of the B-spline volume is 2 × 2 × 1; its lattice is 4 × 4 × 2. a The original object. b The initial B-spline volume
in end interpolation. Its control points are evenly distributed. There is warp in the deformation. c The uniform B-spline volume with
its identity Jacobian. d The end-interpolation B-spline volume converted from (c). The object in both (c) and (d) has no warp when it
is deformed by the initial B-spline volume.

transformation from (x, y, z) to (u, v, w). It is a fast
way since only the linear transformation is computed. However, the object may be warped after it
is deformed by the initial B-spline volume. An example of this is given in Fig. 2b, which shows thedeformation of the object in Fig. 2a by the initial
B-spline volume. The reason for this is that the initial B-spline volume in general is non-linear. The
previous accurate FFD algorithms adopted this mapping technique [10, 11]. The second way is calledbackward mapping, where for each point (x, y, z) on
the object, a triple (u, v, w) is computed such that
R(u, v, w) = (x, y, z). In general, there is no analytical solution for this problem and numerical method
is adopted. By using the backward mapping, there
is no warp when the object is deformed by the initial B-spline volume, but it is a time-consuming procedure, and the local coordinates are non-linearly
distributed in general. It is not suitable for the accurate FFD algorithms. Could we design an initial
B-spline volume that has the advantages of both ofthese two kinds of mapping? The answer is yes.
We will give a constructive method in the following
section.
3.1.2 Definition of uniform B-spline volume
Let the initial B-spline volume R(u, v, w) be notedasR(u, v, w)=nu −1 nv −1 nw −1
i=0

j=0 k=0

Rijk Ni,ku (u)N j,kv (v)Nk,kw (w),(1)

where n u , n v and n w are the numbers of the control points, ku , kv and kw are the degrees, Ni,ku (u),N j,kv (v) and Nk,kw (w) are the B-spline functions
along the u, v and w directions respectively. Rijk are
the control points, namely, the lattice in general. The
object bounding box aligned to the coordinates axis
is noted as (x min , ymin , z min ) and (x max , ymax, z max ).
For simplicity, the x-, y- and z-coordinates are
mapped as u, v and w respectively. Here the uniform
B-spline volume means that both the knot vectors
and the control points are uniformly distributed inthree parametric directions. The knot vectors {u i },{v j } and {wk } are defined as follows:u i = x min + (i − ku ) ∗ (x max − x min )/(n u − ku )(i = 0, 1, . . . , ku + n u ) ;
v j = ymin + ( j − kv ) ∗ (ymax − ymin)/(n v − kv )( j = 0, 1, . . . , kv + n v ) ;






wk = z min + (k − kw ) ∗ (z max − z min )/(n w − kw )(k = 0, 1, . . . , kw + n w ) .
From the definition, we know that the valid parametric domain is identical with the object boundingbox; that is, [x min , x max ] × [ymin, ymax ] × [z min , z max ].
Next, the control point Rijk = (x ijk, yijk, z ijk ) is defined by the following equation: 
x ijk = x min + i ∗ (x max − x min )/(n u − 1)
(i = 0, 1, . . . , n u − 1) ;



 
yijk = ymin + j ∗ (ymax − ymin)/(n v − 1)
( j = 0, 1, . . . , n v − 1) ;







 z ijk = z min + k ∗ (z max − z min )/(n w − 1)
(k = 0, 1, . . . , n w − 1) .

J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

a

b

c

497

d

Fig. 3a–d. The control point configurations for the end-interpolation B-spline volume of degree 3 × 3 × 1. The numbers of the control
points are noted under each subfigure. a 4 × 4 × 2; b 5 × 5 × 2; c 6 × 6 × 2; d 7 × 7 × 2

In general, the extension of R(u, v, w) with theabove control points {(x ijk, yijk, z ijk )} is smaller
than the object bounding box. If the object is deformed by R(u, v, w), there will be no warp introduced except for constrictions. To make the extension of R(u, v, w) identical with the object bounding box, the control points should be scaled. Let(x min, ymin, z min ) and (x max, ymax, z max ) be noted as
the corner points of the R(u, v, w) extension, whichcorrespond to R(x min , ymin, z min ) and R(x max , ymax ,z max ). Three scalars can be computed by the following equations:sx = (x max − x min )/(x max− x min);s y = (ymax − ymin)/(ymax − ymin);sz = (z max − z min )/(z max − z min ) .
After x ijk, yijkand z ijk are scaled by using scalars
sx , s y and sz respectively, the extension of the new
B-spline volume is identical with the object bounding box. Furthermore, its Jacobian is an identity
matrix. This means that when the object is deformed by such an initial B-spline volume, the object
will remain unchanged. One example is shown in
the Fig. 2c. It is the deformation of the object in
Fig. 2a by a uniform B-spline volume as defined
above.

3.1.3 Definition of end-interpolation B-spline
volumeIn practical applications, an end-interpolation
B-spline volume is preferable, since its convex hullof control points is tighter than that of the uniform
B-spline volume. It can be accomplished by using
the knot-insertion algorithm for the end knots [9].

An example is shown in Fig. 2d, which is a conversion result of Fig. 2c. Besides being free ofwarp, another advantage is that the object point coordinates are identical with its local coordinates.
Thus, linear-transformation computational costs are
saved.
Since the above end-interpolation B-spline volume
results from a uniform one, an interesting result
can be derived according to the knot-insertion algorithm. Let us illustrate the result with the 1D case.
Supposing R(u) is an end-interpolation B-spline
curve converted from a uniform one. It has n control points {Ri }n−1i=0 and its degree is k. Let R0 R1
be noted as vector u. If n = k + 1, it is a Bézier
curve. The control points are evenly distributed.
Else, the results are described by Eq. (2). Thus theend-interpolation B-spline volume can be defined ina straightforward way. Some examples are shownin Fig. 3.

R0 R1 = Rn−2 Rn−1 = u R1 R2 = Rn−3 Rn−2 = 2u
..
..
.
(2)
. =
.
= ..




 Rk−1 Rk = Rn−k−1 Rn−k = ku
Rs Rs+1 = ku (other s).
With the above steps, we defined an initial B-spline
volume. It has the following properties. Its parametric domain is identical with the object bounding box.
It has end interpolation, and it is linear. Its Jacobian is
an identity matrix. When the object is embedded into
its parametric domain by using either forward mapping or backward mapping, the deformation results
are the same. Each object point coordinate is identical to its local coordinate.

498

J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

P0

P0

Q0

P4

P1

P1

P2

a

P0

Q0

P4

P3

Q2

P4

P1

P2

b

P3

P2

Q1

P3
Q1

Q3

c

Fig. 4. a Original polygon and its cycle-linked list. b Knot plane crosses a polygon vertex P0 . The vertex P0 is marked as theintersection and linked with another intersection point Q 1 . c The final linked list after all intersection points are evaluated

3.2 Subdivision of a polygonal object in theparametric space
Before introducing the subdivision algorithm, wefirst define the knot box. A knot box is a 3D
knot interval for the B-spline volume. For example, [u i , u i+1 ] × [v j , v j+1 ] × [wk , wk+1 ] is a knot
box for Eq. (1). It can be indexed by (i, j, k) forku ≤ i < n u , kv ≤ j < n v , kw ≤ k < n w . The B-spline
volume limited on a knot box is a Bézier volume. The purpose of the subdivision is to makeeach polygon in the object lie in a knot box after subdivision. The polygons after subdivision are
called sub-polygons throughout the rest of paper.
The proposed subdivision algorithm is described asfollows:for ( each parametric direction u,v,w ){for ( each polygon in the object ){form the cycle-linked list for the vertex of thepolygon;find indices i0 and i1 such that knot plane ui( i0 < i < i1 ) crosses the polygon;
for ( each knot plane ui ( i0 < i < i1 ) ){compute the intersections between eachknot plane and polygon;
insert the intersections intoa cycle-linked list;}
extract the sub-polygons fromthe cycle-linked list;}}

3.2.1 Computing intersections between the knot
planes and the polygon
We assume that the object is composed of convex polygons. Let {Pi }n−1i=0 be an n-sided polygon.
It may be from the original object or from previous subdivision steps. First a cycle-linked list is built
for the polygon vertex. An example is shown in
Fig. 4a. Then the knot indices i 0 and i 1 are foundsuch that the plane determined by knot u i , vi or wi(i 0 < i < i 1 ) crosses the current polygon. For eachknot plane crossing the polygon, there are two intersections on the polygon. The two intersections
should be linked to each other for traversal purposes in the next step. If the knot plane crossesa polygon vertex, then the crossed polygon vertex
is re-marked as an intersection point, and this islinked with the other intersection point, as shown in
Fig. 4b. When the sub-polygons are extracted fromthe linked cycle list, each polygon vertex will be visited once and each intersection will be visited twice.
Thus each point will be assigned with a Counter:
Counter = 1 for the polygon vertex; Counter = 2
for the intersection.
Here, if a knot plane crosses a polygon edge or one
intersection is coincident with another, then there isno intersection between the knot plane and the polygon. Examples of crossing an edge or special vertex
are shown in Fig. 5a,b. If we do not check such singular cases carefully, degenerate sub-polygons will
be produced. The degenerate polygons can be neglected in our previous methods, since they only
cause some wasteful computations [10, 11]; but, inJ. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

Q0

P0

Q1

P3

P1

P0

P1

499

Q0 P 0

P1

P2

a

Q0
P2

b

Q1

P3

P2

Q1

P3

c

Fig. 5a–d. Singular cases in subdivision, with the knot plan drawn as a bold line: a if the knot plane were not regarded as crossing edge P0 P3 because of a numerical error, then a 0-width polygon P0 Q 0 Q 1 P3 would be generated; b if the knot plane were
not regarded as crossing an extreme vertex P3 , a degenerate triangle Q 0 Q 1 P3 would be produced; c if the knot plane were notregarded as crossing vertex P0 , the polygon P0 Q 0 Q 1 P3 will possess a 0-length edge P0 Q 0

the method proposed here, they cause fatal errors
in the following process steps. This is the reason
why we check carefully for any singular cases in thesubdivision.
3.2.2 Extracting the sub-polygons from the linked
cycle list
For convenience of extracting sub-polygons, the data
structure of the polygon vertex and the intersection is
proposed as follows:class node {double Position[3];// vertex or intersection coordinates
short int Counter;// "1" if vertex ; "2" if intersection
class node *NextNode;// pointer to the next points:// vertex or intersection;
class node *Intersection;// NULL if object vertex}

The algorithm for extracting sub-polygons begins
from a node whose Counter is 1. The key stepis to determine the next node from the linked cycle list. There will be some degenerate sub-polygons
generated in the old subdivision algorithm [10]. The
traversal algorithm to extract a sub-polygon is described in the following. When all sub-polygons
are extracted from the linked list, all Counters
are 0. According to the algorithm, there are three
sub-polygons extracted in Fig. 4c. The three subpolygons in correct vertex order are: P1 P2 Q 1 Q 0 ,Q 0 Q 1 Q 3 Q 2 and Q 3 P3 P4 Q 2 .

ExtractSubPolygon(){pHead = GetFirstNode();pHead->Counter--;if(pHead->NextNode->Counter == 0) {pNode = pHead->Intersection;
status = INTERSECTION_EDGE;}
else {
pNode = pHead->NextNode;status = POLYGON_EDGE;}
while(pNode != pHead) {pNode->Counter--;if(status == POLYGON_EDGE) {if(pNode->Intersection == NULL)
pNode = pNode->NextNode;else {
pNode = pNode->Intersection;
status = INTERSECTION_EDGE;}}
else {
pNode = pNode->NextNode;status = POLYGON_EDGE}}}

Compared with the subdivision algorithm in [10],the proposed algorithm has two improvements. The
new algorithm is more robust than the old one. There
are no degenerate sub-polygons produced, which are“0-area” or “0-edge” sub-polygons. Another improvement is that the new algorithm produces fewer
sub-polygons than before since no-triangulation
step is necessary. Some examples are shown in

500

J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

Figs. 7–13, and a comparison of the sub-polygon
numbers is shown in Table 2 in Sect. 4.

3.3 Marking coplanar sub-polygons anddetermining surface degrees
3.3.1 B-spline deformation of a quadrilateral rather
than a triangle
After subdivision, each sub-polygon lies in a knot
box. In the previous results [10, 11], all sub-polygons
are triangles. The accurate deformation of a triangle
is a triangular Bézier surface, whose degree equals
the degree of the B-spline volume. Unfortunately,
most graphics APIs (such as OpenGL and Direct3D)
and graphics hardware cannot support the display of
a triangular Bézier surface directly. The triangularBézier surface needs to be tessellated into triangularmeshes or converted into three tensor-product Bézier
surfaces in advance. Another inconvenience is thatthe triangular Bézier surface cannot be integrated
into the current systems because most current geometric modelling and computer animation systems
do not support the data structure of triangular Bézier
surfaces directly. To overcome the above problems,
it is necessary to consider tensor-product Bézier surfaces as an accurate deformation result.
According to the results in [8], the B-spline deformation of a quadrilateral in a knot box is a tensorproduct Bézier surface, whose degree is the tensor
product of the degree of the B-spline volume. The
brute-force solution is that the sub-polygons are tessellated into quadrilaterals, and then the quadrilaterals are deformed as tensor-product Bézier surfaces.
This is an inefficient solution. Firstly, there are many
triangles among the sub-polygons in general. Eachtriangle will be tessellated as three quadrilaterals;that is, there will be three tensor-product Bézier surfaces generated. This is a heavy burden in termsof both computation and storage space. Secondly,
some sub-polygons in a knot box may be coplanar,and it is inefficient to deform them one by one. Infact they are parts of a “large” tensor-product Bézier
surface. If such a “large” surface can be computed,
trimmed surfaces can be adopted for representationof the coplanar sub-polygons deformation. Thus,
both computational cost and storage space are saved.
3.3.2 Marking coplanar sub-polygons in a knot box
To mark coplanar sub-polygons, the data structure
Normalis defined as follows. The coplanar sub-polygons in a knot box will be marked and then classified into the Normal structure.
class Normal{
short i,j,k;// The index of the knot box;
double n[3], d;// Plane equation of coplanar polygons;
double s[3], t[3];// Local frames for the bounding box;
short ns, nt;// Degree of resulting Bezier surface.
RECTANGLE box;// 2D bounding box of sub-polygons in s-t plane.
POLYGON *List;// Pointer to the coplanar sub-polygon list}

Each sub-polygon has a plane equation. It can be
described by a normalized normal (n x , n y , n z ) andscalar d as follows:
xn x + yn y + zn z + d = 0.
If a sub-polygon is degenerate, there will be noplane equation. In practical implementations, thiswill cause the algorithm to become unstable, since its
normal will be the zero-vector. Thus, a robust subdivision algorithm is necessary in the previous step.
If two sub-polygons are coplanar, their normalized
normals (n x , n y , n z ) and scalars d should be equal.
Because of numerical error, equal comparisons are
implemented as the following equations:n1 · n2 − 1 < εnd1 − d2  < εd .
In our implementation, εn = 2.2 × 10−8 and εd =1.1 × 10−6 . By using these two thresholds, if two
polygons in the original object are coplanar, thenthe sub-polygons generated from them will be also
coplanar. No visual artefact will occur because of thenumerical error. By using the above two criteria, thecoplanar sub-polygons in a knot box will be classified into the same Normal structure. In general,
there are many Normal structures in a knot box,
whose number is smaller than the number of subpolygons in the knot box.
3.3.3 Determining the degree of the tensor-productBézier surface according to the normal
In this section, the degree of the “large” tensorproduct Bézier surface is determined for each
Normal structure in a knot box. According to the

J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

2D bounding box
w

knot box

t
n

v
s

u
sub-polygons

Fig. 6. The local coordinate system (s, t, n) and the 2D
bounding box of the sub-polygons in the s-t plane

theoretical results in [8], the degree of the tensorproduct Bézier surface, which is the accurate deformation of a quadrilateral by the R(u, v, w) in
Eq. (1), is (ku + kv + kw ) × (ku + kv + kw ). In our algorithm, the degrees of large surfaces are always
lower than this value, and it could be known for
special cases. For example, if a quadrilateral lies in
the knot plane u = ũ, with its normal (n x , n y , n z ) =(1, 0, 0), then the degree of its corresponding tensorproduct Bézier surface is kv × kw or kw × kv . In fact,
the surface is the part of the iso-parameter surface u = ũ.
How can we reduce the surface degrees for generalcases? Let us illustrate this with Fig. 6. Supposing
the normal n = (n x , n y , n z ) has no zero component,
or n x · n y · n z  = 0, and the degree of the B-spline

501

volume is satisfied: ku ≤ kv ≤ kw . Then a normalized vector s orthogonal to the normal n can be defined ass = (sx , s y , 0) =(−n y , n x , 0).
(−n y , n x , 0)

Together with another normalized vector t = n × s,
the three vectors (s, t, n) define a local coordinate
system. The sub-polygons in the Normal structurelie in the s-t plane of the local coordinate system. In
the s-t plane, we can compute a 2D bounding box
for sub-polygons in the Normalstructure, which isshown as the yellow rectangle in Fig. 6. We can provethat the degree of the Bézier surface is (ku + kv ) ×(ku + kv + kw ), which is the accurate deformation of
the rectangle. The reason is that along the s-directionin the local coordinate system, the degree of the surface is independent of the parametric w-direction of
B-spline volume.
For a given B-spline volume, the accurate deformation results can be classified by seven cases. The
local coordinate system (s, t, n) and surface degree
are listed in Table 1. This covers all possible configurations of sub-polygons in a Normalstructure.
Till now, we have determined the degree of Bézier
surfaces and the rectangle (2D bounding box of thecoplanar sub-polygons in the s-t plane) for each
Normalstructure. The Bézier surface accompanying the sub-polygons forms the trimmed-surface representation of the accurate B-spline deformation.

3.4 Sampling and interpolation
So far, the preprocessing mechanisms have been
introduced. These can be implemented once. The

Table 1. Determining the local coordinate system and the degrees of the Bézier surface according to the normal surface
Normal nst

degrees(±1, 0, 0)(0, ±1, 0)(0, 0, ±1)(0, n y , n z )(n x , 0, n z )(n x , n y , 0)(0, 1, 0)(0, 0, 1)(1, 0, 0)(1, 0, 0)(0, 1, 0)(0, 0, 1)(0, 0, ±1)(±1, 0, 0)(0, ±1, 0)(0, n z , −n y )(n z , 0, −n x )(n y , −n x , 0)

kv × kw
kw × ku
ku × kv
ku × (kv + kw )kv × (ku + kw )kw × (ku + kv )(n x , n y , n z )(0, −n z , n y )(0, −n z , n y ) if (ku ≥ kv , ku ≥ kw )(n z , 0, −n x ) if (k ≥ k , k ≥ k )v
u vw
(n z , 0, −n x )(−n y , n x , 0)
if (kw ≥ ku , kw ≥ kv )(−n y , n x , 0)

(kv + kw ) × (ku + kv + kw )n×s

(ku + kw ) × (ku + kv + kw )(ku + kv ) × (ku + kv + kw )502

J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

work described here, which is to compute the control
points of the tensor-product Bézier surface, is processed by an interactive method.
3.4.1 Solving the functional composition of theBernstein polynomial
There are two methods of computing the control
points of the Bézier surface introduced in Sect. 3.3.3.
One solution is through the shifting operators and thegeneralized de Casteljau algorithm, which is similar to the work in [10]. It is a numerically stable but
time-consuming algorithm. Firstly, the generalized
de Casteljau algorithm is time consuming in itself.
Secondly, the B-spline volume must be converted
into piecewise continuous Bézier volumes during theinteraction, and this leads to a heavy computationalburden. Another disadvantage is that the source code
to solve the functional composition is complex when
the resultant surface is of tensor-product form. The
second solution is through polynomial interpolation,
which was used for our previous accurate FFD [11].
This solution is much faster than the first one, is suitable for real-time interaction and the source code isvery simple. We will choose it to solve the functional
composition between the rectangle and the B-spline
volume.
From a knowledge of numerical analysis, we knowthat a tensor-product Bézier surface R(s, t) of degree n s × n t can be interpolated if (n s + 1) × (n t + 1)
sampling points are given on the surface, where thesampling points are distributed in the same way as
the surface control points. The interpolation can beaccomplished by solving the linear equation system.
There are two interpolation approaches. One approach is that the sampling points are arranged as in
the 1D case. This is similar to the work for interpolating the triangular Bézier surface [11]. It is inefficient,since the equation size of tensor-product surface islarge and many of the submatrices are 0. Anotherapproach is to convert the tensor-product surfaceinterpolation problem into two stages of curveinterpolation problems. This is more efficient [9], so
we will choose the latter solution.
3.4.2 Polynomial interpolation for tensor-productBézier surface
First, the sampling points are computed on the surface. The procedure is to choose sampling points
uniformly distributed in the parametric domain of
the Bézier surface; that is, R(i/n s , j/n t ) for i =0, 1, . . . , n s and j = 0, 1, . . . , n t . In our context,
such sampling points are computed in two steps: thefirst step is to compute the sampling points in therectangle (2D bounding box in the Sect. 3.3.3) andthen compute the R(i/n s , j/n t ) from the B-spline
volume by using the de Boor–Cox algorithm. We
need not convert the B-spline volume into piecewiseBézier volumes as in [10, 11], and this results in theneed for many fewer computations.
There is a subtlety in the evaluation of the sampling
points on the B-spline volume. Because of the choiceof the local coordinate system (s, t, n) in Sect. 3.3.3,
the 2D bounding box in a Normalstructure may
not lie totally in the corresponding knot box. If wecompute the points on the B-spline volume by using de Boor–Cox algorithm directly, it may givea wrong interpolation result. The reason lies in thatthe functional composition between a rectangle anda B-spline volume limited on a knot box is equivalent to a composition between a rectangle anda Bézier volume. If the rectangle lies partly outside of the Bézier volume, we also compute thesampling points using this Bézier volume for thesake of consistency. In our algorithm, the de Boor–
Cox algorithm has a slight modification: the sampling points are always evaluated in the same knot
box index by (i,j,k) in the Normalstructure.
This can be viewed as the generalized de Boor–Cox
algorithm.
The two-step polynomial interpolations should bestarted from a parametric direction whose degree islower than a corresponding one. Without loss of generality, suppose n s ≤ n t . Then, for j = 0, 1, . . . , n t ,we should solve the following linear equations:
 B




0,n s (0)B0,n s n1s..
.

B0,n sns
ns

B1,ns (0)B1,ns n1s..
.
B1,ns nn ss

· · · Bn s ,ns (0)   R0, j· · · Bn s ,ns n1s  R
1, j


..
..

 ... .
.
n s , jR· · · Bn s ,ns n s
ns



R

j
0
ns , ntR 1 , j  ns nt 


=
.
..


.




R

(3)

j
ns
ns , nt

The left-hand (n s + 1) × (n s + 1) matrix in Eq. (3) isconstant for a given degree n s , which is independentof the underlying curve and is only related to the de-J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

503

Fig. 7. The degree of the B-spline volume is 2 × 2 × 2. The number of control points is 4 × 4 × 4

gree. Let this be noted as Bk for degree k. Obviously,
this can be inverted. The inverse matrix of Bk canbe precomputed. In our implementation, B−1
k is evaluated for k ≤ 12 in advance. The reason for taking
degree 12 is that we seldom use a B-spline volume
ij arewhose degree is greater than 4 × 4 × 4. After R
obtained, the control points Rij can be obtained for
i = 0, 1, . . . , n s by:Ti,nt T . (4)i,0 R
i,1 · · · R
Ri,0 Ri,1 · · · Ri,nt = B−1
Rnt
Combining with the sub-polygons, the trimmed
tensor-product Bézier surfaces are the accurate deformation result.

4 Implementation, results anddiscussions
4.1 Implementation and comparison results
We have implemented the proposed algorithm on
a PC with an AMD Athlon 1 GHz CPU. OpenGL
was adopted as graphics API. The trimmed NURBS

surfaces can be displayed directly by using the functions in the OpenGL Utility Library. In each of
Figs. 7–13, there are 8 subfigures, and these are as
follows:• The original object with the initial B-spline volume lattice.
• The traditional free-form deformation result without additional sampling.
• The free-form deformation result where the object is uniformly subdivided.
• The previous object subdivision method [10]. The
subdivision results in triangular meshes.
• The previous accurate deformation result [11].
The results are triangular Bézier surfaces that correspond with the triangles in the 4th subfigure• The object subdivision result by using the proposed method. The coplanar sub-polygons ina knot box are displayed by using the same color,which will be deformed as the trimmed surfaces
from a “large” Bézier surface.
• The “large” Bézier surface before being trimmed.
• The final result, the trimmed Bézier surfaces.

504

J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

89Fig. 8. The degree of the B-spline volume is 2 × 2 × 1. The number of control points is 4 × 4 × 2
Fig. 9. The degree of the B-spline volume is 2 × 2 × 1. The number of control points is 5 × 5 × 2

J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

505

Fig. 10. The degree of the B-spline volume is 2 × 2 × 2. The number of control points is 6 × 8 × 6

Of course, the final result in the 8th subfigure is thesame as the 5th subfigure except for their mathematical expression. A comparison of the numbers of
sub-polygons, the numbers of Bézier triangles andtensor-product Bézier surfaces and the run-time between the proposed method and our previous method[11] is given in Table 2.

4.2 Discussion
From implementation comparisons, we derived thatthe proposed algorithm is faster than the algorithmin [11]. The run-time acceleration results from twofacts. Firstly, the number of tensor-product Bézier
surfaces is much less than that of the triangular

Table 2. Implementation comparisonsFig.

Degree of
B-spline
volumeLattice

Triangles
in original
objectSub-polygons
number
Old:NewBézierpatches
Old:New

Run-time (s)comparison
Old:New

7
8
9
10
11
12
13

2×2×2
2×2×1
2×2×1
2×2×2
2×2×2
2×2×2
2×2×2

4×4×4
4×4×2
5×5×2
6×8×6
6×9×3
6×5×6
7×3×6

12
56
200
260
178
96
956

60 : 36
142 : 96
518 : 344
1436 : 840
814 : 438
840 : 390
1980 : 1424

60 : 24
142 : 47
518 : 157
1636 : 318
814 : 194
840 : 207
1980 : 1154

0.071/0.000
0.070 : 0.010 = 7.0
0.261 : 0.05 = 5.22
1.442 : 0.060 = 24.03
0.781 : 0.030 = 26.03
0.751 : 0.051 = 14.73
1.722 : 0.671 = 2.57

506

J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

Fig. 11. The degree of the B-spline volume is 2 × 2 × 2. The number of control points is 6 × 8 × 3

Bézier surfaces in [11], and secondly, the resultant
tensor-product Bézier surface has a low degree. Of
course, the run-time acceleration times are also influenced by the initial object configuration. If theobject has many polygons that are parallel to thecoordinate planes or if it has many coplanar polygons, then run-time acceleration will be evident, as
in Figs. 10 and 11. Otherwise, acceleration will be

J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

507

Fig. 12. The degree of the B-spline volume is 2 × 2 × 2. The number of control points is 6 × 5 × 6

less evident, as in Fig. 13. The storage comparison
is approximately equal to the run-time comparison, since the polynomial interpolation method isadopted for both algorithms. The number of sampling points is equal to the number of control pointsof the Bézier surfaces or Bézier triangles since matrix multiplication can be omitted, compared withthe evaluation of the sampling points on the Bspline volume. Thus we can say that the storage
cost of the proposed algorithm is also lower thanthat of the previous algorithm [11]. The third advantage of the proposed algorithm is that the deformation result is a trimmed tensor-product Bézier,where the trimmed surface description is consistent with that in the industrial standard, STEP. Thus,
the algorithm can be integrated into current geometric modelling systems and computer animation
systems without any difficulty. Finally, the sourcecode to solve the functional composition between
the rectangle and the B-spline volume is very simple compared with the generalized de Casteljau
algorithm.
Real-time interaction is always beneficial in computer graphics. For a simple object, the proposed
algorithm is fast enough to achieve real-time editing. For interactively deforming a complex object, however, it may not be such a good choice,since it is time consuming compared with the tra-ditional FFD method. For such cases, FFD acting
directly on the vertex should be adopted for theinteraction. If the deformation is far away fromthe real deformation result, uniform subdivision
for a polygonal object can be conducted to solve
the sampling problem. An accurate FFD of thepolygonal object is recommended for the final deformation result. We chose this strategy in our
implementation.

5 Conclusion
In this paper, a new accurate B-spline free-form
deformation algorithm for a polygonal object is proposed. First, a new initial B-spline volume definitionis given. Its Jacobian is an identity matrix. Then theobject is subdivided in the parametric space of the
B-spline volume. The sub-polygons so obtained arefree of degenerate sub-polygons. After subdivision,
the coplanar sub-polygons are combined, and thedegree of the resulting surface is then determined. Finally the tensor-product Bézier surface is evaluatedthrough polynomial interpolation, and the deformation results in the trimmed tensor-product Bézier
surfaces.
The algorithm has advantages in both run-time andstorage over the previous one [11]. The algorithm

508

J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

Fig. 13. The degree of the B-spline volume is 2 × 2 × 2. The number of control points is 7 × 3 × 6

can also be integrated into geometric modellingand computer animation systems without difficulty,
since the deformation results are the trimmed tensorproduct Bézier surfaces, whose representation isin accordance with that of the industrial standard,
STEP. The display of the resultant surface can be
accelerated if the graphics API is supported byhardware.
Acknowledgements. The work is supported by National Natural Science
Foundation of China (No.69903008) and the Scientific Foundation ofChina for Innovative Research Groups (No.60021201). The authors
also thank Dr. Neil Dodgson from the Computer Laboratory at the University of Cambridge, who supplied us with the reference [21]. We also
thank the reviewers who gave helpful suggestions on improving thequality of the paper.

References
1. Bechmann D, Bertrand Y, Thery S (1997) Continuous
free form deformation. Comput Networks ISDN Syst
29(14):1715–1725
2. Bechmann D (1998) Multidimensional free-form deformation tools. In: de Sousa A, Hopgood B (eds) State of the art
report (Eurographics ’98). Blackwell, Lisbon, pp 102–110
3. Chadwick JE, Haumann DR, Parent RE (1989) Layered
construction for deformable animated characters. ACM
Comput Graph (SIGGRAPH ’89) 23(3):243–252
4. Coquillart S (1990) Extended free-form deformation:a sculpturing tool for 3D geometric modeling. ACM Comput Graph (SIGGRAPH ’90) 24(4):187–193
5. Coquillart S, Jancene P (1991) Animated free-form deformations: an interactive animation technique. ACM Comput
Graph (SIGGRAPH ’91) 25(4):23–26

J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

6. Davis OR, Burton RP (1991) Free-form deformation
as an interactive modeling tool. J Imaging Technology
17(4):181–187
7. DeRose T (1988) Compositing Bézier simplex. ACM Trans
Graph 7(3):198–221
8. DeRose T, Goldman R, Hagen H, Mann S (1993) Functional composition algorithms via blossoming. ACM Trans
Graph 12(3):113–135
9. Farin G (1990) Curves and surfaces for computer-aided geometric design, 2nd edn. Academic Press, New York10. Feng J, Heng P, Wong T (1998) Accurate B-spline freeform deformation of polygonal objects. J Graph Tools
3(3):11–27
11. Feng J, Peng Q (2000) Accelerating accurate B-spline
free-form deformation through polynomial interpolation.
J Graph Tools 5(1):1–8
12. Gain J, Dodgson N (1999) Adaptive refinement and decimation under free-form deformation. Eurographics UK ’99,
Cambridge. http://people.cs.uct.ac.za/∼jgain/publications/EUUK99.ps.gz. Cited 17 July 2002
13. Griessmair J, Purgathofer W (1989) Deformation of solids
with trivariate B-Splines. In: Hansmann W, Hopgood FRA,
Strasser W (eds) Proceedings of Eurographics’89. Elsevier
Science Publisher, North-Holland, pp 137–148.
14. Gudukbay U, Ozguc (1990) Free-form solid modeling using deformation. Comput Graph 14(3/4):491–500
15. Hsu W, Hughes J, Kaufman H (1992) Direct manipulation on free-form deformation. ACM Comput Graph (SIGGRAPH’92) 26(2):177–184
16. Karla P, Mangli A, Thalmann NM, Thalmann D (1992)
Simulation of facial muscle actions based on rational freeform deformation. Comput Graph Forum (Eurographics’92) 2(3):59–69

509

17. Lamousin H, Waggenspack W (1994) NURBS-based freeform deformation. IEEE Comput Graph Appl 14(9):59–65
18. Liu W, Mann S (1997) An optimal algorithm for expanding the composition of polynomials. ACM Trans Graph
16(2):155–178
19. MacCracken R, Joy K (1996) Free-form deformation withlattice of arbitrary topology. In: Rushmeier H (ed) SIGGRAPH ’96 Conference Proceedings. Addison-Wesley,
Reading, Mass., pp 181–183
20. Moccozet L, Thalmann NM (1997) Dirichlet free-form deformations and their application to hand simulation. In:
Thalman NM, Thalmann D (eds) Proceedings of ComputerAnimation’97. IEEE Computer Society Press, Geneva, pp
93–102
21. Nimscheck UM (1995) Rendering for free-form deformations Technical report 381, Computer Laboratory, University of Cambridge, Oct. 1995
22. Parry S (1996) Free-form deformations in a constructive
solid geometry modeling system. PhD thesis, Department
of Civil Engineering, Brigham Young University, Utah,
Apr. 1986.
23. Sederberg T, Parry S (1986) Free-form deformation of solid
geometric models. ACM Comput Graph (SIGGRAPH’96)
20(4):537–541
24. Vergeest SM (1991) CAD surface data exchange using
STEP. Comput-Aided Des 23(1):269–281

Photographs of the authors and their biographies are given on
the next page.

510

J. Feng et al.: B-spline free-form deformation of polygonal object as trimmed Bézier surfaces

J IEQING F ENG is an associate professor at the State Key
Lab. of CAD&CG, Zhejiang
University, People’s Republicof China. He received his BSc
degree in applied mathematics
from the National University of
Defense Technology in 1992,
PhD in computer graphics from
Zhejiang University in 1997. His
research interests include space
deformation, computer-aided geometric design and computer
animation.

X IAOGANG J IN is a professor of the State Key Lab. of
CAD&CG, Zhejiang University,
People’s Republic of China. He
received his BSc degree in computer science in 1989, MSc and
PhD degrees in applied mathematics in 1992 and 1995 respectively, all from Zhejiang University. His research interests include implicit surface modeling,
space deformation, computer animation and realistic image synthesis.

T OMOYUKI N ISHITA isa professor in the Department
of Information Science at University of Tokyo, Japan, since1998. He taught at Fukuyama
University from 1979 to 1998.
He was an associate researcher
in the Engineering Computer
Graphics Laboratory at Brigham
Young University from 1988
to 1989. His research interests
center on computer graphics, including lighting models, hiddensurface removal, and antialiasing. Nishita received his BE, ME
and PhD in engineering in 1971, 1973, and 1985 respectively,from Hiroshima University.

Q UNSHENG P ENG is a Professor of computer graphics at
Zhejiang University. His research interests include realistic
image synthesis, computer animation, scientific data visualization, virtual reality, featuresmodeling and parametric design.
Professor Peng graduated from
Beijing Mechanical College in1970 and received a PhD fromthe Department of Computing
Studies, University of East Anglia, in 1983. He is currently the
Director of the State Key Lab. of
CAD & CG at Zhejiang University and serves as a member ofthe editorial boards of several Chinese journals.