Deep learning versus $\ell^1$-minimization for compressed sensing photoacoustic tomography

TL;DR

This paper compares deep learning methods and $ ext{l}^1$-minimization for compressed sensing in photoacoustic tomography, showing deep learning excels with deterministic measurements while $ ext{l}^1$-minimization is better with Bernoulli measurements.

Contribution

It provides a comparative analysis of deep learning and $ ext{l}^1$-minimization techniques for CS in PAT, highlighting their relative performance under different measurement schemes.

Findings

Deep learning yields more accurate results with deterministic measurements.

$ ext{l}^1$-minimization outperforms deep learning with Bernoulli measurements.

Nullspace network outperforms residual network in MSE.

Abstract

We investigate compressed sensing (CS) techniques for reducing the number of measurements in photoacoustic tomography (PAT). High resolution imaging from CS data requires particular image reconstruction algorithms. The most established reconstruction techniques for that purpose use sparsity and $\ell^1$-minimization. Recently, deep learning appeared as a new paradigm for CS and other inverse problems. In this paper, we compare a recently invented joint $\ell^1$-minimization algorithm with two deep learning methods, namely a residual network and an approximate nullspace network. We present numerical results showing that all developed techniques perform well for deterministic sparse measurements as well as for random Bernoulli measurements. For the deterministic sampling, deep learning shows more accurate results, whereas for Bernoulli measurements the $\ell^1$-minimization algorithm performs best. Comparing the implemented deep learning approaches, we show that the nullspace network uniformly outperforms the residual network in terms of the mean squared error (MSE).

Figures

Tables

Table 1. Table 1: Performance averaged over 50 blood vessel images.

	SME	PSNR	SSIM
Sparse measurements
FBP	$15.1\text{×}{10}^{-4}$	28.6	$4.83\text{×}{10}^{-1}$
${\ell }^1$-minimization	$3.35\text{×}{10}^{-4}$	35.0	$8.50\text{×}{10}^{-1}$
residual network	$3.08\text{×}{10}^{-4}$	35.6	$9.30\text{×}{10}^{-1}$
nullspace network	$2.22\text{×}{10}^{-4}$	37.0	$9.17\text{×}{10}^{-1}$
Bernoulli measurements
FBP	$20.2\text{×}{10}^{-4}$	27.3	$4.18\text{×}{10}^{-1}$
${\ell }^1$-minimization	$1.89\text{×}{10}^{-4}$	37.5	$9.06\text{×}{10}^{-1}$
residual network	$6.32\text{×}{10}^{-4}$	32.6	$8.89\text{×}{10}^{-1}$
nullspace network	$2.21\text{×}{10}^{-4}$	36.9	$8.89\text{×}{10}^{-1}$

Table 2. Table 2: Performance for the Shepp-Logan phantom.

	SME	PSNR	SSIM
Sparse measurements
${\ell }^1$-minimization	$6.73\text{×}{10}^{-4}$	31.7	$8.04\text{×}{10}^{-1}$
residual network	$6.32\text{×}{10}^{-4}$	32.0	$8.90\text{×}{10}^{-1}$
nullspace network	$5.29\text{×}{10}^{-4}$	32.8	$8.59\text{×}{10}^{-1}$
Bernoulli measurements
${\ell }^1$-minimization	$6.03\text{×}{10}^{-4}$	32.2	$8.19\text{×}{10}^{-1}$
residual network	$19.2\text{×}{10}^{-4}$	27.2	$7.63\text{×}{10}^{-1}$
nullspace network	$6.92\text{×}{10}^{-4}$	31.6	$7.67\text{×}{10}^{-1}$

Equations

$$\mathbf{g}=\mathcal{A}\mathbf{f}+ϵ,$$

$$\begin{gathered} {\partial }_t^{2}p(\mathbf{r},t)-c^2\Delta p(\mathbf{r},t) \\ ={\delta }^{\prime }(t)p_0(\mathbf{r})\text{ for }(\mathbf{r},t)\in {\mathbb{R}}^2\times {\mathbb{R}}_{+}. \end{gathered}$$

$$p_0(\mathbf{r})=-\frac{1}{\pi R}{\int }_{\partial B_R}{\int }_{|\mathbf{r}-z|}^{\infty }\frac{({\partial }_{t}tp)(\mathbf{s},t)}{\sqrt{t^2-|\mathbf{r}-\mathbf{s}{|}^2}}\mathrm{d}t\mathrm{d}S(\mathbf{s}).$$

$$p\left({\mathbf{s}}_k,t_{\ell }\right)\text{ for }(k,\ell )\in \left\{1,\dots ,M\right\}\times \left\{1,\dots ,Q\right\},$$

$${\mathbf{s}}_k$$

$$t_{\ell }$$

$$M\ge 2R_0{\lambda }_0$$

$$\mathbf{g}(j,\ell )\triangleq \sum _{k=1}^M\mathbf{S}[j,k]p({\mathbf{r}}_k,t_{\ell })\text{ for }j\in \left\{1,\dots ,m\right\},$$

$$\mathbf{g}=\mathcal{A}\mathbf{f}\text{ with }\mathcal{A}\triangleq \mathcal{S}\circ \mathcal{W}\in {\mathbb{R}}^{mQ\times n}.$$

$$(1-\delta )\Vert \mathbf{z}{\Vert }^2\le \Vert \mathcal{A}\mathbf{z}{\Vert }^2\le (1+\delta )\Vert \mathbf{z}{\Vert }^2,$$

$$\exists \mathbf{\eta }\in {\mathbb{R}}^{mQ}:{\mathcal{A}}^{\mathsf{T}}\mathbf{\eta }\in \mathrm{sign}(\mathbf{h})$$

$$\left|({\mathcal{A}}^{\mathsf{T}}\mathbf{\eta }{)}_i\right|<1\text{ for }i\notin \mathrm{supp}(\mathbf{h}),$$

$${\mathbf{h}}_{\beta }^{\delta }\in \underset{\mathbf{z}}{\mathrm{argmin}}\frac{1}{2}\Vert \mathcal{A}\mathbf{z}-{\mathbf{g}}^{\delta }{\Vert }_2^2+\beta \Vert \mathbf{z}{\Vert }_1$$

$$\begin{gathered} {\partial }_t^{2}q(\mathbf{r},t)-c^2\Delta q(\mathbf{r},t) \\ ={\delta }^{\prime }(t)c^2\Delta \mathbf{f}(\mathbf{r})\text{ for }(\mathbf{r},t)\in {\mathbb{R}}^2\times {\mathbb{R}}_{+}. \end{gathered}$$

$$\forall \mathbf{f}\in {\mathbb{R}}^n:{\mathcal{D}}_t^2\mathcal{A}[\mathbf{f}]=\mathcal{A}[c^2{\mathcal{L}}_{\mathbf{r}}\mathbf{f}],$$

$$\underset{\mathbf{h}}{\min }\Vert \mathbf{h}{\Vert }_1\text{such that }\mathcal{A}\mathbf{h}={\mathcal{D}}_t^2\mathbf{g}.$$

$$\underset{(\mathbf{f},\mathbf{h})}{\min }\Vert \mathbf{h}{\Vert }_1+I_C(\mathbf{f})$$

$$\text{such that }\left[\begin{array}{c}\mathcal{A}\mathbf{f},\mathcal{A}\mathbf{h},{\mathcal{L}}_{\mathbf{r}}\mathbf{f}-\mathbf{h}/c^2\end{array}\right]=\left[\begin{array}{c}\mathbf{g},{\mathcal{D}}_t^2\mathbf{g},0\end{array}\right].$$

$$\begin{gathered} \frac{1}{2}\Vert \mathcal{A}\mathbf{f}-\mathbf{g}{\Vert }_2^2+\frac{1}{2}\Vert \mathcal{A}\mathbf{h}-{\mathcal{D}}_t^2\mathbf{g}{\Vert }_2^2+\frac{\alpha }{2}\Vert {\mathcal{L}}_{\mathbf{r}}\mathbf{f}-\mathbf{h}/c^2{\Vert }_2^2 \\ +\beta \Vert \mathbf{h}{\Vert }_1+I_C(\mathbf{f})\to \underset{(\mathbf{f},\mathbf{h})}{\min }. \end{gathered}$$

$$\begin{gathered} \Phi (\mathbf{f},\mathbf{h})\triangleq \frac{1}{2}\Vert \mathcal{A}\mathbf{f}-\mathbf{g}{\Vert }_2^2 \\ +\frac{1}{2}\Vert \mathcal{A}\mathbf{h}-{\mathcal{D}}_t^2\mathbf{g}{\Vert }_2^2+\frac{\alpha }{2}\Vert {\mathcal{L}}_{\mathbf{r}}\mathbf{f}-\mathbf{h}/c^2{\Vert }_2^2 \end{gathered}$$

$${\mathrm{prox}}_F(\mathbf{f})\triangleq \mathrm{argmin}\left\{F(\mathbf{z})+\frac{1}{2}\Vert \mathbf{f}-\mathbf{z}{\Vert }_2^2\mid \mathbf{z}\in {\mathbb{R}}^n\right\}.$$

$${\nabla }_{\mathbf{f}}\Phi (\mathbf{f},\mathbf{h})$$

$${\nabla }_{\mathbf{h}}\Phi (\mathbf{f},\mathbf{h})$$

$$\mathrm{prox}(\mathbf{f},\mathbf{h})$$

$${\mathrm{prox}}_{I_C}(\mathbf{f}{)}_i$$

$${\mathrm{prox}}_{\beta \Vert \cdot {\Vert }_1}(\mathbf{h}{)}_i$$

$${\mathbf{f}}^{k+1}$$

$${\mathbf{h}}^{k+1}$$

$${\mathcal{R}}_{\theta }={\mathcal{N}}_{\theta }\circ {\mathcal{A}}^{♯}:{\mathbb{R}}^{mQ}\to {\mathbb{R}}^n.$$

$${\mathcal{A}}^{♯}:{\mathbb{R}}^{mQ}\to {\mathbb{R}}^n$$

$${\mathcal{N}}_{\theta }:{\mathbb{R}}^n\to {\mathbb{R}}^n.$$

$${\mathcal{R}}_{\theta }^{\mathrm{res}}=(\mathrm{Id}+{\mathcal{U}}_{\theta }){\mathcal{A}}^{♯},$$

$${\mathcal{R}}_{\theta }^{\mathrm{null}}=(\mathrm{Id}+{\mathcal{P}}_{\mathrm{Ker}(\mathcal{A})}{\mathcal{U}}_{\theta }){\mathcal{A}}^{♯}.$$

$${\mathbf{f}}^{(0)}$$

Full text

Keywords: Compressed sensing, sparsity, $\ell^{1}$-minimization, deep learning, residual learning, nullspace network

1 Introduction

Compressed sensing (CS) allows to reduce the number of measurements in photoacoustic tomography (PAT) while preserving high spatial resolution. A reduced number of measurements can increase the measurement speed and reduce system costs [1, 2, 3, 4, 5] . However, CS PAT image reconstruction requires special algorithms to achive high resolution imaging. In this work, we compare $\ell^{1}$-minimization and deep learning algorithms for 2D PAT. Among others, the two-dimensional case arises in PAT with integrating line detectors [6, 7].

In the case that a sufficiently large number of detectors is used, according to Shannon’s sampling theory, implementations of full data methods yield almost artifact free reconstructions [8]. As the fabrication of an array of detectors is demanding, experiments using integrating line detectors are often carried out using a single line detector, scanned on circular paths using scanning stages [9, 10], which is very time consuming. Recently, systems using arrays of $64$ parallel line detectors have been demonstrated [11, 12]. To keep production costs low and to allow fast imaging, the number of measurements will typically be kept much smaller than advised by Shannon’s sampling theory and one has to deal with highly under-sampled data.

After discretization, image reconstruction in CS PAT consists in solving the inverse problem

[TABLE]

where $\boldsymbol{f}\in\mathbb{R}^{n}$ is the discrete photoacoustic (PA) source to be reconstructed, $\boldsymbol{g}\in\mathbb{R}^{mQ}$ are the given CS data, $\epsilon$ is the noise in the data and $\mathcal{A}\colon\mathbb{R}^{n}\to\mathbb{R}^{mQ}$ is the forward matrix. The forward matrix is the product of the PAT full data problem and the compressed sensing measurement matrix. Making CS measurements in PAT implies that $mQ\ll n$ and therefore, even in the case of exact data, solving (1.1) requires particular reconstruction algorithms.

1.1 CS PAT recovery algorithms

Standard CS reconstruction techniques for (1.1) are based on sparse recovery via $\ell^{1}$-minimization. These algorithms rely on sparsity of the unknowns in a suitable basis or dictionary and special incoherence of the forward matrix. See [2, 3, 4, 5] for different CS approaches in PAT. To guarantee sparsity of the unknowns, in [1] a new sparsification and corresponding joint $\ell^{1}$-minimization have been derived. Recently, deep learning appeared as a new reconstruction paradigm for CS and other inverse problems. Deep learning approaches for PAT can be found in [13, 14, 15, 16, 17, 18, 19].

In this work, we compare the performance of the joint $\ell^{1}$-minimization algorithm of [1] with deep learning approaches for CS PAT image reconstruction. For the latter we use the residual network [20, 21, 13] and the nullspace network [22, 23]. The nullspace network includes a certain data consistency layer and even has been shown to be a regularization method in [23]. Our results show that the nullspace network uniformly outperforms the residual network for CS PAT in terms of the mean squared error (MSE).

1.2 Outline

In Section 2, we present the required background from CS PAT. The sparsification strategy and the joint $\ell^{1}$-minimization algorithm are summarized in Section 3. The employed deep learning image reconstruction strategies using a residual network and an (approximate) nullspace network are described in Section 4. In Section 5 we present reconstruction results for sparse measurements and Bernoulli measurements. The paper ends with a discussion in Section 6.

2 Compressed photoacoustic tomography

2.1 Photoacoustic tomography

As illustrated in Figure 1.1, PAT is based on generating an acoustic wave inside some investigated object using short optical pulses. Let us denote by $p_{0}\colon\mathbb{R}^{d}\to\mathbb{R}$ the initial pressure distribution which provides diagnostic information about the patient and which is the quantity of interest in PAT [24, 6, 25]. For keeping the presentation simple and focusing on the main ideas we only consider the case of $d=2$. Among others, the two-dimensional case arises in PAT with so called integrating line detectors [6, 7]. Further, we restrict ourselves to the case of a circular measurement geometry, where the acoustic measurements are made on a circle surrounding the investigated object.

In two spatial dimensions, the induced pressure in PAT satisfies the 2D wave equation

[TABLE]

Here $\mathbf{r}\in\mathbb{R}^{2}$ is the spatial location, $t\in\mathbb{R}$ the time variable, $\Delta_{\mathbf{r}}$ the spatial Laplacian, $c$ the speed of sound, and $p_{0}(\mathbf{r})$ the PA source that is assumed to vanish outside the disc $B_{R}\triangleq\{x\in\mathbb{R}^{2}\mid\|x\|<R\}$ and has to be recovered. The wave equation (2.1) is augmented with $p(\mathbf{r},t)=0$ on $\left\{t<0\right\}$. The acoustic pressure is then uniquely defined and referred to as the causal solution of (2.1).

PAT in a circular measurement geometry consist in recovering the function $p_{0}$ from measurements of $p(\mathbf{s},t)$ on $\partial B_{R}\times(0,\infty)$. In the case of full data, exact and stable PAT image reconstruction is possible [26, 27] and several efficient methods for recovering ${\boldsymbol{f}}$ are available. As an example, we mention the FBP formula derived in [28],

[TABLE]

Note the inversion operator in (2.2) is also the adjoint of the forward operator, see [28].

2.2 Discretization

In practical applications, the acoustic pressure can only be measured with a finite number of acoustic detectors. The standard sampling scheme for PAT in circular geometry assumes uniformly sampled values

[TABLE]

with

[TABLE]

The number $M$ of detector positions in (2.3) is directly related to the resolution of the final reconstruction. Namely,

[TABLE]

equally spaced transducers are required to stably recover any PA source $p_{0}$ that has maximal essential wavelength $\lambda_{0}$ and is supported in a disc $B_{R_{0}}\subseteq B_{R}$; see [8]. Image reconstruction in this case can be performed by discretizing the inversion formula (2.2). The sampling condition (2.6) requires a very high sampling rate, especially when the PA source contains narrow features, such as blood vessels or sharp interfaces.

Note that temporal samples can easily be collected at a high sampling rate compared to the spatial sampling, where each sample requires a separate sensor. It is therefore beneficial to keep $M$ as small as possible. Consequently, full sampling in PAT is costly and time consuming and strategies for reducing the number of detector locations are desirable.

2.3 Compressive measurements in PAT

To reduce the number of measurements we use CS measurements. Instead of collecting $M$ individually sampled signals as in (2.3), we take general linear measurements

[TABLE]

with $m\ll M$. Several choices for the measurement matrix $\mathbf{S}$ are possible and have been used for CS PAT [2, 3, 4]. In this work, we take $\mathbf{S}$ as deterministic sparse subsampling matrix or Bernoulli random matrix; see Subsection 5.1.

Let us denote by $\mathcal{W}\in\mathbb{R}^{MQ\times n}$ the discretized solution operator of the wave equation and by $\mathcal{S}\triangleq\mathbf{S}\otimes\mathbf{I}\in\mathbb{R}^{mQ\times MQ}$ the Kronecker (or tensor) product between the CS measurement matrix $\mathbf{S}$ and the identity matrix $\mathbf{I}$. Then the CS data (2.7) written as column vector $\boldsymbol{g}\in\mathbb{R}^{mQ}$ are given by

[TABLE]

In the case of CS measurements we have $mQ\ll n$ and therefore (2.8) is highly underdetermined and image reconstruction requires special reconstruction algorithms.

3 Joint $\ell^{1}$-minimization for CS PAT

Standard CS image reconstruction is based on $\ell^{1}$ minimization and sparsity of the unknowns to be recovered. In [1] we introduced a sparse recovery strategy that we will use in the present paper and recall below.

3.1 Background from $\ell^{1}$-minimization

An element ${\boldsymbol{h}}\in\mathbb{R}^{n}$ is called $s$-sparse if it contains at most $s$ nonzero elements. If we are given measurements $\mathcal{A}{\boldsymbol{h}}=\boldsymbol{g}$ where ${\boldsymbol{h}}\in\mathbb{R}^{n}$ and $\boldsymbol{g}\in\mathbb{R}^{mQ}$ with $mQ\ll n$, then stable recovery of ${\boldsymbol{h}}$ from $\boldsymbol{g}$ via $\ell^{1}$-minimization can be guaranteed if ${\boldsymbol{h}}$ is sparse and the matrix $\mathcal{A}$ satisfies the restricted isometry property of order $2s$. The latter property means that for all $2s$-sparse vectors $\boldsymbol{z}\in\mathbb{R}^{n}$ we have

[TABLE]

for an RIP constant $\delta<1/\sqrt{2}$; see [29].

Bernoulli random matrices satisfy the RIP with high probability [30] whereas the subsampling matrix clearly does not satisfy the RIP. In the case of CS PAT, the forward matrix is given by $\mathcal{A}=(\mathbf{S}\otimes\mathbf{I})\circ\mathcal{W}$. It is not known whether $\mathcal{A}$ satisfies the RIP for either the Bernoulli of the subsampling matrix. In such situations one may use the following stable reconstruction result from inverse problems theory.

Theorem 3.1** ($\ell^{1}$-minimization).**

Let $\mathcal{A}\in\mathbb{R}^{mQ\times n}$ and ${\boldsymbol{h}}\in\mathbb{R}^{n}$ Assume

[TABLE]

where $\mathrm{sign}({\boldsymbol{h}})$ is the set valued signum function and $\operatorname{supp}({\boldsymbol{h}})$ the set of all nonzero entries of ${\boldsymbol{h}}$, and that the restriction of $\mathcal{A}$ to the subspace spanned by $e_{i}$ for $i\in\operatorname{supp}({\boldsymbol{h}})$ is injective. Then for any $\boldsymbol{g}^{\delta}\in\mathbb{R}^{mQ}$ with $\|\mathcal{A}{\boldsymbol{h}}-\boldsymbol{g}^{\delta}\|_{2}\leq\delta$, any minimizer of the $\ell^{1}$-Tikhonov functional

[TABLE]

satisfies $\|{\boldsymbol{h}}_{\beta}^{\delta}-{\boldsymbol{h}}\|_{2}=\mathcal{O}(\delta)$ provided $\beta\asymp\delta$. In particular, ${\boldsymbol{h}}$ is the unique $\|\,\cdot\,\|_{1}$-minimizing solution of $\mathcal{A}\boldsymbol{z}=\boldsymbol{g}$.

Proof.

See [31]. ∎

In [32, 31] it is shown that the RIP implies the conditions in Theorem 3.1. Moreover, the smaller $\operatorname{supp}({\boldsymbol{h}})$, the easier the conditions in Theorems are satisfied. Therefore, sufficient sparsity of the unknowns is a crucial condition for the success of $\ell^{1}$-minimization.

3.2 Sparsification strategy

The used CSPAT approach in [1] is based on following theorem which allows bringing sparsity into play.

Theorem 3.2**.**

Let $p_{0}$ be a given PA source vanishing outside $B_{R}$, and let $p$ denote the causal solution of (2.1). Then $\partial_{t}^{2}p$ is the causal solution of

[TABLE]

In particular, up to discretization error, we have

[TABLE]

where $\mathcal{A}=(\mathbf{S}\otimes\mathbf{I})\circ\mathcal{W}$ denotes the discrete CS PAT forward operator defined by (2.8), ${\mathcal{L}}_{\mathbf{r}}$ is the discretized Laplacian, and ${\mathcal{D}}_{t}$ the discretized temporal derivate.

Proof.

See [1]. ∎

Typical phantoms consist of smoothly varying parts and rapid changes at interfaces. For such PA sources, the modified source $c^{2}{\mathcal{L}}_{\mathbf{r}}{\boldsymbol{f}}$ is sparse or at least compressible. The theory of CS therefore predicts that the modified source can be recovered by solving via $\ell^{1}$-minimization

[TABLE]

Having obtained an approximate minimizer ${\boldsymbol{h}}$ by either solving (3.7) or its relaxed version, one can recover the original PA source ${\boldsymbol{f}}$ by subsequently solving the Poisson equation ${\mathcal{L}}_{\mathbf{r}}{\boldsymbol{f}}={\boldsymbol{h}}/c^{2}$ with zero boundary conditions. Using the above two-stage procedure, we observed disturbing low frequency artifacts in the reconstruction. Therefore, in [1] we introduced a different joint $\ell^{1}$-minimization approach based on Theorem 3.2 that jointly recovers ${\boldsymbol{f}}$ and $c^{2}{\mathcal{L}}_{\mathbf{r}}{\boldsymbol{f}}$.

3.3 Joint $\ell^{1}$-minimization framework

The modified data ${\mathcal{D}}_{t}^{2}\boldsymbol{g}$ is well suited to recover singularities of ${\boldsymbol{f}}$, but hardly contains low-frequency components of ${\boldsymbol{f}}$. On the other hand, the low frequency information is contained in the original data, which is still available to us. This motivates the following joint $\ell^{1}$-minimization problem

[TABLE]

Here $I_{C}$ is the indicator function of $C\triangleq[0,\infty)^{n}$, defined by $I_{C}({\boldsymbol{f}})=0$ if ${\boldsymbol{f}}\in C$ and $I_{C}({\boldsymbol{f}})=\infty$ otherwise, and guarantees non-negativity.

Theorem 3.3**.**

Assume that ${\boldsymbol{f}}\in\mathbb{R}^{n}$ is non-negative, that the measurement matrix $\mathcal{A}$ and the modified PA source ${\boldsymbol{h}}=c^{2}{\mathcal{L}}_{\mathbf{r}}{\boldsymbol{f}}$ satisfy Equations (3.2), (3.3), and denote $\boldsymbol{g}=\mathcal{A}{\boldsymbol{f}}$. Then, the pair $[{\boldsymbol{f}},c^{2}{\mathcal{L}}_{\mathbf{r}}{\boldsymbol{f}}]$ can be recovered as the unique solution of the joint $\ell^{1}$-minimization problem (3.8).

Proof.

See [1]. ∎

In the case the data is only approximately sparse or noisy, we propose, instead of (3.8), to solve the $\ell^{2}$-relaxed version

[TABLE]

Here $\alpha>0$ is a tuning and $\beta>0$ a regularization parameter.

3.4 Numerical minimization

We will solve (3.9) using a proximal forward-backward splitting method [33], which is well suited for minimizing the sum of a smooth and a non-smooth but convex part. In the case of (3.9) we take the smooth part as

[TABLE]

and the non-smooth part as $\Psi(f,h)\triangleq\beta\|{\boldsymbol{h}}\|_{1}+I_{C}(f)$.

The proximal gradient algorithm then alternately performs an explicit gradient step for $\Phi$ and an implicit proximal step for $\Psi$. For the proximal step, the proximity operator of a function must be computed. The proximity operator of a given convex function $F\colon\mathbb{R}^{n}\to\mathbb{R}$ is defined by [33]

[TABLE]

The regularizers we are considering here have the advantage, that their proximity operators can be computed explicitly and do not cause a significant computational overhead. The gradient $[\nabla_{\boldsymbol{f}}\Phi,\nabla_{\boldsymbol{h}}\Phi]$ of the smooth part can easily be computed to be

[TABLE]

The proximal operator of the non-smooth part is given by

[TABLE]

With this, the proximal gradient algorithm is given by

[TABLE]

where $({\boldsymbol{f}}^{k},{\boldsymbol{h}}^{k})$ is the $k$-th iterate and $\mu$ the step size. We initialize the proximal gradient algorithm with ${\boldsymbol{f}}^{0}={\boldsymbol{h}}^{0}=0$.

4 Deep learning for CS PAT

As an alternative to the joint $\ell^{1}$-minimization algorithm we use deep learning or CS image reconstruction. We thereby use a trained residual network as well as a corresponding (approximate) nullspace network, which offers improved data consistence.

4.1 Image reconstruction by deep learning

Deep learning is a recent paradigm to solve inverse problems of the form (1.1). In this case, image reconstruction is performed by an explicit reconstruction function

[TABLE]

The reconstruction operator $\mathcal{R}_{\theta}$ is the composition of a backprojection operator and a convolutional neural network

[TABLE]

The backprojection $\mathcal{A}^{\sharp}$ performs an initial reconstruction that is subsequently improved by the CNN $\mathcal{N}_{\theta}$. In this work, we use the filtered backprojection (FBP) algorithm [28] for $\mathcal{A}^{\sharp}$, which is a discretization of the inversion formula (2.2). For the CNN $\mathcal{N}_{\theta}$ we use the residual network (see Subsection 4.2) and the nullspace network (see Subsection 4.3).

The CNN is taken from a parameterized family, where parameterization $\theta\in\Theta\mapsto\mathcal{N}_{\theta}$ is determined by the network architecture. For adjusting the parameters, one assumes a family of training data $((\boldsymbol{b}_{k},\boldsymbol{f}_{k}))_{k=1}^{N}$ is given where any training example consist of artifact-free output image $\boldsymbol{f}_{k}$ and a corresponding input image $\boldsymbol{b}_{k}=\mathcal{A}^{\sharp}\mathcal{A}(\boldsymbol{f}_{k})$. The free parameters $\theta$ are chosen in such a way, that the overall error of the network for predicting $\boldsymbol{f}_{k}$ from $\boldsymbol{b}_{k}$ is minimized. The minimization procedure used in this paper is described in Subsection (5.2).

4.2 Residual network

The architecture of the CNN is a crucial step for the performance of tomographic image reconstruction with deep learning. A common architecture in that context is the following residual network

[TABLE]

where $\mathcal{U}_{\theta}$ is the Unet, originally introduced in [34] for biomedical image segmentation. The residual network 4.4 has successfully been used for various tomographic image reconstruction tasks [13, 20, 21] including PAT.

Using $\operatorname{Id}+\mathcal{U}_{\theta}$ instead of $\mathcal{U}_{\theta}$ affects that actually the residual images $\boldsymbol{f}+\boldsymbol{b}$ are learned by the Unet. The residual images often have a simpler structure than the original outputs $\boldsymbol{f}$. As argued in [21], learning the residuals and adding them to the inputs after the last layer is more effective than directly training for the outputs. The resulting deep neural network architecture is shown in Figure 4.1.

4.3 Nullspace network

Especially when applying $\mathcal{R}_{\theta}^{\rm res}$ to objects very different from the training set, the residual network (4.4) lacks data consistency, in the sense that $\mathcal{R}_{\theta}^{\rm res}\boldsymbol{g}$ is not necessarily a solution of the given equation $\mathcal{A}{\boldsymbol{f}}=\boldsymbol{g}$. To overcome this limitation, as an alternative we use the nullspace network [23],

[TABLE]

One strength of the nullspace network is that the term $\mathcal{P}_{\operatorname{Ker}(\mathcal{A})}\mathcal{U}_{\theta}$ only adds information that is consistent with the given data. For example, if $\mathcal{A}^{\sharp}=\mathcal{A}^{+}$ equals the pseudoinverse, then $\mathcal{R}_{\theta}^{\rm null}\boldsymbol{g}$ even is fully data consistent as implied by the following theorem.

Theorem 4.1**.**

Let $\boldsymbol{g}=\mathcal{A}({\boldsymbol{f}}^{\star})$ be in the range of the forward operator, write $L(\mathcal{A},\boldsymbol{g})$ for the set of solutions of the equation $\mathcal{A}{\boldsymbol{f}}=\boldsymbol{g}$ and take $\mathcal{A}^{\sharp}=\mathcal{A}^{+}$ as the pseudoinverse.

$\mathcal{R}_{\theta}^{\rm null}(\boldsymbol{g})$* is a solution of $\mathcal{A}{\boldsymbol{f}}=\boldsymbol{g}$.* 2. 2.

We have $\mathcal{R}_{\theta}^{\rm null}(\boldsymbol{g})=\mathcal{P}_{L(\mathcal{A},\boldsymbol{g})}\mathcal{R}_{\theta}^{\rm res}(\boldsymbol{g})$. 3. 3.

Consider the iteration

[TABLE]

with step size $0<s<\|\mathcal{A}\|^{-2}$. Then:

1. (a)

$\|{\boldsymbol{f}}^{\star}-{\boldsymbol{f}}^{(k)}\|$* is monotonically decreasing* 2. (b)

$\lim_{k\to\infty}{\boldsymbol{f}}^{(k)}=\mathcal{R}_{\theta}^{\rm null}(\boldsymbol{g})$** 3. (c)

$\|{\boldsymbol{f}}^{\star}-{\boldsymbol{f}}^{(k)}\|\leq\|{\boldsymbol{f}}^{\star}-\mathcal{R}_{\theta}^{\rm res}(\boldsymbol{g})\|$.

Proof.

Will be presented elsewhere. ∎

Theorem 4.1 implies that iteration (4.6), (4.7) defines a sequence

[TABLE]

that monotonically converges to $\mathcal{R}_{\theta}^{\rm null}(\boldsymbol{g})$. It implies that the nullspace network as well as the approximate nullspace network $\mathcal{R}_{\theta}^{\mathrm{null},(k)}(\boldsymbol{g})$ have a smaller reconstruction error than the residual network. Moreover, according to Theorem 4.1 the nullspace network yields a solution of the equation $\mathcal{A}{\boldsymbol{f}}=\boldsymbol{g}$ even for elements very different from the training data.

5 Numerical results

In this section we numerically compare the joint $\ell^{1}$-minimization approach with the residual network and the nullspace network. We also compare the results with plain FBP. We use Keras [35] with TensorFlow [36] to train and evaluate the CNN. The FBP, the $\ell^{1}$-minimization algorithm and the iterative update (4.7) is implemented in MATLAB. We ran all our experiments on a computer using an Intel i7-6850K and an NVIDIA 1080Ti. The phantoms as well as the FBP reconstruction from fully sampled data are shown in Figure 4.2. Note that we use limited view data which implies that the reconstructions in Figure 4.2 contain some artefacts.

5.1 Measurement setup

The entries of any discrete PA source ${\boldsymbol{f}}\in\mathbb{R}^{n}$ with $n=256^{2}$ correspond to discrete samples of the continuous source at a $256\times 256$ Cartesian grid covering the square $[$-5\text{,}\mathrm{\SIUnitSymbolMicro m}$,$9\text{,}\mathrm{\SIUnitSymbolMicro m}$]\times[$-12.5\text{,}\mathrm{\SIUnitSymbolMicro m}$,$1.5\text{,}\mathrm{\SIUnitSymbolMicro m}$]$. The full wave data $\boldsymbol{g}\in\mathbb{R}^{MP}$ corresponds to $P=747$ equidistant temporal samples in $[0,T]$ with $T=$4.9749\text{\times}{10}^{-2}\text{,}\mathrm{\SIUnitSymbolMicro s} and $M=240$ equidistant sensor locations on the circle of radius $40\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and polar angles in the interval $[$35\text{\,}\mathrm{\SIUnitSymbolDegree}$,$324\text{\,}\mathrm{\SIUnitSymbolDegree}$]$. The sound speed is taken as $c=$1.4907\text{\times}{10}^{3}\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}. The wave equation is evaluated by discretizing the solution formula of the wave equation, and the inversion formula (2.2) is discretized using the standard FBP procedure described in [28, 37]. Recall that the continuous setting the inversion integral in (2.2) equals the adjoint of the forward operator. Therefore the above procedure gives a pair

[TABLE]

of forward operator and unmatched adjoint.

We consider $m=60$ spatial measurements which corresponds to a compression factor of four. For the sampling matrices $\mathbf{S}\in\mathbb{R}^{m\times M}$ we use the following instances:

• •

Deterministic sparse subsampling matrix with entries

[TABLE]

• •

Random Bernoulli matrix where each entry is taken independently as $\pm 1/\sqrt{m}$ with equal probability.

The Bernoulli matrix satisfies then RIP with high probability, whereas the sparse subsampling matrix doesn’t. Therefore, we expect the $\ell^{1}$-minimization approach to work better for Bernoulli measurements. On the other hand, in the subsampling case the artefacts have more structure which therefore is expected to be better for the deep learning approaches. Our findings below confirm these conjectures.

5.2 Construction of reconstruction networks

For the residual and the nullspace network we use the backprojection layer

[TABLE]

and the same trained CNN. For that purpose, we construct $N=5000$ training examples $(\boldsymbol{b}_{k},\boldsymbol{f}_{k})_{k=1}^{N}$ where $\boldsymbol{f}_{k}$ are taken as projection images from three dimensional lung blood vessel data as described in [17]. All images $\boldsymbol{f}_{k}$ are normalized to have maximal intensity one. The corresponding input images are computed by $\boldsymbol{b}_{k}=\mathcal{A}^{\sharp}\mathcal{A}\boldsymbol{f}_{k}$. The CNN is constructed by minimizing the mean absolute error

[TABLE]

using stochastic gradient descent with batch size 1 and momentum 0.9. We trained for 200 epochs and used a decaying learning parameter between $0.005$ to $0.0025$.

Having computed the minimizer of (5.3) we use the trained residual network $\mathcal{R}_{\theta}^{\rm res}$ as well as the corresponding approximate nullspace network $\mathcal{R}_{\theta}^{\mathrm{null},(10)}$ for image reconstruction.

5.3 Blood vessel phantoms

First we investigate the performance on 50 blood vessel phantoms that are not contained in the training set. We consider sparse sampling as well as Bernoulli measurements. For the joint recovery approach, we use 70 iterations of the iterative thresholding procedure with coupling parameter $\alpha=0.001$, regularization parameter $\beta=0.005$ and step size $\mu=0.125$. For the (approximate) nullspace network $\mathcal{R}_{\theta}^{\mathrm{null},(10)}$ we use 10 iterations to approximately compute the projection. Results for one of the vessel phantoms are visualized in Figure 5.1. To quantitatively evaluate the results we computed the MSE (mean square error), the PSNR (peak signal to noise ratio) and the SSIM (structural similarity index [38]) averaged over all 50 blood vessel phantoms. The reconstruction errors are summarized in Table 1 where the best results are framed.

From Table 1 we see that the hybrid as well as the deep learning based methods significantly outperform the FBP reconstruction. Moreover, the deep learning approach even outperforms the joint recovery approach for the sparse sampling. The nullspace network in all cases decreases the MSE (increases the PSNR) compared to the residual network.

5.4 Shepp-Logan type phantom

Next we investigate the performance on a Shepp-Logan type phantom that contains structures completely different from the training data. For the joint recovery approach, we use 50 iterations of the iterative thresholding procedure with $\alpha=0.001$ regularization parameter $\beta=0.005$ and step size $\mu=0.1$. For the nullspace network we use $\mathcal{R}_{\theta}^{\mathrm{null},(10)}$. Results are shown in Figure 5.2. Table 2 shows the MSE, the PSNR and the SSIM for the head phantom, where the best results are again framed.

As the considered Shepp-Logan type phantom is very different from the training data it is not surprisingly the standard residual network does not perform that well for the Bernoulli measurements. Surprisingly the residual network still works well in the sparse data case. The nullspace network yields significantly improved results compared for the residual network, especially for the Bernoulli case. In the Bernoulli case, the $\ell^{1}$-minimization approach performs best, however only slightly better than the approximate nullspace network.

6 Conclusion

In this paper we compared $\ell^{1}$-minimization with deep learning for CS PAT image reconstruction. The two approaches have been tested on blood vessel data (test data not contained in the training set that consists of similar objects) as well as a Shepp-Logan type phantom (with structures very different from the training data). For the CS PAT measurements we considered deterministic subsampling as well as random Bernoulli measurements. For the used reconstruction networks, we considered the Unet with residual connection and an approximate nullspace network which contains an additional data consistency layer.

In terms of reconstruction quality, our findings can be summarized as follows:

Sparse recovery and deep learning both significantly outperforms filtered backprojection for both measurement matrices. If the training data are not accurate for the object to be reconstructed, for the deep learning approach this conclusion only holds for the null-space network. 2. 2.

In the case of the sparse measurement matrix, the deep learning approach outperforms $\ell^{1}$-minimization. In the case of Bernoulli measurement, the $\ell^{1}$-minimization algorithms yields better performance. 3. 3.

The nullspace network contains a data consistence layer and yields good results even for phantoms very different from the training data. Even for the test data similar to the training data it yields an improved PSNR compared to the residual network

According to the above results we can recommend the $\ell^{1}$-minimization algorithm in the case of random measurements and the nullspace network in the case of sparse measurements. We point out that application of the CNN only takes fractions of second (actually, less than $0.01$ seconds) in Keras whereas the joint recovery approach requires around 2 minutes for 50 iterations in Matlab. Note that this comparison is not completely fair and with a recent GPU implementation in PyTorch we have been able the reduce the computation time to about one second for 50 iterations. Nevertheless, the deep learning based methods are still significantly faster. Therefore, especially the nullspace network is very promising for high quality real-time CS PAT imaging.

Acknowledgement

The work of M.H and S.A. has been supported by the Austrian Science Fund (FWF), project P 30747-N32.