Image Anal Stereol 2013;32:101-105 doi:10.5566/ias.v32.p101-105 Original Research Paper EXACT SIMULATION OF A BOOLEAN MODEL Christian Lantuejoulb 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) ABSTRACT A Boolean model is a union of independent objects (compact random subsets) located at Poisson points. Two algorithms are proposed for simulating a Boolean model with non uniformly bounded objects in a bounded domain. The first one applies only to stationary models. It generates the objects prior to their Poisson locations. Two examples illustrate its applicability. The second algorithm applies to stationary and non-stationary models. It generates the Poisson points prior to the objects. Its practical difficulties of implementation are discussed. Both algorithms are based on importance sampling techniques, and the generated objects are weighted. Keywords: Boolean model, importance sampling, Minkowsky functionals, Steiner formula. 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 Rd with Lebesgue measure vd. The origin of Rd 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 e Rd and A c Rd, txA stands for the translation of A w.r.t. vector ox. The dilation of A by another subset B is defined as sBA = {X e Rd : tXB n A = 0} . DEFINITION AND MAIN PROPERTIES The basic ingredients of a Boolean model are (i) a Poisson point process P with an intensity function d = (dx,x e Rd) that is assumed to be locally integrable. (ii) a family of nonempty and mutually independent compact random subsets (Ax, x e Rd). Ax 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 Ax can be considered (Matheron, 1975). It assigns each compact subset K of Rd (in short K e K) the probability that it intersects with Ax, TX(K )= P{AX n K = 0} K e K. (1) Definition 1 A Boolean model is a union of objects located at Poisson points, U Ax xeP 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 0 (K) = ( 0X Tx (K) dx , K e K (Z) . (2) In particular, the avoiding functional of X n Z, defined for each K e K (Z) as Qx nZ (K) = P{(X n Z) n k = 0} = p{x n k = 0} is equal to Qx nZ (K) = exp(—0 (K)) K e K (Z) . (3) STATIONARY CASE ALGORITHM In this case, the intensity function is constant, say e, and all objects have the same hitting functional, up to a translation (4) Tx(K) = To(T-xK) . In what follows, it is convenient to write To(K) = I dF(A) 1AnK=0 , JK where F (symbolically) denotes the distribution of the parameters of Ao. Then we have Tx(K) = [ dF (A) 1TxAnK=0 . K Applying this formula to Eq. 2 and permuting integrals, the mean number of objects hitting Z becomes 0(Z) = e [ dF(A)vd(SaZ) = eE{vd(SaZ)}, (5) K which gives the following expression for the avoiding functional of X Z. Qxnz(K) = exp ^-0 jf^dF(A)vd {SaK) Let us write it slightly differently Qx nz (K)= exp (-, (z) f d.z (A) V^) , which involves a weighted version Fz of F dF (A)vd (Saz) dFz (A) E {vd (Saz )} (7) An interpretation of Eq. (6) is as follows: X n Z is the union of a Poisson number (mean 0(Z)) of independent objects of the form txA where A is distributed like FZ and x is a uniform point over SAZ. Hence the simulation algorithm: Algorithm 1 (i) setX = 0; (ii) generate n ~ P (0(Z)); (iii) if n = 0, return X n Z and stop; (iv) generate A ~ dFZ; (v) generate x ~ U (SAZ); (vi) put X = X U txA, n = n — 1 and goto (iii). The main difficulty with this algorithm is step (iv): How to simulate the weighted distribution dFZ? 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, p). Then Steiner formula applies and shows that the mean number of objects hitting Z depends only on the expected Minkowski functionals of the objects: *(z) = d £if)E{Wk(A)}pk k=0 (8) Moreover, dFz can be simulated as a mixture of distributions of objects weighted by their Minkowski functionals: yd id dFz (A) = dF (A)_ 7k=°}k )Wk (A)pk (6) yd=o (d)E {Wk (A)}pk y Pk dFk(A) , k=0 with and dFk(A) dF (A)Wk(A) E {Wk (A)} pk (dk)E{Wk (A)}p k EU (d)E{Wi(A)}p1 k = 0,..., d, (9) k = 0,..., d. (10) Using this mixture of weighted distributions, the algorithm for simulating dFZ becomes Algorithm 2 (i) generate k ~ p; (ii) generate A ~ dFk(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 0 (z ) = e e{ ®d (p + r)d} ed!®d Á (ap)k ad L " k=0 k! where rnd = nd/2/r(d/2 + 1) is the d-volume of the unit ball. The next step is the simulation of dFk(A). Actually, it can also be written dFk(r) because the only random element of A is its radius. If A = B(o, r), then Wk (A) = rndrd-k and E{Wk(A)} = ad(d — k)!/ad—k. Plugging these values into Eq. 9, we obtain dFk (r) = a exp(-ar)(ar) (d -k)! dk This is a gamma distribution with parameter d — k + 1 and scale factor a. A simple way to simulate it is to consider — ln(ui • • • ud—k+\)/a where ui,...,ud—k+\ 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 d = 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%. Fig. 1. Simulation of a Boolean model of discs with exponential radii. 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 A, see Fig. 2. Fig. 2. Realization of a Poisson line process. The polygons delimited by the lines are Poisson polygons. These polygons have been extensively studied in Miles (1969) and Matheron (1975). In particular their expected Minkowski functionals are E{Wo(A)} = E{Wx(A)} = A E{W2(A)} = n . Using these expected values, Eq. (8) gives d 0 *(Z) = ^ (1 + nAp)2 . Not so simple is the simulation of dFZ(A) = p0F0(A) + piFi (A) + p2F2(A). The explicit values for the weights are 1 2nAp po = (1 + nAp)2 , p1 = (1 + nAp)2 , _ n2A2p2 p2 = (1 + nAp)2 . With probability p0, a polygon must be simulated from F0. 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 p1 a polygon must be simulated from F\. To do this, a polygon is first generated from F0, 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 p2 a polygon must be simulated from F2. The algorithm proposed by Miles (1974) consists of taking the intersection between a polygon generated from Fo 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 (y) = y sin y/n. Fig. 3. Weighted polygons generated from distributions F0 (left), F1 (middle) and F2 (right). To illustrate the construction, a simulation of a Boolean model of Poisson polygons is depicted in Fig. 4. With a Poisson intensity d and a Poisson line intensity A respectively set to 10 and 1.625, the background proportion is close to 30%. provided that û (Z) > 0 (if û (Z) = 0, then X n Z = 0 a.s.). On the other hand, Eq. 2 can be used to derive the following expansion û (K) QXTX(Z) Tx(K) û(Z) JRd û(Z) Tx(Z) dx, which shows that $ (K)/$ (Z) is the hitting functional of some object of Z, say AZ. More precisely, AZ is located according to the pdf f (x) = PxTx(Z) û (Z) x dx (at each point x) (ii) its population of objects (A'x) satisfies A'x D Ax (iii) there exists an algorithm to conditionally simulate Ax given Ax. Then the idea is firstly to generate the objects of X' hitting Z. Then each generated object A'x is saved with probability Qx/Q'x and discarded with the complementary probability 1 — Qx/Qx. Finally, each remaining object A'x is replaced by a new object Ax 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. 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 ApplProb 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.