Development of a Stationary 3D Photoacoustic Imaging System Using Sparse Single-Element Transducers: Phantom Study

Abstract

Photoacoustic imaging (PAI) is an emerging label-free and non-invasive modality for imaging biological tissues. PAI has been implemented in different configurations, one of which is photoacoustic computed tomography (PACT) with a potential wide range of applications, including brain and breast imaging. Hemispherical Array PACT (HA-PACT) is a variation of PACT that has solved the limited detection-view problem. Here, we designed an HA-PACT system consisting of 50 single element transducers. For implementation, we initially performed a simulation study, with parameters close to those in practice, to determine the relationship between the number of transducers and the quality of the reconstructed image. We then used the greatest number of transducers possible on the hemisphere and imaged copper wire phantoms coated with a light absorbing material to evaluate the performance of the system. Several practical issues such as light illumination, arrangement of the transducers, and an image reconstruction algorithm have been comprehensively studied.

1. Introduction

Photoacoustic (PA) imaging is a fast-developing hybrid imaging modality that combines optical absorption and ultrasound detection techniques to visualize the hemodynamic changes in biological tissues [1,2,3,4,5,6,7]. In contrast to all-optical imaging methods [8,9,10], it uses light for illumination and ultrasound for detection. In PA imaging (PAI), a nanosecond pulsed laser deposits energy onto a tissue, causing a local temperature increase of light absorbing molecules and subsequent thermal expansion through the thermoacoustic effect [1,6,7,11]. This thermal expansion causes a localized pressure rise, which propagates from the sample to be imaged by an ultrasound transducer (or an array of ultrasound transducers), similar to traditional ultrasound detection. Finally, depending on the configuration of the PAI system, a two-dimensional (2D) or a three-dimensional (3D) image is reconstructed from the absorption distribution across the tissue [12,13,14,15,16,17].

Among different implementations of PAI, photoacoustic computed tomography (PACT) has characteristics that help its clinical translatability. Some of these characteristics are: (i) deep penetration, compared to all-optical imaging methods such as Functional Near-Infrared Spectroscopy (fNIRS), where in both illumination and detection paths, travelling photons are involved, and compared to the use of acoustic detection in PACT, which has orders of magnitude smaller signal attenuation; (ii) high frame rate, compared to magnetic resonance imaging (MRI), PACT can generate faster frame-rate volumetric images. It is worth noting, in the current study, the channel data was averaged many times and worsened the temporal resolution. Ideally, PACTs are developed using high-sensitivity transducers and fast repetition-rate (rep-rate) lasers; the frame rate is defined by the rep-rate of the laser or at least one-tenth of that if 10 iterations of acquired signals are averaged at each channel. The ability of PACT systems to render 3D volumetric images, i.e., 3D-PACT, relies on covering angles to which the acoustic waves are traveling to [18]. Such an imaging configuration will produce adequate channel data for image reconstruction because reconstruction algorithms are developed based on the assumption that the entire imaging target is surrounded by transducers [19,20]. 3D-PACT has been studied for volumetric mammography, neuroimaging in small animal, and human neonates [7,20,21,22,23,24,25,26,27].

As for ultrasound detection implementation, one way to perform 3D imaging is to use a 2D ultrasound matrix array [28]. However, these arrays consist of many ultrasonic transducer elements, significantly increasing the system’s complexity and cost [28]. A cost-effective alternative method for 3D imaging is to use a linear array and stack multiple cross-sectional images; this method causes motion artifacts and usually it is not time efficient [29,30]. Apart from 2D arrays, expensive arc-shaped [31,32] or circular ring arrays [33,34,35] along with some mechanical translations have also been used for 3D-PACT; these acoustic detection configurations have not been optimized for the best temporal resolution and/or spatial resolution.

As for light illumination, an overhead illumination using a single large optical fiber or a full-field illumination using a fiber bundle, have been studied [15,36,37]. We have previously shown in [38,39] that full-field illumination provides a more homogenous and deeper illumination than overhead illumination; using full-field illumination, we have access to a larger illumination area, and according to the literature [40], the larger the area of illumination, the more light we can couple to the tissue. This issue is of paramount importance especially in the areas of brain and breast imaging where there is a need for a large amount of fluence inside the tissue, which is to be maintained below the tissue damage threshold.

Different 3D photoacoustic imaging systems with their detection and illumination characteristics as well as their applications are listed in

Table 1. Different 3D PA signal detection configuration.

Authors	Configuration (Element No.)	Illumination Configuration	Application	Ref.
Ephrat et al.	Spherically sparse array (15)	Center illumination from bottom guided by mirror	Imaging moving phantoms	[41]
Xiang et al.	Spherical array (192)	Cylindrically-shaped guided by fiber bundle	Monitoring needle-based drug delivery, hemodynamic changes and temperature variation	[42]
Dean-Ben et al.	Spherical array (256)	Custom-made optical fiber bundle through central opening	Measuring hemodynamic and oxygen parameters	[26]
Tang et al.	2D rectangular array (3 × 64)	Center illumination from top using liquid light guide (LLG)	Measuring hemodynamic responses in the primary visual cortex	[43]
Xia et al.	Linear array scanning (512)	Center illumination guided by fiber bundle	Hair phantom and ex-vivo mouse embryo imaging	[44]
Gateau et al.	Linear array rotate/translate scanning (128)	Side illumination using fiber bundle	In-vitro	[45]
Wang et al.	2D matrix array (50 × 50)	Top illumination using fiber bundle	Mapping of the sentinel lymph node in rat model	[20]
Wygant et al.	2D CMUT array (synthetic aperture: 64 × 64)	Side illumination using diffuser	Phantom study	[46]

. In the 3D-PACT systems given in Table 1, it can be seen that either a costly transducer array is used that has not solved the problem of limited detection view completely, or spatially equidistant low-cost single element transducers have been used with non-uniform light illumination techniques.

In this study, we have tested the feasibility of a stationary 3D hemispherical array PACT (HA-PACT) imaging system using 50 single element transducers with a full-field illumination scheme produced by a distributed optical fiber bundle. It is anticipated that the 3D HA-PACT mitigates the limited detection view problem by using the hemispherical ultrasound detection and improves the penetration depth by increasing the signal-to-noise level of the detected signal utilizing a full-field illumination scheme.

For implementation, initially we performed a simulation study with parameters close to those in practice to determine the relationship between the number of transducer elements and the quality of the reconstructed image. We then used the greatest number of transducers possible on the hemisphere surface and evaluated the performance of the resultant imaging system with two different illumination configurations, i.e., overhead and full-field illuminations, to image different imaging targets.

2. Materials and Methods

In this section, initially, a different number of elements in the hemispherical array and corresponding reconstructed images are simulated using k-wave toolbox in Matlab to justify the use of maximum packing of transducers into the hemisphere. This was done to determine if a reduced number of transducers, below the maximal packing method, would achieve satisfactory image quality. Next, the experimental setup and phantom preparation are explained, followed by a description of an image quality assessment protocol. The performance of different image reconstruction algorithms was evaluated on the simulation data. The phantom study was mainly designed to assess the performance parameters of the 3D HA-PACT and optimize the light illumination and the image reconstruction algorithm.

2.1. Simulation Study

To understand the relation between the quality of the reconstructed images and number of transducers in the 3D HA-PACT, we have modeled a hemispherical cap in k-wave toolbox [47]. The diameter of the hemisphere was chosen to be 12 cm based on the consideration of neonatal brain imaging as a potential application of the proposed system. Typically, the neonatal head circumference is between 20 cm to 35 cm [48]; we chose the maximum diameter and considered room for a thin layer of coupling medium inside the hemisphere. Based on the diameter of the transducer and the thickness of the designed hemisphere, 50 transducers are the maximum number that can be placed. Any further number of transducers would result in hollow space overlap in the hemisphere model. Additionally, we have simulated 14, 26, 38, and 50 equidistant transducer elements on the hemisphere and measured the quality of the reconstructed image. We modeled directivity, frequency range, size, and the arrangement of transducers in our simulations. In all the above configurations, transducers with a diameter of 12.7 mm, frequency of 3.5 MHz, 63% −6 dB fractional bandwidth, with an equidistant fashion from each other, center to center, were used. The reason that the transducers used in the simulations were 3.5 MHz (lower than what was used in the experiment, i.e., 5 MHz and 63% bandwidth) was because we did not have sufficient computational power to run the simulations with spatial sampling below the Nyquist limit. The difference between the simulation and experimental setup will not affect the justification of using the maximum number of transducer on the hemisphere, because the presented results show the relative difference in different numbers of transducers; using a higher frequency transducer will only improve the resolution of the resultant image. To model the element size of the transducers, each point transducer was expanded to a neighborhood of 4000 omnidirectional point transducers to mimic the size of 12.7 mm diameter. Transducer directivity was realized using spatial averaging on the surface of each single element transducer consistent of the 4000 point transducers. To place the transducers equidistant, we used the algorithm presented in [49]. The algorithm for the equidistant placement of transducers on a surface of the hemisphere was derived by first plotting a 2D curve of n points where the coordinates are determined using: $\left(x,y\right)=\left(s\cdot a,\pi /2\times \mathrm{signum}\left(s\right)\times \left(1-\sqrt{1-\left|s\right|}\right)\right)$ where s can be ranged from $-1+1/\left(n-1\right)$ to $-1-1/\left(n-1\right)$ in steps of $\left(2-2/\left(n-1\right)\right)/\left(n-1\right)$, $a=0.1+1.2n$, and signum is a function that determines the sign of s. These 2D cartesian points are then converted into 3D Cartesian points using $\left(x^{\prime },y^{\prime },z^{\prime }\right)=\left\{\cos \left(x\right)\cos \left(y\right),\sin \left(x\right)\cos \left(y\right),\cos \left(y\right)\sin \left(y\right)\right\}$.

The arrangement of the highest number (50) of equidistant elements on the hemisphere is shown in Figure 1a.

For the generation of simulation data in a hemisphere with a diameter of 12 cm, a k-wave volume of 12 × 12 × 12 cm³ was created using 750 × 750 × 750 voxels with 160 µm width in every direction. Our simulation provides two spatial samples per wavelength and meets Nyquist’s limit. Since k-wave is a pseudospectral method, for a completely accurate simulation it requires at least 10 grid points per wavelength, which is beyond our computational power. We, therefore, inspected the contrast to noise ratio (CNR) of reconstructions of a simple triangle phantom with differing grid-sizes (i.e., 300 × 300 × 300, 400 × 400 × 400, 500 × 500 × 500, 750 × 750 × 750 voxels). CNR was calculated based on the formula provided in [50], as the absolute value of the difference between the average object and background pixels is divided by the square root of the sum of the variances of the object and background pixels. Areas of greater signal with visually well-defined boundaries of the imaging target were chosen to be object pixels. Background area was chosen near the boundary of the reconstruction where the object was known to not exist. From least-dense to most-dense grid-size, the CNR values were found to be: 1.07, 1.38, 1.44, and 1.45. From the result of negligible difference in CNR with further increasing grid-size and simplicity in our simulated environment (no scatterers) as well as maximal use of our computational power, our use of 750 × 750 × 750 voxels grid-size is justified. For the simulation study, an imaging target consisting of three thin orthogonal rectangular legs with lengths of 1cm and widths of 2 mm each, was defined in k-Wave simulation [47] as shown in Figure 1a. The simulated object is similar to the imaging target that we used in the experiments. The initial pressure was set to 10 Pa and sampling frequency was set to 50 MHz, same as what was used in the experiment.

Image reconstruction is a process through which transducer signals are transformed into an image. Our reconstruction volume for both simulation and experimental data was 3 × 3 × 3 cm³ and made up of 400 × 400 × 400 voxels with 75 µm width in every direction. Different beamformers treat transducer signals differently to generate an image [51]; some of them are simple to implement such as delay-and-sum (DAS), some are very accurate and exact such as universal back projection (UBP), and some are more robust to noise and artifacts such as delay-multiply-and-sum (DMAS). DAS is the most commonly used beamforming method. It is popular due to its simplicity and the fact that it can be used to reconstruct both US and PA images which is ideal for commercial products [52,53,54]. DAS is a blind beamforming method which treats all the detected pressure waves the same way. Mathematically, DAS is the simplified version of the UBP algorithm at which only the first term is used to backproject (see equation 20 in [55]); UBP algorithm is the exact time-domain back-projection method, as illustrated in [55]. The different frequencies accompanying UBP in Figure 1c are the cut-off frequencies of the filter used in UBP. This filter is utilized to mitigate the effects of the amplified noise by the ramp filter existing in UBP formula. In DAS, sidelobes of the radiation pattern of the US transducers affect the image quality. To suppress the effects of sidelobes and off-axis signals, coherence factor (CF) (defined as the ratio of the coherent and incoherent summation of the recorded signals [56]) can be used along with DAS, we call that DAS-CF. CF is an index to show focusing quality and works based on the ratio of mainlobe energy to the total energy.

While DAS, UBP and filtered back projection follow a general delay and sum procedure, delay-multiply-and-sum (DMAS) is another beamforming method outperforming DAS in terms of resolution and contrast that use a correlation process between the received signals [56]. It should be noticed that all these algorithms are based on an ideal point transducer element with an infinite angular view. Figure 1c demonstrates the performance of different reconstruction algorithms including DAS, DAS-CF, DMAS, and UBP with multiple cutoff frequencies in terms of CNR of the simulated imaging target consisting of three thin orthogonal rectangular legs with 50 transducer elements.

DMAS [56] algorithm was chosen in this study for image reconstruction because it provided images with the highest CNR. DMAS algorithm has several variations, we implemented the fast DMAS following the method described in [57]. DMAS is conceived from the traditional DAS at which the detected signals are delayed proportional to the distance between the imaging target and the position of the elements. Summation of the delayed signal intensities is then performed, and an image is reconstructed. The equation for DAS, to determine each voxel intensity, is as follows:

$$y_{DAS}(k)={\displaystyle \sum _{i=1}^M{x}_i}(k-{\Delta }_i)$$

(1)

where $y_{DAS}(k)$ is the output of the beamformer, $k$ is the time index, and $M$ is the number of elements. $x_i\left(k\right)$ and ${\Delta }_i$ are the detected signals and the corresponding time delay for detector $i$, respectively. The time delay for each transducer element is equal to the distance between the transducer and the focused imaging pixel (this imaging pixel is specified by time index k) and is represented as:

$${\Delta }_i=d_i/c$$

(2)

DMAS consists of time-shifting radio-frequency (RF) signals of each transducer element to align them in phase (delaying procedure), followed by multiplying and summing them. DMAS can be represented as:

$$\begin{array}{l}{y}_{DMAS}(k)={\displaystyle \sum _{i=1}^{M-1}{\displaystyle \sum _{j=i+1}^M{x}_{id}}}(k)x_{jd}(k)=\\ \underset{\mathrm{first}\mathrm{term}}{\underbrace{[x_{1d}(k)x_{2d}(k)+x_{1d}(k)x_{3d}(k)+\dots +x_{1d}(k)x_{Md}(k))]}}+\underset{\mathrm{second}\mathrm{term}}{\underbrace{[x_{2d}(k)x_{3d}(k)+x_{2d}(k)x_{4d}(k)+\dots +x_{2d}(k)x_{Md}(k)]}}+\\ \dots +\underset{(\mathrm{M}-2)\mathrm{th}\mathrm{term}}{\underbrace{[x_{(M-2)d}(k)x_{(M-1)d}(k)+x_{(M-2)d}(k)x_{Md}(k)]}}+\underset{(\mathrm{M}-1)\mathrm{th}\mathrm{term}}{\underbrace{[x_{(M-1)d}(k)x_{Md}(k)]}}\end{array}$$

(3)

where $x_{id}(k)$ is the delayed signal corresponding to the $i\mathrm{th}$ element of the array. The issue of the squared dimension of the output of Equation (3) is also solved using the method introduced in [58].

2.2. Experimental Setup

A Q-switched Nd: YAG laser (Changchun Industries Optoelectronics Tech. Co., Changchun, Jilin Province, China) with 10 Hz, 8.5 ns pulse width at 532 nm wavelength has been used for illumination purposes. For full-field illumination, 40, 5-mm diameter poly-methyl methacrylate (PMMA) (epef-10, Ever Heng Optical Co., Shenzhen, China) optical fibers were bundled into a custom-made aluminum fiber bundle holder and aligned to the laser beam using 1′′ round top-hat 20° diffuser (EDI-C20-MD, Thorlabs, Newton, NJ, USA). The overhead illumination was done using a single fiber with a 12-mm diameter; the bundle was removed and a single optical fiber was coupled to the laser.

The integral part of the system is the 3D-printed hemisphere cap to hold the distal end of the optical fibers and transducer-amplifier units. Based on the simulation results, we developed this hemispherical dome-shaped cap (Fusion 360, v2.0.5688-Year: 2019, San Rafael, CA, USA) with a total of 90 equidistant cylindrically hollow spaces (Figure 1b): 50 for the transducers and 40 for optical fibers. The experimental setup is shown in Figure 2a. The cap has a diameter of 12 cm and is made of polyactic acid (PLA). The cap has three holes at the base to hold it in upright position using three horizontal optical rods. These rods are connected to three vertical rods which are fixed to the optical table. Out of 90 hollow spaces, 40 equally spaced smaller holes around the hemisphere are dedicated to fit in the distal end of the 5 mm PMMA optical fibers (length: 1.75 m, NA: 0.70) for homogeneous illumination (see Figure 2a(ii)). These holes were left empty when the single fiber (length: 1 m, NA: 0.35) overhead illumination was used, where the fiber was held in place directly above the center of the hemisphere (see Figure 2a(iii)). The other 50 equidistantly located hollow spaces (diameter: 16 mm) throughout the cap are dedicated to hold 5 MHz single element ultrasound transducers (ISL-0504-GP, Technisonic Research Inc., Fairfield, CT, USA). These transducers are directional and would only provide decent imaging at the center due to the acoustic field overlap they create, as compared to that of focused transducers which is less homogenized but covers a larger field of view. The transducers were connected to low noise 24 dB amplifiers (ZFL-500LN, Mini Circuits RF/Microwave Components, Brooklyn, NY, USA) to amplify the raw signals. SMA cables were used to establish a connection between transducer-amplifier units and the data acquisition (DAQ) system (NI PXIe-1078, National Instruments, Austin, TX, USA). The NI system contains seven eight-channel 14-bit data acquisition cards (NI PXIe-5170R). All 50 transducer-amplifier units were sampled simultaneously with the sampling frequency of 50 MHz. Two-hundred frames of data were collected and averaged for SNR improvement at each transducer location with a total acquisition time of 20s, considering that the laser repetition rate was 10 Hz. The trigger management for the laser and DAQ unit synchronization was done in Labview. System specification details are provided in

Table 2. 3D HA-PACT System specification.

Parameter	Value
Hemisphere radius	6 cm
Transducer active aperture/number	12.7 mm/50
Optical fiber diameter/number—full-field illumination	5 mm/40
Optical fiber numerical aperture—full-field illumination	0.70
Optical fiber diameter/number—overhead illumination	12 mm/1
Optical fiber numerical aperture—overhead illumination	0.35

.

To prepare the phantom, the hemisphere was held in place with the imaging target held by a 90° rod holder (RA90, ThorLabs, Newton, NJ, USA) inside the space of the hemisphere. Two objects were constructed for imaging, (1) triangle: an equilateral triangular shape with 2-cm leg length (Figure 2b(i)), and (2) four leg: pyramid shape (without base) of 1.5-cm leg length using 2-mm diameter copper wire soldered together (Figure 2b(ii)). All copper rods were covered in black tape as an absorbing coating. The inner surface of the hemisphere was covered with transparent (both optically and acoustically) saran wrap (uBoxes, Miramar, FL, USA). We used a multibeam swept-source optical coherence tomography system (SS-OCT) (VivoSight, Michelson Diagnostic™ Inc., Maidstone, UK) and measured the thickness of the saran wrap to be 40 µm [59]. The result is shown in Figure 2a(i). Next, 8% gelatin solution in distilled water was poured into the hemisphere and cured at room temperature. The transducer-sensing surfaces were coated in ultrasound gel before placing into assigned positions in contact with the saran wrap.

Finally, based on the transducer diameter and central frequency, the target object in the hemisphere is within the near-field. Imaging within the near-field results in rapid oscillation in received pressure with distance away from the transducer. We have accounted for this oscillation by recording the received pressure profile by imaging a thin lead (with the diameter of 500 microns) from a distance of 4.5 to 7.5 cm in 100-µm increments and recording the PA signal. We used a digital moving stage (Applied Scientific Instrumentation, Eugene, OR, USA) which travelled at 1 mm/s during acquisition. A complementary weighting profile was then constructed from the near-field profile for compensation. This experiment was repeated 10 times to reach sufficient precision. This profile was interpolated and applied to the raw data before reconstruction. It is possible to remove this step by utilizing highly focused transducers that have very short focal length.

3. Results

3.1. Simulation Results

The simulated object was reconstructed and the reconstructed images using a different number of equidistant transducers, i.e., 14, 26, 38, and 50, are shown in Figure 1d. The quality of the reconstructed images as a function of number of elements was quantitatively assessed by evaluating the CNR of the images; the results are given in Figure 1e. Additionally, we present a square root function fit to demonstrate the change in CNR with an increasing number of elements. As defined in Section 2.1, CNR was calculated using the square root of the sum of the variances of the object pixels and background pixels. Therefore, it was assumed that CNR was proportional to the square root of the number of elements in the system. Results of Figure 1e promoted the use of the maximum number of elements that can fit on the surface of the 12-cm diameter hemisphere.

For light illumination, we evaluated two configurations: (i) single-fiber overhead illumination, in which a single large optical fiber was used on top of the sample and illuminated the entire sample, and (ii) full-field illumination using a bundled optical fiber, in which 40 optical fibers were uniformly distributed on the residual space between the transducer holes in the designed hemisphere using the same algorithm described in [49]. This method was designed based on the fact that the larger the area of illumination, the more light we will be able to couple to the tissue [40]. The 3D projection model of the hemisphere showing the transducer and optical fiber locations is shown in Figure 1b. The fluence map of the two configurations are visualized using 3DSlicer software (Brigham and Women’s Hospital, Harvard Medical School, USA) and presented in Figure 1f. The simulations were performed using Zemax. The parameters considered in the Zemax simulations are as follows: 1,000,000 rays, 10 mJ input energy for both configurations; NA for the full field and overhead illuminations were 0.70 and 0.35 respectively; and diameters for the full-field fibers and overhead fiber were 5 mm and 12 mm respectively. As shown in Figure 1f(i,ii), the full-field illumination configuration provides a more homogenous illumination with higher optical energy compared to single-fiber illumination. We quantitatively compared these fluence maps within a sphere of 2-cm diameter right above the center of the cap, where the imaging target is located at. The average fluence and standard deviation for the full-field illumination within this sphere were 0.038 and 0.003, while those were 0.078 and 0.01 for overhead illumination, which suggests a less homogeneous illumination by the overhead configuration and consequently lower quality images.

3.2. Experimental Results

Experimental setup of the data acquisition is shown in Figure 2a, demonstrating our two methods of illumination: (ii) full-field illumination and (iii) single-fiber overhead illumination. The fast DMAS reconstruction algorithm was used for image reconstruction with the near-field compensation curve applied prior (please see the excplanation in Section 2.2). We demonstrated in Figure 2b that full-field illumination (v,vi) has significantly improved the reconstruction compared to overhead illumination (Figure 2b(vii,viii). The simulated target object has also been reconstructed for comparison purposes as shown in Figure 2b(iii,iv). The main reason was that a more homogeneous fluence map was generated by full-field illumination compared to overhead illumination (see the explanation provided in Section 3.1). CNR of the reconstructed images using both methods on two imaging targets was quantified; they are given in Figure 2c. As expected, the full-field illumination yielded improved image quality due to the even deposition of optical energy around the imaging target, thus producing a uniform initial induced pressure. The absolute CNRs were calculated for experimental images. For the triangle-shaped phantom, CNR of the image generated by the full-field illumination is 3.2 times greater than that generated by the overhead illumination. For the pyramid-shaped phantom, the improvement is 1.6 times. In terms of spatial resolution, the developed system was able to reconstruct objects to be as close as 0.9 mm to each other.

3.3. System Calibration

Manual placement of the transducers in the hemisphere resulted in small deviations in the distance between each sensing surface and the target. Post-processing was performed to shift the signals of each channel in time domain. Using an optimization algorithm and an appropriate cost function, an automatic calibration can be performed. In this study, we used the continuous sequential (CS) optimization algorithm [60]. According to the Nyquist theorem, the step size should be smaller than half the wavelength. Considering the central frequency of the transducers, i.e., 5 MHz, its bandwidth and the speed of sound of 1500 m/s, an optimization step size of 100 µm is adequate which is equivalent to three samples with a 50 MS/s DAQ sampling frequency. A flowchart of the continuous sequential algorithm is shown in Figure 3a. In this optimization algorithm (see Figure 3a), each channel was shifted with an increment of 100 µm to a maximum displacement of 1.5 mm forward and 1.5 mm backward, creating 31 possible sets of shifted data per channel. These data were then reconstructed using the DMAS algorithm. Due to the very large search space, i.e., 31⁵⁰ combinations and lack of a reference image, an automatic optimization process will take a very long time. Therefore, the optimality of the channel data shift was determined by the user through visually assessing reduction in background signal and alignment of imaged object features. The best reconstructed image out of 31 for each of the 50 channels was chosen sequentially and the corresponding signal shift was stored. The reconstructed images before and after optimization are shown in Figure 3b(i,ii), and the quantification of the CNR for each volume is shown in Figure 3c. As seen in Figure 3c, CNR was improved by 30% after optimization.

4. Discussion

Photoacoustic imaging has shown great promise in preclinical studies [11,15,24,61,62,63,64,65,66]. It has also shown great potential for clinical translatability. Such translation, especially for neonatal brain and breast imaging, requires an imaging system with a large detection view, fast acquisition and volumetric imaging capability [7,20,21,22,23,24,25]. Towards clinical translation, we have tested the feasibility of a stationary 3D HA-PACT imaging system with a full-field illumination scheme. Through quantitative analysis of simulation results, we have demonstrated that using a greater number of transducers in the HA-PACT system results in higher quality images; this was expected as the greater number of view angles will provide the image reconstruction algorithm with improved spatial averaging, which aids in artifact removal [56]. Therefore, based on the diameter of transducers and the thickness of the designed hemispherical structure, 50 single element transducers are the maximum amount that can be placed on a 12-cm diameter hemispherical dome. Although the physical increase in the number of transducers is limited by the space on the hemisphere, we ran a series of simulations (Figure 4a,b) to show that using a greater number of transducers with smaller sizes improves the CNR as shown in Figure 4c. Such configuration can be implemented using either smaller transducers or the same size transducers with a rotating scheme.

We tested two illumination configurations, full-field illumination and single-fiber overhead illumination. Our experimental results (shown in Figure 2) demonstrated that the proposed homogeneous illumination outperforms the traditional overhead illumination. This was mainly because the overhead illumination does not provide uniform deposition of optical energy, especially on the underside of the objects where they are directly in the field of view of the transducer elements. This inadequate illumination induced a deteriorated signal at the propagation angle to the elements and, therefore, worsened the reconstructed image quality. According to the literature [40], the larger the area of illumination, the more light we can couple to the tissue. This justifies the use of the entire area of the hemisphere and placement of fibers in an equidistant fashion on it.

Table 3. Estimated cost of the proposed system.

Parameter	Specifications	Cost ($)
Transducers	50 single elements/512 elements	~15K/~65K
Amplifiers	Low-noise 24 dB/customized	~7.5K/~25K
DAQ	56 channel/512 channel (or 64 channel with MUX)	~20K/~50K

provides an estimated cost to develop the proposed system. The overall cost of the proposed hemispherical system is approximately $42K. Other PACT systems such as the one developed by Upputuri et al. cost $15K using the rotation of one single element. Another system based on ring array presented in [67] that costs over $100K.

Our preliminary design was a proof of concept to investigate the feasibility of 3D imaging using single element transducers, in which we did not use the most optimum transducers in terms of size, type and sensitivity. The laser repetition rate was not optimal either. Our end goal is to develop a real-time 3D imaging system that is capable of hemodynamic imaging. One way to improve the proposed design is to rotate the hemisphere; such design will maintain inexpensive hardware, but virtually increase the number of view angles and hence dramatically improve the quality of the reconstructed image [68,69]. Utilizing a laser with a rep-rate of 200 Hz with 20 rotation steps and 10 frame acquisition per step for averaging will allow for 1 Hz frame rate full-volume imaging which will satisfy the hemodynamic imaging requirement [33].

We demonstrated (see Figure 1c) that the implementation of a more sophisticated reconstruction algorithm improved the quality of the reconstructed image [50]. This is because reconstruction methods such as DAS, UBP, and filtered back projection follow a general delay-and-sum procedure. These algorithms are based on an ideal point transducer element with an infinite angular view. Whereas more sophisticated algorithms such as the ones presented in [70,71,72,73,74] consider a more realistic acoustic field for the transducers in the image reconstruction process. Further, it is possible to utilize a speed of sound map in the reconstruction of the imaging targets, especially in biological tissue, and use iterative image reconstruction algorithms [69] for further improvement of the reconstructed images.

With regards to optimization, utilizing optimization methods such as the simulated annealing algorithm or Genetic algorithm along with an effective cost function will allow an automated and more effective calibration [60,75,76,77,78,79,80,81].

5. Conclusions

In this paper, we have developed a stationary 3D Hemispherical Array Photoacoustic Imaging System using 50 single element transducers. The performance of the system on several phantoms, using two light illumination configurations, i.e., overhead and full-field illumination, with different image reconstruction algorithms, was evaluated. CNR of the images were used for quantitative analysis of the results. The DMAS reconstruction algorithm with full-field illumination produced the images with the best quality and quantitative results. In the future, utilizing techniques to increase the number of view angles for ultrasonic detection, a more sophisticated image reconstruction algorithm considering the speed of sound map of the biological tissue, and a fluence compensation algorithm will be explored.