Non-convex non-local flows for saliency detection

Image processing by variational methods is an established field in applied mathematics and computer vision which aims to model typical low and mid level image reconstruction and restoration processes such as denoising, deconvolution, impainting, segmentation and registration of digital degraded images. Since the mathematical approach of Tikhonov and Arsenin on ill-posed inverse problems [3] and the applied work of Rudin, Osher and Fatemi [21] through their celebrated ROF model, there have been a standing effort to deal with new image processing tasks through variational methods. Currently there is a growing interest in image processing and computer vision applications for visual saliency models, able to focus on perceptually relevant information within digital images.

Semantic segmentation, object detection, object proposals, image clustering, retrieval and cognitive saliency applications, such as image captioning and high-level image understanding, are just a few examples of saliency based models. Saliency is also of interest as a means for improving computational efficiency and for increasing robustness in clustering and thresholding.

Despite of the lack of a general consensus on a proper mathematical definition of saliency, it has a clear biologically perceptive meaning: it models the mechanism of human attention, that is, of finding relevant objects within an image.

Recently, there has been a burst of research on saliency due to its wide application in leading medical disciplines such as neuroscience and cardiology. In fact, when considering medical images like those acquired in magnetic resonance imaging (MRI) or positron emission tomography (PET), the automatic obtention of saliency maps is useful for pathology detection, disease classification [22], location and segmentation of brain strokes, gliomas, myocardium detection for PET images, tumors quantification in FLAIR MRI [23], etc.

The role of saliency in models is application dependent, and several different techniques and approaches have been introduced to construct saliency maps. They vary from low dimensional manifold features minimization [27] to non-local sparse minimization [25], graphs techniques [9], partial differential equations (PDE) [12], superpixel [14], learning [13], or neural networks based approaches [4] (MIT-Benchmark).

With the aim to explore the applications and algorithms of non-smooth, non-local, non-convex optimization of saliency models, we present in this article a new variational model applied to Fluid Attenuated Inversion Recovery (FLAIR) MR images for accurate location of tumor and edema.

Models for saliency detection try to transform a given image, $f$, defined in the pixel domain, $\Omega$, into a constant-wise image, $u$, whose level sets correspond to salient regions of the original image. They are usually formulated through the inter-relation among three energies: fidelity, regularization, and saliency, being the latter the mechanism promoting the classification of pixels into two or more classes. There is a general agreement in considering the fidelity term as determined by the $L^{2}$ norm, that is

so that departure from the original state is penalized in the minimization procedure.

For regularization, an edge preserving energy should be preferred. The use of the Total Variation energy,

defined on the space of Bounded Variation functions has been ubiquitous since its introduction as an edge preserving restoration model [21]. Indeed, the Total Variation energy allows discontinuous functions to be solutions of the corresponding minimization problem, with discontinuities representing edges. This is in contrast to Sobolev norms, which enforce continuity across level lines and thus introduce image blurring. Only recently, the non-local version of the Total Variation energy and, in general, of the energy associated to the $p$-Laplacian, for $p>1$, has been considered in restoration modeling. In saliency modeling, only the range $p>2$ seems to have been treated [12]. One of the main advantages of introducing these non-local energies is the lack of the hard regularizing effect influencing their local counterparts.

For the saliency term, a phase-transition Ginzburg-Landau model consisting in a double-well absorption-reaction term is often found in the literature [8]. The resulting energy is a functional of the type

whose minimization drives the solution towards the discrete set of values $\{-1,1\}$, facilitating in this way the labeling process. However, due to the vanishing slope of $g(u)=(1-|u|^{2})^{2}$ at $u=\pm 1$, the resulting algorithm has a slow convergence to the minimizer.

In this article, we introduce new regularizing and saliency terms that enhance the convergence of the classification algorithm while keeping a good quality compromise. As regularizing term, we propose the non-local energy associated to the $p$-Laplacian for any $p>0$, see Eq. (3). Although there is a lack of a sound mathematical theory for the concave range $0<p<1$, computational evidence of the ability of these energies to produce sparse gradient solutions has been shown [11]. In any case, our analysis also includes the energies arising for $p\geq 1$, enjoying a well established mathematical theory.

With respect to the saliency term, we consider the combination of two effects. One is captured by the concave energy

with $\delta>0$ constant, which fastly drives the minimization procedure so that $H(u)\to-\infty$. To counteract this tendency and remain in the meaningful interval $u\in I=[0,1]$, being $\{0,1\}$ the saliency labels employed in this article, we introduce an obstacle which penalizes the minimization when the solution lies outside $I$. The modeling of such obstacle is given in terms of the indicator function

and the resulting saliency term is then defined as a weighted sum of the operators $H(u)$ and $\mathcal{I}_{I}(u)$.

Finally, the Euler-Lagrange equation for the addition of these three energies leads to a formulation of a multi-valued non-local reaction-diffusion problem which is later rendered to a single-valued equation by means of the Yosida’s approximants. The introduction of a gradient descent and the discretization of the corresponding evolution equation are the final tools we use to produce our saliency detection algorithm.

The main contributions contained in this paper may be summarized as follows.

• •

  The non-local $p$-Laplacian convex model proposed in [12], valid for $p>2$, is extended to include the range $0<p\leq 2$. Notably, the non-differentiable non-convex fluxes related to the range $0<p<1$ are investigated.

• •

  An absorption-reaction term of convex-concave energy nature is considered for fast saliency detection. In absence of a limiting mechanism, the reactive part could drive the solution to take values outside the image range $[0,1]$. Thus, a penalty term is introduced to keep the solution within such range.

• •

  The resulting multivalued constrained variational formulation is approximated and shown to be stable.

• •

  A fast algorithm overcoming the computational drawbacks of non-local methods is presented.

• •

  A 3D generalization of the model is considered providing promising results specially valuable in medical image processing and tumor detection.

The paper is organized as follows. In Section 2, we introduce the variational mathematical framework of our models. Starting with the local equations as guide for the modelling exercise, we focus on the non-local diffusive terms, explicited in the form of $p$-Laplacian, for $p>1$ and extend it to the range $0<p\leq 1$ through a differentiable family of fluxes to cover the resulting non-local non-convex hyper-Laplacian operators. Then, we introduce a multi-valued concave saliency detection term which defines an obstacle problem for the non-local diffusion models. In Section 3 we deduce the corresponding Euler-Lagrange equations. A gradient descent approximation is finally used to solve the elliptic non-local problems until stabilization of the associated evolution problems. In Section 4, we give the discretization schemes used for the actual computation of the solutions of the non-local diffusion problems. In Section 5 a simplified computational approach is described in order to reduce the time execution with a view to a massive implementation on the proposed data-set. Section 6 contains the numerical experiments on the proposed model and present the simulations performed on FLAIR sequences of MR images obtained from the BRATS2015 dataset [16]. Finally, in Section 7, we give our conclusions.

Gathering the fidelity, the regularizing and the saliency energies, we define a bilateral constrained obstacle problem associated to the following energy

where $\alpha>0$ is a parameter modulating the relationship between regularization and saliency promotion. Observe that there is no use in multiplying $\mathcal{I}_{I}(u)$ by the constant $1/\alpha$, so we omit it for clarity.

The Euler-Lagrange equation corresponding to (5) together with the use of a gradient descent method leads to the consideration of non-local multi-valued evolution problems.

For the sake of presentation, we have introduced the following notation in (6): we rewrote the Fréchet differential of $H(u)/\alpha+\lambda F(u)$ as $au(x)-b(x)$, with

and defined the non-local hyper-Laplacian ($0<p<1$) and $1$-Laplacian ($p=1$) diffusion operators

with differential kernels

Notice that, while for the local diffusion problem we must explicitly impose the homogeneous Neumann boundary conditions, which are the most common boundary conditions for image processing tasks, for the non-local diffusion problem this is no longer necessary since these conditions are implicitly imposed by the non-local diffusion operator [2].

In the previous section we have shown that solutions of the multivalued problem (6) may be approximated by the introduction of Yosida’s approximants, leading to the single-valued problems $P_{\epsilon,r}(u)$ in (10) that depend on the Yosida’s approximation parameter $r$. We also proved that in the limit $r\to 0$ the corresponding solutions lie in the relevant range of values for image processing tasks, this is, in the interval $[0,1]$.

In this section we provide a fully discrete algorithm to numerically approximate the limit solution of (10) when $r\to 0$. First, we introduce a time semi-implicit Euler discretization of the evolution equation in (10), that we show to retain the stability property of its continuous counterpart, stated in Theorem 1.

The resulting space dependent PDE is discretized by finite differences. Since the problem is nonlinear and, in addition, we want to pass to the limit $r\to 0$, we introduce an iterative algorithm which renders the problem to a linear form and, at the same time, replaces the fixed parameter $r$ by a decreasing sequence $r_{j}\to 0$.

In the previous sections we have deduced, through a series of approximations, a discrete algorithm to compute approximated solutions of the obstacle problem $P_{\epsilon}(u)$ in (6). We have shown that our scheme is stable with respect to the approximating parameter $r$, producing solutions of problems $P_{r}(u)$ that, in the limit $r\to 0$, lie effectively in the image value range $[0,1]$, apart from producing the required edge preserving saliency detection on images.

In this section, by introducing some hard nonlinearities (truncations) to replace one of the iterative loops of (19), we provide a simplified algorithm for solving a problem closely related to (6). In addition, we use an approximation technique, based on the discretization of the image range, to compute the non-local diffusion term. These modifications allow for a fast computation of what we demonstrate to be fair approximations to the solutions of the original problem (6).

Considering the time discrete problem (15), we introduce two changes which greatly alleviate the computational burden:

1. 1.

   Compute the non-local diffusion term fully explicitly, and

2. 2.

   Replace the obstacle term by a hard truncation.

Thus, we replace problem $P_{\epsilon,r}(u^{n})$ in (15) by the following which can be deduced from problem $P_{\epsilon}(u)$ in (6) using the two above strategies.

Set $u^{0}=f$. For $n=0,\ldots,N-1$, and for $k=0,\ldots,LM-1$, compute

and truncate $u^{n+1}[k]$

We shall show that there are very small differences between the solutions of the explicit truncated problem $P_{\epsilon,0}(u^{n})$ and the solutions of the $P_{\epsilon,r}(u^{n})$ problems for $r$ sufficiently small. Nevertheless, the numerical scheme is greatly improved and much more efficient because costly iteration in $r$-loop is avoid as it can be seen in section (6). We finally describe the efficient approach of [26] (see also [7], and and [6] for a related approach) that we use for computing the sum in (21), corresponding to the non-local diffusion term, by discretizing also the range of image values. Let $q=\{q_{1},\ldots,q_{Q}\}$ be a quantization partition, with $0=q_{1}<q_{2}<\ldots<q_{Q-1}<q_{Q}=1$, where $Q$ is the number of quantization levels. Let $v:\Omega\to[0,1]$ be a quantized function, that is, taking values on $q$. For each $i=1,\ldots,Q$, we introduce the discrete convolution operator

where we recall that $w[k,m]=w(x_{k}-x_{m})$. We then have

Notice that, for any $v$ taking values in $q$, the computation of each $K_{\epsilon,p}^{i}(v)$ may be carried out in parallel by means of fast convolution algorithms, e.g. the fast Fourier transform.