Group Title: Department of Computer and Information Science and Engineering Technical Reports
Title: Snake pedals : compact and versatile geometric models with physics-based control
CITATION PDF VIEWER THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00095426/00001
 Material Information
Title: Snake pedals : compact and versatile geometric models with physics-based control
Series Title: Department of Computer and Information Science and Engineering Technical Reports
Physical Description: Book
Language: English
Creator: Vemuri, Baba C.
Guo, Yanlin
Affiliation: University of Florida
Princeton University
Publisher: Department of Computer and Information Science and Engineering, University of Florida
Place of Publication: Gainesville, Fla.
Copyright Date: 1998
 Record Information
Bibliographic ID: UF00095426
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.

Downloads

This item has the following downloads:

1998279 ( PDF )


Full Text









Snake Pedals: Compact and Versatile Geometric Models

with Physics-Based Control *t


Baba C. Vemuri and Yanlin Guo3





Abstract

In this paper, we introduce a novel geometric shape modeling scheme which allows
for representation of global and local shape characteristics of an object. Geometric
models are traditionally well suited for representing global shapes without local detail.
However, in this paper we propose a powerful geometric shape modeling scheme which
allows for the representation of global shapes with local detail and permits model
shaping as well as topological changes via physics-based control.
1 Ir proposed modeling scheme consists of representing shapes by pedal curves and
surfaces pedal curves/surfaces are the loci of the foot of perpendiculars to the tan-
gents of a fixed curve/surface from a fixed point called the pedal point. By varying
the location of the pedal point, one can synthesize a large class of shapes which exhibit
both local and global deformations. We introduce physics-based control for shaping
these geometric models by letting the pedal point vary and use a snake to represent
the position of this varying pedal point. I Ir model dubbed as a ---I .i. pI ..I' allows
for interactive manipulation via forces applied to the snake. We develop a fast numer-
ical algorithm for shape recovery from image data using this novel geometric shape
modeling scheme. I I algorithm involves the Levenber-:- i.l.i p l.ii it (LM) method in
the outer loop for solving the global parameters and the Alternating Direction Implicit
(ADI) method in the inner loop for solving the local parameters of the model. I IrI
combination for the global and local scheme leads to an efficient numerical solution to
the model fitting problem. We demonstrate the applicability of this modeling scheme
and efficiency of the model fitting methods via examples of shape synthesis and shape
estimation from real image data.



Index Terms Geometric i'-.dels, Snakes, Pedal Curves/Surfaces, Alternating Direction
Implicit method, Levenber- '. i, rquardt method, Schur Complement.
*This research was supported in part by the NSF grant IIS-9811042 and the NIH grants R01-RR13197
and RO1-LM05944.
submitted to the IEEE Trans. on PAMI
1B. C. Vemuri is with the Department of Computer and Information Sciences, University of Florida,
Gainesville, FL 32611 (email: 11 !!L,, ,-. [t!! 1 tl I
Y. Guo is with the Sarnoff Labs in Princeton (email: '..- .. !!








1 Introduction


.i modeling shapes is a fundamental constituent of computer vision as well as computer graph-

ics. Tli, 1, are numerous shape modeling techniques in computer vision and computer graph-

ics literature. In the context of shape modeling, there are two primary tasks of importance in

computer vision namely, shape recovery from image data and shape recognition. Tli, shape

recovery task imposes the requirement that the shape model be capable of representing fine

detail and hence requires that the model have large number of degrees of freedom. On the

other hand, shape recognition requires that the model representation be compact in terms

of the number of parameters so as to facilitate, easy storage of a large database of models

and easy matching for retrieval. Tli -, two requirements are potentially conflicting and de-

signing shape models that can possibly '-.1i-" these requirements has been the focus of

many researchers in the recent past in the computer vision community [1, 2, 3, 4, 5].

Geometric models (lumped parameter models) can represent a shape with very fetw param-

eters and hence are good for object recognition. On the other hand, physics-based models

(distributed parameter models) require a large number of parameters for their representation

and are well suited for shape recovery and reconstruction. In [5], Terzopoulos and :i i, .1,-.

developed a hybrid dynamic model a deformable superquadric model which combines the

descriptive power of geometric and physics-based models via superposition of a membrane

spline on a core superquadric. This hybrid modeling scheme allows for modeling shapes at

both ends of the shape spectrum (one end of which is occupied by lumped parameter mod-

els and the other by distributed models). Superimposing the free form deformation (FFD)

on a core superquadric allows representation of a variety of global shapes with local detail

as was reported in [1, 6, 7, 8]. This modeling scheme requires only a few parameters to

represent global shape deformations. In [9], Vemuri and Radisavljevic developed a multires-

olution, hybrid modeling scheme which represents the displacement function of a deformable

superquadric model in a wavelet basis. Pentland and Williams [4] reported a novel shape









modeling scheme which is based on modal analysis. Tli, number of parameters required for

this representation depends on the desired level of detail.

In [2, 10], a hybrid hyperquadric model was developed which adds an exponential term

to a hyperquadric primitive. Tlii addition of the exponential term to the hyperquadric

equation allows for local control of the shape generated. This modeling scheme allows the

synthesis of an arbitrary number of concavities. It differs from most of the schemes described

above because the models are geometric but allow for global as well as local deformations.

Ti, model representation is quite compact in terms of the number of parameters needed to

describe the shape. However, the computation of local concavities is nontrivial and human

interaction is needed in the model fitting.

Recently, O'Donnel et al. [11] developed a novel solid shape modeling scheme. This model

includes built-in offsets from a core base model. Local deformation over the scaled-offset

model is allowed via displacement vectors attached to the default model. '.ldel fitting to

data is achieved in a dynamic framework as in [5]. This modeling scheme is well suited for

tasks such as cardiac motion recovery from 1.I-i-, d `.li; images.

An interesting and topologically adaptable deformable model was developed in [3] by De-

Carlo and .1 ii,\,i. via the use of parameterized blending functions. T111-, models can

determine the number and position of the holes in an unknown object and accordingly

change the topology by use of appropriate error of fit criteria. ;i'li, [ I l-. techniques based

on the concept of level-sets have been reported on literature [12, 13]. Another modeling

scheme called "topologically adaptable -ii..!-" was introduced in [14] by :-.1 iii, 11 ,y and

Terzopoulos. This novel version of snakes allows for automatic change in topology if and

when required during shape recovery. In all these models, a large number of parameters are

required to describe the recovered shape.

.Idel fitting to desired features in the image data can be posed as a nonlinear minimiza-

tion problem and solved numerically. In the process of numerically fitting the model to the








data, Vemuri and Radisavljevic [9] used a multiresolution gradient descent (GD) method,

while Terzopoulos and '.1% i1.x1-. [5] used a first order Euler scheme to fit their dynamically

deformable superquadric model to the data. Both of these methods are slow since only a

small time step can be taken in each iteration in order to get a stable solution. A relatively

faster and more stable integration scheme based on the Runge-Kutta method was developed

by DeCarlo and :'.1 1..1-, [3]. In [15], Cohen introduced a two-step iterative scheme based

on auxiliary variables. This scheme is very general and can be applied to most optimization

problems. Ti, auxiliary variables may be interpreted as representing intermediary recon-

struction steps. :'.iii1ii.-,11i ii, of several non-convex potentials were transformed to convex

minimization problems. However, the transformation was shown primarily for potentials

involving a function of distance to the closest data and not image-based potentials.


1.1 Overview of Our Model and Numerical Algorithms


In this paper, we propose a powerful geometric shape modeling scheme which permits the

representation of global shapes with local detail using only a -, iii-global" characterization

and relatively simple and efficient numerical methods. Ti, proposed modeling scheme con-

sists of representing shapes by pedal curves and surfaces -pedal curves/surfaces are the loci

of the foot of perpendiculars to the tangents of a fixed curve/surface from a fixed point called

the pedal point. By varying the location of the pedal point, one can synthesize a large class

of shapes which exhibit both local and global deformations. We introduce physics-based

control for shaping these geometric models by letting the pedal point vary and use a snake

to represent the position of this varying pedal point. Ti, model dubbed as a -iil.I pedal"

allows for interactive manipulation via forces applied to the snake. TIl key contribution

of this paper is that our modeling scheme is inherently geometric but allows for local and

global deformations as would a physics-based model. Also, the model is compact, capable

of global deformations such as bending, twisting, and tapering without the introduction of

additional parameters. In addition, the modeling scheme allows for synthesis of objects of








varying topology.

To fit our model to a sparse set of data points in 2D/3D as well as to image data for shape

recovery, we develop a new nonlinear optimization scheme which contains the Levenberg-

:.irquardt (1 '.I) method in the outer loop for estimating global parameters of the model

and the alternating direction implicit (ADI) method in the inner loop for estimating the

local parameters of the model. Such a combination of a global and a local scheme is quite

stable and yields a locally optimal solution very quickly.


1.2 Organization of the Paper


Ti, rest of the paper is organized as follows. In Section 2 we give the definitions of pedal

curves/surfaces and show examples of the same in 2D and 3D. Following which we develop

the proposed models namely the -Cil; pe II!- and illustrate the ideas by synthesizing a

variety of shapes. In Section 3, we develop the mathematics of fitting the -i,,IIl; pedals" to

real image data and present the experimental results in Section 4 followed by conclusions in

Section 5. Ti, appendix contains details on the derivation of the stiffness matrix in a finite

element (FE) framework and its representation as a Kronecker product.


2 "Snake Pedals": The Proposed Model


Tlii definitions of pedal curves and surfaces presented here maybe found in any book on

elementary differential geometry text and we refer the reader to [16].


2.1 Pedal Curves


Let a be a planar curve, the pedal curve of a is defined as the locus of points on the foot

f of the perpendicular from a fixed point p called the pedal point to a variable tangent of

the planar curve a.

Suppose that the planar curve a is a parameterized curve with domain parameter t. Let 3





















Figure 1: f is on the pedal curve of a with respect to the pedal point p.


be the pedal curve of a with respect to the fixed pedal point p, and let a(t) = g, 3 (t) = f

at a particular t, as shown in Fig. 1. T l, projection of a (t) p in the direction Ja '(t)
must be ( (t) p, where a' (t) is the tangent line of plane curve a (t), J : R2 R i2 is a
linear map given by



J(,y) =(-y,). (1)

Note that applying the operation J on a vector results in another new vector and J can be

geometrically interpreted as a rotation by 7r/2 in a counterclockwise direction. We can thus
define a pedal curve as follows [16]:

Definition : T.,, pedal curve of a regular curve a : (c, d) R2 with respect to a fixed

(pedal) point p S 2 is given by



pedal(t) = p + ( p) J (t). (2)
| Ja 1(t)
In Fig. 2, we present examples depicting the pedal curves of an ellipse for different positions

of the pedal point as shown. Note that the pedal curve is capable of exhibiting global as
well as local deformations. Tli, global shape of a pedal curve of an ellipse is approximated

by a new ellipse, whose long axis is almost the same as that of the original ellipse, while



















(a) (b) (c) (d)

Figure 2: Examples of pedal curves (in green) of an ellipse (in red) for different pedal point
positions. Pedal points (in red) are shown by a dot in each case.


the short axis is much larger than that of the original ellipse. Ti, new ellipse is pinched in

the direction of the short axis to generate the pedal curve. Ti, location and extent of this

pinching are determined by the position of the pedal point. As shown in Fig. 2, if the pedal

point is located in the center of the original ellipse (Fig. 2 (a)), the pinching occurs around

the center part of the pedal curve, while as the pedal point moves to an upper position (Fig.

2 (b)), the pinching in the pedal curve moves up. Similarly, as the pedal point moves to the

right of the ellipse origin (Fig. 2 (c)), more pinching occurs on the right side in comparison

to the left of the origin. This pinching operation constitutes the local deformation of the

pedal curve, but the type of the local deformation that can be exhibited by a pedal curve is

by no means limited to the pinching type deformation. '. .1,re complex local deformations can

be realized in the pedal curve as shown subsequently. In summary, by moving the position

of the pedal point, it is possible to synthesize a variety of local deformations. Ti!, curve a(t)

will henceforth be referred to as the generator for the pedal curve 0 (t) and the process of

generating a pedal curve will be referred to as the pedaling operation.


2.2 Snake Pedals


1 .ire general shapes may be synthesized by letting the pedal point be different for each

point of the generator when we apply the pedaling operation. We can let the pedal points





















(a) (b)

Figure 3: (a) Tli process of generating a snake pedal with an ellipse as a generator. (b)
-i~1. pedal" controlled by the snake using an ellipse generator.


be specified by another curve p(t), which can be represented by a standard snake/B-snake

[17, 18] and then apply the pedaling operation to each point on the generator ac =a (ti)

with respect to corresponding pedal point pi = p(ti), where a and pi are points on the

generator a (t) and curve p(t) respectively. In this paper, we will focus on parameterized

generator curves. Ti, pedaling operation generates a new curve that we dub a snake pedal

x(t) as shown in Fig. 3.

If the generator is an ellipse as shown in Fig. 3 (a), we can represent it in a parametric

form by

L cosO sinO a cost m1 in 1
-sinO cosO b sint m2

where a, b are aspect ratio parameters, 0 is the rotation angle between the world coordinates

and the object centered coordinates, m = (ml, m2)' is the centroid of the generator in the

world coordinates. We collect the generator parameters into a vector g = (a, b, 0, mi, m2)'

and the vector g is referred to as the global shape parameter vector.

For an ellipse, the term Ja'(t) in the Eq. (2) of the pedal curve can be expressed as

Ja (t) n[ ]= -al sin t sin 0 a2 cOs t cOs 0 (4)
n2 --al sin t cos 0 + a2 cos t sin 0









Note that .[a'l = n + n = a, sin2t a cos2t, which is independent of the rotation angle

0. We will see that .*a'I is a very important quantity in the model fitting process.

Ti! pedal point can belong to a curve, which is not restricted to be a parameterized function.

We can sample the curve at the same rate as the sampling rate for the generator so that

there is a one-to-one correspondence between this curve and the generator. We use a snake

as this controlling curve and get a sequence of position vectors, p; = (xi, ;, )', by sampling

the snake. We collect all the snake position vectors in another vector, p = (p', p ,..., p)

where 1[ is the number of sampling points. Note that the snake may be replaced with a

B-snake to achieve compactness in the representation. Vector p is referred to as the local

shape parameter vector. A snake pedal curve is completely determined by the global and

local shape parameter vectors.

In Fig. 4, we depict some examples of snake pedals, which are generated using the snakes

and an ellipse as the generator. Note the variety of local deformations that can be generated

using this modeling technique. We remind the reader that the snake pedal itself is a geometric

model and that it is not directly responsive to the application of external forces unlike the

standard snake models [17]. Also, it is easy to synthesize snake pedals whose topology is

different from that of the generator as will be seen in examples depicted subsequently. TI

uniqueness of our modeling scheme lies in the ability to generate snake like local deformations

in a non reactive geometric shape namely the pedal curve/surface.

Tli, pedal curve definition can be modified slightly by subtracting the second term from

the first term in Eq. (2). This modification allows us to synthesize a large class of shapes

which appear very different from the shapes produced using Eq. (2). Til, key feature

of using this modification is that the snake pedal curve allows for more local deformations

including shrinkage and expansion of the snake pedal. Tli, shrinkage and expansion behavior

is controlled by the location of the snake. When the snake is located inside the generator,

the snake pedal exhibits shrinkage behavior, while expansion is exhibited when the snake is

outside the generator. Tli, following equation describes this modification:



















(a) (b) (c)

Figure 4: Examples of -i., ;I pel .l!-" using an ellipse generator and a snake.




(ca(t) p) J '(t)(5)
pedal(t) = p J )(tM. (5)

Fig. 5 (a)-(c) depict examples of shapes generated using this modification. Note the amount

of local deformations allowed in the snake pedal curves using this modification. We compare

the original and i i lio1i I" snake pedals using the same generators and same snakes in Fig.

5 (d)-(e).

Topology changes maybe realized by applying the pedaling operation to a single generator

and several pedal points (or snakes). Fig. 6 depicts an example of the topology change

realized in 2D. When the snake is located inside the generator, a single closed snake pedal

curve is generated, when the snake is split into two parts with both parts located outside

the generator, two separated closed snake pedal curves are obtained. In the later case, only

part of the generator is used to produce the snake pedal curves. Note that any number

of closed snake pedals can be generated from one generator using this procedure if the

position of snakes and parts of the generator are properly chosen. T11 -, quantities may be

automatically determined in data fitting applications using an appropriate technique.


















(b)


Figure 5: (a) (c) Examples of modified -nil.l;' pedals" using an ellipse generator and a
snake. (d)-(e) Comparisons of original(red) and modified (blue) -il.;h, pe l.i!- using the
same generators and same snakes in both cases.


J I

I
JI




(a)


I I
I I




0

(b)


I









(C)


Figure 6: Examples of topology change. (a) depicts a single snake pedal curve (solid)
obtained from an ellipse generator (dashed) and a single pedal point. (b) same as (a) except
two distinct pedal points located outside the generator are used. (c) same as (a) except
several pedal points on a snake are used. (d) same as (c) except a discontinuous snake is
used.


(C)









2.3 Pedal Surfaces


A pedal surface is the surface analog of the pedal curve. It is the locus of the points on the

foot of the perpendicular from a fixed pedal point to a variable tangent plane of the surface.

1, i'e precisely, it is defined as follows [16]:

Definition : T.,, pedal surface of a parameterized surface or a patch a : U -+ R3 with

respect to a point p E S3 is defined by


(a(u, v) p) (a, x a0,)
pedal(u, v) = p + ( )a~U x av, (6)
au X v

where U is an open subset of S2, aQ(u, v) is a parameterized generator surface with parameters

u and v, a, and av are tangent vectors in the u, v directions.

Fig. 7 depicts examples of pedal surfaces synthesized using an ellipsoidal generator surface

and a fixed pedal point located inside the generator surface in each example. Note that

change in topology is achievable by simply choosing an appropriate location of the pedal

point as shown in Fig. 7. However, this type of topology change is different from what we

described in Fig. 6.

As in the 2D case, we can let the pedal point vary for each point on the generator surface.

Tliil-. we have the snake pedal surfaces in 3D whose shape can be controlled by snakes which

are either curves or surfaces in 3D. For examples of the 3D snake pedals, we refer the reader

to the next section.

In comparison to the deformable superquadric model, the snake pedal model can exhibit a

larger amount of deformation and can easily exhibit a large amount of bending, twisting and

.,l, '.,I without the explicitly using parameters to represent these operations in the model

representation. Ti, i, fi 'e, it is well suited for modeling and recovery of shapes with complex

geometries, examples of which are Rabundant in the field of medical image analysis.











T1 Y


(a) (b) (c)




(d) (e) (f)

Figure 7: Examples of pedal surfaces with fixed pedal points located inside ellipsoidal gen-
erator surfaces.


3 Model Fitting

In this section we will develop an efficient technique for fitting our -~I.1Ii pedal" model to
a sparse set of data points in 2D/3D as well as to image data for shape recovery. Ti, model
fitting is posed as a nonlinear optimization problem and can be solved numerically. Til, 1,
are two types of parameters which need to be determined in this optimization process. One
is the global shape parameters namely, generator parameters g = (a, b, 0, mi, mi2)T (in 2D),
and another is the local shape parameters, i.e., snake position p = {p }1. A fast numerical
algorithm is proposed to estimate these two types of parameters. Ti, algorithm involves
employing the Levenber-:'. i.Irquardt (i_ -) technique to solve for the global parameters in
the outer loop and the Alternating Direction Implicit (ADI) method [19, 20, 21, 22] to solve
for the local parameters in the inner loop (for fixed global parameters) as shown in Fig. 8.
We present the details of each of these methods in the following sections.


OL











While ( Energy > e )


One iteration of

solving global parameters


Solve local parameters

completely





Figure 8: Big picture of numerical schemes for model fitting.

3.1 LM Method for Global Parameter Estimation


Ti, estimation of the global shape parameters in the outer loop is a nonlinear optimization

problem, and the Levenber,--.i, 1rquardt (1 .i1) method [23] is utilized for achieving the opti-

mization. T1i, I1-'. 1 method essentially combines the '; "ton method and the gradient descent

(GD) method in solving the nonlinear function minimization problem. When the current

solution is far from the true solution, the i -.1 method adopts the gradient descent (GD)

algorithm. On the other hand, when the current solution is close to the true solution, the

nonlinear function can be approximated by a quadratic form and the I -.1 method switches

to the '; "ton method in this case. For a more detailed description of this algorithm, we

refer the reader to [23].

We now adopt the I -.1 method to our model fitting problem. To fit the snake pedal to a

given data set of m points {Di}l, in 3D, we minimize the sum of the squared distances from

each data point Di to the corresponding closest point xi on the pedal curve/surface. Let h

be the function that returns the closest point on the snake pedal from a given data point








Di i.e., xi = hi(g, pi, ti). This can be achieved by employing the technique of migration of

the point of influence as in [5]. Let {pi}1, be a set of \[ discrete points along the snake

p(t), {xi}l'l be the corresponding points on the snake pedal curve. Let D = {D } and

p = {pf}I1, then the quantity to be minimized for the model fitting can be written as


1m
ED(D, g,p) = 3(D x)2, (7)
Si=1
where xi = hi(g, pi, ti) and the magnitude of 3 is related to the uncertainty of data mea-

surements. T!i, quantity ED(D, g, p) can be considered as the external energy imposed on

the model in the fitting process.

In 1 -.1 the gradient of ED(D, g, p) with respect to generator parameters g = (gi, g,..)T

(a, b, 0, mi, mn2)T has the following components:


OED Thi (g, pi, ti)
-- (Di hi(g, pi, ti) ). (8)

If we assume that there is a one-to-one correspondence between the points on the generator

and the pedal snake then the partial derivative OED can be written as:


OED ihk (, Pktk)
-E = --P(Dk hk(g, Pk, tk) ). (9)
(0 1. 0.,.

Taking an additional partial derivative gives xL which is needed in the I-.1 iteration.

For the derivative of the energy required in the i. algorithm, we use the chain rule to get
OEDo = (OED)(h ). -. v1ire specifically, if we denote each component hi of h as hi = (hil, hi2)T'

then


[Oh il Oh1 Oh1 Oh Oh 1
OED [ (Dil- hil) (Di2 Oa Ob 00 0ni1 0'
Oa Ob 90 Omn1 O,1,
(10)








After the estimation of the global parameters, we solve the local shape parameters in the

inner loop on the parameter estimation algorithm using the ADI method.


3.2 ADI Method for the Local Parameter Estimation

When solving for the local parameters in the inner loop for fixed global parameters, the

quantity to be minimized has a quadratic form and its minimization requires the solution to

a linear system, namely, Kx = b, where K is the stiffness matrix, x is the local parameter

vector (set x = p), and b = '" which is the external force described later. Because

of the special structure of the matrix K in this linear system, we can use the alternating

direction implicit(ADI) method for the local parameter estimation. We can prove that the

ADI method takes only O(N) time to converge with soft data constraints where N is the

number of nodes in the rectangular finite element discretization of the snake in our model.

Note that in soft data constraints, the model is not required to interpolate the data but is

expected to approximate it.

We now present the structure of the stiffness matrix K. In the snake pedal model, the

controlling snake p is a controlled continuity spline under the influence of image forces or

other external constraint forces. In our model, the spline deformation energy imposed on

the snake is a membrane energy plus thin plate energy given in a continuous form by


EO(p)+( OP+ a22 P2 a2 (+( a2(1
E(p) W ( )2 ( 2( )_ 2( ) du dv, 11)

where wi and .,_ are the tension and rigidity factors respectively.

T!, external energy used can be the same as the energy ED(D, g, p) discussed in Section

3.1. In this case, the external force imposed on the snake is given by


8ED 8ED 9hi 9h,
O- O E3[ (Dil hi1) (Di2 hi2) O, (12)
Opi Ohi Opi pi








/South Pole




Total number of nodes :

n Xn+2
V
Denote : N = nX n




North Pole
u

Figure 9: :'.i del tessellation in intrinsic coordinates.



where

Ohi n 2 n + n n + n* (
OhP P1 1 2h 1 2 2

1_ 1 2 2

with ni, n2 defined as before. Notice that OED/Ohi is the external force in a standard snake,

but in our case, the external force is OED/Opi, which is related to the force in a standard

snake via a Jacobian matrix 0hi/Op, in Eq. (13).

In order to fit the snake pedal curve to salient features in an image, we use the following

equation as the energy to be minimized:



ED = -(G, VI) -|(G, VI(x) -, (14)

where x = x(x(t),y(t)) describes the snake pedal curve, G, is a Gaussian with standard

derivation a, and denotes the convolution operator. '. .re sophiscated boundary finding

energies can be synthesized and we refer the reader to [5]. By minimizing the above energy,

the model is made to fit to points of high gradient in the image. Ti, procedure to minimize

above energy is similar to that of fitting of snake pedal curves/surfaces to data points.

We then discretize the model in model coordinates using finite elements. Geometrically, the








snake pedal surfaces are closed surfaces in space whose intrinsic coordinates u = (u, v) are

defined on a domain Q. Ti, positions of points on the model in an inertial frame of reference

are given by a vector-valued function pedal(u) = (Ul),(u), x 3(u)). Til tessellation of

u = (u, v) in intrinsic coordinates is shown in Fig. 9. We use a combination of quadrilateral

and triangular elements in the discretization. In the inner quadrilateral mesh, we have

(n x n) = N nodes and 2 additional nodes are used in representing the north and south pole
leading to a total of (N + 2) nodes. Due to the presence of the north and south poles, after

the discretization of the model in intrinsic coordinates u using finite elements, the stiffness

matrix K C R(N+2)(N+2) in our problem can be written as a (2 x 2) block matrix:



K K K 12 (15)
K K21 K22
where K11 C RNxN contains coefficients which corresponds to the inner quadrilateral mesh

and K22 C -2x2 is a diagonal block matrix containing the coefficients corresponding to the

north and south poles, and K12 = KT1 contains the coefficients of interaction between the

poles and other nodes.

To solve the linear system Kx = b for the local parameter estimation, we employ the

following Schur complement formula [22]


K-lb Kll + K1K12S-1K21K11 -K1 K12S-1
K _b K S1KK11 S- 1 S1 b, (16)
-S-1K21K-l S_
where

S = K22 K21K1K12. (17)

From Eq. (16), it is evident that to obtain the solution of K-lb, we need the solution of

several Kllbs. Further analysis shows that the component K11 of the stiffness matrix K

can be expressed as the sum of two Kronecker products:



K11 = wi(AB+B A)+ ,A A









= ADC+CA, (18)


where C = wiB + (.r_/2)A with A and B being (n x n) tridiagonal and symmetric positive

definite (SPD) matrices of the form



1 -1 2 1
-1 2 -1 1 4 1
A = .. .. .. ~~nx"; B = .. .. .. C Rnxn. (19)
-1 2 -1 1 4 1
-1 1 1 2

Tli, Kronecker product of C and A, denoted by C A, is the matrix of order N x N

(where N = n x n) and given by


c11A c12A ... clnA
c21A c22A ... C2nA
CA= "' (20)


CniA Cn2A ... c,, A

When computing with Kronecker products, it is computationally convenient for us to rep-

resent vectors using matrices. When solving (C A)x, we represent the vector x of order

n2 by the matrix X = {xkt} of order n x n defined by



Xkl = Xk+n(t-1). (21)

It can be shown that this representation can lead to efficient procedures for computing

(C A)x and solving (C A)x = b, respectively. Tli, following lemma is used for this

purpose.


LEMMA 1

Let C = {ckl}, A = {akl} and X = {xki} be matrices of order n x n, then the n x n matrix

representation of n2 x 1 vector (C A)x is given by

19








(C (D A)x (C(AX)T)T.


Using this representation, the linear system KZ1x


b can be written as:


CXAT + AXCT = D,


(22)


where D, X are the n x n representation of the n2 vector b and x. Since A and C are also

SPD matrices, Eq. (22) can simply be written as


CXA + AXC = D.


(23)


After some simple linear algebra operations, it follows that above equation is equivalent to


A'X + XA' = D',


(24)


with A' = A-1C and D' = A-1DA-1

Eq. (24) is the standard Lyapunov matrix equation [24]. Ti1, ADI method can be applied

to solve this matrix equation and has the following steps in each iteration j = 1, 2,..., J,


(A' + pjI)Xj-i

Xj(A' + pjI)


(25)

(26)


D' Xj_i(A' pjI)

D'- (A'- pjI)Xj-


For each iteration in ADI, we need to solve two equations, namely, Eq. (25) and (26), and

therefore Xj. is introduced as an intermediate variable. Since A' is a SPD matrix, we
j2
can use Cl .1 -l.y decomposition for the tridiagonal SPD matrix (A' + pjI) to compute the

update solution. This only takes O(N) time per iteration. Tli, parameters pj are called

ADI parameters and J is the number of iterations needed. T!i ~ are chosen according to the

method described in [25].








4 Model Fitting Results


In this section, we present a set of five experiments, two of these experiments are in the

medical imaging domain and demonstrate model fitting to sparse 3D data points placed

by an expert neuroscientist along the boundaries of a hippocampus and a -- i-. in selected

slices of an .i: brain scan. In both the experiments, the spatial resolution of the images

was (256,256) and a region of interest was identified interactively. Sparse set of points are

placed on the boundary of the shape of interest because these shapes of interest do not

have sufficient gradients along large segments of the boundary. This requires that prior

anatomical knowledge be introduced into the model fitting process. What we describe here

in this paper is the stage of the fitting process requiring prior knowledge from an expert

in the form of a sparse set of boundary points on the shape of interest. This is used in

conjunction with image gradients as external forces in the model fitting process. Once the

model is successfully fitted, its parameters over a large number of fits can be learnt for later

use. Tli, third example depicts the snake pedal fitting in 3D to a Vulcun (Spock from Star

Wars) head. And the last two experiments demonstrate the power of the model in fitting

to 3D synthetic data from a bent shape (hair pin) and a twisted shape respectively. Once

again, we would like to stress the fact that the model does not make use of any explicit

1, iiliii- /twisting transformations to recover the underlying shapes from these data.

Tli, medical data example in Fig. 10 is organized as follows: left to right, the images depict

a slice of an i, i: brain scan in which the shape of interest -hippocampus or a -- i i-. has

been identified by a neuroscientist via sparsely placed points (only in selected slices of the

.1i: scan) on the shape boundary. TI, next image shows the collection of these 3D points

(in red) and the initialized snake pedal model (in green) followed by images depicting the

intermediate stage of fitting and the final fitted model respectively. Tli, model fitting process

takes 258s and 291s CPU time on SGI-O2 respectively. As evident, the model achieves a

visually accurate fit to the data and the model initialization is quite far from the final true















(a) (b) )()


(e) (f)


Figure 10: Firl in results of snake pedal surfaces to a hippocampus (d) and a.-- 1 i-. (h) with
initializations shown in (b) and (f), the intermediate stages are shown in (c) and (g). (a)
and (e) are slices of, 1I i1 scans for the hippocampus and -- i ii-

shape in each of these examples. In addition, we have verified that these fits are in agreement

with xwhat the expert neuroscientist would achieve via manual segmentation.

Ti, non medical data examples in Fig. 11 are organized as follows: left to right, the

images depict a 3D synthetic data set (in red) along with the model initialization (in green)

superimposed which is followed by images of intermediate stage of fitting and the final model

fit. Tlii CPU time for the model to converge to the desired shape -spock, bent shape and

twisted shape -is 267s, ;l..-. and '"','- on SGI-02 respectively. Note that in the synthetic

data examples, in one case, the data exhibits a large amount of bending and in another, a

large amount of twist. What we would like to stress here, is the fact that the model did not

require any explicit bending or twisting transformations (unlike the deformable superquadric)

to achieve the fits. As mentioned earlier, modeling schemes that require explicit incorporation

of bending, twisting transforms are highly nonlinear and are known to introduce numerical

instabilities into the fitting algorithms. Ti, fitting times are quite large in these examples


(c) (d)








than in the previous two because of the large amount of bending and twisting present in the

desired shapes and also because the initialization was quite far from the final solution.

We performed a comparison between our model fitting and the deformable superquadric

model fits to both synthetic and medical data using a fixed CPU time (that required in

achieving the fits with our model to the .- il-. and the twisted shape data shown in Figs.

10 (e)-(f) and 11 (g) respectively). For these fixed CPU times, the deformable superquadric

fits did not converge to the desired solution as is evident from Fig. 12 (b)-(e). Figs. 12

(b)-(c) depict model fitting to data shown in (a) using the deformable superquadric and the

snake pedal models respectively. A similar arrangement is depicted in Fig. 12 (e)-(f) for

data shown in (d). In both deformable superquadric model and snake pedal model fitting,

the numerical techniques used were the same. As evident from these comparisons, the snake

pedal model was able to achieve visually accurate fits while the deformable superquadric

yielded results that are far from satisfactory (within the fixed CPU time).


5 Conclusion


In this paper, we introduced a novel geometric shape modeling scheme which allows for

representation of global and local shape characteristics of an object. Geometric models are

traditionally well suited for representing global shapes but not the local details. However, we

proposed a powerful geometric shape modeling scheme which allows for the representation of

global shapes with local detail and permits model shaping as well as topological changes via

physics-based control. Our modeling scheme consists of representing shapes by pedal curves

and surfaces. Physics-based control was introduced for shaping these geometric models by

using a dynamic spline to represent the position of the pedal point. Ti, model dubbed as a

-i~i.1: pedal" allows for interactive manipulation via forces applied to the snake.

We also developed a novel technique for fitting our model to a sparse set of data points in

2D/3D as well as to image data for shape recovery. Tli, model fitting is posed as a nonlinear































(d) (e)


(g) (h) (i)

Figure 11: Filling results of snake pedal surfaces to (c) spock, (f) a bent shape and (i)
a twisted shape. Initializations are shown in (a), (d) and (g) respectively followed by the
intermediate stages in (b), (e) and (h).
























(d) (e)


Figure 12: Comparison of deformable superquadric and snake pedal model fitting perfor-
mance. (a) and (d) are model initializations, (b) and (e) exhibit model fits with deformable
superquadric models, (c) and (f) are model fits with snake pedal models. CPU time for fits
in (b) & (c) and (e) & (f) respectively are the same.


optimization problem and can be solved numerically. Tli, nonlinear optimization scheme

developed in this paper contains the Levenbers- i .iArquardt (i -.) method in the outer loop

and the Alternating Direction Implicit (ADI) method in the inner loop, and the combination

of the global and local scheme is quite stable and yields a locally optimal solution in very

few iterations.

We demonstrated the applicability of our modeling scheme and numerical algorithms via

model fitting to real image data and range data. Future work will focus introducing dynamics

into the model and experimenting with more image data.


6 Appendix


In Section 3.2, we described the internal deformation energy in an active/deformable model.

To obtain the expression for the stiffness matrix K, we employ the finite element method

(FEM) to discretize the model in material coordinates u = (u, v). As shown in Fig. 9, we










v I- a "1

S(-1,1) (1,1)
W1 T

(u v,) )

(-1,-1) 2 (1,-1)


u


Figure 13: Bilinear quadrilateral element, where (u, v) are material coordinates and (w, y)
are the reference coordinates.

use bilinear quadrilateral elements in conjunction with linear triangular elements. TilI later
is used for representing the regions around the north and south poles. We will first discuss
the bilinear quadrilateral elements, and the linear triangular elements can be handled in the
similar way.


6.1 Bilinear Quadrilateral Elements

Bilinear quadrilateral elements are quite common in FEM literature [26]. In this section,

we present the basic transformation equations between material and element coordinates.
Denote the finite element nodal shape functions by Ni where i = 1, 2, ..., no, and no is the
number of nodes associated with each element Ej. For the bilinear quadrilateral elements,
we have

N(w,) = (1 + ww) (1 + ), (27)
4
where (wi, .' ) are the reference coordinates of node i as shown in Fig. (13).

T11, reference coordinates (w, ) are related to the material coordinates u = (u, v) by

2 2
w = (u Uc) = (v ), (28)
a b








where (uc, v,) are the coordinates of the element center, a and b are the length and width of

the element. In our implementation, we use a = b. Differentiating Eq. (27) with respect to

u and v gives


ON8 ONZD w N1 8
Du Dw Du 9 D u
-ONi Ni D w +Ni 9D =
Ov Ow Ov 9 Dyv
When integrating function f(u, v) over element Ej,

ence system:


I IE f(u,v)dudv


1
2aW(1 + )
2a


(29)


(1 + wwi (30)
2v

we perform the integration in the refer-


.4


(31)


6.2 Derivation of the Local Stiffness Matrix


As described in section 3.2, we use the membrane plus thin plate energy expressed in Eq.(11)

as the internal deformation energy imposed on the snake,

EO !+ Dp+ p2 2 + D a 2p
E,(p)= wI(( )2( )2 (( + 2( )2 2 ) du dv. (32)
Ou OV DuOv

By using the procedure described in [5], we can obtain the local stiffness matrix. For the

membrane energy part, if we denote the stiffness matrix for the element Ej by K c SR xnoo,

the entries of Km can be represented as

/ DN1,DNm D8NN D,
(k'I)wm (= Ni ( + D D ) du dv, (33)
M BE u Ov Ou 0V

where m = 1,2,...,no. Since we use a quadrilateral element, no = 4 in this case. Using

the formula in Eq. (29), (30) and (31), we obtain the following local stiffness matrix for the

membrane energy


K = (w )
m 6 a


1

-2 6 bL
2


-2
1
1
2


(34)























V



Figure 14: A small 3 x 4 mesh in (u, v) space.


Adopting the similar scheme, we obtain the following stiffness matrix Kj for the thin plate

energy part in element E,


1 1
KT g s
t ab 1
S-1

T!l, summation of Kj, and Kj gives the full


-1 1 -1
1 -1 1
-1 1 -1
1 -1 1

stiffness matrix for element Ei


(35)


(36)


K' = K J+ Ki.


Note that both KN and KI are symmetric positive definite (SPD) matrixes, therefore their

summation, Kj, is also a SPD matrix.


6.3 Kronecker Tensor Product Representation of the Assembled
Stiffness Matrix


Ti, matrix Kj discussed earlier is only the local stiffness matrix for element Ej, which

describes the stress-strain relationship between all the nodes within element Ej. For

example, k2 = 4 means applying one unit displayment to node 3 results in 4 unit force on

node 2. We need to assemble the local matrix over all the elements to obtain the entire


E1 E2 E3
5 6 7

E4 E5 E6
9 10 11








wib
6a
0


wib
-x
6a
0


0

02


02


---x
wib
6a
00


wib
-x
6a
0


0


Figure 15: '.1I iii I oI iii molecules
@--@--@


w1b 4x8 -
x a--


@--@--@


Figure 16: Assembled membrane molecular representation.

stiffness matrix K. To achieve this, we introduce a novel method via the Kronecker tensor
product operation. A small 3 x 4 mesh in (u, v) space as shown in Fig. 14 is used to illustrate
the idea. We will only consider the first part of the KJ (the first term in the summation of
Eq. (34)) for explanation, the second part of KJ and Kj can be handled in the similar way.

In Fig. (14), we assign a number to each node, ranging from 1 to 12, and denote each
element by El, E2, ..., E. We choose node 6 for our illustration, and investigate the stress
-strain realtionship between node 6 and every other node in the entire mesh. From Eq.
(34), observe that applying one unit displacement to node 6 will result in W (b) x 1 unit,
(w )b x (-1) unit, l(b) x (-2) unit force on node 2, node 1, and node 5, respectively. This








concept can be visualized by the nodal molecule representation, as shown in Fig. (15) (a). We

can use a similar representation for elements E2, E4 and r -, resulting in the nodal molecules

as shown in Figs. (15) (b), (c), and (d), respectively. Note that there is no interaction

between node 6 and nodes in elements E3 and Eg. From the nodal molecular representation,

when applying one unit displacement to node 6, we can obtain the force resulting on every

other node in the entire 3 x 4 mesh by summing up all the molecular values at the same

molecular location. Tli, result is shown in Fig. (16).

Ti, assembled molecular representation can be expressed as a Kronecker product formula


(- [ -1 2 -1 ] 4 (37)
6a I

Ti, assembled stiffness matrix that corresponds to the first part of the membrane energy

(first term in the summation of Eq. (34)) is



Karti) = AX By, (38)

where
2 1 1 2 -1 0 -1 -
-1 2 -1 0
A = 1 4 1 and By 1 (39)
0 1 2 -1 0 -1 2

We now prove the above conclusion by expanding the Kronecker tensor product expression

in (38). According to the Kronecker product definition in Eq. (20),











2 1 2 -1 0 -1
K(partl) ( b) 4 1 2 -1 0
m 6 a 1 0 -1 2 -1
0 2 -1 0 -1 2
2By By 0
a= ( ) By 4By By
0 By 2By

-4 -2 0 -2 2 -1 0 -1
-2 4 -2 0 -1 2 -1 0
0 -2 4 -2 0 -1 2 -1
-2 0 -2 4 -1 0 -1 2
2 -1 0 -1 8 -4 0 -4 2 -1 0 -1
(wb -1 2 -1 0 -4 8 -4 0 -1 2 -1 0
6a 0 -1 2 -1 0 -4 8 -4 0 -1 2 -1
-1 0 -1 2 -4 0 -4 8 -1 0 -1 2
2 -1 0 -1 4 -2 0 -2
-1 2 -1 0 -2 4 -2 0
0 -1 2 -1 0 -2 4 -2
-1 0 -1 2 -2 0 -2 4
(40)

In this matrix, each entry (Y ,nortl)),l represents the stress-strain relation between node I

and node rn, where I, r = 1, 2,..., 12. For example, the entries in the 6th column (. (n'rt1))16

represent the force caused in node I by the unit displacement in node 6. IWe can see from

column 6 that one unit displacement in node 6 contributes (w1 b) x 8 unit force to itself, and
6 a
(w1 b) x (-4) unit force to node 5 and 7, (wV ) x 2 unit force to node 2 and 10, (w 6) x (-1)
S6 a/ v \/ "' 6 a' '" 6 a '
unit force to node 1, 3, 5 and 7, and zero force to node 4, 8 and 12. This is exactly what

Fig. (16) conveys.

For the second part of the membrane energy, performing a similar procedure, we have


KI = (6w ) B, A, (41)

where
0 4 1 0 1 -
1 4 1 0
A,= -1 2 -1 and B= 0 1 4 1 ()
-1 1 1 0 1 4









Since a = b in our implementation, the entire stiffness matrix corresponding to the membrane

energy is


K = K (arti) + K..

= (W) (A, 0 By + B 0 Ay).
Similarly, we have the stiffness matrix corresponding to the thin plate energy


K1t ( ) A, (D Ay.
ab


In general, for an m x n mesh, denote

2 1
1 4 1
A X ... .. ...
1 4
1

4 1
1 4 1
Ay 1 .. . ..
1 4
1 1


(43)


(44)


(45)


E Smxm
1
2

1


1
4


(46)


1 -1
-1 2 -1
Bx ... ... ... R SMxm (
-1 2 -1
-1 1

2 -1 -1
-1 2 -1
By = ... ... ... (
-1 2 -1
-1 -1 2

then the assembled stiffness matrix corresponding to the energy expressed in Eq.(11) is:


K Km + Kt = wl(Ax 0 By + Bx Ay) + ,,'_(A ( A,), (


47)





48)


49)


where w1 and .,_ are tension and rigidity parameters. Note that the difference between the

first column/row and last column/row in A, and Ay or B, and By is due to the different

boundary conditions namely, the natural boundary and the periodic boundary conditions








respectively. When employing the same boundary conditions and let m = n, A = Ax = Ay,

B = Bx = By, Eq. (49) becomes


K = wi(A B + B A) + ,,'_(A A). (50)


This is the same equation as (18).


References


[1] E. Bardinet, N. Ayache, and L.D. Cohen, "Fitting of iso-surface using superquadrics

and free-form deformations," in Proc. of IEEE Workshop on Biomedical Image Analysis

WBIA '94, Seattle, 1994, pp. 184-193.

[2] I. Cohen and L.D. Cohen, "Hyperquadric model for 2d and 3d data fitting," in Proc.

of ICPR, Piscataway, NJ, 1994, pp. 403-405.

[3] D. DeCarlo and D. :'. i1 .~i\- "Blended deformable models," IEEE Trans. Pattern

Analysis and Machine Intelligence, vol. 18, pp. 443-448, 1996.

[4] A. Pentland and J. Williams, "Good vibrations: :., del dynamics for graphics and

animation," Computer Graphics, vol. 23, no. 3, pp. 215-222, 1989.

[5] D. Terzopoulos and D. '.1 ii, \.-, "Dynamic 3d models with local and global deforma-

tions : Deformable superquadrics," IEEE Pattern Analysis and Machine Intelligence,

vol. 13, no. 7, pp. 703-714, 1991.

[6] '. Hebert, K. Ikeuchi, and H. Delingette, "A spherical representation for recognition

of free-form surfaces," IEEE Pattern Analysis and Machine Intelligence, vol. 17, no. 7,

pp. 681-690, 1995.

[7] W.:'.1 Hsu, "Direct manipulation of free-form deformation," Computer Graphics, vol.

26, no. 2, pp. 177-184, 1992.









[8] T.W. Sederberg and S.R. Parry, i I -t I,, deformation of solid geometric models," in

ACM SIGGRAPH, Dallas, TX, 1986, pp. 151-160.

[9] B.C. Vemuri and A. Radisavljevic, il liresolution stochastic hybrid shape models

with fractal priors," ACM Trans. on Graphics, pp. 177-207, 1994.

[10] S. Han, D.B. Goldgof, and K.W. Bowyer, "Using hyperquadrics for shape recovery from

range data," in IEEE Proceedings of the 3rd International Conference on Computer

Vision, Berlin, 1993, pp. 492-496.

[11] T. O'Donnel, T. Boult, and A. Gupta, "Global models with parametric offsets as applied

to cardiac," in IEEE Proc. of CVPR, 1996, pp. 293-299.

[12] V. Caselles, F. Catte, T. Coll, and F. Dibos, "A geometric model for active contours in

image processing," A,..o, cische Mathematik, vol. 66, pp. 1-31, 1993.

[13] R. i.illili J.A. Sethian, and B.C. Vemuri, 'lIlpe modeling with front propagation

: A level set approach," IEEE Trans. Pattern Analysis and Machine Intelligence, vol.

17, no. 2, pp. 158-175, 1995.

[14] T. '.1I ii,, i,,y and D. Terzopoulos, "Topologically adaptable snakes," in IEEE ICCV,

1995, pp. 840-845.

[15] L.D. Cohen, "Auxiliary variables and two-step iterative algorithms in computer vision

problems," Journal of Mathematical Imaging and Vision, vol. 6, pp. 59-83, 1996.

[16] A. Gray, Modern differential geometry of curves and surfaces, CRC Press, Boca Raton,

1993.

[17] :. Kass, A. Witkin, and D. Terzopoulos, "Snakes : A active contour models," Int.

Journal of Computer Vision, vol. 1, pp. 321-331, 1987.

[18] P.S. i-.i,'c, S. i.i, 1 and G. -.i, ii1,ii "B-snakes: implementation and application to

stereo," in IU Workshop, 1990, pp. 720-726.









[19] D. Calvetti and L. Reichel, "Application of adi iterative methods to the restoration of

noisy images," SIAM J. Matrix Anal. Appl., vol. 17, no. 1, pp. 165-186, 1996.

[20] W.R. Dyksen, "Tensor product generalized adi methods for separable elliptic problems,"

Siam J. i~,.,,. c. Anal., vol. 24, no. 1, pp. 59-76, 1987.

[21] R.E. Lynch, J.R. Rich, and D.H. T!,ii_,,-1 "Tensor product analysis of alternating

direction implicit methods," J. Soc. Indust. Appl. Math, vol. 13, no. 4, pp. 995-1007,

1965.

[22] D.V. Ouellette, "Schur complements and statistics," Linear Algebra and its Applica-

tions, vol. 36, pp. 187-2'1i., 1981.

[23] W.H. Press, S.S. Teukolsky, W.T., Vetterling, and B.P Flannery, a',,,, rical Recipes in

C, Cambridge University Press, 1992.

[24] A. Lu and E.L. Wachspress, "Solution of lyapunov equations by alternating direction

implicit iteration," Computers Math. Applic., vol. 21, no. 9, pp. 43-58, 1991.

[25] S. H. Lai and B. C. Vemuri, "Physically-based adaptive preconditioning for early vision,"

IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 19, no. 6, pp.

.Q"I ,1 7, 1997.


[26] O.C. Zienkiewicz and R.L. Taylor, T7. Finite Element Method, :'. Graw-Hill, V:- -York,

NY, 1989.




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs