SECOND MOMENT MEASURE AND K-FUNCTION FOR PLANAR STIT TESSELLATIONS

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 logm).


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, i.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 paircorrelation 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 2 the Euclidean plane and R 2 the σ -algebra of Borel sets in R 2 .By ℓ 2 , ℓ 1 , ℓ + we denote the Lebesgue measures on R 2 , R and [0, ∞), respectively.Further, we use the notations ψ for a planar fibre system and Ψ for a random fibre process in the sense of Stoyan et al. (1995), P Ψ denotes the distribution of Ψ.A fibre process Ψ is said to be macroscopically homogeneousi.e., spatially stationary -if P Ψ is invariant under all translations of R 2 .It is isotropic if its distribution is invariant under rotations.The total fibre length of the fibre system ψ in B ∈ R 2 is written as ψ(B).For a stationary fibre process the mean total fibre length per unit area (length intensity) is L A .
The second moment measure µ (2) If Ψ is stationary, the reduced second moment measure K Ψ can be defined by µ (2) where 1 B denotes the indicator function of B.
In the case of a stationary and isotropic fibre process it is sufficient to consider the reduced second moment function K Ψ (called K-function) given by where B(o, r) is the circle with radius r centered in o.
The product L A K Ψ (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 secondorder process if L A K Ψ (B) < ∞ for all bounded B ∈ R 2 .In the following we only consider second-order processes.And furthermore, we assume that also the intersection point processes Ψ ∩ g on all lines g are of second order.
If the K-function is differentiable, the paircorrelation function g exists with We consider two particular cases of fibre systems.
a) Let Ψ be a stationary and isotropic Poisson line process with length intensity L A .Then the Kfunction is given by (cf.(Stoyan et al., 1995)) and the pair-correlation function is b) For line segments we choose the parametrization (x, α, s) where x ∈ R 2 is the centre, α ∈ [0, π) the normal direction (i.e., the angle between the normal of the segment and the positive x-axis) and s ∈ (0, ∞) the length.Thus a segment process Ψ can be considered as a point process on the space R 2 × [0, π) × (0, ∞).Notice that the normal direction will be used in the parametrization of both the segments and the lines.Now let Ψ be a stationary and isotropic Boolean segment process with intensity N 1 of the centre point (germ) process and with length distribution L of the typical segment.Then the length intensity is L A = N 1 s where s denotes the mean length of the typical segment.For the K-function we obtain (cf.Stoyan et al., 1995) The formulae for L A K Ψ (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 P Ψ , 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 G the set of all lines in R 2 and by y−z the distance between two points y and z in R 2 .Further, let dg be the element of the motion invariant measure on the set of all lines with 1{g∩B(o, 1) = / 0}dg = 2π.
For a fibre system ψ (here considered as a subset of R 2 ) and g ∈ G we consider ψ ∩ g which consists of isolated intersection points and of linear segments of positive length of ψ on g: where S(ψ ∩g) denotes the set of all intersection points and L(ψ ∩ g) the union of all straight line segments on g.It is clear, that only for countably many g we have L(ψ ∩ g) = / 0.
An intersection point y ∈ S(ψ ∩ g) is marked by the fibre tangent angle w(y, g).This is the angle between the tangent of ψ 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 ∈ G is parametrized as g = {y(t) : t ∈ R} such that y(t 1 ) − y(t 2 ) = |t 1 −t 2 | for all t 1 , t 2 ∈ R.This yields a parametrization of the union of segments L(ψ ∩ g) as L g = {t ∈ R : y(t) ∈ L(ψ ∩ g)}.
Then the following stereological formula for the second moment measure holds.Let Ψ be a stationary planar fibre process.Then for all measurable nonnegative functions f on R 2 × R 2 we have where u = w(y, g), 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.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 globali.e., in the whole plane R 2 -construction was given in Mecke et al. (2008a).This construction can be understood as a process of sequential cell division at random times.
Let Λ = L A ℓ 1 × R be a measure on R × [0, π) and R is a probability measure on [0, π), the directional distribution.In order to generate a tessellation it is assumed that R is not concentrated on a single direction.A point (p, α) ∈ R × [0, π) represents a line g(p, α) ∈ 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 α is the angle between its normal and the positive x-axis.For a set At time t = 0 the construction starts with a compact convex polygonal window W ⊂ R 2 (e.g., a rectangle).After a random life-time t 1 that is exponentially distributed with parameter Λ( ).Thus two new polygons are born, and the birth time t 1 is attributed to them.Then, sequentially, all the extant polygons W ′ 1 , . . .,W ′ k , say, with the respective birth times t ′ 1 , . . .,t ′ k are divided in the following way.The life-time of W ′ i is a random variable that is exponentially distributed with parameter Λ([W ′ i ]), i.e., the parameter depends on the size of W ′ i .In the particular case that R is the uniform distribution on [0, π) it is proportional to the perimeter of W ′ i .At the end of its life-time the polygon is divided by a chord that comes from a random line that has the distribution Λ( Thus W ′ i dies and two new polygons are born.The state at a fixed time a > 0 is a tessellation in W and it is denoted by Φ(a,W ).
A crucial feature of the construction 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 construction: Let (τ j , γ j ), j = 1, 2, . .., be a sequence of independent and identically distributed (i.i.d.) random variables, where also τ j and γ j are independent.The τ j are exponentially distributed positive numbers with parameter Λ([W ]).The γ j are random lines with the distribution Λ( ).Any of these random variables can be used at most once during the construction.The division of an extant polygon W ′ i ⊂ W is performed according to the rejection rule: For a pair (τ j , γ j ) throw -after the time τ j elapsed - the line γ j onto the window W .If it hits W ′ i then use the generated chord to divide this polygon.If γ j ∩ W ′ i = / 0 then reject γ j and make a new trial with another pair (τ j ′ , γ j ′ ).Also in the case of rejection τ j has to be added to the lifetime of W ′ i .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 construction 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 construction 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 Φ of the whole R 2 such that the symbol D = stands for the identity of distributions of random variables.
The distribution of the tessellation Φ does not depend on W , and one can show that this formula holds for all compact and convex 2-dimensional sets W ⊂ R 2 .Moreover, the existence of a stationary tessellation is sufficient that the intensities and the K-function are well defined.
Generalizations to higher dimensions are given in Nagel and Weiß (2005); Mecke et al. (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 2 .
The operation of iteration (also referred to as nesting) for tessellations is defined as follows.Denote by Y = {Y 1 ,Y 2 , . ..} a sequence of independent and identically distributed (i.i.d.) stationary random tessellations.Further assume that Y 0 is a stationary random tessellation which is independent of Y .It is useful to consider the set C(Y 0 ) = {p 1 , p 2 , . ..} of the cells (which are convex polytopes) of Y 0 .The iteration of the tessellation Y 0 and the sequence Y is defined as This definition means that a cell p k of the so called 'frame' tessellation Y 0 is -independently of all other cells -subdivided by the cells p ki , i = 1, 2, . . . of the tessellation Y k which intersect the interior of p k .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 Y 0 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 Y 0 .This is now defined more precisely.For a real number r > 0 the tessellation rY is generated by transforming all points (x, y) ∈ Y into (rx, ry).Accordingly, rY means that this transformation is applied to all tessellations of the sequence Y .Let Y 0 be a stationary random tessellation and Y 1 , Y 2 , . . .a sequence of sequences of tessellations such that all the occurring tessellations (including Y 0 ) are i.i.d.Then the sequence I 2 (Y 0 ), I 3 (Y 0 ), . . . of rescaled iterations is defined in Nagel and Weiß (2003;2005) as Here m is the rescaling factor which is chosen such that the results of iteration do not degenerate for m → ∞.
Definition: A stationary random tessellation Y is said to be stable with respect to iteration (STIT) if 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) Stability with respect to iteration: The stationary tessellation Φ given in (3) is stable with respect to iteration.

Characteristic entities:
The segments (edges of cells) in the tessellation Φ form a stationary segment process with L A (which appeared in the definition of Λ) as the mean total edge length per unit area and the directional distribution R, i.e., the distribution of the normal direction in a randomly chosen point on the segments.By L A and R the distribution of stationary STIT tessellations is already uniquely determined.
Poisson typical cell: Now we consider the interior of the typical cell of Φ.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 Φ with directional distribution R and intensity L A has the same distribution as the interior of the typical cell of a stationary Poisson line tessellation with the same R and L A , see Nagel and Weiß (2003).
Length distribution of the typical I-segment: For a stationary and isotropic STIT tessellation Φ with intensity L A the density of the length of the typical Isegment is for x > 0 see Mecke et al. (2007).
This length distribution is a mixture of exponential distributions 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-function and the pair-correlation function of a STIT tessellation.
If Ψ is any stationary segment process (as a point process on the parameter space introduced above) its intensity measure Λ Ψ can be factorized where N 1 is the mean number of segment centers per unit area and ρ is the joint distribution of direction and length of the typical segment.Notice that L , the length distribution of the typical segment, is a marginal distribution of ρ.The directional distribution R is the length weighted directional distribution derived from ρ.
Theorem 1 Let Φ be a stationary STIT tessellation with intensity L A,Φ and joint distribution ρ Φ of length and direction of the typical I-segment.Further, let Ψ be a stationary Boolean segment process with intensity L A,Ψ and joint distribution ρ Ψ of length and direction of the typical segment.For the respective second order measures µ (2) Φ and µ (2) Ψ we obtain Ψ .
Corollary: With Eq. 1 and Eq. 5 the K-function of a stationary and isotropic planar STIT tessellation Φ with s = π/L A as the mean length of the typical Isegment, see Nagel and Weiß (2006), is Using an exponential integral Ei(x) and the Euler-Mascheroni constant γ we obtain Then the pair-correlation function is Proof of the Theorem: The proof will be based on Eq. 2 and we will show that the respective addends for Φ und Ψ are identical.
1st addend: For any g ∈ G with normal direction γ g (i.e., the angle between the normal and the positive x-axis) we consider the following four marked section point processes on g and show that they are identically distributed.The marks will always be the section angles in [0, π) between g and the intersecting lines or segments respectively.Notice that the segment and the line processes have the same intensity L A and directional distribution R.
i) The Boolean model Ψ with L A and intensity measure Λ Ψ as given in Eq. 6 generates on g a stationary Poisson point process with intensity P L = L A | sin(α − γ g )|R(dα) and with independent marks.The mark distribution H is given by for β ∈ (0, π], see Stoyan et al. (1995), p. 289.
ii) Let Γ be a Poisson line process (as a point process on the parameter space introduced above) with the intensity measure Λ = L A ℓ 1 × R. It can also be considered as a planar fibre process (according to Stoyan et al., 1995, p. 280), and thus its intersection with g generates a stationary Poisson point process with intensity P L = L A | sin(α − γ g )|R(dα) and with independent marks.The mark distribution H coincides with that one in (i).
iii) Let Π be a Poisson point process on G × (0, ∞) (i.e., lines with birth times) with the intensity measure Λ × ℓ + where Λ is the same as in (ii).For all a > 0 the line process Π a = {h ∈ G : (h,t) ∈ Π, t < a} has the length intensity aL A .Hence the intersection with g generates a stationary Poisson point process with intensity aP L = aL A | sin(α − γ g )|R(dα) 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 W ∩ g = / 0. Then for a STIT tessellation Φ the marked intersection point process within W , i.e., W ∩ g ∩ Φ, can be generated by the algorithm that is described above, and using the measure Λ = L A ℓ 1 × R. Its distribution on W ∩ g coincides for a = 1 with that one from (iii) (for a proof see below).Since the construction of STIT tessellations fulfills a consistency property (see Theorem 1 in Nagel and Weiß, 2005) we can conclude that the marked section point process on g generated by the STIT Φ is identically distributed as that one from (iii).
For a stationary STIT tessellation Φ it is known from Nagel and Weiß (2003), Lemma 5, that for all g ∈ G the intersection process Φ ∩ g is a stationary Poisson point process on g.Thus it is easy to see that all the (unmarked) intersection point processes on g 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 Γ and Π 1 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 ∈ G be a line with W ∩ g = / 0.
We consider (iv).First note, that if the fixed line g intersects an extant polygon W ′ i then the time until the chord W ′ i ∩ g is intersected by a random line γ j is exponentially distributed with the parameter Λ([W ′ i ∩ g]).This follows from the calculation of the capacity functional for the constructed tessellation in Nagel and Weiß (2005).Now, for a fixed time t > 0 let W ′ 1 , . . .,W ′ m be all those cells of Φ(t,W ) with . Assume that the intervals J ′ t,i , i = 1, . . ., m, are ordered from left to right (if g is vertical then bottom-up).Let α ′ i be the section angle (with g) in the point that separates J ′ t,i and J ′ t,i+1 , i = 1, . . ., m − 1. (We can neglect the case that g goes through a vertex of a polygon.)Denote by (X t ) t≥0 the random process which at time t is in the state ).The state X 0 at time t = 0 is (W ∩ g) a.s.Further, (X t ) t≥0 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 W ∩ g the process (X t ) t≥0 is piecewise constant and it jumps into a new state when a new intersection point appears.For any fixed t > 0 the distribution of the time until the next appearance of a section point is exponentially distributed with parameter Λ( ) which is the distribution of the minimum of m independent exponentially distributed random variables with the respective parameters Λ([J ′ t,i ]), i = 1, . . ., m.The increment or jump when the process changes its state is determined by a marked intersection point (x, α).The probability that x appears in (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 J ′ t,i with the respective weights.Thus x is uniformly distributed on W ∩ g.The product form of the measure Λ = L A ℓ 1 × R implies that the intersection angle α is independent of x, and the distribution of α is given by its distribution function given in Eq. 8. Now we study (iii) and observe that the process Π induces a process (Z t ) t≥0 with states 0} and α 1 , . . .α n are the section angles (again left to right or bottom-up, respectively) and J t,1 , . . ., J t,n are the ordered intervals in W ∩ g generated by the n − 1 sections with lines of Π.The state Z 0 at time t = 0 is (W ∩ g) a.s.Since Π is a Poisson point process the process (Z t ) t≥0 has the Markov property.For the bounded segment W ∩ g the process (Z t ) t≥0 is piecewise constant and it jumps into a new state when a new intersection point appears.For any fixed t > 0 the distribution of the time until the next appearance of a section point is exponentially distributed with parameter Λ([W ∩ g]).
Due to the product form of the intensity measure Λ, the distribution of the appearing marked intersection point (x, α) given by the independence of x and α, is the uniform distribution of x on W ∩ g and angle distribution as in Eq. 8.
Thus it is shown that for all g ∈ G , all t > 0 and all bounded intervals W ∩ g on g the distributions of X t and Z t 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 ∈ G the marked section point processes on g either induced by Φ or by Ψ respectively, are identically distributed.
2nd addend: The second summand of Eq. 2 is formed of L(ϕ ∩ g) or L(ψ ∩ g), respectively.For Φ as well as for Ψ and for g ∈ G we have either L((•) ∩ g) = / 0 or L((•) ∩ g) consists of a.s.exactly one I-segment of ϕ or of one segment of ψ, respectively.Hence we can rewrite the sum ∑ g : L(ψ∩g) = / 0 in Eq. 2 as the sum over all I-segments s of ψ and analogously for ϕ and its I-segments.With the intensity measure Λ Ψ of Ψ, the stationarity of Ψ and the notation L s = {t ∈ R : y(t) ∈ s} we obtain Since L A,Φ = L A,Ψ and ρ Φ = ρ Ψ we have N 1,Φ = N 1,Ψ and hence Λ Φ = Λ Ψ and thus the equality of the second summand in Eq. 2 for a STIT tessellation Φ 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 useful 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 function 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 function after a convolution of the edge with a kernel function κ, 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 function g and length intensity L A the covariance function cov Ψ is given by cov 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 g in Eq. 7, the covariance function cov Similar to the convolution of functions we define the convolution µ * f : R 2 → R ∪ {−∞, ∞} of a measure µ on R 2 with a function f that is nonnegative and measurable w. r. t. µ, As above, we identify the edge system of a tessellation with a random (length) measure Ψ.Let κ : R 2 → R be a nonnegative function with compact support and satisfying R 2 κ(x)dx = 1.We define the random function f Ψ = Ψ * κ − L A , associated with a random fibre system Ψ.As mentioned above, Ψ * κ can be interpreted as a smoothing (of the contrasts in an image) of the edge system.
The function Ψ * κ is not necessarily integrable but almost surely integrable on every compact subset of R 2 .Finally, from E[Ψ * κ](x) < ∞ for all x ∈ R 2 , it follows that the random function f Ψ is almost surely locally integrable and the reduced covariance measure of f Ψ has a density cov f , the covariance function of f Ψ .It follows that Now let W be a compact window with nonempty interior.We are smoothing the fibre process Ψ and observe it in W .This yields the random function which is almost surely integrable and thus its Fourier transform exists.
By f = F f we denote the 2-dimensional Fourier transform of an integrable function f : R 2 → C, The corresponding cotransform may be denoted by F , Now we can formulate a Wiener-Khintchine type theorem.
Theorem 2 Let Ψ be a stationary random fibre process with a locally finite first moment measure and an existing covariance function cov Ψ .Let κ : R 2 → R be a bounded nonnegative function of compact support and satisfying R 2 κ(x)dx = 1.Furthermore, let W be a compact window of nonempty interior.For the windowed random function f Ψ,W it follows that for ξ where c W is the window function A proof is given in Unverzagt (2005), Section 2.2.
From Eq. 9 it immediately follows that a smoothed version (κ * κ * ) * cov Ψ of the density cov Ψ can be estimated via frequency space.Assume that W is compact and the origin belongs to its interior, then for all x in the interior of W .This means, 2π F | fΨ,W | 2 (x)/c W (x) is an unbiased estimator for the expression on the left-hand side.
Lemma: For a stationary planar fibre process Ψ, a compact window W with nonempty interior and a family of bounded nonnegative kernel functions κ ε : , is an asymptotically unbiased estimator for the covariance function cov Ψ of Ψ as ε → 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 function (or the covariance function cov Ψ , respectively) of Ψ.Let L 2 = aZ 2 be a square lattice with lattice spacing a > 0 and the unit cell C = [0, a] 2 , where [0, a] is the closed segment between 0 and a.A digitization (or sample in the statistical sense) Ψ L 2 of Ψ on the lattice L 2 may be defined as the set of the lattice points x with This Ψ L 2 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 Ψ L 2 .In the following we assume that only the sampling Ψ L 2 is known but not Ψ itself.Our aim is now to estimate the pair-correlation function of the unknown process Ψ from the data Ψ L 2 observed in a window W .
We refer to Ψ L 2 as the set of the foreground pixels and the complementary set Ψ c L 2 = L 2 \ Ψ L 2 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 ∩ L 2 .Pictograms of these pixel configurations ξ ℓ are shown in Table 1, where the full 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 ξ 15−ℓ = (C ∩ L 2 ) \ ξ ℓ , ℓ = 0, . . ., 15.
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 w ℓ .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 2 with positive area fraction is considered, where the intersection of the random set with L 2 is suggested as an appropriate digitization mode and a corresponding discretization of a Crofton integral formula is applied.)Let now W be a rectangular window with edges parallel to the co-ordinate axes.By W ′ = W ⊖ Č we denote the reduced window, where ⊖ is the Minkowski subtraction and Č is the unit cell reflected at the origin, Č = −C.We assume that the window is much larger than the unit cell such that #(L 2 ∩ W ′ ) ≫ 0. Then the length density L A of Ψ can be estimated using the length weights w ℓ given in Table 1.Let h = (h ℓ ) be the vector of the number of pixel configurations in with the vector w = (w ℓ ) and hw is the scalar product of the vectors h and w.If Ψ is stationary and isotropic, then the estimator LA is asymptotically unbiased for L A as a ↓ 0 (multigrid convergent, see (Ohser and Schladitz, 2009, Section 5.2)).In order to estimate the covariance function cov Ψ by applying a discrete Fourier transform, a discretization of the estimator on the right-hand side of Eq. 10 is needed.The length weights w ℓ given in Table 1 are used in order to obtain an appropriate representation of the function f Ψ,W .Consider the function f L 2 : L 2 → R mapping each lattice point to the local contribution aw ℓ for the (random) length measure of Ψ.We determine the index ℓ of the local pixel configuration (Ψ L 2 − x) ∩ (C ∩ L 2 ) of Ψ L 2 at the lattice point x and assign an appropriate length weight aw ℓ to f 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 cov Ψ 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 L A (small a for large L A ).
The use of the method described above is demonstrated in the following example.First, we observe in Fig. 2 that the difference between the paircorrelation 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 L A = 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 ∆g = g − g of estimates g of pair-correlation functions and the true pair-correlation function g of the isotropic STIT.We remark that the length density L A (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 Tnodes which is a more significant criterion in order to discriminate from STITs.However, ILPLs often serve as a benchmark for comparison studies, and their paircorrelation functions are known explicitly.The function f L 2 differs from f Ψ * κ, i. e. we can not find a kernel function κ such that f L 2 (x) = [ f Ψ * κ](x), x ∈ L 2 .As a consequence we can not show multigrid convergence for the estimator of the paircorrelation 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.

Remarks
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 2W : f Ψ,2W (x) = f Ψ,W (x) if x ∈ W and f Ψ,2W (x) = 0 if x ∈ 2W \W , i.e., the original image is padded with zeros.This increases the number of pixels to 4m.Still the complexity belongs to O(m log m) which is a considerable gain compared to the usual estimation of the covariance function by a convolution which is of complexity O(m 2 ).
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 cov Ψ , also known as the scattering intensity of Ψ, see Ohser and Schladitz (2009, Section 6.3).

:
For m pixels of an image of Ψ the covariance cov Ψ can be computed by the use of the Fast Fourier Transform (FFT) with a complexity in O(m log m), because the FFT has a complexity of O(m log m) and the window function can be computed via c W = F |F 1 W | 2 .

Table 1 .
The local pixel configurations of L 2 , the corresponding pictograms and the length weights w ℓ .