1. Introduction

Photoacoustic (PA) imaging (PAI) has developed rapidly in the past two decades. In PA image acquisition, an object is excited by a light source, and ultrasound (US) waves emerge because of the PA effect. Then, one can recover the light absorber distribution by collecting the US waves with an US transducer and using a reconstruction algorithm. PAI has achieved promising results in microvascular imaging [1–3] and endoscopic imaging [4], to name a few. PA microscopy (PAM) is one implementation of PAI and can be classified into optical-resolution PAM (OR-PAM) and acoustic-resolution PAM (AR-PAM), depending on the focusing approach. Overall, OR-PAM enjoys high resolution but suffers from low imaging depth (∼1 mm) while AR-PAM can achieve deep imaging depth by taking advantage of diffused light and good acoustic penetration [5].

In AR-PAM, its lateral resolution depends on the center frequency and numerical aperture (NA) of a focused US transducer. High lateral resolution requires that the transducer has high center frequency and high NA, which in turn limit penetration and depth of focus (DOF) of AR-PAM. The former is because high-frequency US waves attenuate dramatically in biological tissue, and the latter is because a high-NA transducer has limited DOF. Besides, the ultimate resolution of AR-PAM is restricted by acoustic diffraction. Fortunately, previous works have demonstrated image enhancement in lateral resolution, DOF, and signal-to-noise ratio (SNR) by using algorithms, which circumvents the issues associated with the high-frequency and high-NA transducer.

Synthetic aperture focusing technique (SAFT) is often used in AR-PAM reconstruction for enhancing lateral resolution in the out-of-focus region. SAFT was first developed for US imaging to realize beamforming [6] and then was applied to PAI based on a point-like acoustic detector [7]. Later, Li et al. proposed the virtual-point-detector (VPD) concept to apply SAFT to AR-PAM using a focused transducer [8]. Afterward, SAFT was extended to two dimensions by Deng et al. [9]. Considering the line pattern of blood vessels in AR-PAM images, adaptive SAFT (A-SAFT) was further investigated [10]. To improve resolution in the focal region, spatial-dependent contributions were weighted using the spatial impulse response (SIR) in SAFT [11]. More recently, delay-multiply-and-sum-based SAFT [12], two-dimensional (2D) spherical SAFT (S-SAFT) [13], and directional SAFT (D-SAFT) [14] have been further studied to improve lateral resolution and SNR in the out-of-focus region. Lateral resolution in the out-of-focus region can be effectively improved by SAFT but is still ultimately limited to the in-focus lateral resolution. Therefore, deconvolution-based algorithms are used to further improve lateral resolution in PAI [13,15,16] and US imaging [17,18]. Specifically, in AR-PAM, the point spread function (PSF) is dominated by the finite US detection beam size of a transducer. Deconvolution of acquired PA images with PSF could recover the raw PA distribution. Previously, Cai et al. [13,16] used Richardson-Lucy (R-L) deconvolution algorithm to improve the AR-PAM lateral and axial resolutions in both the focal and out-of-focus regions. However, the nature of the line pattern of blood vessels was not considered, and thus, the recovered lines suffer the discontinuity issue. An alternative model-based (MB) deconvolution method assumes the acquired PA image is the result of convolution of a PA source (an absorber) and a dictionary matrix. Then, the PA source distribution could be recovered by using an optimization approach. It has been demonstrated that the MB deconvolution method can improve resolution in both PAI and US imaging [19–21]. Besides, prior knowledge of sparsity was considered as a regularization term in the model to further enhance the resolution [20,21]. As a matter of fact, the nature of the line pattern of blood vessels in AR-PAM images has only one-dimensional (1D) sparsity (in contrast to 2D sparsity in [20,21]). If it is not considered, results will have the discontinuity issue, as also encountered in the R-L deconvolution algorithm mentioned above.

To overcome the above challenges, we propose a novel two-step algorithm including a Fourier accumulation SAFT (FA-SAFT) and a directional MB (D-MB) deconvolution method to improve resolution, SNR, and image fidelity in both the focal and out-of-focus regions for AR-PAM [22]. A preliminary version of this work has been reported [22]. In this work, substantial progress was made on further refinement of the algorithms and their comprehensive evaluation on phantom and in vivo images with more complex structures. Our contributions can be summarized in two aspects. (i) FA-SAFT uses a sharpening filter to improve the high-frequency components of PA images in Fourier domain (k-space). Such a sharpening filter is adaptive to different patterns of the acquired PA images, leading to enhanced resolution and SNR. (ii) D-MB deconvolution takes 1D sparsity of the line pattern into account, which further enhances resolution and keeps image fidelity. The results by our algorithm are promising. By imaging tungsten wires, high resolution over large DOF and good image fidelity are experimentally achieved using our algorithm. Further, in vivo imaging is also demonstrated to show the superior capability of our algorithm. Compared with existing algorithms (D-SAFT, R-L deconvolution, etc.), our algorithm enables great imaging performance in resolution, SNR, and fidelity for both phantom and in vivo imaging demonstrations. The effect of SNR on the performance of our algorithm is also investigated and discussed.

2. Related techniques

We first describe some existing techniques related to our algorithm. The conventional SAFT algorithms are introduced in Sections 2.1. Then, adaptive-based SAFT algorithms including A-SAFT and D-SAFT are introduced in Sections 2.2. Finally, R-L deconvolution and MB deconvolution are introduced in Sections 2.3.

2.1 Conventional SAFT

In AR-PAM, linear scanning is conducted in one dimension or two dimensions. We first illustrate SAFT for the 1D case, as shown in Fig. 1. In SAFT, signal synthesis is performed based on a delay-and-sum algorithm along the scanning line, which forms a synthetic aperture, and can be expressed as:

 

Fig. 1. Schematic of SAFT for (a) a point detector and (b) a focused transducer (VPD) in one dimension. Scanning is conducted along the X direction.

Download Full Size | PDF

Due to the spatial response of the focused transducer, the previous study further proposed improved SAFT algorithm by incorporating SIR to weigh the spatial-dependent contributions of $\textrm{R}{\textrm{F}_i}(t )$ [11]. Equation (1) is thus modified as:

Furthermore, to suppress the out-of-phase summation in Eq. (1), coherence factor (CF) [7] is applied in SAFT and is defined as:

2.2 A-SAFT and D-SAFT

AR-PAM is mostly used for vascular imaging. Considering the directional characteristic of the line pattern of blood vessels, A-SAFT and D-SAFT were proposed by Deng et al. [10] and Jeon et al. [14], respectively, aiming at SAFT for vascular images. A-SAFT is based on 1D SAFT with the signal synthesis directions adaptively changed with the directions of vessel line branches. However, the performance of A-SAFT depends on the successful recognition of the vessel branches’ directions before applying A-SAFT. In other words, if one branch is not observed or recovered properly before A-SAFT is applied, it will not be recovered properly after A-SAFT.

To deal with this problem, D-SAFT was designed to recover vessel line branches in a directional way [14]. In D-SAFT, 1D SAFT along a series of directions is first applied, as described in the following equation:

2.3 R-L deconvolution and MB deconvolution

As mentioned previously, acoustic diffraction still limits lateral resolution to no better than in-focus lateral resolution. To further enhance resolution, one can use deconvolution algorithm. Overall, R-L deconvolution and MB deconvolution are two typical methods. The underlying assumption is that an AR-PAM imaging system is a linear and spatial shift-invariant system. Therefore, the acquired PA image, $\textrm{RF}({x,y} )$, by AR-PAM can be expressed as:

3. Methods

Overall, our algorithm has two parts, FA-SAFT and D-MB deconvolution, as shown in the green and blue parts in Fig. 2, respectively. They are introduced in Sections 3.1 and 3.2 in the following.

 

Fig. 2. Procedures of the proposed algorithm.

Download Full Size | PDF

3.1 FA-SAFT

Inspired from the work done by Delbracio et al. [28], we propose to improve resolution and SNR with FA-SAFT. FA-SAFT is modified from D-SAFT. Compared with D-SAFT, FA-SAFT exploits the Fourier bursts among 1D SAFTs along different directions to adaptively enhance the spectrum (k-space), and thus, FA-SAFT can produce improved results.

Note that using CF in SAFT (same for D-SAFT and FA-SAFT) can already increase SNR along with its spatial resolution. The key to optimizing image enhancement is to identify the optimal synthesis direction among 1D SAFTs along different directions. Although 1D SAFTs along different directions are utilized in both D-SAFT and FA-SAFT, FA-SAFT adopts an improved approach to combine the 1D SAFTs from all directions, leading to further improved resolution and SNR. Detailed comparison is elaborated in Section 5. FA-SAFT includes the following steps. First, 1D SAFT is applied to a series of directions using Eqs. (6) and (7). Secondly, the signal is transformed into k-space using Eq. (8). Thirdly, ${\textrm{K}_\theta }({u,v,\omega } )$ is accumulated according to the following equations:

Note that the proposed FA-SAFT approach can be further divided into 2 types, depending on whether the directional term ${\textrm{W}_{{\theta _n}}}({u,v,\omega } )$ is applied in Eq. (18) or not. The 2 types are termed FA-SAFT-dir0 and FA-SAFT-dir1, which are the cases of without applying the directional term and with applying it in Eq. (18), respectively. Then, the results after FA-SAFT will be used for deconvolution.

3.2 D-MB deconvolution

As the line pattern only has 1D sparsity, it stands to reason to apply the 1D deconvolution. Therefore, it is crucial to recognize the branches in different directions so that the 1D deconvolution along different directions can be applied accordingly. We use the hamming window methods as before (Eq. (10)) to decompose the branches into two directions that are perpendicular to each other and apply the 1D MB deconvolution with sparsity (Eq. (17)) for the branches. The algorithm is thus termed D-MB deconvolution. The detailed procedure is described in the following equation:

Note that ${\theta _n}$ indicates the deconvolution direction, which lies in the lateral plane and means the direction along which the 1D deconvolution is conducted. The final processed result after D-MB deconvolution is obtained by:

FA-SAFT can be considered to be independent of D-MB deconvolution. That is, FA-SAFT and D-MB deconvolution (with well-estimated PSF) can be applied alone. Therefore, $\mathrm{N^{\prime}}$ in (18) is independent of $\textrm{N}^{\prime\prime} $ in (24), and $\mathrm{N^{\prime}}$ can be different from $\textrm{N}^{\prime\prime} $. Note that $\textrm{N}^{\prime\prime} $ = 2 is used because we found that too large $\textrm{N}^{\prime\prime} $ leads to too straight branches and thus loses fidelity.

4. Results

We demonstrated the feasibility and advantages of our algorithm through imaging experiments of tungsten wire phantom, leaf skeleton phantom, and mouse dorsal and ear blood vessels in vivo. These experiments were conducted using a focused US transducer with center frequency of 50 MHz, NA of 0.44, focal length of 6.7 mm, and DOF of ∼300 µm [13]. The results are described in the following subsections. The number of iterations for R-L deconvolution (used for comparison) is 15 for tungsten wire experiments and 10 for leaf skeleton experiment and in vivo mouse experiments.

4.1 Tungsten wire phantom imaging

We designed two tungsten wire phantom imaging experiments to elaborate the effect of our algorithm. First, we imaged two crossed tungsten wires with an angle of ∼90° to quantify full width at half maximum (FWHM) and DOF. Secondly, we imaged two closely-located tungsten wires to calibrate the minimum resolvable distance.

4.1.1 Crossed tungsten wires

Two tungsten wires with diameter of 20 µm were evenly placed at four different OFDs of 0 mm, 0.3 mm, 0.6 mm and 0.9 mm (i.e., the axial distance of 6.7 mm, 7.0 mm, 7.3 mm and 7.6 mm, respectively, from the surface of the transducer), with one wire overlapping on the other. The four OFDs are denoted as D1 − D4, respectively. The lateral MAP images of raw data (i.e., $\textrm{RF}({x,y} )$ without further processing) are shown in Supplementary Fig. S1. Supplementary materials are available online. In the images in Supplementary Fig. S1, lateral resolution deteriorates as the $\textrm{abs}$ (OFD) of the wires increases. For convenient comparison of several algorithms, we selected the same regions with the same size of the images for display (indicated by the white boxes in Supplementary Fig. S1).

The workflow of our algorithm is shown in Fig. 3. PAM images of D2 are used for illustration. The $\gamma $ used in FA-SAFT is set to 0.2 (Eq. (19)). We also compared FA-SAFT using different $\gamma $ (see Supplementary Fig. S2). Larger $\gamma $ generates smaller features size (or FWHM), but too large $\gamma $ leads to artifacts. As a result, the optimal $\gamma $ is set to 0.2. For Fig. 3, $\mathrm{N^{\prime}} = 2$ and $\textrm{M} = 1$ are used for illustration. Note that $\mathrm{N^{\prime}} = 16$ (as recommended previously [14]) and $\textrm{M} = 4$ are used for other figures unless otherwise specified. $\textrm{M} = 4$ is chosen to balance the processed effects and the computation time (see Supplementary Fig. S3). Supplementary Fig. S3 shows that the initial phase shift ${\varphi _0}$ does not affect the D-MB results for the cases of $\textrm{M} = 4$ and $\textrm{M} = 8$. Experimentally measured PSF of 65 µm is used when applying the deconvolution algorithm [13]. Note that since FA-SAFT is to restore the lateral resolution in the out-of-focus region to the in-focus lateral resolution, the same lateral PSF can be assumed for AR-PAM images at different OFDs when applying deconvolution [13].

 

Fig. 3. Illustration of the workflow of our algorithm. The numbers inside circles here correspond to the 6 steps in Fig. 2. $\mathrm{N^{\prime}} = 2$ and $\textrm{M} = 1$ are used here for illustration. (a) MAP image of raw data at D2 (OFD of 0.3 mm). (b,c) By 1D SAFT along the horizontal and vertical directions, respectively. (d − f) (e) and (f) are the Fourier transform of (b) and (c), respectively. (d) is the merged result from (e) and (f). That is, (d) is ${\textrm{K}_{\textrm{merged} - \textrm{FA}}}({u,v} )$. (g) By FA-SAFT (i.e., $\textrm{R}{\textrm{F}_{\textrm{FA} - \textrm{SAFT}}}({x,y} )$). (h,i) The images in spatial domain after directional decomposition of (g). (j − l) (k) and (l) are 1D MB deconvolution results of (h) and (i), respectively. (j) is the final result (i.e., $\textrm{R}{\textrm{F}_{\textrm{Final}}}({x,y} )$). IFT, inverse Fourier transform.

Download Full Size | PDF

We then studied the effects of the directional term (Eq. (18) with and without ${\textrm{W}_{{\theta _n}}}({u,v,\omega } )$) and sharpening term (Eq. (18) with and without $\textrm{D}({u,v,\omega } )$) in FA-SAFT. PAM images of D2 are used for illustration. The results are shown in Fig. 4. Figures 4(a)−4(d) show the results as follows: without ${\textrm{W}_{{\theta _n}}}({u,v,\omega } )$ and without $\textrm{D}({u,v,\omega } )$ for Fig. 4(a) (termed pure SAFT); with ${\textrm{W}_{{\theta _n}}}({u,v,\omega } )$ and without $\textrm{D}({u,v,\omega } )$ for Fig. 4(b) (i.e., D-SAFT); without ${\textrm{W}_{{\theta _n}}}({u,v,\omega } )$ and with $\textrm{D}({u,v,\omega } )$ for Fig. 4(c) (i.e., FA-SAFT-dir0); with ${\textrm{W}_{{\theta _n}}}({u,v,\omega } )$ and with $\textrm{D}({u,v,\omega } )$ for Fig. 4(d) (i.e., FA-SAFT-dir1). We found that Fig. 4(d) achieves better resolution than Fig. 4(b). We also compared the performance of R-L deconvolution and D-MB deconvolution based on Figs. 4(b) and 4(d). Figures 4(e) and 4(f) are the results by R-L deconvolution from Figs. 4(b) and 4(d), respectively, and Figs. 4(g) and 4(h) are the results by D-MB deconvolution from Fig. 4(b) and 4(d), respectively. The results show that Fig. 4(f) (finer feature size and fewer artifacts) is better than Fig. 4(e), showing the advantage of the proposed FA-SAFT over D-SAFT for the following deconvolution. Similarly, Fig. 4(h) shows slightly smaller feature size than Fig. 4(g). Further, Fig. 4(h) shows smaller feature size, higher SNR, and fewer artifacts than Fig. 4(f), indicating the benefit of D-MB deconvolution over R-L deconvolution. Moreover, we compared the deconvolution results based on Figs. 4(c) and 4(d) (see Supplementary Fig. S4). The results show that Supplementary Fig. S4f has comparable resolution to Supplementary Fig. S4c (similarly, Supplementary Fig. S4e comparable to Supplementary Fig. S4b), meaning that both FA-SAFT-dir0 and FA-SAFT-dir1 suit for the following deconvolution. Overall, for SAFT, FA-SAFT outperforms D-SAFT. As for deconvolution, D-MB deconvolution presents better results than R-L deconvolution.

 

Fig. 4. Comparison of different SAFT and the following deconvolution methods. (a) By pure SAFT. (b) By D-SAFT. (c) By FA-SAFT-dir0. (d) By FA-SAFT-dir1. (e,f) R-L deconvolution from (b) and (d), respectively. (g,h) D-MB deconvolution from (b) and (d), respectively.

Download Full Size | PDF

Then, we compared the results by FA-SAFT-dir1 and D-SAFT in the out-of-focus region (D2 − D4). Besides, R-L deconvolution and D-MB deconvolution were applied and compared for the four OFDs (D1 − D4). The results are shown in Fig. 5, which are displayed in MAP images. Figure 5(a) shows the MAP images of raw data at D1 − D4. Figures 5(b) and 5(c) are the results by D-SAFT and FA-SAFT-dir1, respectively. Figures 5(d)−5(f) show the results by R-L deconvolution, 2D MB deconvolution (Eq. (17)), and D-MB deconvolution from Fig. 5(c), respectively. The FWHM and SNR of the images in Fig. 5 are summarized in Table 1. Note that “R-L”, “MB”, and “D-MB” in Table 1 represent the deconvolution results from Fig. 5(c) (FA-SAFT-dir1). Overall, the results in Fig. 5 and Table 1 for the results at D1 − D4 are consistent with Fig. 4. First, FA-SAFT-dir1 performs better than D-SAFT. The former has smaller FWHM and higher SNR. Secondly, compared with R-L deconvolution, D-MB deconvolution achieves better results with smaller FWHM, higher SNR, better pattern continuity, and fewer artifacts. Compared with MB deconvolution, D-MB deconvolution keeps the previous advantages except for larger FWHM at D2 and D3. Finally, we obtained FWHMs of 26 − 31 µm within DOF of 1.8 mm [= 0.9 mm×2], when considering below and above the focal plane, by D-MB deconvolution from FA-SAFT images.

 

Fig. 5. Results of the first tungsten wire experiment. (a) MAP images of raw data at D1 − D4. (b) By D-SAFT. (c) By FA-SAFT-dir1. (d) R-L deconvolution from (c). (e) 2D MB deconvolution from (c). (f) D-MB deconvolution from (c).

Download Full Size | PDF

Table 1. FWHM (µm) and SNR (dB)

View Table

4.1.2 Closely-located tungsten wires

Two wires with diameter of 20 µm were closely placed with a narrow angle between them. The phantom was imaged at OFDs of 0 mm (i.e., in-focus) and −0.45 mm (i.e., out-of-focus), respectively. The results are shown in Fig. 6. PAM images of raw data at the focal and out-of-focus regions are shown in Figs. 6(a) and 6(d), respectively. For the in-focus image, we directly applied R-L deconvolution (Fig. 6(b)) and D-MB deconvolution (Fig. 6(c)), respectively. For the out-of-focus image, we first applied FA-SAFT-dir1 (Fig. 6(e)). Then, Fig. 6(e) was further processed by R-L deconvolution (Fig. 6(f)) and D-MB deconvolution (Fig. 6(g)), respectively.

 

Fig. 6. Results of the second tungsten wire experiment. (a) PAM image of raw data at OFD of 0 mm. (b,c) R-L deconvolution and D-MB deconvolution from (a), respectively. (d) PAM image of raw data at OFD of −0.45 mm. (e) By FA-SAFT-dir1 from (d). (f,g) R-L deconvolution and D-MB deconvolution from (e), respectively. (h,i) 1D profiles of (a − c) along the two lines indicated by #1 and #2 in (a), respectively. (j,k) 1D profiles of (d − g) along the two lines indicated by #3 and #4 in (d), respectively. #1-#4: The position below which the two tungsten wires start to be distinguishable in (c), (a), (g), and (e), respectively. That is, (h − k) show the resolution of (c), (a), (g), and (e), respectively.

Download Full Size | PDF

The resolution is determined by checking the minimum distinguishable separation (i.e., minimum resolvable distance) of the two tungsten wires (≥6 dB for the contrast between the minimum of the two peaks and the valley). Figures 6(h)−6(k) show 1D profiles along the lines #1−#4 in Fig. 6, respectively (#1 and #2 at the same positions for Figs. 6(a)−6(c); #3 and #4 at the same positions for Figs. 6(d)−6(g)). Note that lines #1−#4 in Fig. 6 are chosen based on the above criterion to determine the resolution of Figs. 6(c), 6(a), 6(g), and 6(e), respectively. Similarly, the resolution of Figs. 6(b), 6(d), and 6(f) can be determined using the same criterion (related 1D profiles not shown). As a result, the measured resolutions of Figs. 6(a)−6(g) are 129 µm, 76 µm, 49 µm, 406 µm, 108 µm, 65 µm, and 46 µm, respectively. For OFD of 0 mm, the resolution is enhanced from 129 µm to 49 µm (∼2.6 times) before and after applying D-MB deconvolution (Fig. 6(a) vs. Fig. 6(c)). For OFD of −0.45 mm, the resolution is enhanced from 108 µm to 46 µm (∼2.3 times) before and after applying D-MB deconvolution (Fig. 6(e) vs. Fig. 6(g)). Besides, the resolution achieved by D-MB deconvolution is also better than that by R-L deconvolution (∼1.4 − 1.6 times improvement).

4.2 Leaf skeleton phantom imaging

To investigate the ability of our algorithm for more complex structures, we imaged a Banyan leaf skeleton phantom and applied our algorithm. Before imaging, the leaf was immersed in carbon ink to increase light absorption of the phantom. The results are shown in Fig. 7. PAM images were acquired at both the focal (Fig. 7(a)) and out-of-focus (Fig. 7(d)) regions. For the in-focus image, we applied the deconvolution algorithms directly (Figs. 7(b) and 7(c)). For the out-of-focus image, different SAFT and the following deconvolution were applied (Figs. 7(e)−7(i)).

 

Fig. 7. Results of leaf skeleton experiment. (a) PAM image of raw data at OFD of 0 mm (in-focus). (b,c) R-L deconvolution and D-MB deconvolution from (a), respectively. (d) PAM image of raw data at OFD of 0.25 mm (out-of-focus). (e − g) By D-SAFT, FA-SAFT-dir1, and FA-SAFT-dir0, respectively. (h,i) R-L deconvolution and D-MB deconvolution from (g), respectively. (j − l) 1D profiles of the images (d − i) along the three lines indicated by #1, #2, and #3 in (d), respectively.

Download Full Size | PDF

First, Fig. 7(c) presents much better image fidelity and smaller feature size than Fig. 7(b), demonstrating the great performance of D-MB deconvolution over R-L deconvolution. Secondly, Figs. 7(e)−7(g) are the results by D-SAFT, FA-SAFT-dir1, and FA-SAFT-dir0, respectively. The latter two images show less noise than the former. Similar to Supplementary Fig. S4, FA-SAFT-dir0 (Fig. 7(g)) is comparable to FA-SAFT-dir1 (Fig. 7(f)). Therefore, we chose to apply the deconvolution algorithms to the image by FA-SAFT-dir0 (Fig. 7(g)). The results are shown in Figs. 7(h) (R-L deconvolution) and 7(i) (D-MB deconvolution). Similar to Figs. 7(b) and 7(c), Fig. 7(i) shows better image fidelity and smaller feature size than Fig. 7(h). We also applied the deconvolution algorithms to the image by FA-SAFT-dir1 (see Supplementary Fig. S5). The comparable deconvolution results from FA-SAFT-dir1 and FA-SAFT-dir0 in Supplementary Fig. S5 are similar to those in Supplementary Fig. S4.

For FWHM comparison, we plot 1D profiles along the three lines with different directions in Figs. 7(d)−7(i). Note that the three lines are at the same positions in Figs. 7(d)−7(i). For convenience, the three lines are only indicated in Fig. 7(d) and denoted as #1, #2, and #3. The 1D profiles are shown in Fig. 7(j)−7(l) for lines #1, #2, and #3, respectively. FA-SAFT-dir0 and FA-SAFT-dir1 achieve better FWHM than D-SAFT. Besides, as expected, D-MB deconvolution produced smaller FWHM compared with R-L deconvolution. In addition, we rotated Fig. 7(a) by 45° to test the robustness of D-MB deconvolution to different orientations of the line pattern (see Supplementary Fig. S6).

4.3 In vivo mouse imaging

We further demonstrated our algorithm applied to in vivo mouse imaging. This animal experiment was conducted in conformity with the laboratory animal protocol approved by Laboratory Animal Care Committee of Shanghai Jiao Tong University. Mouse dorsal subcutaneous vessels of a BALB/c mouse (∼20 g, Slac Laboratory Animals) with hair removed were first imaged by our AR-PAM system.

Similar to Figs. 7(d)−7(i), the images of raw data (at OFD of ∼0.4 mm), by D-SAFT, by FA-SAFT-dir1, by FA-SAFT-dir0, R-L deconvolution after FA-SFAT-dir0, and D-MB deconvolution after FA-SAFT-dir0 are shown in Figs. 8(a)−8(f), respectively. Besides, for resolution comparison, we also plot 1D profiles in Figs. 8(g)−8(i) along the three lines in Figs. 8(a)−8(f). The results are overall similar to Fig. 7. First, FA-SAFT-dir0 (Fig. 8(d)) and FA-SAFT-dir1 (Fig. 8(c)) present lower noise than D-SAFT (Fig. 8(b)). Interestingly, FA-SAFT-dir0 (Fig. 8(d)) shows slightly higher SNR than FA-SAFT-dir1 (Fig. 8(c)). Therefore, the deconvolution based on Fig. 8(d) was conducted, and the results are shown in Figs. 8(e) and 8(f) by R-L deconvolution and D-MB deconvolution, respectively. Secondly, as expected, D-MB deconvolution (Fig. 8(f)) achieves better results than R-L deconvolution in terms of image fidelity and SNR. Similarly, the deconvolution applied to Fig. 8(c) was also conducted for comparison (see Supplementary Fig. S7). Overall, Supplementary Fig. S7f shows slightly higher SNR than the corresponding Supplementary Fig. S7c (similarly for Supplementary Figs. S7e vs. S7b), which is consistent with the results before applying deconvolution (Supplementary Figs. S7d vs. S7a). Besides, by checking Figs. 8(g)−8(i), FA-SAFT-dir1 and FA-SAFT-dir0 present higher resolution than D-SAFT, and D-MB deconvolution enables the smallest FWHM. Figures 8(j)−8(o) show zoom-in images of the dashed box regions (indicated in Fig. 8(a)) in Figs. 8(a)−8(f), respectively. Note that the relatively weak PA amplitude of some small vessels in the dashed box region in Fig. 8(f) does not mean that these small vessels disappear after applying deconvolution, but is due to the adopted color map for display. As can be seen in Fig. 8(o), these small vessels can be clearly visualized with good SNR.

 

Fig. 8. Results of in vivo mouse experiment (mouse dorsal blood vessels). (a) PAM image of raw data at OFD of ∼0.4 mm (out-of-focus). (b − d) By D-SAFT, FA-SAFT-dir1, and FA-SAFT-dir0, respectively. (e,f) R-L deconvolution and D-MB deconvolution from (d), respectively. (g − i) 1D profiles of the images (a − f) along the three lines indicated by #1, #2, and #3 in (a), respectively. (j − o) Zoom-in images of the dashed box region (indicated in (a)) in (a − f), respectively. Note that the box indicates the same regions for (a − f).

Download Full Size | PDF

For further validation, we tested our algorithm on two other in vivo AR-PAM data: (i) Another image of mouse dorsal blood vessels in vivo [13] (see Supplementary Fig. S8); (ii) mouse ear blood vessels in vivo imaged at the focal and out-of-focus regions (see Supplementary Figs. S9 and S10). Supplementary Fig. S9 shows the images of raw data and processed by SAFT and deconvolution, and Supplementary Fig. S10 shows 1D profiles and zoom-in images indicated in Supplementary Fig. S9 for better comparison and display. The results are promising. For Supplementary Fig. S8, overall, FA-SAFT performs better than D-SAFT (notably high resolution in Supplementary Fig. S8d over Fig. S8b), and D-MB deconvolution produces better results than R-L deconvolution (especially pattern continuity in Supplementary Fig. S8f over Fig. S8e). Similarly, for Supplementary Fig. S9, FA-SAFT enables higher resolution than D-SAFT (FA-SAFT vs. D-SAFT: Supplementary Figs. S9f vs. S9e). Besides, compared with R-L deconvolution, D-MB deconvolution renders finer feature size (D-MB vs. R-L: Supplementary Figs. S9c vs. S9b; Supplementary Figs. S9i vs. S9h) and better pattern continuity (D-MB vs. R-L: Supplementary Figs. S10f vs. S10e; Supplementary Figs. S10h vs. S10g). Notably, the high resolution by D-MB deconvolution is clearly evidenced by Supplementary Figs. S10a − S10d, where D-MB can better resolve two closely-located vessels. Similar to Fig. 8, the relatively weak PA amplitude of some small vessels in the dashed boxes in Supplementary Fig. S9 is due to the adopted color map for display. These small vessels can be clearly visualized with good SNR in Supplementary Figs. S10e − S10h.

5. Discussion

We proved that our algorithm can improve the lateral resolution, SNR, and pattern fidelity in both phantom and in vivo imaging experiments. The great performance can be attributed to two parts, FA-SAFT and D-MB deconvolution. More details are discussed in the following.

5.1 FA-SAFT

Note that besides SAFT, a compressed sensing-based approach can also be used to enhance images in the out-of-focus region in AR-PAM [29].

FA-SAFT produced higher resolution and sharper edges than D-SAFT. The resolution improvement in FA-SAFT can be explained from the sharpening term, which is newly introduced in FA-SAFT and acts as an adaptive high-pass filter. For better illustration, we compared the processed results and their spectra (k-space) by high-pass filtering (HPF), D-SAFT (i.e., $\gamma $ = 0), and FA-SAFT-dir1 of the crossed tungsten wire images over D2 − D4 (Section 4.1.1) (see Supplementary Figs. S2 and S11). In Supplementary Fig. S2, we compared the results by HPF using different cut-off frequencies and by FA-SAFT-dir1 using different $\gamma $. We found that HPF with higher cut-off frequency produces smaller feature size, but HPF with too high cut-off frequency of ≥3 units (larger units correspond to higher cut-off frequency) results in evident artifacts. Thus, the optimal cut-off frequency is set to 2 units. As mentioned previously, the optimal $\gamma $ is set to 0.2. Therefore, we then compared the optimized results by HPF with cut-off frequency of 2 units and by FA-SAFT-dir1 with $\gamma $ = 0.2. As can be seen, FA-SAFT-dir1 with $\gamma $ = 0.2 shows well-suppressed background (or high SNR).

In Supplementary Fig. S11, we first checked the case of D2. Compared with HPF, FA-SAFT-dir1 did not filter out all the low-frequency components but adjusted their weights adaptively, as indicated by the red arrows in Supplementary Fig. S11. Besides, compared with [HPF] − [D-SAFT] (“−” means difference), [FA-SAFT-dir1] − [D-SAFT] can more selectively enhance the images along the directions of the tungsten wires (i.e., increased weights in the spectra of FA-SAFT-dir1 along the two diagonal directions corresponding to the main signals of tungsten wires), showing that the sharpening term in FA-SAFT-dir1 can better reflect directional features of the data compared with HPF. In other words, the sharpening term in FA-SAFT-dir1 can increase the weights of high-frequency components adaptively. Here the “adaptive” means that the sharpening term changes accordingly with the 1D SAFT results as in Eq. (19) and enhances the spectrum directionally. Specifically, for example, as indicated by the green arrows in Supplementary Fig. S11, the spectrum not along the directions of the tungsten wires are more suppressed in [FA-SAFT-dir1] − [D-SAFT] than in [HPF] − [D-SAFT]. As can be seen in Supplementary Fig. S11, the same characteristics are observed in D3 and D4.

Therefore, Supplementary Fig. S11 would explain that compared with HPF, FA-SAFT-dir1 can enhance resolution (having similar high-frequency components) and without sacrificing SNR (enhancing the spectrum directionally, as detailed in the previous paragraph). Similarly, Supplementary Fig. S11 would also explain compared with D-SAFT, FA-SAFT-dir1 enables higher resolution (keeping more high-frequency components, as shown in the columns of “FA-SAFT-dir1” vs. “D-SAFT”) and higher SNR (enhancing the spectrum directionally, as shown in the column of [FA-SAFT-dir1] − [D-SAFT]).

Another variable in FA-SAFT is the number of directions for 1D SAFT ($\mathrm{N^{\prime}}$ in Eq. (18)). As discussed in D-SAFT [14], using more directions (larger $\mathrm{N^{\prime}}$) is useful to improve the resolution of branches in different orientations. However, when $\mathrm{N^{\prime}}$ is large enough ($\ge 16$), the image quality may not improve anymore. Therefore, $\mathrm{N^{\prime}}$ is also set to 16 in our demonstrations.

Overall, FA-SAFT-dir0 produces comparable or slightly better results than FA-SAFT-dir1 (e.g., slightly better for Fig. 8 and Supplementary Fig. S8). Besides, either FA-SAFT-dir0 or FA-SAFT-dir1 performs better than D-SAFT for the demonstrations in this work (Figs. 4 − 8). This may imply that the sharpening term plays a more important role than the directional term in Eq. (18). On the other hand, the following deconvolution (either R-L or D-MB) exhibits similar trends (e.g., similar results after deconvolution are obtained given similar image quality of FA-SAFT-dir0 and FA-SAFT-dir1).

SAFT in AR-PAM has been continuously improved to deal with more complicated patterns, especially for microvascular patterns, to meet the requirements of in vivo imaging applications. By considering 1D sparsity of the line pattern, A-SAFT and D-SAFT were proposed [10,14]. However, the former requires successful recognition of diverse orientations of vessel branches prior to applying A-SAFT, while the latter can better recover the vessel branches by utilizing the weighted 1D SAFTs. Based on D-SAFT, we proposed FA-SAFT to further improve resolution and SNR by introducing the sharpening term, as detailed previously. Deconvolution can perform better given the original input images with better image quality. Taking advantage of the high-quality image by FA-SAFT, D-MB deconvolution using FA-SAFT images outperforms that using D-SAFT images.

5.2 D-MB deconvolution

After FA-SAFT, deconvolution algorithms further improve the resolution and SNR. R-L deconvolution algorithm was initially designed for astronomy applications, where the targets are usually point-like objects. Compared to R-L deconvolution and MB deconvolution, the proposed D-MB deconvolution is more suitable for the line pattern because 1D sparsity of the line pattern is considered in D-MB deconvolution. To our knowledge, the directional deconvolution method for handling the line pattern is proposed for the first time. Besides the line pattern, D-MB deconvolution with $\textrm{M}$≥4 can well preserve the point features (see Supplementary Fig. S12). Another parameter in D-MB deconvolution is the sparsity regularization term, which contributes to the pattern centralization (i.e., resolution enhancement) and noise reduction to an extent (see Supplementary Fig. S13). The sparsity regularization in D-MB deconvolution may also contribute to the artifacts reduction. As can be seen, noise and artifacts are much reduced in Supplementary Fig. S13f ($\lambda $ = 1.0) compared with Supplementary Fig. S13a ($\lambda $ = 0; i.e., without the sparsity regularization). Finally, our algorithm (FA-SAFT and D-MB deconvolution) achieved lateral FWHMs of 26 − 31 µm when imaging 20 µm tungsten wires (Fig. 5), which is better than measured PSF of 65 µm. Besides, lateral resolution of 46 − 49 µm was obtained when resolving two closely-located tungsten wires (Fig. 6), which is better than original resolution of 129 − 406 µm. Note that the “sparsity” here does not mean that the pattern of the input PA image has to be sparse in order to apply the D-MB deconvolution. It works well in Figs. 7, 8, and S8, which proves that it is applicable to the AR-PAM images with a certain degree of pattern density.

To improve 3D resolution for the line pattern of blood vessels in AR-PAM images, the proposed D-MB deconvolution can be applied to 2D lateral slices at each axial position, and then 1D R-L or 1D MB deconvolution can be applied along the axial direction [13]. Considering that the line pattern of blood vessels mainly appears over the 2D lateral slices, it would be reasonable to apply 1D deconvolution along the axial direction without dealing with the directional characteristic. Besides, such a method for 3D deconvolution can save computation time and reduce memory usage [13].

5.3 Other Discussion

As AR-PAM is a 3D imaging modality, we evaluated the axial resolution changes when using our algorithm (see Supplementary Fig. S14 and Supplementary Table S1). We found that axial resolution after FA-SAFT and D-MB deconvolution is comparable or slightly improved. SAFT is typically not expected to improve the axial resolution [13]. On the other hand, in Supplementary Fig. S14, D-MB deconvolution is applied to 2D lateral slices at each axial position (multiple-times deconvolution for many slices). Because deconvolution is not applied along the axial direction, the resolution enhancement along the axial direction is not obvious. Note that except Supplementary Fig. S14, D-MB deconvolution is performed over the MAP images (e.g., Figs. 7(i) and 8(f)).

The robustness of the algorithm to different SNR is also an important factor worthy to discuss. Indeed, sufficient SNR is essential for deconvolution. Fortunately, in our algorithm, the first step, FA-SAFT, helps improve SNR for out-of-focus images prior to applying the second step, D-MB deconvolution. In Supplementary Fig. S15 and Supplementary Table S2, the effect of SNR on the resolution enhancement by FA-SAFT and D-MB deconvolution is evaluated. Note that the resolution in Supplementary Table S2 is determined using a criterion similar to Fig. 6. When artifacts become too evident to determine resolution, the resolution is regarded as “not applicable (N/A).” As can be seen, when the PAM image of raw data has peak SNR (PSNR) down to 20 dB (Fig. 6(d) with noise added), FA-SAFT and D-MB deconvolution can still enhance resolution from N/A to 54 µm. On the other hand, in Supplementary Fig. S16 and Supplementary Table S3, the effect of SNR on the resolution enhancement by D-MB deconvolution alone is investigated. The determination of resolution in Supplementary Table S3 is the same as Supplementary Table S2 mentioned above. As can be seen, when the PAM image by FA-SAFT (prior to applying deconvolution) has PSNR down to 21 dB (Fig. 6(e) with noise added), D-MB deconvolution can still enhance resolution from N/A to 54 µm, while R-L deconvolution suffers artifacts and the resolution is still N/A. Overall, in Supplementary Fig. S16, D-MB deconvolution produces better image quality (i.e., fewer artifacts) than R-L deconvolution.

Another important factor is the performance of the proposed algorithm when applied to images with penetration depth in tissue. We conducted additional in vivo AR-PAM imaging of subcutaneous blood vessels around the lower abdomen region of a mouse. As can be seen in Supplementary Fig. S17, the proposed algorithm demonstrates good performance for the in vivo AR-PAM images with penetration depth of ∼0.32 mm. Besides, although the skin surface is not flat and the blood vessels distribute over a large axial range, our algorithm applied for a large axial range of 1.35 mm (OFD from −0.525 mm to 0.825 mm) shows satisfactory results. We further and intentionally sought deeper blood vessels in AR-PAM imaging experiments. As can be seen in Supplementary Fig. S18 (notably the blood vessel V3), the proposed algorithm still demonstrates the effectiveness of resolution enhancement for the in vivo AR-PAM image of deep vessels up to ∼1.3 mm below the skin surface.

In the experiments described above, the penetration depth of ∼1.3 mm was demonstrated. AR-PAM is able for deeper penetration in tissue. In this case, however, SNR degrades, and low SNR adversely affects the performance of the proposed algorithm, as described previously (Supplementary Fig. S15, Table S2, Fig. S16, and Table S3). To advance the penetration of in vivo AR-PAM when applying our algorithm, one could consider signal averaging and/or contrast agents to secure sufficient SNR, which would ensure in vivo applicability of our algorithm.

In Section 4.1, the tungsten wires were immersed in non-scattering water media and placed at different OFDs (e.g., 0 − 0.9 mm), and the performance of the algorithms at different OFDs can be conveniently compared (Table 1 and Figs. 5 and 6) by excluding the scattering effect. Similar experiments were conducted also in non-scattering water media in a previous study of resolution enhancement in AR-PAM [30]. Similar to the above discussion, strong scattering leads to low SNR, which degrades the performance of the proposed algorithm. Therefore, for scattering media, approaches to boost SNR are needed to ensure the performance of our algorithm.

The computation time of our algorithm is another important aspect. For FA-SAFT, the computation time is mainly determined by the number of directions for 1D SAFT ($\mathrm{N^{\prime}}$) and is proportional to $\mathrm{N^{\prime}}$, which is similar to the case of D-SAFT (Supplementary Table S4). That is, the inclusion of the sharpening term (specifically, for FA-SAFT-dir1) in the calculation will increase the computation time little (the time increase of ∼0.2 seconds in Supplementary Table S4 for $\mathrm{N^{\prime}}$ = 16). Therefore, once the same $\mathrm{N^{\prime}}$ is chosen for D-SAFT and FA-SAFT, they will have similar computation time (Supplementary Table S4). For D-MB deconvolution, more computation time is required compared with R-L deconvolution (Supplementary Table S5). Fortunately, the computation time of D-MB deconvolution is less than that of FA-SAFT, and thus, the overall computation time of our algorithm is still at the same time scale as FA-SAFT (FA-SAFT of ∼72 s ($\mathrm{N^{\prime}}$ = 16) and D-MB deconvolution of ∼8 s ($\textrm{M}$ = 4) in Supplementary Table S4 and Table S5, respectively).

6. Conclusions

Tailored for AR-PAM imaging applications of microvasculature, we proposed a novel algorithm consisting of FA-SAFT and D-MB deconvolution for AR-PAM to achieve evident image enhancement in both the focal and out-of-focus regions. Specifically, FA-SAFT improved resolution and SNR, and D-MB deconvolution further enhanced resolution and SNR while keeping image fidelity. As a result, compared with other SAFT and/or deconvolution algorithms, our algorithm achieved great performance in terms of resolution, SNR, and image fidelity. Notably, we demonstrated promising results for in vivo images of mouse blood vessels by our algorithm given the images of raw data with sufficient SNR. Our work opens up new opportunities for microvascular AR-PAM imaging applications. The proposed algorithms may benefit other imaging modalities to handle the line pattern, e.g., FA-SAFT for ultrasound imaging and D-MB deconvolution for optical coherence tomography angiography.