NONPARAMETRIC SHAPE PRIORS FOR ACTIVE CONTOUR-BASED

IMAGE SEGMENTATION

Junmo Kim, Müjdat Çetin, and Alan S. Willsky

Laboratory for Information and Decision Systems
Massachusetts Institute of Technology
77 Massachusetts Ave., Cambridge, MA 02139, USA

ABSTRACT techniques have been applied to segmentation problems in-

When segmenting images of low quality or with missing volving low SNR or occluded images successfully, especially
data, statistical prior information about the shapes of the ob- when the shape variability is small. However, there are two
jects to be segmented can significantly aid the segmentation major shortcomings of such techniques. First, these meth-
process. However, defining probability densities in the space ods treat the signed distance functions as elements of a lin-
of shapes is an open and challenging problem. In this paper, ear vector space, and perform operations such as averaging.
we propose a nonparametric shape prior model for image Yet, the space of signed distance functions is a nonlinear
segmentation problems. In particular, given example train- manifold and is not closed under linear operations. For ex-
ing shapes, we estimate the underlying shape distribution by ample, the average of signed distance functions, which is
extending a Parzen density estimator to the space of shapes. commonly used to obtain a mean shape, is not necessarily a
Such density estimates are expressed in terms of distances signed distance function. Therefore, the use of linear analy-
between shapes, and we propose two distance metrics that sis tools such as PCA gives rise to an inconsistent framework
could be used in this framework. We then incorporate the for shape modeling [3]. Second, these technique can handle
learned shape prior distribution into a maximum a posteriori only unimodal, Gaussian-like shape densities. In particular,
estimation framework for segmentation. This results in an these methods cannot deal with “multimodal” shape densi-
optimization problem, which we solve using active contours. ties, which involve multiple classes of shapes (e.g. a shape
We demonstrate the effectiveness of the resulting algorithm density of handwritten digits, composed of multiple digits).
in segmenting images that involve low-quality data and oc- In our work, we propose a framework for constructing
clusions. The proposed framework is especially powerful in nonparametric shape densities from example training shapes.
handling “multimodal” shape densities, involving multiple In particular, we assume that the training shapes are drawn
classes of objects. from an unknown shape distribution, and we estimate the
underlying shape distribution by extending a Parzen density
1. INTRODUCTION estimator to the space of shapes. Such density estimates are
expressed in terms of distances between shapes. We propose
We consider image segmentation problems that involve lim- two specific distance metrics to be used for nonparametric
ited and low-quality data. Such segmentation problems are density estimation, although other metrics could be used in
ill-posed and require the incorporation of prior information our framework as well. We then formulate the segmenta-
about the objects to be segmented. When human experts tion problem as maximum a posteriori (MAP) estimation,
segment images, they clearly make use of such prior infor- in which we use the learned nonparametric shape density
mation. For example, a radiologist can usually manually as the prior. This leads to an optimization problem for the
segment an organ (e.g. the prostate) in a magnetic reso- segmenting curve, for which we develop and use an active
nance image, although the boundaries are mostly invisible contour-based iterative algorithm. We present experimental
to a layperson. Clearly, radiologists augment the observed results of segmenting low-quality and occluded images. We
data with their expertise, in other words with statistical prior also demonstrate how the proposed algorithm can success-
information, about the shape of the organ. Existing auto- fully incorporate shape densities involving multiple object
matic segmentation methods (either explicitly or implicitly) classes.
enforce only very simple constraints about the underlying
shapes. For example, many active contour-based methods
(which is the framework we also use in our work) involve 2. NONPARAMETRIC SHAPE DENSITIES
a curve length penalty, which translates to the assumption
In this section, we address the problem of estimating an un-
that shorter curves are statistically more likely than longer
known shape probability density. Given n example train-
ones. However, in many applications, more structured prior
ing curves C1 , . . . , Cn , we first align1 them to obtain the
information about the shapes is available. Yet the challenge
is how to construct probabilistic descriptions in the space aligned training set C̃1 , . . . , C̃n . Ideally, this operation re-
of shapes, and then incorporate such statistical information moves the pose variability in the training data, and what
into the segmentation process. remains is just shape variability. The problem then is to
Early work on this problem involved landmark-based rep- estimate pC̃ (C̃) from which the training samples C̃1 , . . . , C̃n
resentations of shapes, and the construction of typical shapes are drawn. This density is a probability density over an in-
and typical variability based on a set of training shapes via finite dimensional space. We would like to leave the shape
principal component analysis (PCA) [1] . The use of land- of this density unconstrained, therefore we adopt a nonpara-
marks has the drawback that the performance of shape anal- metric density estimation route. Assuming that we are given
ysis depends on the quality of those landmarks. Recently, a distance metric dC (·, ·) in the space of curves C, we can form
there has been increasing interest in using level set-based rep- a Parzen density estimator as follows:
resentations for shape priors [2, 3], which avoid landmarks.
In [2] and [3], PCA of the signed distance functions of train- 1 We use the alignment algorithm of [3], but other techniques

ing data is used to capture the variability of shapes. These can be used as well.
 Tφ̃ D

1X
n
p̂C̃ (C̃) = k(dC (C̃, C̃i ), σ) (1) D
n i=1
||φ̃ − φ̃j ||L2
φ˜j φ˜i
where k(·, σ) denotes a Gaussian kernel with kernel size σ,
1 x2
i.e. k(x, σ) = √2πσ 2
exp(− 2σ 2 ).

Conceptually, the nonparametric density estimate in (1)

can be used with a variety of distance metrics, yet the key
issue is what kind of metrics make sense. This is related
φ̃
to the question of what it means for two shapes to be simi-
lar. In the following sections, we consider two specific met-
rics, namely the template metric2 and the Euclidean (or L2 )
distance between signed distance functions. The template Figure 1: Notional illustration of the space of signed distance
metric is given by the area of the set-symmetric difference functions D, in a scenario where example shapes form two
between interior regions of two shapes, and can be expressed clusters, corresponding to two classes of shapes.
as a norm of difference between two binary maps represent- However, computing a geodesic distance in an infinite di-
ing the shapes. On the other hand, the L2 metric we use is a mensional manifold is a challenging problem. There is some
norm of the difference between two signed distance functions. previous work on computing geodesic distances in the space
From a practical standpoint, the key difference between these of curves [7, 8], but there is little work when the shape is
two metrics is that the template metric puts equal weight represented by signed distance functions.
on pixels, whereas the L2 distance between signed distance Instead, we now consider the Parzen density estimate
functions puts variable weight. with the L2 distance in L:
1X
n
2.1 Template Metric
p̂φ̃ (φ̃) = k(dL2 (φ̃, φ̃i ), σ). (6)
We now consider the Parzen density estimate in (1) with n i=1
a specific metric, namely the template metric dT (C̃, C̃i ) =
Area(Rinside C̃
Rinside C̃i ) [5], where
denotes set sym- Now we discuss how the density estimate in (6) might be a
metric difference. The density estimate with the template good approximation to the one with geodesic distance in (5).
metric is given by If the example shapes are of small variation then the part of
the manifold supporting the example shapes can be approx-
1X
n
p̂C̃ (C̃) = k(dT (C̃, C̃i ), σ). (2) imately flat provided that the manifold does not have too
n i=1 much curvature. This is why methods involving linear anal-
ysis tools such as PCA of signed distance functions [2, 3]
In [6], we also consider a variant of this metric, namely its work reasonably well in the case of small shape variation.
square-root, as an alternative distance measure. Due to the same phenomenon, the L2 distance used in the
Parzen density estimate can be close to the geodesic distance.
2.2 Euclidean Distance between Signed Distance Now consider the case where the example shapes belong to
Functions multiple classes and form multiple clusters as illustrated in
Fig. 1. In this case, the part of the manifold supporting the
We now consider representing each shape C̃i , i ∈ {1, ..., n}, samples is no longer flat, and PCA would no longer be a
by its corresponding signed distance function φ̃i , and con- valid approach. In contrast, the density estimate with L2
structing a shape density estimate p̂C̃ (C̃) = p̂φ̃ (φ̃).3 We ob- distance can still be a good approximation of (5) for the fol-
serve that the space of signed distance functions D is a subset lowing reasons. For two shapes that are in the same cluster
of an infinite dimensional Hilbert space L
{φ|φ : Ω → R}, (e.g. φ̃ and φ̃i in Fig. 1), hence are of small variation, the L2
where Ω denotes the image domain. We can define an inner norm will be a good approximation of the geodesic distance.
product and an induced L2 distance in L as follows: On the other hand, for shapes that are in different clusters
Z (e.g. φ̃ and φ̃j in Fig. 1), there will be an error in approxima-
φ1 , φ2  = φ1 (x)φ2 (x)dx, (3) tion of the geodesic distance, but the overall error in density
Ω estimation will be small as long as the kernel size σ is small
p
dL2 (φ1 , φ2 ) = φ1 − φ2 , φ1 − φ2 . (4) compared to the distance dL2 (φ̃, φ̃j ).

Since the space D is embedded in a Hilbert space, a natural 3. SHAPE-BASED SEGMENTATION

metric d(φ1 , φ2 ) for this space would be a minimum geodesic Now we combine the nonparametric shape prior with a data
distance, i.e. the distance of the shortest path from φ1 to φ2 term, and formulate the segmentation problem as MAP es-
lying in D. If one could compute such minimum geodesic timation. This leads to the following energy functional to be
distances dgeodesic (·, ·), the corresponding Parzen density minimized for segmentation4 :
estimate based on aligned training samples {φ̃i } would be E(C) = − log p(data|C) − log pC (C). (7)
1X
n
p̂φ̃ (φ̃) = k(dgeodesic (φ̃, φ̃i ), σ). (5) The data term we use here is based on [9], and the second
n i=1 term is based on the nonparametric shape priors introduced
in Section 2. We minimize this functional iteratively by gra-
2 We note that the work in [4] (which came to our attention
dient descent, for which we need the gradient flow for the
while preparing this manuscript) has commonalities with our tem- curve C or the corresponding signed distance function φ.
plate metric-based approach.
3 φ and φ̃ denote the signed distance functions before and after 4 In the rest of the paper, we drop the hat for simplicity in

alignment, respectively. density estimate p̂.

1. Evolve C without the shape prior for time t ∈ [0, t0 ]
2. For the curve C|t=t0 , compute the pose p|t=t0 by aligning
C|t=t0 with respect to {C̃i }
3. Iterate until convergence:
(a) fix p and
i. compute C̃ = T [p]C
(a) (b) (c)
ii. compute ∂∂tC̃ for the shape prior log pC̃ (C̃)
iii. compute ∂C
∂t
from ∂∂tC̃ by ∂C
∂t
= T −1 [p] ∂∂tC̃ Figure 2: Segmentation of an occluded aircraft image. (a)
(b) update C by both the data and the shape force Test image. (b) Result without a shape prior. (c) Result
with the nonparametric shape prior.
(c) fix C and update p using an alignment technique
the velocity field in (9) that directly impact the shape evolu-
Algorithm 1: Iterative algorithm for updating the pose p
and the curve C. tion are those defined at the points on the boundary C̃, i.e.
where the φ̃ = 0. So we do not modify those components.
The gradient flow for the data term is computed as is done The key property we exploit is the following: if we start
in [9]. So we only focus on the gradient flow ∂C ∂t
or ∂φ
∂t
for the with a signed distance function and evolve it with a velocity
shape term log pC (C), where t denotes an artificial iteration field that is constant along the direction normal to the cor-
time. responding curve C̃, the evolving level set function remains
Note that any given curve C is not necessarily aligned a signed distance function [11]. So we construct a modified
with the training data, whereas the shape densities we pro- velocity field by copying the velocities on the boundary from
posed in Section 2 are based on aligned curves C̃ and {C̃i }. (9), and then extending these values in the direction normal
As a result, in order to compute ∂C ∂t
, we first compute ∂∂tC̃ , to the boundary.
∂C ∂ C̃
and then compute ∂t from ∂t . To this end, let C̃ = T [p]C,
where T [p] is a similarity transformation with the pose pa- 4. EXPERIMENTAL RESULTS
rameter p capturing translation, rotation, and scale. We Now we present experimental results demonstrating our seg-
update C and p iteratively as described in Algorithm 1. We mentation method based on nonparametric shape priors. We
discuss Step 3-(a)-ii below for the two distance metrics we first show results for segmenting occluded objects. Here the
consider. shape prior involves a single class of object shapes. Next, we
When we use the template metric-based shape prior present experimental results on a more challenging problem,
p̂C̃ (C̃) of (2) in (7), we obtain the following gradient flow [6]: where the prior involves multiple classes of object shapes.
In particular, we consider the problem of segmenting hand-
∂ C̃ 1 1 1 X  written digits (with low quality or missing data), where the
= k(dT (C̃, C̃i ), σ)dT (C̃, C̃i )(1 − 2H(φi ))N (8)
∂t pC̃ (C̃) n σ 2 i prior density involves all digits, and the algorithm does not
know which digit class the test data belong to. For the sake
where N  is the outward normal direction with respect to of brevity we present results with only the L2 metric; the
the curve, and H(·) is the Heaviside function, i.e. H(φ)=1 results with the template metric are similar.
if φ ≥ 0 and H(φ) = 0 if φ < 0. Note that the gradient
flow is composed of a weighted average of several directions, 4.1 Segmentation of Occluded Objects
where the i-th direction is an optimal (gradient) direction
In this section, we demonstrate our shape-based segmen-
that decreases the distance between the i-th training shape
tation algorithm on synthetic aircraft images. As example
and the evolving shape. We implement (8) using level set
shapes of this class, we have a set of 11 binary images (not
methods.
shown here), whose boundaries provide the training curves
Next, we consider the Euclidean distance-based shape
C1 , . . . , Cn (n = 11). We present segmentation results on
prior pφ̃ (φ̃) of (6). Let us first consider evolving φ̃ in the the image of an aircraft whose shape is not included in the
space L along the gradient of log pφ̃ (φ̃) w.r.t. φ̃. The result- training set. Fig. 2(a) shows the noisy aircraft test image
ing gradient flow would be [6]: with an occluded left-wing. Fig. 2(b) shows the segmenta-
tion of this image without a shape prior (the shape prior is
∂ φ̃ 1 1 1X
= k(dL2 (φ̃, φ̃i ), σ)(φ̃i − φ̃). (9) replaced by a curve length penalty term). In Fig. 2(c) we
∂t pφ̃ (φ̃) σ 2 n i show the result of our proposed technique, where the shape
prior information allows us to recover the wing of the air-
This flow involves a weighted average of {φ̃i − φ̃}n i=1 , where
craft. We also note that our algorithm does not have access
φ̃i − φ̃ is the direction toward the i-th training shape φ̃i . Note to any information about which parts of the image contain
occlusions.
that the weight k(dL2 (φ̃, φ̃i ), σ) for the velocity component
φ̃i − φ̃ increases as φ̃ gets closer to φ̃i . 4.2 Segmentation of Multiple Classes of Objects
One issue with the evolution in (9) is that the evolving
level set function does not necessarily remain on the man- We now consider the case where the shape prior density in-
ifold of signed distance functions D. The evolution of the volves multiple classes of shapes. This is a scenario which
zero level set (i.e. the curve) in such a case can be less stable cannot be readily handled by most existing shape-based seg-
than the case where the evolving level set function remains mentation techniques. In particular, we consider the prob-
a signed distance function [10]. Hence, it is desirable to con- lem of segmenting handwritten digits, where there are 10
strain the evolving level set function to stay on the manifold. classes of digits, i.e. 0, 1, . . . , 9. We use a training set of 100
Based on these observations, we modify the evolution sample images with 10 segmented images of each digit.
equation (9). The goal here is to extract relevant informa- Let us consider the low-SNR test images (not included in
tion for shape evolution from the velocity field in (9) and to the training set) in Fig. 3(a). Segmentations without a shape
construct a new velocity field such that the resulting trajec- prior (and with a curve length penalty instead) are shown
tory is on D. First we observe that the only components of in Fig. 3(b). The results of our shape-based segmentation
 (a) (b) (c) (a) (b) (c)

Figure 3: Segmentation of low-SNR digit images. (a) Test Figure 4: Segmentation of digit images with missing data.
images. (b) Results without a shape prior. (c) Results with (a) Test images. (b) Results without a shape prior. (c)
the proposed shape prior. Results with the proposed shape prior.
method are shown in Fig. 3(c), which appear to be quite ac- in Proc. IEEE Conf. on Computer Vision and Pattern
curate despite the low-quality data. Finally, we consider an Recognition, 2000, pp. 316–323.
example involving test images with missing data, as shown [3] A. Tsai, A. Yezzi, Jr., W. Wells, C. Tempany, D. Tucker,
in Fig. 4(a). Our shape-based segmentation results shown A. Fan, W. E. L. Grimson, and A. S. Willsky, “A shape-
in Fig. 4(c) provide accurate segmentations despite the data based apprach to the segmentation of medical imagery
limitations. using level sets,” IEEE Trans. on Medical Imaging, vol.
5. CONCLUSION 22, no. 2, pp. 137–154, Feb. 2003.
We have considered the problem of estimating shape prior [4] D. Cremers, S. Osher, and S. Soatto, “Kernel density
densities from example shapes and proposed a shape-based estimation and intrinsic alignment for knowledge-driven
segmentation method. In particular, we have developed a segmentation: Teaching level sets to walk,” in Pattern
framework for estimating shape priors from training shapes Recognition, Aug. 2004, vol. 3175 of Lecture Notes in
nonparametrically. Based on such nonparametric shape pri- Computer Science, pp. 36–44.
ors, we have formulated the shape-based segmentation prob- [5] D. Mumford, “Mathematical theories of shape: Do they
lem as a MAP estimation problem. Evaluation of the non- model perception?,” in SPIE Proceedings Vol. 1570 Ge-
parametric shape prior for a candidate curve for segmenta- ometric Methods in Computer Vision, 1991, pp. 2–10.
tion is given in terms of distances between the candidate [6] J. Kim, Nonparametric Statistical Methods for Image
curve and the training curves. We considered the template Segmentation and Shape Analysis, Ph.D. thesis, Mas-
metric and the L2 distance between signed distance func- sachusetts Institute of Technology, Feb. 2005.
tions, but other metrics can also be used for nonparametric [7] P. W. Minchor and D. Mumford, “Riemannian geome-
shape priors. We have derived curve evolution equations tries on spaces of plane curves,” J. Eur. Math. Soc.
based on the nonparametric shape priors. We have pre- (JEMS), to appear.
sented experimental results of segmenting partially-occluded
images. We have considered the case in which the training [8] E. Klassen, A. Srivastava, W. Mio, and S. H. Joshi,
shapes form multiple clusters, and demonstrated that our “Analysis of planar shapes using geodesic paths on
nonparametric shape priors model such shape distributions shape spaces,” IEEE Trans. on Pattern Analysis and
successfully without requiring prior knowledge on the num- Machine Intelligence, vol. 26, no. 3, pp. 372–383, Mar.
ber of clusters. 2004.
[9] T. Chan and L. Vese, “Active contours without edges,”
REFERENCES IEEE Trans. on Image Processing, vol. 10, no. 2, pp.
266–277, February 2001.
[1] T. F. Cootes, C. J. Taylor, D. H. Cooper, and J. Gra-
ham, “Active shape models— their training and appli- [10] S. Osher and R. Fedkiw, Level Set Methods and Dy-
cation,” Computer Vision and Image Understanding, namic Implicit Surfaces, Springer, 2003.
vol. 61, no. 1, pp. 38–59, 1995. [11] H.-K. Zhao, T. Chan, B. Merriman, and S. Osher, “A
[2] M. E. Leventon, W. E. L. Grimson, and O. Faugeras, variational level set approach to multiphase motion,”
“Statistical shape influence in geodesic active contours,” Journal of Comp. Phys., vol. 127, pp. 179–195, 1996.