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 groups1,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 Burgholzer6 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.
Correction for Ultrasonic Attenuation in OAT Using SVD
For soft biological tissue, the attenuation coefficient has a power-law dependence on the frequency given by:11Ref. 11).
Using the Kramers-Kronig relations, for power-law absorption with and , one obtains the dependence of phase velocity on frequency due to frequency-dependent attenuation as:9
For , the dependence of phase velocity on frequency due to frequency-dependent attenuation is given by:
So for an attenuating medium, Eq. (1) becomes:
This equation can be solved using standard Green’s function methods as:12
Angular spectrum expansion of the measured pressure signals
The angular spectrum expansion of the Green’s function is given by:138) and consider the pressure measurements on the plane and we get:
We assume that the photoacoustic object lies in the plane . On taking the Fourier transform on both sides and reducing the resulting expression, one obtains the angular spectrum expansion of measured photoacoustic pressure as:
SVD of integral operator relating pressure to optical-absorption coefficient
Define a two-dimensional (2-D) spatial wave vector as:
Define , . Thus, one can obtain an integral equation:8
We can obtain a simpler expression for by performing the integral over . Equation (18) can be written as:
This can be reduced to:
Thus the desired image of the absorbed optical energy is obtained by:
The computation of the inverse of the matrix 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 . 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 is:
1. Take the Fourier transform of the measured pressure (at ) with respect to time and the 2-D spatial components to obtain .
2. Divide by the scaling factor .
3. Compute the matrix by discretizing and using Eq. (15).
4. Compute the pseudoinverse of matrix .
5. Compute using Eq. (17).
6. Take the 2-D inverse Fourier transform of to obtain .
Calculation of eigenvalues of M matrix with nonzero attenuation
To investigate the posedness of the inverse problem we examined the properties of the matrix and its singular values for a medium with nonzero attenuation. We computed this matrix for a typical tissue with attenuation and . 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 ranged from 0 to . This corresponds to the transducer elements spaced 0.05 cm apart. We performed an SVD of the matrix . The eigenvalues of the matrix were obtained from the singular values of matrix .
The variation of eigenvalues of for various spatial frequencies “” is shown in Fig. 1. Figure 2 shows the variation of the eigenvectors of matrix with depth for a specific eigenvalue. From these plots one concludes that:
1. Smaller values of 2-D spatial-frequency wave vector “” 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.
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 method13 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 voxels was used () 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 voxels was used () with the same temporal sampling interval and 512 time steps. In both cases the sensor data was recorded in the plane , corresponding to sensor positions for the point-source simulations and sensor positions for the 3-D phantom. The attenuation coefficient was set to , a value of was used for the power-law frequency dependence, and the speed of sound was set to . For each investigated case, several reconstructions were performed. First, the SVD-based method was used by obtaining the inverse of the matrix using truncated SVD (as described in Sec. 2.1.2). The magnitude of the eigenvalues of for a specific was used to regularize the pseudoinverse of for that . 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.)
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 and related it to a certain percentage of total eigenvalue energy at that . Specifically, the truncation parameter for a given matrix for a given was chosen to be:Figure 3 shows how the reconstructed SVD-based images for a point source, at two different depths, vary with . The regularization parameter that gives the best image also depends on the value of attenuation and how much noise is present in the data, and it has to be adjusted when noise level or attenuation changes.
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 , 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.
In the following results, a value of was used for the SVD-based images. All the images were reconstructed using attenuated-pressure data generated using the -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.
Random Gaussian noise was added to the attenuated-pressure data, generated using the -Wave toolbox, to obtain a SNR of 45 dB. A value of was used for regularization in the SVD-based method for the slice at 0.25-cm depth. A value of 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.
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.
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.
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|
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.
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 . 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 literature17 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.
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 which is used to recover the absorbed optical energy. We found that smaller values of the 2-D spatial frequency vector 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 via SVD.
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.