Baddeley, A.J. (Adrian)

Lieshout, M.N.M. (Marie-Colette) van

English

CWI's Institutional Repository

241 MARKOV CHAIN MONTE CARLO METHODS FOR CLUSTERING OF IMAGE FEATURES MN M van Lieshout University of Warwick United Kingdom 1 Introduction The identification of centres of clustering is of iuterest in many areas of applications, for instance .. dge detector output has to be grouped into mean-ingful curves. In this paper we argue that stochas-r 1c geometry models are helpful both in providing ruodels for clustering and as a prior distribution to combat overestimation of the number of clusters and to improve robustness. The general idea in connection with object recognition was proposed by Baddeley and Van Lieshout [1]. See also Van Lieshout [13]. Inde-pendently, in an epidemiological context, a differ-Pllt Gibbs sampler technique for detection of cluster centres in a Cox process was developed by Lawson [11). Earlier attempts include [3, 5, 14, 15). For 1·xt.ensions and more examples see Lawson et al. [12). 2 Cluster processes The data consist of a set of 'events' y = {y1 , ... , Ym} ~ T, where T ~ IR-2 (say) is the win-dow of observation and it is required to determine r lie locations of an unspecified number of cluster centres x = {x1 , •. • ,xn} ~ U, n:::: 0. In the sim-plest case, U = T or a somewhat larger set to ac-,·nunt for edge effects, but note that U may be a different space as in Section 5 below. In probabil-ity theory, a set of points such as x or y are real-isat.ions of a random point process. An arbitrary dummy centre xo is introduced both for technical reasons, and to allow for events not belonging to anv cluster . • The mathematical model we adopt is an inde-pendent cluster model (see eg. [19]) where, condi-tional on x = { x 0, . ... Xn}, the observations result from the superposition of independent, finite point processes y = LJzx,· i=O For simplicity, we confine ourselves here to the case where each set of 'daughters' Zu, u E U, is a real-Image Processing And Its Applications, 4-8 July 1995 Conference Publication No. 410, ©IEE 1995. A J Baddeley University of WA Australia isation of an inhomogenous Poisson process on T with intensity function h( · I tt) : T --+ [O, oo). By the superposition property, the combined offspring form a Poisson process with intensity function n ,\(-1 x) = Lh(·- x;). i=O In other words, given a fixed, finite refer-ence measure µ on T, representing determinis-tic but spatially varying factors, the number of events n(y) in y is Poisson distributed with mean fr >.(t I x) dµ(t) and conditionally on n(y) = m the events have joint density D.i=1 >.(yj Ix) fm(Y Ix)= (JT ,\(t Ix) dµ(t))"' with respect to µm. 3 Detection of cluster centres We interpret the identification of centres of clus-tering as a statistical estimation problem. The unknown configuration of centres is regarded as a model parameter, and upon observation of y, has likelihood l(x I y) = f(y I x). A maximum likeli-hood estimator solves x = argmax xf(y Ix), but note that there is no guarantee these <"quations have a solution, nor that a solution is unique. In-deed, a maximum likelihood estimator of x may run into difficulties similar to those encountered in the context of object recognition [l]. If his smooth and almost flat near its maximum, and the dat.a pattern is 'dense', the ma.ximum likelihood estimate tends to contain multiple responses to each true cluster. To overcome these problems we propose a prior distribution to penalise scenes that contain too many 'similar' centres. A suit.able choice is a Markov spatial process [2, 17). For brevity, we con-sider only pairwise interactions models, defined by their density (1) (with respect to a unit rate Poisson model). Here, a- is the normalising constant, f3 > 0 a model pa-rameter. The product ranges over all pairs of sim-ilar objects x; ,..., Xj and g(-,-) : U x U - [O, oo) is the interact.ion function. The case g = 1 is a Poisson process with rate f3; for g(-, . ) < l, config-urations with many similar centres are unlikely. An important property of (1), and of Markov models in general, is that the likelihood ratio p(xU{u})_aIT (. ) (x) - fJ g x,, u P x;-u depends only on those Xi E x that are similar to u, signifying that all interaction is 'local'. In the statistical physics interpretation, - log p(xU { u}) + logp(x) is the energy required to add.a new point u to an existing configuration x; in probabilistic terms >.(u;x) = p(xU {u})/p(x) is the Papan-gelou conditional intensity at u given the rest of the pattern x on U \ {u}, see (6]. Roughly speaking, >.(u; x)du is the conditional probability of a point in the infinitesimal region du centred at u given the configuration agrees with x outside this region. After collection of the data, the posterior density of x is p(x I y) ex f(y I x)p(x) by an application of Bayes' formula. However, due to the normalising constant, this distribution tends to be rather intractable and we have to resort to Markov chain Monte Carlo methods [4] for opti-misation and sampling. The basic idea is to build a Markov chain (or in continuous time, a Markov process) with the target distribution p(x J y) as equilibrium and transitions that are 'easy' to per-form. · The classical approach in a point process context is via spatial birth-and-death processes [2, 16, 18] but other methods such as Metropolis-Hastings [9], simulated tempering or annealing can be used as well. Furthermore, functionals of the posterior dis-tribution such as the distribution of the number of clusters, the probability that there is no cluster in a particular region, and the first-order intensity of cluster locations can be estimated. All these techniques iteratively update the clus-ter centre configuration by addition and deletion, with transition criteria based on the posterior like-lihood ratios. For instance, a spatial birth-and-death process is a continuous time Markov process for which the only transitions are the birth of a new object (instantaneous transition from x to x U { u}) or the death of an existing one (transition from x to x \ { x;}). Transitions are governed by death rate D(-., ·) and birth rate B(., ·) as follows. • the probability of a death x - x \ { x;} during a time interval (t, t + h), h - 0, is D(x \ {x;}, x;)h + o(h); 242 • the probability of a birth x - x U { u} during time (t, t + h), where u lies in a given subset F ~ U, is B(x, F)h + o(h); • the probability of more than one transition during (t, t + h) is o(h). We will assume that B(x,·) has a density b(x,·) with respect to some finite reference measure v on U, so that intuitively b(x, u) is the transition rate for a birth x - x U { u}. Good choices for the transition rates are b( ) =p(xU{u}Jy) x,u p(xJy) and, writing n(x) for the number of non-dummy centres in x, D(x\{x;}, x;) = 1/n(x), which results in detailed balance between births and deaths: b(x, u)p(x J y) = D(x, u)p(x U {u} I y). To avoid explosion, i.e. an infinite number of transitions occurring in finite time, the rates have to satisfy certain assumptions [16]. Typically in this context we need an inhibitory prior g( ·, ·) 5 1 and an upper bound on the daughter intensity h(-) [13]. The Metropolis-Hastings method is a two-step discrete time chain, which proposes to change the current state x to a new candidate state x', ran-domly sampled from a probability density q(x, ·). The proposal x' is accepted with probability A( ') _ . {l p(x' I y) q(x', x)} x, x - mm , ( I ) ( , ) . p X y q X,X It is easily verified that transitions out of state x are belanced by transitions into x: p(x' J y)q(x',x)A(x',x) = p(x I y)q(x,x')A(x,x'). The transition densities q(-, ·) are built as fol-lows: • with probability q(x) generate a new point u from a density b(x, u) with respect to 11; • otherwise (with probability 1 - q(x)) delete a point x; E x at random; The simplest choices for b(-, -) and q( ·) are 1 1 q(-) = 2' b(·,·):: v(U)' but this may tend to generate too many proposals with low acceptance probabilities. Another possi-bility resembling spatial birth-and-death processes is to set b(x,u)= /(yJxU{u})p(xU{u}) 1 f(y I x) p(x) B(x) and where x _ B(x) q( ) - n(x) + B(x)' B(x) = f p(x u {u} I y) dv(u) lu p(x!y) or some other density that tends to select proposals. u where p(xu {u} I y) is large. Under mild conditions, the resulting Markov chains converge to p(x I y) [9]. From a computational point of view, to (say) add u E U to x we need to evaluate f(y I xUx{u}) = H(u) IT [i + nh(YJ I u) ] , J(y I ! J=l I:i=O h(YJ I Xi) where H(u) =exp {-fr h(t I u) dµ(t) }, and p(xu {u}) p(x) (2) (3) Ratios (2) are straightforward to calculate and can be compared to the Hough transform [10] in the sense that each point YJ votes with variable strength for a cluster centre at point u. Ratios of the form ( 3) are easy to compute if p( ·) is a nearest-neighbour Markov point process (1). In particular, the intractable normalising constant cancels out. Parameters in f and p can be estimated in ad-vance or during iteration. Alternatively, a Gibbs sampler can be specified but note that this ap-proach requires an extra set of prior distributions, one for each parameter. 4 Offspring labelling So far, we concentrated on estimating the pro-cess x of cluster centres, but it might well be of interest to label the observed events by the cluster they belong to. This opens the possibility to esti-mate functionals such as the probability that two data points are siblings, or the distribution func-tion of offspring displacements. In order to be able to do inference on cluster membership, we need the full distribution of w = {(xo,Z:r0 ), ••• ,(xn),Zx,.)}, parents {x1, ... ,xn} marked by their offspring Z:r,, i = 0, ... , n. Recall that (in contrast to e.g. the analysis of fi-nite mixture models [7]) the number of points in the pattern is not fixed in advance. This is quite cru-cial, as it implies the impossibility to build a Gibbs sampler, alternatingly sampling the marks (or cen-tres) given data and centers (or marks). We pro-pose a two-step Metropolis-Hastings algorithm (13], where the transitions consist of births and deaths of centres, followed by an adjustment of the marks, 243 by drawing from the conditional distribution of off-spring. It is important to note that, when a new configuration x' is generated from the current x and offspring Z, x' does not have to conform to Z ( eg. x' could have one point fewer than x). We then update the daughter sets to Z' and require confor-mity with x'. Specifically, we have the following algorithm. Given an initial paren~ configuration x, a labelling Zx = (Z:r,)~~~) and data y, alternate between a Metropolis-Hastings and an adjustment step as follows. 1. given w = {(xo,Z:r0 ), ... ,(xn,Zz,J} and y, . perform a Metropolis-Hastings step for x only, independent of the marks Zx, yielding the new configuration x'; 2. given x', y reallocate the offspring by a Gibbs step, sampling Z~, from the conditional dis-tribution P(Z~, I x', y). Replace x by x', Zx by Z~, and return to step 1. The detailed balance equations reduce to those for the centres only (see Section 3) and hence con-. vergence of the Markov process to p( w I y) follows under similar conditions. In step 2, an assignment of data to clusters can be interpreted as an ordered partition </> : y - x of the data. The labels t/J( ·) ascribing each obser-vation to a cluster centre are independent, with probabilities h(YJ I Xq,(y;)) I:i h(yj I Xi). Hence, simulation of the mark process is easy. The independence of labels provides a justifica-tion for the nearest-parent classifier [11]. Given a set of sites y = {Y1, ... , Ym} and a set of parents x = {xo, ... , xn}, suppose the task is to assign to each Yi a label <P(j) E {O, ... , n} corresponding to centre xr/>(fr By independence, a maximum likeli-hood classifier is ~(j) = argmax i=O, ... ,nh(y; I Xi), (4) which assigns each point to the closest centre if h(y; I Xi) is a decreasing function of llY; - Xiii To conclude this section, if one considers un-ordered partitions of the data instead of unordered ones, a Gibbs sampler [4] can be developed. In-tuitively, these unordered partition elements corre-spond to 'sibling sets' generated by a common, if unspecified, parent. Then, by allowing parents not having any offspring at all, the number of cluster centres can be changed, but the conditional distri-butions involved are rather awkward [13]. 5 Applications The simplest examples arise in spatial statistics where T <;;;; U are bounded regions in the plane. Typical applications include forestry, where the task is to reconstruct the ancestors of the current generation of trees in a wood or the epidemiological analysis of rare diseases. As an illustration, Figure 1 shows the locations of 62 redwood seedlings in a square of side approx-imately 23 m. The data was extracted by Ripley [18] from a larger data set in Strauss [20]. A biolog-ical explanation [20] for the apparent clustering is the presence of stumps known to exist in the plot, but whose position has not been recorded. Our goal is to reconstruct these centres of clustering. ·. I . . · . . . . ... . . . .· • I . . . Figure 1: Positions of 62 redwood seedlings in a unit square (Ripley 197(). Previous analyses of this data set include Strauss [20] who attempted to fit a pairwise interaction Markov model, Ripley [18], Diggle [8] and Lawson [11]. From a biological point of view, the above works tend to find an implausibly large number of stumps (clusters); 25 in [8] and 16 in [11], mak-ing a case for a Bayesian analysis with a prior dis-tribution penalising configurations with too many stumps close together. The analysis below is taken from Van Lieshout [13]. Following [8, 11] we assume the number of daughters per parent is Poisson and seedlings follow a radially symmetric Gaussian distribution around their ancestor. In contrast to the aforementioned papers, a pairwise interaction Markov prior ( 1) with strict inhibition g :: I < 1 is introduced. We used an interaction distance 1· = .084 [8], that is u ~ v if and only if u and v are less than r apart. Using a spatial birth-and-death sampler, a real-isation from the posterior distribution of stumps is shown in Figure 2. Figure 3 displays the estimated 244 posterior intensity, roughly speaking the posterior 'probability' that a point u belongs to x. For rea-sons of clarity, here black corresponds to high val-ues. > ' . : ~ . . . ' . .. . . . Figure 2: Realisation from the posterior distribu-tion taken after 2 time units (black), for a gaussian model with µ = 6.5, CT = .05 and a prior with log,B = log/ = -10, r = .084. The data is dis-played in grey . . :: .·.:--> . Figure 3: Posterior intensity of redwood seedlings (Ripley) estimated over 50 time units, for a gaus-sian model with µ = 6.5, CT = .05 and a prior witli log ,B = log/ = -10, r = .084. For an analysis of the full data [20] see (13]. As another example [1], suppose that the data again consist of a point pattern in a bounded region T <;;;; IR.2, but that the points are believed to lif' close to a curve and the objective is to estimate the curves. This includes the image analysis task of joining a dot pattern into a curvilinear boundary, but also eg. the identification of ancient roads or trade routes given information about the location of archaeological finds such as pottery or coins [19, p. 139], or the analysis of earthquake occurrences in relation to geographical fault patterns. Our final example is the problem of identifying large scale edges in a scene using the output of a low-level edge detector. The 'data' y consist of a pattern of line segments and the objective is to cluster them around a small number of larger line segments [14]. Let T denote the set of possible outputs of the low-level edge detector. For example these may be line segments restricted to have unit length (= 1 pixel width) and orientation which is a multiple of 45 degrees. The space U of objects we are looking for in this case also contains line segments, but of unrestricted length and orientation. Hence, y is a superposition of conditionally in-dependent line segment processes Zx, associated with each true line segment Xi- Typically the ex-pected number of segments in Zx, will depend on the length of Xi. The benefits of a prior model for x include the ability to encourage long lines and continuity between lines, and to penalise lines that cross one another. References (1] A. J. Baddeley and M. N. !VI. van Lieshout. Stochastic geometry models in high-level vi-sion. K. V. Mardia and G. K. Kanji (Eds.) Statistics and Images, Volume 1. Journal of Applied Statistics, 20:233-258, 1993. [2] A. J. Baddeley and J. M~ller. Nearest-neighbour Markov point processes and random sets. International Statistical Review, 57:89-121, 1989. [3] M. Baudin. Note on the determination of clus-ter centres from a realizat.ion of a multidimen-sional Poisson cluster process. Journal of Ap-plied Probability, 20: 136-143, 1983. (4] J. Besag, P.J. Green, D. Higdon and K. Mengersen. Bayesian computation and stochastic systems. Statistical Science, 1995. To appear. (5) J. Besag and J. Newell. The detection of clus-ters in rare diseases. Journal of the Royal Sta-tistical Society, Series A, 154:143-155, 1991. [6) D.J. Daley and D. Vere-Jones. An introduc-tion to the theory of point processes. Springer Verlag, New York, 1988. [7] J. Diebolt and C.P. Robert. Estimation of finite mixture distributions through Bayesian 245 sampling. Journal of the Royal Statistical So-ciety, Series B, 56:363-375, 1994. (8) P. J. Diggle. Statistical analysis of spatial point patterns. Academic Press, London, 1983. (9] C. J. Geyer and J. M!llller. Simulation proce-dures and likelihood inference for spatial point processes. Research report 260, Mathematical Institute, University of Aarhus, january 1993. [10] [11) (12] [13] J. Illingworth and J. Kittler. A survey of the Hough transform. Computer Vision, Graphics and Image Processing, 44:87-116, 1988. A. Lawson. Discussion contribution. Jour-nal of the Royal Statistical Society, Series B, 55:61-62, 1993. A. B. Lawson, M. N. M. van Lieshout and A. J. Baddeley. Markov chain Monte Carlo methods for clustering of spatial point pro-cesses. In preparation, 1995. M. N. M. van Lieshout. Stochastic geometry models in image analysis and spatial statis-tics, PhD thesis Vrije Universiteit Amsterdam, 1994. To appear as monograph in CWI-tract series, 1995. (14] P. Nacken. A metric for line segments. IEEE Transactions on Pattern Analysis and Ma-chine Intelligence, 12:1312-1318, 1993. (15) [16] S. Openshaw, M.G. Charlton, C. Wymer, and A.W. Craft. A mark 1 geographical analysis machine for the automated analysis of point data sets. International Journal on Geograph-ical Information Systems, 1:335-358, 1987. C. J. Preston. Spatial birth-and-death pro-cesses. Bulletin of the International Statistical Institute, 46:371 - 391, 1977. [17] B. D. Ripley and F. P. Kelly. Markov point processes. Journal of the London Mathemati-cal Society, 15:188-192, 1977. [18] B. D. Ripley. Modelling spatial patterns (with discussion). Journal of the Royal Statistical Society, Series B, 39:172 - 212, 1977. [19] D. Stoyan, W. S. Kendall and J. Mecke. Stochastic Geometry and its Applications. John Wiley and Sons, Chichester, 1987. (20) D.J. Strauss. A model for clustering. Biometrika, 62:467-475, 1975. 

Markov chain Monte Carlo methods for clustering of image features

Abstract

Similar works

Full text

Available Versions

CWI's Institutional Repository