Combining Contrast Invariant L1 Data Fidelities with Nonlinear Spectral Image Decomposition

TL;DR

This paper introduces a novel combination of contrast-invariant L1 data fidelities with nonlinear spectral decompositions, enhancing multi-scale image analysis and segmentation, especially in biomedical imaging.

Contribution

It presents the first integration of L1-based multi-scale methods with nonlinear spectral decompositions for improved image segmentation.

Findings

L1-TV promotes sparsity in spectral responses.

Contrast invariance improves multi-scale decomposition.

Enhanced segmentation results on biomedical images.

Abstract

This paper focuses on multi-scale approaches for variational methods and corresponding gradient flows. Recently, for convex regularization functionals such as total variation, new theory and algorithms for nonlinear eigenvalue problems via nonlinear spectral decompositions have been developed. Those methods open new directions for advanced image filtering. However, for an effective use in image segmentation and shape decomposition, a clear interpretation of the spectral response regarding size and intensity scales is needed but lacking in current approaches. In this context, $L^1$ data fidelities are particularly helpful due to their interesting multi-scale properties such as contrast invariance. Hence, the novelty of this work is the combination of $L^1$-based multi-scale methods with nonlinear spectral decompositions. We compare $L^1$ with $L^2$ scale-space methods in view of spectral image representation and decomposition. We show that the contrast invariant multi-scale behavior of $L^1-TV$ promotes sparsity in the spectral response providing more informative decompositions. We provide a numerical method and analyze synthetic and biomedical images at which decomposition leads to improved segmentation.

Equations

$$\underset{u}{\min }\frac{1}{p}{\Vert u-f\Vert }_{L^p}^p+\alpha J(u)$$

$$\text{vspace}-.5em0\in \partial {\Vert {\partial }_{t}u(t,x)\Vert }_{L^p}+\partial J(u(t,x)).$$

$$\text{vspace}-.5em{\partial }_{t}u=-q(u)\text{with }q\in \partial J(u).$$

$$\lambda u\in \partial J(u).$$

$$\underset{x\in \partial C}{\text{ess sup}}{\kappa }_C(x)\le \frac{Per(C)}{|C|}\text{(C is Cheeger in itself)}$$

$$u(t_{\alpha },x)=\{\begin{array}{ll}\left(c-\frac{2}{r}{t}_{\alpha }\right)\cdot \text{mathbbm}{1}_{B_r(0)}(x)& \text{mbox}if0\le t_{\alpha }<\frac{cr}{2}\\ 0& \text{mbox}otherwise.\end{array}$$

$$\varphi (t,x)=u_{tt}(t,x)\cdot t\text{and}S(t)={\Vert \varphi (t,x)\Vert }_{L^1}.$$

$$f(x)={\int }_0^{\infty }\varphi (t,x)dt+\bar{f}$$

$$f_H(x)={\int }_0^{\infty }H(t)\varphi (t,x)dt+H(\infty )\bar{f}$$

$$\underset{u}{\min }{\Vert u-f\Vert }_{L^1}+t_{\alpha }TV(u).$$

$$u(t_{\alpha },x)=\{\begin{array}{ll}f& \text{mbox}if0\le t_{\alpha }<\frac{r}{2}\\ c^{\prime }\cdot f\text{ with }{c}^{\prime }\in [0,1]& \text{mbox}if{t}_{\alpha }=\frac{r}{2}\\ 0& \text{mbox}otherwise.\end{array}$$

$$\varphi (t,x)=-u_t(t,x)\text{ fulfilling }f(x)={\int }_0^{\infty }\varphi (t,x)dt+\hat{c}.$$

$$S^2(t)=⟨\varphi (t,x),f(x)⟩.$$

$${\Vert f\Vert }^2=⟨f,f⟩={\int }_0^{\infty }⟨\varphi (t,x),f(x)⟩\mathrm{d}t={\int }_0^{\infty }{S}^2(t)\mathrm{d}t.$$

$$f_H(x)={\int }_0^{\infty }H(t)\varphi (t,x)dt+H(\infty )\hat{c}$$

$$f_{\text{seg},H}(x)=\left({\int }_0^{\infty }H(t)\varphi (t,x)dt>0\right).$$

$$⟨\nabla u,g⟩+{\Vert u-f\Vert }_{L^1}-t_{\alpha }{\delta }_P(g)⟶\underset{u}{\min }\underset{g}{\max }$$

Full text

1 Introduction

In imaging science, the solution of inverse problems is often addressed by the modeling and analysis of variational methods of the form:

[TABLE]

where $f$ denotes a noisy signal, $u$ is a desired image function defined on $\Omega$ in $\mathbb{R}^{2}$. The data fidelity measures the residual in the $L^{p}$ norm and $J$ is a regularization functional that has a weighting parameter $\alpha\geq 0$. In this paper, we concentrate on image denoising methods with convex and one-homogeneous regularization functionals $J$ which can address image decomposition and segmentation adequately. More specifically, we focus on total variation regularization $J(u)=TV(u)$ and norms with $p=1$ versus $p=2$. One option to generate a scale-space method is a gradient flow, based on the functional (1) with initial condition $u(0,x)=f(x)$ and subdifferential inclusions forming the doubly nonlinear evolution equation:

[TABLE]

In the case of $p=2$, this simplifies to well-known gradient flows of the form:

[TABLE]

The regularization parameter $\alpha$ is now hidden in the time dependency. Corresponding inverse scale space gradient flows can be constructed via Bregman distances [1, 2]. The analysis of linear eigenvalue problems and spectral decompositions, e.g. via the Fourier transform, is a well-known and widely used theory in the fields of signal, image and graph-based data processing. Due to the continuing success of nonlinear regularization functionals in imaging, there is strong interest in generalizing spectral theory to the nonlinear case. The general idea is to examine solutions to the nonlinear eigenvalue problem:

[TABLE]

In [3, 4], Gilboa introduced the idea of nonlinear spectral decompositions for the TV transform. By transferring solutions (eigenfunctions) of (4) to sparse peaks in a spectral domain, the idea of advanced filters, suppressing or enhancing image components similar to the Fourier transform, came about. This concept was studied for one-homogenous functionals [5, 6, 7] and scale-space flows of the form (3) with $p=2$. In this way, a decomposition of the input signal $f$ into significant components is possible while an exact representation of $f$ can still be guaranteed.

For $p=1$, the gradient flow in (2) is interesting and more challenging than for $p=2$. Another way to obtain a forward scale-space is to construct a sequence of variational problems of the form (1) with increasing regularization parameter $t_{\alpha}$ replacing the fixed $\alpha$. Here, the scale parameter $t_{\alpha}$ corresponds to the time variable $t$ used in (2). In this paper, we will focus on this type of model particularly in view of $p=1$. From pioneering works on nonlinear $L^{1}$ filtering [8, 9, 10, 11], it is well known that such variational reconstruction methods share interesting multi-scale properties including contrast invariance. For this reason, $L^{1}$ data fidelities have been successfully used for advanced image reconstruction techniques [12, 13, 14], vector field estimation [15] and image decompositions regarding texture [10, 16]. The special multi-scale behavior becomes clear in Figure 1. In the $L^{1}-TV$ case shown on the right, the contrast invariance leads to plateaus in the scale-space graph indicating an abrupt disappearance of TV-eigenshapes regarding the input image on the left.

Motivated by [9, 11, 17] and nonlinear spectral methods [7, 18], the main goal of this work is to study $L^{1}$ versus $L^{2}$ in view of the sparsity of nonlinear spectral decompositions. Does $L^{1}$ imply sparsity and hence a more informative spectral response? Can we expect more reliable image decompositions via backtransformation facilitating improved image segmentation? What happens if complex shapes, e.g. nonconvex or compositions of eigenfunctions, are involved?

2 Modeling

In the following section we first give a short overview of the eigenshapes of the total variation functional and calibrable sets before we introduce the spectral framework for $L^{2}-TV$ denoising in more detail. Finally, we show why a combination of the spectral framework with $L^{1}-TV$ denoising seems promising and how the spectral framework can be adapted in this case.

2.1 Geometry: Eigenfunctions and Calibrable Sets

In the introduction in Figure 1 we have already seen an input image composed of TV eigenfunctions. A general geometric description of TV eigenfunctions is given in [19, 20, 21] in terms of calibrable sets, more precisely by convex sets which are Cheeger in themself. An indicator function $\chi_{C}(x)$ of a convex and connected set $C$ with finite perimeter $Per(C)$, for which $C$ admits

[TABLE]

where $|C|$ denotes the area and $\kappa_{C}$ the curvature of $\partial C\in C^{1,1}$, is a solution of (4) with unit norm and therefore an eigenfunction of TV. With this geometric interpretation of eigenfunctions for TV as Cheeger sets, the role of perimeter and volume is significant for contrast invariant image decompositions. In the convex case, Duval and collaborators [17] proved that exact solutions of the $L^{1}-TV$ problem are given by an morphological opening followed by a simple test over the perimeter-area ratio. This fact was first published but not proved in [22]. For more complex shapes, formed by compositions of eigenfunctions, a better understanding of the $L^{1}$ scale-space flows is therefore a very promising direction.

2.2 $L^{2}-TV$ based Spectral Analysis

The nonlinear spectral analysis framework was first introduced by Gilboa in [3, 4] for the total variation regularization functional. A forward scale-space was constructed via the TV flow (3). Here, the first iterate is the original data $f$ which is then smoothed in every time-step such that increasingly fine scales are removed. This concept was later generalized to more general one-homogenous regularization functionals $J(u)$ in [5, 6, 7]. Moreover, it was shown that the forward scale-space can also be constructed via a variational approach by iteratively solving the ROF model (1) for $p=2$ and $J(u)=TV(u)$ using increasing regularization parameters $t_{\alpha}$. The idea of the spectral filtering approach is to transform the signal $f$ into a spectral domain, where both filtering of certain scales and decomposition of $f$ into significant signals is possible. Thus, eigenfunctions will be mapped onto a peak in the spectral domain. The shape of the eigenfunction is determined by the chosen regularization functional, for $J(u)=TV(u)$ the most prominent eigenfunction is a disc with radius $r$ surrounding the origin. In [23], Strong and Chan analyzed how the solution of the ROF model behaves for increasing $t_{\alpha}$ if $f$ is an eigenfunction. Thus, let $f(x)=c\cdot\mathbbm{1}_{B_{r}(0)}(x)$ be a disc of constant height $c$ and radius $r$ surrounding the origin and with a background of zero. Then the solution of the $L^{2}-TV$ model is given as:

[TABLE]

That means that even in the noise-free case, the reconstructed solution $u(t_{\alpha},x)$ never reaches the true value $f$ and for increasing regularization parameter $t_{\alpha}$ the disc height decreases. Solutions are corrupted by a systematic contrast loss that is dependent on the regularization strength but also on the radius $r$ and the height $c$ of input data $f$. To transform such eigenfunctions to a peak in the spectral domain, Gilboa defined the spectral transform function $\phi(t,x)$ and the spectral response function $S(t)$ as:

[TABLE]

The definition of $\phi$ allows, under certain conditions, that the original signal $f$ can be reconstructed via:

[TABLE]

where $\bar{f}$ is the mean of $f$. Filtered versions of $f$ can be constructed by applying:

[TABLE]

where $H(t)$ is the filter function.

However, a disadvantage of the $L^{2}$ based spectral framework is ambiguity with respect to size and intensity scales. The method is not able to clearly differentiate size and intensity scales since the timepoint $t_{d}$ at which a disc disappears and a peak occurs is determined by both values together.

2.3 $L^{1}-TV$ based Spectral Analysis

In the following section we want to combine the idea of nonlinear spectral TV analysis with the $L^{1}$ denoising model:

[TABLE]

As mentioned earlier, this model shows very interesting multi-scale decomposition behavior and seems therefore very suitable to be combined with the spectral approach. In [9], Chan and Esedoglu showed that the behavior of solutions of the $L^{1}-TV$ model if $f$ is the TV eigenfunction $f(x)=c\cdot\mathbbm{1}_{B_{r}(0)}(x)$ is significantly different from the $L^{2}$ case. Here, the solution is given by:

[TABLE]

An interesting observation is that the solution $u(t_{\alpha},x)$ is not dependent on the height $c$ of the disc but only on the radius $r$. That means that the $L^{1}-TV$ denoising approach is contrast-invariant and therefore highly suitable for decomposing a signal $f$ based on size scales. Note that the solution of (7) is not unique; for $f$ defined as above and $t_{\alpha}=\frac{r}{2}$ there exist, for example, an infinite number of solutions.

Since the gradient flow (2) is more challenging for $p=1$ rather than $p=2$ (due to non-smoothness of the data-term), we construct a forward scale-space by taking a variational approach; in other words, we solve (7) for increasing regularization parameters $t_{\alpha_{1}}<t_{\alpha_{2}}<...<t_{\alpha_{N}}$. To transform an eigenfunction of $TV(u)$ onto a peak in the spectral domain, a suitable spectral transformation function is now given by:

[TABLE]

Here, the first time derivative of $u(t,x)$ is in a distributional sense and the constant $\hat{c}$ is the median of $f$. Since the disc’s height does not decrease in every time step but remains the constant $c$ until it immediately decreases to 0, already the first time-derivative leads to a delta peak in the spectral transform function defined via:

[TABLE]

This definition was first introduced by Burger and collaborators in [6] and leads to an analogue to Parseval’s identity:

[TABLE]

The signal can be filtered based on size via:

[TABLE]

where $H(t)$ is again the filter function. A segmentation of objects in a certain size range (nearly independent of their intensity; only objects with intensity $\hat{c}$ cannot be found) is given by:

[TABLE]

Applications where this intensity-independent segmentation approach is especially helpful are described in section 4.

3 Numerical Approach

The following section introduces the numerical realization of the spectral $L^{1}-TV$ framework presented previously. The main component is the numerical solution of the denoising problem. To find solutions of this purely primal nonlinear minimization problem (7), we make use of the first order primal dual algorithm proposed by Chambolle and Pock [24]. As the name already suggest, the minimization scheme works with a primal dual version of (7) given by:

[TABLE]

where $P=\left\{g:\|g\|_{\infty}\leq 1\right\}$ and $\delta_{P}(g)$ equals [math] if $g\in P$, and equals $\infty$ if $g\notin P$. We define $K(u)=\nabla u$, $F(u)=\left\|u\right\|_{L^{1}}$ and $G(u)=\left\|u-f\right\|_{L^{1}}$.

The minimization algorithm proposed by Chambolle and Pock consists of three update steps: the first step is a dual update using the resolvent operator of $F^{\ast}$ and the second is a primal update using the resolvent operator of $G$. These are followed by a simple weighting step between the previous two primal iterates. See [24] for more details. The resolvent operators for $G$ and $F^{\ast}$ are presented in [24, chpt. 6.2]. To construct a forward scale space, we solve the $L^{1}-TV$ denoising model with increasing regularization parameter $t_{\alpha_{1}}<t_{\alpha_{2}}<...<t_{\alpha_{N}}$ and compute the spectral transform function $\phi(t,x)$ via backward-differences and the response function $S(t)$ based on these solutions. Both the resulting primal-dual algorithm to minimize (7) and the computation of the spectral functions are embodied in Algorithm 1. Note that to receive a very high degree of convergence we needed a very large number of iterations resulting in a long computational time. With fewer iterations, the algorithm did not converge completely and eigenfunctions lost contrast altough, according to (8), this should not be the case.

4 Results

In this section, we describe the main properties and advantages of the $L^{1}-TV$ based spectral approach introduced above. We discuss both several synthetic experiments and real experiments from biological cell imaging and retina imaging. For some examples, the results are compared to results from the $L^{2}$ based spectral framework to illustrate the differences between both models. For all results, we applied the framework summarized in Algorithm 1 and cluster the spectral response function afterwards. In our experiments, we set $\tau=0.2,\sigma=0.625,\theta=1$ and $maxIts=50.000$. For synthetic datasets we used $N=20$ linearly spaced $t_{\alpha_{i}}$ and increased $N$ to 50 for the experimental datsets. After solving the variational model with increasing regularization parameter, we manually clustered $S$ but comparable results can be achieved with common histogram thresholding methods such as Otsu’s method, the Triangle method or methods designed for more than two classes. The different classes in the histogram are always visualized by different colors. Reconstruction of the filtered signal was performed via (12) where $H(t)$ was defined as an indicator function of the filtered time interval and the resulting signal was color-coded with the same color as in the spectral response and multiplied with the original gray values so that intensity changes in the input signal remain detectable.

Detection of Varying Eigenshapes: Size versus intensity. In this set of experiments, we focus on circular objects that are all eigenfunction of the total variation functional. The aim of these experiments is to investigate which scales can be detected and reconstructed using our method and where the differences to an $L^{2}$ based approach are. Figure 3 shows an example with four discs in front of a uniform background, all showing the same intensity level. The corresponding spectral response function (b) clearly shows the four peaks each corresponding to one disc, as can be seen in the color-coded reconstruction in (c). Since $L^{1}$ is a purely size-based approach, it decomposes these eigenfunctions clearly. The opposite case is shown in Figure 3. Here, the input data is again composed of four discs but now with only one size and various intensity levels. Although an $L^{2}$ based approach would again show four peaks in the spectral response function due to its size/intensity ambiguity mentioned previously, we now see only one peak in the response (b). Filtering only the signal belonging to this red peak returns all four discs; see (c). This example clearly shows the contrast invariance of the $L^{1}$ data fidelity term which is a major difference between $p=2$ and $p=1$.

A direct comparison of both spectral frameworks is shown in Figures 4 and 5. Figure 4 shows as input data (a) two discs with different sizes and intensities. From a visual perspective, this would be clearly identified as two different scales but the $L^{2}-TV$ denoising approach is not able to distinguish between the objects. The scale that this approach uses is always a mixture of the size scale (small for the disc on top and large for the other one) with the intensity scale (large for the one on top and small for the second disc) and therefore it can occur that they both end up with the same “medium” scale. In the spectral domain, they are represented by one peak (b) and can therefore not be reconstructed separately. However, this is different for the contrast invariant $L^{1}-TV$ approach since it is purely size based. The spectral response function shows two easily separable peaks (d) that can be reconstructed one by one (e). A similar behavior can be seen for two discs on top of each other (see Figure 5(a)). Although these two discs are represented by two separate peaks in both approaches, the $L^{2}$ based approach mixes them if they are on top of each other. In this case, the peaks become less sparse and apart from each other (b) and in the reconstruction of the small scales (c) we see some artifacts of the larger disc (bluish ring around red disc). This is not the case for $L^{1}$. Both peaks are clearly separable (e) and the reconstruction gives a clear separation of both discs. The larger red circle in (f) is just a discretization artifact.

Segmentation of Experimental Cell Data. In Figure 6, we present a dataset that was also used in [18]. The experimental dataset (a) shows a fluorescent microscopy image of Circulating Tumor Cells. Here, the goal is to reliably segment all cells although they differ much in size and intensity. A multi-scale segmentation approach was presented in [18], see (d), but due to the intensity dependency of this approach, the method was not able to detect those cells that are very dim (highlighted with red boxes). Since our new $L^{1}$ based spectral approach ignores intensity differences, the method also finds the very dim cells. When reconstructing the orange part of the spectral response function (b) and using the thresholding formula (13), we obtain a segmentation that contains all four cells highlighted with a red box.

Experimental Data of Network Structures. In Figure 7 (a), a manually segmented blood vessel network taken from the STARE dataset [25] is shown. A problem often occurring in retinal blood vessel segmentation is that small vessels are also very dim and therefore even more challenging to detect. In (b), we added an intensity bias to the original data to test whether this influences our segmentation/clustering approach or not. In (c)-(f), the spectral response functions and reconstructions for $L^{2}$ denoising are shown. We see that there are fewer clear peaks in the spectral response function, especially in the case of intensity biased input. Another problem that we observe is that the vessels are not removed while retaining their original shape but are reshaped to more circular objects. This leads to a mixing of all scales in the reconstruction and the dim appearance of the visualization. For $L^{1}$ based denoising, the spectral responses are much sparser in (h) and (j) and therefore easier to cluster. The reconstructions (g) and (j) are also much easier to interpret since the vessels are not reshaped but removed in one step based on their diameter. Thus, for both networks (even for the one with an intensity bias) we see a clear clustering of blood vessels based on size. The difference between both methods is clearly shown in the magnified box. While in the top row the largest scale is just a turquoise shadow around the thickest blood vessels, the blood vessel is clearly reconstructed in the bottom row.

5 Conclusion and Outlook

In this paper, we have described the study of contrast-invariant $L^{1}$ data fidelities for variational multi-scale methods in combination with nonlinear spectral image analysis. We have shown that the contrast invariance results in an improved sparsity of spectral responses. In comparison to standard $L^{2}$, this allows a more informative spectral image representation to be obtained. We presented a model, an efficient algorithm and numerical results. In the particular case of experimental data sets that have complex shapes and strong intensity variations of objects or the background, our method outperforms the current standard method for nonlinear spectral decomposition. For future studies, it will be important to extend the ideas to nonlocal graph-based problems such as $L^{1}$ with nonlocal $TV$, and further analyze the relationship of the doubly nonlinear scale-space flow to the proposed scale-space procedure.

Acknowledgements

LZ, LT, CB acknowledge support by the EUFP7 program #305341 CTCTrap and the IMI EU program #115749 CANCER-ID. CB acknowledges support for his TT position by the NWO via Veni grant 613.009.032 within the NDNS+ cluster.