Computational Statistics & Data Analysis 51 (2006) 859 – 871

www.elsevier.com/locate/csda

Statistical reconstruction of random point patterns

A. Tscheschela,∗ , D. Stoyanb
a Department of Mathematics, Hong Kong Baptist University, Kowloon Tong, Hong Kong
b Institute for Stochastics, TU Bergakademie Freiberg, 09599 Freiberg, Germany

Received 13 April 2005; received in revised form 12 September 2005; accepted 12 September 2005
Available online 6 October 2005

Abstract
A general reconstruction method is described which simulates point patterns possessing prescribed summary characteristics, which
are free of explicit model conditions. The characteristics are for instance the intensity, the L-function, the spherical contact distribution
function and the kth nearest neighbour distance distributions. The use of the statistical reconstruction method is demonstrated on
both a theoretical and practical example.
© 2005 Elsevier B.V. All rights reserved.

Keywords: Conditional simulation; kth Nearest neighbour distance distribution; L-Function; Spherical contact distribution; Reconstruction; Point
process

1. Introduction

Simulation of spatial point processes is an important task in spatial statistics, see MZller and Waagepetersen (2003).
It is needed for visualization, testing statistical procedures for point processes, model tests and parameter estimation.
Very often, as also in the present paper, it is assumed that the point process considered is stationary and isotropic, i.e.
that its distribution is motion invariant. This paper considers the planar case and rectangular simulation windows W,
which are not serious limitations.
The usual way leading to simulations in applied point process statistics is as follows. In some window Wobs a point
pattern
is observed which ‘looks stationary’. The window of observation Wobs can be of arbitrary shape and in
particular different from the simulation window W. The statistician assumes that
is a sample of a stationary point
process and estimates some non-parametric summary characteristics (s.c.) such as intensity, K- or L-function, pair
correlation function, kth nearest neighbour distance distribution function (d.f.) and spherical contact d.f. (also called
empty space function), see Cressie (1993), Diggle (2003), Stoyan and Stoyan (1994) and Stoyan et al. (1995). These
s.c.’s help to find a suitable model, for example a Cox process, a cluster process or a Gibbs process. Then with estimated
parameters the models are fitted to the data and tested. If accepted, then simulation can start. This approach needs a lot
of probabilistic and statistical knowledge, but it has the advantage to lead to elegant models with a small number of
parameters and supports building physical or biological models.
The approach of this paper, which is called reconstruction, is different. Its principal aim is simulation and not
derivation of a point process model with a small number of interpretable parameters. Reconstruction can be used

∗ Corresponding author. Tel.: +852 3411 2526; fax: +852 3411 5811.
E-mail address: [email protected] (A. Tscheschel).

0167-9473/$ - see front matter © 2005 Elsevier B.V. All rights reserved.
doi:10.1016/j.csda.2005.09.007
860 A. Tscheschel, D. Stoyan / Computational Statistics & Data Analysis 51 (2006) 859 – 871

to generate irregular point patterns which have s.c.’s close to estimates obtained statistically from the observation
window. Such patterns may serve for visualization (patterns in larger windows with shapes different from that of Wobs ),
for various kinds of non-parametric statistics (investigation of the behaviour of non-parametric s.c., determination of
the window size necessary for prescribed precision of intensity or s.c. estimation) and for constructing neighbours for
points close to the window’s edge.
When such a ‘reconstructed’ pattern is obtained, some measure of quality may be useful. It may consist in comparing
s.c.’s which were not used in the reconstruction, for the original and the simulated patterns. Of course, behind this
method there is also some point process model, but it is not needed here explicitly. The corresponding model parameters
are the numerical values of the s.c.’s. The reconstruction is based on some optimization method similar to simulated
annealing. So its use can be seen as a further application of the modern optimization heuristics to modelling problems
as discussed in Winker and Gilli (2004).
This method has the potential to become popular in spatial statistics. Variants of this methods have some tradition
in statistical physics, where they have been used for simulating random closed sets, typically as models for porous
media; see the exposition in Torquato (2001) and the references given there. In the diploma dissertation Tscheschel
(2001), written under the supervision of the second author and Gareth O. Roberts, this method was first-time used for
point processes. The paper Pommerening (2006), which uses another variant of the method in a forestry context, led
the authors to write this paper.
The application of the reconstruction method is demonstrated here by means of an example which is not typical
for real statistical applications, but which enables a very good quality test and may convince the readers of this paper
of the power of the approach. Instead of taking some empirical point pattern and reconstructing it, where never full
certainty can be reached (since the model behind the data, if any, is unknown), here a well-known point process is used,
a particular Neyman–Scott cluster process. Its theoretical s.c.’s play the role of the starting s.c.’s. In the reconstructions,
however, the knowledge of this particular underlying point process has not been employed—only some of its s.c.’s
have been used.
As the authors believe, the experience with the various s.c. as presented in this paper may be useful also for other
cases, since the results are very plausible. To satisfy readers who are more practically oriented, also the famous amacrine
cell pattern studied e.g. in Diggle (2003) is reconstructed. The method can be easily extended to conditional simulation,
where the positions of the points in some region are fixed while points in the other regions have to be reconstructed.
Finally, also the case of marked point processes is a straightforward application of the reconstruction method.

2. Theory

Let  denote an arbitrary stationary planar point process. In practical applications, the point process  itself is
usually not given in a constructive way but  may be accurately described by a number of s.c.’s. As mentioned above,
the usual way to simulate samples from  is based on an appropriate theoretical model, which is however in general
quite difficult to find. This task is even the more difficult the more precise information about  is available because
then a statistical test will reject each simple model.
In the reconstruction method presented in this paper, one does not need to specify a specific theoretical model for
. Instead, the reconstruction method only assumes knowledge of some prescribed s.c.’s, which are usually estimates
from observed realizations of . They are assumed to be the theoretical s.c.’s of the unknown point process . Those
s.c.’s should represent sufficient information about . In particular, if  is anisotropic then some s.c.’s should describe
the anisotropy in some way.
Let the prescribed
s.c.’s be some numerical s.c.’s ni for i = 1, 2, . . . , I and some functional s.c.’s fj (r) for j =
1, 2, . . . , J and r ∈ 0, Rj , where the Rj are suitably chosen limits. Examples for numerical s.c.’s are the intensity
 or so-called structural indexes as widely used in forestry, see von Gadow and Hui (2002). Examples for functional
s.c.’s are the K- or L-function, the kth nearest neighbour distance d.f.’s Dk and the spherical contact d.f. Hs .
Now let ñi (
) and f˜j (r;
) denote statistical estimates of ni and fj (r) by some appropriate method for the observed
point pattern
in Wobs . The basic idea of the reconstruction algorithm is to simulate a point pattern  in the sampling
window W which has estimated s.c.’s n̂i () and fˆj (r; ) close to the prescribed s.c.’s ñi (
) and f˜j (r;
), i.e.

n̂i () ≈ ñi (
) for i = 1, . . . , I

 A. Tscheschel, D. Stoyan / Computational Statistics & Data Analysis 51 (2006) 859 – 871 861

and


fˆj (r; ) ≈ f˜j (r;
) for r ∈ 0, Rj and j = 1, . . . , J .

The deviations of the estimated s.c.’s of a simulated pattern  in W from the prescribed s.c.’s are in the following
combined within one common deviation E() given by
I
 J

E() = Eni () + Efj () (1)
i=1 j =1

with

2
Eni () = ñi (
) − n̂i () for i = 1, 2, . . . , I (2)

and
 Rj  2
Efj () = f˜j (r;
) − fˆj (r; ) dr for j = 1, 2, . . . , J . (3)
0

Now the reconstruction method can be formulated as an optimization problem: The configuration space is the set of
all finite point patterns  within the simulation window W. The objective function is the deviation E() which is to
minimize. Of course, this optimization problem is not easy to solve because the configuration space is of high dimension
and the computation of the objective function E() may be time-consuming. Furthermore, the configuration space is
expected to contain not only a single global minimum but rather a high number of local minima  with small deviation
E(). However, the latter is not a disadvantage but rather an advantage. The optimization routine will randomly select
one of these local minima yielding different results in each simulation, which depends on the start configuration. In the
present paper, the following algorithm is used:
Start with an arbitrary pattern 0 in W, for example with a sample of a binomial point process. This pattern has
n points where n is the rounded value of  · 2 (W ) where  is the prescribed intensity and 2 denotes the Lebesgue
measure. The number of points n will not be changed during the reconstruction; thus it is guaranteed that the final
pattern has the right intensity.
In the following reconstruction steps always one point of the present pattern is randomly chosen and a candidate
point is created which will perhaps replace the chosen point. Let s be the point pattern after the step s, and let s be the
pattern without the chosen point and with a candidate uniformly placed in the window W. ‘Randomly chosen’ means
that every of the n points of s has the same chance 1/n to be selected. Then the successor s+1 of s is determined
by
   
s if E s
E s ,
s+1 = (4)
s else.
The iteration step from s to s+1 is carried out until either a maximum number of iterations steps have been carried
out or at least m iteration steps have been performed and
 
E s−m − E s < 

for some small  > 0 andbig m. 

Let denote
Es = E s − E s . It is possible tomodify the  main simulation step so that with some small
probability p (
Es ) also candidates s are accepted if E s > E s . Then the simulation is analogous to that of a
Gibbs process by means of the Metropolis–Hastings algorithm, see MZller and Waagepetersen (2003). If
Es > 0 then
a typical choice for p (
Es ) is of the form


Es
p (
Es ) = exp −
Ts
for some chosen Ts > 0 which can be interpreted as temperature at time s. Thus, this modified simulation could also
be called low temperature sampling.
862 A. Tscheschel, D. Stoyan / Computational Statistics & Data Analysis 51 (2006) 859 – 871

If the temperatures Ts are decreasing for increasing s, then the simulation follows the spirit of the simulated annealing
algorithm, Kirkpatrick et al. (1983). The simulated annealing algorithm has, in general, better convergence properties
than the above described improvements-only algorithm. The authors’ experience however showed that the deviation
E() of the final pattern  is always very small if W is big enough and the prescribed s.c.’s are not contradictory.
A combination of low temperature sampling and simulated annealing would be theoretically preferable against the
improvements-only method. In practical applications, however, the temperatures Ts introduce a set of parameters which
is difficult to choose.
Note that the results of the reconstruction algorithm can qualitatively differ if deficient s.c.’s are used for the
reconstruction. It is well-known that there are different point process models which have the same intensity and further
equal s.c.’s, for example L(r), see Baddeley and Silverman (1984).
In the algorithm described above, the estimators n̂i () and fˆj (r; ) were not yet specified in particular. However
these estimators should be chosen carefully as they directly determine the deviation E(). Since usually the point
patterns are simulated using periodic boundary conditions, it is natural to use also periodic boundary conditions in the
estimation as explained in the following section.

3. Summary characteristics and their estimation

This paper uses four different s.c.’s, which are widely used in practice, see for instance Cressie (1993), Diggle (2003),
Stoyan and Stoyan (1994), Stoyan et al. (1995) and Schladitz and Baddeley (2000). The first one is the L-function—a
normalization of the K-function. It is given in the planar case by

K(r)
L(r) = for r 0.

The kth nearest neighbour distance d.f.’s Dk (r) form a sequence of further s.c.’s. Dk (r) is, for k = 1, 2, 3, . . ., the d.f.
of the distance from the ‘typical’ point to its kth nearest neighbour, see also Stoyan and Stoyan (1994, p. 267).
For any point process with intensity , the K-function and the kth nearest neighbour distance d.f.’s Dk (r) are related
by
∞

K(r) = Dk (r) for r 0. (5)
k=1

So the sequence of all functions Dk (r) contains more information than the K- or L-function.
The third s.c. used in this paper is the spherical contact d.f. Hs (r). It is given by

Hs (r) = 1 − Prob( ∩ b(o, r) = ∅) for r 0

or equivalently by

Hs (r) = E2 b(o, r) ∩ [0, 1]2 for r 0,

where b(o, r) denotes the closed ball around the origin o with radius r. The symbol  denotes the Minkowski—addition
of two sets and is defined by

AB = {x + y : x ∈ A, y ∈ B} for A, B ⊂ R2 .

As an example for a third-order quantity, the characteristic T (r) in Schladitz and Baddeley (2000) is considered as a
fourth s.c. in this paper.
The estimation of the functions K(r), Dk (r), Hs (r) and T (r) is a crucial part of the reconstruction method. First,
these functions need to be estimated accurately from the given data. There are many publications dealing with that
subject.
For the practical example in Section 5, both the kth nearest neighbour distance d.f. Dk (r) and the spherical contact
d.f. Hs (r) were estimated using the border method, see for instance Baddeley (1999, p. 54). K(r) was estimated using
 A. Tscheschel, D. Stoyan / Computational Statistics & Data Analysis 51 (2006) 859 – 871 863

Ripley’s edge correction, see for instance Diggle (2003, p. 50). T (r) was estimated using the translation correction,
see Schladitz and Baddeley (2000).
Estimates are needed furthermore in each step of the reconstruction algorithm, where all used s.c.’s have to be
determined for the candidate pattern s . Let xd and yd denote the width and height of the rectangular sampling window
W. Then, by applying periodic boundary conditions, the following natural estimators were used for r < 21 min (xd , yd ):
xd · yd  
K̂(r; ) = 1( x − y < r), (6)
n(n − 1)
x∈
y∈
y =x

⎛ ⎞
1  ⎜  ⎟
D̂k (r; ) = 1⎜
⎝ k
1( x − y < r)⎟
⎠, (7)
n
x∈ y∈
y =x

xd2 · yd2
T̂ (r; ) =
n(n − 1)(n − 2)
  
× 1(max{ x − y , x − z , y − z } < r), (8)
x∈ y∈ z∈
z =x
y =x z =y

and
2 (W ∩ [b(o, r)])
Ĥs (r; ) = . (9)
xd · y d
The Euclidean distance · is used in (6)–(8) in the spirit of periodic boundary conditions. Periodic boundaries have
also to be applied when computing the area measure 2 of the dilated point pattern b(o, r). Its analytical calculation
is not very difficult for such systems of discs if one applies Gauss’ formula as demonstrated by Brodatzki and Mecke
(2002).
Although all three estimators given in (6)–(9) can be computed relatively fast, their estimation can still be made more
efficient using the fact that the candidate pattern s only differs in one point from s . This is particularly easy for the
K-function and the kth nearest neighbour distance d.f.’s Dk (r), where instead of a double sum only a single sum related
to the point changed needs to be calculated anew. For the estimation of the spherical contact d.f. Hs (r), the Lebesgue
measure 2 needs only to be computed anew in the surrounding area of the point deleted and added.
The L-function was estimated using an estimate of the K-function by

K̂(r; )
L̂(r; ) = for r 0.


4. Reconstructing a Neyman–Scott process

The following presents an example for the application of the reconstruction algorithm for a well-known model, for
which many s.c.’s are analytically known so that a precise quality control is possible. The aim is to demonstrate the
power of the reconstruction method and to show which s.c.’s should be used in the reconstruction method.
The example follows Tscheschel (2001) in reconstructing particular cluster processes in a 2 × 2 square, namely a
special class of Neyman–Scott processes, see Stoyan et al. (1995). Such processes are Poisson cluster processes: there
is a homogeneous Poisson process of parent points with intensity p . Around each of the parent points a cluster of
independent daughter points is scattered. The number of points per cluster is random with mean c and the probability
to have k points is pk . Then the intensity of the cluster process is

 = p · c. (10)
864 A. Tscheschel, D. Stoyan / Computational Statistics & Data Analysis 51 (2006) 859 – 871

Let the daughter points be independently and uniformly distributed in a disc of radius R around the parent point. In this
case the d.f. F (r) of the distance between two randomly chosen points of the same cluster is given by
⎛
1 ⎝ 3 r r
F (r) = 4R arcsin + 4Rr 2 arccos
2R 3 2R 2R
 ⎞
r 2
⎠
−r 2R 2 + r 2 1−
4R 2

for 0
r < 2R.

If the parent point is at the origin o, the cluster is called the representative cluster and denoted by Xo .
In all following practical examples, the cluster radius R = 0.03 was used.
Using the notation above, the K-function of a Neyman–Scott process has the general form
∞
F (r) 
K(r) = r 2 + k(k − 1)pk /c for r 0, (11)

k=2

see Stoyan et al. (1995, p. 159).

The intensity of the cluster process is now assumed to be  = 100 and the value of
∞

k(k − 1)pk /c
k=2

as 3. Three possible cases are

(a) p = 25 with p4 = 1 and pk = 0 for k = 4,

(b) p = 70 with p1 = 20
21 , p10 = 21 and pk = 0 for k ∈
1
/ {1, 10}, and
3k
(c) p = 100
3 with pk = k! e−3 for k = 0, 1, . . . .

Case (a) corresponds to the fixed number 4 of cluster points in every cluster, while in the case (b) there is one point per
cluster mostly but with small probability there are 10. The case (c) is a particular Matérn cluster process (Stoyan et al.,
1995) with a mean number of three cluster points. All these cluster processes share the same intensity and the same K-
respectively L-function.
Let (a) and (b) denote the point processes corresponding to the cases (a) and (b). Both processes will be subject
of the further investigations. Fig. 1 shows simulated patterns of (a) and (b) .
Our aim was to reconstruct point patterns that are statistically close to the processes (a) , respectively, (b) . In the
following, three experiments allow a comparison of the effects of different s.c.’s used during the reconstruction. In
these experiments, one or two of the following three s.c.’s are considered: the L-function, the spherical contact d.f.
Hs (r) and the kth nearest neighbour distance d.f.’s Dk (r). As a measure of the quality of the final reconstructions
also a third-order quantity, namely the characteristic T (r) in Schladitz and Baddeley (2000), was considered. In the
implementation, the integral in (3) is replaced by a sum as


100
2
EL () = r L(i · r) − L̂(i · r; ) , (12)
i=1


100
2
EHs () = r Hs (i · r) − Ĥs (i · r; ) , (13)
i=1


41 
100
2
EDk () = r Dk (i · r) − D̂k (i · r; ) (14)
k=1 i=1
 A. Tscheschel, D. Stoyan / Computational Statistics & Data Analysis 51 (2006) 859 – 871 865

(a) (b)

Fig. 1. Samples of two different Neyman–Scott processes in a 2 × 2 square with the same intensity and L-function as explained in the text. In the left
hand case all clusters consist of 4 points, while on the right hand there are usually one-point clusters and with small probability 10-point clusters:
(a) p = 25 with p4 = 1; (b) p = 70 with p1 = 20 1
21 and p10 = 21 .

and


100
2
ET () = r T (i · r) − T̂ (i · r; ) (15)
i=1

with r = 0.0025.
The theoretical formula for L(r) for the two processes (a) and (b) could easily be obtained from Eq. (11). The
functions Dk (r) were calculated numerically, using the fundamental formula for the Palm distribution of Poisson cluster
processes, see Stoyan et al. (1995, p. 158). Statistical estimates were used for the theoretical function for T (r). The
formula for Hs (r) is known from the theory of Boolean models (Stoyan et al., 1995) and it holds
  
Hs (r) = 1 − exp −p · E2 Xo b(o, r) for r 0. (16)

In the case of the process (a) we have

   4 
 R+r
r,R (t)
E2 Xo b(o, r) = 2 t 1− 1− dt for r 0 (17)
0 R 2

while the process (b) yields

   10 
 R+r
20 1 r,R (t)
E2 Xo b(o, r) = · r 2 + · 2 t 1− 1− dt (18)
21 21 0 R 2

for r 0, where R is the cluster radius and

r,R (t) = 2 (b(o, R) ∩ b(t, r)) for r, R, t 0

and t being any vector of length t. The formula for r,R (t) is for instance given in Stoyan and Stoyan (1994, p. 365).
All following simulations are carried out using n = 400 points in the window W = [0, 2] × [0, 2], i.e. a square ofside
length 2. The random start configuration 0 used for all reconstructions is shown in Fig. 2(a). The deviations E 0
defined in Eqs. (12)–(15) for the start configuration 0 are shown in Tables 1 and 2 in the column denoted by 0.
866 A. Tscheschel, D. Stoyan / Computational Statistics & Data Analysis 51 (2006) 859 – 871

(a) (b)

Fig. 2. (a) The start configuration consisting of 400 points randomly scattered in the 2 × 2 square. (b) Result of the reconstruction in the 2 × 2 square
based only on intensity and L-function. Compare with Fig. 1.

Table 1
Deviations of four different s.c.’s with respect to the point process (a) for the start configuration (0) and the final configurations of three experiments
(I–III) as explained in the text

0 I II III

EL () 3.2 × 10−4 3.6 × 10−9 6.7 × 10−9 6.0 × 10−8

EHs () 8.7 × 10−3 1.6 × 10−3 2.0 × 10−8 8.6 × 10−6
EDk () 2.1 × 10−1 2.5 × 10−2 1.0 × 10−2 1.6 × 10−4
ET () 3.8 × 10−6 2.8 × 10−8 2.8 × 10−8 8.9 × 10−10

The deviation of the s.c.’s used in the respective reconstruction is shown in bold.

Table 2
As Table 1 but with respect to the point process (b)

0 I II III
−9 −9
EL () 3.2 × 10−4 3.6 × 10 8.4 × 10 5.8 × 10−8
EHs () 6.4 × 10−4 1.0 × 10−3 4.9 × 10−9 1.1 × 10−6
EDk () 1.7 × 10−1 2.3 × 10−2 6.5 × 10−3 1.9 × 10−4
ET () 5.5 × 10−6 8.7 × 10−8 1.5 × 10−8 7.4 × 10−9

In the first reconstruction experiment, only the L-function is used as s.c., beside the intensity , i.e. it is I = 1 and
J = 1 with

f1 (r) = L(r).

Clearly, the L-function does not separate the two processes (a) and (b) . However, it is interesting to observe to which
geometry the algorithm leads.
The final pattern after s = 106 steps is shown in Fig. 2(b). It looks like a cluster process. Of course, the number of
points per cluster is random, and perhaps one could say that the pattern is somehow ‘between’ the patterns of Fig. 1.
The deviations E s of the final pattern s of the four considered s.c.’s are shown in Tables 1 and 2 in the column
denoted by I. The deviation EL s became very small, so the algorithm worked well.
 A. Tscheschel, D. Stoyan / Computational Statistics & Data Analysis 51 (2006) 859 – 871 867

(a) (b)

Fig. 3. Results of the reconstructions in the 2 × 2 square based on the intensity, the L-function and the spherical contact d.f. Compare with Fig. 1.

In the second experiment, the functional E() was refined. Because the L-function does not separate the point
processes (a) and (b) , a further functional s.c. f2 (r) was introduced, namely the spherical contact d.f. Hs (r), i.e. it
is I = 1 and J = 2 with
f1 (r) = L(r) and f2 (r) = Hs (r).
Fig. 3 shows the result of reconstruction with two functional s.c.’s obtained again after s = 106 iteration steps. The
clusters in Fig. 3(a) are now similar to that shown in Fig. 1(a), where every cluster has 4 points, but, of course, it cannot
be expected that every cluster in the final reconstructed pattern has exactly 4 members. It is hard anyway to identify
single clusters because some may overlap. On the other side, in Fig. 3(b), the ‘clusters’ consist mostly of only one point
but some clusters contain
 a high number of points—similar to the process shown in Fig. 1(b).
The deviations E s of the final  pattern s of
 the four considered s.c.’s are shown in Tables 1 and 2 in the column
denoted by II. The deviations EL s and EHs s of both s.c.’s used during the reconstruction became very small 
 pattern. Thus the reconstruction algorithm worked very successful in both cases. The deviations EDk s
in both final
and ET s also became relatively small at the same time. This indicates that the two s.c.’s L(r) and Hs (r) already
contain a lot of statistical information about the underlying point process.
In the third experiment, instead of L(r) and Hs (r) a sequence of kth nearest neighbour distance d.f.’s Dk (r) were used.
The reconstruction algorithm was used now to reconstruct the k-nearest neighbour distance d.f.’s D1 (r), . . . , D41 (r),
i.e. it is I = 1 and J = 41 with
fk (r) = Dk (r) for k = 1, 2, . . . , 41.
For each of the processes (a) and (b) , one final pattern of the reconstruction process after s = 106 steps is shown in
Fig. 4. The precision of the reconstruction is better than in the simultaneous use of L-function and spherical contact d.f.
The reason is that the sequence of functions Dk (r) includes not only information about the L-function as mentioned
above but also precise information about the geometry of the clusters. For example, for the point process (a) which
contains four-point clusters only, it is
Dk (r) = 1 for r 2R and k = 1, 2, 3.
This leads to reconstructions in which for almost every point the next three neighbours are within distance 2R.
Fig. 5 shows D1 (r), D2 (r), . . . , D41 (r) for both cluster processes. The amount of structural information about the
point processes contained
 in the sequence D1 (r), D2 (r), . . . becomes obvious.
The deviations E s of the final pattern s of the four considered s.c.’s are shown in Tables 1 and 2 in the column

denoted by III. The deviation EDk s is for both cases (a) and (b) about 1000 times smaller than the deviation
868 A. Tscheschel, D. Stoyan / Computational Statistics & Data Analysis 51 (2006) 859 – 871

(a) (b)

Fig. 4. Results of the reconstructions in the 2 × 2 square based on the intensity and the kth nearest neighbour distance d.f.’s D1 (r), . . . , D41 (r).
Compare with Fig. 1.

1 1

0.8 0.8

0.6 0.6
Dk(r)
Dk(r)

0.4 0.4

0.2 0.2

0 0
0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25
(a) r (b) r

Fig. 5. Theoretical kth nearest neighbour distance d.f.’s D1 (r), . . . , D41 (r) for the point processes (a) and (b) shown in Fig. 1. D1 (r) is the
respectively top most curve; D41 (r) the respectively lowest curve.

EDk 0 of the start configuration. Thus also the reconstruction in the third experiment worked well. The other three
s.c.’s L(r), Hs (r) and T (r) serve as measures of quality as they were not used during the reconstruction. Although

the deviation EHs s is larger than in the second experiment, it is quite small compared to the deviation EHs 0
 It is not surprising that the deviation EL s became very small because of the relation (5).
of the start configuration.
Also the deviation ET s is remarkably small. This indicates that also the third-order behaviour is well hit by the
simulated patterns, at least in the sense of T (r).
An important criterion for practical use is the speed of such reconstructions. On a recent desktop computer with a
single 2.4 GHz processor, about 22 iteration steps per second were performed in the reconstructions using L(r) and
Hs (r), whereas about 110,000 iteration steps per second were performed when using the L-function only. In the third
experiment, in which 41 nearest neighbour distance d.f.’s needed to be estimated in each step, about 14,000 iteration
steps per second were carried out. These numbers of steps per second are practically independent of the size of the
observation window and thus of the number of simulated points.

5. Reconstructing an empirical point pattern

In practical applications one has to reconstruct empirical patterns, and instead of the theoretical functions estimates
from the given pattern can be used.
 A. Tscheschel, D. Stoyan / Computational Statistics & Data Analysis 51 (2006) 859 – 871 869

(a) (b)

Fig. 6. (a) Amacrine cell point pattern: both ‘on’ and ‘off’ cells. (b) Reconstruction of the Amacrine cell point pattern using the first 17 nearest
neighbour distance d.f.’s Dk (r).

In order to demonstrate also such an application, Fig. 6(b) shows a reconstruction of the famous amacrine cell
point pattern shown in Fig. 6(a), which is thoroughly discussed in Diggle (2003) and Diggle et al. (2005). This is a
bivariate pattern of ‘on’ and ‘off’ cells, but it is considered here as a univariate pattern, since this paper restricts itself
to non-marked point processes. While that pattern
of 294 points was given in a rectangular window Wobs of side
lengths 1.601 and 1, the size of the reconstruction window W is now 2 × 2 and the point number is n = 734.
Statistical estimates L̃(r;
), D̃k (r;
), H̃s (r;
) and T̃ (r;
) for L(r), Dk (r), Hs (r) and T (r) obtained from the
observed point pattern
in Wobs are interpreted as prescribed values. The reconstruction is based on the kth nearest
neighbour distance d.f.’s D1 (k), . . . , D17 (k), while the other s.c.’s L(r), Hs (r) and T (r) serve as measures of quality.
Thus it is I = 1 and J = 17 with

fk (r) = Dk (r) for k = 1, 2, . . . , 17.

The integral in (3) is replaced by a sum as


100
2
EL () = r L̃(i · r;
) − L̂(i · r; ) ,
i=1

100
2
EHs () = r H̃s (i · r;
) − Ĥs (i · r; ) ,
i=1

17 
100
2
EDk () = r D̃k (i · r;
) − D̂k (i · r; ) ,
k=1 i=1

and

100
2
ET () = r T̃ (i · r;
) − T̂ (i · r; )
i=1

with r = 0.0015.
The deviations EL (), EDk (), EHs () and ET () of both the random start configuration 0 and the final config-
uration s obtained after s = 106 steps are shown in Table 3. The reconstruction method worked very well for both
Dk (r) and the other three s.c.’s L(r), Hs (r) and T (r).
870 A. Tscheschel, D. Stoyan / Computational Statistics & Data Analysis 51 (2006) 859 – 871

Table 3
As Table 1 but with respect to the amacrine point pattern

0 I

EL () 1.8 × 10−5 6.1 × 10−8

EHs () 4.7 × 10−4 7.5 × 10−8
EDk () 3.0 × 10−2 2.2 × 10−5
ET () 7.6 × 10−9 1.7 × 10−11

(a) (b)

Fig. 7. Reconstruction of the Amacrine cell point pattern conditioned on the observed point pattern in the respectively outlined subwindow. The
subwindow used is (a) the whole observed window and (b) a disc-shaped subset of the observed window shown in Fig. 6(a).

In order to show also the potential of the reconstruction method for conditional simulation, the amacrine pattern is
conditionally reconstructed in two ways. Fig. 7(a) shows the result of a reconstruction of the amacrine pattern in a
window W of size 2 × 2 given the original 294 points in the observed window Wobs of size 1.601 × 1, i.e. the points
within Wobs are fixed while the remaining 440 points were reconstructed in W \Wobs . The fixed points in Wobs interact
with points in W \Wobs .
The conditional reconstruction is not restricted to rectangular observation windows. Thus, Fig. 7(b) shows another
conditional simulation, in which now only the points within a disc of radius 0.5 inside Wobs are fixed.

6. Discussion

The reconstruction method is an efficient way to produce simulations of any stationary planar point process if only
some of its s.c.’s are given. This method is particularly precise and efficient when using the kth nearest neighbour
distance d.f.’s Dk (r). The functions Dk (r) can be estimated very fast and reconstructions can be performed within
seconds or some few minutes. The spherical contact d.f. Hs (r) can deal as a measure of quality for the final pattern.
Except of k = 1, the kth nearest neighbour distance d.f.’s Dk (r) have been considered only rarely for statistically
describing point patterns until now. This paper encourages their use also in other applications since they obviously
contain valuable information.
The potential of the method is higher than shown in this paper. It can be applied also to reconstruct marked point
processes, where marked versions of the kth nearest neighbour distance d.f.’s Dk (r) may be useful. A further application
is conditional simulation. By means of this method, points outside the original observation window can be simulated.
These artificial points can be used for the estimation of local characteristics such as nearest neighbour distances of
single points, see Pretzsch (2002).
 A. Tscheschel, D. Stoyan / Computational Statistics & Data Analysis 51 (2006) 859 – 871 871

Acknowledgements

We thank the anonymous referees for helpful comments and suggestions. The first author thanks for financial support
by grants from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project Numbers
HKBU2048/02P and HKBU200503).

References

Baddeley, A.J., 1999. Spatial sampling and censoring. In: Barndorff-Nielsen, O.E., Kendall, W.S., van Lieshout, M.N.M. (Eds.), Stochastic Geometry:
Likelihood and Computation. Chapman & Hall/CRC, London, pp. 37–78.
Baddeley, A.J., Silverman, B.W., 1984. A cautionary example on the use of second-order methods for analyzing point patterns. Biometrics 40,
1089–1094.
Brodatzki, U., Mecke, K., 2002. Simulating stochastic geometries: morphology of overlapping grains. Comput. Phys. Comm. 147, 218–221.
Cressie, N.A.C., 1993. Statistics for Spatial Data. Review ed. Wiley, New York.
Diggle, P.J., 2003. Statistical Analysis of Spatial Point Patterns. Second ed. Arnold, Paris.
Diggle, P.J., Eglen, S.J., Troy, J.B., 2005. Modelling the bivariate spatial distribution of amacrine cells. In: Baddeley, A.J., Gregori, P., Mateu, J.,
Stoica, R., Stoyan, D. (Eds.), Case Studies in Spatial Point Process Modeling. Lecture Notes in Statistics, Vol. 185, Springer, Berlin, pp. 215–234.
Kirkpatrick, S., Gelatt Jr., C.D., Vecchi, M.P., 1983. Optimization by simulated annealing. Science 220, 671–680.
MZller, J., Waagepetersen, R.P., 2003. Statistical Inference and Simulation for Spatial Point Processes. Chapman & Hall/CRC, Boca Raton.
Pommerening, A., 2006. Evaluating structural indices by reversing forest structural analysis. Forest Ecology and Management (in print).
Pretzsch, H., 2002. Grundlagen der Waldwachstumsforschung. Blackwell, Berlin.
Schladitz, K., Baddeley, A.J., 2000. A third order point process characteristic. Scand. J. Statist. 27, 657–671.
Stoyan, D., Stoyan, H., 1994. Fractals, Random Shapes and Point Fields. Wiley, Chichester.
Stoyan, D., Kendall, W.S., Mecke, J., 1995. Stochastic Geometry and its Applications. Second ed. Wiley, Chichester.
Torquato, S., 2001. Random Heterogeneous Materials. Springer, New York.
Tscheschel, A., 2001. Reconstruction of random porous media. Diploma dissertation, TU Bergakademie Freiberg, Germany.
von Gadow, K., Hui, G.Y., 2002. Characterizing forest spatial structure and diversity. In: Björk, L. (Ed.), Sustainable Forestry in Temperate Regions.
Proceedings of the SUFOR International Workshop, Lund University, pp. 20–30.
Winker, P., Gilli, M., 2004. Applications of optimization heuristics to estimation and modelling problems. Comput. Statist. Data Anal. 47, 211–223.