Detecting anomalies in fibre systems using 3-dimensional image data

TL;DR

This paper introduces a novel method for detecting anomalies in 3D fibre images by classifying local regions based on fibre direction attributes and applying statistical tests, validated on simulated and real data.

Contribution

It proposes a new spatial SAEM algorithm and a change point test for anomaly detection in 3D fibre systems, combining clustering and significance testing.

Findings

Effective anomaly detection in simulated 3D images

Successful application to real fibre reinforced polymer data

New spatial SAEM algorithm improves classification accuracy

Abstract

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 statistics. The proposed methodology is first validated on simulated images. Finally, it is applied to a 3D image of a fibre reinforced polymer.

Tables

Table 1. Table 1: The plug-in entropy estimator for the uniform directional distribution on 𝕊 2 superscript 𝕊 2 \mathbb{S}^{2}

Sample size	Mean	Variance	MSE
500000	2.476	$2.848\times {10}^{-6}$	0.003
375000	2.456	$5.034\times {10}^{-6}$	0.006
250000	2.418	$6.293\times {10}^{-6}$	0.013
125000	2.309	$7.245\times {10}^{-6}$	0.049
62500	2.099	$2.739\times {10}^{-5}$	0.187
12500	0.981	$8.917\times {10}^{-5}$	2.403
6250	0.359	$1.251\times {10}^{-4}$	4.718
1250	0.008	$4.695\times {10}^{-4}$	6.36

Table 2. Table 2: The Dobrushin estimator ( 6 ) for the uniform directional distribution on the sphere.

Sample size	Mean	Variance	MSE
125	2.51	0.02	0.02
64	2.50	0.03	0.02

Table 3. Table 3: Critical values y 0.05 subscript 𝑦 0.05 y_{0.05} based on inequality ( 19 ) for different values of m 𝑚 m and σ 𝜎 \sigma .

${\sigma }^2$	$m=10$	$m=9$	$m=8$	$m=7$	$m=6$
1	1.0757	1.0492	0.8793	0.7198	0.5711
4	2.1513	2.0985	1.7587	1.4394	1.1423
8	3.0424	2.9677	2.4871	2.0357	1.6154

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

Attribute	Sample var.	Test statistics	$p-$value
$\tilde{x}$	0.04360	0.0344	1.00
$\tilde{y}$	0.03743	0.0130	1.00
$\tilde{z}$	0.03749	0.0146	1.00
$\tilde{E}$	0.08984	0.0942	1.00

Table 5. Table 5: Change-point test for mean local directions of layered RSA data.

Attribute	Sample var.	Test statistics	$p-$value
$\tilde{x}$	0.10592	0.44036	$4.6\times {10}^{-30}$
$\tilde{y}$	0.10948	0.43163	$2.8\times {10}^{-23}$
$\tilde{z}$	0.06151	0.18764	0.301
$\tilde{E}$	0.3583	1.07030	0.00

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

Attribute	Sample var.	Test statistics	$p-$value
$\tilde{x}$	0.04589	0.15995	1.00
$\tilde{y}$	0.06795	0.44733	$2.1\times {10}^{-10}$
$\tilde{z}$	0.07982	0.43383	$1.3\times {10}^{-6}$
$E$	0.30126	0.46811	$3.96\times {10}^{-8}$

Equations

$$\varphi (B)=\sum _{{\gamma }^{(i)}\in \varphi }h({\gamma }^{(i)}\cap B)$$

$$\Psi (B\times L)={\int }_B\mathbb{I}\{w(x)\in L\}\Phi (dx)$$

$$\frac{1}{\lambda |A|}\mathbb{E}{\int }_A\mathbb{I}\{w(x)\in \cdot \}\Phi (dx)$$

$$\ne \frac{1}{\lambda |W\setminus A|}\mathbb{E}{\int }_{W\setminus A}\mathbb{I}\{w(x)\in \cdot \}\Phi (dx).$$

$$\begin{array}{c}{M}_{x,\vec{l}}=\frac{1}{|S_{\vec{l}}|}\sum _{\vec{i}\in S_{\vec{l}}}{\hat{w}}_{i_1},M_{y,\vec{l}}=\frac{1}{|S_{\vec{l}}|}\sum _{\vec{i}\in S_{\vec{l}}}{\hat{w}}_{i_2},\\ M_{z,l}=\frac{1}{|S_{\vec{l}}|}\sum _{\vec{i}\in S_{\vec{l}}}{\hat{w}}_{i_3}.\end{array}$$

$$E_X=-{\int }_{\Omega }\ln (\mathbf{P}(X(\omega )))\mathbf{P}(d\omega )$$

$$E_f:=E_X=-{\int }_{\mathbb{X}}\ln (f(x))f(x)\sigma (dx).$$

$${\hat{f}}_B(y)=\frac{1}{|B|}\sum _{\mathbf{i}\in B\cap J}\frac{|\sin \rho (y,X_{\mathbf{i}})|}{h^2\rho (y,X_{\mathbf{i}})}K\left(\frac{\rho (y,X_{\mathbf{i}})}{h}\right),$$

$${\hat{E}}_{f,\mathbf{l}}=-\frac{1}{|S_{\mathbf{l}}|}\sum _{\mathbf{i}\in S_{\mathbf{l}}}\ln {\hat{f}}_{B+\mathbf{i}}(X_{\mathbf{i}}),$$

$$E_f=-\underset{{\mathbb{S}}^2}{\int }\frac{1}{4\pi }\ln \frac{1}{4\pi }\sigma (dx)=\ln 4\pi \approx \mathrm{2.53.}$$

$$\sigma (B_{\delta }(x))\sim c{\delta }^d,\text{ as }\delta \to 0+,$$

$$\hat{E}=\frac{d}{N}\sum _{i=1}^N\ln {\rho }_i+\ln (c(N-1))+\gamma ,$$

$$\begin{array}{c}\hat{E}=\frac{d}{{\sum }_{i=1}^N\mathbb{I}\{{\rho }_i>{\rho }_0\}}\sum _{i=1}^N\mathbb{I}\{{\rho }_i>{\rho }_0\}\ln {\rho }_i+\ln c\\ +\ln \left(\sum _{i=1}^N\mathbb{I}\{{\rho }_i>{\rho }_0\}-1\right)+\gamma \end{array}$$

$$s_k={\tilde{\xi }}_k+h\mathbb{I}\{k\in I_{{\theta }_0}\},k\in W,$$

$${\Theta }_0=\{\theta \in \Theta :{\gamma }_0|W|\le |I_{\theta }|\le {\gamma }_1|W|\}.$$

$${\Theta }_1=\{\theta \in \Theta :|I_{\theta }|<{\gamma }_0|W|,\text{ or }|I_{\theta }^c|<(1-{\gamma }_1)|W|\}.$$

$$Z(\theta )=\frac{1}{|I_{\theta }|}\sum _{k\in I_{\theta }}{s}_k-\frac{1}{|I_{\theta }^c|}\sum _{k\in I_{\theta }^c}{s}_k$$

$$=\frac{1}{|I_{\theta }|}\sum _{k\in I_{\theta }}({\xi }_k+\mu +h\mathbb{I}\{k\in I_{{\theta }_0}\})$$

$$-\frac{1}{|I_{{\theta }_0}^c|}\sum _{k\in I_{\theta }^c}({\xi }_k+\mu +h\mathbb{I}\{k\in I_{{\theta }_0}\})$$

$$=\frac{1}{|I_{\theta }|}\sum _{k\in I_{\theta }}{\xi }_k-\frac{1}{|I_{\theta }^c|}\sum _{k\in I_{\theta }^c}{\xi }_k+h\left(\frac{|I_{\theta }\cap I_{{\theta }_0}|}{|I_{\theta }|}-\frac{|I_{\theta }^c\cap I_{{\theta }_0}|}{|I_{\theta }^c|}\right).$$

$$\begin{array}{c}\eta (\theta ):=Z(\theta )-\mathbb{E}Z(\theta )\\ =\frac{1}{|I_{\theta }|}\sum _{k\in I_{\theta }}{\xi }_k-\frac{1}{|I_{\theta }^c|}\sum _{k\in I_{\theta }^c}{\xi }_k,\theta \in \Theta \end{array}$$

$$\begin{array}{c}{T}_W(S)=\underset{\theta \in {\Theta }_0}{\max }|Z(\theta )|\\ =\underset{\theta \in {\Theta }_0}{\max }\left|\frac{1}{|I_{\theta }|}\sum _{k\in I_{\theta }}{s}_k-\frac{1}{|I_{\theta }^c|}\sum _{k\in I_{\theta }^c}{s}_k\right|.\end{array}$$

$${\mathbf{P}}_{H_0^{\prime }}(\underset{\theta \in {\Theta }_0}{\max }|Z(\theta )|\ge y_{\alpha })$$

$$={\mathbf{P}}_{H_0^{\prime }}\left(\underset{\theta \in {\Theta }_0}{\max }\left|\eta (\theta )+h\left(\frac{|I_{\theta }\cap I_{{\theta }_0}|}{|I_{\theta }|}-\frac{|I_{\theta }^c\cap I_{{\theta }_0}|}{|I_{\theta }^c|}\right)\right|\ge y_{\alpha }\right)$$

$$=\mathbf{P}(|\underset{\theta \in {\Theta }_0}{\max }|\eta (\theta )|\ge y_{\alpha }).$$

$$\mathbb{E}|{\xi }_k^p|\le \frac{p!}{2}{H}^{p-2}{\sigma }^2,p=2,3,\dots$$

$$\mathbf{P}\left(\left|\sum _{k\in W}{\xi }_k{b}_k\right|\ge y\right)\le 2\exp \left(-\frac{y^2}{4m^D{\sigma }^2\Vert b{\Vert }_2^2}\right)$$

$$\times \mathbb{I}\left\{0<y\le \frac{{\sigma }^2\Vert b{\Vert }_2^2}{H\Vert b{\Vert }_{\infty }}\right\}$$

$$+2\exp \left(-\frac{y}{2H{m}^D\Vert b{\Vert }_{\infty }}+\frac{{\sigma }^2\Vert b{\Vert }_2^2}{4H^2{m}^D\Vert b{\Vert }_{\infty }^2}\right)$$

$$\times \mathbb{I}\left\{y>\frac{{\sigma }^2\Vert b{\Vert }_2^2}{H\Vert b{\Vert }_{\infty }}\right\},$$

$$\mathbf{P}\left(\left|\sum _{k\in W}{\xi }_k{b}_k\right|\ge y\right)=\mathbf{P}\left(\sum _{k\in W}{\xi }_k{b}_k\ge y\right)$$

$$+\mathbf{P}\left(\sum _{k\in W}{\xi }_k(-b_k)\ge y\right)\le \frac{\mathbb{E}\exp (u{\sum }_{k\in W}{\xi }_k{b}_k)}{\exp (uy)}$$

$$+\frac{\mathbb{E}\exp (u{\sum }_{k\in W}{\xi }_k(-b_k))}{\exp (uy)}.$$

$$\mathbb{E}\exp \left(u\sum _{k\in W}{\xi }_k{b}_k\right)$$

$$=\mathbb{E}\left(\prod _{l\in \{1,\dots ,m{\}}^D}\exp \left(u\sum _{k\in W(l)}{\xi }_k{b}_k\right)\right)$$

$$\le \prod _{l\in \{1,\dots ,m{\}}^D}{\left(\mathbb{E}\exp \left(m^{D}u\sum _{k\in W(l)}{\xi }_k{b}_k\right)\right)}^{1/m^D}$$

$$=\prod _{l\in \{1,\dots ,m{\}}^D}{\left(\prod _{k\in W(l)}\mathbb{E}{e}^{m^{D}u{\xi }_k{b}_k}\right)}^{1/m^D}$$

Full text

1 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 BrassNik93 ; BrodDarkh00 ; Carlstein_ed94 ; ChenGupta12 ; TimeS1 ; Wu05 . Changes to be detected may concern the mean, variance, correlation, spectral density, etc. of the (stationary) sequence $\{X_{k},k\geq 0\}$. 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 $m-$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 paper Laurent ). 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 $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.

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 $\gamma$ is a simple curve $\{\gamma(t)$ $=(\gamma_{1}(t),\gamma_{2}(t),\gamma_{3}(t)),$ $t\in[0,1]\}$ of finite length satisfying the following assumptions:

• •

$\{\gamma(t),t\in[0,1]\}$ is a $C^{1}$-smooth function.

• •

$\|\gamma^{\prime}(t)\|_{3}^{2}>0$ for all $t\in[0,1],$ where $\|\gamma^{\prime}(t)\|_{3}^{2}=|\gamma_{1}^{\prime}(t)|^{2}+|\gamma_{2}^{\prime}(t)|^{2}+|\gamma_{3}^{\prime}(t)|^{2}.$

• •

A fibre does not intersect itself.

The collection of fibres forms a fibre system $\phi$ if it is a union of at most countably many fibres $\gamma^{(i)},$ such that any compact set is intersected by only a finite number of fibres, and $\gamma^{(i)}((0,1))\cap\gamma^{(j)}((0,1))=\varnothing,$ if $i\neq j,$ i.e., the distinct fibres may have only end-points in common. The length measure corresponding to the fibre system $\phi$ (and denoted by the same symbol) is defined by

[TABLE]

for bounded Borel sets $B\in\mathcal{B}(\mathbb{R}^{3}),$ where $h(\gamma\cap B)=\int_{0}^{1}\mathbb{I}\{\gamma(t)\in B\}\sqrt{|\gamma^{\prime}(t)|^{2}}dt$ is the length of fibre $\gamma$ in window $B.$ Then $\phi(B)$ is the total length of fibre pieces in the window $B.$

Definition 1

A fibre process $\Phi$ is a random element with values in the set $\mathbb{D}$ of all fibre systems $\phi$ with $\sigma$-algebra $\mathcal{D}$ generated by sets of the form $\{\phi\in\mathbb{D}:\phi(B)<x\}$ for all bounded Borel sets $B$ and real numbers $x$. The distribution $P$ of a fibre process is a probability measure on $[\mathbb{D},\mathcal{D}].$ The fibre process $\Phi$ is said to be stationary if it has the same distribution as the translated fibre process $\Phi_{x}=\Phi+x$ for all $x\in\mathbb{R}^{3}.$

For classification needs we consider an abstract fibre characteristic $w$. Let $(E,\mathcal{E},\sigma)$ be a measurable space where $E$ is a (compact) Riemannian manifold equipped with a metric $\rho.$ Let $w(x)\in E$ be some characteristic of a fibre at point $x\in\mathbb{R}^{3},$ assuming that exactly one fibre of $\Phi$ passes through $x.$ Then a weighted random measure $\Psi$ can be defined by

[TABLE]

for bounded $B\in\mathcal{B}(\mathbb{R}^{3})$ and $L\in\mathcal{E}.$ Thus, $\Psi(B\times L)$ is the total length of all fibre pieces of $\Phi$ 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 $\mathbb{S}^{2}$), their length or curvature (both with $E=\mathbb{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 $\Phi$ is stationary then the intensity measure of $\Psi$ can be written as $\mathbb{E}\Psi(B\times L)=\lambda|B|f(L),$ where $\lambda>0$ is called the intensity of $\Psi,$ $|\cdot|$ is the Lebesgue measure in $\mathbb{R}^{3}$ and $f$ is a probability measure on $\mathbb{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 $\oplus$ and $\ominus$ be the dilation (erosion, resp.) operation on images as introduced e.g. in stoyan . Assume that we observe a dilated version $\Xi=\Phi\oplus B_{r}$ of $\Phi$ within a window $W=[a_{1},b_{1}]\times[a_{2},b_{2}]\times[a_{3},b_{3}],$ $a_{i}<b_{i},i=1,2,3,$ 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 $\varepsilon>0$ such that $\Xi$ is morphologically closed w.r.t. $B_{\varepsilon},$ i.e., $(\Xi\oplus B_{\varepsilon})\ominus B_{\varepsilon}=\Xi.$ This condition ensures that the local fibre direction is uniquely defined in each point within $\Xi.$

We would like to test the hypothesis

• $H_{0}:$ $\Phi$ is stationary with intensity $\lambda$ and directional distribution $f$ vs.
• $H_{1}:$ There exists a compact set $A\subset W$ with $|A|>0$ and $|W\setminus A|>0$ such that

[TABLE]

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$.

3 Data and clustering criteria

We assume that the dilated fibre system $\Xi\cap W$ is observed as a 3D greyscale image. Several methods for estimating the local fibre direction $w(x)$ in each fibre pixel $x\in\Xi$ 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 $\sigma$ required by the method is chosen as $\sigma=\hat{r}$, where $\hat{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 $\widetilde{W}_{\vec{i}}$ (see Fig. 1) of the same size, whose edge length $\Delta$ equals three times the fibre diameter. The principal axis $\hat{w}_{\vec{i}}$ of local directions (e.g. ImageAnalStereol1489 ) in each $\widetilde{W}_{\vec{i}}$, here referred to as “average local direction”, is then computed using the function SubfieldFibreDirections in MAVI.

Let $J_{T}=\{\mathbf{i}=(i_{1},i_{2},i_{3}),i_{1}=\overline{1,n_{1}},i_{2}=\overline{1,n_{2}},i_{3}=\overline{1,n_{3}}\}$ be the regular grid of cubes $\widetilde{W}_{\vec{i}}.$ Some of the cubes $\widetilde{W}_{\vec{i}}$ may not contain enough fibre voxels to obtain a reliable estimate of the local fibre direction $\hat{w}_{\vec{i}}$. Let $J\subset J_{T}$ be the subset of indices of cubes which allow for such an estimation. For each point $\vec{i}\in J$ denote by $X_{\vec{i}}=(x_{\vec{i}},y_{\vec{i}},z_{\vec{i}})^{T}$ the average local direction estimated from $\widetilde{W}_{\mathbf{i}}.$ We assume that our fibres are non-oriented and can be then transformed such that $X_{\vec{i}}\in\mathbb{S}^{2}_{+},$ where $\mathbb{S}^{2}_{+}=\{(x,y,z)\in\mathbb{R}^{3}:x^{2}+y^{2}+z^{2}=1,z\in[0,1]\}$ is a hemisphere. The size of this sample is $N=|J|\leq n_{1}n_{2}n_{3}.$

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 $\widetilde{W}_{\vec{i}}$ (having edge length $\Delta$) to a larger cube $W_{\vec{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\Delta$ and the corresponding grid of the larger cubes is denoted by $J_{W}=\{\vec{l}=(l_{1},l_{2},l_{3}):l_{1}=\overline{1,m_{1}},l_{2}=\overline{1,m_{2}},l_{3}=\overline{1,m_{3}}\}$. The set of indexes of non-empty cubes $\widetilde{W}_{\vec{i}}$ within $W_{\vec{l}}$ is denoted by $S_{\vec{l}}=\{(i_{1},i_{2},i_{3})\in J,\widetilde{W}_{(i_{1},i_{2},i_{3})}\subset W_{\vec{l}}\}$. For each window $W_{\vec{l}},\vec{l}\in J_{W}$, we estimate the entropy and the mean of the local directions, based on the estimates $\{\hat{w}_{\vec{i}},\vec{i}\in S_{\vec{l}}\}$ as described below.

3.1 Mean of local fibres direction

The vector $M_{\vec{l}}=(M_{x,\vec{l}},M_{y,\vec{l}},M_{z,\vec{l}})^{T}$ is calculated for each window $W_{\vec{l}},\vec{l}\in J_{W}$ as the coordinate-wise sample mean of local directions (MLD):

[TABLE]

Note that $|S_{\vec{l}}|\leq M^{3}$ and normally $|S_{\vec{l}}|\approx M^{3}.$

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 $\mathbf{P}$ be a probability distribution of a random element $X$ on an abstract measurable phase-space $(\mathbb{X},\sigma).$ The value

[TABLE]

is called the Shannon (differential) entropy of $X$. If $\mathbf{P}$ is absolutely continuous with respect to some measure $\sigma$ then there exists the Radon-Nikodym derivative (or density) $f=\frac{d\mathbf{P}}{d\sigma},$ and the entropy of $X$ has the following form

[TABLE]

In what follows, we assume that the random variable $X$ is absolutely continuous with density $f:\mathbb{X}\to\mathbb{R}_{+}.$ In our problem setting, $X$ corresponds to the local fibre direction, $\mathbb{X}$ is the sphere $\mathbb{S}^{2}$ and $\sigma$ is the spherical surface area measure on $\mathbb{S}^{2}.$ Since our fibres are non-oriented ($X\in\mathbb{S}^{2}_{+}$), we may consider even local direction densities $f$ on the whole sphere $\mathbb{S}^{2}$ where appropriate, instead of a density $f$ defined on $\mathbb{S}^{2}_{+}.$ However, choosing another classifying characteristic $w$ will lead to a different measurable space $(\mathbb{X},\sigma).$

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 $\mathbb{R}^{D},$ 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 $\mathbb{S}^{2}$ with even densities $f:$ $f(x)=f(-x),x\in\mathbb{S}^{2}.$ Its general definition on compact manifolds $\mathbb{X}$ can be found e.g. in Alonzo .

For a directional distribution density $f,$ take the kernel estimator $\widehat{f}_{B}(\cdot)$ on a window $B\subseteq J_{T}$ of the form

[TABLE]

where $h>0$ denotes the bandwidth, $K:\mathbb{R}_{+}\to\mathbb{R}$ is a kernel function and $\rho:\mathbb{S}^{2}\times\mathbb{S}^{2}\to\mathbb{R}_{+}$ is a geodesic metric given by $\rho(\mathbf{x},\mathbf{y})=\arccos\langle\mathbf{x},\mathbf{y}\rangle,\mathbf{x},\mathbf{y}\in\mathbb{S}^{2},$ where $\langle\cdot,\cdot\rangle$ is the Euclidean scalar product in $\mathbb{R}^{3}$.

Then the plug-in estimator of $E_{f}$ in the window $W_{\mathbf{l}}$ is given by

[TABLE]

where $B\subset J_{T}$ is the sub-window and $B+\mathbf{i}:=\{\mathbf{j}+\mathbf{i},\mathbf{j}\in B\}$ denotes the translation of $B.$

For homogeneous marked Poisson point processes, the plug-in estimator $\widehat{E}_{f}$ 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, $\widehat{E}_{f}$ needs a large number of points in sub-windows $B^{\prime}$ 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 $\mathbb{S}^{2},$ i.e., the density is $f(x)=\frac{1}{4\pi},x\in\mathbb{S}^{2}.$ We generate a sample from this distribution and estimate its entropy $E_{f}$ using the plug-in estimator (4) with $|B|=|\mathbb{S}^{2}|^{1/9}$ (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

[TABLE]

Obviously, the bias of $\widehat{E}_{f}$ 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 $J_{T}$ 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 $W_{\mathbf{l}}.$ In other words, in order to test the hypotheses $H_{0}$ vs. $H_{1}$ 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 $m=1$ the minimal sample size $|W|$ must be 1000.

3.3.2 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 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 $d$-dimensional Euclidean space which is flat. In our setting, the phase space $(\mathbb{S}^{2},\sigma)$ is a manifold of positive constant curvature with geodesic metric $\rho.$ Therefore, we take a version of Dobrushin estimator for the case of an $d-$dimensional compact Riemannian manifold $\mathbb{X}$ with geodesic metric $\rho$ and Hausdorff measure $\sigma$.

For defining this estimator, the following results will be useful. Denote by $B_{\delta}(x)$ the ball in $(\mathbb{X},\rho)$ with radius $\delta>0$ and center $x,$ i.e., $B_{\delta}(x)=\{y\in\mathbb{X}:\rho(x,y)\leq\delta\}.$ Since a $\rho-$ball and a Euclidean $d$-dimensional ball are bi-Lipschitz equivalent, $d$ coincides with the Hausdorff dimension of the manifold $\mathbb{X},$ see (Fal, , Corollary 2.4). Furthermore, for $\sigma-$almost all points $x\in\mathbb{X}$ the Lebesgue density theorem is true, i.e.,

[TABLE]

where $d=\dim_{H}X$ is the Hausdorff dimension of $\mathbb{X}$ and $c=2^{d}D_{X}>0,$ where $D_{X}$ is the Hausdorff density of $\mathbb{X}$, see (Fal, , Proposition 4.1,5.1) and (rogers, , Theorem 30).

Definition 2

Let $(\xi_{1},\ldots,\xi_{N})$ be a sample of i.i.d. $\mathbb{X}-$ valued random elements with continuous density function $f:\mathbb{X}\to\mathbb{R}_{+}.$ Denote by $\rho_{i}$ the distance to the nearest neighbor of $\xi_{i},i=\overline{1,N},$ i.e., $\rho_{i}=\min_{j=\overline{1,N},j\neq i}\{\rho(\xi_{i},\xi_{j})\}.$ Define the statistic

[TABLE]

where $c$ and $d$ are defined by (5) and $\gamma=-\int_{0}^{\infty}(\ln y)e^{-y}dy\approx 0.5772$ is Euler’s constant. The statistic $\widehat{E}$ 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 $\xi_{i}$ are used instead of geodesic distances $\rho_{i}.$ The $L_{2}$-consistency of $\widehat{E}$ is proven in (Penrose, , Theorem 2.4) for i.i.d. samples $(\xi_{1},\ldots,\xi_{N})$ as above with bounded density $f$ 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.

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 $\xi_{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 $\Delta$ of small observation windows $\widetilde{W}.$ Such fibres yield several almost equal values of average local directions $X_{\mathbf{i}}.$ This leads to very small values of a distance to the nearest neighbor $\rho_{i}$ and, consequently, to the large negative bias of $\hat{E}$ which is computed using $\log\rho_{i}.$ Trying to eliminate this bias, we propose to use the following version of (6)

[TABLE]

with penalty value $\rho_{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 $\mathbb{S}^{2}$. We have computed the Dobrushin statistic and compared it with the exact entropy value $\ln(4\pi)\approx 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.

4 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\in W\}.$ Here we follow the ideas from Brodsky , where change-point problems for mixing random fields on general parametric (disorder) regions were considered. The field $\{s_{k},k\in 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_{\theta}\subset 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 BucchiaWendler17 .

4.1 Random fields with inhomogeneities in mean

Let $\{\tilde{\xi}_{k},k\in\mathbb{Z}^{3}\}$ be an integrable stationary real-valued random field with $\mu=\mathbb{E}\tilde{\xi}_{k}.$ Denote by $\xi_{k}$ the centered field $\xi_{k}=\tilde{\xi}_{k}-\mu,k=(k_{1},k_{2},k_{3})\in\mathbb{Z}^{3}.$ Moreover, we assume that $\{\xi_{k},k\in\mathbb{Z}^{3}\}$ is $m-$dependent, i.e, $\xi_{k}$ and $\xi_{l}$ are independent if $\max_{i=1,2,3}|k_{i}-l_{i}|>m,k=(k_{1},k_{2},k_{3})\in\mathbb{Z}^{3},l=(l_{1},l_{2},l_{3})\in\mathbb{Z}^{3}.$ Let $\Theta$ be a finite parameter space. For every $\theta\in\Theta$ we define subsets of anomalies $I_{\theta}\subset\mathbb{Z}^{3}$ completely determined by a parameter $\theta.$ Then for some $\theta_{0}\in\Theta$ we observe

[TABLE]

where $W=[1,M_{1}]\times[1,M_{2}]\times[1,M_{3}]\cap\mathbb{Z}^{3},$ and $h\in\mathbb{R}.$ Assume that $I_{\theta}\subset W$ for every $\theta\in\Theta.$ Denote $I_{\theta}^{c}=W\setminus I_{\theta}.$

Let $\Theta_{0}$ correspond to the values of $\theta$ for anomalies which we consider as significant, i.e., they are neither extremely small nor represent the majority of the data. Formally, for $\gamma_{0},\gamma_{1}\in(0,1),$ $\gamma_{0}<\gamma_{1},$ we let

[TABLE]

Then $\Theta_{1}=\Theta\setminus\Theta_{0}$ corresponds to extremely small or large anomalies, i.e.,

[TABLE]

4.2 Testing the change of expectation

Now we can formulate the change-point hypotheses for the random field $\{s_{k},k\in W\}$ with respect to its expectation as follows.

• $H_{0}^{\prime}:$ $\mathbb{E}s_{k}=\mu$ for every $k\in W$ (i.e. $h=0$) vs.
• $H_{1}^{\prime}:$ There exists $\theta_{0}\in\Theta_{0}$ such that $\mathbb{E}s_{k}=\mu+h,k\in I_{\theta_{0}},~{}h\neq 0,$ and $\mathbb{E}s_{k}=\mu,k\in I_{\theta^{c}_{0}}.$

Consider the following change-in-mean statistics for the sample $S=\{s_{k},k\in W\}:$

[TABLE]

Denote by

[TABLE]

the centered field $Z(\theta).$ In order to test $H^{\prime}_{0}$ vs. $H^{\prime}_{1}$ we use the test statistics

[TABLE]

We reject $H_{0}^{\prime}$ if $T_{W}(S)$ exceeds the critical value $y_{\alpha}.$ Let us find such $y_{\alpha}>0$ via the probability of the 1st-type error $\mathbf{P}_{H_{0}^{\prime}}(\max_{\theta\in\Theta_{0}}|Z(\theta)|\geq y_{\alpha})\leq\alpha.$ It holds

[TABLE]

Thus, we find the bounds for tail probabilities of the random variable $\max_{\theta\in\Theta_{0}}|\eta(\theta)|.$ To do so, we use the ideas from Heinrich to get the following bounds for $m-$dependent random fields. For the sake of generality, our results are formulated in $\mathbb{Z}^{D},D\in\mathbb{N}.$

Theorem 4.1

Let $\{\xi_{k},k\in\mathbb{Z}^{D}\}$ be a stationary real-valued $m-$dependent centered random field and $\{b_{k},k\in\mathbb{Z}^{D}\}$ be real numbers. Assume that there exist $H,\sigma>0$ such that

[TABLE]

Then for any $W\subset\mathbb{Z}^{D},$ $|W|<\infty$ we have

[TABLE]

where $\|b\|_{\infty}=\max_{k\in W}|b_{k}|$ and $\|b\|_{2}^{2}=\sum_{k\in W}b^{2}_{k}.$

Proof

Using the Markov inequality we have for any $u>0$ that

[TABLE]

Denote by $W(l)=\{l+mi\in W|i\in\mathbb{Z}^{D}\}$ for $l\in\{1,\ldots,m\}^{D}.$ It follows from Hölder’s inequality and $m-$dependence that

[TABLE]

From Taylor’s expansion we have for $u\in\left[0,\frac{1}{2Hm^{D}|b_{k}|)^{-1}}\right]$

[TABLE]

Combining (14) and (15) we continue for the first term in (13) with the following bound for $0\leq u\leq$ $(2Hm^{D}\max_{k\in W}|b_{k}|)^{-1}:$

[TABLE]

The minimum of (16) is achieved for $u=y/(2m^{D}\sigma^{2}\|b\|_{2}^{2}).$ Moreover, bound (16) is valid for the second term in (13), too. Therefore, for $0\leq y\leq\frac{\sigma^{2}\|b\|_{2}^{2}}{H\|b\|_{\infty}}$ we have

[TABLE]

For $y>\frac{\sigma^{2}\|b\|_{2}^{2}}{H\|b\|_{\infty}}$ we put $u=(2Hm^{D}\|b\|_{\infty})^{-1}$ in (16) and obtain

[TABLE]

This completes the proof.

We apply Theorem 4.1 to $\eta(\theta),\theta\in\Theta_{0},D=3.$ First, we rewrite $\eta(\theta)$ as

[TABLE]

where

[TABLE]

Corollary 1

Let $|I_{\theta}|\leq|I_{\theta}^{c}|$ for $\theta\in\Theta_{0}$ and $|\xi_{k}|\leq M_{0},k\in W$ a.s., then under the conditions of Theorem 4.1 we have that

[TABLE]

Proof

From the definition of $b_{k}(\theta)$ we have

[TABLE]

and

[TABLE]

Since $|\xi_{k}|\leq M_{0}$ then $\mathbb{E}|\xi_{k}^{p}|\leq\mathbb{E}(\xi^{2}_{k}|\xi_{k}|^{p-2})\leq M_{0}^{p-2}\mathbb{E}\xi^{2}_{k}$ $\leq M_{0}^{p-2}\sigma^{2}.$ Therefore, we put $H=M_{0}$ in (12) and obtain the statement of the corollary.

Since $\mathbf{P}\left(\max_{\theta\in\Theta_{0}}|\eta(\theta)|\geq y\right)\leq\sum_{\theta\in\Theta_{0}}\mathbf{P}\left(|\eta(\theta)|\geq y\right),$ we have the following corollary.

Corollary 2

Let $0<\gamma_{0}<\gamma_{1}<1/2$ and assume that the conditions of Theorem 4.1 and Corollary 1 are satisfied. Then

[TABLE]

Hence, we reject $H_{0}^{\prime}$ if the test statistic $T_{W}(S)\geq y_{\alpha},$ where critical value $y_{\alpha}$ is the minimum positive number such that

[TABLE]

Remark 3

In the case of Gaussian random field $\{\xi_{k},k\in\mathbb{Z}^{D}\},$ we can obtain the same statements of Corollaries (1) and (2) with $M_{0}=\sigma.$ Indeed, $\mathbb{E}|\xi_{0}|^{p}\leq\sigma^{p-2}\sigma^{2}\mathbb{E}|\zeta|^{p-2},$ $p\geq 2,$ where $\zeta\sim N(0,1).$ So we put $H=M_{0}$ in (12).

Simplifying the result of Corollary 2, we get that (17) and (18) is bounded by

[TABLE]

Particularly, if $y<\frac{\sigma^{2}}{M_{0}(1-\gamma_{0})}$ then

[TABLE]

and if $y>\frac{\sigma^{2}}{M_{0}(1-\gamma_{1})}$ then

[TABLE]

4.3 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 $\mathbf{P}_{H_{0}^{\prime}}(T_{W}(S)\geq y_{\alpha})$ with respect to different values of $\sigma^{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 $|\Theta|$ up to a feasible quantity for computational reasons (see bounds (20) and (21)). Let $W=[1,M_{1}]\times[1,M_{2}]\times[1,M_{3}]\cap\mathbb{N}^{3}.$ We fix $\gamma_{0}=0.05$ and $\gamma_{1}=0.5,$ as the anomaly should not cover the majority of the window. In this paper we restrict $I_{\theta}$ to be a single rectangular parallelepiped of the form $[1+\Delta_{0}i_{1},1+\Delta_{0}i_{1}+\Delta_{1}l_{1}]\times[1+\Delta_{0}i_{2},1+\Delta_{0}i_{2}+\Delta_{1}l_{2}]\times[1+\Delta_{0}i_{3},1+\Delta_{0}i_{3}+\Delta_{1}l_{3}]$. Then the parametric set of significant anomaly regions is given by

[TABLE]

The offset parameters $\Delta_{0}$ and $\Delta_{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 $\xi$. 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\in W\}$ and assess $m$ as the range when this empirical covariance is sufficiently close to zero. From relation (21) with $\frac{\sigma^{2}}{M_{0}(1-\gamma_{1})}<y<M$ we obtain the following approximate bound for an admissible $m:$

[TABLE]

For example, for $|\Theta_{0}|=10^{4},$ $\alpha=0.05,$ $\gamma_{0}=0.05,$ one gets

[TABLE]

Let us now compare the empirical probability of the error of the 1-st type with the bounds (19) for $\mathbf{P}_{H_{0}^{\prime}}(T_{W}(\cdot)\geq y_{\alpha}).$ We generate 300 realizations of a Gaussian centered $m$-dependent ($m=10$) random field $\{Y_{k},k\in W\}$ with $W=[1,80]^{3}\cap\mathbb{N}^{3}$ (which is matched to the considered data sets) and $Y_{k}\sim N(0,1).$ The dependence is modelled as follows: random variables $Y_{1+mk},k\geq 0$ are independent, and $Y_{1+mk}=Y_{l+mk},k\in\mathbb{N}^{3}$ for all $l\in\{1,\ldots,m\}^{3}.$ We take $\Delta_{0}=\Delta_{1}=8,$ $\gamma_{0}=0.05,$ and $\gamma_{1}=0.5.$ In this case, $|\Theta_{0}|=11954.$ Based on the simulated sample of values of the test statistics $T_{W}(Y)$ we compute the empirical critical value $\hat{y}_{\alpha}=0.6396$ for $\alpha=0.05.$ From comparison of $\hat{y}_{\alpha}$ with critical values $y_{\alpha}$ based on inequality (19) with $M_{0}=\sigma$ (presented in Table 3), we see that even under the exact value of $\sigma^{2}=1,$ critical values $y_{\alpha}$ are quite conservative. For example, $y_{\alpha}=0.7198$ for $m=7$ is still greater than $\hat{y}_{\alpha}=0.6396$ generated for $m=10.$

Therefore, we can use critical values from inequality (19) with $m$ smaller than its real value.

5 Cluster based anomaly detection

For processes $\Xi$ 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\in 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 $\Xi$ at $x\in W$ or one of its coordinates (introduced in Section 3 as $\hat{w}_{i}$).

Assume that the anomaly test presented in Section 4 rejected the hypothesis $H_{0}^{\prime}$ (and hence $H_{0}$), i.e., we have an evidence of an anomalous fibre behaviour in the rectangular subregion $I_{\theta_{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, $\bigcup_{\vec{l}\in J_{W}}W_{\vec{l}}$, is a union of scanning windows $W_{\vec{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_{\vec{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.

5.1 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. EM 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. SEM ; EM ).

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. Celeux92 .

Assume that the observable distribution has a density of the form

[TABLE]

where $\varphi(x,\delta_{i})$ is a multivariate Gaussian density with unknown parameter $\delta_{i}=(\mu_{i},\Sigma_{i}),$ $i=1,2$ and $\beta\in[0,1].$ Here $\mu_{i}$ is the mean and $\Sigma_{i}$ is the covariance matrix of Gaussian component $i=1,2.$ The combined unknown parameter is $\vec{\delta}=(\beta,\delta_{1},\delta_{2}).$ We call $\varphi(\cdot,\delta_{1})$ and $\varphi(\cdot,\delta_{2})$ the first and the second component of the mixture, respectively. For each observation $x_{\mathbf{l}},\mathbf{l}\in J_{W},$ we define a new variable $y_{\mathbf{l}}=\mathbb{I}\{x_{\mathbf{l}}$ belongs to the first component$\}$. Therefore, we have two samples: observable $\mathbf{x}=\{x_{\mathbf{l}},\mathbf{l}\in J_{W}\}$ and unobservable $\mathbf{y}=\{y_{\mathbf{l}},\mathbf{l}\in J_{W}\}.$ Then the log-likelihood function equals

[TABLE]

where $\nu_{1}=\sum_{\mathbf{l}\in J_{W}}y_{\mathbf{l}}$ denotes the number of observations belonging to the first mixture component and $\nu_{2}=m_{1}m_{2}m_{3}-\nu_{1}$ observations belong to the second one.

Assume that we know the a posteriori probability $q_{\mathbf{l}}^{(k-1)},\mathbf{l}\in J_{W},$ that $x_{\mathbf{l}}$ belongs to the first component $\varphi(\cdot,\delta_{1}),$ where $k-1$ is the iteration number. Let us describe the EM part. During the M-Step we obtain new estimates of the parameters $\hat{\delta}^{(k)}_{1}=(\mu^{(k)}_{1},\Sigma^{(k)}_{1}),$ $\hat{\delta}^{(k)}_{2}=(\mu^{(k)}_{2},$ $\Sigma^{(k)}_{2}),\hat{\beta}^{(k)}$ by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In the E-step we compute the new probabilities based on (25)-(28) as

[TABLE]

In the SEM-part we act in a different way. In the S-step we generate independent Bernoulli-distributed random variables $y_{\mathbf{l}}^{(k)}\in\{0,1\}$ with probabilities$\mathbf{P}(y_{\mathbf{l}}^{k}=1)=q_{\mathbf{l}}^{(k-1)},$ $\mathbf{l}\in J_{W}$.

During the M-Step we get $\nu_{1}^{(k)}=\sum_{\mathbf{l}\in J_{W}}y_{\mathbf{l}}^{(k)}$ and the estimates $\hat{\delta}^{(k)}_{1}=(\mu^{(k)}_{1},\Sigma^{(k)}_{1}),$ $\hat{\delta}^{(k)}_{2}=(\mu^{(k)}_{2},\Sigma^{(k)}_{2}),$ $\hat{\beta}^{(k)}$ by

[TABLE]

[TABLE]

[TABLE]

In the E-Step, we compute the updated probabilities $q_{\mathbf{l}}^{(k,SEM)}$ based on (30)-(32) by relation (29).

The essential idea of the SAEM algorithm is to mix $q_{\mathbf{l}}^{(k,EM)}$ and $q_{\mathbf{l}}^{(k,SEM)}$ in iteration step $k$ as

[TABLE]

which gives the a posteriori probabilities for the next $(k+1)$th iteration. Here in (33), $\{\lambda_{k},k\geq 1\}$ is a sequence of positive real numbers $\lambda_{k}\in(0,1)$ decreasing to zero. We stop the SAEM algorithm after the $k-$th iteration if $\sum_{\mathbf{l}\in J_{W}}|q_{\mathbf{l}}^{(k-1)}-q_{\mathbf{l}}^{(k)}|\leq\varepsilon,$ where $\varepsilon$ is some threshold. In the following, the choice of parameters is a result of experimental tuning to our image data yielding good practical results. Particularly, we use $\lambda_{k}=\frac{50}{50+k^{2}},k\geq 1$ and $\varepsilon=0.0001$ in our computations.

When the SAEM algorithm stops in the $k_{0}$-th iteration, we obtain the values $\{q^{(k_{0})}_{\mathbf{l}},\mathbf{l}\in J_{W}\}$ which indicate that $x_{\mathbf{l}}$ belongs to the first component if $q^{(k_{0})}_{\mathbf{l}}>1/2$ and to the second one, otherwise.

Applying the above SAEM algorithm to our image data yields diffuseness in the resulting clusters (see Figure 7). To avoid this, we propose a smoothing modification (Spatial SAEM), which takes the spatial location of the sample data into account. Let us describe the new Spatial step.

Let SAEM stop after $k_{0}$ iterations. For each sample entry $x_{\mathbf{l}},\mathbf{l}=(l_{1},l_{2},l_{3})\in J_{W}$ we define the coordinate $v_{\mathbf{l}}=\left((l_{1}-1)M\Delta,(l_{2}-1)M\Delta,(l_{3}-1)M\Delta\right),$ that is a vertex of the cube $W_{\mathbf{l}}.$ In each further iteration, i.e. for $k>k_{0}$, Bernoulli random variables $y^{(k)}_{\mathbf{l}}\in\{0,1\},\mathbf{l}\in J_{W}$ with success probability $q^{(k_{0})}_{\mathbf{l}},\mathbf{l}\in J_{W}$ are simulated. Now $y^{(k)}_{\mathbf{l}}$ classifies $x_{\mathbf{l}}$ in such a way that $y^{(k)}_{\mathbf{l}}=0$ indicates that $x_{\mathbf{l}}$ belongs to the first component and the second one, otherwise. Then, we compute the number $a^{(k)}_{\mathbf{l}}$ of neighbors of $x_{\mathbf{l}}$ belonging to the same cluster as $x_{\mathbf{l}}$. Neighborhood is defined in terms of the $r-$neighborhood of $v_{\mathbf{l}}$ such that

[TABLE]

$r>0.$

If all $a^{(k)}_{\mathbf{l}}$ are greater than or equal to a certain threshold $a$ (for our image data, $a=3$ is used) we call the classification $y^{(k)}_{\mathbf{l}},\mathbf{l}\in J_{W}$, admissible and move to the next iteration. Otherwise, for sample entries $x_{\mathbf{l}}$ with $a^{(k)}_{\mathbf{l}}$ being less than $a,$ we change $y^{(k)}_{\mathbf{l}}$, hence, the class of $x_{\mathbf{l}}$. If the new set $y^{(k)}_{\mathbf{l}},\mathbf{l}\in J_{W}$, is admissible, then we pass to iteration $k+1.$ If no, we resimulate $y^{(k)}_{\mathbf{l}},\mathbf{l}\in J_{W}$, until $a^{(k)}_{\mathbf{l}}\geq a$ for all $\mathbf{l}\in J_{W}.$

The smoothing procedure stops when $K$ admissible classifications have been generated ($K=1000$ is used). The final a posteriori probabilities $q_{\mathbf{l}}$ in the space of all admissible classifications are computed over the sample of $\{y^{(k)}_{\mathbf{l}},\mathbf{l}\in J_{W}\},k_{0}\leq k\leq k_{0}+K$ by

[TABLE]

We also get estimates of components’ weights $(\hat{\beta},1-\hat{\beta})$ by (28). If $\hat{\beta}\geq 0.5$ we say that the second component corresponds to the “anomaly”. The observation window $W_{\mathbf{l}}$ thus belongs to the zone of homogeneous material if $q_{\mathbf{l}}\geq 0.5$ and to the anomaly zone if $q_{\mathbf{l}}<0.5.$ If $\hat{\beta}<0.5$ the roles of the components are swapped.

6 Application to 3D image data of fibre materials

6.1 Simulated data

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. andra2014geometric ; Red11 . Figure 3 shows simulated RSA fibre data in an image of $2000\times 2000\times 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 Franke . In the two outer layers, the preferred direction is the $x$-direction and the concentration parameter is $\beta=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 ($\beta=0.5$).

Additionally, we investigate a homogeneous RSA data set where no anomalies should be detected. The data set consists of an image of $2000\times 2000\times 2100$ voxels. Here, the concentration parameter of the fibre direction distribution is $\beta=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_{\mathbf{k}},\mathbf{k}\in J$. In order to avoid cancelling effect of averaging, we build for each coordinate $x,y,z,$ the samples $\tilde{x}_{\mathbf{k}},\tilde{y}_{\mathbf{k}},\tilde{z}_{\mathbf{k}},$ such that their entries lie in the hemispheres $\mathbb{S}_{x}^{2}=\{(x,y,z)\in\mathbb{S},x\geq 0\},\mathbb{S}_{y}^{2}=\{(x,y,z)\in\mathbb{S},y\geq 0\},\mathbb{S}_{z}^{2}=\{(x,y,z)\in\mathbb{S},z\geq 0\},$ respectively, i.e., $\tilde{x}_{\mathbf{k}}=|x|_{\mathbf{k}},$ $\tilde{y}_{\mathbf{k}}=|y|_{\mathbf{k}},$ $\tilde{z}_{\mathbf{k}}=|z|_{\mathbf{k}},$

The sample of the estimated entropy values $\hat{E}_{\mathbf{k}}$ of directional distribution of fibres $X_{\mathbf{i}},\mathbf{i}\in J$ in the windows $W_{\mathbf{k}},\mathbf{k}\in J_{W}$ is build by estimator (7) over the transformed directions $\hat{X}_{\mathbf{i}}\in\mathbb{S}_{+}^{2},\mathbf{i}\in J.$

We apply the results of Section 4 consequently to the random fields $s_{\mathbf{k}}=$ $\tilde{x}_{\mathbf{k}},\tilde{y}_{\mathbf{k}},$ or $\tilde{z}_{\mathbf{k}},$ and $\hat{E}_{\mathbf{k}}.$ Therefore, we have the following 4 pairs of hypotheses of $(H_{0}^{\prime},H_{1}^{\prime})$-type.

• •

$H_{0}^{x}:$ $\mathbb{E}\tilde{x}_{\mathbf{k}}=\mu_{x}$ for every $\mathbf{k}\in W$ vs.

• •

$H_{1}^{x}:$ $\exists\theta_{0}\in\Theta_{0}$ such that $\mathbb{E}\tilde{x}_{\mathbf{k}}=\mu_{x}+h_{x},\mathbf{k}\in I_{\theta_{0}}$ and $\mathbb{E}\tilde{x}_{\mathbf{k}}=\mu_{x},\mathbf{k}\in I_{\theta^{c}_{0}},$ $h_{x}\neq 0$;

• •

$H_{0}^{y}:$ $\mathbb{E}\tilde{y}_{\mathbf{k}}=\mu_{y}$ for every $\mathbf{k}\in W$ vs.

• •

$H_{1}^{y}:$ $\exists\theta_{0}\in\Theta_{0}$ such that $\mathbb{E}\tilde{y}_{\mathbf{k}}=\mu_{y}+h_{y},\mathbf{k}\in I_{\theta_{0}}$ and $\mathbb{E}\tilde{y}_{\mathbf{k}}=\mu_{y},\mathbf{k}\in I_{\theta^{c}_{0}},$ $h_{y}\neq 0$;

• •

$H_{0}^{z}:$ $\mathbb{E}\tilde{z}_{\mathbf{k}}=\mu_{z}$ for every $\mathbf{k}\in W$ vs.

• •

$H_{1}^{z}:$ $\exists\theta_{0}\in\Theta_{0}$ such that $\mathbb{E}\tilde{z}_{\mathbf{k}}=\mu_{z}+h_{z},\mathbf{k}\in I_{\theta_{0}}$ and $\mathbb{E}\tilde{x}_{\mathbf{k}}=\mu_{z},\mathbf{k}\in I_{\theta^{c}_{0}},$ $h_{z}\neq 0$;

• •

$H_{0}^{E}:$ $\mathbb{E}\hat{E}_{\mathbf{k}}=\mu_{E}$ for every $\mathbf{k}\in J_{W}$ vs.

• •

$H_{1}^{E}:$ $\exists\theta_{0}\in\Theta_{0}$ such that $\mathbb{E}\hat{E}_{\mathbf{k}}=\mu_{E}+h_{E},\mathbf{k}\in I_{\theta_{0}}$ and $\mathbb{E}\hat{E}_{\mathbf{k}}=\mu_{E},\mathbf{k}\in I_{\theta^{c}_{0}},$ $h_{E}\neq 0.$

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 $\frac{1}{4}\alpha.$

Before running the algorithms we need to choose the right size of scanning windows. From the initial layered and homogeneous RSA images with $2000\times 2000\times 2100$ voxels we obtain $83\times 83\times 87$ small windows $\widetilde{W}_{\vec{i}}$ with $24\times 24\times 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 $\tilde{X}_{\mathbf{k}}$ to be $m-$dependent with $m=5$ and $\sigma^{2}=0.2,M_{0}=\frac{1}{2}.$ For mean local directions $\tilde{x},\tilde{y},$ and $\tilde{z}$, the parametric set $\Theta_{0}$ is constructed with $\Delta_{0}=\Delta_{1}=8,$ $\gamma_{0}=0.05,\gamma_{1}=0.5,$ and $L_{M}=22$ in (22), $|\Theta_{0}|=39395.$

We point out that the samples $(\tilde{x},\tilde{y},\tilde{z})$ and $\hat{E}$ have different sizes due to the construction described in Section 3. Therefore, the parameters $m,\sigma^{2}$ and the parameter set $\Theta_{0}$ in (22) for $\hat{E}$ differ from the ones for $\tilde{x},\tilde{y},\tilde{z}.$

For the sample of estimated entropy $\hat{E}_{\mathbf{k}},\mathbf{k}\in J_{W},$ we have $m=1,\sigma^{2}=0.5$ and the parametric set $\Theta_{0}$ is constructed with $\Delta_{0}=\Delta_{1}=2,$ $\gamma_{0}=0.05,\gamma_{1}=0.5,$ and $L_{M}=4$ in (22), $|\Theta_{0}|=16536.$ Entropy values approximately have a normal distribution, so we put $M_{0}=\sigma$ in (19).

The computed statistics (given in (11)) $T_{W}(\tilde{x}),T_{W}(\tilde{y}),$ $T_{W}(\tilde{z}),$ $T_{W}(\hat{E})$ 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 change-point 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.

It follows from (Penrose, , 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 of fibre locations and directions, the estimated local directional entropy $\widehat{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\sigma$-rule with $\sigma^{2}$ being the sample variance of $\widehat{E}$, compare ruiz2016entropy , to find anomaly regions in both Figures 3 and 4.

One can see that all centered entropy values lie in the interval $[-3\sigma;3\sigma],$ so the $3\sigma-$method does not distinguish between homogeneous (Figure 4) and inhomogeneous (Figure 3) images. We conclude that the $3\sigma$-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 $\widehat{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\times 5\times 5$ small windows $\widetilde{W}_{\vec{k}}$ was selected, i.e. $M=5.$ Additionally, we put $r=\Delta 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 everitetal11 ; ClustHandbook16 ; WKClus18 for an overview. In addition to the results reported here, we tried the recent AWC algorithm AWC . However, the spatial SAEM approach yields better results, cf. ConfMat . 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 $\mathbb{R}^{3}.$ Therefore, the estimates in the M-step (25)-(27), (30),(31) have to be modified, cf. Franke . This task goes, however, beyond the scope of the present article.

6.2 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 IVWFibres .

We apply the change point analysis from Section 4 to real data with $970\times 1469\times 1217$ voxels and the estimated radius of 3 voxels. We obtain $64\times 97\times 80$ small windows $\widetilde{W}_{\vec{i}}$ with $15\times 15\times 15$ voxels.

For mean local directions, the parametric set $\Theta_{0}$ is constructed with $\Delta_{0}=\Delta_{1}=8,$ $\gamma_{0}=0.05,\gamma_{1}=0.5$ and $L_{M}=22$ in (22), which gives $|\Theta_{0}|=33004.$ To choose the suitable value of $m$ for $m-$dependence we need additional investigation. Under the hypotheses $H_{0}^{x},H_{0}^{y},H_{0}^{z}$ the random fields $\tilde{x}$,$\tilde{y}$, and $\tilde{z}$ are assumed to be stationary, so that we can estimate their covariance functions. We use the standard approach and estimate e.g. $\rho_{x}(\mathbf{h})=\mathbf{Cov}(\tilde{x}_{\mathbf{1}},\tilde{x}_{\mathbf{1}+\mathbf{h}})$ as

[TABLE]

where $K=\{(k_{1},k_{2},k_{3})\in\mathbb{N}^{3},1\leq k_{1}\leq M_{1}-h_{1},1\leq k_{2}\leq M_{2}-h_{2},1\leq k_{3}\leq M_{3}-h_{3}\},$ and $\bar{x}_{0}$ and $\bar{x}_{\mathbf{h}}$ are the sample means of $\tilde{x}$ over the index ranges $K$ and $K+\mathbf{h}$ respectively. To visualize $\hat{\rho}_{x}$ we compute its maximum values in the following way:

[TABLE]

The values of $\hat{\rho}_{x,max}$(black color), $\hat{\rho}_{y,max}$(red color), $\hat{\rho}_{z,max}$(blue color) are given in Figure 13. We choose the value of $m$ in such a way that $\hat{\rho}_{x,max}(i)\leq\varepsilon_{0},$ $\hat{\rho}_{y,max}(i)\leq\varepsilon_{0},$ $\hat{\rho}_{z,max}(i)\leq\varepsilon_{0},$ for all $i\geq m,$ where $\varepsilon_{0}$ is a threshold. Here we use $\varepsilon_{0}=0.04$ and obtain $m=9.$ We have $M_{0}=0.5$ and assume that $\sigma^{2}=0.2.$

Moreover, due to simulation experiments in Section 4, we compute critical values for the change-point statistics $T_{W}(\cdot)$ 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 $\Theta_{0}$ is constructed with $\Delta_{0}=\Delta_{1}=2,$ $\gamma_{0}=0.05,\gamma_{1}=0.5$ and $L_{M}=4$ in (22), which gives $|\Theta_{0}|=12366.$ We assume that $\sigma^{2}=0.5$ and $M_{0}=\sigma.$ The result of our change point analysis is presented in Table 6.

Our change point test detects the evidence of anomaly regions in real fibre data at significance level $\alpha=8.4\times 10^{-10}.$

Similarly to the case of RSA data, the detection of anomalies by the $3\sigma-$rule gives meaningless results. The Spatial SAEM algorithm works much better. Its results are presented in Figures 14, 15, and 16, where the color labelling of points is the same as for the simulated data.

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.

Acknowledgements.

We are grateful to Dr. Stefanie Schwaar from Fraunhofer ITWM for valuable discussions and to Jan Niedermeyer for simulating RSA data.