Negativity artifacts in back-projection based photoacoustic tomography

The generation and propagation of photoacoustic signals $p(\mathbf{r}, t)$ under the condition of thermal confinement in lossless homogenous media are governed by the following photoacoustic wave equation:

$$\begin{equation} \nabla^{2} p(\mathbf{r}, t)-\frac{1}{c^{2}} \frac{\partial^{2}}{\partial t^{2}} p(\mathbf{r}, t) = -\frac{\beta}{C_{p}} \frac{\partial}{\partial t} H(\mathbf{r}, t), \end{equation} \tag{ 1 }$$

where $H(\mathbf{r}, t)$ is the heating source, representing the energy deposited in the tissue per unit volume per unit time, c is the speed of sound, β is the isobaric thermal volume expansion coefficient and C_p denotes the specific heat capacity at constant pressure. Under the condition of stress confinement, the duration of the laser pulse is much shorter than the time it takes sound to travel across the heated region. The heating function can be decomposed as $H(\mathbf{r}, t) \approx A(\mathbf{r}) \delta(t)$, where $A(\mathbf{r})$ is the energy deposited in the tissue per unit volume and δ(t) is the Dirac delta function. As a result, the initial acoustic pressure $p_{0}(\mathbf{r})$ at the position $\mathbf{r}$ can be written as:

$$\begin{equation} p_{0}(\mathbf{r}) = \Gamma A(\mathbf{r}), \end{equation} \tag{ 2 }$$

where Γ is the Gr$\ddot{\mathrm{u}}$eneisen coefficient, a dimensionless constant that represents the efficiency of the conversion of heat energy to pressure. By solving equation (1) using Green's function, the acoustic pressure $p(\mathbf{r}, t)$ recorded by a detector at location $\mathbf{r}_{\mathrm{d}}$, as depicted in figure 1, can be written as:

$$\begin{equation} p(\mathbf{r}_{\mathrm{d}}, t) = \frac{1}{4 \pi c^{2}} \frac{\partial}{\partial t} \int_{V} \frac{p_{0}\left(\mathbf{r}_{\mathrm{s}}\right)}{\left|\mathbf{r}_{\mathrm{s}}-\mathbf{r}_{\mathrm{d}}\right|} \delta\left(t-\frac{\left|\mathbf{r}_{\mathrm{s}}-\mathbf{r}_{\mathrm{d}}\right|}{c}\right) \textrm{d} V, \end{equation} \tag{ 3 }$$

where dV is a 3D volume element at $\mathbf{r}_{\mathrm{s}}$. This equation mathematically represents the spherical Radon transform and suggests that the initial acoustic pressure $p_{0}\left(\mathbf{r}_{\mathrm{s}}\right)$ in a small tissue volume at $\mathbf{r}_{\mathrm{s}}$ linearly contributes to the recorded acoustic pressure $p(\mathbf{r}_{\mathrm{d}}, t)$. However, the contribution falls as the inverse of the distance between the source and the detector, i.e. as 1/$\left|\mathbf{r}_{\mathrm{s}}-\mathbf{r}_{\mathrm{d}}\right|$.

Zoom In Zoom Out Reset image size

Image reconstruction in PAT aims to recover the initial acoustic pressure distribution $p_{0}\left(\mathbf{r}_{\mathrm{s}}\right)$ from measured photoacoustic signals $p(\mathbf{r}_{\mathrm{d}}, t)$ by reverting the spherical radon transform in equation (3). To achieve this, several analytical inversion approaches have been proposed [17–21]. One of the most well known is the BP algorithm proposed by Xu et al [18], which reconstructs the initial pressure distribution $p_{0}\left(\mathbf{r}_{\mathrm{s}}\right)$ through the following inversion formula:

$$\begin{equation} p_{0}\left(\mathbf{r}_{\mathrm{s}}\right) = \int_{\Omega} b(\mathbf{r}_{\mathrm{d}}, t) \frac{{\textrm{d}} \Omega}{\Omega}, \end{equation} \tag{ 4 }$$

where the BP term is

$$\begin{equation}b(\mathbf{r}_{\mathrm{d}}, t) = 2\left[p(\mathbf{r}_{\mathrm{d}}, t)-t \frac{\partial p(\mathbf{r}_{\mathrm{d}}, t)}{\partial t}\right] \delta\left(t-\frac{\left|\mathbf{r}_{\mathrm{s}}-\mathbf{r}_{\mathrm{d}}\right|}{c}\right),\end{equation} \tag{ 5 }$$

Ω is the entire solid angle covered by the detection surface (Ω = 2π for the infinite planar geometry and Ω = 4π for the spherical and cylindrical geometries) and dΩ is the solid angle subtended by the element dσ of the detection surface and can be written as

$$\begin{equation} \textrm{d} \Omega = \frac{\textrm{d} \sigma}{\left|\mathbf{r}_{\mathrm{s}}-\mathbf{r}_{\mathrm{d}}\right|^{2}}\left(\mathbf{n}_{\mathrm{d}} \cdot \frac{\mathbf{r}_{\mathrm{s}}-\mathbf{r}_{\mathrm{d}}}{\left|\mathbf{r}_{\mathrm{s}}-\mathbf{r}_{\mathrm{d}}\right|}\right), \end{equation} \tag{ 6 }$$

where $\mathbf{n}_{\mathrm{d}}$ denotes the unit normal vector of the detector surface pointing to the region of interest. It has been proven that the impulse response of the BP reconstruction algorithm for unbounded planar, unbounded cylindrical and closed spherical imaging geometries is a Dirac delta function and the point spread function (PSF) takes the form [18] of:

$$\begin{equation} \mathrm{PSF}\left(\mathbf{r}_{\mathrm{s}}^{^{\prime}}\right) = \delta\left(\mathbf{r}_{\mathrm{s}}^{^{\prime}}-\mathbf{r}_{\mathrm{s}}\right), \end{equation} \tag{ 7 }$$

where $\mathbf{r}_{\mathrm{s}}^{^{\prime}}$ and $\mathbf{r}_{\mathrm{s}}$ represent the locations of the point of interest and the point source, respectively. This equation suggests that, under ideal circumstances, the BP algorithm can achieve theoretically exact reconstruction of the point source without causing negativity artifacts.

The conclusion that the impulse response of the BP algorithm is a delta function is only valid under ideal circumstances, such as unbounded planar, unbounded cylindrical and closed spherical imaging geometries and infinite detector bandwidth. Such conditions, however, can never be fully met in real experiments. The deviation from these ideal conditions contaminates the PSF and eventually degrades the image quality. The presence of negativity artifacts at the edge of photoacoustic sources is a typical manifestation of image degradation.

To understand the formation process of negativity artifacts in BP-based image reconstruction, a group of simulations was performed. In this example, the photoacoustic source is a unit-intensity spherical absorber with a diameter of 5 mm located at the origin and the detection geometry is a spherical surface with a diameter of 40 mm, as shown in figure 2(a). Figure 2(b) is the cross section of figure 2(a) in the x–y plane showing the spherical absorber reduces to a disk and the spherical detection geometry becomes a circle. Figures 2(c) and (d) show the recorded approximately N-shaped signal $p(\mathbf{r}_{\mathrm{d}}, t)$ by a bandwidth-limited ((0, 4) MHz) detector located on the circle and the BP signal $b(\mathbf{r}_{\mathrm{d}}, t)$ used for image reconstruction, respectively. Figures 2(e)–(h) are the reconstructed images of the disk using 4, 16, 64, and 256 detectors uniformly distributed on the detection circle, respectively. Although the recovered disk images more and more resemble the original one as the number of detectors increases, they always contain negative values or negativity artifacts at the edge of the source, which does not disappear even with many detectors. The negativity artifacts do not arise from the negative components of the N-shaped signals but from the pair of negative wings of the BP signals in figure 2(d).

Zoom In Zoom Out Reset image size

The image reconstruction results shown in this example do not support the exact reconstruction theory presented in equation (7). The fundamental cause is the aforementioned deviation from the ideal imaging conditions. In this simulation, the detector bandwidth is not infinite but has been restricted to (0, 4) MHz. Limited bandwidth causes the detectors only to respond to a certain frequency range of the broadband photoacoustic signals and eventually introduces negativity artifacts. Moreover, the detection surface only encloses the target in a two-dimensional (2D) plane (2π radian view angle) rather than 3D space (4π steradian view angle). Since photoacoustic signals propagate in 3D space, the vast majority of the signals from the target are lost. In the following section, we quantitatively investigate the impact of the detector bandwidth and view angle on the formation of the negativity artifacts. It is worth mentioning that limited bandwidth and limited view angle will result in positive artifacts and negative artifacts simultaneously. However, negative artifacts are usually more prominent and are the focus of this work.

To investigate the impact of bandwidth on negativity artifacts, we assume that the signal detection geometry is a closed sphere and the bandwidth is the only impacting factor. The study starts with the bandwidth-limited PSF. In [25], Xu et al developed a PSF formula incorporating the effect of detector bandwidth, which can be written as:

$$\begin{equation} \mathrm{PSF}\left(\mathbf{r}_{\mathrm{s}}^{^{\prime}}\right) = \frac{1}{2 \pi^{2}} \int_{0}^{+\infty} H(k) j_{0}\left(k\left|\mathbf{r}_{\mathrm{s}}^{^{\prime}}-\mathbf{r}_{\mathrm{s}}\right|\right) k^{2} \textrm{d} k, \end{equation} \tag{ 8 }$$

where k = 2πf/c is the wavenumber, f is the ultrasound frequency, j₀ is the zeroth-order spherical Bessel function of the first kind, H(k) is the frequency response of the ultrasound detector and $\left|\mathbf{r}_{\mathrm{s}}^{^{\prime}}-\mathbf{r}_{\mathrm{s}}\right|$ is the distance from the point source $\mathbf{r}_{\mathrm{s}}$ to the point of interest $\mathbf{r}_{\mathrm{s}}^{^{\prime}}$. The PSF expressed in equation (8) applies to the planar, cylindrical and spherical detection geometries.

Specifically, if the detector has a flat frequency response over the entire frequency band, i.e. H(k) = 1 when $k \in(0,+\infty)$, equation (8) reduces to a Dirac delta function, the same as equation (7). However, this case may never happen in real experiments because ultrasound detectors cannot respond beyond a maximum frequency, which we call the upper cutoff frequency k₂. Assume that H(k) = 1 within the band (0, k₂) and H(k) = 0; otherwise, equation (8) can be simplified as:

$$\begin{equation} \mathrm{PSF}\left(\mathbf{r}_{\mathrm{s}}^{^{\prime}}\right) = \frac{k_{2}^{3}}{2 \pi} \frac{j_{1}\left(k_{2}\left|\mathbf{r}_{\mathrm{s}}^{^{\prime}}-\mathbf{r}_{\mathrm{s}}\right|\right)}{k_{2}\left|\mathbf{r}_{\mathrm{s}}^{^{\prime}}-\mathbf{r}_{\mathrm{s}}\right|}, \end{equation} \tag{ 9 }$$

where j₁ is the first-order spherical Bessel function of the first kind. Real ultrasound detectors not only have an upper cutoff frequency k₂ but also a lower cutoff frequency k₁. Assume that H(k) = 1 within the bandwidth $\left(k_{1}, k_{2}\right)$ and H(k) = 0; otherwise, equation (8) can be further written as:

$$\begin{equation} \mathrm{PSF}\left(\mathbf{r}_{\mathrm{s}}^{^{\prime}}\right) = \frac{k_{2}^{3}}{2 \pi^{2}} \frac{j_{1}\left(k_{2}\left|\mathbf{r}_{\mathrm{s}}^{^{\prime}}-\mathbf{r}_{\mathrm{s}}\right|\right)}{k_{2}\left|\mathbf{r}_{\mathrm{s}}^{^{\prime}}-\mathbf{r}_{\mathrm{s}}\right|}-\frac{k_{1}^{3}}{2 \pi^{2}} \frac{j_{1}\left(k_{1}\left|\mathbf{r}_{\mathrm{s}}^{^{\prime}}-\mathbf{r}_{\mathrm{s}}\right|\right)}{k_{1}\left|\mathbf{r}_{\mathrm{s}}^{^{\prime}}-\mathbf{r}_{\mathrm{s}}\right|}. \end{equation} \tag{ 10 }$$

Equations (9) and (10) show that the PSFs of a bandwidth-limited detector no longer converges to a Dirac delta function but are broadened as a result of the spherical Bessel function. Since the spherical Bessel functions have negative side lobes, the PSF of a bandwidth-limited detector contains negativity artifacts in principle. As an example, figures 3(a) and (b) show the PSFs of two ultrasound detectors with bandwidths of (0, 5) MHz and (4, 10) MHz, respectively, while figure 3(c) illustrates their center profiles. Ringing negativity artifacts due to limited bandwidth are visible in the PSFs, which eventually result in negative values in the final reconstructed images.

Zoom In Zoom Out Reset image size

Furthermore, we studied the impact of bandwidth on the negativity artifacts in photoacoustic images reconstructed by BP. The numerical phantom used in this example is a unit-intensity spherical absorber located at the origin and the detection geometry is a closed spherical surface, similar to the configuration in figure 2(a). In this case, the diameters of the spherical absorber and the spherical detection surface are 2 and 30 mm, respectively. Signals from the absorber are first captured by 32 768 bandwidth-limited detectors (center frequency f_c = 3 MHz) evenly distributed on the spherical detection surface in 3D space and then used for image reconstruction. Figures 4(a)–(d) show the recorded signal by a detector with a fractional bandwidth of 100%, its frequency spectrum, the BP signal and the final reconstructed image, respectively. From these results, we may draw several conclusions. First, due to limited detector bandwidth, the recorded time-domain signal (figure 4(a)) loses a significant portion of the frequency contents of the original signal (figure 4(b)) and is no longer N-shaped. Second, the distorted time-domain signal yields an inaccurate BP signal (figure 4(c)), which differs greatly from the theoretical BP signal, especially at the edge. Third, when back-projected, the inaccurate BP signal produces a distorted initial pressure image (figure 4(d)) with many negativity artifacts at the edge. As a comparison, figures 4(e)–(h) show another group of simulation with similar parameter settings except that the fractional bandwidth of the detector is expanded to 500%, i.e. bandwidth ranging from 0 to 15 MHz. Bandwidth expansion greatly eases the problem of frequency contents loss and minimizes the distortion of the recorded signal, which in turn yields approximately perfect images having significantly fewer negativity artifacts. To quantify the impact of bandwidth on the amplitude of the negativity artifacts in this example, the fractional bandwidth of the detectors was tuned from 100% to 1000%. Figure 5 shows the evaluation results, which suggests that the amplitude of the negativity artifacts decreases with the increase of the bandwidth, although the improvement is not very significant when the bandwidth is broad enough. Limited detector bandwidth is a fundamental cause of negativity artifacts in BP-based photoacoustic image reconstruction. Bandwidth expansion is helpful to mitigate the problem.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Besides detector bandwidth, the limited view angle is another fundamental factor that gives rise to negativity artifacts. Theoretically, photoacoustic signals generated from the source will not be confined in a particular plane or direction but will propagate in 3D space. Collecting signals from all view angles helps perfect image reconstruction. However, real experiments may never fulfill this requirement, which unavoidably causes negativity artifacts. To investigate the impact of detector view angle on negativity artifacts, we assume that the detector bandwidth is infinite and the view angle is the only impacting factor in the following discussion. Since the derivation of analytical PSF formulas for imaging geometries with limited covering angle is challenging, numerical simulation is adopted to study the impact of detector view angle on the negativity artifacts.

A multi-sphere phantom designed to reflect the distribution of negativity artifacts under different detection view angles is shown in figure 6(a). It consists of nine spherical absorbers with unit intensity. Eight of the absorbers with diameters uniformly varying from 1 to 2 mm are evenly distributed on a circle with a radius of 5 mm while the ninth one, having a diameter of 1.6 mm, is seated at the origin. Figures 6(b) and (c) give the cross sections of the phantom in the x–z and x–y planes, respectively.

Zoom In Zoom Out Reset image size

To study the impact of detector view angle on the characteristics of negativity artifacts, the multi-sphere phantom is imaged by three open spherical detection segments, as shown in figures 7(a)–(c). The three spherical detection segments have the same geometric radius (15 mm) and equal detector distribution density (32 768 point detectors evenly distributed over the entire spherical surface) but different surface areas or solid view angles (π, 1.5π and 2π steradian for figures 7(a)–(c), respectively). The frequency response of each detector was set to approximately full bandwidth to minimize its impact on the results. The imaging results of the x–z cross section of the phantom at view angles of π, 1.5π and 2π are shown in figures 7(d)–(f), respectively. It can be seen that negativity artifacts appear in all three x–z cross-section images as a result of the limited angle of view but their distribution and amplitude at 2π are significantly less than those at smaller view angles. This indicates that a larger view angle can produce better image quality and suppress the negativity artifacts more efficiently. The imaging results of the x–y cross sections at view angles of π, 1.5π and 2π are shown in figures 7(g)–(i). The distribution and amplitude of the negativity artifacts look different from those in the x–z cross sections but obey the same rule. The difference lies in that the artifacts, in this case, are much fewer. This is because the x–y cross section of the photoacoustic source is approximately parallel to the detection surfaces, which satisfies the visibility condition [24, 26] in PAT and can achieve better signal detection and image reconstruction. Nevertheless, the conclusion that a larger view angle improves the image quality and produces fewer artifacts still holds. Note that we do not perform simulations for view angles greater than 2π ($\Omega\gt2 \pi$) because their artifacts are exactly opposite in sign to those at 4π − Ω under the condition of exact image reconstruction (equation (7)).

Zoom In Zoom Out Reset image size

The above simulation demonstrates that the view angle of the detection geometry needs to be expanded to reduce the negativity artifacts for the spherical detection geometry. To study how detection geometry affects the distribution and amplitude of the negativity artifacts, we performed another set of simulations using planar, cylindrical and spherical geometries with finite surface areas, a practical situation in experiments. Figures 8(a)–(c) show the arrangements of the multi-sphere phantom and the three detection geometries, where the planar surface is located 12 mm away from the x–y plane, and the cylindrical surface, as well as the spherical detection surface, is located at the origin. For comparison, all three geometries have the same number of point detectors (32 768) and an approximately equal detection area (figure 8(a) 150 mm × 150 mm; figure 8(b) base radius: 35 mm, height: 100 mm; figure 8(c) radius: 42 mm). The planar and cylindrical surfaces are thus bounded while the spherical surface is closed with respect to the photoacoustic source. The frequency response of each detector was set to approximately full bandwidth to minimize its impact on the results. Figures 8(d)–(f) present the reconstruction results of the cross section of the phantom in the x–z plane while figures 8(g)–(i) show the results in the x–y plane. From these results, we may conclude as follows. First, the spherical detection geometry produces invisible negativity artifacts and the best image quality among the three, as shown in figures 8(f) and (i). This is easy to understand because the spherical detection surface is closed and strictly satisfies the exact image reconstruction condition required by BP. Second, the bounded planar detection surface yields significantly fewer negativity artifacts in the x–y plane (figure 8(g)) compared with the result produced by the cylindrical surface (figure 8(h)). This is because the photoacoustic signals emanating from the source in the x–y plane can be largely captured by the planar detector (the normal of the source intersects with the detector surface) and, therefore, the source can be reasonably reconstructed [24, 26]. However, this is not true for the cylindrical detection surface. Third, the bounded cylindrical detection surface yields fewer negativity artifacts in the x–z plane (figure 8(e)) compared with the result produced by the planar surface (figure 8(d)). This is again due to the amount of signal collected by the detection surface and can be explained using the visibility condition in [24, 26]. This example illustrates that view angle plays an important role in the generation of negativity artifacts in different detection geometries.

Zoom In Zoom Out Reset image size

In all of the numerical experiments presented so far, we used a fixed detector density, that is, 2¹⁵ (32 768) point detectors uniformly distributed over the entire spherical surface to perform the investigation. In principle, the number of detectors should be infinite because the signal-collecting surface is continuous and the detector is assumed to be a point. Decreasing the number of detectors is equivalent to decreasinge the solid view angle, which may produce negativity artifacts. Figure 9 presents a simulation reflecting the relationship between the amplitude of the negativity artifacts and the number of detectors under the same parameter settings as figure 8(c). When the number of detectors is small, significant negativity artifacts are produced, which, however, can be greatly improved by increasing the number of detectors. The number of detectors, 32 768, used in this study produces negligible negativity artifacts in this case and is, therefore, sufficient for this study.

Zoom In Zoom Out Reset image size