*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.^{1} 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.^{2} 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^{1}^{,}^{3}^{–}4.^{5} 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^{6} 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.,^{7}^{,}^{8} 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^{2}^{,}^{9}^{,}^{10} and is given by:

## Eq. (1)

$${\nabla}^{2}p(\mathbf{r},t)-\frac{{\partial}^{2}p(\mathbf{r},t)}{{c}_{0}^{2}\partial {t}^{2}}+L(t)*p(\mathbf{r},t)=-\frac{\beta}{{C}_{p}}\frac{\partial}{\partial t}H(\mathbf{r},t),$$For soft biological tissue, the attenuation coefficient has a power-law dependence on the frequency given by:^{11}

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

## Eq. (5)

$$\frac{1}{c(\omega )}=\frac{1}{{c}_{0}}+{\alpha}_{0}\text{\hspace{0.17em}}\mathrm{tan}(\pi n/2)({\omega}^{n-1}-{\omega}_{0}^{n-1}),$$For $n=1$, the dependence of phase velocity $c(\omega )$ on frequency due to frequency-dependent attenuation is given by:

## Eq. (6)

$$\frac{1}{c(\omega )}=\frac{1}{{c}_{0}}-\frac{2}{\pi}{\alpha}_{0}\text{\hspace{0.17em}}ln\left|\frac{\omega}{{\omega}_{0}}\right|.$$So for an attenuating medium, Eq. (1) becomes:

## Eq. (7)

$${\u25bf}^{2}\tilde{p}(\mathbf{r},\omega )+{k}^{2}(\omega )\tilde{p}(\mathbf{r},\omega \text{})=-\frac{i\omega \beta}{{C}_{p}}A(\mathbf{r})\tilde{I}(\omega ),$$This equation can be solved using standard Green’s function methods as:^{12}

## Eq. (8)

$$\tilde{p}(\mathbf{r},\omega )=-i\eta \omega \tilde{I}(\omega )\iiint {\mathrm{d}}^{3}{r}^{\prime}A({\mathbf{r}}^{\prime})G(\mathbf{r}-{\mathbf{r}}^{\prime}),$$## 2.1.1.

#### Angular spectrum expansion of the measured pressure signals

The angular spectrum expansion of the Green’s function is given by:^{13}

## Eq. (9)

$$G(r)=\frac{\mathrm{exp}(i{k}^{2}(\omega )r)}{4\pi r}=\frac{1}{8{\pi}^{2}}\iint \frac{1}{{[{\alpha}_{x}^{2}+{\alpha}_{y}^{2}-{k}^{2}(\omega )]}^{1/2}}\phantom{\rule{0ex}{0ex}}\times \mathrm{exp}[-|z|{[{\alpha}_{x}^{2}+{\alpha}_{y}^{2}-{k}^{2}(\omega )]}^{1/2}+ix{\alpha}_{x}+iy{\alpha}_{y}]\mathrm{d}{\alpha}_{x}\mathrm{d}{\alpha}_{y}.$$## Eq. (10)

$$\tilde{p}(x,y,0,\omega )=\frac{-i\eta \omega \tilde{I}(\omega )}{8{\pi}^{2}}\iint {\mathrm{d}}^{3}{r}^{\prime}A({x}^{\prime},{y}^{\prime},{z}^{\prime})\phantom{\rule{0ex}{0ex}}\iint \mathrm{d}{\alpha}_{x}\mathrm{d}{\alpha}_{y}\frac{1}{{[{\alpha}_{x}^{2}+{\alpha}_{y}^{2}-{k}^{2}(\omega )]}^{1/2}}\phantom{\rule{0ex}{0ex}}\times \mathrm{exp}[-|{z}^{\prime}|{{[\alpha}_{x}^{2}+{\alpha}_{y}^{2}-{k}^{2}(\omega )]}^{1/2}\phantom{\rule{0ex}{0ex}}+i(x-{x}^{\prime}){\alpha}_{x}+i(y-{y}^{\prime}){\alpha}_{y}].$$We assume that the photoacoustic object lies in the plane ${z}^{\prime}\ge 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},\omega )=\frac{-i\eta \omega \stackrel{~}{I}(\omega )}{2{[{k}_{x}^{2}+{k}_{y}^{2}-{k}^{2}(\omega )]}^{1/2}}\int \mathrm{d}{z}^{\prime}\stackrel{~}{A}({k}_{x},{k}_{y},{z}^{\prime})\phantom{\rule{0ex}{0ex}}\mathrm{exp}[-{z}^{\prime}{[{k}_{x}^{2}+{k}_{y}^{2}-{k}^{2}(\omega )]}^{1/2}].$$## 2.1.2.

#### SVD of integral operator relating pressure to optical-absorption coefficient

We follow Markel’s^{7}^{,}^{8} 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 $\mathbf{q}$ as:

We then discretize the temporal frequency, $\omega $ as ${\omega}_{n}$. For notational simplicity, we use $z$ instead of ${z}^{\prime}$ in the following equations. One can then write the pressure as:## Eq. (13)

$$\tilde{p}(\mathbf{q},{\omega}_{n})=\frac{-i\eta {\omega}_{n}\tilde{I}({\omega}_{n})}{2{[{q}^{2}-{k}^{2}({\omega}_{n})]}^{1/2}}\int \mathrm{d}z\tilde{A}(\mathbf{q},z)\phantom{\rule{0ex}{0ex}}\mathrm{exp}[-z{[{q}^{2}-{k}^{2}({\omega}_{n})]}^{1/2}],$$Define ${P}_{n}(\mathbf{q})\equiv \stackrel{~}{p}(\mathbf{q},{\omega}_{n})$, ${k}_{n}\equiv k({\omega}_{n})$. Thus, one can obtain an integral equation:

## Eq. (14)

$${P}_{n}(\mathbf{q})={S}_{n}(\mathbf{q})\int \mathrm{d}z\stackrel{~}{A}(\mathbf{q},z){K}_{n}(\mathbf{q},z),$$## Eq. (16)

$${S}_{n}(\mathbf{q})\equiv \frac{-i\eta {\omega}_{n}\stackrel{~}{I}({\omega}_{n})}{2{({q}^{2}-{k}_{n}^{2})}^{1/2}},$$^{8}

## Eq. (17)

$$\stackrel{~}{A}(\mathbf{q},z)=\sum _{m,n}{K}_{m}^{*}(\mathbf{q},z){M}_{mn}^{-1}(\mathbf{q})\frac{{P}_{n}(\mathbf{q})}{{S}_{n}(\mathbf{q})},$$## Eq. (18)

$${M}_{mn}(\mathbf{q})\equiv \int \mathrm{d}z{K}_{m}(\mathbf{q},z){K}_{n}^{*}(\mathbf{q},z),$$We can obtain a simpler expression for ${M}_{mn}(\mathbf{q})$ by performing the integral over $z$. Equation (18) can be written as:

## Eq. (19)

$${M}_{mn}(\mathbf{q})={\int}_{0}^{\infty}\mathrm{d}z\text{\hspace{0.17em}}\mathrm{exp}\text{\hspace{0.17em}}[-z\{{Q}_{m}(\mathbf{q})+{Q}_{n}^{*}(\mathbf{q})\}],$$This can be reduced to:

## Eq. (20)

$${M}_{mn}(\mathbf{q})=\left[\frac{1}{{Q}_{m}(\mathbf{q})+{Q}_{n}^{*}(\mathbf{q})}\right].$$Thus the desired image of the absorbed optical energy is obtained by:

## Eq. (21)

$$A(\mathbf{r})=\frac{1}{(2\pi {)}^{2}}\int {\mathrm{d}}^{2}q\text{\hspace{0.17em}}\mathrm{exp}(i\mathbf{q}.\mathbf{\rho})\sum _{m,n}{K}_{m}^{*}(\mathbf{q},z){M}_{mn}^{-1}(\mathbf{q})\frac{{P}_{n}(\mathbf{q})}{{S}_{n}(\mathbf{q})},$$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(\mathbf{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 $\stackrel{~}{p}({k}_{x},{k}_{y},\omega )$.

2. Divide $\stackrel{~}{p}({k}_{x},{k}_{y},\omega )$ by the scaling factor ${S}_{n}(\mathbf{q})$.

3. Compute the matrix ${K}_{n}(\mathbf{q},z)$ by discretizing $\mathbf{q}$ and using Eq. (15).

4. Compute the pseudoinverse of matrix ${M}_{mn}(\mathbf{q})$.

5. Compute $\stackrel{~}{A}(\mathbf{q},z)$ using Eq. (17).

6. Take the 2-D inverse Fourier transform of $\stackrel{~}{A}(\mathbf{q},z)$ to obtain $A(\mathbf{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 ${\alpha}_{0}=0.1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{dB}\text{\hspace{0.17em}}\mathrm{MHz}{}^{-1}\text{\hspace{0.17em}}{\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{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.

## 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,^{14}^{,}^{15} which uses the k-space method^{13} 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\text{\hspace{0.17em}}\times 16\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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 ${\alpha}_{0}=3.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{dB}\text{\hspace{0.17em}}{\mathrm{MHz}}^{-n}\text{\hspace{0.17em}}{\mathrm{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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}/\mathrm{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}^{\prime}$. For comparison, the images were also reconstructed using TR.^{1}^{,}^{5} 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.^{5} 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 $\mathbf{q}$ and related it to a certain percentage of total eigenvalue energy at that $\mathbf{q}$. Specifically, the truncation parameter $\u03f5$ for a given matrix $M$ for a given $q$ was chosen to be:

where $0<f\le 1$. The eigenvalue energy, for a given $\mathbf{q}$, for $n$ eigenvalues was defined as:## Eq. (23)

$$\mathrm{eigenvalue}\text{\hspace{0.17em}}\mathrm{energy}=\sum _{i=1}^{n}\mathrm{eigenvalue}(i{)}^{2}.$$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,^{16} 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.

Figure 5 shows the line used for the y-profile shown in Fig. 6, in the slices $z=0.25$ and $z=1.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{cm}$ of the reconstructed $A(\mathbf{r})$ for noiseless 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.

Figure 8 shows the $y$-profile along the line shown in Fig. 5, in the slices $z=0.25$ and $z=1.0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{cm}$ of the reconstructed absorbed optical energy for noisy 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.

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^{17} 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,^{18} 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 $\mathbf{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.