Photoacoustic image reconstruction in an attenuating medium using singular-value decomposition

Dimple Modgil; Patrick J. La Rivière; Bradley E. Treeby

doi:10.1117/1.JBO.17.6.061204

7 May 2012 Photoacoustic image reconstruction in an attenuating medium using singular-value decomposition

Dimple Modgil, Patrick J. La Rivière, Bradley E. Treeby

Author Affiliations +

Journal of Biomedical Optics, Vol. 17, Issue 6, 061204 (May 2012). https://doi.org/10.1117/1.JBO.17.6.061204

Abstract

Attenuation effects can be significant in photoacoustic tomography since the generated pressure signals are broadband, and ignoring them may lead to image artifacts and blurring. La Rivière et al. [Opt. Lett.31(6), pp. 781-783, (2006)] had previously derived a method for modeling the attenuation effect and correcting for it in the image reconstruction. This was done by relating the ideal, unattenuated pressure signals to the attenuated pressure signals via an integral operator. We derive an integral operator relating the attenuated pressure signals to the absorbed optical energy for a planar measurement geometry. The matrix operator relating the two quantities is a function of the temporal frequency, attenuation coefficient and the two-dimensional spatial frequency. We perform singular-value decomposition (SVD) of this integral operator to study the problem further. We find that the smallest singular values correspond to wavelet-like eigenvectors in which most of the energy is concentrated at times corresponding to greater depths in tissue. This allows us to characterize the ill-posedness of recovering the absorbed optical energy distribution at different depths in an attenuating medium. This integral equation can be inverted using standard SVD methods, and the initial pressure distribution can be recovered. We conduct simulations and derive an algorithm for image reconstruction using SVD for a planar measurement geometry. We also study the noise and resolution properties of this image-reconstruction method.

1. Introduction

The majority of the image-reconstruction algorithms developed for photoacoustic tomography (PAT) or optoacoustic tomography (OAT) have assumed a lossless acoustic medium. The effect of frequency-dependent attenuation on acoustic waves can be significant, since the ultrasonic waves generated by pulsed lasers can be extremely broadband.¹ Reconstructed images may exhibit distortion and artifacts if these effects are not taken into account. Several researchers have looked at the effect of acoustic attenuation in OAT. Previous work on dispersive acoustic media done by La Rivière et al.² focused on incorporating the frequency-dependent attenuation effects into the OAT model. They related the attenuated pressure in a lossy medium to the ideal pressure in a lossless medium via an integral operator. More recently, several groups¹^,³^–4.⁵ have used the time-reversal (TR) approach to address attenuation in PAT. This method is based on the use of an acoustic forward model in which the absorption operators are reversed in sign. Roitner and Burgholzer⁶ have looked at a complex frequency-approach regularization method to correct for acoustic attenuation in PAT.

In this work, we will use an approach similar to that by Markel et al.,⁷^,⁸ in optical-diffusion tomography, to derive an inversion formula for the absorbed optical absorption energy using singular-value decomposition (SVD). This formula is applicable in a planar detector geometry. It provides insight into the conditioning of the inverse problem and offers a promising method for image reconstruction in an attenuating medium. This expression directly relates the measured attenuated pressure to the absorbed optical energy. We also study the noise and resolution properties of this SVD-based algorithm.

2. Correction for Ultrasonic Attenuation in OAT Using SVD

2.1.

Methods

In an attenuating medium, the optoacoustic wave equation includes a loss term²^,⁹^,¹⁰ and is given by:

Eq. (1)

\nabla^{2} p (r, t) - \frac{\partial^{2} p (r, t)}{c_{0}^{2} \partial t^{2}} + L (t) * p (r, t) = - \frac{β}{C_{p}} \frac{\partial}{\partial t} H (r, t),

where

L (t) = 1 / (2 π) \int d ω {k (ω)^{2} - (ω^{2} / c_{0}^{2})} e^{- i ω t},

and

Eq. (2)

k (ω) = \frac{ω}{c (ω)} + i α (ω),

where

ω

is the temporal frequency,

α (ω)

is the linear attenuation coefficient,

k (ω)

is the complex wave number,

p (r, t)

is the measured pressure at point

r

and time

t

,

β

is the coefficient of volume expansion,

C_{p}

is the specific heat and

c_{0}

is the speed of sound.

H (r, t)

is the heating function that denotes the energy deposited per unit time per unit volume by the illuminating optical pulse and is given by:

Eq. (3)

H (r, t) = A (r) I (t),

where

A (r)

is the absorbed optical energy and

I (t)

is the temporal profile of the optical pulse.

For soft biological tissue, the attenuation coefficient has a power-law dependence on the frequency given by:¹¹

Eq. (4)

α (ω) = α_{0} | ω |^{n} .

In tissues, typical values for ultrasonic attenuation are:

n \approx 1

and

α_{0} \approx 0.1 dB {MHz}^{- 1} {cm}^{- 1}

(Ref. 11).

Using the Kramers-Kronig relations, for power-law absorption with $0 < n < 3$ and $n \neq 1$ , one obtains the dependence of phase velocity $c (ω)$ on frequency due to frequency-dependent attenuation as:⁹

Eq. (5)

\frac{1}{c (ω)} = \frac{1}{c_{0}} + α_{0} \tan (π n / 2) (ω^{n - 1} - ω_{0}^{n - 1}),

where

c_{0}

is the speed of sound at a reference frequency of

ω_{0}

.

For $n = 1$ , the dependence of phase velocity $c (ω)$ on frequency due to frequency-dependent attenuation is given by:

Eq. (6)

\frac{1}{c (ω)} = \frac{1}{c_{0}} - \frac{2}{π} α_{0} l n | \frac{ω}{ω_{0}} | .

So for an attenuating medium, Eq. (1) becomes:

Eq. (7)

▿^{2} \tilde{p} (r, ω) + k^{2} (ω) \tilde{p} (r, ω) = - \frac{i ω β}{C_{p}} A (r) \tilde{I} (ω),

where

\tilde{p} (r, ω)

is the temporal Fourier transform of measured pressure and

\tilde{I} (ω)

is the temporal Fourier transform of the illumination pulse.

This equation can be solved using standard Green’s function methods as:¹²

Eq. (8)

\tilde{p} (r, ω) = - i η ω \tilde{I} (ω) ∭ d^{3} r^{'} A (r^{'}) G (r - r^{'}),

where

η \equiv β / C_{p}

and

G (r - r^{'})

is the Green’s function.

2.1.1.

Angular spectrum expansion of the measured pressure signals

The angular spectrum expansion of the Green’s function is given by:¹³

Eq. (9)

G (r) = \frac{\exp (i k^{2} (ω) r)}{4 π r} = \frac{1}{8 π^{2}} \iint \frac{1}{{[α_{x}^{2} + α_{y}^{2} - k^{2} (ω)]}^{1 / 2}} \times \exp [- | z | {[α_{x}^{2} + α_{y}^{2} - k^{2} (ω)]}^{1 / 2} + i x α_{x} + i y α_{y}] d α_{x} d α_{y} .

Substitute this in Eq. (8) and consider the pressure measurements on the plane

z = 0

and we get:

Eq. (10)

\tilde{p} (x, y, 0, ω) = \frac{- i η ω \tilde{I} (ω)}{8 π^{2}} \iint d^{3} r^{'} A (x^{'}, y^{'}, z^{'}) \iint d α_{x} d α_{y} \frac{1}{{[α_{x}^{2} + α_{y}^{2} - k^{2} (ω)]}^{1 / 2}} \times \exp [- | z^{'} | {{[α}_{x}^{2} + α_{y}^{2} - k^{2} (ω)]}^{1 / 2} + i (x - x^{'}) α_{x} + i (y - y^{'}) α_{y}] .

We assume that the photoacoustic object lies in the plane $z^{'} \geq 0$ . On taking the Fourier transform on both sides and reducing the resulting expression, one obtains the angular spectrum expansion of measured photoacoustic pressure as:

Eq. (11)

\tilde{p} (k_{x}, k_{y}, ω) = \frac{- i η ω \tilde{I} (ω)}{2 {[k_{x}^{2} + k_{y}^{2} - k^{2} (ω)]}^{1 / 2}} \int d z^{'} \tilde{A} (k_{x}, k_{y}, z^{'}) \exp [- z^{'} {[k_{x}^{2} + k_{y}^{2} - k^{2} (ω)]}^{1 / 2}] .

2.1.2.

SVD of integral operator relating pressure to optical-absorption coefficient

We follow Markel’s⁷^,⁸ approach to optical-diffusion tomography to obtain an integral operator relating the measured attenuated pressure to the absorbed-optical-energy density.

Define a two-dimensional (2-D) spatial wave vector $q$ as:

Eq. (12)

q \equiv (k_{x}, k_{y}, 0) .

We then discretize the temporal frequency,

ω

as

ω_{n}

. For notational simplicity, we use

z

instead of

z^{'}

in the following equations. One can then write the pressure as:

Eq. (13)

\tilde{p} (q, ω_{n}) = \frac{- i η ω_{n} \tilde{I} (ω_{n})}{2 {[q^{2} - k^{2} (ω_{n})]}^{1 / 2}} \int d z \tilde{A} (q, z) \exp [- z {[q^{2} - k^{2} (ω_{n})]}^{1 / 2}],

where

q \equiv {(k_{x}^{2} + k_{y}^{2})}^{1 / 2}

.

Define $P_{n} (q) \equiv \tilde{p} (q, ω_{n})$ , $k_{n} \equiv k (ω_{n})$ . Thus, one can obtain an integral equation:

Eq. (14)

P_{n} (q) = S_{n} (q) \int d z \tilde{A} (q, z) K_{n} (q, z),

where

Eq. (15)

K_{n} (q, z) \equiv \exp [- z {(q^{2} - k_{n}^{2})}^{1 / 2}],

Eq. (16)

S_{n} (q) \equiv \frac{- i η ω_{n} \tilde{I} (ω_{n})}{2 {(q^{2} - k_{n}^{2})}^{1 / 2}},

This equation can be inverted as:⁸

Eq. (17)

\tilde{A} (q, z) = \sum_{m, n} K_{m}^{*} (q, z) M_{m n}^{- 1} (q) \frac{P_{n} (q)}{S_{n} (q)},

where

Eq. (18)

M_{m n} (q) \equiv \int d z K_{m} (q, z) K_{n}^{*} (q, z),

the operator

K^{*}

denotes the complex conjugate of matrix

K

, and

M^{- 1}

is the inverse of matrix

M

.

We can obtain a simpler expression for $M_{m n} (q)$ by performing the integral over $z$ . Equation (18) can be written as:

Eq. (19)

M_{m n} (q) = \int_{0}^{\infty} d z \exp [- z {Q_{m} (q) + Q_{n}^{*} (q)}],

where

Q_{n}^{*} (q) \equiv {(q^{2} - k_{n}^{2})}^{1 / 2}

.

This can be reduced to:

Eq. (20)

M_{m n} (q) = [\frac{1}{Q_{m} (q) + Q_{n}^{*} (q)}] .

Thus the desired image of the absorbed optical energy is obtained by:

Eq. (21)

A (r) = \frac{1}{(2 π)^{2}} \int d^{2} q \exp (i q . ρ) \sum_{m, n} K_{m}^{*} (q, z) M_{m n}^{- 1} (q) \frac{P_{n} (q)}{S_{n} (q)},

where

ρ

is the 2-D spatial vector.

The computation of the inverse of the matrix $M$ is the key step in the procedure to recover the absorbed optical energy. Its pseudoinverse can be computed by performing the SVD of the matrix $K$ . Some of the eigenvalues of this matrix may be zero or very small. One will need to use regularization to circumvent this problem.

Thus, the procedure for recovering $A (r)$ is:

1. Take the Fourier transform of the measured pressure (at $z = 0$ ) with respect to time and the 2-D spatial components to obtain $\tilde{p} (k_{x}, k_{y}, ω)$ .
2. Divide $\tilde{p} (k_{x}, k_{y}, ω)$ by the scaling factor $S_{n} (q)$ .
3. Compute the matrix $K_{n} (q, z)$ by discretizing $q$ and using Eq. (15).
4. Compute the pseudoinverse of matrix $M_{m n} (q)$ .
5. Compute $\tilde{A} (q, z)$ using Eq. (17).
6. Take the 2-D inverse Fourier transform of $\tilde{A} (q, z)$ to obtain $A (r)$ .

2.1.3.

Calculation of eigenvalues of M matrix with nonzero attenuation

To investigate the posedness of the inverse problem we examined the properties of the matrix $K$ and its singular values for a medium with nonzero attenuation. We computed this matrix for a typical tissue with attenuation $α_{0} = 0.1 dB MHz^{- 1} {cm}^{- 1}$ and $n = 1$ . A set of 250 discrete temporal frequencies and 64 discrete spatial frequencies was used. Temporal frequencies ranged from 0 to 5 MHz and spatial frequencies $q$ ranged from 0 to $123.7 {cm}^{- 1}$ . This corresponds to the transducer elements spaced 0.05 cm apart. We performed an SVD of the matrix $K$ . The eigenvalues of the matrix $M$ were obtained from the singular values of matrix $K$ .

The variation of eigenvalues of $M$ for various spatial frequencies “ $q$ ” is shown in Fig. 1. Figure 2 shows the variation of the eigenvectors of matrix $M$ with depth for a specific eigenvalue. From these plots one concludes that:

1. Smaller values of 2-D spatial-frequency wave vector “ $q$ ” are recovered much better.
2. The eigenvectors that are nonzero at greater depths correspond to smaller eigenvalues.

The behavior of eigenvectors and eigenvalues indicates that shallower objects can be recovered much better than deeper objects because they correspond to eigenvectors with much higher eigenvalues. We also observe that eigenvectors that correspond to smaller eigenvalues also have much higher frequency components. This implies that resolution will degrade as you go deeper and that the higher frequencies cannot be recovered as accurately.

Fig. 1

Plot of eigenvalues of matrix $M$ for different spatial frequencies $q$ .

Fig. 2

Variation of select eigenvectors of matrix $M$ with depth (a) 10th eigenvalue, (b) 125th eigenvalue, (c) 230th eigenvalue.

2.2.

Image-Reconstruction Details

To investigate the performance of the SVD-based algorithm, several reconstructions were performed. The sensor data was simulated using the k-Wave toolbox,¹⁴^,¹⁵ which uses the k-space method¹³ to model the propagation of optoacoustic waves in lossy media. Both a point source placed at different depths and a three-dimensional (3-D) phantom were used to define the initial pressure distribution. For the point source, a grid size of $64 \times 64 \times 32$ voxels was used ( $32 \times 32 \times 16 mm$ ) supporting temporal frequencies up to 1.5 MHz. The temporal sampling interval was set to 100 ns and the simulations performed using 256 time steps. For the 3-D phantom, a grid size of $128 \times 128 \times 32$ voxels was used ( $64 \times 64 \times 16 mm$ ) with the same temporal sampling interval and 512 time steps. In both cases the sensor data was recorded in the plane $z = 0$ , corresponding to $64 \times 64$ sensor positions for the point-source simulations and $128 \times 128$ sensor positions for the 3-D phantom. The attenuation coefficient was set to $α_{0} = 3.0 dB {MHz}^{- n} {cm}^{- 1}$ , a value of $n = 1.5$ was used for the power-law frequency dependence, and the speed of sound was set to $1500 m / s$ . For each investigated case, several reconstructions were performed. First, the SVD-based method was used by obtaining the inverse of the matrix $M$ using truncated SVD (as described in Sec. 2.1.2). The magnitude of the eigenvalues of $M$ for a specific $q$ was used to regularize the pseudoinverse of $M$ for that $q^{'}$ . For comparison, the images were also reconstructed using TR.¹^,⁵ The TR reconstructions were performed both with and without compensation for acoustic absorption. When absorption compensation was included, the reconstructions were regularized by filtering the absorption and dispersion terms within the governing equations using a Tukey window. This restricts the range of frequencies that are allowed to grow during the reconstruction. The filter cut-off frequency (which acts as the regularization parameter for the TR algorithm) was chosen by examining the spectrum of the recorded sensor data.⁵ To separate the effect of limited-view artifacts (which are common to all reconstruction algorithms for this sensor geometry), an ideal reconstruction was also generated using TR with the absorption parameters set to zero in both the forward and inverse problems. (This is referred to as the ideal-TR-based reconstruction in the following sections.)

3. Results

3.1.

Sensitivity of SVD-Based Reconstruction to Regularization Parameter

We found that the quality of the reconstructed images using our SVD-based method was very sensitive to the regularization parameter used. The resulting magnitude and sharpness of the reconstructed images varied with this parameter. We chose the regularization parameter to depend on $q$ and related it to a certain percentage of total eigenvalue energy at that $q$ . Specifically, the truncation parameter $ϵ$ for a given matrix $M$ for a given $q$ was chosen to be:

Eq. (22)

ϵ = total (eigenvalue (q) energy) \times f,

where

0 < f \leq 1

. The eigenvalue energy, for a given

q

, for

n

eigenvalues was defined as:

Eq. (23)

eigenvalue energy = \sum_{i = 1}^{n} eigenvalue (i)^{2} .

Figure 3 shows how the reconstructed SVD-based images for a point source, at two different depths, vary with

f

. The regularization parameter that gives the best image also depends on the value of attenuation

α_{0}

and how much noise is present in the data, and it has to be adjusted when noise level or attenuation changes.

Fig. 3

Variation of SVD-based images with the regularization parameter as a percentage of total energy of the eigenvalues.

In general, there is no clear rule to help us decide what is a good value for the regularization parameter. It depends on how one evaluates image quality and what is the detection task. Image quality assessment depends both on the task and the observer,¹⁶ and a given technique may not perform well for all possible applications. One has to consider various image-quality aspects like spatial resolution, contrast, signal-to-noise ratio (SNR), artifacts, suppression of noise, signal detectability, quantitative accuracy, etc. We found that even in the absence of random noise in the data, we still had to use regularization as defined above to obtain a good-quality image that matches well with the true phantom. The need for regularization in this case may be arising from the very nature of the matrix $K$ , since it is exponentially damped. It could also arise from the effect of rounding and discretization errors. Further investigation is needed to understand and address this issue.

3.2.

Noiseless Data

In the following results, a value of $f = 0.99996$ was used for the SVD-based images. All the images were reconstructed using attenuated-pressure data generated using the $k$ -Wave toolbox. The cropped reconstructed images using SVD-based reconstruction and the TR algorithm are shown below in Fig. 4. The reconstructed images using SVD bear a close resemblance to the ideal image.

Fig. 4

Reconstructed $z$ -slice at a depth of 0.25 cm of a 3-D phantom, left to right: true phantom, SVD-based image, ideal-TR-based image, corrected-TR-based image, uncorrected-TR-based image.

Figure 5 shows the line used for the y-profile shown in Fig. 6, in the slices $z = 0.25$ and $z = 1.0 cm$ of the reconstructed $A (r)$ for noiseless data.

Fig. 5

True phantom showing the line used for y-profile.

Fig. 6

Line profile through phantom at depths of 0.25 and 1.0 cm for noiseless attenuated pressure data.

3.3.

Noisy Data

Random Gaussian noise was added to the attenuated-pressure data, generated using the $k$ -Wave toolbox, to obtain a SNR of 45 dB. A value of $f = 0.99996$ was used for regularization in the SVD-based method for the slice at 0.25-cm depth. A value of $f = 0.9999$ was used for regularization in the SVD-based method for the slice at 1.0 cm depth so as to reduce noise amplification. Figure 7 shows the cropped reconstructed phantom images using noisy attenuated pressure.

Fig. 7

Reconstructed $z$ -plane at a depth of 0.25 cm reconstructed using noisy attenuated pressure, left to right: true phantom, SVD-based image, ideal-TR-based image, corrected-TR-based image, uncorrected TR-based image.

Figure 8 shows the $y$ -profile along the line shown in Fig. 5, in the slices $z = 0.25$ and $z = 1.0 cm$ of the reconstructed absorbed optical energy for noisy data.

Fig. 8

Line profile through phantom at depths of 0.25 and 1.0 cm for noisy attenuated-pressure data.

The reconstructed images based on the noisy simulatedpressures obtained via the SVD-based method are in good agreement with the ideal image. The SVD-based image-reconstruction algorithm offers comparable images even with noisy data, although it is more sensitive to the presence of noise.

3.4.

Depth Resolution

Figure 9 shows the profile through a point source placed at various depths reconstructed using SVD-based attenuated pressure. The plot labeled “True” refers to the original point source. The label “TR-based” refers to corrected TR-based reconstructed image, while the label “Uncorrected TR-based” refers to uncorrected TR-based reconstructed image.

Fig. 9

Line profile of a point source placed at a depth of 0.15 and 0.65 cm for images reconstructed using: TR method corrected for attenuation, TR method, and SVD-based method.

Table 1 lists the full width at half maxima (FWHM) for the point-source images constructed using the different algorithms.

Table 1

FWHM in cm for a point-source image at different depths constructed using SVD-based and TR-based algorithms.

Depths	0.15 cm	0.35 cm	0.65 cm
SVD-based	0.0500	0.0500	0.0624
Corrected TR-based	0.0515	0.0542	0.0620
Uncorrected TR-based	0.0546	0.0603	0.0845

The SVD-based image reconstruction algorithm shows very good depth resolution, especially for depths up to 0.35 cm. At greater depths this algorithm offers worse resolution when compared to images reconstructed using the TR-based method with or without attenuation correction. However, the intensity of the image reconstructed using the SVD-method at greater depths relative to shallower depths is much higher than that obtained from the TR-based method with attenuation correction. As mentioned before, the depth resolution of the images reconstructed via the SVD-based method are also dependent on the choice of the regularization parameter. In practice, we may be able to attain better depth resolution using the SVD-based method by adjusting this parameter.

4. Discussion

Our SVD-based method provides insight into the ill-conditioning of the inverse problem of image reconstruction in PAT in a lossy medium. However, there are several important issues to consider in using this approach. One needs to address the choice of regularization method used to obtain the matrix inverse of matrix $M$ . We found that there was a need for regularization even in the absence of noise. We used a truncated SVD method for regularization. The truncation of eigenvalues was based on a certain percentage of total energy contained in all the eigenvalues. Our choice of this type of truncation was based on obtaining good resolution at greater depths. However, there are many other methods for regularizing, and there are several excellent references in the literature¹⁷ that address this issue extensively. A systematic study of the choice of regularization parameter for our SVD-based image-reconstruction method remains a subject for further study.

Both the TR-based and the SVD-based image-reconstruction techniques use regularization to prevent high-frequency noise present in the sensor data being amplified during the reconstruction. However, there is a distinct difference between these techniques in how this regularization is applied. In the TR-based approach, the range of frequencies that is allowed to grow is restricted by filtering the absorption parameters within the governing equations. This restriction applies to signals originating from all depths within the medium. However, in practice, some of the high frequency components from shallow features within the image may be above the noise floor, while the same frequency components from deeper features may not. This method of regularization does not allow for the selective compensation of the attenuation at these frequencies. Conversely, in the SVD-based approach, the regularization is based on the magnitude of the eigenvalues where the corresponding eigenvectors have features at different depths. This allows for the compensation of the attenuation of high-frequency components originating from shallow features, even if the corresponding frequency components from deeper features are not above the noise floor. Finally, it is useful to note that, while the SVD-based method is computationally intensive, this may be justified in cases where the quantitative accuracy of the reconstructed pressure distribution is of importance. For example, in quantitative photoacoustics where the goal is to derive accurate chromophore concentrations,¹⁸ errors in the reconstructed pressure due to acoustic attenuation will ultimately result in errors in the extracted quantitative information.

5. Conclusions

We derived an operator relating the attenuated pressure to the absorbed optical energy to reconstruct images in PAT in an attenuating medium for a planar geometry. We derived an SVD-based algorithm based on this operator to reconstruct the absorbed optical energy in PAT in an attenuating medium. We looked at the eigenvalues of the matrix $M$ which is used to recover the absorbed optical energy. We found that smaller values of the 2-D spatial frequency vector $q$ are recovered much better than the larger values. We also observed that the smaller eigenvalues typically corresponded to eigenvectors with most of their energy at greater depths. The behavior of eigenvectors and eigenvalues indicated that shallower objects can be recovered much better than deeper objects because they correspond to eigenvectors with much higher eigenvalues. The SVD-based approach to the attenuation problem in PAT offers a promising method for direct image reconstruction in a planar geometry. It also provides a way to study the conditioning of this inverse problem. This SVD-based algorithm shows good depth resolution and noise stability. However, we found that the image quality of the reconstructed images was very sensitive to the choice of regularization parameter used for obtaining the inverse of $M$ via SVD.

Acknowledgments

Dimple Modgil would like to thank Mark A. Anastasio for helpful discussions related to this work. This work was supported in part by a breast cancer predoctoral research award (W81XWH-08-1-0331) to Dimple Modgil.

References

1.

B. E. Treebyet al., “Acoustic attenuation compensation in photoacoustic tomography: application to high-resolution 3d imaging of vascular networks in mice,” Proc. SPIE, 7899 78992Y (2011). http://dx.doi.org/10.1117/12.874530 Google Scholar

2.

P. J. La RivièreJ. ZhangM. A. Anastasio, “Image reconstruction in optoacoustic tomography for dispersive acoustic media,” Opt. Lett., 31 (6), 781 –783 (2006). http://dx.doi.org/10.1364/OL.31.000781 OPLEDP 0146-9592 Google Scholar

3.

P. Burgholzeret al., “Compensation of acoustic attenuation for high-resolution photoacoustic imaging with line detectors,” Proc. SPIE, 6437 643724 (2007). http://dx.doi.org/10.1117/12.700723 Google Scholar

4.

P. Burgholzeret al., “Information changes and time reversal for diffusion-related periodic fields,” Proc. SPIE, 7177 717723 (2009). http://dx.doi.org/10.1117/12.809074 Google Scholar

5.

B. E. TreebyE. Z. ZhangB. T. Cox, “Photoacoustic tomography in absorbing acoustic media using time reversal,” Inv. Prob., 26 (11), 115003 (2010). http://dx.doi.org/10.1088/0266-5611/26/11/115003 INPEEY 0266-5611 Google Scholar

6.

H. RoitnerP. Burgholzer, “Efficient modeling and compensation of ultrasound attenuation losses in photoacoustic imaging,” Inv. Prob., 27 (1), 015003 (2011). http://dx.doi.org/10.1088/0266-5611/27/1/015003 INPEEY 0266-5611 Google Scholar

7.

V. A. MarkelJ. C. Schotland, “Inverse problem in optical diffusion tomography. I. Fourier-Laplace inversion formulas,” J. Opt. Soc. Am., 18 1336 –1347 (2001). http://dx.doi.org/10.1364/JOSAA.18.001336 JOSAAH 0030-3941 Google Scholar

8.

V. A. MarkelV. MitalJ. C. Schotland, “Inverse problem in optical diffusion tomography. III. inversion formulas and singular-value decomposition,” J. Opt. Soc. Am. A, 20 890 –902 (2006). http://dx.doi.org/10.1364/JOSAA.20.000890 JOAOD6 0740-3232 Google Scholar

9.

T. L. Szabo, “Time-domain wave equations for lossy media obeying a frequency power law,” J. Acoust. Soc. Am., 96 (1), 491 –500 (1994). http://dx.doi.org/10.1121/1.410434 JASMAN 0001-4966 Google Scholar

10.

N. V. SushilovR. S. C. Cobbold, “Frequency-domain wave equation and its time-domain solutions in attenuating media,” J. Acoust. Soc. Am., 115 1431 –1436 (2004). http://dx.doi.org/10.1121/1.1675817 JASMAN 0001-4966 Google Scholar

11.

F. A. Duck, Physical Properties of Tissue: A Comprehensive Reference, Academic, London (1990). Google Scholar

12.

Y. XuD. FengL. V. Wang, “Exact frequency-domain reconstruction for thermoacoustic tomography, I: planar geometry,” IEEE Trans. Med. Imag., 21 823 –828 (2002). http://dx.doi.org/10.1109/TMI.2002.801172 ITMID4 0278-0062 Google Scholar

13.

C. Matson, “A diffraction tomographic model of the forward problem using diffuse photon density waves,” Opt. Exp., 1 6 –11 (1997). http://dx.doi.org/10.1364/OE.1.000006 OPEXFF 1094-4087 Google Scholar

14.

B. E. TreebyB. 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

15.

B. E. TreebyB. T. Cox, “Modeling power law absorption and dispersion for acoustic propagation using the fractional laplacian,” J. Acoust. Soc. Am., 127 (5), 2741 –2748 (2010). http://dx.doi.org/10.1121/1.3377056 JASMAN 0001-4966 Google Scholar

16.

H. H. BarrettK. J. Myers, Foundations of Image Science, Wiley Interscience, Hoboken, NJ (2004). Google Scholar

17.

P. C. Hansen, Rank-Deficient and Discrete Ill-Posed Problems, SIAM, Philadelphia (1998). Google Scholar

18.

J. Lauferet al., “Quantitative determination of chromophore concentrations from 2D photoacoustic images using a nonlinear model-based inversion scheme,” Appl. Opt., 49 (8), 1219 –1233 (2010). http://dx.doi.org/10.1364/AO.49.001219 APOPAI 0003-6935 Google Scholar

Citation Download Citation

Dimple Modgil, Patrick J. La Rivière, and Bradley E. Treeby "Photoacoustic image reconstruction in an attenuating medium using singular-value decomposition," Journal of Biomedical Optics 17(6), 061204 (7 May 2012). https://doi.org/10.1117/1.JBO.17.6.061204

Published: 7 May 2012

Access the abstract

JOURNAL ARTICLE
9 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 21 scholarly publications and 1 patent.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Signal attenuation

Reconstruction algorithms

Image restoration

Photoacoustic spectroscopy

Lanthanum

Inverse optics

Absorption

1.

Introduction

2.

Correction for Ultrasonic Attenuation in OAT Using SVD

2.1.

Methods

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

2.1.1.

Angular spectrum expansion of the measured pressure signals

Eq. (9)

Eq. (10)

Eq. (11)

2.1.2.

SVD of integral operator relating pressure to optical-absorption coefficient

Eq. (12)

Eq. (13)

Eq. (14)

Eq. (15)

Eq. (16)

Eq. (17)

Eq. (18)

Eq. (19)

Eq. (20)

Eq. (21)

2.1.3.

Calculation of eigenvalues of M matrix with nonzero attenuation

Fig. 1

Fig. 2

2.2.

Image-Reconstruction Details

3.

Results

3.1.

Sensitivity of SVD-Based Reconstruction to Regularization Parameter

Eq. (22)

Eq. (23)

Fig. 3

3.2.

Noiseless Data

Fig. 4

Fig. 5

Fig. 6

3.3.

Noisy Data

Fig. 7

Fig. 8

3.4.

Depth Resolution

Fig. 9

Table 1

4.

Discussion

5.

Conclusions

Acknowledgments

References

Show All Keywords

Keywords/Phrases

Search In:

Publication Years