1.

Introduction

The capabilities of biomedical photoacoustic (PA) imaging were remarkably enhanced by significant advances in computer technology and the use of confocal and toroidal piezoelectric transducer arrays.^1–7 The complexity of design of PA and laser ultrasonic (LU) imaging systems requires the development of criteria and techniques for their assessment at the design stage and calibration after manufacturing. The main characteristics of any imaging system are the size of sensitivity region (the region in space in which a signal source can be reliably recognized), spatial resolution (the minimum distance between distinguishable objects in an image), and depth of field (the size of the region within which high spatial resolution is achieved). In designing biomedical PA and LU imaging systems, it is important to maximize the size of the sensitivity region and depth of field as well as to improve the spatial resolution.

Real-time diagnostics of biological objects and nondestructive testing often require two-dimensional (2-D) visualization of the object under study. Reconstruction and analysis of three-dimensional (3-D) images in real time is a separate complex task, requiring considerable computing power. Real-time PA and LU image reconstruction is usually based on the backprojection algorithm,⁸ which can be effectively parallelized on graphics processors.⁹ (Real-time operation mode precludes the use of computationally intensive model-based algorithms, which can compensate for distortions caused by the finite size of receivers.^10,11) In 2-D imaging, the array should be designed to receive signals from sources in the image plane. To form a sufficiently “thin” image plane, the array is focused in the direction perpendicular thereto.¹² Focusing the array in the image plane itself allows the field of view to be increased and spatial resolution to be improved^13,14 at the expense of the smaller sensitivity region. Indeed, the smaller the size of an ellipsoidal Gaussian PA source, the wider its directivity pattern in the corresponding direction;¹³ that is why the greater the angular aperture of the array for a given point of the image plane, the higher the spatial resolution at this point.¹⁴ In the case of a planar array, the field of view is limited to the angle of total internal reflection of ultrasound at the boundary between the immersion fluid and the array material. Toroidal arrays are focused in both the image plane and the direction perpendicular thereto. Thus, real-time 2-D imaging systems equipped with such arrays appear promising for investigation of biological objects,^5–7,15 femtosecond laser filaments,¹⁶ etc., where high spatial resolution in a limited spatial region is required.

The limited view issue has been widely discussed in PA imaging. In Ref. 17, the “visibility” condition is formulated. If there exists a point $P$ on the sharp boundary $L$ of the object under study for which the normal to $L$ does not pass through any receiver, then the image of $L$ will be mandatorily blurred away at $P$ (“invisible”). In scenarios where enclosing the object in the detection surface and achieving full view are possible (for instance, small animal imaging^18,19 or imaging of femtosecond laser filaments¹⁶) using full-ring detection is advantageous because it eliminates partial-view artifacts and mitigates image blurring. However, in many realistic PA imaging scenarios, limited view effects are unavoidable.²⁰ Beyond increasing acoustic coverage, multiple images acquired in various conditions can be superposed to enhance image quality. The object under study or the sensor array can be rotated, and images acquired for various rotation angles can be averaged.²¹ Multiple images of random sparse distributions of artificially created small PA sources within “invisible” structures can be nonlinearly combined.²⁰ However, approaches suggesting acquisition of multiple images can mitigate partial-view artifacts at the expense of longer imaging time which may be problematic for moving objects. In this study, we assume the single-impulse acquisition scenario for real-time PA imaging.

The size of the sensitivity region, depth of field, and spatial resolution of arrays are determined by their geometry, sizes of receivers, and number of receivers. To obtain the spatial resolution in a given point in space, the Rayleigh integral¹² is usually calculated and point spread function (PSF) is reconstructed. In Refs. 12, 2223.–24, for arrays focused only on the image plane, the “thickness” of the latter was estimated using sensitivity maps constructed for a single receiver in a plane perpendicular to the image plane using the amplitudes of the received pressure signals. For such arrays, the relationship among spatial resolution in focus and field of view, number of receivers, and their sizes was studied thoroughly. Also, the dependence of the focal resolution on the angular aperture of the arrays focused only in the image plane was addressed.^14,25 The main disadvantage of the latter type of arrays is a “thick” image plane. In Ref. 7, the sensitivity fields for spherical and toroidal arrays in the image plane were calculated by rotating and summing single-element sensitivity fields, and sensitivity was defined as the maximum pressure signal amplitude from a given position. Although the aforementioned techniques have been used to build experimental sensor arrays and should be considered for quick performance analysis, the ever-growing power of modern computers makes more computationally intensive methods applicable. As the final result of the imaging process is an image, we suggest that it is the image and its features that need to be quantified and used for assessing and comparison.

In this study, it is proposed that the size of the sensitivity region, depth of field, and spatial resolution of toroidal arrays for real-time PA imaging should be estimated using sensitivity maps and spatial resolution maps in the image plane, constructed for the sensor array as a whole using the backprojection algorithm. The technique for construction of such maps is described in detail and the results of the computational research are discussed. The numerical estimates of the size of the sensitivity region and spatial resolution are given for arrays with different apertures and receivers of various sizes. The relationship between the size of the sensitivity region and the spatial resolution is studied for real-time image visualization systems with toroidal arrays. Sensitivity maps and spatial resolution maps might help in optimization of array geometry and the number and size of receiving elements.

2.

Construction of Sensitivity Maps and Spatial Resolution Maps

In PA imaging, the pressure field $p(\overrightarrow{r},t)$ in an acoustically homogeneous medium is the result of absorption of a short laser pulse. Due to the small duration of the laser pulse, the thermal energy density released at point $\overrightarrow{r}$ may be considered proportional to the Dirac delta function $\delta (t)$ of time $(t)$. Neglecting viscosity and thermal conductivity of the medium, wave equation for $p(\overrightarrow{r},t)$ in this case can be written as^1,26

Here $p_0(\overrightarrow{r})=p(\overrightarrow{r},t=0)$ is the initial pressure distribution in the medium, $\mathrm{\Delta }$ is the Laplace operator, and $c$ is the speed of sound. In the case of PA point source with amplitude $P_0$ [$p_0(\overrightarrow{r}|{\overrightarrow{r}}_0)=P_0\delta (\overrightarrow{r}-{\overrightarrow{r}}_0)$] placed at point ${\overrightarrow{r}}_0$, the solution $p^{pt}(\overrightarrow{r}|{\overrightarrow{r}}_0,t)$ of Eq. (1) has the following form:

In order to construct sensitivity maps and spatial resolution maps correctly, the finite size of the receiver with area $\mathrm{\Delta }S$ and its frequency response $\tilde{H}(\omega )$ should be taken into account. The Fourier transform of pressure ${\tilde{p}}_d^{pt}({\overrightarrow{r}}_S|{\overrightarrow{r}}_0,\omega )$ measured by the receiver located at the point defined by vector ${\overrightarrow{r}}_S$ is equal to the product of the Fourier transform of pressure [Eq. (2)]

In Eq. (4), the receiving element has the shape of a segment of a torus, a cylinder, or a sphere. The Fourier transform of pressure ${\tilde{p}}_d^{pt}({\overrightarrow{r}}_S|{\overrightarrow{r}}_0,\omega )$ in our calculation was approximated by the sum

Fig. 1

The typical example of toroidal array.

Estimation of the sensitivity region and spatial resolution of the imaging system is based on the analysis of $\mathrm{PSF}(\overrightarrow{r}|{\overrightarrow{r}}_0)$ (point spread function), which is defined as an image of PA point source reconstructed by the system. $\mathrm{PSF}(\overrightarrow{r}|{\overrightarrow{r}}_0)$ functionally depends on the image spatial coordinates $\overrightarrow{r}$ and parametrically depends on the coordinates of the point source ${\overrightarrow{r}}_0$ and it is considered to be the fundamental characteristic feature of an image in theoretical models of PA and LU imaging formation. $\mathrm{PSF}(\overrightarrow{r}|{\overrightarrow{r}}_0)$ is also influenced by a number of factors. Loss of high-frequency components of the pressure signal due to acoustic attenuation in the medium and limited bandwidth of receivers lead to image blurring.^10,16,27 PSF is greatly influenced by the receiver size as the signal is averaged over its surface.^10,27 If the number of receivers is rather small, the image exhibits severe artifacts that can be interpreted as nonexistent sources.¹⁰

PSF essentially depends on the reconstruction algorithm. If the pressure signal from the PA point source $p^{pt}({\overrightarrow{r}}_S|{\overrightarrow{r}}_0,t)$ is known at every point determined by ${\overrightarrow{r}}_S$ on the detection surface $S$ that subtends the solid angle $\mathrm{\Omega }$, then the initial pressure distribution $p_0(\overrightarrow{r}|{\overrightarrow{r}}_0)$ in the medium can be reconstructed using the backprojection algorithm⁸

Here ${\overrightarrow{n}}_S$ is the normal to surface $S$ at the point specified by vector ${\overrightarrow{r}}_S$. The reconstructed initial pressure distribution $p_0(\overrightarrow{r}|{\overrightarrow{r}}_0)$ in Eq. (6) is the PSF for an ideal system with infinite number of point-like receivers with infinite bandwidth. In our study, $\mathrm{PSF}(\overrightarrow{r}|{\overrightarrow{r}}_0)$ was reconstructed from $N$ pressure signals $p_d^{pt}({\overrightarrow{r}}_{Sj}|{\overrightarrow{r}}_0,t)$ measured by each of $N$ identical receivers located on the surface $S$ at the points defined by vectors ${\overrightarrow{r}}_{Sj}$ ($j=\mathrm{1,2},\dots ,N$) using the formula⁸

Note that, unlike $p_0$, PSF can take negative values due to the finite number of receivers.

Functional dependence of PSF on image spatial coordinates $\overrightarrow{r}$, when ${\overrightarrow{r}}_0$ is fixed, enables its maximum amplitude

Although 3-D maps provide complete information required to assess and compare different array configurations, construction of 3-D maps is a computationally intensive consumable task. We focus on a current study investigation of 2-D slices of 3-D sensitivity and spatial resolution maps in the image plane ($y_0=0$): $I(x_0,z_0)$, $A(x_0,z_0)$, and $L(x_0,z_0)$ for toroidal arrays and their dependence on array parameters. All calculations were done in MATLAB software package (MathWorks) on AMD FX-9590 CPU @ 4.7 GHz (AMD) and accelerated on NVIDIA GeForce GTX 780 GPU (NVIDIA).

Figure 2 shows a typical example of $\mathrm{PSF}(x,z|x_0,z_0)$ with corresponding maximum amplitude $I(x_0,z_0)$, axial $A(x_0,z_0)$, and lateral $L(x_0,z_0)$ spatial resolutions (a) and the typical 3-D sensitivity map $I({\overrightarrow{r}}_0)$ (b) for toroidal array with $R=60\mathrm{mm}$, $f=40\mathrm{mm}$, $N=16$, $d=1\mathrm{mm}$, $\phi =45\mathrm{deg}$, and $\theta =23\mathrm{deg}$ as shown in Fig. 1, where $D$ denotes the FWHM in the lateral $z$-direction of the 2-D sensitivity map in the image plane $xz$.

Fig. 2

The typical example of $\mathrm{PSF}(x,z|x_0,z_0)$ for fixed point $(x_0,z_0)$ in the image plane ($y_0=0$), axial $A(x_0,z_0)$, and lateral $L(x_0,z_0)$ spatial resolutions in this case (a) and typical 3-D sensitivity map $I({\overrightarrow{r}}_0)$ (b) for toroidal array as shown in Fig. 1 with $R=60\mathrm{mm}$, $f=40\mathrm{mm}$, $N=16$, $d=1\mathrm{mm}$, $\phi =45\mathrm{deg}$, and $\theta =23\mathrm{deg}$.

2-D sensitivity and spatial resolution maps were constructed for cylindrical ($R=\infty$), spherical ($R=f$), and toroidal ($R\ne f$) arrays with different apertures and with up to 64 identical receivers of various sizes (see Fig. 1), where $R$ is the curvature radius in the image plane and $f$ is the curvature radius in the direction perpendicular thereto. Centers of curvature (foci of individual receivers) in the direction perpendicular to the image plane formed the focal line. Toroidal arrays also had two angular apertures: $\phi$ (in the image plane) and $\theta$ (in the direction perpendicular to the image plane). Each receiver had the same width $d$ in the image plane. 2-D sensitivity maps and spatial resolution maps were constructed in the $xz$ image plane when the point $(x_0,y_0,z_0)$ belonged to the square region ($32\times 32{\mathrm{mm}}^2$) determined by four points $(f-R\pm 16\mathrm{mm},0,0\pm 16\mathrm{mm})$. A single PA point source was placed in turn at every point $(x_0^i,z_0^j)$ of a uniform $129\times 129$ rectangular grid, and pressure signals were calculated at all the receivers of the array. Pressure on the surface of every receiver was recorded in $20\times 200$ points. Then, using Eq. (7), 2-D image of $\mathrm{PSF}(x,z|x_0^i,z_0^j)$ was reconstructed for points $(x,z)$ belonging to the square region ($3\times 3{\mathrm{mm}}^2$) determined by four points $(x_0^i\pm 1.5\mathrm{mm},z_0^j\pm 1.5\mathrm{mm})$ with a resolution of $300\times 300\text{pixels}$.

Having reconstructed $\mathrm{PSF}(x,z|x_0^i,z_0^j)$, its maximum amplitude, $I(x_0^i,z_0^j)$, was determined according to Eq. (9). To find axial and lateral resolutions, which were defined as FWHM of PSF in two mutually perpendicular directions, PSF was thresholded: in all points of the image, the values of PSF smaller than $0.5I(x_0^i,z_0^j)$ were replaced with zeros and the values of PSF greater than $0.5I(x_0^i,z_0^j)$ were replaced with ones. The resulting binary image was approximated by an ellipse with the same second moments as the calculated image using the MATLAB regionprops function. The minor and major axes of the ellipse were assumed to be equal to axial $A(x_0^i,z_0^j)$ and lateral $L(x_0^i,z_0^j)$ resolutions, respectively.

3.

Results and Discussion

Figure 3 shows normalized sensitivity maps and spatial resolution maps for spherical arrays with $N=32$ receiving elements, $R=f=40\mathrm{mm}$, $\theta =23\mathrm{deg}$, and various receiver widths $d$ (0.5 and 1.0 mm) and angular apertures in the image plane $\phi$ (45 deg and 90 deg). Here and below in the maps of axial $A(x_0,z_0)$ and lateral $L(x_0,z_0)$ resolution regions with high resolution are shown in red (60 to $80\mu \mathrm{m}$ for axial resolution and 200 to $400\mu \mathrm{m}$ for lateral resolution); axial resolution less than $150\mu \mathrm{m}$ is shown in blue. All maps of the same type (in each row) are shown in the same unified color scale.

Fig. 3

2-D normalized sensitivity maps $I(x_0,z_0)$, axial resolution maps $A(x_0,z_0)$, and lateral resolution maps $L(x_0,z_0)$ for spherical arrays with $R=f=40\mathrm{mm}$, $N=32$, $\theta =23\mathrm{deg}$ under different $d$ and $\phi$.

We have found that sensitivity maps and spatial resolution maps depend only weakly on the number of receivers for $N\ge 16$ (see Fig. 4). Earlier,¹² it was found that focal resolution did not depend on $N$ for $N>8$. The map calculation time is roughly proportional to $N$. Notwithstanding visual similarities of the sensitivity maps $I(x_0,z_0)$, the size of the high-sensitivity region essentially depends on $N$. Real experimental images are usually thresholded at some level to reduce clutter. The level of backprojection artifacts is proportional to $1/N$, and the signal-to-noise ratio of experimental images is proportional to $\sqrt{N}$. Consequently, the threshold level can be decreased for high enough $N$, and the size of high-sensitivity region increases with $N$ while its shape remains virtually unchanged.

Fig. 4

2-D sensitivity maps $I(x_0,z_0)$, axial resolution maps $A(x_0,z_0)$, and lateral resolution maps $L(x_0,z_0)$ for spherical arrays with $R=f=40\mathrm{mm}$, $\theta =23\mathrm{deg}$, $\phi =90\mathrm{deg}$, $d=0.5\mathrm{mm}$, and various number of receivers $N$.

The size of high-sensitivity regions, spatial resolution, and depth of field are essentially dependent on the parameters of the array. We have found that high axial resolution in the high-sensitivity region is achieved for all the spherical array configurations presented. To compare various array configurations, we denote the FWHM of the sensitivity region by $D$ (see Fig. 1) and the FWHM of PSF (lateral resolution) at the point with maximum sensitivity by $\tilde{L}$. For spherical arrays, $D$ and $\tilde{L}$ depend on the receiver width $d$ and angular aperture $\phi$ in the image plane (see Fig. 5).

Fig. 5

The dependencies of $D$ and $\tilde{L}$ on receiver width $d$ when $\phi =45\mathrm{deg}$ (a) and angular aperture $\phi$ when $d=1\mathrm{mm}$ (b) for spherical arrays with $R=f=40\mathrm{mm}$, $N=16$. Squares are calculations, dashed lines represent approximations.

Arrays with cylindrical geometry provide the widest region of high sensitivity and relatively poor lateral resolution while spherical arrays have smaller sensitivity regions and significantly higher lateral resolution. Figure 6 shows that toroidal arrays ($R\ne f$) with different curvature radii provide an opportunity to optimize the size of high-sensitivity region and the value of lateral resolution.

Fig. 6

2-D sensitivity maps $I(x_0,z_0)$, axial resolution maps $A(x_0,z_0)$, and lateral resolution maps $L(x_0,z_0)$ for toroidal arrays with $f=40\mathrm{mm}$, $N=16$, $\theta =23\mathrm{deg}$, $\phi =45\mathrm{deg}$, $d=1\mathrm{mm}$, and various curvature radius in the image plane $R$.

As curvature radius in the image plane $R$ increases, the size of sensitivity region $D$ increases. If $R\to \infty$ (toroidal array approaches cylindrical shape), when $D$ strives for $D_{\infty }$, but the lateral resolution $\tilde{L}$ deteriorates up to ${\tilde{L}}_{\infty }$. Spherical arrays ($R=f$) provide the smallest high-sensitivity region $D(R=f)$ and the best lateral resolution $\tilde{L}(R=f)$. The calculated dependencies $D(R)$ and $\tilde{L}(R)$ are shown in Fig. 7. They can be approximated by the following formulas:

Fig. 7

(a) The dependencies of FWHM of sensitivity region $D(R)$ (top) and (b) lateral resolution at the point of maximal sensitivity $\tilde{L}(R)$ (bottom) on the curvature radius of the array $R$ for arrays with $N=16$, $\theta =23\mathrm{deg}$, $f=40\mathrm{mm}$, and various receiver widths $d$ and array angular apertures $\phi$.

4.

Conclusions

A toroidal array quantitative assessment and comparison procedure is proposed, which is required for design of real-time PA and LU imaging systems. Such imaging systems will make it possible to study the structure of biological objects, rocks, composite materials, and femtosecond laser filamentation with high spatial resolution in real time. The proposed approach, based on analysis of sensitivity maps and spatial resolution maps constructed using the backprojection algorithm, allows comprehensive assessment of the parameters of arrays to choose the most suitable configuration for a given problem.

Importantly, as toroidal detection geometries remain difficult to manufacture, the calculation procedure of sensitivity maps can be applied in their calibration process to improve imaging quality. Individual receivers can have different inherent sensitivities and geometrical sizes after manufacturing. Due to these imperfections, sensitivity maps constructed from experimental data may differ from numerically calculated ones. However, weight factors may be introduced in the backprojection algorithm for individual receivers to partially compensate for their imperfections. These weighting factors may be found by minimizing the difference between numerical (ideal) maps and experimental maps—a typical nonlinear optimization problem with well-developed solution algorithms.