Vector extrapolation methods for accelerating iterative reconstruction methods in limited-data photoacoustic tomography

Navchetan Awasthi; Sandeep Kumar Kalva; Manojit Pramanik; Phaneendra K. Yalavarthy

doi:10.1117/1.JBO.23.7.071204

5 February 2018 Vector extrapolation methods for accelerating iterative reconstruction methods in limited-data photoacoustic tomography

Navchetan Awasthi, Sandeep Kumar Kalva, Manojit Pramanik, Phaneendra K. Yalavarthy

Author Affiliations +

Journal of Biomedical Optics, Vol. 23, Issue 7, 071204 (February 2018). https://doi.org/10.1117/1.JBO.23.7.071204

Abstract

As limited data photoacoustic tomographic image reconstruction problem is known to be ill-posed, the iterative reconstruction methods were proven to be effective in terms of providing good quality initial pressure distribution. Often, these iterative methods require a large number of iterations to converge to a solution, in turn making the image reconstruction procedure computationally inefficient. In this work, two variants of vector polynomial extrapolation techniques were deployed to accelerate two standard iterative photoacoustic image reconstruction algorithms, including regularized steepest descent and total variation regularization methods. It is shown using numerical and experimental phantom cases that these extrapolation methods that are proposed in this work can provide significant acceleration (as high as 4.7 times) along with added advantage of improving reconstructed image quality.

1. Introduction

Photoacoustic (PA) imaging is a noninvasive imaging modality, which gives the advantage of optical contrast at ultrasonic resolution.¹^–⁶ Application of PA imaging ranges from blood vessel imaging, tumor imaging, and biomarkers to gene activities.¹ It consists of irradiating the tissue with a nanosecond laser pulse, which gets absorbed by the tissue and generates ultrasonic waves due to thermoelastic expansion. The ultrasonic detectors placed at the boundary of the domain detect these ultrasonic waves.² The collected data gets utilized in an inverse framework resulting in the distribution of initial pressure rise.

Several image reconstruction algorithms exist for the reconstruction of PA tomographic images, such as filtered back projection, time-reversal, and Fourier transform methods.⁷^,⁸ They require sufficiently large amount of data for accurate reconstruction and may not be useful for quantitative comparisons.⁹^–¹¹ Acquisition of large amounts of data is often linked with large data-acquisition time and/or increase in the instrumentation cost. Moreover, the setups used for recording the acoustic signals cover an aperture, which may not enclose the complete object, resulting in incomplete or limited data scenarios.¹²^–¹⁴ Obtaining quantitatively accurate reconstruction using limited data has been preferred in the recent past.⁶^–⁸^,¹⁵ In limited data cases, many iterative image reconstruction algorithms were proposed, which improve the quantitative accuracy of the reconstructed images.⁹^–¹¹^,¹⁵ The direct methods (noniterative type) often have limitation in terms of not being robust to noise, and typically in the noisy data cases, the iterative type algorithms are preferred as they exhibit better stability.⁷^,⁸ Moreover, iterative algorithms often provide inherent regularization in terms of number of iterations, especially for the ill-posed problems like the one at hand. In terms of computation, typically, iterative algorithms require matrix–vector computations [ $O (n^{2})$ operations per iteration with $n$ being the number of unknowns] and direct reconstruction algorithms that utilize regularization require matrix–matrix computations [ $O (n^{3})$ operations]. Above all, the memory requirement for the iterative algorithms is at least one order less, thus enabling the solving of large problems.⁷^,⁸ Even though these iterative reconstruction methods are preferred for PA image reconstruction, often these methods require prohibitively large computational time to converge especially when the initial guess is far away from the expected one, making them less appealing in the real-time.⁷^,⁸

In this work, we address this problem by deploying a vector polynomial extrapolation method, which can accelerate any iterative scheme, thus, making iterative reconstruction methods more appealing. Vector polynomial extrapolation approach relies on finding the limit or antilimit of the initial solutions of an iterative method. Here, we have used these methods to accelerate the most widely used iterative reconstruction methods, namely regularized steepest descent (RSD)¹⁶ and total variation (TV) regularization.¹⁷^,¹⁸ Both numerical and experimental studies involving regular and irregular objects were performed to show the efficacy of the proposed method in comparison to the original iterative methods.

2. Photoacoustic Image Reconstruction

PA wave equation can be written as follows:⁶

Eq. (1)

\nabla^{2} P (x, t) - \frac{1}{c^{2}} \frac{\partial^{2} P (x, t)}{\partial t^{2}} = \frac{- β}{C_{p}} \frac{\partial H (x, t)}{\partial t},

where

P (x, t)

is the pressure at position

x

and time

t

,

c

is the speed of sound in the medium,

β

is the thermal expansion coefficient,

C_{p}

is the specific heat capacity, and

H (x, t)

represents the energy deposited per unit time per unit volume. The PA data collected at the boundary of the imaging domain can be obtained by solving Eq. (1) for time

t

. The reconstruction problem requires an estimate of the initial pressure [

P (x, t)

] at

t = 0

inside the imaging domain, given the observations at the boundary at time

t

, making the reconstruction problem equivalent to initial value problem.

The PA wave propagation model can be written in terms of linear system of equations and known to provide a simple yet effective way of solving the PA wave equation (given in Refs. 19 and 20). In this approach, the forward model for PA imaging can be written as follows:⁷

Eq. (2)

A x = b,

where

A

is the system matrix containing impulse responses of all pixels in the imaging region as columns,

x

is the unknown vector representing the initial pressure rise, and

b

is the measurement vector. The time varying data representing the impulse response recorded corresponding to each pixel is stacked as a long vector accounting for the column of the system matrix. Thus, the number of columns in the system matrix (

A

) is equal to number of pixels in the imaging domain. The computation time required for building

A

is improved by calculating the impulse response for only one corner pixel and using the shifting and attenuation properties of the PA signal to build the whole system matrix.¹⁹^,²⁰ The simple linear back-projected (LBP) image

x_{bp}

can be obtained via²⁰^,²¹

Eq. (3)

x_{bp} = A^{T} b,

where

T

represents the transpose of the matrix. Earlier works have also introduced filtered back projection as well as delay-and-sum,²² which are also noniterative as well as analytic types of reconstruction methods. Since all of them are noniterative, they are computationally efficient but only provide qualitative results along with requirement of having complete data.⁷ This

x_{bp}

is typically used as an initial guess (

x_{0}

) for all iterative methods that are discussed in this work, which is a standard practice.²²

Typically, for limited data cases, a least-squares framework can be utilized along with Tikhonov regularization, resulting in the cost function:⁸^,²²

Eq. (4)

Ω = {‖ A x - b ‖}_{2}^{2} + α {‖ x ‖}_{2}^{2},

where

α

is the regularization parameter. The

Ω

is minimized with respect to

x

and as the regularization functional is smooth in nature, it promotes smooth solution. The direct solution for this equation can be written as follows:

Eq. (5)

x = {(A^{T} A + α I)}^{- 1} A^{T} b .

The direct solution is prone to numerical instability as in most cases $A$ tends to be sparse and/or contains small values. In these scenarios, it is very common to use iterative technique to minimize the cost function given in Eq. (4).

2.1.

Regularized Steepest Descent Method

The simplest of them is the RSD, given in Algorithm 1. First, in step 2, the residual ( $r_{n}$ ) is computed, then the direction of ascent ( $L_{n}^{α}$ ) is computed in step 3, and finally, the step size ( $k_{n}^{α}$ ) is computed in step 4. The solution vector $x_{n}$ is updated in step 5. For all methods discussed in the work, the iterations are stopped when there is less than 1% change in the relative residual [for 40-dB signal-to-noise ratio (SNR) case]. The convergence criteria were varied similarly for different noise levels.

Algorithm 1

Algorithm for RSD method.

Input: $A, b, x_{0}, α$

Output: Solution vector: $x$

1 Repeat 2 to 5 till convergence is satisfied

2

r_{n} = A (x_{n}) - b

3

L_{n}^{α} = A^{T} r_{n} + α x_{n}

4

k_{n}^{α} = \frac{{‖ L_{n}^{α} ‖}^{2}}{{‖ A L_{n}^{α} ‖}^{2} + α {‖ L_{n}^{α} ‖}^{2}}

5

x_{n + 1} = x_{n} - k_{n}^{α} L_{n}^{α}

The choice of the regularization parameter ( $α$ ) is critical as it effects the reconstructed image quality (trade-off between the residual and the minimum norm solution). Usually, $α$ is decreased at every iteration ( $α_{0} > α_{1} > α_{2} > \dots > α_{n}$ ) for the RSD methods, which is also known as adaptive regularization strategy.¹⁶

2.2.

Total Variation Regularization Method

TV regularization is another popular method used for solving ill-posed linear inverse problems. In this, the cost-function is written as follows:¹⁷^,¹⁸

Eq. (6)

Ω = ‖ A x - b ‖_{2}^{2} + λ ‖ x ‖_{TV},

where

λ

is the regularization parameter, which balances the data model misfit and variation in

x

.

‖ . ‖_{TV}

represents the isotropic TV in this case, which is utilizing the same form, as in Ref. 23.

A popular method, split augmented Lagrangian shrinkage algorithm (SALSA),²⁴ was utilized here for the TV regularization, which utilizes the alternating direction method of multipliers framework given in Ref. 25, detailed steps are given in Algorithm 2. The step 2 of the algorithm has a closed form solution as given below:

Eq. (7)

x_{k + 1} = {(A^{T} A + λ I)}^{- 1} [A^{T} b + λ (v_{k} + d_{k})] .

Algorithm 2

Algorithm of SALSA.

Input: $A, b, λ$ , $k = 0, v_{0}$ and $d_{0}$

Output: Solution vector: $x$

1 Repeat 2-5 till convergence is satisfied

2

x_{k + 1} = \arg \min_{x} {‖ A x - b ‖}_{2}^{2} + λ {‖ x - v_{k} - d_{k} ‖}_{2}^{2}

3

v_{k + 1} = \arg \min_{v} τ ϕ (v) + (λ / 2) {‖ x_{k + 1} - v - d_{k} ‖}_{2}^{2}

4

d_{k + 1} = d_{k} - (x_{k + 1} - v_{k + 1})

5

k = k + 1

The step 3 of the algorithm has a solution of the form:

Eq. (8)

v_{k + 1} = ψ_{τ ϕ / λ} (x_{k + 1} - d_{k}),

where

ψ

denotes the soft thresholding function and

ϕ

denotes the TV functional

{‖ x ‖}_{TV}

. In case of TV, where it does not have a closed form solution as in Eq. (8) and needs an iterative solution, the convergence can be guaranteed by Theorem 1 in Ref. 24.

2.3.

Proposed Vector Extrapolation Methods

The iterative methods refine the current iterate and in an asymptotic sense, this sequence of iterates produced by these algorithms converge to the exact solution. Typically, these required iterations are in the range of few hundred and ideally, one would like to converge to the solution with few tens of iterations. The acceleration of convergence to achieve this can be obtained by fast cyclic iterative algorithms,²⁶ which use vector extrapolation technique. Here, each cycle of the algorithm generates a finite number of consecutive iterates by the algorithms discussed above in Algorithm 1 or Algorithm 2. These are used to obtain an extrapolated estimate of the solution. The estimated iterates are used as the next guess for the iterative algorithms. This process is repeated until convergence is achieved.²⁷^,²⁸ Among the available methods, minimal polynomial extrapolation (MPE)²⁷ and reduced rank extrapolation (RRE)²⁶^,²⁸ are utilized in this work, which are known to be the state-of-the-art extrapolation methods.

2.3.1.

Minimal polynomial extrapolation method

Polynomial extrapolation methods compute the approximation as a weighted sum of the iterates, where the weights are determined by the coefficients of the minimal polynomial $P (λ)$ of $A$ with respect to $u_{0}$ , where $u_{0}$ is defined as follows:

Eq. (9)

u_{0} = \nabla x_{1} = x_{1} - x_{0},

where

x_{0}, x_{1}, \dots, x_{n}

are the initial iterates, i.e., the unique monic polynomial will satisfy the following property:

Eq. (10)

P (A) u_{0} = 0 .

Suppose $k$ is the degree of the minimal polynomial $P (λ)$ , where $k \leq N$ . Here, $N$ denotes the dimension of the solution vector. Then, minimal polynomial can be written as

Eq. (11)

P (λ) = \sum_{j = 0}^{k} c_{j} λ^{j},

where

c_{j}

are unknown coefficients (to be determined). Combining Eqs. (10) and (11) results in

Eq. (12)

\sum_{j = 0}^{k} c_{j} u_{j} = 0,

which can be written as follows:

Eq. (13)

U_{k - 1} c = - u_{k},

where

c = {(c_{0}, c_{1}, \dots, c_{k - 1})}^{T}

and

c_{k} = 1

. Using

QR

factorization of

U_{k - 1}

, where

Q_{k - 1}^{*} Q_{k - 1} = I_{k - 1 \times k - 1}

:

Eq. (14)

Q_{k - 1} R_{k - 1} c = - u_{k} .

Multiplying both sides by

Q_{k - 1}^{*}

, it becomes

Eq. (15)

Q_{k - 1}^{*} Q_{k - 1} R_{k - 1} c = - Q_{k - 1}^{*} u_{k},

which reduces to

Eq. (16)

R_{k - 1} c = - Q_{k - 1}^{*} u_{k},

making

Q_{k - 1}^{*} u_{k} = (r_{0 k}, r_{1 k}, \dots, r_{k - 1, k}) = ρ_{k}

, simplifies to

Eq. (17)

R_{k - 1} c = - ρ_{k} .

Solution to the above linear system of equations will give $c$ , which are the unknown coefficients of the minimal polynomial. When $U$ is rank deficient, Eq. (13) reduces to

Eq. (18)

c = - U_{k - 1}^{+} u_{k},

where

U_{k - 1}^{+}

represents the pseudoinverse of

U_{k - 1}

.

Suppose, the error at every iteration, $e_{n}$ , is defined as

Eq. (19)

e_{n} = x_{n} - s,

where

s

denotes the exact solution of the system and

x_{n}

is the current iterate. From the error propagation theory for linear systems, it is well known that

Eq. (20)

e_{n} = A^{n} e_{0} .

Since $P (λ)$ is the minimal polynomial of $A$ . Therefore, $P (A) e_{0} = 0$ . Combining Eqs. (19) and (20) lead to

Eq. (21)

\sum_{j = 0}^{k} c_{j} A^{j} (x_{o} - s) = \sum_{j = 0}^{k} c_{j} (x_{j} - s) = 0,

solving for

s

results in

Eq. (22)

s = \frac{\sum_{j = 0}^{k} c_{j} x_{j}}{\sum_{j = 0}^{k} c_{j}} .

To simplify above equation, write $γ_{j} = \frac{c_{j}}{\sum_{i = 0}^{k} c_{i}}$ , where $j = 0,1, \dots, k$ with $\sum_{j = 0}^{k} γ_{j} = 1$ . Multiplying both sides of Eq. (13) by ${(\sum_{j = 0}^{k} c_{j})}^{- 1}$ converts it to

Eq. (23)

U_{k} γ = 0 .

Simplifying Eq. (22), writing it in terms of $γ_{j}$ leads to

Eq. (24)

s = \sum_{j = 0}^{k} γ_{j} x_{j} .

Since $x_{j} = x_{0} + \sum_{i = 0}^{j - 1} u_{i}$ , Eq. (24) can be rewritten as follows:

Eq. (25)

s = x_{0} + \sum_{i = 0}^{k - 1} ζ_{i} u_{i},

where

Eq. (26)

ζ_{i} = \sum_{j = i + 1}^{k} γ_{j},

which in the matrix form can be written as

Eq. (27)

s = x_{0} + U_{k - 1} ζ .

In MPE, we will define matrix $U_{k - 1} \in R^{n * k}$ as comprising of vectors $u_{0}, u_{1}, u_{2}, \dots, u_{k - 1}$ , where $u_{i} \in R^{n}$ . Let $c$ be the least square solution of Eq. (13). Start with $c_{k} = 1$ and compute $γ_{0}, γ_{1}, \dots, γ_{k}$ using

Eq. (28)

γ_{j}^{MPE} = \frac{c_{j}}{\sum_{i = 0}^{k} c_{i}},

where

j = 0,1, \dots, k

. Thus, the MPE estimate becomes

Eq. (29)

s_{k}^{MPE} = \sum_{j = 0}^{k} γ_{j}^{MPE} x_{j} .

The computation of the solution ( $s$ ) thus requires only an estimate of $γ_{j}$ ‘s and followed by simple vector–vector multiplication.

2.3.2.

Reduced rank extrapolation method

In RRE, the definition of $U$ changes to $U_{k} \in R^{n * (k + 1)}$ as comprising of vectors $u_{0}, u_{1}, u_{2}, \dots, u_{k}$ , where $u_{i} \in R^{n}$ . Similar to MPE, the RRE estimate becomes

Eq. (30)

s_{k}^{RRE} = \sum_{j = 0}^{k} γ_{j}^{RRE} x_{j} .

A unified algorithm for implementing MPE and RRE is given in Algorithm 3.²⁹ It is based on the modified QR factorization.

Algorithm 3

Algorithm listing important steps for implementing MPE and RRE.

Input: $x_{n}, x_{n + 1}, x_{n + 2}, \dots, x_{n + k + 1}$

Output: Solution vector: $s_{n, k}$

1

u_{i} = Δ x_{i} = x_{i + 1} - x_{i}

where

i = n, n + 1, \dots, n + k

2

U_{j} = [u_{n} | u_{n + 1} | \dots | u_{n + j}]

where

j = 0, 1, \dots

3

U_{k} = Q_{k} R_{k}

, QR factorization of

U_{k}

,

Q_{k}

is unitary and

R_{k}

is upper triangular

4 (MPE)

•

R_{k - 1} c^{'} = - ρ_{k}; ρ_{k} = {[r_{0 k}, r_{1 k}, \dots, r_{k - 1, k}]}^{T}, c^{'} = {[c_{o}, c_{1}, \dots, c_{k - 1}]}^{T}

•

c_{k} = 1; α = \sum_{i = 0}^{k} c_{i}

•

γ_{i} = c_{i} / α

where

i = 0,1, \dots, k

(RRE)

•

R_{k}^{*} R_{k} d = e, d = {[d_{o}, d_{1}, \dots, d_{k}]}^{T}, e = {[1, 1, \dots, 1]}^{T} \in R^{k + 1}

•

λ = {(\sum_{i = 0}^{k} d_{i})}^{- 1}

•

γ = λ d

(Common steps)

ζ = 1 - γ_{0}; ζ_{i} = ζ_{i - 1} - γ_{j}

where

j = 1, \dots, k - 1

s_{n, k} = x_{n} + Q_{k - 1} (R_{k - 1} ζ)

Algorithm 3 is called for $L$ cycles with $k$ representing the order of extrapolation. First, one of the above algorithms (Algorithm 1 or Algorithm 2) is called for $k$ steps. The generated $k$ iterates are extrapolated to get the solution ( $s$ ). This completes one cycle of the algorithm. The solution obtained becomes an input to Algorithm 1 or Algorithm 2 as an initial guess and the process is repeated till the convergence is obtained. These methods combining vector extrapolation with RSD, via Algorithm 1, will be called as minimal polynomial extrapolated regularized steepest descent (MPERSD) and reduced rank extrapolated regularized steepest descent (RRERSD). The ones where TV regularization, via Algorithm 2, will be called as minimal polynomial extrapolated total variation (MPETV) and reduced rank extrapolated total variation (RRETV).

3. Figures of Merit

The figures of merit utilized for effectively comparing the reconstruction methods introduced in this work were Pearson correlation (PC), contrast-to-noise ratio (CNR), and SNR.³⁰

3.1.

Pearson Correlation

PC is defined as follows:

Eq. (31)

PC (x^{target}, \tilde{x}) = \frac{COV (x^{target}, \tilde{x})}{σ (x^{target}) σ (\tilde{x})},

where

x^{target}

is the expected initial pressure distribution and

\tilde{x}

is the reconstructed initial pressure distribution. Here,

σ

denotes standard deviation, and COV is the covariance. It measures the correlation between the target and the reconstructed image, ranging from

- 1

to 1. The higher the value, the closer is the reconstructed image with the expected image.

3.2.

Contrast to Noise Ratio

CNR is defined as follows:

Eq. (32)

CNR = \frac{μ_{RoI} - μ_{back}}{{(δ_{RoI}^{2} a_{RoI} + δ_{back}^{2} a_{back})}^{1 / 2}},

where

μ

and

δ

represent the mean and the variance corresponding to the region of interest (RoI) and the background (back) in the reconstructed initial pressure. The

a_{RoI} = \frac{A_{RoI}}{A_{tot}}

and

a_{back} = \frac{A_{back}}{A_{tot}}

represent the area ratio.

A_{RoI}

represents the number of pixels with nonzero initial pressure distribution (here, it is the region in the target image with pixel values being 1) and

A_{back}

represents the number of pixels with zero initial pressure rise (region containing pixel values as zero in the target image).

A_{tot}

is the sum of areas of the region of interest and the background.

a_{RoI}

and

a_{back}

are the noise weights for the region of interest and the background, respectively. Higher value of CNR represents the contrast recovery of the reconstruction technique.

3.3.

Signal-to-Noise Ratio

In cases, especially experimental and in-vivo, it may not be possible to know the target initial pressure distribution to calculate both PC and CNR. In these cases, the figure of merit that has been used is SNR of the reconstructed image. It is expressed as given below:

Eq. (33)

SNR (dB) = 20 \times \log_{10} (\frac{S}{n}),

with

S

denoting the peak-to-peak initial pressure value and

n

corresponds to standard deviation. The higher the SNR represents the lesser noise in the reconstructed image and thus better performance of the reconstruction method.

4. Numerical and Experimental Phantom Studies

The acceleration achieved by the proposed vector extrapolation methods was tested using both numerical and experimental phantom cases. The measurement of actual initial pressure is challenging in real experiments. It leads to difficulty in comparing the accuracy of different algorithms. Thus, numerical phantoms are used for comparing the accuracy of different reconstructions.

The two-dimensional numerical phantoms have a size of $401 \times 401$ spanning the imaging region of $20.1 mm \times 20.1 mm$ . The schematic diagram showing the imaging domain along with detectors placement is given in Fig. 1. The experimental data were generated on a high dimensional grid ( $401 \times 401$ ) and the reconstruction was performed on a lower dimensional grid ( $201 \times 201$ ). The data generated using the high dimensional grid was added with a white Gaussian noise, resulting in different SNR levels ranging from 20 to 60 dB. An open source tool box $k$ -WAVE in MATLAB was used for generating the data.³¹ The computational imaging grid has a size of $501 \times 501 pixels$ ( $0.1 mm / pixel$ ) and the detectors were placed on a circle of radius 22 mm (Fig. 1). The sampling frequency for the data collection is 20 MHz with number of time samples for each detector being 500. The hundred detectors were considered as point detectors having center frequency of 2.25 MHz and 70% bandwidth. Thus, the system matrix ( $A$ ) has a dimension of $50,000 \times 40,401$ . The speed of sound was assumed to be $1500 m / s$ and the medium was assumed to be homogeneous with no absorption or dispersion of sound. In all numerical experiments conducted here, 100 cycles of second order extrapolation was utilized.

Fig. 1

Schematic diagram showing the PA data acquisition geometry with hundred acoustic detectors (shown by dots) around the imaging domain. The computational imaging grid size is $50 mm \times 50 mm$ .

Initially, two numerical phantoms with a maximum initial pressure of 1 kPa [blood vessel and Derenzo given in Figs. 3(a) and 4(a), respectively] for targets were considered for proving the effectiveness of the proposed method. The numerical blood vessel phantom [Fig. 3(a)] consisted of thick and thin blood vessels mimicking the blood vessel structure. The modified Derenzo phantom [Fig. 4(a)] consisted of circular objects with varying diameter grouped according to size.

Fig. 2

(a) Schematic of the experimental setup used for PA data acquisition with number of detector positions being 100 for the horse hair and tube phantom described in this work. RD, rotating disc; P1, P2, P3, P4, uncoated right-angled prisms; L1, planoconcave lens; R/A/F, receiver, amplifier, and filter for the PA signal; DAQ, data acquisition card; UST, ultrasound transducer; and SM, stepper motor. (b) Photograph of triangular-shaped horse hair phantom. (c) Photograph of circular-shaped tube phantom filled with black Indian ink.

Fig. 3

(a) The numerical blood vessel (BV) phantom that was considered in this work. The reconstructed initial pressure rise using hundred detectors data (schematic of data collection geometry is shown in Fig. 1) with SNR of 40 dB using (b) RSD, (c) MPERSD, (d) RRERSD, (e) LBP, (f) TV, (g) MPETV, and (h) RRETV methods. The corresponding figures of merit, PC and CNR, are given in Fig. 5. The computational time corresponding to these methods is presented in fourth column of Tables 1 and 2.

Table 1

Computational time (in s) for the results presented in Figs. 5(a) and 5(b) corresponding to RSD, MPERSD, and RRERSD methods.

Phantom	Method	SNR=20 dB	SNR=40 dB (Figs. 3 and 4)	SNR=60 dB
BV [Figs. 5(a) and 5(b)]	RSD	3.07	89.61	1334.22
	MPERSD	2.47	86.16	279.77
	RRERSD	2.54	73.02	566.48
Derenzo [Figs. 5(a) and 5(b)]	RSD	9.32	115.76	1662.89
	MPERSD	6.08	87.61	558.53
	RRERSD	6.26	104.66	426.65

Table 2

Computational time (in s) for the results presented in Figs. 5(c) and 5(d) corresponding to TV, MPETV, and RRETV methods.

Phantom	Method	SNR=20 dB	SNR=40 dB (Figs. 3 and 4)	SNR=60 dB
BV [Figs. 5(c) and 5(d)]	TV	6.11	5.92	39.11
	MPETV	7.21	6.99	13.04
	RRETV	6.99	6.97	15.79
Derenzo [Figs. 5(c) and 5(d)]	TV	5.70	5.47	41.17
	MPETV	6.93	6.92	12.78
	RRETV	6.88	9.19	15.19

The schematic of the experimental setup for the collection of the PA data is shown in Fig. 2(a) [similar to Fig. 1(e) of Ref. 32]. A Q-switched Nd:YAG laser was used for delivering 532-nm wavelength laser of 5-ns duration at 10-Hz repetition rate. Four right-angle uncoated prisms (PS911, Thorlabs) and one uncoated planoconcave lens L1 (LC1715, Thorlabs) were used to deliver the laser pulses to the sample. The laser energy density on the phantom was $\sim 9 mJ / {cm}^{2}$ ( $< 20 mJ / {cm}^{2}$ : ANSI safety limit³³). A triangular shaped horse hair phantom was utilized for imaging. The side-length and diameter of hair are $\sim 10$ and 0.15 mm, respectively. The hair phantom was glued to the pipette tips adhered on acrylic slab.³⁴ The PA data were collected continuously around the hair phantom in full 360 deg using a 2.25-MHz flat ultrasound transducer (Olympus NDT, V306-SU) with 13 mm diameter active area and 70% nominal bandwidth. Another experimental phantom was also used to evaluate the effectiveness of the proposed method. It was circular in shape and made using low density polyethylene tubes (5-mm inner diameter), which were filled with black Indian ink. The tubes were placed at 0 and 15 mm from the scanning center and affixed at the bottom of the acrylic slab. The acquired PA signals were first amplified and filtered using a pulse amplifier (Olympus-NDT, 5072PR), and then recorded using a data acquisition (DAQ) card (GaGe, CompuScope 4227) using a single channel with a sampling frequency of 25 MHz inside a desktop (Intel Xeon 3.7 GHz 64-bit processor, 16 GB RAM, running windows 10 operating system). Synchronization of data acquisition with laser illumination was achieved using a sync signal from laser. The reconstructed PA imaging region has a size of $40 mm \times 40 mm$ containing $200 \times 200 pixels$ .

Fig. 4

(a) The numerical Derenzo phantom that was considered in this work. The reconstructed initial pressure rises using hundred detectors data (schematic of data collection geometry is shown in Fig. 1) with SNR of 40 dB using (b) RSD, (c) MPERSD, (d) RRERSD, (e) LBP, (f) TV, (g) MPETV, and (h) RRETV methods. The corresponding figures of merit, PC and CNR, are given in Fig. 5. The computational time corresponding to these methods are presented in fourth column of Tables 1 and 2.

The data collected have 2400 A-lines averaged over six times resulting in 400 detected signals collected by the ultrasound transducers acquiring the data continuously around the hair phantom in full 360 deg for an acquisition time of 240 s with a rotational speed of 1.5 deg/s. This averaging of A-lines reduced the energy fluctuations during multiple firings of the laser. The ultrasound transducers and the phantom are immersed in water to enable ultrasound coupling. The PA data were acquired with sampling frequency of 25 MHz (1024 samples) and subsampled at half the rate at 12.5 MHz to result in 512 time samples. The system matrix has a dimension of $51,200 \times 40,000$ (51,200: number of detectors is 100 with each collecting 512 time samples, 40,000: dimensions of the imaging domain being $200 \times 200$ ). Determining the actual initial pressure rise (target values) in these experiments is not plausible, so only a comparison between the proposed and standard methods was presented here with SNR being the figure of merit for quantitative comparison. Note that a Linux workstation with 32 cores of Intel Xeon processor having a speed of 3.1 GHz with 128 GB RAM was used for all computations performed in this work.

5. Results and Discussion

The reconstructed initial pressure distribution for the numerical blood vessel phantom [Fig. 3(a)] data having SNR of 40 dB using RSD method is shown in Fig. 3(b). The reconstructions using MPERSD and RRERSD methods are shown in Figs. 3(c) and 3(d), respectively. The LBP result is shown in Fig. 3(e). Even though the required computation time for obtaining LBP results is only 0.5 s, it can be observed that the LBP performance is poor compared to other methods in terms of both image quality (streak artifacts are more observable in LBP image) as well as reconstructed dynamic range of initial pressure (the maximum value is only going up to 0.1 against expected 1). The figures of merit, PC and CNR, for these results along with the results (not shown) for SNR of data being 20 and 60 dB are shown in Figs. 5(a) and 5(b). The corresponding computational time (in seconds) is given in Table 1. The reconstruction times presented here are for smallest order of extrapolation (quadratic). From Table 1, it is evident that the time taken by proposed extrapolation methods (MPERSD and RRERSD) was always less as compared to the standard (RSD) method. From Figs. 5(a) and 5(b), it can be seen that the PC and CNR values are nearly same for both standard and proposed methods for SNR values of 40 and 60 dB. For 20 dB, the standard method performs slightly better in terms of PC and CNR. Note that the proposed extrapolation methods generally perform well when the initial iterates are robust to noise and in the case of 20-dB initial iterates may not be well regularized solutions, thus causing the convergence to a local minima. In terms of speedup (top rows corresponding to blood vessel phantom results of Table 1), for SNR of 60 dB, it is 4.7 for MPERSD and 2.3 for RRERSD. Clearly, for low noise cases (SNR of 60 dB), there is a significant reduction in the computation time. For SNR of 40 dB, the speedup is 1.2 for RRERSD, which shows an improvement over the standard RSD method. For the case of SNR of 20 dB, the overall computation time was reduced using extrapolation methods, but these are not significant.

Fig. 5

Figures of merit (a) PC coefficient and (b) CNR for RSD, MPERSD, and RRERSD methods for the blood vessel (BV) and Derenzo phantoms for SNR in the data being 20, 40, and 60 dB (given in the parenthesis of each caption of bar graph). The computational times corresponding to these results are presented in Table 1. Similar results of (c) PC and (d) CNR for TV, MPETV, and RRETV. The computational times corresponding to these results are presented in Table 2. The reconstructed initial pressure distributions corresponding to SNR of 40 dB was shown in Figs. 3 and 4.

The reconstructed initial pressure distribution for the blood vessel phantom [Fig. 3(a)] using TV is shown in Fig. 3(f). The reconstructions using MPETV and RRETV are shown in Figs. 3(g) and 3(h). The figures of merit for these are shown in Figs. 5(c) and 5(d), along with results corresponding to varying SNR values of 20 and 60 dB. The corresponding computational time (in seconds) is reported in Table 2. From Table 2, it is evident that the time taken by MPETV and RRETV is less as compared to the standard TV method for the 60-dB SNR case. From Figs. 5(c) and 5(d), it is evident that the PC and CNR values are similar for the methods discussed in this work. For the case of SNR being 60 dB, the speedup obtained is 2.9 for MPETV and 2.4 for RRETV.

Similarly for the Derenzo phantom, the reconstructed initial pressure distribution is shown in Fig. 4(b) for the RSD method. The reconstructions for MPERSD and RRERSD are shown in Figs. 4(c) and 4(d) correspondingly. The computational times reported in Table 1 show that the speedup obtained for 60-dB SNR case is 2.9 for MPERSD and 3.89 for RRERSD. The same trend of LBP result having poor performance as observed in blood vessel phantom case was also observed here [result is presented as Fig. 4(e)]. So, for rest of experimental phantom cases, the LBP method was not utilized in the comparison. The reconstruction result pertaining to TV for Derenzo phantom case is given in Fig. 4(f). The reconstructions for MPETV and RRETV are shown in Figs. 4(g) and 4(h). Similarly, from Table 2, the speedup obtained for 60-dB SNR is 3.2 for MPETV and 2.7 for RRETV. There could not be any speed-up achieved by the proposed methods (MPETV and RRETV) for SNR values of 20 and 40 dB in these numerical phantom cases. TV has the property of preserving the features of an image. It finds the solution, which has the total minimal variation but is not biased toward a smooth or sharp solution. Thus, the disk boundaries are smooth in case of TV regularization. The RSD is equivalent of solving least-squares, which promotes only smooth solutions, does not have any properties of minimizing the variation and hence, the disk boundaries have stronger pixel values.

The reconstruction results for the experimental horse hair phantom are shown in Figs. 6(a) and 6(d) for the standard methods, RSD and TV correspondingly. The reconstruction results pertaining to MPERSD and RRERSD methods are shown in Figs. 6(b) and 6(c) correspondingly. The reconstructions for the MPETV and RRETV methods are shown in Figs. 6(e) and 6(f) correspondingly. From the results, it is evident that the TV-based methods provide better performance compared to RSD methods. Table 3 provides computational time needed for computing these results, and it is evident that the extrapolation methods that are proposed here provide significant speedup, specifically, it is 1.2 for MPERSD and 1.6 for RRERSD, whereas it is 3 for MPETV and 3.1 for RRETV. The extrapolated methods are clearly efficient as compared to the standard methods in terms of computation and also provide better SNR in terms of reconstruction results.

Table 3

Computational time (in s) required for obtaining the results presented in Figs. 6 and 7 corresponding to experimental horse hair and the tube phantom.

Phantom	RSD	MPERSD	RRERSD	TV	MPETV	RRETV
Horse hair	150.88	127.19	94.93	70.47	23.16	22.70
Tube	287.45	192.55	153.02	54.82	29.46	30.12

Fig. 6

The reconstructed initial pressure distribution with experimental horse hair phantom data using (a) RSD, (b) MPERSD, (c) RRERSD, (d) TV, (e) MPETV, and (f) RRETV methods. The figure of merit, SNR (in dB), corresponding to each reconstructed image is given below. The computational time required for obtaining these results is presented as second row in Table 3.

The axial resolution that was calculated from the reconstructed horse hair phantom image (estimated FWHM in Fig. 6) was found to be $301 μ m$ .³⁵ This resolution matches closely with the theoretical expected spatial resolution ( $\sim 333 μ m$ ) for a 2.25-MHz transducer. As stated earlier, the diameter of the hair phantom was $150 μ m$ . To resolve such object, one can use higher frequency transducer of 5 MHz³⁶ and as the aim of this part of the work is to show the real-time utility of the proposed vector extrapolation methods, this was not attempted.

The reconstruction results for the experimental tube phantom are shown in Figs. 7(a) and 7(d) for the standard methods, RSD and TV correspondingly. The reconstruction results pertaining to MPERSD and RRERSD methods are shown in Figs. 7(b) and 7(c) correspondingly. The reconstructions for the MPETV and RRETV methods are shown in Figs. 7(e) and 7(f) correspondingly. Here, as the expected objects are more circular in nature, the improvement in terms of SNR is not significant with TV-based methods showing a marginal improvement compared to its counterparts. Table 3 provides computational time needed for obtaining these results, and it is evident that the extrapolation methods that are proposed here provide significant speedup, specifically, it is 1.5 for MPERSD and 1.9 for RRERSD, whereas it is 1.9 for MPETV and 1.8 for RRETV.

Fig. 7

The reconstructed initial pressure distribution with experimental tube phantom data using (a) RSD, (b) MPERSD, (c) RRERSD, (d) TV, (e) MPETV, and (f) RRETV methods. The figure of merit, SNR (in dB), corresponding to each reconstructed image is given below. The computational time required for obtaining these results is presented as last row in Table 3.

The postprocessing of reconstructed images using Hilbert transform can improve the visualization of initial pressure distributions (currently, they are bipolar, although the expected is unipolar).³⁷^–³⁹ As the main aim of this work is to show the speed up obtained by vector extrapolation schemes, this postprocessing step was not included to show the raw performance of proposed methods. The performance of the reconstruction algorithm can also be characterized by the dynamic range of the reconstructed initial pressure [for example, compare Figs. 4(e) and 4(f)], normalizing these plots either using Hilbert or other transformation methods will make the representations inaccurate.

The work presented here is geared toward accelerating the image reconstruction procedure in case of limited data. The definition of limited data here is relative to the state-of-the-art setups. To be precise, the current clinical PA tomography scanners with advances in instrumentation acquire large amounts of data with number of detectors being as high as 512 and number of time samples for each detector being 2048.⁴⁰ Note that in this work, the maximum number of detectors that were utilized is 100 with number of time samples being 512, thus making the data available here relatively limited.

It is important to note that the vector extrapolation methods presented in this work are known to be providing global optima as they estimate the limit (or antilimit) of the initial solutions of an iterative method. Thus, the performance of these methods largely depends on the initial solutions and often regularized solutions (like the ones obtained using TV) can improve the performance of these methods (which is also observed here). Ideally, in noise-free cases, both extrapolation-based iterative methods and standard iterative methods should converge to the same solution as long as the stopping criterion is same for both methods (which is the case here). Unfortunately, noise-free cases are far away from real-time scenarios, so, this exercise was not taken up and interested readers can refer to Refs. 26–29 for detailed mathematical framework.

The important step in the PA tomography is image reconstruction and developing computationally efficient methods that can be deployable in real-time is desirable to translate PA imaging into a clinical imaging modality. Especially, in cases of limited data, most reconstruction algorithms present in the literature even though provide better quantitatively accurate results, often they are prohibitively expensive in terms of computation. In this work, we have deployed extrapolation techniques to accelerate two standard reconstruction methods and shown using example cases that these methods can be very effective in reducing the total computation time required especially in experimental scenarios. It is important to note that the extrapolation techniques that are presented in this work are more of universal in nature and can be employed in any iterative reconstruction method. The results presented here are limited in nature as the focus is more toward showing the utility of these extrapolation methods, one could expect the similar speed-up in other cases as well. The corresponding algorithms along with necessary code are given as an open-source for enthusiastic users.⁴¹

6. Conclusions

The vector extrapolation methods have been shown to improve any iterative procedure, where there will be a reduction in the residual error as the iterations progress (convergence is guaranteed). Application of these extrapolation methods in the frame work of two iterative methods that are commonly used in PA imaging was presented here. Specifically, the initial iterates that were calculated using RSD or TV were given as input to these vector extrapolation methods and was shown to provide significant speed-up (as high as 4.7 times) especially in numerical cases. The extrapolation techniques are important in achieving the limit of the function without iterating through a large number of steps as frequently performed by many iterative methods. It is shown using both numerical and experimental phantom cases, the extrapolation methods (two of such variants were presented in this work), have a potential in terms of speeding up the reconstruction algorithm along with improving the reconstructed image quality.

Disclosures

No conflicts of interest, financial or otherwise, are declared by the authors.

References

1.

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

2.

M. Pramanik et al., “Design and evaluation of a novel breast cancer detection system combining both thermoacoustic (TA) and photoacoustic (PA) tomography,” Med. Phys., 35 (6), 2218 –2223 (2008). http://dx.doi.org/10.1118/1.2911157 MPHYA6 0094-2405 Google Scholar

3.

L. V. Wang, “Tutorial on photoacoustic microscopy and computed tomography,” IEEE J. Sel. Top. Quantum Electron., 14 (1), 171 –179 (2008). http://dx.doi.org/10.1109/JSTQE.2007.913398 IJSQEN 1077-260X Google Scholar

4.

P. K. Upputuri and M. Pramanik, “Recent advances toward preclinical and clinical translation of photoacoustic tomography: a review,” J. Biomed. Opt., 22 (4), 041006 (2017). http://dx.doi.org/10.1117/1.JBO.22.4.041006 JBOPFO 1083-3668 Google Scholar

5.

L. Li et al., “Single-impulse panoramic photoacoustic computed tomography of small-animal whole-body dynamics at high spatiotemporal resolution,” Nat. Biomed. Eng., 1 (5), 0071 (2017). http://dx.doi.org/10.1038/s41551-017-0071 Google Scholar

6.

Y. Zhou, J. Yao and L. V. Wang, “Tutorial on photoacoustic tomography,” J. Biomed. Opt., 21 (6), 061007 (2016). http://dx.doi.org/10.1117/1.JBO.21.6.061007 JBOPFO 1083-3668 Google Scholar

7.

G. Paltauf et al., “Iterative reconstruction algorithm for optoacoustic imaging,” J. Acoust. Soc. Am., 112 (4), 1536 –1544 (2002). http://dx.doi.org/10.1121/1.1501898 Google Scholar

8.

K. Wang et al., “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography,” Phys. Med. Biol., 57 (17), 5399 –5423 (2012). http://dx.doi.org/10.1088/0031-9155/57/17/5399 PHMBA7 0031-9155 Google Scholar

9.

X. L. Dean-Ben, V. Ntziachristos and D. Razansky, “Acceleration of optoacoustic model-based reconstruction using angular image discretization,” IEEE Trans. Med. Imaging, 31 (5), 1154 –1162 (2012). http://dx.doi.org/10.1109/TMI.2012.2187460 ITMID4 0278-0062 Google Scholar

10.

X. L. Dean-Ben et al., “Accurate model-based reconstruction algorithm for three-dimensional optoacoustic tomography,” IEEE Trans. Med. Imaging, 31 (10), 1922 –1928 (2012). http://dx.doi.org/10.1109/TMI.2012.2208471 ITMID4 0278-0062 Google Scholar

11.

A. Buehler et al., “Model-based optoacoustic inversions with incomplete projection data,” Med. Phys., 38 (3), 1694 –1704 (2011). http://dx.doi.org/10.1118/1.3556916 MPHYA6 0094-2405 Google Scholar

12.

Y. Xu et al., “Reconstructions in limited-view thermoacoustic tomography,” Med. Phys., 31 (4), 724 –733 (2004). http://dx.doi.org/10.1118/1.1644531 MPHYA6 0094-2405 Google Scholar

13.

S. Arridge et al., “Accelerated high-resolution photoacoustic tomography via compressed sensing,” Phys. Med. Biol., 61 (24), 8908 –8940 (2016). http://dx.doi.org/10.1088/1361-6560/61/24/8908 PHMBA7 0031-9155 Google Scholar

14.

S. R. Arridge et al., “On the adjoint operator in photoacoustic tomography,” Inverse Prob., 32 (11), 115012 (2016). http://dx.doi.org/10.1088/0266-5611/32/11/115012 INPEEY 0266-5611 Google Scholar

15.

S. Gutta et al., “Modeling errors compensation with total least squares for limited data photoacoustic tomography,” IEEE J. Sel. Top. Quantum Electron., (2017). http://dx.doi.org/10.1109/JSTQE.2017.2772886 Google Scholar

16.

M. S. Zhdanov, Geophysical Inverse Theory and Regularization Problems, 36 Elsevier, Amsterdam, Netherlands (2002). Google Scholar

17.

M. Tao, J. Yang and B. He, “Alternating direction algorithms for total variation deconvolution in image reconstruction,” Nanjing, Jiangsu, China (2009). Google Scholar

18.

Y. Wang et al., “A new alternating minimization algorithm for total variation image reconstruction,” SIAM J. Imaging. Sci., 1 (3), 248 –272 (2008). http://dx.doi.org/10.1137/080724265 Google Scholar

19.

C. B. Shaw et al., “Least squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography,” J. Biomed. Opt., 18 (8), 080501 (2013). http://dx.doi.org/10.1117/1.JBO.18.8.080501 JBOPFO 1083-3668 Google Scholar

20.

J. Prakash et al., “Basis pursuit deconvolution for improving model-based reconstructed images in photoacoustic tomography,” Biomed. Opt. Express, 5 (5), 1363 –1377 (2014). http://dx.doi.org/10.1364/BOE.5.001363 BOEICL 2156-7085 Google Scholar

21.

G. L. Zeng, Medical Image Reconstruction: A Conceptual Tutorial, Springer, New York (2010). Google Scholar

22.

A. Rosenthal, V. Ntziachristos and D. Razansky, “Acoustic inversion in optoacoustic tomography: a review,” Curr. Med. Imaging Rev., 9 (4), 318 –336 (2013). http://dx.doi.org/10.2174/15734056113096660006 Google Scholar

23.

A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vision, 20 (1), 89 –97 (2004). https://doi.org/10.1023/B:JMIV.0000011325.36760.1e Google Scholar

24.

M. V. Afonso, J. M. Bioucas-Dias and M. A. Figueiredo, “Fast image recovery using variable splitting and constrained optimization,” IEEE Trans. Image Process., 19 (9), 2345 –2356 (2010). http://dx.doi.org/10.1109/TIP.2010.2047910 IIPRE4 1057-7149 Google Scholar

25.

J. Eckstein and D. P. Bertsekas, “On the Douglas–Rachford splitting method and the proximal point algorithm for maximal monotone operators,” Math. Program., 55 (1), 293 –318 (1992). http://dx.doi.org/10.1007/BF01581204 MHPGA4 1436-4646 Google Scholar

26.

M. Mešina, “Convergence acceleration for the iterative solution of the equations

x = ax + f

,” Comput. Methods Appl. Mech. Eng., 10 (2), 165 –173 (1977). http://dx.doi.org/10.1016/0045-7825(77)90004-4 CMMECC 0045-7825 Google Scholar

27.

S. Cabay and L. Jackson, “A polynomial extrapolation method for finding limits and antilimits of vector sequences,” SIAM J. Numer. Anal., 13 (5), 734 –752 (1976). http://dx.doi.org/10.1137/0713060 SJNAEQ 0036-1429 Google Scholar

28.

D. A. Smith, W. F. Ford and A. Sidi, “Extrapolation methods for vector sequences,” SIAM Rev., 29 (2), 199 –233 (1987). http://dx.doi.org/10.1137/1029042 SIREAD 0036-1445 Google Scholar

29.

A. Sidi, “Minimal polynomial and reduced rank extrapolation methods are related,” Adv. Comput. Math., 43 (1), 151 –170 (2017). http://dx.doi.org/10.1007/s10444-016-9481-0 ACMHEX 1019-7168 Google Scholar

30.

J. Prakash et al., “Sparse recovery methods hold promise for diffuse optical tomographic image reconstruction,” IEEE J. Sel. Top. Quantum Electron., 20 (2), 74 –82 (2014). http://dx.doi.org/10.1109/JSTQE.2013.2278218 IJSQEN 1077-260X Google Scholar

31.

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

32.

S. K. Kalva and M. Pramanik, “Experimental validation of tangential resolution improvement in photoacoustic tomography using modified delay-and-sum reconstruction algorithm,” J. Biomed. Opt., 21 (8), 086011 (2016). http://dx.doi.org/10.1117/1.JBO.21.8.086011 JBOPFO 1083-3668 Google Scholar

33.

American Institute, “American national standard for safe use of lasers (ansi z136. 1-2000),” (2000). Google Scholar

34.

P. K. Upputuri and M. Pramanik, “Pulsed laser diode based optoacoustic imaging of biological tissues,” Biomed. Phys. Eng. Express, 1 (4), 045010 (2015). http://dx.doi.org/10.1088/2057-1976/1/4/045010 Google Scholar

35.

S. K. Kalva and M. Pramanik, “Use of acoustic reflector to make a compact photoacoustic tomography system,” J. Biomed. Opt., 22 (2), 026009 (2017). http://dx.doi.org/10.1117/1.JBO.22.2.026009 JBOPFO 1083-3668 Google Scholar

36.

P. K. Upputuri and M. Pramanik, “Performance characterization of low-cost, high-speed, portable pulsed laser diode photoacoustic tomography (PLD-PAT) system,” Biomed. Opt. Express, 6 (10), 4118 –4129 (2015). http://dx.doi.org/10.1364/BOE.6.004118 BOEICL 2156-7085 Google Scholar

37.

G. Li et al., “Multiview Hilbert transformation for full-view photoacoustic computed tomography using a linear array,” J. Biomed. Opt., 20 (6), 066010 (2015). http://dx.doi.org/10.1117/1.JBO.20.6.066010 JBOPFO 1083-3668 Google Scholar

38.

L. Li et al., “Multiview Hilbert transformation in full-ring transducer array-based photoacoustic computed tomography,” J. Biomed. Opt., 22 (7), 076017 (2017). http://dx.doi.org/10.1117/1.JBO.22.7.076017 JBOPFO 1083-3668 Google Scholar

39.

Q. Sheng et al., “A constrained variable projection reconstruction method for photoacoustic computed tomography without accurate knowledge of transducer responses,” IEEE Trans. Med. Imaging, 34 (12), 2443 –2458 (2015). http://dx.doi.org/10.1109/TMI.2015.2437356 ITMID4 0278-0062 Google Scholar

40.

R. A. Kruger et al., “Dedicated 3D photoacoustic breast imaging,” Med. Phys., 40 (11), 113301 (2013). http://dx.doi.org/10.1118/1.4824317 MPHYA6 0094-2405 Google Scholar

41.

N. Awasthi et al., “Vector extrapolation based methods for accelerating iterative reconstruction methods in limited-data photoacoustic tomography,” (2017) https://sites.google.com/site/sercmig/home/patextrapolation September ). 2017). Google Scholar

Biography

Navchetan Awasthi received his MTech in computational science from Indian Institute of Science (IISc), Bangalore, in 2016. He also received his BTech in electronics and communication engineering from the National Institute of Technology (NIT), Jalandhar, in 2011. He is currently a PhD student in the Department of Computational and Data Sciences at Indian Institute of Science, Bangalore. His research interests include inverse problems in biomedical optics and medical image reconstruction.

Sandeep Kumar Kalva received his master’s degree in biomedical engineering from the Indian Institute of Technology (IIT), Hyderabad, India, in 2015. He is currently a PhD student at the School of Chemical and Biomedical Engineering, Nanyang Technological University, Singapore. His research area is on functional photoacoustic imaging for various clinical applications.

Manojit Pramanik received his PhD in biomedical engineering from Washington University in St. Louis, Missouri. Currently, he is an assistant professor at the School of Chemical and Biomedical Engineering, Nanyang Technological University, Singapore. His research interests include the development of photoacoustic/thermoacoustic imaging systems, image reconstruction methods, clinical application areas, such as breast cancer imaging, molecular imaging, contrast agent development, and Monte Carlo simulation of light propagation in biological tissue.

Phaneendra K. Yalavarthy received his MSc degree in engineering from Indian Institute of Science, Bangalore, and his PhD degree in biomedical computation from Dartmouth College, Hanover, New Hampshire, in 2007. He is an associate professor with the Department of Computational and Data Sciences at Indian Institute of Science, Bangalore. His research interests include medical image computing, medical image analysis, and biomedical optics. He is a senior member of IEEE and serves as an associate editor of IEEE Transactions on Medical Imaging.

Citation Download Citation

Navchetan Awasthi, Sandeep Kumar Kalva, Manojit Pramanik, and Phaneendra K. Yalavarthy "Vector extrapolation methods for accelerating iterative reconstruction methods in limited-data photoacoustic tomography," Journal of Biomedical Optics 23(7), 071204 (5 February 2018). https://doi.org/10.1117/1.JBO.23.7.071204

Received: 29 September 2017; Accepted: 8 January 2018; Published: 5 February 2018

Access the abstract

JOURNAL ARTICLE
11 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 13 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Signal to noise ratio

Reconstruction algorithms

Data acquisition

Sensors

Image restoration

Blood vessels

Iterative methods

1.

Introduction

2.

Photoacoustic Image Reconstruction

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

Eq. (5)

2.1.

Regularized Steepest Descent Method

Algorithm 1

2.2.

Total Variation Regularization Method

Eq. (6)

Eq. (7)

Algorithm 2

Eq. (8)

2.3.

Proposed Vector Extrapolation Methods

2.3.1.

Minimal polynomial extrapolation method

Eq. (9)

Eq. (10)

Eq. (11)

Eq. (12)

Eq. (13)

Eq. (14)

Eq. (15)

Eq. (16)

Eq. (17)

Eq. (18)

Eq. (19)

Eq. (20)

Eq. (21)

Eq. (22)

Eq. (23)

Eq. (24)

Eq. (25)

Eq. (26)

Eq. (27)

Eq. (28)

Eq. (29)

2.3.2.

Reduced rank extrapolation method

Eq. (30)

Algorithm 3

3.

Figures of Merit

3.1.

Pearson Correlation

Eq. (31)

3.2.

Contrast to Noise Ratio

Eq. (32)

3.3.

Signal-to-Noise Ratio

Eq. (33)

4.

Numerical and Experimental Phantom Studies

Fig. 1

Fig. 2

Fig. 3

Table 1

Table 2

Fig. 4

5.

Results and Discussion

Fig. 5

Table 3

Fig. 6

Fig. 7

6.

Conclusions

Disclosures

References

Biography

Show All Keywords

Keywords/Phrases

Search In: