2.1

3D fusion based on joint weighted evaluation optimization

The multi-focus 3D photoacoustic data was firstly decomposed into 8 wavelet coefficients via 3D stationary wavelet transform (3D-SWT)^[7]. The coefficients were then divided into several equal-sized blocks and STD is used to detect whether a block locates within the depth of field automatically^[8]. The blocks detected as focused were preserved in the fused coefficients thereafter and large-volumetric and high-resolution photoacoustic data can be achieved by inverse stationary wavelet transform. As shown in Figure 1, the multi-focus photoacoustic data P1 and P2 with size of H×W×L were obtained through the virtual OR-PAM^[9]. P₁ and P₂ were decomposed into 8 wavelet coefficients LLL, LLH, LHL, LHH, HHL, HHH, HLH and HLL (where H represents high-frequency filtering, L represents low-frequency filtering) using 3D-SWT, respectively. The wavelet coefficients were divided into several sub-blocks with equal size of h×w×l based on spatial location, and the fusion result varies with the block size. Oversized block may include both focused and defocused region and the fused data may contain out-of-focus area. Undersized block may not be able to provide sufficient data points for accurate focus detection. Therefore, DE algorithm was employed to iterate the size of the block region to obtain the optimized block size and corresponding fusion result adaptively^[10]. A weighted evaluation function L is proposed to evaluate the maximum amplitude projection (MAP) image of fused 3D data which combines three quantification parameters, namely information entropy (EN)^[11], average gradient (AVG)^[12], and structural similarity index (SSIM)^[13] as shown in formula (1)

Figure 1.

Flowchart of 3D photoacoustic information fusion based on joint weighted evaluation optimization. UT, ultrasonic transducer; P1, P2, multi-focus 3D photoacoustic data; SWT, stationary wavelet transform; ISWT, inverse stationary wavelet transform.

Where L_AVG, L_EN and L_SSIM are the quantitative evaluations of the MAP image of the fused 3D data. λ_AVG,λ_EN and λ_SSIM are the pre-defined weights for the three quantification parameters. Larger L of the fused data shows a better fusion quality. To find the optimal block size of sub-coefficients, the optimization problem was solved through DE as shown in expression (2).

Based on the optimized block size given by DE, the STD of corresponding sub-blocks of coefficients decomposed from multi-focus data are compared to detect the focus regions. The blocks with larger STD, which contain more details are kept as sub-blocks of fused coefficients. Blocks with smaller STD are discarded. The fused 3D data with extended DoF can be achieved by performing ISWT on the fused coefficients. Figure 1 takes the wavelet coefficient LLL as an example to illustrate the whole process of focus detection and block size optimization. The other seven wavelet coefficients are processed in the same way.

2.2

Multi-focus 3D data acquisition using virtual OR-PAM

A virtual Gaussian-beam photoacoustic microscopy based on k-Wave^[9] was adopted to obtain multi-focus 3D photoacoustic data in our experiment. The Gaussian beam in the virtual OR-PAM platform is generated by an objective lens with a numerical aperture of NA = 0.2, and the light wavelength was set as 532 nm. The number of three-dimensional grids is Nx ×Ny×Nz = 80×80×80 and the pixel size is dx = dy = 2 μm, dz = 3 μm. Water was used as the medium around the sample, and the speed of sound was set to 1500 m/s. A ultrasonic detector with a center frequency of 75 MHz and a bandwidth of 67% was used to collect photoacoustic signals. The sample (such as vessel) was set up as required and placed in the 3D grid. Diagonally positioned fiber with a diameter of about 2 μm was imaged and fused to measure the performance of proposed 3D fusion algorithm. A set of simulated vessels was used to further validate the effectiveness of proposed method. The photoacoustic signal was obtained by raster scanning and reconstruction. 3D photoacoustic microscopy of different focal positions can be acquired by changing the position of the optical focus. All simulations were performed with a 64-bit Windows 10, Intel(R) Core(TM) i7-10750H CPU @ 2.60GHz desktop running windows operating system.

3.1

3D fusion based on joint weighted evaluation optimization

The multi-focus vertically tilted fibers and fusion results are shown in Figure 2. Figures 2(a) and 2(b) are a set of multifocus vertically tilted fibers obtained by the virtual OR-PAM. Figure 2(a) is the MAP image of the tilted fiber with optical focus positioning at z = 30 in the 3D grid. Figure 2(a) is the MAP image of the tilted fiber with optical focus positioning at z = 40 (the optical focus is moved down by 30 μm) in the 3D grid and Figure 2(c) is the MAP image of fusion result of the tilted fibers through proposed method. The depth-coding MAP images of Figures 2(a)-2(c) are displayed in the lower right corner, repectively. The white arrows indicate the focal plane of Focus 1, and the yellow arrows indicate the focal plane of Focus 2. NPA is the normalized photoacoustic amplitude. It can be observed from Figures 2(a) and 2(b) that due to the limited DoF of OR-PAM, high-resolution and large-volumetric titled fiber cannot be captured during a single imaging. The upper and lower halves of the fused fiber preserve the in-focus portion of source data and the fused image has a larger DoF, as shown in Figure 2(c). Figures 2(d)-2(f) are the B images corresponding to the dotted lines in Figures 2(a)-2(c). Figure 2(g) shows full width at half maximum (FWHM) of the profile of the fiber before and after fusion. A smaller FWHM represents a better lateral resolution. The black circle corresponds to Focus 1, the blue star corresponds to Focus 2, and the red circle corresponds to Fusion. The lateral resolution of the focused part (0-70 μm along z direction of Focus 1 and 70-120 along z direction μm of Focus 2) is small, while lateral resolution of the defocused part (70-120 μm along z direction of Focus 1 and 0-70 μm along z direction of Focus 2) is relatively large, as shown in Figure 2(c). The DoF of OR-PAM is defined as the depth range in which the FWHM of the fiber becomes twice the nadir. The DoF of Focus 1 is estimated to be about 37.3 μm where the FWHM varies from 1.1 μm to 2.2 μm. The DoF of Focus 2 is estimated to be about 33.0 μm where the FWHM varies from 0.99 μm to 1.98 μm. Through our proposed method, the DoF of OR-PAM is extended to 67.2 μm where FWHM varies from 0.91 μm to 1.82 μm. The DoF assessments of fibers before and after fusion quantitively show that the DoF of OR-PAM can be doubled without the sacrifice of lateral resolution through our proposed method. To observe the 3D structure of fibers more intuitively, the above three sets of photoacoustic data were imported into the Amira software for 3D visualization as shown in Figures 2(h)-2(j).

Figure 2.

(a) and (b) are the MAP images projected from the tilted fibers obtained by the virtual OR-PAM when the focus positions are z = 30 and z = 40, respectively. (c) is the MAP image of the fusion data obtained through proposed method. Depth-coding MAP images of (a)-(c) are presented in the lower right corner, respectively. NPA is the normalized photoacoustic amplitude. (d)-(f) are the B images corresponding to the white dotted line in (a)-(c). (g) is the variation of lateral resolution along with the depth before and after fusion. (h)-(j) are the 3D visualizations rendered by Amira software corresponding to (a)-(c).

3.2

Large volumetric imaging of vessel

The multi-focus vessel data and fusion results are presented in Figure 3. Figures 3(a) and 3(b) are a set of MAP images projected from multi-focus simulation vessels obtained by the virtual OR-PAM. The positions of optical focus set in Figures 3(a) and 3(b) are z = 40 and z = 60 (The optical focus is moved down by 60 μm). NPA is the normalized photoacoustic amplitude. The lower right corner of Figures 3(a)-3(c) are the close-up MAP images of the areas indicated by the white dashed rectangles in the figures, respectively. From Figures 3(a) and 3(b), it can be seen that due to the limited DoF of OR-PAM, the whole vessel cannot be captured through a single imaging. Figure 3(c) is the MAP image of the fusion results obtained using the proposed method. To observe the detailed structure of the vessel, the above three sets of photoacoustic data were imported into the Amira for 3D visualization as shown in Figures 3(d)-3(f). Figures 3(g)-3(i) are the normalized amplitude distributions along the dashed lines in the partially enlarged views of Figures 3(a)-3(c) correspondingly. The FWHM of the signal intensity distribution within focused regions with a better lateral resolution shown in Figure 3(g) is measured to be 3.4 μm, while the FWHM of the signal intensity distribution positioning out of focused regions with a worse lateral resolution in Figure 3(g) is measured to be 9 μm. The corresponding signal intensity distribution of the fused vessel (Figure 3(i)) preserves the better lateral resolution of Figure 3(g), which further verifies the effectiveness of the proposed method in processing complex photoacoustic information.

Figure 3.

(a) and (b) are the MAP images projected from the simulated vessels obtained by the virtual OR-PAM when the focus positions are z = 40 and z = 60, respectively. (c) is the MAP image of the fusion data obtained through proposed method. The lower right corners of the (a)-(c) present the partially enlarged view of the white dotted rectangles in the figure. NPA is the normalized photoacoustic amplitude. (d)-(f) are the 3D visualizations rendered by Amira software corresponding to (a)- (c). (g)-(i) are the normalized intensity distributions of the white dotted lines in the enlarged images of (a)-(c).