## 1.

## Introduction

Optoacoustic imaging is a hybrid imaging technique that has attracted a lot of attention in recent years.^{1, 2} It is based on the photoacoustic/optoacoustic effect, which refers to acoustic wave generation upon 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 subsequent contraction 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 (3-D) inhomogeneous wave equation^{3}:

## Eq. 1

$\frac{{\partial}^{2}p(\mathbf{r},t)}{\partial {t}^{2}}-{c}^{2}{\nabla}^{2}p(\mathbf{r},t)=\frac{\beta}{{C}_{p}}\frac{\partial}{\partial t}H(\mathbf{r},t),$Optoacoustic tomography (OAT) is inherently a 3-D inverse problem. The sound waves generated by a 3-D distribution of optoacoustic sources are spherical waves radiated into the volume surrounding the sources. These 3-D optoacoustic signals can be detected using isotropic ultrasound detectors arrayed on a 2-D measurement aperture, and a 3-D image can be reconstructed using these signals. However, detection of these signals using a 1-D linear array of transducers and reconstruction of a 2-D image slice is sometimes more practical and cost-effective, especially in a clinical setting. The problem can be reduced to 2-D by making one of the following assumptions, using the terminology in the paper by Kostli and Beard:^{4}

• Two-dimensional source distribution

• Two-dimensional source directivity

• Two-dimensional detector directivity.

Two-dimensional source distribution implies that the source is truly 2-D and that it lies in the detection plane. This is not a very realistic scenario for biomedical applications since there are generally 3-D sources present in the human body. The second assumption implies that even though the source is 3-D, it is constant in the third direction. This kind of source will be highly directional, and the signals received by the detector will be only from sources in the detection plane. The third assumption is relevant for highly directional detectors that are insensitive to signals coming from outside the detection plane. Thus, a 2-D cross-sectional slice of the unknown source is reconstructed, and the image reconstruction problem is reduced to 2-D. All three assumptions imply that detected acoustic signals are from only in-plane sources. However, these assumptions are not exactly equivalent. The first two assumptions impose constraints on the source geometry that may not exist. This limitation can affect the accuracy of the resulting 2-D reconstructed image that may be especially important for quantitative optoacoustic imaging. In this paper, we consider the scenario based on the third assumption, where the problem is reduced to 2-D due to detector’s directivity. This is achieved by using a linear array of anisotropic transducers that have tall, narrow elements that are essentially insensitive to out-of-plane acoustic waves.

There are several algorithms that have been proposed for 3-D image reconstruction in OAT for measurements over a planar aperture. These include the Fourier-based algorithm,^{4, 5, 6} synthetic aperture (SA) algorithm,^{7, 8} synthetic aperture plus coherent weighting algorithm,^{9} and universal back-projection algorithm.^{10} Some of these algorithms have also been applied to image reconstruction in 2-D OAT. The Fourier-based algorithm is theoretically exact in 3-D for continuously sampled data on an infinite measurement aperture but not necessarily in 2-D. This algorithm implicitly uses the preceding second assumption when applied in 2-D—namely, that the object does not vary in the third direction. Synthetic aperture and coherent weighting algorithms,^{7, 8, 9} on the other hand, are approximate reconstruction algorithms. There exists an algorithm in reflection mode tomography, proposed by Norton,^{11} that is theoretically exact in 2-D for continuously sampled data on an infinite measurement aperture. We have applied this algorithm to OAT for the first time (to the best of our knowledge). We have implemented this algorithm for a planar geometry using a linear transducer array with tall, narrow elements.

Of course, no algorithm is theoretically exact for sampled data acquired on a finite interval, so this paper compares these algorithms for that practically relevant regime by examining image contrast, resolution, and noise properties. The studies are performed using simulated optoacoustic pressure data. We go beyond the standard image quality metrics by computing noise texture measures like local noise power spectra (LNPS) and resolution measures like local modulation transfer function (LMTF). These noise and resolution measures are used to obtain the local noise equivalent quanta (LNEQ) metric that is known to predict signal detectability under certain conditions.^{12}

## 2.

## Methods

The optoacoustic pressure signals,
$p(\mathbf{r},t)$
, for an impulse optical illumination, are related to the absorbed optical energy
$A\left(\mathbf{r}\right)$
^{13} as

## 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-\frac{|\mathbf{r}-{\mathbf{r}}^{\prime}|}{c})}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime}|},$Let us consider a line of transducers along the $x$ axis (i.e., at $y=0$ , $z=0$ ). Let $A(x,z)$ represent an effective 2-D slice of the absorbed optical energy function in the half-plane ( $y=0$ , $z>0$ ). The pressure signals in 2-D reduce to (as derived in the Appendix { label needed for app[@id='x0'] })

## Eq. 3

$p(x,z=0,t)=\frac{\eta c}{4\pi}\int \int \mathrm{d}{x}^{\prime}\mathrm{d}{z}^{\prime}A({x}^{\prime},{z}^{\prime})\frac{\partial}{\partial t}\delta \{ct-{[{(x-{x}^{\prime})}^{2}+{z}^{\prime 2}]}^{1\u22152}\}.$## Eq. 4

$g(x,t)\equiv \frac{4\pi}{\eta c}\int p(x,{t}^{\prime})\mathrm{d}{t}^{\prime}=\int \int \mathrm{d}{x}^{\prime}\mathrm{d}{z}^{\prime}A({x}^{\prime},{z}^{\prime})\delta \{ct-{[{(x-{x}^{\prime})}^{2}+{z}^{\prime 2}]}^{1\u22152}\}.$We will consider three algorithms for our study: an approach based on Norton’s algorithm, the Fourier-based approach, and the synthetic aperture algorithm. We will not consider the synthetic aperture plus coherent weighting algorithm since it is nonlinear due to the presence of the coherence factor, and nonlinear algorithms cannot be meaningfully characterized using the LMTF, LNPS, and LNEQ functions.

## 2.1.

### Application of Norton’s Algorithm for Reflection Mode Tomography to OAT

In this section, we derive an exact expression for optoacoustic image reconstruction in 2-D closely following derivation of an algorithm by Norton for reflection mode tomography.^{11}

Letting $r\equiv ct$ and using the identity, $\delta (r-a)=2r\delta ({r}^{2}-{a}^{2})$ , one can write Eq. 4 as:

## Eq. 5

$g(x,t)=2r\int \int \mathrm{d}{x}^{\prime}\mathrm{d}{z}^{\prime}A({x}^{\prime},{z}^{\prime})\delta [{r}^{2}-{(x-{x}^{\prime})}^{2}-{z}^{\prime 2}].$## Eq. 6

$g(x,\rho )=\sqrt{\rho}\int \int \mathrm{d}{x}^{\prime}\mathrm{d}{\zeta}^{\prime}\frac{A({x}^{\prime},\sqrt{{\zeta}^{\prime}})}{\sqrt{{\zeta}^{\prime}}}\delta [\rho -{(x-{x}^{\prime})}^{2}-{\zeta}^{\prime}].$## Eq. 9

${g}^{\prime}(x,\rho )=\int \mathrm{d}{x}^{\prime}\mathrm{d}{\zeta}^{\prime}{A}^{\prime}(x,\zeta )\delta [\rho -{\zeta}^{\prime}-{(x-{x}^{\prime})}^{2}],$This convolution relation can in principle be solved for ${A}^{\prime}$ by taking a 2-D Fourier transform on both sides of Eq. 10 with respect to $x$ and $\rho $ . The 2-D convolution in Fourier space becomes multiplication of the 2-D Fourier transforms due to the convolution–multiplication theorem. This can then be explicitly solved for the Fourier transform of ${A}^{\prime}$ , and on taking the inverse 2-D Fourier transform, one obtains ${A}^{\prime}(x,\rho )$ . The transformation of $g(x,t)$ into ${g}^{\prime}(x,\rho )$ by Eq. 7 can be regarded as a 1-D geometric distortion of the $x-r$ plane in the $r$ direction. The double convolution relation, Eq. 10, implies that the time-integrated pressure signal is equivalent to a linear, space-variant mapping of the absorbed optical energy from points in the $x-z$ object plane to hyperbolas in the $x-r$ measurement plane (see Ref. 11, Sec. 2, for more details).

Another approach that is more direct to solving Eq. 10 is to seek a solution such that:

where we need to determine $R(x,\rho )$ . This can be done using the method outlined in Norton’s paper^{11}and is also derived in Sec. 6. The final result is:

## Eq. 12

$A(x,z)=2z{\nu}_{c}^{2}\int \int g({x}^{\prime},r){R}_{1}\left\{{\nu}_{c}[{z}^{2}-{r}^{2}+{(x-{x}^{\prime})}^{2}]\right\}\mathrm{d}r\mathrm{d}{x}^{\prime},$The preceding relation is an exact equation relating the absorbed optical energy in the medium to a filtered back-projection of time-integrated pressure signals.

In the case when $z\u2aa2{\nu}_{c}^{-1\u22152}$ , this can be approximated as:

## Eq. 13

$A(x,z)=z{\nu}_{c}^{3\u22152}\int \int \frac{g({x}^{\prime},r)}{r}{R}_{1}\left({\sqrt{\nu}}_{c}\{{[{z}^{2}+{(x-{x}^{\prime})}^{2}]}^{1\u22152}-r\}\right)\mathrm{d}r\mathrm{d}{x}^{\prime}.$## Eq. 14

$A(x,z)=z{\nu}_{c}^{3\u22152}\int G\{{x}^{\prime},{[{z}^{2}+{(x-{x}^{\prime})}^{2}]}^{1\u22152}\}\mathrm{d}{x}^{\prime}.$Equation 13 is similar to the filtered back-projection (FBP) expression used for image reconstruction in computer tomography (CT). The function
${R}_{1}$
is the same as the Fourier transform of the truncated ramp filter used in the FBP expression.^{14} Eq. 14 can be seen as back-projection of the filtered function
$G$
. So this algorithm is equivalent to 1-D filtration at each transducer position
${x}^{\prime}$
followed by back-projection on circular arcs. Here,
${\nu}_{c}^{-1\u22152}$
can be regarded as a measure of the resolution of signal
$g(x,r)$
in the temporal direction since
${\nu}_{c}$
is the bandwidth of the function
${g}^{\prime}(x,\rho )$
with respect to the square of the temporal variable
$r$
. We found that the exact Norton’s algorithm was extremely sensitive to the choice of cutoff frequency and did not give us good results. Hence, we implemented the approximate Eq. 13 as the Norton-based algorithm.

So $A(x,z)$ can be obtained via the following steps in the approximate Norton-based algorithm:

1. Convolve the 2-D time-integrated pressure signal $g$ with 1-D filter ${R}_{1}$ with respect to the temporal variable $r$ .

2. Map the result onto a circular grid for a given $(x,z)$ .

3. Sum the resulting expression over all the transducers.

4. Multiply the result by the distance from the transducer axis, $z$ , and other constant factors.

## 2.2.

### Fourier-Based Algorithm

This algorithm has been derived by Kostli and Beard^{4} and Xu
^{5} It relates the Fourier transform of the absorbed optical energy function to the Fourier transform of the measured optoacoustic pressures. The relation in 2-D is given by:^{4}

## Eq. 15

$A\{{k}_{x},{k}_{z}={[{\left(\frac{\omega}{c}\right)}^{2}-{k}_{x}^{2}]}^{1\u22152}\}=\frac{2c{({\omega}^{2}-{c}^{2}{k}_{x}^{2})}^{1\u22152}}{\omega}P({k}_{x},\omega ),$## Eq. 16

$P({k}_{x},\omega )=\int \int p(x,t)\mathrm{exp}(-i{k}_{x}x)\mathrm{cos}\left(\omega t\right)\mathrm{d}x\mathrm{d}t,$Here, $A(x,z)$ can be obtained via the following steps:

1. Take the real part Fourier transform of pressure, $p(x,t)$ with respect to time.

2. Take the Fourier transform of the result with respect to $x$ . This gives us $P({k}_{x},\omega )$ .

3. Scale $P({k}_{x},\omega )$ via Eq. 15 to obtain $A\{{k}_{x},{k}_{z}=[{(\omega \u2215c)}^{2}-{k}_{x}^{2}]\}$ .

4. Map $A\{k,{k}_{z}={[{(\omega \u2215c)}^{2}-{k}_{x}^{2}]}^{1\u22152}\}$ to a Cartesian grid via interpolation to obtain $A({k}_{x},{k}_{z})$ . We employed bilinear interpolation. Note that only values for which ${k}_{x}^{2}\u2a7d{(\omega \u2215c)}^{2}$ are used in the mapping. Higher values would produce imaginary values of ${k}_{z}$ . (Physically, they correspond to rapidly decaying evanescent wave components.)

5. Inverse Fourier transform $A({k}_{x},{k}_{z})$ to obtain $A(x,z)$ .

## 2.3.

### Synthetic Aperture–Based Algorithm

This algorithm relates the signal intensity at each image point
$A(x,z)$
to the sum of signals from the transducer at different positions delayed with the time it takes the signal to travel from the transducer position to the point.^{8}

## Eq. 17

$A(x,z)=\int g\{{x}^{\prime},\frac{{[{(x-{x}^{\prime})}^{2}+{z}^{2}]}^{1\u22152}}{c}\}\mathrm{d}{x}^{\prime}.$## 2.4.

### Details of the Simulation

All the simulations were performed using 128 transducer positions spaced $0.1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ apart and 128 time samples of width $67\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ . Time-integrated pressure data were simulated by integrating the optoacoustic source over circles of radii $ct$ , with the transducer placed at the center of these circles. Simulated 2-D pressure data were used for the Fourier-based algorithm, while time-integrated simulated pressure data were used for the synthetic aperture and Norton-based algorithms. In order to focus simply on the acquisition geometry and the inherent differences between the algorithms, we did not simulate a low-pass or bandpass transducer response.

Simulated pressure data were generated for two different phantoms—a circular phantom and a phantom consisting of a line of rectangles. An exact analytical expression was used to generate the time-integrated pressure data for the circular and point source phantoms A discrete numerical method was used to generate the time-integrated pressure data for the line phantom. This method is based on an implementation of a variation of Siddon’s algorithm^{15} for computing the intersection lengths of an arc specified by the coordinates of a source and receiver with a pixel. The circular phantom was of radius
$1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
with its center placed at a distance of
$2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
from the transducer axis. The line of rectangles was placed at a distance of
$2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
from the transducer axis with each rectangle being
$0.5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\times 0.3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
wide. The images were constructed on a
$128\times 128$
grid.

To study the resolution, simulated pressure data were generated for a point source of size
$0.1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
(same as pixel width) placed at a distance of
$1.0\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
from the transducer axis. A zoomed-in image of a point source of was reconstructed using the three algorithms with a zoom factor of 10. The images were reconstructed on a
$64\times 64$
grid. We used these point source images to compute a local impulse response (LIR) function. The LIR is a generalization of the point-spread function applicable when the image acquisition and reconstruction processes are not shift-invariant, as is the case here. We then computed the 2-D Fourier transforms of the LIRs to obtain what Barrett and Myers have called a local modulation transfer function (LMTF).^{12} A standard modulation transfer function (MTF) is the absolute value of the Fourier transform of the system’s point-spread function and predicts the ratio of output modulation to input modulation as a function of spatial frequency. The localized modulation transfer function (LMTF) is the generalization of MTF to linear, shift-variant systems. Higher and broader LMTF indicates better resolution properties.

Random Gaussian noise with mean 0 and a standard deviation of 1.0 was used for noisy pressure signals for the noise studies. Noisy images were constructed on a $64\times 64$ grid using a zoom factor of 10. The noise studies were performed for 500 realizations. LNPS is a generalization of the noise power spectra (NPS) concept that can be used for linear systems without the assumption of shift invariance that does not hold for finite transducer aperture systems in OAT. LNPS was computed by first generating a set of 500 realizations of reconstructed images for each algorithm corresponding to pure Gaussian noise pressure. For each set of these reconstructed images, the mean image was computed. The mean image was then subtracted from the other 500 images. LNPS was then obtained by taking the average squared modulus of the Fourier transform (FT) of the subtracted images:

LNEQ is defined as the ratio of the square of the LMTF to the LNPS:

LNEQ is a kind of frequency-dependent signal-to-noise ratio generalized to linear, shift-variant systems. Higher LNEQ implies higher signal detectability performance for the so-called ideal observer in the task when both signal and background are known exactly.^{12}

## 3.

## Results

## 3.1.

### Phantom Images

In all the phantom images shown in this section, the transducer axis is along the bottom of the images. Figure 1 shows non-zoomed-in images produced by the three algorithms for a circular phantom of radius $1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ placed at a distance of $2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ from the transducer axis. The image reconstructed via the Fourier based algorithm has sharp edges, but it is nonuniform and has lower contrast. The synthetic aperture image is quite blurred. The image reconstructed via Norton’s approach is quite sharp and fairly uniform. Note that we observe two images in the Fourier-based algorithm due to the implicit symmetry assumption in the reconstruction. Figure 2 shows non-zoomed-in images of a line of small rectangles of size $0.5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\times 0.3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ placed at a distance of $2.0\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ from the transducer axis. The circular arc artifacts are more visible in the synthetic aperture algorithm than the Norton-based algorithm due to the additional filtration step that is performed in the Norton-based algorithm. In addition, the rectangles themselves are much sharper and more filled-in in the Norton-based algorithm compared to the other two algorithms.

## 3.2.

### Spatial Resolution

The images of zoomed-in point sources are shown in Fig. 3, where the transducer axis is along the bottom of the images. The LIR plots for the three algorithms are shown in Fig. 4 . These show that the Fourier-based algorithm shows the best resolution perpendicular to the transducer array, while Norton’s algorithm shows the best resolution parallel to the transducer array. The main difference between Norton’s algorithm and the SA algorithm is filtering. This results in a much narrower LIR for Norton’s algorithm than for the SA algorithm. In general, the lateral resolution for the Norton-based and SA algorithms is much better than the depth resolution (perpendicular to the transducer axis). The full width at half maxima (FWHM) results are shown in Table 1 .

## Table 1

FWHM values for a point source of size 0.1mm with pixel width=0.1mm .

FWHM | Fourier-based | Norton-based | Synthetic aperture |
---|---|---|---|

Depth (mm) | 0.154 | 0.200 | 0.471 |

Lateral (mm) | 0.161 | 0.151 | 0.189 |

Figure 5 shows the LMTF images for the three algorithms. The LMTF images exhibit an asymmetry due to the finite transducer length. The reciprocal relationship between LIR and LMTF is exhibited in these images. LMTF for SA algorithm is the narrowest, since LIR for SA is the broadest. Figure 6 shows the LMTF plots. LMTF for the Fourier-based algorithm is best in the direction perpendicular to the transducer axis, especially for smaller frequencies. Norton’s algorithm produces the best LMTF profile in the lateral direction, which is expected since it had the smallest lateral resolution.

## 3.3.

### Noise Texture

The noise texture images are simply images of a uniform background reconstructed from pressure signals to which zero-mean Gaussian noise has been added. Figure 7 shows the noise texture in the unzoomed images reconstructed via the three algorithms. The noise in the Fourier-and Norton-based algorithms seems uniformly speckled, while the smeared out noise texture in the SA algorithm seems to exhibit some long-range correlations. Such “blobbiness” in the noise can impede detectability of signals of size comparable to the blob size. LNPS describes the frequency content of the background texture in the region surrounding the location of the signal in the image.^{12} Figure 8 shows the zoomed-in (250%) images of LNPS for the three algorithms. LNPS images for the Norton-based and synthetic aperture algorithms are fairly symmetric. But this is not the case for the Fourier-based algorithm. Note that the input to the Fourier-based algorithm was Gaussian noise pressure, while the input to the other two algorithms was time-integrated Gaussian noise pressure, which introduces noise correlations that can affect the form of the LNPS. Figure 9 shows the LNPS plots in the two directions. The plot for the SA algorithm was omitted because it was several orders of magnitude higher than the other two algorithms and could not be fit into the same plot. The shapes for LNPS are very different for the Norton-based and Fourier-based algorithms.

## 3.4.

### Signal Detectability/LNEQ

Figure 10 shows the zoomed-in (250%) images of LNEQ for the three algorithms. Figure 11 shows the LNEQ plots. In general, LNEQ images are somewhat difficult to interpret visually, as they represent a kind of generalized frequency–domain detectability transfer function, but bright areas correspond to frequency components that are more likely to be detected against a background of correlated noise, as captured by the LNPS. To calculate the so-called ideal observer signal-to-noise ratio for a small low-contrast signal, one would calculate the squared magnitude of the Fourier transform of the signal, multiply it by the LNEQ shown in the images, and integrate over all frequencies. The plots in Fig. 11 give a better sense of the relative magnitude of the LNEQ for the different algorithms. Higher values are better, and the Norton-based algorithm produces the highest LNEQ in both directions followed by the Fourier-based algorithm. This indicates superior ideal observer signal detectability in images reconstructed by the use of the Norton-based algorithm. For all algorithms, the LNEQ becomes small at $5\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ in the depth direction and $7\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ in the lateral direction.

## 4.

## Discussion

The 2-D geometry that we consider in this paper reduces the inherently 3-D optoacoustic image reconstruction to 2-D due to transducer’s 2-D directivity. This results in the measured 2-D optoacoustic signal being related to the time derivative of the integral of the absorbed optical energy over circular arcs centered at the transducer as given by Eq. 3. The three algorithms considered here for evaluation have some intrinsic differences. The SA algorithm is approximate to start with and involves only the delay and sum technique based on back-projection over circular arcs. The Norton-based algorithm improves on the SA algorithm by providing a filtration step that filters the time-integrated pressure data before using the back- projection over circular arcs technique. The Fourier-based algorithm is not exact in the 2-D geometry that we are considering since it assumes that the optoacoustic source is constant in the third dimension. These intrinsic differences in the algorithms explain the differences in sharpness and quality of the reconstructed images.

Norton-based and SA algorithms use time-integrated pressure as input signal, while the Fourier-based algorithm uses pressure signal as input. Integration of noisy data introduces noise correlations that can affect the LNPS. This was reflected in the blobbiness of the noise texture images of the SA algorithm. The additional filtration step in the Norton-based algorithm aids in the removal of such blobbiness and gives a much better noise texture. The Fourier-based algorithm does not show such blobbiness in noise texture.

## 5.

## Conclusions

We explored three different ways in which a 2-D image can be reconstructed in OAT. We implemented and evaluated three algorithms for 2-D image reconstruction in OAT Fourier-based, Norton-based, and synthetic aperture algorithms. We found that the 2-D Fourier-based algorithm offers better resolution and LMTF in the depth direction, while Norton’s algorithm offers the best lateral resolution. However, we found that in reconstructions of phantoms, the images produced by the Norton-based algorithm looked the sharpest and most uniform. The LNEQ data suggests that the Norton-based algorithm has the best signal detectability.

## Appendices

## { label needed for app[@id='x0'] }

### Appendix

#### 2-D Relationship between Pressure and Absorbed Optical Energy Function

Starting with Eq. 2,

## 18.

## Eq. 19

$p(x,z,t)=\frac{\eta c}{4\pi}\frac{\partial}{\partial t}{\int}_{|\rho -{\rho}^{\prime}|=ct}\left(ct\right)A\left({\rho}^{\prime}\right)\mathrm{d}{\varphi}^{\prime},$Using the integral property of delta functions, this can be written in polar coordinates as:

## Eq. 20

$p(x,z,t)=\frac{\eta c}{4\pi}\frac{\partial}{\partial t}\int {\int}_{|\rho -{\rho}^{\prime}|=ct}|\rho -{\rho}^{\prime}|d\left(|\rho -{\rho}^{\prime}|\right)\delta (ct-|\rho -{\rho}^{\prime}|)A\left({\rho}^{\prime}\right)\mathrm{d}{\varphi}^{\prime}.$This can be written equivalently in 2-D Cartesian coordinates as:

#### Derivation of Norton-Based Algorithm

In this section, we shall derive Eq. 12 following the method detailed in Norton’s paper.^{11} Define the Fourier transform with respect to
$r$
as:

## Eq. 22

${\stackrel{\u0303}{g}}^{\prime}(x,\nu )={\stackrel{\u0303}{A}}^{\prime}(x,\nu )\ast \mathrm{exp}\left(i2\pi {x}^{2}\nu \right).$On convolving both sides with $\mathrm{exp}(-i2\pi {x}^{2}\nu )$ , we get

## Eq. 23

${\stackrel{\u0303}{g}}^{\prime}(x,\nu )\ast \mathrm{exp}(-i2\pi {x}^{2}\nu )={\stackrel{\u0303}{A}}^{\prime}(x,\nu )\ast \mathrm{exp}\left(i2\pi {x}^{2}\nu \right)\ast \mathrm{exp}(-i2\pi {x}^{2}\nu ),$## Eq. 24

${\stackrel{\u0303}{g}}^{\prime}(x,\nu )\ast \mathrm{exp}(-i2\pi {x}^{2}\nu )={\stackrel{\u0303}{A}}^{\prime}(x,\nu )\ast [\frac{\delta \left(x\right)}{2\left|\nu \right|}+\frac{\delta \left(\nu \right)}{2\left|x\right|}]=\frac{1}{2\left|\nu \right|}{\stackrel{\u0303}{A}}^{\prime}(x,\nu )+\frac{1}{2}\delta \left(\nu \right)[{\stackrel{\u0303}{A}}^{\prime}(x,\nu )\ast \frac{1}{\left|x\right|}].$## Eq. 25

${A}^{\prime}(x,\rho )={g}^{\prime}(x,\rho )\ast \ast F{T}^{-1}{\left[2\left|\nu \right|\mathrm{exp}(-i2\pi {x}^{2}\nu )\right]}_{\rho}.$Comparing Eq. 25 with Eq. 11, we see that

## 26.

## Eq. 27

$R(x,\rho )=2F{T}^{-1}{\left[\left|\nu \right|\mathrm{rect}\left(\frac{\nu}{2{\nu}_{c}}\right)\right]}_{\rho +{x}^{2}}.$Substituting for ${A}^{\prime}(x,\rho )$ , ${g}^{\prime}(x,\rho )$ and for $\zeta $ and $\rho $ , one obtains

which is the desired result.## Acknowledgments

This work was supported in part by the University of Chicago’s SPORE grant for breast cancer research, a DOD breast cancer predoctoral fellowship (W81XWH-08-1-0331), and an American Cancer Society Research Scholar award. D.M. would like to thank Phillip A. Vargas for implementing the forward model in 2-D for OAT and for helpful discussions related to it.