2.1.

Proposed HD-UNet Architecture

Since its advent, U-Net-based CNN has been widely used in complex imaging-related tasks and it comprises contraction and expansion layers with skip connections resembling a symmetrical U-shape.⁴⁷ However, for improved accuracy and performance in U-Net, extensional techniques are needed.⁴⁸ A modified version of U-Net, called fully dense U-Net (FD-U-Net), was first proposed for artifact removal and was then attuned for various PAI applications.^49,50 The FD-U-Net incorporates dense blocks in both the contracting and expansive layers to enable the learning of additional feature maps from the knowledge gained by previous layers through concatenation. Furthermore, the dense blocks increase the network’s depth without incrementing the number of layers. An enhanced version of FD-U-Net termed dense dilated U-Net (DD-U-Net) was then proposed for correcting the artifacts in three-dimensional (3D) PAT systems.⁵¹ The dense dilated blocks employed in the DD-U-Net uses atrous convolutions along with standard convolutions in the dense blocks to increase the receptive field to extract additional information. Furthermore, the incorporation of atrous convolutions in the dense blocks allows the CNN to learn multiscale features in an exponential means.⁵² However, a significant limitation with the DD-U-Net is the memory constraint due to a large number of training parameters and gridding artifacts if dilated convolutions of large receptive fields are employed.⁵¹ Thus, for improving the frame rate of the PAT systems, we developed a DL architecture HD-UNet by leveraging the benefits of both dilated convolution and standard convolution.

The proposed network incorporates dilated dense blocks in the encoding path followed by the standard dense blocks in the decoding path along with a residual block as the bridge. The schematic of the HD-UNet architecture is shown in Fig. 1. Depending on the layer level ($l$), the dense block employed in the encoding path intends to learn $f_l$ feature maps from the input feature map $f_i$ by iteratively learning $k_l$ features maps at each step. The standard dense blocks employed in the expanding path learns $f_l=2^{l-1}\times f_i$ at the growth factor of $k_l=2^{l-1}\times 8$ using $3\times 3$ convolutions of dilation rate 1. On the encoding path, the dilated convolutions implemented also learn $k_l$ features and it can be represented as $k_l={[\frac{k_l}{2}]}_s+{[\frac{k_l}{2}]}_d$, where ${[\frac{k_l}{2}]}_d$ refers to features from convolutions with dilation rates 1 (standard convolution) and ${\frac{k_l}{2}}_s$ are from convolutions with dilation rate 2 (dilated convolution). The dilation rate is limited to 2 for reducing the gridding artifacts. The down sampling operation in the encoding path is carried out by a $1\times 1$ convolution block followed by a $3\times 3$ convolution block of stride 2 and the up-sampling operation in the decoding path is performed by a transposed $3\times 3$ convolution block of stride 2. Skip connections were also implemented at each level to prevent the loss of any spatial information. Two $3\times 3$ convolution blocks were employed at the end of decoding to generate the resultant image. Each convolution block used in the model consists of batch normalization preceded by convolution and rectified linear unit (RELU) activation [$\mathrm{RELU}(x)=\max \{x,0\}$, where $x$ is the input to the neuron]. The proposed HD-UNet accepts an input image $X$ of size $128\times 128\text{pixels}$ and generates an output image $Y$ of size $128\times 128\text{pixels}$.

Fig. 1

Schematic of the proposed HD-UNet architecture incorporating dense blocks with dilated convolution in contracting path and dense blocks with standard convolutions in expanding path. L1, L2, L3, and L4 refers to the different layer levels. $X$ and $Y$ are the input and output image of size $128\times 128\text{pixels}$.

2.2.

Network Optimization and Implementation

The HD-UNet was implemented in Python 3.9 using the Tensorflow (V2.7) DL library.⁵³ The optimization of the network was performed on an Nvidia Tesla V100-32 GB GPU using the nodes of the Gekko cluster, High-Performance Computing Centre, Nanyang Technological University, Singapore. Adam optimizer with a call-back monitor reducing the learning rate by a factor of 0.5 on instances of no improvement in the monitored metrics was used. The initial learning rate was set to be 0.001. The loss function employed in the model is a composite of two-loss functions with weights $k_1$ and $k_2$ and the equation is expressed as

The $L_{\mathrm{FMAE}}$ is the Fourier mean absolute error loss (FMAE) and is applied to enforce the pixel-wise similarity between the ground truth $Y_g$ and predicted image $Y_p$, which is given by

The weights $k_1$ and $k_2$ used for optimizing the network are 1 and 0.001, respectively. The weights were chosen in such a way that the pixel-wise MAE serves as the primary loss. As the FMAE can contribute to instability in training, a smaller weighing factor was chosen. In total, the model was trained for 100 epochs with a batch size of two, and its performance was evaluated after the training.

2.3.

Simulated Photoacoustic Datasets for Training

DL is a data-based optimization approach, and its performance relies on the quality of the training data. In general, the training dataset used for optimizing the model comprises an input image and ground truth image. Although it is viable to generate a large amount of input data experimentally, the ground truth experimental data are sometimes difficult to obtain. Thus, for optimizing the HD-UNet, we used the simulated dataset with $k$-Wave MATLAB toolbox⁴⁶ for training the model. For generating the simulated data to improve the frame rate in the 1-UST-PAT system, three numerical phantoms such as: five-point targets (in which the position, orientation, and source strength of the point sources were varied randomly), triangles (in which the position, orientation, and size of the triangles were varied randomly), and vessel shapes-mimicking the cerebral venous sinuses of the rodent brain (in which the orientation, magnitude, and position were varied randomly) were used [Figs. 2(a)–2(c)]. A computational grid of $82\times 82\mathrm{mm}$ ($0.2\mathrm{mm}/\text{pixel}$) and a perfectly matched bounding layer were used for the simulation [Fig. 2(d)]. The imaging region was constrained to 40 mm. The SNR was maintained at 40 dB and 40 ns step size with 1500-time steps was used. The medium chosen was acoustically homogeneous and the speed of sound used is $1500\mathrm{m}/\mathrm{s}$. For generating the input PA data, the number of detector positions (sensor points) was randomly varied between 10 and 50 at steps of 10 ($\sim 1$ to 5 s scanning time), a large aperture unfocused UST (13 mm active area) of central frequency 2.25 MHz with 70% nominal bandwidth was used as the detector. For the ground truth data generation, 4800 detector positions with an ideal point detector of central frequency 2.25 MHz and 70% nominal bandwidth were considered.

Fig. 2

(a) Five-point sources numerical phantom. (b) Triangular numerical phantom. (c) Numerical blood vessel phantom. (d) Schematic of the $k$-wave simulation geometry of 1-UST-PAT system, (e) simulation geometry of the 8-UST-PAT system. (f) Five-point source phantom made up of pencil leads. (g) Triangular phantom made up of horsehair. (h) Photograph of the rat brain area used for in vivo imaging. (i) Schematic of the 1-UST-PAT system employed for phantom imaging. (j) Photograph of the 8-UST-PLD-PAT system used for in vivo brain imaging. AMP, amplifier; SM, stepper motor; DAQ, data acquisition card; PC, personal computer; UST, ultrasound transducer; AM, anesthesia machine; AH, animal holder; GG, ground glass; P1, P2, and P3 are uncoated prisms; S is the imaging sample.

For improving the frame rate of the 8-UST-PAT system, vessel shapes resembling the rodent cerebral sinuses were used. A computational grid consisting of $82\times 82\mathrm{mm}$ ($0.1\mathrm{mm}/\text{pixel}$) and a perfectly matched bounding layer were considered. The schematic of the computational grid is shown in Fig. 2(e). The SNR was maintained at (10 to 20 dB). 1500-time steps with a step size of 40 ns were used for recording the A-lines. The medium used is acoustically homogeneous and the speed of sound was maintained at $1500\mathrm{m}/\mathrm{s}$. For the input data generation, eight large aperture unfocused UST (5 MHz central frequency with 70% nominal bandwidth) with 240 detector positions (30 detector locations per UST) were used. For the generation of ground truth data, an ideal point detector (5 MHz central frequency and 70% nominal bandwidth) with 1600 detector positions was considered. In both cases, conventional delay-and-sum beamformer was employed to reconstruct the PA data into cross-sectional PAT images of size $128\times 128\text{pixel}$. The reconstructed PAT images were then normalized by rescaling it in the new range of 0 to 1 without the loss of bipolar information using the equation,

2.4.

Experimental Phantom Data

Experimental phantom imaging was performed for evaluating the performance of the optimized HD-UNet. For obtaining the experimental phantom data, two types of phantoms namely, five-point targets (made up of pencil leads) [Fig. 2(f)] and triangular phantom (made up of horsehair) [Fig. 2(g)] were used. 1-UST-PAT system was employed for phantom PAT imaging.⁵⁴ The schematic of the 1-UST-PAT system is shown in Fig. 2(i). A Q-switched Nd:YAG laser delivering 532 nm laser pulses with 10 pulses per second at a pulse width of 5 ns was employed as the excitation source. The emergent laser beam was homogenized using an optical diffuser and the laser energy density was maintained at $\sim 6\mathrm{mJ}/{\mathrm{cm}}^2$. An unfocused UST of 2.25 MHz (Olympus-NDT, V306-SU) central frequency (70% nominal bandwidth) was used to acquire the PA signals. An ultrasound pulse-receiver (Olympus-NDT, 5072PR) with a gain of 48 dB is used to amplify the PA signals. The amplified PA signals were then stored inside the desktop computer using a data acquisition card (DAQ) [GaGe, compuscope 4227]. Conventional delay-and-sum beamformer was used to reconstruct the cross-sectional PAT images from the PA data.

2.5.

In Vivo Experimental Data

The performance of the proposed HD-UNet was also evaluated on the in vivo PAT imaging. Sprague Dawley rats ($\sim 95\mathrm{gm}$) obtained from InVivos Pte. Ltd., Singapore were utilized for imaging [Fig. 2(h)]. The rats were anesthetized by the intraperitoneal administration of ketamine ($100\mathrm{mg}/\mathrm{mL}$) and xylazine ($20\mathrm{mg}/\mathrm{mL}$) mixture. The hair on the rat head was then removed using depilatory cream, and the ocular gel was applied before imaging. A layer of ultrasound was applied to the scalp and a constant supply of anesthesia ($1.0\mathrm{L}/\min$ oxygen and 0.75% isoflurane) was maintained during imaging. All the animal experiments were performed as per the guidelines of the Institutional Animal Care and Use Committee, Nanyang Technological University, Singapore (Protocol No.: A0331). The in vivo PAT imaging was performed using the 8-UST-PLD-PAT system.²³ The image of the 8-UST-PLD-PAT system is shown in Fig. 2(j). A PLD capable of delivering $\sim 816\mathrm{nm}$ wavelength at 2000 pulses per second with a pulse width of $\sim 107\mathrm{ns}$ and 3.4 mJ per pulse energy was used as the excitation source. An optical diffuser was employed to homogenize the emergent rectangular beam from the PLD, and the laser energy density was maintained at $\sim 0.17\mathrm{mJ}/{\mathrm{cm}}^2$, below the safety limits of American National Standards Institute (ANSI).⁵⁵ Eight unfocused USTs (5 MHz central frequency with 70% nominal bandwidth) fitted with 45 deg acoustic reflectors (Olympus-NDT, F102) were employed to acquire the PA data. The PA signals were then amplified using a 48 dB low signal noise amplifier (Mini-circuits, ZFL-500LNBNC, two of them in series each with 24 dB gain) and stored inside a computer (IntelXeon, 3.7 GHz 64-bit processor, 16 GB RAM) using a DAQ card (Spectrum, M2i.4932-Exp). Conventional delay-and-sum beamformer was used to reconstruct the cross-sectional PAT brain images.

For in vivo imaging, the maximum permissible exposure (MPE) is limited by the ANSI laser safety standards.⁵⁵ For the wavelength in the range of 700 to 1050 nm, the maximum per pulse energy density on the skin surface should not exceed $20\times {10}^{2(\lambda -700)/1000}\mathrm{mJ}/{\mathrm{cm}}^2$. For 816 nm wavelength, the MPE is $\sim 34.12\mathrm{mJ}/{\mathrm{cm}}^2$. For the illumination period of 1.5 s ($t=1.5\mathrm{s}$), the MPE safety limit is $1.1\times {10}^{2(\lambda -700)/1000}\times t^{0.25}\mathrm{J}/{\mathrm{cm}}^2$ ($=2.07\mathrm{J}/{\mathrm{cm}}^2$). Thus, for the scan time of 1.5 s, the MPE per pulse is $0.69\mathrm{mJ}/{\mathrm{cm}}^2$. Similarly, for a scan time of 0.3 s, the MPE safety limit is $1.39\mathrm{J}/{\mathrm{cm}}^2$, and per pulse, it is $2.31\mathrm{mJ}/{\mathrm{cm}}^2$. As the per pulse energy was maintained at $\sim 0.17\mathrm{mJ}/{\mathrm{cm}}^2$ during the in vivo imaging (the PLD used in our study produces a max energy of $\sim 3.4\mathrm{mJ}$ per pulse, and it was illuminating $\sim 20{\mathrm{cm}}^2$ area), the per pulse energy does not exceed the ANSI safety limit for the scan time of 0.3 and 1.5 s.

3.1.

Performance Comparison

After the optimization, $k$-fold cross-validation ($k=10$) was implemented to evaluate the performance of the proposed HD-UNet, and it was compared with the performances of other DL architectures such as the FD-UNet, 2D-DD-UNet (an adapted version of 3D-DD-UNet⁵¹), and U-Net, using a variety of loss metrics such as Pearson correlation coefficient (PCC), structural similarity index measure (SSIM), peak signal-to-noise ratio (PSNR), and MAE. For the performance evaluation, the original dataset was randomly split 10 times in to into training, validation, and testing. For each dataset, the FD-Unet, 2D-DD-UNet, and U-Net were optimized for 100 epochs and its performance was evaluated on the testing dataset. On evaluation, the HD-UNet exhibited superior performance over the FD-UNet, 2D-DD-UNet, and U-Net over all the metrics (Table 1) and it signifies the generalizability of the developed HD-UNet.⁵⁶

Table 1

k-Fold cross-validation (k=10) to compare the performance of the HD-UNet (mean ± standard deviation). The best values are shown in bold.

3.2.

Performance of HD-UNet on Simulated Phantoms

The reconstructed PAT images of three numerical phantoms (nine-point target phantom, triangular phantom, and vessel phantom) are shown in Fig. 3. Figure 3(a) shows the PAT image of the nine-point target phantom simulated using 1-UST-PAT configuration for a scan time of 5 s ($\sim 50$ A-lines). Figure 3(b) depicts the PAT image of the nine-point target phantom reconstructed using the HD-UNet and Fig. 3(c) shows the ground truth image [1-UST-PAT configuration, 8 min scan time ($\sim 4800$ A-lines)]. From Fig. 3(a) it can be noted that point targets are not clearly visible and were also marred by the presence of artifacts arising from the sparse data acquisition at higher scanning speeds. When the HD-UNet has been applied, the artifacts were corrected and the point targets were very well reconstructed [Fig. 3(b)]. Furthermore, the improvement in the tangential resolution over the points can also be noted especially at the farthest point (marked with small yellow arrows) and is very close to the ground truth image [Fig. 3(c)]. As the HD-UNet was not trained on the nine-point targets, the ability of the network to improve the quality nine-point target PAT image signifies its ability on unknown phantom data. Figures 3(d)–3(f) show the PAT images of the triangular phantom obtained using 1-UST-PAT geometry with a scan time of 5 s ($\sim 50$ A-lines), with the HD-UNet, and the expected ground truth. The potential of the HD-UNet to preserve the target shape along with improvement in the artifacts (indicated with small yellow arrows) can be observed by comparing Figs. 3(d) and 3(e). The vessel phantom simulated using 8-UST-PAT configuration for a scan time of 0.3 s (30 A-lines per transducer) is shown in Fig. 3(g). Figure 3(h) depicts the PAT image of the vessel phantom reconstructed using the HD-UNet, and Fig. 3(i) shows the ground truth image [8-UST-PAT configuration, 1.5 s scan time ($\sim 1200$ A-lines: 150 A-lines per transducer)]. From the comparison of Figs. 3(h) and 3(i), it can be noted that the HD-UNet can produce the image very close to the ground truth even when the scan time was five times less than the ground truth. The improvement noticed over the cerebral venous and veins can be envisaged through visual comparison of the areas indicated by small yellow arrows in Figs. 3(g)–3(i). The higher PCC values of the HD-UNet PAT reconstructed images signify that it plays a very good role by preserving the shape of the target along with improvement in image quality.

Fig. 3

(a)–(c) PAT images of nine-point target numerical phantom (1-UST-PAT configuration): (a) simulated for a scan time of 5 s, (b) reconstructed with HD-UNet, and (c) simulated for a scan time of 8 min (ideal image, or the ground truth). (d)–(f) PAT images of triangular phantom (1-UST-PAT configuration): (d) simulated for a scan time of 5 s, (e) reconstructed with HD-UNet, and (f) simulated for a scan time of 8 min (ideal image, or the ground truth). (g)–(i) PAT images of numerical vessel phantom using (8-UST-PAT configuration): (g) simulated for a scan time of 0.3 s, (h) reconstructed with HD-UNet, and (i) simulated for a scan time of 1.5 s (ideal image, or the ground truth).

3.3.

Performance of HD-UNet on Experimental Phantom Images

Experimental phantom imaging was performed on the 1-UST-PAT system to evaluate the performance of the HD-UNet at higher imaging speeds. As discussed before, two types of phantoms were utilized for the imaging. Figure 4(a) depicts the five-point target phantom PAT image obtained in a scan time of 5 s. The HD-UNet reconstructed image of the point target phantom is shown in Fig. 4(b). Figures 4(c)–4(e) show the reconstructed PAT image using the FD-Unet, 2D-DD-UNet, and U-Net. Figure 4(f) depicts the PAT image of the point target phantom obtained for a scan time of 8 min. From Fig. 4(f), it can be noted that even though the scan time is high (8 min) the shape of the point targets is not well preserved when its distance from the scanning center increases. When the HD-UNet was employed the shape of the point targets are well preserved along with the removal of artifacts [Fig. 4(b)]. Furthermore, the ability of the HD-UNet to preserve the target shape along with improvement in the quality of the image can be visualized by comparing the PAT image of triangular phantom obtained with a scan time of 5 s [Fig. 4(g)] and the reconstructed image using HD-UNet [Fig. 4(h)]. It can be noted that the edges of the triangular phantom (marked with yellow arrows), which are murkier at higher imaging speeds can be visualized when the HD-UNet was applied, and the resultant image quality is better than that of the PAT image of triangular phantom imaged at a scan time of 8 min [Fig. 4(i)].

Fig. 4

(a)–(f) Reconstructed PAT images of five-point target phantom (1-UST-PAT system): (a) obtained with a scan time of 5 s, (b) reconstructed with HD-UNet, (c) reconstructed with FD-UNet, (d) reconstructed with 2D-DD-UNet, (e) reconstructed with U-Net, and (f) obtained with a scan time of 8 min. (g)–(i) Reconstructed PAT images of triangular horsehair phantom (1-UST-PAT system): (g) obtained with a scan time of 5 s, (h) reconstructed with HD-UNet, and (i) obtained with a scan time of 8 min.

The 8-UST-PLD-PAT system described before is used for generating the in vivo brain images. Figure 5(a) shows the in vivo brain images obtained at a high frame rate of 0.3 s. It can be noted from the images that the cerebral venous sinuses such as the transverse sinuses (TS) are imperceptible (shown with small yellow arrows) due to low SNR at higher scanning speeds, which is a major hindrance for the analysis of various in vivo morphological studies such as intracranial hypotension, cerebral hemorrhage, etc. Furthermore, the presence of white artifacts arising from the limited bandwidth detection also limits the visual analysis of sagittal sinus (SS). When the HD-UNet was applied for reconstruction, the cerebral venous sinuses were perceptible (marked with small yellow arrows) along with the removal of artifacts [Fig. 5(b)]. It can also be noted that the HD-UNet also improves the tangential resolution without compromising the image quality in comparison with the conventionally reconstructed in vivo brain images obtained using a scan time of 1.5 s [Fig. 5(c)].

Fig. 5

Reconstructed in vivo PAT brain images (8-UST-PLD-PAT system): (a) obtained with a scan time of 0.3 s, (b) reconstructed with HD-UNet, and (c) obtained with a scan time of 1.5 s. TS, transverse sinus; SS, sagittal sinus.