Boundary conditions in photoacoustic tomography and image reconstruction

Lihong V. Wang; Xinmai Yang

doi:10.1117/1.2709861

1 January 2007 Boundary conditions in photoacoustic tomography and image reconstruction

Lihong V. Wang, Xinmai Yang

Author Affiliations +

Journal of Biomedical Optics, Vol. 12, Issue 1, 014027 (January 2007). https://doi.org/10.1117/1.2709861

Abstract

Recently, the field of photoacoustic tomography has experienced considerable growth. Although several commercially available pure optical imaging modalities, including confocal microscopy, two-photon microscopy, and optical coherence tomography, have been highly successful, none of these technologies can penetrate beyond ~1 mm into scattering biological tissues because all of them are based on ballistic and quasiballistic photons. Consequently, heretofore there has been a void in high-resolution optical imaging beyond this depth limit. Photoacoustic tomography has filled this void by combining high ultrasonic resolution and strong optical contrast in a single modality. However, it has been assumed in reconstruction of photoacoustic tomography until now that ultrasound propagates in a boundary-free infinite medium. We present the boundary conditions that must be considered in certain imaging configurations; the associated inverse solutions for image reconstruction are provided and validated by numerical simulation and experiment. Partial planar, cylindrical, and spherical detection configurations with a planar boundary are covered, where the boundary can be either hard or soft. Analogously to the method of images of sources, which is commonly used in forward problems, the ultrasonic detectors are imaged about the boundary to satisfy the boundary condition in the inverse problem.

1. Introduction

The motivation for photoacoustic tomography (PAT) is to combine the spatial resolution of ultrasound with the contrast of optical absorption for deep imaging of biological tissues in the diffusion regime.¹ In this regime, pure optical imaging suffers from strong light scattering by biological tissues. As a result, pure optical imaging is plagued with either shallow imaging depth or low spatial resolution. Ultrasound scattering, however, is two to three orders of magnitude weaker than optical scattering in biological tissues. Consequently, pure ultrasonic imaging provides better resolution than pure optical imaging at depths greater than approximately the optical transport mean free path $(\sim 1 mm)$ . However, because ultrasound imaging is based on the detection of mechanical properties in biological tissues, it provides only weak contrast for detecting early stage tumors. Moreover, compared with ultrasound contrast, optical absorption is highly sensitive to biochemical molecules, such as oxygenated and deoxygenated hemoglobins, and can thus provide functional imaging.

PAT overcomes the limitations of both the pure optical and the pure ultrasonic imaging modalities. In this hybrid technique, the contrast is related to the optical properties of the sample, yet the resolution is not limited by optical diffusion. The image resolution, as well as the maximum imaging depth, is scaleable with the ultrasonic bandwidth as long as the maximum depth is within the reach of diffuse photons. As the ultrasonic center frequency and bandwidth increase, the spatial resolution improves at the expense of maximum imaging depth. For example, photoacoustic signals of a 1-MHz bandwidth can provide a $\sim 1.5$ -mm spatial resolution since the speed of sound in soft tissues is $\sim 1.5 mm ∕ μ s$ , whereas photoacoustic signals of a 50-MHz bandwidth can provide a $\sim 30 - μ m$ spatial resolution. In addition, PAT images are devoid of speckle artifacts, which are conspicuously present in ultrasound images.

PAT can be implemented in either direct or reconstruction form. A focused ultrasonic transducer is used in the former, whereas an unfocused ultrasonic transducer is used and an algorithm is required in the latter. Until now, however, a boundary-free infinite acoustic medium has been assumed in all of the exact reconstruction algorithms ^{2, 3, 4, 5, 6, 7} for PAT. In this paper, the boundary conditions that must be taken into account in certain imaging configurations are presented and the associated inverse solutions for image reconstruction are provided. Partial planar, cylindrical, and spherical detection configurations with a planar boundary are considered, where the boundary can be either hard or soft. In the sections to follow, the photoacoustic equation, the forward problem in an infinite medium, and the inverse problem in an infinite medium are summarized for completeness; then, the boundary conditions for a finite medium and the inverse problem in a finite medium are presented.

2. Photoacoustic Equation

When a pressure rise is initiated by laser heating, acoustic waves develop and begin to propagate.⁸ The acoustic wave generation and propagation in an inviscid medium is described by

Eq. 1

(\nabla^{2} - \frac{1}{v_{s}^{2}} \frac{\partial^{2}}{\partial t^{2}}) p (r, t) = - \frac{B}{κ v_{s}^{2}} \frac{\partial^{2} T (r, t)}{\partial t^{2}},

where

p (r, t)

is the acoustic pressure at location

r

and time

t

,

T

is the temperature rise,

v_{s}

is the speed of sound,

κ

is the isothermal compressibility, and

β

is the thermal coefficient of volume expansion. The isothermal compressibility

κ

can be expressed as

Eq. 2

κ = \frac{C_{P}}{ρ v_{s}^{2} C_{V}},

where

ρ

is the mass density; and

C_{P}

and

C_{V}

are the specific heat capacities at constant pressure and volume, respectively. The left-hand side of Eq. 1 describes the acoustic wave propagation, while the right-hand side represents the source term.

In thermal confinement, the thermal equation becomes

Eq. 3

ρ C_{V} \frac{\partial T (r, t)}{\partial t} = H (r, t),

where

H (r, t)

is the heating function defined as the thermal energy produced in tissue by the laser light per unit time and volume. The heating function is related to the optical energy deposition

Q_{e}

by

Q_{e} = \int H d t

and to the optical fluence rate

Φ

by

H = μ_{a} Φ

if all of the absorbed optical energy is converted to heat. Substituting Eq. 3 into Eq. 1, we obtain the following less general photoacoustic equation, which is valid in thermal confinement only:

Eq. 4

(\nabla^{2} - \frac{1}{v_{s}^{2}} \frac{\partial^{2}}{\partial t^{2}}) p (r, t) = - \frac{β}{C_{P}} \frac{\partial H}{\partial t} .

Here, the source term is expressed as the first derivative of the heating function

H

with respect to time

t

. This means that constant heating does not produce a pressure wave; only a change in heating does.

3. Forward Solution in an Infinite Medium

The general photoacoustic equation [Eq. 1] can be solved by the Green’s function approach. Green’s function is defined as the response to an impulse source term represented by the product of the spatial and the temporal Dirac delta functions:

Eq. 5

(\nabla^{2} - \frac{1}{v_{s}^{2}} \frac{\partial^{2}}{\partial t^{2}}) G (r, t; r^{'}, t^{'}) = - δ (r - r^{'}) δ (t - t^{'}),

where

r^{'}

and

t^{'}

denote the source location and time, respectively. In an infinite medium, where no boundaries are involved, the Green’s function can be expressed as

Eq. 6

G (r, t; r^{'}, t^{'}) = \frac{δ [t - t^{'} - ∣ r - r^{'} ∣ ∕ v_{s}]}{4 π ∣ r - r^{'} ∣},

which represents an impulse outgoing spherical wave.

Based on the Green’s function approach, the pressure due to an arbitrary source in an infinite medium can be represented by

Eq. 7

p (r, t) = \int_{0}^{t^{+}} d t^{'} \int d r^{'} G (r, t; r^{'}, t^{'}) \frac{β}{κ v_{s}^{2}} \frac{\partial^{2} T (r^{'}, t^{'})}{\partial {t^{'}}^{2}} .

Substituting Eq. 6 into Eq. 7 leads to

Eq. 8

p (r, t) = \frac{β}{4 π κ v_{s}^{2}} \int d r^{'} \frac{1}{∣ r - r^{'} ∣} {\frac{\partial^{2} T (r^{'}, t^{'})}{\partial {t^{'}}^{2}} ∣}_{t^{'} = t - \frac{∣ r - r^{'} ∣}{v_{s}}} .

In thermal confinement, Eq. 3 holds. Thus, Eq. 8 becomes

Eq. 9

p (r, t) = \frac{β}{4 π C_{P}} \frac{\partial}{\partial t} \int d r^{'} \frac{1}{∣ r - r^{'} ∣} H (r^{'}, t - \frac{∣ r - r^{'} ∣}{v_{s}}) .

If the heating function can be decomposed into the product of the spatial and temporal components as $H (r^{'}, t^{'}) = H_{s} (r^{'}) H_{t} (t^{'})$ , then Eq. 9 can be further simplified to

Eq. 10

p (r, t) = \frac{β}{4 π C_{P}} \frac{\partial}{\partial t} \int d r^{'} \frac{H_{s} (r^{'})}{∣ r - r^{'} ∣} H_{t} (t - \frac{∣ r - r^{'} ∣}{v_{s}}) .

The physical meaning of Green’s function must be interpreted carefully. The spatial delta function simply represents a point absorber that provides a point acoustic source. However, because the source term of the photoacoustic equation is proportional to the first time derivative of the heating function in thermal confinement, or the second time derivative of the temperature in general, a temporal delta function in the source term can be translated into a step function in the heating function or a ramp function in the temperature rise. Therefore, Green’s function of the photoacoustic equation represents the response to step heating, rather than impulse heating, of a point absorber.

For impulse heating, i.e., $H_{t} (t^{'}) = δ (t^{'})$ , Eq. 10 can be used to compute the delta-heating response of an arbitrary absorbing object:

Eq. 11

p (r, t) = \frac{\partial}{\partial t} [\frac{β}{4 π C_{P}} \frac{1}{v_{s} t} \int d r^{'} H_{s} (r^{'}) δ (t - \frac{∣ r - r^{'} ∣}{v_{s}})],

where the quantity within the square brackets is the step-heating response. Since the initial pressure response on delta heating can be expressed as

p_{0} (r, t) = Γ H_{s} (r^{'})

, where the Grueneisen coefficient is defined by

Γ = v_{s}^{2} β ∕ C_{P}

, Eq. 11 can be rewritten as

Eq. 12

p (r, t) = \frac{1}{4 π v_{s}^{2}} \frac{\partial}{\partial t} [\frac{1}{v_{s} t} \int d r^{'} p_{0} (r^{'}) δ (t - \frac{∣ r - r^{'} ∣}{v_{s}})],

which can be used to derive the response to delta heating of an arbitrary object. In PAT, Eq. 12 is inverted for image reconstruction.

4. Image Reconstruction in an Infinite Medium

A universal reconstruction algorithm for an enclosing planar, cylindrical, and spherical detection surface in an acoustically homogeneous infinite medium, is summarized here.⁷ The detection surface $S$ that encloses the source $p_{0} (r^{'})$ consists of double parallel planes or a single cylindrical or spherical surface.

The following Fourier transformation pair on variable $\bar{t} = v_{s} t$ is used to convert the pressure between the time and frequency domains:

Eq. 13

{\begin{matrix} \tilde{F} (k) = \int_{- \infty}^{+ \infty} F (\bar{t}) \exp (i k \bar{t}) d \bar{t} \\ F (\bar{t}) = \frac{1}{2 π} \int_{- \infty}^{+ \infty} \tilde{F} (k) \exp (- i k \bar{t}) d k \end{matrix},

where

k = ω ∕ v_{s}

with

ω

being the angular frequency.

The initial pressure can be reconstructed by the following back-projection formula:

Eq. 14

p_{0}^{(b)} (r) = \frac{1}{π} \int_{- \infty}^{+ \infty} d k \int_{S} d S {\tilde{p}}^{*} (r_{0}, k) [{\hat{n}}_{0}^{s} ∙ \nabla_{0} {\tilde{G}}_{k}^{(out)} (r, r_{0})] .

Here,

\tilde{p} (r_{0}, k)

represents the spectrum of the pressure

p (r_{0}, \bar{t})

detected at

r_{0}

on

S

;

*

denotes complex conjugate;

\nabla_{0}

denotes the gradient over variable

r_{0}

;

{\hat{n}}_{0}^{s}

denotes the normal vector of surface

S

pointing inward; the term in the square brackets represents dipole radiation;

{\tilde{G}}_{k}^{(out)} (r, r_{0})

is a Green’s function representing a monochromatic outgoing spherical wave:

Eq. 15

{\tilde{G}}_{k}^{(out)} (r, r_{0}) = \exp (i k ∣ r^{'} - r_{0} ∣) ∕ (4 π ∣ r^{'} - r_{0} ∣) .

Equation 14 is actually a form of time reversal⁹ in the frequency domain and can be rewritten in the following form since the reconstructed pressure is real:

Eq. 16

p_{0}^{(b)} (r) = \frac{1}{π} \int_{S} d S \int_{- \infty}^{+ \infty} d k \tilde{p} (r_{0}, k) [{\hat{n}}_{0}^{s} ∙ \nabla_{0} {\tilde{G}}_{k}^{(in)} (r, r_{0})],

where

{\tilde{G}}_{k}^{(in)} (r, r_{0})

is a Green’s function corresponding to a monochromatic incoming spherical wave:

Eq. 17

{\tilde{G}}_{k}^{(in)} (r, r_{0}) = \exp (- i k ∣ r - r_{0} ∣) ∕ (4 π ∣ r - r_{0} ∣) .

5. Boundary Conditions for a Finite Medium

When an acoustic wave encounters an interface between two media of different acoustic impedances, the boundary conditions require that both the pressure and the normal medium velocity (velocity component perpendicular to the interface) be continuous, where the acoustic characteristic impedance $(Z)$ is defined by $Z = ρ v_{s}$ . The pressure continuity (boundary condition I) can be considered as a consequence of the conservation of energy and the normal velocity continuity (boundary condition II) as a consequence of the conservation of mass.

In the case where a plane wave propagating in medium 1 of impedance $Z_{1}$ is normally incident on a planar interface with medium 2 of $Z_{2}$ , the pressure reflection coefficient can be calculated by

Eq. 18

r_{n} = \frac{Z_{2} - Z_{1}}{Z_{2} + Z_{1}} .

In the case of oblique incidence with an angle of incidence

θ_{i}

, the angle of transmission

θ_{t}

is determined by Snell’s law:

Eq. 19

\sin θ_{i} ∕ v_{s 1} = \sin θ_{t} ∕ v_{s 2},

where

v_{s 1}

and

v_{s 2}

are the speeds of sound in the two media, respectively. The pressure reflection coefficient is generalized to

Eq. 20

r_{0} = \frac{Z_{2} \cos θ_{i} - Z_{1} \cos θ_{t}}{Z_{2} \cos θ_{i} + Z_{1} \cos θ_{t}} .

Of course, when

θ_{i}

is zero,

r_{0}

reduces to

r_{n}

. Two special types of boundaries are of special interest in photoacoustic tomography.

At a hard boundary, the acoustic impedance in the incident medium is much less than that in the transmitted medium, i.e., $Z_{1} ≪ Z_{2}$ , leading to $r_{a} = 1$ . In this case, boundary condition II is zero (zero normal velocity or zero normal pressure gradient), which is mathematically referred to as a homogeneous Neumann boundary condition, and the incident pressure is reflected in equal amplitude with no sign change. Acoustic wave reflection in soft tissue by an interface with metal can be modeled approximately by a hard boundary. For example, $r_{n} = 0.94$ for a water-steel interface, where the acoustic impedance of water is close to that of soft biological tissue.

At a soft boundary, the acoustic impedance in the incident medium is much greater than that in the transmitted medium, i.e., $Z_{1} ≫ Z_{2}$ , leading to $r_{a} = - 1$ . In this case, boundary condition I is zero (zero pressure), which is mathematically referred to as a homogeneous Dirichlet boundary condition, and the incident pressure is reflected in equal amplitude with a sign change. Acoustic wave reflection in soft tissue by an interface with air can be modeled well by a soft boundary. For example, $r_{n} = - 0.9994$ for a water-air interface.

6. Image Reconstruction in a Finite Medium

When boundaries are present, the boundary conditions must be incorporated into the inverse problem. First, we consider planar detection geometry with a planar hard boundary, where the detection plane is parallel to the boundary [Fig. 1a ]. In the following, we examine the forward problem for a point photoacoustic source in a finite medium with the boundary condition included. As the problem in question is linear, the conclusion will be applicable to an arbitrary photoacoustic source.

Fig. 1

Schematic for image reconstruction in planar detection geometry when a planar (a) hard or (b) soft boundary is present in PAT. Column 1 shows the original configuration, column 2 shows the added image source and image detectors that are symmetric to the original ones about the boundary, and column 3 shows the boundary-free infinite medium with the added image source and image detectors.

The method of images can be used to satisfy the boundary condition for the forward problem. For a hard boundary, the image point source has the same sign as the original point source. Once this image source is in place, the boundary condition is satisfied; as a result, the boundary can be removed. Therefore, we now have two point sources embedded in an infinite medium free of boundaries. Based on the Green’s function that provides the solution for a single source in an infinite medium, as shown in Sec. 3, the solution for these two point sources in an infinite medium can be derived as a superposition of the solutions for each of the two point sources alone.

If we image the original detectors for the inverse problem as well, we notice immediately that the photoacoustic signals that are received by the image detectors are the same as those received by the corresponding original detectors for reasons of symmetry. Therefore, in the “imaged” infinite medium, we can use the same inverse reconstruction algorithm presented in Sec. 4 to reconstruct the original point source as well as the image point source. Of course, there is no need to reconstruct the latter in practice. This conclusion applies to an arbitrary source equally.

In summary, we can use the following algorithm to reconstruct the image:

1. Image the original detectors about the planar hard boundary from $r_{0}$ to $r_{0}^{*}$ .
2. Assign the received photoacoustic signals from the original detectors to the corresponding image detectors with no sign change: $\tilde{p} (r_{0}^{*}, k) = \tilde{p} (r_{0}, k)$ . Of course, we also have $p (r_{0}^{*}, t) = p (r_{0}, t)$ .
3. Apply the reconstruction algorithm in Eq. 14 or Eq. 16 for an infinite medium to this “imaged” infinite medium, where the detection surface includes both the original and the image detectors.

We can easily extend this algorithm to the case of a planar detection configuration with a soft boundary [Fig. 1b]. In this case, the soft boundary alters the sign of the pressure on reflection. As a result, the image point source has a sign that is opposite to that of the original point source. If we also image the original detectors, the photoacoustic signals that are received by the image detectors are the same in magnitude but are opposite in sign to those received by the corresponding original detectors for reasons of antisymmetry. Therefore, a sign change is required in step 2, whereas steps 1 and 3 remain unaffected.

Second, we consider a double-semiplanar or semicylindrical detection configuration with a planar boundary (Fig. 2 ) where the detection surfaces are truncated at the boundary. For brevity, only hard boundaries are considered and only the counterpart of column 2 in Fig. 1a is drawn. The double-semiplanar detection consists of two parallel half planes that are perpendicular to the boundary, whereas the semicylindrical detection consists of a half cylinder that is perpendicular to the boundary. Once the detection surfaces are imaged about the boundary to form complete planes, or a complete cylinder, steps 2 and 3 in the preceding algorithm can be executed.

Fig. 2

Schematic for image reconstruction in double-semiplanar or semicylindrical detection geometry when a planar hard boundary is present in PAT.

Third, we similarly consider a semicylindrical or hemispherical detection configuration with a planar boundary (Fig. 3 ), where the detection surface is truncated at the boundary. Both the axis of the semicylinder and the center of the hemisphere lie on the boundary. Once the detection surfaces are imaged about the boundary to form a complete cylinder or sphere, steps 2 and 3 in the preceding algorithm can be implemented.

Fig. 3

Schematic for image reconstruction in semicylindrical or hemispherical detection geometry when a planar hard boundary is present in PAT.

Fourth, we generalize the configuration in Fig. 2 to dual parallel boundaries (Fig. 4 ), where the detection surfaces are finite in length because they are truncated by both the boundaries. Because each boundary mirrors the images, in addition to the original ones, of the other boundary alternately, an infinite series of images of the original source and detectors are formed. However, when step 3 is executed, one can judiciously truncate the series as the distal image detectors make increasingly diminishing contributions to the image in the original space. Of course, if a boundary is of the soft kind, the sign is reversed whenever a source or detector is imaged about this boundary. This configuration implies that the two parallel boundaries can turn a finite ultrasound array into an infinite array at no additional cost, thereby improving the image resolution due to the full view¹⁰ of detection and boosting the signal-noise ratio as well due to the reflected signals from the boundaries.

Fig. 4

Schematic for image reconstruction in double-finite-planar or finite-cylindrical detection geometry when two parallel planar hard boundaries are present in PAT.

Finally, we similarly generalize the configuration in Fig. 3 to dual orthogonal boundaries (Fig. 5 ), where the detection surface is only a quarter of a complete surface because it is truncated by both the boundaries. In this case, three images of the original source and detectors are formed because each boundary mirrors the images of the other boundary. This configuration can be further extended to half a detection plane that is parallel with one boundary and truncated by the other.

Fig. 5

Schematic for image reconstruction in quarter-cylindrical or quarter-spherical detection geometry when two orthogonal planar hard boundaries are present in PAT.

The planar detection configuration deserves additional comment. Unlike the cylindrical and spherical detection configurations, a nonenclosing single-detection plane is sufficient to provide an exact solution. Therefore, we can remove one of the planar detection surfaces in Figs. 2 and 4 without sacrificing the exactness of the final solution. More strikingly, we can even ignore the boundary conditions in Fig. 1 and simply reconstruct the image in the original space with only the original detectors. However, a division by 2 is in order in the algorithm of step 3 because a single detection plane has a solid angle of $2 π$ rather than $4 π$ .

The aforementioned configurations may apply in breast imaging when a breast is compressed between a rigid compression plate and an ultrasonic detection array. The configurations in Figs. 4 and 5 can potentially be used for the imaging of small animals or extremities. In these applications, the new algorithm improves the image quality significantly for the following reason. Since the original actual detection has a limited view (incomplete detection aperture), it is expected that the directly reconstructed image using the conventional reconstruction algorithm will produce significant artifacts. The inclusion of image detectors provides a closed detection surface, making the view (detection aperture) complete. The equivalent full-view data consequently produce better images.

7. Theoretical Simulation

We simulated 2-D planar imaging of a medium that has either a hard or a soft boundary ( $Z = 1$ or $- 1$ ). The buried object is a triangle that consists of three optically absorbing lines [see Fig. 6a ], which represent three typical orientations of line objects—parallel, perpendicular, and oblique to a boundary. We can write the detected signal from the object as follows:

Eq. 21

p_{m} = \frac{1}{r} \int_{L} I (t - r ∕ c_{0}) d l,

where

I (t)

denotes the impulse response of the ultrasonic detector to a point object,

L

denotes the object,

c_{0}

denotes the speed of sound in the propagation medium, and

r

denotes the distance between the detector and a small differential element

d l

on the object.

Fig. 6

(a) Image reconstructed with simulated data from a $2 π$ (full-circular) scan around the medium that has no boundary and (b) 1-D distribution of image amplitude along the dashed line in (a).

In our numerical simulations, we choose¹¹

Eq. 22

I (t) = A t^{- 3} \cos (2 π f_{0} t) \exp (- K f_{0} t),

where

A

is a constant, and

f_{0}

is the center frequency. With

K = 3.833

, this function gives a bandwidth of 80%. In our simulation,

f_{0}

is set to

6 MHz

. In addition, an additive white Gaussian noise is imposed to the signal.

Figure 6a shows an image reconstructed with simulated data from a $2 π$ (full-circular) scan in a boundary-free medium. Figure 6b shows the one-dimensional (1-D) distribution of image amplitude along the dashed line in Fig. 6a. Figure 7 shows the simulation result when a hard planar boundary $(Z = 1)$ is placed at $x = 0$ . Figure 7a shows the image reconstructed with simulated data from a $π$ (semicircular) scan using the conventional algorithm. Figure 7b shows the 1-D distribution of image amplitude along the dashed line in Fig. 7a. Compared with the result for the full-circular scan shown in Fig. 6b, an extra half cycle is present at the end of every pulse, which deteriorates the spatial resolution by causing a “shadow” in the reconstructed image. Figure 7c shows the image reconstructed with the same data using the new algorithm proposed in this paper. In both Figs. 7a and 7c, the mirrored images of the buried triangular object are clearly visible. In Fig. 7a, shadows and streaks appear owing to the limited view; in addition, the mirror image of the horizontal line is not recovered. By implementing the new reconstruction algorithm described in Sec. 6, we eliminate the preceding artifacts, as shown in Fig. 7c. The 1-D distribution of image amplitude along the dashed line in Fig. 7c is presented in Fig. 7d, which shows that the extra half cycle disappears. Figure 8 shows results similar to those in Fig. 7. The medium, however, has a soft boundary $(Z = - 1)$ . As a result, the mirror image has an opposite polarity.

Fig. 7

Simulation with a hard planar boundary placed at $x = 0$ : (a) image reconstructed with simulated data from a $π$ (semicircular) scan using the conventional algorithm, (b) 1-D distribution of image amplitude along the dashed lines in (a), (c) image reconstructed with the same data using the new algorithm, and (d) 1-D distribution of image amplitude along the dashed lines in (c).

Fig. 8

Simulation with a soft planar boundary placed at $x = 0$ : (a) image reconstructed with simulated data from a $π$ (semicircular) scan using the conventional algorithm, (b) 1-D distribution of image amplitude along the dashed lines in (a), (c) image reconstructed with the same data using the new algorithm, and (d) 1-D distribution of image amplitude along the dashed lines in (c).

It is also obvious that the new algorithm increases the image amplitude, which improves the SNR of the images. For example, the SNRs for the peaks on the left sides of Figs. 7b and 8b are 32 and $35 dB$ , respectively; by contrast, the corresponding SNRs in Figs. 7d and 8d are 35 and $38 dB$ , respectively. Therefore, a 3-dB improvement is achieved in both cases.

8. Experimental Results

In the experiment, a virtual ultrasound point detector developed in our lab was used for photoacoustic detection.¹² The virtual point detector has a center frequency around $6 MHz$ . The buried object was made of three human hairs with a $\sim 40 - μ m$ diameter and was triangularly shaped as well. A glass plate served as a hard planar boundary, whereas a piece of white plastic foam served as a soft planar boundary. The boundary was placed next to one of the right edges of the buried object. Water was used as the acoustic coupling medium. The reconstructed data are presented in Figs. 9, 10, 11 . Figure 9a shows the reconstructed image from a $2 π$ scan of the buried object in a boundary-free medium. The 1-D distribution of image amplitude along the dashed line in Fig. 9a is plotted in Fig. 9b. Like Figs. 7 and 8, Figs. 10 and 11 show the images reconstructed from data acquired under the hard and soft boundary conditions, respectively. Figure 10a shows the image reconstructed from a $π$ scan using the conventional reconstruction algorithm. The experimental results resemble the simulation results closely. Again, the image amplitude is increased by the new algorithm. For example, the SNR of the left peak is increased from 25 to $30 dB$ in the hard boundary case and from 26 to $30 dB$ in the soft boundary case.

Fig. 9

(a) Image reconstructed with experimental data from a $2 π$ (full-circular) scan around the medium that has no boundary, and (b) 1-D distribution of image amplitude along the dashed line in (a).

Fig. 10

Experiment with a hard planar boundary placed near $x = 0$ : (a) image reconstructed with simulated data from a $π$ (semicircular) scan using the conventional algorithm, (b) 1-D distribution of image amplitude along the dashed lines in (a), (c) image reconstructed with the same data using the new algorithm, and (d) 1-D distribution of image amplitude along the dashed lines in (c).

Fig. 11

Experiment with a soft planar boundary placed near $x = 0$ : (a) image reconstructed with simulated data from a $π$ (semicircular) scan using the conventional algorithm, (b) 1-D distribution of image amplitude along the dashed lines in (a), (c) image reconstructed with the same data using the new algorithm, and (d) 1-D distribution of image amplitude along the dashed lines in (c).

9. Summary

In summary, the boundary conditions that must be considered in certain imaging configurations in PAT and the associated inverse solutions for image reconstruction were presented and validated. Partial planar, cylindrical, and spherical detection configurations with a planar boundary were considered, where the boundary can be either hard or soft. Analogously to the method of images of sources used for forward problems, the ultrasonic detectors were imaged about the boundary to satisfy the boundary condition for the inverse problem. The final surfaces that include both the original and the imaged detection surfaces form a completely enclosing surface in an infinite medium. The exact algorithms that were developed for infinite media can be applied subsequently. Consequently, this new algorithm improves the quality of reconstructed images, compared with the limited-view reconstruction algorithms. The conclusions that are drawn here are equally applicable to thermoacoustic tomography, where radio frequency, rather than optical, excitation is employed.¹

Acknowledgments

We thank Geng Ku for experimental assistance. This work is supported by the National Institutes of Health grants R01 NS46214 and R01 EB000712.

References

1.

A. A. Oraevsky, and L. V. Wang, Proc. SPIE, 6086 (2006). 0277-786X Google Scholar

2.

K. P. Kostli, D. Frauchiger, J. J. Niederhauser, G. Paltauf, H. P. Weber, and M. Frenz, “Optoacoustic imaging using a three-dimensional reconstruction algorithm,” IEEE J. Sel. Top. Quantum Electron., 7 918 –923 (2001). https://doi.org/10.1109/2944.983294 1077-260X Google Scholar

3.

Y. Xu, D. Z. 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

4.

Y. Xu, M.-H. Xu, and L. V. Wang, “Exact frequency-domain reconstruction for thermoacoustic tomography. II: Cylindrical geometry,” IEEE Trans. Med. Imaging, 21 829 –833 (2002). https://doi.org/10.1109/TMI.2002.801171 0278-0062 Google Scholar

5.

M.-H. 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

6.

D. Finch, S. K. Patch, Rakesh, “Determining a function from its mean values over a family of spheres,” SIAM J. Math. Anal., 35 1213 –1240 (2003). https://doi.org/10.1137/S0036141002417814 0036-1410 Google Scholar

7.

M.-H. 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

8.

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

9.

Y. Xu and L. V. Wang, “Time reversal and its application to tomography with diffracting sources,” Phys. Rev. Lett., 92 033902 (2004). https://doi.org/10.1103/PhysRevLett.92.033902 0031-9007 Google Scholar

10.

Y. Xu, L. V. Wang, G. Ambartsoumian, and P. Kuchment, “Reconstructions in limited-view thermoacoustic tomography,” Med. Phys., 31 724 –733 (2004). https://doi.org/10.1118/1.1644531 0094-2405 Google Scholar

11.

J. L. San Emeterio and L. G. Ullate, “Diffraction impulse response of rectangular transducers,” J. Acoust. Soc. Am., 92 651 –662 (1992). https://doi.org/10.1121/1.403990 0001-4966 Google Scholar

12.

X. Yang, M. L. Li, and L. V. Wang, “Photoacoustic tomography with a virtual point detector,” Proc. SPIE, 6437 643744 (2007). 0277-786X Google Scholar

Citation Download Citation

Lihong V. Wang and Xinmai Yang "Boundary conditions in photoacoustic tomography and image reconstruction," Journal of Biomedical Optics 12(1), 014027 (1 January 2007). https://doi.org/10.1117/1.2709861

Published: 1 January 2007

Access the abstract

JOURNAL ARTICLE
10 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 38 scholarly publications and 3 patents.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Reconstruction algorithms

Sensors

Photoacoustic tomography

Image restoration

Image sensors

Acoustics

Photoacoustic spectroscopy

1.

Introduction

2.

Photoacoustic Equation

Eq. 1

Eq. 2

Eq. 3

Eq. 4

3.

Forward Solution in an Infinite Medium

Eq. 5

Eq. 6

Eq. 7

Eq. 8

Eq. 9

Eq. 10

Eq. 11

Eq. 12

4.

Image Reconstruction in an Infinite Medium

Eq. 13

Eq. 14

Eq. 15

Eq. 16

Eq. 17

5.

Boundary Conditions for a Finite Medium

Eq. 18

Eq. 19

Eq. 20

6.

Image Reconstruction in a Finite Medium

Fig. 1

Fig. 2

Fig. 3

Fig. 4

Fig. 5

7.

Theoretical Simulation

Eq. 21

Fig. 6

Eq. 22

Fig. 7

Fig. 8

8.

Experimental Results

Fig. 9

Fig. 10

Fig. 11

9.

Summary

Acknowledgments

References

Show All Keywords

Keywords/Phrases

Search In:

Publication Years