Fast spatiotemporal image reconstruction based on low-rank matrix estimation for dynamic photoacoustic computed tomography

Kun Wang; Jun Xia; Changhui Li; Lihong V. Wang; Mark A. Anastasio

doi:10.1117/1.JBO.19.5.056007

9 May 2014 Fast spatiotemporal image reconstruction based on low-rank matrix estimation for dynamic photoacoustic computed tomography

Kun Wang, Jun Xia, Changhui Li, Lihong V. Wang, Mark A. Anastasio

Author Affiliations +

Journal of Biomedical Optics, Vol. 19, Issue 5, 056007 (May 2014). https://doi.org/10.1117/1.JBO.19.5.056007

Abstract

In order to monitor dynamic physiological events in near-real time, a variety of photoacoustic computed tomography (PACT) systems have been developed that can rapidly acquire data. Previously reported studies of dynamic PACT have employed conventional static methods to reconstruct a temporally ordered sequence of images on a frame-by-frame basis. Frame-by-frame image reconstruction (FBFIR) methods fail to exploit correlations between data frames and are known to be statistically and computationally suboptimal. In this study, a low-rank matrix estimation-based spatiotemporal image reconstruction (LRME-STIR) method is investigated for dynamic PACT applications. The LRME-STIR method is based on the observation that, in many PACT applications, the number of frames is much greater than the rank of the ideal noiseless data matrix. Using both computer-simulated and experimentally measured photoacoustic data, the performance of the LRME-STIR method is compared with that of conventional FBFIR method followed by image-domain filtering. The results demonstrate that the LRME-STIR method is not only computationally more efficient but also produces more accurate dynamic PACT images than a conventional FBFIR method followed by image-domain filtering.

1. Introduction

Photoacoustic computed tomography (PACT), also known as optoacoustic tomography, is an emerging imaging modality that holds great promise for a wide range of biomedical imaging applications.¹^–³ PACT is a hybrid imaging modality that combines the high spatial resolution of ultrasound imaging and the high soft tissue contrast of optical imaging.⁴^–⁶ In PACT, a short laser pulse is employed to illuminate an object and internal acoustic wavefields are produced via the thermoacoustic effect. The acoustic wavefields propagate out of the object and are measured by using ultrasonic transducers.⁴^,⁵^,⁷ From the measured wavefield data, an image that depicts the absorbed optical energy density distribution within the object, hereafter referred to as the object function, is produced by using a reconstruction algorithm. The vast majority of PACT image reconstruction algorithms developed to date assume static imaging conditions, in which the sought-after object function is independent of time.

In PACT studies that involve dynamic physiological processes, the object function changes with time. The goal of dynamic PACT⁸^–¹¹ is to reconstruct a sequence of object function estimates that correspond to a collection of time points. These temporal samples of the object function will be referred to as object frames. When the spatial support of the object function is two-dimensional (2-D) or three-dimensional (3-D), the problem can be interpreted as 3-D or four-dimensional (4-D) tomographic imaging, where time is counted as a dimension. At each temporal sample in a dynamic PACT study, a static PACT data set is recorded. This data set will be referred to as a data frame. It is assumed that the object function remains static during acquisition of each data frame, which can be approximately satisfied if the temporal resolution of the imaging system is sufficiently high. An estimate of each object frame can be reconstructed from the corresponding data frame by using a conventional static PACT reconstruction algorithm. This image reconstruction approach will be referred to as a frame-by-frame image reconstruction (FBFIR) method and has been utilized in the vast majority of previous studies of dynamic PACT.⁹^–¹² As demonstrated in mature imaging modalities,¹³^–¹⁹ a limitation of FBFIR methods is that they fail to exploit correlations between data frames and are, therefore, known to be statistically and computationally suboptimal. For example, because each object frame is solely determined by its corresponding data frame, in order to avoid strong reconstruction artifacts, each data frame has to be densely sampled. This places limitations on the temporal resolution of the system. The fact that object frames are independently reconstructed in FBFIR methods also causes them to be computationally burdensome, especially if iterative image reconstruction algorithms are employed for 4-D dynamic PACT applications.¹¹^,²⁰^,²¹ Moreover, because they fail to exploit temporal information contained in the measured data frames, FBFIR methods do not optimally mitigate the effects of measurement noise.¹³

Several image reconstruction methods have been proposed for dynamic PACT.⁹^–¹² Gamelin et al. proposed to synthesize a densely sampled data frame by estimating the pressure data at the locations between transducers.²² It was demonstrated that the use of the method could improve the temporal resolution of a PACT imaging system by permitting sufficiently accurate image reconstruction from data frames of reduced sizes. However, being an FBFIR method, it is subject to the limitations discussed above. In a different study, an image-domain Karhunen Loève (KL) filter [i.e., principal component analysis (PCA) filter] and an independent component analysis filter were applied to the images reconstructed by using an FBFIR method.²³

Unlike FBFIR methods, spatiotemporal image reconstruction (STIR) methods for dynamic tomography modalities jointly reconstruct the sequence of object frame estimates by using all of the data frames. These algorithms exploit statistical correlations or deterministic linear dependency between either data or image frames.²⁴^,²⁵ In this way, STIR methods can circumvent the limitations of FBFIR methods that are described above. For example, statistical correlations can be exploited by Bayesian estimation approaches¹⁴^,²⁶ or data-domain KL filtering¹³^,¹⁶^,²⁷ to improve reconstructed image quality and/or reduce the computational burden as compared to FBFIR methods. Recently studied matrix completion methods exploit the deterministic linear dependency of data matrices¹⁵^,¹⁷^,¹⁸^,²⁸^,²⁹ to achieve the same goal. A variety of advanced STIR methods have been developed for use with established dynamic tomography modalities.¹³^,¹⁷^,²⁹^,³⁰ When utilized in conjunction with optimized data-acquisition protocols, STIR methods can improve the temporal resolution of an imaging system, so that rapid dynamic physiological events can be visualized.¹⁵^,²⁹ Despite these advantages, the development and investigation of STIR methods for dynamic PACT remains largely unexplored.

In this study, a low-rank matrix estimation-based STIR (LRME-STIR) method is investigated for dynamic PACT applications. The LRME-STIR method is based on the observation that, in many PACT applications, the number of frames is much greater than the rank of the ideal noiseless data matrix. Our work is inspired by the successful application of similar ideas in dynamic magnetic resonance imaging.¹⁵^,¹⁷ The remainder of the article is organized as follows. The static and dynamic PACT image reconstruction problems are reviewed in Sec. 2. Section 3 describes the LRME-STIR method for PACT. A description of the numerical and experimental studies is provided in Sec. 4. Section 5 contains the results of these studies, and the paper concludes with a discussion in Sec. 6.

2. Background

2.1.

Static Image Reconstruction in Conventional PACT

A static PACT imaging system can be accurately described by a continuous-to-discrete (C-D) imaging model as⁵^,⁷^,²⁰^,³¹^,³²

Eq. (1)

{[g]}_{j I + i} = h^{e} (t) *_{t} \frac{1}{Ω_{j}} \int_{Ω_{j}} d r^{'} \frac{β}{4 π C_{p}} \int_{V} d r A^{s} (r) \frac{d}{d t} \frac{δ (t - \frac{| r^{'} - r |}{c_{0}})}{| r^{'} - r |} |_{t = i Δ T}, \binom{i = 0, 1, \dots, I - 1}{j = 0, 1, \dots, J - 1},

where

h^{e} (t)

is the electrical impulse response (EIR) of the transducer,³²^,³³

*_{t}

denotes the temporal convolution operation,

δ (t)

is the one-dimensional Dirac delta function, and

β

,

c_{0}

, and

C_{p}

denote the thermal coefficient of volume expansion, the (constant) speed-of-sound, and the specific heat capacity of the medium at constant pressure, respectively. The static object function,

A^{s} (r)

, is assumed to be bounded and contained within the volume

V

. The quantities

r

and

r^{'}

specify the spatial coordinates within

V

and on the measurement surface, respectively. The vector

g \in R^{I J}

represents a lexicographically ordered collection of sampled values of the electrical signals produced by the ultrasonic transducers, where

J

and

I

denote the number of transducers employed in the imaging system and the number of temporal samples recorded by each transducer, respectively. The notation

{[g]}_{j I + i}

will be utilized to denote the

(j I + i)

’th element of

g

. Here, the integer-valued indices

j

and

i

describe the transducer location and temporal sample, respectively. The temporal sampling interval is

Δ_{t}

, and

Ω_{j}

denotes the detection area of the

j

’th transducer, which is assumed to be a subset of the measurement surface.

The sought-after object function is related to the optical absorption coefficient, $μ_{a}^{s} (r)$ , and optical fluence, $Φ^{s} (r)$ , within the object as $A^{s} (r) = μ_{a}^{s} (r) Φ^{s} (r)$ . Note that the superscript $s$ indicates that these quantities are static. While quantification of $μ_{a}^{s} (r)$ has been actively investigated,³⁴^,³⁵ an estimate of $A^{s} (r)$ represents the quantity produced in the majority of current PACT studies.

Based on Eq. (1), iterative image reconstruction algorithms have been developed for estimation of $A^{s} (r)$ .²⁰^,³²^,³⁶ When the transducer size is small and/or the object is located near the center of a relatively large measurement geometry, the surface integral over $Ω_{j}$ can be neglected. In these cases, a variety of analytic formulae³⁷^–⁴⁰ can also be employed for image reconstruction after deconvolving the EIR from the recorded signals.⁴¹ Linear static image reconstruction algorithms will be described as

Eq. (2)

{\hat{A}}^{s} = B g,

where the matrix

B \in R^{N \times I J}

represents a discrete image reconstruction operator and

{\hat{A}}^{s} \in R^{N}

is the reconstructed digital image arranged in a lexicographical order with

N

being the number of pixels or voxels.

2.2.

Dynamic PACT and Frame-by-Frame Image Reconstruction

In dynamic PACT, the optical absorption coefficient $μ_{a} (r, t^{s})$ and the optical fluence $Φ (r, t^{s})$ are time-dependent. Note that the optical fluence $Φ (r, t^{s})$ varies over time because of the temporal variation of $μ_{a} (r, t^{s})$ .⁴ The time-dependent object function is accordingly defined as $A (r, t^{s}) \equiv Φ (r, t^{s}) μ_{a} (r, t^{s})$ . Here and throughout the manuscript, $t^{s}$ (or the corresponding index $k$ defined below) will be referred to as a slow-time (i.e., the time of each frame) coordinate, while the time coordinate $t$ (or the corresponding index $i$ ) in Eq. (1) will be referred to as a fast-time (i.e., the arrival time of acoustic signals) coordinate.

The $k$ ’th frame of the dynamic object is defined as ${A_{k} (r) \equiv Φ (r, t^{s}) μ_{a} (r, t^{s}) |}_{t^{s} = k Δ_{t^{s}}}$ , for $k = 0, 1, \dots, K - 1$ , where $K$ denotes the number of temporal samples, indexed by $k$ , with temporal sampling interval $Δ_{t^{s}}$ . Replacing $A^{s} (r)$ in Eq. (1) by $A_{k} (r)$ yields the $k$ ’th data frame denoted by $g_{k} \in R^{I J}$ , for $k = 0, 1, \dots, K - 1$ .

The goal of dynamic PACT is to estimate the collection of object frames ${A_{k} (r)}_{k = 0}^{K - 1}$ from the measured data frames of ${g_{k}}_{k = 0}^{K - 1}$ . To accomplish this, a linear FBFIR method operates as

Eq. (3)

{\hat{A}}_{k} = {B g}_{k}, for k = 0, 1, \dots, K - 1,

where

{\hat{A}}_{k} \in R^{N}

represents a finite-dimensional estimate of the

k

‘th object frame

A_{k} (r)

. Note that

B

has to be applied

K

times to obtain the sought-after sequence of image estimates

{{\hat{A}}_{k}}_{k = 0}^{K - 1}

. We define the matrices

\hat{A} \in R^{N \times K}

and

G \in R^{I J \times K}

as

Eq. (4)

\hat{A} = [{\hat{A}}_{0} | {\hat{A}}_{1} | \dots | {\hat{A}}_{K - 1}], and G = [g_{0} | g_{1} | \dots | g_{K - 1}],

where

G

will be referred to as the data matrix. In terms of these matrices, Eq. (3) can be expressed as

Eq. (5)

\hat{A} = B G .

3. Low-Rank Matrix Estimation-Based STIR for Dynamic PACT

3.1.

Singular Value Decomposition-Based STIR

Consider the singular-value decomposition (SVD)²⁴ of the data matrix $G$ :

Eq. (6)

G = \sum_{k = 0}^{R - 1} μ_{k} v_{k} u_{k}^{†},

where

R

is the rank of

G

, the superscript

†

denotes the matrix adjoint, and

{u_{k}, v_{k}, μ_{k}^{2}}_{k = 0}^{R - 1}

is the associated singular system that satisfies

Eq. (7)

G^{†} {G u}_{k} = μ_{k}^{2} u_{k} and {G G}^{†} v_{k} = μ_{k}^{2} v_{k} .

The orthonormal vectors

{u_{k}}

and

{v_{k}}

are related as

Eq. (8)

v_{k} = \frac{1}{μ_{k}} {G u}_{k}, for k = 0, 1, \dots, R - 1 .

Here and throughout the manuscript, the singular values will be labeled in descending order, i.e.,

μ_{0} \geq μ_{1} \geq \dots \geq μ_{R - 1} > 0

.

On substitution from Eq. (6) into Eq. (5), one obtains

Eq. (9)

\hat{A} = \sum_{k = 0}^{R - 1} B ({G u}_{k}) u_{k}^{†} .

Equation (9) describes an STIR method that is mathematically equivalent to the FBFIR method in Eq. (3) or (5), but, as investigated below, can be computationally more efficient. Note that, by using resolution of the identity on the canonical basis for the slow-time coordinate, Eq. (5) can be expressed as

Eq. (10)

\hat{A} = \sum_{k = 0}^{K - 1} B ({G e}_{k}) e_{k}^{†},

where

{e_{k}}_{k = 0}^{K - 1}

denotes the canonical basis of

R^{K}

. Equations (9) and (10) reveal that the reconstruction operator

B

needs to be applied

R

and

K

times for the STIR and FBFIR methods, respectively. Since

R \leq K

, this indicates that the STIR method is computationally more efficient than an FBFIR method by a factor of

K / R

.

For many PACT applications, the data matrix $G$ is likely to be of low rank with $R ≪ K$ . Based on a discrete-to-discrete approximation⁷^,²⁴ of the C-D model in Eq. (1), it has been demonstrated that for a linear imaging system, the rank of $G$ is determined by the rank of the object function’s finite-dimensional representation,¹⁵^,⁴² denoted by $A \in R^{N \times K}$ . When pixels are employed, $A$ can be defined as⁷^,²⁴

Eq. (11)

{[A]}_{n, k} = A_{k} (r_{n}), \begin{array}{l} n = 0, 1, \dots, N - 1 \\ k = 0, 1, \dots, K - 1 \end{array},

with

r_{n}

denoting the location of the

n

’th pixel. Each row of

A

will be referred to as a pixel time activity curve (TAC). Intuitively, in the same anatomical structure, the pixel TACs are likely to be linearly dependent. For certain applications, even different organs may possess linearly dependent TACs. If the object to be imaged contains only a small number of organs with independent TACs, the rank of

A

will be small compared with the number of slow-time frames, as will the rank of

G

.¹⁵^,⁴² In functional brain imaging,⁹ for example, the pixel TACs can often be grouped into three clusters that correspond to the cortical arteries, the veins, and other static background, suggesting a rank of 3 for

G

. This idealized scenario neglects several factors, including movement, vessel expansion/contraction, noise, and other uncertainties. These issues are discussed in Sec. 6.

3.2.

Low-Rank Regularized Data Matrix Denoising

The derivation of the STIR formula in Eq. (9) did not consider data noise. Because of measurement noise, the measured data matrix, denoted by $\underset{̲}{G}$ , is likely to be of full rank, i.e., $\underset{̲}{R} = K$ , where $\underset{̲}{R}$ is the rank of $\underset{̲}{G}$ . Simply replacing $G$ in Eq. (9) with a full-rank data matrix $\underset{̲}{G}$ will expropriate the computational advantages of the STIR method. Furthermore, directly applying the STIR formula with noisy measurements may result in image artifacts, just as would occur in FBFIR. These artifacts can be mitigated by incorporating regularization in the static image reconstruction operator $B$ .²⁰^,³⁶ Alternatively, an explicit denoising approach can be employed to estimate the noise-free $G$ , from which $A$ can be reconstructed by using algorithms developed for idealized noise-free measurements,⁴¹ such as Eq. (9). In this study, the latter strategy is employed.

Hereafter, random quantities are underlined. The measured data matrix $\underset{̲}{G}$ can be expressed as

Eq. (12)

\underset{̲}{G} = G + \underset{̲}{N},

where

\underset{̲}{N}

is a random noise matrix of dimensions

I J \times K

and

G

is the noiseless data matrix. Estimation of

G

from

\underset{̲}{G}

is a classic image denoising problem. A variety of image denoising algorithms can be employed for this task, including recently studied sparsity-regularized denoising methods.⁴¹^,⁴³^–⁴⁵ In this study, denoising of

\underset{̲}{G}

was performed by solving the following optimization problem:¹⁵^,⁴⁶^,⁴⁷

Eq. (13)

\hat{G} = \arg \min_{G} \frac{1}{2} {∥ \underset{̲}{G} - G ∥}_{F}^{2} + β Rank (G),

where

Rank (G)

denotes the rank of

G

,

β

is a regularization parameter, and

{∥ \cdot ∥}_{F}

denotes the Frobenius norm. The Frobenius norm is defined as the square root of the sum of the absolute squares of the matrix’s elements and can be viewed as an extension of the

ℓ_{2}

vector norm to matrices. The penalty term

Rank (G)

will promote solutions

\hat{G}

that are of low rank. Equation (13), which will be referred to as the LRME problem, possesses a closed-form solution that can be calculated via the singular value hard thresholding (SVHT) estimator as⁴⁸

Eq. (14)

\hat{G} = \sum_{k = 0}^{\hat{R} - 1} {\underset{̲}{μ}}_{k} {\underset{̲}{v}}_{k} {\underset{̲}{u}}_{k}^{†},

where

{{\underset{̲}{u}}_{k}, {\underset{̲}{v}}_{k}, {\underset{̲}{μ}}_{k}^{2}}

denotes the singular system of

\underset{̲}{G}

and

\hat{R}

is the maximum of

{k | {\underset{̲}{μ}}_{k} > β}

.

Note that the SVD of $\underset{̲}{G}$ , in general, can be efficiently calculated because the number of slow-time frames is often $< 100$ for many dynamic PACT applications.⁹^–¹¹ In addition, each row vector of $\hat{G}$ can be interpreted as a linear combination of the row vectors of $\underset{̲}{G}$ , as shown in Appendix A.

Multiple optimal estimators of the regularization parameter $β$ in Eq. (13), or the equivalent $\hat{R}$ in Eq. (14), have been proposed based on various criteria, including minimizing the asymptomatic mean square error⁴⁷ and minimizing the estimation risk.⁴⁸ Most of these studies assume a random white noise model with a known noise level. In practice, there are multiple sources of noise, including acoustic noise that can be object-dependent⁴⁹ and, therefore, difficult to model. Also, these methods only ensure that $\hat{G}$ is an optimal estimate of $G$ , which does not necessarily yield an optimal estimate ${{\hat{A}}_{k}}$ . For example, we have tested the optimal threshold value proposed in Ref. 47. Though not included here, our computer-simulation studies suggest that the method nearly always provides an optimal $\hat{G}$ that minimizes ${∥ \hat{G} - G ∥}_{F}^{2}$ . However, the images reconstructed from the optimal $\hat{G}$ possess larger errors than do the corresponding results presented in this article. Therefore, in this study, the regularization parameter was empirically selected for each data set, respectively, as described in Sec. 5.

3.3.

LRME-Based STIR Method

Because both are conducted in the singular system of data matrices, the low-rank regularized data matrix denoising and the STIR method can be naturally combined as a two-step LRME-STIR method. Because we employed the SVHT estimator in Eq. (14), we obtained not only a low-rank estimate $\hat{G}$ but also its singular system as ${{\underset{̲}{u}}_{k}, {\underset{̲}{v}}_{k}, {\underset{̲}{μ}}_{k}^{2}}_{k = 0}^{\hat{R} - 1}$ with $\hat{R} ≪ K$ . The implementation procedure is summarized as follows: (1) Compute the SVD of the data matrix $\underset{̲}{G}$ . This can be solved by many available numerical libraries.⁵⁰ The resultant ${{\underset{̲}{u}}_{k}, {\underset{̲}{v}}_{k}, {\underset{̲}{μ}}_{k}^{2}}_{k = 0}^{K - 1}$ are stored in memory. (2) Estimate $\hat{G}$ by using the SVHT estimator, i.e., Eq. (14). Many approaches have been proposed to estimate the optimal threshold value $β$ or, equivalently, the optimal choice of $\hat{R}$ .⁴⁷^,⁴⁸ A brief discussion will be provided in Sec. 6. (3) Apply a static image reconstruction algorithm to the first $\hat{R}$ singular components

Eq. (15)

{\tilde{A}}_{k} = {\underset{̲}{μ}}_{k} B {\underset{̲}{v}}_{k}, for k = 0, 1, \dots, \hat{R} - 1,

where

{{\tilde{A}}_{k}}_{k = 0}^{\hat{R} - 1}

provide an estimate of the object function in the coordinate system specified by

{{\underset{̲}{v}}_{k}}_{k = 0}^{\hat{R} - 1}

. (4) Project

{{\tilde{A}}_{k}}_{k = 0}^{\hat{R} - 1}

onto the canonical coordinate system as

Eq. (16)

\hat{A} = \sum_{k = 0}^{\hat{R} - 1} {\tilde{A}}_{k} {\underset{̲}{u}}_{k}^{†} .

4. Numerical Studies

Computer-simulation and experimental data studies were conducted to demonstrate the use of the LRME-STIR method for dynamic PACT reconstruction. The performance of the LRME-STIR was compared with an FBFIR method followed by image-domain filtering.

4.1.

Description of Computer-Simulation Studies

4.1.1.

Numerical phantom

We employed a numerical phantom consisting of $K = 90$ slow-time object frames, denoted by ${A_{k} (r)}_{k = 0}^{89}$ . The $k$ ’th frame of the phantom was defined as

Eq. (17)

A_{k} (r) = \sum_{n = 0}^{6} A_{n}^{s} (r) {[f_{n}]}_{k},

where each

A_{n}^{s} (r)

corresponded to a circular absorber for

n = 0, 1, \dots, 6

, and the vector

f_{n}

described the slow-time activity of the

n

‘th circular absorber. The vector

f_{n}

was of dimensions

90 \times 1

, with the

k

‘th element denoted by

{[f_{n}]}_{k}

. The pixel TAC of the numerical phantom was a superposition of the time activities of the seven absorbers as described in Eq. (17). We designed one element of

{f_{n}}_{n = 0}^{6}

as a linear combination of the other six elements. The slow-time sampling interval was set to be 1.6 s, which is consistent with an existing dynamic PACT imaging system.⁹ Finite-dimensional representations of the object frames

A_{k} (r)

were created according to Eq. (11) with

{r_{n}}_{n = 0}^{N - 1}

distributed on a uniform Cartesian grid of spacing 0.05 mm. These images were of dimension

440 \times 440

and, when lexiographically ordered, represented the columns of

A

. These images were displayed continuously and recorded as a video shown in the left panel of Video 1.

4.1.2.

Simulated measurement data

In the computer-simulation studies, we assumed a ring-shaped transducer array with a radius of 25 mm. The array contained $J = 512$ elements, which were uniformly distributed on the ring. Each data frame $g_{k}$ was analytically calculated by using Eq. (1) with a fast-time sampling rate of 40 MHz. We assumed idealized point-like transducers, i.e., $h^{e} (t) = δ (t)$ and $Ω_{j} = 0$ . For each transducer, $I = 650$ fast-time samples were calculated. Accordingly, the simulated noise-free data matrix $G$ had dimensions of $332,800 \times 90$ . The noise matrix $\underset{̲}{N}$ contained white Gaussian entries with variances specified as 20, 30, and 40% of ${∥ G ∥}_{F}^{2} / (I J K)$ . Summation of $\underset{̲}{N}$ and $G$ resulted in the noisy data matrix $\underset{̲}{G}$ according to Eq. (12).

4.1.3.

Image reconstruction

The STIR formula, i.e., Eq. (9), was validated by using simulated noise-free data. An FBFIR method was also implemented for comparison. For both reconstruction methods, the static image reconstruction operator $B$ was defined to be a discrete version of a 2-D filtered backprojection (FBP) formula⁵¹ whose description in continuous form is given by

Eq. (18)

{\hat{A}}_{k} (r) = \frac{C_{p}}{π c_{0} β R^{'}} \int_{Ω} d r^{'} \int_{0}^{\frac{2 R^{'}}{c_{0}}} d t \frac{\partial}{\partial t} t \times \frac{\partial}{\partial t} \log | t^{2} c_{0}^{2} - {| r^{'} - r |}^{2} | \int_{0}^{t} d τ g_{k} (r^{'}, τ),

where

R^{'} = 25

mm was the scanning radius,

c_{0} = 1.5 mm / μ s

,

Ω

was the measurement circle,

C_{p} / β

was set to be 1.0 in arbitrary units, and the measured pressure

g (r^{'}, t)

for the

k

’th frame was a function of location

r^{'} \in Ω

and fast-time

t

. Since we assumed idealized point-like transducers, samples of

g (r^{'}, t)

were equivalent to the simulated measurement data. Spatial samples of

A (r)

that resided on the same Cartesian grid employed in Eq. (11) were computed by discretizing Eq. (18).

The LRME-STIR method described in Sec. 3.3 was evaluated by using simulated noisy measurement data. We also implemented the FBFIR method followed by two types of image-domain filtering, namely, a Hann window low-pass filter and a PCA filter. These image-domain filters were applied to every pixel TAC of the images reconstructed by using the FBFIR method. Implementations of these filters are described in Appendices B and C. The FBFIR method followed by either the Hann window low-pass filter or the PCA filter will be referred to as the FBFIR-Hann method or the FBFIR-PCA method, respectively. Note that the performance of each method was affected by the choice of regularization parameter $β$ , the cutoff frequency $f_{c}$ , and the cutoff order $K_{c}$ , respectively. In order to compare optimal performances of these methods, we swept the regularization parameter values over wide ranges until the Euclidean distance between the reconstructed images and the phantom was minimized.

The accuracy of reconstructed images was quantified by the mean squared error (MSE)

Eq. (19)

MSE = \frac{1}{N K} {∥ A - \hat{A} ∥}_{F}^{2} .

4.2.

Experimental Studies

The LRME-STIR method was also evaluated with experimental data. The experimental data, also referred to as raw data, were acquired in a previous study of the wash-in process of Evans blue dye in a mouse brain.⁹ The small-animal photoacoustic imaging system contained a 512-element, i.e., $J = 512$ , full-ring transducer array and is described in detail elsewhere.¹²^,⁵² In order to form a data frame, each element acquired $I = 1300$ fast-time samples at the sampling rate of 40 MHz. Each transducer element in the array had a central frequency of 5 MHz and was focused in elevation. The combined foci of all elements created a 2-cm-diameter central imaging region with 100 μm in-plane resolution.⁸^,⁵³ Light illumination was provided by an optical parametric oscillator laser, tunable from 400 to 680 nm. The mouse was anesthetized using isoflurane and mounted in the system on a lab-made holder. During the experiment, Evans blue was slowly injected through the tail vein over 70 s, and the animal was simultaneously imaged at a frame rate of $0.625 s / frame$ . Accordingly, the slow-time sampling interval was 1.6 s and a total of $K = 90$ frames were acquired. Additional details regarding the experimental procedure can be found in Ref. 9.

For both the FBFIR and the STIR methods, the static image reconstruction operator $B$ was defined to be a discrete version of the simple backprojection operator.

Eq. (20)

A (r) \approx \int_{Ω} {d r^{'} p (r^{'}, t) |}_{t = \frac{| r - r^{'} |}{c_{0}}} .

This operator has been employed for image reconstruction in experimental PACT studies and provides only qualitative images.⁵⁴ A variety of FBP formulae have been reported.⁹^,¹⁰^,⁵² However, they are based on idealized imaging models that limit their ability to improve image quality. In our study, Eq. (20) was employed to reconstruct images that were sampled at locations on a uniform

360 \times 240

Cartesian grid with a grid spacing of 0.05 mm. Accordingly,

N = 360 \times 240

.

Because the imaged blood vessels in the mouse brain were shallow and highly absorbing, the raw data set possessed a high signal-to-noise ratio (SNR). In order to emulate imaging conditions or other hardware configurations that would produce lower SNR data, we added computer-simulated Gaussian white noise to the raw data, from which images were reconstructed by using the LRME-STIR, FBFIR-Hann, and FBFIR-PCA methods. The artificially synthesized data sets will be referred to as the high-noise data. We lacked knowledge of ground truth and, therefore, employed the images reconstructed from the original raw data by using the FBFIR method as a reference. The variance of the simulated noise was 300% of ${∥ G ∥}_{F}^{2} / (I J K)$ . The accuracy of reconstructed images was quantified by computing the MSE of the reconstructed images and the reference images.

5. Results

5.1.

Images Reconstructed from Simulated Noise-Free Data

The rank of the data matrix $G$ was six because ${f_{n}}_{n = 0}^{6}$ contained six linearly independent vectors.¹⁵ Therefore, the FBP algorithm was applied six times when implementing the STIR method. As shown in Video 1, the reconstructed images are nearly identical to the numerical phantom, with $MSE = 5.08 \times 10^{- 4}$ . Several isolated slow-time frames of the reconstructed images are displayed in Fig. 1. The pixel TACs of the reconstructed images align closely with those of the numerical phantom, as shown in Fig. 2. Though not shown here, the images reconstructed by using the FBFIR method are very similar to those reconstructed by the STIR method. However, the FBFIR method required applying the FBP algorithm 90 times. This idealized noise-free simulation corroborates the mathematical equivalence of the FBFIR and the STIR methods and demonstrates that the STIR method is a computationally efficient alternative to conventional FBFIR when the data matrix is of low rank.

Fig. 1

Slow-time frames of the images reconstructed by the use of spatiotemporal image reconstruction (STIR) from simulated noise-free data (Video 1, QuickTime 164 KB) [URL: http://dx.doi.org/10.1117/1.JBO.19.5.056007.1].

Fig. 2

Pixel time activity curves (TACs) of the phantom (solid) and of the images reconstructed by using STIR from noise-free data (circle), where blue, red, green, and black correspond to the pixels marked by A, B, C, and D in Fig. 1, respectively.

5.2.

Images Reconstructed from Simulated Noisy Data

As shown in Figs. 3 Fig. 4 Fig. 5 to 6, images reconstructed by using the LRME-STIR method are more accurate than those reconstructed via the FBFIR-Hann method and the FBFIR-PCA method. Figure 3 displays the 60’th slow-time frames of the numerical phantom and of the images reconstructed by all three methods. It is interesting to note that no obvious artifacts due to noise are visibly obvious in the isolated slow-time frames. Likely, the image noise was effectively mitigated during backprojection due to the densely sampled measurement data employed. However, the noise becomes obvious along the slow-time axis in the reconstructed images as shown in Figs. 4 and 5. Particularly, streak artifacts can be observed in the image produced by the FBFIR-PCA method shown in Fig. 4(c). Also, as expected, the Hann-window low-pass filter removes high-frequency components of both noise and signals, as shown in the TACs of pixel A in Fig. 5(a). In contrast, the LRME-STIR method mitigates noise with minimal degradation of signals as shown in Fig. 5(c). The improved accuracy from the LRME-STIR method can also be observed in Video 2. It is interesting to observe that the FBFIR-PCA, in general, performs worst among the three algorithms (see Fig. 5). This observation is likely due to the noise correlation introduced by the FBP algorithm, which degrades the performance of the PCA-based filtering.²⁷ The superior performance of the LRME-STIR method is consistent across varying noise levels, as shown in Fig. 6. In addition, the LRME-STIR method is computationally much more efficient, since the FBP algorithm needs to be applied only six times (for the 20% noisy data), as opposed to the FBFIR methods that require applying the FBP algorithm 90 times.

Fig. 3

Sixtieth slow-time frames of (a) the phantom, and the images reconstructed by using (b) the frame-by-frame image reconstruction (FBFIR)-Hann, (c) the FBFIR-PCA, and (d) the low-rank matrix estimation-based STIR (LRME-STIR) methods, respectively, from the noisy data contaminated with 20% Gaussian white noise (Video 2, QuickTime 258 KB) [URL: http://dx.doi.org/10.1117/1.JBO.19.5.056007.2]. The gray-scale window is $[- 0.1, 1.1]$ .

Fig. 4

Estimates of ${A_{k} (x, y = 0)}_{k = 0}^{89}$ by using (b) the FBFIR-Hann method, (c) the FBFIR-PCA method, and (d) the LRME-STIR method from the noisy data contaminated with 20% Gaussian white noise, respectively. The original phantom is plotted in panel (a). The gray-scale window is $[- 0.1, 1.1]$ .

Fig. 5

Pixel TACs of the phantom (solid) and of the images reconstructed by using (a) FBFIR-Hann, (b) FBFIR-PCA, and (c) LRME-STIR, respectively, from the noisy data contaminated with 20% Gaussian white noise. Blue, red, green, and black correspond to the pixels marked by A, B, C, and D in Fig. 1, respectively.

Fig. 6

Plots of the mean squared errors (MSE’s) of reconstructed images versus the noise levels.

5.3.

Images Reconstructed from Experimental Data

As shown in Fig. 7, an L-shaped singular value spectrum is observed for the measured data matrix $\underset{̲}{G}$ , revealing that its energy is concentrated in the first few singular components. This suggests that the noise-free data matrix $G$ is approximately of low rank. Similar to the computer-simulation studies, the isolated slow-time frames of the images reconstructed by using the LRME-STIR and the FBFIR methods are, in general, very similar (see Fig. 8). Here, $β$ was selected so that the resultant denoised data matrix $\hat{G}$ was of rank $\hat{R} = 4$ , corresponding to the elbow point in Fig. 7. This heuristic estimator of $β$ has been widely employed in similar studies.¹³^,¹⁷^,²⁷ Estimates of ${A_{k} (x = 0, y)}_{k = 0}^{89}$ reconstructed by using the FBFIR and the LRME-STIR methods are displayed in Fig. 9. Images reconstructed by using the LRME-STIR method appear to be less noisy than those reconstructed by using the FBFIR method (revealed clearly in Video 3). It is evident that the LRME-STIR method accurately captured the dye wash-in process, as shown in Fig. 10, which represents the dynamic process of interest. In addition, the computation required in the LRME-STIR method was only $4 / 90$ of that required in the FBFIR method.

Fig. 7

The singular value spectra of the measured raw data $\underset{̲}{P}$ .

Fig. 8

Slow-time frames of the images reconstructed from the experimental raw data by using the FBFIR (top row) and the LRME-STIR (bottom row) methods, respectively (Video 3, QuickTime 262 KB) [URL: http://dx.doi.org/10.1117/1.JBO.19.5.056007.3]. The gray-scale window is $[- 2.2, 5.3]$ .

Fig. 9

Estimates of ${A_{k} (x = 0, y)}_{k = 0}^{89}$ reconstructed by using (a) the FBFIR method and (b) the LRME-STIR method from the experimental raw data, respectively. The gray-scale window is $[- 2.2, 5.3]$ .

Fig. 10

Pixel TACs of the images reconstructed from experimental raw data by using the FBFIR (solid) and the LRME‐STIR (circles) methods, respectively. Blue, red, green, and black correspond to the pixels marked by A, B, C, and D in Fig. 8, respectively.

In the images reconstructed from the emulated high-noise-level data, the isolated slow-time frames in Fig. 11 reveal little difference among the FBFIR-Hann, the FBFIR-PCA, and the LRME-STIR methods, as expected. However, Fig. 12(b) and the slow-time TACs in Fig. 13 show severe distortions by using the FBFIR-Hann method. The $MSE$ of the images reconstructed by using the FBFIR-Hann, the FBFIR-PCA, and the LRME-STIR methods are $1.58 \times 10^{- 2}$ , $1.43 \times 10^{- 2}$ , and $1.42 \times 10^{- 2}$ , respectively. Unlike the results in the computer-simulation studies, the images reconstructed by using the FBFIR-PCA and the LRME-STIR methods appear to be similar. This observation is likely due to the fact that the simple backprojection does not correlate the noise. Nonetheless, the computation was reduced by a factor of $90 / 4$ when the LRME-STIR method was employed.

Fig. 11

Slow-time frames of the images reconstructed from the high-noise data set (Video 4, QuickTime 262 KB) [URL: http://dx.doi.org/10.1117/1.JBO.19.5.056007.4]. Columns from left to right correspond to the FBFIR‐Hann, the FBFIR‐PCA, and the LRME-STIR methods, respectively. The gray-scale window is $[- 2.2, 5.3]$ .

Fig. 12

Estimates of ${A_{k} (x = 0, y)}_{k = 0}^{89}$ reconstructed by using (a) the FBFIR method, (b) the FBFIR-Hann, (c) the FBFIR-PCA, and (d) the LRME-STIR methods from the experimental high-noise data set, respectively. The gray-scale window is $[- 2.2, 5.3]$ .

Fig. 13

Pixel TACs of the reference images (solid) and the images reconstructed by using (a) the FBFIR-Hann, (b) the FBFIR-PCA, and (c) the LRME-STIR methods, respectively, from the experimental data contaminated with computer-simulated noise. Blue, red, green, and black correspond to the pixels marked by A, B, C, and D in Fig. 8, respectively.

6. Conclusions and Discussion

In this study, we proposed and investigated an LRME-STIR method for dynamic PACT image reconstruction. The method employs a data denoising step followed by image reconstruction conducted in the singular system of the data matrix. In the ideal noise-free scenario, the method is mathematically equivalent to, but can be computationally far more efficient than, the conventional FBFIR method. In practice, compared with the conventional FBFIR method, the LRME-STIR method can improve not only the computational efficiency but also the quantitative accuracy due to the low-rank regularization employed in the data denoising step. Although our studies employed a 2-D imaging geometry, our conclusions are equally applicable to the general 3-D case.

The LRME-STIR method exploits the fact that, for many dynamic PACT applications, the data matrix can be approximately described as a low-rank matrix whose rank is typically much smaller than the number of slow-time frames. The low-rank assumption has been explored for other dynamic imaging applications.¹⁵^,¹⁷ However, if a small moving structure is of interest, the data matrix may not be effectively low-rank. Instead, the data matrix can be decomposed into a low-rank component and a sparse component by using robust PCA.¹⁸^,²⁸^,²⁹ Application of robust PCA to dynamic PACT remains an interesting topic for future studies.

Under certain conditions, the implementation of the LRME-STIR method is identical to that of conventional data-domain KL filtering.¹³^,²⁵^,²⁷ More specifically, this is true when (1) the sample covariance matrix, i.e., ${\underset{̲}{G}}^{†} \underset{̲}{G}$ , is employed to estimate the covariance matrix for the data-domain KL filtering and (2) the low-rank assumption is formulated as the optimization problem in Eq. (13) for the LRME-STIR. Without satisfying these conditions, the implementation, as well as the performance, of the LRME-STIR method will be different from that of the data-domain KL filtering. It is well-known that a variety of other algorithms can be employed to estimate the covariance matrix.⁵⁵ Also, based on the same low-rank matrix assumption, the data matrix denoising problem can be formulated in different ways.⁴⁸ For example, a commonly employed nuclear norm minimization formulation is expressed as¹⁸^,²⁹^,⁵⁶

Eq. (21)

\hat{G} = \arg \min_{G} \frac{1}{2} {∥ \underset{̲}{G} - G ∥}_{F}^{2} + β {∥ G ∥}_{*},

where the nuclear norm of

G

, denoted by

{∥ G ∥}_{*}

, functions as a convex relaxation of the

Rank (G)

in Eq. (13). It is currently unclear which formulation is optimal for PACT applications. However, it is clear that the rationales behind the LRME-STIR method and the data-domain KL filtering are substantially different.

The computational advantage of the LRME-STIR method will be more significant if advanced iterative image reconstruction algorithms are employed to implement the linear image reconstruction operator $B$ in Eq. (9).²⁰ Combining the LRME-STIR method with iterative image reconstruction algorithms remains as an important topic for future studies.

Appendices

Appendix A:

Interpretation of the Row Vectors of $\hat{G}$ as Linear Combinations of the Row Vectors of $\underset{̲}{G}$

Assume $\underset{̲}{G}$ is of full rank, i.e, $Rank (\underset{̲}{G}) = K$ , with $K \leq (I J)$ . Each right singular vector of $\underset{̲}{G}$ can be written as a linear combination of the transpose of the row vectors of $\underset{̲}{G}$ , i.e.,

Eq. (22)

u_{k} = \sum_{m = 0}^{I J - 1} c_{k, m} q_{m}, for k = 0, 1, \dots, K - 1,

where

q_{m}

denotes the transpose of the

m

’th row vector of

\underset{̲}{G}

and

c_{k, m}

denotes the corresponding coefficient. On substitution of

u_{k}

from Eq. (22) into Eq. (8), one obtains

Eq. (23)

{\underset{̲}{v}}_{k} = \frac{1}{{\underset{̲}{μ}}_{k}} [\begin{matrix} \sum_{m = 0}^{I J - 1} c_{k, m} q_{0}^{†} q_{m} \\ \sum_{m = 0}^{I J - 1} c_{k, m} q_{1}^{†} q_{m} \\ ⋮ \\ \sum_{m = 0}^{I J - 1} c_{k, m} q_{I J - 1}^{†} q_{m} \end{matrix}], for k = 0, 1, \dots, K - 1 .

Substituting Eq. (23) into Eq. (14) results in

Eq. (24)

\hat{G} = [\begin{matrix} \sum_{m^{'} = 0}^{I J - 1} (\sum_{k = 0}^{\hat{R} - 1} \sum_{m = 0}^{I J - 1} c_{k, m} c_{k, m^{'}} q_{0}^{†} q_{m}) q_{m^{'}}^{†} \\ \sum_{m^{'} = 0}^{I J - 1} (\sum_{k = 0}^{\hat{R} - 1} \sum_{m = 0}^{I J - 1} c_{k, m} c_{k, m^{'}} q_{1}^{†} q_{m}) q_{m^{'}}^{†} \\ ⋮ \\ \sum_{m^{'} = 0}^{I J - 1} (\sum_{k = 0}^{\hat{R} - 1} \sum_{m = 0}^{I J - 1} c_{k, m} c_{k, m^{'}} q_{I J - 1}^{†} q_{m}) q_{m^{'}}^{†} \end{matrix}] .

Note that the quantities in all parentheses in Eq. (24) are scalars. Therefore, each row vector of

\hat{G}

can be interpreted as a linear combination of the row vectors of

\underset{̲}{G}

.

Appendix B:

Implementation of Image-Domain Hann Window Filter

The image-domain Hann window low-pass filter is defined as

Eq. (25)

W (f) = \frac{1}{2} [1 - \cos (π \frac{f_{c} - f}{f_{c}})],

where

f

is the temporal frequency coordinate corresponding to the slow-time coordinate and

f_{c}

is a cutoff frequency functioning as a regularization parameter. Since the slow-time sampling interval was 1.6 s, we swept the value of

f_{c}

from 0.02 to 0.3 Hz, with a step size of 0.01 Hz.

Appendix C:

Implementation of Image-Domain Principal Component Analysis Filter

Implementation of the principal component analysis filter is summarized as follows. First, we estimate the covariance matrix of the random slow-time activity vector by⁵⁷

Eq. (26)

S = \frac{1}{N - 1} {(\underset{̲}{\hat{A}} - \bar{A})}^{†} (\underset{̲}{\hat{A}} - \bar{A}),

where

S

is the sample covariance matrix of dimension

K \times K

,

N

is the number of pixels in each slow-time frame, and

\underset{̲}{\hat{A}}

denotes the image matrix reconstructed by using the FBFIR method from noisy data. Note that Eq. (26) treats each row of

\underset{̲}{\hat{A}}

as a realization of the random slow-time activity vector. The

N

realizations become zero-mean after subtraction of the mean image matrix

\bar{A}

, which is defined as

Eq. (27)

\bar{A} = [\begin{matrix} {\bar{A}}_{0} & {\bar{A}}_{1} & \dots & {\bar{A}}_{K - 1} \\ {\bar{A}}_{0} & {\bar{A}}_{1} & \dots & {\bar{A}}_{K - 1} \\ ⋮ & ⋮ & \dots & ⋮ \\ {\bar{A}}_{0} & {\bar{A}}_{1} & \dots & {\bar{A}}_{K - 1} \end{matrix}],

where

Eq. (28)

{\bar{A}}_{k} = \frac{1}{N} \sum_{n = 0}^{N - 1} {[\underset{̲}{\hat{A}}]}_{n, k}, for k = 0, 1, \dots, K - 1 .

Accordingly, both

\underset{̲}{\hat{A}}

and

\bar{A}

are of dimensions

N \times K

. Second, we calculate the eigenvectors of

S

. The eigenvectors, denoted by

{η_{k}}_{k = 0}^{K - 1}

, are sorted so that the corresponding eigenvalues decrease as

k

increases. Finally, we filter the high-order components.

Eq. (29)

\hat{A} = \sum_{k = 0}^{K_{c} - 1} (\underset{̲}{\hat{A}} - \bar{A}) η_{k} η_{k}^{†} + \bar{A},

where

K_{c}

is the cutoff order functioning as a regularization parameter.

Acknowledgments

This work was supported in part by NIH awards EB016963, EB010049, and CA1744601. K.W. would also like to thank Dr. Robert W. Schoonover for useful discussions.

References

1.

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

2.

S. Yanget al., “Functional imaging of cerebrovascular activities in small animals using high-resolution photoacoustic tomography,” Med. Phys., 34 (8), 3294 –3301 (2007). http://dx.doi.org/10.1118/1.2757088 MPHYA6 0094-2405 Google Scholar

3.

H.-P. Brechtet al., “Whole-body three-dimensional optoacoustic tomography system for small animals,” J. Biomed. Opt., 14 (6), 064007 (2009). http://dx.doi.org/10.1117/1.3259361 JBOPFO 1083-3668 Google Scholar

4.

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

5.

A. A. OraevskyA. A. Karabutov, “Optoacoustic tomography,” Biomedical Photonics Handbook, 34/1 –34/34 CRC Press LLC, Boca Raton, FL (2003). Google Scholar

6.

R. KrugerD. ReineckeG. Kruger, “Thermoacoustic computed tomography—technical considerations,” Med. Phys., 26 (9), 1832 –1837 (1999). http://dx.doi.org/10.1118/1.598688 MPHYA6 0094-2405 Google Scholar

7.

K. WangM. A. Anastasio, “Photoacoustic and thermoacoustic tomography: image formation principles,” Handbook of Mathematical Methods in Imaging, 781 –815 Springer, New York, NY (2011). Google Scholar

8.

J. Gamelinet al., “A real-time photoacoustic tomography system for small animals,” Opt. Express, 17 (13), 10489 –10498 (2009). http://dx.doi.org/10.1364/OE.17.010489 OPEXFF 1094-4087 Google Scholar

9.

C. Liet al., “Real-time photoacoustic tomography of cortical hemodynamics in small animals,” J. Biomed. Opt., 15 (1), 010509 (2010). http://dx.doi.org/10.1117/1.3302807 JBOPFO 1083-3668 Google Scholar

10.

A. Buehleret al., “Video rate optoacoustic tomography of mouse kidney perfusion,” Opt. Lett., 35 (14), 2475 –2477 (2010). http://dx.doi.org/10.1364/OL.35.002475 OPLEDP 0146-9592 Google Scholar

11.

L. Xianget al., “4-D photoacoustic tomography,” Sci. Rep., 3 (1113), 1 –8 (2013). http://dx.doi.org/10.1038/srep01113 SRCEC3 2045-2322 Google Scholar

12.

M. R. Chatniet al., “Tumor glucose metabolism imaged in vivo in small animals with whole-body photoacoustic computed tomography,” J. Biomed. Opt., 17 (7), 076012 (2012). http://dx.doi.org/10.1117/1.JBO.17.7.076012 JBOPFO 1083-3668 Google Scholar

13.

M. WernickE. InfusinoM. Milosevic, “Fast spatio-temporal image reconstruction for dynamic PET,” IEEE Trans. Med. Imaging, 18 (3), 185 –195 (1999). http://dx.doi.org/10.1109/42.764885 ITMID4 0278-0062 Google Scholar

14.

D. S. LalushB. M. W. Tsui, “Block-iterative techniques for fast 4D reconstruction using a priori motion models in gated cardiac SPECT,” Phys. Med. Biol., 43 (4), 875 –886 (1998). http://dx.doi.org/10.1088/0031-9155/43/4/015 PHMBA7 0031-9155 Google Scholar

15.

J. HaldarZ.-P. Liang, “Spatiotemporal imaging with partially separable functions: a matrix recovery approach,” in IEEE Int. Symp. on Biomedical Imaging: From Nano to Macro, 716 –719 (2010). Google Scholar

16.

H. Pedersenet al., “k-t PCA: temporally constrained k-t BLAST reconstruction using principal component analysis,” Magn. Reson. Med., 62 (3), 706 –716 (2009). http://dx.doi.org/10.1002/mrm.v62:3 MRMEEN 0740-3194 Google Scholar

17.

C. Brinegaret al., “Improving temporal resolution of pulmonary perfusion imaging in rats using the partially separable functions model,” Magn. Reson. Med., 64 (4), 1162 –1170 (2010). http://dx.doi.org/10.1002/mrm.v64:4 MRMEEN 0740-3194 Google Scholar

18.

S. Lingalaet al., “Accelerated dynamic MRI exploiting sparsity and low-rank structure: k-t SLR,” IEEE Trans. Med. Imaging, 30 (5), 1042 –1054 (2011). http://dx.doi.org/10.1109/TMI.2010.2100850 ITMID4 0278-0062 Google Scholar

19.

J. TsaoP. BoesigerK. P. Pruessmann, “k-t blast and k-t sense: dynamic MRI with high frame rate exploiting spatiotemporal correlations,” Magn. Reson. Med., 50 (5), 1031 –1042 (2003). http://dx.doi.org/10.1002/(ISSN)1522-2594 MRMEEN 0740-3194 Google Scholar

20.

K. Wanget 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

21.

X. Dean-Benet 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

22.

J. K. GamelinA. AguirreQ. Zhu, “Fast, limited-data photoacoustic imaging for multiplexed systems using a frequency-domain estimation technique,” Med. Phys., 38 (3), 1503 –1518 (2011). http://dx.doi.org/10.1118/1.3533669 MPHYA6 0094-2405 Google Scholar

23.

J. Glatzet al., “Blind source unmixing in multi-spectral optoacoustic tomography,” Opt. Express, 19 (4), 3175 –3184 (2011). http://dx.doi.org/10.1364/OE.19.003175 OPEXFF 1094-4087 Google Scholar

24.

H. BarrettK. Myers, Foundations of Image Science, John Wiley and Sons, New York, NY (2004). Google Scholar

25.

M. N. WernickJ. N. Aarsvold, Emission Tomography, the Fundamentals of PET and SPECT, Elsevier Academic Press, San Diego, California (2004). Google Scholar

26.

M. Jinet al., “4D reconstruction for low-dose cardiac gated SPECT,” Med. Phys., 40 (2), 022501 (2013). http://dx.doi.org/10.1118/1.4773868 MPHYA6 0094-2405 Google Scholar

27.

Z. Chenet al., “Temporal processing of dynamic positron emission tomography via principal component analysis in the sinogram domain,” IEEE Trans. Nucl. Sci., 51 (5), 2612 –2619 (2004). http://dx.doi.org/10.1109/TNS.2004.834816 IETNAE 0018-9499 Google Scholar

28.

E. J. Candèset al., “Robust principal component analysis,” J. ACM, 58 (3), 11:1 –11:37 (2011). http://dx.doi.org/10.1145/1970392 JOACF6 0004-5411 Google Scholar

29.

H. Gaoet al., “Robust principal component analysis-based four-dimensional computed tomography,” Phys. Med. Biol., 56 (11), 3181 (2011). http://dx.doi.org/10.1088/0031-9155/56/11/002 PHMBA7 0031-9155 Google Scholar

30.

W. Qiet al., “A quantitative study of motion estimation methods on 4D cardiac gated SPECT reconstruction,” Med. Phys., 39 (8), 5182 –5193 (2012). http://dx.doi.org/10.1118/1.4738377 MPHYA6 0094-2405 Google Scholar

31.

L. V. WangH.-I. Wu, Biomedical Optics, Principles and Imaging, Wiley, Hoboken, N.J. (2007). Google Scholar

32.

K. Wanget al., “An imaging model incorporating ultrasonic transducer properties for three-dimensional optoacoustic tomography,” IEEE Trans. Med. Imaging, 30 (2), 203 –214 (2011). http://dx.doi.org/10.1109/TMI.2010.2072514 ITMID4 0278-0062 Google Scholar

33.

A. Conjusteauet al., “Measurement of the spectral directivity of optoacoustic and ultrasonic transducers with a laser ultrasonic source,” Rev. Sci. Instrum., 80 (9), 093708 (2009). http://dx.doi.org/10.1063/1.3227836 RSINAK 0034-6748 Google Scholar

34.

P. C. Beardet al., “Quantitative photoacoustic imaging: measurement of absolute chromophore concentrations for physiological and molecular imaging,” Photoacoustic Imaging and Spectroscopy, 121 –143 CRC Press, Boca Raton, FL (2009). Google Scholar

35.

B. Coxet al., “Quantitative spectroscopic photoacoustic imaging: a review,” J. Biomed. Opt., 17 (6), 061202 (2012). http://dx.doi.org/10.1117/1.JBO.17.6.061202 JBOPFO 1083-3668 Google Scholar

36.

A. RosenthalV. NtziachristosD. Razansky, “Model-based optoacoustic inversion with arbitrary-shape detectors,” Med. Phys., 38 (7), 4285 –4295 (2011). http://dx.doi.org/10.1118/1.3589141 MPHYA6 0094-2405 Google Scholar

37.

D. FinchS. PatchRakesh, “Determining a function from its mean values over a family of spheres,” SIAM J. Math. Anal., 35 (5), 1213 –1240 (2004). http://dx.doi.org/10.1137/S0036141002417814 SJMAAH 0036-1410 Google Scholar

38.

M. XuL. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E, 71 (71), 016706 (2005). http://dx.doi.org/10.1103/PhysRevE.71.016706 PLEEE8 1063-651X Google Scholar

39.

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 INPEEY 0266-5611 Google Scholar

40.

P. KuchmentL. Kunyansky, “Mathematics of photoacoustic and thermoacoustic tomography,” Handbook of Mathematical Methods in Imaging, 817 –865 Springer, New York (2011). Google Scholar

41.

K. Wanget al., “Sparsity regularized data-space restoration in optoacoustic tomography,” Proc. SPIE, 8223 822322 (2012). http://dx.doi.org/10.1117/12.909690 PSISDG 0277-786X Google Scholar

42.

Z.-P. Liang, “Spatiotemporal imaging with partially separable functions,” in 4th IEEE Int. Symp. on Biomedical Imaging: From Nano to Macro, 2007, 988 –991 (2007). Google Scholar

43.

D. Donoho, “De-noising by soft-thresholding,” IEEE Trans. Inf. Theory, 41 (3), 613 –627 (1995). http://dx.doi.org/10.1109/18.382009 IETTAW 0018-9448 Google Scholar

44.

A. BuadesB. CollJ. Morel, “A review of image denoising algorithms, with a new one,” Multiscale Model. Simul., 4 (2), 490 –530 (2005). http://dx.doi.org/10.1137/040616024 MMSUBT 1540-3459 Google Scholar

45.

J.-L. StarckE. CandesD. Donoho, “The curvelet transform for image denoising,” IEEE Trans. Image Process., 11 (6), 670 –684 (2002). http://dx.doi.org/10.1109/TIP.2002.1014998 IIPRE4 1057-7149 Google Scholar

46.

E. CandèsB. Recht, “Exact matrix completion via convex optimization,” Found. Comput. Math., 9 (6), 717 –772 (2009). http://dx.doi.org/10.1007/s10208-009-9045-5 1615-3375 Google Scholar

47.

D. L. DonohoM. Gavish, “The optimal hard threshold for singular values is 4/√(3),” (2014) https://statistics.stanford.edu/sites/default/files/2013-04.pdf March ). 2014). Google Scholar

48.

E. Candeset al., “Unbiased risk estimates for singular value thresholding and spectral estimators,” IEEE Trans. Signal Process., 61 (19), 4643 –4657 (2013). http://dx.doi.org/10.1109/TSP.2013.2270464 Google Scholar

49.

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

50.

B. Smithet al., Matrix Eigensystem Routines. EISPACK Guide, Volume 6 of Lecture Notes in Computer Science, Springer-Verlag, New York (1976). Google Scholar

51.

D. FinchM. HaltmeierRakesh, “Inversion of spherical means and the wave equation in even dimensions,” SIAM J. Appl. Math., 68 (2), 392 –412 (2007). http://dx.doi.org/10.1137/070682137 SMJMAP 0036-1399 Google Scholar

52.

J. Xiaet al., “Whole-body ring-shaped confocal photoacoustic computed tomography of small animals in vivo,” J. Biomed. Opt., 17 (5), 050506 (2012). http://dx.doi.org/10.1117/1.JBO.17.5.050506 JBOPFO 1083-3668 Google Scholar

53.

J. Xiaet al., “Three-dimensional photoacoustic tomography based on the focal-line concept,” J. Biomed. Opt., 16 (9), 090505 (2011). http://dx.doi.org/10.1117/1.3625576 JBOPFO 1083-3668 Google Scholar

54.

C. Huanget al., “Aberration correction for transcranial photoacoustic tomography of primates employing adjunct image data,” J. Biomed. Opt., 17 (6), 066016 (2012). http://dx.doi.org/10.1117/1.JBO.17.6.066016 JBOPFO 1083-3668 Google Scholar

55.

Statistical Digital Signal Processing and Modeling, 1st ed.John Wiley & Sons Inc., New York, NY (1996). Google Scholar

56.

A. MajumdarR. K. Ward, “Causal dynamic MRI reconstruction via nuclear norm minimization,” Magn. Reson. Imaging, 30 (10), 1483 –1494 (2012). http://dx.doi.org/10.1016/j.mri.2012.04.012 MRIMDQ 0730-725X Google Scholar

57.

R. Vershynin, “How close is the sample covariance matrix to the actual covariance matrix?,” J. Theor. Probab., 25 (3), 655 –686 (2012). http://dx.doi.org/10.1007/s10959-010-0338-z JTPREO 0894-9840 Google Scholar

Biography

Kun Wang received his PhD degree in biomedical engineering from Washington University in St. Louis (WUSTL) in 2012. Currently he is a research scientist in the group of Prof. Mark A. Anastasio in WUSTL, where his research interests include tomographic reconstruction in various imaging modalities. He has firstly demonstrated the feasibility of iterative image reconstruction algorithm for PACT in practice. He has published 12 peer-reviewed articles and two book chapters.

Jun Xia earned his PhD degree at the University of Toronto and is currently a postdoctoral fellow at Washington University in St. Louis, under the mentorship of Dr. Lihong V. Wang. His research interests are the development of novel biomedical imaging techniques including photoacoustic imaging and ultrasonic imaging. He has published more than 20 peer-reviewed journal articles in photoacoustic and photothermal research.

Changhui Li is an associate professor of the Department of Biomedical at Peking University. He received his BS and MS degrees at Peking University in 1997 and 2000, respectively. He got his PhD in physics in 2006 from the Department of Physics at Texas A&M University. His research focuses on biomedical optics, functional and molecular imaging, and computational electromagnetics.

Lihong V. Wang earned his PhD degree at Rice University, Houston, Texas. He currently holds the Gene K. Beare Distinguished Professorship of Biomedical Engineering at Washington University in St. Louis. He has published 342 peer-reviewed journal articles and delivered 370 keynote, plenary, or invited talks. His Google Scholar h-index and citations have reached 81 and over 26,000, respectively.

Mark A. Anastasio earned his PhD degree at the University of Chicago in 2001. He is currently a professor of biomedical engineering at Washington University in St. Louis (WUSTL). His research interests include tomographic image reconstruction, imaging physics, and the development of novel computed biomedical imaging systems. He is an internationally recognized expert in the fields of diffraction tomography, x-ray phase-contrast x-ray imaging, and photoacoustic computed tomography.

Citation Download Citation

Kun Wang, Jun Xia, Changhui Li, Lihong V. Wang, and Mark A. Anastasio "Fast spatiotemporal image reconstruction based on low-rank matrix estimation for dynamic photoacoustic computed tomography," Journal of Biomedical Optics 19(5), 056007 (9 May 2014). https://doi.org/10.1117/1.JBO.19.5.056007

Published: 9 May 2014

Access the abstract

JOURNAL ARTICLE
12 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 26 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Photoacoustic tomography

Image restoration

Reconstruction algorithms

Principal component analysis

Image filtering

Transducers

Biological research

1.

Introduction

2.

Background

2.1.

Static Image Reconstruction in Conventional PACT

Eq. (1)

Eq. (2)

2.2.

Dynamic PACT and Frame-by-Frame Image Reconstruction

Eq. (3)

Eq. (4)

Eq. (5)

3.

Low-Rank Matrix Estimation-Based STIR for Dynamic PACT

3.1.

Singular Value Decomposition-Based STIR

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

Eq. (11)

3.2.

Low-Rank Regularized Data Matrix Denoising

Eq. (12)

Eq. (13)

Eq. (14)

3.3.

LRME-Based STIR Method

Eq. (15)

Eq. (16)

4.

Numerical Studies

4.1.

Description of Computer-Simulation Studies

4.1.1.

Numerical phantom

Eq. (17)

4.1.2.

Simulated measurement data

4.1.3.

Image reconstruction

Eq. (18)

Eq. (19)

4.2.

Experimental Studies

Eq. (20)

5.

Results

5.1.

Images Reconstructed from Simulated Noise-Free Data

Fig. 1

Fig. 2

5.2.

Images Reconstructed from Simulated Noisy Data

Fig. 3

Fig. 4

Fig. 5

Fig. 6

5.3.

Images Reconstructed from Experimental Data

Fig. 7

Fig. 8

Fig. 9

Fig. 10

Fig. 11

Fig. 12

Fig. 13

6.

Conclusions and Discussion

Eq. (21)

Appendices

Appendix A:

Interpretation of the Row Vectors of G^ as Linear Combinations of the Row Vectors of G̲

Eq. (22)

Eq. (23)

Eq. (24)

Appendix B:

Implementation of Image-Domain Hann Window Filter

Interpretation of the Row Vectors of $\hat{G}$ as Linear Combinations of the Row Vectors of $\underset{̲}{G}$