From Global to Local, a Continuum of Shape
Models with Fractal Priors *
B. C. Vemuri1 and A. Radisavljevic2
1Department of Computer & Information Sciences
2Department of Electrical Engineering
University of Florida, Gainesville, FL 32611
Abstract
3D Shape modeling has been a very prominent part of Computer Vision over the
past decade. Several shape modeling techniques have been proposed in literature, some
are local (distributed parameter) while others are global (lumped parameter) in terms
of the parameters required to describe the shape. Hybrid models that combine both
ends of this parameter spectrum have been in vogue only recently. However, they do
not allow a smooth transition between the two extremes of this parameter spectrum.
In this paper, we introduce a new shape modeling scheme that can t,.f, n,,'
smoothly from local to /glhal models or viceversa. The modeling scheme utilizes a
hybrid primitive called the deformable superquadrics constructed in an orthonormal
wavelet bases. These multiresolution bases provide the power to continuously trans
form from local to global shape deformations and thereby allow for a continuum of
shape models from those with local to those with global shape descriptive power
to be created. The characteristic of continuously transforming from local to global
shape deformations allows us to generate fractal surfaces of arbitrary degree that can
be useful in describing natural detail.
We embed these multiresolution shape models in a probabilistic framework and
use them for segmenting anatomical structures in the human brain from MRI data.
A salient feature of our modeling scheme is that it can naturally allow for the incor
poration of prior statistics of a rich variety of shapes. This stems from the fact that,
unlike other modeling schemes, in our modeling, we require relatively few parameters
to describe a large class of shapes completely.
*This work was supported in part by the University of Florida Brain Institute
1 Introduction
Modeling shapes is an important and integral part of Computer Vision as well as Computer
Graphics. In Computer Vision, modeling shapes is important for shape reconstruction and
shape recognition from sensed data, while in Computer Graphics, the application is for shape
synthesis. Many shape modeling schemes have been proposed in the literature. In Computer
Vision, the motivation for modeling shapes has been primarily driven by either shape re
construction or shape recognition tasks. Shape reconstruction from sensed data requires a
broad geometric coverage. The models must recover detailed structure from noisy data using
only the weakest of the possible assumptions about the observed shapes. Generalized spline
models with continuity constraints are well suited for fulfilling the goals of shape reconstruc
tion (see [16, 4]). Generalized splines are the key ingredient of the dynamic shape modeling
paradigm introduced by Terzopoulos et. al., [18]. Incorporating dynamics into shape mod
eling with generalized splines, allows for creating realistic animation for computer graphics
applications and also for tracking moving objects in computer vision applications. With the
advent of the dynamic shape modeling paradigm, there was a flurry of research activity in
this area, with numerous application motivated modifications to the modeling primitives,
and the external forces derived from data [20, 21, 22, 8, 5]. However, the distributed na
ture of these models each point of the spline surface potentially contributes three spatial
degrees of freedom makes them unsuitable for shape recognition tasks. Shape recognition
a highlevel vision task requires that the shape models be characterized by a small set
of parameters. Lumped parameter models are thus better suited for such tasks. Some ex
amples of lumped parameter models include, superquadrics [3, 12, 2] and the parametrically
deformable models [23, 9].
Recently, a hybrid modeling scheme dubbed, deformablee superquadrics" was introduced
by Terzopoulos et al., [19]. These models have the advantage of combining the descriptive
power of lumped and distributed parameter models and thus, simultaneously satisfying the
requirements of shape reconstruction and shape recognition. However, these hybrid models
do not exhibit a smooth transition in the number of parameters required for their description.
For example, the two extremes of the spectrum are occupied by lumped models at one end
and distributed models at the other. The former requiring, few parameters and the later
requiring a large number of parameters for their description. The modeling scheme described
in [19] involves a superposition of the lumped and distributed parameters to completely
characterize the plethora of generated shapes. Although, the hybrid models can be used to
generate both local and global shape deformations, they do not depict a smooth transition
between these types of deformations. In addition, there is no smooth transition in the number
of parameters that can capture descriptions of the range of generated shapes.
In this paper we introduce a new multiresolution dw.,ii;'. modeling scheme that uses
an orthonormal wavelet basis [11, 7]. Our models have the unique property of being able
to ,ii..,,'/ lilscale up or down the range (from local to global) of possible deformations and
the number of parameters required to characterize them. This is achieved by adding to the
global (lumped) parameters of the deformable superquadrics, a set of wavelet coefficients
of the membrane spline at any resolution. Earlier work reported in Pentland [13] on the
use of wavelets in the context of shape modeling was limited to surface reconstruction. He
used wavelets basis to approximate the Eigen vectors of the stiffness matrix assembled in the
finite element solution to the surface reconstruction problem using the thinplatemembrane
splines. The main idea of his work was to speed up the solution to the surface reconstruction
problem by a change of basis, from regular nodal basis to the natural basis represented by
the Eigen vectors of the stiffness matrix. Wavelets basis were used as an approximation
to the natural basis i.e., the Eigen vectors of the stiffness matrix. In the wavelet basis,
the stiffness matrix is approximately diagonal and this approximate diagonalization leads to
savings in storage space as well as computation time for the solution. In our research, the
key focus is on development of a general shape modeling scheme which when embedded in
a probabilistic framework can be utilized for easily incorporating a larger class of shapes as
prior information. In addition, we also emphasize on the efficient computation of the model
fitting process to the data.
The thinplatemembrane regularizers were also used for a very different purpose by Szeliski
et. al., [15] namely, generation of constrained fractals which are useful in representing natu
ral images. They derive the fractal characteristic of a membrane and thinplate regularizers
and show that the membrane splines give raise to fractals of degree two while the thinplate
splines yield degree four fractals. They approximate a fractal of inbetween degree by using
a blend of the membrane and thinplate models. The modeling scheme we introduce in this
paper allows us to generate fractals of inbetween degree without having to resort to approx
imations with a blend of the plate membrane models. We achieve this by modulating the
frequency characteristics of the wavelet basis with a polynomial of appropriate degree. In
order to generate stochastic fractal surfaces, we convert the energies into probability distri
butions using ideas from statistical mechanics [10] and employ a stochastic multiresolution
algorithm. We will present some examples of the synthesized fractal surfaces in a subsequent
section.
The rest of the paper is organized as follows: In section 2 we briefly describe the geometry
of the modeling scheme and the multiresolution basis for expressing the geometry, followed
by a new way to augment the global shape of the deformable superquadrics primitive for
achieving the smooth transition from global to local shape descriptions. Section 3 contains
discussion on how to embed our models in an probabilistic framework and the derivation of
the mechanics of the modeling scheme. In section 4 we derive the internal energy of the prior
model in the wavelet basis and present the fractal characteristics of this prior model with
the aid of synthesized fractal surfaces. Section 5 contains a description of, the unsupervised
learning technique that we use for incorporating specific a priori knowledge into our modeling
framework and, implementation results on real MRI data followed by conclusions in section
6.
2 Multiresolution Geometry of the Modeling Scheme
In this section, we will briefly describe the geometry and the construction of a multiresolution
wavelet basis for the modeling scheme. To dramatically increase the library of shapes that
can be modeled and to achieve a smooth transition from global to local shape descriptions,
we will present a new method that will augment the global parameters of the deformable
superquadrics.
2.1 Geometry of the Deformable Superquadrics
The deformable superquadrics are closed surfaces in space whose intrinsic coordinates u =
(u, v) are defined on a domain Q. The positions of points on the model in an inertial frame
of reference are given by a vectorvalued function x(u) = (xI(u), x2(u), x3(u))T. In a model
centered coordinate frame, the position vector x becomes
x = c + Rp (1)
where c is the location of the center of the model, and the rotation matrix R specifies
the orientation of the modelcentered coordinate frame. Hence, p(u) denotes the canonical
positions of points on the model relative to the model frame. We further express p as the sum
of a reference shape s(u) and a displacement d(u) namely p = s + d. For the parameterized
reference shape s(u) we use the superquadrics with the bending deformation (as in [2]), due
to it's attractive global shape characterization. The reference shape s is given by,
s = B(s) (2)
Where B is a bending deformation, a function of two parameters: the radius of curvature
k and angle of the bending plane a (see [2] for details). The parametric equation of a
superquadric s is given by
/aCJC1 S 2
S= a2CU6SU2 (3)
a3 Se1
Where, r/2 < u < v/2 7 < v < w, and S = sgn(sin w) sin wl and C' = sgn(cos w)cos wl .
Here, 0 < al, a2, a3 < 1 are aspect ratio parameters, (1, 62 > 0 are "squai. i.  parameters.
This defines the geometry of our underlying reference shape namely, the superquadric. We
collect all these shape parameters along with the two bending parameters into a vector q,.
We will denote c by the vector qc and use qo to denote the vector of rotational coordinates
of the model in a quaternion representation (see [19]).
An important aspect of our reference shape which differs from previous work of Terzopou
los et al., [19] and Bajscy et al., [2], is that we have constructed a uniform tessalation of
the superquadric reference shape. Superquadric shapes exhibit anomalous behavior when
the squareness parameters 1, 62 depart from unity, causing grid lines to become unevenly
spaced leading to distortion of the deformational properties of the superimposed membrane
surface and inefficient use of computational resources. In addition, numerical instabilities
may be caused in the computation of derivatives due to the large grid spacing. A Uniform
tessalation of the superquadric can be used to overcome these problems. To achieve the
uniform tessalation, we first construct uniformly spaced grid lines in the parameter direction
u and apply the same technique to achieve uniform spacing along the parameter direction
v. Thus, we seek a function of the parameter u (same argument applies to the parameter v)
which when substituted into the equation of the superquadric i.e., equation 3, would yield
a constant differential arc length increment for a constant differential increment in u. For a
horizontal elliptical contour,
ds= (x/)2+ (yIdu (4)
Setting ds/du = constant = (circumference/27) and using (x = alcos'(f(u)), y =
a2sin(f(u)), z = 0) results in an ordinary differential equation for f(u)
df I tan(f) (
= K (5)
du e/a# l tan(f) 41 cos(f) 2 +a sin(f) 12
f(O) = 0 as the initial condition (6)
Where K is a constant. This ordinary differential equation can be solved numerically. In
figure 1, a superquadric with its natural nonuniform tessalation is shown side by side to the
same superquadric ,.., l,,, tessellated using the method described above.
In the following sections, we show that the displacement vector function d = )Sqd is
completely determined by a set of wavelet coefficients collected into the vector qd, with
S corresponding to the discrete wavelet transform and ) containing nodal interpolation
functions. Thus, the state of the model is fully expressed by the vector q =q q q q T q .
2.2 Multiresolution Analysis and Wavelet Bases
Multiresolution analysis has been quite popular in computer vision research for the past
decade. Most of the research in the past has dealt with speeding up computations in solu
Figure 1: Tessellations of a Superquadric, (a) nonuniform (b) uniform
tions of various lowlevel vision problems using multiresolution techniques [17]. With the
advent of the wavelet transform (WT), multiresolution analysis of functions has become
more attractive due to the fact that wavelet decomposition of a function (signal) into dif
ferent frequency components allows for studying each component at a resolution matched
to its scale. In the following, we will very briefly describe the construction of a wavelet
based multiresolution approximation in 1D and its extension to a 2D function specifically
the membrane displacement function d(u, v).
2.2.1 Detail Function from Multiresolution Approximations
A multiresolution approximation of a scalar valued function d(u) e L2 (R) at a resolution j
is obtained by projecting d(u) into a vector subspace Vj C L2(1) where integer (k) translates
Pjk(u) = 2/j (2 u k) of a dilated scaling function 0(u) form the orthonormal basis of
Vj. The approximation dj(u) is expressed as :
d((u)= ^JlkJU}k
k
Where the set of coefficients ack = (d(u), Ojk) for all integer k corresponds to the discrete
characterization of this approximation.
An important property of multiresolution approximations is that all the necessary infor
mation to compute the function at a smaller resolution j + 1 is contained in the function at
Figure 2: Division of the frequency domain for a) multiresolution and b) wavelet represen
tation.
resolution j. In other words the subspaces Vj satisfy the containment hierarchy,
.....V3 C V2 C V1 C VO (8)
The scaling function 0 has the effect of a lowpass filter with the cutoff frequency decreasing
proportionally to 2j (see figure 2(a)). The key idea of the wavelet basis is to span the
difference spaces Wj+1 between the two approximations Vj and Vj+1 such that Wj+1 I Vj+i.
A wavelet decomposition therefore consists of splitting the function at resolution j into the
detail function in Wj+1 and the next lower (j+1) resolution approximation in Vj+i. This
leads directly to a scale wise, orthogonal decomposition of L2(R). To represent the function
d(u) of a finite resolution we can assume d(u) e Vo, where Vo is given by,
..... W3 ( w2 ( w1 = (9)
with a denoting the direct sum operator. The hierarchicalorthonormal basis for d(u) is given
by a family of dilations and translations of a wavelet function 0(u) : Ojk(u) = j(2ju 
k). Thus, d(u) can be be written in this basis as,
d(u) = ;' "(u) (10)
j,cc cc
Where the discrete wavelet coefficients are computed as Sjk = Kd(u), Yjk).
As in the case of the multiresolution approximation, wherein a low pass filtering interpre
tation was given to the scaling function, a filtering interpretation can also be given to the
WT. At each scale j the WT corresponds to a logarithmically spaced bandpass filter (see
Figure 2(b)). Our work in this paper is primarily based on this filter bank interpretation of
the WT as discussed in Rioul et. al. [14].
2.2.2 Extension to 2D and Fast Implementation
Extension to 2D of the multiresolution approximation using wavelet basis can be easily
achieved using separable wavelets obtained from the products of onedimensional wavelets
and scaling functions [11]. For the multiresolution approximation of the membrane d(u, v)
we use the tensor products of p and y given by ( see Mallat [11]),
'"(UV) = 9(U)jk0 (v)jl
'(UV) = p(U)j (V)j, (11)
Note that, since the displacement d is a vector valued function, we have to repeat the tensor
product approximations shown above for each component of the displacement vector. The
superscript H (L) represents high (low) pass filtering in each of the directions u and v in
the parameter domain Q, and the indices k and 1 represent the integer shifts (as in the 1D
case) in u and v directions respectively. We collect the coefficients of the inner product of
the basis and the function d(u, v) i.e., 6 6&{H and 6HH into a vector qd as illustrated in
Figure 3. This displacement vector qd contains a complete hierarchical description of the
membrane surface. Note that qd has the same number of elements as the original nodal
discretization at j = 0. For our WT implementation we adopted the specific 0(u) and y(u)
functions as given in Mallat [11], which are based on the cubic spline interpolants. If we
denote the original nodal discretization at resolution j = 0 by a vector 0o then, we have,
ao =S qd Reconstruction
qd = ST ao Decomposition (12)
d(u,v) = D(u, v) ao Interpolation
where 1 represents a vector of interpolation functions (as in the 1D case, see equation 7)
and multiplication by the orthogonal matrix S corresponds to performing the WT. A fast
Figure 3: a) Organization of the displacement vector qd, b) separation of the "global" and
"local" parts of qd
implementation of the WT is achieved using Quadrature Mirror Filters (QMFs). A self
explanatory block diagram of the filterbank interpretation of the QMF is presented in figure
4 (see Mallat [11] for details). Figure 4 shows an efficient way to implement the WT, with
a computational cost of order O(N) (where N is the number of grid nodes).
Finally, vector qd is expressed in an orthonormal wavelet basis. Along with other pa
rameters of the model it is subjected to the iterative minimization process as discussed
subsequently.
2.3 From Global to Local Shapes with Wavelet Coefficients
In figure 3, we explicitly depict the hierarchical decomposition of the displacement vector qd
in an orthonormal wavelet basis. The coefficients 6j toward the upper left corner are fewer
then those toward the bottomright corner, i.e. the corresponding grid nodes are spaced
farther apart and the basis functions are coarser. Few parameters needed in the description
and ability to depict the global shape gives the low frequency components of the wavelet
representation a "semiglobal" character. Depending on the application we can draw the
boundary at an arbitrary jo to separate the global and local part of vector qd,
qd = [qda qd]
qdg = [ ,"....,6J0 1 (13)
qdl = [j0oi, ***, 1]
Figure 4: a) Decomposition and b) reconstruction of a 2dimensional array adopted from
Mallat[89]
where qd contains the top (coarser) portion of the multiresolution pyramid. One of the
useful aspects of including qd into the global description along with superquadric parameters
is the greater variety of shapes we can generate, such as coarse dips, bumps etc. (see figure
5). For a comparison of possible global deformations generated via tweaking of superquadric
parameters and the coarse level wavelet coefficients of the displacement function, we illustrate
the deformed shapes adjacent to one another. Figure 5 depicts the original reference shape
on the top left hand corner, global deformation resulting from tweaking a superquadric
parameter is shown on the top right hand corner. The bottom two shapes are generated via
global deformations caused by tweaking the coarse level wavelet coefficients.
When we build our prior model with mean values and covariances on q, and q, the model
assumes this enhanced shape as it's new rest state and it deforms away from this shape under
the influence of sensed data. The need for incorporating more specific information such as
this, into the prior model is especially present in difficult segmentation problems encountered
in Medical Imaging applications eg. segmentation of anatomical shapes from MR images.
We will present two such examples in this paper.
Figure 5: a) Global shape deformations generated via tweaking global parameters of the
superquadric model, b) global deformation generated by tweaking of coarse scale wavelet
coefficients
3 Probabilistic Framework and Model Mechanics
In the last section, we have described our multiresolution deformable superquadric model
in a wavelet basis. The multiresolution wavelet basis provides a computationally efficient
solution to the problem of fitting the deformable superquadric model to data sets and also
provides us a way to enhance the library of shapes that can be described using a few global
parameters. In this section, we will present a way to incorporate specific prior knowledge
about the shape to be extracted from the data and then derive the mechanics of our modeling
scheme.
We cast the model fitting process in a probabilistic framework and incorporate prior dis
tributions on the vector of geometric parameters that we are estimating. The combination
of prior and sensor models results in a Bayesian model, since the posterior distribution of the
parameters we are trying to estimate conditioned on the data can be computed using Baye's
rule. In this paper, we use the multiresolution physicallybased shape models described in
earlier sections as the prior models. We collect the statistics on the global parameters of the
model via a training process described in a subsequent section and for the local parameters,
we convert the continuous strain energy that governs the deformation of the deformable su
perquadric model away from its natural shape into a probability distribution over expected
shapes, with lower energy shapes being the more likely. This is achieved via a technique of
statistical mechanics conversion of energies into probabilities using Boltzmann, or Gibbs,
distribution [10]. The link between physicallybased deformablee models) models and priors
is conveniently established using a Gibbs distribution of the form
1
p(q)= ep(E,(q)), (14)
Zp
where E,(q) is the discretized version of the internal smoothness energy of the model, and
Z, (called the partition function) is a normalizing constant. What was originally an elastic
energy restoring a model towards a rest state now becomes a probability distribution over
expected shapes, with lower energy shapes being more likely. To complete the formulation
of the estimation problem, we combine this prior model with a simple sensor model based
on linear measurements with Gaussian noise
1
p(D/q) = ep(ED(D, q)), (15)
ZD
where ED(D, q) can be either an edgebased potential energy synthesized from the MR
images or the potential energy in the springs attached from 3D data points to the surface
of the deformable superquadric constraining it to conform to the observed data. Combining
the prior and sensor models (from equations 14 and 15 respectively) using Baye's rule, we
obtain the posterior distribution
p(D/q)p(q) 1
p(q/D) ep(D= (E(q)), (16)
p(D) Z
where
E(q) = E(q) + ED(D, q). (17)
Computing the Maximum A Posteriori (MAP) estimate, i.e., the value of q that maximizes
the conditional probability p(q/D), provides the same result as finding the minimum energy
configuration of the physicallybased model. However, the advantage of the probabilistic
framework being that, it allows for explicitly incorporating prior model and sensor model
characteristics in addition to providing an uncertainty in the estimates.
The internal strain energy Ep of the model in our case is given by a quadratic form,
Ep = 1/2[(q q) TK(q q)]. Where, q is the rest state of the model and the matrix K]
corresponds to the stiffness of deformation. We will discuss the internal energy Ep further in
a subsequent section. The two types of sensor model energies for ED(D, q) discussed earlier
are expressed as
f i /3(D x(q, u)2 spring energy
ED(D, q)= /3 P(x(q, u)) image potential (
In the above equations ui determines the closest point x to the given measured point Di.
We adjust the locations of these grid coordinates as the shape of the model evolves (see [19],
for more on migrating point of influence). The magnitude of parameter 3 is related to the
uncertainty of sensor measurements. P(x(q,u)) denotes an edge based potential derived
from the MR image data. Maximization of of the posteriori amounts to minimizing the
total energy E = Ep + ED for which we employ a gradient descent method that amounts to
iteratively solving,
qk+l = qk AE(q) (19)
Oq
For the derivative of the energy on the right hand side in the above equation, we use the
chain rule to get a= Oa ax where the contribution of each parameter group is
ax ac a(Rp) Ras RO(4ISqd)
q 9qc,' sq0 9q,' 09qd
= [I B RJ iS] (21)
= L (22)
Here, we adopt the notation as in [19]. Partial derivatives of the prior model and data
energies defined above can be now written as
=E)
OE = K (q tj) (23)
aq
and
OED Eq /(Di x(q, u1)L spring energy
Oq q = V P(x(q, u)) L image potential (24)
By substituting equations 23 and 24 into the equation for gradient descent 19 we obtain
an expression for iterative updating of the state vector q
(q(k+) (k) +A (fk) K(q(k) ) (25)
where A is a diagonal matrix containing step sizes and consequently effecting the speed
of convergence. Close similarity of this evolution equation to the one derived in [19] reveals
that equation 25 can be modified slightly and effectively used to generate motion sequences
corresponding to first order dynamics Cql + K(q q) = fq. In addition to producing the
statistically most likely state q, this framework provides the uncertainty in the estimate
via the associated covariance matrices. In our work, we employ the above gradient descent
procedure in a multiresolution framework for computational efficiency as well as alleviation
from entrapment in local minima.
4 Prior Model in Wavelet Basis
In the previous section we described a probabilistic framework for embedding our model. We
now derive the explicit form of the internal energy of the model in the wavelet basis. This
energy is then converted to a Gibbsian probability distribution using the technique discussed
in section 3.
4.1 Internal Energy of the Prior Model
In addition to generic smoothness assumptions, a prior model should incorporate the infor
mation about the possible shapes that an observed object can assume. This can be very
important in applications involving segmentation of known natural shapes from "unfriendly"
(low contrast and noisy) images. We express the probability of a given configuration q
through a Gibbs distribution (see section 3), which assigns a low probability to configura
tions with large energy Ei and vice versa. The configuration with Ei = 0 represents the
rest state of the model. The internal energy of the membrane spline (which augments the
reference superquadric shape s of the model), in the continuous form is given by,
Ea = / wod2 + w+ a) du (26)
With no additional prior information about the shapes being modeled, the underlying
global shape of the model namely, the superquadric, can be interpreted as the mean value
or the rest state of the membrane surface expressed by the displacement function d. If we
rewrite the deformation energy in terms of vector p = s + d we have an equivalent expression
to equation 26,
/((O 8(p s) 2 + p S)
Ea(p S. +) du (27)
Where the role of the global (reference) shape s as the rest state becomes apparent. We
extend our underlying global shape s with a low resolution mean displacement function
d = dt0 corresponding to parameters qdg as defined in equation 13 and shown in figure 3.
This is done by altering the deformation energy equation 26 to,
E = / wo(d d)2 + w (d d) + ((d d) ) du (28)
.a1 u av
The discrete version of this internal strain energy can be written as,
1
ET = (qd qd)T Kd (qd qd) (29)
2
Where, Kd is the stiffness matrix that determines the elastic properties of the model and
qT = [qdg, 0,...,0] incorporates the statistics of the possible shapes that can be encoun
tered in the segmentation problem. These statistics are gathered through training process
described in section 5. It is evident that the number of parameters in qd, is low, i.e. the
representation of the model remains compact and thus making the training process feasible.
We include other global shape parameters of the model in our prior statistics by extending
the above energy as follows,
E = q(q q)T C (q, q) +
(qb q)T b b 1 (qb q) +
2(qd qd)T Kd (qd qd) (30)
Where s denotes the superquadric shape without the bending while, Cs represent the co
variance matrices which are accumulated in the iterative training phase. In section 4.2 we
show the derivation of stiffness matrix K in the wavelet basis and comment on its structure.
4.2 The Stiffness Matrix in Wavelet Basis
The stiffness matrix for the membrane operator is a very large, sparse and banded matrix in
the ordinary nodal basis. By expressing this in an orthonormal wavelet basis, we can achieve
near diagonalization of this matrix. We will show how the filter bank interpretation of WT
(see figure 2) relates directly to the diagonal entries of the stiffness matrix K. This approach
in turn gives us a straightforward and inexpensive way to compute K.
The standard way to obtain stiffness matrix in wavelet basis is essentially a similarity
transformation. By combining equations 12, 28 and 29 we get,
KO = woST (f 4 du) S (31)
K1 = wiS Q/< l) du) S (32)
K = K + K1 (33)
Here Ko controls the local magnitude and K1 the local variation of the (membrane) displace
ment function. Due to orthogonality of the interpolation functions in ) and orthogonality
of S, equation 31 simplifies to K = wol. However, the part K1 of the stiffness matrix in
equation 33 needs to be computed by a similarity transformation K1 = STjI S, where K1
represents the membrane stiffness matrix in nodal basis (see [13]). The drawback of comput
ing K1 using this similarity transformation is the large computational cost incurred (O(N2),
where N is the number of nodal variables) which can be avoided by using the property that
differential operators in wavelet basis remain sparse (see [1]) and have strong diagonal dom
inance (see [13]). This allows us to compute the diagonal elements of the stiffness matrix
directly by taking inner products of the corresponding finite element basis functions (an
O(N) computational cost),
K3k = W1 2 Tk() 2d (34)
oo
Figure 6: (a) Scaling of the diagonal entries in the stiffness matrix K of the membrane
model with spectral density, a2 and an arbitrary density w~, (b) Increasing energy with
increasing resolution
which is done in frequency domain by employing Parseval's theorem [6]. Extension to 2
D case is straightforward and achieved by using the tensor product of the basis functions.
The factor jw results from the Fourier transform of the first derivatives representing the
membrane regularizer. Note that both K and K1 have diagonal structure and that only the
relative magnitude of their diagonal entries determine the smoothing effects of each matrix.
The filter bank interpretation of WT (Figure 2) explains how coefficients corresponding
to lower resolutions j contribute less energy than those corresponding to higher resolutions.
Figure 6 illustrates two aspects of the modeling, (a) the membrane displacement function
which has a spectral density characterized by w2 and is a fractal process of degree two,
can be used to modulate the ideal bandpass characteristics of the wavelet basis for creating
the appropriate stiffness matrix. In addition, one can synthesize fractals of arbitrary degree
by using a polynomial function w~ to modulate the bandpass characteristic of the wavelet
basis for creation of the appropriate stiffness matrix; (b) The concept of increasing energy
contribution with increasing resolution via a 1D example. For instance, if we assume an
ideal bandpass characteristic for Tj(w) between wj and 2wj we would have the following
relationship between elements of K at j and (j + 1) resolutions.
Kik = w Ijw l2d = 8I1(j+i)k (35)
W]A
Figure 7: Diagonal entries of K as an elevation map, for a (32X32) grid computed via a
similarity transformation
Therefore, for the membrane regularizer, we notice that the energy at a resolution (j + 1) is
eight times that at the next lower resolution j.
For a verification of the above discussed sclaing behaviour, we computed the diagonal
entries of the K matrix for a (32 X 32) grid of nodes using a similarity transformation
technique discussed earlier. Figure 7 depicts the diagonal entries as an elevation map,
where the elevation corresponds to the magnitute of the entries. The entries are organized
as in a standard representation of images in wavelet basis (see figure 3 for the standard
organization).
4.3 Prior Model of General Fractal Degree
Splines and stocahstic fractals are complementary techniques for modeling freeform shapes.
Splines are well suited for modeling smooth man made objects and are easily controlled.
Fractals on the other hand are well suited for modeling irregular shapes and are difficult to
constrain. Szeliski et. al. [15] describe a way to constrain fractals by establishing a formal
connection between splines and fractals using Fourier analysis of energy minimizing splines.
Our prior model corresponding to the membrane spline is a Markov random field with a
Gibbsian distribution. Its spectrum can be determined using Fourier analysis techniques. In
Figure 8: Surface fitting result a) MAP Estimate and b) Representative sample from the a
posteriori distribution, corresponding to a fractal surface (membrane model)
their work on constrained fractals Szeliski et. al. [15] show how prior models corresponding
to a membrane and thinplate regularizers posses a spectrum w2 and 4 respectively, which
is fractal i.e., it is selfaffine over scale. These spectra characterize special cases of fractals
of general form given by,
1
S(w) oc (36)
In [15] fractals of an inbetween degree 2 < 3 < 4 were approximated by blending the
membrane and thinplate terms. In contrast, our multiresolution wavelet based modeling
scheme discussed earlier allows us to construct fractals of any degree directly by utilizing
the natural scale space characteristics of wavelet basis. Elements of K governing the prior
distribution are computed using either equation 34 or 35 with w2 replaced by 2A. Thus,
yielding a stiffness matrix corresponding to regularizer whose spectral distribution is of order
3. A random sample from the corresponding Gibbsian distribution of this prior model will
be a fractal surface of degree 3.
Figure 8 shows a surface fitting result to 3D point data set with (a) MAP estimate and (b)
representative sample from the Gibbs distribution p(qD) reflecting the underlying fractal
prior. Fractal surfaces are very useful for modeling natural detail. In our application of
Medical Imaging, our goal is to realistically model the shape of the anatomical structures
Figure 9: Supervised learning to incorporate information into the prior model
found in the human body, specifically the brain. Experiments toward achieving this goal are
currently underway.
5 Supervised Learning and Implementation Results
In our work, prior knowledge about the object shapes that are to be segmented from the
MRI data is collected via a supervised learning technique. A training sequence of several
iterations provides samples of correctly segmented shapes for specific classes of objects in
our experiments a i/i, i, and the hippocampus from the human brain. Parameter sets q(i)
corresponding to these shapes represent the training sequence for the global part (q, + qdq)
of our prior model (see figure 9). With each iteration i we improve estimates of the mean
vector q and covariance matrix C through a simple equation for accumulating measurements
as indicated in the Figure 9. Since we update only the statistics of the global parameters
of the model, the local deformations remains constrained only by smoothness assumptions
captured in the stiffness matrix K (refer to equation 30). Therefore, most of the flexibility of
the locally deformable surface remains unchanged. We recommend human visual supervision
of the surface fitting in the training stage to assure that only good segmentation results enter
the statistics of the prior model. An extension to multiclass case is straightforward, and
can be achieved by defining multiple prior distributions pj(q) for each class j.
This method clearly demonstrates the useful aspect of global models which offer a compact
description of the overall features of the object. Deformable superquadrics with the global
shape extended by some offset deformation (as presented in previous sections) allow a much
wider range of shapes to be representable in the prior model, all at the expense of few
additional parameters. This makes our model very attractive for shape recognition as well
as shape reconstruction tasks.
We now present two examples depicting segmentation of a gyrus and a hippocampus from
MRI data. The data sets consisted of 64 slices each, with each slice of 1.25mm thickness. We
applied a simple 3D edge detection technique to create a gradient force field which was used
to attract the inialized model toward the edges of the gyrus and the hippocampus in the
first and second examples respectively. Figure 10, illustrates our results displayed with slices
of the reconstructed shape superimposed on the corresponding slices from the original data.
The region of interest has been boxed in the original image and our results are displayed
only for this region of interest. Left to right, the images are, a slice from the original data,
a slice from the initialized 3D model, the corresponding slice from the final segmentation.
Further experimentation with our modeling scheme is currently under way and we expect
to better our results by using more sophisticated 3D edge detectors [5] for defining the
external forces to draw the model into the desired shape.
6 Conclusion
In this paper, we introduced a new shape modeling scheme that can transform ,i..,... I,/
from local to global models or viceversa. The modeling scheme utilizes the hybrid primitive
called the deformable superquadric constructed in an orthonormal wavelet bases. These bases
provide the power to continuously transform from local to global shape deformations and
Figure 10: Segmentation results for (a) a gyrus and (b) a hippocampus from MRI data;
Figure shows slices from 3D, of data, of the initialization, and the reconstruction respectively.
thereby allow for a continuum of shape models from those with local to those with global
shape descriptive power to be created.
When embedded in a probabilistic framework, our modeling scheme allows us to incor
porate more detailed prior information than when using lumped parameter models, while
keeping the number of parameters required to describe the shape low. The same level of de
tail may also be incorporated by using some existing modeling techniques [9, 19] but, at the
expense of increasing the number of parameters tremendously. From a computational point
of view, our modeling scheme is very attractive for segmentation purposes as multiresolution
relaxation is a natural by product of implementing the gradient relaxation method in the
wavelet basis. Also, entrapment in local minima during the minimization of the nonlinear
energy function can be alleviated. Our future efforts will be directed toward testing with
more MRI data and using better 3D edge detectors to define external forces.
References
[1] B. K. Alpert, [1992], "Wavelets and other bases for fast numerical linear algebra," in
Wavelets A Tutorial in Ti .. ry and Applications, C. K. Chui (ed.), pp. 181216.
[2] R. Bajscy and R. Solina, [1987], "Threedimensional object representation revisited,"
in IEEE First Conference on Computer Vision, London, England, pp. 231240.
[3] A. H. Barr, "Superquadrics and anglepreserving transformations," IEEE Comp. Graph
ics and Applications, Jan. 1981, pp. 1123.
[4] A. Blake and A. Zisserman [1987], Visual Reconstruction, MIT Press, Cambridge, MA.
[5] L. D. Cohen and I. Cohen, [1992], "Deformable models for 3D medical images using finite
elements and balloons," IEEE Conference on Computer Vision and Pattern Recognition,
June, Urbana, Champagne, Ill.
[6] R. N. Bracewell, [1978], The Fourier Transform and its Applications, McGrawHill Book
Company, NY.
[7] I. Daubechies [l'"], "Orthonormal bases of compactly supported wavelets," Commu
nications of Pure and Applied Mathematics, XLI: pp. 909996.
[8] H. Delingette, M. Hebert, and K. Ikeuchi, [1991], "Shape represnetation and image
segmentation using deformable surfaces," IEEE Conference on Computer Vision and
Pattern Recognition, Maul, Hawaii, pp. 467472.
[9] L. H. Staib and J. S. Duncan, [1991], "Boundary finding with parametrically deformable
contour models," IEEE Transactions on Pattern A,,,,/..', and Machine Intelligence, Vol.
14, No. 11, pp. 10611075.
[10] S. Geman, and D. Geman, [1984], "Stochastic relaxation, Gibbs distribution, and
Bayesian restoration of images," IEEE Transactions on Pattern A,,,1/..', and Machine
Intelligence, Vol. 6(6), pp. 721724.
[11] S. G. Mallat, [1989], "A theory of multiresolution signal decomposition: the wavelet
representation," IEEE Transcations on Pattern A,,,,il.' and Machine Intelligence, Vol.
11(7), pp. 6741,', ;
[12] A. P. Pentland, [1','i.], "Parts: Structured descriptions of shape," in Proceedings of
AAAI, Philadelphia, PA, pp. 695701.
[13] A. P. Pentland, [1990], "Fast surface estimation using wavelet basis," M.I.T. Media lab
Vision and Modeling Group, Technical Report No. 142, M.I.T., Cambridge, MA.
[14] O. Rioul and M. Vetterli, [1991], "Wavelets and signal processing," IEEE SP magazine,
October, pp. 1438.
[15] R. Szeliski and D. Terzopoulos, [1989], "Constrained fractals," in Proceedings of ACM
SIGGRAPH, August, pp. 5160.
[16] D. Terzopoulos [1' '1.], "Regularization of inverse visual problems involving discontinu
ities," IEEE Trans. on Pattern A,.,,lii.' and Machine Intelligence, PAMI8, 413424.
[17] D. Terzopoulos, [1','1.] "Image analysis using multigrid relaxation methods," IEEE
Trans. on Pattern A,.,,lii..' and Machine Intelligence, PAMI8, No. 2, pp. 129139.
[18] D. Terzopoulos, A. Witkin and M. Kass, [1'""], "Constraints on deformable models:
Recovering 3D shape and nonrigid motion," Artificial Intelligence, 36, pp. 91123.
[19] D. Terzopoulos and D. Metaxas, [1991], "Dynamic 3D models with local and global
deformations: Deformable superquadrics," IEEE Transactions on Pattern A,,l.' and
machine Intelligence, Vol. 13(7), pp. 703714.
[20] B. C. Vemuri and R. Malladi [1991], "Deformable models: Canonical parameters for
invariant surface representation and multipleview integration," in Proc. of the IEEE
Conference on Computer Vision and Pattern Recognition, Maul, Hawaii, June 36.
[21] B. C. Vemuri and R. Malladi, "Constructing intrinsic parameters with active models for
invariant surface reconstruction," IEEE Transactions on Pattern A,,,1i1.', and Machine
Intelligence, to appear.
[22] Y. F. Wang, J. F. Wang, [1990], "Surface reconstruction using deformable models with
interior and boundary constraints," in Proceedings of IEEE Conference on Computer
Vision, Osaka, Japan, pp. 300303.
[23] A. L. Yuille, D. Cohen, and P. W. Hallinan, [1989], 1'., lire extraction from faces
using deformable templates," IEEE Conference on Computer Vision, San Diego, CA.,
pp. 104109.
