Confidence Regions for Means of Random Sets using Oriented Distance Functions

# Confidence Regions for Means of Random Sets using Oriented Distance Functions

Hanna K. Jankowski
York University
Larissa I. Stanberry
University of Washington

February 25, 2011
###### Abstract

Image analysis frequently deals with shape estimation and image reconstruction. The objects of interest in these problems may be thought of as random sets, and one is interested in finding a representative, or expected, set. We consider a definition of set expectation using oriented distance functions and study the properties of the associated empirical set. Conditions are given such that the empirical average is consistent, and a method to calculate a confidence region for the expected set is introduced. The proposed method is applied to both real and simulated data examples.

Keywords: Random Closed Set, Simultaneous Confidence Interval, Image Reconstruction, Shape Estimation.

## 1 Introduction

Image analysis frequently deals with the problem of shape estimation. Many instances of this may be found, for example, in medical imaging, where shape analysis of brain structures is often used to differentiate between different populations of subjects. Examples include the study of the hippocampus for schizophrenic patients and corpus callosum for adults with fetal alcohol exposure, as well as other neuroanatomical structures (Styner et al., 2004; Bookstein et al., 2002b, a; Styner et al., 2003; Levitt et al., 2009).

To analyze the observed shapes, it is natural to think of them as realizations of random sets. The problem of shape inference then translates into finding the average set and describing its variability.

As the space of closed sets is nonlinear, there is no natural way to define the mean of a set. Indeed, many different definitions of the expected set exist, and therefore one must first select a definition of the mean to work with. Quoting Molchanov (2005), “the definition of the expectation depends on what the objective is, which features of random sets are important to average, and which are possible to neglect”. Here, we focus on the definition of set expectation given in Jankowski and Stanberry (2010), which is based on the oriented distance function (ODF) 111To differentiate it from others, we will refer to the Jankowski and Stanberry (2010) definition as the ODF definition. However, the definition Baddeley and Molchanov (1998) may also be based on the ODF.. The definition is akin to those considered in Baddeley and Molchanov (1998); Lewis et al. (1999), which also rely on the distance function. The ODF definition, however, has several desirable theoretical properties and was found to outperform other definitions in image applications; for a detailed comparison see Jankowski and Stanberry (2010). A thorough review of other existing definitions appears in Molchanov (2005).

In this paper, we study the empirical estimators of the expected set and expected boundary using the ODF definition. We give conditions for consistency of these estimators. We also study their variability using the concept of confidence regions.

Seri and Choirat (2004) present two methods for the calculation of the confidence region of the Aumann expectation. To our best knowledge, this is the only other time confidence regions were studied in the context of random sets. The methodology was developed for convex sets, and therefore the Aumann expectation is a natural choice, as it always yields a convex answer.

Molchanov (1998, Theorem 3.1) gives a central limit theorem for the Hausdorff distance of sublevel sets. Baddeley and Molchanov (1998) use this result to obtain a central limit theorem for their expectation, and this may also be done for the ODF definition. Such central limit theorems could potentially be used to calculate confidence regions. The regions may be found as a dilation of the empirical set estimator. However, this approach requires the estimation of a complex functional of the derivative of the expected distance transform, which renders the method impractical. Furthermore, a dilation approach would provide a uniform confidence region: the distance between the boundary of the mean set and the boundary of the confidence set would be equivalent at all points. Thus, the dilation method would mask important information about the local variability of the estimators.

In this paper, we propose a new and simple approach to calculate confidence regions for both the mean set and its boundary. The method works for both convex and non-convex sets. The resulting confidence regions are conservative in that they cover the expected set with at least probability. We show that the confidence regions satisfy certain natural equivariance properties, which are analogous to those of confidence intervals for real–valued parameters. Moreover, the confidence regions provide a simple visual representation of the variability of the estimators and are able to detect local changes in variability. The method can also be used in Bayesian inference to study the behaviour of the posterior sample.

Our definition of the confidence region is based on the quantiles of the supremum of a Gaussian random field. We consider several examples where this quantile is calculated easily. When these quantiles are not available analytically, we propose a bootstrap method to provide approximate confidence regions. The bootstrap approach also avoids making an assumption of an underlying distribution of the observed sets, which could be quite difficult for practitioners.

The outline of this paper is as follows. In Section 2, we review the definitions of set and boundary expectations given in Jankowski and Stanberry (2010). Consistency is studied in Section 3 and the confidence regions are described in Section 4. Section 5 gives several examples of our approach, including a simulation study. We consider the toy image reconstruction example discussed by Baddeley and Molchanov (1998), and analyse the sand grains data first discussed in Stoyan and Molchanov (1997). Lastly, we consider a medical imaging example, and apply our methods to a boundary reconstruction problem in a mammogram image. Proofs and technical details appear in the Appendix, which is available online as supplementary material.

### 1.1 Notation and Assumptions

Throughout, we let denote the domain on which the sets are observed. Unless otherwise stated, we assume that is the working domain and write, for example, without stating that explicitly. We also assume that is a subset of , and denote the Euclidean norm of as .

We write for the closed ball of radius centered at . For a set , we write and to denote its interior, closure, complement and boundary. Unless noted otherwise, set operations are calculated relative to the domain . That is, and so forth. Deterministic sets are denoted using capital letters , while bold upper-case lettering, , is used for random sets.

The notation is used to denote the space of continuous functions endowed with the uniform topology. That is, in if , for all compact subsets We write to say that converges weakly to . Throughout the paper, when handling weak convergence of stochastic processes or random fields, we assume that these take values in .

## 2 Random Closed Set and Its Expectation

Let be the family of closed sets of and let denote the family of all compact subsets of . For a probability triple , a random closed set (RCS) is the mapping such that for every compact set

 {ω:A(ω)∩K≠∅}∈A.

For more background on random closed sets, we refer to Matheron (1975); Ayala et al. (1991).

### 2.1 Definition of Expectation

For a nonempty set the distance function is defined as for all Given the distance function of a closed set , the original set may be recovered via . Also, the Hausdorff distance may be calculated using the distance function as

for any sets . We refer to Delfour and Zolésio (2001) for more mathematical properties of the distance function. The distance function has also been used in the context of image thresholding, see for example Friel and Molchanov (1999); Molchanov and Terán (2003).

The oriented distance function (ODF) of any set such that is defined as for all For a closed set, the set and its boundary may now be recovered by and . Note that, given only the information at a fixed point the ODF is more informative than the distance function. Given the value of we know the value of but the converse statement is not true.

For an RCS with a.s. non-empty boundary we define the random function , and denote its pointwise mean as . The mean set and mean boundary are then defined as follows.

###### Definition 2.1.

Suppose that is a random closed set such that almost surely and assume that for some . Then

 E[A] = {x:E[bA(x)]≤0}, Γ[A] = {x:E[bA(x)]=0}.

Furthermore, we define the expectation of the complement as

###### Remark 2.2.

If is closed and then is open, however, the oriented distance function of continues to be well–defined, and indeed we have that . Thus, .

###### Example 1 (disc with random radius).

Suppose that , for some non–negative, integrable, real-valued random variable . Then and hence and . That is, the expected set is a disc with radius , with boundary equal to the expected boundary.

###### Example 2 (random singleton).

Suppose that for some -valued random variable , then

In Example 1, Definition 2.1 yields a natural answer, while the result of Example 2 seems counterintuitive. As mentioned previously, the choice of expectation depends on the problem at hand. The ODF definition was motivated by problems in imaging, where random singletons are regarded as noise. Example 2 illustrates a natural de-noising property of the expectation; see also the discussion in Jankowski and Stanberry (2010).

The function provides additional information about the mean and its boundary. Recall that is Lipschitz with constant one almost surely (Delfour and Zolésio, 1994). Therefore, is also Lipschitz, which also implies that and are closed sets. Furthermore, if the function is smooth in a neighbourhood of , then is also smooth near

###### Proposition 2.3.

Suppose that and fix , and suppose that the function is in a neighbourhood of , and that Then, there exists a (possibly different) neighbourhood of such that is .

We next consider the estimation of and

###### Definition 2.4.

Suppose that we observe random sets, and let

 ¯bn(x)=1nn∑i=1bAi(x).

Then the empirical mean set and the empirical mean boundary are defined as

 ¯An={x:¯bn(x)≤0}  and  ¯¯¯¯Γn={x:¯bn(x)=0}.

Similarly to the mean ODF, the function is Lipschitz, and therefore the empirical sets and are both closed.

### 2.2 Consistency and Fluctuations of the Empirical ODF

Suppose that we have independent and identically distributed (IID) RCSs . Then the following results are valid for the average ODF.

###### Theorem 2.5.

Suppose that are IID, and that for some , Then for all compact subsets

 limnsupx∈K|¯bn(x)−E[bA(x)]|=0,

almost surely.

###### Theorem 2.6.

Suppose that are IID, and that for some . Then

 Zn(x)≡√n(¯bn(x)−E[bA(x)])⇒Z(x),

where is a mean zero Gaussian random field with covariance

 cov(Z(x),Z(y))=E[bA(x)bA(y)]−E[bA(x)]E[bA(y)]

for

The next result shows that has very smooth sample paths. Recall that a function is Hölder of order if it satisfies for some positive finite constant and , for all in the domain of .

###### Proposition 2.7.

For any

 var(Z(x)−Z(y)) ≤ |x−y|2, (2.1) |cov(Z(y)−Z(x),Z(y′)−Z(x′))| ≤ 2|y−x||y′−x′|. (2.2)

Moreover, the sample paths of are Hölder of order , for any .

## 3 Consistency

Here, we study consistency of the estimators and , assuming that almost surely in By Theorem 2.5, these results apply to IID random closed sets. Following Molchanov (1998), we say that converges strongly to , if the Hausdorff distance almost surely, for any compact set . The following theorem follows from Molchanov (1998, Theorem 2.1) and Cuevas et al. (2006, Theorem 1).

###### Theorem 3.1.

Suppose that

 limnsupx∈D|¯bn(x)−E[bA(x)]|=0

almost surely, and that is well-defined. Suppose also that the expected ODF, , satisfies

 {x:E[bA(x)]≤0} = ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯{x:E[bA(x)]<0}. (3.1)

Then converges strongly to . If also satisfies

 {x:E[bA(x)]≥0} = ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯{x:E[bA(x)]>0}, (3.2)

then converges strongly to .

Condition (3.1) says that the expected ODF is not allowed to have a local minimum on , while (3.2) excludes mean ODFs which have a local maximum on (again, this need not be unique). Alternatively, since is a continuous, condition (3.1) says that is a topologically regular closed set, while condition (3.2) says that is a topologically regular open set.

###### Remark 3.2.

It is possible that violates the conditions (3.1) and/or (3.2) and consistency still holds. For example, consider the RCS almost surely. Then IID sampling trivially produces a consistent estimate, but fails to satisfy (3.1).

###### Example 3 (half plane).

For , consider the RCS , where is a real-valued random variable with finite mean Then , and The mean ODF satisfies both conditions (3.1) and (3.2), and therefore and are consistent estimators of and . Indeed, we may easily check that which converges to zero almost surely.

The following result provides some further insight into the consistency conditions.

###### Proposition 3.3.

Conditions (3.1) and (3.2) may be re–written in terms of mean set properties.

1. Condition (3.1) holds iff and iff

2. Condition (3.2) holds iff , and iff

###### Example 4 (set and its boundary).

Suppose that is either or with equal probability. Then , while . On the other hand, if is seen with probability , then if , we have that ; If , then and The case provides a setting where neither (3.1) nor (3.2) are satisfied.

We observe independent sets , where each is either or with equal probability. Let denote the proportion of the random sets equal to . Then

 ¯bn(x)=ˆpnb[0,1](x)+(1−ˆpn)b{0,1}(x),

and it follows that whenever , , and for , while . Clearly, convergence to the expected set or the expected boundary can never be achieved.

###### Example 5 (missing centre).

Suppose that is either a disc or an annulus in ; that is,

 A = {{x:|x|≤1}with probability % p,{x:0.5≤|x|≤1}otherwise.

Then the expected set is an annulus for , and a disc for . For , the expected ODF satisfies both (3.1) and (3.2). For we have

 E[bA(x)] = {|x|−1for |x|≥0.75,−|x|/3otherwise,

and hence while . Since has a local maximum at , it fails to satisfy (3.2), and therefore this point may be omitted by the estimators .

Suppose that the RCS is either a rectangle or a union of two squares with equal probability. Specifically, define

 A1 = {x:0≤x1≤3,0≤x2≤1}, A2 = {x:0≤x1≤1,0≤x2≤1}∪{x:2≤x1≤3,0≤x2≤1}.

Then with probability 0.5, and otherwise Thus, half of the time the has its “middle” removed. The resulting mean set and mean boundary are shown in Figure 1. Here, the expected ODF has no local maxima, or minima, along the boundary, and therefore both (3.1) and (3.2) are satisfied.

## 4 Confidence Regions

We now construct confidence regions, or confidence supersets, for both and To do this, we assume that the sets are observed on a compact window . We also assume that in which holds, for example, under IID sampling by Theorem 2.6.

###### Definition 4.1.

Let and be numbers such that , and . Then, a confidence region for is

 {x∈W:¯bn(x)≤q1/√n} (4.1)

and a confidence region for is

 {x∈W:|¯bn(x)|≤q2/√n}. (4.2)

By Proposition 2.7, the limiting field is continuous, and therefore both and are well-defined. Further, Proposition 2.7 gives a uniform upper bound on the variability of the increments of the random field. Understanding the path properties of the process , such as smoothness, provides information about the variability of the quantiles and and therefore also on the tightness of the confidence sets. We also make the following comments about the new definition.

• The confidence region is conservative, in that it covers the set or with a probability of at least One reason why the method is conservative is our use of the supremum of the fluctuation field to find the cut–off quantile values. However, the field is very smooth and highly correlated by Proposition 2.7. Therefore, we expect that the proposed method, although conservative, yields reasonable answers, especially for sets that have been co–registered apriori. We explore the question of over–coverage via simulations in Section 5.2.

• The confidence region is “immune” to the consistency conditions (3.1) and (3.2); see Example 8.

• If the exact distribution of is unknown (as is often the case in practice), the quantiles may be approximated using a bootstrap approach. For computational reasons, it is often easier to calculate asymmetric quantiles in (4.2).

• The confidence region gives no information as to the geometry of the random sets. That is, the shape of the confidence region may be similar for observed random sets which are stars, as well as for those which are squares.

###### Example 7 (disc in R2 with random centre).

The random set is a disc with radius one centred at where is Uniform, and suppose that we observe 100 IID random sets from this model. The expected set is shown in Figure 2. Moreover, since where , it follows that is contained inside the disc of radius one centered at . Confidence regions were formed for both and using re-sampling techniques to estimate the quantiles of and where the window These are illustrated in Figure 2.

###### Example 8 (confidence regions for the set and its boundary).

Let and suppose that is either or with equal probability. Suppose also that we observe a simple random sample of size from this model. Recall that , and satisfies neither (3.1) nor (3.2). As before, let denote the proportion of times that the set is observed. If , then . If , then with

The fluctuation field is given by

 Zn(x) = √n(ˆpn−0.5)(b[0,1](x)−b{0,1}(x)) ⇒ Z(b[0,1](x)−b{0,1}(x)),

where is a univariate normal random variable with mean zero and variance . The largest difference for occurs at , and hence, for any window such that , we have

 supx∈WZ(x)=max{−Z,0},   and %  supx∈W|Z(x)|=|Z|.

Therefore, the exact quantiles are and , and the confidence region for is given by Now, for any , . Therefore the confidence region misses a part of if and only if . For large , this happens with a probability of roughly . On the other hand, the Hausdorff distance whenever This example illustrates that the confidence regions are not affected by the violation of the consistency conditions.

###### Example 9 (confidence set for disc with random radius).

Suppose that is a disc with random radius with and . Then the expected set is a circle with radius . Also, the 95% confidence interval for is a circle with radius , while the 95% confidence set for is the band .

###### Example 10 (pacman in R2).

Define the pacman with radius , , to be a disc with radius centred at the origin with its upper left quadrant removed. That is,

 A(r)={x:|x|≤r}∩{{x:x1≤0}∪{x:x2≤0}}.

Figure 3 (right, dashed) shows the contour of for . Suppose that where is a uniform random variable on

The expected set is a smoothed version of , as seen in Figure 3 (left). Recall that the smoothness of the boundary depends on the smoothness of . As the latter is an integral, which tends to have more smoothness than the original function, we expect that in general the expectation is as or more smooth than the original realizations of the boundary. This is exactly what one sees in Figure 3: the mean boundary is smoothed out in the regions where there is movement; however, since the origin is a fixed point of the random boundary, no additional smoothness is introduced here by the average.

Figure 3 also shows bootstrapped confidence sets for both and The accuracy of the estimate and the apparent centering of the confidence intervals around the mean set may be explained upon closer inspection of the sample. The ODF of the pacman is similar to that of the circle; indeed, they are identical in the lower left quadrant. Therefore, the behaviour of the estimators and of the confidence regions is not unlike that of estimators and confidence intervals for the real-valued . For our sample of , we observed which explains the accuracy of the estimator . On the other hand, the confidence region still shows the large variability of .

RCS with simplified ODF structure.

Consider an RCS and suppose that there exist functions and random variables such that

 bA(x)=k∑j=1hj(x)ηj

for all (see also Section 3 in Jankowski and Stanberry (2010)). For example, the ball centered at with random radius takes this form. Here, , and hence and

For such sets, both the empirical and ODF average have a similar simplified form. That is, and Furthermore, the fluctuation field for these sets is particularly straightforward. Suppose that we observe IID samples of the random vector , and assume that for all . Then

 Zn(x)⇒k∑j=1Zjhj(x),

where is a multivariate normal random variable with mean zero and variance matrix given by

###### Example 11 (Random half-plane).

Suppose that Uniform and let That is, is the plane above the line which goes through the origin and has angle with the -axis. Here, and hence and Some calculations also reveal that and

The confidence region for is a strip centred on the line of width where is the quantile of

### 4.1 Equivariance Properties

The following proposition gives equivariance properties of the confidence regions under dilation and rigid motion. The result corresponds to the classical scaling results for the mean and standard deviation of univariate data.

###### Proposition 4.2.

Consider a random closed set Let and denote the confidence regions for and respectively.

1. Fix , and let with Then the confidence regions for and are and respectively.

2. Fix a rigid motion , and let with Then the confidence regions for and are and respectively.

### 4.2 A Modified Approach

It is not hard to see that the variability of and depends on the local fluctuations of the field around (see also Molchanov (1998, Theorem 3.1)). One natural way to incorporate this idea into calculation of the confidence region is described below.

1. Calculate the confidence regions for as described in the previous section. Let denote this region.

2. Find the upper quantile of , and call this .

3. The confidence region for is then calculated as

The modification is designed to decrease the size of the quantile in (4.1). In the first step we reduce the size of the domain of the supremum to a set which is likely to contain the boundary. Note also that the set , and therefore also depend on the sample size . Thus, the larger the sample size, the more effective the modification. A similar approach yields a modified confidence region for We find that this modification yields a slight improvement in the coverage probabilities for some settings, although visually the confidence regions remain quite similar. Notably, iterating (i)–(iii) does not improve the size of the region or the coverage probabilities.

## 5 Examples

### 5.1 Implementation

There exist several efficient algorithms to calculate the distance function of any set, which allows for easy implementation of our methods (Breu et al., 1995; Freidman et al., 1977; Rosenfeld and Pfaltz, 1966). For many of the examples presented here, the oriented distance function may be calculated exactly. When this was not possible, our calculations were implemented in MATLAB (MathWorks), where the bwdist command was used to compute the oriented distance function of a set. Here, the examples need to be discretized to pixels (in ) or voxels (in ). This discretization introduces an additional source of error; see Serra (1984) for a thorough treatment of the induced difficulties.

To minimize the effect of discretization, in the simulations described below, we selected a fine gird, which was calibrated to give accurate results. Suppose that , and that we observe discs centered at the origin with random radius Uniform[0,1]. Here, the confidence regions for the mean boundary are exact (modulo the sample size approximations). Setting , we obtained empirical coverage probabilities of 95.10, 95.02, 95.30, 95.30, and 94.96 for the 95% confidence regions, for grids with side length , respectively (the standard error due to bootstrap sampling was 0.0031). Finally, we selected for our simulations.

### 5.2 Simulation Study

We next simulate coverage probabilities for several examples of random closed set models. The particular examples we consider are given below.

1. The set and its boundary when , considered in Example 4. Recall that in this example netiher the set estimator nor the boundary estimator is consistent. Here

2. The pacman with random radius Uniform[0,2] (see Example 10). Here Although the pacman RCS is not decomposable, it still exhibits similar behaviour.

3. The random set where and are both random discs. The two cases we consider are as follows.

1. is centered at the origin and has radius Uniform[0,2]. is centered at the point and has radius Uniform[0,1]. Here

2. is centered at the origin and has radius Uniform[0,1]. is centered at the point and has radius Uniform[0,1]. Here The mean set and some sample sets are shown in Figure 4.

4. A random ellipse with boundary parameterized as . Let be two independent random variables with distribution Uniform The two cases we consider are as follows. Throughout, we assume that

1. .

2. .

In the above examples, the image was discretized and the quantiles were estimated using Monte Carlo methods (with repetitions and ), except for example (A), where the quantile is known exactly. The empirical coverage was estimated using Monte Carlo simulations in each case. It is important to note that there are four sources of error: sample size, discretization, the bootstrap estimation of the quantile, and the Monte Carlo simulations themselves. The maximal standard error due to Monte Carlo of the empirical coverage probabilities is , but it is not easy to quantify the standard error due to the other three sources. In particular, we note that the quantile was estimated once per example, and the same quantile was used in each Monte Carlo simulation (otherwise, the simulations would have been prohibitive), which increases the bias of our results.

The results of the simulation study show that our methods, though conservative, achieve good coverage probabilities. The greatest over-coverage is seen in Example (C2), which is the most difficult of the examples because the observed sets have not been co-registered; see Figure 4 for some sample sets. In this example, it is possible that a slight improvement of the coverage probability could be seen by a modification to the confidence region discussed in Section 4.2. In general, to minimise over-coverage, we recommend choosing as small as possible in practice.

### 5.3 A Toy Image Reconstruction Example

Image averaging arises in various situations, for example, when multiple images of the same scene are observed or when the acquired images represent objects of the same class and the goal is to determine the average object (shape) that can be described as typical. Here we consider an example of image averaging studied in Baddeley and Molchanov (1998). The data set consists of 15 independent samples of a reconstructed newspaper image (Figure 5, left), and is available from http://school.maths.uwa.edu.au/homepages/adrian. For details on how the data was generated we refer to Baddeley and Molchanov (1998).

The empirical average of the 15 observed images, is shown in Figure 5 (right). The average describes the “typical” image reconstruction, and as such, may be thought of as an estimator of the true text image. Next, we compute 95% confidence regions for the ODF-average reconstruction based on 5K bootstrap samples. We may think of these regions as a measure of the variability of , and . Figure 6 (top) shows the confidence region for with the boundary of the true image overlayed in black. The confidence set contains all of the true text image and appears tight, although there are a number of spurious bounds induced by noise. The confidence set for is also shown in Figure 6 (right). The boundary of the lower confidence set consists of only a few closed contours scattered throughout the text. The lack of tightness in the confidence set is explained by the small sample size and a relatively thin font width. Hypothetically increasing the sample size would produce tighter confidence sets for both the expected set and its boundary. For example, the bootstrap confidence set for the boundary based on 50 and 100 samples is tighter as compared to the one based on 15 samples (see Figure 7). We expect that actual confidence regions for and would exhibit even less noise than the hypothetical regions shown in Figure 7. It should be noted that the confidence regions in Figure 7 were created using the original window , and the picture is a close-up of the result.

Note that these images were generated under a Bayesian framework, where the goal is to reconstruct the original image by selecting an appropriate parameter from the posterior distribution. Here, this parameter is not computable directly, and is instead estimated by and . Although a frequentist concept, the confidence regions allow us to determine the variability of these estimators. As Figure 7 shows, this variability depends on the sample size.

### 5.4 Analysis of Sand Grains

Next, we apply the proposed methods to the sand grains data previously described in Stoyan and Molchanov (1997) and Kent et al. (2000). The sand particles were collected from the shores of the Baltic Sea and banks of the Zelenchuk River in Ossetia. The grains were photographed on the same scale and the data resembles two dimensional projections represented by binary images. Both images may be found online at www.math.yorku.ca
/hkj/Research/SandGrains
(please note that the online images are not to scale). The observed sand grains generally have a smooth rounded shape, but are not necessarily covex.

The grains from the two regions appear to differ both in shape and size with sea grains being more spherical and larger as compared to river grains that are more oblong and smaller in diameter. To summarize the data, we find the average grains and their average boundaries for each group. We also construct confidence regions to describe the variability of the grain shapes.

To begin with, the particles were realigned using the generalised Procrustes analysis as implemented in the shapes package in R (R Development Core Team, 2009). To apply the Procrustes analysis, we use the digitized data as described by Kent et al. (2000). The data was digitized so that each sand grain was represented by 50 vertices approximately equally spaced on the boundary. After digitization, the arc-length between the vertices is around 10-20 pixels with grain particles represented by high-resolution images of size . The realignment was done with and without scaling. Using the scaling, we essentially remove the size effect and can examine differences in average shapes. Alternatively, average particles based on Procrustes analysis without scaling reflect differences both in size and shapes of the particles. The median (IQR) centroid sizes are 1481 (1396, 1665) and 2076 (1867, 2376) for the river and sea grains, respectively, indicating the sea particles to be bigger as compared to the river ones.

Figure 8 shows the confidence regions for the average particle (top row) and the expected boundary (bottom row) for the river (left column) and the sea (right column) sands using scaled realignment. White contours show the empirical mean boundary. The confidence regions are based on 5K bootstrap samples. The average river grain is more oblong as compared to the sea grain. The variability within the two groups appears to be rather similar. Figure 9 shows the confidence regions for the average particle (top row) and the expected boundary (bottom row) for the river (left column) and the sea (right column) sands, based on realignment without scaling. White contours show the empirical mean boundary. The confidence regions are based on 5K bootstrap samples. The average sea grain is more spherical in shape as compared to the river average. It is also considerably larger in size. The variability within the two groups again appears to be rather similar, however, the unscaled images appear more variable than the scaled ones. Overall, the discrepancies in shape and size between the two averages reflect the differences between the raw data sets. Note that the boundaries of the average sets in both figures are rather smooth.

There is also a marked difference between the scaled and unscaled river sand grain averages. Image results from the Procrustes analysis re–alignment with and without re–scaling are quite different. This is to be expected, as the scaling, location, and centering re–alignments in Procrustes analysis are highly interdependent. It would be of interest to compare other re–alignment methods, such as those proposed by Stoyan and Molchanov (1997), but this is beyond the scope of this work.

Whether or not scaled Procrustes realignment is applied, the averages of the particles show a clear difference between the two groups. However, this difference is at this time only visual. An interesting and important problem is to develop quantifiable methodology to test for presence and locations of differences between the mean shapes.

### 5.5 Application to Medical Imaging

We next consider an example of boundary reconstruction in mammography, where the skin-air contour is used to determine the radiographic density of the tissue and to estimate breast asymmetry. Both measures are known to be associated with the risk of developing breast cancer (Scutt et al., 1997; Ding et al., 2008). In Stanberry and Besag (2009), B-spline curves were used to reconstruct a smooth connected boundary of an object in a noisy image, and a Bayesian approach was applied to estimate the tissue boundary in mammograms.

The boundary reconstruction was performed on a binary image which was obtained after filtering and thresholding the original greyscale mammogram. Let be the compact domain of the mammogram image, and let denote the random set describing the breast tissue, or foreground, of the image under the prior belief distribution. Next, let denote the noisy binary mammogram image observed. The skin-air boundary estimate is reconstructed as , the expected boundary of from the posterior distribution given the observed data . The posterior distribution of the random set is too difficult to compute, and was approximated via Markov chain Monte Carlo. Hence, the skin-air boundary estimate was also approximated as from a sub-sample of observed random sets generated from the MCMC simulation. Further details can be found in Stanberry and Besag (2009). Here, we apply the proposed method to construct a confidence region for the posterior mean boundary. We emphasize that the confidence region is not a credible set, but rather it describes the variability of as an estimator of

Figure 10 (left) shows a typical digitized mammogram image, characterised by a low contrast-to-noise ratio. A probability integral transform improves the contrast by increasing the dynamic range of image intensities (Figure 10, centre). The solid white line in Figure 10 (centre and insets on right) shows the reconstructed boundary Note that what appears to be a nipple is, in fact, a duct system leading to the nipple, so that the estimators correctly follow the skin line.

The 95% confidence set (dashed) for the true boundary in Figure 10 is obtained using a bootstrap resampling of size 1000. The confidence set is tight and fits the image well. It also shows that the reconstructed boundary is more variable toward the inside of the breast tissue. The background of the black and white mammogram image has considerably more noise than the foreground. Consequently, the posterior boundary samples show more variability toward the inside of the tissue. More details can be seen in insets in Figure 10 (right).

Note also that to apply our methods, we have assumed that the observed ODFs are independent and identically distributed, whilst the boundary reconstruction is based on Markov chain Monte Carlo sampling from the posterior. To ensure the independence of the curve samples, we construct the confidence set for the boundary using 100 samples from the posterior, which were acquired every 250th sweep after a burn-in period of 1000 sweeps.

## 6 Discussion

In this paper, we studied consistency of set and boundary averages under random sampling. We also presented a method for the construction of confidence regions for the mean set and mean boundary. The confidence regions, though conservative, have appealing equivariance properties and are straightforward to implement. Simulations indicate that they achieve good coverage probabilities. Unlike previous developments, our methods are applicable to both convex and non-convex sets and allow for differences in local variability.

As there exists no notion of the standard deviation of a random set, one can also use the confidence regions as an informal assessment of the variability of the mean set estimator. This relationship is strengthened by the aforementioned equivariance properties, which mimic the scaling properties of the standard deviation and confidence region in the univariate setting.

In Section 5 we considered several empirical examples. The observed sets in these cases are non-convex, and therefore methods based on the Aumann expectation would not work well, although they may yield reasonable approximations for the sand grains example. In Section 5.5 we applied the proposed methods to a boundary reconstruction problem in a mammogram image. There, the confidence region technique was used to assess the variability of the MCMC estimation of the posterior mean. The proposed method was able to detect an increase in the variability of the sampling toward the inside of the breast tissue. A dilation method, such as one based on the results of Molchanov (1998), would not be able to detect this difference.

## Acknowledgements

The first author would like to thank Tom Salisbury from York University for several useful discussions. Both authors thank Ian Dryden from the University of South Carolina and Dietrich Stoyan from the Freiberg University of Mining and Technology for the sand grains data and for advice on the Procrusters transformations.

## References

• Ayala et al. (1991) Ayala, G., Ferrandiz, J. and Montes, F. (1991). Random set and coverage measure. Adv. Appl. Prob. 23 972–974.
• Baddeley and Molchanov (1998) Baddeley, A. and Molchanov, I. (1998). Averaging of random sets based on their distance functions. J. Math. Imaging Vision 8 79–92.
• Bookstein et al. (2002a) Bookstein, F. L., Sampson, P. D., Connor, P. D. and Streissguth, A. P. (2002a). Midline corpus callosum is a neuroanatomical focus of fetal alcohol damage. The Anatomical record 269 162–174.
• Bookstein et al. (2002b) Bookstein, F. L., Streissguth, A. P., Sampson, P. D., Connor, P. D. and Barr, H. M. (2002b). Corpus callosum shape and neuropsychological deficits in adult males with heavy fetal alcohol exposure. NeuroImage 15 233–251.
• Breu et al. (1995) Breu, H., Gil, J., Kirkpatrick, D. and Werman, M. (1995). Linear time euclidean distance transform algorithms. IEEE Trans. Pattern Anal. Mach. Intell. 17 529–533.
• Cuevas et al. (2006) Cuevas, A., González-Manteiga, W. and Rodríguez-Casal, A. (2006). Plug-in estimation of general level sets. Aust. N. Z. J. Stat. 48 7–19.
• Delfour and Zolésio (1994) Delfour, M. C. and Zolésio, J.-P. (1994). Shape analysis via oriented distance functions. J. Funct. Anal. 123 129–201.
• Delfour and Zolésio (2001) Delfour, M. C. and Zolésio, J.-P. (2001). Shapes and geometries, vol. 4 of Advances in Design and Control. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA. Analysis, differential calculus, and optimization.
• Ding et al. (2008) Ding, J., Warren, R., Warsi, I., Day, N., Thompson, D., Brady, M., Tromans, C., Highnam, R. and Easton, D. (2008). Evaluating the effectiveness of using standard mammogram form to predict breast cancer risk: Case-control study. Cancer Epidemiology Biomarkers and Prevention 17 1074–1081.
• Freidman et al. (1977) Freidman, J., Bentley, J. and Finkel, R. (1977). An algorithm for finding best matches in logrithmic expected time. ACM Trans. Pattern Anal. Softw. 3 209–226.
• Friel and Molchanov (1999) Friel, N. and Molchanov, I. (1999). A new thresholding technique based on random sets. Pattern Recognition 32 1507–1517.
• Jankowski and Stanberry (2010) Jankowski, H. and Stanberry, L. (2010). Expectations of random sets and their boundaries using oriented distance functions. J. Math. Imaging Vision 36 291–303.
• Kent et al. (2000) Kent, J. T., Dryden, I. L. and Anderson, C. R. (2000). Using circulant symmetry to model featureless objects. Biometrika 87 527–544.
• Kunita (1990) Kunita, H. (1990). Stochastic flows and stochastic differential equations, vol. 24 of Cambridge Studies in Advanced Mathematics. Cambridge University Press, Cambridge.
• Levitt et al. (2009) Levitt, J., Styner, M., M., N., S., B., Koo, M., Voglmaier, M., Dickey, C., Niznikiewicz, M., Kikinis, R., McCarley, R. and Shenton, M. (2009). Shape abnormalities of caudate nucleus in schizotypal personality disorder. Schizophr Res 110 127–139.
• Lewis et al. (1999) Lewis, T., Owens, R. and Baddeley, A. (1999). Averaging feature maps. Pattern Recognition 32 1615–1630.
• Matheron (1975) Matheron, G. (1975). Random sets and integral geometry. John Wiley & Sons, New York-London-Sydney. With a foreword by Geoffrey S. Watson, Wiley Series in Probability and Mathematical Statistics.
• Molchanov (2005) Molchanov, I. (2005). Theory of random sets. Probability and its Applications (New York), Springer-Verlag London Ltd., London.
• Molchanov and Terán (2003) Molchanov, I. and Terán, P. (2003). Distance transforms for real-valued functions. J.Math.Anal.Appl. 278 472–484.
• Molchanov (1998) Molchanov, I. S. (1998). A limit theorem for solutions of inequalities. Scandinavian Journal of Statistics 25 235–242.
• R Development Core Team (2009) R Development Core Team (2009). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0.
URL http://www.R-project.org
• Rosenfeld and Pfaltz (1966) Rosenfeld, A. and Pfaltz, J. (1966). Sequential operations in digital picture processing. J. ACM 13 471–494.
• Schwartz (1981) Schwartz, L. (1981). Cours d’analyse. 1. 2nd ed. Hermann, Paris.
• Scutt et al. (1997) Scutt, D., Manning, J., Whitehouse, G., Leinster, S. and Massey, C. (1997). The relationship between breast asymmetry, breast size and the occurrence of breast cancer. Br. J. Radiol 70 1017–1021.
• Seri and Choirat (2004) Seri, R. and Choirat, C. (2004). Confidence sets for the Aumann mean of a random closed set. In Lecture Notes in Computer Science, vol. 3045, Proceedings of ICCSA 2004, International Conference in Computational Science and its Applications. 298–307.
• Serra (1984) Serra, J. (1984). Image analysis and mathematical morphology. Academic Press Inc. [Harcourt Brace Jovanovich Publishers], London. English version revised by Noel Cressie.
• Spivak (1965) Spivak, M. (1965). Calculus on manifolds. A modern approach to classical theorems of advanced calculus. W. A. Benjamin, Inc., New York-Amsterdam.
• Stanberry and Besag (2009) Stanberry, L. and Besag, J. (2009). Boundary reconstruction in binary images using splines. Preprint.
• Stoyan and Molchanov (1997) Stoyan, D. and Molchanov, I. S. (1997). Set-valued means of random particles. J. Math. Imaging Vision 7 111–121.
• Styner et al. (2003) Styner, M., Gerig, G., Lieberman, J., Jones, D. and Weinberger, D. (2003). Statistical shape analysis of neuroanatomical structures based on medial models. Medical Image Analysis 7 207 – 220. Functional Imaging and Modeling of the Heart.
• Styner et al. (2004) Styner, M., Lieberman, J. A., Pantazis, D. and Gerig, G. (2004). Boundary and medial shape analysis of the hippocampus in schizophrenia. Medical Image Analysis 8 197–203.
• Winkler (1964) Winkler, W. (1964). Stetigkeitseigenschaften der Realisierungen Gauss’scher zufälliger Felder. In Trans. Third Prague Conf. Information Theory, Statist. Decision Functions, Random Processes (Liblice, 1962). Publ. House Czech. Acad. Sci., Prague, 831–839.

## Appendix

Recall that the boundary of a set is in a neighbourhood if there exists a bijective map , which is and whose inverse is also , which maps the boundary into the set That is, the boundary is at if locally it is a manifold.

###### Proof of Proposition 2.3.

We first note that the condition is necessary. For example, if has a local minimum at then the boundary of contains the isolated point , and no smoothness properties may be carried from the expected oriented distance function to the expected boundary.

Write with Since , there exists a such that Let and define by The boundary is now described via the set and we may apply the implicit function theorem (Theorem 1-12, page 41 of Spivak (1965) and Theorem 31 p. 299 of Schwartz (1981)). The differentiability condition on the Jacobian in the implicit function theorem holds since

 ∂xjH(x(j),xj)=∂jE[bA(x0)]≠0.

It follows that there exists a function , a neighbourhood of , , and a neighbourhood of , such that is and describes the boundary, , near . ∎

###### Proof of Theorem 2.5.

Pointwise convergence follows immediately by the law of large numbers. Since both and are Lipschitz functions (Jankowski and Stanberry, 2010), we also obtain uniform convergence over compact sets. ∎

The next result may be found, for example, in Kunita (1990, Theorem 1.4.7). We repeat it here for convenience.

###### Theorem F.1 (Kolmogorov’s tightness criterion.).

For a compact set , let be a sequence of continuous random fields with values in . Assume that there exist positive constants and with such that

 E[|Yn(x)−Yn(y)|γ] ≤ C(d∑i=1|xi−yi|αi)  % for all x,y∈D, E[|Yn(x)|γ] ≤ C,    for all x∈D,

holds for any . Then is tight in .

###### Lemma F.2.

There exists a constant , depending only on , such that

 E[|Zn(x)−Zn(y)|2d] ≤ C(d)|x−y|2d,

for any and

###### Proof.

The case is immediate. Next, consider ,

 E[|Zn(x)−Zn(y)|4] = n−2n∑i,j,k,l=1E[b∗ib∗