1.

Introduction

In photoacoustic (also known as thermoacoustic or optoacoustic) imaging, a sample is illuminated by a short pulse of electromagnetic radiation, such as a laser or a microwave pulse. The electromagnetic radiation is absorbed to different extents in different regions of the sample, and the resulting weak heating and thermal expansion launch an acoustic wave. This acoustic pressure signal is measured by detectors outside the sample and can be used for reconstruction of the distribution of the absorbed electromagnetic energy inside the sample. Therefore, photoacoustic imaging is a technique that combines the advantages of optical methods (good optical absorption contrast between different types of tissues) and ultrasound (high spatial resolution due to low ultrasound scattering). Utilizing acoustic detectors with large aperture that detect signals simultaneously from all areas of the sample requires tomographic methods for image reconstruction. This is known as photoacoustic tomography (good overviews are provided in Refs. 1, 2, and 3).

Mathematically, photoacoustic wave propagation can be modeled as an initial value problem for the three-dimensional (3-D) wave equation. Within this mathematical framework, image reconstruction algorithms have been developed for detectors so small that they can be modeled as points. Physical implementations of detector arrays use piezoelectric sensors, which have the drawback of decreasing signal-to-noise ratio with decreasing size. Therefore, piezoelectric detector elements are used with a size in the range of several millimeters, leading to some blurring artifacts if reconstruction algorithms for point sensors are applied to measured data. An alternative is to use linear or planar detectors that integrate the acoustic pressure over one or two spatial dimensions. Because of pressure integration over a finite area with known shape, noise is reduced and signal-to-noise ratio is enhanced. Reconstruction algorithms have been adapted to handle such special kinds of detectors. For integrating line detectors, the inversion of the 3-D wave equation is usually split into the following two steps:^4,5 first, for a given orientation of the line detectors, a two-dimensional (2-D) wave equation is inverted to recover linear projections of the initial pressure distribution. Second, multiple projections in different directions are used to reconstruct the original 3-D pressure distribution via a 2-D inverse Radon transform.

In this work, we concentrate on the first step of this procedure, i.e., the inversion of the 2-D wave equation. However, also in these modified algorithms, the finite size in the remaining dimensions is usually neglected and artifacts arise in reconstructed images. Our aim is to explain the particular form of these artifacts and to study additional signal processing algorithms that can be applied at the image postprocessing stage to reduce these artifacts. In Ref. 6, two deconvolution methods were proposed for planar measurement geometry: Wiener deconvolution and piecewise polynomial truncated singular value decomposition (PP-TSVD) deconvolution (see Sec. 3). The latter method was actually implemented. For cylindrical geometry, it was mentioned⁷ that the same methods could be used, but it has not been demonstrated on numerical or experimental data. One of the aims in this paper is to test the assertion in Ref. 7 and to compare the two methods in this context. The implementation as well as the results are different from what was reported for the planar case, which justifies our careful re-examination of these methods for cylindrical geometry.

Matrix-model based reconstruction algorithms⁸ can build arbitrary shapes and spatial impulse responses of the detectors directly into the model and are well suited to reduce finite detector artifacts, but require high numerical effort for the large matrices, even the application of graphics processing units for parallel computing.⁹ In our work, we show that these artifacts can be effectively reduced in a matter of a second of additional computation time, at least for some regular geometrical detector shapes, after using standard reconstruction methods like filtered backprojection¹⁰ or time reversal¹¹ that require, in the 2-D case, about a minute on a standard PC.

After a short review of photoacoustic tomography, in Sec. 2 we address some issues that arise for image reconstruction if finite-sized detectors are used. Furthermore, we describe the theory behind deconvolution algorithms for compensation of spin-blur artifacts that appear in connection with detectors that are parts of cylindrical surfaces with finite angular aperture, or parts of their approximating tangential planes if the aperture is not too large. The accuracy of this approximation is also discussed in this section.

In Sec. 3, computer-simulation studies are performed to test the two proposed deconvolution algorithms on simulated data from different source profiles located at various distances from imaging center and recorded with detectors of various apertures. Section 4 describes real measurements on a cylindrical plastisol phantom with six thin holes filled with OrangeG as an absorbing liquid. The holes were located on a spiral emanating from the center of the cylinder. The acoustic signals were recorded by a piezoelectric detector [Polyvinylidenfluorid (PVDF) foil] having the form of a narrow long strip. The time reversal reconstruction shows the expected blurring, and we demonstrate that the blurring can be almost completely eliminated by applying the proposed deblurring algorithms. A preliminary version of this work appeared in Ref. 12.

2.

Background

2.1.

Photoacoustic Tomography with Finite-Sized Detectors

In photoacoustic tomography with line detectors, ultrasound wave propagation of the pressure integrated over the $z$-direction, $p(x,y,t)$, is described by an inhomogeneous 2-D wave equation with a very short initial forcing (see, for example, Ref. 4).

$$\frac{{\partial }^{2}p}{\partial t^2}-c^2\mathrm{\Delta }p=p_0(\mathbf{r})·\frac{\mathrm{d}}{\mathrm{d}t}\delta (t).$$

Here $c$ is the sound speed, $p_0(\mathbf{r})$ is the (integrated) initial pressure produced by thermoelastic expansion, and $\delta (t)$ is the Dirac delta function. This is equivalent to an initial value problem for the homogeneous 2-D wave equation. ⁵

$$\frac{{\partial }^{2}p}{\partial t^2}-c^2\mathrm{\Delta }p=0p(\mathbf{r},0)=p_0(\mathbf{r})\frac{\partial }{\partial t}p(\mathbf{r},0)=0.$$

The goal in photoacoustic tomography is to recover the spatial distribution of absorbed energy density inside the sample, which is proportional to $p_0(\mathbf{r})$, from acoustic pressure signals $p({\mathbf{r}}_S,t)$ measured outside the sample at the surface $S$ using a detector D (see Fig. 1). Ideal line detectors (in real 3-D space) correspond to ideal point detectors for the 2-D wave equation.

Fig. 1

Illustration of photoacoustic tomography. D represents a detector at position ${\mathbf{r}}_S$ scanning across a closed surface $S$ enclosing the initial source $p_0(\mathbf{r})$. $\mathit{R}$ is the volume enclosed by $S$.

A simple and efficient way to solve the 2-D photoacoustic inverse problem (or analogously the 3-D problem) for point detectors and to recover the initial pressure distribution is a time reversal algorithm:¹¹ assume a time $T_0>0$ to be large enough that the pressure satisfies $p(\mathbf{r},T_0)\simeq 0$ within the whole volume $R$ enclosed by the surface $S$. As the initial pressure distribution outside the surface $S$ is zero, for the 3-D problem, $T_0$ can be chosen, e.g., as the maximal diameter of $R$, divided by $c$. A 2-D wave propagation does not satisfy a strict form of the Huygens principle, which implies that the solution of the wave equation never exactly vanishes after a certain time, but exhibits an infinite, algebraically decaying tail. However, using the Abel transform and its inverse, one can see that the long tail in the 2-D measured signal contains no further information about the initial pressure distribution (see Ref. 12). So it is justified also in 2-D to choose such a time $T_0$ when all pressure transients above a small level of size $\epsilon$ (e.g., $1/1000$ of the detected pressure maximum) have passed all detection points. At that time $T_0$, the time evolution of the pressure field is reversed; we start to “rewind the film.” Mathematically, this can be accomplished by setting the pressure values on the surface $S$ equal to $p({\mathbf{r}}_{\mathbf{S}},T_0-t)$, the measured pressure values in time-reversed order. A finite difference or Fourier transform scheme for the wave equation is executed on a spatial computational grid with these time-reversed values as boundary conditions. After another time period $T_0$, we obtain an approximation to the initial pressure $p_0(\mathbf{r})$. We have recently proposed a pseudospectral variant of time reversal that yields highly accurate results in 2-D and 3-D.¹³ This combines an exact time discretization with a Fourier ($\mathbf{k}$-space) approximation to the spatial derivatives $[p_n(\mathbf{k})=p(\mathbf{k},n\mathrm{\Delta }t)]$.

Eq. (1)

$$p_{n+1}(\mathbf{k})=[2–4{\sin }^2(ck\mathrm{\Delta }t/2)]·p_n(\mathbf{k})-p_{n-1}(\mathbf{k}).$$

The time needed for the Fourier transforms (enabling high accuracy in the spatial derivatives) is offset by the fact that the Courant-Friedrichs-Lewy condition does not apply here and, thus, larger timesteps can be taken.

In this work, we will assume the more realistic scenario that instead of point-wise pressure data we have pressure values integrated over a spatially extended aperture. In Fig. 2, we illustrate our strategy of accounting for such a finite detector size in the simple case of a planar detector array: according to source-receiver reciprocity for acoustic waves, measuring signals from a point-like source with a planar detector array is equivalent to measuring signals from an extended source with an array of point-like detectors. Hence, we may (almost) recover the imaging result from a point-like detector array by deconvolution of the reconstructed image that was obtained by measuring with an array of finite-size detectors.

Fig. 2

For a planar integrating detector array, by source receiver reciprocity, we can use the well-known reconstruction formulas from point-like detectors followed by a spatial deconvolution to approximately reconstruct the initial pressure distribution, which is a point source in this figure.

A similar compensation strategy can also be applied for cylindrical or spherical detection geometry: using the translational and rotational invariance of the wave equation, Haltmeier and Zangerl show in Ref. 14 that the reconstructed image becomes blurred with a certain blurring kernel if a finite-sized detector is used instead of a point (3-D) or line (2-D) detector. This has been proven in 2-D for a detector extending to a cylinder of finite radius, which applies to optical detection by a laser beam¹⁵ or an elongated optic fiber sensor.¹⁶ In our case, a piezoelectric detector is used that extends tangentially to a finite strip (Sec. 4.1). The theory of Ref. 14 can also be adapted to such a situation as we show in the Appendix.

A different kind of finite-size issue is raised by the fact that infinite line detectors in reality can only be approximated by finite lines. In our experiment (Sec. 4.1), we have photoacoustic sources that are located in an $xy$-plane, detected at a radial distance of $r_S=16\mathrm{mm}$ from the coordinate origin. The detectors are rectangular strips whose long axis points in the $z$-direction. To illustrate the influence of the $z$-extent of the detectors in this setup, Fig. 3 shows simulated signals from a small ball source with 1-mm diameter, measured by an integrating strip detector growing in $z$-length. For the $3\times 3\text{-}\mathrm{mm}$ detector, one gets an N-shaped 3-D signal similar to a signal observed by a point detector. For the 10-mm-long detector, the signal initially behaves like the 2-D signal, but shows a clear deviation from the 2-D signal after 11 μs. For the $3\times 25\text{-}\mathrm{mm}$ detector, one gets the 2-D pressure signal up to a measurement time of 13 μs; then a small negative peak appears. This peak, which is due to the arrival of the wave at the ends of the strip, could be shifted to longer times by taking an even longer strip; but it is apparent that the signal and, hence, the reconstruction at this $z$-length, are already sufficiently close to the ideal 2-D situation. A more detailed discussion of the required length of such integrating detectors can be found in Ref. 17. This simulation also demonstrates the gain in signal amplitude by a factor of $\sim 1.5$ with increasing integration length of the detector.

Fig. 3

Simulated signals from a small ball with 1-mm diameter as a source measured at a distance of 16 mm with an integrating detector of constant $x$-width (3 mm) and growing in $z$-length: $3\times 3\mathrm{mm}$ (solid), $3\times 10\mathrm{mm}$ (dash), $3\times 25\mathrm{mm}$ (dash-dot); for comparison infinite strip detector signal (dotted); the gain in positive peak amplitude by a factor of $\sim 1.5$ mainly occurs between $z$-lengths of 3 and 6 mm.

2.2.

Algorithms for Compensation of Spin Blur

As described in the previous subsection, a narrow strip detector of sufficient length measures signals governed by the 2-D wave equation. By rotating the detector along a curve, image reconstruction can, therefore, be performed by inversion of the 2-D wave equation. Due to the width of the strip detectors, application of reconstruction methods designed for ideal line detectors will result in a blurred reconstruction.

Suppose that instead of the pointwise pressure values $p({\mathbf{r}}_S,t)$, integrals $p^{\alpha }({\mathbf{r}}_S,t)$ of the acoustic pressure over an arc centered at points ${\mathbf{r}}_S$ on a circular recording curve are measured. Here the superscript indicates that the detector covers an angle $2\alpha$ on the recording surface (compare Fig. 5). Then applying reconstruction methods, such as time reversal or backprojection, the integrated pressure values $p^{\alpha }({\mathbf{r}}_S,t)$ yield as reconstruction (see the Appendix).

$$p_0^{\alpha }(r,\phi )=\underset{-\alpha }{\overset{\alpha }{\int }}{p}_0(r,\phi +{\alpha }^{\prime })\mathrm{d}{\alpha }^{\prime }=\underset{0}{\overset{2\pi }{\int }}{p}_0(r,{\phi }^{\prime })k_{\alpha }(\phi -{\phi }^{\prime })\mathrm{d}{\phi }^{\prime }.$$

Hence, the obtained reconstruction $p_0^{\alpha }(r,\phi )$ [expressed in polar coordinates $\mathbf{r}=(r,\phi )$] is a blurred version of the unblurred initial pressure distribution $p_0(r,\phi )$, obtained by convolving $p_0(r,\phi )$ with a purely angular blurring kernel given by the box function

$$k_{\alpha }(\phi -{\phi }^{\prime })=\{\begin{array}{ll}1& \text{if}|\phi -{\phi }^{\prime }|\le \alpha \\ 0& \text{if}|\phi -{\phi }^{\prime }|>\alpha \end{array}.$$

We emphasize that from a mathematical point of view, there is no restriction of our proposed deconvolution methods to a box function as blurring kernel. The algorithms can be applied with minor modifications to a more general kernel, i.e., a general signal impulse response function describing a variable sensitivity within the detector foil.

Since the spatial extent of a circular arc with aperture $2\alpha$ is proportional to its radius, the blurring effect increases linearly for increasing distance to the center of the reconstruction domain and is equal to the detector width close to detection circle. We refer to the purely angular blurring proportional to the distance of the center of rotation as spin blurring because the same type of blurring arises in photographic imaging when a camera is spinning around the lens axis during exposure. Figure 4 shows an initial pressure distribution consisting of four points (smoothed by box function convolution) located on a spiral emanating from the detection center (left), and the time reversal reconstruction using signals from an arc-shaped detector with $2\alpha =45\mathrm{deg}$ in Cartesian (center) and polar (right) coordinates.

Fig. 4

Initial pressure (a) and two-dimensional (2-D) time reversal reconstruction in Cartesian (b) and polar (c) coordinates for a circular integrating detector of aperture $2\alpha =45\mathrm{deg}$ (bold white arc), which is rotated in 80 equidistant positions around the sample (white dots). The best resolution is near the rotation center and blurring increases toward the detector circle.

While in Cartesian coordinates, pixels are coupled in a very intricate way by the spin, in polar coordinates, the blurring is constant in $\phi$ and independent of $r$ (Fig. 4, right) and the problem simplifies considerably, with pixels at a given radius decoupled from all other pixels. We convert the $n\times n$ image, which is a discrete representation of $p_0^{\alpha }(\mathrm{r})$, given in Cartesian pixels, to a polar coordinate image ${\mathbf{Y}}_{\mathrm{pol}}$, using spline interpolation to obtain pixel values at $[n/2]$ values of r and $2n$ values of $\phi$. Each column of ${\mathbf{Y}}_{\mathrm{pol}}$ corresponds to a particular radius, and therefore, the blur of pixels in that column is independent of pixels in every other column. Thus, we have $[n/2]$ independent blurring problems, each of size $2n$ and all involving the same blurring kernel.

Let us, for the moment, omit the radial variable and write $y(\phi )$ for the blurred initial pressure $p_0^{\alpha }(r,\phi )$ and $x_0(\phi )$ for the desired unblurred version $p_0(r,\phi )$. We then face the problem of estimating the signal $\phi \mapsto x_0(\phi )$ from the data

$$y(\phi )=(x_0*k_{\alpha })(\phi )+n(\phi )\text{for}\phi \in [0,2\pi ].$$

Here $n(·)$ is the deterministic or stochastic noise and $k_{\alpha }(·)$ is the blurring kernel. The blurred initial pressure is itself the result of a reconstruction algorithm, such as time reversal, and therefore, the noise will be nonstationary and correlated even if the noise in the original data is white and Gaussian.

Deconvolution is a typical inverse problem, where the exact solution may either not exist or the noise is severely amplified. In order to make the reconstruction process well-posed, regularization techniques have to be applied, where approximate solutions are constructed that are stable with respect to data perturbations. The most common technique for solving a linear inverse problem is Tikhonov regularization. There an estimate for the unknown is defined as the unique minimizer of the Tikhonov functional

$$x\mapsto {\Vert x*k_{\alpha }-y\Vert }^2+\lambda {\Vert x\Vert }^2.$$

The minimizer of the strictly convex Tikhonov functional can be found by setting its gradient to zero. Recalling that the convolution operation is diagonalized by the Fourier transform, the unique minimizer of the Tikhonov functional is given by $x=y*g_{\alpha }$, where the filter $g_{\alpha }$ has the Fourier representation

$$F[g_{\alpha }](\omega )=\frac{\overline{F[k_{\alpha }](\omega )}}{{|F[k_{\alpha }](\omega )|}^2+\lambda }.$$

Here $F[h](\omega )=\stackrel{2\pi }{\underset{0}{\int }}h(\phi )\exp (-i\omega \phi )\mathrm{d}\phi$ is the Fourier transform of some function $h$ and the bar indicates complex conjugation. For signal and image deconvolution with additive random noise, this estimator is also known as Wiener deconvolution. The angular extent of the blurring kernel is given by $2\alpha =2\arctan (w/2r_S)$, where $w$ is the detector width and $r_S$ is the radius of the detection circle (Fig. 5). In the context of angular spin deblurring, Wiener deconvolution works best when combined with reconstruction algorithms that already use polar coordinates, such as the Fourier domain algorithms derived in Refs. 18 and 19. However, it may also be applied separately on images reconstructed by any other algorithm.

Fig. 5

Schematic of curved detector versus flat detector.

For our numerical implementation, we use a slightly different formulation:^20,21 the $[n/2]$ blurring problems are written in discretized form as $\mathbf{K}{\mathrm{x}}_{\mathrm{pol}}\cong {\mathrm{y}}_{\mathrm{pol}}$, where ${\mathrm{y}}_{\mathrm{pol}}$ is one column of ${\mathbf{Y}}_{\mathrm{pol}}$. The blur matrix $\mathbf{K}$ is a $2n\times 2n$ circulant band matrix with band size proportional to the blurring angle and has an SVD $\mathbf{K}=\mathbf{U}\mathrm{\Sigma }{\mathbf{V}}^{\mathrm{T}}$, where $\mathbf{U}$ and $\mathbf{V}$ are orthogonal matrices and $\mathrm{\Sigma }$ is diagonal. We again define regularized solutions as minimizers on ${\mathbb{R}}^{2n}$ of the Tikhonov function $x\mapsto {\Vert \mathbf{K}\mathrm{x}-{\mathrm{y}}_{\mathrm{pol}}\Vert }^2+\lambda {\Vert \mathrm{x}\Vert }^2$. Using the SVD of the matrix $\mathbf{K}$, the minimizer of the Tikhonov function is given by

$${\mathrm{x}}_{\mathrm{pol}}={({\mathbf{K}}^{\mathrm{T}}\mathbf{K}+\lambda \mathbf{I})}^{-1}{\mathbf{K}}^{\mathrm{T}}{\mathrm{y}}_{\mathrm{pol}}=\mathbf{V}{({\mathrm{\Sigma }}^2+\lambda \mathbf{I})}^{-1}{\mathrm{\Sigma }\mathbf{U}}^{\mathbf{T}}{\mathrm{y}}_{\mathrm{pol}}.$$

Since any circulant $N\times N$ matrix has complex eigenvectors

$${\mathrm{v}}_m=[1,\exp (\frac{2\pi \mathrm{im}}{N}),\cdots ,\exp (\frac{2\pi \mathrm{im}(N-1)}{N}\left)\right]\text{for}m=0,\dots ,N-1,$$

the SVD in this case essentially amounts to a discrete Fourier transform. A small difference to Wiener deconvolution as described above is that here discretization to a finite vector is done before Tikhonov regularization. In the context of general least squares problems, the first approach is often referred to as “solving the normal equations by eigen-decomposition” as opposed to “solving via SVD.”²² The matrix formulation is also convenient to apply the standard generalized cross-validation (GCV) algorithm²³ to determine an optimal regularization parameter λ for each image and to diminish the influence of noise. The columns ${\mathrm{x}}_{\mathrm{pol}}$ give us a matrix ${\mathbf{X}}_{\mathrm{pol}}$, which we transform back to Cartesian pixels $\mathbf{X}$, again using spline interpolation. Pixels in $\mathbf{X}$ outside the circular domain are set to zero. The image $\mathbf{X}$ is then a discrete approximation of the unblurred initial pressure $p_0(x,y)$ in Cartesian coordinates. We note that with minor modifications, the squared ${\ell }^2$-norm ${\Vert x\Vert }^2$ in the Tikhonov functional can be replaced by the squared ${\ell }^2$-norm of a solution derivative. However, replacing ${\Vert x\Vert }^2$ by a nonquadratic penalty makes the minimization of the resulting Tikhonov functional more difficult.

In Ref. 24, Hansen described such a variant of TSVD that seeks to regularize with the ${\ell }^1$-norm of a solution derivative (in our case, the first derivative) instead of regularizing with the ${\ell }^2$-norm of the solution. Recall that the ${\ell }^1$-norm of a vector is defined by the sum of the absolute values of all components. The ${\ell }^1$-norm regularization produces a vector of delta functions (sparse spike train) for the first derivative, and therefore, the solutions will be piecewise constant functions with at most $k$ discontinuities, where $k$ is the truncation index in the TSVD. If higher solution derivatives are used for regularization, piecewise polynomial solutions are obtained, hence the name PP-TSVD. Obviously, this deblurring method will work best for images that are themselves better described by piecewise constant functions (having sharp edges) than by smooth functions. The method was reported in the context of photoacoustic imaging with planar detection geometry in Ref. 6, where it was demonstrated on a phantom of two cylinders of constant optical absorbance. In Ref. 7, it is mentioned that it could also be transferred to a cylindrical detection geometry.

In Ref. 25, Hansen discusses PP-TSVD for deconvolution of 2-D images in the case of a separable 2-D convolution kernel. In our case, the kernel (when expressed in polar coordinates) is not only separable, but also independent of $r$, so that the radial part after discretization involves only the identity matrix. The extension to 2-D as described by Hansen (using the Kronecker product of the two kernel factor matrices), in this case, reads as follows: apply the deconvolution method using the same blurring matrix $\mathbf{K}$ as for Wiener deblurring; for each column representing a fixed radius, execute a one-dimensional PP-TSVD, where the algorithm is taken from Ref. 26. The optimal choice of the regularization parameter $k$ was not considered in Ref. 25. In the section on reconstructions from computer-simulated data, we will determine an optimal $k$ empirically by comparing reconstruction quality for a sample image and a series of $k$-values.

2.3.

Approximation of a Curved Detector by a Flat Detector

For the measurements described in Sec. 4.1, a $3\times 25\text{-}\mathrm{mm}$ PVDF-foil detector was used, which is a good approximation of an extended integrating line detector (Fig. 3). However, this detector is not curved but flat. In this section, we investigate the difference between curved and flat aperture analytically.

Omitting the $z$-coordinate, the pressure data from the experiment are integrated over a line of length $w$ tangent to a detection circle of radius $r_S$ (Fig. 5). Denoting by $l=l({\alpha }^{\prime })=r_S\tan ({\alpha }^{\prime })$ the distance of a point on the line segment from its center, using a first-order Taylor approximation of $p(r,\phi ,t)$ in the radial variable (around the center $r=r_S$), and applying the approximation $\sqrt{1+\tan {({\alpha }^{\prime })}^2}=1+[({\alpha }^{\prime 2})/2]+O({\alpha }^{\prime 4})$ for small ${\alpha }^{\prime }$, we obtain

$$\begin{gathered} \underset{-w/2}{\overset{w/2}{\int }}p(r,\phi ,t)\mathrm{d}l=\underset{-\alpha }{\overset{\alpha }{\int }}p[r_S\sqrt{1+\tan {({\alpha }^{\prime })}^2},{\phi }_0+{\alpha }^{\prime },t][1+\tan {({\alpha }^{\prime })}^2]r_S\mathrm{d}{\alpha }^{\prime } \\ =\underset{-\alpha }{\overset{\alpha }{\int }}p(r_S,{\phi }_0+{\alpha }^{\prime },t)r_S\mathrm{d}{\alpha }^{\prime }+\underset{-\alpha }{\overset{\alpha }{\int }}p(r_S,{\phi }_0+{\alpha }^{\prime },t){\alpha }^{\prime 2}{r}_S\mathrm{d}{\alpha }^{\prime } \\ +\underset{-\alpha }{\overset{\alpha }{\int }}\frac{\partial p}{\partial r}(r_S,{\phi }_0+{\alpha }^{\prime },t)\frac{{\alpha }^{\prime 2}}{2}{r}_S^2\mathrm{d}{\alpha }^{\prime }+O({\alpha }^4). \end{gathered}$$

This implies that the difference

$$\underset{-w/2}{\overset{w/2}{\int }}p(r,\phi ,t)\mathrm{d}l-\underset{-\alpha }{\overset{\alpha }{\int }}p(r_S,{\phi }_0+{\alpha }^{\prime },t)r_S\mathrm{d}{\alpha }^{\prime },$$

between the pressure integrals over curved detector and flat detector is of order $O({\alpha }^3)$. As a consequence, the measured signals for curved and flat detector shapes agree up to second order in $\alpha$. Hence, for small apertures $2\alpha$ (in our experiment, $2\alpha =10.7\mathrm{deg}$), the measured signals will differ little for the two detector types and, therefore, also the blurred reconstructions (which follows, for example, from the stability of the time reversal algorithm).

However, if the aperture gets larger ($>20\mathrm{deg}$), a blur that was of constant width in the radial coordinate will split up at both ends for the flat detector. In fact, the ends look like a pair of tongs grasping two balls with radii proportional to the largest deviation of the flat detector from the detection circle [Fig. 6(b)].

Fig. 6

Reconstruction of a ball of 1-mm diameter located 12 mm from detection center ($r_S=16\mathrm{mm}$) from a flat strip detector with aperture 20 deg (a) and 40 deg (b). The bold white bar indicates a flat detector that is tangent to the detection circle (dashed white). In both cases, we applied a reconstruction algorithm that is theoretically exact for point-like data.

3.

Computer-Simulation Studies

For testing the two deconvolution methods described in the previous section, we use a phantom similar to the four-point phantom shown in Fig. 4. The points are refined in their structure and described by spatial profiles that are either hollow cylinders (HC) or Gaussians (G). The HC profile type was slightly smoothed by convolving with a box function to make it look more realistic and to make it more interesting for PP-TSVD. We want to examine reconstruction quality in dependence of source distance from center and of detector aperture. Quality will be assessed by visual inspection of 2-D profiles; after visual inspection, it will become apparent which quantitative measures are additionally useful.

In all numerical studies, the solution of the 2-D wave equation has been computed by numerically solving the wave equation. The image reconstruction (prior to deconvolution) has been performed by time reversal. The computational domain for the simulations is $2\times 2\mathrm{mm}$, the detection circle radius $r_S=0.8\mathrm{mm}$, and the spatial grid size $\mathrm{\Delta }x$ is $2/500\mathrm{mm}$ for the forward computation and $2/441\mathrm{mm}$ for time-reversed computation. Also, the timestep size $\mathrm{\Delta }t$ as well as the ending and starting times $T_0$ in the forward and time-reversed computations are chosen somewhat different to avoid an inverse crime. Both the forward computations (using $p(\mathbf{k},n\mathrm{\Delta }t)=p(\mathbf{k},0)·\cos (ckn\mathrm{\Delta }t)$, $\mathrm{\Delta }t=1.33\mathrm{ns}$, sound speed $c=1500\mathrm{m}/\mathrm{s}$) and the time-reversed computations [using Eq. (1) and $\mathrm{\Delta }t=1.68\mathrm{ns}$] are performed in $\mathbf{k}$-space. Pressure signals are recorded at 720 equidistant points around the detection circle. A curved integrating detector of aperture $2\alpha$ degrees is simulated by summing the recorded signals of consecutive measurement points lying within the designated aperture and dividing by the number of these points (in order to make signal sizes comparable to the point detector case, although physically a summation and not an averaging is done). The detector aperture $2\alpha$ (in degrees) is chosen in a way that $2\alpha /360·720=4\alpha$ is an integer, which is the case for selecting $2\alpha =10\mathrm{deg}$ or $2\alpha =20\mathrm{deg}$. After image reconstruction, we present profiles through the center of a point object of the phantom. These profiles are drawn over abscissas representing a pixel sequence in either the radial or tangential direction. Note that the profiles in the radial and tangential direction are equal for the initial pressure, but will be different for the reconstructions. For demonstrating the effect of source distance from center, we display the profiles for the first and the fourth point on the spiral (counted from center). The effect of detector aperture is shown by displaying the results for $2\alpha =10\mathrm{deg}$ and $2\alpha =20\mathrm{deg}$. In summary, our tests with simulated data include 32 simulations (two choices for methods, points, profiles, apertures, and directions).

At this point, we want to determine an optimal parameter $k$ for PP-TSVD deconvolution. The effect of $k$ on deconvolution is examined using $k=80$, $k=200$ $k=n=442$ (the reconstructed image is $n\times n\text{pixles}$), and $k=850$. In Fig. 7, a close-up of the third and fourth point reconstructions for the HC profile is provided. For the annular $xy$-profiles of the HC sources, the reconstruction quality is distinctly better along a cross-section in radial direction than along a cross-section in tangential direction. This feature can be observed in all deblurred reconstructions, including those using Wiener deconvolution. However, with PP-TSVD the shape of the reconstruction is more sharp-edged for $k=80$, and is smoother for the higher-$k$ reconstructions. On the other hand, computing time increases about linearly with $k$ (ranging from 30 min for $k=80$ to $>3\mathrm{h}$ for $k=850$); moreover, the background noise level increases somewhat with $k$, becoming intolerable when $k$ approaches $2n$. Inspecting Fig. 7, we decided to use $k=n=442$ in our simulations. In any case, choosing $k=n$, i.e., truncating at half matrix size, seems to be reasonable for the other images as well in our scenario if no additional noise is added.

Fig. 7

Comparison of reconstructions for third and fourth point [hollow cylinders (HC) profile] for varying piecewise polynomial truncated singular value decomposition (PP-TSVD) regularization parameter: (a) $k=80$, (b) $k=200$, (c) $k=442$, and (d) $k=850$.

For the visual inspection, we display the results on four 4-subplot figures. Each four-subplot figure shows four curves corresponding to initial pressure (solid), reconstruction without deconvolution (dash-dotted), reconstruction after Wiener deconvolution (dashed), and reconstruction deconvolved with PP-TSVD (dotted). The four subplots in each figure show the profiles through the two chosen points in the two directions. The amplitude range for display is always [$-0.1$ 1.1], assuming that the initial pressure is normalized to [0,1]. The following general observations pertain to reconstructed amplitudes: with point detectors, our algorithm returns amplitude values in the range [$-0.12$, 0.90] for our phantoms. The values returned with the integrating (here averaging) detectors lie in the range [$-0.08$, 0.78] because they are smeared over a larger area. Deconvolution provides partial amplitude compensation: values lie in the range [$-0.13$, 0.83]. Note that these numbers apply only to the point closest to the center that is reconstructed best; for the other three points, the amplitude degradation is even higher.

The inspection of Figs. 8 and 9 yields two surprising insights: both the deconvolution method and the detector aperture hardly have any influence on the deblurred results. There is some amplitude decrease between the first and fourth point as mentioned before, and the difference of reconstruction quality in radial and tangential directions is very pronounced and affects both amplitude and full width at half maximum (FWHM) of the profiles. A distinctive feature is also the much stronger smoothing of the hollow interior of the HC profile for the fourth point, compared to the first point. In Figs. 10 and 11, the same simulations were performed for Gaussian (G) initial profiles. The two insights gained with the HC profiles are confirmed once more: there is virtually no difference among (deblurred) results for the two deconvolution methods and two detector apertures. The amplitude decrease from first point to fourth point is as expected and the difference between radial and tangential profiles is again pronounced, but here it mainly affects the width of the Gaussian profile.

Fig. 8

Detector aperture $2\alpha =10\mathrm{deg}$: four subplots, each showing profiles of initial pressure HC (blue, solid), reconstruction before deblurring (green, dash-dot), reconstruction after Wiener deblurring (red, dashed), and reconstruction after PP-TSVD deblurring (black, dotted), all for first point $r$-direction (a), first point $\phi$-direction (b), fourth point $r$-direction (c), and fourth point $\phi$-direction (d).

Fig. 9

Detector aperture $2\alpha =20\mathrm{deg}$: four subplots, each showing profiles of initial pressure HC (blue, solid), reconstruction before deblurring (green, dash-dot), reconstruction after Wiener deblurring (red, dashed), and reconstruction after PP-TSVD deblurring (black, dotted), all for first point $r$-direction (a), first point $\phi$-direction (b), fourth point $r$-direction (c), and fourth point $\phi$-direction (d).

Fig. 10

Detector aperture $2\alpha =10\mathrm{deg}$: four subplots, each showing profiles of Gaussian initial pressure (blue, solid), reconstruction before deblurring (green, dash-dot), reconstruction after Wiener deblurring (red, dashed), and reconstruction after PP-TSVD deblurring (black, dotted), all for first point $r$-direction (a), first point $\phi$-direction (b), fourth point $r$-direction (c), and fourth point $\phi$-direction (d).

Fig. 11

Detector aperture $2\alpha =20\mathrm{deg}$: four subplots, each showing profiles of Gaussian initial pressure (blue, solid), reconstruction before deblurring (green, dash-dot), reconstruction after Wiener deblurring (red, dashed), and reconstruction after PP-TSVD deblurring (black, dotted), all for first point $\mathbf{r}$-direction (a), first point $\phi$-direction (b), fourth point $\mathbf{r}$-direction (c), and fourth point $\phi$-direction (d).

As a quantitative measure we list in Table 1 the FWHM of all the Gaussian-like profiles in Figs. 10 and 11. The measures are given in pixel units $\mathrm{\Delta }x$ (spatial step size of the time-reversed grid). Linear interpolation was done between the integer pixel grid, so that accuracy is $\pm 0.05\mathrm{\Delta }x$. As can already be seen from visual inspection, the deconvolved reconstruction significantly improves resolution compared to the blurred image in terms of reducing the FWHM in $\phi$-direction of isolated objects. However, there are other aspects of degraded image quality, such as the filling up of holes, which can be seen in Figs. 8(c), 8(d), 9(c), and 9(d), that cannot be avoided even with deconvolution and that are strongly dependent on source distance from center.

Table 1

Full width at half maximum (FWHM) in pixel units of the profiles displayed in Figs. 10 and 11.

An explanation of the similarity between PP-TSVD and Wiener deconvolution is that the high number (442) of possible discontinuities produces quasi-smooth solutions in the PP-TSVD case also (note that our plots show only 27 or 33 pixel values) so that the different regularizers produce very similar results. The constancy in results for varying aperture could be confirmed for much higher apertures also. The noticeable amplitude degradation for an increased distance of the source from the center may result from the fact that a uniform grid spacing in polar coordinates transforms to a grid spacing in Cartesian coordinates that decreases in density with distance from center. Since our algorithms work separately on each radius, using a denser angular grid for large radii might reduce this effect. Such an investigation is beyond the scope of our present work.

Computation time for Wiener deblurring is 1 s; for PP-TSVD deblurring with an appropriate $k$, it is at least an hour if the code pptsvd.m from Ref. 26 is used. The documentation for this algorithm (see Ref. 25, Sec. 5) reveals that it was written especially for large-sized problems. For our small-sized problems, the following measures would make the algorithm faster: (1) instead of MATLAB® function svds, use svd to perform SVD; (2) optimize the function l1c to solve a discrete ${\ell }^1$ linear optimization. We estimate that with these measures, the algorithm for $442\times 442\text{pixel}$ reconstruction in this section and for $251\times 251\text{pixel}$ reconstruction in Sec. 4.2 could be made as fast as a few seconds. Still, the Wiener deconvolution is faster and provides the same image quality.

So we can currently recommend the Wiener deblurring for problems of the type presented in this work. It must be stressed, however, that image-specific regularization parameter optimization via GCV or a similar algorithm is essential for this excellent quality of Wiener deblurring. With these provisions, we were not able to observe any instabilities or oscillations as was reported for planar geometry in Ref. 6. Further numerical studies shown below demonstrate that these conclusions are also valid in the case of noisy data.

It is well known that the performance of deblurring algorithms often strongly depends on the noise statistics in the data. For our algorithms, we investigated this issue numerically by adding Gaussian white noise to the numerically computed solution of the 2-D wave equation and applying time reversal to the noisy data. The results of Wiener deblurring and the PP-TSVD deblurring algorithm on the intermediate images are shown in Fig. 12. For Wiener deconvolution, the regularization parameter has again been computed in a data-driven way by applying the GCV algorithm. Note that the regularization parameter computed by GCV increases with increasing noise level. The result is $\lambda =9.2\times {10}^{-4}$ for the simulated data without noise, $\lambda =2.5\times {10}^{-3}$ for 3% noise, and $\lambda =2.7\times {10}^{-3}$ for 10% noise. Such a behavior is expected for an ill-conditioned problem since the regularization parameter acts as trade-off between accuracy for exact data and stability with respect to noise. For the PP-TSVD algorithm, the parameter $1/k$ acts as regularization parameter, and therefore, $k$ should be chosen smaller for higher noise levels. Optimal values for $k$ depending on the noise level have been found empirically by performing similar experiments as in the case of noise-free data. For the result shown in Fig. 12, we found $k=300$ to be optimal for 3% noise and $k=120$ for 10% noise.

Fig. 12

Reconstruction results using detector aperture $2\alpha =10\mathrm{deg}$ after adding 3 and 10% Gaussian white noise, respectively, to the solution of the 2-D wave equation. (Here by p% noise we mean that the standard deviation of the added random variable is p% of the maximal data value.) The regularization parameters for the Wiener deconvolution (red, dashed) have been determined by the generalized cross-validation. For PP-TSVD deblurring (black, dotted), we have taken $k=300$ in the case of 3% noise and $k=120$ in the case of 10% noise, compared to $k=442$ (half dimension of the matrix) without noise in Figs. 10(c) and 10(d).

It can be seen that both deblurring algorithms still provide a comparable image quality. Since time reversal is a linear algorithm, the noise in the blurred intermediate images is still Gaussian but neither uncorrelated nor stationary. However, because of the rotational invariance of the image reconstruction problem, the noise is stationary in the angular direction (i.e., long the angular direction, the variance of the noise is constant and the correlation depends only on angular spacing). This may be exploited to further improve the performance of the Wiener deconvolution algorithm. A precise analysis of the noise statistics, however, is beyond the scope of this paper.

4.

Experimental Results

4.1.

Experimental Setup

Top and front views of the measurement setup are shown in Fig. 13. A homemade piezoelectric ultrasound detector was used to record the acoustic signals. The detector was made of a strip of PVDF-foil with a thickness of 28 μm, with electrodes deposited on either side attached on a polymer substrate (backing material). The width $w$ and the length $L$ of the stripe (sensitive area) were 3 and 25 mm, respectively. The signals were recorded by connecting the detector with an active probe to the digital storage oscilloscope. To improve the SNR, the signals were recorded with bandwidth limitation (dc up to 20 MHz) by the oscilloscope. This is the only limitation to the bandwidth of the ultrasonic transducer.

Fig. 13

Top and front view of the experimental setup.

The cylindrical phantom sample ($D=20\mathrm{mm}$) was made of plastisol with six channels oriented along $z$-direction, distributed along a spiral curve in the cross-section area (see top view in Fig. 13). The sound speed of plastisol is $\sim 1400\mathrm{m}/\mathrm{s}$; further acoustic properties can be found in Ref. 27. The channels, each with a diameter of 1 mm, were filled with a light-absorbing dye solution (OrangeG, absorption coefficient $\sim 100{\mathrm{cm}}^{-1}$).

For photoacoustic excitation, the sample was illuminated by 10-ns laser pulses from a frequency-doubled, Q-switched Nd:YAG laser with a repetition rate of 10 Hz. The output beam had a wavelength of 532 nm and was split into two beams to illuminate the sample from two opposite sides parallel to the $y$-direction. Furthermore, the beams were focused with cylindrical lenses to narrow the extension of the illumination along $z$-direction ($\sim 1.5\mathrm{mm}$), forming a line-shaped illumination pattern on the surface of the sample (see front view in Fig. 13). The rotation axis of the sample was positioned at a normal distance $r_S=16\mathrm{mm}$ to the detector and was oriented parallel to the $z$-direction. The $z$-position of the line-shaped illumination was fixed to the height of the center of the $z$-extension of the detector. Signals were acquired while rotating the sample over five full rotations with an angular increment of 0.9 deg. The five times averaged temporal signals are shown for the 400 orientations over 360 deg in Fig. 14. The dye solution in the holes is strongly absorbing, so the absorbed energy density decreases rapidly inside the holes. This can already be seen from the sinogram in Fig. 14 (signals are generated only near the boundaries of the absorbers), but more clearly in the reconstructed images in Sec. 4.2.

Fig. 14

Measured pressure signals (sinogram) as a function of the rotation angle ($x$ axis) and the time; on the $y$ axis, the time has been multiplied by the sound speed in water ($c=1500\mathrm{m}/\mathrm{s}$). Already in the sinogram it can be seen that the laser light is strongly absorbed by the dye solution in the holes.

4.2.

Reconstruction and Deblurring Results

The sinogram in Fig. 14 shows data measured over 400 equidistant angular positions using a sampling interval of $\mathrm{\Delta }t=10\mathrm{ns}$. Note that the $y$ axis displays $r=ct$ (time multiplied by sound speed in water). Beginning at $r_0=7\mathrm{mm}$, 1210 samples were recorded, ending at $\sim r_1=25\mathrm{mm}$. The computational domain was $33.6\times 33.6\mathrm{mm}$, and the spatial grid size was $33.6/400\mathrm{mm}$. A 2-D time-reversal reconstruction was performed using Eq. (1). The acoustic properties of the plastisol²⁷ were considered to be similar enough to those of water so that the two media were not modeled separately.

For closer inspection of the reconstructed sources, which are located in the inner half of the detection area ($r_S=16\mathrm{mm}$), only the inner $251\times 251\text{pixel}$ square was cut out of the original $401\times 401\text{pixel}$ reconstructed image; then the deblurring algorithms were applied to this inner square. Computation time is 0.5 s for Wiener deconvolution and $\sim 10\min$ for PP-TSVD deconvolution with pptsvd.m. The choice for the regularization parameter $k=n$ is also a reasonable one here.

In Figs. 15 and 16, we see cross-sections through the centers of the six holes in the plastisol phantom in the radial and tangential directions, respectively. Again, the similarity between the results from the two deblurring methods is remarkable, although there are more differences from the simulated data (see fifth hole in Fig. 15, for instance). Amplitude compensation by deblurring is somewhat better in the radial direction than in the tangential direction, whereas FWHM compensation and, thus, resolution is similar for the two directions. In the reconstructions from our experimental data, we observe some unusually large negative pixel values. Reasons that this occurs may be as follows: modeling the plastisol medium the same as water, modeling the temporal impulse response of the piezoelectric detector as a delta function, disturbances or feedback during electronic acquisition. However, looking at the 2-D reconstructed images (Fig. 17, where negative pixel values were set to zero), we conclude that these imperfections do not seriously affect the overall quality of the information about the optical absorbance of the dye solution in the holes. The reconstructions show that most of the absorption occurs near the boundaries of the holes, which is in agreement with the absorption properties of Orange G given in Sec. 4.1 and with the images of the original phantom drawn in Fig. 13.

Fig. 15

Slices through the six reconstructed holes in $(r,z)$-direction with hole distance from center increasing from left to right: before deblurring (blue, solid), after Wiener deblurring (red, dashed), and after PP-TSVD deblurring (black, dotted).

Fig. 16

Slices through the six reconstructed holes in $(\phi ,z)$-direction with hole distance from center increasing from left to right: before deblurring (blue, solid), after Wiener deblurring (red, dashed), and after PP-TSVD deblurring (black, dotted).

Fig. 17

Reconstructed $xy$-profiles of the entire six-hole phantom: before deblurring (left), after Wiener deblurring (center), and after PP-TSVD deblurring (right).

5.

Conclusion

The mathematical model of photoacoustic wave propagation is the 3-D linear wave equation. If line detectors are used that extend to infinity in one direction (or are at least sufficiently long), the reconstruction process typically involves the inversion of a 2-D wave equation for the pressure projections onto the remaining dimensions, followed by a 2-D inverse Radon transform. Considering only the first problem, we are confronted with a cylindrical geometry. Then the ideal 2-D model still is an abstraction from reality if (piezoelectric) detectors have sizeable extent in angular coordinates.

If the detector size is increased from an (idealized) line detector to a finite width, the resolution starts to deteriorate in parts of the image. Xu and Wang have first quantified this effect of the detector size on the spatial resolution in Ref. 28. Haltmeier and Zangerl have explained these effects as a consequence of the rotational and translational invariance of the wave equation.¹⁴ They have described the effect of the finite detector size in terms of a blurring kernel that has to be convolved with the image obtained from point or line detectors. This suggests the use of deconvolution methods to reverse the effect.

We have demonstrated with computer simulations and with an example on experimental data that both the Wiener and the PP-TSVD deconvolution methods yield an excellent compensation of the blurring. The computational cost of Wiener deblurring is very small ($<1\mathrm{s}$); the cost for PP-TSVD deblurring is significantly larger, but still in the range of a few seconds (assuming image sizes of a few hundred pixels squared). So there is every reason to apply Wiener deblurring, together with appropriate optimization of the regularization parameter, for any blurring kernel that is explicitly known. We did not observe any of the problems with Wiener deblurring that were reported in Ref. 6. If a 3-D object is to be imaged, this deconvolution can be applied at the first step that involves the inversion of the 2-D wave equation. After correction of the 2-D projections, the final step using the inverse Radon transform can be applied as usual for obtaining a full 3-D image.

Matrix-model based reconstruction algorithms^8,9 (or Refs. 29 and 30, using a different model) have the capability to incorporate arbitrary detector shapes and spatial impulse responses into the imaging model and, thus, into the reconstruction algorithm itself. However, such approaches require solving optimization algorithms working iteratively with very large-sized matrices and entail high computational effort. The conclusion of Ref. 9 states that “current implementations of (3-D) iterative image reconstruction algorithms still require several days to process a densely sampled data set.” In Ref. 29, we read “the run time for constructing and inverting the (2-D) model matrix for a $81\times 81\text{grid}$ was approximately one hour.” While this cannot be avoided for more complicated ultrasonic transducer shapes, we demonstrated that effective reduction of imaging artifacts in 2-D caused by the finite size of some regular-shaped detectors is possible in the framework of traditional reconstruction algorithms (like time-reversal or filtered backprojection) and takes only a fraction of the time needed with matrix-model algorithms.