Multiresolution stochastic hybrid shape models with fractal priors.
3D shape modeling has received enormous attention in computer graphics and 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. We introduce a new shape-modeling scheme that can transform smoothly from local to global models or vice versa. The modeling scheme utilizes a hybrid primitive called the deformable superquadric constructed in an orthonormal wavelet basis. The multiresolution wavelet basis provides the power to continuously transform from local to global shape deformations and thereby allow for a continuum of shape modelsfrom those with local to those with global shape descriptive powerto be created. The multiresolution wavelet basis allows us to generate fractal surfaces of arbitrary order that can be useful in describing natural detail. We embed these multiresolution shape models in a probabilistic framework and use them for recovery of anatomical structures in the human brain from MRI data. A salient feature of our modeling scheme is that it can naturally allow for the incorporation 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.
Introduction
Modeling shapes is an important and integral part of computer graphics as well as computer vi-
sion. In computer graphics, modeling shapes is an important ingredient of shape synthesis while
in computer vision, it is needed for shape reconstruction and shape recognition from sensed data.
Many shape modeling schemes have been proposed in the graphics [31, 4, 28, 26, 33, 13, 37, 21]
as well as vision literature [6, 43, 29, 2, 41, 9, 44, 45, 42, 14]. In this paper, we will be concerned
with shape modeling with a view toward facilitation of shape recovery from 3D data for computer
vision/medical imaging applications. In computer vision, the motivation for modeling shapes has
been primarily driven by either shape reconstruction 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 reconstruction [38, 43, 8]. Generalized splines are the key ingredient of the dynamic shape
modeling paradigm introduced by Terzopoulos et. al., [40, 41]. Incorporating dynamics into shape
modeling 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; several methods were also developed to improve the numerical stability of the
algorithms by using finite element discretization [44, 45, 50, 17, 13, 14]. However, the distributed
nature 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 high-level
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 examples of lumped parameter models
include, generalized cylinders [7], superquadrics [4, 29, 2] and the parametrically deformable
models [52, 34].
Recently, a hybrid modeling scheme dubbed, "deformable superquadrics" was introduced by
et al., [42]. 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. Also, the lumped parameter models allow only for global deformations while
the distributed parameter models allow for local deformations. The modeling scheme described in
[42] involves a superposition of the lumped and distributed parameters to completely characterize
the plethora of shapes occupying the spectrum. Although, these 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 multi-resolution hybrid modeling scheme [47] that uses an
orthonormal wavelet basis [23, 16]. Our hybrid models have the unique property of being able to
smoothly scale up or down the range (from local to global) of possible deformations and the number
of parameters required to characterize them. These properties are inherited by virtue of the
multiresolution wavelet representation. Depending on the application, one may choose a set of
wavelet coefficients at a particular coarse-level from the multi-resolution representation of the hybrid
model and augment the global (lumped) parameters of the model. This augmented set of
"global" (lumped) parameters define a larger class of shapes than those characterized by purely
lumped parameters. We will use this enhanced class of shapes characterized by a few global parameters
as prior information in our probabilistic shape recovery algorithm. In section 5, we describe a
supervised learning method that may be used in selecting an appropriate level in the multiresolution
representation for choosing the set of wavelet coefficients in the aforementioned augmentation
process. Tweaking these coarse-level wavelet coefficients will produce global deformations that are
very different from those produced by tweaking the global (lumped) parameters of the model (see
the bottom four shapes in figure 6). The number of nodes per level in the multiresolution algorithm
used in our work is a power of 2 but the number of coefficients chosen in the augmentation process
need not be a power of two. Thus, giving the power to smoothly/"continuously" scale up or down
the range of possible deformations.
A question that arises in the context of constructing hybrid models in a multiresolution wavelet
basis is, can we achieve the same results via standard finite element subdivision techniques? The
answer is no. Firstly, single level subdivisions of the elements will not achieve the same results
as a multiresolution wavelet representation since, there is more to a wavelet transform than sub-sampling
or up-sampling operations (used to achieve projection from one subspace to another in
the multiresolution spaces) specifically, the high-pass and low-pass filtering operations (see section
2.2). Also, the single level finite element subdivision can not lead to a smooth transition from
local to global deformations as there is no scope for descending/ascending through the space of
deformations. This however may be achieved via the use of a hierarchical basis finite element
scheme. What then are the advantages of a wavelet basis in this context? There are several
advantages and some of them are, (a) the number of coefficients used in a wavelet representation
are the same as that in the original nodal representation of the function i.e., no additional nodes
are added in creating the hierarchical wavelet-based description of a function. This leads to a
compact representation technique. (b) The wavelet transform is a perfectly invertible transform
and therefore, there is no loss of information in this representation. (c) The model stiffness matrix
in the wavelet basis is nearly diagonal which leads to a very fast numerical solution to the model
fitting problem. (d) The wavelet-based representation allows us to synthesize fractal surfaces of
arbitrary order without having to resort to blending approximations. These are but a few of the
list of advantages of the wavelet-based representation.
Some of the reasons for the nonequivalence of the hierarchical basis finite element representation
i.e., the multi-level splitting of the finite element spaces as in [51] and the multiresolution wavelet-based
representation are, using a hierarchical basis for the finite element spaces improves the
condition number (over the one obtained with the usual nodal bases) of the model stiffness matrix.
Whereas, a multiresolution orthonormal wavelet basis does not change the condition number of
the stiffness matrix since, the wavelet transform is performed via a similarity transformation (a
transformation that leaves the Eigen values and hence the condition number unchanged). Even if
we use an orthogonal basis in the hierarchical finite element method, it would still not be equivalent
to the multiresolution wavelet representation because, not all orthogonal bases are wavelets. Finally,
unlike in the construction of a wavelet basis representation, there is no high-pass filtering operation
in the construction of a hierarchical basis finite element representation.
Unlike the windowed Fourier transform, where the analyzing sinusoid function (kernel) exists on
an interval of fixed length, the analyzing wavelet's (``mother wavelet'') time support reduces with
increasing frequency. Wavelet analysis therefore works best for functions where low frequency components
are stationary over longer time intervals while high frequency components are stationary
over shorter time intervals. Most of the images/signals encountered in practice actually exhibit
these characteristics which accounts for the increasing popularity of wavelet analysis. For our purposes
of modeling the deformable superquadrics, the choice of wavelet basis not only facilitates a
hierarchical representation of the hybrid model surface but, also provides good visual appearance due
to this natural trade-off between time and frequency resolutions. Thus, the multiresolution hybrid
models can be used to synthesize a variety of shapes that look natural in the context of interactive
sculpting. We present some examples of these shapes in section 2.
Earlier work reported in Pentland [30] on the use of wavelets in the context of shape modeling
was limited to computational issues in surface fitting applications. 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 thin-plate-membrane splines. The main goal 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. Wavelet 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 [33], Sclaroff et al., report techniques for generalizing implicit function representation of shapes
by use of modal deformations and displacement maps. The displacement map function is defined
by using a bilinear interpolation between point samples and is stored as an image array. When
smoother interpolation is desired, they propose to employ a spline or pyramid-based multiresolution
method [12]. Wavelets or quadrature mirror filter pyramids are used for data compression and hence
savings in storage of the displacement maps. Also, the modal deformations are assumed to be the
vibration modes of the object which simplifies the simulation of the dynamics of the object. The
multiresolution pyramid method discussed in [12] and proposed for use in [33] has many limitations
compared to a wavelet-based multiresolution representation as pointed out in Mallat [23]. Note
that the wavelet-based pyramid is used in their work only for the data compression purposes. The
implications of this wavelet-based representation on improving the expressive power of the modeling
scheme are however left unexplored.
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. The reason for being able to incorporate a larger class of shapes as
priors in shape extraction applications is due to the enhanced class of shapes that we can now
characterize (in our modeling scheme) using a small set of global parameters. Another salient
feature of our modeling scheme is that, it allows for smooth transitions in the spectrum of model
deformations and parameters. In addition, we also emphasize on the efficient computation of the
model fitting process to the data. Specifically, we give a closed form expression for the stiffness
matrix computation in the wavelet basis (see section 4.2) and employ a multiresolution relaxation
process in solving the model fitting problem.
Thin-plate-membrane splines were also used for a very different purpose by Szeliski et. al.,
namely, generation of constrained fractals which are useful in representing natural images.
They derive the fractal characteristics of the membrane and thin-plate splines and show that the
membrane splines give raise to fractals of while the thin-plate splines yield fractals of
A fractal of in-between fi (order) was approximated by using a blend of the membrane and thin-plate
models. The modeling scheme we introduce in this paper allows us to generate fractals of
in-between order without having to resort to approximations with a blend of the plate membrane
models. We achieve this by modulating the frequency characteristics of the wavelet basis with
a rational function (which is a reciprocal polynomial) of appropriate order. To synthesize fractal
surfaces, we convert the energies into probability distributions using ideas from statistical mechanics
[20] and employ a multi-resolution stochastic relaxation algorithm. We will present some examples
of the synthesized fractal surfaces in a subsequent section.
Like any deformable surfaces, our multiresolution hybrid models can also be interactively manipulated
to mold into desired shapes via the use of externally applied forces. In this paper, external
forces will be derived from volumetric medical image data (MRI brain scans) and the model will
be forced into desired anatomical shapes in the MRI data of a human brain. The shape recovery
from data is achieved in two stages namely, a learning phase and an operational phase . The shape
recovery process involves the use of prior information about the anatomical shape of interest in the
MRI brain scan. This prior information is obtained via a supervised learning procedure wherein,
statistics on the augmented global parameters of the multiresolution hybrid model are collected. In
doing so, the model is subject to external image based forces and also forces derived from the interactive
placement of data points by an expert (neuroanatomist) on the boundary of the anatomical
shape of interest. The expert anatomist has the freedom to interactively manipulate the model
until a "satisfactory" fit to the data is obtained. The learning of a specific shape is carried out over
a fairly large data set and statistics for the augmented global parameters of the hybrid model are
collected. This information is utilized in the initialization of the model during the operational phase
of the shape recovery. We give a recursive algorithm for collecting the statistics on the augmented
global parameters of the model in section 5.
The rest of the paper is organized as follows: In section 2 we briefly describe the geometry
of the modeling scheme and the multi-resolution basis for expressing the geometry, followed by a
new way to augment the global shape of the deformable superquadric primitive for achieving a
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 model mechanics. In
section 4 we derive the internal energy of the prior model in a 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 model fitting to MRI data
followed by conclusions in section 6.
Multi-resolution Geometry of the Modeling Scheme
In this section, we will briefly describe the geometry and the construction of a multi-resolution
wavelet basis for the modeling scheme. For a detailed description of the deformable superquadric
geometry, we refer the reader to Terzopoulos [42]. To dramatically increase the library of shapes
that can be characterized by a small set of global parameters, and to achieve a smooth transition
from global to local shape descriptions, we will present a new technique that will augment the
global parameters of the deformable superquadric.
2.1 Geometry of the Deformable Superquadrics
The deformable superquadrics are closed surfaces in space whose intrinsic coordinates
are defined on a
domain\Omega\Gamma The positions of points on the model in an inertial frame of reference are
given by a vector-valued function In a model-centered coordinate
frame, the position vector x becomes
where c is the location of the center of the model, and the rotation matrix R specifies the orientation
of the model-centered 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, 1). For the parameterized reference shape
s(u) we use the superquadric with a bending deformation (as in Bajscy et al., [2]), due to it's
attractive global shape characterization. The reference shape s is given by,
frame
Model
Deformable
y
x
z
s
d
c
x
Z
Y

Figure

1: Deformable superquadric geometry (adopted from [42])
Where B is a bending deformation, a function of two parameters: the radius of curvature k and
angle of the bending plane ff (see [2] for details). The parametric equation of a superquadric ~ s is
given by
a 3 S ffl 1
uC A (3)
Where, \Gamma-=2 - u -=2 \Gamma- v -, and S ffl
are aspect ratio parameters, ffl 1 are "squareness" parameters and
a - 0 is a scale parameter. 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 s . We will denote c by the vector q c and use q ' to denote the vector of rotational
coordinates of the model in a quaternion representation (see [42]).
An important aspect of our superquadric reference shape which differs from previous definitions
in the work of Terzopoulos et al., [42] and Bajscy et al., [2] is that, we have constructed a uniform
tessalation of the superquadric. Superquadric shapes exhibit anomalous behavior when the square-
ness parameters ffl 1 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
(a) (b)

Figure

2: Tessalations of a Superquadric, (a) nonuniform (b) uniform
spacing along the parameter direction v. Thus, we seek a function of the parameter u (same reasoning
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
Setting using
results in an ordinary differential equation for f(u)
df
du
a 2
as the initial condition (6)
Where K is a constant. This ordinary differential equation can be solved numerically. In figure 2, a
superquadric with its natural nonuniform tessalation is shown side by side to the same superquadric
uniformly tessalated using the method described above. Very fast retessalation of the superquadric
can be achieved using the Runge-Kutta method [15]. Specifically, we use an order three Runge-Kutta
method. The retessalation needs to be done whenever the parameters of the superquadric
are subject to large (greater than a preselected threshold) changes.
In the following sections, we show that the displacement vector function completely
determined by a set of wavelet coefficients collected into the vector q d , with S corresponding to
the discrete wavelet transform and \Phi containing nodal interpolation functions. Thus, the state of
the model is fully expressed by the vector
d
2.2 Multiresolution Analysis and Wavelet Bases
Multiresolution analysis has been quite popular in signal processing and computer vision research
for the past decade. It has been used quite extensively for speeding up computations in solutions to
problems in computational vision [39] and in computer graphics for generating stochastic fractals
[36]. With the advent of the wavelet transform (WT), multi-resolution analysis of functions has
become more attractive due to the fact that wavelet decomposition of a function into different
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 multi-resolution
approximation in 1D and its extension to a 2D function specifically the membrane displacement
function d(u; v).
2.2.1 Detail Function from Multi-resolution Approximations
A multi-resolution approximation of a scalar valued function d(u) 2 L 2 (!) at a resolution j
is obtained by projecting d(u) into a vector subspace
of a dilated scaling function OE(u) form the orthonormal basis of V j . The
approximation -
ff jk OE jk (u) (7)
Where the set of coefficients ff corresponds to the discrete characterization
of this approximation.
An important property of multi-resolution approximations is that all the necessary information
to compute the function at a smaller resolution contained in the function at resolution j.
In other words the subspaces V j satisfy the containment hierarchy,
The scaling function OE has the effect of a low-pass filter with the cutoff frequency decreasing
proportionally to 2 \Gammaj (see figure 3(a)). The key idea of the wavelet basis is to span the difference
spaces W j+1 between the two approximations V j and V j+1 such that W j+1 ? V j+1 . A wavelet
decomposition therefore consists of splitting the function at resolution j into the detail function in
W j+1 and the next lower (j+1) resolution approximation in V j+1 . This leads directly to a scale
wise, orthogonal decomposition of L 2 (!). To represent the function d(u) of a finite resolution we
(a) (b)
f
low-pass filters
band-pass filters

Figure

3: Division of the frequency domain for (a) multi-resolution and (b) wavelet representation.
can assume
with \Phi denoting the direct sum operator. The hierarchical orthonormal basis for d(u) is given by a
family of dilations and translations of a wavelet function
d(u) can be written in this basis as,
Where the discrete wavelet coefficients are computed as
As in the case of the multi-resolution approximation, wherein a low pass filtering interpretation
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 band-pass filter (see Figure 3(b)). Our
work in this paper is primarily based on this filter bank interpretation of the WT as discussed in
Rioul et. al. [32].
2.2.2 Extension to 2-D and Fast Implementation
Extension to 2D of the multi-resolution approximation using wavelet basis can be easily achieved
using separable wavelets obtained from the products of one-dimensional wavelets and scaling functions
[23]. For the multi-resolution approximation of the membrane d(u; v) we use the tensor
products of OE and / given by ( see Mallat [23]),B @
(j;k;l)2Z
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)
d
d
d
d
d
d
d
d
d
3
HH
HH
HH
a 3
LHd
d
HH
d
d
d LH
a
d dl

Figure

4: a) Organization of the displacement vector q d , b) separation of the "global" and "local"
parts of q d
represents high (low) pass filtering in each of the directions u and v in the parameter
domain\Omega\Gamma and
the indices k and l 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., ffi LH
jkl
and
jkl into a vector q d as illustrated in Figure 4. This displacement vector q d contains a complete
hierarchical description of the membrane surface. Note that q d has the same number of elements
as the original nodal discretization at our WT implementation we adopted the specific
OE(u) and /(u) functions as given in Mallat [23], which are based on the cubic spline interpolants.
If we denote the original nodal discretization at resolution by a vector ff 0 then, we have,
Reconstruction
where \Phi 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 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 5 (see Mallat [23] for details).

Figure

5 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 q d is expressed in an orthonormal wavelet basis. Along with other parameters of
columns
rows
convolve with filter X
keep one column out of two
keep one row out of two
G
G
G ~
~
~
~
~
~ 11122212a
a
d
d
d
HH
convolve with filter X
G
G
G
. 421put one row of zeros between each row
put one column of zeros between each column
multiply by 4
HH
a
d
d
d
a
columns

Figure

5: a) Decomposition and b) reconstruction of a 2-dimensional function, adopted from [23]
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 4, we explicitly depict the hierarchical decomposition of the displacement vector q d in an
orthonormal wavelet basis. The coefficients toward the upper left corner are fewer than those
toward the bottom-right 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 "semi-global"
character. Depending on the application we can draw the boundary at an arbitrary j 0 to separate
the global and local part of the vector q d i.e.,
where q dg contains the top (coarser) portion of the multi-resolution pyramid. Since the separation of
the wavelet decomposition into "global" and local parts depends on the application problem, it can
be achieved via a supervised learning procedure that can be used to cluster the "global" parameters
which well describe the mean shapes of interest. Given a good population size for each shape of
interest, based on the clustering of the parameters, we can choose the frequency j 0 to separate the
global and local parts of the vector q d . Our future research will focus on implementation of the
parameter clustering technique for choosing the separation frequency j 0 .
One of the useful aspects of including q dg 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 6). 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 6(a) depicts the original reference shape, global
deformation resulting from tweaking a superquadric parameter is shown in (b). The bottom four
shapes (c) through (f) are generated via global deformations caused by tweaking the coarse level
wavelet coefficients. Note that, tweaking these coarse-level wavelet coefficients produces global deformations
that are very different from those produced by tweaking the global (lumped) parameters
of the model.
When we build our prior model with mean values and covariances on q s and q dg , 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 shape recovery problems encountered in Medical
Imaging applications eg. recovery of anatomical shapes from MR images. We will present two such
examples in this paper.
3 Probabilistic Framework and Model Mechanics
In the last section, we described our multiresolution deformable superquadric model in a wavelet
basis. This basis provides a computationally efficient solution to the model fitting problem 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 recovered from the data and then derive the mechanics of our modeling scheme.
We cast the model fitting process in a probabilistic framework [35] and incorporate prior distributions
on the vector of geometric parameters that are being estimated. 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 physically-based 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 (supervised learning) described in a subsequent section and for the local parameters, we
convert the continuous strain energy - which governs the deformation of the model away from
its natural shape - into a probability distribution over expected shapes. This is achieved via a
(a) (b)
(c) (d)

Figure

deformations generated via tweaking global parameters of the
superquadric model, (c), (d), (e) & (f) global deformations generated by tweaking of coarse scale
wavelet coefficients
technique of statistical mechanics - conversion of energies into probabilities using Boltzmann, or
Gibbs, distribution [20]. The link between physically-based (deformable models) models and priors
is conveniently established using a Gibbs distribution of the form [36]
is the discretized version of the internal smoothness energy of the model, and Z p (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 [25]
exp(\GammaE D (D; q)); (15)
where ED (D; q) can be either an edge-based potential energy synthesized from image data 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
where
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 physically-based model [49]. 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 E p of the model in our case is given by a quadratic form,
is the rest state of the model and the matrix K p corresponds to
the stiffness of deformation. We will discuss the internal energy E p further in a subsequent section.
The two types of sensor model energies for ED (D; q) discussed earlier are expressed as
ED (D;
spring energy
R\Omega fi jP (x(q; u))j image potential (18)
In the above equations u i determines the closest point x to the given measured point D i . We
adjust the locations of these grid coordinates as the shape of the model evolves (see [42, 46], for
more on migrating point of influence). The magnitude of parameter fi 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 a posteriori amounts to minimizing the total energy
for which we employ a gradient descent method that amounts to iteratively solving,
where   is a diagonal matrix containing step sizes and consequently effecting the speed of conver-
gence. For the derivative of the energy on the right hand side in the above equation, we use the
chain rule to get @E
@x
@x
, where the contribution of each parameter group - using the notation
in [42] - is
@x
R@s
@q d
Partial derivatives of the prior model and data energies defined above can be now written as
and
@ED
spring energy
R\Omega fir jP (x(q; u))j 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
Close similarity of this evolution equation to the one derived in [42] reveals that equation 25 can
be modified slightly and effectively used to generate motion sequences corresponding to first order
dynamics C -
7). The damping matrix C determines the velocity
dependent dissipation of the kinetic energy which governs the exponentially decaying motion of the
model toward the equilibrium state. For simplicity, the equations can be decoupled by assuming
the C matrix to be diagonal and constant over time. In addition to producing the statistically most
likely state -
q, the Bayesian framework provides the uncertainty in the estimate via the associated
covariance matrices.
(a) (b) (c) (d)

Figure

7: Evolution of a deformable superquadric, (a) initialization, (b),(c), and (d) configurations
of the model after 10, 20, 40, iterations, respectively.
(a) (b) (c)

Figure

8: Coarse to fine model fitting, from (a) through (c) resolution of the membrane increases.
While the above iterative equation will eventually converge to the correct estimate, in practice
it may be unacceptably slow. To accelerate convergence we use coarse-to-fine search on a multi-resolution
pyramid corresponding to the wavelet basis. The problem is first solved on a coarser
mesh, then this solution serves as a starting point for the next finer mesh. Figure 8 shows the
evolution of the model at three successive resolutions. This was achieved by incorporating wavelet
coefficients of increasingly finer levels in equation 25. Besides improving computational efficiency,
this framework leads to significant smoothing of the energy potential at coarser resolutions which
consequently reduces the problem of 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 earlier.
4.1 Internal Energy of the Prior Model
In addition to generic smoothness assumptions, a prior model should incorporate the information
about the possible shapes that an observed object can assume. This can be very important in
applications involving recovery 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 configurations with large energy E i and vice versa.
The configuration with 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 (see [42]),
Z
du
where, the functions ! 0 (u) and ! 1 (u) control the local magnitude and variation of the deformation
respectively. In this paper, we will assume these functions to be known scalar parameters.
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 we have an equivalent expression to equation 26,
Z
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 -
corresponding
to parameters q dg as defined in equation 13 and shown in figure 4. This is done by altering the
deformation energy equation 26 to,
Z
d)
d)
d)
The discrete version of this internal strain energy can be written as,
Where, K d is the stiffness matrix that determines the elastic properties of the model and -
incorporates the statistics of the possible shapes that can be encountered in the
shape recovery problem. These statistics are gathered through a training process described in
section 5. It is evident that the number of parameters in -
q dg 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,
Where ~ s denotes the superquadric shape without the bending while, Cs represent the covariance
matrices which are accumulated in the iterative training phase. In section 4.2, we give the derivation
for the 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
relates directly to the diagonal entries of the stiffness matrix K. This approach in turn gives us
a straight-forward and inexpensive way to compute K.
The standard way to obtain stiffness matrix in wavelet basis is essentially a similarity trans-
formation. By combining equations 12, 28 and 29 and for simplicity, depicting the 1D case we
get,
'Z
@ \Phi
du
Here K 0
d controls the local magnitude and K 1
d the local variation of the (membrane) displacement
function. Due to orthogonality of the interpolation functions in \Phi and orthogonality of S, equation
simplifies to K 0
I . However, the part K 1
d of the stiffness matrix in equation 33 needs to
be computed by a similarity transformation K 1
represents the membrane
membrane

Figure

9: (a) Scaling of the diagonal entries in the stiffness matrix K of the membrane model with
spectral density, ! \Gamma2 and an arbitrary density ! \Gammafi , (b) Increasing energy with increasing resolution
stiffness matrix in nodal basis (see [30]). The drawback of computing K 1
d using this similarity
transformation is the large computational cost incurred (O(N 3 ), 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 dominance (see [30]). 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),
which is done in frequency domain by employing Parseval's theorem [11]. Extension to 2-D case
is straightforward and achieved by using the tensor product of the basis functions. The factor j!
in equation 34 results from the Fourier transform of the first derivatives representing the membrane
regularizer. Note that both K 0
d
and K 1
d
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 explains how coefficients corresponding to lower resolutions
contribute less energy than those corresponding to higher resolutions. Figure 9 illustrates two
aspects of the modeling, (a) the membrane displacement function which has a spectral density
characterized by ! \Gamma2 and is a fractal process of 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 order (fi) by using a function ! \Gammafi 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 band-pass characteristic for \Psi j (!) between ! j and 2! j we would have the following
K K
HH
HH

Figure

10: Diagonal entries of K as an elevation map, for a (32X32) grid computed via a similarity
transformation
relationship between elements of K 1 at j and (j
Z 2!j
Therefore, for the membrane regularizer, we notice that the energy at a resolution j is eight times
that at the next lower resolution (j 1).
For a verification of the above discussed scaling behavior, we computed the diagonal entries of
the K 1
d matrix for a (32 X 32) grid of nodes using a similarity transformation technique discussed
depicts the diagonal entries as an elevation map, where the elevation corresponds
to the magnitude of the entries. The entries are organized as in a standard representation of images
in wavelet basis (see figure 4 for the standard organization).
4.3 Prior Model of General Fractal Degree
Splines and stochastic fractals are complementary techniques for modeling free-form shapes [36].
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. There
is abundant literature on fractal generation algorithms and we refer the reader to relevant literature
in [19, 48, 3]. However, these algorithms do not describe techniques to constrain fractals. Szeliski
et. al. [36] describe a novel technique for constraining fractals by establishing a formal connection
between splines and fractals using Fourier analysis of energy minimizing splines.

Figure

Model fitting result (a) MAP Estimate and (b) Representative sample from the a
posteriori distribution, corresponding to a fractal surface (membrane model)
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 [36]. In their work
on constrained fractals Szeliski et. al. [36] show how prior models corresponding to a membrane
and thin-plate regularizers possess a spectral density proportional to ! \Gamma2 and ! \Gamma4 respectively,
which are fractal i.e., they are self-affine over scale. These spectra characterize special cases of
fractals of general form possessing a density given by,
The spectral density of a stochastic fractal obeys the inverse power law in equation 36 and characterizes
a fractal Brownian function with fractal dimension is
[48]. Here, E denotes the dimension of the Euclidean space.
In [36] fractals of an in-between fi namely, 2 - fi - 4 were approximated by blending the
membrane and thin-plate terms. In contrast, our multi-resolution wavelet based modeling scheme
discussed earlier allows us to construct fractals of any fi 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 replaced by ! fi . Thus, yielding a stiffness matrix corresponding to
a regularizer whose spectral distribution is a rational function of order fi. A random sample from
the corresponding Gibbsian distribution of this prior model will be a fractal surface of order fi.

Figure

11 shows a surface fitting result to 3D point data set with (a) MAP estimate and (b)
representative sample from the Gibbs distribution p(qjD) depicting the underlying fractal prior.
For the purpose of visually comparing surfaces of various fractal orders fi, we choose to normalize
the cumulative energy contributed by stiffness matrix K. Normalization was achieved by setting
the Frobenious norm of K to be constant for various fi. As a result, the only factor governing the
behavior of fractal surfaces in figure 12 is the relative distribution of energy along the frequency axis
(a) (b) (c)

Figure

12: Fractal surfaces of varying degree fi; (a),(b) & (c) -
according to curve ! \Gammafi . As expected, a higher value of fi emphasizes low frequency content and
deemphasizes higher frequencies. These surfaces were obtained by taking a representative sample
from the associated Gibbsian distributions that describe the surface probabilistically. The explicit
evaluation of the Gibbs distribution is very difficult due to the large computational requirements
(summation overall possible states) in evaluating the partition function Z. As in [36], this difficult
computation can be circumvented by generating random samples from the Gibbs distribution using
a simple technique called the Gibbs Sampler [20]. In the Gibbs Sampler technique, at each iteration,
a new random state is selected from the Gibbs distribution associated with the local energy function
of the variable being updated. The convergence of this updating scheme to a random sample
from the overall distribution is guaranteed. However, the rate of convergence may be very slow.
We employ the multiresolution wavelet-based pyramid to accelerate the convergence by using the
results at a coarse level as the starting point of the relaxation for a finer level.
The implementation of the coarse-to-fine relaxation in a wavelet basis is quite different from
that discussed in [36]. Specifically, after transformation of the nodal basis to the wavelet basis, the
solution is obtained at a coarse resolution, to this solution, we add fine detail (detail function) and
then recompute the new solution and so on. White noise is added to the surface at each iteration in
a given level. The retransformation to nodal basis is done only after the entire solution is obtained
at the finest level in the wavelet basis. Strictly speaking, if the model stiffness is perfectly diagonal
then, the multiresolution solution can be obtained directly in one step. In our model, the stiffness
matrix (inverse covariance) is for all practical purposes diagonal and therefore, neglecting the close
to zero off diagonal terms, the white noise was added to the final solution surface to produce a
representative sample (fractal surface) from the corresponding Gibbs distribution.
Fractal surfaces are very useful for modeling natural detail. In our application of medical imag-
ing, our goal is to create realistic shape models of the anatomical structures found in the human
body, specifically the brain. Future experiments will focus on achieving this goal.
5 Supervised Learning and Model Fitting Results
From an application point of view, our goal is to apply the new multi-resolution shape modeling
scheme described in previous sections for recovering the 3D shape of cortical and subcortical structures
such as the gyrus and hippocampus, from brain MRI data. The shapes will be accurately
quantified by a parameter vector that can be used in computing shape measures such as volume.
Volumetric analysis of structures contributes to understanding of various neurological disorders
such as dyslexia, specific language impairement etc. (Leonard, et al., [22]). Nontrivial shapes, low
contrast images and high noise levels in this application domain prompted us to take advantage of
our statistical modeling scheme by incorporating object-specific information into our prior model.
To facilitate construction of such specialized prior models we divide shape recovery process into
a training phase and an operational phase. In the training phase human supervision via interaction
and manipulation is needed to provide good training samples, while in the operational phase
the model-fitting process becomes fully autonomous and relies heavily on statistics accumulated in
prior model.
Shape recovery from poor quality data can be formidable. Any approach that does not make
use of prior knowledge of the shape of interest is doomed to fail due to the underconstrained nature
of the problem. Imperfect data can be augmented with a priori information that a model of the
structure of interest provides. For example, when dealing with specific anatomical shapes, one
can use models of these shapes to guide the recovery algorithms. In the earlier sections, we have
described a new modeling scheme that captures global as well as semi-global shape using a small set
of parameters. This key attractive feature of the modeling scheme makes it conducive to perform
statistical data collection on global model parameters which are then used (in combination with
the Gibbsian probability distribution of the local part of the model) as prior information in the
model fitting process.
As discussed in section 3, Bayesian modeling requires knowledge of the a priori probability p(q)
of the model q and the probability of q conditioned on the sensed data i.e., p(q j D). Unfortunately,
in surface modeling applications we rarely if ever have this kind of complete knowledge about the
probabilistic structure of the problem. Typically, we have a number of representative samples of
what we want to model. In an approach known as Bayesian learning [18], these samples are used
Set of training data samples
Prior Model p( q )
and covariances of the global model
obtained from sensor
Surface fitting by maximizing
current prior
model p( q ) is used.
Iterative improvment of the mean

Figure

13: Supervised learning to incorporate information into the prior model
to estimate the unknown probabilities which in turn are used as if they were the true values. It
is important to distinguish between supervised learning and unsupervised learning. The distinction
being, in supervised learning, we know the class (label) of each sample in advance, whereas in
unsupervised learning, the natural clustering of samples is used to define the classes. In our work,
prior knowledge about the object shapes that are to be recovered from the MRI data is collected
via a supervised learning technique. A training sequence of several iterations provides samples of
correctly reconstructed shapes for specific classes of objects - in our experiments a gyrus and the
hippocampus from the human brain. Parameter sets q(i) corresponding to these shapes represent
the training sequence for the global part (q s our prior model (see figure 13). 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 13. Since we update only the
statistics of the global parameters of the model, the local deformations remain 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 model fitting in the training stage to assure that only good model fitting
results enter the statistics of the prior model. An extension to multi-class case is straightforward,
and can be achieved by defining multiple prior distributions p j (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 the model fitting process to a gyrus and a hippocampus
shape, given MRI data of a human brain. The data set consisted of 64 slices, with each slice of
1.25mm thickness. Figure 14, illustrates our results for a gyrus. From left to right and top to
bottom, the displayed images are, (a) a sagittal section of the brain showing a region of interest
(containing a gyrus) indicated in a box, (b) a slice of the initialized 3D shape model superimposed
on the corresponding slice from the original data, (c) superimposition of a slice from the final recovered
shape on a corresponding slice from the original data, and (d) the 3D recovered shape.
The initialization shown here has incorporated statistical information (mean and variance) of the
gyrus from a small sample data set. The initialized model was subject to externally applied forces
synthesized from a 3D edge based potential obtained by computing edges in 3D using the recursive
filtering discussed in [27]. Figure 15 depicts the reconstruction results for a hippocampus, arranged
left to right top to bottom in sequence. Figure 15(a) is a slice of the raw MRI data showing
the region of interest (ROI), (b) is a slice from the 3D model initialization superimposed on the
ROI, (c) depicts a slice from the model initialization superimposed on the edge based potential field
computed using the recursive filtering method in [27], (d) is a slice from the recovered 3D shape
superimposed on the ROI, and (e) depicts the recovered 3D shape of the hippocampus.
Further experimentation with our modeling scheme is currently under way and we expect to
better our results by using, (a) larger data sets for collecting the statistics, (b) better 3D edge
finding methods, and (c) adaptive clustering for automatically choosing the cut off frequency j 0 in
separating the global and local parameters of the model.

Figure

14: 3D Shape recovery results for a gyrus.
6 Conclusion
In this paper, we introduced a new shape modeling scheme that can transform smoothly from
local to global models or vice-versa. The modeling scheme utilizes a hybrid primitive called the
deformable superquadric constructed in an orthonormal wavelet basis. These basis provide the
power to continuously transform 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. In addition, they also allow to smoothly scale up and down the spectrum of models
with lumped parameter models at one end and distributed parameter models at the other. The
multi-resolution wavelet basis allows us to generate fractal surfaces of arbitrary order which can be

Figure

15: 3D Shape recovery results for a hippocampus
useful in describing natural detail.
When embedded in a probabilistic framework, our modeling scheme allows us to incorporate
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 detail may also be
incorporated by using some existing modeling techniques [34, 42] 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 multi-resolution relaxation is a natural by product
of implementing the gradient relaxation method in the wavelet basis.

Acknowledgments

The authors would like to thank the reviewers of this paper for their useful and insightful comments
and Mr. S. H. Lai for the fruitful technical discussions and his help in preparing the figures in
LATEX.



--R

"Wavelets and other bases for fast numerical linear algebra,"
"Three-dimensional object representation revisited,"

"Superquadrics and angle-preserving transformations,"
"The Beta-spline: A local representation based on shape parameters and fundamental geometric measures,"
"Three-dimensional object recognition,"
"Visual Perception by computer,"

"On surface reconstruction methods,"
"Visual surface reconstruction using sparse depth data,"

"A multiresolution spline with application to image mosaics,"
"Deformable curve and surface finite-elements for free-form shape design,"
"Deformable models for 3D medical images using finite elements and balloons,"

"Orthonormal bases of compactly supported wavelets,"
"Shape represnetation and image segmentation using deformable surfaces,"

"Computer rendering of stochastic models,"
"Stochastic relaxation, Gibbs distribution, and Bayesian restoration of images,"
"Surface reconstruction from unorganized points,"
"Angelman and Prader- Willi Syndrome: A Magnetic resonance imaging study of differences in cerebral structure."
"A theory of multi-resolution signal decomposition: the wavelet repre- sentation,"
"Singularity detection and processing with wavelets,"


"3D edge detection using recursive filtering: Application to scanner images,"
"Volumetric description of range data using the Blobby Model,"
"Parts: Structured descriptions of shape,"
"Fast solutions to physical equilibrium and interpolation problems,"
"Representations of rigid solid objects,"
"Wavelets and signal processing,"
"Generalized implicit functions for computer graphics,"
"Boundary finding with parametrically deformable contour models,"

"Constrained fractals,"
"Surface modeling with oriented particle systems,"
"Regularization of inverse visual problems involving discontinuities,"
"Image analysis using multigrid relaxation methods,"
"Elastically deformable models,"
"Constraints on deformable models: Recovering 3D shape and nonrigid motion,"
"Dynamic 3D models with local and global defor- mations: Deformable superquadrics,"
"Curvature-based representation of objects from range data,"
"Deformable models: Canonical parameters for invariant surface representation and multiple-view integration,"
"Constructing intrinsic parameters with active models for invariant surface reconstruction,"
"Invariant surface and motion estimation from sparse range data,"
"From global to local, a continuum of shape models with fractal priors,"
"Random fractal forgeries,"

"Surface reconstruction using deformable models with interior and boundary constraints,"
"On the multi-level splitting of finite element spaces,"
"Feature extraction from faces using deformable templates,"
--TR
A multiresolution spline with application to image mosaics
Three-dimensional object recognition
Regularization of inverse visual problems involving discontinuities
On the multi-level splitting of finite element spaces
Image analysis using multigrid relaxation methods
Curvature-based representation of objects from range data
Visual reconstruction
Elastically deformable models
Constraints on deformable models: recovering 3D shape and nongrid motion
A Theory for Multiresolution Signal Decomposition
On Three-Dimensional Surface Reconstruction Methods
Dynamic 3D Models with Local and Global Deformations
Geometrically deformed models
Volumetric shape description of range data using MYAMPERSANDldquo;Blobby ModelMYAMPERSANDrdquo;
Generalized implicit functions for computer graphics
Deformable curve and surface finite-elements for free-form shape design
Surface reconstruction from unorganized points
Surface modeling with oriented particle systems
Boundary Finding with Parametrically Deformable Models
Wavelets and other bases for fast numerical linear algebra
Representations for Rigid Solids: Theory, Methods, and Systems
Computer rendering of stochastic models
Constructing Intrinsic Parameters with Active Models for Invariant Surface Reconstruction
The beta-spline
Bayesian modeling of uncertainty in low-level vision

--CTR
Baba C. Vemuri , Yanlin Guo, Snake Pedals: Compact and Versatile Geometric Models with Physics-Based Control, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.22 n.5, p.445-459, May 2000
Baba C. Vemuri , Chhandomay Mandal , Shang-Hong Lai, A Fast Gibbs Sampler for Synthesizing Constrained Fractals, IEEE Transactions on Visualization and Computer Graphics, v.3 n.4, p.337-351, October 1997
Pak Chung Wong , R. Daniel Bergeron, Authenticity Analysis of Wavelet Approximations in Visualization, Proceedings of the 6th conference on Visualization '95, p.184, October 29-November 03, 1995
Baba C. Vemuri , Chhandomay Mandal, A fast Gibbs sampler for synthesizing constrained fractals, Proceedings of the 7th conference on Visualization '96, p.29-35, October 28-29, 1996, San Francisco, California, United States
Baba C. Vemuri , Yanlin Guo , Zhizhou Wang, Deformable Pedal Curves and Surfaces: Hybrid Geometric Active Models for Shape Recovery, International Journal of Computer Vision, v.44 n.2, p.137-155, September 2001
Douglas DeCarlo , Dimitri Metaxas, Blended Deformable Models, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.18 n.4, p.443-448, April 1996
Gene C.-H. Chuang , C.-C. Jay Kuo, Cartoon Animation and Morphing with Wavelet Curve Descriptor, Multidimensional Systems and Signal Processing, v.8 n.4, p.423-447, Oct. 1997
Hong Qin , Chhandomay Mandal , Baba C. Vemuri, Dynamic Catmull-Clark Subdivision Surfaces, IEEE Transactions on Visualization and Computer Graphics, v.4 n.3, p.215-229, July 1998
Chhandomay Mandal , Hong Qin , Baba C. Vemuri, Dynamic Modeling of Butterfly Subdivision Surfaces, IEEE Transactions on Visualization and Computer Graphics, v.6 n.3, p.265-287, July 2000
Dimitris N. Metaxas , Ioannis A. Kakadiaris, Elastically Adaptive Deformable Models, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.24 n.10, p.1310-1321, October 2002
surface fitting;multiresolution representation;stiffness matrix;fractal surfaces;superquadrics;orthonormal wavelet basis;deformable surfaces;bayesian estimation