Multiview Attenuation Estimation and Correction

TL;DR

This paper introduces a novel method for estimating and correcting attenuation coefficients using multiview observations, applicable to devices like lidar and microscopes, based on solving a nonlinear inverse problem with new algorithms.

Contribution

It presents a new computational approach for attenuation estimation and correction, including robust solvers and efficient algorithms, applicable to various imaging modalities.

Findings

Validated on simulated data demonstrating effectiveness

Developed algorithms for attenuation correction in imaging devices

Proposed efficient methods for solving related nonlinear equations

Abstract

Measuring attenuation coefficients is a fundamental problem that can be solved with diverse techniques such as X-ray or optical tomography and lidar. We propose a novel approach based on the observation of a sample from a few different angles. This principle can be used in existing devices such as lidar or various types of fluorescence microscopes. It is based on the resolution of a nonlinear inverse problem. We propose a specific computational approach to solve it and show the well-foundedness of the approach on simulated data. Some of the tools developed are of independent interest. In particular we propose an efficient method to correct attenuation defects, new robust solvers for the lidar equation as well as new efficient algorithms to compute the Lambert W function and the proximal operator of the logsumexp function in dimension 2.

Equations

$$u_1(x)=\beta (x)\exp \left(-{\int }_0^x\alpha (t)dt\right)\text{mbox}forx\in [0,1]$$

$$u_2(x)=\beta (x)\exp \left(-{\int }_x^1\alpha (t)dt\right)\text{mbox}forx\in [0,1].$$

$$v(x)={\int }_0^x\alpha (t)dt-{\int }_x^1\alpha (t)dt.$$

$$\alpha (x)=\frac{1}{2}\frac{\partial }{{\partial }_x}v(x)$$

$$\beta (x)=\frac{u_1(x)}{\exp \left(-{\int }_0^x\alpha (t)dt\right)}.$$

$$u_1(x)=\mathcal{P}\left(\frac{C}{x^2}\beta (x)\exp \left(-2{\int }_0^x\alpha (t)dt\right)\right),$$

$$u_1(x,y,z)=\mathcal{P}\left(C\beta (x,y,z)\exp \left(-(A_t^{+}{\alpha }_t+A_e^{+}{\alpha }_e)(x,y,z)\right)\right)$$

$$u_2(x,y,z)=\mathcal{P}\left(C\beta (x,y,z)\exp \left(-(A_t^{-}{\alpha }_t+A_e^{-}{\alpha }_e)(x,y,z)\right)\right).$$

$$(A_1^t{\alpha }_t)(x,y,z)={\int }_0^z{\alpha }_t(x,y,t)dt$$

$$(A_2^t{\alpha }_t)(x,y,z)={\int }_1^z{\alpha }_t(x,y,t)dt$$

$$(A_1^e{\alpha }_e)(x,y,z)={\int }_0^z{\alpha }_e(x,y,t)dt$$

$$(A_2^e{\alpha }_e)(x,y,z)={\int }_1^z{\alpha }_e(x,y,t)dt.$$

$$u_1(x,y,z)=\mathcal{P}\left(\beta (x,y,z)\exp \left(-(A_1\alpha )(x,y,z)\right)\right)$$

$$u_2(x,y,z)=\mathcal{P}\left(\beta (x,y,z)\exp \left(-(A_2\alpha )(x,y,z)\right)\right),$$

$$u_{1,+}=\mathcal{P}\left(C\beta \exp (-A_{(1,t,+)}{\alpha }_t-A_{(1,e,+)}{\alpha }_e)\right)$$

$$u_{1,-}=\mathcal{P}\left(C\beta \exp (-A_{(1,t,-)}{\alpha }_t-A_{(1,e,-)}{\alpha }_e)\right)$$

$$u_{2,+}=\mathcal{P}\left(C\beta \exp (-A_{(2,t,+)}{\alpha }_t-A_{(2,e,+)}{\alpha }_e)\right)$$

$$u_{2,-}=\mathcal{P}\left(C\beta \exp (-A_{(2,t,-)}{\alpha }_t-A_{(2,e,-)}{\alpha }_e)\right).$$

$$u_i=\mathcal{P}\left(\beta \exp (-{\mathbf{A}}_i\mathbf{\alpha })\right),$$

$$\left\{\begin{array}[]{ll}u_{1}&=\mathcal{P}\left(\beta\exp(-A_{1}\alpha)\right)\\ u_{2}&=\mathcal{P}\left(\beta\exp(-A_{2}\alpha)\right).\end{array}\right.$$

$$(A_{1}u)[i]=\sum _{j=1}^i{u}_1[j].$$

$$A_1=\left(\begin{array}{cccccc}1& 0& 0& 0& \dots & 0\\ 1& 1& 0& 0& \dots & 0\\ 1& 1& 1& 0& \dots & 0\\ \dots & \dots & \dots & \ddots & \dots & 0\\ 1& 1& 1& 1& \dots & 1\end{array}\right)$$

$$\underset{\alpha \in {\mathbb{R}}^n,\beta \in {\mathbb{R}}^n}{\max }\mathrm{p}(\alpha ,\beta |u_1,u_2).$$

$$\underset{\alpha \in {\mathbb{R}}^n,\beta \in {\mathbb{R}}^n}{\min }-\log (\mathrm{p}(u_1,u_2|\alpha ,\beta ))-\log (\mathrm{p}(\alpha ,\beta )).$$

$${\lambda }_1=\beta \exp (-(A_1\alpha ))\text{mbox}and{\lambda }_2=\beta \exp (-(A_2\alpha )).$$

$$\mathbb{P}(X=k)=\frac{{\lambda }^k{e}^{-\lambda }}{k!},$$

$$\mathrm{p}(u_1,u_2|\alpha ,\beta )=\sum _{i=1}^n{\lambda }_1[i]+{\lambda }_2[i]-u_1[i]\log ({\lambda }_1[i])-u_2[i]\log ({\lambda }_2[i])+C,$$

$$\mathrm{p}(\alpha )\propto \exp (-R_{\alpha }(\alpha ))\text{mbox}and\mathrm{p}(\beta )\propto \exp (-R_{\beta }(\beta )),$$

$$R_{\alpha }(\alpha )=+\infty \text{ if }\exists i\in \{1,\dots ,n\},\alpha [i]<0$$

$$R_{\beta }(\beta )=+\infty \text{ if }\exists i\in \{1,\dots ,n\},\beta [i]<0.$$

Full text

1 Introduction

The ability to analyze the composition of gases in the atmosphere, the organization of a biological tissue, or the state of organs in the human body has invaluable scientific and societal repercussions. These seemingly unrelated issues can be solved thanks to a common principle: rays traveling through the sample are attenuated and this attenuation provides an indirect measurement of absorption coefficients. This is the basis of various devices such as X-ray and optical projection tomography [1, 2, 3] or lidar [4]. The aim of this paper is to provide an alternative approach based on the observation of the sample from a few different angles.

1.1 The basic principle

Let us provide a flavor of the proposed idea in an idealized 1D system. Assume that two measured signals $u_{1}$ and $u_{2}$ are formed according to the following model:

[TABLE]

and

[TABLE]

The function $\beta:[0,1]\to\mathbb{R}_{+}$ will be referred to as a density throughout the paper. It may represent different physical quantities such as backscatter coefficients in lidar or fluorophore densities in microscopy. The function $\alpha:[0,1]\to\mathbb{R}_{+}$ will be referred to as the attenuation and may represent absorption or extinction coefficients. The signals $u_{1}$ and $u_{2}$ can be interpreted as measurements of the same scene under opposite directions. Equations (1) and (2) coincide with the Beer-Lambert law that is a simple model to describe attenuation of light in absorbing media. The question tackled in this paper is: can we recover both $\alpha$ and $\beta$ from the knowledge of $u_{1}$ and $u_{2}$?

Under a positivity assumption $\beta(x)>0$ for all $x\in[0,1]$, the answer is straightforwardly positive. Setting $v(x)=\log\left(\frac{u_{2}(x)}{u_{1}(x)}\right)$, equations (1) and (2) yield:

[TABLE]

Therefore

[TABLE]

and

[TABLE]

Unfortunately, formulas (4) and (5) only have a theoretical interest: they cannot be used in practice since computing the derivative of a log of a ratio is extremely unstable from a numerical point of view. We will therefore design a numerical procedure based on a Bayesian estimator to retrieve the density $\alpha$ and attenuation coefficient $\beta$ in a stable and efficient manner. It is particularly relevant when the data suffer from Poisson noise.

1.2 Contributions

This paper contains various contributions listed below.

• •

We show that it is possible to retrieve attenuation coefficients from multiview measurements in different systems such as lidar, confocal or SPIM microscopes. To the best of our knowledge, this fact was only known in lidar until now [5, 6, 7]111We became aware of these works while finishing this work..

Fig. 1 summarizes the proposed idea. The attenuation, which is usually considered as a nuisance in confocal microscopy is exploited to measure attenuation. The algorithm successfully retrieves estimates of the density and attenuation from two attenuated and noisy images. Let us mention that some researchers already proposed to measure absorption and correct attenuation by combining optical projection tomography and SPIM imaging [8]. The principle outlined here shows that much simpler optical setups (a traditional confocal microscope) theoretically allows estimating the same quantities.

• •

We propose novel Bayesian estimators for the density $\alpha$ and the attenuation $\beta$. To the best of our knowledge, this is the first attempt to provide a clear statistical framework to recover these quantities.

• •

The proposed estimators are solutions of a nonconvex problem. We develop a specific nonlinear programming solver based on a warm-start convex initialization. This leads us to develop an efficient algorithm to compute the proximal operator of the logsumexp function in dimension 2 as well as a new algorithm to compute the Lambert W function.

• •

The proposed estimators also seem to be novel for the standard mono-view inverse problem in lidar and for correcting attenuation defects under a Poisson noise assumption.

• •

We perform a numerical validation of the proposed ideas on synthetic data, showing the well-foundedness of the approach. The validation of the method on specific devices is left as an outlook for future works.

2 Applications

In this section, we show various applications where the methodology proposed in this paper can be applied. We finish by describing precisely the discrete model considered in our numerical experiments.

2.1 Lidar

In lidar, an object (atmosphere, gas,…) is illuminated with a laser beam. Particles within the object reflect light. The time to return of the reflected light is then measured with a scanner. The received signal $u_{1}(x)$ is the backscattered mean power at altitude $x$ for a specific wavelength. The density $\beta$ corresponds to the backscattered coefficient, while $\alpha$ is called extinction coefficient. The equation relating $u_{1}$ to $\alpha$ and $\beta$ is:

[TABLE]

where $C$ is independent of $x$. The notation $\mathcal{P}(z)$ stands for a Poisson distributed random variable of parameter $z$. The Poisson distribution is a rather good noise model in lidar, since measurements describe a number of detected photons. The term $\frac{C}{x^{2}}\beta(x)$ appears in the lidar equation (6) instead of simply $\beta$. The algorithm developed later will allow retrieving $\frac{C}{x^{2}}\beta(x)$ instead of $\beta$. This is not a problem since there is a direct known relationship between both.

Remark 1**.**

In Raman lidar, the coefficient $\beta$ corresponds to the molecular density of the atmosphere, while $\alpha$ is the sum of extinction coefficients at different wavelengths. The theory developed herein also applied to this setting.

When the backscatter coefficient $\beta$ has a known analytical relationship with the extinction coefficient $\alpha$, direct inversion is possible. A popular method is Klett’s formula [9] for instance. Alternative formula exist [10] when the backscatter coefficient is known. The recent trend consists in using iterative methods coming from the field of inverse problems [11, 12, 13], leading to improved robustness. All these approaches crucially depend on a precise knowledge of the backscatter coefficient. This is a strong hypothesis that is often rough or unreasonable in practice.

To overcome this issue, a few authors proposed to use two opposite lidars and to retrieve the attenuation coefficient using equation (4) [5, 6, 7]. The stability to noise was ensured by linear filtering of the input and output data. The algorithms proposed later will provide a more robust and statistically motivated approach.

2.2 Confocal microscopy

The principle proposed herein can also be applied to fluorescence microscopy, in particular confocal and multi-SPIM microscopy. To the best of our knowledge, this idea is novel. We expect it to be relevant when absorption dominates scattering and diffusion. A particular case of interest should be optically cleared samples that are commonly used in optical projection tomography [2].

In confocal microscopy, a laser beam is focused at a specific point in 3D space and excites fluorophores. The emitted light is then collected on a camera. A 3D image can be created by scanning the whole sample volume with the focal spot, see Fig. 2.

In order to apply the proposed principle, two images from opposite sides have to be taken. This can be done by either rotating the sample or using $4\pi$ microscopes [14, 15]. The image formation model can then be written as follows:

[TABLE]

and

[TABLE]

The coefficients $\alpha_{t}$ and $\alpha_{e}$ refer to the attenuation coefficients for the transmitted and emitted light respectively. They may differ since the excitation light is typically red, while the emitted light is usually green. In order to retrieve them, we will assume a linear relationship of type $\alpha=\alpha_{e}=c\alpha_{t}$ between both for a certain constant $c$. The symbols $A_{t}^{+}$, $A_{e}^{+}$, $A_{t}^{-}$ and $A_{t}^{+}$ are integral operators describing the optical path for the transmitted and emitted light respectively. In the numerical experiments of this paper, we will simply use the following model:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The proposed algorithm also applies to more general linear operators integrating the attenuation coefficients along cones. With the mentioned hypotheses and setting $C=1$, equations (7) and (8) simplify to:

[TABLE]

and

[TABLE]

where $A_{1}=A^{t}_{1}+cA^{e}_{1}$ and $A^{2}=A^{t}_{2}+cA^{e}_{2}$.

2.3 Multiview SPIM

SPIM is an acronym for Selective Plane Illumination Microscopy [16]. An alternative name is Light Sheet Fluorescence Microscopy (LSFM). Its principle is explained in Fig. 3, left: a laser beam passing through a cylindrical lens focuses in only one direction, getting the shape of a light sheet. This light sheet excites fluorophores in a sample. The light sheet coincides with the focal plane of a microscope, allowing imaging 2D sections of the object. A 3D image can be reconstructed by shifting the sample or the light sheet and stacking the obtained images.

Many different variants of this microscope exist. Here, we are interested in the multiview SPIM [17, 18, 19, 20], which simultaneously produces 4 images generated with different optical paths, see Fig. 3, right and Fig. 4. Similarly to the confocal microscope, let $\alpha_{t}$ and $\alpha_{e}$ denote the attenuation maps for the transmitted and emitted light respectively. The four images can be written as:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The eight operators $A_{\cdot,\cdot,\cdot}$ describe integrals along the different directions for the emitted and transmitted waves. It is relatively easy to show, using a reasoning similar to that of paragraph 1.1, that the four images allow retrieving the attenuation coefficients $\alpha_{e}$ and $\alpha_{t}$, without assuming linear relationships between them. This is an advantage over the simpler confocal microscope.

2.4 Summary

We showed three different applications where measuring attenuation could be attacked with a similar methodology. In all mentioned applications, we have access to a set of $m$ signals $(u_{i})_{1\leq i\leq m}$ with the following expression:

[TABLE]

where $\boldsymbol{\alpha}$ is a concatenation of attenuation coefficients at different wavelengths and $\boldsymbol{A}$ is a concatenation of integral operators. In this paper, we will focus on the simplest setting where only two views are available. We developed the algorithms in such a way that their extension to an arbitrary number of views be rather straightforward. In addition, we will assume that the two views are opposite in the numerical experiments. This is a requirement in 1D, but not when dealing with 2D or 3D images.

3 MAP estimator and numerical evaluation

3.1 The discretized model

The discrete model considered in this paper reads:

[TABLE]

The signals $u_{1}$, $u_{2}$, $\beta$ and $\alpha$ are assumed to be nonnegative and belong $\mathbb{R}^{n}$, where $n=n_{1}\ldots n_{d}$ denotes the number of observations and $d$ is the space dimension. The value of a vector $u_{1}$ at location $i=(i_{1},\ldots,i_{d})$ will be denoted either $u_{1}[i]$ or $u_{1}[i_{1},\ldots,i_{d}]$. The matrices $A_{1}$ and $A_{2}$ in $\mathbb{R}^{n\times n}$ are discretization of linear integral operators. In our numerical experiments, the product $A_{1}u_{1}$ represents the cumulative sum of $u_{1}$ along one direction and the product $A_{2}u_{2}$ represents the cumulative sum of $u_{2}$ in the opposite direction. For instance, for a 1D signal, we set:

[TABLE]

Therefore, matrix $A_{1}$ has the following lower triangular shape:

[TABLE]

We are now ready to design a Bayesian estimator of $\alpha$ and $\beta$ from model (20).

3.2 A Bayesian estimator

The Maximum A Posteriori (MAP) estimators $\hat{\alpha}$ and $\hat{\beta}$ of $\alpha$ and $\beta$ are defined as the maximizers of the conditional probability density:

[TABLE]

By using the Bayes rule and a negative log-likelihood, this is equivalent to finding the minimizers of:

[TABLE]

Let us evaluate $\mathrm{p}(u_{1},u_{2}|\alpha,\beta)$. To this end, set

[TABLE]

Since the distribution of a Poisson distributed random variable with parameter $\lambda$ has the following probability mass function:

[TABLE]

we get:

[TABLE]

where $C$ is a value that does not depend on $\alpha$ and $\beta$. Next, we assume that $\alpha$ and $\beta$ are independent random vectors with probability distribution functions of type:

[TABLE]

where $R_{\alpha}:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ and $R_{\beta}:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ are regularizers describing properties of the density and attenuation maps. In particular, we choose them so as to impose nonnegativity of the estimators:

[TABLE]

and

[TABLE]

Overall, the optimization problem characterizing the MAP estimates reads:

[TABLE]

where

[TABLE]

Remark 2**.**

With an arbitrary number $m$ of views, it suffices to replace $\sum_{j=1}^{2}$ by $\sum_{j=1}^{m}$ in the above expression. For $m=1$ view, the problem $\min_{\alpha}F(\alpha,\beta)$ allows recovering the attenuation knowing the density: this is an inverse problem met in lidar. To the best of our knowledge, the proposed formulation is novel for this problem. The problem $\min_{\beta}F(\alpha,\beta)$ corresponds to correcting the attenuation on the density map. This is also a frequently met problem [21, 22, 23, 24] and the proposed approach also seems novel.

Remark 3**.**

The hypothesis of independence between $\alpha$ and $\beta$ can probably be improved in some applications. For instance, in lidar, it is well known that extinction and backscatter coefficients are strongly related. Similarly, strongly absorbing parts of biological specimens are likely to have a specific density of fluorophores. We prefer stating this independence property, since our aim is to derive a generic algorithm.

3.3 The overall optimization algorithm

If $R_{\alpha}$ and $R_{\beta}$ are convex functions, then $F$ is convex in each variable separately. Unfortunately, $F$ is usually nonconvex on the product space $\mathbb{R}_{+}^{n}\times\mathbb{R}_{+}^{n}$ since the 2D function $f(x,y)=\exp(-x)y+x-\log(y)$ is nonconvex. Finding global minimizers in reasonable times therefore seems complicated, except for specific regularizers $R_{\alpha}$ and $R_{\beta}$ that could compensate for the nonconvexity of the data term. The main observation in this paragraph is that it is possible to find a global minimizer of $F$ when $R_{\alpha}$ is a standard convex regularizer and $R_{\beta}$ is just the indicator function of the positive orthant. Set

[TABLE]

With this specific choice, the optimality conditions of problem (31) with respect to variable $\beta$ read:

[TABLE]

This expression can be seen as a simple estimator of $\beta$ knowing $u_{1}$, $u_{2}$ and $\alpha$. By replacing this expression in (32), we obtain an optimization problem in variable $\alpha$ only:

[TABLE]

Proposition 1**.**

Problem (35) is convex for a convex regularizer $R_{\alpha}$.

Proof.

The term $u_{j}[i](A_{j}\alpha)[i]$ is linear, hence convex. The term $\log\left(\sum_{j=1}^{2}\exp(-(A_{j}\alpha)[i])\right)$ is the composition of the convex logsumexp function with a linear operator, hence it is convex. ∎

The above observation motivates using the minimizer of (35) as an initial guess. Then, it is possible to use an arbitrary convex regularizer $R_{\beta}$ and to minimize (31) using an alternate minimization between $\alpha$ and $\beta$. This idea is captured in Algorithm 1, it ensures a monotonic decay of the cost function: $F(\alpha_{k+1},\beta_{k+1})\leq F(\alpha_{k},\beta_{k})$. However the algorithm does not necessarily converge to a stationary point. This is not a critical issue since we will see later that only 1 iteration is usually preferable.

3.4 Minimizing the initial convex program: warm start

We now delve into the numerical resolution of the warm start initialization (35). First, we need to choose a convex regularizer $R_{\alpha}$. In this paper, we propose to simply use the total variation [25], which is well known to preserve sharp edges. Its expression is given by:

[TABLE]

where $\nabla:\mathbb{R}^{n}\to\mathbb{R}^{dn}$ is a discretization of the gradient and $\lambda_{\alpha}\geq 0$ is a regularization parameter. We will use the standard discretization proposed in [26] in our numerical experiments.

Problem (35) is convex, but rather hard to minimize for various reasons listed below. First, the vectors $\alpha$ and $\beta$ may be very high dimensional, preventing the use of an arbitrary black-box method. Second, the regularizer $R_{\alpha}$ is non differentiable. Third, the operators $A_{i}$ have a spectral norm depending on the dimension $n$, preventing the use of gradient based methods since the Lipschitz constant of the gradient would be too high, see Proposition 2. Last, the proximal operator associated to the logsumexp function has no simple analytical formula.

Proposition 2**.**

Matrix $A_{1}$ in (22) satisfies $\|A_{1}\|_{2\to 2}\gtrsim n$, where $\|\cdot\|_{2\to 2}$ stands for the spectral norm.

Proof.

[TABLE]

∎

A large number of splitting methods have been developed to solve problems of type (35), and we refer to the excellent review paper [27] for an overview. Among them, the Simultaneous Direction Method of Multipliers (SDMM), a variant of the ADMM [28, 29] is particularily adapted to the structure of our problem. This algorithm allows solving problems of type:

[TABLE]

where functions $g_{i}:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ are convex closed and the operators $L_{i}:\mathbb{R}^{n}\to\mathbb{R}^{m_{i}}$ are linear and such that $Q=\sum_{i=1}^{n}L_{i}^{T}L_{i}$ is an invertible matrix. The SDMM then takes the algorithmic form described in Algorithm 2.

To cast problem (35) into form (37), we use the following choices. We set

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The numbers $c_{1},c_{2},c_{3}$ are positive constants allowing to accelerate the algorithm’s convergence by balancing the relative importance of each term. This can also be seen as a simple diagonal preconditioner. In our numerical experiments, we set $c_{1}=1$ and tune $c_{2}$ and $c_{3}$ manually to accelerate convergence.

In order to apply Algorithm 2, we need to compute the proximal operators of each function $g_{i}$, defined by:

[TABLE]

The proximal operators of $g_{2}$ and $g_{3}$ have closed form solutions found in nearly all recent total variation minimization solvers. We refer to [29] for instance. Unfortunately, the proximal operator of $g_{1}$ has no closed-form expression. In order to compute it, we propose using a non trivial Newton based algorithm described in section 6.1. Finally, we need to evaluate matrix-vector products with $Q^{-1}$. This can be achieved using either a LU factorization or a conjugate gradient. In our codes, we simply use a conjugate gradient algorithm.

To conclude this paragraph, we illustrate the results obtained by the described procedure in Figure 5.

3.5 Recovering the attenuation knowing the density

In this paragraph, we focus on the resolution of:

[TABLE]

where $L_{1}$, $L_{2}$ and $L_{3}$ are defined as in the previous section, $h_{2}=g_{2}$, $h_{3}=g_{3}$ and

[TABLE]

As in the previous section, the SDMM is a good candidate to solve this problem. The only difficulty is to compute the proximal operator $\mathrm{prox}_{\gamma h_{1}}$. This amounts to solving $2n$ 1D optimization problems of the following type:

[TABLE]

where $u$, $z_{0}$ and $\beta$ are real numbers.

Proposition 3**.**

The solution $z^{*}$ of problem (44) is given by:

[TABLE]

where $a=z_{0}-\gamma u$ and $W$ is the so-called Lambert W function, i.e. the reciprocal function of $x\mapsto x\exp(x)$.

Proof.

The optimality conditions for this problem read:

[TABLE]

so that

[TABLE]

with $a=z_{0}-\gamma u$. Now, we can write $z^{*}=a+w$, hence $w=\beta\gamma\exp(-a-w)$, which is still equivalent to:

[TABLE]

or $w=W(\beta\gamma\exp(-a))$. ∎

Remark 4**.**

The Lambert W function has been studied thoroughly. The standard approach to compute it consists of using Halley’s method (a Househölder method with cubic convergence rate) initialized with the first terms of an asymptotic expansion. This method was proposed in the excellent review paper [30] and seems to be the default solver in MAPLE and MATLAB. Unfortunately, this method fails for our problem since the number $\beta\gamma\exp(-a)$ can exceed the maximum number available in double precision. We therefore develop a specific method in Appendix 6.2.

3.6 Recovering the density knowing the attenuation

In this paragraph, we focus on the resolution of $\min_{\beta\in\mathbb{R}^{n}}F(\alpha_{k},\beta)$. This amounts to simultaneously correcting the attenuation and denoising the resulting image. This is a rather simple inverse problem, but it seems original due to the noise statistics. A Poisson distributed variable multiplied by a positive constant different from 1 is not Poisson anymore. This makes the proposed algorithm similar, but different from existing approaches developed for Poisson noise in [31, 32] for instance.

A simple idea to regularize the problem is to use the total variation again, i.e. to set $R_{\beta}(\beta)=\lambda_{\beta}\sum_{i=1}^{n}\|(\nabla\beta)[i]\|_{2}$. Once again the resulting problem can be solved with the SDMM. Let us detail this procedure. Define

[TABLE]

The problem then reads:

[TABLE]

with $L_{1}=c_{1}I_{n}$,

[TABLE]

$L_{2}=c_{2}\nabla$ and

[TABLE]

The proximal operators of $j_{2}$ is standard and we do not detail it here. The proximal operator of $j_{1}$ is provided below:

Proposition 4**.**

We have:

[TABLE]

Proof.

It suffices to write the first order optimality conditions of $\min_{z\geq 0}1/2\|z-z_{0}\|_{2}^{2}+a/c_{1}z-u\log(z/c_{1})$. This shows that $z$ is the root of a second order polynomial. Its only positive root is given in (52). ∎

We show a typical result of total variation minimization in Figure (6). Pamareter $\lambda_{\beta}$ was chosen manually so as to maximize the SNR of the result.

4 Additional comments

4.1 Parameter selection

Data terms

The two data term parameters are $\lambda_{\alpha}$ and $\lambda_{\beta}$. They specify the regularity of the attenuation and the density respectively. In all our experiments, we optimized them by trial and error. We observed experimentally, that similar results are obtained within a relatively large range, making a manual optimization quite easy. In addition, for a given measurement device, the same parameter is likely to be always the same, decreasing the interest of an automatized procedure such as SURE.

Algorithms parameters

The optimization algorithms are based on the SDMM and their convergence rates depend a lot on the paramaters $\gamma$, $c_{1}$, $c_{2}$, $c_{3}$ and $c_{4}$. They may converge to a satisfactory solution rapidly (about 50 iterations) or slowly (more than 10000 iterations) depending on these choices. Unfortunately, we found no systematic method to choose them and also used a trial and error strategy in our numerical experiments. Our numerical experiments suggest that these parameters are suitable for a wide range of data (image size, maximum attenuation, image dynamics), so that the tedious tuning can be done once for all for a given application.

4.2 Number of outer iterations

Algorithm 1 depends on a number of outer iterations denoted $Nit$. One may wonder how many iterations are needed to obtain a satisfactory solution. It turns out that 1 iteration leads to the best results in terms of SNR: iterating more tends to degrade the solution, even though the cost function continues to decrease. The precise reason behind this is still unknown, but it is likely that this is total variation bias. Low contrasted solutions are favored at the expense of SNR. This seems to be one more example, where the Maximum A Posteriori principle should be taken with caution [33, 34, 35]. It is useful to get a rough idea of a functional to minimize (here it allowed designing the warm start estimate), but should not be considered as the best possible estimator.

This remark being stated, the proposed algorithm can be seen as a simple two step procedure:

• •

Find the warm start estimate (35).

• •

Correct the density (see section 3.6).

It is possible to evaluate the attenuation once again (see section 3.5), but we observed that this influenced the result very little.

4.3 Computing times

All the experiments of the paper were performed on a laptop with an Intel i7 processor with 4-cores. The codes were written mostly in Matlab (natively parallel), with some parts written in C with OpenMP support.

The complexity of the proposed algorithms scale roughly linearly with the number of pixels $n$, as shown in Fig. 7. We observed that the number of iterations of the SDMM to reach a given relative accuracy remains the same whatever the size $n$, while the cost per iteration scales linearly with it (at least for the cumulative sum integral operators considered herein).

As can be seen on Fig. 7 the algorithm takes around 48 seconds for a $256\times 256$ image. Out of these, 45 seconds are spent to recover the attenuation, while the 3 remaining are dedicated to correct the density.

All codes can be easily parallelized on a GPU. A speed-up of 100 can be expected on such an architecture, making the proposed methods suitable for large 2D or 3D images.

4.4 Influence of attenuation and signal-to-noise ratio

Two parameters strongly influence the ability to recover the attenuation and density: the signals dynamics (or signal-to-noise ratio) and the attenuation dynamics.

As the signal-to-noise ratio decreases, it becomes impossible to recover fine details. For instance, the fine stripes are not recovered in Fig. 1(e), but they are recovered for signals with a much higher amplitude. In Fig. 8, it can indeed be verified that a high dynamics of $10^{5}$ allows recovering most of the stripes. This experiment shows that highly sensitive EMCCD cameras should be preferred over more standard devices for this specific application.

The attenuation amplitude also plays a key role: if it is too low, then no attenuation can be detected. On the contrary, if it is too high, then the signals $u_{1}$ and $u_{2}$ will vanish too rapidly, making it impossible to evaluate the attenuation. This is illustrated in Fig. 8. It is remarkable that the algorithm manages to recover the attenuation partially for very low signal-to-noise ratio. In Fig. 8(c), we observe that the attenuation is partially recovered with no more than 30 expected photons per pixel!

4.5 Toolbox

A Matlab toolbox containing all the main algorithms described in the paper is provided on the website of the authors https://www.math.univ-toulouse.fr/~weiss/ and on GitHub https://github.com/pierre-weiss/MAEC. The Lambert W function and the proximal operator of logsumexp have been implemented with C-mex files with OpenMP support for multicore acceleration. Demonstration scripts are available for testing.

5 Conclusion & outlook

We proposed a robust and efficient approach to recover attenuation and correct density from multiview measurements. This principle was already known in the field of lidar and solved with simple filtering approaches. The algorithms proposed herein are based on a clear and versatile statistical framework. The approach seems promising in various devices such as lidar or some fluorescence microscopes. It is likely that its scope be much wider and we therefore provide a free Matlab toolbox on our website.

As a prospective, we plan to confront our algorithms with real data coming from lidar and microscopy. The total variation based algorithm to correct attenuation defects is somewhat disappointing since it is unable to recover fine textures. A promising and unexplored issue is to extend the nonlocal means algorithms [36] to solve this problem. To conclude, let us mention a serious limitation of the proposed approach: it is not so common to find a couple optical system-sample, where attenuation dominates scattering. We do not know at the present time how many applications can reasonably be modeled by equation (19). This question is central to precisely understand the strengths and limits of the proposed approach.

6 Appendices

6.1 Proximal operator of logsumexp in dimension 2

In this section, we propose a fast and accurate numerical algorithm based on Newton’s method to solve the following problem:

[TABLE]

where $z=\begin{pmatrix}z_{1}\\ z_{2}\end{pmatrix}$ and $x=\begin{pmatrix}x_{1}\\ x_{2}\end{pmatrix}$ are vectors in $\mathbb{R}^{2n}$. This problem may seem innocuous at first sight, but turns out to be quite a numerical challenge. The first observation is that it can be decomposed as $n$ independent problems of dimension 2 since:

[TABLE]

To simplify the notation, we will skip the index $i$ in what follows. The following proposition shows that our problem is equivalent to finding the proximal operator associated to the “logsumexp” function.

Proposition 5**.**

Define the logsumexp function $\mathrm{lse}(x_{1},x_{2})=\log\left(\sum_{j=1}^{2}\exp(x_{j})\right)$. The solution of problem (53) coincides with the opposite of the proximal operator of $\mathrm{lse}$:

[TABLE]

where $a=\gamma(u_{1}+u_{2})$ and $y_{j}=\gamma u_{j}-z_{j}$.

Proof.

The first order optimality conditions for problem (53) read

[TABLE]

By letting $a=\gamma(u_{1}+u_{2})$ and $y_{j}=\gamma u_{j}-z_{j}$, this equation becomes

[TABLE]

It now suffices to make the change of variable $x^{\prime}_{i}=-x_{i}$ to retrieve the optimality conditions of problem (55)

[TABLE]

∎

Remark 5**.**

To the best of our knowledge, this is the first attempt to find a fast algorithm to evaluate the prox of logsumexp. This function is important in many regards. In particular, it is a $C^{\infty}$ approximation of the maximum value of a vector. In addition, its Fenchel conjugate coincides with the Shannon entropy restricted to the unit simplex. We refer to [37, §3.2] for some details. The algorithm that follows has potential applications outside the scope of this paper.

We now design a fast and accurate minimization algorithm for problem (55) or equivalently, a root finding algorithm for problem (58). This algorithm differs depending on whether $y_{1}\geq y_{2}$ or $y_{2}\geq y_{1}$. We focus on the case $y_{1}\geq y_{2}$. The case $y_{2}\geq y_{1}$ can be handled by symmetry.

Let $\lambda=\frac{\exp(x^{\prime}_{1})}{\exp(x^{\prime}_{1})+\exp(x^{\prime}_{2})}$ and notice that $\frac{\exp(x^{\prime}_{2})}{\exp(x^{\prime}_{1})+\exp(x^{\prime}_{2})}=1-\lambda$. Therefore (58) becomes:

[TABLE]

Hence

[TABLE]

Taking the logarithm on each side yields 222Applying the logarithm is important for numerical purposes. When $y_{2}-y_{1}-a$ is very small, the exponential cannot be computed accurately in double precision.:

[TABLE]

We are now facing the problem of finding the root $\lambda^{*}$ of the following function:

[TABLE]

There are two important advantages for this approach compared to the direct resolution of (58). First, we have to solve a 1D problem instead of a 2D problem. More importantly, we directly constrain $x^{\prime}$ to be of form $x^{\prime}=y-a\delta$, where $\delta$ lives on the 2D simplex.

Let us collect a few properties of function $f$. First, we have:

[TABLE]

Therefore, $f$ is increasing on $(0,1)$. To use convergence results of Newton’s algorithm, we need to compute $f^{\prime\prime}$ as well:

[TABLE]

Proposition 6**.**

If $y_{1}\geq y_{2}$, then $x_{1}^{\prime}\geq x_{2}^{\prime}$ and

[TABLE]

Proof.

The first statement can be proven by contradiction. Assume that $x_{2}^{\prime}>x_{1}^{\prime}$, then equation (58) indicates that $y_{2}>y_{1}$.

For the second statement, it suffices to evaluate $f$ at the extremities of the interval since $f^{\prime}>0$. We get $f(1/2)=y_{2}-y_{1}\leq 0$ and $f\left(\frac{1}{1+\exp(y_{2}-y_{1})}\right)=-a+\frac{2a}{1+\exp(y_{2}-y_{1})}\geq 0$. ∎

Proposition 7**.**

Set $\lambda_{0}=\frac{1}{1+\exp(y_{2}-y_{1})}$. Then, the following Newton’s method

[TABLE]

converges to the root $\lambda^{*}$ of $f$, with a locally quadratic rate.

Proof.

First notice that $f^{\prime\prime}(\lambda)\geq 0$ on the interval $[1/2,1)$. Hence $f^{\prime\prime}$ is also positive on $I=[\lambda^{*},\lambda_{0}]$. This ensures that

[TABLE]

We prove this assertion by recurrence. Notice that $\lambda_{0}\geq\lambda^{*}$ by Proposition 6. Now, assume that $\lambda_{k}\geq\lambda^{*}$, then

[TABLE]

Hence, $\lambda_{k}-\lambda^{*}\geq\frac{f(\lambda_{k})}{f^{\prime}(\lambda_{k})}$ and $\lambda_{k+1}\geq\lambda^{*}$. In addition $\frac{f(\lambda_{k})}{f^{\prime}(\lambda_{k})}\geq 0$ on $I$, so that $\lambda_{k+1}\geq\lambda_{k}$.

The sequence $(\lambda_{k})_{k\in\mathbb{N}}$ is monotonically decreasing and bounded below, therefore it converges to some value $\lambda^{\prime}\geq\lambda^{*}$. Necessarily $\lambda^{\prime}=\lambda^{*}$, since for $\lambda^{\prime}>\lambda^{*}$, $\frac{f(\lambda^{\prime})}{f^{\prime}(\lambda^{\prime})}>0$.

To prove the locally quadratic convergence rate, we just invoke the celebrated Newton-Kantorovich’s theorem [38, 39], that ensures local quadratic convergence if $f^{\prime\prime}$ is bounded in a neighborhood of the minimizer. ∎

Finally, let us mention that computing $\lambda_{0}$ on a computer is a tricky due to underflow problems: in double precision the command $1+\exp(y_{2}-y_{1})$ will return $1$ for $y_{2}-y_{1}<-37\simeq\log(10^{-16})$. This may cause the algorithm to fail since $f$ and its derivatives are undefined at $\lambda=1$. In practice we therefore set $\lambda_{0}=1/(1+\exp(y_{2}-y_{1}))-10^{-16}$. Similarly, by bound (65), we get $\lambda^{*}=1$ up to machine precision whenever $y_{2}-y_{2}-a<\log(10^{-16})$. Algorithm 3 summarizes all the ideas described in this paragraph.

An attentive reader may have remarked that the convergence of Newton’s algorithm depends only on the difference $y(1)-y(2)$ and $a$. A shift of $y(1)$ and $y(2)$ by the same value does not change Newton’s iteration. In Fig. 9, we show that the algorithm behaves very well for a wide range of parameters. For $y(1)-y(2)$ and $a$ varying in the interval $[2^{-10},2^{20}]$, the algorithm never requires more than $18$ iterations to reach machine precision and needs $2.8$ iterations in average.

6.2 Computing the Lambert W function for large numbers

In this section, we propose a numerical algorithm to solve the following equation:

[TABLE]

The solution is given by $x^{*}=W(\exp(z))$, where $W$ is the Lambert W function.

When $z$ is sufficiently small, Halley’s method proposed in [30] can be applied. However, this method fails whenever $z$ is too large due to numerical instabilities. For instance, the command $\log(\exp(z))$ returns $+\infty$ for $z>710$. In that case, it makes sense to apply a logarithm on each side of equation (69) and solve the following problem instead:

[TABLE]

The two first terms of the asymptotic expansion of $W$ yield $W(\exp(z))\simeq z-\log(z)$ for large $z$. This value is a good starting point for a root finding algorithm. Here, we simply propose to use Newton’s algorithm for large $z$. The complete method is detailed in Algorithm 4. It never takes more than 5 iterations to reach machine precision. The initialization $w=\sqrt{5.43z+2}-1$ is obtained using a Padé approximation and the threshold $0.12$ has been determined experimentally in [30].

Acknowledgment

V. Debarnot is supported by Plan CANCER, MIMMOSA project. The authors wish to thank Emilio Gualda, Jan Huisken, Philipp Keller, Théo Liu, Jürgen Mayer and Anne Sentenac for interesting discussions and feedbacks on the model.