Fibre composites, e.g., fibre reinforced polymers or high performance concrete, are an important class of functional materials. Physical properties of a fibre composite such as elasticity or crack propagation are influenced by its microstructure characteristics including the fibre volume fraction, the size or the direction distribution of the fibres. Therefore, an understanding of the relations between the fibre geometry and macroscopic properties is crucial for the optimisation of materials for certain applications. During the last years, micro computed tomography (μCT) has proven to be a powerful tool for the analysis of the three-dimensional microstructure of materials.
In the compression moulding process of glass fibre reinforced polymers, the fibres order themselves inside the raw material as a result of mechanical pressure. During this process, deviations from the requested direction may occur, creating undesirable fibre clusters and/or deformations. These inhomogeneities are characterized by abrupt changes in the direction of the fibres, and their detection is studied in this paper.
The problem of detecting change points in random sequences, (multivariate) time series, panel and regression data has a long history, see the books BrassNik93 ; BrodDarkh00 ; Carlstein_ed94 ; ChenGupta12 ; TimeS1 ; Wu05
. Changes to be detected may concern the mean, variance, correlation, spectral density, etc. of the (stationary) sequence. This kind of change detection has been considered by various authors starting with Page . Sen and Srivastava Sen considered tests for a change in the mean of a Gaussian model. An overview can also be found in Brodsky . The CUSUM procedure, Bayesian approaches as well as maximum likelihood estimation are often used. Scan statistics come also into play naturally, see e.g. BrodDarkh99 ; BrodDarkh00 .
First approaches to change point analysis for random fields (or measures) have been developed in the papers Buc14 ; BucHeus15 ; CaoWors99 ; Chamb02 ; HanubiaMnatsakanov96 ; PiterJar11 ; Kaplan90 ; Kaplan92 ; Lai08 ; MuellSong94 ; Ninomiya04 ; Sharietal16 ; SiegYakir08 ; SiegWors95 , see also the review in (BrodDarkh99, , Section 2, D) and (BrodDarkh00, , Chapter 6). The involved methods include M-estimation, minimax methods for risks, the geometric tube method, some nonparametric and Bayesian techniques. However, much is still to be done in this relatively new area of research.
In this paper, we develop a change-point test for dependent random fields. In the spirit of the book BrodDarkh00 , it uses inequalities for tail probabilities of suitable test statistics. It is applied to the mean and the entropy of the local directional distribution of fibres observed in a 3D image of a fibre composite obtained by micro computed tomography. Characteristics are estimated in a moving scanning window that runs over the observed material sample, cf. ruiz2016entropy ; Alonzo . Our main task is to detect areas with anomalous spreading of the fibres. Even though we focus on anomalies in fibres’ directions, our method will work with any local characteristic of fibres with values in a (compact) Riemannian manifold such as fibre length or mean curvature.
If an anomaly is present, its location is detected using a new spatial modification of the Stochastic Approximation Expectation Maximization (SAEM) algorithm (see EM
for a review of Expectation Maximization (EM) algorithms for the separation of components in a mixture of Gaussian distributions as well as a recent paperLaurent ). It allows for spatial clustering of the whole fibre material into a “normal” and an “anomaly” zones.
The paper is organized as follows. In Section 2, we introduce the stochastic model of a fibre process. In Section 3, we describe the procedure of generating the sample data, introduce the mean of local directions as well as their entropy. There, we compare two methods for entropy estimation: plug-in and nearest neighbor statistics. In Section 4, we consider the detection of anomalies as a change-point problem for the corresponding dependent random fields. In Section 5, we localize the anomalous region of fibres solving a clustering problem for multivariate random fields. For this purpose, we propose a new spatial modification of SAEM algorithm, which decreases the diffuseness of clusters. In Section 6, we apply our methods to 3D images of simulated (Section 6.1) and real (Section 6.2) fibre materials and compare their performance.
2 Problem setting
In this section, we give some basic definitions and results for fibre processes. For more details, see, for example, the book stoyan . In 3-dimensional Euclidean space, a fibre is a simple curve of finite length satisfying the following assumptions:
is a -smooth function.
for all where
A fibre does not intersect itself.
The collection of fibres forms a fibre system if it is a union of at most countably many fibres such that any compact set is intersected by only a finite number of fibres, and if i.e., the distinct fibres may have only end-points in common. The length measure corresponding to the fibre system (and denoted by the same symbol) is defined by
for bounded Borel sets where is the length of fibre in window Then is the total length of fibre pieces in the window
A fibre process is a random element with values in the set of all fibre systems with -algebra generated by sets of the form for all bounded Borel sets and real numbers . The distribution of a fibre process is a probability measure on The fibre process is said to be stationary if it has the same distribution as the translated fibre process for all
For classification needs we consider an abstract fibre characteristic . Let be a measurable space where is a (compact) Riemannian manifold equipped with a metric Let be some characteristic of a fibre at point assuming that exactly one fibre of passes through Then a weighted random measure can be defined by
for bounded and Thus, is the total length of all fibre pieces of in such that their characteristic lies in range
As classifying characteristics we can for instance choose the fibres’ local direction (with being the sphere ), their length or curvature (both with ). In this article we focus on local directions of fibres, but the results can easily be applied to other choices of
If the fibre process is stationary then the intensity measure of can be written as where is called the intensity of is the Lebesgue measure in and is a probability measure on which is called the directional distribution of fibres. The distribution is the fibre direction distribution in the typical fibre point, hence length-weighted. In what follows, is either the cardinality of a finite set or the Lebesgue measure of , if is uncountable and measurable.
Let and be the dilation (erosion, resp.) operation on images as introduced e.g. in stoyan . Assume that we observe a dilated version of within a window where is the ball of radius centered at the origin. In our setting, we assume that the fibres’ length is significantly larger than their diameter Moreover, we assume that there is such that is morphologically closed w.r.t. i.e., This condition ensures that the local fibre direction is uniquely defined in each point within
We would like to test the hypothesis
is stationary with intensity and directional distribution vs.
There exists a compact set with and such that
If holds true, the region is called an anomaly region. In the following, we discuss how to test the hypothesis and how to detect the anomaly region .
3 Data and clustering criteria
We assume that the dilated fibre system is observed as a 3D greyscale image. Several methods for estimating the local fibre direction in each fibre pixel are discussed in ImageAnalStereol1489 . We use the approach based on the Hessian matrix that is implemented in the 3D image analysis software tool MAVI mavi . The smoothing parameter required by the method is chosen as , where is an estimate of the (constant) thickness of the typical fibre. In the simulated samples, it is known. For the real data, it is obtained manually from the images.
We divide the observation window into small cubes (see Fig. 1) of the same size, whose edge length equals three times the fibre diameter. The principal axis of local directions (e.g. ImageAnalStereol1489 ) in each , here referred to as “average local direction”, is then computed using the function SubfieldFibreDirections in MAVI.
Let be the regular grid of cubes Some of the cubes may not contain enough fibre voxels to obtain a reliable estimate of the local fibre direction . Let be the subset of indices of cubes which allow for such an estimation. For each point denote by the average local direction estimated from We assume that our fibres are non-oriented and can be then transformed such that where is a hemisphere. The size of this sample is
Our main task is to determine the anomaly regions or, in other words, to classify the set of points into two clusters corresponding to the “homogeneous material” and the “anomaly” (one of these clusters can be empty). To do so, we combine of the small cubes (having edge length ) to a larger cube , such that the 3D image is divided into larger non-empty cubic observation windows (see Figure 1). The larger cubes have side length and the corresponding grid of the larger cubes is denoted by . The set of indexes of non-empty cubes within is denoted by . For each window , we estimate the entropy and the mean of the local directions, based on the estimates as described below.
3.1 Mean of local fibres direction
The vectoris calculated for each window as the coordinate-wise sample mean of local directions (MLD):
Note that and normally
3.2 The entropy of the directional distribution
The entropy of a random variable is a certain measure of the diversity/concentration of its range. Let
be a probability distribution of a random elementon an abstract measurable phase-space The value
is called the Shannon (differential) entropy of . If is absolutely continuous with respect to some measure then there exists the Radon-Nikodym derivative (or density) and the entropy of has the following form
In what follows, we assume that the random variable is absolutely continuous with density In our problem setting, corresponds to the local fibre direction, is the sphere and is the spherical surface area measure on Since our fibres are non-oriented (), we may consider even local direction densities on the whole sphere where appropriate, instead of a density defined on However, choosing another classifying characteristic will lead to a different measurable space
3.3 Entropy estimation
In the literature, there is a large number of papers devoted to non-parametric entropy estimation for i.i.d. random vectors in see e.g. the review in BDGM and references in BD . We will dwell upon two important estimates: the plug-in and the nearest neighbor ones.
3.3.1 Plug-in estimator of entropy
For simplicity, define the plug-in estimator for directional distributions on the sphere with even densities Its general definition on compact manifolds can be found e.g. in Alonzo .
For a directional distribution density take the kernel estimator on a window of the form
where denotes the bandwidth, is a kernel function and is a geodesic metric given by where is the Euclidean scalar product in .
Then the plug-in estimator of in the window is given by
where is the sub-window and denotes the translation of
For homogeneous marked Poisson point processes, the plug-in estimator as above is considered in Alonzo . See also ruiz2016entropy for the context of Boolean models of line segments. We also made an attempt to apply this method to our 3D image data. But we met difficulties which basically come from the relatively small amount of data available. Namely, needs a large number of points in sub-windows during the estimation of together with a large number of such sub-windows. Let us illustrate these difficulties on a simple example.
Consider the uniform distribution on the spherei.e., the density is We generate a sample from this distribution and estimate its entropy using the plug-in estimator (4) with (as in ruiz2016entropy ). The results are presented in Table 1. Moreover, we run the procedure 100 times and compare the obtained values with the exact value of the entropy
Obviously, the bias of is too large with less than 62000 entries, which is in accordance with M stating the impracticability of plug-in entropy estimates for samples in higher dimensions. There are 430741 entries in for the real data (Figure 11) and 463537 entries in RSA fibre data (Figure 3), that allows us to subdivide the images only into 4 non-intersecting regions with more than 100000 cubes In other words, in order to test the hypotheses vs. with test statistics based on estimated entropy (4) we have a sample of size 4, which is too small, compare Section 4, inequality (24). There, for the minimal sample size must be 1000.
3.3.2 Nearest neighbor estimator of entropy
In order to overcome the above difficulties we apply another estimator of introduced in the paper by Kozachenko and Leonenko Leon . We call this estimator “Dobrushin estimator” because its main idea is due to Dobrushin Dobr . The estimator from Leon cannot be applied directly, because it is designed for random vectors in a -dimensional Euclidean space which is flat. In our setting, the phase space is a manifold of positive constant curvature with geodesic metric Therefore, we take a version of Dobrushin estimator for the case of an dimensional compact Riemannian manifold with geodesic metric and Hausdorff measure .
For defining this estimator, the following results will be useful. Denote by the ball in with radius and center i.e., Since a ball and a Euclidean -dimensional ball are bi-Lipschitz equivalent, coincides with the Hausdorff dimension of the manifold see (Fal, , Corollary 2.4). Furthermore, for almost all points the Lebesgue density theorem is true, i.e.,
Let be a sample of i.i.d. valued random elements with continuous density function Denote by the distance to the nearest neighbor of i.e., Define the statistic
where and are defined by (5) and is Euler’s constant. The statistic is called nearest neighbor (Dobrushin) estimator of the entropy.
It coincides with the nearest-neighbour entropy estimate given in (Penrose, , p. 2169) with the only difference that in Penrose Euclidean distances between are used instead of geodesic distances The -consistency of is proven in (Penrose, , Theorem 2.4) for i.i.d. samples as above with bounded density of compact support.
In fact, a large class of parametric distributions on a sphere, including the Fisher, the Watson or the Angular Central Gaussian distribution, has bounded densities with compact support.
In many problems of probability theory, limit theorems for independent observations remain true for weakly dependent data. Since the fibre materials are weakly dependent (the fibers have a finite length), we can assume that the entropy and mean local directions are weakly dependent as well. The proof of consistency of (6) for weakly dependent is non-trivial and goes beyond the scope of this paper.
Our data sets consist of straight fibres which are longer than the edge length of small observation windows Such fibres yield several almost equal values of average local directions This leads to very small values of a distance to the nearest neighbor and, consequently, to the large negative bias of which is computed using Trying to eliminate this bias, we propose to use the following version of (6)
with penalty value found by computational tuning.
In order to test the accuracy of the Dobrushin estimator, we have generated 100 samples from the uniform directional distribution on . We have computed the Dobrushin statistic and compared it with the exact entropy value The results are presented in Table 2.
Based on these results, we conclude that the Dobrushin estimator (6) is quite accurate for small sample sizes. Even for a sample with 64 entries the entropy is estimated much better than by the plug-in method.
4 Change point detection in random fields
To test the hypothesis against we check the existence of anomaly regions in a realization of an dependent geometric random field Here we follow the ideas from Brodsky , where change-point problems for mixing random fields on general parametric (disorder) regions were considered. The field will be chosen in a way such that the hypothesis implies that it is stationary, whereas means the presence of a region with different mean value of Later in our application to fibre materials in Section 4.3, we assume the anomaly region to be a box BucchiaWendler17 .
4.1 Random fields with inhomogeneities in mean
Let be an integrable stationary real-valued random field with Denote by the centered field Moreover, we assume that is dependent, i.e, and are independent if Let be a finite parameter space. For every we define subsets of anomalies completely determined by a parameter Then for some we observe
where and Assume that for every Denote
Let correspond to the values of for anomalies which we consider as significant, i.e., they are neither extremely small nor represent the majority of the data. Formally, for we let
Then corresponds to extremely small or large anomalies, i.e.,
4.2 Testing the change of expectation
Now we can formulate the change-point hypotheses for the random field with respect to its expectation as follows.
for every (i.e. ) vs.
There exists such that and
Consider the following change-in-mean statistics for the sample
the centered field In order to test vs. we use the test statistics
We reject if exceeds the critical value Let us find such via the probability of the 1st-type error It holds
Thus, we find the bounds for tail probabilities of the random variable To do so, we use the ideas from Heinrich to get the following bounds for dependent random fields. For the sake of generality, our results are formulated in
Let be a stationary real-valued dependent centered random field and be real numbers. Assume that there exist such that
Then for any we have
Using the Markov inequality we have for any that
Denote by for It follows from Hölder’s inequality and dependence that
From Taylor’s expansion we have for
For we put in (16) and obtain
This completes the proof.
We apply Theorem 4.1 to First, we rewrite as
Let for and a.s., then under the conditions of Theorem 4.1 we have that
From the definition of we have
Since then Therefore, we put in (12) and obtain the statement of the corollary.
Since we have the following corollary.