The deformable surface that we use is based on the Generalised Biquadratic B-Spline (GBBS) developed by Loop and De Rose [ 10 , 11 ] for computer graphics. It is a difficult task to compute the surface due mainly to the difficulty of dealing with multi-indices of Bezier simplices.
The GBBS is a powerful and elegant generalisation of the Biquadratic
B-Spline. It automatically maintains
continuity. The GBBS is defined by a set of
M
3D control points
together with connectivity information. The connectivity defines a mesh
topology which is restricted to 4-sided faces, (see for example figure
1(b), 3(d)). Each vertex gives rise to an
n
-sided surface patch where
n
is the number of 4-sided faces using that vertex. Patch
m
depends only on the 2
n
+1 element control vector
consisting of
and the set of all vertices on adjacent faces, i.e. with
k
in the neighbourhood of
m
,
N
(
m
).
In this paper we introduce a new matrix-based scheme to compute the
surface based on notation similar to that of [
2
]. The principal steps in computing the surface are as follows. The
control vector
is converted to the Sabin net vector
by
. The Sabin net vector is converted to a vector of Bezier control points
by another matrix
. This is combined with the Bezier polynomials
to compute the point. Thus we obtain the surface patch
as a mapping from points
p
=(
u
,
v
) contained in a regular
n
-gon Domain
to a 3D surface
The whole surface
S
is the union of the patches
,
. The control vector for patch
m
,
can be obtained from the vector of all control points
by a connectivity matrix
.
The remainder of this section is of a technical nature. We show how to recast the generalised B-spline into the standard form used in computer vision. Readers interested only in how to use the surface may skip to the next section.
B-Splines are intimately related to Bezier curves [
14
,
5
]. When discussing Bezier curves it is useful to replace the usual
single parameter
u
with two parameters
and
and a constraint that
. A depth
d
Bezier curve
is defined in terms of
d
+1 control points
and the Bernstein-Bezier polynomials
as follows
The Bezier curve admits an elegant generalisation called a B-form that
maps a [(
k
+1)-variate]
k
-dimensional parameter space onto any number of scalars. Firstly we must
define multivariate Bernstein-Bezier polynomials. For these we will need
a notation for multi-indices
. The symbol
denotes a multi-index whose components are all zero except for the
j
component which is 1. It is useful to define a modulus of a multi-index
as
. The k-variate depth
d
Bernstein-Bezier polynomials are a simple generalisation of equation (
2
).
The Loop and De Rose scheme is based on S-patches, which are n-sided
smooth patches which map a point
p
=(
u
,
v
) inside
n
-sided domain polygon
D
to a 3D surface. Firstly we form
n
barycentric variables
defined as follows. Define the
n
vertices of the regular
n
-gon as
. Define the fractional areas
as the area of a triangle enclosed by points
p
,
, and
divided by the total area of
D
. Now form
n
new variables
by
Then form normalised variables
The S-patch is now simply defined in terms of the variables
and the Bezier control points
. It is a mapping
where
is a
n
-sided domain polygon and
Note that the n -sided patch uses a k +1 variate Bezier polynomial where n = k +1. The fact that a different dimension multi-index is required for each different sized patch is a bookkeeping headache for implementation.
The Bezier control points for a depth 5
n
-sided S-patch are computed from the Sabin net control points,
v
,
, and
. Letting
, then up to symmetries the Bezier control points are for
i
=1..
n
We now observe that this set of equations is a linear transform from the
Sabin net control points to the Bezier control points. Thus the
coefficients from the above equations may be used to populate the matrix
. (We note that some minor additional steps are necessary since this
leaves some matrix elements unspecified. Consult page 353 of [
11
] for details.) The Sabin net control points are obtained from the GBBS
control points by the recipe that
v
is the central vertex, the
are centroids of each face, and the
are centroids of each edge. This linear mapping can easily be converted
to the matrix
. If we rewrite equation
6
in matrix form as
then we finally obtain
as required.
The constructive procedure specified above is lengthy and complex. We
firstly deal with the indexing problem. The multi-indices
need to be converted to a one dimensional indexing scheme. Each index
may be interpreted as a number in base
d
+1, thus descending order in this value defines a unique mapping to a 1D
array. Enumeration of constant depth 1D array position to multi-index is
straightforwardly done through a lookup table, but the reverse lookup
table from multi-index to 1D array position is inconveniently large
since it is of size
, (recall
d
=5,
n
= no. of sides). However the need for this table can be avoided.
Next we must actually compute the sum specified in equation ( 6 ). Instead of the naive implementation we use the technique of de Casteljau depth reduction [ 5 ]. This is a procedure by which a set of Bezier control points are reduced iteratively in depth to d =0 which corresponds to a point. (This is in fact the basis for the geometric interpretation of Bezier curves!) For the Bezier curve of equation ( 2 ) we note that
where
. This procedure generalises to a Bezier simplex. Thus we don't ever
need to compute the Bezier polynomials, we simply depth reduce the
control points recursively until
d
=0. The remaining trick to convenient computation is a lookup table with
the 1D array positions for depth reduction. In this way multi-indices
need never actually be used. This scheme is based on the suggestions of
Shoemake [
15
].
We have shown that we can compute points on the GBBS in the form
. We now address the speed issues. There is no doubt that the
constructive procedure described above is between one and two orders of
magnitude slower than the equivalent tensor product B-Spline form. This
is partly because the number of Bezier control points (
) grows rapidly with
n
. The solution is to note that the vector
is only 2
n
+1 elements long and needs to be computed only once for each point
p
on the patch with fixed (
u
,
v
) coordinates.
When it is necessary to render the patch we convert the patches by
subdivision into triangles and the same points
p
are always used. Thus lookup tables may be created in advance containing
, and the computation for any point on the patch is reduced to 2
n
+1 multiply and add operations.
We minimise a cost
which balances a smoothness term
and a data force
. The smoothness force may be written
.
In our previous work the data force attached directly to the control points rather than the surface itself. This is unsatisfactory. We now show how to obtain a force which attaches directly to the surface itself.
We suppose that the force arises from a volumetric potential energy term
defined everywhere in space. All points on the surface should seek to
minimise this potential, i.e. we wish to minimise
The integral is approximated as a discrete sum.
With the new matrix form of GBBS the data force becomes easier to evaluate.
The total force
is made up of the sum over sample points
and patches
m
.
The smoothness force that we use is very similar to that obtained by
for curves. The force applied to each control point is equal to the
centroid of its edge neighbours less itself. A reaction force is applied
to the neighbours so that the overall force is always zero. Thus we
solve the dynamic equation
After convergence
and a local minimum of
E
is achieved.
Andrew Stoddart
