Computationally efficient error estimate for evaluation of regularization in photoacoustic tomography

Manish Bhatt; Atithi Acharya; Phaneendra K. Yalavarthy

doi:10.1117/1.JBO.21.10.106002

20 October 2016 Computationally efficient error estimate for evaluation of regularization in photoacoustic tomography

Manish Bhatt, Atithi Acharya, Phaneendra K. Yalavarthy

Author Affiliations +

Journal of Biomedical Optics, Vol. 21, Issue 10, 106002 (October 2016). https://doi.org/10.1117/1.JBO.21.10.106002

Abstract

The model-based image reconstruction techniques for photoacoustic (PA) tomography require an explicit regularization. An error estimate (η²) minimization-based approach was proposed and developed for the determination of a regularization parameter for PA imaging. The regularization was used within Lanczos bidiagonalization framework, which provides the advantage of dimensionality reduction for a large system of equations. It was shown that the proposed method is computationally faster than the state-of-the-art techniques and provides similar performance in terms of quantitative accuracy in reconstructed images. It was also shown that the error estimate (η²) can also be utilized in determining a suitable regularization parameter for other popular techniques such as Tikhonov, exponential, and nonsmooth (ℓ₁ and total variation norm based) regularization methods.

1. Introduction

Photoacoustic tomography (PAT) is a rapidly growing noninvasive, in vivo hybrid biomedical imaging technique combining both optics and acoustics. It has applications in breast cancer imaging, brain imaging, temperature monitoring, sentinel lymph node imaging, and so on.¹^–⁶ It holds a great potential for anatomical, functional, and molecular imaging in preclinical and clinical medicine.⁶^–⁹ In PAT, to-be-imaged biological tissue is illuminated with a nanosecond laser pulse, which generates internal acoustic wavefields through the thermoacoustic effect.⁶ The initial amplitude of induced acoustic wavefield (pressure rise) is proportional to the absorbed optical energy density in the tissue. These acoustic wavefields propagate out of the tissue and are measured by the use of an array of ultrasonic transducers. From the measured acoustic data, an image reconstruction algorithm, analytical, or model-based in nature,¹⁰^–¹⁴ is employed to obtain an estimate of the image that depicts the initial pressure rise or the spatially variant optical energy density distribution within the imaged tissue.

Many analytical image reconstruction algorithms have been proposed in the literature.¹⁵ These algorithms, in the form of filtered back projection, are based on spherical Radon transform. They require a large amount of data to accurately estimate the initial pressure distribution. Large data acquisition needs expensive ultrasonic transducer arrays along with increased data collection time when a mechanical scanning is employed. Moreover, the analytical algorithms do not provide much required quantitative accuracy in limited data cases.¹⁶

Model-based image reconstruction algorithms¹⁵^,¹⁶ can improve image quality in these limited data cases and have a larger utility for practical scenarios. Thus, development and investigation of model-based algorithms for PAT image reconstruction is an important research topic of current interest. In limited-data case, the regularized least-squares QR (LSQR) based image reconstruction scheme was earlier deployed by our group¹⁷ and proven to be highly efficient and quantitatively more accurate in providing regularized solution.¹⁷^,¹⁸ Traditional regularization schemes, such as Tikhonov regularization, become relatively expensive in terms of computation with a requirement of having an explicit appropriate regularization parameter.

Photoacoustic (PA) image reconstruction problem deals with estimation of initial condition of an observed effect. This is a discrete ill-posed problem, where the system matrix is ill conditioned in nature. The computation of numerical solution is challenging due to the possibility of severe error propagation (caused by measurement inaccuracies, transmission errors of the data, and discretization errors). The regularized solutions, such as Tikhonov regularization, reduce the sensitivity of the computed solution to errors in the data. These regularization methods require a suitable value of regularization parameter to be chosen. In this study, we investigate a class of error estimates $(η_{2})$ for the approximate solution and propose to apply this error estimate to compute a suitable value of the regularization parameter. This study focuses on the widely popular Tikhonov regularization framework, which uses Lanczos bidiagonalization to provide dimensionality reduction to the system matrix in case of limited-view datasets.¹⁷^,¹⁸ The applications of the proposed error estimates to various direct and iterative regularization methods were also discussed, and the effectiveness of this analysis is proven by extensive numerical experimentation. It is proven that the class of error estimate that is proposed here is very effective and computationally efficient for evaluation of regularization in PA imaging.

2. Forward Problem in Photoacoustic Tomography: System Matrix Approach

A system matrix-based approach was utilized,¹⁸ where the PA system is represented by a Toeplitz matrix of a time-varying causal system. This approach was extensively discussed in Refs. 17 and 18, herein it is only briefly reviewed.

The image representing the imaging domain, having a dimension of $n \times n$ , is vectorized by stacking all columns in lexicographic manner, converting it into a $n^{2} \times 1$ size vector. This vector was represented as $x$ , the unknown to be estimated in the PA inverse problem. The time varying PA data collected at the ultrasonic transducers was also stacked into a $m \times 1$ dimensional vector, represented as the measurement vector $b$ . The PA system matrix $(A)$ has a dimension of $m \times n^{2}$ , each column representing impulse response (IR) of a corresponding pixel in the image, as explained in Refs. 17 and 18. The k-Wave toolbox,¹⁹ which simulates the acoustic wave propagation in two dimension, was utilized to construct the PA system matrix. This simulation is computationally expensive to perform for all pixels. To reduce the computational time for building the system matrix, IR is recorded for only one pixel (corner pixel), and IRs for the rest of the pixels were determined using the shifting and attenuation properties of PA signal. The IRs for all pixels are computed and filled in column-by-column in the system matrix form.

It was assumed that the medium has homogeneous ultrasound properties. The sound of speed was modeled to be $1500 m / s$ (approximate value for typical soft tissues). The ultrasound detectors were considered to be point detectors with a frequency response of 2.25 MHz as center frequency and 70% bandwidth. These detectors were placed equidistantly around a 22-mm radius circle. A perfectly matched layer was used to satisfy the boundary condition. The simulation geometry had a k-Wave computational grid of size $501 \times 501 pixel$ size, each pixel measuring 0.1 mm ( $50 mm \times 50 mm$ ). The imaging domain had $201 \times 201 pixels$ $(n^{2} = 40,401)$ . The time step for data collection was 50 ns, with a total of 500 steps. For a total of 60 detectors’ geometry, $m = 60 \times 500$ . Thus the size of the system matrix $(m \times n^{2})$ becomes $30,000 \times 40,401$ .

To summarize, the forward model of PA imaging can be expressed as

Eq. (1)

A x = b,

where

A

is the system matrix (collection of IRs of each pixel),

x

is initial pressure value (unknown) for each pixel in the imaging domain, and

b

is the vector containing the recorded measurement data from ultrasound transducers.

3. Error Estimate for a Regularized Solution

The error ( $e$ ) in the measurement vector ( $b$ ), also referred as noise, can be caused by measurement inaccuracies and discretization errors. Let $\hat{b} \in R^{m}$ (unavailable) be the noise-free vector associated with $b$ :

Eq. (2)

b = \hat{b} + e,

and let

\hat{x} \in R^{n^{2}}

be the minimal-norm least squares solution of the following problem,

Eq. (3)

\underset{x \in R^{n^{2}}}{\arg \min} {‖ A x - b ‖}_{2}^{2} .

Here,

{‖ \cdot ‖}_{2}

represents

ℓ_{2}

norm. The aim of the PA inverse problem is to estimate the initial pressure distribution,

x

, from the above equation.¹⁷^,¹⁸ In order to compute a useful approximation, the minimization problem is replaced by a nearby problem, which is less sensitive to noise in

b

, also known as Tikhonov regularization. The cost function minimization for this problem becomes

Eq. (4)

x_{λ} = \underset{x \in R^{n^{2}}}{\arg \min} {‖ A x - b ‖}_{2}^{2} + λ {‖ R x ‖}_{2}^{2},

where

λ > 0

is a regularization parameter and

R

is a regularization operator. For the conventional Tikhonov regularization

R = I

; other possible choices of

R

are discussed in Ref. 20. The minimization of cost function in Eq. (4) leads to the solution

Eq. (5)

x_{λ} = {(A^{T} A + λ I)}^{- 1} A^{T} b,

for any

λ > 0

. The regularization parameter (

λ

) determines the sensitivity of the solution

x_{λ}

to the noise/error

e

in

b

, and how close

x_{λ}

would be to,

\hat{x}

, the desired solution of Eq. (3). Thus, determining a suitable choice of

λ

is a very important part of finding solution. This work addresses this problem and proposes a method based on error estimates to find the regularization parameter.

The initial pressure rise ( $x$ ) in the imaging domain is unknown in experiments. The associated residual vector is defined to be $r := \hat{b} - A x$ . Brezinski et al.²¹ derived the family of error estimates for the norm of the difference, i.e., $d = x - \hat{x}$ , given as²²

Eq. (6)

{‖ d ‖}_{2}^{2} \approx η_{ν}^{2} := d_{o}^{ν - 1} d_{1}^{5 - 2 ν} d_{2}^{ν - 3}, ν \in R,

where

Eq. (7)

d_{o} := {‖ r ‖}_{2}^{2}, d_{1} := {‖ A^{T} r ‖}_{2}^{2}, d_{2} := {‖ {AA}^{T} r ‖}_{2}^{2} .

Proceeding from Eq. (6), the error estimate for $ν = 2$ can be defined as

Eq. (8)

η_{2} = \frac{{‖ r ‖}_{2} {‖ A^{T} r ‖}_{2}}{{‖ {AA}^{T} r ‖}_{2}} .

Since the noise-free measurement vector ( $\hat{b}$ ), and the true solution ( $x$ ) is not available, the residual vector ( $r$ ) in Eq. (7) is written as

Eq. (9)

r = r_{λ} := b - A x_{λ} .

The regularization parameter ( $λ$ ) is determined by the minimization of the expression in Eq. (8):

Eq. (10)

λ = \underset{λ \in R > 0}{\arg \min} η_{2} (λ) = \underset{λ \in R > 0}{\arg \min} \frac{{‖ r_{λ} ‖}_{2}^{2} {‖ A^{T} r_{λ} ‖}_{2}^{2}}{{‖ {AA}^{T} r_{λ} ‖}_{2}^{2}} .

Note that in this work, only an approximation to $η_{2} (λ)$ was utilized (achieved through step 2 of Algorithm 1), thus removing the requirement of computing $r_{λ}$ to find $λ$ . However, Shaw et al.¹⁷ utilized the ${‖ r_{λ} ‖}_{2}^{2}$ , which requires estimation of $x_{λ}$ , which involves matrix–matrix operations. If the singular value decomposition (SVD) of the system matrix ( $A$ ) is available, the computation of residual $r$ [Eq. (9)] and the estimates [Eq. (6)] is inexpensive. However, it may not be the case and as $A$ is large, these computations become expensive in the Tikhonov regularization framework. The cost efficiency of the system can be accomplished by reducing the matrix $A$ to a small lower bidiagonal matrix by applying few steps of Lanczos bidiagonalization.¹⁷^,¹⁸ This framework for PA image reconstruction with the application of proposed regularization method is described in next section.

Algorithm 1

Algorithm for evaluation of regularization parameter and number of Lanczos bidiagonalization steps.

Input: $A; b; [λ_{\min}, λ_{\max}]; k_{\max}$ = max number of Lanczos bidiagonalization steps

Output: Required number of Lanczos iterations: $k$ ; regularization parameter: $λ$

1 Make a grid of

q

log-equispaced values of

λ

in the interval [

λ_{\min}

,

λ_{\max}

]. Define

S_{1} := {λ_{1}, λ_{2}, \dots, λ_{q}}

and

S_{2} := Ø

.

2 For each

λ_{j}

in

S_{1}

, compute the lower bound

η_{2}^{l} (λ_{j})

, the upper bounds

η_{2}^{u} (λ_{j})

, and their average

{\overline{η}}_{2} (λ_{j})

. (Refer Theorem 3.1 of Ref. 22).

3 Find the smallest index

p

, such that

Eq. (11)

| η_{2}^{u} (λ_{j}) - η_{2}^{l} (λ_{j}) | < β \cdot {\bar{η}}_{2} (λ_{j}); j = p, p + 1, \dots, # S_{1},

where

β

is user specified tolerance (typically,

10^{- 4}

) and

# S_{1}

is the cardinality of the set

S_{1}

.

4 ifFor above $p$ ’th index, Eq. (11) = truethen

5 Move

{λ_{p}, λ_{p + 1}, \dots, λ_{# N}}

from

S_{1}

to

S_{2}

6 End

7 Increase

k := k + 1

, i.e., perform another Lanczos bidiagonalization step and repeat the algorithm steps 1 to 5 until

{\bar{η}}_{2} (λ_{p}) > {\bar{η}}_{2} (λ_{p + 1})

for some pair

{λ_{p}, λ_{p + 1}} \subset S_{2}

, this gives local minima

λ^{*}

of the function

{\overline{η}}_{2} (λ)

and

k

is the corresponding number of Lanczos steps.

8 Add points around

λ^{*}

, by bisection in log-scale, to refine the set

S_{2}

.

9 Repeat step 3. Update value of

λ^{*}

.

10 if $λ^{*} - λ_{adjacent} < tolerance (10^{- 4})$ then

11

Required λ = λ^{*}

;

12 Else

13 Repeat steps 8 to 11;

14 End

3.1.

Regularized LSQR Method

The equations of Lanczos bidiagonalization can be written as²³

Eq. (12)

G_{k + 1} (β_{0} e_{1}) = b,

Eq. (13)

{AH}_{k} = G_{k + 1} B_{k},

Eq. (14)

A^{T} G_{k + 1} = H_{k} B_{k}^{T} + α_{k + 1} h_{k + 1} e_{k + 1}^{T},

where

G_{k}

and

H_{k}

represent the left and right orthogonal Lanczos matrices,

β_{0}

represents the

l_{2}

norm of

b

,

e_{1}

as a unit vector of dimension

k \times 1

, and

B_{k}

represents the lower bidiagonal matrix having

α_{k}

in the main diagonal. The cost function can now be written as¹⁸^,²⁴

Ω = {‖ β_{0} e_{1} - B_{k} x^{(k)} ‖}_{2}^{2} + λ {‖ x^{(k)} ‖}_{2}^{2},

where

x^{(k)}

is the dimensionality reduced version of

x

, such that

k < n^{2}

. The updated equation for the solution estimate is given as

Eq. (15)

x_{est}^{(k)} = {(B_{k}^{T} B_{k} + λ I)}^{- 1} β_{0} B_{k}^{T} e_{1} .

Once

x_{est}^{(k)}

is estimated, then the initial pressure distribution can be obtained using the relation

p_{o} = x = H_{k} x_{est}^{(k)}

.

Determining the required number of Lanczos bidiagonalization steps ( $k$ ) and the regularization parameter ( $λ$ ) using the error estimate ( $η_{2}$ ) is given in Algorithm 1. We start the process with a grid of $q$ log-equispaced values of $λ$ in the interval [ $λ_{\min}$ , $λ_{\max}$ ]:

Eq. (16)

0 < λ_{1} < λ_{2} < \dots < λ_{q} .

These grid points are divided into two sets, $S_{1}$ and $S_{2}$ , as discussed in the Algorithm 1. The algorithm has two phases, the first phase computes a nearby estimate of the suitable regularization parameter $λ^{*}$ . The initialization starts with $S_{1} := {λ_{1}, λ_{2}, \dots, λ_{q}}$ and $S_{2} := Ø$ . The upper and lower bounds of the error estimate, $η_{2} (λ)$ , can be computed with the help of Gauss–Radau quadrature rules, as described in Theorem 3.1 of Ref. 22. A smallest $λ_{p}$ within the set $S_{1}$ was searched such that the absolute difference of upper and lower bound is smaller than a set threshold [See Eq. (11)]. Once this $λ_{p}$ is found, all parameters in set $S_{1}$ after $p$ ’th index will be moved to set $S_{2}$ . This indicates that ${\bar{η}}_{2} (λ_{j})$ was considered as an approximation of $η_{2} (λ_{j})$ for $j = p, p + 1, \dots, # S_{1}$ .²² Now the number of Lanczos bidiagonalization steps was increased, and the same procedure was performed until one finds an averaged error estimate [ ${\bar{η}}_{2} (λ_{p + 1})$ ], which is lower than its own value in the previous step [ ${\bar{η}}_{2} (λ_{p})$ ]. This means that a local minima $λ^{*}$ of the function $η_{2} (λ)$ was found at an internal point of $S_{2}$ (this step could be seen in Fig. 2). Once $λ^{*}$ is found, no more Lanczos bidiagonalization steps were required. This completes the first phase of the algorithm.

The second phase improves the accuracy of $λ^{*}$ . This was achieved by adding the points around $λ^{*}$ in set $S_{2}$ . The points were added by bisection in log-scale and the value of $λ^{*}$ was being updated such that Eq. (11) is valid. This step was performed until the minimum distance between adjacent grid points around $λ^{*}$ in $S_{2}$ reaches a user defined tolerance ( $10^{- 4}$ ). Thus, refined $λ = λ^{*}$ is the suitable regularization parameter.

Other popular techniques to determine regularization parameter are L-curve, generalized cross-validation (GCV) method, and LSQR method. Shaw et al.¹⁷ showed that LSQR was faster and quantitatively more accurate than L-curve and GCV methods. The LSQR method is based on the minimization of residual ${‖ b - A x ‖}_{2}^{2}$ . The detailed description of regularized LSQR method is omitted here, interested readers can refer to Algorithm 1 in Ref. 17, which minimizes the squared $ℓ_{2}$ norm of residual vector ( $r$ ). The main aim of the present study is to show that the proposed method is computationally more efficient than other popular techniques, without compromising the reconstructed image quality.

It will also be shown in the next subsections that the proposed error estimate ( $η_{2}$ ) is generic in nature and can be utilized (for example, instead of $r$ in LSQR-based method proposed in Ref. 17) for other popular model-based reconstruction techniques, with examples of traditional Tikhonov regularization, recently proposed exponential filtering,²⁴ and $ℓ_{1}$ - or TV-norm based regularization.

3.2.

Traditional Tikhonov Regularization

The Tikhonov-based $ℓ_{2}$ norm smoothness constraint is widely utilized for solving the reconstruction problem. As seen in Eq. (5), the Tikhonov regularization solution can be written as

Eq. (17)

(A^{T} A + λ I) x_{λ} = A^{T} b .

The SVD of the system matrix, i.e., its reduction to diagonal form, can be represented as

Eq. (18)

A = U Σ V^{T},

where

U

and

V

are left and right orthogonal matrices and

Σ

is a diagonal matrix containing singular (spectral) values on its diagonal. Using the SVD of

A

[Eq. (18)], the solution to Eq. (17) for this case can be rewritten as²⁴

Eq. (19)

x_{λ} = {(V Σ^{T} Σ V^{T} + λ I)}^{- 1} V Σ^{T} U^{T} b,

Eq. (20)

x_{λ} = V Σ_{F}^{†} U^{T} b,

where the filtered singular values

Σ_{F}^{†}

is given as

Eq. (21)

Σ_{F}^{†} = diag (\frac{φ_{i}}{σ_{i}}),

where

σ_{i}

is

i

’th singular value of

A

, and

φ_{i}

are the filter factors defined as²⁴

Eq. (22)

φ_{i} = \frac{σ_{i}^{2}}{σ_{i}^{2} + λ} .

Rewriting Eq. (17) in terms of residual [Eq. (9)] makes

Eq. (23)

A^{T} r_{λ} = λ x_{λ},

where

r_{λ}

is the residue expressed in Eq. (9).

x_{λ}

can be evaluated using Eq. (20). Substituting

A^{T} r_{λ}

with

λ x_{λ}

in Eq. (8), the expression becomes

Eq. (24)

η_{2} (λ) = \frac{{‖ r ‖}_{2} {‖ x_{λ} ‖}_{2}}{{‖ A x_{λ} ‖}_{2}} .

For the case of linear systems treated by Tikhonov method, above representation of the error estimate becomes numerically much more effective than Eq. (8). This estimate can be utilized in the framework of Algorithm 1 (especially in steps 8 to 14) to determine the required regularization parameter.

3.3.

Exponential Filtering or Showalter’s Method

The exponential filtering method was recently proposed for PA image reconstruction.²⁴ This method considers the exponential filtering of singular values of PA system matrix ( $A$ ). This filtering was applied in the Tikhonov framework with solution being determined using Eq. (20) with a difference that the filter factors ( $φ_{i}$ ‘s) in this case are written as

Eq. (25)

φ_{i} = 1 - \exp (- σ_{i}^{2} / λ) .

Similar to Eq. (17), the solution can be uniquely determined by the consistent linear system:

Eq. (26)

(A^{T} A + λ R^{T} R) x_{λ} = A^{T} b,

which implies

Eq. (27)

A^{T} r_{λ} = λ E x_{λ},

where

E = R^{T} R

is diagonal regularization matrix. Comparing Eqs. (20) and (26), the

E

in this case can be constructed as

Eq. (28)

E = R^{T} R = \frac{1}{λ} V Σ^{T} F Σ V^{T},

where

F

is a diagonal filtering matrix, whose

i

’th diagonal element can be constructed from filter factors defined as

Eq. (29)

F_{i} = \frac{1}{\exp (σ_{i}^{2} / λ) - 1} .

Thus, $i$ ’th diagonal element of $E$ can be written as

Eq. (30)

E_{i} = \frac{σ_{i}^{2} / λ}{\exp (σ_{i}^{2} / λ) - 1} .

Substituting Eq. (27) in Eq. (8), the expression for error estimate becomes

Eq. (31)

η_{2} (λ) = \frac{{‖ r ‖}_{2} {‖ E x_{λ} ‖}_{2}}{{‖ AE x_{λ} ‖}_{2}} .

Now similar to the traditional Tikhonov method, one can use this expression in the framework of Algorithm 1 to determine the required regularization parameter (

λ

).

Local regularization parameter $λ^{*} = λ_{p}$ at which $η_{2} (λ_{p})$ is minimum. Follow steps 8 to 14 of Algorithm 1 to find required $λ$ .

3.4.

$ℓ_{1}$ Norm or Total Variation-Norm Based Regularization Method

The $ℓ_{1}$ and total variation (TV) norms are sparsity promoting norms, which have been previously utilized in various image reconstruction applications.²⁵^,²⁶ In this scheme, the cost function to be minimized becomes as

Eq. (32)

\underset{x \in R^{n^{2}}}{\arg \min} {‖ A x - b ‖}_{2}^{2} + λ {‖ X ‖}_{*},

where

{‖ X ‖}_{*}

is either

ℓ_{1}

norm or TV term.

X

is the matrix form of vectorized image

x

. For

ℓ_{1}

-norm case:

Eq. (33)

{‖ X ‖}_{*} = {‖ x ‖}_{1} .

The TV term is widely written as the

ℓ_{1}

norm of gradient of image (

x

). Hence, for the TV-norm case:

Eq. (34)

{‖ X ‖}_{*} = Σ_{i = 1}^{n^{2}} {‖ D_{i} x ‖}_{1},

where

D_{i}

denotes the local finite difference operator (with boundary condition) at pixel

i

and

Σ_{i} {‖ D_{i} x ‖}_{1}

is a discretization of the TV of

x

.²⁶ The solution

x

can be computed by minimizing Eq. (32).

The detailed description of these nonsmooth regularization schemes is beyond the scope of this study and these details can be found in Refs. 25 26.–27. The exact algorithm that was utilized for finding the regularization parameter for $ℓ_{1}$ - or TV-norm regularization schemes was provided in Algorithm 2. Note that the initial step is similar to step 1 of Algorithm 1, i.e., the first step was to start with the $i$ grid values chosen in interval [ $\log (λ_{\min})$ , $\log (λ_{\max})$ ].

Algorithm 2

Algorithm for evaluation of regularization parameter for ℓ1- or TV-norm based regularization method.

1 Input: $A; b; [λ_{\min}, λ_{\max}]; x_{o}$ ;.

2 Output: local regularization parameter: $λ^{*}$ .

1.

{while}^{1}

η_{2} (λ_{p}) < η_{2} (λ_{p + 1})

(i.e., step 7 of Algorithm 1)

2.

{while}^{2}

\frac{Ω_{k} - Ω_{k + 1}}{Ω_{k} + Ω_{k + 1}} \geq Tol (10^{- 4})

(minimization of cost function

Ω

)

• Compute objective function for current iteration (

k

):

Eq. (35)

Ω_{k} = {‖ A x_{k} - b ‖}_{2}^{2} + λ_{p} {‖ X_{k} ‖}_{*}

(

X_{k}

is the matrix form (dimension:

n \times n

) of vectorized image

x_{k}

(dimension:

n^{2} \times 1

), with

{‖ . ‖}_{*}

being either

ℓ_{1}

[Eq. (33)] or TV [Eq. (34)])

• i.

x_{k} = x_{k - 1} + A^{T} (b - A x_{k - 1})

ii. Form the matrix

X_{k}

by reshaping

x_{k}

iii. SVD:

X_{k} {= U}^{'} {SV}^{' T}

iv. Soft threshold of the singular values

\hat{S} = soft (diag (S), λ_{p} / 2)

v.

X_{k + 1} {= U}^{'} {\hat{S} V}^{' T}

. Form

x_{k + 1}

by vectorizing

X_{k + 1}

vi. Update

k = k + 1

• Compute objective function for next iterate:

Ω_{k + 1} = {‖ A x_{k + 1} - b ‖}_{2}^{2} + λ_{p} {‖ X_{k + 1} ‖}_{*}

•

{End while}^{2}

3.

r_{λ} = b - A x_{k + 1}

. Compute

η_{2} (λ_{p}) = \frac{{‖ r_{λ} ‖}_{2} {‖ A^{T} r_{λ} ‖}_{2}}{{‖ {AA}^{T} r_{λ} ‖}_{2}}

4. Update

p = p + 1

5.

{End while}^{1}

Then, the regularized solution (which minimized cost function) and the error estimate ( $η_{2}$ ) was computed. For this step, the algorithm contained two loops. The outer loop started from a higher value of $λ$ , and the inner loop solved Eq. (32). The outer loop progressively reduced the value of $λ$ until the condition in step 7 of Algorithm 1 is met. The error estimate $η_{2} (λ)$ was evaluated using the following expression:

Eq. (36)

η_{2} (λ) = \frac{{‖ r_{λ} ‖}_{2} {‖ A^{T} r_{λ} ‖}_{2}}{{‖ {AA}^{T} r_{λ} ‖}_{2}} .

The suitable $λ^{*}$ was selected as the one, which minimized $η_{2} (λ)$ . Then, to accurately refine the value of $λ$ , we iteratively added points to the grid around the minima $λ^{*}$ , as described in steps 8 to 14 of Algorithm 1.

4. Numerical Experiments

Numerical experiments are ideal for comparing the quantitative accuracies and efficiency of the various methods as measuring actual initial pressure rise in experimental phantom is very challenging. In this work, three numerical phantoms were chosen to justify effectiveness of the proposed model. Since PA imaging is widely used for visualizing internal blood vessel structures, first is a numerical blood vessel network phantom with an initial pressure rise of 1 kPa. Another numerical phantom is a variation of the Derenzo phantom, which consists of small and large size circular targets distributions over the imaging region. Both these phantoms dimensions were $201 \times 201$ .

A bigger dimension, $301 \times 301$ , the blood vessel phantom was considered as a third phantom. In this case, 100 ultrasound detectors were utilized, thus making $A$ as big as $50,000 \times 90,601$ . The size of vector $b$ became $90,601 \times 1$ . These phantoms are shown in Figs. 1(a), 1(e), and 1(i). The SNR of simulated data was kept at 40 dB to closely mimic the experimental data.

The main focus of this study was to compare the proposed LSQR method with Shaw et al.’s¹⁷ regularized LSQR method. This comparison was performed via running numerical simulations for the three phantoms introduced here. The study also showed that the proposed error estimate can be utilized with any regularization framework and to support the claim image reconstruction routines for Tikhonov regularization, exponential filtering, $ℓ_{1}$ -norm based regularization, and TV regularization frameworks were performed. A Linux workstation with dual eight-core Intel Xeon processors having a speed of 3.10 GHz with 128 GB RAM was used for all computations performed in this work.

Fig. 1

The numerical phantoms used in this study. (a) Blood vessel network (dimension: $201 \times 201$ ), (e) Derenzo phantom (dimension: $201 \times 201$ ), and (i) blood vessel phantom (dimension: $301 \times 301$ ). Reconstructed images of these phantoms using (b,f,j) regularized LSQR method¹⁷ and (c,g,k) proposed method. (d,h,l) A 1-D cross-sectional plot for the reconstructed images along the dotted line in target phantoms. The corresponding computation time and quantitative metric are listed in Table 1.

5. Results and Discussion

The minimization of the error estimate ( $η_{2}$ ) has been illustrated in Fig. 2. The circles show where ${\bar{η}}_{2} (λ)$ was evaluated [Fig. 2(b)]. Adjacent circles were connected for visual clarity and the lines have no significance. The residue in the figure was minimized to find $λ$ in Ref. 17. The local minima of residue curve [Fig. 2(a)] was situated close to the local minima of $η_{2}$ curve [Fig. 2(b)]. Once regularization is determined, the desired solution can be computed using Eq. (15).

For a 60 detectors’ geometry PA system matrix [Figs. 1(a) and 1(e)], the required number of Lanczos bidiagonalization steps ( $k$ ) turns out to be 44. This makes the size of matrix $B_{k}^{T}$ as $45 \times 44$ (see Ref. 20), making computation of ${(B^{T} B + λ I)}^{- 1}$ , which has a dimension of $45 \times 45$ , inexpensive. The number of Lanczos bidiagonalization steps ( $k$ ) required for the 100-detector case [Fig. 1(i)] was 77. The comparison of computation time between utilization of $201 \times 201$ dimension phantoms versus $301 \times 301$ dimension phantom is listed in Table 1.

Fig. 2

The minimization of (a) residue ( ${‖ b - A x ‖}_{2}^{2}$ ) for $k = 44$ and (b) error estimate ( $η_{2}$ ) [Eq. (8)] through Algorithm 1 for the phantom, as shown in Fig. 1(a).

Table 1

Performance comparison of two methods, regularized LSQR versus proposed methods for the results presented in Fig. 1.

Phantom	Regularized LSQR method17			Proposed
Phantom	Time (s)	λ	PC	Time (s)	λ	PC
Blood vessel network	$\sim 412$	0.0423	0.59	$\sim 90$	0.0594	0.58
Derenzo phantom	$\sim 424$	0.0158	0.69	$\sim 96$	0.0205	0.76
Blood vessel phantom	$\sim 1871$	0.0890	0.534	$\sim 526$	0.0672	0.529

The reconstruction results obtained for the blood vessel network using the regularized LSQR method¹⁷ and the proposed method are shown in Figs. 1(b) and 1(c). A one-dimensional cross-sectional plot for the reconstructed image along the red dotted line in the target image is shown in Fig. 1(d). This plot shows that the performance of the proposed method in terms of quantitative accuracy of the reconstructed image is similar to the state-of-the-art regularized LSQR method. Figures 1(f)–1(h) show the similar results for the Derenzo phantom, which testify the ability to reconstruct the objects of different shapes and sizes. The results for the larger phantom ( $301 \times 301$ dimension, 100 detector points) are shown in Figs. 1(i)–1(l).

To check the performance of the image quality, a figure of merit, such as Pearson correlation (PC), was used for the comparison. PC is a quantitative metric that measures the degree of correlation between the target image and the reconstructed image. The PC is defined as¹⁸

Eq. (37)

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

where

x

is the expected initial pressure distribution,

x^{recon}

is the reconstructed initial pressure distribution, COV denotes the covariance, and

σ

denotes the standard deviation. Higher value of PC corresponds to the accurate detectability of the target (i.e., spatial fidelity).¹⁸^,²⁸

The proposed method is computationally much faster than previously proposed methods. Table 1 shows the comparison results, where it can be seen that there is little difference in PC, while the computation time is much lower for the proposed method ( $\sim 4.5 \times$ speed up) compared to the regularized LSQR method proposed by Shaw et al.¹⁷ Note that the one-time system matrix building time ( $\sim 181 s$ ) was excluded from the computation time, as listed in Tables 1 and 2.

Table 2

The computation time, regularization parameter, and PC for the results presented in Fig. 3.

Method	Computation	Blood vessel network		Derenzo phantom
Method	Time (s)	λ	PC	λ	PC
Traditional Tikhonov regularization	$\sim 221$	0.0035	0.58	0.0004	0.76
Exponential filtering²⁴	$\sim 210$	0.0077	0.57	0.0008	0.75
$l_{1}$ -norm based regularization²⁵	$\sim 617$	0.0010	0.63	0.0010	0.59
TV-norm based regularization²⁶	$\sim 581$	0.0012	0.64	0.0011	0.65

The proposed regularization method was also tested for other popular model-based reconstruction schemes, i.e., traditional Tikhonov regularization, exponential filtering of SVD-based²⁴ regularization, $ℓ_{1}$ -norm based regularization,²⁵ and TV-norm based regularization²⁶ technique to evaluate the regularization parameter and reconstruct the PA images. Figure 3 presents the reconstruction results for blood vessel network and Derenzo phantom obtained through these methods. Figures 3(a) and 3(b) show the reconstructed images obtained using Tikhonov regularization, and Figs. 3(c) and 3(d) show images obtained using recently proposed exponential filtering method.²⁴ The reconstruction results for nonsmooth regularization, such as $ℓ_{1}$ - and TV-norm based methods, are shown in Figs. 3(e)–3(f) and 3(g)–3(h), respectively. The $ℓ_{1}$ - and TV-norm based algorithm becomes expensive in computation for more iterations [Eq. (35)]. It is to be noted that the results presented for $ℓ_{1}$ - and TV-based reconstructions²⁵ are for 10 iterations, with $ℓ_{1}$ -norm based method being slightly more expensive in terms of computation (refer to Table 2). However, the reconstruction accuracy can be improved with a larger number of iterations. The algorithm to determine regularization parameter for nonsmooth regularization methods was described in Algorithm 2. As it can be seen, the proposed error estimate ( $η_{2}$ ) is capable of determining the required regularization parameter both in the smooth and non-smooth regularization frameworks. The regularization parameter determined using the proposed $η_{2}$ estimate and the image quality metrics, PC, has been listed in Table 2.

Fig. 3

Reconstructed images using (a and b) traditional Tikhonov regularization,²⁴ (c and d) exponential filtering of SVD-based reconstruction,²⁴ (e and f) $ℓ_{1}$ -norm based regularization,²⁵ and (g and h) TV-norm based regularization²⁶ for the two numerical phantoms [target phantoms are given in Figs. 1(a) and 1(e)], where the regularization parameter was determined using proposed error estimate minimization. The quantitative metrics are given in Table 2.

The regularization schemes that are utilized in this work include smooth and nonsmooth regularizations. Among the smooth regularization (Tikhonov) schemes, we have implemented the proposed scheme both in dimensionality reduction frame work (regularized LSQR) as well as traditional Tikhonov frame work. From Tables 1 and 2, it is clear that the regularized LSQR method is computationally more efficient and observed reconstruction performance (in terms of PC) for both methods is same as expected. It should be noted that the Lanczos bidiagnolization was not utilized for nonsmooth regularization schemes, including the presented $ℓ_{1}$ - and TV-norm based ones. These (nonsmooth regularization schemes) were implemented in traditional fashion (same as traditional Tikhonov) as the main aim of this work was to show that the proposed error estimates have utility for any regularization frame work (traditional/reduced dimension combined with smooth/nonsmooth).

The assessment of image quality in emission tomography as well as biological imaging has been previously performed through figure of merits, such as PC.¹⁸^,²⁸ Even though the task-based image quality metrics²⁹ can better assess the reconstructed image quality, the use of residual errors (least-squares metric) in comparison of reconstruction algorithms performance (including the parameters used in them) has been a common practice across tomographic problems.³⁰ The estimate of these residual errors has been utilized in this work to find the appropriate reconstruction parameters; here, it is being regularization parameter. The future scope of this study is to validate the performance of proposed algorithm through experimental studies, which can further assure the efficiency and usability of proposed method.

6. Conclusions

The LSQR model has been extensively used and shown to be effective¹⁵ in full-view data cases ( $m ≫ n^{2}$ ), where there is no requirement of regularization. However, in limited data ( $m ≪ n^{2}$ ) cases, Tikhonov framework has to be used for regularizing the solution.¹⁷ Determining a suitable regularization parameter is important as the solution is sensitive to this choice. This study addressed the problem of determining the regularization parameter automatically based on error estimate ( $η_{2}$ ) analysis without any prior information about the noise. Utilizing this with Lanczos bidiagonalization framework adds the advantage of reducing the dimensionality of the problem. This makes the proposed model highly computationally efficient and promising as compared to the other popular state-of-the-art techniques. The proposed method is generic in nature and can be used in any regularization framework to determine the regularization parameter. Here, it was also shown that the proposed method is 4.5 times faster compared to existing state-of-the-art methods.

Acknowledgments

Authors would like to thank the anonymous reviewers for their thoughtful comments, which improved this work. This work was supported by Department of Biotechnology (DBT) Innovative Young Biotechnologist Award (IYBA) (Grant No.: BT/07/IYBA/2013-13), the Department of Biotechnology (DBT) Rapid Grant for Young investigator (RGYI) (Grant No.: BT/PR6494/GDB/27/415/2012), and DBT Bioengineering Grant (Grant No.: BT/PR7994/ MED/32/284/2013). The work of M.B. was supported by National Mathematics Initiative (NMI).

References

1.

M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum., 77 41101 (2006). http://dx.doi.org/10.1063/1.2195024 RSINAK 0034-6748 Google Scholar

2.

R. Li et al., “Assessing breast tumor margin by multispectral photoacoustic tomography,” Biomed. Opt. Express, 6 (4), 1273 –1281 (2015). http://dx.doi.org/10.1364/BOE.6.001273 BOEICL 2156-7085 Google Scholar

3.

G. Diot, A. Dima and V. Ntziachristos, “Multispectral opto-acoustic tomography of exercised muscle oxygenation,” Opt. Lett., 40 (7), 1496 –1499 (2015). http://dx.doi.org/10.1364/OL.40.001496 OPLEDP 0146-9592 Google Scholar

4.

B. Z. Chen et al., “Photoacoustic imaging of cerebral hypoperfusion during acupuncture,” Biomed. Opt. Express, 6 (9), 3225 –3234 (2015). http://dx.doi.org/10.1364/BOE.6.003225 BOEICL 2156-7085 Google Scholar

5.

X. Li et al., “High resolution functional photoacoustic tomography of breast cancer,” Med. Phys., 42 (9), 5321 –5328 (2015). http://dx.doi.org/10.1118/1.4928598 Google Scholar

6.

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 Google Scholar

7.

T. N. Erpelding et al., “Sentinel lymph nodes in the rat: noninvasive photoacoustic and US imaging with a clinical US system,” Radiology, 256 (1), 102 –110 (2010). http://dx.doi.org/10.1148/radiol.10091772 RADLAX 0033-8419 Google Scholar

8.

M. Nasiriavanaki et al., “High-resolution photoacoustic tomography of resting-state functional connectivity in the mouse brain,” Proc. Natl. Acad. Sci. U. S. A., 111 (1), 21 –26 (2014). http://dx.doi.org/10.1073/pnas.1311868111 PNASA6 0027-8424 Google Scholar

9.

P. van Es et al., “Initial results of finger imaging using photoacoustic computed tomography,” J. Biomed. Opt., 19 (6), 060501 (2014). http://dx.doi.org/10.1117/1.JBO.19.6.060501 JBOPFO 1083-3668 Google Scholar

10.

L. A. Kunyansky, “Explicit inversion formulae for the spherical mean Radon transform,” Inverse Probl., 23 (1), 373 –383 (2007). http://dx.doi.org/10.1088/0266-5611/23/1/021 Google Scholar

11.

Y. Xu, D. Feng and L. V. Wang, “Exact frequency-domain reconstruction for thermoacoustic tomography: I. Planar geometry,” IEEE Trans. Med. Imaging, 21 (7), 823 –828 (2002). http://dx.doi.org/10.1109/TMI.2002.801172 ITMID4 0278-0062 Google Scholar

12.

Y. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E, 71 16706 (2005). http://dx.doi.org/10.1103/PhysRevE.71.016706 Google Scholar

13.

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 JASMAN 0001-4966 Google Scholar

14.

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

15.

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

16.

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

17.

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

18.

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

19.

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

20.

P. C. Hansen, Rank Deficient and Discrete ill-Posed Problems, Society for Industrial and Applied Mathematics, Philadelphia (1998). Google Scholar

21.

C. Brezinski, G. Rodriguez and S. Seatzu, “Error estimates for the regularization of least squares problems,” Numer. Algorithms, 51 (1), 61 –76 (2009). http://dx.doi.org/10.1007/s11075-008-9243-2 NUALEG 1017-1398 Google Scholar

22.

L. Reichel, G. Rodriguez and S. Seatzu, “Error estimates for large-scale ill-posed problems,” Numer. Algorithms, 51 (3), 341 –361 (2009). http://dx.doi.org/10.1007/s11075-008-9244-1 NUALEG 1017-1398 Google Scholar

23.

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

24.

M. Bhatt, S. Gutta and P. K. Yalavarthy, “Exponential filtering of singular values improves photoacoustic image reconstruction,” J. Opt. Soc. Am. A, 33 (9), 1785 –1792 (2016). http://dx.doi.org/10.1364/JOSAA.33.001785 JOAOD6 0740-3232 Google Scholar

25.

A. Majumdar and R. K. Ward, “Nuclear norm-regularized SENSE reconstruction,” J. Magn. Reson. Imaging, 30 (2), 213 –221 (2012). http://dx.doi.org/10.1016/j.mri.2011.09.014 Google Scholar

26.

Y. H. Xiao and J. F. Yang, “Alternating algorithms for total variation image reconstruction from random projections,” Inverse Probl. Imaging, 6 (3), 547 –563 (2012). http://dx.doi.org/10.3934/ipi Google Scholar

27.

A. Majumdar and R. K. Ward, “An algorithm for sparse MRI reconstruction by Schatten p-norm minimization,” Magn. Reson. Imaging, 29 (3), 408 –417 (2011). http://dx.doi.org/10.1016/j.mri.2010.09.001 MRIMDQ 0730-725X Google Scholar

28.

J. Kuntz et al., “Constrained reconstructions for 4-D intervention guidance,” Phys. Med. Biol., 58 (10), 3283 –3300 (2013). http://dx.doi.org/10.1088/0031-9155/58/10/3283 PHMBA7 0031-9155 Google Scholar

29.

A. Petschke and P. J. La Rivire, “Comparison of photoacoustic image reconstruction algorithms using the channelized hotelling observer,” J. Biomed. Opt., 18 (2), 026009 (2013). http://dx.doi.org/10.1117/1.JBO.18.2.026009 JBOPFO 1083-3668 Google Scholar

30.

W. S. Phillips and M. C. Fehler, “Traveltime tomography: a comparison of popular methods,” Geophysics, 56 (10), 1639 –1649 (1991). http://dx.doi.org/10.1190/1.1442974 GPYSA7 0016-8033 Google Scholar

Biography

Manish Bhatt received his BTech degree from the National Institute of Technology, Hamirpur, Himachal Pradesh, India, in 2011. Currently, he is a graduate student at the National Mathematics Initiative, Indian Institute of Science (IISc), Bangalore, India.

Atithi Acharya is currently an undergraduate student in the Department of Physics, Indian Institute of Science (IISc), Bangalore, India. His research interests include medical physics, biomedical optics, and scientific computing.

Phaneendra K. Yalavarthy received his MSc degree in engineering from the Indian Institute of Science, Bangalore, India, and his PhD degree in biomedical computation from Dartmouth College, Hanover, New Hampshire, in 2007. He is an associate professor in the Department of Computational and Data Sciences, Indian Institute of Science, Bangalore. His research interests include medical image computing, medical image analysis, and biomedical optics.

Citation Download Citation

Manish Bhatt, Atithi Acharya, and Phaneendra K. Yalavarthy "Computationally efficient error estimate for evaluation of regularization in photoacoustic tomography," Journal of Biomedical Optics 21(10), 106002 (20 October 2016). https://doi.org/10.1117/1.JBO.21.10.106002

Published: 20 October 2016

Access the abstract

JOURNAL ARTICLE
9 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 10 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Error analysis

Reconstruction algorithms

Biological research

Image restoration

Photoacoustic tomography

Blood vessels

Data modeling

1.

Introduction

2.

Forward Problem in Photoacoustic Tomography: System Matrix Approach

Eq. (1)

3.

Error Estimate for a Regularized Solution

Eq. (2)

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

Algorithm 1

Eq. (11)

3.1.

Regularized LSQR Method

Eq. (12)

Eq. (13)

Eq. (14)

Eq. (15)

Eq. (16)

3.2.

Traditional Tikhonov Regularization

Eq. (17)

Eq. (18)

Eq. (19)

Eq. (20)

Eq. (21)

Eq. (22)

Eq. (23)

Eq. (24)

3.3.

Exponential Filtering or Showalter’s Method

Eq. (25)

Eq. (26)

Eq. (27)

Eq. (28)

Eq. (29)

Eq. (30)

Eq. (31)

3.4.

ℓ1 Norm or Total Variation-Norm Based Regularization Method

Eq. (32)

Eq. (33)

Eq. (34)

Algorithm 2

Eq. (35)

Eq. (36)

4.

Numerical Experiments

Fig. 1

5.

Results and Discussion

Fig. 2

Table 1

Eq. (37)

Table 2

Fig. 3

6.

Conclusions

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years

$ℓ_{1}$ Norm or Total Variation-Norm Based Regularization Method