Variational Multi-Phase Segmentation using High-Dimensional Local Features

TL;DR

This paper introduces a versatile multi-phase image segmentation method leveraging high-dimensional local features, effective in noisy crystal images and texture segmentation, using PCA and convex optimization techniques.

Contribution

It presents a novel framework that employs high-dimensional local features and PCA for robust multi-phase segmentation, adaptable beyond crystal images.

Findings

Competitive performance on Prague benchmark

Robust to extreme noise in crystal images

No prior knowledge of crystal structure needed

Abstract

We propose a novel method for multi-phase segmentation of images based on high-dimensional local feature vectors. While the method was developed for the segmentation of extremely noisy crystal images based on localized Fourier transforms, the resulting framework is not tied to specific feature descriptors. For instance, using local spectral histograms as features, it allows for robust texture segmentation. The segmentation itself is based on the multi-phase Mumford-Shah model. Initializing the high-dimensional mean features directly is computationally too demanding and ill-posed in practice. This is resolved by projecting the features onto a low-dimensional space using principle component analysis. The resulting objective functional is minimized using a convexification and the Chambolle-Pock algorithm. Numerical results are presented, illustrating that the algorithm is very competitive in texture segmentation with state-of-the-art performance on the Prague benchmark and provides new possibilities in crystal segmentation, being robust to extreme noise and requiring no prior knowledge of the crystal structure.

Tables

Table 1. Table 1: Gray-scale texture segmentation comparison on the Outex_US_00000 test suite with known number of segments. FSEG  [ 33 ] uses the ICPR2014 code of FSEG with fixed number of segments ( segn = 5 segn 5 \mathrm{segn=5} ). FSEG ∗ and Algorithm  2 both combine spectral histograms with 11 11 11 bins, window sizes s = 15 𝑠 15 s=15 and s = 30 𝑠 30 s=30 (stacked with weights 0.8 and 0.2 respectively), and Gabor filters of kernel sizes σ = 5 , 7 , 9 𝜎 5 7 9 \sigma=5,7,9 and orientations θ = 0 , 1 2 ​ π , 1 4 ​ π , − 1 4 ​ π 𝜃 0 1 2 𝜋 1 4 𝜋 1 4 𝜋 \theta=0,\frac{1}{2}\pi,\frac{1}{4}\pi,-\frac{1}{4}\pi . FSEG ∗ - TxtMerge TxtMerge \mathrm{TxtMerge} is the same as FSEG ∗ but runs without FSEG’s TxtMerge TxtMerge \mathrm{TxtMerge} post-processing. Bold face highlights the best, a star the second-best result in each column.

Method	CS $\uparrow$	O $\downarrow$	C $\downarrow$	CA $\uparrow$	CO $\uparrow$	CC $\uparrow$
FSEG∗-$\mathrm{TxtMerge}$	85.20 $\pm$ 22.7	*4.76 $\pm$ 3.9	6.05 $\pm$ 5.5	*82.49 $\pm$ 13.4	*89.23 $\pm$ 9.6	*88.20 $\pm$ 11.8
Algorithm 2	*80.60 $\pm$ 28.3	5.82 $\pm$ 8.3	5.62 $\pm$ 6.8	82.73 $\pm$ 15.5	89.24 $\pm$ 11.3	88.84 $\pm$ 13.3
FSEG∗	68.60 $\pm$ 24.3	3.61 $\pm$ 3.2	*6.04 $\pm$ 8.8	70.67 $\pm$ 17.8	80.40 $\pm$ 13.1	73.51 $\pm$ 18.6
Clustering	60.00 $\pm$ 26.8	13.36 $\pm$ 8.1	15.52 $\pm$ 12.0	66.41 $\pm$ 12.5	77.54 $\pm$ 10.0	77.51 $\pm$ 11.2
FSEG [33]	45.80 $\pm$ 26.0	19.29 $\pm$ 34.2	17.55 $\pm$ 19.1	50.65 $\pm$ 20.6	64.95 $\pm$ 16.3	52.88 $\pm$ 21.9

Table 2. Table 2: Color texture segmentation on the Prague ICPR2014 contest dataset with unknown number of segments. Bold face highlights the best, a star the second-best value in each column, and † indicates that no corresponding publication could be found at the time of writing.

Method	CS $\uparrow$	OS $\downarrow$	US $\downarrow$	ME $\downarrow$	NE $\downarrow$	O $\downarrow$	C $\downarrow$	CA $\uparrow$	CO $\uparrow$	CC $\uparrow$
VRA-PMCFA${}^{\text{†}}$	75.32	*11.95	*9.65	4.57	*4.63	4.51	8.87	83.50	88.16	90.73
Algorithm 2$+\mathrm{TxtMerge}$	*72.27	18.33	9.41	4.19	3.92	*7.25	6.44	*81.13	*85.96	91.24
FSEG [33]	69.18	14.69	13.64	5.13	*4.63	9.25	12.55	78.22	84.44	87.38
SegTexCol${}^{\text{†}}$	61.19	1.92	27.02	9.33	9.05	15.17	12.12	71.69	81.16	76.34
MW3AR8 [16]	53.66	51.40	14.21	*5.54	6.33	19.86	84.27	70.15	75.41	89.36
RS${}^{\text{†}}$	46.02	13.96	30.01	12.01	11.77	35.11	29.91	58.75	68.89	69.30
Deep Brain [13]	36.20	41.87	53.87	7.38	9.06	47.53	99.56	49.97	62.62	70.08

Equations

$$\underset{({\Omega }_l{)}_{l=1}^k}{\min }\sum _{l=1}^k\left\{{\int }_{{\Omega }_l}{f}_l\mathrm{d}x+\lambda \text{Per}({\Omega }_l,\Omega )\right\}.$$

$$E[u]:=\sum _{l=1}^k\left\{{\int }_{\Omega }{f}_l{u}_l\mathrm{d}x+\lambda |u_l{|}_{\text{TV}(\Omega )}\right\}$$

$$\mathcal{U}:=\left\{u\in \text{BV}(\Omega {)}^k:u\ge 0\wedge \sum _{l=1}^k{u}_l=1\text{ a.e. in }\Omega \right\}$$

$$|u{|}_{\text{TV}(\Omega )}:=\underset{\genfrac{}{}{0ex}{}{p\in C_c^1(\Omega ,{\mathbb{R}}^2)}{\Vert p{\Vert }_{\infty }\le 1}}{\sup }{\int }_{\Omega }u\mathrm{div}p$$

$$\underset{\mathbf{u}\in {\mathcal{U}}_h}{\min }\underset{\genfrac{}{}{0ex}{}{{\hat{p}}_{li}}{\Vert {\hat{p}}_{li}\Vert \le 1}}{\max }\sum _{l=1}^k\sum _{i=1}^n\left\{f_{li}{u}_{li}+\lambda ⟨(K{\mathbf{u}}_l{)}_i,{\hat{p}}_{li}⟩\right\},$$

$$R_1(\mathbf{p})$$

$$R_2(\mathbf{u})$$

$$\mathcal{F}[\mathbf{g}](x_i):=\mathcal{F}(\mathbf{g}(W_s(x_i))).$$

$$\underset{l=1,\dots ,k}{\max }\underset{x,x^{\prime }\in {\Omega }_l^{\ast }}{\max }\Vert \mathcal{F}[\mathbf{g}](x)-\mathcal{F}[\mathbf{g}](x^{\prime })\Vert$$

$$\underset{\begin{array}{c}l,l^{\prime }=1,\dots ,k\\ l\ne l^{\prime }\end{array}}{\min }\underset{x\in {\Omega }_l^{\ast },x^{\prime }\in {\Omega }_{l^{\prime }}^{\ast }}{\min }\Vert \mathcal{F}[\mathbf{g}](x)-\mathcal{F}[\mathbf{g}](x^{\prime })\Vert$$

$$\boldsymbol{c}=\bigg{(}\frac{1}{|\Omega_{l}^{*}|}\sum_{x\in\Omega_{l}^{*}}\mathcal{F}[\boldsymbol{g}](x)\bigg{)}_{l=1}^{k}\in\mathbb{R}^{m\times k},$$

$$f_{li}:=\Vert \mathcal{F}[\mathbf{g}](x_i)-{\mathbf{c}}_l{\Vert }^2.$$

$${\mathbf{c}}_l^{\{t+1\}}=\sum _{i=1}^n\mathcal{F}[\mathbf{g}](x_i)u_{li}^{\{t\}}/\sum _{i=1}^n{u}_{li}^{\{t\}}.$$

$$\Vert \mathcal{F}[\mathbf{g}](x_i)-{\mathbf{c}}_l{\Vert }^2=\Vert {\mathbf{\alpha }}_i^{(m)}-{\mathbf{\gamma }}_l^{(m)}{\Vert }^2$$

$$U^{(m)}{\mathbf{\gamma }}_l^{(m)}={\mathbf{c}}_l-\mathbf{\mu }[\mathbf{g}]=U^{(m)}\frac{{\sum }_{i=1}^n{\mathbf{\alpha }}_i^{(m)}{u}_{il}}{{\sum }_{i=1}^n{u}_{il}}$$

$$\big{|}f_{li}-\|\boldsymbol{\alpha}^{(r)}_{i}-\boldsymbol{\gamma}^{(r)}_{l}\|^{2}\big{|}\leq 2\textstyle\sum_{j=r+1}^{m}\lambda_{j}$$

$${\sum }_{j=r+1}^m{\lambda }_j=\Vert A{\Vert }_F^2-{\sum }_{j=1}^r{\lambda }_j.$$

$${\left({\mathcal{F}}_{\text{SH}}[\mathbf{g}](x)\right)}_{ij}=\frac{\sum _{x_k\in W_s(x)}{\int }_{z_{ij}}^{z_{i,j+1}}\delta (z-g_{ik}))\mathrm{d}z}{|W_s(x)|}$$

$$L_{a_1,a_2}=\{n_1{a}_1+n_2{a}_2:n_1,n_2\in \mathbb{Z}\},$$

$$f_l(x):=\frac{1}{N}{\sum }_{k=1}^{N}d(g(x),g(x+M({\alpha }_l)q_k)),$$

$${\mathcal{F}}_{\text{FFT}}[\mathbf{g}](x)=(|\text{FFT}[\mathbf{g}(W_s(x)){]}_{ij}|{)}_{ij}.$$

Full text

1 Introduction

Image segmentation, i.e. the task of decomposing an image into disjoint regions that are roughly homogeneous in a suitable sense, is one of the fundamental image processing problems. If three or more regions are sought, one speaks of multi-phase segmentation. This problem has been studied thoroughly in the literature and entirely different concepts have been put forward as the basis for image segmentation, such as fuzzy region competition [21], contour detection [2], random walks [9], markov random fields [22], just to name a few. Due to the variety of proposed methods, providing a comprehensive list is beyond the scope of this article, but we refer the interested reader to [31]. Then there is, of course, the class of variational approaches based on the famous Mumford-Shah energy [25].

The most straight-forward application in multi-phase segmentation is to divide images into regions based on their gray or color intensities [8]. A more complex task is to segment images based on their local structure. This has applications in texture segmentation [28], as well as many medical applications, such as the segmentation of blood vessels [14]. Algorithms for structure classification and segmentation usually extract local features from the image, which analyze important properties of the structures of interest, such as the image intensity, position and orientation of edges, or the local frequency spectrum [29]. In the case of texture segmentation, Gabor filters are arguably the most popular source of feature discrimination [32], often combined with other filters in so-called local spectral histograms [23]. Other methods rely on linear transforms, such as the short-time Fourier transform [3], wavelet transforms [7], or, more recently, the Stockwell transform [11]. While the part of this paper on texture segmentation uses well-proven spectral histograms to recognize regions, it differs from established methods by their integration into a variational framework, allowing to control the regions’ connectedness.

Dealing with complex structures, such as textures, often implies high-dimensionality of the parameters describing the problem. However, in image segmentation, one is mostly interested in classifying structure into a few categories, potentially allowing for lower-dimensional representations. Dimension reduction of high-dimensional data is an immensely broad topic [30] and finds applications in many different areas of research. There exist different approaches, but the most widely used techniques are arguably clustering [27] and principal component analysis (PCA) [19]. The latter two are connected in the sense that the relaxed solution of $k$-means, one of the most popular clustering algorithms, is given by principal components [10]. PCA has been investigated in the context of variational image segmentation before, both as a means for dimension reduction [26] and to increase the contrast of color-texture indicators in natural images [17].

In materials science, an important application of structure-based segmentation is the analysis of crystals. Available methods are based on variational minimization of Mumford-Shah energies that require the local stencil of a reference crystal as prior knowledge [4, 5, 12].

1.1 Key contributions

• •

a widely applicable framework for image segmentation by structure is discussed, including a novel combination of PCA of high-dimensional features, Mumford-Shah and a robust initialization strategy, which allows for a broad choice of feature descriptors

• •

the framework is shown to work very well, even for extremely noisy data, in crystal segmentation, where it generalizes existing methods in the sense that no a-priori information about the crystals is required

2 Variational multi-phase segmentation

In this section, we briefly recall the Mumford-Shah model [25] for multi-phase segmentation based on suitable indicator functions. Furthermore, we recall a convexification approach that enables an efficient numerical minimization of the model. Let $\Omega=[0,1]^{2}$. The task is to divide $\Omega$ into pairwise disjoint regions $\Omega_{l}$, $l=1,\dots,k$ based on given indicator functions $f_{1},\dots,f_{k}:\Omega\rightarrow\mathbb{R}_{\geq 0}$. $f_{l}(x)$ can be interpreted as the cost of putting a point $x\in\Omega$ into the set $\Omega_{l}$. For instance, if an image $g:\Omega\rightarrow\mathbb{R}$ is supposed to be segmented based on its gray values, possible indicator functions are $f_{l}(x)=(g(x)-c_{l})^{2}$. Here, $c_{l}$ is the average gray value of $g$ in the $l$-th region.

A segmentation of $\Omega$ based on the indicator functions that guarantees a certain regularity of the segments can be achieved by minimizing the Mumford-Shah energy [25]:

[TABLE]

Here, $\text{Per}(\Omega_{l},\Omega)$ denotes the perimeter of the set $\Omega_{l}$ in $\Omega$ [1]. Roughly speaking, the perimeter is the length of the boundary of $\Omega_{l}$, not counting the parts of the boundary that are also on the boundary of $\Omega$. This problem is hard to address numerically since the unknown variables are sets. In particular, its discrete counterpart, known as Pott’s model, is NP-hard. Thus, various convex relaxation approaches have been proposed in the past. For the sake of simplicity, we use one of the most straightforward approaches, given in [34]. Let us stress that our framework does not rely on this particular choice, but can also be combined with more sophisticated convexification approaches. Let

[TABLE]

where $u\in\mathcal{U}$ is a vector valued labeling function and

[TABLE]

is the admissible set. Here,

[TABLE]

denotes the total variation and $\text{BV}(\Omega)$ is the space of functions of bounded variation, i.e. the space of Lebesgue integrable functions with finite total variation. Then, the convex relaxation of (1) is to minimize $E$ over the set $\mathcal{U}$. The minimizer $u^{*}$ can be interpreted as a soft segmentation and can be converted into a hard segmentation by setting $\Omega_{l}:=\{x\in\Omega:u^{*}_{l}(x)\geq u^{*}_{j}(x)\,\forall j\neq l\}$.

In order to address this minimization numerically, a discretization of the energy (2) and the admissible set (3) is required. Let $X=(x_{i})_{i=1}^{n}\in\mathbb{R}^{2\times n}$ be a regular 2D pixel grid. We use piecewise constant approximations $f_{li}=f_{l}(x_{i})$ and $u_{li}=u_{l}(x_{i})$ for $i=1,\dots,n$. The corresponding column vectors and matrices of all pixel values are denoted by a boldface letter, e.g. $\boldsymbol{u}=(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{k})\in\mathbb{R}^{n\times k}$. Furthermore, we denote with $K:\mathbb{R}^{n}\rightarrow\mathbb{R}^{2\times n}$ the discrete gradient operator corresponding to the grid $X$ and forward differences. Using this operator to discretize the total variation (4), the minimization of the discretized energy (2) can be posed as the following discrete saddle point problem:

[TABLE]

where $\mathcal{U}_{h}$ is the discrete counterpart of $\mathcal{U}$ and $\hat{p}_{li}=(p_{li1},p_{li2})$ discretizes $p$ from (4) for $u_{l}$ at node $x_{i}$. Problems of this form can be solved with the Chambolle-Pock algorithm [6], summarized in Algorithm 1.

The required resolvent operators are given by

[TABLE]

Here, $\boldsymbol{f}=(\boldsymbol{f}_{1},\ldots,\boldsymbol{f}_{k})$ and $\pi_{\mathcal{U}_{h}}(\boldsymbol{u})$ denotes the orthogonal projection of $\boldsymbol{u}$ onto the set $\mathcal{U}_{h}$. This projection can be calculated with $\mathcal{O}(k)$ operations using an iterative algorithm described in [24]. In this work, all numerical experiments use the parameters $\sigma=\tau=\frac{1}{8},\theta=0.7,\epsilon=0.001,t_{\text{max}}=10000$. The regularization parameter is chosen as $\lambda=0.01$ (Table 2), $\lambda=0.005$ (Table 1, Figure 1) and $\lambda=25$ (Figure 2).

3 Local features for structure characterization

3.1 Description and relation to Mumford-Shah

Our aim is to provide a method to segment images into regions of different structure based on the information from local features. In the discrete setting, local features corresponding to a pixel $x_{i}$ are encoded in the values of an input image $\boldsymbol{g}$ in a $(2s+1)\times(2s+1)$ window $W_{s}(x_{i})$ centered at $x_{i}$. Here, $s\in\mathbb{N}$ determines the scale that is still considered to be local. From these values, features are extracted by an operator of the form $\mathcal{F}:\mathbb{R}^{(2s+1)^{2}}\rightarrow\mathbb{R}^{m}$, that should fulfill certain properties, which we will detail later. Applying $\mathcal{F}$ to the $(2s+1)\times(2s+1)$ matrix $\boldsymbol{g}(W_{s}(x_{i}))$ containing the image pixel values in the window $W_{s}(x_{i})$ gives the feature vector corresponding to the pixel $x_{i}$:

[TABLE]

Let $\Omega_{l}^{*}\subset X$, $l=1,\dots,k$ denote the sought discrete regions, i.e. the true sets of pixels belonging to the different structure regions. Then, a suitable feature extractor (as defined in (8)) for discriminating regions of different structures can be characterized by the following two properties:

[TABLE]

i.e. local features should vary as little as possible within each region and offer as much contrast as possible between different regions. Examples for robust feature extractors for texture and crystal segmentation will be discussed in Section 4.

Given a suitable feature extractor and the true mean features within the different structure regions

[TABLE]

the following indicator can be used for segmentation in (5):

[TABLE]

In practice, the mean values $\boldsymbol{c}$ are of course unknown. However, given some approximate guess $\boldsymbol{c}^{\{t\}}$ for the mean values, Algorithm 1 can be applied, resulting in a segmentation $\boldsymbol{u}^{\{t\}}$. Then, the following update rule can be applied to refine the mean features

[TABLE]

This way, given some initial guess $\boldsymbol{c}^{\{0\}}$, both the segmentation and the mean features can be refined in an alternating fashion. Note that we use curly brackets instead of round ones for the index here, to differentiate between the iterations within Algorithm 1 and these outer iterations. Unfortunately, the result of this alternating minimization strategy depends heavily on the initial guess $\boldsymbol{c}^{\{0\}}$. In the literature, it is often suggested to approximate $\boldsymbol{c}^{\{0\}}$ via clustering, which is equivalent to minimizing (1) for $\lambda=0$ with (12) as indicator and with respect to both the regions $\Omega_{l}$ and $\boldsymbol{c}$. This clustering problem is NP-hard itself, but efficient iterative solvers, such as $k$-means, are available and have proven to work well in the case of low-dimensional indicator functions (e.g. in color segmentation) [8]. However, robust feature extractors suitable for structure discrimination tend to be high-dimensional ($m$ greater than $100$ or even $1000$). In this case, clustering becomes unfeasible in practice, because the available solvers are likely to get stuck in undesired local minima when applied in such high dimensions.

3.2 Dimension reduction and decorrelation

Clustering of high-dimensional data is a well studied problem in the literature [27]. It has been noted that often many of the dimensions are irrelevant for the core information expressed by a given data set and that they might mask the essential clusters due to noise. Therefore, several approaches for subspace clustering have been proposed to address this problem [20]. In our context, dimension reduction and decorrelation via principal component analysis (PCA) should work well: given a feature extractor $\mathcal{F}$ fulfilling (9) & (10), $\left(\mathcal{F}[\boldsymbol{g}](x_{i})\right)_{x\in\Omega_{l}^{*}}$ is of low variance for any $l\in\{1,\dots,k\}$ and, compared to this, for $l\neq l^{\prime}$ the set $\left(\mathcal{F}[\boldsymbol{g}](x),\mathcal{F}[\boldsymbol{g}](x^{\prime})\right)_{x\in\Omega_{l}^{*},x^{\prime}\in\Omega_{l^{\prime}}^{*}}$ is of high variance.

Performing PCA on the matrix of mean-centralized features $A=(\mathcal{F}[\boldsymbol{g}](x_{1})-\boldsymbol{\mu}[\boldsymbol{g}],\dots,\mathcal{F}[\boldsymbol{g}](x_{n})-\boldsymbol{\mu}[\boldsymbol{g}])\in\mathbb{R}^{m\times n}$ with $\boldsymbol{\mu}[\boldsymbol{g}]=\frac{1}{n}\sum_{i=1}^{n}\mathcal{F}[\boldsymbol{g}](x_{i})$, results in a lower-dimensional coefficient representation $\boldsymbol{\alpha}^{(r)}=(U^{(r)})^{T}A\in\mathbb{R}^{r\times n}$, where $U^{(r)}\in\mathbb{R}^{m\times r}$ is the matrix of eigenvectors belonging to the largest $r$ eigenvalues of $AA^{T}$. Clustering the coefficients $\boldsymbol{\alpha}^{(r)}$ into $k$ clusters gives a coefficient representation $\boldsymbol{\gamma}^{(r)}\in\mathbb{R}^{r\times k}$, which results in the initial guess $\boldsymbol{c}^{\{0\}}=U^{(r)}\boldsymbol{\gamma}^{(r)}+\boldsymbol{\mu}[\boldsymbol{g}]$. Since we need $\boldsymbol{c}\in\mathbb{R}^{m\times k}$, $r=k$ is a natural choice for dimension reduction.

In [33], it was noted that the clustering can get stuck in local minima due to effects caused by the inhomogeneity of the features across the boundary between two regions. Unlike purely point-wise indicators ($s=0$), local feature extractors cause points within about half the window size of a region boundary (in 2D space) to spread between the two mean features corresponding to the regions adjacent to the boundary (in coefficient space). In order to prevent the $k-$means minimizer from getting stuck in-between such two clusters, Yuan et al. proposed to disregard such boundary points when clustering by thresholding an edgeness indicator, given by finite differences of the features on the scale of the window size [33]. As this approach is only based on the assumption of homogeneity of the features within each structure region, it can be used for general feature extractors and allows us to adopt this technique for the initial clustering.

While the above resembles a robust method to retrieve an initial value for $\boldsymbol{c}$ in the full dimensional feature space, the dimension reduction we now have at hand also suggests itself to reduce the noise of the high dimensional feature vectors and to increase their inter-region contrast within the subsequent variational segmentation framework. First of all, let us point out that as $U^{(m)}$ forms an orthonormal basis of $\mathbb{R}^{m}$, we can express the indicator (12) and thus the fidelity term in (5) in terms of $\boldsymbol{\gamma}^{(m)}$ and $\boldsymbol{\alpha}^{(m)}$:

[TABLE]

Furthermore, definition (13) can be rewritten as

[TABLE]

i.e. the mean values $\boldsymbol{\gamma}$ can be updated using the coefficients $\boldsymbol{\alpha}$ instead of the feature vectors $\mathcal{F}[\boldsymbol{g}](x_{i})$. Reducing the dimension to $r<m$ introduces an error, which can be bounded by the eigenvalues $\lambda_{1}\geq\dots\geq\lambda_{m}$ of $AA^{T}$:

[TABLE]

This inequality can be deduced by applying the triangle inequality to the difference of the left- and right-hand side in (14), splitting $\boldsymbol{\alpha}^{(m)}_{i},\boldsymbol{\gamma}^{(m)}_{l}$ into $\boldsymbol{\alpha}^{(r)}_{i},\boldsymbol{\gamma}^{(r)}_{l}$ and remaining parts $\bar{\boldsymbol{\alpha}}_{i},\bar{\boldsymbol{\gamma}}_{l}$, as well as representing $\bar{\boldsymbol{\alpha}}_{i},\bar{\boldsymbol{\gamma}}_{l}$ as convex combinations of the columns of $(U^{(m)})^{T}A$ with non-zero coefficients in columns corresponding to $\lambda_{r+1},\dots,\lambda_{m}$.

Note that the error in (16) can be estimated without calculating all $m$ eigenvalues, which may become computationally expensive when the dimension $m$ becomes large:

[TABLE]

Here, $\|A\|_{F}=\sqrt{\sum_{i=1}^{m}\sum_{j=1}^{n}|A_{ij}|^{2}}$ denotes the Frobenius norm. This way, the error in the fidelity term can still be monitored, when the eigenvectors of $AA^{T}$ are calculated iteratively, e.g. using a deflation type of strategy. We use $r=k$ within the entire framework. Algorithm 2 summarizes the proposed method. All numerical experiments use $\hat{t}_{\text{max}}=3$. The edgeness threshold on the finite differences before clustering in Algorithm 2 is chosen as $\delta=0.5$ (Table 2), $\delta=0.25$ (Table 1, Figure 1) and $\delta=1.0$ (Figure 2).

3.3 Properties and advantages of the method

PCA has been used, for instance, as a concept for dimension reduction of PET data and subsequent variational segmentation [26], as well as a tool for increasing the contrast in the region descriptors for natural color-texture images in a variational segmentation framework [17]. Moreover, Yuan et al. [33] utilized the related concept of singular value decomposition to compute a low-rank factorization of a local spectral histogram based feature matrix and estimate subsequent template features via clustering. We want to stress that the initialization step in the proposed method shares the idea of dimension reduction and clustering of features, albeit differing slightly in the details, and is, in this regard, similar to [33]. However, our work embeds these ideas into a variational segmentation framework, which grants the following two main advantages:

First, the proposed method can be applied to a very general class of feature extractors, since it only relies on the natural properties (9), (10). In particular, in contrast to [33], it does not rely on the assumption that the feature vectors are linear combinations of the mean features in each region (this assumption and its consequences are discussed later in this section). Among others, the generality of our framework allows the usage of globally coupling, convolution based linear transforms. Functions of this type, such as the short-time Fourier transform [3], the Stockwell-Transform [11], or different types of wavelet transforms [7], have been studied for texture segmentation and shown their performance.

Second, the dimension reduction of the fidelity term helps to increase the degree to which it fulfills (9), (10). In particular, incorporating the PCA not only in the initial clustering, but also throughout the entire variational minimization, helps to suppress noise in the fidelity term. In Section 4.2, we will demonstrate how effective this strategy performs in the case of extremely noisy crystal images, using the Fourier transform as the feature extractor.

Unlike [23, 33], the general applicability of the proposed framework is tied to the need for a regularization of the segment boundaries, which is covered by the Mumford-Shah model. This need arises from an unexpected behavior of the indicator functions near segment boundaries. Due to the window size, the feature extractor sees a mixture of different segments near boundaries. For general feature extractors, this means that feature vectors near boundaries are not necessarily a linear combination of the cluster centers corresponding to the adjacent segments. In case $k>2$, it might happen that the feature vector at a boundary between two regions $\Omega_{1}^{*},\Omega_{2}^{*}$ is nearer to the mean feature vector of a third region $\Omega_{3}^{*}$ than it is to $\Omega_{1}^{*}$ or $\Omega_{2}^{*}$ itself. This means that the indicator $\boldsymbol{f}$ cannot necessarily identify the correct segment within a distance of $s$ to the sought segments. Note that this effect does not arise for $k=2$. As mentioned above, the perimeter regularization within the Mumford-Shah model addresses this problem for practical purposes. For input data where the regularization alone is not sufficient, we suggest to combine feature extractors of different window sizes.

Beyond this, the proposed method is an extension of [23, 33] in the sense that the decoupling of the coefficient representation from the segmentation allows for an outer iteration to refine the mean features, whereas in [33] the clusters are solely computed from the feature matrix.

Let us point out that, since the method is based on local windows, it has the common limitations inherent to such methods. The feature scale is tied to the window size $s$, so the method can only reliably detect regions that are at least somewhat larger than the window $W_{s}(x)$. Furthermore, special care has to be taken close to the boundary, where the window $W_{s}(x)$ leaves the support of the image. Please note that the proposed method enforces the region boundaries to approach the image boundary orthogonally, which is due to the natural boundary conditions in the Euler-Lagrange equation of (2). This effect can be reduced by introducing ghost cells at the image boundary with a zero extension of all indicators, but it is still noticeable (cf. Figure 1).

4 Applications and numerical results

4.1 Texture segmentation

Apart from plain gray value or color intensities, among the most thoroughly studied types of structures in image segmentation are textures [18]. In the image processing sense, a texture essentially consists in a more or less strictly repetitive pattern of the spatial arrangement of the gray or color values in an image. Thus, indicators for texture segmentation need to take into account image information from a whole neighborhood, at least on the scale of the spatial distance between repetitions. There are two main classes of operators that have been proposed in the literature, namely 1) local spectral histograms combined with a suitable bank of filters and 2) localized linear transforms, both of which fall into the class of feature extractors described earlier. In the context of texture segmentation, we limit our analysis to the first class, while the second class will be utilized for crystal segmentation in the next section.

Local spectral histograms are defined as follows: first a bank of $p$ filters is selected and applied to the image, resulting in a sequence of filtered images $\boldsymbol{g}_{1},\dots,\boldsymbol{g}_{p}$. Then, the feature extractor is defined by

[TABLE]

Here, $z_{ij}$, $i=1,\dots,p$, $j=1,\dots,q+1$ define the binning of the histograms and are often chosen such that $z_{i,1}=\min_{x}g_{i}(x),z_{i,q}=\max_{x}g_{i}(x)$ and equidistant in between. Thus, the dimension of the extracted feature at every pixel $x$ is $m=pq$. The most popular filter used in this context is arguably the Gabor filter. Other commonly used filters are Gaussian filters, Laplacian of Gaussian filters, or just the intensity filter (i.e. the identity). For a thorough description of spectral histograms of filtered images and their application to texture segmentation, we refer to [23].

In the following, we compare our approach to the Outex_US_00000 test suite of the Outex texture database (http://www.outex.oulu.fi) and the Prague ICPR2014 contest [15] (http://mosaic.utia.cas.cz/icpr2014/).

On the Outex database, we provide a thorough comparison between the proposed method and Factorization-Based Texture Segmentation (FSEG) [33]. We chose to focus on FSEG here, because it 1) ranks best in the ICPR2014 contest among methods with available code and 2) is similar to our framework. Table 1 quantifies the mean segmentation performance and its standard deviation over all 100 texture mosaics from the Outex_US_00000 test suite. Three different versions of FSEG, described in the caption of Table 1, are used for this comparison. Note that running FSEG without its $\mathrm{TxtMerge}$ post-processing (first row) makes sense in this case, since the number of segments is known. However, this also disables filling holes in the segments, as seen in the sixth column of Figure 1. While FSEG∗-$\mathrm{TxtMerge}$ performs best in correct segmentation (CS), and FSEG∗ achieves the smallest omission error (O), our method ranks highest in all remaining measures. Note that in Table 1 we also compare our method to plain clustering in order to evaluate the benefit of 1) the improved initialization strategy via PCA and 2) the subsequent variational optimization including region boundary smoothing. Indeed, the proposed method performs significantly better than plain clustering. Finally, a visual inspection of Figure 1 indicates that the proposed method provides a good compromise between fidelity of region boundaries and reduction of artifacts (holes, missing regions).

Next, we compare the proposed method to results from the Prague ICPR2014 contest [15]. Here, we use the same feature extractors as above (on the lightness channel), except that the kernel size $\sigma=9$ is omitted. Therefore, we add an intensity filter on all three channels (in Lab color space) to each spectral histogram. The number of segments is estimated as $k=\min\{k^{\prime}:\frac{1}{n}\sum_{j=k^{\prime}+1}^{m}\lambda_{j}<\omega\}$ with $\omega=0.05$. Since this estimate is not precise (even for an optimal choice of $\omega$), we additionally employ FSEG’s $\mathrm{TxtMerge}$ post-processing. Table 2 quantifies the mean segmentation quality over all 80 colored texture mosaics from the Prague ICPR2014 contest dataset. While our method produces larger over-segmentation (OS) than the other best-ranked methods, indicating stronger over-estimation of the number of segments, it performs best for under-segmentation (US), indicating a good coverage of all ground truth segments, reflecting the good initialization strategy. Moreover, our method performs second-best for correct segmentation (CS), omission error (O), class accuracy (CA) and correct assignment (CO). Most notably, according to all other presented measures and in total half of them, our method performs best among all competitors.

Note that VRA-PMCFA resolves fine boundary features but produces labeling noise, whereas our method smoothes region boundaries in favor of suppressing labeling noise. Thus, it depends on the application which of the two methods is likely to be more suitable.

4.2 Unsupervised crystal segmentation

A fundamental research topic in materials science is the analysis and modeling of crystals. Modern transmission electron microscopes (TEM) allow for imaging at atomic scale, which makes the crystal grid visible. In a perfect setting, the crystal is given by a Bravais lattice

[TABLE]

where $a_{1},a_{2}\in\mathbb{R}^{2}$ are the lattice vectors defining the orientation and spacing of the crystal. However, crystals of interest usually exhibit a more complicated behavior, like discontinuous orientation changes along so-called grain boundaries. The fully automatic analysis of grain geometries in TEM images is subject of ongoing research [12].

Available variational approaches for grain segmentation [4, 5, 12] are built on the assumption that all grains can be characterized through transformations of a local stencil $q_{1},\dots,q_{N}\in\mathbb{R}^{2}$, corresponding to a reference crystal given by all linear combinations of $a_{1},a_{2}$ with coefficients in $\{-1,0,1\}$. Then, the Mumford-Shah model with an indicator function of the following type can be used [4, 5]:

[TABLE]

Here, $d$ denotes a suitable intensity distance function and $M(\alpha_{l})\in\mathbb{R}^{2\times 2}$ is an orthogonal matrix, rotating the stencil by the angle of the the $l$-th grain relative to the reference.

The need for a-priori knowledge of the reference crystal structure inherent to indicators like (20) is a severe limitation of available methods. As we will show, this limitation can be overcome by using our proposed framework with the modulus of the 2D-FFT as the local feature extractor:

[TABLE]

Let us assume that the window $W_{s}(x)$ is large enough to cover at least one period of the crystal in either direction at any location. Then, the modulus of the Fourier transform $\mathcal{F}_{\text{FFT}}[\boldsymbol{g}](x)$ automatically encodes the local stencil $\left(M(\alpha_{l})q_{k}\right)_{k=1}^{N}$ within the positions of Bragg reflections.

Assuming that the window size $s$ matches the period of the crystal and the unit cell is a square, i.e. the discrete signals $\boldsymbol{g}(W_{s}(x))$ are exactly periodic, the translation of the window across the image causes a phase shift in frequency domain. This phase shift is canceled by the absolute value in the modulus, making the feature extractor $\mathcal{F}_{\text{FFT}}[\boldsymbol{g}](x)$ translation invariant inside crystal regions with fixed lattice parameters. Though in practice this assumption is not met, artifacts in frequency domain caused by window boundary effects are easily handled by the perimeter regularization, as long as the window size $s$ is chosen reasonably large. Furthermore, these are also reduced through the proposed dimension reduction of the fidelity term (14). Note that crystal images are usually far from periodic at the boundary ($s$ pixels in orthogonal direction) and thus $\mathcal{F}_{\text{FFT}}[\boldsymbol{g}](x)$ cannot be reasonably defined there. Here, we simply extend the segmentation constantly to cover the boundary region.

Figure 2 shows segmentation results obtained by the proposed method with a 2D-FFT modulus based feature extractor of sizes $s=15$ (rows 1-3) and $s=20$ (last row). The crystals in the first three rows consist of regions differing only in crystal orientation. From visual inspection, the results for the noise-free images are exact up to inter-atomic distance. In the first row, despite the noisy grain (third column) begin hardly recognizable due to high noise power (Gaussian noise with a standard deviation of 100% of the maximum noise-free image intensity), the segmentation deviates little from that of the noise-free grain. Similar results are observed for the three-phase scenario (second row). A lower noise power (66%) was chosen, because otherwise the small bottom grain could not be detected, likely due to its small size compared to the window size. The multi-phase segmentation also works very well for five regions (third row) under the presence of very strong noise (100%). Furthermore, as seen in the bottom row of Figure 2, the proposed Fourier-based segmentation is feasible and robust to large amounts of noise (100%), even if the individual grains have entirely different crystal lattices. This is a type of material of practical relevance to material scientists that cannot be handled by the stencil based methods [4, 5, 12].

5 Conclusions

We have discussed a variational framework for multi-phase image segmentation based on structural information from high-dimensional local features. The framework imposes no special constraints on the used indicator functions, except that they are suitable for structure discrimination in the sense that they should be roughly homogeneous inside the structures of interest and provide some contrast across the different regions of interest. A robust initialization strategy for the segmentation algorithm was presented in this context, based on dimension reduction and decorrelation via PCA, as well as edgeness detection and clustering. Numerical results for two applications were presented. For texture segmentation, the proposed framework provides very competitive results, including state-of-the-art performance on the Prague benchmark. Using the 2D-FFT as feature extractor, robust and unsupervised crystal segmentation can be achieved, including segmentation of crystals with entirely different structure from extremely noisy data and without a-priori information about the crystals. We would like to point out that the proposed method can also be applied directly to high-dimensional data, for instance in spectroscopy.

The source code of the proposed method, including executables reproducing all presented results, is available at http://nmevenkamp.github.io/pcams.

Acknowledgments

The authors thank Paul Voyles for providing simulated STEM images of a two-phase crystal.