|
1.IntroductionOptoacoustic imaging is a hybrid imaging technique that has attracted a lot of attention in recent years. It is based on the thermoacoustic effect, which refers to acoustic wave generation on absorption of pulsed optical energy by a medium. A slight rise in temperature of the medium due to the absorption of the incident electromagnetic wave results in thermoelastic expansion. This thermoelastic expansion and then contraction due to the pulsed electromagnetic waves leads to the generation of acoustic waves. Under the constraints of thermal and stress confinement, this thermal expansion leads to a rise in pressure, , that satisfies the three-dimensional inhomogeneous wave equation1 where is the speed of sound; , the heating function, is the thermal energy deposited by the electromagnetic radiation per unit time per unit volume; is the isobaric volume expansion coefficient; and is the specific heat of the medium. The heating function can be expressed as the product of a spatially varying absorbed optical energy, , and a time-dependent optical illumination function, . If the speed of sound is constant, , the measured pressure signal can be related to the optical absorption function, assuming delta pulse illumination, aswhere . Equation 2 states that the time integral of acoustic pressure at a point and time is given by the integral of the optical absorption function over a spherical surface of radius centered at . A simple but inexact way to reconstruct is to spatially resolve the optoacoustic waves and to backproject time-integrated pressure signals over hemispheres of radii .However, the speed of sound in tissues can vary from (in fat) to about (for skin) for ultrasound waves in the range.2 Using a constant speed of sound in image reconstruction can result in image distortion and poor image quality. There has been some work done previously to address the effect of speed-of-sound variations in optoacoustic tomography (OAT). Xu3 and Xu and Wang4 looked at variable speed-of-sound reconstruction in breast thermoacoustic tomography (TAT). They concluded that there is minor amplitude distortion in breast TAT and that variation in travel time is typically a first-order perturbation in a weakly acoustically heterogeneous medium. Zhang and Anastasio5 derived a heuristic method for reconstructing acoustic speed and optical absorption distributions in a stepwise manner. Their work was also based on the first-order travel-time effect. Kolkman 6 devised a method to determine the speed of sound in tissue using two concentric ring-based photoacoustic sensors based on the first-order travel-time effect. Manohar 7 also devised a new and improved method to determine the speed of sound using a photoacoustic imager. Thus, overall previous work on the speed-of-sound variation that is based on the generalized radon transform (GRT) model has only looked at the first-order effect of variable speed of sound.4, 5, 6 In the GRT model, the pressure is given by where is the travel time for the sound wave to travel between points and . In the first-order models, it was assumed that sound rays continue to travel in a straight line in the presence of the acoustic heterogeneity and the travel time between points and is given as a line integral over the slowness mapwhere is the straight line joining the points and , is the differential length element along that line, and is the variable speed of sound. The effect of variable speed of sound on signal amplitude was neglected. Previous work based on the first-order GRT model does not incorporate the effect of ray bending, which may be significant when the speed of sound varies by 10% or more.There are some other approaches that make less restrictive assumptions about the speed of sound, although they usually require the assumption of a closed detector surface, which is not always achievable in practice. Agranovsky and Kuchment8 have recently derived an analytic reconstruction formula for OAT for arbitrary detector geometry, as long as the point detectors are placed on a closed surface. It accommodates a variable speed of sound. Their analytic formula leads to reconstruction in terms of an eigenfunction expansion. Zhen and Jiang9 devised an iterative reconstruction algorithm, using a finite element method that incorporates both attenuation and variable speed of sound. Recently, Hristova 10 have used the time-reversal method to incorporate variable speed of sound in optoacoustic image reconstruction. Their method is exact and performs well when the speed-of-sound map is known. However, it has the usual limitations of the time-reversal algorithm that the detection surface needs to be a closed surface, enclosing the object to avoid artifacts due to incomplete data. In this work, we investigate a higher-order geometrical acoustic approximation that incorporates a first-order effect on the signal amplitude and higher-order effects on the travel times. We use this higher-order approximation to construct the system matrix. We then iteratively reconstruct the images using third-, second-, first-, and zeroth-order travel-time effects. We show that the higher-order GA approximation offers much better image quality and accuracy when the speed of sound varies by 5% or more. 2.MethodsWe treat the variation in the speed of sound as a perturbation to the uniform background sound speed, , where characterizes the magnitude of the perturbation.We denote the unperturbed Green’s function by . It satisfies the background homogeneous Helmholtz equation The Green’s function that satisfies the perturbed inhomogeneous Helmholtz equation is given by 2.1.Derivation of GA ApproximationThe GA approximation ignores diffraction effects. It is valid in the short wavelength regime when the size of the inhomogeneity is much greater than the wavelength. In a scattering medium, this approximation is valid when the speed of sound does not change significantly over one wavelength. Under this approximation, the Green’s function can be written as where is the wave number.This can be written in the zeroth-order approximation in wavelength as where is known as the eikonal function. On substituting this into Eq. 7 and equating equal powers in we get in the limitEquation 10 is called the eikonal equation, and it implies where denotes the arc length along the path between the points and as shown in Fig. 1 . Note that the curve connecting the points and is not, in general, a straight line. The ray trajectory is determined by the path that minimizes the acoustic path length (Fermat’s principle) or equivalently, that minimizes the travel time. For an acoustically homogeneous medium, this is a straight line. Previous work by Snieder and Aldridge11 concludes that to first order in perturbation, this trajectory can be chosen to be the path along the reference ray that satisfies the eikonal equation (assuming the variation in speed of sound is slowly varying).Using standard Green’s function techniques12, 13 and Eq. 9, the pressure is given by Taking the Fourier transform with respect to , we obtain the generalized radon transform (GRT) Equation 11 can be solved to obtain where is the constant of integration. Because this equation must also be satisfied if , we get2.2.Higher-Order Geometrical Acoustics ApproximationWe can thus improve the first-order GA model that has been used previously by incorporating higher-order effects on the amplitude and travel times. We can use Eq. 18 for the amplitude of the Green’s function. We can incorporate the effect of ray bending by considering the higher-order perturbations in the eikonal. The assumption for first-order GA is that the speed of sound is slowly varying so that the time of travel can be obtained using linear rays. However, if this assumption is not true, it can lead to higher-order perturbations in travel times. Higher-order travel-time perturbations contribute to ray bending. The ray bends toward the region that has a higher refractive index (or lower speed of sound). Assume that the reference speed of sound is perturbed as given by Eq. 5. Following the methodology of Snieder and Aldridge,11 this leads to a perturbation in travel time as Let denote the reference ray associated with the reference eikonal . Let the reference ray be parametrized by the variable so that , where is the unit vector along the reference ray. If we substitute for and in the eikonal equation and collect terms with equal powers of the perturbation and solve, we get andThus, one can calculate the perturbed time of flight using these equations if one knows the reference ray. Note that this perturbation theory approach only works if the nonlinear perturbations are small and there is no multipathing.11 Thus, we can obtain a higher-order estimate of the pressure than that given by Eq. 24 by using Eq. 18 for amplitude and Eq. 19 and keeping up to second-order terms in Eq. 15. One can obtain even higher-order perturbations to travel times as described by Snieder and Aldridge11 These effects are expected to be smaller than the second-order travel-time effects. The equations for third- and fourth-order effects are given in the Appendix{ label needed for app[@id='x0'] }. As an aside, note that if we keep only the zeroth-order term in , we get from Eq. 18 Using this expression for in Eq. 15, we will obtain the following expression for pressure: This expression along with the first-order travel time given by Eq. 4 is what has been used implicitly in the previous work addressing speed-of-sound variation based on the GRT model. 2.3.Limits of GA Approximation2.3.1.Requirement for validity of GAThe GA approximation is valid under certain assumptions about the rate of variation of the speed of sound and the size of acoustic inhomogeneity. It requires that the scale of variation of the speed of sound should be much greater than the acoustic wavelength. Let the speed of sound be given by where .GA requires that14 This implies that If denotes the scale of variation of speed of sound, then this implies that For a typical ultrasound wavelength of , this limit implies For a 10% variation, i.e., , we get 2.3.2.Regime where ray bending may become importantIn the GRT model, the travel times are used to determine where the optoacoustic signal at a certain transducer position at a certain time came from. If the travel times calculated using the higher-order GA approximation differ from those calculated using the first-order approximation by greater than the temporal sampling interval, then this would lead to incorrect calculation of . Second-order travel-time effects incorporate the leading effect of ray bending. Thus, ray bending becomes important if the second-order correction to the travel time is of the order of the temporal sampling interval . This implies that But, Eq. 22 for implies that where is the length of the inhomogeneous region.This implies that, in order for the second-order effect to be significant, we must have and for , , Eq. 33 implies thatand for a 10% variation (i.e., ), we getThus, even a inclusion of 10% variation above the background speed of sound could cause ray bending to be significant. 2.4.Computer-Simulation StudiesWe performed simulations to study the effect of the higher-order GA approximation on travel times and the resulting reconstructed images. We wanted to validate that our perturbative approach to travel-time calculation matches the analytically computed travel times for a given speed of sound map. We explored some speed-of-sound maps where ray bending may become important. We wanted to study the effect of ray bending on the computed travel times and the resultant reconstructed images. We computed the travel times in 2-D for two different speed-of-sound maps, which are shown in Fig. 2 . The first speed-of-sound map is a continuous, slowly varying region as shown in Fig. 2. The speed of sound in this case is given by where is a constant that controls the variation in speed of sound. The travel times for this speed-of-sound map can be calculated analytically as described in Ref. 15. We wanted to see how the perturbed travel times compare to the analytical travel times for this slowly varying speed-of-sound map.The second speed-of-sound map contains an elliptical acoustic inhomogeneity as shown in Fig. 2, which may be more common clinically. The blurred elliptical region in this map has the variable speed of sound with respect to the background. The background speed of sound was set to . We were not able to compute the travel times analytically in this case for comparison. The travel times were computed assuming that the transducers are placed on a line to the left of the phantom. These travel times were calculated using the higher-order GA approximation in 2-D using the Eqs. 20, 21, 22, 44, 45. The reference ray path was assumed to be a straight line. The line integrals over the reference ray path were evaluated numerically. 2.4.1.Travel-time calculationsTravel times were calculated using up to fourth-order correction for different pixels for a specific transducer position. This was done in 2-D for a grid for the two specific speed-of-sound maps as shown in Fig. 2. The travel times were calculated for a 5, 10, and 15% speed-of-sound variation. The pixel size was set to . The time-sampling interval was set to . The transducers were assumed to lie on a line ( , ) to the left of the phantom, apart. 2.4.2.Image reconstructionThe images were then reconstructed iteratively using the least-squares (LS) method and conjugate gradient method. Iterative methods require the computation of system matrices. The PAT reconstruction problem can be written as16 where is the time-integrated pressure data, is the object function that represents the absorbed optical energy distribution, and is the system matrix. The function is defined asWe consider the 2-D geometry where the inherently 3-D OAT problem is reduced to 2-D due to the transducer’s directivity. In this case,17 for a transducer along the line , using the higher-order GA approximation, can be written as The system matrix is computed by calculating the contribution of pressure signal from each pixel for a given time and transducer position. The contribution from each pixel is calculated on the basis of travel time from that pixel to a given transducer position. The system matrices were constructed in 2-D incorporating zeroth, first-, second-, third-, and fourth-order travel-time effects. The fourth-order travel-time system matrix was used for the forward model using the phantom in Fig. 2. The LS method minimizes the difference between the observed data and that obtained by projecting the object data via the system matrix. This can be written symbolically as16 Equation 39 has the solution where is the pseudo-inverse of . The iterative LS method is based on the following additive update formula:16 where is the object estimate, is the iteration number, is the step size, and is the step direction. The conjugate gradient method was used to find the new step direction. In this algorithm, the step directions are chosen so that they are conjugate to each other. This method, as outlined in Chapter 21 of Ref. 16, was used by us to perform the iterative reconstruction.The images were then reconstructed iteratively using third-, second-, first-, and zeroth-order matrices. The image reconstruction was done on a grid with parameters as specified above. Noisy data were generated by adding Gaussian noise that was 1% of the maximum value of the noiseless fourth-order generated forward data. 3.ResultsThe results below are for the specific speed-of-sound maps as shown in Figs. 2 and 2. The speed of sound in the ellipse varies with respect to the background. The line of transducers is placed to the left of the phantom along the line . Some of these results related to the elliptical speed-of-sound map are slightly different from those in Ref. 18. This is due to some minor improvements in the implementation of the travel-time calculations related to how accurately the slowness map is integrated over the line joining the two points. 3.1.Relative Travel TimesIn order to compare the calculated travel times with actual data, we used a speed-of-sound map that would allow for an analytical calculation of travel times. We specifically looked at the speed-of-sound map shown in Fig. 2. The relative travel times were calculated using where refers to a specific-order perturbative travel time and is the analytically calculated travel time.The relative travel times for the transducer at for phantom pixels along the line are shown in Fig. 3 . From Fig. 3, we can see that even though the relative travel-time differences are small between the perturbed times and the analytical times, the second-order and higher travel times agree much more closely with the analytical times than the first-order travel time. We also observe that the higher-order travel times converge to the second-order travel times. Figure 4 depicts the change to relative travel times for pixels located along the line with the transducer placed at in pixel coordinates for the speed-of-sound map shown in Fig. 2. The relative travel times were calculated using where refers to a specific-order perturbative travel time and is the zeroth-order travel time. One observes that the higher-order travel-time corrections become perceivable when the speed of sound varies by 10% or more. One can also see that the travel-time perturbations seem to converge when one goes up to the fourth-order term. This suggests that one should use up to fourth-order travel-time corrections in the forward model to most accurately represent the signal.3.2.Iteratively Reconstructed ImagesIterative image reconstruction was performed using both the speed-of-sound maps. Because the variation in speed of sound is so gradual in the circularly varying speed-of-sound map, no perceivable differences in iteratively reconstructed images could be seen between using first- and second-order system matrices. The results reported here are for the elliptical speed-of-sound map. Each group of reconstructed images use the same gray scale. 3.2.1.Noiseless dataFigure 5 shows the iteratively reconstructed, noiseless images incorporating up to third-, second-, first-, and zeroth-order travel-time effects for the elliptical speed-of-sound map. The root-mean-square (rms) error in the reconstructed image intensities is given in Table 1 . Table 1rms error in noiseless reconstructed images.
For a 1% speed-of-sound variation, there is almost no difference in images reconstructed using first-order correction and higher-order corrections. For a 5% speed-of-sound variation, the rms error in the reconstructed images are about three times smaller for higher-order reconstructions than the first-order one. For a 10% speed-of-sound variation, the images reconstructed using higher-order corrections have about a five time smaller rms error than that using the first-order one. Figure 6 shows a line profile through the noiseless, reconstructed images at . One can see that the difference between higher and first-order approaches is not large for a 5% speed-of-sound variation. However, the difference is more significant for a 10% speed-of-sound variation, especially for objects lying behind the inhomogeneity. 3.2.2.Noisy dataFigure 7 shows the iteratively reconstructed, noisy images incorporating up to third-, second-, first-, and zeroth-order travel-time effects for the elliptical speed-of-sound map. Figure 8 shows a line profile through the noisy reconstructed images at . One can see that the difference between higher- and first-order approaches is again small for a 5% speed-of-sound variation. However, the difference is more significant for a 10% speed-of-sound variation especially for objects lying behind the inhomogeneity. 3.2.3.Incorrect speed-of-sound mapTo study the effect of an incorrect speed-of-sound map on higher-order GA approximation, we constructed the simulated data using the elliptical speed-of-sound map with 5.33% variation in speed of sound using the fourth-order correction. Images were iteratively reconstructing using third- and lower-order corrections assuming that the speed of variation was only 5%. Figures 9 and 10 show the reconstructed images using noiseless and noisy data, respectively. Table 2 lists the rms error in the reconstructed images using noiseless data. Table 2rms error in noiseless, reconstructed images with a 6% incorrect speed-of-sound map.
In this case, the higher-order approaches still give slightly better reconstructed images when compared to the first-order approach. However, the residual artifacts in the higher-order corrections suggest that the approach is sensitive to the accuracy of the speed-of-sound map. 4.DiscussionThe higher-order GA approximation offers better results than a first-order one for a realistic speed-of-sound map that has no sharp edges with a maximum variation of . This method is of course approximate and assumes that the speed-of-sound map is known. This method is sensitive to the accurate knowledge of the speed-of-sound map both in terms of the magnitude and the location of acoustic heterogeneity. It however does not make any assumptions about the transducer geometry. There will be certain speed-of-sound maps in practice that may result in shadow regions or caustic regions or regions where the rays get trapped. In these cases, one will not be able to correctly calculate the travel time. However, these regions will also be a problem for the GRT model using the first-order GA approximation. Hristova 10 have recently presented an application of the time-reversal algorithm to solving the photoacoustic image reconstruction problem with variable speed of sound. This is an exact approach and offers promising results but requires the detector surface to enclose the object in order to reconstruct images free from artifacts. Our higher-order GA approach could be used in those practical scenarios when the detector surface is a plane or does not completely enclose the object. It would also be useful for modeling variable speed of sound in iterative image reconstruction algorithms that are trying to model other physical factors as well in OAT. Even though our simulations were performed in 2-D, these results should also be valid in 3-D. This is because we reduced the 3-D problem to 2-D using the 2-D directivity of the transducer, which means that the acoustic propagation is still modeled by the 3-D wave equation. We just assumed the transducer is insensitive to out-of-plane waves. In this case, the isochronous, perturbed spherical surfaces become isochronous, perturbed circles in 2-D. 5.ConclusionsWe derived a higher-order geometrical acoustics approximation to the GRT model in PAT. This incorporates the first-order correction to the pressure amplitude and higher-order correction to the travel times. We found that some differences can be seen in the travel times between higher- and first-order corrections when the speed of sound varies by 5% or more. These differences in travel times translated into the reconstructed images as well. Images that were iteratively reconstructed using the higher-order corrections were more accurate than those constructed using the first-order corrections. Similar results were obtained using noisy data, although the results were much better for the noiseless data. However, the differences between first- and higher-order models were not as pronounced when using an incorrect speed-of-sound map with uncertainty in the magnitude of the acoustic inhomogeneity. We will conduct further studies using actual pressure data in an inhomogeneous medium in 3-D to quantify the effect of this higher-order GA approximation. AcknowledgmentsThe authors thank the anonymous reviewers for very helpful and valuable comments. D.M. also thanks Peter Burgholzer, Ben Cox, Kun Wang, and Brad Treeby for helpful discussions related to this work. This work was supported in part by a DoD breast cancer predoctoral fellowship (No. W81XWH-08-1-0331) to D.M. and an American Cancer Society Research Scholar award to P.L.R. ReferencesR. A. Kruger,
P. Liu,
Y. R. Fang, and
C. R. Appledorn,
“Photoacoustic ultrasound (PAUS)—reconstruction tomography,”
Med. Phys., 22 1605
–1609
(1995). https://doi.org/10.1118/1.597429 0094-2405 Google Scholar
R. S. Cobbold, Foundations of Biomedical Ultrasound, Oxford University Press, New York
(2007). Google Scholar
Y. Xu,
“Reconstruction in tomography with diffracting sources,”
Department of Biomedical Engineering, Texas A&M University,
(2003). Google Scholar
Y. Xu and
L. V. Wang,
“Effects of acoustic heterogeneity in breast thermoacoustic tomography,”
IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 50 1134
–1146
(2003). https://doi.org/10.1109/TUFFC.2003.1235325 0885-3010 Google Scholar
J. Zhang and
M. A. Anastasio,
“Reconstruction of speed-of-sound and electromagnetic absorption distributions in photoacoustic tomography,”
Proc. SPIE, 6086 608619
(2006). https://doi.org/10.1117/12.647665 0277-786X Google Scholar
T. G. M. Kolkman,
W. Steenbergen, and
T. G. van Leeuwen,
“Reflection mode photoacoustic measurement of speed of sound,”
Opt. Express, 15 3291
–3300
(2007). https://doi.org/10.1364/OE.15.003291 1094-4087 Google Scholar
S. Manohar,
R. G. H. Willemink, and
T. G. van Leeuwen,
“Speed-of-sound imaging in a photoacoustic imager,”
Proc. SPIE, 6437 64370R
(2007). https://doi.org/10.1117/12.700078 0277-786X Google Scholar
M. Agranovsky and
P. Kuchment,
“Uniqueness of reconstruction and an inversion procedure for thermoacoustic and photoacoustic tomography with variable speed of sound,”
Inverse Probl., 23 2089
–2012
(2007). https://doi.org/10.1088/0266-5611/23/5/016 0266-5611 Google Scholar
Y. Zhen and
H. Jiang,
“Three-dimensional finite-element-based photoacoustic tomography: reconstruction algorithm and simulations,”
Med. Phys., 34 538
–546
(2007). https://doi.org/10.1118/1.2409234 0094-2405 Google Scholar
Y. Hristova,
P. Kutchment, and
L. Nguyen,
“Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media,”
Inverse Probl., 24 1
–25
(2008). https://doi.org/10.1088/0266-5611/24/5/055006 0266-5611 Google Scholar
R. Snieder and
D. F. Aldridge,
“Perturbation theory for travel times,”
J. Acoust. Soc. Am., 98 1565
–1569
(1995). https://doi.org/10.1121/1.413422 0001-4966 Google Scholar
M. Xu and
L. V. Wang,
“Photoacoustic imaging in biomedicine,”
Rev. Sci. Instrum., 77 041101
(2006). https://doi.org/10.1063/1.2195024 0034-6748 Google Scholar
Y. Xu,
D. Feng, and
L. V. Wang,
“Exact frequency-domain reconstruction for thermoacoustic tomography-I: Planar geometry,”
IEEE Trans. Med. Imaging, 21 823
–828
(2002). https://doi.org/10.1109/TMI.2002.801172 0278-0062 Google Scholar
A. D. Wheelon, Electromagnetic Scintillation I. Geometrical Optics, Cambridge University Press, Cambridge, England
(2001). Google Scholar
M. Born and
E. Wolf, Principles of Optics, Cambridge University Press, Cambridge, England
(1980). Google Scholar
M. N. Wernick and
J. Aarsvold, Emission Tomography, Elsevier Academic Press, San Diego
(2004). Google Scholar
D. Modgil and
P. La Riviere,
“Implementation and comparison of reconstruction algorithms for two-dimensional optoacoustic tomography using a linear array,”
J. Biomed. Opt., 14 044023
(2009). https://doi.org/10.1117/1.3194293 1083-3668 Google Scholar
D. Modgil,
M. Anastasio,
K. Wang, and
P. La Riviere,
“Image reconstruction in photoacoustic tomography with variable speed of sound using a higher order geometrical acoustics approximation,”
Proc. SPIE, 7177 71771A
(2009). https://doi.org/10.1117/12.809001 0277-786X Google Scholar
Appendices{ label needed for app[@id='x0'] }Appendix: Third- and Fourth-Order Travel-Time CorrectionsThese results are based on the derivation of Snieder and Aldridge.11 The third-order travel-time correction is given by The fourth-order travel time is given bywhere and are the third-and fourth-order travel times. |