Parameter estimation for a deformable template model.
In recent years, a number of statistical models have been proposed for the purposes of high-level image analysis tasks such as object recognition. However, in general, these models remain hard to use in practice, partly as a result of their complexity, partly through lack of software. In this paper we concentrate on a particular deformable template model which has proved potentially useful for locating and labelling cells in microscope slides Rue and Hurn (1999). This model requires the specification of a number of rather non-intuitive parameters which control the shape variability of the deformed templates. Our goal is to arrange the estimation of these parameters in such a way that the microscope user's expertise is exploited to provide the necessary training data graphically by identifying a number of cells displayed on a computer screen, but that no additional statistical input is required. In this paper we use maximum likelihood estimation incorporating the error structure in the generation of our training data.
INTRODUCTION
The use of high-level models in image analysis has increased in popularity during the last few years as
the problem of interpretation and understanding of images has become more important. One specific
task for which such models have been proposed is object recognition, locating objects in a scene and
providing estimates of their type, position and outline. The confocal microscopy image in Figure 1 is
one such application, showing an optical section through cartilage growth plate in which two types of
cell are present; one stage of the analysis of these data is to identify and distinguish cells on the basis
of their shape and size. Manual segmentation of each cell would be a tedious and time consuming task,
especially as this image forms just one of a large number of similar images, and so there is a great need
for some automated procedure. A number of different approaches, both statistical and non-statistical,
exist for tackling the task. In this paper we concentrate on a high-level statistical model which represents
each cell by a polygon outline model embedding these deformable templates into a marked point
process to handle the unknown number of objects. A point process of stochastic geometrical objects
is thus constructed which is conditioned on the observed data, i.e. the data image, and simulation techniques
are then used to generate samples from this posterior distribution. These samples are used to
estimate quantities such as the number of objects, their type, outline and position.
Although there has been much progress in understanding such models recently, they remain somewhat
inaccessible in part due to the difficulty in specifying many of the parameters associated. Further work
is needed before microscope users could routinely use these techniques for automatic extraction of in-
formation. A significant step towards this goal would be a way of automating parameter estimation,
both for the parameters describing the stochastic variations in shape for the objects and for the parameters
in the data model. The aim of this paper is to provide such a solution in a way that the only
interaction between the microscope user and the estimation procedures is purely graphical; the end-user
supplies computer-drawn templates of the cells and the training data for parameter estimation
then take the form of polygons fitted on the screen to a display of image data.
The plan of the paper is as follows: In Section 2 we review the high-level deformable template model
and we propose a likelihood model suitable for confocal microscopy applications. In Section 3 we
briefly discuss inference from the model using reversible jump Markov chain Monte Carlo methods.
We then discuss parameter estimation in Section 4, treating the image degradation parameters in a fully
Bayesian way and using likelihood methods for the remaining parameters. We concentrate on deriving
the maximum likelihood estimates of the deformation parameters under the assumption that the training
data are recorded with measurement error due to the imperfect handtracking of individual cells.
We conduct studies of this latter estimation problem in Section 5, using artificial data for which comparisons
can be made with known values, and also using the cartilage data set.
Adopting a Bayesian approach, we denote the scene by x and the recorded pixellated data by y, the posterior
of interest being given by -(xjy) / -(yjx) -(x). Since we are concerned with the components,
i.e. the cells, of the scene, x will be specified as a "high-level" variable, composed of an unknown
length list of the constituent components of the scene. This approach has been pioneered by Grenander
and his co-authors, and recent sources for references can be found in Grenander & Miller (1994),
Jain, Zhong & Lakshmanan (1996), Baddeley & Van Lieshout (1993) and Blake & Isard (1998).
2.1 THE OBJECT PRIOR -(x)
The prior model for a random configuration of cells is built around a stochastic template model for a
single cell embedded in a mixture model of different cells types which is used as the mark distribution
of a marked point process model.
2.1.1 THE TEMPLATE MODEL FOR A SINGLE CELL
The prototypical cell is described by an n-sided template defined by a set of vectors
which give the edges of the polygon (Figure 2a). For example, if one type of cell is characteristically
circular, then these edges describe a polygonal approximation to a circle. The closure of the polygonal
template is equivalent to the condition
Described in this way, the template does not
have any location information, and so we will consider its first vertex to be located at the origin, the
second to be located at g 0 , the second at and so on (Figure 2b). In order to accommodate
scaling and rotational effects, the template may be globally scaled by a scalar R and rotated through
an angle ff (Figure 2c). To model natural shape variability occuring between cells of the same type,
each edge g i is subject to a stochastic deformation which incorporates the scaling R, the rotation ff
plus an additional edge-specific Gaussian deformation in length and direction (Figure 2d). This final
edge-specific effect describes the change in length and direction between the undeformed g i and the
new edge prior to global scaling and rotation. Writing the deformed edge as s Rff
i is the
representing these changes, we thus haveR
sin ff cos ff
s Rff
Writing t (0)
determines that
s Rff
\GammaR sin
Specifying the distribution of r i and ' i to be the angular and radial components of a bivariate Gaussian
with zero correlation, then conditional on fixed values of R and ff, t (0)
are independent
Gaussians with mean zero. Ignoring for a moment the constraint that the deformed template must be
closed, i.e.
suggest a first order cyclic Markov
structure on the ft (0)
and the ft (1)
independently with each having an n-dimensional Gaussian distribution
with mean 0 and circulant inverse covariance matrix incorporating the scaling
where all other entries are zero.
Define the length 2n vector
considering only the
unconstrained non-closed polygon case
Imposing the closure constraint on the deformed template will destroy the simple structure of (2); we
will derive the distribution under the closure constraint in Section 4.2.2.
2.1.2 THE MARKED POINT PROCESS FOR MULTIPLE OBJECTS
It is possible to allow objects of different types by forming a mixture model using different basic templates
with different deformation parameters, where the weights in the mixture model represent the
relative frequencies of occurrence; refer to Rue & Hurn (1999) for details. This mixture distribution
is then used as the mark distribution of a marked point process model to model an unknown number
of objects (Baddeley & Van Lieshout, 1993). Within this framework, each object x i comprises a point
which gives location (unspecified by the deformation model) and a set of marks which then specify
its other attributes, in this case type and outline. The points lie in a window   related to the data observation
coordinates, and the marks in the space M associated with the object types and shapes. A
configuration of cells is described as a finite unordered set an object
process, i.e. a marked point process on   \Theta M with a Poisson object process as the basic reference
process. Under the Poisson process, conditional on the number of objects k, the objects are uniformly
and identically distributed. The joint density of fx is defined by a density f(x;
to the Poisson object process. For example, to model pairwise interactions between objects which are
defined to be neighbours by some relation -, the function
Y
could be used; Van Lieshout (1995) discusses various interaction models. To model a situation where
objects are not allowed to overlap, such as confocal microscopy where the images are optical sections,
all objects are defined to be neighbours and the interaction function h(x taken to be zero if
objects x i and x j overlap and one otherwise. This model is known as the hard core object process. We
take the point, denoted c, to be the location of the first vertex for each object. The first mark is then
a categorical variable defining the object type, and the remaining marks are the deformations t of the
standard template for this type. We refer to Rue & Syversveen (1998) for more details.
2.2 A POISSON LIKELIHOOD MODEL -(yjx) FOR CONFOCAL MICROSCOPY
Confocal microscopy can produce relatively high-quality sectional images for certain specimen types.
The technique works by scanning the focal point of a laser through the stained specimen, the fluorescence
emitted from this focal point is then recorded, and the image is built up as a grid of fluorescent
dye concentration measurements. The likelihood model for the data y approximates the microscope
measuring fluorescence photon counts. We begin by assuming that cells all have a typical fluorescence
level -, while the background emits no fluorescence. The signal is measured on an 8 bit scale (i.e. on
the integers 0 to 255) and in order to minimise clipping, the signal is operator adjusted at the experimental
stage by a scale factor a and a shift b, so that the mean background value is now b, while the mean
foreground value is b (and both are sufficiently far from the ends of the range of measurement).
In addition, the recording process involves additive instrument white noise. We will approximate the
record at pixel p by a Gaussian random variable with mean equal to the black-level b plus the expected
scaled Poisson count for p, and with variance equal to the variance of the instrument noise oe 2 again
plus the scaled expected Poisson count variance,
where I [p2cell] is the indicator or whether the pixel is labelled as cell or background. The fy p g are
assumed to be conditionally independent given the configuration x.
To do inference for the high-level object model, we have to resort to Markov chain Monte Carlo techniques
(refer to Gilks, Richardson & Spiegelhalter (1996) for an introduction and list of references.)
The key idea is to construct a Markov chain which has the posterior as the equilibrium density, and to
estimate expected values by empirical average along the trajectory of a realization of the Markov chain.
This makes it possible to access both estimated expected values and their uncertainty represented by
the estimated variance of the estimated marginal density. To construct this Markov chain, we need a set
of basic transition rules that will take us around in the state-space, hopefully moving efficiently around
and thereby creating good estimates with low variance for the empirical averages. The basic transition
types for our problem are as follows:
ffl Change a part of the template boundary of a randomly chosen object.
ffl Change the type of an object.
ffl Propose the birth of a new object.
ffl Propose the death of an object.
The details of the MCMC algorithm for these transitions are given in Rue & Husby (1998), Rue &
Syversveen (1998), where the birth and death moves rely on applications of the Green (1995) algorithm
for variable dimension distributions. In practice the basic transitions listed above do not ensure
sufficient mixing of the Markov chain, the main problems being transitions between the various object
types and difficulties when there is uncertainty over the whether one "object" is actually two objects
lying close together. Rue & Hurn (1999) approach this problem by extending the template model to allow
it to have a variable number of sides, demonstrating that this allows smoother transitions between
the object types, and easing the transitions which split and merge objects.
The model described in the previous sections relies on the specification of the parameters controlling
the marked point process prior and the degradation model. Some of these parameters, particularly those
involved in the deformations, are hard to interpret in any intuitive way. In this section we propose a parameter
specification approach which will hopefully allow models of this type to be more widely used
in practice. The parameters of the degradation process are handled in a Bayesian way. The remaining
parameters will be estimated using maximum likelihood methods.
4.1 BAYESIAN TREATMENT OF THE DEGRADATION PARAMETERS
The unknown parameters of the measurement process are - the typical fluorescence level, a and b the
operator controlled measurement adjustments, and oe 2 the recording noise. Since estimation of these
parameters from a section of image data would be dependent on a good initial segmentation, we propose
instead to treat them as hyperparameters. Our choice of prior distributions is intended mainly to
impose the constraints that the variables are positive, with 255. The noise variance
rather than the precision is given a - 2 (1) prior distribution since it is expected to be of the order of a
few squared pixel units. The mean fluorescence is given a fairly diffuse - 2 (255) prior reflecting that it
is positive variable expected to take a reasonably large value, although this will vary considerably between
specimens. As the data carry considerable information regarding the parameters, it is intended
that these priors will not exert much influence; refer to Hurn (1998) for more discussion in a similar
setting,
a
The joint posterior distribution of interest is then
Markov chain Monte Carlo methods are needed for the sampling which generally will be done one
parameter at a time requiring the conditional posterior distributions
Y
Y
pixels p2cell
Y
pixels p2cell
Y
4.2 ESTIMATION OF THE PRIOR MODEL PARAMETERS
It is not computationally feasible to treat the remaining parameters in a similar Bayesian manner since
it would require the intractable normalising constants which depend on these parameters. The alternative
we suggest is to use likelihood-based methods which work with training data generated by the
microscopist tracing cell outlines on a number of slices of image data using the computer mouse. This
method has two benefits, the first is that the microscopist's expertise is being used in the data collection,
and the second is that there is no further user-interaction in the estimation procedure.
4.2.1 THE MARKED POINT PROCESS PARAMETERS p i AND fl
The use of a hard-core interaction model for the marked point process means that the normalising
constant will be too complex for maximum likelihood estimation. However, previous experience in
object recognition problems suggests that the results are not particularly sensitive to the specification
of these parameters. We propose to use the naive estimators, that is the observed proportions of the
different cell types in the training data for the fp i g, and the observed frequency per unit area for fl.
These estimators will be least accurate in situations where there is significant packing observed in the
training data. If this is the case, then an alternative is to simulate from the prior model (having first
estimated the deformation parameters) searching over a range of values for fp i g and fl to match the
observed and the simulated frequencies; this approach could be coded in advance. We also note some
recent work by Baddeley & Turner (1998) who use maximum pseudo-likelihood methods to fit various
point processes, although such an approach here would be considerably more complex.
4.2.2 THE DEFORMATION PARAMETERS fi AND ffi AND TRACKING ERRORS
We derive the likelihood function for fi and ffi using data which are a number of polygons the vertices
of which are observed with error. Although the inclusion of the observation error complicates this
derivation, it is a realistic assumption given the intended manner of collecting training data. We begin
by transforming from the deformation model for the polygon edges to the model for the corresponding
vertex locations. Recall that the first vertex defines the location of the entire polygon; if the first vertex
is at location then the second vertex v 1 is located at
s Rff
Note that there are vertices in the non-closed polygon. Considering the x and y components
separately, the vertices can be
x
y
\Gammag xg y
\Gammag x
\Gammag xg y
\Gammag x
\Gammag x
t (1)t (1)t (1)
where I is the 2n \Theta 2 matrix whose first column consists of n 1's followed by n 0's, and whose second
column is n 0's followed by n 1's, and v Rff is the vector of vertex x and y positions of the undeformed
template with first vertex at zero rotated through ff and scaled by R. We will write (3) in the form
We assume that we record the n vertices of the closed polygon with error
~
and we further assume that these errors are normally distributed with
Since the true location of the polygon is unobservable, we begin by finding the unconstrained distribution
of v given R; ff; \Phi; ffl 0 and the observable ~
which by (3) and (4) is the distribution of
Denoting the marginal variance of (ffl x
22 , the unconstrained distribution is
To find the distribution of (v for an arbitrary v   , we reorder
the components of v from x then y components to the vertex pairs, rewriting (5) as
~
where the partitioning of the mean and variance correspond to partitioning into the sets (v
and v n . Note that E(v n j(R; ff; \Phi; ~
by closure of the undeformed template v Rff . Denote
the partitioned inverse of the variance matrix
We will also rearrange the ordering of the fffl i g so that
By manipulation of the conditional density -(v
show in
the

Appendix

that
The particular v   in which we are interested is v 0 , i.e. ~ which case
To find the distribution of the vertices under the closure constraint, we use the known marginal distribution
together with (8), integrating out ffl 0 to show in the Appendix that
I \Phi 22
I T
is the first 2n \Gamma 2 components of the reordered I~
I is the
of the reordered I .
The distribution of (~
under the same conditioning uses (4). Since (ffl
may not be independent of ffl 0 , the resulting distribution is
I \Phi 22
I
I \Phi T
I
The distribution of ~ v 0 arises as the convolution of the distribution of v 0 with the distribution of ffl 0 .
Assuming that v 0 is uniformly distributed in the observation window, and that the variance of ffl 0 is
small in comparison to the window size, ~
v 0 is taken to be also uniformly distributed in the window.
Finally, under an assumption of independence of the polygon shapes, the likelihood for m cells will
be the product of the individual likelihoods.
Finally, since in many applications it will not be reasonable to assume that the scale and rotation of
each cell is known, the model allows each cell i to have an associated scaling R i and rotation ff i . We
will concentrate on the case where there is a common unknown R for all cells (as would be used in
an example where size was a distinguishing characteristic in cell recognition). Treating scale and rotation
as additional nuisance parameters adds parameters R; ff , to the optimisation in order
to obtain the maximum-likelihood estimates. However in the Appendix, we show that the rotation
optimisation may be done exactly for any set of fi; ffi; R; \Phi, thus reducing the computational burden.
5 EXAMPLES OF THE ESTIMATION OF THE DEFORMATION PARAMETERS
5.1 AN ARTIFICIAL EXAMPLE
We begin with an artificial example to demonstrate the properties of the estimation procedure since
it allows us to consider the effect of various model misspecifications. The basic template for our experiments
is an eight-sided regular approximation to a circle of radius 1, and the deformation model
takes the parameters it is analytically intractable to find the maximum
likelihood estimates of fi and ffi, Matlab (The MathWorks Inc, 1994) has been used for the numerical
optimisation; it can also be used for direct entry of vertex locations using the mouse.
5.1.1 THE EFFECT OF MODEL MISSPECIFICATION
In order to assess the importance of deriving the maximum likelihood estimates under the assumption
of measurement error in recording the vertices, simulated realisations from the deformed
template model have been generated using independent measurement errors with variance OE 2, and estimates
of fi and ffi have then been found from these data when misspecifying that OE 2

Table

6 shows
the mean and median estimates resulting from 250 repetitions of this procedure to assess the variability
of the estimates. Not surprisingly, as the true variance increases, the estimates become severely biased;
fi and - ffi both decrease since the observed shapes demonstrate greater variability and less smoothness
than expected.
Similarly, it is possible to fit a model with measurement error but for which the scale R and rotation ff
are assumed known but are incorrectly specified. The effect here is that the observed polygons do have
the expected shapes but that in the attempt to fit the incorrectly specified template, OE 2
1 is overestimated,
while the balance of -
fi to -
is perturbed to the point of degeneracy latter effect is due
to the extremely high correlations required between deformations which are needed in order to mimic
a rotation or a rescale.
5.1.2 ESTIMATION UNDER DIFFERENT ERROR MODELS
Having seen the importance of including measurement error in the modelling, we now investigate different
possible models for the tracking errors. It seems unlikely that the errors will in practise be independent
and so we consider the following three structures:
Independent errors,
First-order Markov structure for the variance,
where \Phi 1 has a circulant
structure with OE 1 down the diagonal, and OE 2 on the second diagonal
MODEL 3 First-order Markov structure for the inverse variance, \Phi
where \Phi 1 is
although note that OE 1 and OE 2 have different interpretations under model 2 and 3.
Our training data are a set of m simulated cell outlines, where m takes the values 25, 50 or 100. In
the simplest noise model, we have also found estimates when the scale and rotation are taken to be unknown
(the nuisance parameter R takes the value 1). Table 2 summarises the results of the numerical
experiments, while Figure 3 shows the pairwise scatterplots of the four estimated parameters for the
case under the independent noise model 1 with unknown scale and rotation. Both mean and
median estimates have been quoted in all cases to highlight the skewed nature of the sampling distribu-
tions, a feature which can also be seen in Figure 3 (marginalising by eye). Figure 3 also demonstrates
the high correlation structure within the estimates, particularly between fi and ffi and between the deformation
parameters and the observation error variance. Such relationships might be explained by noting
that in the unconstrained model whose inverse variance is given by (1) the variance of any particular
given the remaining t (0)
components is 1=fi while the correlation of two consecutive t (0)
given the
remaining components is \Gammaffi=fi. The near non-identifiability due to the partitioning of the variability
of the observed vertex locations into that due to the deformations and that due to the observational error
is one of the reasons why this estimation problem is particularly difficult.
5.2 ESTIMATION FOR IMAGING

Figure

4 shows the two template types together with 100 cell outlines traced by one of the authors from
a series of images of cell cartilage. For display purposes, the outlines have been relocated on a regular
grid, the scalings and rotations have not been altered. Using these training data, the proposed estimation
technique has been applied assuming unknown scale R and rotation and a common measurement
error for both template types. The same three measurement noise models have been considered. Table
3 gives the resulting parameter estimates together with the quantity \Gamma2\Theta the natural logarithm of
the likelihood at those parameter estimates; this quantity will be used for determining which of the
three models is most appropriate using the Akaike Information Criterion. Using this criterion, noise
model 3 is chosen. Figure 5 shows some simulations from this fitted model; the top row shows 50
realisation from each of the two deformed template models, the bottom row shows the same realisations
convolved with realisations from the fitted noise model. Note that the rotations here have been
set to the value zero for displaying these outlines, which consequently appear less ordered than those
in

Figure

4 as a result of the natural structure of the cartilage.
6 DISCUSSION
User-friendly techniques for parameter estimation in deformable template models are an important
issue if such models are to find more wide-spread acceptance in applied work. Recent work by Kent,
Dryden & Anderson (1998) has also considered the question of estimators for deformation parameters,
although valid only for circular objects observed without measurement error. We have demonstrated
that it is possible to treat all the parameters necessary for a Bayesian object recognition approach in
either a fully Bayesian way or requiring only graphical intervention by the microscope user. We hope
that approaches such as these will facilitate further work in the area of statistical image analysis in
collaboration with the end-users of such techniques.



--R

Practical maximum pseudolikelihood for spatial point patterns
Stochastic geometry models in high-level vision
Active Contours
Markov Chain Monte Carlo in Practice
Reversible jump MCMC computation and Bayesian model determination
Hands: a Pattern Theoretic Study of Biological Shapes
Representations of knowledge in complex systems (with discus- sion)
Confocal fluorescence microscopy of leaf cells: An application of Bayesian image anal- ysis
Object matching using deformable templates
Using circulant symmetry to model featureless objects
Bayesian object identification
of partly destroyed objects using deformable templates
Bayesian object recognition with Baddeley's delta loss
Markov point processes and their applications in high-level imaging (with discussion)
--TR
object recognition;Markov chain Monte Carlo;bayesian image analysis;deformable templates;parameter estimation;confocal microscopy