# Exact simulation of a Boolean model.

INTRODUCTION

The Boolean model is certainly one of the most currently used random set models in mathematical morphology, stochastic geometry and spatial statistics. It is defined as the union of a family of independent random compact subsets (denoted for short as "objects") located at the points of a locally finite Poisson process. It is stationary if the objects are identically distributed (up to their location) and the Poisson process is homogeneous, and non-stationary otherwise.

Despite its widespread use, it seems that little attention has been paid to the following problem: How to perform exact simulations of a Boolean model in a bounded domain?

The solution to that problem is not straightforward, unless the objects are uniformly bounded. The difficulty lies in that the intersection of a Boolean model and a bounded domain is also a Boolean model, but its parameters are different. The more remote the object to the domain, the less chance they have to hit it. on the other hand, the larger the objects, the more chance they have to hit the domain.

After a brief reminder on the Boolean model, this problem is investigated. Although most emphasis is placed on the stationary case because of its possible connections with stereology, the general case is also treated in a second part of the paper. Two examples serve to illustrate the algorithms and their implementations.

Here is the set of notation that will be used throughout the paper. Unless specified, the workspace is the d-dimensional Euclidean space [R.sup.d] with Lebesgue measure [v.sub.d]. The origin of [R.sup.d] is denoted by o. More generally, points are denoted by lower case letters, subsets by capital letters and family of subsets by calligraphic letters. If x [member of] [R.sup.d] and A [subset] [R.sup.d], [[tau].sub.x]A stands for the translation of A w.r.t. vector [??]x. The dilation of A by another subset B is defined as

[[delta].sub.B]A = {x [member of] [R.sup.d] : [[tau].sub.x]B [intersection] A [not equal to] [empty set]}.

DEFINITION AND MAIN PROPERTIES

The basic ingredients of a Boolean model are

(i) a Poisson point process P with an intensity function [theta] = ([[theta].sub.x], x [member of] [R.sup.d]) that is assumed to be locally integrable.

(ii) a family of nonempty and mutually independent compact random subsets ([A.sub.x], x [member of] [R.sup.d]). [A.sub.x] is called the object located at x. If it takes simple shapes, then its statistical properties can be specified by elementary descriptors (e.g., distribution of its radius for a disk). For more intricate shapes, the hitting functional of [A.sub.x] can be considered (Matheron, 1975). It assigns each compact subset K of [R.sup.d] (in short K [member of] K) the probability that it intersects with [A.sub.x],

[T.sub.x](K) = P{[A.sub.x] [intersection] K [not equal to] [empty set]}, K [member of] K. (1)

Definition 1 A Boolean model is a union of objects located at Poisson points,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Our objective is to simulate the Boolean model in a compact domain, say Z. Insofar as this mainly consists of simulating the objects of X that intersect with Z, it is crucial to assume their number to be almost surely finite. Under this assumption, it can be shown (Lantuejoul, 2002) that the number of objects hitting any compact subset K of Z is Poisson distributed with mean

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (2)

In particular, the avoiding functional of X [intersection] Z, defined for each K [member of] K(Z) as

[Q.sub.X[intersection]Z](K) = P{(X [intersection] Z) [intersection] K = [empty set]} = P{X [intersection] K = [empty set]}

is equal to

[Q.sub.X[intersection]Z](K) = exp(-v(K)) K [member of] K(Z). (3)

STATIONARY CASE

ALGORITHM

In this case, the intensity function is constant, say [theta], and all objects have the same hitting functional, up to a translation

[T.sub.x](K) = [T.sub.o]([[tau].sub.-x]K). (4)

In what follows, it is convenient to write

[T.sub.o](K) = [[integral].sub.K] dF(A)[1.sub.A[intersection]K[not equal to][empty set]],

where F (symbolically) denotes the distribution of the parameters of [A.sub.o]. Then we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Applying this formula to Eq. 2 and permuting integrals, the mean number of objects hitting Z becomes

v(Z) = [theta] [[integral].sub.K] dF(A)[v.sub.d]([[delta].sub.A]Z) = [theta]E{[v.sub.d]([[delta].sub.A]Z)), (5)

which gives the following expression for the avoiding functional of X [intersection] Z.

[Q.sub.X[intersection]Z](K) = exp(-[theta] [[integral].sub.K] dF(A)[v.sub.d]([[delta].sub.A]K)).

Let us write it slightly differently

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (6)

which involves a weighted version [F.sup.Z] of F

d[F.sup.Z](A) = dF(A)[v.sub.d]([[delta].sub.A]Z)/E{[v.sub.d]([[delta].sub.A]Z). (7)

An interpretation of Eq. (6) is as follows: X [intersection] Z is the union of a Poisson number (mean v(Z)) of independent objects of the form [[tau].sub.[??]]A where A is distributed like [F.sup.Z] and [??] is a uniform point over [[delta].sub.A]Z. Hence the simulation algorithm:

Algorithm 1

(i) set X = [empty set];

(ii) generate n ~ P(v(Z));

(iii) if n = 0, return X [intersection] Z and stop;

(iv) generate A ~ d[F.sup.Z];

(v) generate x ~ U ([[delta].sub.A]Z);

(vi) put X = X [union] [[tau].sub.x]A, n = n - 1 and goto (iii).

The main difficulty with this algorithm is step (iv): How to simulate the weighted distribution d[F.sup.Z]? The next section shows that interesting simplifications arise when the objects are convex.

CONVEX OBJECTS

Such algorithmic simplifications actually arise only when the simulations are performed within a ball-shaped domain. Accordingly, it is advantageous to firstly extend the simulation field to a ball, then perform the simulations in the extended domain, and finally restrict the simulations produced to the actual simulation field. The choice of the ball is unimportant, as long as it encloses the simulation field. One possibility is the ball circumscribed to the simulation field.

In what follows, the simulation field Z is assumed to be a ball, say B(o, [rho]). Then Steiner formula applies and shows that the mean number of objects hitting Z depends only on the expected Minkowski functionals of the objects:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (8)

Moreover, d[F.sup.Z] can be simulated as a mixture of distributions of objects weighted by their Minkowski functionals:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

with

d[F.sub.k](A) = dF(A)[W.sub.k](A)/E{[W.sub.k](A)} k = 0, ..., d, (9)

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (10)

Using this mixture of weighted distributions, the algorithm for simulating d[F.sup.Z] becomes

Algorithm 2

(i) generate k ~ p;

(ii) generate A ~ d[F.sub.k](A);

(iii) return A and stop.

The following two examples show how this algorithm can be implemented.

Example 1

The objects are balls and their radii follow independent exponential distributions with mean 1/a

P{R > r} = exp(-ar).

Starting from Eq. 5, the mean number of objects hitting Z is

v(Z) = [theta]E{[[omega].sub.d][([rho] + R).sup.d]} = [theta]d![[omega].sub.d]/[a.sup.d] [d.summation over (k=0)] [(a[rho]).sup.k]/k!,

where [[omega].sub.d] = [[pi].sup.d/2]/[GAMMA](d/2 + 1) is the d-volume of the unit ball.

The next step is the simulation of d[F.sub.k](A). Actually, it can also be written d[F.sub.k](r) because the only random element of A is its radius. If A = B(o, r), then [W.sub.k](A) = [[omega].sub.d][r.sup.d-k] and E{[W.sub.k](A)} = [[omega].sub.d](d - k)!/[a.sup.d-k]. Plugging these values into Eq. 9, we obtain

d[F.sub.k](r) = a exp(-ar)[(ar).sup.d-k]/(d - k)!.

This is a gamma distribution with parameter d - k + 1 and scale factor a. A simple way to simulate it is to consider -ln([u.sub.1] ... [u.sub.d-k+1])/a where [u.sub.1], ..., [u.sub.d-k+1] are independent uniform values on [0,1].

As an illustration, Fig. 1 shows a simulation of a Boolean model of discs with exponential radii. The displayed simulation field is a 7 x 5 rectangle. The Poisson intensity is [theta] = 10 and the parameter a of the exponential distribution has been chosen equal to 7.22 so as to provide a background proportion of 30%.

Example 2

Here is a somewhat more elaborate example, even if the Boolean model considered is only two-dimensional. The objects are typical Poisson polygons derived from a network of Poisson lines with intensity [lambda], see Fig. 2.

These polygons have been extensively studied in Miles (1969) and Matheron (1975). In particular their expected Minkowski functionals are

E{[W.sub.0](A)} = [1/[pi][[lambda].sup.2]] E{[W.sub.1](A)} = [1/[lambda]] E{[W.sub.2](A)} = [pi].

Using these expected values, Eq. (8) gives

v(Z) = [theta]/[pi][[lambda].sup.2] [(1 + [pi][lambda][rho]).sup.2].

Not so simple is the simulation of d[F.sup.Z](A) = [p.sub.0][F.sub.0](A) + [p.sub.1][F.sub.1](A) + [p.sub.2][F.sub.2](A). The explicit values for the weights are

[p.sub.0] = 1/[(1 + [pi][lambda][rho]).sup.2], [p.sub.1] = 2[pi][lambda][rho]/[(1 + [pi][lambda][rho]).sup.2],

[p.sub.2] = [[pi].sup.2][[lambda].sup.2][[rho].sup.2]/[(1 + [pi][lambda][rho]).sup.2].

With probability [p.sub.0], a polygon must be simulated from [F.sub.0]. The standard procedure to do this consists of generating Poisson lines sequentially by increasing distance from the origin. The procedure is continued until the generation of additional lines no longer affects the polygon containing the origin.

With probability [p.sub.1] a polygon must be simulated from [F.sub.1]. To do this, a polygon is first generated from [F.sub.0], and then split by a uniformly oriented line through the origin. It remains to select at random one of the two polygons thus delimited (see Fig. 3).

With probability [p.sub.2] a polygon must be simulated from [F.sub.2]. The algorithm proposed by Miles (1974) consists of taking the intersection between a polygon generated from [F.sub.0] and a cone delimited by two random rays emanating from the origin (see Fig. 3). Both rays are uniformly oriented and separated by an angle with p.d.f. f([phi]) = [phi] sin [phi]/[pi].

To illustrate the construction, a simulation of a Boolean model of Poisson polygons is depicted in Fig. 4. With a Poisson intensity [theta] and a Poisson line intensity [lambda] respectively set to 10 and 1.625, the background proportion is close to 30%.

NON STATIONARY CASE

ALGORITHM

Let us start again with Eq. 3 that provides the avoiding functional of X [intersection] Z, and let us write it like

[Q.sub.X[intersection]Z](K) = exp(-v(Z) [v(K)/v(Z)]),

provided that v(Z) > 0 (if v(Z) = 0, then X [intersection] Z = [empty set] a.s.). On the other hand, Eq. 2 can be used to derive the following expansion

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

which shows that v(K)/v(Z) is the hitting functional of some object of Z, say [A.sup.Z]. More precisely, [A.sup.Z] is located according to the pdf

f(x) = [[theta].sub.x][T.sub.x](Z)/v(Z), x [member of] [R.sup.d]. (11)

Moreover, the conditional hitting functional of [A.sup.Z] given its location x is equal to

[T.sub.x](K)/[T.sub.x](Z) = P{([A.sub.x][intersection]Z) [intersection] K [not equal to] [empty set]| [A.sub.x][intersection]Z [intersection] Z [not equal to] [empty set]. (12)

[A.sup.Z] is called a typical object hitting Z.

The interest of typical objects lies in the following property (Lantuejoul, 2002):

Proposition 1 X [intersection] Z is the union of a Poisson number (mean v(Z)) of typical objects hitting Z.

Indeed, we can readily check that the avoiding functional of such a union of typical objects coincides with that of X [intersection] Z:

[[infinity].summation over (n=0)] exp(-v(Z)) [[v.sup.n](Z)/n!] [(1 - [v(K)/v(Z)]).sup.n] = exp(-v(K)).

From this proposition, the following algorithm is derived for simulating a Boolean model, even non stationary, in the domain Z.

Algorithm 3

(i) set X = [empty set];

(ii) generate n ~ P(v(Z));

(iii) if n = 0, return X and stop;

(iv) generate a typical object [A.sup.Z] hitting Z;

(v) put X = X [union] [A.sup.Z], n = n - 1 and goto (iii).

PRACTICAL IMPLEMENTATION

Algorithm 3 calls for several remarks. Step (ii) assumes that v(Z) is explicitly known, but this is not always the case because the integral of [[theta].sub.x][T.sub.x](Z) may not be mathematically tractable. Moreover, step (iv) requires simulating the pdf (Eq. 11) that specifies the location of the typical objects hitting Z. This distribution may have a complicated expression. Step (iv) also involves generating conditionally typical objects given their locations. It is possible to generate objects [A.sub.x] sequentially till one is produced that hits Z. This algorithm can be easily implemented but lacks efficiency.

In order to alleviate these difficulties, rejection and coupling methods can be resorted to. A typical approach is to consider another Boolean model X' that dominates X, in the sense that its ingredients satisfy the following three properties:

(i) its Poisson intensity function [theta]' satisfies [[theta]'.sub.x] [greater than or equal to] [[theta].sub.x] (at each point x)

(ii) its population of objects ([A'.sub.x]) satisfies [A'.sub.x] [contains] [A.sub.x]

(iii) there exists an algorithm to conditionally simulate [A.sub.x] given [A'.sub.x].

Then the idea is firstly to generate the objects of X' hitting Z. Then each generated object [A'.sub.x] is saved with probability [[theta].sub.x]/[[theta]'.sub.x] and discarded with the complementary probability 1 - [[theta].sub.x]/[[theta]'.sub.x]. Finally, each remaining object [A'.sub.x] is replaced by a new object [A.sub.x] generated using the conditional algorithm mentioned in (iii). A simulation of X in Z is given by the union of these new objects in Z.

Hence, a Boolean model can effectively be generated whenever it is dominated by another Boolean model that is numerically tractable and can be simulated. There is no general rule for building such a dominating model. This must be done on a case by case basis.

doi: 10.5566/ias.v32.p101-105

ACKNOWLEDGEMENTS

The topic of this paper was presented at the S4G Conference, June 25-28, 2012 in Prague, Czech Republic. The author thanks P. Calka for fruitful discussions about the simulation of weighted Poisson polygons.

REFERENCES

Lantuejoul C (2002). Geostatistical simulations: Models and algorithms. Berlin: Springer.

Matheron G (1975). Random sets and integral geometry. New York: Wiley.

Miles RE (1969). Poisson flats in Euclidean spaces. Adv Appl Prob 1:211-37.

Miles RE (1974). A synopsis of Poisson flats in Euclidean spaces. In: Harding EF, Kendall DG, eds. Stochastic geometry. New York: Wiley. pp. 202-27.

CHRISTIAN LANTUEJOUL ([mail])

MinesParisTech, 35 rue Saint-Honore, 77305 Fontainebleau, France

e-mail: christian.lantuejoul@mines-paristech.fr

(Received February 27, 2013; revised April 4, 2013; accepted May 5, 2013)
COPYRIGHT 2013 Slovenian Society for Stereology and Quantitative Image Analysis
No portion of this article can be reproduced without the express written permission from the copyright holder.