Implementation and comparison of reconstruction algorithms for two-dimensional optoacoustic tomography using a linear array

Dimple Modgil; Patrick J. La Riviere

doi:10.1117/1.3194293

1 July 2009 Implementation and comparison of reconstruction algorithms for two-dimensional optoacoustic tomography using a linear array

Dimple Modgil, Patrick J. La Riviere

Author Affiliations +

Journal of Biomedical Optics, Vol. 14, Issue 4, 044023 (July 2009). https://doi.org/10.1117/1.3194293

Abstract

Our goal is to compare and contrast various image reconstruction algorithms for optoacoustic tomography (OAT) assuming a finite linear aperture of the kind that arises when using a linear-array transducer. Because such transducers generally have tall, narrow elements, they are essentially insensitive to out-of-plane acoustic waves, and the usually 3-D OAT problem reduces to a 2-D problem. Algorithms developed for the 3-D problem may not perform optimally in 2-D. We have implemented and evaluated a number of previously described OAT algorithms, including an exact (in 3-D) Fourier-based algorithm and a synthetic-aperture-based algorithm. We have also implemented a 2-D algorithm developed by Norton for reflection mode tomography that has not, to the best of our knowledge, been applied to OAT before. Our simulation studies of resolution, contrast, noise properties, and signal detectability measures suggest that Norton's approach-based algorithm has the best contrast, resolution, and signal detectability.

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 (r, t)$ , that satisfies the three-dimensional (3-D) inhomogeneous wave equation³:

Eq. 1

\frac{\partial^{2} p (r, t)}{\partial t^{2}} - c^{2} \nabla^{2} p (r, t) = \frac{β}{C_{p}} \frac{\partial}{\partial t} H (r, t),

where

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

C_{p}

is the specific heat of the medium. The heating function can be expressed as the product of a spatially varying absorbed optical energy in the medium,

A (r)

, and a time-dependent optical illumination function,

I (t)

. The absorbed optical energy in the medium is the product of the optical absorption function and the optical fluence.

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:⁴

• 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,⁹ and universal back-projection algorithm.¹⁰ 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,¹¹ 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.¹²

2. Methods

The optoacoustic pressure signals, $p (r, t)$ , for an impulse optical illumination, are related to the absorbed optical energy $A (r)$ ¹³ as

Eq. 2

p (r, t) = η \int d^{3} r^{'} A (r^{'}) \frac{\partial}{\partial t} \frac{δ (t - \frac{| r - r^{'} |}{c})}{4 π | r - r^{'} |},

where

η = β ∕ C_{P}

. Equation 2 states that the time integral of acoustic pressure at a point

r

and time

t

is given by the integral of the absorbed optical energy function over a spherical surface of radius

| r - r^{'} | = c t

centered at

r

. A simple but inexact way to reconstruct

A (r)

is to spatially resolve the optoacoustic waves by using the speed of sound and to back-project them over hemispheres. In 2-D geometry, this reduces to spatially resolving the optoacoustic waves by using the speed of sound and back-projecting over semicircles to obtain a 2-D slice of

A (r)

. This is the method followed by the synthetic aperture algorithm.

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{η c}{4 π} \int \int d x^{'} d z^{'} A (x^{'}, z^{'}) \frac{\partial}{\partial t} δ {c t - {[{(x - x^{'})}^{2} + z^{' 2}]}^{1 ∕ 2}} .

Define

g (x, t)

as:

Eq. 4

g (x, t) \equiv \frac{4 π}{η c} \int p (x, t^{'}) d t^{'} = \int \int d x^{'} d z^{'} A (x^{'}, z^{'}) δ {c t - {[{(x - x^{'})}^{2} + z^{' 2}]}^{1 ∕ 2}} .

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

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

Eq. 5

g (x, t) = 2 r \int \int d x^{'} d z^{'} A (x^{'}, z^{'}) δ [r^{2} - {(x - x^{'})}^{2} - z^{' 2}] .

Defining new variables,

ρ = r^{2}

,

ζ = z^{2}

, and substituting in Eq. 5 yields

Eq. 6

g (x, ρ) = \sqrt{ρ} \int \int d x^{'} d ζ^{'} \frac{A (x^{'}, \sqrt{ζ^{'}})}{\sqrt{ζ^{'}}} δ [ρ - {(x - x^{'})}^{2} - ζ^{'}] .

Last, setting

Eq. 7

g^{'} (x, ρ) = \frac{g (x, \sqrt{ρ})}{\sqrt{ρ}},

and

Eq. 8

A^{'} (x, ζ) = \frac{A (x, \sqrt{ζ})}{\sqrt{ζ}},

Eq. 6 becomes

Eq. 9

g^{'} (x, ρ) = \int d x^{'} d ζ^{'} A^{'} (x, ζ) δ [ρ - ζ^{'} - {(x - x^{'})}^{2}],

which can be written as a 2-D convolution:

Eq. 10

g^{'} (x, ρ) = A^{'} (x, ρ) * * δ (ρ - x^{2}) .

This convolution relation can in principle be solved for $A^{'}$ by taking a 2-D Fourier transform on both sides of Eq. 10 with respect to $x$ and $ρ$ . 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^{'}$ , and on taking the inverse 2-D Fourier transform, one obtains $A^{'} (x, ρ)$ . The transformation of $g (x, t)$ into $g^{'} (x, ρ)$ 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:

Eq. 11

A^{'} (x, ρ) = g^{'} (x, ρ) * * R (x, ρ),

where we need to determine

R (x, ρ)

. This can be done using the method outlined in Norton’s paper¹¹ and is also derived in Sec. 6. The final result is:

Eq. 12

A (x, z) = 2 z ν_{c}^{2} \int \int g (x^{'}, r) R_{1} {ν_{c} [z^{2} - r^{2} + {(x - x^{'})}^{2}]} d r d x^{'},

where

R_{1} (u) = 4 sinc (2 u) - 2 {sinc}^{2} (u)

, and

ν_{c}

is the cutoff frequency that dictates the band-limit of the measured signals.

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 ⪢ ν_{c}^{- 1 ∕ 2}$ , this can be approximated as:

Eq. 13

A (x, z) = z ν_{c}^{3 ∕ 2} \int \int \frac{g (x^{'}, r)}{r} R_{1} ({\sqrt{ν}}_{c} {{[z^{2} + {(x - x^{'})}^{2}]}^{1 ∕ 2} - r}) d r d x^{'} .

Defining

G (x^{'}, r) \equiv g (x^{'}, r) ∕ r * R_{1} [\sqrt{ν_{c}} r]

and substituting in Eq. 13 yields:

Eq. 14

A (x, z) = z ν_{c}^{3 ∕ 2} \int G {x^{'}, {[z^{2} + {(x - x^{'})}^{2}]}^{1 ∕ 2}} d x^{'} .

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.¹⁴ 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^{'}$ followed by back-projection on circular arcs. Here, $ν_{c}^{- 1 ∕ 2}$ can be regarded as a measure of the resolution of signal $g (x, r)$ in the temporal direction since $ν_{c}$ is the bandwidth of the function $g^{'} (x, ρ)$ 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⁴ and Xu ⁵ 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:⁴

Eq. 15

A {k_{x}, k_{z} = {[{(\frac{ω}{c})}^{2} - k_{x}^{2}]}^{1 ∕ 2}} = \frac{2 c {(ω^{2} - c^{2} k_{x}^{2})}^{1 ∕ 2}}{ω} P (k_{x}, ω),

where

Eq. 16

P (k_{x}, ω) = \int \int p (x, t) \exp (- i k_{x} x) \cos (ω t) d x d t,

and

ω = c {(k_{x}^{2} + k_{z}^{2})}^{1 ∕ 2}

. Note that our notation is different from that in Ref. 4.

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}, ω)$ .
3. Scale $P (k_{x}, ω)$ via Eq. 15 to obtain $A {k_{x}, k_{z} = [{(ω ∕ c)}^{2} - k_{x}^{2}]}$ .
4. Map $A {k, k_{z} = {[{(ω ∕ c)}^{2} - k_{x}^{2}]}^{1 ∕ 2}}$ 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} ⩽ {(ω ∕ c)}^{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.⁸

Eq. 17

A (x, z) = \int g {x^{'}, \frac{{[{(x - x^{'})}^{2} + z^{2}]}^{1 ∕ 2}}{c}} d x^{'} .

So the step for obtaining

A (x, z)

is.

• For each point $(x, z)$ in the source, sum over time-integrated pressure samples corresponding to transducer positions $x^{'}$ at different times $t$ such that ${[{(x - x^{'})}^{2} + z^{2}]}^{1 ∕ 2} = c t$ .

2.4.

Details of the Simulation

All the simulations were performed using 128 transducer positions spaced $0.1 mm$ apart and 128 time samples of width $67 ns$ . Time-integrated pressure data were simulated by integrating the optoacoustic source over circles of radii $c t$ , 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¹⁵ 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 mm$ with its center placed at a distance of $2 mm$ from the transducer axis. The line of rectangles was placed at a distance of $2 mm$ from the transducer axis with each rectangle being $0.5 mm \times 0.3 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 mm$ (same as pixel width) placed at a distance of $1.0 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).¹² 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:

L N P S = \frac{1}{N} \sum_{i = 1}^{N} {| F T [noisyImage (i) - meanImage] |}^{2},

where

N

is the number of realizations.

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

L N E Q = \frac{{(L M T F)}^{2}}{L N P S} .

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.¹²

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 mm$ placed at a distance of $2 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 mm \times 0.3 mm$ placed at a distance of $2.0 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.

Fig. 1

Circular phantom images: Fourier-based; Norton-based; synthetic aperture.

Fig. 2

Images of a line of rectangles: Fourier-based; Norton-based; synthetic aperture.

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 .

Fig. 3

Zoomed-in point source images: Fourier-based; Norton-based; synthetic aperture.

Fig. 4

LIR plots: perpendicular to the transducer axis; parallel to the transducer axis.

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.

Fig. 5

Zoomed-in LMTF images (250%): Fourier-based; Norton-based; synthetic aperture.

Fig. 6

LMTF plots: perpendicular to the transducer axis; parallel to the transducer axis.

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.¹² 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.

Fig. 7

Noisy images: Fourier-based;Norton-based; synthetic aperture.

Fig. 8

Zoomed-in LNPS images (250%): Fourier-based; Norton-based; synthetic aperture.

Fig. 9

LNPS plots: perpendicular to the transducer axis; parallel to the transducer axis.

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 {mm}^{- 1}$ in the depth direction and $7 {mm}^{- 1}$ in the lateral direction.

Fig. 10

LNEQ images: Fourier-based; Norton-based; synthetic aperture.

Fig. 11

LNEQ plots: perpendicular to the transducer axis; parallel to the transducer axis.

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.

p (r, t) = η \int d^{3} r^{'} A (r^{'}) \frac{\partial}{\partial t} \frac{δ (t - \frac{| r - r^{'} |}{c})}{4 π | r - r^{'} |},

= \frac{η c}{4 π} \frac{\partial}{\partial t} \int_{| Δ r = c t |} \frac{A (r^{'})}{c t} d S^{'},

where

Δ r = | r - r^{'} |

,

d S^{'}

is the differential surface area, and the integral is over a spherical surface centered on

r

and of radius

Δ r

. This can be written as:

p (r, t) = \frac{η c}{4 π} \frac{\partial}{\partial t} \int_{| Δ r | = c t} {(c t)}^{2} \frac{A (r^{'})}{c t} d Ω^{'},

where

d Ω^{'} = \sin θ^{'} d θ^{'} d ϕ^{'}

is the solid angle. This equation in 2-D (

x - z

plane) reduces to:

Eq. 19

p (x, z, t) = \frac{η c}{4 π} \frac{\partial}{\partial t} \int_{| ρ - ρ^{'} | = c t} (c t) A (ρ^{'}) d ϕ^{'},

where

ρ = {(x^{2} + z^{2})}^{1 ∕ 2}

is the polar radial coordinate.

Using the integral property of delta functions, this can be written in polar coordinates as:

Eq. 20

p (x, z, t) = \frac{η c}{4 π} \frac{\partial}{\partial t} \int \int_{| ρ - ρ^{'} | = c t} | ρ - ρ^{'} | d (| ρ - ρ^{'} |) δ (c t - | ρ - ρ^{'} |) A (ρ^{'}) d ϕ^{'} .

This can be written equivalently in 2-D Cartesian coordinates as:

p (x, z, t) = \frac{η c}{4 π} \frac{\partial}{\partial t} \int \int d x^{'} d z^{'} A (x^{'}, z^{'}) δ {c t - {[{(x - x^{'})}^{2} + z^{' 2}]}^{1 ∕ 2}} .

Derivation of Norton-Based Algorithm

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

Eq. 21

\tilde{f} (x, ν) = \int f (x, r) \exp (i 2 π ν r) d r .

Taking the Fourier transform of Eq. 10 on both sides with respect to

ρ = r^{2}

, we get

Eq. 22

{\tilde{g}}^{'} (x, ν) = {\tilde{A}}^{'} (x, ν) * \exp (i 2 π x^{2} ν) .

On convolving both sides with $\exp (- i 2 π x^{2} ν)$ , we get

Eq. 23

{\tilde{g}}^{'} (x, ν) * \exp (- i 2 π x^{2} ν) = {\tilde{A}}^{'} (x, ν) * \exp (i 2 π x^{2} ν) * \exp (- i 2 π x^{2} ν),

where the convolution is with respect to

x

. Using the identity

\exp (i 2 π x^{2} ν) * \exp (- i 2 π x^{2} ν) = δ (2 ν x),

= [(\frac{δ (x)}{2 | ν |} + \frac{δ (ν)}{2 | x |})],

Equation 23 becomes

Eq. 24

{\tilde{g}}^{'} (x, ν) * \exp (- i 2 π x^{2} ν) = {\tilde{A}}^{'} (x, ν) * [\frac{δ (x)}{2 | ν |} + \frac{δ (ν)}{2 | x |}] = \frac{1}{2 | ν |} {\tilde{A}}^{'} (x, ν) + \frac{1}{2} δ (ν) [{\tilde{A}}^{'} (x, ν) * \frac{1}{| x |}] .

Solving for

{\tilde{A}}^{'} (x, ν)

gives

{\tilde{A}}^{'} (x, ν) = 2 | ν | {\tilde{g}}^{'} (x, ν) * \exp (- i 2 π x^{2} ν) - | ν | δ (ν) [A^{'} (x, ν) * \frac{1}{| x |}],

Using the identity

| ν | δ (ν) = 0

to eliminate the second term on the right and taking the inverse Fourier transform

(F T^{- 1})

with respect to

ν

on both sides, one finds

Eq. 25

A^{'} (x, ρ) = g^{'} (x, ρ) * * F T^{- 1} {[2 | ν | \exp (- i 2 π x^{2} ν)]}_{ρ} .

Comparing Eq. 25 with Eq. 11, we see that

26.

R (x, ρ) = F T^{- 1} {[2 | ν | \exp (- i 2 π x^{2} ν)]}_{ρ},

= 2 F T^{- 1} {[| ν |]}_{ρ + x^{2}} .

If the data

g^{'} (x, ρ)

is band-limited by a cutoff frequency

ν_{c}

, then the preceeding convolution relation is unchanged if we impose the same band-limit on

R (x, ρ)

. So one gets:

Eq. 27

R (x, ρ) = 2 F T^{- 1} {[| ν | rect (\frac{ν}{2 ν_{c}})]}_{ρ + x^{2}} .

One can write out the convolution in Eq. 25 as

A^{'} (x, ρ) = ν_{c}^{2} \int \int g^{'} (x^{'}, ρ) R_{1} {ν_{c} [ζ - ρ + {(x - x^{'})}^{2}]} d ρ d x^{'},

where

R_{1} (u) = 4 sinc (2 u) - 2 {sinc}^{2} (u)

.

Substituting for $A^{'} (x, ρ)$ , $g^{'} (x, ρ)$ and for $ζ$ and $ρ$ , one obtains

A (x, z) = 2 π ν_{c}^{2} \int \int g (x^{'}, r) R_{1} {ν_{c} [z^{2} - r^{2} + {(x - x^{'})}^{2}]} d r d x^{'},

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.

References

1.

R. A. Kruger, D. R. Reinecke, and G. A. Kruger, “Thermoacoustic computed tomography—technical considerations,” Med. Phys., 26 1832 –1837 (1999). https://doi.org/10.1118/1.598688 0094-2405 Google Scholar

2.

M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum., 77 041101-1 –041101-22 (2006). https://doi.org/10.1063/1.2195024 0034-6748 Google Scholar

3.

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

4.

K. P. Kostli and P. C. Beard, “Two-dimensional photoacoustic imaging by use of Fourier-transform image reconstruction and a detector with an anisotropic response,” Appl. Opt., 42 1899 –1908 (2003). https://doi.org/10.1364/AO.42.001899 0003-6935 Google Scholar

5.

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

6.

K. P. Kostli, M. Frenz, H. Bebie, and H. P. Weber, “Temporal backward projection of optoacoustic pressure transients using Fourier transform methods,” Phys. Med. Biol., 46 1863 –1872 (2001). https://doi.org/10.1088/0031-9155/46/7/309 0031-9155 Google Scholar

7.

C. G. A. Hoelen and F. F. M. de Mul, “Image reconstruction for photoacoustic scanning of tissue structures,” Appl. Opt., 39 5872 –5883 (2000). https://doi.org/10.1364/AO.39.005872 0003-6935 Google Scholar

8.

D. Feng, Y. Xu, G. Ku, and L. V. Wang, “Microwave-induced thermoacoustic tomography: reconstruction by synthetic aperture,” Med. Phys., 28 2427 –2431 (2001). https://doi.org/10.1118/1.1418015 0094-2405 Google Scholar

9.

C. K. Liao, M. L. Li, and P. C. Li, “Optoacoustic imaging with synthetic aperture focusing and coherence weighting,” Opt. Lett., 29 2506 –2508 (2004). https://doi.org/10.1364/OL.29.002506 0146-9592 Google Scholar

10.

M. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E, 71 016706 (2005). https://doi.org/10.1103/PhysRevE.71.016706 1063-651X Google Scholar

11.

S. J. Norton, “Reconstruction of a reflectivity field from line integrals over circular paths,” J. Acoust. Soc. Am., 67 853 –863 (1980). https://doi.org/10.1121/1.383964 0001-4966 Google Scholar

12.

H. H. Barrett and K. J. Myers, Foundations of Image Science, 1224 –1225 Wiley Interscience, Hoboken, NJ (2004). Google Scholar

13.

M. Xu and L. V. Wang, “Time domain reconstruction for thermoacoustic tomography in a spherical geometry,” IEEE Trans. Med. Imaging, 21 814 –822 (2002). https://doi.org/10.1109/TMI.2002.801176 0278-0062 Google Scholar

14.

A. C. Kak and M. Slaney, Principles of Computerized Tomography Imaging, 71 –72 IEEE Press, New York (1988). Google Scholar

15.

R. L. Siddon, “Fast calculation of the exact radiological path for a three-dimensional CT array,” Med. Phys., 12 252 –255 (1985). https://doi.org/10.1118/1.595715 0094-2405 Google Scholar

Citation Download Citation

Dimple Modgil and Patrick J. La Riviere "Implementation and comparison of reconstruction algorithms for two-dimensional optoacoustic tomography using a linear array," Journal of Biomedical Optics 14(4), 044023 (1 July 2009). https://doi.org/10.1117/1.3194293

Published: 1 July 2009

Access the abstract

JOURNAL ARTICLE
9 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 4 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Reconstruction algorithms

Signal detection

Transducers

Fourier transforms

Algorithm development

Tomography

Image restoration

1.

Introduction

Eq. 1

2.

Methods

Eq. 2

Eq. 3

Eq. 4

2.1.

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

Eq. 5

Eq. 6

Eq. 7

Eq. 8

Eq. 9

Eq. 10

Eq. 11

Eq. 12

Eq. 13

Eq. 14

2.2.

Fourier-Based Algorithm

Eq. 15

Eq. 16

2.3.

Synthetic Aperture–Based Algorithm

Eq. 17

2.4.

Details of the Simulation

3.

Results

3.1.

Phantom Images

Fig. 1

Fig. 2

3.2.

Spatial Resolution

Fig. 3

Fig. 4

Table 1

Fig. 5

Fig. 6

3.3.

Noise Texture

Fig. 7

Fig. 8

Fig. 9

3.4.

Signal Detectability/LNEQ

Fig. 10

Fig. 11

4.

Discussion

5.

Conclusions

Appendices

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

Appendix

2-D Relationship between Pressure and Absorbed Optical Energy Function

18.

Eq. 19

Eq. 20

Derivation of Norton-Based Algorithm

Eq. 21

Eq. 22

Eq. 23

Eq. 24

Eq. 25

26.

Eq. 27

Acknowledgments

References

Show All Keywords

Keywords/Phrases

Search In:

Publication Years