## 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(\mathbf{r},t)$
, that satisfies the three-dimensional inhomogeneous wave equation^{1}

## Eq. 1

$$\frac{{\partial}^{2}p(\mathbf{r},t)}{\partial {t}^{2}}-c{\left(\mathbf{r}\right)}^{2}{\nabla}^{2}p(\mathbf{r},t)=\frac{\beta}{{C}_{p}}\frac{\partial}{\partial t}H(\mathbf{r},t),$$## Eq. 2

$$p(\mathbf{r},t)=\eta \int {d}^{3}{r}^{\prime}A\left({\mathbf{r}}^{\prime}\right)\frac{\partial}{\partial t}\frac{\delta [t-(|\mathbf{r}-{\mathbf{r}}^{\prime}|\u2215c)]}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime}|},$$However, the speed of sound in tissues can vary from
$1350\phantom{\rule{0.3em}{0ex}}\mathrm{m}\u2215\mathrm{s}$
(in fat) to about
$1700\phantom{\rule{0.3em}{0ex}}\mathrm{m}\u2215\mathrm{s}$
(for skin) for ultrasound waves in the
$1\u201310\phantom{\rule{0.3em}{0ex}}\mathrm{MHz}$
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). Xu^{3} and Xu and Wang^{4} 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 Anastasio^{5} 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(\mathbf{r},t)=\eta \int {d}^{3}{r}^{\prime}A\left({\mathbf{r}}^{\prime}\right)\frac{\partial}{\partial t}\frac{\delta [t-\tau (\mathbf{r},{\mathbf{r}}^{\prime})]}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime}|},$$## Eq. 4

$$\tau (\mathbf{r},{\mathbf{r}}^{\prime})\approx {\tau}_{1}(\mathbf{r},{\mathbf{r}}^{\prime}){\int}_{{\mathbf{r}}_{\mathbf{0}}}ds\u2215c\left(s\right),$$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 Kuchment^{8} 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 Jiang^{9} 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, ${c}_{0}$ ,

where $\u03f5$ characterizes the magnitude of the perturbation.We denote the unperturbed Green’s function by ${G}_{0}$ . It satisfies the background homogeneous Helmholtz equation

## Eq. 6

$$({\nabla}^{2}+\frac{{\omega}^{2}}{{c}_{0}^{2}}){G}_{0}(\mathbf{r},{\mathbf{r}}^{\prime},\omega )=-\delta (\mathbf{r}-{\mathbf{r}}^{\prime}).$$The Green’s function that satisfies the perturbed inhomogeneous Helmholtz equation is given by $G(\mathbf{r},{\mathbf{r}}^{\prime},\omega )$

## Eq. 7

$$({\nabla}^{2}+\frac{{\omega}^{2}}{{c}^{2}\left(\mathbf{r}\right)})G(\mathbf{r},{\mathbf{r}}^{\prime},\omega )=-\delta (\mathbf{r}-{\mathbf{r}}^{\prime}).$$## 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(\mathbf{r},{\mathbf{r}}^{\prime})$ can be written as

## Eq. 8

$$G(\mathbf{r},{\mathbf{r}}^{\prime},\omega )=({g}_{0}+\frac{{g}_{1}}{ik}+\frac{{g}_{2}}{{\left(ik\right)}^{2}}+\cdots )\mathrm{exp}\left[i\omega \tau (\mathbf{r},{\mathbf{r}}^{\prime})\right],$$This can be written in the zeroth-order approximation in wavelength $\lambda $ as

## Eq. 9

$$G(\mathbf{r},{\mathbf{r}}^{\prime},\omega )={g}_{0}(\mathbf{r},{\mathbf{r}}^{\prime})\mathrm{exp}\left[i\omega \tau (\mathbf{r},{\mathbf{r}}^{\prime})\right],$$## Eq. 10

$$O\left({k}^{2}\right)\phantom{\rule{1em}{0ex}}{\left[{\nabla}_{\mathbf{r}}\tau (\mathbf{r},{\mathbf{r}}^{\prime})\right]}^{2}={c}^{-2}\left(\mathbf{r}\right),$$## Eq. 11

$$O\left({k}^{1}\right)\phantom{\rule{1em}{0ex}}2{\nabla}_{\mathbf{r}}\tau \cdot {\nabla}_{\mathbf{r}}{g}_{0}+{g}_{0}{\nabla}_{\mathbf{r}}^{2}\tau =0,$$## Eq. 12

$$O\left({k}^{0}\right)\phantom{\rule{1em}{0ex}}{\nabla}_{\mathbf{r}}^{2}{g}_{0}(\mathbf{r},{\mathbf{r}}^{\prime})=-\delta (\mathbf{r}-{\mathbf{r}}^{\prime}).$$Equation 10 is called the eikonal equation, and it implies

## Eq. 13

$$\tau (\mathbf{r},{\mathbf{r}}^{\prime})={\int}_{{\mathbf{r}}^{\prime}}^{\mathbf{r}}ds\u2215c\left(s\right),$$^{11}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 techniques^{12, 13} and Eq. 9, the pressure is given by

## Eq. 14

$$p(\mathbf{r},\omega )\simeq -i\omega \eta \int {d}^{3}{r}^{\prime}A\left({r}^{\prime}\right){g}_{0}(\mathbf{r},{\mathbf{r}}^{\prime})\mathrm{exp}\left[i\omega \tau (\mathbf{r},{\mathbf{r}}^{\prime})\right].$$Taking the Fourier transform with respect to $\omega $ , we obtain the generalized radon transform (GRT)

## Eq. 15

$$p(\mathbf{r},t)\simeq \eta \int {d}^{3}{r}^{\prime}A\left({r}^{\prime}\right){g}_{0}(\mathbf{r},{\mathbf{r}}^{\prime})\frac{d}{dt}\delta [t-\tau (\mathbf{r},{\mathbf{r}}^{\prime})].$$Equation 11 can be solved to obtain

## Eq. 16

$${g}_{0}(\mathbf{r},{\mathbf{r}}^{\prime})=C\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-\frac{1}{2}{\int}_{{\mathbf{r}}^{\prime}}^{\mathbf{r}}c\left(s\right)\left({\nabla}^{2}\tau \right)ds),$$Equation 16 can be solved to obtain^{14}

## 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 ${c}_{0}$ 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

$$\tau (\mathbf{r},{\mathbf{r}}^{\prime})={\tau}_{0}(\mathbf{r},{\mathbf{r}}^{\prime})+\u03f5{\tau}_{1}(\mathbf{r},{\mathbf{r}}^{\prime})+{\u03f5}^{2}{\tau}_{2}(\mathbf{r},{\mathbf{r}}^{\prime})+\cdots .$$Let ${\mathbf{r}}_{0}$ denote the reference ray associated with the reference eikonal ${\tau}_{0}(\mathbf{r},{\mathbf{r}}^{\prime})$ . Let the reference ray be parametrized by the variable $s$ so that ${\mathbf{r}}_{0}=s{\widehat{\mathbf{t}}}_{0}$ , where ${\widehat{\mathbf{t}}}_{0}$ is the unit vector along the reference ray. If we substitute for $\tau (\mathbf{r},{\mathbf{r}}^{\prime})$ and $c\left(\mathbf{r}\right)$ in the eikonal equation and collect terms with equal powers of the perturbation $\u03f5$ and solve, we get

## Eq. 20

$${\tau}_{0}(\mathbf{r},{\mathbf{r}}^{\prime})={\int}_{{\mathbf{r}}_{\mathbf{0}}}ds\u2215{c}_{0},$$## Eq. 21

$${\tau}_{1}(\mathbf{r},{\mathbf{r}}^{\prime})=-{\int}_{{\mathbf{r}}_{\mathbf{0}}}\frac{{c}_{1}\left(s\right)}{{c}_{0}^{2}}ds,$$## Eq. 22

$${\tau}_{2}(\mathbf{r},{\mathbf{r}}^{\prime})={\int}_{{\mathbf{r}}_{\mathbf{0}}}\frac{{c}_{0}}{2}[\frac{{c}_{1}^{2}\left(s\right)}{{c}_{0}^{4}}-{|\nabla {\tau}_{1}(\mathbf{r},{\mathbf{r}}^{\prime})|}^{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 Aldridge^{11} 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\left(\mathbf{r}\right)$ , we get from Eq. 18

## Eq. 23

$${g}_{0}(\mathbf{r},{\mathbf{r}}^{\prime})=\frac{1}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime}|}.$$Using this expression for ${g}_{0}(\mathbf{r},{\mathbf{r}}^{\prime})$ in Eq. 15, we will obtain the following expression for pressure:

## Eq. 24

$$p(\mathbf{r},t)\simeq \eta \int {d}^{3}{r}^{\prime}\frac{A\left({\mathbf{r}}^{\prime}\right)}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime}|}\frac{d}{dt}\delta [t-\tau (\mathbf{r},{\mathbf{r}}^{\prime})].$$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

where $\sigma \left(r\right)\equiv \delta \left(\mathbf{r}\right)\u2215{c}_{0}$ .GA requires that^{14}

This implies that

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

For a typical ultrasound wavelength of $\lambda =0.3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , this limit implies

For a 10% variation, i.e., $\sigma =0.1$ , we get

## 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\left(\mathbf{r}\right)$ . 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 $\Delta t$ . This implies that

But, Eq. 22 for ${\tau}_{2}$ implies that

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

and for $\Delta t=10\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ , ${c}_{0}=1.5e5\phantom{\rule{0.3em}{0ex}}\mathrm{cm}\u2215\mathrm{s}$ , Eq. 33 implies thatand for a 10% variation (i.e., $\sigma =0.1$ ), we getThus, even a $3\text{-}\mathrm{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

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.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\times 1.4\phantom{\rule{0.3em}{0ex}}\mathrm{mm})$ in this map has the variable speed of sound with respect to the background. The background speed of sound was set to $1500\phantom{\rule{0.3em}{0ex}}\mathrm{m}\u2215\mathrm{s}$ . 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\times 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.01\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ . The time-sampling interval was set to $33.33\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ . The transducers were assumed to lie on a line ( $z=0$ , $x\u2a7e0$ ) to the left of the phantom, $0.01\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ 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 as^{16}

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)=\int \int d{x}^{\prime}d{z}^{\prime}A({x}^{\prime},{z}^{\prime})\sqrt{\frac{c(x,0)}{c({x}^{\prime},{z}^{\prime})}}\delta [t-\tau (x,{x}^{\prime},{z}^{\prime})]\frac{{c}_{0}t}{\sqrt{{(x-{x}^{\prime})}^{2}+{z}^{\prime 2}}}.$$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 as^{16}

Equation 39 has the solution

where ${H}^{+}$ is the pseudo-inverse of $H$ . The iterative LS method is based on the following additive update formula:^{16}where $\widehat{f}$ is the object estimate, $n$ is the iteration number, $\alpha $ is the step size, and $\Delta 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\times 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

$${\tau}_{\text{relative},\text{order}}=\frac{({\tau}_{\text{analytical}}-{\tau}_{\text{order}})}{{\tau}_{\text{analytical}}},$$The relative travel times for the transducer at $(0,2.3\phantom{\rule{0.3em}{0ex}}\mathrm{mm})$ for phantom pixels along the line $(4.5\phantom{\rule{0.3em}{0ex}}\mathrm{mm},0)$ 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 $(0,1.0\phantom{\rule{0.3em}{0ex}}\mathrm{mm})$ with the transducer placed at $(0,2.3\phantom{\rule{0.3em}{0ex}}\mathrm{mm})$ in pixel coordinates for the speed-of-sound map shown in Fig. 2. The relative travel times were calculated using

## Eq. 43

$${\tau}_{\text{relative},\text{order}}=\frac{({\tau}_{0}-{\tau}_{\text{order}})}{{\tau}_{0}},$$## 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 .

## Table 1

rms error in noiseless reconstructed images.

Speed of sound variation(%) | Third-order | Second-order | First-order | Zeroth-order |
---|---|---|---|---|

0 | 0.07202 | 0.07202 | 0.07202 | 0.07375 |

1 | 0.07207 | 0.07207 | 0.08195 | 0.25713 |

5 | 0.07514 | 0.07516 | 0.2343 | 0.69133 |

10 | 0.10248 | 0.10353 | 0.52213 | 0.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.0\phantom{\rule{0.3em}{0ex}}\mathrm{mm})$ . 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 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.0\phantom{\rule{0.3em}{0ex}}\mathrm{mm})$ . 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 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.

## Table 2

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

Third-order | Second-order | First-order | Zeroth-order |
---|---|---|---|

0.26809 | 0.26809 | 0.31771 | 0.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
$\sim 10\mathrm{\%}$
. 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 $\sim 6\mathrm{\%}$ 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

## 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