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,8 which can be effectively parallelized on graphics processors.9 (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.12 Focusing the array in the image plane itself allows the field of view to be increased and spatial resolution to be improved13,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;13 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.14 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,16 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 on the sharp boundary of the object under study for which the normal to does not pass through any receiver, then the image of will be mandatorily blurred away at (“invisible”). In scenarios where enclosing the object in the detection surface and achieving full view are possible (for instance, small animal imaging18,19 or imaging of femtosecond laser filaments16) 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.20 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.21 Multiple images of random sparse distributions of artificially created small PA sources within “invisible” structures can be nonlinearly combined.20 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 integral12 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.
Construction of Sensitivity Maps and Spatial Resolution Maps
In PA imaging, the pressure field 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 may be considered proportional to the Dirac delta function of time . Neglecting viscosity and thermal conductivity of the medium, wave equation for in this case can be written as1,26
Here is the initial pressure distribution in the medium, is the Laplace operator, and is the speed of sound. In the case of PA point source with amplitude  placed at point , the solution 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 and its frequency response should be taken into account. The Fourier transform of pressure measured by the receiver located at the point defined by vector 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 in our calculation was approximated by the sumFig. 1). The centers of all receiving elements were situated on the circular arc with the center of curvature at the origin for toroidal and spherical arrays and on the straight-line segment for cylindrical arrays. In our calculation, , where , the cutoff frequency.
Estimation of the sensitivity region and spatial resolution of the imaging system is based on the analysis of (point spread function), which is defined as an image of PA point source reconstructed by the system. functionally depends on the image spatial coordinates and parametrically depends on the coordinates of the point source and it is considered to be the fundamental characteristic feature of an image in theoretical models of PA and LU imaging formation. 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.10
PSF essentially depends on the reconstruction algorithm. If the pressure signal from the PA point source is known at every point determined by on the detection surface that subtends the solid angle , then the initial pressure distribution in the medium can be reconstructed using the backprojection algorithm8
Here is the normal to surface at the point specified by vector . The reconstructed initial pressure distribution in Eq. (6) is the PSF for an ideal system with infinite number of point-like receivers with infinite bandwidth. In our study, was reconstructed from pressure signals measured by each of identical receivers located on the surface at the points defined by vectors () using the formula8
Note that, unlike , PSF can take negative values due to the finite number of receivers.
Functional dependence of PSF on image spatial coordinates , when 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 (): , , and 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 with corresponding maximum amplitude , axial , and lateral spatial resolutions (a) and the typical 3-D sensitivity map (b) for toroidal array with , , , , , and as shown in Fig. 1, where denotes the FWHM in the lateral -direction of the 2-D sensitivity map in the image plane .
2-D sensitivity and spatial resolution maps were constructed for cylindrical (), spherical (), and toroidal () arrays with different apertures and with up to 64 identical receivers of various sizes (see Fig. 1), where is the curvature radius in the image plane and 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: (in the image plane) and (in the direction perpendicular to the image plane). Each receiver had the same width in the image plane. 2-D sensitivity maps and spatial resolution maps were constructed in the image plane when the point belonged to the square region () determined by four points . A single PA point source was placed in turn at every point of a uniform rectangular grid, and pressure signals were calculated at all the receivers of the array. Pressure on the surface of every receiver was recorded in points. Then, using Eq. (7), 2-D image of was reconstructed for points belonging to the square region () determined by four points with a resolution of .
Having reconstructed , its maximum amplitude, , 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 were replaced with zeros and the values of PSF greater than 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 and lateral resolutions, respectively.
Results and Discussion
Figure 3 shows normalized sensitivity maps and spatial resolution maps for spherical arrays with receiving elements, , , and various receiver widths (0.5 and 1.0 mm) and angular apertures in the image plane (45 deg and 90 deg). Here and below in the maps of axial and lateral resolution regions with high resolution are shown in red (60 to for axial resolution and 200 to for lateral resolution); axial resolution less than is shown in blue. All maps of the same type (in each row) are shown in the same unified color scale.
We have found that sensitivity maps and spatial resolution maps depend only weakly on the number of receivers for (see Fig. 4). Earlier,12 it was found that focal resolution did not depend on for . The map calculation time is roughly proportional to . Notwithstanding visual similarities of the sensitivity maps , the size of the high-sensitivity region essentially depends on . Real experimental images are usually thresholded at some level to reduce clutter. The level of backprojection artifacts is proportional to , and the signal-to-noise ratio of experimental images is proportional to . Consequently, the threshold level can be decreased for high enough , and the size of high-sensitivity region increases with while its shape remains virtually unchanged.
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 (see Fig. 1) and the FWHM of PSF (lateral resolution) at the point with maximum sensitivity by . For spherical arrays, and depend on the receiver width and angular aperture in the image plane (see Fig. 5).
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 () with different curvature radii provide an opportunity to optimize the size of high-sensitivity region and the value of lateral resolution.
As curvature radius in the image plane increases, the size of sensitivity region increases. If (toroidal array approaches cylindrical shape), when strives for , but the lateral resolution deteriorates up to . Spherical arrays () provide the smallest high-sensitivity region and the best lateral resolution . The calculated dependencies and are shown in Fig. 7. They can be approximated by the following formulas:Fig. 7).
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.
The authors have no relevant financial interests in the paper and no other potential conflicts of interest to disclose.
This work was carried out with financial support from the Ministry of Education and Science of the Russian Federation in the framework of the Increase Competitiveness Program of NUST “MISiS” (No. K1-2015-025) and the Russian Science Foundation (Grant No. 16-17-10181).
Anton S. Bychkov graduated Lomonosov Moscow State University, Russia, in 2016. At present, he is a PhD student in the Laser Optoacoustic Laboratory at Lomonosov MSU and an engineer in the Laser Ultrasonic Nondestructive Testing Laboratory at the National University of Science and Technology MISiS, Russia. His current research interests include combined real-time laser photoacoustic and laser-induced ultrasound imaging of biological objects, nondestructive testing of industrial products, and experimental research automation.
Elena B. Cherepetskaya is a professor at the National University of Science and Technology MISiS, Moscow, Russia. Her research has focused on optoacoustic spectroscopy of heterogeneous media, including rocks, concrete, and carbon composites on application of laser ultrasound in NDT and material evaluation.
Alexander A. Karabutov is professor at the International Laser Center of M.V. Lomonosov Moscow State University (ILC MSU), head of the laboratory of laser optoacoustics of ILC MSU and leading scientist of the Laser Ultrasonic Nondestructive Testing Laboratory at the National University of Science and Technology MISiS. Since 1980 he has been developing theoretical foundations and experimental research in the field of nonlinear acoustics and fundamentals of laser-ultrasonic spectroscopy.
Vladimir A. Makarov is a doctor of science, professor of physics, chair of general physics and wave processes, and the director of the International Laser Center in M. V. Lomonosov Moscow State University. His research interests are in the field of laser physics, nonlinear optics, interaction of laser radiation with matter, laser applications in life sciences, and fundamentals of laser-ultrasonic spectroscopy. He is an author of more than 200 refereed scientific publications, including some books.