Image Anal Stereol 2010;29:121-131
Original Research Paper
SECOND MOMENT MEASURE AND K-FUNCTION FOR PLANAR STIT TESSELLATIONS
Viola Weiss1, Joachim Ohser2 and Werner Nagel3
1 Fachhochschule Jena, FB Grundlagenwissenschaften, Carl-Zeiss-Promenade 2, D-07745 Jena, Germany; 2Hochschule Darmstadt, FB MN, SchOfferstraße 3, D-64295 Darmstadt, Germany; 3 Institut fiir Stochastik, Friedrich-Schiller-Universitat Jena, D-07737 Jena, Germany e-mail: nagel@minet.uni-jena.de
(Accepted May 12, 2010)
ABSTRACT
For STIT tessellations - stationary tessellations that are stable under the operation iteration of tessellations -the second-order measure of the edge system is studied. A result is that this measure coincides with that one of a Boolean segment process. In the isotropic case an explicit formula for the pair-correlation function is given. An estimator for the covariance function of the edge length measure is derived and adapted to digitized images of tessellations. For m pixels of an image the algorithm is of complexity O (m log m).
Keywords: covariance function, estimation, K-function, pair-correlation function, random tessellations, second moment measure, STIT tessellations, stochastic geometry .
INTRODUCTION
The model of STIT tessellations - stationary tessellations that are stable under the operation iteration of tessellations - was introduced in Nagel and Weiß (2005). This stochastic stability is an essential property which also allows to derive many theoretical results. In particular, formulae were published for mean values (first-order moments) and also for certain distributions which concern nodes, edges, and cells. A next step is to find quantitative expressions which describe the mutual arrangement of cells. This points to the study of second-order entities. In the present paper we study the second moment measure for the edge system of planar random STIT tessellations. This can also be expressed by the K-function, the pair-correlation function or the covariance function, respectively. An explicit formula for the pair-correlation function in the isotropic case is derived. Finally, an asymptotically unbiased estimator for the covariance function is given.
The initial problem was whether this K-function reflects the structure of the tessellation, in particular the mutual arrangement of the cells. The result of the present paper shows that the K-function of a STIT tessellation coincides with the K-function of a Boolean model of segments with appropriately chosen parameters. This means that in this case the K-function even does not express whether the random segments form a tessellation or not.
In stochastic geometry second-order quantities, ^.e., the second moment measure, the K-function or the pair-correlation function respectively, were
considered to represent essential information about a random geometric structure, in particular the arrangement of geometric objects. But already Baddeley and Silverman (1984) provided examples of point processes with the same K-function but with evidently different point patterns - thus showing that the K-function does not necessarily comprise sufficient information. our result contributes an example from the class of segment processes and tessellations.
Firstly we recall the definitions of the second moment measure, the K-function and the pair-correlation function for random segment processes in the plane, and a theorem is cited which relates the second moment measure to marked section point processes when the random segment process is intersected with a line. The marks of the section points are the section angles. Then the STIT tessellations are explained and their basic properties are summarized. This is followed by a proof of the theorem, that the second moment measures of STIT tessellations and of Boolean segment processes coincide. This yields the K-function and the pair-correlation function for the segment system of STIT tessellations. In the last section, an asymptotically unbiased estimator for the covariance function of the edge system is derived. It is based on a convolution with a kernel function (i.e., a "smoothing") of the edge system. This approach is not specific to STIT tessellations and can be applied to arbitrary fibre processes with existing pair-correlation function. Finally, it is adapted to digitized images of tessellations.
SECOND MOMENT MEASURE AND K-FUNCTION FOR PLANAR FIBRE PROCESSES
Let R denote the set of real numbers, R^ the Euchdean plane and the (7-algebra of Borel sets in R^. By £2, ii, we denote the Lebesgue measures on R^, R and [0,oo), respectively.
Further, we use the notations Y for a planar fibre system and for a random fibre process in the sense of Stoyan et al. (1995), Ry denotes the distribution of
A fibre process is said to be macroscopically homogeneous - i.e., spatially stationary - if is invariant under all translations of R^. It is isotropic if its distribution is invariant under rotations. The total fibre length of the fibre system i/a in 5 G is written as For a stationary fibre process the mean total fibre length per unit area (length intensity) is La.
The second moment measure of a fibre process is described by
= lxi/{Bi)-xi/{B2)Ri>{dxi/),
{orBi,B2e^2-
If is stationary, the reduced second moment measure J^ can be defined by
xB2) = LI j J lB,{x)lB,{x+h) J^(dA)dx,
where 15 denotes the indicator fianction of B.
In the case of a stationary and isotropic fibre process it is sufficient to consider the reduced second moment function Kyy (called K-fianction) given by
= J^(5(o,r)), for r>0,
where r) is the circle with radius reentered in o.
The product Z,/iÄ'vj/(r) can be interpreted as the mean total length of fibres of in a circle with radius r, when the center of the circle is located in the typical fibre point. Such a typical point can be understood as a randomly chosen point inside a bounded window, and its distribution is the uniform distribution concentrated on the fibre system.
The stationary fibre process is called a second-order process if LaJ(^{B) < »o for all bounded B G 9^2- In the following we only consider second-order processes. And fiarthermore, we assume that also the intersection point processes n ^on all lines ^are of second order.
If the K-fianction is differentiable, the pair-correlation fianction ^exists with
är) =
1 di('vp(r) iTir dr
for r > 0 .
We consider two particular cases of fibre systems.
a) Let *F be a stationary and isotropic Poisson line process with length intensity La. Then the K-fianction is given by cf Stoyan et af (1995)
LAK^{r) = 2r+ Laki^ for r > 0 and the pair-correlation fianction is
g{r) = l +
LATir
for r>0.
b) For line segments we choose the parametrization {x,a,s) where x G R^ is the centre, a G [0,:/r) the normal direction (i.e., the angle between the normal of the segment and the positive x-axis) and 5 G (0,°°) the length. Thus a segment process *F can be considered as a point process on the space R^ x [0, n) x (0,°°). Notice that the normal direction will be used in the parametrization of both the segments and the lines.
Now let *F be a stationary and isotropic Boolean segment process with intensity N\ of the centre point (germ) process and with length distribution ^ of the typical segment. Then the length intensity is La = N\s where 5 denotes the mean length of the typical segment. For the K-fianction we obtain {cf Stoyan et af, 1995)
\o
LAK^{r) = š-'
00
+J (Irs-ids)
Laki^. (1)
/
The formulae for LAK^{r) for both particular cases can be interpreted as follows. The first summand arises from the segment on which the typical point (the center of the circle) lies. The second one is formed from the remainder of the process which has again the distribution due to the independence properties of the models.
A STEREOLOGICAL FORMULA FOR THE SECOND MOMENT MEASURE OF PLANAR FIBRE PROCESSES
There are at least three stereological methods for determination of second-order quantities of planar fibre processes, see Weiß and Nagel (1994). In the present paper we make use of one of these methods which was originally presented by Schwandtke (1988), where the fibre process is intersected by a line and the intersection process with the corresponding intersection angles is observed.
Denote by the set of all lines in R^ and by 11/- z| | the distance between two points / and zin R^. Further, let d^ be the element of the motion invariant measure on the set of all lines with fl{gnB(o, l)^@}dg=27i.
For a fibre system y/ (here considered as a subset of R^) and ^^ we consider i/Afl^ which consists of isolated intersection points and of linear segments of positive length of i/A on ^
xi/ng=S{xi/ng)uL{xi/ng),
where	denotes the set of all intersection points
and L{xi/rig) the union of all straight line segments on g. It is clear, that only for countably many g we have
An intersection point 7G S{xi/rig) is marked by the fibre tangent angle w{y,g). This is the angle between the tangent of Xj/ in y and g. If the tangent is not uniquely defined, then put w{y,g) = 0. Notice that the intersection angle between the tangent and g is the same as the angle between the respective normal directions of the fibre and of the line.
Assume that any line ^ G ^^ is parametrized as {y{t): t eR} suchthat	|| = ki - ^2! for all
ti, /2 G R- This yields a parametrization of the union of segments	as Lg= {teR: y{t) G L{\^r\g)}.
Then the following stereological formula for the second moment measure holds. Let be a stationary planar fibre process. Then for all measurable nonnegative fianctions / on R^ x R^ we have

= Ll
I
y,zeS{\ifr\g) y+z
J Jf{yy+z)J(r^{dz)dy \\y-4
sm usm V
f{yz) dgPi>id\i/)
+ / X [ [ f{y{ti),y{t2))dtidt2Ri>{dxif), (2) where u= w{yg), v= w(z,g).
SECOND MOMENT MEASURE AND K-FUNCTION FOR PLANAR STIT TESSELLATIONS
A new model for random tessellations, the so-called STIT tessellations, was introduced in Nagel and Weiß (2005). The simulation of a planar STIT tessellation in Fig. 1 suggests that STIT tessellations are potential models for crack or fissure structures.
Fig. 1. Simulation of an isotropic STIT tessellation.
The characteristic and eponymous property is the stability of their distribution under the operation iteration.
We will give here only a short overview, for more details see the cited papers.
Construction: A construction of STIT tessellations in bounded windows was described in all details in Nagel and Weiß (2005) and a global - i.e., in the whole plane R^ - construction was given in Mecke et al. (2008a). This construction can be understood as a process of sequential cell division at random times.
Let A = LaIi X ^ be a measure on R x [0, n) and ^ is a probability measure on [0,:/r), the directional distribution. In order to generate a tessellation it is assumed that ^ is not concentrated on a single direction. A point {p, a) G R x [0, n) represents a line g{p,a) G where p is its signed distance from the origin (positive iff the intersection of the line with its perpendicular through the origin lies in the upper halfplane) and a is the angle between its normal and the positive x-axis. For a set 4 C R^ denote [4] = {(p, a) G
At time ? = 0 the construction starts with a compact convex polygonal window C R^ {e-g-, a rectangle). After a random life-time ti that is exponentially distributed with parameter A([I4^) a random line ji is thrown onto W with the distribution	- n
[W\). Thus two new polygons are bom, and the birth time t\ is attributed to them. Then, sequentially, all the extant polygons	say, with the respective
birth times /j,..., are divided in the following way. The life-time of W- is a random variable that is exponentially distributed with parameter A([I^']), i.e., the parameter depends on the size of W-. In the particular case that M is the uniform distribution on [0,:7r) it is proportional to the perimeter of W^. At the end of its life-time the polygon is divided by a chord that comes from a random line that has the distribution A([I^'])-iA(- n [I^']). Thus W; dies and two new polygons are bom. The state at a fixed time a > 0 is a tessellation in W and it is denoted by
A cmcial feature of the constmction is that the parameters of the exponentially distributed life-time depend on the (random) size of the relevant cells, such that smaller cells have a longer expected life than larger ones. This implies a certain dependence between the random lifetimes of different cells. Therefore we used the following rejection method for a formal description of the constmction: Let J= 1,2,..., be a sequence of independent and identically distributed (i.i.d.) random variables, where also Tj and jj are independent. The Tj are exponentially distributed positive numbers with parameter A([W^). The jj are random lines with the distribution	n [I^). Any of these
random variables can be used at most once during the constmction. The division of an extant polygon
C is performed according to the rejection mle:
For a pair (tj, Jj) throw - after the time Tj elapsed -the line Jj onto the window W. If it hits W- then use the generated chord to divide this polygon. If /jfi W- =% then reject jj and make a new trial with another pair {tj^Yj). Also in the case of rejection Tj has to be added to the lifetime of Wf. This yields the lifetime distribution of the cells described above.
The chords - if they do not end on the boundary of W - that appear during the constmction are called I-segments. Later occurring I-segments have their endpoints in the relative interior of two already extant I-segments or chords, respectively.
Even if this constmction is related to a fixed and bounded window W it yields a distribution that is consistent in the following sense.
Existence: There exists a stationary tessellation O of the whole R^ such that
D
0(1,140 = on l^.
(3)
D
the symbol = stands for the identity of distributions of random variables.
The distribution of the tessellation O does not depend on W, and one can show that this formula holds for all compact and convex 2-dimensional sets W C R^. Moreover, the existence of a stationary tessellation is sufticient that the intensities and the K-fianction are well defined.
Generalizations to higher dimensions are given in Nagel and Weiß (2005); Mecke etal. (2008b).
Now we recall some important properties of STIT tessellations, the proofs were given in earlier papers (Nagel and Weiß, 2003; 2005; 2006; Mecke et al, 2007; Nagel and Weiß, 2008; Mecke et al, 2010). A random tessellation Y is understood as the set of all its edge points, i.e., as a random subset of R^.
The operation of iteration (also referred to as nesting) for tessellations is defined as follows. Denote by = {Y\,Y2,...} a sequence of independent and identically distributed (i.i.d.) stationary random tessellations. Further assume that Jo is a stationary random tessellation which is independent of . It is usefial to consider the set C{Yq) = {pi, /32, • • •} of the cells (which are convex polytopes) of Yq . The iteration of the tessellation Yq and the sequence is defined as
=	(4)
k>\
This definition means that a cell pk of the so called 'frame' tessellation Jo is - independently of all other cells - subdivided by the cells pi^, i= 1,2,... of the tessellation which intersect the interior of pk-The result of an iteration of stationary tessellations
is a stationary tessellation. A strict formalization and a proof of homogeneity are given in Mecke et al. (2008b).
According to the general concept of stochastic stability w.r.t. an operation, a random tessellation Yq is called stable w.r.t. iteration (STIT) if the distribution of the iterated tessellation, multiplied by a rescaling factor, is the same as that one of Yq. This is now defined more precisely. For a real number r > 0 the tessellation rY is generated by transforming all points {x,y) G Y into {rx,ry). Accordingly, iW means that this transformation is applied to all tessellations of the sequence '3^. Let Jo be a stationary random tessellation and	■ ■ ^ sequence of sequences of tessellations
such that all the occurring tessellations (including Yo) are i.i.d. Then the sequence /2(}o),-6(Jo),• • • of rescaled iterations is defined in Nagel and Weiß (2003; 2005)as
l2{Yo) = I{2Yo,2fy,),

m) = /(/(3	= Il^-^m), 3^2 ] ,
\
UYo) = I
m
m-l

/
in = 3,A,...
Here m is the rescaling factor which is chosen such that the results of iteration do not degenerate for
Deßnition: A stationary random tessellation Y is said to be stable with respect to iteration (STIT) if
Y=I^{Y) forallm = 2,3,... ,
i.e., if its distribution is not changed by repeated rescaled iteration with sequences of tessellations with the same distribution.
It can be shown (see Mecke et al, 2010) that, equivalents, Y is STIT if 7 = /2 (7).
Stability with respect to iteration: The stationary tessellation O given in (3) is stable with respect to iteration.
Characteristic entities: The segments (edges of cells) in the tessellation O form a stationary segment process with La (which appeared in the definition of A) as the mean total edge length per unit area and the directional distribution i.e., the distribution of the normal direction in a randomly chosen point on the segments. By La and M the distribution of stationary STIT tessellations is already uniquely determined.
Poisson typical cell: Now we consider the interior of the typical cell of O. That means more intuitively the single isolated cell neglecting additional nodes or edges emanating outside on their boundary. One can
show, that the interior of the typical cell of O with directional distribution ^ and intensity La has the same distribution as the interior of the typical cell of a stationary Poisson line tessellation with the same M and La, see Nagel and Weiß (2003).
Length distribution of the typical I-segment: For a stationary and isotropic STIT tessellation O with intensity La the density of the length of the typical I-segment is for x > 0


Llx^
2n
n
2 \
X LaX^ 12 3
A^-J
(5)
see Mecke et al. (2007).
This length distribution is a mixture of exponential distributions
p{x) = J -se 0
-isx
2s
d5.
Results for the non-isotropic case are derived in Mecke (2009). Further distributions are given in Mecke et al. (2007; 2010).
Mean values of STIT tessellation are calculated in Nagel and Weiß (2006; 2008). It can be shown, that a section of STIT tessellation with a lower-dimensional plane is again a STIT tessellation, stereological formulae are given in Mecke et al. (2009).
Now we will consider the second moment measure, the K-fianction and the pair-correlation fianction of a STIT tessellation.
If is any stationary segment process (as a point process on the parameter space introduced above) its intensity measure Axp can be factorized
Avp = Nil2 X p ,
(6)
where N\ is the mean number of segment centers per unit area and p is the joint distribution of direction and length of typical segment. Notice that the length distribution of the typical segment, is a marginal distribution of p. The directional distribution M is the length weighted directional distribution derived irom p.
If is additionally isotropic, then p(d(a,5)) = lda^(d5).
Theorem 1 Let ^ be a stationary STIT tessellation with intensity Z,/i o and Joint distribution po of length and direction of the typical I-segment. Further, let be a stationary Boolean segment process with intensity La^ and Joint distribution pxp of length and direction
of the typical segment. For the respective second order (2) (2) measures and we obtain
L^^qt =	and po = p^
(2) (2) Mo =MV-
Corollary: With Eq. 1 and Eq. 5 the K-flinction of a stationary and isotropic planar STIT tessellation O with 5 = tt/La as the mean length of the typical I-segment, see Nagel and Weiß (2006), is


A=1
k-k\
Using an exponential integral Ei(x) and the Euler-Mascheroni constant y we obtain
LAK^{r) =	I^r) -Ei(-fl^r)) + W.
Then the pair-correlation function is 1 dÄo(r) . 1
g{r) =
2nr dr
= 1
2Ly
l-e-l^^Ar] (7)
Proof of the Theorem: The proof will be based on Eq. 2 and we will show that the respective addends for O und are identical.
1st addend: For any ^ G ^^ with normal direction Yg {i.e., the angle between the normal and the positive x-axis) we consider the following four marked section point processes on ^and show that they are identically distributed. The marks will always be the section angles in [0, n) between g and the intersecting lines or segments respectively. Notice that the segment and the line processes have the same intensity La and directional distribution
i)	The Boolean model with La and intensity measure Axp as given in Eq. 6 generates on ^ a stationary Poisson point process with intensity Pl = La sin(a -	and with independent marks. The mark distribution H is given by
Pi//((0,j8]) =	|sm(a-r^)|^(da), (8)
for ß e (0,7i], see Stoyan et af (1995), p. 289.
ii)	Let r be a Poisson line process (as a point process on the parameter space introduced above) with the intensity measure A = LaIi x It can also be considered as a planar fibre process (according to Stoyan et af, 1995, p. 280), and thus its intersection with g generates a stationary Poisson point process with intensity Pl = LaJ\ sin(a — Yg) l^(da) and with independent marks. The mark distribution //coincides with that one in (i).
iii)	Let 11 be a Poisson point process on ^^ x (0,oo) (i.e., lines with birth times) with the intensity measure Ax£^ where A is the same as in (ii). For all a > 0 the line process 11^ = {A G ^ : {h, t) G n, t < a} has the length intensity aLA. Hence the intersection with g generates a stationary Poisson point process with intensity aPi = aLA f I sin(a — 7g)|^(da) and with independent marks, and the mark distribution is the same as in (i) and (ii).
iv)	Choose a compact convex polygonal window W with	Then for a STIT tessellation O the marked intersection point process within W, i.e., H^n^nO, can be generated by the algorithm that is described above, and using the measure A = LaI\ x Its distribution on H^n^ coincides for a = 1 with that one from (iii) (for a proof see below). Since the construction of STIT tessellations fialfills a consistency property (see Theorem 1 in Nagel and Weiß, 2005) we can conclude that the marked section point process on ^generated by the STIT O is identically distributed as that one from (iii).
For a stationary STIT tessellation O it is known from Nagel and Weiß (2003), Lemma 5, that for all ^ G ^^ the intersection process O n ^ is a stationary Poisson point process on g. Thus it is easy to see that all the (unmarked) intersection point processes on ^are stationary Poisson point processes. Therefore the focus of the proof is to show that the intersection angles are independent and identically distributed.
The identity of distributions of the marked section point processes in (i) and (ii) is obvious. The equivalence for (ii) and (iii) for a = 1 is straightforward since the distributions of the line processes T and 111 are identical. The identity of the distributions for (iii) with a = 1 and (iv) within any compact convex window W can be shown as follows. Let g e^^ he a line with
wng^0.
We consider (iv). First note, that if the fixed line g intersects an extant polygon Wj then the time until the chord fi ^ is intersected by a random line Yj is exponentially distributed with the parameter A([Wf n ). This follows from the calculation of the capacity fimctional for the constructed tessellation in Nagel and Weiß (2005).
Now, for a fixed time ? > 0 let	he all
those cells ofO(/', W) with 7; ^ = Wlng^ 0,... = n ^ / 0. Assume that the intervals Jf pi = l,...,m, are ordered from left to right (if g is vertical then bottom-up). Let aj be the section angle (with g) in the point that separates TJy and J, 1= \,...,m— 1. (We
can neglect the case that g goes through a vertex of a polygon.) Denote by {Xt)t>o the random process which at time t is in the state	• • • >J-
The state Xq at time t = 0 is {Wf^g) a.s. Further, {Xt)t>o is a Markov process due to the independence and the exponentially distributed lifetimes of the cells used in the algorithm. For a fixed bounded segment n^ the process {Xt)t>o is piecewise constant and it jumps into a new state when a new intersection point appears. For any fixed ? > 0 the distribution of the time until the next appearance of a section point is exponentially distributed with parameter A([I4^n g\) = E?LiA([)j,J) which is the distribution of the minimum of w independent exponentially distributed random variables with the respective parameters i= l,...,m. The increment or jump when the process changes its state is determined by a marked intersection point (x, a). The probability that x appears m IS	= UJ/\Wn
where | • | denotes the length of a segment. (We emphasize that this formula holds for arbitrary directional distributions.) Hence the distribution of x is a mixture of m uniform distributions on the intervals Jfj with the respective weights. Thus x is uniformly distributed on WCig. The product form of the measure A = La£i X ^ implies that the intersection angle a is independent of x, and the distribution of a is given by its distribution fianction given in Eq. 8.
Now we study (iii) and observe that the process n induces a process {Zt)t>o with states {Jt,h0ChJt,2,---,0Cn-hJt,n) where n-\	:
{h,s) G n,5 < t,gr\hr\W ^ and	are
the section angles (again left to right or bottom-up, respectively) and 7f,i,..., Jt,n are the ordered intervals in n ^ generated by the n — I sections with lines of n. The state Zq at time t = 0 is ( Wrig) a.s. Since n is a Poisson point process the process {Zt)t>o has the Markov property. For the bounded segment H^n^the process {Zt)t>o is piecewise constant and it jumps into a new state when a new intersection point appears. For any fixed ? > 0 the distribution of the time until the next appearance of a section point is exponentially distributed with parameter Due to the product form of the intensity measure A, the distribution of the appearing marked intersection point {x,a) given by the independence of x and a, is the uniform distribution of x on Wf^g and angle distribution as in Eq. 8.
Thus it is shown that for all ^ G ^ , all ? > 0 and all bounded intervals H^n^on ^the distributions of Xt and Zt and hence of the corresponding marked section point processes are identical. Together with the consistency result (Eq. 3) for the construction of STIT
tessellations in bounded windows the identity of the distributions in (iii) and (iv) is shown.
Summarizing, we can conclude that for all ^ G ^^ the marked section point processes on ^either induced by O or by respectively, are identically distributed.
2nd addend: The second summand of Eq. 2 is formed of L{(pf^g) ox L{•^ff^g), respectively. For O as well as for and for ^G ^^ we have either L{{■)f^g) = 0 or /,((•) n^) consists of a.s. exactly one I-segment of (p or of one segment of y/, respectively. Hence we can rewrite the sum X in Eq. 2 as the sum over
all I-segments s of \j/ and analogously for (p and its I-segments. With the intensity measure Avp of the stationarity of and the notation Lg = {t eR : y{t) e s} we obtain
I I [ [ f{y{ti),At2))dt,dt2R,{dxif) I. f f f{Ati),Ät2))dhdt2P^{dxif)
= 111 fim,Ät2))dhdt2A^id5).
L, L,
Since L/i^o = La^ and po = p^ v^Q have A'l o = Ni^^ and hence Ao = Axp and thus the equality of the second summand in Eq. 2 for a STIT tessellation O and a Boolean segment process respectively.
Interpretation of the result: Theorem 1 shows, that the second moment measure of the considered segment processes does not indicate whether the segments are arranged "completely randomly", i.e., independently, or in a "highly dependent" way such that they yield a tessellation. In particular, the second moment measure does not depend on the type of crossings or nodes that are generated by the segments.
Nevertheless, it is an open problem whether the second moment measure is applicable to discriminate between different tessellations (e.g., STIT, Voronoi, Poisson line tessellations).
ESTIMATION OF THE COVARIANCE FUNCTION
In the previous sections we studied the second (reduced) moment measure of STIT tessellations with the help of section points and angles on fixed lines. While this was a usefial tool for a proof of our result it is not the method of choice in image analysis. In particular, the measurement of section angles in
digitized images is vague and problematic. The aim of estimation the covariance by image analysis methods is to compare the theoretical covariance function of a STIT tessellation with estimates from images of a fibre system in order to get a (necessary) criterion for the hypothesis that the observed fibre system forms a STIT tessellation.
In the present section we describe a feasible estimator for the covariance fianction of the edge system of tessellations. This estimator will be based on the Fourier transform of the smoothed edge system inside an observation window W. Here, a smoothed edge means a Sanction after a convolution of the edge with a kernel Sanction K, i.e., a smoothing of the contrast in an image. In the following an estimator for digitized images is suggested. This approach holds for a general fibre process and it is not specific for STIT tessellations.
We commence with some notation and a review of some facts for general stationary fibre processes (the edge system of a tessellation can be considered as a particular fibre process). A more detailed presentation is given in Ohser and Schladitz (2009), Section 6.4.
For a stationary and isotropic planar fibre process with pair-correlation Sanction ^and length intensity La the covariance fianction cowp is given by
covvi/(x) = La

This covariance is the density of the covariance measure of the length measure that is induced by We remark that for a stationary and isotropic planar STIT with ^in Eq. 7, the covariance fianction cowp is not integrable, ^^ covvj/(x)dx »o as r^ »o.
The convolution /i * /2 of two integrable fianctions /i,/2 : M^ M is defined by [/1 * fz^x) = /r2 /1 (7) f2{x—y)dy,xe M^. Similar to the convolution of fianctions we define the convolution ß * f: M.^ ^ RU {—of a measure ß on R^ with a fianction / that is nonnegative and measurable w. r. t. ß.
[ß*f]{x)= f{x-y)ß{dy),
xG:
As above, we identify the edge system of a tessellation with a random (length) measure Let (C : R^ R be a nonnegative function with compact support and satisfying /^2 K{x)dx= 1. We define the random function /j/ = * (c — La, associated with a random fibre system *F. As mentioned above, *F * (c can be interpreted as a smoothing (of the contrasts in an image) of the edge system.
The function *F * (C is not necessarily integrable but almost surely integrable on every compact subset of
R2. Finally, from E[*F * k]{x) < 00 for all x G R^, it follows that the random function fy is almost surely locally integrable and the reduced covariance measure of % has a density cov /, the covariance function of fy. It follows that
COV/ = (ic* K*) * COVvp
with K* = - K.
Now let Whe a compact window with nonempty interior. We are smoothing the fibre process *F and observe it in W. This yields the random function
= {[^'*k]{x)-La)Iw{x), xer\
which is almost surely integrable and thus its Fourier transform exists.
By f=ß^fwe denote the 2-dimensional Fourier transform of an integrable function /: R^ C,
The corresponding cotransform may be denoted by 1
[J^f]{x) =
2n


Now we can formulate a Wiener-Khintchine type theorem.
Theorem 2 Let be a stationary random Ghre process with a locally ßnite ßrst moment measure and an existing covariance function covvp. Let (C : R^ R be a bounded nonnegative function of compact support and satisfying /^2 fc(x)dx = 1. Furthermore, let W be a compact window of nonempty interior For the windowed random function /j/it follows that for
2;rE|	= ^{cw {{K* /C*) *covvp))(^), (9)
where cw is the window function cw =
K proof is given in Unverzagt (2005), Section 2.2.
From Eq. 9 it immediately follows that a smoothed version {k* K*) * covxp of the density covxp can be estimated via frequency space. Assume that W is compact and the origin belongs to its interior, then
(k* K*) * covxp] (x) = IKE
cw{x)
(10)
for _all x in the interior of W. This means, 271,:^(x)/cw{x) is an unbiased estimator for the expression on the left-hand side.
Let (fCe)e>o be a family of bounded nonnegative kernel functions (Cg : R^ i-^ R with /^2 Ke{x)dx= 1 and K'g(x) = 0 for ||x|| > e. Then, if the density covvp is continuous in x, it follows that
lim((K-e
ejO
* (Cg) * COVvp ) = COVvp ,
pointwise.
Lemma: For a stationary planar fibre process a compact window W with nonempty interior and a family of bounded nonnegative kernel fianctions (Cg : R with /]g,2 Ke{x)dx= 1 and Ke{x) = 0 for HxH >
cw{x)
with fy^ w,s (^) = (* K-fi] (^) - La) 1 w(x), a- G R^, is an asymptotically unbiased estimator for the covariance fianction covvp of as e ^ 0.
AN ESTIMATOR FOR DIGITIZED IMAGES
In practical applications, realizations of fibre processes are usually observed on lattices and we apply methods of image analysis in order to estimate the pair-correlation fianction (or the covariance fianction covxp, respectively) of Let L^ = aZ^ be a square lattice with lattice spacing a > 0 and the unit cell C= [0,a]^, where [0,a] is the closed segment between 0 and a. A digitization (or sample in the statistical sense) ^ii of on the lattice L^ may be defined as the set of the lattice points x with
This is only a simplified model for the much more complicated sampling of real fibre processes on a lattice. Furthermore we remark that, in general, the process can not be reconstructed from ^fi. In the following we assume that only the sampling *I'l2 is known but not itself Our aim is now to estimate the pair-correlation fianction of the unknown process fi-om the data *I'l2 observed in a window W.
We refer to *I'l2 as the set of the foreground pixels and the complementary set = L^ \ ^12 is called the background. Following the approach in Ohser et al. (2009) we consider local pixel configurations (^0, • • •, <^15 defined as subsets of the set of the vertices of the unit cell, C Cn L^. Pictograms of these pixel configurations are shown in Table 1, where the fiall discs mark the foreground pixels and the circles mark the background. In our setting, the indexing of the pixel configurations is chosen such that the
complementary pixel configurations are = (Cn
The total length will be estimated by summing up the contributions from the local pixel configurations in the digitized image. The local contribution to the estimation of the length of of such pixel configurations is given by weights wi. An appropriate choice of these weights is a non-trivial problem. Here, pragmatically, we use the weights which we calculated in Ohser et al. (2009) for another mode of digitization of sets. (In this article the boundary length estimation for a random set on R^ with positive area fraction is considered, where the intersection of the random set with Jj is suggested as an appropriate digitization mode and a corresponding discretization of a Crofton integral formula is applied.)
Table 1. The local pixel conßgurations ofh^, the corresponding pictograms and the length weights w^.
I
1
2
3
4
5
6 7
Wl
i
O
0.335 190 0.335 190 0.474030 0.335 190 0.474 030 0.392 699 0.335 190
9 10 11 12
13
14
15
Wl
0.335 190 0.392 699 0.474030 0.335 190 0.474 030 0.335 190 0.335 190 O
Let now be a rectangular window with edges parallel to the co-ordinate axes. By W = 14^0 Č we denote the reduced window, where 0 is the Minkowski subtraction and C is the unit cell reflected at the origin, C = —C. We assume that the window is much larger than the unit cell such that #(L2 n W^') > 0. Then the length density La of can be estimated using the length weights wi given in Table 1. Let h = {hi) be the vector of the number of pixel configurations in ^L^nH/',
xel7r\W'
= 0,..., 15. Then La can be estimated using
hw
La =
a#{h^nW') '
with the vector w= (wg) and hw is the scalar product of the vectors h and w. If *F is stationary and isotropic, then the estimator La is asymptotically unbiased for La as a I 0 (multigrid convergent, see Ohser and Schladitz, 2009, Section 5.2).
gir)
4-
2-
0.5
r/lA
Fig. 2. The pair-correlation function g{r) of the isotropic STIT (thick line) compared with the pair-correlation function of the stationary ILPL.
0.2-
^gir)
0.1-
^Oo
Oo
Oo
°°°°°°ooooo
ooo
0-1-1-1-1-1-1-1-1-1-r-
0	0.5	,,
rjlA
Fig. 3. Differences Ag= g— g of estimates g of pair-correlation functions and the true pair-correlation function g of the isotropic STIT; full discs: g from realizations of a STIT; circles: g from realizations of an IPLP. The solid curve shows the difference of the true pair-correlation functions for the ILPL and the STIT, respectively.
In order to estimate the covariance function covxp by applying a discrete Fourier transform, a discretization of the estimator on the right-hand side of Eq. 10 is needed. The length weights wi given in Table 1 are used in order to obtain an appropriate representation of the function f^^w- Consider the function 4,2 : L^ i-^ R mapping each lattice point to the local contribution awi for the (random) length measure of We determine the index I of the local pixel configuration {^'j^i-x) n (CnL^) of at the lattice point Xand assign an appropriate length weight awi to
42« =
a j^wtli^td -c -x)-la,
i=0
for X G L2 n W. The set {{x, 42(x)) : xeh^nW'} forms a (random) grey-value image with real-valued pixels. This can be understood as the result of a nonlinear filtering (in contrast to the linear filter in the continuous model), applied to Finally, the function covvp can be estimated from this image via the inverse space {i.e., the space of the Fourier transform) using formula (10) where the continuous Fourier transform is replaced with a discrete Fourier transform. The choice of the lattice spacing a depends on the length density La (small a for large La).
The use of the method described above is demonstrated in the following example. First, we observe in Fig. 2 that the difference between the pair-correlation functions of isotropic STITs and stationary and isotropic Poisson line processes (IPLP) is very small. Thus, due to estimation errors, it seems to be difficult to see these differences also in estimates of the pair-correlation functions. However, for large samples (i.e., large windows or averaging over many realizations) and small lattice spacings, it is possible to discriminate between both processes. Estimates g of pair-correlation functions and the pair-correlation functions g of the isotropic STIT (full discs) and an ILPL (circles), respectively, are compared in Fig. 3. The data for g were computed for realizations with La = 32, observed in a square window of edge length 1 and on a square lattice of spacing a = 1/256. The estimates were averaged over 16 realizations. In order to display the very small estimation errors. Fig. 3 shows the differences ls.g = g— g oi estimates g of pair-correlation functions and the true pair-correlation function g of the isotropic STIT. We remark that the length density La (used for the scaling of the x-axis in the diagram) was estimated from the same data. Fig. 3 shows that estimates of the pair-correlation function can be used in order to discriminate between different kinds of tessellations even if the differences between the pair-correlation functions are small (as for STIT and ILPL). Clearly, an ILPL does not have T-nodes which is a more significant criterion in order to discriminate from STITs. However, ILPLs often serve as a benchmark for comparison studies, and their pair-correlation functions are known explicitly.
Remarks: For m pixels of an image of the covariance cowp can be computed by the use of the Fast Fourier Transform (FFT) with a complexity in iff{m\o%m), because the FFT has a complexity of (^(mlogm)_and the window function can be computed via cw =	■
The function ^2 differs from />j/ * (C, i. e. we can not find a kemel function K such that 42 (^) = [ ^ * k:]{x), x eh^. As a consequence we can not show
multigrid convergence for the estimator of the pair-correlation function based on the procedure peresented in the previous section. Nevertheless, the examples in Fig. 3 show that the pair-correlation functions for a STIT tessellation and a Poisson line process can be estimated with high accuracy.
Unfortunately, the assumption of periodicity in the discrete Fourier transform causes an overlapping effect (edge effect), see Koch et al. (2003). This effect can be eliminated by expanding the function f^^w to the window 214^: f^^wi^ =	H x G W and
h'ßwix) = 0 if xG 2W\ W, i.e., the original image is padded with zeros. This increases the number of pixels to 4w. Still the complexity belongs to «^(mlogm) which is a considerable gain compared to the usual estimation of the covariance function by a convolution which is of complexity ffin?).
Furthermore, notice that the Bartlett spectrum of does not exist. The Bartlett spectrum of is the quantity in the inverse space {i.e., the space of the Fourier transforms) associated with cowp, also known as the scattering intensity of see Ohser and Schladitz (2009, Section 6.3).
REFERENCES
Baddeley AJ, Silverman BW (1984). A cautionaiy example on the use of second-order methods for analyzing point patterns. Biometrics 40:1089-93.
Koch K, Ohser J, Schladitz K (2003). Spectral theoiy for random closed sets and estimating the covariance via frequency space. Adv Appl Prob 35:603-13. Mecke J (2009). Joint distribution of direction and length of the typical I-segment in a homogeneous random planar tessellation stable under iteration. J Contemp Math Anal 44:45-53.
Mecke J, Nagel W, Weiß V (2007). Length distributions of edges in planar stationary and isotropic stit tessellations. J Contemp Math Anal 42:28-43. Mecke J, Nagel W, Weiß V (2008a). A global construction of homogeneous random planar tessellations that are stable under iteration. Stochastics Int J Prob 80:51-67.
Mecke J, Nagel W, Weiß V (2008b). The iteration of random tessellations and a construction of a homogeneous process of cell divisions. Adv Appl Prob SGSA 40:4959.
Mecke J, Nagel W, Weiß V (2009). Homogeneous STIT tessellations - a stereological aspect. In: Capasso V, Aletti G, Micheletti A, eds., Stereology and Image Analysis. ECS 10 The 10th European Congress of ISS, MIRIAM. ECS 10, Milano, Bologna: Societa Editrice Esculapio Bologna.
Mecke J, Nagel W, Weiß V (2010). Some distributions for I-segments of planar random homogeneous STIT tessellations. Math Nachr (in press).
Nagel W, Weiß V (2003). Limits of sequences of stationaiy planar tessellations. Adv Appl Prob SGSA 35:123-38.
Nagel W, Weiß V (2005). Crack tessellations: characterization of stationaiy random tessellations stable with respect to iteration. Adv Appl Prob SGSA 37:859-83.
Nagel W, Weiß V (2006). STIT tessellations in tiie plane. Rend Circ Math Palermo II Suppl. 77:441-58.
Nagel W, Weiß V (2008). Mean values for homogeneous stit tessellations in 3D. Image Anal Stereol 27:29-37.
Ohser J, Schladitz K (2009). 3D images of materials structures - Processing and analysis. Weinheim, Berlin: Wiley VCH.
Ohser J, Schladitz K, Nagel W (2009). Miles formulae for Boolean models observed on lattices. Image Anal Stereol 28:77-92.
Schwandtke A (1988). Second-order quantities for stationaiy weighted fibre processes. Math Nachr 139:321-34.
Stoyan D, Kendall WS, Mecke J (1995). Stochastic and integral geometiy. Chichester: J. Wiley & Sons, 2nd ed.
Unverzagt K (2005). Spectral analysis of random closed sets: the surface measure associated with a random closed set. Master's thesis, Universität Siegen.
Weiß V, Nagel W (1994). Second-order stereology for planar fibre processes. Adv Appl Probab 26:906-18.