- Email: [email protected]
Spatial pattern analysis of second-phase particles in composite materials
Materials Science and Engineering A356 (2003) 245 /257 www.elsevier.com/locate/msea
Spatial pattern analysis of second-phase particles in composite materials J.D. Scalon a,*, N.R.J. Fieller b, E.C. Stillman b, H.V. Atkinson c a
Department of Biomedical Engineering, Federal University of Sa˜o Joa˜o Del Rey, Prac¸a Dom Helve´cio 64, Sa˜o Joa˜o Del Rey MG 36301-160, Brazil b Department of Probability and Statistics, University of Sheffield, P.O. Box 597, Sheffield S10 2UN, UK c Department of Engineering Materials, University of Sheffield, P.O. Box 600, Sheffield S1 4DU, UK Received 30 July 2002; received in revised form 12 February 2003
Abstract The spatial distribution of second-phases is one of the most important factors affecting the mechanical properties of composite materials and therefore, the quantitative characterisation of such distributions is of prime importance in materials science. In this paper, we present a suitable procedure for analysing spatial distribution of particle centres held in three planar sections of an aluminium silicon carbide (Al/SiCp) composite material. We suggest that the spatial pattern analysis can be performed by following five consecutive steps: visualising spatial distribution, testing against complete spatial randomness, choosing and fitting appropriate stochastic models and testing goodness-of-fit. These steps show how statistical methods that rely on functional pattern descriptors and Monte Carlo simulations can be useful for characterising spatial distribution of the reinforcing particles in composite materials. # 2003 Elsevier Science B.V. All rights reserved. Keywords: Particle distribution; Reinforced aluminium; Spatial analysis; Simulation
1. Introduction It is well known that the mechanical properties (e.g. fracture toughness) of particulate composite materials depend not only on the shape and the volume fraction, but also on the spatial distribution of the reinforcing particles held in the matrix [1]. For example, there is evidence that surface cracking associated with clustered particles appear to be more prevalent than those associated with less clustered composite materials [1,2]. Thus, the quantitative characterisation of the spatial distribution of particles is of interest not only for a better understanding of the relationship between inclusions and mechanical behaviour, but also for better control of the production of composite materials. Consequently, there is an interest among material researchers in developing methods to provide a more rigorous quantitative analysis of the spatial distribution
* Corresponding author. Tel.: /55-32-3379-2505; fax: /55-323379-2328. E-mail address: [email protected] (J.D. Scalon).
of second-phase particles in two-dimensional distributed multi-phase systems. During the last few years, statistical methods have been appearing, to a certain extent, in materials literature to provide analysis of spatial distribution of particles in composite materials. Methods, such as quadrat count [3 /6], nearest neighbour distances [4,7,8], Dirichlet tessellation [2,4,8 /14], spatial pattern descriptors [15 /17] and correlation function [18] have been applied, basically, to test departure from randomness of the spatial distribution of the reinforcing particles. We have observed that some methods of spatial statistics were completely ignored by the researchers, for example, the fitting of a stochastic model for spatial distribution of particles and the use of techniques based on computer intensive methods for providing hypothesis tests. Thus, the aim of the present paper is to further explore the use of methods of spatial statistics (empirical distribution functions of pattern descriptors, such as F -, G - and K -functions and Monte Carlo simulation) for characterising the spatial distribution of reinforcing particles held in planar sections of composite materials.
0921-5093/03/$ - see front matter # 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0921-5093(03)00138-2
246
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
We suggest that one can provide a suitable quantitative analysis for characterising the spatial distribution of the particles by following five consecutive steps: visualising spatial distribution, testing against randomness, choosing an appropriate stochastic model, fitting the model to the pattern and testing goodness-of-fit. The main aim of this approach is to provide assessments of deviation from randomness towards either regularity or clustering of the particles and to suggest some alternative models for the spatial distribution of the reinforcing particles.
2. Description of experimental data We have applied the statistical analysis to three metallographic samples (S1 /S3) of an aluminium silicon carbide composite material produced by the Department of Engineering Materials, University of Sheffield, UK, where the silicon carbide second-phase particles had a volume fraction equal to 11%. The metallographic sample areas were :/192 /288 mm. The three metallographic samples were placed on a computer controlled optical microscope stage analyser (Polyvar), which allowed fully automatic adjustment, focusing, positioning and scanning of the samples. The overall microscope magnification used was 600 times, yielding a pixel size of 0.375 mm. Thus, the Polyvar produced three digital images with area frames equal to 512 /767 square pixels. The two-dimensional digital images were analysed by using image-processing techniques to extract the coordinates (centres) of each particle within the images. The actual images were 512/767 pixels, out of which we used only the particles located in the left top square of 512 /512 pixels, yielding 419 particles in sample S1, 507 particles in sample S2 and 597 particles in sample S3. The 512 /512 images were transformed into patterns with unit square area to facilitate the analysis. Then, we have provided the quantitative analysis of the spatial distributions of the reinforcing particles by following the suggested consecutive steps.
3. Visualising spatial distribution of the reinforcing particles The first requirement in any data analysis is the ability to see the data being analysed. Visualising spatial patterns means mapping the particles held in the matrix of the composite material. A natural tool for visualising a spatial distribution of particles is to plot their location in the matrix as a dot map. This map provides an idea of the distribution of particles and shows any obvious patterns presented in the sample. The dot maps of the samples presented in Fig. 1(a), Fig. 2(a) and Fig. 3(a) give a visual impression of the distribution of particles
Fig. 1. Dot map of locations (a) and contour plot (b) of the nonparametric kernel estimates of intensity using quartic kernel estimator and t /0.14 of 419 particles in sample S1. Contour values shown are the quantitative measures of local intensity of particles.
along the matrix. One must observe that it is very difficult to judge the nature of the patterns based only on these figures, that is, whether there is a tendency for aggregation or regularity. A more advanced statistical tool is the use of a spatially smooth estimate of the way in which the intensity of particles (the expected number of particles per unit area) is varying in the matrix. Diggle [19] describes a sophisticated approach for obtaining a spatially smooth intensity of particles based on a kernel estimator. More formally, let (s1, s2, . . . , sn ) be the spatial locations of n particles in a bounded study region jAj. The intensity, l (s), at s is estimated by lˆt (s)
n X 1 (s si ) ; k 2 t i1 t
(1)
where t /0 is the bandwidth parameter that determines the amount of smoothing. Following Diggle’s sugges-
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
Fig. 2. Dot map of locations (a) and contour plot (b) of the nonparametric kernel estimates of intensity using quartic kernel estimator and t /0.14 of 507 particles in sample S2. Contour values shown are the quantitative measures of local intensity of particles.
tion [19], we have adopted t /0.68n0.2. The k (.) is a bivariate probability density function, known as the kernel, which is symmetric about the origin. A simple quartic function can be used as an alternative for k (.) and then the estimate of intensity lˆt (s) may be simply expressed as n X 3 h2i 2 ˆ lt (s) 1 ; (2) 2 t2 i1 pt where hi is the distance between the point s and the observed particle location si . The smoothed intensity estimates lˆt (s) can be plotted in a variety of ways, e.g. perspective and contour. These
247
Fig. 3. Dot map of locations (a) and contour plot (b) of the nonparametric kernel estimates of intensity using quartic kernel estimator and t /0.13 of 597 particles in sample S3. Contour values shown are the quantitative measures of local intensity of particles.
plots can help identify peaks in the resulting surface corresponding to possible locations of clusters. Contour plots of our samples are presented in Fig. 1(b), Fig. 2(b) and Fig. 3(b). We can observe that these plots allow identify locally high concentrations of particles (the one with the local intensity is much bigger than the number of particles in the sample) and reduced local concentrations (where the local intensity is much lower than the number of particles in the sample) in the samples. These figures show evidence of the presence of both high and low local concentrations of particles in both samples and therefore, there is evidence that the distribution of particles along the matrix is not random.
248
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
One can observe that the kernel estimate is essentially a sum of ‘bumps’ placed at observations. The kernel function determines the shape of the bumps, while the bandwidth parameter determines their width. When estimating intensity, the most important thing is the choice of an appropriate bandwidth value, while the choice of kernel function is of secondary importance since any reasonable kernel (e.g. Gaussian and quartic) will give close to optimal results [19]. If the bandwidth tends to zero, the intensity is estimated as a series of ‘spikes’ in the individual observations, while if the bandwidth becomes large, all detail is obscured. In practice, we can try different values of bandwidth in order to explore the variation of intensity at different scales. Two examples are presented (Fig. 4) where values of t /0.06 and t/0.16 reveal that the higher values ‘smooth’ the distribution rather too much, while the lower value of 0.06 gives a rather too ‘spiky’ impression
Fig. 4. Contour plots of the nonparametric kernel estimates of intensity, in sample S1, using quartic kernel estimator for different values of the bandwidth parameter (t ). Contour values shown are the quantitative measures of local intensity of particles.
of the spatial distribution of the particles in sample S1. Diggle [19] suggests using t /0.68n0.2 for estimating the intensity on a unit square or other more efficient methods that attempt to optimise the value of bandwidth parameter based on estimates of the mean square error. We observe that the maps presented in Figs. 1/4 may be sufficient to answer some questions about the spatial distribution of second-phase particles in samples of composite material. However, they did not allow a conclusive answer about the nature (regular or clustered) of the pattern and therefore, it is appropriate to test whether the particle distribution can be considered random.
4. Testing against complete spatial randomness A large number of statistical methods have been proposed for testing departure against the null hypothesis of completely spatial randomness (CSR) in point patterns. The first group of techniques (quadrat count methods) is based on the determination of the number of particles in small test regions in order to assess spatial particle distribution. The observed frequency distribution of the number of particles per quadrat is compared to a theoretical Poisson distribution. If the particles are randomly distributed, the observed distribution will coincide closely with the theoretical distribution. In general, the x2-test for goodness-of-fit is applied to assess the degree of agreement [3 /6]. One must observe that the choice of the region size is a crucial problem in quadrat methods and therefore they depend on the arbitrary choice of both the quadrat area size and quadrat shape. The second group of techniques is based on two approaches of nearest neighbour distances that require information about the locations of a sample of particle centres held in the composite material. The first approach is based on the distance between a randomly chosen particle and its nearest neighbour. The second approach consists in choosing a point at random and then the nearest distance between the point and particle in the pattern is determined. The aim, in both approaches, is to compare the observed distribution of the nearest neighbour distances and the theoretical distribution of distances under randomness (Poisson distribution). If the particles are randomly distributed, the observed distribution will coincide closely with the theoretical distribution [4,7,8]. The drawback of these techniques is that the reduction of complex point patterns to a single nearest neighbour summary statistic, calculated from locations of a sample of particles, results in a considerable loss of information. Thus, this sort of methods is not recommended when the information
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
about the location of all second-phase particle centres held in the matrix of composite material is available. The third and most widely used of all statistical methods for the analysis of spatial distribution of particles in composite materials is based on Dirichlet (Delaunay or Voronoi) tessellation. This method uses a geometrical construction for partitioning a field into a unique set of convex polyhedra (polygons), each of which is associated with and contains one of the particles. The polygons form a continuous cellular structure. Analysis of the characteristics of these polyhedra (e.g. area) provides information about the local environment of the particles [2,4,8 /12]. It is observed that Dirichlet-based methods are not effective in distinguishing cluster, regular and random patterns for high volume (area) fractions of particles [13]. Otherwise, Dirichlet-based methods may neglect some significant areas of space containing no events [14]. We note that the empirical distribution function for some variables related to Dirichlet tessellation (e.g. area of the polygon) are not known for either Poisson processes or a number of potentially useful classes of spatial models. This fact can limit the use of Dirichlet tessellation for either conducting formal significant tests against randomness or estimating parameters of spatial models. One of the least frequently used group of all statistical methods for the analysis of spatial distributions of particles in composite materials is based on spatial pattern descriptors, such as F -, G - and K -functions. These methods require a completely enumerated point pattern, that is, one needs information about the location of all particles held in the composite material. These functions have been used on a few works for the analysis of spatial distribution of inclusions in composite material [15 /17], but the full potential has not been totally explored. Thus, in this paper we are using tests based on these pattern descriptors. For practical purposes, each second-phase particle is treated as a point defined by its co-ordinates. Thus, our spatial data can be assumed to be a map of all particle locations in an essential planar region and therefore the analysis can be founded on the theory of spatial point patterns. Typically, an observed point pattern can be thought of as the outcome (a realisation) of a spatial point process. Because spatial point pattern data are usually a single realisation of a spatial point process, the additional assumptions of stationarity and sometimes isotropy, are often made to reduce the parameter space and to allow parameter estimation [20]. One can never hope to totally characterise the spatial point process, but one can investigate the behaviour of functional pattern descriptors, proposed by Ripley [21] and described in detail by Diggle [20] and Cressie [22], that represent important aspects of the process.
249
4.1. The K -function A convenient way of characterising spatial point patterns is provided by the K -function. The theoretical K -function of a stationary spatial point process is defined by K (t)/l 1E (number of further events within distance t of an arbitrary event). Under the hypothesis of CSR, the K -function is given by K (t)/ pt2. An unbiased estimator for K (t) that corrects edge effects is given by ˆ K(t)
n X n ½A½ X It (uij )
n2
i1 j1
wij
;
(3)
where uij is the distance between events i and j (i "/j), It (uij ) is an indicator function that takes the value 1 when uij is less than distance t, l is the intensity, n is the number of events in the analysed pattern of area jA j and wij is a weighting factor that represents the proportion of the circumference of the circle around event i , passing through the point j that lies within jAj. The K -function represent a natural and valuable tool for the description of a spatial point process, but it does not provide a complete description of the probability structure of a spatial point process. It is possible to find different processes that have identical K -functions [23]. For a complete description of a spatial point pattern it would be necessary to include other functional pattern descriptors, such as those based on nearest neighbour distances. 4.2. The F - and G -functions The F- function is based on the distance x , that is, the distance from an arbitrary chosen point (not event) to the nearest event. It is defined by the distribution function F (x )/P {distance from an arbitrary point to the nearest event is at most x }. The theoretical F function under the null hypothesis of CSR is given by F (x)1exp(lpx2 ) x 0:
(4)
The simplest suitable edge-corrected estimator of the F -function is provided by the equation n X
Ix (xi ; ri ) i1 ˆ F (x) n X Ix (ri )
x 0;
(5)
i1
where m is the number of sample points in the analysed pattern, xi denotes the distance from the ith chosen point to the nearest of the n events in the analysed pattern, Ix (ri ) is an indicator function that takes the value 1 when xi is /x , ri is the distance from each event to the nearest point on the boundary of jAj and Ix (xi ,ri )
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
250
is an indicator function that takes the value 1 when xi is /x and ri is ‘/x. We can use m sample points in a regular grid k /k to estimate the values of the F pffiffiffi function, where k $ n:/ The G -function is based on the distance y, that is, the distance between two nearest events. It is given by the expression G (y) /P {distance from an arbitrary event to the nearest other event is at most y }. The theoretical G function under the null hypothesis of CSR and an edgecorrected estimator for G (y) is provided by equations analogous to Eq. (4) and Eq. (5), respectively, substituting the distance x by the distance y. We suggest the use of F - and G -functions to test departure against CSR because they can suggest a way of fitting alternative models in case of rejection of the hypothesis of CSR. We also suggest the use of the K function for estimating parameters of proposed models because it provides an effective summary of spatial dependence over a wide range of scales [20,21].
4.3. Formal method for testing against complete spatial randomness To test departure against CSR, based on the F function as defined in Section 4.2 (it could be any empirical distribution function), ideally one should compute the difference Fˆ(x)E(Fˆ(x)) and compare that with the S.D. of Fˆ(x): Because of the difficulties in calculating E(Fˆ(x)) and var(Fˆ(x)); Diggle [20] describes a Monte Carlo based approach to provide the test. This approach examines the dof between the observed empirical distribution function and the expected empirical distribution function under CSR. It is given by x0
FIi
g fFˆ (x)[1exp(lpx )]g dx: 2
i
2
(6)
0
The statistic FIi does not have known distributions. The following Monte Carlo based method can be used to perform the test against the hypothesis of CSR and estimate the critical values. Let Fˆ1 (x) be the F -function of an observed point pattern with n events as defined in Eq. (5) and Fˆ2 (x); :::; Fˆs (x) the F -functions from s simulations of CSR patterns with n events. Calculate the statistic FIi for the observed and simulated patterns. Then, the value FIi for the observed pattern is compared with values FI2 ; . . . ; FIs for the simulated patterns. If FIi ranks among the largest of FI2 ; . . . ; FIs ; it indicates departure from CSR. Suppose F1 /F(j ) for some j f1; . . . ; sg then reject the hypotheses of CSR if Pvalue ((s1j)=2)5a:/ To carry out these tests, we used 99 simulations (s/ 100) from CSR to estimate the values of the F -function. The integral was, approximately, calculated by a Rie-
Table 1 P -values of the tests against the hypothesis of CSR based on the integrated squared difference of the F - and G -functions Test
F -function G -function
Sample S1
S2
S3
0.82 0.01
0.18 0.04
0.25 0.01
P -values were attached based on 99 Monte Carlo simulations (s/ 100) of a homogeneous Poisson process.
mann sum at 100 intervals between 0 and 0.05. The results of the tests against CSR are presented in Table 1. The results presented in Table 1 show that the test based on the G -function leads to an emphatic rejection of the hypotheses of CSR in all samples (P B/0.05), while the test based on the F -function suggests the acceptance of CSR in all samples. One can note that both functions provided contradictory answers about the nature of the pattern. We can also observe that the test based on the G -function does not give us information about the cause of the rejection of CSR and therefore one cannot decide whether the rejection is toward regular or cluster. The cause of the rejection of CSR can be explained by the graphical approach presented below. 4.4. Graphical method for testing against complete spatial randomness Let Fˆ1 (x) be the F-function of an observed point pattern with n events. Now, generate s-1 sets of n events by simulation of the random mechanism assumed to operate under null hypothesis of CSR. For these s -1 sets of simulated data, the upper and lower simulation envelopes, respectively, are: U(x)max Fˆi (x) i 2; . . . ; s; L(x) min Fˆi (x) i 2; . . . ; s:
(7) (8)
Now, plot Fˆ1 (x); U (x ), L (x ) as the ordinate (Fhat) against distances x . A plot of Fˆ1 (x) inside the simulation envelopes would not suggest rejection of the hypotheses of CSR, that is, events do not interact with each other and the expected number of events per area (intensity) is constant. The nature of the observed patterns depends on the type of function. If Fˆ1 (x) lies below the lower simulation envelope, it indicates deficiency of small nearest neighbour distance between point and event and that is typical of aggregated patterns. Otherwise, if Fˆ1 (x) lies above the upper simulation envelope, this typifies a regular pattern and indicates excess of small nearest neighbour distance between point and event. If Gˆ 1 (y) lies below the lower simulation envelope, it indicates deficiency of small nearest neighbour distance between events and that is typical of regular patterns.
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
251
Otherwise, if Gˆ 1 (y) lies above the upper simulation envelope, this typifies aggregated patterns and indicates excess of small nearest neighbour distance between events. We can see in Figs. 5 /7 that the empirical F -functions lie between the simulation envelopes and therefore the hypotheses of CSR are accepted in all samples. The most interesting feature of the G -function plots is the complete absence of small nearest neighbour distances between particles that are typical of regular patterns and therefore the hypothesis of CSR is rejected in all samples. This regular mechanism can correspond to the obvious fact that particles very close together can only coexist if their diameters are small. This premise is absolutely true in the material context, since there is non-overlapping particles in the planar sections of composite materials.
Fig. 6. Empirical F - and G -functions (dotted lines) with lower and upper envelopes from 99 simulations of a Poisson point process with intensity 507 (solid lines) in sample S2.
Fig. 5. Empirical F - and G -functions (dotted lines) with lower and upper envelopes from 99 simulations of a Poisson point process with intensity 419 (solid lines) in sample S1.
The results of the present section show the importance in using both the F - and G -functions to provide tests against CSR. It is clear that each of these two functions has a better power in detecting certain features of the pattern and therefore, the use of both functions makes a more sensitive analysis than would be obtained using either function alone. For example, it is well known [20 / 22] that analysis based on the G -function rather than analysis based on the F -function is thought to be more sensitive to regular departure from CSR, while the F function is more sensitive for detecting cluster. This behaviour explains the differences observed between the results provided by the F - and G -functions.
252
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
ical properties (e.g. toughness) of the composite material under investigation. Secondly, it could be used as a basis for comparing the variability of the mechanical properties in a sequence of planar sections by studying the variation of parameters of the spatial models of either the same composite material or different composite materials. Finally, on the basis of given parameters estimated from observed data, the model can be used for the prediction of spatial distribution of particle in a particular composite material. 5.1. The proposed models
Fig. 7. Empirical F - and G -functions (dotted lines) with lower and upper envelopes from 99 simulations of a Poisson point process with intensity 597 (solid lines) in sample S3.
5. Choosing and fitting stochastic models to spatial distribution of particles The analysis presented above was valuable in suggesting the kind of spatial structure that may be present in the data. Actually, the previous analysis indicated broad regularity at small scale. If we wish to explain the particular nature of regularity, then we need formal spatial point process models to allow us to do this. Modelling spatial distribution of reinforcing particles held in matrix of composite materials may be of interest for materials research from, at least, three viewpoints. Firstly, it could be fundamental in simulation work in order to relate spatial particle distribution and mechan-
In the present paper, we consider just the most simple models (simple sequential inhibition and Strauss) that can be appropriate to model particles regularly distributed over planar sections of composite material. These models are especially appropriate for modelling patterns where there are forces encouraging the events to be distributed more regularly. The simplest class of spatial point process models whose realisation exhibit regularity is the class of hardcore models. Such models arise most naturally by the imposition of a minimum permissible inhibition distance (8 ) between any two events. This may or may not simply reflect the physical size of the particles. Diggle [20] considers a simple sequential inhibition (SSI) process under which events are generated sequentially over a finite region A . At each stage, the event is added to the pattern only if it is at a distance of at least 8 from every other event already in the region A , otherwise it is ignored. The procedure terminates when a prespecified number of events have been placed in the region A , or it is impossible to continue. Using the following packing version can generate the Strauss processes, introduced by Kelly and Ripley [24]. The first event is placed uniformly in A. Subsequent events are generated uniformly in A , and accepted with probability gp , where p is the number of existing events closer than D to the possible new event. This process is repeated until a fixed number of events are tried or a given number is accepted. 5.2. Parameter estimation Since we have chosen some model to spatial distribution of second-phase inclusions, we need methods for estimating the parameters for them. Diggle [20] suggests a method based on the K- function to obtain leastsquares estimates of parameters in spatial point process models. He argues that an analysis using the K -function is attractive because its mathematical form is known, either explicitly or as an integral, for a number of potentially useful classes of models. Even when the mathematical form of the K -function is unknown, it can be easily estimated by using Monte Carlo simulation.
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
ˆ be the empirical distribution More formally, let K(t) function obtained from an observed point pattern and K(t; u) the corresponding theoretical K -function of a point process model with vector parameter u . Then, a (modified) least-squares estimator for u is obtained by minimising H(u ) in Eq. (9) by using the criterion based on the integrated squared distance between the values over the range of values 0/t0. That is,
253
samples. For example, the average particle size (graphical scale) is 0.0077 in sample S1. This means that the SSI model may allow overlapping particles and therefore, it is not a reasonable model for the spatial distribution of particles in our samples of composite material. We can conclude that the inhibition mechanism observed in our samples of composite material reflects the fact that the particles cannot overlap.
t0
g
c
c
ˆ H(u) [fK(t)g fK(t; u)g ]2 dt;
(9) 5.4. Fitting the Strauss model to the data
0
where c and t0 are tuning constants. The choice of the tuning constants is not thought to be critical. If the pattern has unit square, the limit of integration (0 B/t0 B/0.25) and the power transformation (c/0.5) are reasonable choices for regular patterns. For aggregated patterns, the choice would be (0 B/t0 B/0.25) and (c/0.25). In practice, one can try different combinations of tuning constants in order to provide desirable estimation properties [20]. For some point process models (e.g. Strauss), the value of K (t;u ) is known algebraically. For other models (e.g. SSI) where K (t;u) is unknown and cannot be evaluated either explicitly or numerically, Diggle [20] suggested that K (t;u ) might be estimated and replaced by s X
˜ K(t; u)
We have also applied the estimation method described above to estimate the inhibition distance (D ) and the interaction between neighbours value (g ) of the Strauss model in the three samples. The values of the tuning constants in Eq. (9) were t0 /0.05 and c /0.25. We have again used the estimator for K -function given by Eq. (3). The approximate theoretical K -function for the Strauss process is known [25] as: K(t; u)$
s
;
(10)
where Kˆ j (t) is calculated from s simulated realisations of the model. 5.3. Fitting the simple sequential inhibition model to the data We have applied the estimation method described above to estimate the inhibition distance (8 ) of the SSI in our samples. We have also used the form of the estimator for the K -function given in Eq. (3). The values of the tuning constants in Eq. (9) were d0 /0.05 and c/ 0.5. Since the K(d; u) is unknown for the SSI model, we replaced it in Eq. (9) by Eq. (10) calculated from 50 simulated realisations of the model evaluated. The estimates of 8 (in the transformed distances of the unit square) of the SSI were: 0.007421 (sample S1), 0.005729 (sample S2) and 0.006158 (sample S3). These values indicate the minimum distance allowed between one particular particle and other particles in a unit square planar section of composite materials. The inhibition distance infers that no two particles can be found closer than the estimate of 8 . We have observed that the estimates of 8 were smaller than the average particle size (radius of the equivalent disc) for all
0 Bt5D : t D
(11)
We have substituted K(t; u)/in Eq. (9) by Eq. (11). The estimates of the inhibition distance (D ) and the interaction between neighbours’ value (g) of the Strauss model in the three samples are presented in Table 2. The values of gˆ in Table 2 indicate the fraction of particles allowed to remain within the radius of inhibiˆ in the unit square planar sections of the tion (/D) composite material. For example, in sample S1, the interaction value implies that only 0.047% of the particles can be found closer than 0.0084. The same conclusion can be drawn from samples S2 and S3. We must observe that the values of Dˆ for all samples were greater than the observed minimum distance between two particles. For example, in sample S1, we found an inhibition distance equal to 0.008397, while the observed minimum distance between two particles in this sample is 0.0073. These results indicate that the Strauss model can provide a more reasonable model for the spatial distribution of particles since it may not allow overlapping particles.
Kˆ j (t)
j1
gpt2 pt2 (1g)pD2
Table 2 ˆ and the interaction value (/g) Estimates of the inhibition distance (/D) ˆ of the Strauss model in the three samples of the particulate composite material Parameter
Dˆ gˆ
/ /
Sample S1
S2
S3
0.008397 0.000473
0.006578 0.001152
0.007541 0.000338
254
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
statistic uses an integrated square distance between the distance values over the range of values 0 /x0. It is given by x0
fIi
g fFˆ (x) F¯ (x)g dx; 2
i
i
(12)
0
where Fˆ1 (x) is the F -function of an observed point pattern with n events as defined in Section 4.2 (it could be any pattern descriptor), Fˆ2 (x); . . . ; Fˆs (x) are the F functions from s simulations of the fitted model with n events and F¯i (x) is given by
Fig. 8. Goodness-of-fit plots based on the F - and G -functions (sample S1). The dotted lines are the observed data. The solid lines are the upper and lower envelopes from 99 simulations of fitted simple sequential inhibition process with 8 /0.007421.
6. Assessment of the fitted models Since we have the estimates of the parameters for both models, we have to check whether the models fit well to those samples. Thus, the present section is reserved for checking the goodness-of-fit of the models. There are two approaches to assessing the goodness-of-fit of the models: hypothesis test and graphical approaches. These methods are similar to testing departure against CSR. 6.1. Goodness-of-fit test The goodness-of-fit test, described in Diggle [20], examines the degree of agreement between the observed empirical distribution function and the expected empirical distribution function under the fitted model. The
Fig. 9. Goodness-of-fit plots based on the F - and G -functions (sample S2). The dotted lines are the observed data. The solid lines are the upper and lower envelopes from 99 simulations of fitted simple sequential inhibition process with 8 /0.005729.
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257 s X
Fˆj (x)
j2 F¯i (x) s1
:
(13)
The following Monte Carlo approach is applied to assess the goodness-of-fit. The value fI1 for the observed pattern is compared with values fI2,. . ., fIs from s simulations of the fitted model. Then rank fI1 from 1 to s and let j be the rank of fI1. If fI1 ranks among the largest of fI2,. . ., fIs , it indicates departure from CSR.
255
The significance level of the goodness-of-fit test is estimated to be Pvalue /(s/1/j )/s. If P-value is small then a poor fit of the observed pattern is inferred. We used 99 simulations (s /100) from the fitted models and 50 intervals from 0 to 0.05 to carry out these goodness-of-fit tests for both SSI and Strauss models. In sample S1, the P -values were: SSI (Ffunction /0.49, G- function /0.73) and Strauss (Ffunction /0.72, G- function/0.73). In sample S2, the P -values were: SSI (F-function /0.73, G- function / 0.53) and Strauss (F- function/0.89, G function / 0.37). In sample S3, the P -values were: SSI (Ffunction /0.71, G- function /0.05) and Strauss (Ffunction /0.89, G- function /0.01). These results, based on F -function, suggest that both the SSI and Strauss models can be reasonable choices to explain the spatial distribution of reinforcing particles in all samples The G -function rejects both models in sample S3. For a conclusive decision about the best model for the spatial distribution of particles in our samples, we provide further analyses in Section 6.2. 6.2. The graphical approach for testing of goodness-of-fit Let Fˆ1 (x) be the F-function of an observed point pattern with n events (it could be any function). Now, generate s -1 sets of n events by simulation of the fitted model. We define the upper and lower simulation envelopes for these s-1 sets of simulated data as being: U(x)maxFˆi (x) i 2; . . . ; s; L(x) minFˆi (x) i 2; . . . ; s
Fig. 10. Goodness-of-fit based on the F - and G -function for sample S3. The dotted line is the observed data. The solid lines are the upper and lower envelopes from 99 simulations of fitted simple sequential inhibition process with 8 /0.006158.
(14) (15)
We can plot Fˆ1 (x); U (x ), L (x ) as the ordinate (Fhat) against distances d as the abscissa. If the plot of Fˆ1 (x) is located inside the lower and upper simulation envelopes, it will not suggest rejection of the fitted model [20]. Figs. 8/10 show the estimated values of the F - and G functions for the observed data (calculated on 50 equally spaced distances between 0 and 0.06) and the upper and lower envelopes from 99 simulations (s/100) of the fitted SSI model as the ordinate against distances as the abscissa. We can see that the observed F -function lies between the simulation envelopes in all samples. On the other hand, the empirical distribution of G -function lies below the lower simulation envelopes in all samples. This behaviour shows that the SSI model did not fit adequately to explain the spatial variation of the particle centres in our samples of the composite material. The graphical method was also applied to test goodness-of-fit of the Strauss model. Figs. 11 /13 show these plots calculated over the same conditions as used before in the plots for the SSI model. We can see that the Strauss model improved the quality of the explanation of the spatial distribution of particles when we compare with the plots from the SSI model. The behaviour of
256
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
We can observe that the estimates of the radius of inhibition parameter were greater than the observed minimum radius of the particles and therefore the Strauss model was capable of detecting the presence of finite size of particles in these three samples. One must observe that the fitted model alone will not allow us to access the physical mechanisms responsible for the spatial distribution of the particles in the matrix. They may be powerful tools in the detection of changes in an experimental piece of work. The purpose of such model fitting is to establish simple (and yet compatible with a real material) definition of the basic characteristics of the spatial distribution of particles, thus providing an objective criterion for comparison between samples or composites obtained from different production processes.
Fig. 11. Goodness-of-fit plots based on F - and G -functions for sample S1. The dotted lines are the observed data. The solid lines are the upper and lower envelopes from 99 simulations of fitted Strauss process with D/0.008397 and g /0.000473.
these plots shows that the estimated values of radius and interaction value appear to fit moderately well to the samples and therefore, the Strauss model can be used to explain moderately well the spatial distribution of particles in our samples. In order to have a simple idea of the fitted Strauss model, we may think of each particle as a centre of a disc with radius D /2 and whose softness is controlled by parameter g. Let us assume that the particles are able to move freely, only subject to inter-particle interaction. If two particles get closer than D , the respective discs will tend to push them apart. Actually, this may be a particularly sensible way of viewing the model in composite materials, since there is evidence of particle migration over the matrix during the process of mixing and solidification of the composite [1].
Fig. 12. Goodness-of-fit plots based on F - and G -functions for sample S2. The dotted lines are the observed data. The solid lines are the upper and lower envelopes from 99 simulations of fitted Strauss process with D/0.006578 and g /0.001152.
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
257
function gives a very useful insight into the estimation of parameters of a stochastic spatial model. The results based on three samples of one aluminium silicon carbide composite material showed a statistically significant departure from CSR towards regularity. We also found that a Strauss process can conveniently model the spatial distribution of the second-phase particles.
Acknowledgements The careful and detailed comments of the referee were very much appreciated. J.D. Scalon was supported by a scholarship from the Ministry of Education (CAPES), Brazil.
References
Fig. 13. Goodness-of-fit based on the F - and G -functions for sample S3. The dotted line is the observed data. The solid lines are the upper and lower envelopes from 99 simulations of fitted Strauss process with D/0.007548 and g /0.000338.
7. Conclusions This paper suggests that one can provide a complete analysis of the spatial distribution of second phase particles held in planar sections of composite materials by carrying out five steps: visualising spatial patterns, testing against CSR, choosing appropriate stochastic models, fitting a model to the data and testing goodnessof-fit. The analysis clearly showed that functional pattern descriptors describe different properties of the spatial point patterns and all should take part in it. The F - and G -functions are powerful tools to summarise the essential properties of the pattern and therefore ideal for testing against both CSR and goodness-of-fit. The K -
[1] D.J. Lloyd, in: N. Hansen, D. Juuljensen, T. Leffers (Eds.), Metal Matrix Composites */Processing, Microstructure and Properties (Proceedings of the Twelfth Risø International Symposium on Materials Science), 1991, pp. 81 /99. [2] J.J. Lewandowski, C. Liu, W.H. Hunt, Mater. Sci. Eng. A107 (1989) (1989) 241. [3] M. Cans, H. Hougardy, Acta Stereol. 5 (1986) 281. [4] M.T. Shehata, J.D. Boyd, in: R. Rungta (Ed.), Inclusions and Their Influence on Material Behaviour, American Society for Metals, Metals Park, OH, 1988, p. 123. [5] T. Weiss, Pract. Metallog. 24 (1987) 507. [6] K. Wiencek, H. Hougardy, Acta Stereol. 16 (1987) 69. [7] I. Saxi, K. Pelikan, M. Besterci, Powder Metall. Int. 19 (1987) 27. [8] W.A. Spitzig, J.F. Kelly, O. Richmond, Metallography 18 (1985) 235. [9] P.J. Wray, O. Richmond, H. Morrison, Metallography 16 (1983) 39. [10] J.B. Parse, J.A. Wert, Model Simul. Mater. Sci. Eng. 1 (1993) 275. [11] R.K. Everett, J.H. Chu, J. Comp. Mater. 27 (1993) 1129. [12] M. Bertram, H. Wendrock, J. Microsc. 181 (1996) 45. [13] S. Ghosh, Z. Nowak, K. Lee, Acta Mater. 45 (1997) 2215. [14] J. Summerscales, D. Green, F.J. Guild, J. Microsc. 169 (1997) 173. [15] P. Louis, A.M. Gokhale, Metall. Mater. Trans. A Phys. Metall. 26 (1995) 1449. [16] J.P. Myles, E.C. Flenley, N.R.J. Fieller, H.V. Atkinson, H. Jones, Phil. Mag. 72 (1995) 515. [17] R. Pyrz, Compos. Sci. Tech. 50 (1995) 197. [18] J.E. Spowart, B. Maruyama, D.B. Miracle, Mater. Sci. Eng. A307 (2001) 51. [19] P.J. Diggle, Appl. Stats. 34 (1985) 138. [20] P.J. Diggle, Statistical Analysis of Spatial Point Patterns, Academic Press, London, 1983, p. 234. [21] B.D. Ripley, J. R. Statist. Soc. B 39 (1977) 212. [22] N. Cressie, Statistics for Spatial Data, Wiley and Sons, New York, 1993, p. 900. [23] A.J. Baddeley, B.W. Silverman, Biometrics 40 (1984) 1089. [24] F.P. Kelly, B.D. Ripley, Biometrika 63 (1976) 357. [25] V. Isham, Proc. R. Soc. Lond. A391 (1984) 39.
Spatial pattern analysis of second-phase particles in composite materials J.D. Scalon a,*, N.R.J. Fieller b, E.C. Stillman b, H.V. Atkinson c a
Department of Biomedical Engineering, Federal University of Sa˜o Joa˜o Del Rey, Prac¸a Dom Helve´cio 64, Sa˜o Joa˜o Del Rey MG 36301-160, Brazil b Department of Probability and Statistics, University of Sheffield, P.O. Box 597, Sheffield S10 2UN, UK c Department of Engineering Materials, University of Sheffield, P.O. Box 600, Sheffield S1 4DU, UK Received 30 July 2002; received in revised form 12 February 2003
Abstract The spatial distribution of second-phases is one of the most important factors affecting the mechanical properties of composite materials and therefore, the quantitative characterisation of such distributions is of prime importance in materials science. In this paper, we present a suitable procedure for analysing spatial distribution of particle centres held in three planar sections of an aluminium silicon carbide (Al/SiCp) composite material. We suggest that the spatial pattern analysis can be performed by following five consecutive steps: visualising spatial distribution, testing against complete spatial randomness, choosing and fitting appropriate stochastic models and testing goodness-of-fit. These steps show how statistical methods that rely on functional pattern descriptors and Monte Carlo simulations can be useful for characterising spatial distribution of the reinforcing particles in composite materials. # 2003 Elsevier Science B.V. All rights reserved. Keywords: Particle distribution; Reinforced aluminium; Spatial analysis; Simulation
1. Introduction It is well known that the mechanical properties (e.g. fracture toughness) of particulate composite materials depend not only on the shape and the volume fraction, but also on the spatial distribution of the reinforcing particles held in the matrix [1]. For example, there is evidence that surface cracking associated with clustered particles appear to be more prevalent than those associated with less clustered composite materials [1,2]. Thus, the quantitative characterisation of the spatial distribution of particles is of interest not only for a better understanding of the relationship between inclusions and mechanical behaviour, but also for better control of the production of composite materials. Consequently, there is an interest among material researchers in developing methods to provide a more rigorous quantitative analysis of the spatial distribution
* Corresponding author. Tel.: /55-32-3379-2505; fax: /55-323379-2328. E-mail address: [email protected] (J.D. Scalon).
of second-phase particles in two-dimensional distributed multi-phase systems. During the last few years, statistical methods have been appearing, to a certain extent, in materials literature to provide analysis of spatial distribution of particles in composite materials. Methods, such as quadrat count [3 /6], nearest neighbour distances [4,7,8], Dirichlet tessellation [2,4,8 /14], spatial pattern descriptors [15 /17] and correlation function [18] have been applied, basically, to test departure from randomness of the spatial distribution of the reinforcing particles. We have observed that some methods of spatial statistics were completely ignored by the researchers, for example, the fitting of a stochastic model for spatial distribution of particles and the use of techniques based on computer intensive methods for providing hypothesis tests. Thus, the aim of the present paper is to further explore the use of methods of spatial statistics (empirical distribution functions of pattern descriptors, such as F -, G - and K -functions and Monte Carlo simulation) for characterising the spatial distribution of reinforcing particles held in planar sections of composite materials.
0921-5093/03/$ - see front matter # 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0921-5093(03)00138-2
246
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
We suggest that one can provide a suitable quantitative analysis for characterising the spatial distribution of the particles by following five consecutive steps: visualising spatial distribution, testing against randomness, choosing an appropriate stochastic model, fitting the model to the pattern and testing goodness-of-fit. The main aim of this approach is to provide assessments of deviation from randomness towards either regularity or clustering of the particles and to suggest some alternative models for the spatial distribution of the reinforcing particles.
2. Description of experimental data We have applied the statistical analysis to three metallographic samples (S1 /S3) of an aluminium silicon carbide composite material produced by the Department of Engineering Materials, University of Sheffield, UK, where the silicon carbide second-phase particles had a volume fraction equal to 11%. The metallographic sample areas were :/192 /288 mm. The three metallographic samples were placed on a computer controlled optical microscope stage analyser (Polyvar), which allowed fully automatic adjustment, focusing, positioning and scanning of the samples. The overall microscope magnification used was 600 times, yielding a pixel size of 0.375 mm. Thus, the Polyvar produced three digital images with area frames equal to 512 /767 square pixels. The two-dimensional digital images were analysed by using image-processing techniques to extract the coordinates (centres) of each particle within the images. The actual images were 512/767 pixels, out of which we used only the particles located in the left top square of 512 /512 pixels, yielding 419 particles in sample S1, 507 particles in sample S2 and 597 particles in sample S3. The 512 /512 images were transformed into patterns with unit square area to facilitate the analysis. Then, we have provided the quantitative analysis of the spatial distributions of the reinforcing particles by following the suggested consecutive steps.
3. Visualising spatial distribution of the reinforcing particles The first requirement in any data analysis is the ability to see the data being analysed. Visualising spatial patterns means mapping the particles held in the matrix of the composite material. A natural tool for visualising a spatial distribution of particles is to plot their location in the matrix as a dot map. This map provides an idea of the distribution of particles and shows any obvious patterns presented in the sample. The dot maps of the samples presented in Fig. 1(a), Fig. 2(a) and Fig. 3(a) give a visual impression of the distribution of particles
Fig. 1. Dot map of locations (a) and contour plot (b) of the nonparametric kernel estimates of intensity using quartic kernel estimator and t /0.14 of 419 particles in sample S1. Contour values shown are the quantitative measures of local intensity of particles.
along the matrix. One must observe that it is very difficult to judge the nature of the patterns based only on these figures, that is, whether there is a tendency for aggregation or regularity. A more advanced statistical tool is the use of a spatially smooth estimate of the way in which the intensity of particles (the expected number of particles per unit area) is varying in the matrix. Diggle [19] describes a sophisticated approach for obtaining a spatially smooth intensity of particles based on a kernel estimator. More formally, let (s1, s2, . . . , sn ) be the spatial locations of n particles in a bounded study region jAj. The intensity, l (s), at s is estimated by lˆt (s)
n X 1 (s si ) ; k 2 t i1 t
(1)
where t /0 is the bandwidth parameter that determines the amount of smoothing. Following Diggle’s sugges-
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
Fig. 2. Dot map of locations (a) and contour plot (b) of the nonparametric kernel estimates of intensity using quartic kernel estimator and t /0.14 of 507 particles in sample S2. Contour values shown are the quantitative measures of local intensity of particles.
tion [19], we have adopted t /0.68n0.2. The k (.) is a bivariate probability density function, known as the kernel, which is symmetric about the origin. A simple quartic function can be used as an alternative for k (.) and then the estimate of intensity lˆt (s) may be simply expressed as n X 3 h2i 2 ˆ lt (s) 1 ; (2) 2 t2 i1 pt where hi is the distance between the point s and the observed particle location si . The smoothed intensity estimates lˆt (s) can be plotted in a variety of ways, e.g. perspective and contour. These
247
Fig. 3. Dot map of locations (a) and contour plot (b) of the nonparametric kernel estimates of intensity using quartic kernel estimator and t /0.13 of 597 particles in sample S3. Contour values shown are the quantitative measures of local intensity of particles.
plots can help identify peaks in the resulting surface corresponding to possible locations of clusters. Contour plots of our samples are presented in Fig. 1(b), Fig. 2(b) and Fig. 3(b). We can observe that these plots allow identify locally high concentrations of particles (the one with the local intensity is much bigger than the number of particles in the sample) and reduced local concentrations (where the local intensity is much lower than the number of particles in the sample) in the samples. These figures show evidence of the presence of both high and low local concentrations of particles in both samples and therefore, there is evidence that the distribution of particles along the matrix is not random.
248
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
One can observe that the kernel estimate is essentially a sum of ‘bumps’ placed at observations. The kernel function determines the shape of the bumps, while the bandwidth parameter determines their width. When estimating intensity, the most important thing is the choice of an appropriate bandwidth value, while the choice of kernel function is of secondary importance since any reasonable kernel (e.g. Gaussian and quartic) will give close to optimal results [19]. If the bandwidth tends to zero, the intensity is estimated as a series of ‘spikes’ in the individual observations, while if the bandwidth becomes large, all detail is obscured. In practice, we can try different values of bandwidth in order to explore the variation of intensity at different scales. Two examples are presented (Fig. 4) where values of t /0.06 and t/0.16 reveal that the higher values ‘smooth’ the distribution rather too much, while the lower value of 0.06 gives a rather too ‘spiky’ impression
Fig. 4. Contour plots of the nonparametric kernel estimates of intensity, in sample S1, using quartic kernel estimator for different values of the bandwidth parameter (t ). Contour values shown are the quantitative measures of local intensity of particles.
of the spatial distribution of the particles in sample S1. Diggle [19] suggests using t /0.68n0.2 for estimating the intensity on a unit square or other more efficient methods that attempt to optimise the value of bandwidth parameter based on estimates of the mean square error. We observe that the maps presented in Figs. 1/4 may be sufficient to answer some questions about the spatial distribution of second-phase particles in samples of composite material. However, they did not allow a conclusive answer about the nature (regular or clustered) of the pattern and therefore, it is appropriate to test whether the particle distribution can be considered random.
4. Testing against complete spatial randomness A large number of statistical methods have been proposed for testing departure against the null hypothesis of completely spatial randomness (CSR) in point patterns. The first group of techniques (quadrat count methods) is based on the determination of the number of particles in small test regions in order to assess spatial particle distribution. The observed frequency distribution of the number of particles per quadrat is compared to a theoretical Poisson distribution. If the particles are randomly distributed, the observed distribution will coincide closely with the theoretical distribution. In general, the x2-test for goodness-of-fit is applied to assess the degree of agreement [3 /6]. One must observe that the choice of the region size is a crucial problem in quadrat methods and therefore they depend on the arbitrary choice of both the quadrat area size and quadrat shape. The second group of techniques is based on two approaches of nearest neighbour distances that require information about the locations of a sample of particle centres held in the composite material. The first approach is based on the distance between a randomly chosen particle and its nearest neighbour. The second approach consists in choosing a point at random and then the nearest distance between the point and particle in the pattern is determined. The aim, in both approaches, is to compare the observed distribution of the nearest neighbour distances and the theoretical distribution of distances under randomness (Poisson distribution). If the particles are randomly distributed, the observed distribution will coincide closely with the theoretical distribution [4,7,8]. The drawback of these techniques is that the reduction of complex point patterns to a single nearest neighbour summary statistic, calculated from locations of a sample of particles, results in a considerable loss of information. Thus, this sort of methods is not recommended when the information
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
about the location of all second-phase particle centres held in the matrix of composite material is available. The third and most widely used of all statistical methods for the analysis of spatial distribution of particles in composite materials is based on Dirichlet (Delaunay or Voronoi) tessellation. This method uses a geometrical construction for partitioning a field into a unique set of convex polyhedra (polygons), each of which is associated with and contains one of the particles. The polygons form a continuous cellular structure. Analysis of the characteristics of these polyhedra (e.g. area) provides information about the local environment of the particles [2,4,8 /12]. It is observed that Dirichlet-based methods are not effective in distinguishing cluster, regular and random patterns for high volume (area) fractions of particles [13]. Otherwise, Dirichlet-based methods may neglect some significant areas of space containing no events [14]. We note that the empirical distribution function for some variables related to Dirichlet tessellation (e.g. area of the polygon) are not known for either Poisson processes or a number of potentially useful classes of spatial models. This fact can limit the use of Dirichlet tessellation for either conducting formal significant tests against randomness or estimating parameters of spatial models. One of the least frequently used group of all statistical methods for the analysis of spatial distributions of particles in composite materials is based on spatial pattern descriptors, such as F -, G - and K -functions. These methods require a completely enumerated point pattern, that is, one needs information about the location of all particles held in the composite material. These functions have been used on a few works for the analysis of spatial distribution of inclusions in composite material [15 /17], but the full potential has not been totally explored. Thus, in this paper we are using tests based on these pattern descriptors. For practical purposes, each second-phase particle is treated as a point defined by its co-ordinates. Thus, our spatial data can be assumed to be a map of all particle locations in an essential planar region and therefore the analysis can be founded on the theory of spatial point patterns. Typically, an observed point pattern can be thought of as the outcome (a realisation) of a spatial point process. Because spatial point pattern data are usually a single realisation of a spatial point process, the additional assumptions of stationarity and sometimes isotropy, are often made to reduce the parameter space and to allow parameter estimation [20]. One can never hope to totally characterise the spatial point process, but one can investigate the behaviour of functional pattern descriptors, proposed by Ripley [21] and described in detail by Diggle [20] and Cressie [22], that represent important aspects of the process.
249
4.1. The K -function A convenient way of characterising spatial point patterns is provided by the K -function. The theoretical K -function of a stationary spatial point process is defined by K (t)/l 1E (number of further events within distance t of an arbitrary event). Under the hypothesis of CSR, the K -function is given by K (t)/ pt2. An unbiased estimator for K (t) that corrects edge effects is given by ˆ K(t)
n X n ½A½ X It (uij )
n2
i1 j1
wij
;
(3)
where uij is the distance between events i and j (i "/j), It (uij ) is an indicator function that takes the value 1 when uij is less than distance t, l is the intensity, n is the number of events in the analysed pattern of area jA j and wij is a weighting factor that represents the proportion of the circumference of the circle around event i , passing through the point j that lies within jAj. The K -function represent a natural and valuable tool for the description of a spatial point process, but it does not provide a complete description of the probability structure of a spatial point process. It is possible to find different processes that have identical K -functions [23]. For a complete description of a spatial point pattern it would be necessary to include other functional pattern descriptors, such as those based on nearest neighbour distances. 4.2. The F - and G -functions The F- function is based on the distance x , that is, the distance from an arbitrary chosen point (not event) to the nearest event. It is defined by the distribution function F (x )/P {distance from an arbitrary point to the nearest event is at most x }. The theoretical F function under the null hypothesis of CSR is given by F (x)1exp(lpx2 ) x 0:
(4)
The simplest suitable edge-corrected estimator of the F -function is provided by the equation n X
Ix (xi ; ri ) i1 ˆ F (x) n X Ix (ri )
x 0;
(5)
i1
where m is the number of sample points in the analysed pattern, xi denotes the distance from the ith chosen point to the nearest of the n events in the analysed pattern, Ix (ri ) is an indicator function that takes the value 1 when xi is /x , ri is the distance from each event to the nearest point on the boundary of jAj and Ix (xi ,ri )
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
250
is an indicator function that takes the value 1 when xi is /x and ri is ‘/x. We can use m sample points in a regular grid k /k to estimate the values of the F pffiffiffi function, where k $ n:/ The G -function is based on the distance y, that is, the distance between two nearest events. It is given by the expression G (y) /P {distance from an arbitrary event to the nearest other event is at most y }. The theoretical G function under the null hypothesis of CSR and an edgecorrected estimator for G (y) is provided by equations analogous to Eq. (4) and Eq. (5), respectively, substituting the distance x by the distance y. We suggest the use of F - and G -functions to test departure against CSR because they can suggest a way of fitting alternative models in case of rejection of the hypothesis of CSR. We also suggest the use of the K function for estimating parameters of proposed models because it provides an effective summary of spatial dependence over a wide range of scales [20,21].
4.3. Formal method for testing against complete spatial randomness To test departure against CSR, based on the F function as defined in Section 4.2 (it could be any empirical distribution function), ideally one should compute the difference Fˆ(x)E(Fˆ(x)) and compare that with the S.D. of Fˆ(x): Because of the difficulties in calculating E(Fˆ(x)) and var(Fˆ(x)); Diggle [20] describes a Monte Carlo based approach to provide the test. This approach examines the dof between the observed empirical distribution function and the expected empirical distribution function under CSR. It is given by x0
FIi
g fFˆ (x)[1exp(lpx )]g dx: 2
i
2
(6)
0
The statistic FIi does not have known distributions. The following Monte Carlo based method can be used to perform the test against the hypothesis of CSR and estimate the critical values. Let Fˆ1 (x) be the F -function of an observed point pattern with n events as defined in Eq. (5) and Fˆ2 (x); :::; Fˆs (x) the F -functions from s simulations of CSR patterns with n events. Calculate the statistic FIi for the observed and simulated patterns. Then, the value FIi for the observed pattern is compared with values FI2 ; . . . ; FIs for the simulated patterns. If FIi ranks among the largest of FI2 ; . . . ; FIs ; it indicates departure from CSR. Suppose F1 /F(j ) for some j f1; . . . ; sg then reject the hypotheses of CSR if Pvalue ((s1j)=2)5a:/ To carry out these tests, we used 99 simulations (s/ 100) from CSR to estimate the values of the F -function. The integral was, approximately, calculated by a Rie-
Table 1 P -values of the tests against the hypothesis of CSR based on the integrated squared difference of the F - and G -functions Test
F -function G -function
Sample S1
S2
S3
0.82 0.01
0.18 0.04
0.25 0.01
P -values were attached based on 99 Monte Carlo simulations (s/ 100) of a homogeneous Poisson process.
mann sum at 100 intervals between 0 and 0.05. The results of the tests against CSR are presented in Table 1. The results presented in Table 1 show that the test based on the G -function leads to an emphatic rejection of the hypotheses of CSR in all samples (P B/0.05), while the test based on the F -function suggests the acceptance of CSR in all samples. One can note that both functions provided contradictory answers about the nature of the pattern. We can also observe that the test based on the G -function does not give us information about the cause of the rejection of CSR and therefore one cannot decide whether the rejection is toward regular or cluster. The cause of the rejection of CSR can be explained by the graphical approach presented below. 4.4. Graphical method for testing against complete spatial randomness Let Fˆ1 (x) be the F-function of an observed point pattern with n events. Now, generate s-1 sets of n events by simulation of the random mechanism assumed to operate under null hypothesis of CSR. For these s -1 sets of simulated data, the upper and lower simulation envelopes, respectively, are: U(x)max Fˆi (x) i 2; . . . ; s; L(x) min Fˆi (x) i 2; . . . ; s:
(7) (8)
Now, plot Fˆ1 (x); U (x ), L (x ) as the ordinate (Fhat) against distances x . A plot of Fˆ1 (x) inside the simulation envelopes would not suggest rejection of the hypotheses of CSR, that is, events do not interact with each other and the expected number of events per area (intensity) is constant. The nature of the observed patterns depends on the type of function. If Fˆ1 (x) lies below the lower simulation envelope, it indicates deficiency of small nearest neighbour distance between point and event and that is typical of aggregated patterns. Otherwise, if Fˆ1 (x) lies above the upper simulation envelope, this typifies a regular pattern and indicates excess of small nearest neighbour distance between point and event. If Gˆ 1 (y) lies below the lower simulation envelope, it indicates deficiency of small nearest neighbour distance between events and that is typical of regular patterns.
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
251
Otherwise, if Gˆ 1 (y) lies above the upper simulation envelope, this typifies aggregated patterns and indicates excess of small nearest neighbour distance between events. We can see in Figs. 5 /7 that the empirical F -functions lie between the simulation envelopes and therefore the hypotheses of CSR are accepted in all samples. The most interesting feature of the G -function plots is the complete absence of small nearest neighbour distances between particles that are typical of regular patterns and therefore the hypothesis of CSR is rejected in all samples. This regular mechanism can correspond to the obvious fact that particles very close together can only coexist if their diameters are small. This premise is absolutely true in the material context, since there is non-overlapping particles in the planar sections of composite materials.
Fig. 6. Empirical F - and G -functions (dotted lines) with lower and upper envelopes from 99 simulations of a Poisson point process with intensity 507 (solid lines) in sample S2.
Fig. 5. Empirical F - and G -functions (dotted lines) with lower and upper envelopes from 99 simulations of a Poisson point process with intensity 419 (solid lines) in sample S1.
The results of the present section show the importance in using both the F - and G -functions to provide tests against CSR. It is clear that each of these two functions has a better power in detecting certain features of the pattern and therefore, the use of both functions makes a more sensitive analysis than would be obtained using either function alone. For example, it is well known [20 / 22] that analysis based on the G -function rather than analysis based on the F -function is thought to be more sensitive to regular departure from CSR, while the F function is more sensitive for detecting cluster. This behaviour explains the differences observed between the results provided by the F - and G -functions.
252
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
ical properties (e.g. toughness) of the composite material under investigation. Secondly, it could be used as a basis for comparing the variability of the mechanical properties in a sequence of planar sections by studying the variation of parameters of the spatial models of either the same composite material or different composite materials. Finally, on the basis of given parameters estimated from observed data, the model can be used for the prediction of spatial distribution of particle in a particular composite material. 5.1. The proposed models
Fig. 7. Empirical F - and G -functions (dotted lines) with lower and upper envelopes from 99 simulations of a Poisson point process with intensity 597 (solid lines) in sample S3.
5. Choosing and fitting stochastic models to spatial distribution of particles The analysis presented above was valuable in suggesting the kind of spatial structure that may be present in the data. Actually, the previous analysis indicated broad regularity at small scale. If we wish to explain the particular nature of regularity, then we need formal spatial point process models to allow us to do this. Modelling spatial distribution of reinforcing particles held in matrix of composite materials may be of interest for materials research from, at least, three viewpoints. Firstly, it could be fundamental in simulation work in order to relate spatial particle distribution and mechan-
In the present paper, we consider just the most simple models (simple sequential inhibition and Strauss) that can be appropriate to model particles regularly distributed over planar sections of composite material. These models are especially appropriate for modelling patterns where there are forces encouraging the events to be distributed more regularly. The simplest class of spatial point process models whose realisation exhibit regularity is the class of hardcore models. Such models arise most naturally by the imposition of a minimum permissible inhibition distance (8 ) between any two events. This may or may not simply reflect the physical size of the particles. Diggle [20] considers a simple sequential inhibition (SSI) process under which events are generated sequentially over a finite region A . At each stage, the event is added to the pattern only if it is at a distance of at least 8 from every other event already in the region A , otherwise it is ignored. The procedure terminates when a prespecified number of events have been placed in the region A , or it is impossible to continue. Using the following packing version can generate the Strauss processes, introduced by Kelly and Ripley [24]. The first event is placed uniformly in A. Subsequent events are generated uniformly in A , and accepted with probability gp , where p is the number of existing events closer than D to the possible new event. This process is repeated until a fixed number of events are tried or a given number is accepted. 5.2. Parameter estimation Since we have chosen some model to spatial distribution of second-phase inclusions, we need methods for estimating the parameters for them. Diggle [20] suggests a method based on the K- function to obtain leastsquares estimates of parameters in spatial point process models. He argues that an analysis using the K -function is attractive because its mathematical form is known, either explicitly or as an integral, for a number of potentially useful classes of models. Even when the mathematical form of the K -function is unknown, it can be easily estimated by using Monte Carlo simulation.
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
ˆ be the empirical distribution More formally, let K(t) function obtained from an observed point pattern and K(t; u) the corresponding theoretical K -function of a point process model with vector parameter u . Then, a (modified) least-squares estimator for u is obtained by minimising H(u ) in Eq. (9) by using the criterion based on the integrated squared distance between the values over the range of values 0/t0. That is,
253
samples. For example, the average particle size (graphical scale) is 0.0077 in sample S1. This means that the SSI model may allow overlapping particles and therefore, it is not a reasonable model for the spatial distribution of particles in our samples of composite material. We can conclude that the inhibition mechanism observed in our samples of composite material reflects the fact that the particles cannot overlap.
t0
g
c
c
ˆ H(u) [fK(t)g fK(t; u)g ]2 dt;
(9) 5.4. Fitting the Strauss model to the data
0
where c and t0 are tuning constants. The choice of the tuning constants is not thought to be critical. If the pattern has unit square, the limit of integration (0 B/t0 B/0.25) and the power transformation (c/0.5) are reasonable choices for regular patterns. For aggregated patterns, the choice would be (0 B/t0 B/0.25) and (c/0.25). In practice, one can try different combinations of tuning constants in order to provide desirable estimation properties [20]. For some point process models (e.g. Strauss), the value of K (t;u ) is known algebraically. For other models (e.g. SSI) where K (t;u) is unknown and cannot be evaluated either explicitly or numerically, Diggle [20] suggested that K (t;u ) might be estimated and replaced by s X
˜ K(t; u)
We have also applied the estimation method described above to estimate the inhibition distance (D ) and the interaction between neighbours value (g ) of the Strauss model in the three samples. The values of the tuning constants in Eq. (9) were t0 /0.05 and c /0.25. We have again used the estimator for K -function given by Eq. (3). The approximate theoretical K -function for the Strauss process is known [25] as: K(t; u)$
s
;
(10)
where Kˆ j (t) is calculated from s simulated realisations of the model. 5.3. Fitting the simple sequential inhibition model to the data We have applied the estimation method described above to estimate the inhibition distance (8 ) of the SSI in our samples. We have also used the form of the estimator for the K -function given in Eq. (3). The values of the tuning constants in Eq. (9) were d0 /0.05 and c/ 0.5. Since the K(d; u) is unknown for the SSI model, we replaced it in Eq. (9) by Eq. (10) calculated from 50 simulated realisations of the model evaluated. The estimates of 8 (in the transformed distances of the unit square) of the SSI were: 0.007421 (sample S1), 0.005729 (sample S2) and 0.006158 (sample S3). These values indicate the minimum distance allowed between one particular particle and other particles in a unit square planar section of composite materials. The inhibition distance infers that no two particles can be found closer than the estimate of 8 . We have observed that the estimates of 8 were smaller than the average particle size (radius of the equivalent disc) for all
0 Bt5D : t D
(11)
We have substituted K(t; u)/in Eq. (9) by Eq. (11). The estimates of the inhibition distance (D ) and the interaction between neighbours’ value (g) of the Strauss model in the three samples are presented in Table 2. The values of gˆ in Table 2 indicate the fraction of particles allowed to remain within the radius of inhibiˆ in the unit square planar sections of the tion (/D) composite material. For example, in sample S1, the interaction value implies that only 0.047% of the particles can be found closer than 0.0084. The same conclusion can be drawn from samples S2 and S3. We must observe that the values of Dˆ for all samples were greater than the observed minimum distance between two particles. For example, in sample S1, we found an inhibition distance equal to 0.008397, while the observed minimum distance between two particles in this sample is 0.0073. These results indicate that the Strauss model can provide a more reasonable model for the spatial distribution of particles since it may not allow overlapping particles.
Kˆ j (t)
j1
gpt2 pt2 (1g)pD2
Table 2 ˆ and the interaction value (/g) Estimates of the inhibition distance (/D) ˆ of the Strauss model in the three samples of the particulate composite material Parameter
Dˆ gˆ
/ /
Sample S1
S2
S3
0.008397 0.000473
0.006578 0.001152
0.007541 0.000338
254
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
statistic uses an integrated square distance between the distance values over the range of values 0 /x0. It is given by x0
fIi
g fFˆ (x) F¯ (x)g dx; 2
i
i
(12)
0
where Fˆ1 (x) is the F -function of an observed point pattern with n events as defined in Section 4.2 (it could be any pattern descriptor), Fˆ2 (x); . . . ; Fˆs (x) are the F functions from s simulations of the fitted model with n events and F¯i (x) is given by
Fig. 8. Goodness-of-fit plots based on the F - and G -functions (sample S1). The dotted lines are the observed data. The solid lines are the upper and lower envelopes from 99 simulations of fitted simple sequential inhibition process with 8 /0.007421.
6. Assessment of the fitted models Since we have the estimates of the parameters for both models, we have to check whether the models fit well to those samples. Thus, the present section is reserved for checking the goodness-of-fit of the models. There are two approaches to assessing the goodness-of-fit of the models: hypothesis test and graphical approaches. These methods are similar to testing departure against CSR. 6.1. Goodness-of-fit test The goodness-of-fit test, described in Diggle [20], examines the degree of agreement between the observed empirical distribution function and the expected empirical distribution function under the fitted model. The
Fig. 9. Goodness-of-fit plots based on the F - and G -functions (sample S2). The dotted lines are the observed data. The solid lines are the upper and lower envelopes from 99 simulations of fitted simple sequential inhibition process with 8 /0.005729.
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257 s X
Fˆj (x)
j2 F¯i (x) s1
:
(13)
The following Monte Carlo approach is applied to assess the goodness-of-fit. The value fI1 for the observed pattern is compared with values fI2,. . ., fIs from s simulations of the fitted model. Then rank fI1 from 1 to s and let j be the rank of fI1. If fI1 ranks among the largest of fI2,. . ., fIs , it indicates departure from CSR.
255
The significance level of the goodness-of-fit test is estimated to be Pvalue /(s/1/j )/s. If P-value is small then a poor fit of the observed pattern is inferred. We used 99 simulations (s /100) from the fitted models and 50 intervals from 0 to 0.05 to carry out these goodness-of-fit tests for both SSI and Strauss models. In sample S1, the P -values were: SSI (Ffunction /0.49, G- function /0.73) and Strauss (Ffunction /0.72, G- function/0.73). In sample S2, the P -values were: SSI (F-function /0.73, G- function / 0.53) and Strauss (F- function/0.89, G function / 0.37). In sample S3, the P -values were: SSI (Ffunction /0.71, G- function /0.05) and Strauss (Ffunction /0.89, G- function /0.01). These results, based on F -function, suggest that both the SSI and Strauss models can be reasonable choices to explain the spatial distribution of reinforcing particles in all samples The G -function rejects both models in sample S3. For a conclusive decision about the best model for the spatial distribution of particles in our samples, we provide further analyses in Section 6.2. 6.2. The graphical approach for testing of goodness-of-fit Let Fˆ1 (x) be the F-function of an observed point pattern with n events (it could be any function). Now, generate s -1 sets of n events by simulation of the fitted model. We define the upper and lower simulation envelopes for these s-1 sets of simulated data as being: U(x)maxFˆi (x) i 2; . . . ; s; L(x) minFˆi (x) i 2; . . . ; s
Fig. 10. Goodness-of-fit based on the F - and G -function for sample S3. The dotted line is the observed data. The solid lines are the upper and lower envelopes from 99 simulations of fitted simple sequential inhibition process with 8 /0.006158.
(14) (15)
We can plot Fˆ1 (x); U (x ), L (x ) as the ordinate (Fhat) against distances d as the abscissa. If the plot of Fˆ1 (x) is located inside the lower and upper simulation envelopes, it will not suggest rejection of the fitted model [20]. Figs. 8/10 show the estimated values of the F - and G functions for the observed data (calculated on 50 equally spaced distances between 0 and 0.06) and the upper and lower envelopes from 99 simulations (s/100) of the fitted SSI model as the ordinate against distances as the abscissa. We can see that the observed F -function lies between the simulation envelopes in all samples. On the other hand, the empirical distribution of G -function lies below the lower simulation envelopes in all samples. This behaviour shows that the SSI model did not fit adequately to explain the spatial variation of the particle centres in our samples of the composite material. The graphical method was also applied to test goodness-of-fit of the Strauss model. Figs. 11 /13 show these plots calculated over the same conditions as used before in the plots for the SSI model. We can see that the Strauss model improved the quality of the explanation of the spatial distribution of particles when we compare with the plots from the SSI model. The behaviour of
256
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
We can observe that the estimates of the radius of inhibition parameter were greater than the observed minimum radius of the particles and therefore the Strauss model was capable of detecting the presence of finite size of particles in these three samples. One must observe that the fitted model alone will not allow us to access the physical mechanisms responsible for the spatial distribution of the particles in the matrix. They may be powerful tools in the detection of changes in an experimental piece of work. The purpose of such model fitting is to establish simple (and yet compatible with a real material) definition of the basic characteristics of the spatial distribution of particles, thus providing an objective criterion for comparison between samples or composites obtained from different production processes.
Fig. 11. Goodness-of-fit plots based on F - and G -functions for sample S1. The dotted lines are the observed data. The solid lines are the upper and lower envelopes from 99 simulations of fitted Strauss process with D/0.008397 and g /0.000473.
these plots shows that the estimated values of radius and interaction value appear to fit moderately well to the samples and therefore, the Strauss model can be used to explain moderately well the spatial distribution of particles in our samples. In order to have a simple idea of the fitted Strauss model, we may think of each particle as a centre of a disc with radius D /2 and whose softness is controlled by parameter g. Let us assume that the particles are able to move freely, only subject to inter-particle interaction. If two particles get closer than D , the respective discs will tend to push them apart. Actually, this may be a particularly sensible way of viewing the model in composite materials, since there is evidence of particle migration over the matrix during the process of mixing and solidification of the composite [1].
Fig. 12. Goodness-of-fit plots based on F - and G -functions for sample S2. The dotted lines are the observed data. The solid lines are the upper and lower envelopes from 99 simulations of fitted Strauss process with D/0.006578 and g /0.001152.
J.D. Scalon et al. / Materials Science and Engineering A356 (2003) 245 /257
257
function gives a very useful insight into the estimation of parameters of a stochastic spatial model. The results based on three samples of one aluminium silicon carbide composite material showed a statistically significant departure from CSR towards regularity. We also found that a Strauss process can conveniently model the spatial distribution of the second-phase particles.
Acknowledgements The careful and detailed comments of the referee were very much appreciated. J.D. Scalon was supported by a scholarship from the Ministry of Education (CAPES), Brazil.
References
Fig. 13. Goodness-of-fit based on the F - and G -functions for sample S3. The dotted line is the observed data. The solid lines are the upper and lower envelopes from 99 simulations of fitted Strauss process with D/0.007548 and g /0.000338.
7. Conclusions This paper suggests that one can provide a complete analysis of the spatial distribution of second phase particles held in planar sections of composite materials by carrying out five steps: visualising spatial patterns, testing against CSR, choosing appropriate stochastic models, fitting a model to the data and testing goodnessof-fit. The analysis clearly showed that functional pattern descriptors describe different properties of the spatial point patterns and all should take part in it. The F - and G -functions are powerful tools to summarise the essential properties of the pattern and therefore ideal for testing against both CSR and goodness-of-fit. The K -
[1] D.J. Lloyd, in: N. Hansen, D. Juuljensen, T. Leffers (Eds.), Metal Matrix Composites */Processing, Microstructure and Properties (Proceedings of the Twelfth Risø International Symposium on Materials Science), 1991, pp. 81 /99. [2] J.J. Lewandowski, C. Liu, W.H. Hunt, Mater. Sci. Eng. A107 (1989) (1989) 241. [3] M. Cans, H. Hougardy, Acta Stereol. 5 (1986) 281. [4] M.T. Shehata, J.D. Boyd, in: R. Rungta (Ed.), Inclusions and Their Influence on Material Behaviour, American Society for Metals, Metals Park, OH, 1988, p. 123. [5] T. Weiss, Pract. Metallog. 24 (1987) 507. [6] K. Wiencek, H. Hougardy, Acta Stereol. 16 (1987) 69. [7] I. Saxi, K. Pelikan, M. Besterci, Powder Metall. Int. 19 (1987) 27. [8] W.A. Spitzig, J.F. Kelly, O. Richmond, Metallography 18 (1985) 235. [9] P.J. Wray, O. Richmond, H. Morrison, Metallography 16 (1983) 39. [10] J.B. Parse, J.A. Wert, Model Simul. Mater. Sci. Eng. 1 (1993) 275. [11] R.K. Everett, J.H. Chu, J. Comp. Mater. 27 (1993) 1129. [12] M. Bertram, H. Wendrock, J. Microsc. 181 (1996) 45. [13] S. Ghosh, Z. Nowak, K. Lee, Acta Mater. 45 (1997) 2215. [14] J. Summerscales, D. Green, F.J. Guild, J. Microsc. 169 (1997) 173. [15] P. Louis, A.M. Gokhale, Metall. Mater. Trans. A Phys. Metall. 26 (1995) 1449. [16] J.P. Myles, E.C. Flenley, N.R.J. Fieller, H.V. Atkinson, H. Jones, Phil. Mag. 72 (1995) 515. [17] R. Pyrz, Compos. Sci. Tech. 50 (1995) 197. [18] J.E. Spowart, B. Maruyama, D.B. Miracle, Mater. Sci. Eng. A307 (2001) 51. [19] P.J. Diggle, Appl. Stats. 34 (1985) 138. [20] P.J. Diggle, Statistical Analysis of Spatial Point Patterns, Academic Press, London, 1983, p. 234. [21] B.D. Ripley, J. R. Statist. Soc. B 39 (1977) 212. [22] N. Cressie, Statistics for Spatial Data, Wiley and Sons, New York, 1993, p. 900. [23] A.J. Baddeley, B.W. Silverman, Biometrics 40 (1984) 1089. [24] F.P. Kelly, B.D. Ripley, Biometrika 63 (1976) 357. [25] V. Isham, Proc. R. Soc. Lond. A391 (1984) 39.