Translator Disclaimer
1 March 2010 Image reconstruction in photoacoustic tomography with variable speed of sound using a higher-order geometrical acoustics approximation
Author Affiliations +
Abstract
Previous research correcting for variable speed of sound in photoacoustic tomography (PAT) based on a generalized radon transform (GRT) model assumes first-order geometrical acoustics (GA) approximation. In the GRT model, the pressure is related to the optical absorption, in an acoustically inhomogeneous medium, through integration over nonspherical isochronous surfaces. Previous research based on the GRT model assumes that the path taken by acoustic rays is linear and neglects amplitude perturbations to the measured pressure. We have derived a higher-order GA expression that takes into account the first-order effect in the amplitude of the measured signal and higher-order perturbation to the travel times. The higher-order perturbation to travel time incorporates the effect of ray bending. Incorrect travel times can lead to image distortion and blurring. These corrections are expected to impact image quality and quantitative PAT. We have previously shown that travel-time corrections in 2-D suggest that perceivable differences in the isochronous surfaces can be seen when the second-order travel-time perturbations are taken into account with a 10% speed-of-sound variation. In this work, we develop iterative image reconstruction algorithms that incorporate this higher-order GA approximation assuming that the speed of sound map is known. We evaluate the effect of higher-order GA approximation on image quality and accuracy.

1.

Introduction

Optoacoustic 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, p(r,t) , that satisfies the three-dimensional inhomogeneous wave equation1

Eq. 1

2p(r,t)t2c(r)22p(r,t)=βCptH(r,t),
where c(r) is the speed of sound; H(r,t) , 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 Cp is the specific heat of the medium. The heating function can be expressed as the product of a spatially varying absorbed optical energy, A(r) , and a time-dependent optical illumination function, I(t) . If the speed of sound is constant, c(r)=c , the measured pressure signal can be related to the optical absorption function, assuming delta pulse illumination, as

Eq. 2

p(r,t)=ηd3rA(r)tδ[t(|rr|c)]4π|rr|,
where ηβCp . Equation 2 states that the time integral of acoustic pressure at a point r and time t is given by the integral of the optical absorption function over a spherical surface of radius |rr|=ct centered at r . A simple but inexact way to reconstruct A(r) is to spatially resolve the optoacoustic waves and to backproject time-integrated pressure signals over hemispheres of radii ct .

However, the speed of sound in tissues can vary from 1350ms (in fat) to about 1700ms (for skin) for ultrasound waves in the 110MHz 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

Eq. 3

p(r,t)=ηd3rA(r)tδ[tτ(r,r)]4π|rr|,
where τ(r,r) is the travel time for the sound wave to travel between points r and r . 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 r and r is given as a line integral over the slowness map

Eq. 4

τ(r,r)τ1(r,r)r0dsc(s),
where r0 is the straight line joining the points r and r , ds is the differential length element along that line, and c(s) 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.

Methods

We treat the variation in the speed of sound as a perturbation to the uniform background sound speed, c0 ,

Eq. 5

c(r)=c0+ϵc1(r),
where ϵ characterizes the magnitude of the perturbation.

We denote the unperturbed Green’s function by G0 . It satisfies the background homogeneous Helmholtz equation

Eq. 6

(2+ω2c02)G0(r,r,ω)=δ(rr).

The Green’s function that satisfies the perturbed inhomogeneous Helmholtz equation is given by G(r,r,ω)

Eq. 7

(2+ω2c2(r))G(r,r,ω)=δ(rr).

2.1.

Derivation of GA Approximation

The 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 G(r,r) can be written as

Eq. 8

G(r,r,ω)=(g0+g1ik+g2(ik)2+)exp[iωτ(r,r)],
where k2πλ is the wave number.

This can be written in the zeroth-order approximation in wavelength λ as

Eq. 9

G(r,r,ω)=g0(r,r)exp[iωτ(r,r)],
where τ is known as the eikonal function. On substituting this into Eq. 7 and equating equal powers in k we get in the limit λ0

Eq. 10

O(k2)[rτ(r,r)]2=c2(r),

Eq. 11

O(k1)2rτrg0+g0r2τ=0,

Eq. 12

O(k0)r2g0(r,r)=δ(rr).

Equation 10 is called the eikonal equation, and it implies

Eq. 13

τ(r,r)=rrdsc(s),
where s denotes the arc length along the path between the points r and r as shown in Fig. 1 . Note that the curve connecting the points r and r 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).

Fig. 1

Path traveled by sound wave between points r and r .

021308_1_006002jbo1.jpg

Using standard Green’s function techniques12, 13 and Eq. 9, the pressure is given by

Eq. 14

p(r,ω)iωηd3rA(r)g0(r,r)exp[iωτ(r,r)].

Taking the Fourier transform with respect to ω , we obtain the generalized radon transform (GRT)

Eq. 15

p(r,t)ηd3rA(r)g0(r,r)ddtδ[tτ(r,r)].

Equation 11 can be solved to obtain

Eq. 16

g0(r,r)=Cexp(12rrc(s)(2τ)ds),
where C is the constant of integration. Because this equation must also be satisfied if c(s)=c0 , we get

Eq. 17

C=14π|rr|.

Equation 16 can be solved to obtain14

Eq. 18

g0(r,r)=14π|rr|c(r)c(r).

2.2.

Higher-Order Geometrical Acoustics Approximation

We 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 c0 is perturbed as given by Eq. 5.

Following the methodology of Snieder and Aldridge,11 this leads to a perturbation in travel time as

Eq. 19

τ(r,r)=τ0(r,r)+ϵτ1(r,r)+ϵ2τ2(r,r)+.

Let r0 denote the reference ray associated with the reference eikonal τ0(r,r) . Let the reference ray be parametrized by the variable s so that r0=st̂0 , where t̂0 is the unit vector along the reference ray. If we substitute for τ(r,r) and c(r) in the eikonal equation and collect terms with equal powers of the perturbation ϵ and solve, we get

Eq. 20

τ0(r,r)=r0dsc0,

Eq. 21

τ1(r,r)=r0c1(s)c02ds,
and

Eq. 22

τ2(r,r)=r0c02[c12(s)c04|τ1(r,r)|2]ds.

Thus, 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 c(r) , we get from Eq. 18

Eq. 23

g0(r,r)=14π|rr|.

Using this expression for g0(r,r) in Eq. 15, we will obtain the following expression for pressure:

Eq. 24

p(r,t)ηd3rA(r)4π|rr|ddtδ[tτ(r,r)].

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 Approximation

2.3.1.

Requirement for validity of GA

The 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

Eq. 25

c(r)=c0+δ(r)

Eq. 26

=c0[1+σ(r)],
where σ(r)δ(r)c0 .

GA requires that14

Eq. 27

2λ|δ(r)|c01.

This implies that

Eq. 28

2λ|σ(r)|1.

If L denotes the scale of variation of speed of sound, then this implies that

σL2λ.

For a typical ultrasound wavelength of λ=0.3mm , this limit implies

Eq. 29

σ1.667Lmm1.

For a 10% variation, i.e., σ=0.1 , we get

Eq. 30

L0.0599mm.

2.3.2.

Regime where ray bending may become important

In the GRT model, the travel times are used to determine where the optoacoustic signal at a certain transducer position at a certain time t 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 A(r) . 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 Δt . This implies that

Eq. 31

τ2O(Δt).

But, Eq. 22 for τ2 implies that

Eq. 32

τ2O(c22c12c04l),
where l is the length of the inhomogeneous region.

This implies that, in order for the second-order effect to be significant, we must have

Eq. 33

σ2l2(Δt)c0
and for Δt=10ns , c0=1.5e5cms , Eq. 33 implies that

Eq. 34

σ2l0.03mm,
and for a 10% variation (i.e., σ=0.1 ), we get

Eq. 35

l3mm.

Thus, even a 3-mm inclusion of 10% variation above the background speed of sound could cause ray bending to be significant.

2.4.

Computer-Simulation Studies

We 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

Eq. 36

c(r)=c0(1+r2a2),
where a 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.

Fig. 2

Speed-of-sound maps and the original phantom used for the simulations: Slowly varying circular speed-of-sound map, elliptical speed-of-sound map, and phantom object representing absorbed optical energy

021308_1_006002jbo2.jpg

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 (2×1.4mm) in this map has the variable speed of sound with respect to the background. The background speed of sound was set to 1500ms . 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 calculations

Travel 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 50×50 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 0.01cm . The time-sampling interval was set to 33.33ns . The transducers were assumed to lie on a line ( z=0 , x0 ) to the left of the phantom, 0.01cm apart.

2.4.2.

Image reconstruction

The 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

Eq. 37

g=Hf,
where g is the time-integrated pressure data, f is the object function that represents the absorbed optical energy distribution, and H is the system matrix. The function g is defined as
g(x,t)4πηc00tp(x,t)dt.

We 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 z=0 , using the higher-order GA approximation, g(x,t) can be written as

Eq. 38

g(x,z=0,t)=dxdzA(x,z)c(x,0)c(x,z)δ[tτ(x,x,z)]c0t(xx)2+z2.

The system matrix H 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

Eq. 39

f̂=argmingHf2.

Equation 39 has the solution

Eq. 40

f̂=H+g,
where H+ is the pseudo-inverse of H . The iterative LS method is based on the following additive update formula:16

Eq. 41

f̂(n+1)=f̂(n)+αΔf(n),
where f̂ is the object estimate, n is the iteration number, α is the step size, and Δf 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 50×50 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.

Results

The 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 x=0 . 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 Times

In 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

Eq. 42

τrelative,order=(τanalyticalτorder)τanalytical,
where τorder refers to a specific-order perturbative travel time and τanalytical is the analytically calculated travel time.

The relative travel times for the transducer at (0,2.3mm) for phantom pixels along the line (4.5mm,0) are shown in Fig. 3 .

Fig. 3

Relative travel times for transducer at (0,2.3mm) for phantom pixels along the line (4.5,0mm) for a slowly varying circular speed-of-sound map: 5, 10, and 15% variation.

021308_1_006002jbo3.jpg

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 (0,1.0mm) with the transducer placed at (0,2.3mm) in pixel coordinates for the speed-of-sound map shown in Fig. 2. The relative travel times were calculated using

Eq. 43

τrelative,order=(τ0τorder)τ0,
where τorder refers to a specific-order perturbative travel time and τ0 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.

Fig. 4

Relative travel times for transducer at (0,2.3mm) for phantom pixels along the line (0,1.0mm) for a specific speed-of-sound map: 5, 10, and 15% variation.

021308_1_006002jbo4.jpg

3.2.

Iteratively Reconstructed Images

Iterative 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 data

Figure 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 .

Fig. 5

Iteratively reconstructed images using fourth-order travel-time correction in the forward model with a 0, 1, 5, and 10% speed-of-sound variation. Images were reconstructed using various-order travel-time corrections.

021308_1_006002jbo5.jpg

Table 1

rms error in noiseless reconstructed images.

Speed of sound variation(%)Third-orderSecond-orderFirst-orderZeroth-order
00.072020.072020.072020.07375
10.072070.072070.081950.25713
50.075140.075160.23430.69133
100.102480.103530.522130.74761

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 (0,2.0mm) . 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.

Fig. 6

Line profile through the phantom at (0,2.0mm) for noiseless data: 5 and 10% variation.

021308_1_006002jbo6.jpg

3.2.2.

Noisy data

Figure 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 (0,2.0mm) . 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.

Fig. 7

Iteratively reconstructed images using fourth-order travel-time correction with noise in the forward model with a 0, 1, 5, and 10% speed-of-sound variation. Images were reconstructed using various-order travel-time corrections.

021308_1_006002jbo7.jpg

Fig. 8

Line profile through the phantom at (0,2.0mm) for noisy data: 5 and 10% variation.

021308_1_006002jbo8.jpg

3.2.3.

Incorrect speed-of-sound map

To 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.

Fig. 9

Images reconstructed using incorrect speed of sound map with noiseless data with the transducer placed to the left of the object along the z -axis.

021308_1_006002jbo9.jpg

Fig. 10

Images reconstructed using incorrect speed-of-sound map with noisy data with the transducer placed to the left of the object along the z -axis.

021308_1_006002jbo10.jpg

Table 2

rms error in noiseless, reconstructed images with a 6% incorrect speed-of-sound map.

Third-orderSecond-orderFirst-orderZeroth-order
0.268090.268090.317710.82745

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.

Discussion

The 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 10% . 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.

Conclusions

We 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 6% 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.

Acknowledgments

The 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.

References

1. 

R. 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

2. 

R. S. Cobbold, Foundations of Biomedical Ultrasound, Oxford University Press, New York (2007). Google Scholar

3. 

Y. Xu, “Reconstruction in tomography with diffracting sources,” Department of Biomedical Engineering, Texas A&M University, (2003). Google Scholar

4. 

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

5. 

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

6. 

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

7. 

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

8. 

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

9. 

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

10. 

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

11. 

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

12. 

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

13. 

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

14. 

A. D. Wheelon, Electromagnetic Scintillation I. Geometrical Optics, Cambridge University Press, Cambridge, England (2001). Google Scholar

15. 

M. Born and E. Wolf, Principles of Optics, Cambridge University Press, Cambridge, England (1980). Google Scholar

16. 

M. N. Wernick and J. Aarsvold, Emission Tomography, Elsevier Academic Press, San Diego (2004). Google Scholar

17. 

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

18. 

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 Corrections

These results are based on the derivation of Snieder and Aldridge.11 The third-order travel-time correction is given by

Eq. 44

τ3(r,r)=c0r0[τ2(r,r)τ1(r,r)]ds.
The fourth-order travel time is given by

Eq. 45

τ4(r,r)=c0r0[τ1(r,r)τ3(r,r)+12|τ2(r,r)|2]ds,
where τ3 and τ4 are the third-and fourth-order travel times.

©(2010) Society of Photo-Optical Instrumentation Engineers (SPIE)
Dimple Modgil, Mark A. Anastasio, and Patrick J. La Rivière "Image reconstruction in photoacoustic tomography with variable speed of sound using a higher-order geometrical acoustics approximation," Journal of Biomedical Optics 15(2), 021308 (1 March 2010). https://doi.org/10.1117/1.3333550
Published: 1 March 2010
JOURNAL ARTICLE
9 PAGES


SHARE
Advertisement
Advertisement
Back to Top