Photo-acoustic tomographic image reconstruction from reduced data using physically inspired regularization

TL;DR

This paper introduces a physically inspired non-convex regularization method for photoacoustic tomography that enables high-quality image reconstruction from reduced data sets, utilizing a custom minimization algorithm for efficiency.

Contribution

It develops a novel non-convex regularization tailored to PAT images and a specialized minimization method, improving reconstruction quality from limited data.

Findings

Effective reconstruction from reduced datasets.

Outperforms existing FISTA-based methods.

Maintains image quality with less data.

Abstract

We propose a model-based image reconstruction method for photoacoustic tomography(PAT) involving a novel form of regularization and demonstrate its ability to recover good quality images from significantly reduced size datasets. The regularization is constructed to suit the physical structure of typical PAT images. We construct it by combining second-order derivatives and intensity into a non-convex form to exploit a structural property of PAT images that we observe: in PAT images, high intensities and high second-order derivatives are jointly sparse. The specific form of regularization constructed here is a modification of the form proposed for fluorescence image restoration. This regularization is combined with a data fidelity cost, and the required image is obtained as the minimizer of this cost. As this regularization is non-convex, the efficiency of the minimization method is crucial in obtaining artifact-free reconstructions. We develop a custom minimization method for efficiently handling this non-convex minimization problem. Further, as non-convex minimization requires a large number of iterations and the PAT forward model in the data-fidelity term has to be applied in the iterations, we propose a computational structure for efficient implementation of the forward model with reduced memory requirements. We evaluate the proposed method on both simulated and real measured data sets and compare them with a recent reconstruction method that is based on a well-known fast iterative shrinkage threshold algorithm (FISTA).

Tables

Table 1. Table 1: SSIM scores of reconstruction with varying input noise levels and different number of transducers for all phantoms. FBM: FISTA-based method of Huang et al. [ 25 ]

Blood Vessel Phantom					Derenzo Phantom				PAT Phantom
SNR	64 trans.		128 trans.		16 trans.		32 trans.		16 trans.		32 trans.
(dB)	FBM	Ours	FBM	Ours	FBM	Ours	FBM	Ours	FBM	Ours	FBM	Ours
20	.640	.758	.705	.762	. 682	.983	.914	.994	.691	.907	.805	.957
30	.723	.762	.737	.763	.725	.997	.982	.999	.774	.976	.976	.999
40	.743	.762	.742	.764	.740	.999	.989	.999	.780	.997	.985	.999

Equations

$$\Big{(}\frac{\partial^{2}}{\partial t^{2}}-c^{2}_{0}\;\nabla^{2}\Big{)}\;p(\textbf{r},t)=\Gamma\;\frac{\partial}{\partial t}H(\textbf{r},t),$$

$$p(\text{r},t)=\frac{\Gamma }{4\pi c_0}\frac{\partial }{\partial t}\underset{|\text{r}-{\text{r}}^{\prime }|=c_{0}t}{\int }\frac{H({\text{r}}^{\prime })}{|\text{r}-{\text{r}}^{\prime }|}d{\text{r}}^{\prime },$$

$$p_0(\text{r})=\Gamma H(\text{r}).$$

$$p(\textbf{r},t)=\mathcal{F}^{-1}\Big{\{}\hat{P}_{o}(\textbf{k})\cos(c_{0}\|\textbf{k}\|t)\Big{\}},$$

$${\text{p}}_m=\text{H}{\text{p}}_0,$$

$${\text{p}}_m=\text{H}{\text{p}}_0+\mathbf{\eta },$$

$${\hat{\text{p}}}_0=\underset{{\text{p}}_0}{\mathrm{argmin}}J({\text{p}}_0)$$

$$J({\text{p}}_0)=||{\text{p}}_m-\text{H}{\text{p}}_0|{|}_2^2+\lambda R({\text{p}}_0).$$

$${R}(\textbf{p}_{0})=\sum_{i}\|\textbf{D}_{o,i}\textbf{p}_{0}\|^{2}_{2}=\sum_{r=1}^{N}\sum_{i}\big{(}(\textbf{D}_{o,i}\textbf{p}_{0})_{r}\big{)}^{2}$$

$$\textbf{p}_{0}=\big{[}\textbf{H}^{T}\textbf{H}+\lambda\;\sum_{i}\textbf{D}_{o,i}^{T}\textbf{D}_{o,i}\big{]}^{-1}\textbf{H}^{T}\textbf{p}_{m}.$$

$${R}_{TV}(\textbf{p}_{0})=\sum_{r=1}^{N}\sqrt{\sum_{i}\big{(}\textbf{D}_{o,i}\textbf{p}_{0}\big{)}_{r}^{2}}$$

$${R}_{TV}(\textbf{p}_{0})=\sum_{r=1}^{N}\sqrt{\epsilon+\sum_{i}\big{(}\textbf{D}_{o,i}\textbf{p}_{0}\big{)}_{r}^{2}}$$

$${R}_{h,1}(\textbf{p}_{0},q)=\sum_{r=1}^{N}\left(\epsilon+\alpha\big{(}\textbf{p}_{0}\big{)}_{r}^{2}+(1-\alpha)\sum_{i}\big{(}\textbf{D}_{o,i}\textbf{p}_{0}\big{)}_{r}^{2}\right)^{q},\;\;s.t.\;\;0<\alpha<1,$$

$${R}_{h,2}(\textbf{p}_{0},q)=\alpha\sum_{r=1}^{N}\left(\epsilon+\big{(}\textbf{p}_{0}\big{)}_{r}^{2}\right)^{q}+(1-\alpha)\sum_{r=1}^{N}\left(\epsilon+\sum_{i}\big{(}\textbf{D}_{o,i}\textbf{p}_{0}\big{)}_{r}^{2}\right)^{q}.$$

$${\hat{\text{p}}}_0=\underset{{\text{p}}_0\ge 0}{\mathrm{argmin}}J({\text{p}}_0,q)$$

$$J({\text{p}}_0,q)=||{\text{p}}_m-\text{H}{\text{p}}_0|{|}_2^2+\lambda R_{h,l}({\mathbf{p}}_0,q)l=1\text{mbox}or2.$$

$$J(\textbf{p}_{0},q)=||\textbf{p}_{m}-\textbf{H}\textbf{p}_{0}||_{2}^{2}+\lambda R_{h,i}(\mathbf{p}_{0},q)+\lambda_{p}\|\mathcal{P}^{-}(\textbf{p}_{0}\big{)}\|_{2}^{2}$$

$$({\mathcal{P}}^{-}({\text{p}}_0){)}_r=\{\begin{array}{ll}0& \text{if}({\text{p}}_0{)}_r\ge 0\\ ({\text{p}}_0{)}_r& \text{if}({\text{p}}_0{)}_r<0.\end{array}$$

$$\nabla J(\mathbf{x},q)={\mathbf{A}}^{(\mathbf{x})}\mathbf{x}-{\mathbf{H}}^T{\mathbf{p}}_m,$$

$${\mathbf{A}}^{(\mathbf{x})}={\text{H}}^T\text{H}+\lambda \alpha {\mathbf{W}}^{(\mathbf{x})}+\lambda (1-\alpha )\sum _i{\mathbf{D}}_{o,i}^T{\mathbf{W}}^{(\mathbf{x})}{\mathbf{D}}_{o,i}+{\lambda }_p{\mathbf{N}}^{(\mathbf{x})}.$$

$$\{{\mathbf{W}}^{(\mathbf{x})}{\}}_{jj}$$

$$\{{\mathbf{N}}^{(\mathbf{x})}{\}}_{jj}$$

$${\mathbf{A}}^{({\mathbf{x}}^{(k)})}={\mathbf{H}}^T\mathbf{H}+\lambda \alpha {\mathbf{W}}^{({\mathbf{x}}^{(k)})}+\lambda (1-\alpha )\sum _i{\mathbf{D}}_{o,i}^T{\mathbf{W}}^{({\mathbf{x}}^{(k)})}{\mathbf{D}}_{0,i}+{\lambda }_p{\mathbf{N}}^{({\mathbf{x}}^{(k)})}.$$

$$\big{[}\textbf{H}^{T}\textbf{H}+\lambda\alpha{\bf I}+\lambda(1-\alpha)\sum_{i}\textbf{D}_{o,i}^{T}\textbf{D}_{o,i}\big{]}\mathbf{y}^{(-1)}=\textbf{H}^{T}\textbf{p}_{m}.$$

$${\mathbf{A}}^{(\mathbf{x})}={\text{H}}^T\text{H}+\lambda \alpha {\bar{\mathbf{W}}}^{(\mathbf{x})}+\lambda (1-\alpha )\sum _i{\mathbf{D}}_{o,i}^T{\hat{\mathbf{W}}}^{(\mathbf{x})}{\mathbf{D}}_{o,i}+{\lambda }_p{\mathbf{N}}^{(\mathbf{x})}.$$

$$\{{\bar{\mathbf{W}}}^{(\mathbf{x})}{\}}_{jj}$$

$$\{{\hat{\mathbf{W}}}^{(\mathbf{x})}{\}}_{jj}$$

$${\mathbf{H}}_i\mathbf{v}=\mathcal{P}{\mathcal{F}}^{-1}\left\{\mathcal{F}(\mathcal{A}\mathbf{v})\cos (c_0|\text{k}|t_i)\right\},$$

$${\mathbf{H}}_i^T\mathbf{u}={\mathcal{A}}^a{\mathcal{F}}^{-1}\left\{\cos (c_0|\text{k}|t_i)\mathcal{F}\left\{{\mathcal{P}}^a\mathbf{u}\right\}\right\},$$

$${\mathbf{H}}_i^T{\mathbf{H}}_i\mathbf{v}={\mathcal{A}}^a{\mathcal{F}}^{-1}\left\{\cos (c_0|\text{k}|t_i)\mathcal{F}\left\{{\mathcal{P}}^a\mathcal{P}{\mathcal{F}}^{-1}\left\{\mathcal{F}(\mathcal{A}\mathbf{v})\cos (c_0|\text{k}|t_i)\right\}\right\}\right\}$$

$${\mathbf{H}}_i^T{\mathbf{H}}_i\mathbf{v}={\mathcal{A}}^a{\mathcal{F}}^{-1}\left\{\cos (c_0|\text{k}|t_i)\mathcal{F}\left\{S(\mathbf{r}){\mathcal{F}}^{-1}\left\{\mathcal{F}(\mathcal{A}\mathbf{v})\cos (c_0|\text{k}|t_i)\right\}\right\}\right\}$$

Full text

K****eywords Data reduction methods ; Image reconstruction in medical imaging

1 Introduction

Photoacoustic tomography (PAT) [1, 2, 3, 4, 5, 6, 7] provides high resolution and high contrast images of deep tissues by imaging the spatial distribution of certain substances that can absorb near infra-red optical energy. Upon shining with a laser pulse, the substance under investigation absorbs the optical energy and undergoes thermoelastic expansion; thus, the spatial distribution of the concentration of the substance gets translated into the distribution of pressure-rise. This initial pressure rise travels outwards as ultrasound waves which are collected by ultrasound transducers placed at the boundary. From the ultrasound signal measured by the transducers as a function of time, a PAT reconstruction method recovers an estimate of the initial pressure-rise by solving the associated inverse problem [8]. Achieving accuracy in image reconstruction is a fundamental challenge in PAT; the reconstructed images suffer from artifacts that depend on the measurement geometry mainly because the reconstruction problem is complex.

There is a class of methods known as analytical inversion methods, which obtain reconstruction (recovery of initial pressure rise) by some transformation on the measured data. Among such methods, filtered back projection (FBP) [9, 10, 11, 12] and delay and sum method [13, 14] is fast and memory-efficient but require a large number of measurements for good quality reconstruction. This leads to increased scan time or expensive instrumentation setups. Time reversal methods [15] are the least demanding class of methods in the category of analytical methods and can be used for any detection geometry [16], acoustic absorption, and dispersion [17], and can accommodate the heterogeneities present in the medium [18]. However, in many applications of PAT, due to the geometrical limitations [19] or by choice to accelerate data acquisition [20], we may have to put restrictions on the spatial and/or temporal sampling of the photoacoustic (PA) signal. In such situations, these methods suffer from blurring and curved-line artifacts.

Model-based reconstruction methods outperform the direct analytic reconstruction methods in the case of limitation in the size of measured data and yield better quality in reconstruction [21, 22, 23, 24]. These methods express the measured pressure samples as a linear transformation on the initial pressure-rise. This transformation is implemented in the form of multiplication of the vector containing image pixels with a very large matrix that represents the PAT forward model (model matrix). This transformation also goes into the data fidelity term, which is minimized along with a regularization functional to achieve reconstruction. The regularization represents a prior belief on the spatial characteristics of the image to be recovered [25, 26, 27, 28]. The minimization is achieved using iterations involving the repeated application of the above-mentioned transformation and its adjoint leading to high computational burden and memory-burden. To reduce memory requirements and computational complexity, several methods have been reported. For example, the memory overhead can be reduced by decreasing the number of measurements [29] or by calculating the matrix-vector products on-the-fly without explicitly storing the model matrix [30]. However, since the same operations need to be performed multiple times, the computational time becomes longer. A fast implementation can be achieved by simplifying the forward model, assuming that the photoacoustic sources lie in a plane [31]. It can also be made faster by decoupling the inversion into smaller subproblems using a discrete wavelet packet decomposition [32]. However, the method uses interpolation, which leads to modeling errors. A more accurate model matrix calculation based on the direct discretization of Poisson-type integral showed better reconstruction accuracy in high-noise and low-SNR imaging conditions [33].

We develop a novel model-based reconstruction method that yields significantly improved reconstruction from datasets of reduced size. Our contributions are threefold. First, we note that the existing model-based methods in PAT use generic regularization forms developed for general images. There has been no reported application of a regularization that can cater to the structural properties of the photoacoustic images. We observe that PAT images have the following characteristics, which are also observed in fluorescence images: in these images, high intensities and high second-order derivatives are jointly sparse. This property was exploited for fluorescence image reconstruction by constructing a regularization that combines the intensity and the second-order derivatives [34]. Here we modify this form such that it is more suitable for the current reconstruction problem. Second, we construct a novel preconditioned gradient method for efficiently minimizing the cost function that combines this regularization and the data-fidelity term. Finally, we develop an efficient method for on-the-fly computation of matrix-vector products derived using a model that uses an exact time propagator to calculate the acoustic field [35, 36]. By taking advantage of the filtering structure present in this formulation, forward and adjoint matrix operations are implemented using filters in a memory-efficient way. We evaluate the proposed method and compare it with a recent reconstruction method based on a well-known fast iterative shrinkage threshold algorithm (FISTA). For this purpose, we use simulated data corresponding to the reduced number of transducers and real measured data.

1.1 Forward model in photoacoustic tomography

The forward problem in PAT accounts for the calculation of the pressure fields in space and time $p(\textbf{r},t)$ from a known photoacoustic source $H(\textbf{r},t)$, which represents the light energy deposited in the medium per unit volume per unit time. The induced pressure waves $p(\textbf{r},t)$ under the condition of thermal and stress confinements obey the following differential equation for an acoustically homogeneous medium [8, 4, 37]

[TABLE]

where $\Gamma$ is a dimensionless parameter called the Gr$\ddot{u}$neisen coefficient, which describes the conversion efficiency of heat to pressure and $c_{0}$ is the speed of sound in the medium. By recognizing that the temporal duration of the laser pulse is shorter than the temporal resolution of the ultrasound detectors in most of the photoacoustic imaging applications, the PA source $H(\textbf{r},t)$ may be approximated by $H(\textbf{r})\delta(t)$, where $H(\textbf{r})$ is the density of deposited energy. Then the solution to the differential equation (1) can be written as [8]

[TABLE]

where the initial pressure field ${p}_{0}(\textbf{r})=p(\textbf{r},t=0)$ can be written as

[TABLE]

The pressure distribution, $p(\textbf{r},t)$, can also be expressed as [36, 35]

[TABLE]

where $\mathcal{F}^{-1}$ denotes the Fourier inversion, and $\hat{P}_{o}(\textbf{k})$ is the Fourier transform of ${p}_{0}(\textbf{r})$ with $\mathbf{k}$ denoting the 2D Fourier frequency.

1.2 Discrete forward model

The discrete representation of the forward problem is the basis of model-based reconstruction algorithms in PAT. Since the measured pressure field is linearly related to the photoacoustic source, the discretization of the forward problem may be written in a matrix form as [22]

[TABLE]

where $\textbf{p}_{m}$ is a $LM\times 1$ vector representing the discrete pressure measurements from $L$ transducers each taking $M$ time-samples and H is the model matrix of size $LM\times N$ with $N$ being the total number of pixels. The initial pressure distribution $p_{0}(\textbf{r})$ is represented spatially using a 2D imaging grid having $N_{x}$ and $N_{y}$ grid points along $x$ and $y$ directions respectively and is denoted by a $N\times 1$ vector $\textbf{p}_{0}$ where $N=N_{x}\times N_{y}$. The model matrix can be calculated by discretizing the integral relation in (2) in the ideal case of homogeneous lossless medium and point detectors. Several methods to improve the accuracy of the model by incorporating transducer responses, heterogeneity in the medium and interpolation techniques for accurate discretization has been reported [25, 38, 31]. Including the measurement noise, the modified imaging model can be written as

[TABLE]

where ${\boldsymbol{\eta}}$ represents measurement noise, which is Gaussian.

1.3 Model based image reconstruction

A reconstruction method in PAT aims to recover the initial pressure distribution $\textbf{p}_{0}$ from the noisy transducer measurement data $\textbf{p}_{m}$. In the limited data case, the PA image reconstruction problem is ill-posed, and hence constraints are imposed on the required solution in the form of regularization. In this case, the image reconstruction problem can be treated as an optimization problem where the solution is obtained by minimizing a cost function. The reconstruction problem can be written as

[TABLE]

where $J(\textbf{p}_{0})$ is the cost function and is given by

[TABLE]

Here ${R}(\textbf{p}_{0})$ is the regularization functional and $\|\cdot\|_{2}$ represents the $l_{2}$ norm. The regularization parameter $\lambda$ controls the amount of regularity in the solution and fidelity to the measured data. The regularization functional typically should be able to smooth the noise and hence contains derivative terms in their formulation. A quadratic and differentiable regularization functional called Tikhonov regularization has been used in limited data cases [8, 39] and is given by

[TABLE]

where $(\cdot)_{r}$ denotes the $r$th component of its vector argument, and $\textbf{D}_{o,i}$ represents the matrix of $i^{th}$ derivative filter of order $o$. For example, $\textbf{D}_{1,i},i=1,2$, are the matrix equivalents of filtering by discrete filters that implement the operators $\frac{\partial}{\partial x}$ and $\frac{\partial}{\partial y}$. Further, $\textbf{D}_{2,i},i=1,2,3$, are the matrix equivalents of filtering by discrete filters that implement the operators $\frac{\partial^{2}}{\partial x^{2}},\frac{\partial^{2}}{\partial y^{2}}$ and $\sqrt{2}\frac{\partial}{\partial x\partial y}$. The resulting minimization of the convex quadratic cost function yields a closed form solution given by

[TABLE]

As it precludes any large derivative values, Tikhonov regularization tends to smooth edges in the reconstructed image.

It should be emphasized that, as of now, only the data-fidelity term requires building a large matrix ${\bf H}$, and the regularization does not require building the matrices. The matrices, $\textbf{D}_{o,i}$s, in the regularization correspond to discrete filters implementing the derivatives $\frac{\partial}{\partial x}$, $\frac{\partial}{\partial y}$, $\frac{\partial^{2}}{\partial x^{2}}$, $\frac{\partial^{2}}{\partial y^{2}}$ and $\sqrt{2}\frac{\partial}{\partial x\partial y}$. These filters can be directly applied to the images without building the matrices, which is more efficient. We use the matrix formulation for derivatives also for achieving notational convenience in describing the minimization method. In the latter part of our development, we will eliminate the need for building the matrix ${\bf H}$ as well.

A widely used non-quadratic regularization is the Total Variation (TV) [40, 20, 41, 25] and it is superior to the quadratic functional in its ability to preserve edges and it is robust to noise. The discrete total variation is given by

[TABLE]

Often a differentiable approximation of the total variation is used which can be written as

[TABLE]

where $\epsilon$ is a small positive number.

2 The proposed reconstruction method

2.1 Proposed Regularization functional

The photoacoustic images have high contrast due to the differential absorption of light in the near-infrared region by chromophores such as hemoglobin. Due to this, high values of initial pressure, $p_{0}(\mathbf{r})$, are sparsely distributed. Further, regions having high derivative values are also sparsely distributed. This pattern is also observed for fluorescence images in the work presented in [34], where the regularization was constructed by adding an intensity term to second-order derivatives. The combined point-wise cost went into a logarithmic function and summed over all pixels. Here we replace the log by a fractional power and write the proposed regularization as

[TABLE]

where the weight $\alpha\in(0,1)$ controls the relative penalization. The advantage of this modification is that it allows an optimization strategy that can efficiently handle non-convex cost function, which will be demonstrated later. Here, we choose $q<0.5$ meaning that the resulting cost functional is non-convex. We also consider a variant of the above form, which is given below:

[TABLE]

For notational convenience in describing the minimization algorithm, we call $q$ the sparsity index.

2.2 The complete cost functional

The initial pressure distribution in PAT is proportional to the fluence distribution and absorption coefficient in the tissue and hence the recovered PAT images should contain only non-negative values in it. Therefore, a non-negativity constraint is imposed on the solution of the optimization problem. The modified optimization problem to be solved is given by

[TABLE]

where

[TABLE]

This constrained optimization problem can be efficiently handled by splitting-based optimization methods. However, such methods have slow convergence and yield reconstructions with artifacts when the regularization is non-convex. Here, we intend to develop a method based on traditional gradient-based iterative schemes, and for this purpose, we introduce an approximate unconstrained formulation. Specifically, we modify cost by adding quadratic penalty term to enforce positivity as given below:

[TABLE]

where

[TABLE]

Note that the problems in the equations (15) and (17) are equivalent only when $\lambda_{p}$ is arbitrarily large. However, setting $\lambda_{p}$ to an arbitrarily large value will lead to numerical instability. Fortunately, through a series of reconstruction trials, we found that setting $\lambda_{p}=10\lambda$ was sufficient to avoid significant negative values.

2.3 Proposed Algorithm

Here, we first develop the minimization algorithm for reconstruction with $R_{h,1}$ and then describe the modifications necessary to use $R_{h,2}$. We adapt the preconditioned gradient search for minimizing the cost function. To this end, we first need to write the expression of the gradient. For notational convenience, we use $\mathbf{x}$ in the place of $\mathbf{p}_{0}$. The gradient expression is given by

[TABLE]

where

[TABLE]

Here $\mathbf{W}^{(\mathbf{x})}$ and $\mathbf{N}^{(\mathbf{x})}$ are diagonal matrices with $j$th diagonal element given by

[TABLE]

The superscript ‘$(\mathbf{x})$’ in the diagonal matrices signify their dependence on $\mathbf{x}$, which makes the operation $\mathbf{A}^{(\mathbf{x})}\mathbf{x}$ a non-linear operation.

The proposed preconditioned gradient search proceeds as follows: given a current estimate of the minimum, $\mathbf{x}^{(k)}$, we update the estimate as $\mathbf{x}^{(k+1)}=\mathbf{x}^{(k)}-\beta_{k}\mathbb{P}(\nabla J({\bf x}^{(k)},q))$, where $\mathbb{P}(\cdot)$ is the preconditioning, and $\beta_{k}$ is the step-size. To describe how the preconditioning is done, let $\mathbf{g}^{(k)}=\nabla J({\bf x}^{(k)},q)$ and $\hat{\mathbf{g}}^{(k)}=\mathbb{P}(\nabla J({\bf x}^{(k)},q))$. Then $\hat{\mathbf{g}}^{(k)}$ is determined by solving the linear system of equations $\mathbf{A}^{(\mathbf{x}^{(k)})}\hat{\mathbf{g}}^{(k)}={\mathbf{g}}^{(k)}$, where

[TABLE]

The idea behind this preconditioning is to mimic damped Newton approach as done in [42]. Here the matrix $\mathbf{A}^{(\mathbf{x}^{(k)})}$, which is the matrix involved in the gradient expression of equations (19), (20), works as an approximation for Hessian of $J(\cdot,q)$. To solve the equation $\mathbf{A}^{(\mathbf{x}^{(k)})}\hat{\mathbf{g}}^{(k)}={\mathbf{g}}^{(k)}$, we use the well-known method of conjugate gradients (CG) and denote the operation by $\hat{\mathbf{g}}^{(k)}=CG(\mathbf{A}^{(\mathbf{x}^{(k)})},{\mathbf{g}}^{(k)},\epsilon_{cg})$ where $\epsilon_{cg}$ denotes the termination tolerance for CG iteration. Next, the step size $\beta_{k}$ is chosen such that $\frac{J({\bf x}^{(k+1)},q)}{J({\bf x}^{(k)},q)}<\epsilon_{s}$ by means of a back-tracking procedure, where $\epsilon_{s}$ is an another user-defined tolerance for back-tracking. Specifically, starting with $\beta_{k}=1$, the required $\beta_{k}$ is determined by series of checks on the condition $\frac{J({\bf x}^{(k+1)},q)}{J({\bf x}^{(k)},q)}<\epsilon_{s}$ with iterated multiplication of $\beta_{k}$ with a factor $\rho\in(0,1)$. Finally, the iterative update on $\mathbf{x}^{(k)}$ is terminated upon the attainment of the condition $\cfrac{||\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}||_{2}}{||\mathbf{x}^{(k)}||_{2}}<\epsilon_{o}$ where $\epsilon_{o}$ is a yet another user-defined tolerance. This iterative method is summarized in the panel Algorithm 1, where the first input $\mathbf{x}^{(0)}$ is an appropriate initialization. Algorithm 1 calls the back-tracking line search method described above. This line search is summarized in the panel Algorithm 2.

Since we choose $q$ to be less than $0.5$, the regularization functional is non-convex. Hence the cost function can have several local minima, and the tolerances used in the inner iterations ($\epsilon_{cg}$ and $\epsilon_{s}$) might influence the actual minimum attained by the overall algorithm. Such a minimum is likely to contain artifacts. To alleviate this problem, we adopt the well-known graduated non-convexity (GNC) approach [43]. Specifically, we introduce an outer loop in reconstruction method, where a series of $n_{s}$ regularized reconstructions are performed with sparsity indices $\{q_{m}:=0.5-m(0.5-q)/n_{s},m=0,1,\ldots,n_{s}\}$. For each sparsity index $q_{m}$, the reconstruction is denoted by $\mathbf{y}^{(m)}$, and it is obtained by calling regularized reconstruction (Algorithm 1) with the initialization set to $\mathbf{y}^{(m-1)}$, which is the reconstruction corresponding to the previous sparsity index $q_{m-1}$. Note that, for $m=0$, $q_{0}=0.5$, and for this value of $q_{0}$, the cost is convex and hence the parameters $\epsilon_{cg}$ and $\epsilon_{s}$ do not affect the reconstruction obtained ($\mathbf{y}^{(0)}$). For subsequent calls for reconstruction with iteration index $m>0$, since the sparsity index $q_{m}$ is typically close to the index of previous reconstruction, $q_{m-1}$, the reconstruction $\mathbf{y}^{(m)}$ is likely to be close to $\mathbf{y}^{(m-1)}$. As a result, the sequence $\{\mathbf{y}^{(m)},m=0,1,\ldots,n_{s}\}$ is mostly determined by the sequence $\{q_{m},m=0,1,\ldots,n_{s}\}$. This will reduce the dependence of the final reconstruction on the parameters $\epsilon_{cg}$ and $\epsilon_{s}$ and will lead a good quality final reconstruction. For $m=0$, we need an efficient initialization for calling the regularized reconstruction. We propose to use the result of minimizing $J({\bf\cdot},q)$ with $q=1$. With $q=1$, $J({\bf\cdot},q)$ is quadratic and its minimum is determined by the following linear system of equations:

[TABLE]

We solve this again by calling the conjugate gradient method. The overall algorithm is given in the Algorithm 3. Note that the possibility of using the GNC method is the main advantage of replacing the logarithm used in the original formulation [34] by the fractional power $q$.

Now, we will specify the modifications required if $R_{h,2}$ is used. It is only required to change the operator $\mathbf{A}^{(\mathbf{x})}$ present in the equation (20). The operator corresponding to the regularization $R_{h,2}$ can be expressed as

[TABLE]

where $\bar{\mathbf{W}}^{(\mathbf{x})}$ and $\hat{\mathbf{W}}^{(\mathbf{x})}$ are diagonal matrices with its elements given by

[TABLE]

[TABLE]

2.4 Matrix-free implementation of the forward model

Finding the gradient $\mathbf{g}^{(k)}$ of the cost function and the search direction $\hat{\mathbf{g}}^{(k)}$ using CG involves the repeated application of $\textbf{H}^{T}\textbf{H}$; this is computationally very expensive and needs large amount of memory. For example, generating H for image size $512\times 512$ corresponding to $128$ transducers with each transducer taking $1024$ samples results in $10^{5}$ elements in $\mathbf{H}$, and thus requires 256 GB RAM [31]. In this work, we derive a formula for matrix free implementation of the proposed algorithm by doing on-the-fly computation of multiplication with $\textbf{H}^{T}\textbf{H}$ without explicitly constructing $\mathbf{H}$.

We first note that $\mathbf{H}$ is of the form $\mathbf{H}=[\mathbf{H}_{1}^{T}\;\mathbf{H}_{2}^{T}\;\cdots\mathbf{H}_{M}^{T}]^{T}$, where $\mathbf{H}_{i}$ represents the operation of obtaining samples from $L$ transducers at $i$th time instant. Let $\mathbf{v}$ be some vector that undergo multiplication by $\mathbf{H}$. Let $v(\mathbf{r})$ be the image obtained by putting elements of $\mathbf{v}$ into image form. Then the operator that is the equivalent of multiplying $\mathbf{v}$ with $\mathbf{H}_{i}$ can be written by using the forward model of the equation (4) as

[TABLE]

where $\mathcal{P}$ is the operation that represents retrieval of samples from transducer locations $\{\mathbf{r}_{s},s=1,\ldots,L\}$ into vector form, and $\mathcal{A}$ represents the operation of assembling the $N\times 1$ vector into an $N_{x}\times N_{y}$ image. Then, multiplication with $\mathbf{H}_{i}^{T}$ can be expressed as

[TABLE]

where $\mathcal{A}^{a}$ and $\mathcal{P}^{a}$ are the adjoints of $\mathcal{A}$ and $\mathcal{P}$ respectively. $\mathcal{A}^{a}$ represents scanning a $N_{x}\times N_{y}$ image into a vector of size $N\times 1$ ($N=N_{x}N_{y}$). Further, assuming that the transducer locations $\{\mathbf{r}_{s},s=1,\ldots,L\}$ are in a subset of image grid points, $\mathcal{P}^{a}$ becomes the operation of embedding an $L\times 1$ vector into an zero-image of size $N_{x}\times N_{y}$. Then, the operator-equivalent of $\mathbf{H}_{i}^{T}\mathbf{H}_{i}$ can be expressed as

[TABLE]

It can be shown that $\mathcal{P}^{a}\mathcal{P}$ becomes equivalent to multiplication by a binary image of $1$s and [math]s, with $1$s corresponding to the transducer locations. Let $S(\mathbf{r})$ be this binary image.

[TABLE]

This results in the following form for $\mathbf{y}=\mathbf{H}^{T}\mathbf{H}\mathbf{v}$

[TABLE]

In actual implementation, we do not build vectors. We used vector-matrix notation in the description given above only for notational convenience. Every vector involved in the algorithm is kept in the image form. Specifically, the above operator’s input will be in image form, and the output will also be in image form. If $v(\mathbf{r})$ is the $N\times N$ image equivalent of the $N^{2}\times 1$ vector $\mathbf{v}$ , and $y(\mathbf{r})$ is the image equivalent of the vector $\mathbf{y}$ the above operator can be represented by

[TABLE]

3 Reconstruction results

3.1 Reconstruction from simulated data

We used three numerical phantoms that are commonly used to evaluate PAT reconstruction methods viz. Blood vessel, Derenzo, and PAT, which is given in Figure 1. All were normalized to the size $128\times 128$ with a corresponding physical size of 12.8 mm $\times$ 12.8 mm as a model for generating the synthetic data.

The synthetic data was generated as per the geometry given in Figure 2, where the dotted circle, whose radius is 12 mm, denotes the trajectory of possible locations for transducers. The required image (image to be reconstructed) was defined to be on a larger grid of size $512\times 512$ with the equivalent physical size of 51.2 mm $\times$ 51.2 mm. This was done to accommodate for the boundary effect caused by Fourier based convolutions represented in Equation (34).

The forward data was generated using Equation (4) and added with Gaussian noise to form simulated data having SNR levels of 20 dB, 30 dB, and 40 dB. Values for $t$ in the model of Eq. (4) was chosen as $\{\delta_{t}i,i=1,\ldots,M\}$ with $1/\delta_{t}=100$ MHz and $M=1600$. The number of transducers to mimic the limited data scenario in our experiments was taken as 16, 32, 64, and 128. Note that most reported methods in the literature typically use not less than 128 transducers. The sound speed in the medium was assumed to be $1.5\;mm/\mu s$, and we considered the medium to be homogeneous with no dispersion or absorption of sound. We choose second-order derivatives for the regularization, i.e., we set $o=2$ in the equation (13) and (14). The parameter $\lambda$ was determined using the model itself, as done in most methods that focus on the design of regularization. Next, for positivity, we set $\lambda_{p}=10\lambda$. Further, we found that $q=0.25$ with $n_{s}=10$ was adequate for all cases. All the tolerances ( $\epsilon,\epsilon_{s},\epsilon_{cg}$, and $\epsilon_{o}$) were set to $10^{-6}$. Further, for all test cases, setting $\alpha=0.5$ gave good quality reconstructions. For the blood vessel phantom, we conducted experiments with the number of transducers set to 64 and 128. For the remaining phantoms, we did experiments with 16 and 32 as the number of transducers. For a given $q$, each reconstruction takes about 38 mins for the proposed method in Dell Precision T7820, whereas the FISTA-based methods take about 30 mins.

Considering the choice of method to compare, we first note that derivative based regularization methods are known to deliver robust and consistent performance in most inverse problems and hence we restrict to this class for choosing the method to compare. These methods can be categorized into Tikhonov filtering methods [21] and total variation methods [20, 25, 26, 27, 29, 39], and the later ones are more robust. Among the total variation methods, the method of Boink et al [27] use total generalized variation (TGV) and others use first order total variation (TV-1). Although TGV is known to be superior in its formulation, it poses numerical difficulties in implementing the reconstruction. This narrows down our choice of method to compare to be within TV-1 methods. All TV-1 methods [20, 25, 26, 29, 39] differ only by speed and we chose the method of Huang et al [25] because of its ease of implementation. We compare the proposed method using structural similarity (SSIM) index [44] as it is more sensitive to reconstruction artifacts.

In the first experiment, we generated measurements corresponding to 128 transducers with an SNR of 20 dB. Figure 3 compares the reconstructed results obtained from this measurement set. The proposed method with both regularization forms yields better reconstruction over the FISTA-based method, although the result with the second form of regularization gives a slightly inferior result. This is reflected by the relative improvement in the SSIM score (0.057 and 0.049). The proposed method was also able to reconstruct the peak amplitude of the initial pressure distribution (1 Pa) more accurately than the FISTA algorithm as evident from the Figures 3.A, 3.B, 3.C and 3.D. Figure 4 shows the scanline based intensity profile of reconstructions in Figure 3. The intensity profiles show that our method closely follows the ground truth compared to the FISTA based method, which is also confirmed by superior SSIM scores. Next, we consider two more noise levels viz. 30 dB and 40 dB. For this phantom, since the first form of proposed regularization outperforms the second, we present the reconstruction results only for the first form for these additional input SNR cases. The relative improvement in our SSIM scores over the FISTA-based method corresponding to 30 dB and 40 dB are 0.026 and 0.022, respectively. The proposed method yielded an SSIM score of 0.7622 for a 20 dB noise level, whereas the FISTA-based produced 0.7052. This indicates that the proposed method retrieves most of the quality that can be obtained from nearly noise-free data from the noisy data. On the other hand, the improvement in reconstruction yielded by the FISTA-based method concerning an increase in the input SNR is more gradual, as indicated by the SSIM scores 0.7052, 0.7370, and 0.7427 corresponding to input SNRs 20 dB, 30 dB, and 40 dB. Further, it also clear that the difference in performance between the proposed and FISTA-based methods is higher for lower input SNR. Next, we repeated the above trials with 64 transducers instead of 128. The relative improvements in the proposed method’s SSIM scores with the first regularization form are 0.118, 0.039, and 0.019 for the corresponding input SNRs 20 dB, 30 dB, and 40 dB.

In the next experiment, we generated measurements using a Derenzo phantom image corresponding to 32 transducers with an SNR of 20 dB. Figure 5 compares the reconstructed results obtained from this measurement set. Here too, the proposed method with both forms of regularization yields a better reconstruction over the FISTA-based method as seen from the Figures 5.A.1, 5.B.1, 5.C.1, and 5.D.1. Figures 5.A.2, 5.B.2, 5.C.2, and 5.D.2 display a cropped region from the images displayed in 5.A.1, 5.B.1, 5.C.1, and 5.D.1 for a closer view. The improvement seen in the displayed images is also reflected by the relative improvement in the SSIM score (0.074 and 0.08). Next, we considered two other noise levels, viz., 30 dB, and 40 dB, for the same phantom with the same number of transducers (32). For these cases, we present the reconstruction result only for the second form since the second form of proposed regularization outperforms the first for the current phantom. The relative improvement in the SSIM scores of the reconstructions from data sets with input SNRs 30 dB and 40 dB is 0.017 and 0.01. This again proves that the proposed method retrieves most of the quality obtained from nearly noise-free data from the noisy data. The performance between the proposed and FISTA is higher for lower input SNR in this experiment also.

Next, we repeated the same set of trials with 16 transducers instead of 32. The proposed method outperformed the FISTA method for all SNRs levels, in this case, also as shown in Figure 6. The quality of the reconstruction for the FISTA method was heavily degraded for all SNRs (20 dB, 30 dB, and 40 dB), and the reconstructed images contained artifacts, as seen from Figure 6.D, Figure 6.E and Figure 6.F respectively. However, the reconstruction quality was not much affected by the proposed method, although the number of transducers used is very less. The proposed method’s improved performance in comparison with the FISTA method is also reflected in the large difference of SSIM scores (0.301, 0.272, and 0.259) for all three input SNR levels: 20 dB, 30 dB, and 40 dB, respectively. These results demonstrate the robustness of our algorithm in the limited data scenario. The above observations were further verified by repeating the same pattern of experiments on the PAT phantom. Reconstruction results for this phantom and the results of all other trials reported in this section are displayed Table 1. For a visual comparison of performance evaluated using PAT phantom, we display the images corresponding to trials with 16 transducers in Figure 7, which confirms the superiority of the proposed method.

3.2 Reconstruction from real measured data

We use a triangular-shaped physical phantom constructed using horsehair to generate real measured data. The details of the experimental setup can be found in Fig.2 of Ref. [45]. A Q-switched Nd: YAG laser operating at 532 nm delivered laser pulses with 5 ns width at a 10 Hz repetition rate onto the sample. One uncoated plano-concave lens (LC1715, Thorlabs) and four right-angle uncoated prisms (PS911, Thorlabs) were utilized to provide a light fluence of 9 mJ/cm2 ($<$ 20 mJ/cm2: ANSI safety limit). The hair phantom having the side-length and diameter of 10 and 0.15 mm, respectively, was attached to the pipette tips adhered to an acrylic slab. For recording the PA data, a 2.25 MHz flat ultrasound transducer (Olympus-NDT, V306-SU) of 13 mm diameter active area and $70\%$ nominal bandwidth was rotated continuously for 360 deg around the sample. A pulse amplifier (Olympus-NDT, 5072PR) first amplified and filtered the acquired PA signals, and then a data acquisition (DAQ) card (GaGe, compuscope 4227) recorded the signals using a sampling frequency of 25 MHz. A sync signal from the laser was used for the synchronization of data acquisition with laser illumination.

The reconstructed PA imaging region has a size of 12.8 mm by 12.8 mm containing 128 by 128 pixels, and data from 400 transducer positions are used to do the reconstruction. Since the actual values of the initial pressure rise are unknown here, we have used the following figure-of-merit as used in the reference [46] to compare different methods.

[TABLE]

where $S$ is the peak intensity value of the reconstructed image and $n$ is the standard deviation of the intensity.

Figure 8 compares reconstructed images for horse hair phantom data. For the proposed method, since we do not have the ground truth for evaluation, we present reconstruction results using both regularization forms. Figure 8.A shows the image reconstructed using the first form of regularization term, and Figure 8.B shows the result obtained by using the second form regularization. Figure 8.C shows the corresponding result obtained from the FISTA based method. The proposed methods were able to give sharp images while giving a 3dB improvement in the FOM values compared with the FISTA based result, as given in the figure.

4 Conclusions

A novel model-based method that can reconstruct high-quality PAT images from significantly reduced measurements was developed. The joint sparse nature of intensity and derivatives in PAT images was exploited to construct the regularization. The reconstruction was formulated as the minimizer of this regularization along with data fidelity cost and a cost that quantifies the deviation from positivity. A novel computational method was developed to carry out the minimization. Further, a new computational formula was derived for implementing the forward model and was used in the data fidelity term to reduce the memory requirement. The algorithm was compared against the FISTA based method of Huang et al. [25] for various levels of reduction in the measured data (16, 32, 64, and 128 transducers) with various levels of measurement noise (20, 30, and 40 dB). We considered both real and simulated data sets for our experiments, and the proposed method yielded superior reconstruction quality in all cases.