Detecting anomalies in fibre systems using 3-dimensional image data

We consider the problem of detecting anomalies in the directional distribution of fibre materials observed in 3D images. We divide the image into a set of scanning windows and classify them into two clusters: homogeneous material and anomaly. Based on a sample of estimated local fibre directions, for each scanning window we compute several classification attributes, namely the coordinate wise means of local fibre directions, the entropy of the directional distribution, and a combination of them. We also propose a new spatial modification of the Stochastic Approximation Expectation-Maximization (SAEM) algorithm. Besides the clustering we also consider testing the significance of anomalies. To this end, we apply a change point technique for random fields and derive the exact inequalities for tail probabilities of a test statistic. The proposed methodology is first validated on simulated images. Finally, it is applied to a 3D image of a fibre reinforced polymer.


Introduction
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 [4,8,14,17,19,52].Changes to be detected may concern the mean, variance, correlation, spectral density, etc. of the (stationary) sequence {X k , k ≥ 0}.This kind of change detection has been considered by various authors starting with [41].Sen and Srivastava [45] considered tests for a change in the mean of a Gaussian model.An overview can also be found in [6].The CUSUM procedure, Bayesian approaches as well as maximum likelihood estimation are often used.Scan statistics come also into play naturally, see e.g.[7,8].
First approaches to change point analysis for random fields (or measures) have been developed in the papers [9,10,13,16,28,31,32,33,36,39,40,46,47,48], see also the review in [7, Section 2, D] and [8,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 m−dependent random fields.In the spirit of the book [8], 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.[1,2].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 [34] for a review of Expectation Maximization (EM) algorithms for the separation of components in a mixture of Gaussian distributions as well as a recent paper [37]).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 changepoint problem for the corresponding m−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.

Problem setting
In this section, we give some basic definitions and results for fibre processes.For more details, see, for example, the book [18].In 3-dimensional Euclidean space, a fibre γ is a simple curve {γ(t) = (γ 1 (t), γ 2 (t), γ 3 (t)), t ∈ [0, 1]} of finite length satisfying the following assumptions: The collection of fibres forms a fibre system φ if it is a union of at most countably many fibres γ (i) , such that any compact set is intersected by only a finite number of fibres, and γ (i) ((0, 1)) ∩ γ (j) ((0, 1)) = ∅, if i = j, 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 is the length of fibre γ in window B. Then φ(B) is the total length of fibre pieces in the window B.
Definition 1 A fibre process Φ is a random element with values in the set D of all fibre systems φ with σ-algebra D generated by sets of the form {φ ∈ D : φ(B) < x} for all bounded Borel sets B and real numbers x.The distribution P of a fibre process is a probability measure on [D, D].The fibre process Φ is said to be stationary if it has the same distribution as the translated fibre process Φ x = Φ + x for all x ∈ R 3 .
For classification needs we consider an abstract fibre characteristic w.Let (E, E, σ) be a measurable space where E is a (compact) Riemannian manifold equipped with a metric ρ.Let w(x) ∈ E be some characteristic of a fibre at point x ∈ R 3 , assuming that exactly one fibre of Φ passes through x.Then a weighted random measure Ψ can be defined by for bounded B ∈ B(R 3 ) and L ∈ E. Thus, Ψ (B × L) is the total length of all fibre pieces of Φ in B such that their characteristic w lies in range L.
As classifying characteristics w we can for instance choose the fibres' local direction (with E being the sphere S 2 ), their length or curvature (both with E = R + ).In this article we focus on local directions of fibres, but the results can easily be applied to other choices of w.
If the fibre process Φ is stationary then the intensity measure of Ψ can be written as EΨ (B × L) = λ|B|f (L), where λ > 0 is called the intensity of Ψ, | • | is the Lebesgue measure in R 3 and f is a probability measure on S 2 which is called the directional distribution of fibres.The distribution f is the fibre direction distribution in the typical fibre point, hence length-weighted.In what follows, |A| is either the cardinality of a finite set A or the Lebesgue measure of A, if A is uncountable and measurable.
Let ⊕ and be the dilation (erosion, resp.)operation on images as introduced e.g. in [18].Assume that we observe a dilated version , where B r is the ball of radius r > 0 centered at the origin.In our setting, we assume that the fibres' length is significantly larger than their diameter 2r.Moreover, we assume that there is ε > 0 such that Ξ is morphologically closed w.r.t.B ε , i.e., (Ξ ⊕ B ε ) B ε = Ξ.This condition ensures that the local fibre direction is uniquely defined in each point within Ξ.
We would like to test the hypothesis H 0 : Φ is stationary with intensity λ and directional distribution f vs. H 1 : There exists a compact set A ⊂ W with |A| > 0 and |W \ A| > 0 such that If H 1 holds true, the region A is called an anomaly region.In the following, we discuss how to test the hypothesis H 0 and how to detect the anomaly region A.

Data and clustering criteria
We assume that the dilated fibre system Ξ ∩W is observed as a 3D greyscale image.Several methods for estimating the local fibre direction w(x) in each fibre pixel x ∈ Ξ are discussed in [51].We use the approach based on the Hessian matrix that is implemented in the 3D image analysis software tool MAVI [26].The smoothing parameter σ required by the method is chosen as σ = r, where r 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 W into small cubes W i (see Fig. 1) of the same size, whose edge length ∆ equals three times the fibre diameter.The principal axis ŵi of local directions (e.g.[51]) in each W i , here referred to as "average local direction", is then computed using the function SubfieldFibreDirections in MAVI.
be the regular grid of cubes W i .Some of the cubes W i may not contain enough fibre voxels to obtain a reliable estimate of the local fibre direction ŵi .Let J ⊂ J T be the subset of indices of cubes which allow for such an estimation.For each point i ∈ J denote by X i = (x i , y i , z i ) T the average local direction estimated from W i .We assume that our fibres are non-oriented and can be then transformed such that where Our main task is to determine the anomaly regions or, in other words, to classify the set of points J into two clusters corresponding to the "homogeneous material" and the "anomaly" (one of these clusters can be empty).To do so, we combine M 3 of the small cubes W i (having edge length ∆) to a larger cube W l , such that the 3D image W is divided into larger non-empty cubic observation windows (see Figure 1).The larger cubes have side length M ∆ and the corresponding grid of the larger cubes is denoted by For each window W l , l ∈ J W , we estimate the entropy and the mean of the local directions, based on the estimates { ŵi , i ∈ S l } as described below.

Mean of local fibres direction
The vector M l = (M x,l , M y,l , M z,l ) T is calculated for each window W l , l ∈ J W as the coordinate-wise sample mean of local directions (MLD):

The entropy of the directional distribution
The entropy of a random variable is a certain measure of the diversity/concentration of its range.Let P be a probability distribution of a random element X on an abstract measurable phase-space (X, σ).The value is called the Shannon (differential) entropy of X.If P is absolutely continuous with respect to some measure σ then there exists the Radon-Nikodym derivative (or density) f = dP dσ , and the entropy of X has the following form In what follows, we assume that the random variable X is absolutely continuous with density f : X → R + .In our problem setting, X corresponds to the local fibre direction, X is the sphere S 2 and σ is the spherical surface area measure on S 2 .Since our fibres are non-oriented (X ∈ S 2 + ), we may consider even local direction densities f on the whole sphere S 2 where appropriate, instead of a density f defined on S 2 + .However, choosing another classifying characteristic w will lead to a different measurable space (X, σ).

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 R D , see e.g. the review in [5] and references in [12].We will dwell upon two important estimates: the plug-in and the nearest neighbor ones.

Plug-in estimator of entropy
For simplicity, define the plug-in estimator for directional distributions on the sphere S 2 with even densities f : Its general definition on compact manifolds X can be found e.g. in [2].
For a directional distribution density f, take the kernel estimator f B (•) on a window B ⊆ J T of the form where h > 0 denotes the bandwidth, K : R + → R is a kernel function and ρ : S 2 × S 2 → R + is a geodesic metric given by ρ(x, y) = arccos x, y , x, y ∈ S 2 , where •, • is the Euclidean scalar product in R 3 .
Then the plug-in estimator of E f in the window W l is given by where B ⊂ J T is the sub-window and B +i := {j+i, j ∈ B} denotes the translation of B.
For homogeneous marked Poisson point processes, the plug-in estimator E f as above is considered in [2].See also [1] 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, E f needs a large number of points in sub-windows B during the estimation of f together with a large number of such sub-windows.Let us illustrate these difficulties on a simple example.
Example 1 Consider the uniform distribution on the sphere S 2 , i.e., the density is f (x) = 1 4π , x ∈ S 2 .We generate a sample from this distribution and estimate its entropy E f using the plug-in estimator (4) with |B| = |S 2 | 1/9 (as in [1]).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 E f is too large with less than 62000 entries, which is in accordance with [38] stating the impracticability of plug-in entropy estimates for samples in higher dimensions.There are 430741 entries in J T for the real data (Figure 11) and 463537 entries in RSA fibre data (Figure 3), that allows us to  4) we have a sample of size 4, which is too small, compare Section 4, inequality (24).There, for m = 1 the minimal sample size |W | must be 1000.

Nearest neighbor estimator of entropy
In order to overcome the above difficulties we apply another estimator of E f introduced in the paper by Kozachenko and Leonenko [35].We call this estimator "Dobrushin estimator" because its main idea is due to Dobrushin [20].The estimator from [35] cannot be applied directly, because it is designed for random vectors in a d-dimensional Euclidean space which is flat.In our setting, the phase space (S 2 , σ) is a manifold of positive constant curvature with geodesic metric ρ.Therefore, we take a version of Dobrushin estimator for the case of an d−dimensional compact Riemannian manifold X with geodesic metric ρ and Hausdorff measure σ.
For defining this estimator, the following results will be useful.Denote by B δ (x) the ball in (X, ρ) with radius δ > 0 and center x, i.e., B δ (x) = {y ∈ X : ρ(x, y) ≤ δ}.Since a ρ−ball and a Euclidean d-dimensional ball are bi-Lipschitz equivalent, d coincides with the Hausdorff dimension of the manifold X, see [24,Corollary 2.4].Furthermore, for σ−almost all points x ∈ X the Lebesgue density theorem is true, i.e., where d = dim H X is the Hausdorff dimension of X and c = 2 d D X > 0, where D X is the Hausdorff density of X, see [24, Proposition 4.1,5.1]and [44,Theorem 30].
Definition 2 Let (ξ 1 , . . ., ξ N ) be a sample of i.i.d.X− valued random elements with continuous density function f : X → R + .Denote by ρ i the distance to the nearest neighbor of ξ i , i = 1, N , i.e., ρ i = min j=1,N ,j =i {ρ(ξ i , ξ j )}.Define the statistic where c and d are defined by ( 5) and γ = − ∞ 0 (ln y)e −y dy ≈ 0.5772 is Euler's constant.The statistic E is called nearest neighbor (Dobrushin) estimator of the entropy.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.
Remark 1 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 ξ i is non-trivial and goes beyond the scope of this paper.
Remark 2 Our data sets consist of straight fibres which are longer than the edge length ∆ of small observation windows W .Such fibres yield several almost equal values of average local directions X i .This leads to very small values of a distance to the nearest neighbor ρ i and, consequently, to the large negative bias of Ê which is computed using log ρ i .Trying to eliminate this bias, we propose to use the following version of ( 6) with penalty value ρ 0 = 0.01 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 S 2 .We have computed the Dobrushin statistic and compared it with the exact entropy value ln(4π) ≈ 2.53.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.

Change point detection in random fields
To test the hypothesis H 0 against H 1 , we check the existence of anomaly regions in a realization of an m−dependent geometric random field {s k , k ∈ W }.Here we follow the ideas from [6], where change-point problems for mixing random fields on general parametric (disorder) regions were considered.The field {s k , k ∈ W } will be chosen in a way such that the hypothesis H 0 implies that it is stationary, whereas H 1 means the presence of a region I θ ⊂ W with different mean value of s k .Later in our application to fibre materials in Section 4.3, we assume the anomaly region to be a box [11].

Random fields with inhomogeneities in mean
Let { ξk , k ∈ Z 3 } be an integrable stationary real-valued random field with µ = E ξk .Denote by ξ k the centered field Let Θ be a finite parameter space.For every θ ∈ Θ we define subsets of anomalies I θ ⊂ Z 3 completely determined by a parameter θ.Then for some θ 0 ∈ Θ we observe where 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 γ 0 , γ 1 ∈ (0, 1), γ 0 < γ 1 , we let Then Θ 1 = Θ \ Θ 0 corresponds to extremely small or large anomalies, i.e.,

Testing the change of expectation
Now we can formulate the change-point hypotheses for the random field {s k , k ∈ W } with respect to its expectation as follows.
H 0 : Es k = µ for every k ∈ W (i.e.h = 0) vs. H 1 : There exists θ 0 ∈ Θ 0 such that Es k = µ + h, k ∈ I θ 0 , h = 0, and Consider the following change-in-mean statistics for the sample S = {s k , k ∈ W } : Denote by the centered field Z(θ).In order to test H 0 vs. H 1 we use the test statistics We reject H 0 if T W (S) exceeds the critical value y α .Let us find such y α > 0 via the probability of the 1st-type error P H 0 (max θ∈Θ 0 |Z(θ)| ≥ y α ) ≤ α.It holds Thus, we find the bounds for tail probabilities of the random variable max θ∈Θ 0 |η(θ)|.
To do so, we use the ideas from [29] to get the following bounds for m−dependent random fields.For the sake of generality, our results are formulated in Z D , D ∈ N.
Theorem 1 Let {ξ k , k ∈ Z D } be a stationary real-valued m−dependent centered random field and {b k , k ∈ Z D } be real numbers.Assume that there exist H, σ > 0 such that Then for any W ⊂ Z D , |W | < ∞ we have Proof Using the Markov inequality we have for any u > 0 that Denote by W (l) = {l + mi ∈ W |i ∈ Z D } for l ∈ {1, . . ., m} D .It follows from Hölder's inequality and m−dependence that From Taylor's expansion we have for u ∈ 0, Combining ( 14) and ( 15) we continue for the first term in (13) with the following bound for 0 The minimum of ( 16) is achieved for u = y/(2m D σ 2 b 2 2 ).Moreover, bound ( 16) is valid for the second term in (13), too.Therefore, for 0 ≤ y ≤ 16) and obtain This completes the proof.
We apply Theorem 1 to η(θ), θ ∈ Θ 0 , D = 3.First, we rewrite η(θ) as s., then under the conditions of Theorem 1 we have that Therefore, we put H = M 0 in (12) and obtain the statement of the corollary.
Simplifying the result of Corollary 2, we get that ( 17) and ( 18) is bounded by Particularly, if y < and if y >

Change-point detection in simulated random fields
In this section, we study the behaviour of the test statistics T W (S) given in (11) and probabilities of 1st-type error P H 0 (T W (S) ≥ y α ) with respect to different values of σ 2 and m.
The form of T W allows to test the existence of the anomaly regions of arbitrary form and arbitrary number of connected components.On the other hand, we need to decrease the value |Θ| up to a feasible quantity for computational reasons (see bounds (20) and ( 21) We fix γ 0 = 0.05 and γ 1 = 0.5, as the anomaly should not cover the majority of the window.In this paper we restrict I θ to be a single rectangular parallelepiped of the form Then the parametric set of significant anomaly regions is given by The offset parameters ∆ 0 and ∆ 1 as well as the parameter L M controlling the minimal edge length of the cuboids have to be chosen by the user.Assuming the m−dependence for our observations, we do not know the exact value of m.Hence, we need to impose some restrictions on the field ξ.First, if we know a-priori the maximum length of a typical fiber we can immediately obtain the bound for m.Second, we can estimate the covariance function of the random field {s k , k ∈ W } and assess m as the range when this empirical covariance is sufficiently close to zero.From relation (21) with < y < M we obtain the following approximate bound for an admissible m : For example, for |Θ 0 | = 10 4 , α = 0.05, γ 0 = 0.05, one gets Let us now compare the empirical probability of the error of the 1-st type with the bounds (19) for P H 0 (T W (•) ≥ y α ).We generate 300 realizations of a Gaussian centered m-dependent (m = 10) random field {Y k , k ∈ W } with W = [1, 80] 3 ∩ N 3 (which is matched to the considered data sets) and Y k ∼ N (0, 1).The dependence is modelled as follows: random variables Y 1+mk , k ≥ 0 are independent, and Y 1+mk = Y l+mk , k ∈ N 3 for all l ∈ {1, . . ., m} 3 .We take ∆ 0 = ∆ 1 = 8, γ 0 = 0.05, and γ 1 = 0.5.In this case, |Θ 0 | = 11954.Based on the simulated sample of values of the test statistics T W (Y ) we compute the empirical critical value ŷα = 0.6396 for α = 0.05.From comparison of ŷα with critical values y α based on inequality (19) with M 0 = σ (presented in Table 3), we see that even under the exact value of σ 2 = 1, critical values y α are quite conservative.For example, y α = 0.7198 for m = 7 is still greater than ŷα = 0.6396 generated for m = 10.Therefore, we can use critical values from inequality (19) with m smaller than its real value.

Cluster based anomaly detection
For processes Ξ of thick fibres introduced in Section 2, the evidence of an anomaly is tested by applying the test of Section 4 to the random field {s k , k ∈ W } of estimated local mean or entropy of the chosen fibre characteristic w.In this paper, w(x) is the average direction vector of the fibres of Ξ at x ∈ W or one of its coordinates (introduced in Section 3 as ŵi ).
Assume that the anomaly test presented in Section 4 rejected the hypothesis H 0 (and hence H 0 ), i.e., we have an evidence of an anomalous fibre behaviour in the rectangular subregion I θ 0 of our image data.Now we are interested in a more accurate estimate of the geometry of this anomal.The search for an anomaly region in a 3D image can be interpreted as a problem of splitting the volume of the image into two disjoint clusters: homogeneous material and anomaly.
In our problem setting, the volume under investigation, l∈J W W l , is a union of scanning windows W l with meaningful local direction information.Each of them yields the clustering attributes mean of local directions (MLD) and entropy.We need to classify all the windows W l as either belonging to the homogeneous material or the anomaly.For this purpose, a spatial version of the Stochastic Expectation Maximization algorithm is used.

Spatial modification of a Stochastic Approximation Expectation Maximization (SAEM) algorithm
We assume that under the alternative H 1 (see page 2), fibres in the material may have two different distributions f 0 and f 1 of local directions.Therefore, the distribution of clustering attributes is a mixture of the distributions f 0 and f 1 .
The Expectation-Maximization algorithm (EM) is commonly used to separate modes in a finite mixture of distributions, cf.[34] for a review.It is an iterative procedure consisting of two steps: Expectation (Estimation) and Maximization.In general, one assumes that the probability law under study is a mixture of k distributions from the same parametric family.In the first step, the hidden parameters of the sample distribution, i.e., the weights of the mixture components, are estimated, while in the second step the resulting parameters are updated by maximizing the likelihood function.
Since the EM algorithm belongs to the so-called "greedy" algorithms, that is, it converges to the first local optimum that has been found, a modification that compensates this deficiency should be used.One way out is a random "shaking" of observations in each iteration.This method is the basis of the Stochastic EM (SEM) algorithm (cf.[27,34]).
The SEM algorithm works relatively fast in comparison with other methods, and its results are non-sensitive to an initial approximation.Random perturbations on the parameter space in the S-step guarantee the convergence to the global maximum of the likelihood function and help to avoid unstable local maxima.On the other hand, the outputs of the SEM algorithm form a Markov chain and the final solution is its stationary distribution.To avoid this additional problem we use a modification called SAEM (Stochastic Approximation of EM) algorithm which brings together advantages of both EM and SEM approaches, e.g.[15].
Assume that the observable distribution has a density of the form where ϕ(x, δ i ) is a multivariate Gaussian density with unknown parameter δ i = (µ i , Σ i ), i = 1, 2 and β ∈ [0, 1].Here µ i is the mean and Σ i is the covariance matrix of Gaussian component i = 1, 2. The combined unknown parameter is δ = (β, δ 1 , δ 2 ).We call ϕ(•, δ 1 ) and ϕ(•, δ 2 ) the first and the second component of the mixture, respectively.For each observation x l , l ∈ J W , we define a new variable y l = I{x l belongs to the first component}.Therefore, we have two samples: observable x = {x l , l ∈ J W } and unobservable y = {y l , l ∈ J W }. Then the loglikelihood function equals ln L(δ, x, y) window W l thus belongs to the zone of homogeneous material if q l ≥ 0.5 and to the anomaly zone if q l < 0.5.If β < 0.5 the roles of the components are swapped.
6 Application to 3D image data of fibre materials First, we illustrate the use of the methods from Sections 4 and 5 on simulated 3D fibre images.We choose a random sequental absorbtion (RSA) model that randomly adds fibres to the existing material, such that they do not intersect each other, cf.[3,43].Figure 3 shows simulated RSA fibre data in an image of 2000 × 2000 × 2100 voxels.The sample exhibits three layers, where fibres differ in their local directional distribution.Each layer has a thickness of 700 voxels and contains 82474 fibres with a constant radius of 4 voxels and length of 100 voxels.Fibre directions are distributed according to a special case of the Angular Central Gaussian distribution described in [25].In the two outer layers, the preferred direction is the x-direction and the concentration parameter is β = 0.1 resulting in a high concentration of the fibres along the main direction.In the middle layer, which is considered the anomaly region, the preferred direction is the y-direction and the fibres are less concentrated (β = 0.5).Here, the concentration parameter of the fibre direction distribution is β = 0.1 in the whole sample.The preferred direction is the x-direction.The fibre radius is 4 voxels, the fibre length is 100 voxels.A visualisation of a realisation of this model is shown in Figure 4. Now we apply the change-point analysis of Section 4 to random fields of mean local directions and entropy estimates for the homogeneous and layered RSA data.
To do so, we transform the data of average local directions X k , k ∈ J.In order to avoid cancelling effect of averaging, we build for each coordinate x, y, z, the samples xk , ỹk , zk , such that their entries lie in the hemispheres S The sample of the estimated entropy values Êk of directional distribution of fibres X i , i ∈ J in the windows W k , k ∈ J W is build by estimator (7) over the transformed directions Xi ∈ S 2 + , i ∈ J.We apply the results of Section 4 consequently to the random fields s k = xk , ỹk , or zk , and Êk .Therefore, we have the following 4 pairs of hypotheses of (H 0 , H 1 )-type. - Since we test only 4 hypotheses simultaneously, we stick to the classical Bonferroni method, e.g.we test each direction and entropy separately with significance level Before running the algorithms we need to choose the right size of scanning windows.From the initial layered and homogeneous RSA images with 2000 × 2000 × 2100 voxels we obtain 83 × 83 × 87 small windows W i with 24 × 24 × 24 voxels each, and 463537 and 460559 nonempty entries, respectively.
Due to the model parameters (fibre length of 100 voxels corresponds to 5 points in W ), we can assume the random field Xk to be m−dependent with m = 5 and σ 2 = 0.2, M 0 = 1 2 .For mean local directions x, ỹ, and z, the parametric set Θ 0 is constructed with ∆ 0 = ∆ 1 = 8, γ 0 = 0.05, γ 1 = 0.5, and L M = 22 in (22), We point out that the samples (x, ỹ, z) and Ê have different sizes due to the construction described in Section 3. Therefore, the parameters m, σ 2 and the parameter set Θ 0 in (22) for Ê differ from the ones for x, ỹ, z.
The computed statistics (given in ( 11)) T W (x), T W (ỹ), T W (z), T W ( Ê) and corresponding p-values from relation (19) are presented in Table 4 for the homogeneous and in Table 5 for the layered RSA data.Thus, there is no evidence to reject H x 0 , H y 0 , H z 0 , H E 0 in the homogeneous case.The described test allows to claim that there is an anomaly region in the layered RSA image data, because we reject H x 0 , H y 0 , and H E 0 , but have no evidence to reject H z 0 .Therefore, the choice of the mean of local directions attribute for the changepoint analysis in our problem with layered fibre image data is reasonable.Moreover, depending on the data (e.g. containing whirlpools of fibres) it may be better to choose entropy or other attributes to test for other types of anomalies.

Attribute
It follows from [ of fibre locations and directions, the estimated local directional entropy E for the homogeneous data seems to have a unimodal distribution, see Figure 6.Assuming that the Gaussian distribution provides a reasonable approximation also in this case, we apply the 3σ-rule with σ 2 being the sample variance of E, compare [1], to find anomaly regions in both Figures 3 and 4.  One can see that all centered entropy values lie in the interval [−3σ; 3σ], so the 3σ−method does not distinguish between homogeneous (Figure 4) and inhomogeneous (Figure 3) images.We conclude that the 3σ-rule for anomaly detection does not work well if the anomaly regions are large enough to produce histograms of the clustering attribute with many modes.
The fact that the distribution of E seems to have two modes might indicate that it is a mixture of two Gaussian distributions.So we apply the Spatial SAEM algorithm from Section 5.1 to separate these modes.By an empirical study, a scanning window W consisting of 5 × 5 × 5 small windows W k was selected, i.e.M = 5.Additionally, we put r = ∆M in (34) by default.
First, let us consider clustering based on the attribute entropy.In the layered data, the Spatial SAEM algorithm finds two clusters and determines the distributional parameters for them.The clustering results are presented in Figurer 7 and 8 .Green labels denote the centers of scanning windows corresponding to the homogeneous fibre material and blue labels mark objects belonging to the anomalous region.As expected, the Spatial SAEM algorithm found no anomaly for homogeneous RSA fibre data.We also ran the Spatial SAEM algorithm with the attributes mean of local direction (MLD) and a vector combining entropy and MLD.The results for the layered data are presented in Figures 9 and 10 .One can see that combination of both attributes gives a more reliable result.Remark 4 The problem of clustering a fibre material into homogeneity and anomaly zones using vector-valued cluster attributes can be solved by a variety of other clustering methods, see the books [23,30,49] for an overview.In addition to the results reported here, we tried the recent AWC algorithm [22].However, the spatial SAEM approach yields better results, cf.[21].Moreover, it does not require a complex parameter tuning and operates fast.
We also tried to use a principal axis of fibre directions as a classification attribute in the described SAEM algorithm.But the results are worse than the ones for the MLD attribute.This effect can be explained by the fact that the distribution of principal axes has its support on a unit sphere and thus cannot be a mixture of Gaussian distributions in R 3 .Therefore, the estimates in the M-step ( 25)-( 27), (30), (31) have to be modified, cf.[25].This task goes, however, beyond the scope of the present article.

Real glass fibre reinforced polymer
Now we apply our anomaly detection approach to a 3D-image of a glass fibre reinforced polymer.The images are provided by the Institute for Composite Materials (IVW) in Kaiserslautern, see Figure 11.For a detailed description of the material we refer to [50].We apply the change point analysis from Section 4 to real data with 970 × 1469 × 1217 voxels and the estimated radius of 3 voxels.We obtain 64 × 97 × 80 small windows W i with 15 × 15 × 15 voxels.
For mean local directions, the parametric set Θ 0 is constructed with ∆ 0 = ∆ 1 = 8, γ 0 = 0.05, γ 1 = 0.5 and L M = 22 in (22), which gives |Θ 0 | = 33004.To choose the suitable value of m for m−dependence we need additional investigation.Under the hypotheses H x 0 , H y 0 , H z 0 the random fields x,ỹ, and z are assumed to be stationary, so that we can estimate their covariance functions.We use the standard approach and estimate e.g.ρ x (h) = Cov(x 1 , x1+h ) as , and x0 and xh are the sample means of x over the index ranges K and K + h respectively.To visualize ρx we compute its maximum values in the following way: ρx,max (i) = max The values of ρx,max (black color), ρy,max (red color), ρz,max (blue color) are given in Figure 13.We choose the value of m in such a way that ρx,max (i) ≤ ε 0 , ρy,max (i) ≤ ε 0 , ρz,max (i) ≤ ε 0 , for all i ≥ m, where ε 0 is a threshold.Here we use ε 0 = 0.04 and obtain m = 9.We have M 0 = 0.5 and assume that σ 2 = 0.   test detects the evidence of anomaly regions in real fibre data at significance level α = 8.4 × 10 −10 .
Similarly to the case of RSA data, the detection of anomalies by the 3σ−rule gives meaningless results.The Spatial SAEM algorithm works much better.Its  We conclude that the Spatial SAEM algorithm with attributes entropy, MLD, and a combination of both produces adequate results.Since the images in Figure 15 and 16 look very similar at first glance, we investigate the results of the Spatial SAEM anomaly detection in more detail.We separate a part of the 3D image (Figure 12) into 9 layers and present the result of clustering for the 1st and 5th layer for the entropy in Figure 17, for MLD in Figure 18, and for the combination of both in Figure 19.
One can observe that the Spatial SAEM algorithm with local entropy attribute detects vortices of fibers in the material as anomaly regions.This is natural since a vortex exhibits a large diversity of fibre directions, and the entropy is a measure of such diversity.The Spatial SAEM anomaly detection using the mean of local fibre directions identifies the central part of the image (cf. Figure 18) as an anomaly region, where the directions of fibres differ from the average throughout the image.Finally, the Spatial SAEM approach using both clustering attributes identifies both vortices of fibres and layers of fibres with principally different main direction, cf. Figure 19.

Fig. 1
Fig. 1 Construction of scanning windows W l .

Fig. 2
Fig. 2 Visualization of simulated material with anomaly zone.

Fig. 5
Fig. 5 Histogram of frequencies of the local entropy of fibre directions and two separated Gaussian probability density functions found by the spatial SAEM algorithm.Layered RSA fibre data, σ = 0.5986.

Fig. 6
Fig. 6 Histogram of frequencies of the local entropy of fibre directions fitted Gaussian probability density.Homogeneous RSA fibre data, σ = 0.2997

Fig. 7
Fig. 7 Anomaly detection in layered RSA fibre data using the local entropy: SAEM algorithm.

Fig. 8
Fig. 8 Anomaly detection in layered RSA fibre data using the local entropy: Spatial SAEM algorithm.

Fig. 9
Fig. 9 Anomaly detection with Spatial SAEM algorithm in the layered RSA fibre data using mean of local directions.

Fig. 10
Fig. 10 Anomaly detection with Spatial SAEM algorithm in the layered RSA fibre data using local entropy and mean of local directions attributes.
2. Moreover, due to simulation experiments in Section 4, we compute critical values for the change-point statistics T W (•) and p−values from inequality (19) with m = 7.For the random field of estimated local entropies we have m = 1 and the parametric set Θ 0 is constructed with ∆ 0 = ∆ 1 = 2, γ 0 = 0.05, γ 1 = 0.5 and L M = 4 in(22), which gives |Θ 0 | = 12366.We assume that σ 2 = 0.5 and M 0 = σ.The result of our change point analysis is presented in

Fig. 12
Fig.123D image of a glass fibre reinforced composite material.The part of the data containing an anomaly (970 × 700 × 660 voxels)

Fig. 15
Fig.15Anomaly detection in the fibre image (Fig.11) using mean of local direction.

Fig. 16 Fig. 17
Fig.16 Anomaly detection in the fibre image (Fig.11) using local entropy and mean of local direction.

Table 1
The plug-in entropy estimator for the uniform directional distribution on S 2 subdivide the images only into 4 non-intersecting regions with more than 100000 cubes W l .In other words, in order to test the hypotheses H 0 vs. H 1 with test statistics based on estimated entropy (

Table 2
[42,Dobrushin estimator(6)for the uniform directional distribution on the sphere.It coincides with the nearest-neighbour entropy estimate given in[42, p. 2169] with the only difference that in[42]Euclidean distances between ξ i are used instead of geodesic distances ρ i .The L 2 -consistency of E is proven in[42, Theorem 2.4]for i.i.d.samples (ξ 1 , . . ., ξ N ) as above with bounded density f of compact support.

Table 3
(19)ical values y 0.05 based on inequality(19)for different values of m and σ.

Table 4
Change-point test for mean local directions of homogeneous RSA data.

Table 5
42, Theorem 2.4.] that the Dobrushin estimator of the entropy of i.i.d.vectors on a C 1 −smooth manifold is asymptotically Gaussian.Although the RSA fibre data do not satisfy the i.i.d.assumption of mutual independence Change-point test for mean local directions of layered RSA data.

Table 6 .
Our change point

Table 6
Change-point test for mean local directions of real data.