Snake Pedals: Compact and Versatile Geometric Models
with PhysicsBased Control *t
Baba C. Vemuri and Yanlin Guo3
Abstract
In this paper, we introduce a novel geometric shape modeling scheme which allows
for representation of global and local shape characteristics of an object. Geometric
models are traditionally well suited for representing global shapes without local detail.
However, in this paper we propose a powerful geometric shape modeling scheme which
allows for the representation of global shapes with local detail and permits model
shaping as well as topological changes via physicsbased control.
1 Ir proposed modeling scheme consists of representing shapes by pedal curves and
surfaces pedal curves/surfaces are the loci of the foot of perpendiculars to the tan
gents of a fixed curve/surface from a fixed point called the pedal point. By varying
the location of the pedal point, one can synthesize a large class of shapes which exhibit
both local and global deformations. We introduce physicsbased control for shaping
these geometric models by letting the pedal point vary and use a snake to represent
the position of this varying pedal point. I Ir model dubbed as a I .i. pI ..I' allows
for interactive manipulation via forces applied to the snake. We develop a fast numer
ical algorithm for shape recovery from image data using this novel geometric shape
modeling scheme. I I algorithm involves the Levenber: i.l.i p l.ii it (LM) method in
the outer loop for solving the global parameters and the Alternating Direction Implicit
(ADI) method in the inner loop for solving the local parameters of the model. I IrI
combination for the global and local scheme leads to an efficient numerical solution to
the model fitting problem. We demonstrate the applicability of this modeling scheme
and efficiency of the model fitting methods via examples of shape synthesis and shape
estimation from real image data.
Index Terms Geometric i'.dels, Snakes, Pedal Curves/Surfaces, Alternating Direction
Implicit method, Levenber '. i, rquardt method, Schur Complement.
*This research was supported in part by the NSF grant IIS9811042 and the NIH grants R01RR13197
and RO1LM05944.
submitted to the IEEE Trans. on PAMI
1B. C. Vemuri is with the Department of Computer and Information Sciences, University of Florida,
Gainesville, FL 32611 (email: 11 !!L,, ,. [t!! 1 tl I
Y. Guo is with the Sarnoff Labs in Princeton (email: '.. .. !!
1 Introduction
.i modeling shapes is a fundamental constituent of computer vision as well as computer graph
ics. Tli, 1, are numerous shape modeling techniques in computer vision and computer graph
ics literature. In the context of shape modeling, there are two primary tasks of importance in
computer vision namely, shape recovery from image data and shape recognition. Tli, shape
recovery task imposes the requirement that the shape model be capable of representing fine
detail and hence requires that the model have large number of degrees of freedom. On the
other hand, shape recognition requires that the model representation be compact in terms
of the number of parameters so as to facilitate, easy storage of a large database of models
and easy matching for retrieval. Tli , two requirements are potentially conflicting and de
signing shape models that can possibly '.1i" these requirements has been the focus of
many researchers in the recent past in the computer vision community [1, 2, 3, 4, 5].
Geometric models (lumped parameter models) can represent a shape with very fetw param
eters and hence are good for object recognition. On the other hand, physicsbased models
(distributed parameter models) require a large number of parameters for their representation
and are well suited for shape recovery and reconstruction. In [5], Terzopoulos and :i i, .1,.
developed a hybrid dynamic model a deformable superquadric model which combines the
descriptive power of geometric and physicsbased models via superposition of a membrane
spline on a core superquadric. This hybrid modeling scheme allows for modeling shapes at
both ends of the shape spectrum (one end of which is occupied by lumped parameter mod
els and the other by distributed models). Superimposing the free form deformation (FFD)
on a core superquadric allows representation of a variety of global shapes with local detail
as was reported in [1, 6, 7, 8]. This modeling scheme requires only a few parameters to
represent global shape deformations. In [9], Vemuri and Radisavljevic developed a multires
olution, hybrid modeling scheme which represents the displacement function of a deformable
superquadric model in a wavelet basis. Pentland and Williams [4] reported a novel shape
modeling scheme which is based on modal analysis. Tli, number of parameters required for
this representation depends on the desired level of detail.
In [2, 10], a hybrid hyperquadric model was developed which adds an exponential term
to a hyperquadric primitive. Tlii addition of the exponential term to the hyperquadric
equation allows for local control of the shape generated. This modeling scheme allows the
synthesis of an arbitrary number of concavities. It differs from most of the schemes described
above because the models are geometric but allow for global as well as local deformations.
Ti, model representation is quite compact in terms of the number of parameters needed to
describe the shape. However, the computation of local concavities is nontrivial and human
interaction is needed in the model fitting.
Recently, O'Donnel et al. [11] developed a novel solid shape modeling scheme. This model
includes builtin offsets from a core base model. Local deformation over the scaledoffset
model is allowed via displacement vectors attached to the default model. '.ldel fitting to
data is achieved in a dynamic framework as in [5]. This modeling scheme is well suited for
tasks such as cardiac motion recovery from 1.Ii, d `.li; images.
An interesting and topologically adaptable deformable model was developed in [3] by De
Carlo and .1 ii,\,i. via the use of parameterized blending functions. T111, models can
determine the number and position of the holes in an unknown object and accordingly
change the topology by use of appropriate error of fit criteria. ;i'li, [ I l. techniques based
on the concept of levelsets have been reported on literature [12, 13]. Another modeling
scheme called "topologically adaptable ii..!" was introduced in [14] by :.1 iii, 11 ,y and
Terzopoulos. This novel version of snakes allows for automatic change in topology if and
when required during shape recovery. In all these models, a large number of parameters are
required to describe the recovered shape.
.Idel fitting to desired features in the image data can be posed as a nonlinear minimiza
tion problem and solved numerically. In the process of numerically fitting the model to the
data, Vemuri and Radisavljevic [9] used a multiresolution gradient descent (GD) method,
while Terzopoulos and '.1% i1.x1. [5] used a first order Euler scheme to fit their dynamically
deformable superquadric model to the data. Both of these methods are slow since only a
small time step can be taken in each iteration in order to get a stable solution. A relatively
faster and more stable integration scheme based on the RungeKutta method was developed
by DeCarlo and :'.1 1..1, [3]. In [15], Cohen introduced a twostep iterative scheme based
on auxiliary variables. This scheme is very general and can be applied to most optimization
problems. Ti, auxiliary variables may be interpreted as representing intermediary recon
struction steps. :'.iii1ii.,11i ii, of several nonconvex potentials were transformed to convex
minimization problems. However, the transformation was shown primarily for potentials
involving a function of distance to the closest data and not imagebased potentials.
1.1 Overview of Our Model and Numerical Algorithms
In this paper, we propose a powerful geometric shape modeling scheme which permits the
representation of global shapes with local detail using only a , iiiglobal" characterization
and relatively simple and efficient numerical methods. Ti, proposed modeling scheme con
sists of representing shapes by pedal curves and surfaces pedal curves/surfaces are the loci
of the foot of perpendiculars to the tangents of a fixed curve/surface from a fixed point called
the pedal point. By varying the location of the pedal point, one can synthesize a large class
of shapes which exhibit both local and global deformations. We introduce physicsbased
control for shaping these geometric models by letting the pedal point vary and use a snake
to represent the position of this varying pedal point. Ti, model dubbed as a iil.I pedal"
allows for interactive manipulation via forces applied to the snake. TIl key contribution
of this paper is that our modeling scheme is inherently geometric but allows for local and
global deformations as would a physicsbased model. Also, the model is compact, capable
of global deformations such as bending, twisting, and tapering without the introduction of
additional parameters. In addition, the modeling scheme allows for synthesis of objects of
varying topology.
To fit our model to a sparse set of data points in 2D/3D as well as to image data for shape
recovery, we develop a new nonlinear optimization scheme which contains the Levenberg
:.irquardt (1 '.I) method in the outer loop for estimating global parameters of the model
and the alternating direction implicit (ADI) method in the inner loop for estimating the
local parameters of the model. Such a combination of a global and a local scheme is quite
stable and yields a locally optimal solution very quickly.
1.2 Organization of the Paper
Ti, rest of the paper is organized as follows. In Section 2 we give the definitions of pedal
curves/surfaces and show examples of the same in 2D and 3D. Following which we develop
the proposed models namely the Cil; pe II! and illustrate the ideas by synthesizing a
variety of shapes. In Section 3, we develop the mathematics of fitting the i,,IIl; pedals" to
real image data and present the experimental results in Section 4 followed by conclusions in
Section 5. Ti, appendix contains details on the derivation of the stiffness matrix in a finite
element (FE) framework and its representation as a Kronecker product.
2 "Snake Pedals": The Proposed Model
Tlii definitions of pedal curves and surfaces presented here maybe found in any book on
elementary differential geometry text and we refer the reader to [16].
2.1 Pedal Curves
Let a be a planar curve, the pedal curve of a is defined as the locus of points on the foot
f of the perpendicular from a fixed point p called the pedal point to a variable tangent of
the planar curve a.
Suppose that the planar curve a is a parameterized curve with domain parameter t. Let 3
Figure 1: f is on the pedal curve of a with respect to the pedal point p.
be the pedal curve of a with respect to the fixed pedal point p, and let a(t) = g, 3 (t) = f
at a particular t, as shown in Fig. 1. T l, projection of a (t) p in the direction Ja '(t)
must be ( (t) p, where a' (t) is the tangent line of plane curve a (t), J : R2 R i2 is a
linear map given by
J(,y) =(y,). (1)
Note that applying the operation J on a vector results in another new vector and J can be
geometrically interpreted as a rotation by 7r/2 in a counterclockwise direction. We can thus
define a pedal curve as follows [16]:
Definition : T.,, pedal curve of a regular curve a : (c, d) R2 with respect to a fixed
(pedal) point p S 2 is given by
pedal(t) = p + ( p) J (t). (2)
 Ja 1(t)
In Fig. 2, we present examples depicting the pedal curves of an ellipse for different positions
of the pedal point as shown. Note that the pedal curve is capable of exhibiting global as
well as local deformations. Tli, global shape of a pedal curve of an ellipse is approximated
by a new ellipse, whose long axis is almost the same as that of the original ellipse, while
(a) (b) (c) (d)
Figure 2: Examples of pedal curves (in green) of an ellipse (in red) for different pedal point
positions. Pedal points (in red) are shown by a dot in each case.
the short axis is much larger than that of the original ellipse. Ti, new ellipse is pinched in
the direction of the short axis to generate the pedal curve. Ti, location and extent of this
pinching are determined by the position of the pedal point. As shown in Fig. 2, if the pedal
point is located in the center of the original ellipse (Fig. 2 (a)), the pinching occurs around
the center part of the pedal curve, while as the pedal point moves to an upper position (Fig.
2 (b)), the pinching in the pedal curve moves up. Similarly, as the pedal point moves to the
right of the ellipse origin (Fig. 2 (c)), more pinching occurs on the right side in comparison
to the left of the origin. This pinching operation constitutes the local deformation of the
pedal curve, but the type of the local deformation that can be exhibited by a pedal curve is
by no means limited to the pinching type deformation. '. .1,re complex local deformations can
be realized in the pedal curve as shown subsequently. In summary, by moving the position
of the pedal point, it is possible to synthesize a variety of local deformations. Ti!, curve a(t)
will henceforth be referred to as the generator for the pedal curve 0 (t) and the process of
generating a pedal curve will be referred to as the pedaling operation.
2.2 Snake Pedals
1 .ire general shapes may be synthesized by letting the pedal point be different for each
point of the generator when we apply the pedaling operation. We can let the pedal points
(a) (b)
Figure 3: (a) Tli process of generating a snake pedal with an ellipse as a generator. (b)
i~1. pedal" controlled by the snake using an ellipse generator.
be specified by another curve p(t), which can be represented by a standard snake/Bsnake
[17, 18] and then apply the pedaling operation to each point on the generator ac =a (ti)
with respect to corresponding pedal point pi = p(ti), where a and pi are points on the
generator a (t) and curve p(t) respectively. In this paper, we will focus on parameterized
generator curves. Ti, pedaling operation generates a new curve that we dub a snake pedal
x(t) as shown in Fig. 3.
If the generator is an ellipse as shown in Fig. 3 (a), we can represent it in a parametric
form by
L cosO sinO a cost m1 in 1
sinO cosO b sint m2
where a, b are aspect ratio parameters, 0 is the rotation angle between the world coordinates
and the object centered coordinates, m = (ml, m2)' is the centroid of the generator in the
world coordinates. We collect the generator parameters into a vector g = (a, b, 0, mi, m2)'
and the vector g is referred to as the global shape parameter vector.
For an ellipse, the term Ja'(t) in the Eq. (2) of the pedal curve can be expressed as
Ja (t) n[ ]= al sin t sin 0 a2 cOs t cOs 0 (4)
n2 al sin t cos 0 + a2 cos t sin 0
Note that .[a'l = n + n = a, sin2t a cos2t, which is independent of the rotation angle
0. We will see that .*a'I is a very important quantity in the model fitting process.
Ti! pedal point can belong to a curve, which is not restricted to be a parameterized function.
We can sample the curve at the same rate as the sampling rate for the generator so that
there is a onetoone correspondence between this curve and the generator. We use a snake
as this controlling curve and get a sequence of position vectors, p; = (xi, ;, )', by sampling
the snake. We collect all the snake position vectors in another vector, p = (p', p ,..., p)
where 1[ is the number of sampling points. Note that the snake may be replaced with a
Bsnake to achieve compactness in the representation. Vector p is referred to as the local
shape parameter vector. A snake pedal curve is completely determined by the global and
local shape parameter vectors.
In Fig. 4, we depict some examples of snake pedals, which are generated using the snakes
and an ellipse as the generator. Note the variety of local deformations that can be generated
using this modeling technique. We remind the reader that the snake pedal itself is a geometric
model and that it is not directly responsive to the application of external forces unlike the
standard snake models [17]. Also, it is easy to synthesize snake pedals whose topology is
different from that of the generator as will be seen in examples depicted subsequently. TI
uniqueness of our modeling scheme lies in the ability to generate snake like local deformations
in a non reactive geometric shape namely the pedal curve/surface.
Tli, pedal curve definition can be modified slightly by subtracting the second term from
the first term in Eq. (2). This modification allows us to synthesize a large class of shapes
which appear very different from the shapes produced using Eq. (2). Til, key feature
of using this modification is that the snake pedal curve allows for more local deformations
including shrinkage and expansion of the snake pedal. Tli, shrinkage and expansion behavior
is controlled by the location of the snake. When the snake is located inside the generator,
the snake pedal exhibits shrinkage behavior, while expansion is exhibited when the snake is
outside the generator. Tli, following equation describes this modification:
(a) (b) (c)
Figure 4: Examples of i., ;I pel .l!" using an ellipse generator and a snake.
(ca(t) p) J '(t)(5)
pedal(t) = p J )(tM. (5)
Fig. 5 (a)(c) depict examples of shapes generated using this modification. Note the amount
of local deformations allowed in the snake pedal curves using this modification. We compare
the original and i i lio1i I" snake pedals using the same generators and same snakes in Fig.
5 (d)(e).
Topology changes maybe realized by applying the pedaling operation to a single generator
and several pedal points (or snakes). Fig. 6 depicts an example of the topology change
realized in 2D. When the snake is located inside the generator, a single closed snake pedal
curve is generated, when the snake is split into two parts with both parts located outside
the generator, two separated closed snake pedal curves are obtained. In the later case, only
part of the generator is used to produce the snake pedal curves. Note that any number
of closed snake pedals can be generated from one generator using this procedure if the
position of snakes and parts of the generator are properly chosen. T11 , quantities may be
automatically determined in data fitting applications using an appropriate technique.
(b)
Figure 5: (a) (c) Examples of modified nil.l;' pedals" using an ellipse generator and a
snake. (d)(e) Comparisons of original(red) and modified (blue) il.;h, pe l.i! using the
same generators and same snakes in both cases.
J I
I
JI
(a)
I I
I I
0
(b)
I
(C)
Figure 6: Examples of topology change. (a) depicts a single snake pedal curve (solid)
obtained from an ellipse generator (dashed) and a single pedal point. (b) same as (a) except
two distinct pedal points located outside the generator are used. (c) same as (a) except
several pedal points on a snake are used. (d) same as (c) except a discontinuous snake is
used.
(C)
2.3 Pedal Surfaces
A pedal surface is the surface analog of the pedal curve. It is the locus of the points on the
foot of the perpendicular from a fixed pedal point to a variable tangent plane of the surface.
1, i'e precisely, it is defined as follows [16]:
Definition : T.,, pedal surface of a parameterized surface or a patch a : U + R3 with
respect to a point p E S3 is defined by
(a(u, v) p) (a, x a0,)
pedal(u, v) = p + ( )a~U x av, (6)
au X v
where U is an open subset of S2, aQ(u, v) is a parameterized generator surface with parameters
u and v, a, and av are tangent vectors in the u, v directions.
Fig. 7 depicts examples of pedal surfaces synthesized using an ellipsoidal generator surface
and a fixed pedal point located inside the generator surface in each example. Note that
change in topology is achievable by simply choosing an appropriate location of the pedal
point as shown in Fig. 7. However, this type of topology change is different from what we
described in Fig. 6.
As in the 2D case, we can let the pedal point vary for each point on the generator surface.
Tliil. we have the snake pedal surfaces in 3D whose shape can be controlled by snakes which
are either curves or surfaces in 3D. For examples of the 3D snake pedals, we refer the reader
to the next section.
In comparison to the deformable superquadric model, the snake pedal model can exhibit a
larger amount of deformation and can easily exhibit a large amount of bending, twisting and
.,l, '.,I without the explicitly using parameters to represent these operations in the model
representation. Ti, i, fi 'e, it is well suited for modeling and recovery of shapes with complex
geometries, examples of which are Rabundant in the field of medical image analysis.
T1 Y
(a) (b) (c)
(d) (e) (f)
Figure 7: Examples of pedal surfaces with fixed pedal points located inside ellipsoidal gen
erator surfaces.
3 Model Fitting
In this section we will develop an efficient technique for fitting our ~I.1Ii pedal" model to
a sparse set of data points in 2D/3D as well as to image data for shape recovery. Ti, model
fitting is posed as a nonlinear optimization problem and can be solved numerically. Til, 1,
are two types of parameters which need to be determined in this optimization process. One
is the global shape parameters namely, generator parameters g = (a, b, 0, mi, mi2)T (in 2D),
and another is the local shape parameters, i.e., snake position p = {p }1. A fast numerical
algorithm is proposed to estimate these two types of parameters. Ti, algorithm involves
employing the Levenber:'. i.Irquardt (i_ ) technique to solve for the global parameters in
the outer loop and the Alternating Direction Implicit (ADI) method [19, 20, 21, 22] to solve
for the local parameters in the inner loop (for fixed global parameters) as shown in Fig. 8.
We present the details of each of these methods in the following sections.
OL
While ( Energy > e )
One iteration of
solving global parameters
Solve local parameters
completely
Figure 8: Big picture of numerical schemes for model fitting.
3.1 LM Method for Global Parameter Estimation
Ti, estimation of the global shape parameters in the outer loop is a nonlinear optimization
problem, and the Levenber,.i, 1rquardt (1 .i1) method [23] is utilized for achieving the opti
mization. T1i, I1'. 1 method essentially combines the '; "ton method and the gradient descent
(GD) method in solving the nonlinear function minimization problem. When the current
solution is far from the true solution, the i .1 method adopts the gradient descent (GD)
algorithm. On the other hand, when the current solution is close to the true solution, the
nonlinear function can be approximated by a quadratic form and the I .1 method switches
to the '; "ton method in this case. For a more detailed description of this algorithm, we
refer the reader to [23].
We now adopt the I .1 method to our model fitting problem. To fit the snake pedal to a
given data set of m points {Di}l, in 3D, we minimize the sum of the squared distances from
each data point Di to the corresponding closest point xi on the pedal curve/surface. Let h
be the function that returns the closest point on the snake pedal from a given data point
Di i.e., xi = hi(g, pi, ti). This can be achieved by employing the technique of migration of
the point of influence as in [5]. Let {pi}1, be a set of \[ discrete points along the snake
p(t), {xi}l'l be the corresponding points on the snake pedal curve. Let D = {D } and
p = {pf}I1, then the quantity to be minimized for the model fitting can be written as
1m
ED(D, g,p) = 3(D x)2, (7)
Si=1
where xi = hi(g, pi, ti) and the magnitude of 3 is related to the uncertainty of data mea
surements. T!i, quantity ED(D, g, p) can be considered as the external energy imposed on
the model in the fitting process.
In 1 .1 the gradient of ED(D, g, p) with respect to generator parameters g = (gi, g,..)T
(a, b, 0, mi, mn2)T has the following components:
OED Thi (g, pi, ti)
 (Di hi(g, pi, ti) ). (8)
If we assume that there is a onetoone correspondence between the points on the generator
and the pedal snake then the partial derivative OED can be written as:
OED ihk (, Pktk)
E = P(Dk hk(g, Pk, tk) ). (9)
(0 1. 0.,.
Taking an additional partial derivative gives xL which is needed in the I.1 iteration.
For the derivative of the energy required in the i. algorithm, we use the chain rule to get
OEDo = (OED)(h ). . v1ire specifically, if we denote each component hi of h as hi = (hil, hi2)T'
then
[Oh il Oh1 Oh1 Oh Oh 1
OED [ (Dil hil) (Di2 Oa Ob 00 0ni1 0'
Oa Ob 90 Omn1 O,1,
(10)
After the estimation of the global parameters, we solve the local shape parameters in the
inner loop on the parameter estimation algorithm using the ADI method.
3.2 ADI Method for the Local Parameter Estimation
When solving for the local parameters in the inner loop for fixed global parameters, the
quantity to be minimized has a quadratic form and its minimization requires the solution to
a linear system, namely, Kx = b, where K is the stiffness matrix, x is the local parameter
vector (set x = p), and b = '" which is the external force described later. Because
of the special structure of the matrix K in this linear system, we can use the alternating
direction implicit(ADI) method for the local parameter estimation. We can prove that the
ADI method takes only O(N) time to converge with soft data constraints where N is the
number of nodes in the rectangular finite element discretization of the snake in our model.
Note that in soft data constraints, the model is not required to interpolate the data but is
expected to approximate it.
We now present the structure of the stiffness matrix K. In the snake pedal model, the
controlling snake p is a controlled continuity spline under the influence of image forces or
other external constraint forces. In our model, the spline deformation energy imposed on
the snake is a membrane energy plus thin plate energy given in a continuous form by
EO(p)+( OP+ a22 P2 a2 (+( a2(1
E(p) W ( )2 ( 2( )_ 2( ) du dv, 11)
where wi and .,_ are the tension and rigidity factors respectively.
T!, external energy used can be the same as the energy ED(D, g, p) discussed in Section
3.1. In this case, the external force imposed on the snake is given by
8ED 8ED 9hi 9h,
O O E3[ (Dil hi1) (Di2 hi2) O, (12)
Opi Ohi Opi pi
/South Pole
Total number of nodes :
n Xn+2
V
Denote : N = nX n
North Pole
u
Figure 9: :'.i del tessellation in intrinsic coordinates.
where
Ohi n 2 n + n n + n* (
OhP P1 1 2h 1 2 2
1_ 1 2 2
with ni, n2 defined as before. Notice that OED/Ohi is the external force in a standard snake,
but in our case, the external force is OED/Opi, which is related to the force in a standard
snake via a Jacobian matrix 0hi/Op, in Eq. (13).
In order to fit the snake pedal curve to salient features in an image, we use the following
equation as the energy to be minimized:
ED = (G, VI) (G, VI(x) , (14)
where x = x(x(t),y(t)) describes the snake pedal curve, G, is a Gaussian with standard
derivation a, and denotes the convolution operator. '. .re sophiscated boundary finding
energies can be synthesized and we refer the reader to [5]. By minimizing the above energy,
the model is made to fit to points of high gradient in the image. Ti, procedure to minimize
above energy is similar to that of fitting of snake pedal curves/surfaces to data points.
We then discretize the model in model coordinates using finite elements. Geometrically, the
snake pedal surfaces are closed surfaces in space whose intrinsic coordinates u = (u, v) are
defined on a domain Q. Ti, positions of points on the model in an inertial frame of reference
are given by a vectorvalued function pedal(u) = (Ul),(u), x 3(u)). Til tessellation of
u = (u, v) in intrinsic coordinates is shown in Fig. 9. We use a combination of quadrilateral
and triangular elements in the discretization. In the inner quadrilateral mesh, we have
(n x n) = N nodes and 2 additional nodes are used in representing the north and south pole
leading to a total of (N + 2) nodes. Due to the presence of the north and south poles, after
the discretization of the model in intrinsic coordinates u using finite elements, the stiffness
matrix K C R(N+2)(N+2) in our problem can be written as a (2 x 2) block matrix:
K K K 12 (15)
K K21 K22
where K11 C RNxN contains coefficients which corresponds to the inner quadrilateral mesh
and K22 C 2x2 is a diagonal block matrix containing the coefficients corresponding to the
north and south poles, and K12 = KT1 contains the coefficients of interaction between the
poles and other nodes.
To solve the linear system Kx = b for the local parameter estimation, we employ the
following Schur complement formula [22]
Klb Kll + K1K12S1K21K11 K1 K12S1
K _b K S1KK11 S 1 S1 b, (16)
S1K21Kl S_
where
S = K22 K21K1K12. (17)
From Eq. (16), it is evident that to obtain the solution of Klb, we need the solution of
several Kllbs. Further analysis shows that the component K11 of the stiffness matrix K
can be expressed as the sum of two Kronecker products:
K11 = wi(AB+B A)+ ,A A
= ADC+CA, (18)
where C = wiB + (.r_/2)A with A and B being (n x n) tridiagonal and symmetric positive
definite (SPD) matrices of the form
1 1 2 1
1 2 1 1 4 1
A = .. .. .. ~~nx"; B = .. .. .. C Rnxn. (19)
1 2 1 1 4 1
1 1 1 2
Tli, Kronecker product of C and A, denoted by C A, is the matrix of order N x N
(where N = n x n) and given by
c11A c12A ... clnA
c21A c22A ... C2nA
CA= "' (20)
CniA Cn2A ... c,, A
When computing with Kronecker products, it is computationally convenient for us to rep
resent vectors using matrices. When solving (C A)x, we represent the vector x of order
n2 by the matrix X = {xkt} of order n x n defined by
Xkl = Xk+n(t1). (21)
It can be shown that this representation can lead to efficient procedures for computing
(C A)x and solving (C A)x = b, respectively. Tli, following lemma is used for this
purpose.
LEMMA 1
Let C = {ckl}, A = {akl} and X = {xki} be matrices of order n x n, then the n x n matrix
representation of n2 x 1 vector (C A)x is given by
19
(C (D A)x (C(AX)T)T.
Using this representation, the linear system KZ1x
b can be written as:
CXAT + AXCT = D,
(22)
where D, X are the n x n representation of the n2 vector b and x. Since A and C are also
SPD matrices, Eq. (22) can simply be written as
CXA + AXC = D.
(23)
After some simple linear algebra operations, it follows that above equation is equivalent to
A'X + XA' = D',
(24)
with A' = A1C and D' = A1DA1
Eq. (24) is the standard Lyapunov matrix equation [24]. Ti1, ADI method can be applied
to solve this matrix equation and has the following steps in each iteration j = 1, 2,..., J,
(A' + pjI)Xji
Xj(A' + pjI)
(25)
(26)
D' Xj_i(A' pjI)
D' (A' pjI)Xj
For each iteration in ADI, we need to solve two equations, namely, Eq. (25) and (26), and
therefore Xj. is introduced as an intermediate variable. Since A' is a SPD matrix, we
j2
can use Cl .1 l.y decomposition for the tridiagonal SPD matrix (A' + pjI) to compute the
update solution. This only takes O(N) time per iteration. Tli, parameters pj are called
ADI parameters and J is the number of iterations needed. T!i ~ are chosen according to the
method described in [25].
4 Model Fitting Results
In this section, we present a set of five experiments, two of these experiments are in the
medical imaging domain and demonstrate model fitting to sparse 3D data points placed
by an expert neuroscientist along the boundaries of a hippocampus and a  i. in selected
slices of an .i: brain scan. In both the experiments, the spatial resolution of the images
was (256,256) and a region of interest was identified interactively. Sparse set of points are
placed on the boundary of the shape of interest because these shapes of interest do not
have sufficient gradients along large segments of the boundary. This requires that prior
anatomical knowledge be introduced into the model fitting process. What we describe here
in this paper is the stage of the fitting process requiring prior knowledge from an expert
in the form of a sparse set of boundary points on the shape of interest. This is used in
conjunction with image gradients as external forces in the model fitting process. Once the
model is successfully fitted, its parameters over a large number of fits can be learnt for later
use. Tli, third example depicts the snake pedal fitting in 3D to a Vulcun (Spock from Star
Wars) head. And the last two experiments demonstrate the power of the model in fitting
to 3D synthetic data from a bent shape (hair pin) and a twisted shape respectively. Once
again, we would like to stress the fact that the model does not make use of any explicit
1, iiliii /twisting transformations to recover the underlying shapes from these data.
Tli, medical data example in Fig. 10 is organized as follows: left to right, the images depict
a slice of an i, i: brain scan in which the shape of interest hippocampus or a  i i. has
been identified by a neuroscientist via sparsely placed points (only in selected slices of the
.1i: scan) on the shape boundary. TI, next image shows the collection of these 3D points
(in red) and the initialized snake pedal model (in green) followed by images depicting the
intermediate stage of fitting and the final fitted model respectively. Tli, model fitting process
takes 258s and 291s CPU time on SGIO2 respectively. As evident, the model achieves a
visually accurate fit to the data and the model initialization is quite far from the final true
(a) (b) )()
(e) (f)
Figure 10: Firl in results of snake pedal surfaces to a hippocampus (d) and a. 1 i. (h) with
initializations shown in (b) and (f), the intermediate stages are shown in (c) and (g). (a)
and (e) are slices of, 1I i1 scans for the hippocampus and  i ii
shape in each of these examples. In addition, we have verified that these fits are in agreement
with xwhat the expert neuroscientist would achieve via manual segmentation.
Ti, non medical data examples in Fig. 11 are organized as follows: left to right, the
images depict a 3D synthetic data set (in red) along with the model initialization (in green)
superimposed which is followed by images of intermediate stage of fitting and the final model
fit. Tlii CPU time for the model to converge to the desired shape spock, bent shape and
twisted shape is 267s, ;l... and '"',' on SGI02 respectively. Note that in the synthetic
data examples, in one case, the data exhibits a large amount of bending and in another, a
large amount of twist. What we would like to stress here, is the fact that the model did not
require any explicit bending or twisting transformations (unlike the deformable superquadric)
to achieve the fits. As mentioned earlier, modeling schemes that require explicit incorporation
of bending, twisting transforms are highly nonlinear and are known to introduce numerical
instabilities into the fitting algorithms. Ti, fitting times are quite large in these examples
(c) (d)
than in the previous two because of the large amount of bending and twisting present in the
desired shapes and also because the initialization was quite far from the final solution.
We performed a comparison between our model fitting and the deformable superquadric
model fits to both synthetic and medical data using a fixed CPU time (that required in
achieving the fits with our model to the . il. and the twisted shape data shown in Figs.
10 (e)(f) and 11 (g) respectively). For these fixed CPU times, the deformable superquadric
fits did not converge to the desired solution as is evident from Fig. 12 (b)(e). Figs. 12
(b)(c) depict model fitting to data shown in (a) using the deformable superquadric and the
snake pedal models respectively. A similar arrangement is depicted in Fig. 12 (e)(f) for
data shown in (d). In both deformable superquadric model and snake pedal model fitting,
the numerical techniques used were the same. As evident from these comparisons, the snake
pedal model was able to achieve visually accurate fits while the deformable superquadric
yielded results that are far from satisfactory (within the fixed CPU time).
5 Conclusion
In this paper, we introduced a novel geometric shape modeling scheme which allows for
representation of global and local shape characteristics of an object. Geometric models are
traditionally well suited for representing global shapes but not the local details. However, we
proposed a powerful geometric shape modeling scheme which allows for the representation of
global shapes with local detail and permits model shaping as well as topological changes via
physicsbased control. Our modeling scheme consists of representing shapes by pedal curves
and surfaces. Physicsbased control was introduced for shaping these geometric models by
using a dynamic spline to represent the position of the pedal point. Ti, model dubbed as a
i~i.1: pedal" allows for interactive manipulation via forces applied to the snake.
We also developed a novel technique for fitting our model to a sparse set of data points in
2D/3D as well as to image data for shape recovery. Tli, model fitting is posed as a nonlinear
(d) (e)
(g) (h) (i)
Figure 11: Filling results of snake pedal surfaces to (c) spock, (f) a bent shape and (i)
a twisted shape. Initializations are shown in (a), (d) and (g) respectively followed by the
intermediate stages in (b), (e) and (h).
(d) (e)
Figure 12: Comparison of deformable superquadric and snake pedal model fitting perfor
mance. (a) and (d) are model initializations, (b) and (e) exhibit model fits with deformable
superquadric models, (c) and (f) are model fits with snake pedal models. CPU time for fits
in (b) & (c) and (e) & (f) respectively are the same.
optimization problem and can be solved numerically. Tli, nonlinear optimization scheme
developed in this paper contains the Levenbers i .iArquardt (i .) method in the outer loop
and the Alternating Direction Implicit (ADI) method in the inner loop, and the combination
of the global and local scheme is quite stable and yields a locally optimal solution in very
few iterations.
We demonstrated the applicability of our modeling scheme and numerical algorithms via
model fitting to real image data and range data. Future work will focus introducing dynamics
into the model and experimenting with more image data.
6 Appendix
In Section 3.2, we described the internal deformation energy in an active/deformable model.
To obtain the expression for the stiffness matrix K, we employ the finite element method
(FEM) to discretize the model in material coordinates u = (u, v). As shown in Fig. 9, we
v I a "1
S(1,1) (1,1)
W1 T
(u v,) )
(1,1) 2 (1,1)
u
Figure 13: Bilinear quadrilateral element, where (u, v) are material coordinates and (w, y)
are the reference coordinates.
use bilinear quadrilateral elements in conjunction with linear triangular elements. TilI later
is used for representing the regions around the north and south poles. We will first discuss
the bilinear quadrilateral elements, and the linear triangular elements can be handled in the
similar way.
6.1 Bilinear Quadrilateral Elements
Bilinear quadrilateral elements are quite common in FEM literature [26]. In this section,
we present the basic transformation equations between material and element coordinates.
Denote the finite element nodal shape functions by Ni where i = 1, 2, ..., no, and no is the
number of nodes associated with each element Ej. For the bilinear quadrilateral elements,
we have
N(w,) = (1 + ww) (1 + ), (27)
4
where (wi, .' ) are the reference coordinates of node i as shown in Fig. (13).
T11, reference coordinates (w, ) are related to the material coordinates u = (u, v) by
2 2
w = (u Uc) = (v ), (28)
a b
where (uc, v,) are the coordinates of the element center, a and b are the length and width of
the element. In our implementation, we use a = b. Differentiating Eq. (27) with respect to
u and v gives
ON8 ONZD w N1 8
Du Dw Du 9 D u
ONi Ni D w +Ni 9D =
Ov Ow Ov 9 Dyv
When integrating function f(u, v) over element Ej,
ence system:
I IE f(u,v)dudv
1
2aW(1 + )
2a
(29)
(1 + wwi (30)
2v
we perform the integration in the refer
.4
(31)
6.2 Derivation of the Local Stiffness Matrix
As described in section 3.2, we use the membrane plus thin plate energy expressed in Eq.(11)
as the internal deformation energy imposed on the snake,
EO !+ Dp+ p2 2 + D a 2p
E,(p)= wI(( )2( )2 (( + 2( )2 2 ) du dv. (32)
Ou OV DuOv
By using the procedure described in [5], we can obtain the local stiffness matrix. For the
membrane energy part, if we denote the stiffness matrix for the element Ej by K c SR xnoo,
the entries of Km can be represented as
/ DN1,DNm D8NN D,
(k'I)wm (= Ni ( + D D ) du dv, (33)
M BE u Ov Ou 0V
where m = 1,2,...,no. Since we use a quadrilateral element, no = 4 in this case. Using
the formula in Eq. (29), (30) and (31), we obtain the following local stiffness matrix for the
membrane energy
K = (w )
m 6 a
1
2 6 bL
2
2
1
1
2
(34)
V
Figure 14: A small 3 x 4 mesh in (u, v) space.
Adopting the similar scheme, we obtain the following stiffness matrix Kj for the thin plate
energy part in element E,
1 1
KT g s
t ab 1
S1
T!l, summation of Kj, and Kj gives the full
1 1 1
1 1 1
1 1 1
1 1 1
stiffness matrix for element Ei
(35)
(36)
K' = K J+ Ki.
Note that both KN and KI are symmetric positive definite (SPD) matrixes, therefore their
summation, Kj, is also a SPD matrix.
6.3 Kronecker Tensor Product Representation of the Assembled
Stiffness Matrix
Ti, matrix Kj discussed earlier is only the local stiffness matrix for element Ej, which
describes the stressstrain relationship between all the nodes within element Ej. For
example, k2 = 4 means applying one unit displayment to node 3 results in 4 unit force on
node 2. We need to assemble the local matrix over all the elements to obtain the entire
E1 E2 E3
5 6 7
E4 E5 E6
9 10 11
wib
6a
0
wib
x
6a
0
0
02
02
x
wib
6a
00
wib
x
6a
0
0
Figure 15: '.1I iii I oI iii molecules
@@@
w1b 4x8 
x a
@@@
Figure 16: Assembled membrane molecular representation.
stiffness matrix K. To achieve this, we introduce a novel method via the Kronecker tensor
product operation. A small 3 x 4 mesh in (u, v) space as shown in Fig. 14 is used to illustrate
the idea. We will only consider the first part of the KJ (the first term in the summation of
Eq. (34)) for explanation, the second part of KJ and Kj can be handled in the similar way.
In Fig. (14), we assign a number to each node, ranging from 1 to 12, and denote each
element by El, E2, ..., E. We choose node 6 for our illustration, and investigate the stress
strain realtionship between node 6 and every other node in the entire mesh. From Eq.
(34), observe that applying one unit displacement to node 6 will result in W (b) x 1 unit,
(w )b x (1) unit, l(b) x (2) unit force on node 2, node 1, and node 5, respectively. This
concept can be visualized by the nodal molecule representation, as shown in Fig. (15) (a). We
can use a similar representation for elements E2, E4 and r , resulting in the nodal molecules
as shown in Figs. (15) (b), (c), and (d), respectively. Note that there is no interaction
between node 6 and nodes in elements E3 and Eg. From the nodal molecular representation,
when applying one unit displacement to node 6, we can obtain the force resulting on every
other node in the entire 3 x 4 mesh by summing up all the molecular values at the same
molecular location. Tli, result is shown in Fig. (16).
Ti, assembled molecular representation can be expressed as a Kronecker product formula
( [ 1 2 1 ] 4 (37)
6a I
Ti, assembled stiffness matrix that corresponds to the first part of the membrane energy
(first term in the summation of Eq. (34)) is
Karti) = AX By, (38)
where
2 1 1 2 1 0 1 
1 2 1 0
A = 1 4 1 and By 1 (39)
0 1 2 1 0 1 2
We now prove the above conclusion by expanding the Kronecker tensor product expression
in (38). According to the Kronecker product definition in Eq. (20),
2 1 2 1 0 1
K(partl) ( b) 4 1 2 1 0
m 6 a 1 0 1 2 1
0 2 1 0 1 2
2By By 0
a= ( ) By 4By By
0 By 2By
4 2 0 2 2 1 0 1
2 4 2 0 1 2 1 0
0 2 4 2 0 1 2 1
2 0 2 4 1 0 1 2
2 1 0 1 8 4 0 4 2 1 0 1
(wb 1 2 1 0 4 8 4 0 1 2 1 0
6a 0 1 2 1 0 4 8 4 0 1 2 1
1 0 1 2 4 0 4 8 1 0 1 2
2 1 0 1 4 2 0 2
1 2 1 0 2 4 2 0
0 1 2 1 0 2 4 2
1 0 1 2 2 0 2 4
(40)
In this matrix, each entry (Y ,nortl)),l represents the stressstrain relation between node I
and node rn, where I, r = 1, 2,..., 12. For example, the entries in the 6th column (. (n'rt1))16
represent the force caused in node I by the unit displacement in node 6. IWe can see from
column 6 that one unit displacement in node 6 contributes (w1 b) x 8 unit force to itself, and
6 a
(w1 b) x (4) unit force to node 5 and 7, (wV ) x 2 unit force to node 2 and 10, (w 6) x (1)
S6 a/ v \/ "' 6 a' '" 6 a '
unit force to node 1, 3, 5 and 7, and zero force to node 4, 8 and 12. This is exactly what
Fig. (16) conveys.
For the second part of the membrane energy, performing a similar procedure, we have
KI = (6w ) B, A, (41)
where
0 4 1 0 1 
1 4 1 0
A,= 1 2 1 and B= 0 1 4 1 ()
1 1 1 0 1 4
Since a = b in our implementation, the entire stiffness matrix corresponding to the membrane
energy is
K = K (arti) + K..
= (W) (A, 0 By + B 0 Ay).
Similarly, we have the stiffness matrix corresponding to the thin plate energy
K1t ( ) A, (D Ay.
ab
In general, for an m x n mesh, denote
2 1
1 4 1
A X ... .. ...
1 4
1
4 1
1 4 1
Ay 1 .. . ..
1 4
1 1
(43)
(44)
(45)
E Smxm
1
2
1
1
4
(46)
1 1
1 2 1
Bx ... ... ... R SMxm (
1 2 1
1 1
2 1 1
1 2 1
By = ... ... ... (
1 2 1
1 1 2
then the assembled stiffness matrix corresponding to the energy expressed in Eq.(11) is:
K Km + Kt = wl(Ax 0 By + Bx Ay) + ,,'_(A ( A,), (
47)
48)
49)
where w1 and .,_ are tension and rigidity parameters. Note that the difference between the
first column/row and last column/row in A, and Ay or B, and By is due to the different
boundary conditions namely, the natural boundary and the periodic boundary conditions
respectively. When employing the same boundary conditions and let m = n, A = Ax = Ay,
B = Bx = By, Eq. (49) becomes
K = wi(A B + B A) + ,,'_(A A). (50)
This is the same equation as (18).
References
[1] E. Bardinet, N. Ayache, and L.D. Cohen, "Fitting of isosurface using superquadrics
and freeform deformations," in Proc. of IEEE Workshop on Biomedical Image Analysis
WBIA '94, Seattle, 1994, pp. 184193.
[2] I. Cohen and L.D. Cohen, "Hyperquadric model for 2d and 3d data fitting," in Proc.
of ICPR, Piscataway, NJ, 1994, pp. 403405.
[3] D. DeCarlo and D. :'. i1 .~i\ "Blended deformable models," IEEE Trans. Pattern
Analysis and Machine Intelligence, vol. 18, pp. 443448, 1996.
[4] A. Pentland and J. Williams, "Good vibrations: :., del dynamics for graphics and
animation," Computer Graphics, vol. 23, no. 3, pp. 215222, 1989.
[5] D. Terzopoulos and D. '.1 ii, \., "Dynamic 3d models with local and global deforma
tions : Deformable superquadrics," IEEE Pattern Analysis and Machine Intelligence,
vol. 13, no. 7, pp. 703714, 1991.
[6] '. Hebert, K. Ikeuchi, and H. Delingette, "A spherical representation for recognition
of freeform surfaces," IEEE Pattern Analysis and Machine Intelligence, vol. 17, no. 7,
pp. 681690, 1995.
[7] W.:'.1 Hsu, "Direct manipulation of freeform deformation," Computer Graphics, vol.
26, no. 2, pp. 177184, 1992.
[8] T.W. Sederberg and S.R. Parry, i I t I,, deformation of solid geometric models," in
ACM SIGGRAPH, Dallas, TX, 1986, pp. 151160.
[9] B.C. Vemuri and A. Radisavljevic, il liresolution stochastic hybrid shape models
with fractal priors," ACM Trans. on Graphics, pp. 177207, 1994.
[10] S. Han, D.B. Goldgof, and K.W. Bowyer, "Using hyperquadrics for shape recovery from
range data," in IEEE Proceedings of the 3rd International Conference on Computer
Vision, Berlin, 1993, pp. 492496.
[11] T. O'Donnel, T. Boult, and A. Gupta, "Global models with parametric offsets as applied
to cardiac," in IEEE Proc. of CVPR, 1996, pp. 293299.
[12] V. Caselles, F. Catte, T. Coll, and F. Dibos, "A geometric model for active contours in
image processing," A,..o, cische Mathematik, vol. 66, pp. 131, 1993.
[13] R. i.illili J.A. Sethian, and B.C. Vemuri, 'lIlpe modeling with front propagation
: A level set approach," IEEE Trans. Pattern Analysis and Machine Intelligence, vol.
17, no. 2, pp. 158175, 1995.
[14] T. '.1I ii,, i,,y and D. Terzopoulos, "Topologically adaptable snakes," in IEEE ICCV,
1995, pp. 840845.
[15] L.D. Cohen, "Auxiliary variables and twostep iterative algorithms in computer vision
problems," Journal of Mathematical Imaging and Vision, vol. 6, pp. 5983, 1996.
[16] A. Gray, Modern differential geometry of curves and surfaces, CRC Press, Boca Raton,
1993.
[17] :. Kass, A. Witkin, and D. Terzopoulos, "Snakes : A active contour models," Int.
Journal of Computer Vision, vol. 1, pp. 321331, 1987.
[18] P.S. i.i,'c, S. i.i, 1 and G. .i, ii1,ii "Bsnakes: implementation and application to
stereo," in IU Workshop, 1990, pp. 720726.
[19] D. Calvetti and L. Reichel, "Application of adi iterative methods to the restoration of
noisy images," SIAM J. Matrix Anal. Appl., vol. 17, no. 1, pp. 165186, 1996.
[20] W.R. Dyksen, "Tensor product generalized adi methods for separable elliptic problems,"
Siam J. i~,.,,. c. Anal., vol. 24, no. 1, pp. 5976, 1987.
[21] R.E. Lynch, J.R. Rich, and D.H. T!,ii_,,1 "Tensor product analysis of alternating
direction implicit methods," J. Soc. Indust. Appl. Math, vol. 13, no. 4, pp. 9951007,
1965.
[22] D.V. Ouellette, "Schur complements and statistics," Linear Algebra and its Applica
tions, vol. 36, pp. 1872'1i., 1981.
[23] W.H. Press, S.S. Teukolsky, W.T., Vetterling, and B.P Flannery, a',,,, rical Recipes in
C, Cambridge University Press, 1992.
[24] A. Lu and E.L. Wachspress, "Solution of lyapunov equations by alternating direction
implicit iteration," Computers Math. Applic., vol. 21, no. 9, pp. 4358, 1991.
[25] S. H. Lai and B. C. Vemuri, "Physicallybased adaptive preconditioning for early vision,"
IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 19, no. 6, pp.
.Q"I ,1 7, 1997.
[26] O.C. Zienkiewicz and R.L. Taylor, T7. Finite Element Method, :'. GrawHill, V: York,
NY, 1989.
