Unsupervised Multi Class Segmentation of 3D Images with Intensity Inhomogeneities

TL;DR

This paper introduces a novel biconvex variational model for segmenting 3D images with intensity inhomogeneities, combining total variation and multiplicative modeling to improve accuracy in challenging imaging conditions.

Contribution

The paper proposes a new biconvex variational model with a specialized optimization method for effective 3D image segmentation under intensity inhomogeneities.

Findings

Demonstrates superior segmentation performance on 3D FIB tomography images.

Shows convergence of the proposed optimization algorithm.

Achieves robust results in handling intensity inhomogeneities.

Abstract

Intensity inhomogeneities in images constitute a considerable challenge in image segmentation. In this paper we propose a novel biconvex variational model to tackle this task. We combine a total variation approach for multi class segmentation with a multiplicative model to handle the inhomogeneities. Our method assumes that the image intensity is the product of a smoothly varying part and a component which resembles important image structures such as edges. Therefore, we penalize in addition to the total variation of the label assignment matrix a quadratic difference term to cope with the smoothly varying factor. A critical point of our biconvex functional is computed by a modified proximal alternating linearized minimization method (PALM). We show that the assumptions for the convergence of the algorithm are fulfilled by our model. Various numerical examples demonstrate the very good performance of our method. Particular attention is paid to the segmentation of 3D FIB tomographical images which was indeed the motivation of our work.

Figures

Equations

$$F(x)=R(x)L(x),$$

$$f(x)=r(x)+l(x),$$

$$u:=(u_k(j){)}_{j\in \mathcal{G},k\in \{1,\dots ,K\}}.$$

$$\sum _{k=1}^K{u}_k(j)=1.$$

$${△}_K:=\{v=(v_k{)}_{k=1}^K\in [0,1{]}^K:\sum _{k=1}^K{v}_k=1\}.$$

$$\hat{k}:=\underset{k=1,\dots ,K}{\mathrm{argmax}}{u}_k(j)$$

$$\sum _{k=1}^K\sum _{j\in \mathcal{G}}(u_k(j){)}^p|F(j)-C_k{|}^2,$$

$${\mathcal{E}}_{data}(u,C,L):=\sum _{k=1}^K\sum _{j\in \mathcal{G}}(u_k(j){)}^p|F(j)-L(j)C_k{|}^2,$$

$$E_{data}(u,c,l):=\sum _{k=1}^K\sum _{j\in \mathcal{G}}(u_k(j){)}^p|f(j)-l(j)-c_k{|}^2,$$

$$\displaystyle\operatorname*{arg\,min}_{u,C,L}\big{\{}{\cal E}_{data}(u,C,L)++\gamma_{1}\Big{\|}\nabla L\Big{\|}^{2}_{2}+\gamma_{2}\Big{\|}\nabla^{2}L\Big{\|}^{2}_{2}\big{\}}$$

$$\mathrm{subject}\mathrm{to}{\left(u_k(j)\right)}_{k=1}^K\in {△}_K,j\in \mathcal{G},$$

$$\displaystyle\operatorname*{arg\,min}_{u,C,L}\Big{\{}E_{data}(u,c,l)$$

$$\displaystyle\qquad\qquad\qquad+\lambda\sum_{k=1}^{K}\sum_{j\in{\cal G}}(u_{k}(j))^{p}\sum_{i\in{\cal N}_{j}}(f(i)-l(i)-c_{k})^{2}\Big{\}}$$

$$\mathrm{subject}\mathrm{to}{\left(u_k(j)\right)}_{k=1}^K\in {△}_K,j\in \mathcal{G},$$

$$0<\sum _{j\in \mathcal{G}}{u}_k(j)<n,k=1,\dots ,K.$$

$$\displaystyle\operatorname*{arg\,min}_{\phi,C,L}\Big{\{}\int\sum_{k=1}^{K}$$

$$\displaystyle+\lambda\int|\nabla H(\phi)|\,dx+\gamma\int(|\nabla\phi|-1)^{2}\,dx\Big{\}}$$

$$e_k(x):=\int K_{\sigma }(y-x)|F(x)-L(y)C_k){|}^{2}dy.$$

$$\displaystyle\operatorname*{arg\,min}_{u,c,l}\Big{\{}E_{data}(u,c,l)+\sum_{k=1}^{K}\lambda_{k}\Big{\|}|\nabla u_{k}|\Big{\|}_{1}+\gamma\Big{\|}\nabla l\Big{\|}_{2}^{2}\Big{\}}$$

$$\mathrm{subject}\mathrm{to}{\left(u_k(j)\right)}_{k=1}^K\in {△}_K,j\in \mathcal{G},$$

$$u_k:={\left(u_k(j)\right)}_{j\in \mathcal{G}}$$

$$\iota_{S}(x)\coloneqq\left\{\begin{array}[]{ll}0&{\rm if}\;x\in S,\\ +\infty&{\rm otherwise},\end{array}\right.$$

$$E(u,c,l)=E_{data}(u,c,l)$$

$$\displaystyle+\gamma\Big{\|}\nabla l\Big{\|}_{2}^{2}+\iota_{\triangle_{K}^{n}}(u).$$

$$E(u^r,c^r,l^r)\to \underset{(u,c,l)}{\inf }E(u,c,l)\mathrm{as}r\to \infty .$$

$$\sum _{j\in \mathcal{G}}{u}_k^{r_i}(j)\to 0\mathrm{as}i\to \infty ,$$

$$\sum _{j\in \mathcal{G}}{u}_k^r(j)\to 0\mathrm{as}i\to \infty .$$

$$\ker \nabla =\{a{1}_n:a\in \mathbb{R}\},$$

$$l^r=l_{\ker }^r{1}_n+l_{{\ker }^{\perp }}^r,$$

$$\mathcal{K}:=\left\{k:\sum _{j\in \mathcal{G}}{u}_k^r(j)\to 0\mathrm{as}r\to \infty \right\}\subset \{1,\dots ,K\}$$

$${\tilde{u}}_k^r$$

$${\tilde{c}}_k^r$$

$${\tilde{l}}^r$$

$$(\sum _{j\in \mathcal{G}}{u}_k^r(j)|f(j)-l^r(j)-c_k^r{|}^2{)}_{r\in \mathbb{N}}$$

$$(|f(j_k^r)-l^r(j_k^r)-c_k^r|{)}_{r\in \mathbb{N}}=(|f(j_k^r)-l_{{\ker }^{\perp }}^r(j_k^r)-{\tilde{c}}_k^r|{)}_{r\in \mathbb{N}}$$

Full text

1 Introduction

Intensity inhomogeneities often occur in real-world images mainly due to different spatial lighting and deficiencies of imaging devices. For example in MRI, imperfections in the radio-frequency coils or problems associated with acquisition sequences cause intensity changes. The motivation for this paper was the task of segmenting 3D images stemming from focused ion beam (FIB) tomography. While classical $X$-ray tomography does often not reach the required material resolution, FIB tomography enables to investigate the 3D morphology of structures on a scale down to several nanometers. The material is successively removed by a focused ion beam and after every section, the surface is displayed by scanning electron microscopy. Several hundred of these serial slices finally form a 3D image. A typical slice of a 3D FIB tomography of aluminum with silicon carbide (SiC) particles (larger black parts) and copper aggregations (small white parts) is shown in Fig. 1(a). The segmentation has to distinguish between the particles, the aggregates and the surrounding aluminum matrix. However, due to the intensity inhomogeneities, a segmentation based on the gray-values gives a very bad result. Fig. 1(b) and 1(c) show segmentation results for a 3D FIB data set using a supervised gray-value based segmentation method without considering the illumination. Therefore, we have to choose one cluster center for each of the three classes, i.e., in total three gray-values that are close to the classes aluminum, SiC and copper, respectively. Unfortunately, the cluster centers cannot be chosen appropriate for all parts of the image so that either too many or not enough particles are detected. Therefore the segmentation of such images has to take the intensity inhomogeneities into account.

There are several techniques for illumination corrections in the literature. These methods could be used in a preprocessing step before applying a segmentation algorithm of choice. In particular in MRI, intensity corrections were proposed by simple homomorphic filtering [12, 21] and polynomial, resp., spline surface fitting approaches [16, 31, 42, 17]. Many spatial illumination correction methods for natural images take hypotheses about the Human Visual System (HVS) into account. In particular, the perceptual work about the Retinex model [22] has found wide acceptance. It states that the HVS does not perceive an absolute lightness but rather a relative, local one. This phenomenon is called lateral inhibition. For example, Fig. 2 shows the experiment of the checkerboard shadow illusion of Adelson [1]. Although the squares A and B have the same gray-value, the perceived intensities are different. In the Retinex model, the light intensity $F$ perceived by the observer or camera is considered to be the product of the reflectance of the objects $R$ in the scene and the amount of source illumination $L$ falling on the objects, see also [18, p. 51],

[TABLE]

where $R(x)\in(0,1)$ and $L(x)\in(0,+\infty)$. While we assume that the reflectance inherits the structures of the objects, e.g. edges, we consider the illumination as spatially smooth, in particular it should not have sharp edges, which represent the image structures. Taking the logarithm in (1) we obtain

[TABLE]

where $f=\log F$, $r=\log R$ and $l=\log L$. Retinex based variational or PDE based approaches for illumination corrections can be found for example in [28, 30, 33, 34].

In MRI the observed intensity is often modeled similarly as in (1) by the product of a structural part $R$ and a so called gain factor $L$. In this paper we also follow the multiplicative intensity model.

In contrast to a two step procedure we consider the simultaneous segmentation and intensity inhomogeneity estimation. This avoids the computational burden of two separate procedures and has moreover the advantage of being able to use intermediate information from the segmentation while performing the update. Besides statistical methods as the computationally extensive EM approach in [43], variational based algorithms were proposed in [2, 3, 25, 26, 29, 36, 45]. We consider the later approaches in more detail in the next section.

Variational segmentation models have shown a very good performance and flexibility in many applications. Level set methods and convex segmentation models which penalize the (nonlocal) discrete total variation (TV) of a relaxed label assignment matrix have been successfully applied [5, 11, 15, 23, 24, 37, 40, 44].

In this paper, we combine the TV based segmentation method with the multiplicative intensity model. This results in a biconvex non smooth functional, which has to be optimized with respect to the label assignment matrix, the cluster centers and the smoothly varying intensity factor. Then we obtain both a segmented image and an estimation of the intensity inhomogeneities. We compute a critical point of the corresponding functional by applying a slight modification of the proximal alternating linearized minimization method (PALM) by Bolte et al. [10].

The paper is organized as follows: In Section 2 we review variational segmentation models which compensate for the varying intensities while segmenting the image. Then we introduce our model. The modified PALM is explained in Section 3. Section 4 contains numerical examples. The paper finishes with conclusions in Section 5.

2 New Model

Let ${\cal G}\coloneqq\{1,\ldots,n_{1}\}\times\ldots\times\{1,\ldots,n_{d}\}$ be a $d$-dimensional image grid. Here we will deal with $d=2$ and $d=3$. Let $n=n_{1}\ldots n_{d}$ be the number of image pixels. We consider images $F\colon{\cal G}\rightarrow\mathbb{R}$ which we want to segment into $K$ classes. We introduce a so-called label assignment matrix

[TABLE]

Ideally we would like to have for a fixed pixel $j\in{\cal G}$ that $u_{k}(j)=1$ if it belongs to class $k$ and $u_{k}(j)=0$ otherwise. However, in the subsequent variational approach this would lead to optimization tasks which are NP hard to solve. Therefore it is common to relax the assumptions on the label assignment matrix and require only that $u_{k}(j)\in[0,1]$ is the probability that pixel $j$ belongs to class $k$. Since every pixel should be assigned to one of the classes, we enforce for every $j\in{\cal G}$ that

[TABLE]

In other words $(u_{k}(j))_{k=1}^{K}$ is an element of the probability simplex

[TABLE]

From the label assignment matrix the segmentation can be obtained by assigning, e.g., the label

[TABLE]

to $j\in{\cal G}$.

A variational model is composed of a data term, which incorporates information about the given image $F$ and penalizing/regularizing terms, which contain prior information on the image. Let $C\coloneqq(C_{1},\ldots,C_{K})^{\scriptscriptstyle\operatorname{T}}$ denote the unknown vector of the centers of the $K$ gray-value clusters. For segmentation problems a typical data term is given by

[TABLE]

where $p\geq 1$. For $p=1$ this data term appears in the $K$ (or $C$) means approach, while $p>1$ is related to the fuzzy $K$-means method [7, 8, 20, 32]. Roughly speaking, the data term indicates the following: If $F(j)$ is close to the center $C_{k}$, then the corresponding summand in (4) remains small for a larger $u_{k}(j)$. Thus the probability of $F(j)$ being in the class $k$ is high. If $F(j)$ is far away from the center $C_{k}$, then the corresponding term in (4) becomes only small if the probability $u_{k}(j)$ that $F(j)$ belongs to class $k$ is small. Clearly other distance measures between $F(j)$ and $C_{k}$ than the squared absolute value can be chosen. Moreover, the approach can be generalized to color images or more general image feature vectors $F$ than gray-values. In this paper we restrict our attention to gray-value images. The above data term does not take intensity inhomogeneities into account. However, as shown in Fig. 1 it is often not possible to find appropriate class centers if the illumination varies. A remedy consists in incorporating the smoothly varying intensity part $L$ into the model and consider

[TABLE]

$p\geq 1$, or its logarithmic counterpart

[TABLE]

where $f=\log F$, $c=\log C$ and $l=\log L$. From (6) we get the illumination by $L=\exp(l)$ and the class centers by $C=\exp(c)$. The functional (5) is triconvex, i.e., fixing two of the variables $u,C,L$ the functional is convex in the third variable. In contrast, the functional (6) is only biconvex in $u$ and $v\coloneqq(c,l)$. This means if we fix $u$, then the functional is convex in $v$ and conversely.

Within the penalizers pixel differences are often used since they indicate smooth or higher frequent image parts as edges. By $\nabla F$ we denote the $d$-dimensional discrete gradient of $F$, where we use forward differences in each direction together with mirror boundary conditions. For the concrete matrix structure of the gradient operator we refer to [40]. Similarly, $\nabla^{2}F$ is the discrete Hessian of $F$. Throughout this paper $\|A\|_{2}$ denotes the square root of the sum of the squared entries of a multidimensional array $A$.

Let us briefly review existing variational models for image segmentation which take intensity inhomogeneities into account. In [36] the following model with data term (5) with $p=2$ was suggested:

[TABLE]

where $\gamma_{i}>0$, $i=1,2$ are regularization parameters.

This model only penalizes non smooth intensity inhomogeneities $L$ via two terms, but does not take the edges of the label assignment matrix $u$ into account. Therefore it is not robust to noise. Moreover it has the drawback that there does not exist a minimizer. This can be simply seen by setting $L\coloneqq F/r$ and $C_{k}\coloneqq r$ with some constant $r>0$ so that the functional becomes \frac{\gamma_{1}}{r}\Big{\|}\nabla L\Big{\|}^{2}_{2}+\frac{\gamma_{2}}{r^{2}}\Big{\|}\nabla^{2}L\|_{2}^{2}. For $r\rightarrow+\infty$ we obtain a minimizing sequence of the functional which does not converge.

Another segmentation model for MRI with data term (6) was proposed by Ahmed et al. [2]:

[TABLE]

Here $\lambda>0$ is a regularization parameter. The second term takes the illumination in the neighborhood ${\cal N}_{j}$ of the $j$-th pixel into account in order to respect the lateral inhibition. However, as in the previous model no penalizer for the label assignment matrix is used, which makes the model sensitive to noise. A similar model which penalizes the simplex constraint was suggested in [29]. Unfortunately, the model formulation in [29] is mathematically not sound and the optimization procedure does not fit to the model.

Li et al. [25, 26] considered a variational level set model, which reads in the continuous setting as

[TABLE]

with

[TABLE]

Here $K_{\sigma}$ is a truncated Gaussian kernel with standard deviation $\sigma$, $H$ a smoothed Heaviside function and $M_{k}$ is a membership function. The model was applied for two-dimensional medical images. Note that the functional is not convex in $\phi$ for fixed $C$ and $L$. A more general level set approach is considered in [45]. Furthermore, a slightly different approach was proposed also in [3].

In this paper we consider the model

[TABLE]

where $p=1$ in the data term. Here $\lambda>0$ and $\gamma>0$ are regularization parameters and

[TABLE]

is the $k$-th label assignment matrix. We have not found this model in the literature. Using the indicator function of sets

[TABLE]

we can express the constraint “subject to $\left(u_{k}(j)\right)_{k=1}^{K}\in\triangle_{K},\;j\in{\cal G}$” by adding $\iota_{\triangle_{K}^{n}}(u)$ to the functional. Hence, the functional in (16) can be written as

[TABLE]

As in (7) we penalize a non smooth illumination $l$, but just by one quadratic first order difference term. The matrix $u_{k}$ is penalized by its discrete total variation (TV) [39], where $|\nabla u_{k}|$ denotes the length of the discrete gradient of $u_{k}$ and the $\ell_{1}$ norm is taken over the grid points $j\in{\cal G}$. This term encourages smooth edges of the segmented objects and makes the model more robust to noise. We allow the choice of different parameters $\lambda_{k}$ to better segment objects of different size. The functional (16) is biconvex in $u$ and $(c,l)$ and possesses a minimizer as the following proposition states.

Proposition 1**.**

The functional $E(u,c,l)$ defined by (18) has a minimizer.

Proof.

It suffices to show that there exists a bounded infimal sequence since $E(u,c,l)$ is lower semi-continuous and bounded from below.

Let $((u,c,l)^{r})_{r\in\mathbb{N}}$, $(u,c,l)^{r}\in\triangle_{K}^{n}\times\mathbb{R}^{K}\times\mathbb{R}^{n}$, be an infimal sequence, i.e.,

[TABLE]

Without loss of generality, we can assume the following: If there exists a subsequence $(u^{r_{i}}_{k})_{i\in\mathbb{N}}$ for some $k\in 1,\ldots,K$ with

[TABLE]

then the whole sequence converges to zero, in other words

[TABLE]

If this was not the case, we could restrict the following analysis to the sequence $((u,c,l)^{r_{i}})_{i\in\mathbb{N}}$ and possibly repeat the procedure for the next $k$ fulfilling (21). Note that it is not possible that (22) holds for all $k\in\{1,\ldots,K\}$ since $u^{r}\in\triangle_{K}^{n}$.

Since the sequence is assumed to be an infimal sequence, we immediately obtain by the constraint on $u^{r}$ that $(u^{r})_{r\in\mathbb{N}}$ is bounded.

Next we consider $l$. Due to mirror extension at the boundary, we have that

[TABLE]

where $\mathbf{1}_{n}$ is the vector consisting of $n$ entries $1$, i.e., the kernel consists of constant images. We decompose $l$ into two parts

[TABLE]

where $l_{\ker^{\perp}}^{r}\in\ker(\nabla)^{\perp}$ and $l_{\ker}^{r}\mathbf{1}_{n}\in\ker(\nabla)$ with $l_{\ker}^{r}\in\mathbb{R}$. Using again that the sequence is an infimal one, we see due to the term $\|\nabla l\|_{2}^{2}$ that $(l^{r}_{\ker^{\perp}})_{r\in\mathbb{N}}$ is bounded.

Now define the (possibly empty) set

[TABLE]

and fix $\hat{k}\not\in\mathcal{K}$. We show that a bounded infimal sequence is given by $((\tilde{u},\tilde{c},\tilde{l})^{r})_{r\in\mathbb{N}}$, where

[TABLE]

We immediately obtain that $(\tilde{u}^{r})_{r\in\mathbb{N}}$ is bounded as $(u^{r})_{r\in\mathbb{N}}$ is bounded. Further, $(\tilde{l}^{r})_{r\in\mathbb{N}}=(l^{r}_{\ker^{\perp}})_{r\in\mathbb{N}}$ is bounded. Next, we show that $(\tilde{c}^{r})_{r\in\mathbb{N}}$ is bounded. From the data term (6), we conclude that

[TABLE]

is bounded for all $k=1,\ldots,K$. Let $k\not\in\mathcal{K}$. By our assumption that (21) implies (22), we know that there also does not exist a subsequence $(r_{i})_{i\in\mathbb{N}}$ with $\sum_{j\in\mathcal{G}}u^{r_{i}}_{k}(j)\rightarrow 0\ \mathrm{as}\ i\rightarrow\infty$. Thus, for $k\not\in\mathcal{K}$, there must exist a series of pixels $(j_{k}^{r})_{r\in\mathbb{N}}$ such that $u_{k}^{r}(j_{k}^{r})\geq\varepsilon>0$ for $r$ sufficiently large. Hence,

[TABLE]

is bounded. Consequently, $(\tilde{c}_{k}^{r})_{r\in\mathbb{N}}$ is bounded for all $k\not\in\mathcal{K}$. Together with $\tilde{c}_{k}^{r}=0$ for $k\in\mathcal{K}$, we obtain that $(\tilde{c}^{r})_{r\in\mathbb{N}}$ is bounded. Thus, the whole sequence $((\tilde{u},\tilde{c},\tilde{l})^{r})_{r\in\mathbb{N}}$ is bounded.

It remains to show that it is an infimal one. Note that adding $l^{r}_{\ker}$ to $c_{k}^{r}$ and subtracting it from $l^{r}(j)$ does not change the objective value.

By construction, the constraint $\tilde{u}^{r}\in\triangle_{K}^{n}$ is still fulfilled. Thus, the only change in the objective value arises for the summands belonging to $k\in\mathcal{K}$ and $\hat{k}$. But since

[TABLE]

and

[TABLE]

the objective value does not change in the limit such that the sequence $((\tilde{u},\tilde{c},\tilde{l})^{r})_{r\in\mathbb{N}}$ is still an infimal sequence. This finishes the proof. ∎

Remark 2**.**

Alternatively to (16) we could deal with the non logarithmic model

[TABLE]

We prefer to work with the biconvex functional (16), although we have obtained similar numerical results by minimizing the above triconvex functional.

3 Minimization Algorithm

In this section we deal with the minimization of our functional (18). First we mention that fixing $(c,l)$ we get

[TABLE]

which is convex, but not strictly convex. Fixing $u$ we obtain the convex functional

[TABLE]

This functional becomes strictly convex if we assume that we have only nonempty classes, i.e.,

[TABLE]

and enforce the mean value of $l$ to be zero. In other words,

[TABLE]

is strictly convex in $c$ and $l$ separately due to the quadratic term. It is furthermore jointly convex in $(c,l)$ by the term $\iota_{\{0\}}(\langle\mathbf{1},l\rangle)$.

For an algorithm which computes in an alternating way

[TABLE]

one can prove similarly as in [6] that the resulting sequence $\{u^{(r)},c^{(r)},l^{(r)}\}_{r}$ has a subsequence which converges to a partial minimizer of (18). For an alternating convex search algorithm for general biconvex problems we refer to [19].

In this paper we want to apply an algorithm which ensures the convergence of the whole sequence $\{u^{(r)},c^{(r)},l^{(r)}\}_{r}$. To this end we need the definition of the proximal mapping. For a proper, lower semi-continuous, convex function $h\colon\mathbb{R}^{N}\rightarrow(-\infty,+\infty]$ and $\lambda>0$, the proximal operator $\mathrm{prox}_{\lambda h}\colon\mathbb{R}^{N}\rightarrow R^{N}$ is defined by

[TABLE]

Indeed, the minimizer of the right-hand side exists and is uniquely determined, see, e.g., [38].

Recently, Bolte, Sabach and Teboulle [10] considered problems of the form

[TABLE]

where $h_{i}\colon\mathbb{R}^{N_{i}}\rightarrow(-\infty,+\infty]$, $i=1,\ldots,m$, are proper lower semi-continuous functions and $H\colon\mathbb{R}^{N_{1}}\times\ldots\times\mathbb{R}^{N_{m}}\rightarrow\mathbb{R}$ is continuously differentiable. Further, $\sum_{i=1}^{m}h_{i}(x_{i})+H(x)$ needs to be a Kurdyka-Łojasiewicz (KL) function. For a discussion of KL functions, we refer to [4, 9, 10]. Here, we only want to emphasize that semi-algebraic functions are KL and our setting fits to the examples for semi-algebraic functions given in [10, Sect. 5]. The authors suggested Algorithm 1 for such problems.

Under certain assumptions it was shown that the sequence $\{x^{(r)}\}_{r}$ converges to a critical point of the functional in (39).

Theorem 3**.**

Let $h_{i}\colon\mathbb{R}^{N_{i}}\rightarrow(-\infty,+\infty]$, $i=1,\ldots,m$ be proper lower semi-continuous functions and $H\colon\mathbb{R}^{N_{1}}\times\ldots\times\mathbb{R}^{N_{m}}\rightarrow\mathbb{R}$ a continuously differentiable function such that $\nabla H$ is Lipschitz continuous on bounded sets. Let $h_{i}$, $i=1,\ldots,m$ and $\sum_{i=1}^{m}h_{i}+H$ be bounded from below. Further, let $\sum_{i=1}^{m}h_{i}(x_{i})+H(x)$ be a KL function. Assume that the sequence of iterates $\{x^{(r)}\}_{r}$ produced by PALM is bounded. Further suppose that for $i=1,\ldots,m$, the partial gradients $\nabla_{x_{i}}H(x)$ are globally Lipschitz for fixed

[TABLE]

with Lipschitz constant

[TABLE]

where

[TABLE]

for all $r\in\mathbb{N}$ and some $\lambda_{i}^{+}\in\mathbb{R}$. Then, for

[TABLE]

the sequence $\{x^{(r)}\}_{r}$ converges to a critical point of the functional in (39).

For minimizing (18) we apply PALM with $m=3$ blocks and the following setting

[TABLE]

More precisely,

[TABLE]

i.e., we have modified the functional (18) by the additional function $h_{3}$, which enforces the logarithmic illumination to have mean value zero, see (37). Note that adding a constant $C\in\mathbb{R}$ to $l(j)$ for all $j\in{\cal G}$, and subtracting this constant from the $c_{k}$, $k=1,\ldots,K$ does not change the value of the functional (18). Therefore we have an ambiguity in the solution consisting of a constant function. By introducing $h_{3}$ we enforce that $l$ has mean zero, i.e., we fix the constant to $C\coloneqq-\langle\mathbf{1},l\rangle/n$ such that the ambiguity is removed.

Of course it would also be possible to use PALM with $m=2$ blocks by handling $(c,l)$ together. However, this leads to joint estimates for the Lipschitz constant in the partial gradient of $H$ with respect to $(c,l)$ which is more restrictive than the Lipschitz constants from the separate gradients with respect to $c$ and $l$, see also proof of Corollary 4.

In PALM we need the following partial gradients of $H$:

[TABLE]

Based on these gradients PALM reads for our setting as follows:

The first proximum $\mathrm{prox}_{\tau_{1,r}h_{1}}$ is not given analytically, but can be computed by several methods. Here we use the primal dual algorithm proposed in [14]. The second step is just a gradient descent step. The last proximum is a projection of $a\coloneqq l^{(r)}-\frac{1}{\tau_{3,r}}\nabla_{l}H(u^{(r+1)},c^{(r+1)},l^{(r)})$ onto the hyperplane $\{l:\langle\mathbf{1},l\rangle=0\}$ and can be computed by subtracting from $a$ its mean value.

Corollary 4**.**

For the functions $h_{i}$, $i=1,2,3$ and $H$ defined in (43) and

[TABLE]

Algorithm 2 converges to a critical point of (18).

Proof.

We have to check that the assumptions of Theorem 3 are fulfilled.

It is easy to check that $H(u,c,l)+h_{1}(u)+h_{2}(c)+h_{3}(l)$ is a semi-algebraic function and therefore a KL function, see the examples in [10, Sect. 5]. The functions $h_{1},h_{2},h_{3}$ are lower semi-continuous and bounded from below. Since $H$ is twice continuously differentiable, its gradient is Lipschitz on bounded sets and by $u\in\triangle_{K}^{n}$ the function is also bounded from below.

Let us consider the partial gradients of $H$ given by (49). Since $\nabla_{u}H(u,c,l)$ does not depend on $u$, we get immediately $L_{1}(c,l)=0$. For the gradient with respect to $c$ it follows with $u\in\triangle_{K}^{n}$ that

[TABLE]

Finally, $\nabla_{l}H(u,c,l)$ is Lipschitz since

[TABLE]

where

[TABLE]

The $d$-dimensional Laplacian $\nabla^{*}\nabla$ has a spectral norm smaller than $4d$. This can be seen as follows. In 1D one simply has the forward difference matrix with mirror boundary conditions

[TABLE]

The eigenvalues of $\nabla^{*}\nabla$ are $4\sin^{2}(\frac{\pi j}{2n})$, $j=0,\ldots,n_{1}-1$, see, e.g., [41]. Hence, $\|\nabla^{*}\nabla\|<4$. In $d$ dimensions, after reshaping the $d$-dimensional image into a vector, the gradient becomes

[TABLE]

where $a_{i}\coloneqq\Pi_{j=1}^{i-1}n_{j}$, $b_{i}\coloneqq\Pi_{j=i+1}^{d}n_{j}$ and $\otimes$ denotes the tensor product of matrices. Then

[TABLE]

Since the eigenvalues of $A\otimes B$ are given by the products of the eigenvalues of $A$ and $B$, one obtains $\|\nabla^{*}\nabla\|\leq\sum_{i=1}^{d}\|D_{n_{i}}^{\scriptscriptstyle\operatorname{T}}D_{n_{i}}\|\leq 4d$. Since $u^{(r)}\in\triangle_{K}^{n}$, we see that $L_{3}(u^{(r)},c^{(r)})\leq 2+8d\gamma$.

Finally, the sequence of iterates produced by the algorithm is bounded by the following arguments. By adding the constraint to $l$, i.e., $\mathbf{1}^{\scriptscriptstyle\operatorname{T}}l=0$, the resulting functional becomes coercive. Together with the decrease property for the PALM iterates [10, Lemma 3.3] one immediately gets the boundedness. ∎

4 Numerical Examples

In this section we demonstrate the performance of our algorithm. The algorithm was implemented in MATLAB with the following general parameter setting:

• $\bullet$

Lipschitz constants in PALM: According to the Lipschitz constants we can choose any $\tau_{1}^{r}>0$. In practice $\tau_{1}^{r}$ needs to be chosen rather small to achieve an acceptable convergence speed, here we set $\tau_{1,r}\coloneqq 10^{-6}$. Further we set $\tau_{2,r}\coloneqq n$. Although this choice does not fulfill the restrictions given by the Lipschitz constants, we observed numerical convergence. In particular compared to a smaller stepsize that fulfills the Lipschitz condition, we got the same results in a faster way. Finally, we set $\tau_{3,r}\coloneqq 2+8d\gamma$ according to the Lipschitz condition.

• $\bullet$

Initialization: The illumination $l$ is initialized as a very smooth version of the input image obtained by convolving the input image with a Gaussian kernel of large standard deviation $\sigma$ displayed in the caption of the images. The codebook $c$ is initialized with a rough initial guess, and $u$ with constant entries $\frac{1}{K}$. We observed that the initialization of $u$ and $c$ is not important. For reasonable initial values the results were almost equal.

• $\bullet$

Iteration number: In all our experiments we performed $2000$ (outer) iterations of PALM. In each outer iteration, we applied $50$ (inner) iterations of the primal dual method for computing the proximal operator in the first step of PALM. If we would apply the inner PDHG iterations on their own, usually more than $50$ iterations would be necessary. Here, we observed that $50$ inner iterations are sufficient since the initialization is already close to the solution after some outer iterations.

• $\bullet$

Regularization parameters: If not stated otherwise, the other parameters, i.e., $\lambda_{k}$, $k=1,\ldots,K$ and $\gamma$ were chosen according to the best visual impression and are stated in the captions of the corresponding figures.

We compare our approach with the following three other methods:

• M1)

Model (13) proposed in [26] for 2D images, where we used the program package of the authors available at http://www.engr.uconn.edu/~cmli/code/.

• M2)

We apply the pure segmentation model

[TABLE]

proposed, e.g., in [40]. This model does not take care of illumination changes. We use this method for comparison to show that it is necessary to include the illumination in the model.

• M3)

We apply a two step procedure. In the first step we use the initialization described above to estimate the illumination and correct the image by dividing it by the illumination. In the second step we apply model M2) to the corrected image.

We start with an artificial 2D image as a ground truth experiment. Next we apply our algorithm to 2D medical images, in particular CT images, since this kind of images was used to test the algorithm M1) in [26]. Finally, we consider 3D FIB tomography data which was the motivation for this work. Since some of the images have pixels with value zero, we add a small constant to avoid problems when taking the logarithm of such images.

Artificial 2D Images

We start with the artificial 2D image with varying illumination in Fig. 3. The figure shows the noise-free case. A three class segmentation by M2) leads to completely wrong results. Our model is able to segment this image exactly and extract the proper illumination $l$. Similar results can be obtained by the method M1) [26].

Next, we add Gaussian noise of standard deviation $s$ and compare the segmentation of our model and M1). For our model we optimized the TV parameter $\lambda=\lambda_{1}=\lambda_{2}$ for each noise-level by a grid search according to the smallest number of wrongly assigned pixels, and took $\gamma=100$ as in the noise-free case. In the model (13) from [26] we set $\sigma=7$ which also provided the best result in the noise-free case. As in our model, the parameter $\lambda$ was optimized by a grid search for each noise-level. Additionally, we optimized the parameter $\gamma$ by a grid search. Fig. 5 shows the percentage of wrongly assigned pixels depending on $s$ for our model and M1). Up to the noise level $s=0.004$ the proposed method clearly outperforms M1). For $s=0.005$ however, one class vanishes due to the noise so that the number of wrongly assigned pixels increases significantly. To avoid this the third center $c_{3}$, which belongs to the vanishing class, can be fixed. Then the proposed method outperforms M1) also for the higher noise levels. Fig. 4 depicts the results for the noise level $s=0.005$.

Medical Images

Next we present results for different medical images. The images shown in Fig. 6 were taken from [25, 27].

The first row shows results for an X-ray image of bones. The result of our model is compared to the result of M1) with the parameters proposed in the paper [26], i.e., $\sigma=4,\gamma=1$ and $\lambda=0.001\cdot 255^{2}$. The white part corresponds to positive values of the resulting function of the level set method M1). Furthermore, we depict the result of the two-step method M3) to show that correcting the illumination only in the preprocessing step before the segmentation is not sufficient. The methods M1) and M3) perform equally well, where each one shows different slight artifacts. The proposed model gives the best segmentation result, in particular the left bone and the upper part of the right bone are segmented correctly.

The second and third row show results for two CT angiography images of vessels. It is clearly necessary to incorporate the illumination since parts of the images differ widely in their brightness. For the first vessel image the parameters proposed in [26] were used. For the second vessel we set $\sigma=10,\gamma=10$. The first segmented vessel by our method shows slightly thicker structures. In the image center the segmentation method M1) therefore gives better results. In the lower part M1) leads to too thin structure. Here our method performs better. M3) leads to the worst results in this example. Many thin structures are not corrected and some additional artifacts are visible. For the second vessel all three methods give good results. However, M1) and M3) show some small artifacts that are not visible with the proposed method.

3D FIB Tomography Images

Next we are interested in the segmentation of 3D images stemming from FIB tomography. We consider 3D data of an aluminum matrix composite. More precisely, the material consists of three phases which we want to segment: the first phase is aluminum, the second one consists of silicon carbide (SiC) particles and the third one are copper aggradations. The aluminum particles are darker than the surrounding aluminum matrix. The copper aggradations are visible as small bright spots.

In Fig. 7 and 8 two 3D data sets of this material with different particle sizes are examined. The necessity to take the illumination into account when segmenting the image in Fig. 7 was already discussed in Fig. 1. In Fig. 7(a)-7(b) and Fig. 8(a)-8(b), respectively, an exemplary slice of the data set is shown together with the slice of the segmented data. In addition a visualization of the segmented 3D data set is provided in Fig. 7(c)-7(d) and Fig. 8(c)-8(d), respectively. Note that due to the very different size of SiC and copper we had to choose a small parameter for the aggradations ($\lambda_{3}$) and a larger one for the particles ($\lambda_{2}$) to segment both correctly. Furthermore we fixed $c_{3}$, which corresponds to the copper aggradations, to avoid that this small segment vanishes similar as in Fig. 4.

5 Conclusions

We proposed a novel biconvex model for segmentation of images with intensity inhomogeneities. that combines a total variation based approach for segmentation with a multiplicative model for the intensity inhomogeneities. This results in a simultaneous segmentation and intensity correction algorithm without preprocessing. For the minimization of the resulting functional we used the PALM algorithm for which we could prove convergence to a critical point. We applied our model very successfully to 3D FIB tomography images and showed that the proposed model leads to good results compared to a state-of-the-art method for various 2D medical images. In future work, we want to address the question how good the relaxed two-class model solves the original binary problem with $u_{k}(j)\in\{0,1\}$, see [35]. The incorporation of the illumination into a segmentation framework with additional linear operators as, e.g., blur, see [13] appears also to be useful. Finally, we want to handle other than gray-valued images.

Acknowledgement. Funding by the German Research Foundation (DFG) within the Research Training Group 1932 ”Stochastic Models for Innovations in the Engineering Sciences”, project area P3, is gratefully acknowledged.

The authors thank K. Schladitz (Fraunhofer ITWM, Kaiserslautern) and F. Balle, T. Beck and S. Schuff (Department of Mechanical and Process Engineering, University of Kaiserslautern) for fruitful discussions. We are thankful to T. Löber (Nano Structuring Center Kaiserslautern) for creating FIB/SEM data of the aluminum matrix composite used in Figure 1, 7 and 8.