Open Access
Add to BibSonomy
Abstract
Optical-resolution photoacoustic microscopy has been validated as an ideal tool for angiographic studies. Quantitative vascular analysis reveals critical information where vessel segmentation plays the key step. The comm-only used Hessian filter method suffers from varying accuracy due to the multi-kernel strategy. In this work, we developed a Hessian filter-assisted, adaptive thresholding vessel segmentation algorithm. Its performance is validated by a digital phantom and in vivo images which demonstrates a superior and consistent accuracy of 0.987 regardless of kernel selection. Subtle vessel change detection is further tested in two longitudinal studies on blood pressure agents. In the antihypotensive case, the proposed method detected a twice larger vasoconstriction over the Hessian filter method. In the antihypertensive case, the proposed method detected a vasodilation of 21.2%, while the Hessian filter method failed in change detection. The proposed algorithm may further push the limit of quantitative imaging on angiographic applications.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
Corrections
12 August 2022: A correction was made to the corresponding author.
1. Introduction
Photoacoustic imaging (PAI), a hybrid technique detecting tissue absorption properties ultrasonically, has been a rapidly developing modality for preclinical and clinical studies in the past decades. PAI enables investigation of structural and functional information at various scales ranging from organelles to organs [1]. As a main embodiment of PAI, optical-resolution photoacoustic microscopy (OR-PAM) is capable of achieving high-contrast, high-resolution and high-sensitivity imaging by confocal and coaxial configuration of the tightly focused optical excitation and high frequency acoustic detection. With tens of years development on this technique, it is progressively mature and widely used in different biomedical studies [2–5]. Regarding the high absorption to the visible light, hemoglobin serves as an ideal endogenous contrast for OR-PAM imaging of vasculatures in vivo. The inherent advantages of OR-PAM and associated endogenous contrasts make it a powerful tool for angiographic studies [5–7].
Vascular features are typically accompanied with physiological processes. Multiple parameters extracted from vascular morphologies are widely used for diagnosis [8]. Diameter distribution change in the brain reflects cerebral blood flow regulations [9]; vessel area fraction (VAF) change is a general sign of vasoconstriction or vasodilation, indicating physiological responses to pressure medicines [10]. In order to acquire these parameters, an accurate vessel segmentation method is required to extract vessel networks from grayscale vascular images.
In previous studies, simple thresholding approaches have been applied to microvessel segmentation, among which the Otsu’s thresholding (OT) shows promising results [11]. Although this single-value thresholding approach is likely to either miss small vessels or overly segment large ones [12,13], it is still widely in use because of its simplicity. Numerous advanced vessel segmentation techniques have been developed, which can mainly be classified into three categories: deep learning methods, tracking methods and filter-based methods. The deep learning methods train artificial neural networks to classify images into vessels and background with pre-labelled data. It avoids heavy loads of calculation during segmentation after training the networks, but a large number of pre-labelled data is necessary during the training [14,15]. The tracking methods trace the vessels iteratively by starting from a pre-selected vessel point, and locating the next until the entire network is extracted. The vessels are segmented during the tracking process. The tracking methods perform well on detecting large vessels, but could fail on small vessels due to their weak signal and small scales [16]. The filter-based methods generally enhance the vessel contrast. It extracts vessel structures by applying a filter, which increases the vessel contrast when the vascular size and orientations are matched by the filter. Hessian filter (HF) is still one of the mostly used techniques since its report 20 years ago, for its high sensitivity to tubular structures and low computational cost [17]. However, the vessel distortions brought along require further explorations on the improvements of this technique [18]. Yang et al. and Zhao et al. developed a Hessian filter-based method which achieved better vessel enhancement on PAM images [19,20]. Khan et al. reported an approach by fusing two segmentation results of different filter scales to mitigate the distortions [21]. Although these studies have pushed forward the limit of the HF method, the manually selected scale ranges would in any way introduce errors to the results.
Different from large vessels in OR-PAM images, small ones generally possess weak signal, which is difficult to extract. Also, user-selected parameters in the HF algorithm, like scale ranges, may further introduce errors due to different user experience. In this work, a robust vessel segmentation algorithm of high accuracy for PAM vascular images is proposed. It applies HF for small vessel extraction and fuses the results with large vessels as well as crossing points extracted by OT. An upward region growing (URG) process is followed to retore the lost connectivity caused by HF. By taking the full diameter at half maximum (FDHM) at each pixel along the vessel skeletons of the binary branches, the distorted vessel sizes are corrected. The performance of the algorithm is validated by a custom-built 3D digital phantom and in vivo mouse cortex as well as ear vasculatures. The algorithm is further tested with two longitudinal studies on vascular responses to a vasoconstricting and a vasodilating drug.
2. Materials and methods
2.1 Algorithm workflow
The overall flowchart and corresponding intermediate example images of the proposed algorithm is shown in Fig. 1. The segmentation algorithm consists of four steps: preprocessing, vessel enhancing, thresholding and segmenting, as shown in Fig. 1(a).
Fig. 1. The flowchart and corresponding intermediate example images of each step in the proposed algorithm. (a) Flowchart of the proposed algorithm. (b) Corresponding intermediate example images after each step; (i)-(v) important intermediate steps. HFE, high frequency emphasis; CLAHE, contrast limit adaptive histogram equalization; RC, radius circular; FDHM, full diameter at half maximum.
2.1.1 Preprocessing
In the preprocessing step, a mean filter and a median filter are applied for denoising, followed by a contrast enhancing process. In a maximum amplitude projection (MAP) vessel image, system noise can be recorded by imaging without targets. The system noise is subtracted and the top 1% pixel values are saturated to increase the contrast. A high-frequency emphasis (HFE) filtering is then performed [20], which emphasizes the high frequency components, referring to small vessels, while maintaining the rest components of the image. The emphasized image is subsequently processed by a contrast limit adaptive histogram equalization (CLAHE) operator, which enhance the local contrast to further differentiate the small vessels from background [22]. Notably, the background noise is also enhanced by CLAHE, as shown in Fig. 1(b), but would be suppressed in the next steps.
2.1.2 Vessel enhancing
The vessel enhancing step is achieved by HF [17]. In this step, the vessels in images are probed by a Hessian matrix and the tubular structures in the image are extracted by calculating the eigenvalues of the matrix. The Hessian matrix at pixel $({x,y} )$ is defined as:
The HF extracts vessels and suppresses non-tubular background signal. An example of the filtered image is shown in Fig. 1(b)(i). However, it also introduces vessel distortions due to the user-selected kernel scales and fails to detect the branching and crossing points because of its sensitivity solely to tubular structures. Two types of distortions are brought into the process, i.e., the vessel size distortion and the disconnectivity of adjacent vessels. As for the vessel size distortion, when the filter scale range covers all vessel sizes, the smaller vessels would be dilated by large filters; when the scale range does not cover large vessels, the filters only enhance the edges of large vessels while suppressing the signal inside, forming ‘hollow’ large vessels [18]. The shape of the second-order partial derivative of the Gaussian kernel inherently determines that it only enhances the pixels of the vessel and suppresses those next to it. When a small vessel of weak signal is connected to a large one, the HF filtering would suppress the connection, leading to disconnectivity [17]. In addition, the branching and crossing points are less likely to be extracted because they are not tubular in shape, to which the HF is not sensitive [23]. These distortions from HF are the key reasons for its inaccuracy in quantification, and will be corrected in the following steps.
2.1.3 Thresholding
In the thresholding step, an adaptive thresholding strategy converts the Hessian filtered image into a binary one and the OT extracts the large vessels as well as branching and crossing points at the same time. The adaptive thresholding is implemented with two Gaussian filters with the sizes of 3 pixels and 11 pixels, both of which convolve with the Hessian filtered image to form two smoothed maps. At each pixel, the smaller value of the two maps serves as the corresponding threshold [24]. The idea behind this strategy is to set a lower threshold for small vessels and a higher one for large vessels, to practically preserves small vessels and avoid over-segmentation on large ones. The OT is simultaneously applied on the HFE filtered image to acquire a binary map of high pixel-value components, referring to large vessels as well as branching and crossing points in OR-PAM MAP images. Small and large vessels, branching and crossing points can all be extracted by fusing the two binary maps, and the ‘hollow’ large vessels are also filled up. Note that, the vascular size distortion introduced by the HF and adaptive thresholding step will be corrected in the following steps, so the choice of the two Gaussian filter sizes is not necessarily strict as long as the target vessels are preserved.
2.1.4 Segmenting
In the segmenting step, the lost connectivity is restored with FDHM vessel restoration followed by URG, and the vessel size distortion is finally corrected by FDHM segmentation.
In the FDHM vessel restoration step, the vascular skeleton and radius map are firstly calculated from the fused binary map after the thresholding step. The skeletonization are implemented by iteratively removing the outmost layer of pixels until unit thickness, as shown in Fig. 1(b). To acquire the radius map, the distance field is calculated [20], and the longest distance between pixels within the vessel and the nearest vascular boundary is considered as the vessel radius. The falsely acquired skeleton branches shorter than the vessel radius are discarded. To restore the vessel branches by FDHM, a circular mask with the size of the radius (RC) calculated above plus 15 µm, which is 3 pixels in this study ($Mask\; Radius = Calculated\; Radius + 3\; pixels$) is employed. It is centered at each point along the skeleton, as shown in Fig. 1(b)(ii). The mask is used to take out the corresponding pixels of the HFE filtered image in the circle area, in which the pixels of larger values than half of the maximum are kept as the vessel segment while the rest are removed, as shown in Fig. 1(b)(iii). The entire vasculature is restored by repeating this process across all the skeleton points and adding all the segments up. The 3 extra pixels added on the RC mask radius ensures the full coverage of the target vessel. The full width at half maximum of a vessel’s crossing profile is widely accepted as its segmentation [25]. The FDHM technique here is inspired by this idea and implemented in a circle mask instead. The vessel size distortions are corrected after this step.
In the URG step, a modified region growing method is used to regain the connectivity between main vessels and adjacent ones. Firstly, the binary image from the last step is mapped onto the HFE filtered image, so that all the pixel values outside the foreground are set to 0, and therefore the vessel branches have clear boundaries. Secondly, the boundary layer pixels of a vessel branch are set as the initial seed points, and the URG only grows upwards to the pixels having larger values than the mean of the seed points. An example image after URG is shown in Fig. 1(b)(iv). When a small vessel branch is originally connected to the adjacent large one, the pixels in between have larger values than the seed points, so it grows to the large one and regain the connectivity. If the small vessel is not originally connected to the large one, the pixels in between have smaller values and thus no growth towards it. It should be noted that, no stop threshold is set for URG in this study, which allows full growth for connectivity. The full growth without a stop threshold also causes vessel dilations, which would be corrected in the next step.
In the FDHM segmentation step, it is basically a repeat of the FDHM vessel restoration step but with the updated binary map from the URG step, and the errors introduced in the URG step are corrected in this step. An example image of the final result is shown in Fig. 1(b)(iv). Until here, the vessel segmentation is finished and all the distortions mentioned above are corrected.
2.2 3D digital phantom and simulations
In order to evaluate the performance of the proposed algorithm, a 3D digital phantom is created to generate vascular images. The digital phantom is created from a 2D binary example vessel image from the open-source k-Wave Matlab Toolbox [26]. The skeleton and radius map are firstly calculated from the image. A ball-shaped vessel volume of the size from the radius map in 3D is centered at each pixel along the skeleton. By adding up all the balls across the skeleton, the 3D digital phantom is created, as shown in Fig S. 1(a).
The simulation parameters of the light and transducer are designed according to our OR-PAM system [27]. The laser wavelength is set to 532 nm and the numerical aperture (NA) is set to 0.1. The central frequency of the transducer is 30 MHz, and the one-way -6 dB bandwidth is 45 MHz, respectively. The unit distance between two adjacent pixels is 5 µm. For simplicity, both the light attenuation and acoustic attenuation are not accounted in the simulations [28].
The simulations run in 2D, where the vessel imaging follows a pattern of raster scanning, which images the phantom vessel slice by slice, as shown in Fig S. 1(b). Specific steps are described as follows: (1) Vessel cross-section slices are extracted from the 3D digital phantom; (2) Each vessel cross-section shifts from the left to the right across the laser focal spot, and is imaged pixel by pixel;(3) All the detected B-scans of the cross-section slices are converted into MAP profiles and then formed the 2D vessel image. The acquired MAP image is shown in Fig. 2(b). Gaussian noise with 0.037 mean and $1.49 \times {10^{ - 5}}$ variance is added to the normalized MAP image. The noise parameters are determined according to the in vivo mouse cortex vessel image in Fig. 3(a), where the mean noise value is 560 and the maximum value is 15120. The variance is measured from the background noise after normalization. The noise mean does not equal to 0 due to the fixed level of system noise. A second digital phantom simulation is implemented to further validate the accuracy and robustness of the proposed method under different vessel structures, as shown in Fig S.2.
Fig. 2. Ground truth and MAP images of a digital vessel phantom and corresponding segmentation results by the OT, HF, and proposed methods with different scale ranges. (a) Ground truth. (b) MAP image. (c) Segmentation result of the OT method. (d-f) Segmentation results of the HF method with scale ranges of [1 1.5], [1 3] and [1 10], respectively; red arrows: under- and over-segmentation. (g-i) Segmentation results of the proposed method with scale ranges of [1 1.5], [1 3] and [1 10], respectively. MAP, maximum amplitude projection; OT, Otsu’s thresholding; HF, Hessian filter.
Fig. 3. Illustration of the intermediate images after major steps of the proposed method and vessel profiles of the original, OT, HF and proposed methods. (a) Original mouse cortex vessel image; yellow arrows: vessel structures experiencing distortions during the segmentation process. (b) Result of the OT method; red arrow: under-segmented vessel. (c) Result of the HF method; red arrow: under-segmented branching and crossing points and over-segmented vessel branch. (d) Fusion image of the OT and HF results; red arrows: vessel disconnectivity caused by HF. (e) Results after FDHM vessel restoration and URG; red arrows: dilated vessel branch caused by URG. (f) Final result after FDHM segmentation; green arrows: corrected vessel structures. OT, Otsu’s thresholding; HF, Hessian filter; FDHM, full diameter at half maximum; URG, upwards region growing.
2.3 Performance evaluation
Accuracy (Acc), Sensitivity (Sn), Specificity (Sp) are three commonly used parameters to evaluate the segmentation performance. Acc provides the overall performance evaluation; Sn reflects the effectiveness of positive pixels detection; Sp specifies the effectiveness of negative pixels detection. The three metrics are defined as follows:
where TP, TN, FP and FN are true positive, true negative, false positive and false negative, which describe correctly segmented positive pixels, correctly segmented negative pixels, falsely detected positive pixels and falsely detected negative pixels, respectively.The 2D projection of the 3D digital phantom serves as the ground truth, as shown in Fig. 2(a). The proposed method is compared with two widely used segmentation techniques: OT and HF. The robustness of the proposed method is further tested with different scale ranges in comparison to HF. Three scale ranges of $= [{1\; 1.5} ],\; [{1\; 3} ]\; \textrm{and}\; [{1\; 10} ]$ are tested, which correspond to the conditions of under, normal, and over-segmentation.
2.4 OR-PAM system
To validate the proposed algorithm on in vivo data, we used a 2nd generation OR-PAM system described previously [27,29]. The laser wavelength is 532 nm and the NA of the objective is 0.1, respectively. The one-way -6 dB bandwidth of the transducer is 45MHz centered at 30MHz. The system provides a lateral resolution of 3 µm. The imaging is implemented in 2D raster scanning by a step motor stage, with a step size of 5 µm. All the parameters applied here is consistent with those in the simulations.
2.5 In vivo vessel segmentation
In vivo experiments were performed on Swiss Webster mice (male; 8–9 weeks old; ∼30 grams in weight). The protocol was approved by the Institutional Animal Care and Use Committee (IACUC) of Tsinghua University, following the ARENA/OLAW IACUC Guidebook [30]. The details of animal experiments can be found in our previous publication [27]. Before imaging, the mouse was anesthetized with 1.5% isoflurane. A 1×1 mm2 region of interest (ROI) was imaged both on the mouse cortex and ear vessels. The vascular images were then segmented by the proposed and HF methods with the scale ranges of [1 1.5], [1 3] and [1 10]. The Gaussian kernels of the three selected ranges can barely, mostly and overly cover all the target vascular sizes. Diameter distribution and VAF are calculated to evaluate the segmentation results. VAF is defined as:
where the vessel area is the area of the segmented foreground and the observed area is the area of ROI.2.6 Longitudinal studies
The performances of the proposed method and HF method are compared in longitudinal studies. Longitudinal observations on the mouse ear vessel changes in response to epinephrine (EP) and sodium nitroprusside (SNP) were implemented. EP is an antihypotensive drug that causes vasoconstriction and SNP is an antihypertensive drug that leads to vasodilation [27,31–33]. Prior to the drug injection, baseline images were recorded. The ROI of the mouse ear is 3×3 mm2. A bolus of 4 µg EP solution and a bolus of 75 µg SNP solution were intravenously administered to mice and continuous observations started immediately after injection for 30 min. The dose of 75 µg SNP in this study is referred to [27]. Meanwhile, the intravenously administered dose of 4 µg EP is adjusted from the subcutaneously injected dose of 5 µg EP [34]. The diameter changes in representative vessels were quantitively analyzed by the proposed method and HF method.
It has always been a challenge to obtain the ground truth (GT) of in vivo vessel diameters. General approaches either apply manually delineated vessels by experienced specialists or use well-established algorithms on grayscale images [16,35]. For OR-PAM MAP images, the vessels are observed in different sizes when various contrasts apply due to the nature of intensity continuity of the images. In this case, it may introduce errors in different extent for manual segmentation. Therefore, a line ray is centered at and rotated around each pixel along the vessel centerlines for 180° with an incremental angle of 5°, and the full width at half maximum (FWHM) of the grayscale ray profiles at each angle are taken as the ray lengths, of which the shortest one is considered as the GT diameter in this study. This approach is used to validate the performances of the proposed and HF method on diameter change detection.
3. Results
3.1 Simulation results
The OT, HF and proposed segmentation results on the digital phantom are demonstrated in Fig. 2. The small vessel ends of weak signal are missed by the OT method, as shown in Fig. 2(c). The HF method with $= [{1\; 1.5} ]$ fails to extract several branching points, as indicated by the red arrows in Fig. 2(d), demonstrating that small scale range causes under-segmentations. The HF method with $= [{1\; 3} ]$ shows dilation on small vessel ends, indicating the vessel enlargement effects from large Hessian kernels, as marked by the red arrows in Fig. 2(e). The HF method with $= [{1\; 10} ]$ results in serious distortions in both size and shape, as shown in Fig. 2(f). The proposed method with the three scale ranges shows close results without apparent distortions, as shown in Fig. 2(g)-(i). The Acc, Sn and Sp values of different methods are listed in Table 1.
Table 1. The accuracy (Acc), sensitivity (Sn) and specificity (Sp) of the OT, HF and proposed methods under different scale ranges.
The Acc, Sn and Sp of the OT method are 0.986, 0.715 and 1 respectively, showing that a portion of the vessel network is not detected but all segmented pixels fall in the true vessel area, which agrees well with the missing vessel ends in Fig. 2 (c). The HF method is highly dependent on the scale range selection. The Acc decreases from 0.985 to 0.984 when the scale range extends from [1 1.5] to [1 3], where the vessel enlargement is more significant than the mitigation of under-segmentation, as suggested by the increase of Sn from 0.766 to 0.920 and decrease of Sp from 0.996 to 0.987, respectively. Acc further decreases to 0.950 with the large scale range of [1 10] due to the serious vessel dilation and distortion. The proposed method with three scale ranges shows nearly identical results of 0.987 in Acc and 0.995 in Sn (0.996 at $= [{1\; 1.5} ]$) in Sp, suggesting its high accuracy and robustness regardless of the scale range selection. It should be noted that even if the OT method has a decent accuracy in this case, the failure on extraction of the complete vessel structures limits its use in complicated cases. The second digital phantom simulation showed similar results (Fig S.2 and Table S1) and validated the accuracy and robustness of the proposed method under different vessel structures.
3.2 In vivo results
To better illustrate the process of the segmentation, the intermediate images of a mouse cortex vasculatures after major steps are presented in Fig. 3. Figure 3(a) shows the original mouse cortex vessel image. The yellow arrows label the representative vessel structures that undergo distortions but are eventually corrected during the process. Figure 3(b) demonstrates the OT result, where evident missing vessels are observed, as indicated by the red arrow. Figure 3(c) shows the HF result. To illustrate the distortion correction effects of the proposed method, a scale range of [1 3] is selected. Under-segmented branching and crossing points, as well as over-segmented vessels can be observed, as labelled by the red arrows. Figure 3(d) is the fused image of the OT and HF results, where the under-segmentations from the HF step are corrected. The vessel disconnectivity caused by HF is labelled by the red arrows in the image. Figure 3(e) shows the image after the URG step, which regains the connectivity, but introduces extra vessel dilation, as indicated by the red arrows. The final result after the FDHM segmentation is shown in Fig. 3(f), where all the distorted structures are corrected, as marked by the green arrows.
The intensity profiles of the original image, the results of the OT, HF and proposed methods along the dashed line in Fig. 3(a) is shown in Fig. 3(g). The small vessel of weak signal in the middle is well extracted by both the HF and proposed methods but missed by the OT method. The large vessel on the left is extracted at the full width at half maximum (FWHM) by the proposed method, but at larger size by the OT and HF methods. Moreover, the vessel on the right is also segmented at FWHM. According to Fig. 3(g), the proposed method applies different thresholds to different vessels, which at best preserves every vessel’s size and shape. In short, the proposed method is capable of providing segmentation results of high accuracy and strong robustness on in vivo data.
The impacts of different scale ranges on the proposed method and the HF method are further evaluated on in vivo vasculature. The segmentations with different scale ranges on the mouse cortex vasculature are shown in Fig. 4. The original vessel image is shown in Fig. 3(a). The results of the proposed method with scale ranges of [1 1.5], [1 3] and [1 10] are shown in Figs. 4(a), 4(c) and 4(e), respectively. The results of the HF method with the three scale ranges are shown in Figs. 4(b), 4(d) and 4(f), respectively. Figure 4(g)-(l) illustrate the corresponding diameter distributions of Fig. 4(a)-(f). According to Figs. 4(g), 4(i) and 4(k), despite the slight differences among the three diameter distributions, the general distribution trends of the proposed method are similar to one another. The diameter distributions of the HF method on the other hand are completely different from one another, demonstrating the heavy impacts of the scale ranges on vessel structures, as shown in Figs. 4(h), 4(j) and 4(l). The vessel sizes under $= [{1\,1.5} ]$ are mostly concentrated at 10-30 µm (Fig. 4(h)), because under-segmentation converts the large vessels into multiple small ones, as marked by the red arrow in Fig. 4(b). With the wider scale range of [1 3], under-segmentation on vessel sizes of 40-60 µm is mitigated (Fig. 4(d)). The scale range of [1 10] dilates all the vessels (Fig. 4(f)), which shifts the entire diameter distribution towards the larger sizes, as shown in Fig. 4(l). The segmentation of mouse ear vessels and corresponding diameter distributions are shown in Fig S.3, where similar results are presented. To be noted, large Gaussian kernel creates or restores vessel continuity [18], which might be the reason for the slight structural differences in Figs. 4(a), 4(c) and 4(e), as indicated by the red arrows. The consistent vessel structures and diameter distributions under different scale ranges in the proposed method ensures the accuracy and robustness of the quantitative results in vivo.
Fig. 4. Vessel segmentation results on mouse cortex vessels of the proposed method and the HF method with different scale ranges, and corresponding vessel diameter distributions. (a) Result of the proposed method, $= [{1\,1.5} ]$; red arrow: structural differences to (c) and (e). (b) Result of the HF method, $= [{1\,1.5} ]$; red arrow: under-segmented vessel structure. (c) Result of the proposed method, $= [{1\,3} ]$; red arrow: structural differences to (a) and (e). (d) Result of the HF method, $= [{1\,3} ]$. (e) Result of the proposed method, $= [{1\,10} ]$; red arrow: structural differences to (a) and (c). (f) Result of the HF method, $= [{1\,3} ]$. (g-l) Corresponding diameter distributions in (a-f). HF, Hessian filter.
The VAFs of the proposed method and HF method on the mouse cortex vasculature are listed in Table 2. The VAFs from the proposed method have differences of ∼0.01 between each two adjacent scale ranges, while those from the HF method show a difference of ∼0.1, around 10 times larger than the proposed method, which further proves the robustness of the proposed method.
Table 2. The VAF of the segmented vessels from the proposed and HF method in different scale ranges
3.3 Longitudinal results
The performance of the proposed method on quantitative analysis of longitudinal applications are validated. Diameter changes on the mouse ear vasculature in response to epinephrine administration are shown in Fig. 5. The original OR-PAM MAP image is shown in Fig. 5(a).
Fig. 5. Quantitative analysis of the vessel diameter changes in response to EP injection. (a) Original OR-PAM MAP image; blue area: target artery segment for diameter analysis. (b)-(d) Segmentation results of the proposed method at BL, +3 min and +15 min; yellow arrow: target artery. (e)-(g) Segmentation results of the HF method at BL, +3 min and +15 min; yellow arrow: target artery. (h) Relative diameter change characterized by the proposed method, the HF method and GT. EP, epinephrine; OR-PAM, optical-resolution photoacoustic microscopy; MAP, maximum amplitude projection; ROI, region of interest; HF, Hessian filter; BL, baseline; GT, ground truth. Error bar: standard deviation.
A subset ROI of vessels is segmented by both the proposed and HF methods with a scale range of [1 5]. The maximum $\sigma $ of 5 is determined by increasing it from 1 till the largest vessel in the image is covered by the result after the HF step. The same scale range is applied to both methods. The diameter mean of an artery, labelled in blue in Fig. 5(a), is calculated at different time points. The segmentation results by the proposed method at baseline (BL), 3 min after injection (+3 min) and 15 min after injection (+15 min) are shown in Fig. 5(b-d). Evident vasoconstriction can be observed at +3 min, as marked by the yellow arrows. The segmentation results by the HF method at BL, +3 min and +15 min are shown in Fig. 5(e-g). The vessel diameter changes characterized by the HF method, the proposed method and the GT are normalized to their baselines and demonstrated in Fig. 5(h). The error bar stands for the standard deviation of the vessel size changes at each pixel along the selected vessel skeleton (centerline for GT). An apparent stronger vasoconstriction of - 50.5% is detected by the proposed method compared to that of - 23.3% in the HF method at +3 min. The vasoconstrictions presented by the proposed and HF methods at +15 min are - 38.8% and - 7.2%, respectively. The GT diameter changes show a vasoconstriction of - 56.3% at +3 min and - 43.1% at +15 min, which is close to the changes detected by the proposed method. The vasoconstrictive response after epinephrine injection followed by a progressive recovery agrees well with previous studies [34]. In this longitudinal study, the proposed method detects a ∼ 2 times vasoconstriction compared to the HF method, and demonstrates a changing trend close to the GT results. The slight difference between the proposed results and GT results could be from different calculation strategies. Vessel radii in the proposed results are calculated by the distance field approach while the GT diameters are directly characterized by the shortest ray length at each pixel along the centerline.
Another longitudinal study of the vessel size changes induced by SNP administration on a mouse ear is implemented. The original OR-PAM image is shown in Fig. 6(a). A subset ROI of vessels is segmented by both the proposed method and the HF method with a scale range of [1 5] in this case. The diameters within a selected artery are averaged in the longitudinal analysis, as labelled in blue in Fig. 6(a). Figure 6(b)-(d) demonstrate the segmentation results by the proposed method at BL, +3 min and +15 min, respectively. Vasodilation can be observed after SNP injection at +3 min, as marked by the yellow arrows. The segmentation results by the HF method at BL, +3 min and +15 min are shown in Fig. 6(e)-(g). The relative diameter changes characterized by the two methods and the GT results are presented in Fig. 6(h). The error bar stands for the standard deviation of the vessel size changes at each pixel along the selected vessel skeleton (centerline for GT).
Fig. 6. Quantitative analysis of the vessel diameter changes in response to SNP injection. (a) Original OR-PAM MAP image; blue area: target artery segment for diameter analysis. (b)-(d) Segmentation results of the proposed method at BL, +3 min and +15 min; yellow arrow: target artery. (e)-(g) Segmentation results of the HF method at BL, +3 min and +15 min; yellow arrow: target artery. (h) Relative diameter change characterized by the proposed, HF methods and GT. SNP, sodium nitroprusside; OR-PAM, optical-resolution photoacoustic microscopy; MAP, maximum amplitude projection; ROI, region of interest; HF, Hessian filter; BL, baseline; GT, ground truth. Error bar: standard deviation.
The vessel diameter presented by the proposed method has an increase of 20.8% at +3 min, and drops down to -2.6% at +15 min. The vessel segmented by the HF method shows a vasodilation of 3.2% at + 3min, and returns to - 0.3% at + 15min. The GT results show a vasodilation of 21.2% at +3 min and 6.6% at +15 min, which is close to the changes detected by the proposed method. In this case, a much stronger vasodilation is detected by the proposed method after the SNP injection, compared to that by the HF method. To be noted, an arterial diameter increase of 3.2% to external stimuli can be negligible as arteries exhibit fluctuations under normal states [36]. Hence, with the results of the HF method with the scale range of [1 5], we cannot conclude whether there is vasodilation in response to SNP.
The results in both longitudinal studies collectively show that the proposed method is more sensitive to subtle vessel changes compared to the HF method. The HF method with a fixed scale range either enlarges the constricted vessels or thins the dilated ones because the Hessian kernels only extract the vessels of matching sizes, which may explain its insensitivity to the subtle vascular changes. On the other hand, segmenting at the half maximum in the proposed method would closely capture the real vascular changes.
4. Conclusion and discussion
In this work, we have developed a Hessian filter-assisted, adaptive thresholding vessel segmentation algorithm. Its performance is quantitively evaluated on a 3D custom-designed digital vessel phantom, demonstrating strong robustness and high accuracy. Further tests on in vivo vascular images show promising consistency in vessel structure preservation regardless of scale range selection. We have also compared the proposed algorithm with the HF method on two longitudinal studies of vascular changes in response to different blood pressure regulating drugs. The results in both studies show much higher sensitivity and accuracy of the proposed algorithm in capturing subtle vessel changes over the HF method.
Different from single-value thresholding methods or adaptive thresholding methods [37–13], the basic idea of the proposed algorithm is to apply different thresholds to each point along the vessel skeleton, in order to better preserve the vascular structures at every point. A fixed threshold applied to different vessels would either overly segment the vessels of larger values or insufficiently segment the vessels of smaller values, both leading to distortions. The strategy proposed here have solved the problem.
The role of the HF step and skeletonization step in the algorithm is to enhance the vessels of weak signal and determine their positions. Without the enhancement, the small vessels can possibly be missed if only grayscale information is applied [16]. Moreover, if two adjacent vessels are accidentally merged by a large scale range, the proposed method can separate them away as long as the half maximum value is able to differentiate them, as demonstrated by the red arrows in Fig S. 3(e) and (f). This also contributes to the algorithm robustness.
The adaptive thresholding following the HF step can better preserve small vessels compared to Otsu’s thresholding. Otsu’s thresholding is a single-value thresholding strategy which automatically sets the threshold that best differentiates all pixels into two groups. The small vessels of weak signals close to background noise are highly likely to be categorized as background, and thus are missed by this method. Meanwhile, the adaptive thresholding applies low thresholds to weaker signals and high thresholds to stronger signals, and therefore can better preserve original vessel structures. Fig S.4 demonstrates the results by the two thresholding strategies after the HF step on the mouse cortex vessel image.
For the HF method, the scale range selection plays a significant role for the quantitative accuracy. Scales ranges need to be adjusted for different scenarios, which cannot provide a general solution for vessel quantitative studies. The resilience to the impacts of different scale ranges in the proposed method allows a general strategy for different applications. The maximum is determined as long as the largest vessel in an image can be segmented in the HF step, and the associated vessel distortions will be corrected automatically in the proposed method.
Because the HF method is widely used in different photoacoustic quantitative studies [39–41], it is employed as a comparison to demonstrate the superiority of the proposed method in this work. However, the HF method is not required for OR-PAM quantitative studies.
The proposed method utilized the Gaussian-shaped distribution of the vessel crossing profiles in MAP images to perform FDHM for 2D vessel segmentation. Because multiple steps of this method can directly be applied in 3D, we believe that it is worthy of exploring its potentials in 3D vessel segmentation in our future work.
During the process of the proposed algorithm, several parameters are selected empirically. In the HFE step, the parameter a controlling the offset is set to 2 and the parameter b controlling the high-frequency emphasis is set to 3. The clipping limit in the CLAHE step is set to 0.005. The sizes of the two Gaussian filters in the adaptive thresholding step are set to 3 and 11 pixels. These parameters together with the scale range collectively determine the extent of vessel enhancement and the manual selection of them is barely evitable in the filter-based method. Therefore, fully automatic segmentation strategies are in particular needs. Novel imaging processing techniques based on neural networks have shown promising potentials [14,15]. We expect that the development of such techniques may further push the limits on fully automatic vessel segmentation.
Funding
National Natural Science Foundation of China (61871251, 62027901, 82071925).
Disclosures
The authors declare that there are no conflicts of interest related to this article.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
Supplemental document
See Supplement 1 for supporting content.
References
1. L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science 335(6075), 1458–1462 (2012). [CrossRef]
2. R. Cao, J. Li, B. Ning, N. Sun, T. Wang, Z. Zuo, and S. Hu, “Functional and oxygen-metabolic photoacoustic microscopy of the awake mouse brain,” NeuroImage 150, 77–87 (2017). [CrossRef]
3. A. Khadria, C. D. Paavola, K. Maslov, F. A. Valenzuela, A. E. Sperry, A. L. Cox, R. Cao, J. Shi, P. L. Brown-Augsburger, E. Lozano, R. L. Blankenship, R. Majumdar, S. A. Bradley, J. M. Beals, S. S. Oladipupo, and L. V. Wang, “Photoacoustic imaging reveals mechanisms of rapid-acting insulin formulations dynamics at the injection site,” Mol. Metab. 62, 101522 (2022). [CrossRef]
4. J. Yao, L. Wang, J.-M. Yang, K. I. Maslov, T. T. W. Wong, L. Li, C.-H. Huang, J. Zou, and L. V. Wang, “High-speed label-free functional photoacoustic microscopy of mouse brain in action,” Nat. Methods 12(5), 407–410 (2015). [CrossRef]
5. J. Yao and L. V. Wang, “Photoacoustic microscopy,” Laser Photonics Rev. 7(5), 758–778 (2013). [CrossRef]
6. Q. Chen, T. Jin, W. Qi, X. Mo, and L. Xi, “Label-free photoacoustic imaging of the cardio-cerebrovascular development in the embryonic zebrafish,” Biomed. Opt. Express 8(4), 2359–2367 (2017). [CrossRef]
7. J.-M. Yang, C. Li, R. Chen, B. Rao, J. Yao, C.-H. Yeh, A. Danielli, K. Maslov, Q. Zhou, K. K. Shung, and L. V. Wang, “Optical-resolution photoacoustic endomicroscopy in vivo,” Biomed. Opt. Express 6(3), 918–932 (2015). [CrossRef]
8. M. C. Zanon Zotin, L. Sveikata, A. Viswanathan, and P. Yilmaz, “Cerebral small vessel disease and vascular cognitive impairment: from diagnosis to management,” Curr. Opin. Neurol. 34(2), 246–257 (2021). [CrossRef]
9. C. L. Jarrett, K. L. Shields, R. M. Broxterman, J. R. Hydren, S. H. Park, J. R. Gifford, and R. S. Richardson, “Imaging transcranial Doppler ultrasound to measure middle cerebral artery blood flow: the importance of measuring vessel diameter,” Am. J. Physiol. Regul. Integr. Comp. Physiol. 319(1), R33–R42 (2020). [CrossRef]
10. B. A. Corliss, C. Mathews, R. Doty, G. Rohde, and S. M. Peirce, “Methods to label, image, and analyze the complex structural architectures of microvascular networks,” Microcirculation 26(5), e12520 (2019). [CrossRef]
11. N. Otsu, “A threshold selection method from gray-level histograms,” IEEE Trans. Syst. Man Cybern. 9(1), 62–66 (1979). [CrossRef]
12. Q. B. Truong and B. R. Lee, “Automatic multi-thresholds selection for image segmentation based on evolutionary approach,” Int. J. Control Autom. Syst. 11(4), 834–844 (2013). [CrossRef]
13. J. Xing, P. Yang, and L. Qingge, “Robust 2D otsu’s algorithm for uneven illumination image segmentation,” Comput. Intell. Neurosci. 2020(3), 1–14 (2020). [CrossRef]
14. N.-K. Chlis, A. Karlas, N.-A. Fasoula, M. Kallmayer, H.-H. Eckstein, F. J. Theis, V. Ntziachristos, and C. Marr, “A sparse deep learning approach for automatic segmentation of human vasculature in multispectral optoacoustic tomography,” Photoacoustics 20, 100203 (2020). [CrossRef]
15. A. Y. Yuan, Y. Gao, L. Peng, L. Zhou, J. Liu, S. Zhu, and W. Song, “Hybrid deep learning network for vascular segmentation in photoacoustic imaging,” Biomed. Opt. Express 11(11), 6445–6457 (2020). [CrossRef]
16. Q. Li, L. Li, T. Yu, Q. Zhao, C. Zhou, and X. Chai, “Vascular tree extraction for photoacoustic microscopy and imaging of cat primary visual cortex,” J. Biophotonics 10(6-7), 780–791 (2017). [CrossRef]
17. A. F. Frangi, W. J. Niessen, K. L. Vincken, and M. A. Viergever, “Multiscale vessel enhancement filtering,” in MICCAI’98 (Springer Berlin Heidelberg, Berlin, Heidelberg, 1998), pp. 130–137.
18. A. Longo, S. Morscher, J. M. Najafababdi, D. Jüstel, C. Zakian, and V. Ntziachristos, “Assessment of hessian-based Frangi vesselness filter in optoacoustic imaging,” Photoacoustics 20, 100200 (2020). [CrossRef]
19. Z. Yang, J. Chen, J. Yao, R. Lin, J. Meng, C. Liu, J. Yang, X. Li, L. Wang, and L. Song, “Multi-parametric quantitative microvascular imaging with optical-resolution photoacoustic microscopy in vivo,” Opt. Express 22(2), 1500 (2014). [CrossRef]
20. Q. Zhao, R. Lin, C. Liu, J. Zhao, G. Si, L. Song, and J. Meng, “Quantitative analysis on in vivo tumor-microvascular images from optical-resolution photoacoustic microscopy,” J. Biophotonics 12, e201800421 (2019).
21. K. BahadarKhan, A. A Khaliq, and M. Shahid, “A morphological Hessian based approach for retinal blood vessels segmentation and denoising using region based Otsu thresholding,” PLoS One 11(7), e0158996 (2016). [CrossRef]
22. M. J. Alwazzan, M. A. Ismael, and A. N. Ahmed, “A hybrid algorithm to enhance colour retinal fundus images using a Wiener filter and CLAHE,” J. Digit. Imaging 34(3), 750–759 (2021). [CrossRef]
23. R. G. Koch, A. Tsamis, A. D’Amore, W. R. Wagner, S. C. Watkins, T. G. Gleason, and D. A. Vorp, “A custom image-based analysis tool for quantifying elastin and collagen micro-architecture in the wall of the human aorta from multi-photon microscopy,” J. Biomech. 47(5), 935–943 (2014). [CrossRef]
24. P. S. Tsai, J. P. Kaufhold, P. Blinder, B. Friedman, P. J. Drew, H. J. Karten, P. D. Lyden, and D. Kleinfeld, “Correlations of neuronal and microvascular densities in murine cortex revealed by direct counting and colocalization of nuclei and vessels,” J. Neurosci. 29(46), 14553–14570 (2009). [CrossRef]
25. J. K. Varma, K. Subramanyan, and J. Durgan, “Full width at half maximum as a measure of vessel diameter in computed tomography angiography,” in Proc.SPIE (2004).
26. B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. 15(2), 021314 (2010). [CrossRef]
27. D. Zhang, R. Li, M. Chen, T. Vu, H. Sheng, W. Yang, U. Hoffmann, J. Luo, and J. Yao, “Photoacoustic imaging of in vivo hemodynamic responses to sodium nitroprusside,” J. Biophotonics 14(7), e202000478 (2021). [CrossRef]
28. X. Song, G. Chen, A. Zhao, X. Liu, and J. Zeng, “Virtual optical-resolution photoacoustic microscopy using the k-Wave method,” Appl. Opt. 60(36), 11241–11246 (2021). [CrossRef]
29. S. Hu, K. Maslov, and L. V. Wang, “Second-generation optical-resolution photoacoustic microscopy with improved sensitivity and speed,” Opt. Lett. 36(7), 1134–1136 (2011). [CrossRef]
30. “Institutional Animal Care and Use Committee Guidebook,” Office of laboratory animal welfare (OLAW), National institutes of health (2002).
31. J. A. Friederich and J. F. Butterworth, “Sodium nitroprusside: twenty years and counting,” Anesth. Analg. 81, 152–162 (1995).
32. R. Nosrati, S. Lin, R. Mohindra, A. Ramadeen, V. Toronov, and P. Dorian, “Study of the effects of epinephrine on cerebral oxygenation and metabolism during cardiac arrest and resuscitation by hyperspectral near-infrared spectroscopy,” Crit. Care Med. 47(4), e349–e357 (2019). [CrossRef]
33. C. D. Mavroudis, T. S. Ko, R. W. Morgan, L. E. Volk, W. P. Landis, B. Smood, R. Xiao, M. Hefti, T. W. Boorady, A. Marquez, M. Karlsson, D. J. Licht, V. M. Nadkarni, R. A. Berg, R. M. Sutton, and T. J. Kilbaugh, “Epinephrine’s effects on cerebrovascular and systemic hemodynamics during cardiopulmonary resuscitation,” Critical Care 24(1), 583 (2020). [CrossRef]
34. B. Lan, W. Liu, Y.-C. Wang, J. Shi, Y. Li, S. Xu, H. Sheng, Q. Zhou, J. Zou, U. Hoffmann, W. Yang, and J. Yao, “High-speed widefield photoacoustic microscopy of small-animal hemodynamics,” Biomed. Opt. Express 9(10), 4689–4701 (2018). [CrossRef]
35. S. Avadiappan, S. Payabvash, M. A. Morrison, A. Jakary, C. P. Hess, and J. M. Lupo, “A fully automated method for segmenting arteries and quantifying vessel radii on magnetic resonance angiography images of varying projection thickness,” Front. Neurosci. 14, 537 (2020). [CrossRef]
36. H. C. Chen, V. Patel, J. Wiek, S. M. Rassam, and E. M. Kohner, “Vessel diameter changes during the cardiac cycle,” Eye 8(1), 97–103 (1994). [CrossRef]
37. J. Brieva, E. Gonzalez, F. Gonzalez, A. Bousse, and J. J. Bellanger, “A level set method for vessel segmentation in coronary angiography,” in Conf Proc IEEE Eng Med Biol Soc (2005), pp. 6348–6351.
38. H. Yazid, H. Arof, and H. Mohd Isa, “Exudates segmentation using inverse surface adaptive thresholding,” Measurement 45(6), 1599–1608 (2012). [CrossRef]
39. H. Zhao, G. Wang, R. Lin, X. Gong, L. Song, T. Li, W. Wang, K. Zhang, X. Qian, H. Zhang, L. Li, Z. Liu, and C. Liu, “Three-dimensional Hessian matrix-based quantitative vascular imaging of rat iris with optical-resolution photoacoustic microscopy in vivo,” J. Biomed. Opt. 23(04), 1 (2018). [CrossRef]
40. T. T. Mai, S. W. Yoo, S. Park, J. Y. Kim, K.-H. Choi, C. Kim, S. Y. Kwon, J.-J. Min, and C. Lee, “In Vivo Quantitative Vasculature Segmentation and Assessment for Photodynamic Therapy Process Monitoring Using Photoacoustic Microscopy,” Sensors 21(18), C1 (2021). [CrossRef]
41. L. Lin, P. Hu, J. Shi, C. M. Appleton, K. Maslov, L. Li, R. Zhang, and L. V. Wang, “Single-breath-hold photoacoustic computed tomography of the breast,” Nat. Commun. 9(1), 2352 (2018). [CrossRef]
