Sparsity-based photoacoustic image reconstruction with a linear array transducer and direct measurement of the forward model

Ruibo Shang; Richard Archibald; Anne Gelb; Geoffrey P. Luke

doi:10.1117/1.JBO.24.3.031015

8 December 2018 Sparsity-based photoacoustic image reconstruction with a linear array transducer and direct measurement of the forward model

Ruibo Shang, Richard Archibald, Anne Gelb, Geoffrey P. Luke

Author Affiliations +

Journal of Biomedical Optics, Vol. 24, Issue 3, 031015 (December 2018). https://doi.org/10.1117/1.JBO.24.3.031015

Abstract

Photoacoustic (PA) imaging is an emerging imaging technique for many clinical applications. One of the challenges posed by clinical translation is that imaging systems often rely on a finite-aperture transducer rather than a full tomography system. This results in imaging artifacts arising from an underdetermined reconstruction of the initial pressure distribution (IPD). Furthermore, clinical applications often require deep imaging, resulting in a low-signal-to-noise ratio for the acquired signal because of strong light attenuation in tissue. Conventional approaches to reconstruct the IPD, such as back projection and time-reversal, do not adequately suppress the artifacts and noise. We propose a sparsity-based optimization approach that improves the reconstruction of IPD in PA imaging with a linear array ultrasound transducer. In simulation studies, the forward model matrix was measured from k-Wave simulations, and the approach was applied to reconstruct simulated point objects and the Shepp–Logan phantom. The results were compared with the conventional back projection, time-reversal approach, frequency-domain reconstruction, and the iterative least-squares approaches. In experimental studies, the forward model of our imaging system is directly measured by scanning a graphite point source through the imaging field of view. Experimental images of graphite inclusions in tissue-mimicking phantoms are reconstructed and compared with the back projection and iterative least-squares approaches. Overall these results show that our proposed optimization approach can leverage the sparsity of the PA images to improve the reconstruction of the IPD and outperform the existing popular reconstruction approaches.

1. Introduction

Photoacoustic (PA) imaging is an emerging biomedical imaging technique that combines optical contrast with ultrasound acquisition and image reconstruction to achieve good resolution and high-contrast deep tissue imaging.¹^–³ In PA imaging, the acoustic signal is generated from tissue optical absorption following irradiation of the tissue by a short laser pulse. Then the signal is detected by an ultrasound transducer to generate the radio-frequency (RF) data. This process can be described mathematically by an optical forward model of light transport and absorption followed by the acoustic forward model of ultrasound wave propagation. Therefore, the acoustic inversion and optical inversion are two independent tasks needed to form a PA image. The acoustic inversion is used to reconstruct the initial pressure distribution (IPD) in the tissue immediately following laser excitation, whereas the optical inversion is used to estimate the optical fluence distribution for a given IPD.⁴ In the acoustic inversion problem, which is the focus of this paper, factors such as acoustic diffraction, spatial variance, artifacts, and weak acoustic signal will all deteriorate the reconstructed IPD.⁵^,⁶ Artifacts in PA imaging commonly come from the finite aperture effects of the ultrasound transducer. The finite aperture effect will degrade the resolution and the signal-to-noise ratio (SNR) of the reconstructed IPD image and can also introduce side-lobe effects in the image.⁷ The weak PA signals are due to the limited penetration depth of the laser light in the tissue, which leads to poor SNR in deeper imaging regions.⁸^–¹⁰

Currently, many research groups are investigating effective approaches to address these issues. On the hardware side, a high-numerical-aperture-based virtual point detector was invented for PA tomography (PAT).¹¹ The virtual point detector detects over a wide acceptance angle so that the sensitivity is higher, and the finite-aperture effect is smaller. This concept, however, cannot be readily applied to a linear array geometry. A synthetic aperture PA-guided ultrasound technique was introduced for artifact reduction in PA imaging, where the backpropagation is applied synthetically based on the ultrasound pulse-echo acquisitions.¹²^,¹³ This helps to reduce clutter but does not address finite aperture or low-SNR effects. On the software (algorithms) side, several model-based reconstruction algorithms have been developed.⁴ These include model-based mean squared error minimization¹⁴ and non-negative optimization¹⁵ reconstruction for full-angle tomography systems, which do not address finite-aperture effects. In addition, a model-based correction of the finite-aperture effect in PAT was demonstrated by introducing a spatiotemporal optimal filter to minimize the reconstruction error in the mean square error sense, which can improve the degraded tangential resolution for PAT with finite-aperture unfocused transducers while maintaining the radial resolution.¹⁶

Sparsity-based optimization is a powerful tool in solving inverse problems for image reconstruction in many fields with the prior knowledge that the image to be reconstructed is sparse in a certain domain.¹⁷^–¹⁹ Researchers have applied the sparsity-based optimization in PA image reconstructions. A sparsity-based acoustic inversion algorithm was developed in cross-sectional multiscale PAT systems using different types of sparsity regularizations.²⁰ The forward model was built based on the image grid and the measurement geometry in a PAT system with the assumption of a lossless dispersion-free homogeneous acoustic medium. $L_{1}$ , total variation (TV), combined TV- $L_{1}$ regularizations were compared in both simulations and experiments. A deconvolution-based PA reconstruction was invented with sparsity regularization in a three-transducer PAT system.²¹ The algorithm is a semianalytical approach, where the sparsity regularization improves the numerical conditioning of the system of equations and reduces the computation time of the deconvolution-based process. A full-wave iterative image reconstruction was developed based on TV regularization in PAT system with acoustically inhomogeneous media.²² The forward and back-projection operators in this algorithm are based on the k-space pseudospectral method to compute numerical solutions to PA wave equation in the time domain with a circular ultrasound transducer geometry for data collection in the experiments. This algorithm is still an analytical approach utilizing the numerical solutions to PA wave equation to build the forward and backward model in a circular detection geometry instead of using a linear array ultrasound transducer. A sparsity-based image reconstruction algorithm with compressed sensing was developed in PAT system with both circular and linear array detecting geometry to effectively accelerate the data acquisition. They used a way to build the forward model where each element of the forward model matrix stores the index of the temporal samples of each measurement as a projection from the initial pressure to the velocity potential measurement.²³ This forward model is also built analytically using the pressure equation with the assumption of a universal delta pulse heating function (i.e., the bandwidth was assumed not to be considered). A sparsity-based reconstruction was developed for super-resolved limited-view PA computed tomography in a scattering medium.²⁴ The reconstructed images of a wire in the experiment were shown to have super resolution compared with the delay- and sum-beamforming approach with a dictionary matrix as the forward model. Using an initial, imprecise dictionary, the algorithm yielded a small cluster of points for the wire. Then the locations of these cluster points were used to generate a superposition of responses as a single-experimental dictionary element for a fixed wire location, and the dictionary elements for other wire locations were generated using appropriate delays and amplitude scaling. This approach is semianalytical since the theoretical pressure equation was used for the dictionary matrix of the cluster points, and the dictionary matrices for other wire locations were only generated with appropriate delays and amplitude scaling instead of direct measurement.

All of these model-based PA image reconstruction algorithms designed the forward model matrix from the pressure field equations, the image grid, and the experimental setup geometry. This was done analytically or semianalytically with little consideration to nonidealities of PA imaging such as the spatial variance of the impulse response, the real attenuation effects in biological tissues in experiments, and variations in element response. Furthermore, most of these model-based algorithms are based on tomography systems with detectors at multiple angles for projections, which often is not feasible for clinical applications. Therefore, in this paper, we propose a sparsity-based-photoacoustic image reconstruction (SPAIR) to optimize the IPD image using the two-step iterative shrinkage/thresholding (TwIST) algorithm²⁵ and a linear array transducer system (rather than tomography systems) with direct measurement of the spatial- and element-variant impulse response to construct the mathematical forward model matrix.

2. Method

2.1.

Mathematical Forward Model

First, we built a mathematical forward model mapping the IPD to the measured RF data. Assuming the time delay for the acoustic signal from the $k$ ’th pixel in the IPD image to reach the $j$ ’th element of the ultrasound transducer is $τ (k, j)$ . Then the imaging forward model can be written as

Eq. (1)

R (t, j) = \sum_{k} P [t - τ (k, j), j, k] \times f (k),

where

f (k)

is the initial pressure in the

k

’th pixel of the vectored reconstructed IPD,

P (t, j, k)

is the spatial- and element-variant impulse response of the

j

’th transducer element from the

k

’th pixel of the vectored reconstructed IPD, and

R (t, j)

is the corresponding acquired RF data.

For simplicity, we can formulate Eq. (1) into a matrix multiplication model as $R = Ψ f$ , where $f$ is the vectored IPD to be reconstructed, $R$ is the vectored measured RF data, and $Ψ$ is the measurement matrix (forward model) mapping $f$ to $R$ .

By taking the spatial- and element-variant properties of the impulse response into consideration in real cases, it might not be reasonable to simply build the forward model matrix analytically or semianalytically based only on the pressure field equation and the experimental geometry, as is commonly done. Instead, we directly measured the impulse response to construct the forward model matrix. First, we assume that the reconstructed IPD image has $N \times N$ pixels in the fixed field of view (FOV). Since $f$ in Eq. (1) is the vectored IPD, we then vectorize the IPD image into $N^{2} \times 1$ . We let the $n$ ’th ( $n \in [1, N^{2}]$ ) point in the vectored IPD image be the sole nonzero pixel, then we measured its impulse response (the RF data) using a linear array transducer and vectorize it as well (suppose the RF data of the impulse response has a size of $M_{1} \times M_{2}$ , where $M_{1}$ is the total number of samples and $M_{2}$ is the element number, then the size of the vectored impulse response is $M_{1} M_{2} \times 1$ ). Finally, we fill the $n$ ’th column of the forward model matrix with this vectored RF data of the impulse response corresponding to the $n$ ’th point in the vectored IPD image and repeat the same process for all $N^{2}$ columns. Therefore, the forward model matrix has a size of $M_{1} M_{2} \times N^{2}$ . Figure 1 shows a subset of the forward model matrix for simulations collected in the k-Wave toolbox.²⁶ The simulated IPD image has $128 \times 128 pixels$ and each RF data of the impulse response has a size of $500 \times 128$ , resulting in a forward model matrix of size $64,000 \times 16,384$ . To save memory, a thresholded sparse matrix representation was used to compress the matrix into $\sim 7 GB$ .

Fig. 1

Structure of the modeled forward model matrix: (a) 2-D plot of subset of the forward model matrix and (b) zoom-in picture of one peak in (a). The waveform in (b) is the impulse response of a single transducer element to a single nonzero pixel.

Figure 1(a) shows a subset of the forward model matrix for the points with pixel numbers from 61st to 70th in the 65th column of the IPD image and the vectored sampling number of the RF data from 39,500 to 44,500. Figure 1(b) shows zoom-in picture of one peak in Fig. 1(a). It corresponds to the impulse response of the point in the 65th column and 65th row from the 84th transducer element.

2.2.

Sparsity-Based Optimization

Since the size of the forward matrix is large, it is infeasible to calculate the pseudoinverse matrix for image reconstruction. Additionally, due to the finite-aperture effects and the weak PA signal, the conventional backprojection or time-reversal approaches will result in artifacts and errors. Therefore, additional information about the IPD image is needed for more accurate reconstruction. Here, we assume that the IPD image is sparse when it is transformed to a certain domain as the additional information. We expect to get a better reconstruction of the IPD image with fewer or even no artifacts and errors in a reasonable time period. Thus with the forward model matrix, we propose to solve the following sparsity-based optimization problem to find the optimized IPD image using the TwIST algorithm:

Eq. (2)

\hat{f} = \underset{f}{argmin} \frac{1}{2} {‖ Ψ f - R ‖}_{2}^{2} + λ {‖ ϕ f ‖}_{1},

where

Ψ

is the forward model matrix mapping the IPD to the collected RF data,

ϕ

is the sparse operator on

f

; (here

ϕ

is chosen to be the identity matrix or TV operator with the knowledge that IPD is usually sparse in spatial or TV domains), and

λ

is the heuristically derived regularization parameter.

{‖ Ψ f - R ‖}_{2}^{2}

is the fidelity term indicating how well the reconstructed IPD is matching to the acquired RF data through

Ψ

. The

{‖ ϕ f ‖}_{1}

term promotes sparsity in the

ϕ

domain. By adjusting

λ

, we can find a good balance between data fidelity term and image sparsity.

We used a specific algorithm TwIST²⁵ to achieve SPAIR. It solves a class of problems in the fields of image reconstruction and other linear inverse cases, combining a linear observation model with a nonquadratic regularizer (TV- and wavelet-based regularization, etc.), exactly the problem we would like to solve in Eq. (2). It has been proven that for a variety class of nonquadratic convex regularizers, TwIST converges to a minimizer of the objective function with a fast convergence rate even for ill-conditioned problems. The iterative process is shown in Eqs. (3) and (4):

Eq. (3)

f_{1} = ϕ_{λ} [f_{0} + Ψ^{T} (R - Ψ f_{0})],

Eq. (4)

f_{t + 1} = (1 - α) f_{t - 1} + (α - β) f_{t} + β ϕ_{λ} [f_{t} + Ψ^{T} (R - Ψ f_{t})],

where

f_{0}

is the initial guess of the IPD image,

f_{t}

is the reconstructed IPD image at the

t

’th iteration,

ϕ_{λ}

is the sparse operator on

f

depending on

λ

,

α

and

β

are the selected constants defined in the TwIST paper, and the other variables have the same definition as Eq. (2).

3. Simulations

3.1.

Simulations on Point Objects

To demonstrate the effectiveness of our proposed approach, we simulated the reconstruction of several point objects, as shown in Fig. 2(a). The simulation was done in the k-Wave toolbox. The simulated ultrasound transducer has 128 elements with a pitch of 0.15 mm and a sampling frequency at 22 MHz. It has a center frequency of 6 MHz and a bandwidth of 4.8 MHz. For all the simulations and experiments, the ultrasound transducer is placed horizontally on top of the object to be imaged. The size of the reconstructed IPD image is $19.2 mm \times 19.2 mm$ with $128 \times 128 pixels$ . Then we used the k-Wave MATLAB toolbox to generate the corresponding RF data as shown in Fig. 2(g) and formulated the imaging forward model $Ψ$ as described in Sec. 2.1. With the RF data and the forward model matrix, we constructed the optimization problem in Eq. (2) to reconstruct the IPD. Here the regularization parameter $λ$ is chosen to be $3.55 \times 10^{- 4}$ , and $ϕ$ is chosen to be the identity matrix with the assumption that this point object is sparse in spatial domain. For comparison, we also used the conventional back projection, time-reversal, frequency-domain, and an iterative $L_{2}$ norm minimization approach (LSQR) for IPD image reconstruction.²⁶^,²⁷ The reconstructed IPD images are shown in Figs. 2(b)–2(f), respectively.

Fig. 2

Simulation of a six-point object: (a) ground truth of the object, (b) the reconstructed IPD image from back projection approach, (c) the reconstructed IPD image from time-reversal approach, (d) the reconstructed IPD image from frequency-domain approach, (e) the reconstructed IPD image from LSQR, (f) the reconstructed IPD image from SPAIR, (g) the corresponding RF data, (h) line plots of the white dotted lines in axial direction in (a) for all the five reconstructed IPD images and the ground truth, and (i) line plots of the white dotted lines in lateral direction in (a) for all the five reconstructed IPD images and the ground truth. The scale bar is 5 mm.

In order to quantitatively compare the reconstruction results from SPAIR and other reconstruction approaches, we used full-width at half-maximum (FWHM) as an indication of the resolution of the reconstructed IPD image. Here we calculated the FWHM of the bottom point in the image both in axial and lateral directions as shown in the white dotted lines in Fig. 2(a). Fig. 2(h) shows the line plots of the point in the axial direction for all the five reconstructed IPD images and the ground truth, respectively. According to FWHM, the axial resolution is 0.16 mm for SPAIR, and 0.17, 0.18, 0.16, 0.16 mm for back projection, time-reversal, frequency-domain, and LSQR approaches, respectively. Similarly, Fig. 2(i) shows the line plots of the point in the lateral direction for all the five reconstructed IPD images and the ground truth. By calculation, the lateral resolution is 0.16 mm for SPAIR, and 0.62, 0.61, 0.62, and 0.47 mm for back projection, time-reversal, frequency-domain, and LSQR approaches, respectively. The image resolution in axial direction is already good enough from back projection, time-reversal, frequency-domain, and LSQR approaches, so there is not much improvement with SPAIR. However, SPAIR improves the lateral resolution by 66% compared with other approaches. The unequal improvement of the resolution in the two directions results from the fact that the lateral resolution is blurred by the finite aperture, whereas the axial resolution is largely unaffected by it.

We also calculated the averaged mean square error (AMSE) defined in

Eq. (5)

AMSE = \frac{\frac{1}{N^{2}} \sum_{i = 1}^{N^{2}} {[I (i) - G (i)]}^{2}}{[\frac{1}{N^{2}} \sum_{i = 1}^{N^{2}} | G (i) |]},

where

I

is the vectored reconstructed IPD,

G

is the vectored ground truth,

N

is the pixel number of the IPD image in one dimension, and

i

is the pixel number from 1 to

N^{2}

.

The AMSE indicates the reconstruction accuracy compared with the ground truth. From calculation, the AMSE for SPAIR is $1.4 \times 10^{- 5}$ while the AMSEs are 0.91, 0.94, 0.85, and 0.57 for back projection, time-reversal, frequency-domain, and LSQR approaches, respectively.

3.2.

Robustness to Noise

In order to test if our proposed approach is robust to noise, we added white Gaussian noise to the RF data with different SNR levels. Then we used back projection, time-reversal, frequency-domain, LSQR, and SPAIR approaches for reconstruction and calculated the AMSE for comparison.

Figure 3(a) shows the AMSE values for different RF SNR values. Here the SNR is defined as the averaged difference between the signal power (dB) and the noise power (dB) from all channels of the ultrasound transducer. For each SNR value, we simulated the corresponding noisy RF dataset and used back projection, time-reversal, frequency-domain approach, LSQR, and SPAIR for image reconstruction. Figures 3(b)–3(f) show the reconstructed IPD images from the five approaches with $- 2 - dB$ noise in the RF data. Then we calculated the corresponding AMSE values using Eq. (5). Figure 3 shows that with the increase of the noise level (decrease of the SNR), the AMSE from back projection, time-reversal, frequency-domain approach, and LSQR increases significantly from $\sim 1$ to $> 30$ while the AMSE from SPAIR remains at a low level below 0.3. This indicates that using SPAIR for image reconstruction is robust to noise while the other approaches are more sensitive.

Fig. 3

(a) The AMSE of the IPD from the RF data with different SNRs for time reversal (blue), back projection (red), frequency-domain approach (yellow), LSQR (purple), and SPAIR (green). The inset shows a zoomed-in view of the proposed method. (b) The reconstructed IPD image from time-reversal approach with $- 2 - dB$ noise in the RF data, (c) the reconstructed IPD image from back projection approach with $- 2 - dB$ noise in the RF data, (d) the reconstructed IPD image from frequency-domain approach with $- 2 - dB$ noise in the RF data, (e) the reconstructed IPD image from LSQR with $- 2 - dB$ noise in the RF data, (f) the reconstructed IPD image from SPAIR with $- 2 - dB$ noise in the RF data, and (g) the corresponding noisy RF data with an SNR of $- 2 - dB$ . The scale bar is 5 mm.

Therefore, Secs. 3.1 and 3.2 indicate that our proposed approach outperforms the other four reconstruction approaches in terms of image resolution, reconstruction accuracy, and noise robustness.

3.3.

Simulations on the Shepp–Logan Phantom

Since most of the biological tissues in PA imaging have complex structures, another k-Wave simulation was done with the Shepp–Logan phantom [Fig. 4(a)]. The imaging parameters and forward model matrix were unchanged. Here the regularization parameter $λ$ is chosen to be $1.66 \times 10^{- 4}$ and $ϕ$ was chosen to be the TV operator here because the phantom is sparse in TV domain (i.e., its edges are sparse). Figure 4(g) shows the corresponding RF data collected in the k-Wave toolbox. Figures 4(b)–4(f) show the reconstructed IPD images from back projection, time-reversal, frequency-domain approach, LSQR, and SPAIR, respectively. Due to the finite aperture, the vertical edges of the phantom cannot be reconstructed in back projection, time-reversal, frequency-domain, and LSQR approaches while better reconstruction is achieved with SPAIR.

Fig. 4

Simulation on the Shepp–Logan phantom using an ultrasound transducer with 6-MHz center frequency and 4.8-MHz bandwidth. (a) Ground truth of the Shepp–Logan phantom, (b) reconstructed IPD image from back projection, (c) reconstructed IPD image from the time-reversal approach, (d) reconstructed IPD image from the frequency-domain approach, (e) reconstructed IPD image from the LSQR, (f) reconstructed IPD image from SPAIR, (g) the corresponding noisy RF data with an SNR of 18 dB, and (h) line plots of the white dotted line in (a) for all the reconstructed IPD images and the ground truth. The scale bar is 5 mm.

In order to further compare the results, the profiles across the white dotted line in Fig. 4(a) are plotted from all the reconstructed IPD images and the ground truth [Fig. 4(h)]. The curve from SPAIR matches the best with the ground truth compared with the curves from other approaches. The mismatch part between the curve from SPAIR and the ground truth results from the initial guess in the TwIST algorithm. In the TwIST, an initial guess about the reconstructed image should be given as the starting point of the optimization. Since we used the result from time-reversal approach as an initial guess, the severe mismatches part between the curve from the k-Wave iteration and the ground truth cannot be fully corrected in the TwIST.

To quantitatively compare the results, the AMSE of the reconstructed IPD images was calculated compared with the ground truth. The AMSE for SPAIR is 0.20 while it is 0.35, 0.31, 0.32, and 0.23 for back projection, time-reversal, frequency-domain, and LSQR approaches, respectively.

To further quantitatively compare the results, SNR of the reconstructed IPD image was calculated corresponding to the mean IPD value in the cyan dotted rectangular region (for the signal level calculation) and the standard deviation of the red dotted rectangular region (for the noise level calculation; here noise is defined as the artifacts on the background) in Fig. 4(a). The SNR for SPAIR is 19.6 dB while it is $- 0.4$ , 7.8, 3.2, and 12.7 dB for back projection, time-reversal, frequency-domain, and LSQR approaches, respectively. Therefore, our proposed approach outperforms the other approaches in terms of AMSE and SNR measurements.

4. Experiments

4.1.

Measurement of the Experimental Impulse Responses

As mentioned in Sec. 2, the impulse responses in PA imaging systems are spatial- and element-variant so that they cannot be reliably calculated analytically or semianalytically based on the pressure field equation in PA imaging and the image geometry. Therefore, we directly measured the impulse responses to construct the forward model matrix in the experiment.

The experimental setup is shown in Fig. 5(a). A 0.2-mm diameter graphite rod was immersed into water and the cross section of it was used as point source for impulse response measurements. In order to simulate the scattering properties in the real biological tissue, $0.035 mg / ml$ of titanium dioxide was added to water as scatterers. The graphite rod was illuminated with 7-ns laser pulses at a 10-Hz repetition rate ( $λ = 770 nm$ ) from an Nd:YAG second harmonic pumped optical parametric oscillator laser system (Phocus Mobile HE, Opotek Inc.). The PA signal from the optically absorbing graphite was acquired simultaneously by the 128 elements of a linear array ultrasound transducer with a pitch size of 0.3 mm (L11-4V, Verasonics) connected to a Verasonics vantage 256 ultrasound imaging system sampling at 22.7 MHz. A 2-D stepper motor was used to move the graphite rod in lateral and axial directions with respect to the ultrasound transducer to measure the impulse response of the point at different locations in the FOV. Ten acquisitions were averaged at each spatial location.

Fig. 5

Measurement of the impulse responses in the experimental setup. (a) The experimental setup, (b) 2-D plot of subset of the forward model matrix, and (c) zoom-in picture of one peak in (b).

Collecting the RF dataset of the impulse response for each point requires $\sim 10 s$ ; therefore, for $128 \times 128$ points, the total acquisition time is on the order of two days. Therefore, we only collected the RF dataset of the impulse response from $64 \times 64$ points in the FOV and used only 64 of the ultrasound transducer elements for RF data collection. Then we constructed the forward model matrix from the measured RF data. Figure 5(b) shows a subset of the experimental forward model matrix for the points in the 17th row with the vectored sampling number from 10,000 to 15,000. Figure 5(c) shows the zoom-in picture of one peak in Fig. 5(b).

4.2.

Experiments on Tissue-Mimicking Phantoms

Finally, an experiment was conducted to verify if our proposed approach could be applied to a real imaging system while still outperforming other conventional reconstruction approaches. Since the RF data from the experiment has noise and back projection performs better than time-reversal, frequency-domain approaches in the simulation with the noisy RF data from Fig. 3(a), only the back projection (noniterative approach) and LSQR (iterative approach) were used to compare with SPAIR. In the experiment, a gelatin phantom containing two 0.2-mm graphite rods was prepared.²⁸ In order to generate the same scattering effect as the impulse response measurement experiment, 7 mg of titanium dioxide was completely dissolved in the phantom for a concentration of $0.035 mg / ml$ .

The phantom is shown in Fig. 6(a)–6(c) for front, top, and side views, respectively. Multiple imaging planes were acquired by moving the transducer with a 2-D stepper motor in the direction of the arrow in Fig. 6(b). Then we used the measured forward model matrix and the TwIST algorithm for IPD image reconstruction. Here the regularization parameter $λ$ is chosen to be $1.11 \times 10^{- 3}$ . Figure 6(d) and (i) show two frames of the 10 averaged RF acquisitions of the phantom at two different scanning locations. Figures 6(e)–6(h), 6(j)–6(m) show the reconstruction results from the back projection, LSQR, SPAIR with TV regularization, and SPAIR with $L_{1}$ regularization, respectively at the two scanning locations. It is apparent that some artifacts and background noise in the reconstructed images have been eliminated by both SPAIR with TV and $L_{1}$ regularizations.

Fig. 6

Experimental results on tissue-mimicking phantoms with two graphite inclusions. (a)–(c) 3-D perspectives of the tissue-mimicking phantom with two graphite inclusions, (d) experimentally collected RF data at the first scanning location in (b), (e) reconstructed IPD image from the back projection approach at the first scanning location in (b), (f) reconstructed IPD image from the LSQR approach at the first scanning location in (b), (g) reconstructed IPD image from SPAIR with TV regularization at the first scanning location in (b), (h) reconstructed IPD image from SPAIR with $L_{1}$ regularization at the first scanning location in (b), (i) experimentally collected RF data at the second scanning location in (b), (j) reconstructed IPD image from the back projection approach at the second scanning location in (b), (k) reconstructed IPD image from the LSQR approach at the second scanning location in (b), (l) reconstructed IPD image from SPAIR with TV regularization at the second scanning location in (b); and (h) reconstructed IPD image from SPAIR with $L_{1}$ regularization at the second scanning location in (b). The scale bar is 5 mm.

To quantitatively compare the results, the SNR was calculated with the peak signal level to be 1 (normalized images) and the standard deviation in the white dotted rectangular region (for noise level calculation). For the first scanning location corresponding to Figs. 6(d)–6(h), the SNR is 43.8 and 60.9 dB for SPAIR with TV and $L_{1}$ regularizations, respectively, 30.4 dB for the back projection, and 31.4 dB for LSQR. For the second scanning location corresponding to Fig. 6(i)–6(m), the SNR is 55.0 dB and infinity (standard deviation of the background is 0) for SPAIR TV and $L_{1}$ regularizations, respectively, 32.0 dB for the back projection and 32.5 dB for LSQR. Thus our proposed method outperforms the back projection and LSQR in terms of the SNR measurement, artifacts, and noise reduction in the experiment.

5. Conclusions and Discussion

In conclusion, we have proposed a sparsity-based optimization approach using the TwIST for reconstruction of the IPD in PA imaging with a linear array ultrasound transducer. The impulse responses were directly measured to construct the forward model matrix to consider the real factors such as spatial- and element-variance of the impulse response. The simulation results show that our proposed approach outperforms conventional back projection, time-reversal, frequency-domain, and LSQR approaches in point object and the Shepp-Logan phantom reconstructions in terms of spatial resolution and SNR improvement. A tissue-mimicking phantom with two graphite inclusions was used for experimental verification. From the experiment, it can be demonstrated that our proposed method can be applied to a real imaging system and still outperforms back projection and LSQR.

In the experiment, although the results from SPAIR with $L_{1}$ regularization perform better than SPAIR with TV regularization, the TV regularization is more versatile. In many real cases, it is difficult to predict the features in the image. If the features are unknown, it is better to use TV since many natural objects are sparse in TV domain.

Although the SPAIR approach showed superior performance in these studies, there are some trade-offs to using the method. Currently, this approach cannot be applied to real-time applications since the postprocessing of the sparsity-based optimization takes $\sim 2.5 h$ with the normal desktop. The main reason for the long processing time is the large size of the forward model matrix, even when sparse representation is used. In the future, we will explore a more efficient way to store and extract the elements of the forward model matrix to make the optimization process faster. Another limitation of this approach is that every time we change geometries or the imaging parameters, we need to measure the impulse responses point by point. In the future, we will explore the use of deep learning approach to predict the impulse responses given geometries or the imaging parameters. In addition, the nonuniform speed of sound will introduce slight errors in terms of time of flight in the forward model. In the future, we will explore methods to incorporate phase in the forward model to deal with the errors from the nonuniform speed of sound. Overall, we envision that the method could complement back projection reconstruction (which can be performed in real time) as a postprocessing step to improve the image quality.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

Acknowledgments

This work was supported by the Department of Defense Breast Cancer Research Program Award No. W81XWH-14-1-0356 (G. L.); Department of Energy ACUMEN Program Award No. ERKJ289 (R. A.); National Science Foundation Division of Mathematical Sciences Award Nos. 1502640 and 1732434; Air Force Office of Scientific Research Award No. FA9550-15-1-0152 (A. G.); and Dartmouth College Pilot Grant (G. L. and A. G.).

References

1.

M. Xu et al., “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum., 77 (4), 041101 (2006). https://doi.org/10.1063/1.2195024 RSINAK 0034-6748 Google Scholar

2.

X. Wang et al., “Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain,” Nat. Biotechnol., 21 (7), 803 –806 (2003). https://doi.org/10.1038/nbt839 NABIF9 1087-0156 Google Scholar

3.

L. V. Wang et al., “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science, 335 (6075), 1458 –1462 (2012). https://doi.org/10.1126/science.1216210 SCIEAS 0036-8075 Google Scholar

4.

A. Rosenthal et al., “Acoustic inversion in optoacoustic tomography: a review,” Curr. Med. Imaging Rev., 9 (4), 318 –336 (2013). https://doi.org/10.2174/15734056113096660006 Google Scholar

5.

S. Vilov et al., “Overcoming the acoustic diffraction limit in photoacoustic imaging by the localization of flowing absorbers,” Opt. Lett., 42 (21), 4379 –4382 (2017). https://doi.org/10.1364/OL.42.004379 OPLEDP 0146-9592 Google Scholar

6.

G. P. Luke et al., “Biomedical applications of photoacoustic imaging with exogenous contrast agents,” Ann. Biomed. Eng., 40 (2), 422 –437 (2012). https://doi.org/10.1007/s10439-011-0449-4 ABMECF 0090-6964 Google Scholar

7.

M. K. A. Singh et al., “Photoacoustic-guided focused ultrasound (PAFUSion) for identifying reflection artifacts in photoacoustic imaging,” Photoacoustics, 3 (4), 123 –131 (2015). https://doi.org/10.1016/j.pacs.2015.09.001 Google Scholar

8.

S. Telenkov et al., “Signal-to-noise analysis of biomedical photoacoustic measurements in time and frequency domains,” Rev. Sci. Instrum., 81 (12), 124901 (2010). https://doi.org/10.1063/1.3505113 RSINAK 0034-6748 Google Scholar

9.

M. F. Beckmann et al., “Optimized SNR simultaneous multispectral photoacoustic imaging with laser diodes,” Opt. Express, 23 (2), 1816 –1828 (2015). https://doi.org/10.1364/OE.23.001816 OPEXFF 1094-4087 Google Scholar

10.

A. Hariri et al., “Development of low-cost photoacoustic imaging systems using very low-energy pulsed laser diodes,” J. Biomed. Opt., 22 (7), 075001 (2017). https://doi.org/10.1117/1.JBO.22.7.075001 JBOPFO 1083-3668 Google Scholar

11.

C. Li et al., “High-numerical-aperture-based virtual point detectors for photoacoustic tomography,” Appl. Phys. Lett., 93 (3), 033902 (2008). https://doi.org/10.1063/1.2963365 APPLAB 0003-6951 Google Scholar

12.

M. K. A. Singh et al., “In Vivo demonstration of reflection artifact reduction in photoacoustic imaging using synthetic aperture photoacoustic-guided focused ultrasound (PAFUSion),” Biomed. Opt. Express, 7 (8), 2955 –2972 (2016). https://doi.org/10.1364/BOE.7.002955 BOEICL 2156-7085 Google Scholar

13.

M. K. A. Singh et al., “Photoacoustic reflection artifacts reduction using photoacoustic-guided focused ultrasound: comparison between plane-wave and element-by-element synthetic backpropagation approach,” Biomed. Opt. Express, 8 (4), 2245 –2260 (2017). https://doi.org/10.1364/BOE.8.002245 BOEICL 2156-7085 Google Scholar

14.

A. Rosenthal et al., “Fast semi-analytical model-based acoustic inversion for quantitative optoacoustic tomography,” IEEE Trans. Med. Imaging, 29 (6), 1275 –1285 (2010). https://doi.org/10.1109/TMI.2010.2044584 ITMID4 0278-0062 Google Scholar

15.

L. Ding et al., “Efficient non-negative constrained model-based inversion in optoacoustic tomography,” Phys. Med. Biol., 60 (17), 6733 –6750 (2015). https://doi.org/10.1088/0031-9155/60/17/6733 PHMBA7 0031-9155 Google Scholar

16.

M. Li et al., “Model-based correction of finite aperture effect in photoacoustic tomography,” Opt. Express, 18 (25), 26285 –26292 (2010). https://doi.org/10.1364/OE.18.026285 OPEXFF 1094-4087 Google Scholar

17.

E. Candes et al., “Sparsity and incoherence in compressive sampling,” Inverse Prob., 23 (3), 969 –985 (2007). https://doi.org/10.1088/0266-5611/23/3/008 INPEEY 0266-5611 Google Scholar

18.

E. J. Candes et al., “An introduction to compressive sampling,” IEEE Signal Process. Mag., 25 (2), 21 –30 (2008). https://doi.org/10.1109/MSP.2007.914731 ISPRE6 1053-5888 Google Scholar

19.

I. U. Haq et al., “Sparse-representation-based denoising of photoacoustic images,” Biomed. Phys. Eng. Express, 3 (4), 045014 (2017). https://doi.org/10.1088/2057-1976/aa7a44 Google Scholar

20.

Y. Han et al., “Sparsity-based acoustic inversion in cross-sectional multiscale optoacoustic imaging,” Med. Phys., 42 (9), 5444 –5452 (2015). https://doi.org/10.1118/1.4928596 MPHYA6 0094-2405 Google Scholar

21.

H. Moradi et al., “Deconvolution based photoacoustic reconstruction with sparsity regularization,” Opt. Express, 25 (3), 2771 –2789 (2017). https://doi.org/10.1364/OE.25.002771 OPEXFF 1094-4087 Google Scholar

22.

C. Huang et al., “Full-wave iterative image reconstruction in photoacoustic tomography with acoustically inhomogeneous media,” IEEE Trans. Med. Imaging, 32 (6), 1097 –1110 (2013). https://doi.org/10.1109/TMI.2013.2254496 ITMID4 0278-0062 Google Scholar

23.

Z. Guo et al., “Compressed sensing in photoacoustic tomography in vivo,” J. Biomed. Opt., 15 (2), 021311 (2010). https://doi.org/10.1117/1.3381187 JBOPFO 1083-3668 Google Scholar

24.

D. M. Egolf et al., “Sparsity-based reconstruction for super-resolved limited-view photoacoustic computed tomography deep in a scattering medium,” Opt. Lett., 43 (10), 2221 –2224 (2018). https://doi.org/10.1364/OL.43.002221 OPLEDP 0146-9592 Google Scholar

25.

J. M. Bioucas-Dias et al., “A new TwIST: two-step iterative shrinkage/thresholding algorithms for image restoration,” IEEE Trans. Image Process., 16 (12), 2992 –3004 (2007). https://doi.org/10.1109/TIP.2007.909319 IIPRE4 1057-7149 Google Scholar

26.

B. E. Treeby et al., “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt., 15 (2), 021314 (2010). https://doi.org/10.1117/1.3360308 JBOPFO 1083-3668 Google Scholar

27.

C. C. Paige et al., “LSQR: an algorithm for sparse linear equations and sparse least squares,” ACM Trans. Math. Software, 8 (1), 43 –71 (1982). https://doi.org/10.1145/355984.355989 ACMSCU 0098-3500 Google Scholar

28.

J. R. Cook et al., “Tissue-mimicking phantoms for photoacoustic and ultrasonic imaging,” Biomed. Opt. Express, 2 (11), 3193 –3206 (2011). https://doi.org/10.1364/BOE.2.003193 BOEICL 2156-7085 Google Scholar

Biography

Ruibo Shang received his bachelor’s degree in optoelectronic science and technology from Nankai University, China, and his master’s degree in electrical engineering from Virginia Tech. Currently, he is pursuing his PhD at Thayer School of Engineering of Dartmouth College. His research focuses on computational photoacoustic and ultrasound imaging.

Richard Archibald was advised by Professor Anne Gelb and received his PhD in mathematics from Arizona State University in 2002. He is an applied mathematician in the Computational and Applied Mathematics Group at Oak Ridge National Laboratory. His research interests lie in data reconstruction and analysis, high-order edge detection, large scale optimization, time integration, and uncertainty quantification.

Anne Gelb received her PhD from the Division of Applied Mathematics at Brown University. She is the John G. Kemeny Parents Professor in the Department of Mathematics at Dartmouth College. Her research is focused on image reconstruction, feature detection, algorithm development for reconstructing magnetic resonance images, and feature extraction for synthetic aperture radar.

Geoffrey P. Luke received his PhD in electrical engineering from the University of Texas at Austin, where he developed innovative ultrasound-based molecular imaging methods and contrast agents for cancer applications. He is an assistant professor at Thayer School of Engineering of Dartmouth College and a member of the Cancer Imaging and Radiobiology Research Program at Norris Cotton Cancer Center. He directs the Functional and Molecular Imaging Research Laboratory. His current research in molecular imaging ranges from basic science to clinical translation and incorporates light, sound, and nanotechnology.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Ruibo Shang, Richard Archibald, Anne Gelb, and Geoffrey P. Luke "Sparsity-based photoacoustic image reconstruction with a linear array transducer and direct measurement of the forward model," Journal of Biomedical Optics 24(3), 031015 (8 December 2018). https://doi.org/10.1117/1.JBO.24.3.031015

Received: 30 June 2018; Accepted: 10 November 2018; Published: 8 December 2018

Access the abstract

JOURNAL ARTICLE
9 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 17 scholarly publications.

Explore citations on Lens.org

KEYWORDS

Signal to noise ratio

Transducers

Image restoration

Data modeling

Ultrasonography

Photoacoustic spectroscopy

3D image reconstruction

1.

Introduction

2.

Method

2.1.

Mathematical Forward Model

Eq. (1)

Fig. 1

2.2.

Sparsity-Based Optimization

Eq. (2)

Eq. (3)

Eq. (4)

3.

Simulations

3.1.

Simulations on Point Objects

Fig. 2

Eq. (5)

3.2.

Robustness to Noise

Fig. 3

3.3.

Simulations on the Shepp–Logan Phantom

Fig. 4

4.

Experiments

4.1.

Measurement of the Experimental Impulse Responses

Fig. 5

4.2.

Experiments on Tissue-Mimicking Phantoms

Fig. 6

5.

Conclusions and Discussion

Disclosures

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years