2.1.

Cylindrical Detector Model

Figure 2 shows the cylindrical detector used for this study. The active sensor area has a concave cylindrical shape, determined by its radius of curvature $R$, its width $x_d$ in $x$-direction, and its height $z_d$ in $z$-direction. Its NA is given by $\mathrm{NA}=z_d/(2R)$. For acquiring a section image, the object is rotated about an axis parallel to the $z$ axis. Several sections can be imaged by translating the object or the sensor in $z$-direction.

Fig. 2

Detector surface for section imaging. The detector is characterized by its radius of curvature $R$, its width in $x$-direction, $x_d$, and its height in $z$-direction, $z_d$. The coordinates relative to the detector are $s$ and $z$ and a point on the detector cross-section is denoted as $(s^{\prime },z^{\prime })$. To acquire an image, the sample is rotated about an axis parallel to the $z$-direction and is shifted by $z_m=m\mathrm{\Delta }z$ to cover the entire region of interest. By shifting the object, a target point at $z_n$ moves across the focus plane at $z=0$.

The width $x_d$ exceeds the size of the object to be imaged. To analyze the properties of such an integrating detector, we define a coordinate system relative to the detector and describe its surface as an arrangement of parallel lines, perpendicular to the $z$ axis. For a single line oriented along $x$ direction, the received signal $P_{\varphi }(y,z,t)$ is the pressure field $p_{\varphi }(x,y,z,t)$ coming from the object, integrated over the coordinate $x$. Subscript $\varphi$ means that the object has been rotated by an angle $\varphi$ relative to a reference orientation.

Pressure wave $p_{\varphi }(x,y,z,t)$ is the solution of the acoustic wave equation with $p_{\varphi }(x,y,z,0)=\mathrm{\Gamma }{w}_{\varphi }(x,y,z)$ as initial condition, where $\mathrm{\Gamma }$ is the Grüneisen parameter and $w_{\varphi }$ is the absorbed energy density in the rotated object. It can be shown that the integrated pressure signal is a solution of the 2-D wave equation with source $\mathrm{\Gamma }{W}_{\varphi }(y,z)$,¹⁴ where $W_{\varphi }(y,z)$ is a 2-D Radon transform of $w_{\varphi }(x,y,z)$ and is given by

Validity of the 2-D solution requires integration limits from minus to plus infinity in Eqs. (1) and (2). It can also be shown, however, that with a finite detector exact integrated signals are measured, provided that the width $x_d$ is sufficiently large.⁴ For a perfectly focusing detector, the source can be described by $W_{\varphi }(s)\delta (z)$ and the detector by a single line in the plane at $z=0$, located at position $y_0$. Here, we use coordinate $s=y–y_0$, which is the distance between the line detector and a point in the source. It is proportional to the arrival time $t=s/c_s$ of the wave emanating from point $y$, where $c_s$ is the speed of sound. The distribution of energy density in a section can then be qualitatively reconstructed by applying the inverse Radon transform to the measured pressure signals.¹² However, it has to be taken into account that the pressure signals from a finite source are bipolar, whereas the energy density is purely positive. It is possible to obtain from $p$ a quantity that is proportional to the Radon transform of $W$ by applying an Abel transform.¹⁵

For a single line detector, $q_{\varphi }(t)$ is proportional to the energy density projected into the plane perpendicular to the line and integrated over a circle with radius $c_{s}t$ around the line position.¹⁵ For the focusing detector, which receives signals only from the slice at $z=0$, the inverse Radon transform of $q_{\varphi }(t)$, measured at a sufficient number of angular positions around the imaged slice, directly reconstructs the energy density distribution in this plane.¹⁶

Due to the imperfections of the real sensor, contributions from sources outside the focus plane, i.e., at $z\ne 0$, also have to be taken into account. Furthermore, for sources at $z=0$, but $s\ne R$, the linear dependence between $t$ and $s$ is not valid. To include all these properties of the real sensor, we describe the forward problem by introducing a system matrix that links the temporal signals $q_{\varphi }(t)$ measured by the detector to the distribution of absorbed energy density $W_{\varphi }(s,z)$. Since the imaging problem focuses on discrete sections, the $z$-position of the sources is discretized in multiples of $\mathrm{\Delta }z$, where $\mathrm{\Delta }z$ is the increment between individual sections.

For the implementation of the iterative image reconstruction algorithm, a discrete representation of the source distribution $w_{\varphi }(x,y,z)$ is required. We use a representation in uniform spherical expansion functions, as it has been done in similar PA imaging algorithms before.^10,17 The resulting system matrix is obtained by calculating the response of the detector to single, homogenously heated spheres residing at grid positions, having a radius on the order of the grid spacing. A response to a single spherical source is calculated by summing the analytical signals received by an infinite line¹⁴ over all line positions (typically 200) forming the detector surface.

The system matrix is first defined in three dimensions as $a_{\mathit{k},\mathit{l},|\mathit{m}-\mathit{n}|}$. It contains the Abel transform and links signals measured at time $t_k$ to sources located at distance $s_l$ in the plane at $z_n=n\mathrm{\Delta }z$ in the object, which has been shifted in $z$-direction by distance $z_m=m\mathrm{\Delta }z$. In our notation, where the sample is rotated and shifted relative to the coordinate system of the sensor, $z_n$ is defined relative to the unshifted sample (dashed contour in Fig. 2). Both $m$ and $n$ are ranging from $-N$ to $+N$. A signal measured by the cylindrical detector at time $t_k=k\mathrm{\Delta }t$, shift position $z_m$, and angular orientation ${\varphi }_j$ is denoted by $q_{j,k,m}\equiv q_{{\varphi }_j}(t_k,z_m)$ and is related to the 2-D initial energy density distribution $W_{j,l,n}\equiv W_{{\varphi }_j}(s_l,z_n)$ by

Due to the symmetry of the sensor, the system matrix depends only on $|m-n|$. According to the discrete representation of the source distribution $w_{\varphi }$, the meaning of $W_{j,l,n}$ is actually the strength of a uniform, spherical voxel located at position $(s_l,z_n)$ in the object rotated by ${\varphi }_j$, which is integrated over an infinite line in direction of the cylindrical detector axis. The elements of the system matrix are thus given by

Here, $h(b,d_{l,|m-n|},t)$ is the analytical signal of a homogenously heated spherical source of radius $b$ (on the order of the spatial resolution), located at $s_l$ in the plane offset by $(m-n)\mathrm{\Delta }z$ from the focus plane of the detector, integrated over a line at position $(s^{\prime },z^{\prime })$ on the detector surface¹⁴ (see Fig. 2). Integration over the detector surface is described by the integral of this signal over the curve $\mathrm{C}(s^{\prime },z^{\prime })$, the circular arc defining the cross-section of the cylindrical detector surface.

The above relation [Eq. (4)] contains one equation for each rotation angle ${\varphi }_j$. The system matrix $a_{k,l,|m-n|}$ is written here in three dimensions. Alternatively, all recorded signals at one angular orientation ${\varphi }_j$ can be written as a vector ${\mathbf{q}}^{(j)}$ with elements $q_u^{(j)}$ and all sources in another vector ${\mathbf{W}}^{(j)}$ with elements $W_v^{(j)}$. The single layers of $a_{k,l,|m-n|}$ (each of them having a different $m-n$), denoted by 2-D matrices ${\mathbf{a}}^{(|m-n|)}$, can then be arranged in a 2-D matrix $\mathbf{A}$ with elements $A_{u,v}$, yielding a more convenient way to write this relation.

The arrangement of submatrices ${\mathbf{a}}^{(|m-n|)}$ is

2.2.

Reconstruction

The goal of image reconstruction is to find for each projection angle a vector ${\mathbf{W}}^{(j)}$ that minimizes the deviation of ${\mathbf{q}}^{(j)}=\mathbf{A}{\mathbf{W}}^{(j)}$ from the measured signals ${\widehat{\mathbf{q}}}^{(j)}$. Once these values are determined, the final step in image reconstruction is to calculate the source distribution in all recorded sections by applying the inverse Radon transform to ${\mathbf{W}}^{(j)}$. For recovering ${\mathbf{W}}^{(j)}$ from ${\widehat{\mathbf{q}}}^{(j)}$, we use the maximum likelihood-expectation maximization (ML-EM) algorithm,^18–20 which is known to provide solutions for problems with incomplete or noisy data. It furthermore preserves positivity of the solution, which is important for reconstructing a physically meaningful ${\mathbf{W}}^{(j)}$. The calculation of the iterations requires an adjoint system matrix, which has the effect of a back-projection operator and is the transpose of the real-valued matrix $\mathbf{A}$. However, we observed a more rapid convergence of the algorithm by using the following procedure for calculating a modified back-projection operator $\mathbf{B}$. In each submatrix ${\mathbf{a}}^{(|m-n|)}$, we seek the maximum value of each column and store its index kmax. This gives the most likely arrival time $t_{\mathrm{kmax}}$ of a wave generated at position $s_l$ in the $|m-n|$’th neighbor plane. “Most likely” in this context means, for instance, that a wave generated by a point source at this position impinges perpendicularly at the sensor surface, giving rise to the maximum detected signal by this source, even if it is not lying on the focal line. Then we build a new matrix as in Eq. (8) that only contains the column wise maxima of ${\mathbf{a}}^{(|m-n|)}$. Transposing this matrix yields the back-projection operator $\mathbf{B}$ with elements $B_{v,u}$, which now links the temporal signals to their most likely source positions. The product $\mathbf{B}$ ${\mathbf{q}}^{(j)}$ then yields a back-projection reconstruction of the data ${\mathbf{q}}^{(j)}$, similar to the focal line concept proposed by Xia et al.²¹ The iterative reconstruction algorithm for obtaining the ($i+1$)’th estimate of element ${\mathrm{W}}_v^{(j)}$ of vector ${\mathbf{W}}^{(j)}$ is given by

To avoid a division by zero, a small scalar $\lambda$ on the order of one or two percent of the maximum of $\mathbf{A}$ ${\mathbf{W}}^{(j),i}$ is added to the denominator on the right-hand side. This iteration is performed for each rotation angle ${\varphi }_j$ separately using the temporal signals of all the recorded sections. It has to be taken into account that the projections ${\mathbf{W}}^{(j)}$ are not independent of each other since they emanate from the same distribution of absorbed energy. Therefore, at each iteration, the projections ${\mathbf{W}}^{(j)}$ are normalized in a way that their integral over distance $s$ gives the same value for each section, i.e., for each $z_n$. This integral can be regarded as a measure of the total energy absorbed in this section and should be equal for all projections. Convergence of the algorithm is monitored by calculating the norm of the residual ${\widehat{\mathbf{q}}}^{(j)}-{\mathbf{AW}}^{(j),i+1}$. After completion of the iterations, the final distribution of the energy density in each of the $z_n$ sections is obtained by applying the inverse Radon transform to the corresponding subsets of ${\mathbf{W}}^{(j)}$.

For comparison, we also used an alternative method for the reconstruction of ${\mathbf{W}}^{(j)}$ from the measured pressure signals ${\widehat{\mathbf{q}}}^{(j)}$ using the lsqr function implemented in MATLAB®. This method, which is suited for large and sparse system matrices, provides a least squares inversion of matrix $\mathbf{A}$.

2.3.

Experiment

The model-based reconstruction was tested with a detector made of a 110-μm-thick piezoelectric polyvinylidene fluoride (PVDF) film (Measurement Specialties, Hampton, Virginia) glued onto the concave cylindrical surface of a polyvinyl chloride backing. Detector parameters were a radius of curvature of $R=30\mathrm{mm}$, a width of $x_d=50\mathrm{mm}$, and a height of $z_d=30\mathrm{mm}$, giving a numerical aperture of $\mathrm{NA}=0.5$. A thin plastic sheet separated the PVDF film from a water bath, where the sample was placed. Illumination was provided from two sides parallel to the cylinder axis of the detector using pulses from a frequency-doubled, Q-switched Nd:YAG laser (Spectron Laser Systems, Rugby, UK) with a pulse duration of 10 ns.

In a first experiment, signals from a small absorber were acquired for comparison with the calculated system matrix. A small, spherical oil droplet with a radius of 150 μm containing black dye was embedded in Agar gel (2% by weight) and was scanned in the $s\text{-}z$-plane perpendicular to the cylinder axis. A range from $z=-3\mathrm{mm}$ to $z=3\mathrm{mm}$ (the position of the focal line of the cylindrical detector is at $z=0$) and from $s=R–6\mathrm{mm}$ to $\mathrm{s}=R+6\mathrm{mm}$ was scanned with increments of $\mathrm{\Delta }z=200\mu \text{m}$ and $\mathrm{\Delta }s=200\mu \text{m}$.

For a tomography experiment, a phantom consisting of clear Agar gel was prepared, containing black oil droplets in two layers $\sim 2\mathrm{mm}$ apart. During data acquisition, the phantom was scanned to 17 $z$-positions with an increment of 200 μm. At each of the $z$-positions, the phantom was rotated by 360 deg with an angular increment of 1.8 deg. For the reconstruction, the increment in $s$-direction was set to 50 μm. Since it is important to obtain signals from all absorbers at each position of the detector, the phantom was illuminated with a wide beam that covered the whole phantom.

In the experiments it turned out that sometimes there remained some negative signal values after the Abel transform. To obtain the purely positive signal necessary for the reconstruction algorithm, the following procedure was used. First, the measured signals were Abel transformed. Then an inverse Radon transform was calculated, yielding a first image of $\mathbf{W}$, which contained both positive and unphysical negative components. In this image, the negative parts were set equal to zero, and to the resulting image, the Radon transform was applied. This gave purely positive signals $\hat{\mathbf{q}}$. Compared to simply removing negative values from the Abel-transformed signals, this procedure preserved more of the information content of the data because it avoided the deletion of features residing on negative parts of the signals.

3.1.

Simulations

Reconstructions from simulated data are shown in Fig. 3. The same phantom was used as in Fig. 1. Again, the central plane with three spheres is shown. For the calculation, signals from five planes were used: the three planes containing the spherical absorbers and two additional planes at a distance of 1 mm below and above the phantom. Figure 3(a) shows the result of a back-projection, applying the operator $\mathbf{B}$ to the data, and Fig. 3(b) shows the result after 20 iterations of the ML-EM algorithm in Eq. (9). A reconstruction using 20 iterations of the lsqr algorithm is shown in Fig. 3(c). The back-projection operator leads to a slight improvement of the image quality compared to a direct reconstruction from the Abel-transformed signals [Fig. 1(c)], but only the model-based algorithms are able to eliminate the artifacts. Comparing ML-EM and lsqr methods, both perform similarly in terms of reducing cross-talk and in-plane distortions. The latter, however, yields some residual negative values, as shown in Fig. 3(e). Here, a horizontal profile through one of the sources is displayed together with the reconstruction that is obtained by applying the inverse Radon transform to Abel-transformed pressure signals. No negative values are visible in the profile of the ML-EM reconstruction shown in Fig. 3(d) due to the embedded positivity constraint of this algorithm. The profiles also show a slightly better suppression of cross-talk with the ML-EM algorithm than with the lsqr method. This is seen at position $x=3\mathrm{mm}$, where the reconstruction from the raw, Abel-transformed signals shows high image values on the order of 0.6. At the same position, the residual peak after lsqr reconstruction reaches a value of $\sim 0.1$, whereas the ML-EM reconstruction yields a negligible value on the order of 0.05. A number of merit for the comparison of the algorithms is the total relative residual, defined as the norm of ${\widehat{\mathbf{q}}}^{(j)}-{\mathbf{AW}}^{(j),i+1}$ divided by the norm of ${\widehat{\mathbf{q}}}^{(j)}$. This residual is calculated for the entire measurement data, including all projection angles ${\varphi }_j$. Assuming noise-less data as in the shown simulations, the values for the total relative residual are 0.04 for lsqr and 0.11 for ML-EM (both after 20 iterations). Here the lsqr algorithm shows better convergence, although other parameters, such as the reduction of cross-talk, indicate better performance of the ML-EM algorithm. We repeated the simulation with additive Gaussian noise with a standard deviation of 5% of the signal maximum and found relative residuals of 0.30 for lsqr and 0.13 for ML-EM. With these more realistic data the better performance of the ML-EM algorithm in the presence of noise is demonstrated.

Fig. 3

Simulation results showing section images of the phantom displayed in Fig. 1(a). (a) Applying back-projection operator $\mathbf{B}$ to the data. (b) Twenty iterations with the maximum likelihood-expectation maximization (ML-EM) algorithm. (c) lsqr reconstruction, all followed by the inverse Radon transform. In (a) to (c) negative values are not displayed. (d) Normalized profiles through the spherical source at $(-3,-3)\mathrm{mm}$, comparing a reconstruction from Abel-transformed raw signals [seen in Fig. 1(c)] with the ML-EM reconstruction. (e) Comparison of profiles through the lsqr reconstruction, again compared to the Abel-transformed raw signals.

3.2.

System Matrix

Experimental and theoretical signals from a small sphere scanned in the $s\text{-}z$-plane are displayed in Fig. 4. The theoretical signals were obtained by integrating the analytical signals of a sphere with 150 μm radius over the detector surface. In Figs. 4(a) and 4(b) the signals for $m=n$ are shown and in Figs. 4(d) and 4(e) the signals for $|m-n|=6$. In the latter case, the distance between the focus plane and the source was $6\mathrm{\Delta }z=1.2\mathrm{mm}$. There is a good agreement between experiment and calculation, showing that the parameters for the simulation of the measurement were properly chosen. This justifies the use of the calculated signal matrix, which was obtained by applying the Abel transform to the columns of the matrix containing the simulated signals. The submatrices with $m=n$ and $|m-n|=6$ are displayed in Figs. 4(c) and 4(f). Matrix ${\mathbf{a}}^{(0)}$ is close to diagonal with a maximum at the focus, which is located at $s=30\mathrm{mm}$, and shows some spreading of values at positions outside the focus. Matrix ${\mathbf{a}}^{(6)}$ has highest values for $s$-positions outside the focus.

Fig. 4

Detected signals, system, and back-projection matrices for $z=0$ (upper row) and $z=1.2\mathrm{mm}$ (lower row). (a) and (d) Experimental signals for a black oil droplet with 0.15 mm radius. (b) and (e) Simulated signals for an absorbing sphere with 0.15 mm radius. (c) System submatrix ${\mathbf{a}}^{(0)}$. (f) System submatrix ${\mathbf{a}}^{(6)}$, both calculated from simulated signals by applying the Abel transform.

3.3.

Phantom Experiment

For image reconstruction, signals from all 17 planes were used, covering the range with the black dye droplets and some extra sections that contained no absorbers. With the spatial increment of 200 μm, the total $z$-range used for reconstruction was 3.2 mm. Reconstructed images in Fig. 5 show one of the sections containing absorbing spherical droplets. Additionally, maximum amplitude projections in $x$-direction reveal the spreading of the imaged objects in direction perpendicular to the sections. A comparison of three reconstruction methods is provided: direct reconstruction applying the inverse Radon transform to Abel-transformed raw data, ML-EM reconstruction, and lsqr inversion. In all images, only positive values are displayed. To quantify the resolution achievable with the different reconstruction methods, Fig. 6(a) compares profiles in $z$-direction through one of the small spheres seen at position $x=3.5\mathrm{mm}$, $y=-3\mathrm{mm}$, $z=1\mathrm{mm}$. This sphere has subresolution size in $z$-direction and is suited for estimating the out-of-plane resolution. In the reconstruction from raw, Abel-transformed pressure signals, the spreading in $z$-direction is largest and amounts to 1.35 mm (FWHM). The two iterative techniques yield better resolution, with 0.86 mm using lsqr and 0.51 mm using ML-EM. Also, a significant improvement of resolution in radial direction becomes evident when plotting the radial intensity profiles of the same sphere, as shown in Fig. 6(b). The radial direction goes through the point (0,0) in the $x-y$-sections displayed in Figs. 5(a) to 5(c).

Fig. 5

Experimental imaging. Sections in an $x-y$-plane [(a) to (c)] and maximum amplitude projections along $x$-direction [(d) to (f)] from Abel-transformed raw signals [(a) and (d)], using ML-EM [(b) and (e)], and using lsqr [(c) and (f)].

Fig. 6

Profiles across the spherical source located at $x=3.5\mathrm{mm}$, $y=-3\mathrm{mm}$, $z=1\mathrm{mm}$ along (a) $z$-direction and (b) radial direction.