Image Anal Stereol 2001;20:71-78 Review Paper REVIEW OF RECENT RESULTS ON EXCURSION SET MODELS Department of Mathematics, University of Queensland, Brisbane, Queensland, 4072, Australia; Cooperative Research Centre for Mining Technology and Equipment, PO Box 883, Kenmore, Queensland, 4069, Australia; 3 Department of Statistics, University of New South Wales, Sydney, New South Wales, 2052, Australia E-mail: rjw@maths.uq.edu.au (Accepted June 7, 2001) ABSTRACT Many images consist of two or more “phases”, where a phase is a collection of homogeneous zones. For example, the phases may represent the presence of different sulphides in an ore sample. Frequently, these phases exhibit very little structure, though all connected components of a given phase may be similar in some sense. As a consequence, random set models are commonly used to model such images. The Boolean model and models derived from the Boolean model are often chosen. An alternative approach to modelling such images is to use the excursion sets of random fields to model each phase. In this paper, the properties of excursion sets will be firstly discussed in terms of modelling binary images. Ways of extending these models to multi-phase images will then be explored. A desirable feature of any model is to be able to fit it to data reasonably well. Different methods for fitting random set models based on excursion sets will be presented and some of the difficulties with these methods will be discussed. Keywords: excursion set, model fitting, random field, texture modelling. R ichard JW ilson12 and D avid JN ott3 INTRODUCTION The work described here was motivated by the need to model segmented images of core samples of sulphide ores, such as that shown in Fig. 1. (A description of the segmentation method used is given in Wilson and Rychlik, 1996). In this sample, there are three minerals present, giving three phases. It was desired to construct statistical models describing the spatial distribution of the phases in such core samples. The initial intentions were to use the fitted models for classification of the ores and to simulate input to numerical models for rock fracture. This paper was presented at the Xth International Congress for Stereology (Melbourne, November 1999) and reviews results found in Nott and Wilson (1997a, 1997b, 2000) and David Nott’s Ph.D. thesis, as well as presenting some additional unpublished insights and comments. In general, the pixels in these images take values in an index set I = {0, ...,k — 1}, where each element in I is associated with one of k different mineral types. Many different approaches have been developed to model the spatial behaviour of such multi-phase images (see Haindl, 1991 for a review). For binary images, random closed sets can be used (see, for example, Stoyan et al, 1995). One way of obtaining a random set is to consider the excursion set of a random field, obtained by truncating a random field at some level. As the images can be obtained at different resolutions, it is assumed that the random field is over a continuous domain to ensure the fitted models at different resolutions are consistent with each other. \s •* **"~ m ? < * 4 * * * ^ \p * WJ '* + * • • ¦ i * ^ Fig. 1. Core sample image of sulphide ore. In the following section, the types of random fields and the properties of their excursion sets will be discussed, including examples. In section 3, various methods for fitting excursion set models to binary images will be presented, with some comments regarding their strengths and weaknesses. In the final section, models for multi-phase images and methods 71 Wilson RJ et al: Excursion set models for fitting these models will be discussed. The details of the methods presented here can be found in Nott (1997) and Nott and Wilson (1997a, 1997b, 2000). EXCURSION SETS Suppose that Y(t);t G M2} is a stationary Gaussian random field with zero mean, variance one and covariance function R(t; 0), where 0 G Md is some vector of unknown parameters. The u-level excursion set of Y, Au(Y), is simply the set of points where Y is not less than u, Au(Y) = {t£R2:Y(t)>u}. The morphology of AuY) depends on u, the type of covariance function used and the value of 0 (see Adler, 1982 for geometrical properties of excursion sets). PROPERTIES As the distribution of a Gaussian random field is completely determined by its mean and covariance function, the distribution of the excursion set must also be considered in terms of these. For convenience, define the binary random field Bu (t) by Bu(t) if t E Au Y), otherwise. If Ru (t) denotes the covariance function (centered about the mean) of Bu (and hence of Au{Y)), it can be shown (see p. 27 of Cramer and Leadbetter, 1967) that Ruit) 1 2k R(t) exp -u2\ 1+z J 1 dz. (1) When u = 0, this expression simplifies to s-1 R0(t) = {2ny arcsin(R(t)). The relationship Eq. 1 is monotone, and hence invertible. Thus, the distribution of the random field Y(t) can be determined from u and Ru(t) (that is, from the second order moments of Bu(t)). In addition, the (joint) distribution of the binary random field over a bounded domain can be expressed in terms ofthat of the random field. Rather than using a Gaussian random field, a %2 random field could be used. For a x12 field, a similar relationship to Eq. 1 exists: Ru(t) n J0 1 /¦!Rt! exp \1z -exp 1 z V1 dz. (2) EXAMPLES By careful choice of the covariance function for the random field Y(t), properties of the image to be modelled can be aligned with corresponding properties of the excursion set, Au{Y). Some of the possibilities are as follows (see Section 2.5.1 of Cressie, 1993 for a more complete list; see Nott and Wilson, 1997b also): R(t)=exp{-01||t||e2}, 2, 01 > 0, 0 < 02 < 2), 1 _ 2 Mi _l_ 1 11t1 if I It II S ft 1 2 8, + 2 ft" Xi HlH ^ ü1> (3) 0, (where d = 1, 01 > 0), otherwise, (4) 1 + - R(t)- (where d = 2, 01 > 0, 02 > 0), R(t) = 012||t||K1(012||t| (5) (6) (where d = 1, and Kv(.) denotes the modified Bessel function of the second kind). In the following, Eqs. 3-6 will be referred to as the (generalised) exponential, spherical, rational quadratic and Whittle covariance functions respectively. Fig. 2 shows simulated excursion sets for these covariance functions using the level u = 0.5. The covariance functions above are isotropic. Anisotropy can be introduced through linear transformations of t; that is, replace t by At where A is a 2 x 2 matrix of parameters. This transformation stretches one axis relative to the other and rotates them. As well, periodicity can be introduced by multiplying the covariance function by cos(a1t1 + a2t2), where a1 and a2 are additional parameters. EXCURSIONS AND THE BOOLEAN MODEL There is a close association between the Boolean model and high level excursion sets of smooth Gaussian random fields. It is well known that, with suitable scaling, the high level maxima of such fields are distributed according to a Poisson random field. As well, the connected components of the excursion set, again after suitable scaling, are approximately distributed as independent ellipses (the grains in the 72 Image Anal Stereol 2001;20:71-78 Boolean model), with their major and minor axes determined by the second order spectral moments of the covariance function (see, for example, Wilson and Adler, 1982). Using the methods of Leadbetter and Hsing (1990), it can be shown that the random measure defined by the area of the excursion set and suitably scaled converges to a compound Poisson random field with the mark distribution given by that of the area of the random ellipses. (e) (f) Fig. 2. Simulated excursion sets with u = 0.5. (a) Exponential model: Q1 = 0.02, 02 = 1 (b) Exponential model: Q1 = 0.02, 02 = 2 (c) Rational quadratic model: 01 = 8, 02 = 0.5 (d) Rational quadratic model: Q1 = 12, 02 = 2 (e) Whittle model: 01 = 0.3 and (f) Spherical model: Q1 = 40. (Reproduced with permission from Elsevier Science - see Nott and Wilson, 1997b.) ESTIMATION FOR TWO PHASE IMAGES There are a number of approaches to fitting excursion set models to images. The first choice to be made is that of a covariance function. This can be done in an ad hoc manner, using experience of the morphological properties of the excursions for each of the various possibilities or by fitting each in turn and determining which achieves the best result. The first requires the development of some exploratory data analysis tools, which are currently being considerd. Upon choosing an appropriate covariance function, there are a number of methods for estimating the parameters. The two most commonly used are (a) to use least squares to fit some functional of the model to the corresponding empirical version from the data, and (b) to use a likelihood based method. Both approaches have advantages and disadvantages, which will be illustrated below. Firstly, using least squares with the covariance function will be discussed. LEAST SQUARES: COVARIANCE FUNCTION Eqs. 1 and 2 suggest a method of estimating 0 using the empirical covariance function of the image. Firstly, the level u is estimated by solving the equation p ˆ=1-0(u ˆ), where p is the observed area fraction, u ˆ is the estimator of u and O denotes the standard normal distribution function. Denote the empirical (centered) covariance function of the image by Ru(t) - this can be obtained using morphological operations. Using Ru(t) instead of Ru(t) in Eq. 1 (or in Eq. 2), invert the integral numerically to obtain R(t), the estimated random field covariance. Least squares can be used to fit a nonlinear regression of either Ru(t; 6) toRu(t) orofR(t; 0) toR{t). Depending on which covariance function is chosen, it may be possible to convert the second of these into a linear regression. For instance, if the exponential covariance function given by Eq. 3 is used, then the problem can be converted to a linear regression of ln(—ln(R ˆ(t)) on ln(01)+ 02ln(||t||). Two problems can arise in this method. Firstly, due to the discretization of the image, the shape of the estimated covariance function near the origin can be distorted, thereby influencing adversely the fitted model. Secondly, depending upon the size and number of images available, the empirical covariance function goes negative, even though the theoretical covariance is strictly positive. This is due to sampling creating negative correlation with “holes” around connected 10 Wilson RJ et al: Excursion set models components. As well as not being able to take logs of the empirical covariance, it increases the slope and so biases the results. The best results seem to be obtained when the first few lags are ignored and only the lags up to about a third of the way to where the empirical covariance first goes negative are considered. It should be noted that the process of inversion using Eq. 1 does not preserve nonnegative definiteness of covariance functions, so that the resulting estimate of R(t) cannot be used directly for texture simulation. In fact, the estimated covariance matrix for the Gaussian vector corresponding to the observation lattice in the excursion set model can have a significant proportion of negative eigenvalues. LEAST SQUARES: SIZE DISTRIBUTIONS Other functionals, which capture important features of the data, can be used instead of the covariance function; for example, size distributions. If possible, an analytic form of the functional should be used; otherwise, Monte Carlo methods can be used to evaluate a numerical version for the parameter values considered in the search for fitted values. For size distributions for excursion sets of Gaussian random fields, the analytic form involves finding probabilities for the maximum of the random field over bounded sets. These can be derived for certain cases only. Numerical calculation of linear size distributions can be done efficiently. For excursion sets, these are related to first passage distributions for certain one dimensional processes and their Palm processes following an upcrossing. There are existing computational methods for rapidly calculating these first passage distributions, making linear size distributions very convenient for exploring excursion set geometry. Here, an example will be given to illustrate what can be achieved (see Nott and Wilson, 1997a for more details) - Fig. 3 has been reproduced with permission of Biometrics and Figs. 4-6 of World Scientific Publishing. Two types of linear size distribution can be defined: number weighted and measure weighted. Heuristically, the first describes the length of a “typical” linear intercept of a random set as a fixed line is traversed, while the second describes the length of an intercept containing a fixed point known to belong to a random set. Number weighted linear size distributions will be used in the example below. As the Gaussian field needs to be “smooth” enough to derive the distribution, it was assumed that the covariance function is exponential and 02 = 2. The example involves data from Diggle (1981), who fitted a Boolean model with random discs. The data (Fig. 3) describes the incidence of heather in Ja¨draa°s, Sweden. An area was divided into square plots and a binary mosaic was constructed by recording ones and zeros, depending on the incidence of heather in each plot. Estimating the level u from the area fraction of the right half gave uR = 0.112. Least squares was used to estimate 0, using lags equally spaced (four pixels apart due to correlations in the empirical distribution) between 8 and 32 and giving 6R = 0.067. Fig. 3. Heather incidence data. (Reproduced wit permission from Biometrics.) Fig. 4. Left half of data (a), simulation from fitted model (b). oo d CD d d CM d o d 10 20 30 Fig. 5. Fittedlinear size distr. for right half (solid) with estimated linear size distr’s for left half in horizontal (dotted line) and vertical (broken line) directions. 74 Image Anal Stereol 2001;20:71-78 Fig. 4 shows the left half of the data with a simulation from the model fitted to the right half, while Fig. 5 shows a plot of the fitted linear size distribution for the right half with the empirical linear size distributions for the left half. The fitted model seems to capture well the geometry of the observed images. Fig. 5 shows that the general shape of the size distribution for the fitted model matches that of the empirical size distribution, though there may be too few small intercepts and too many large intercepts in simulations from the fitted model. 10 20 30 40 50 Fig. 6. Estimated covariance from left half (solid line) and fitted covariance function (broken line) with upper and lower simulation envelopes obtained from 49 simulations from model fitted to right half (dotted). As a check on the adequacy of the model, the empirical covariance for the left half was plotted with the fitted covariance and simulation envelopes obtained from 49 simulations from the model fitted to the right half (Fig. 6). The simulation envelopes from the fitted model to the right half suggest the model is satisfactory. MAXIMUM LIKELIHOOD In theory, the joint distribution of the binary random field can be written down in terms of multiple integrals of the joint distribution of the Gaussian random field. However, this is obviously impracticable. One possible approach is to use the EM algorithm as described in Nott and Wilson (1997a). The EM algorithm (Dempster et al., 1977) is an iterative method for calculating maximum likelihood estimates in missing data problems. In situations where the EM algorithm can be applied, data is observed which is incomplete, while there is unobserved data which can be considered as complete and contains extra information. In the case here, the observed data is whether the random field is above or below a level, while the complete data is the random field itself. The EM principle is to maximise the expected log likelihood of the complete data given the incomplete data iteratively. Even so, the E step (finding the expectation of the log likelihood of the complete data given the incomplete data) is difficult to perform. One way to overcome this is by using Markov chain Monte Carlo methods. This is still time consuming, so the stochastic EM algorithm (in which only one simulation is used to approximate the expected log likelihood) is used in Nott and Wilson (1997b) to estimate the unknown parameters. The function evaluations needed to maximize the log likelihood are difficult for a high dimensional Gaussian vector. To make these more efficient, a modification inspired by methods for simulating random fields and involving embedding the random field values into a larger random vector whose likelihood is easier to calculate can be used. No statistical approximations are introduced by this embedding. See Nott and Wilson (1997b) for more background and details. Despite the various modifications above, the stochastic EM algorithm for excursion sets is still very slow. In the next section, an alternative approach, applied to excursion sets and the Boolean model in Nott and Ryde´n (1999), is described. It is both generally applicable and computationally attractive. PAIRWISE LIKELIHOOD Rather than using the full likelihood, it is possible to use the product of the bivariate or pairwise likelihood as an “approximation” to the full likelihood. This approach has been taken in Hjort and Omre (1994) and Lindsay (1988) for other situations. In practice, it is usually not necessary to take the product in the pairwise likelihood over all distinct pairs of observations, but only over pairs which show significant spatial dependence. Nott and Ryde´n (1999) show consistency and asymptotic normality under fairly general conditions for certain image models (including excursion set models and the Boolean model). Pairwise likelihood is computationally attractive for random set models, since bivariate likelihoods can be computed from the area fraction and random set covariance. Hence, whenever these quantities are easily computed or approximated from simulations then pairwise likelihood can be easily used. Simulation studies suggest that using pairwise likelihood is still reasonably efficient when compared with methods 0 h 75 Wilson RJ et al: Excursion set models based on the full likelihood (Nott and Ryde´n, 1999). An example using pairwise likelihood will be given in the following section. MULTI-PHASE IMAGES In Section “Estimation for two phase images”, various ways of fitting random set models have been discussed. Many images (like Fig. 1) consist of more than two phases and, in this section, methods of constructing and fitting three such models will be discussed. These multi-phase models have quite different properties. It is important to note these and then choose a model appropriate for the observed relationships between the phases. LABELLINGS OF RANDOM SETS In images where there is a “background” phase and the other phases have similar morphology and are not adjacent, it may be possible to assume that the union of the non-background phases constitute a single random set. The connected components of this random set can then be assumed to have been allocated to (or labelled as) the different phases. Estimation for the random set can be carried out using any of the methods of Section “Estimation for two phase images”. This idea of labelling components of a random set has been considered by Molchanov (1984) and others. There are various methods for labelling the connected components. One of the simplest is to assume that there is a collection of covariates (say particle size, perimeter, relationship to other particles, orientation) whose values determine the phase which is most likely for the particle under consideration; that is, the allocation is modelled by a polytomous logistic regression, conditional on the observed random set. Estimation for the parameters in the logistic regression can be carried out using standard methods (see, for example, Hosmer and Lemeshow, 1989). MODELS DEFINED FROM MULTIPLE RANDOM SETS In cases where there is a background phase and the remaining phases can be adjacent, it may be possible to assume that some phases dominate or cover other phases. In this situation, the phases are assumed to be ordered, with the background phase being the “bottom” phase and subsequent phases as covering the lower order phases. Hence the top phase can be observed completely, while the lower phases are masked. The simplest model is to assume the non-background phases are independent random sets, so each can be fitted as in the previous section, taking into account that some parts of the random set (that masked by higher order phases) are missing. These types of models seem to have been first used by Greco et al. (1979) as a model for sinter textures. CONTOURS OF RANDOM FIELDS The third possible model is obtained by considering contours of a random field at different levels. As with the previous model, the phases are again ordered, this time by an inclusion principle; that is, the background phase is again the bottom phase and each subsequent phase is contained within the lower order phases. In this situation, all phases can be observed fully and are not independent of each other. Estimation of the parameters is straightforward using pairwise likelihoods as the bivariate probabilities (which may be across phases and hence contours) simply involve bivariate normal probabilities. ANALYSIS OF A SULPHIDE ORE In this section, the methods described above will be used to fit a multi-phase model based on labellings to a core image of a sulphide ore. This and further examples can be found in Nott and Wilson (2000). (Figs. 1 and 7–10 have been reproduced with permission of Elsevier Science.) (a) (b) Fig. 7. First particle phase (a) and second particle phase (b). EXAMPLE The data of Fig. 1 is to be modelled with a model of the kind described in Section “Labellings of random sets”. Fig. 7 shows the two particle phases of the image. The union of these is to be modelled as an excursion set of a stationary Gaussian random field. The steps involved are as follows: 1. Calculate the empirical covariance function for the excursion set. 76 Image Anal Stereol 2001;20:71-78 2. Transform it to one for the random field via the functional relationship Eq. 1. 3. Choose an appropriate form for the covariance. 4. Estimate the level u from the area fraction (as in Section 3.1). 5. Use the pairwise likelihood to estimate the covariance function parameters. Fig. 8 shows the estimated covariance function for the random field covariance in the east-west and north-south directions. There seems very little evidence for anisotropy and so an isotropic rational quadratic covariance function was chosen. 00 Ö CO Ö d 0M d 0M d 10 20 30 40 50 Fig. 8. Empirical covariance function for random field: E-W(solid), N-S (dotted). The estimate for u was 1.4977. In the pairwise likelihood, the product was taken only over pairs of pixels which were either adjacent or separated by one pixel horizontally or vertically. Support for this choice is given in Nott and Ryden (1999). The estimates for 61 and 02 in the rational quadratic covariance function were 8.3519 and 1.26 respectively. The final stage is to use logistic regression to model the assignation of connected components to the two phases. In this example, the number of particles in the first phase was too few and no covariates were used; that is, each connected component was assumed to be assigned independently (equivalent to fitting only an intercept term in the logistic regression). The estimated probability of being assigned to the first phase was 0.097. The data is shown with three simulations from the fitted texture model in Fig. 9. Fig. 9. (a) Image of sulphide ore; (b)-(d) simulations from fitted model. Fig. 10. Empirical correlation and cross-correlation functions (solid), with means (dashed) and simulation envelope (dotted) from fitted model. (a) First phase: E-W; (b) First phase: N-S; (c) Second phase: E-W; (d) Second phase: N-S; (e) Cross-correlation between phases: E-W; (e) Cross-correlation between phases: N-S. 77 Wilson RJ et al: Excursion set models To check the adequacy of the model, the empirical correlation functions for the phases in the east-west and north-south directions, together with means and upper and lower simulation envelopes for these obtained from 49 simulations from the fitted model, were plotted. Similar plots for the cross-correlations between phases were also obtained. The results are shown in Fig. 10, and suggest that there is no real reason to question the adequacy of the model as the empirical functions lie well within the confidence envelopes. ACKNOWLEDGEMENTS We thank Dr Laurence Marcault who collected the core data while at the W. H. Bryan Mining Geology Research Centre at the University of Queensland. This work was part of a project directed by Dr G. J. Lyman for the C.R.C for Mining Technology and Equipment. REFERENCES Adler RJ (1982). The Geometry of Random Fields. New York: Wiley. Crame´r H, Leadbetter MR (1967). Stationary and Related Stochastic Processes. New York: Wiley. Cressie NAC (1993). Statistics for Spatial Data. New York: Wiley. Dempster AP, Laird NM, Rubin DB (1977). Maximum likelihood from incomplete data via the EM algorithm. J Roy Statist Soc B 39:1-38. Diggle PJ (1981). Binary mosaics and the spatial pattern of heather. Biometrics 37:531-9. Greco A, Jeulin D, Serra J (1979). The use of the texture analyser to study sinter structure. J Microscopy 116:199-211. Haindl M (1991). Texture synthesis. CWI Quarterly 4:305-31. Hjort NL, Omre H (1994). Topics in Spatial Statistics. Scand J Stat 21:289-357. Hosmer DW, Lemeshow S (1989). Applied Logistic Regression. New York: Wiley. Leadbetter MR, Hsing T (1990). Limit theorems for strongly mixing stationary random measures. Stoch Proc Appl 36:231-43. Lindsay B (1988). Composite likelihood methods. In: Prabhu NU, eds. Statistical Inference from Stochastic Processes, Contemp Math 80:221-39. Molchanov I (1984). Labelled random sets. Theor Probab Math Stat 27:113-9. Nott DJ, Ryde´n T (1999). Pairwise likelihood methods for inference in image models. Biometrika 86:661-76. Nott DJ, Wilson RJ (1997a). Size distributions for excursion sets. In Jeulin D, ed., Advances in Theory and Applications of Random Sets. Singapore: World Scientific, 175-96. Nott DJ, Wilson RJ (1997b). Parameter estimation for excursion set texture models. Signal Process 63:199-210. Nott DJ, Wilson RJ (2000). Multi-phase image modelling with excursion sets. Signal Process 80:125-39. Stoyan D, Kendall WS, Mecke J (1995). Stochastic Geometry and its Applications. New York: Wiley. Wilson RJ, Adler RJ (1982). The structure of Gaussian fields near a level crossing. Adv Appl Prob 14:543-65. Wilson RJ, Rychlik I (1996). Segmentation of mineral textures. Proceedings of Image Segmentation Workshop 96. Sydney, 93-6. 78