Many diseases involve either the formation of new blood vessels (e.g., tumor angiogenesis) or the damage of existing ones (e.g., diabetic retinopathy) at the microcirculation level. Optical-resolution photoacoustic microscopy (OR-PAM), capable of imaging microvessels in 3D in vivo down to individual capillaries using endogenous contrast, has the potential to reveal microvascular information critical to the diagnosis and staging of microcirculation-related diseases. In this study, we have developed a dedicated microvascular quantification (MQ) algorithm for OR-PAM to automatically quantify multiple microvascular morphological parameters in parallel, including the vessel diameter distribution, the microvessel density, the vascular tortuosity, and the fractal dimension. The algorithm has been tested on in vivo OR-PAM images of a healthy mouse, demonstrating high accuracy for microvascular segmentation and quantification. The developed MQ algorithm for OR-PAM may greatly facilitate quantitative imaging of tumor angiogenesis and many other microcirculation related diseases in vivo.
© 2014 Optical Society of America
Photoacoustic tomography (PAT) is a rapidly developing biomedical imaging technology that can provide anatomic, functional, and molecular contrasts of intact biological tissue in vivo, at multiple scales from organelles to organs . In PAT, a nanosecond pulsed laser is usually used to illuminate the biological tissue, which generates a localized transient thermoelastic expansion and thus leads to the emission of wideband ultrasonic waves (also termed as photoacoustic waves). The emitted photoacoustic waves are then detected to reconstruct the optical absorption properties of the biological tissue, which correlate with many important physiological parameters, including the total concentration, the oxygen saturation, and the oxygen metabolism of hemoglobin. Optical-resolution photoacoustic microscopy (OR-PAM) is a specific form of PAT that offers optical-diffraction limited transverse spatial resolution, which can be as fine as micrometers or even sub-micrometers [2–5]. OR-PAM uses a tightly focused laser beam for photoacoustic excitation and a focused single-element ultrasonic transducer to record depth-resolved 1D photoacoustic images (A-lines). By simultaneously scanning the laser beam and the ultrasonic transducer in 2D, volumetric (3D) OR-PAM images can be acquired. Using endogenous contrast from hemoglobin, OR-PAM can image microvascular morphology and functions in vivo at high resolution. As a result, it is becoming a powerful tool for studying the microcirculation in many physiological and pathological conditions . For example, recently, OR-PAM has been applied to the longitudinal study of angiogenesis (the imaging of angiogenesis at a series of time points over a period of time—days or weeks—to monitoring the evolution of angiogenesis) [7, 8], which plays a critical role in tumor growth, invasion, and metastasis [9, 10]. Such studies may open up new opportunities to better understand the dynamics of the tumor microenvironment, as well as to assess the efficacy of anti-angiogenic and/or combined tumor therapies in the early stages [11–13].
Four parameters—the vessel diameter distribution, the microvessel density, the vascular tortuosity, and the fractal dimension—have been widely accepted to describe and assess the vascular morphology. Quantitative imaging of these parameters in vivo at different time points can provide us insights into the dynamics of tumor angiogenesis [13–18]. In previous studies, we used a simple thresholding approach, together with cross-section vessel tracking, to segment the microvasculature in OR-PAM images [7, 8]. However, to offer truly accurate quantification for the intricate microvascular images acquired by OR-PAM, a more rigorous vascular segmentation and quantification algorithm is needed. Currently, conventional vascular analysis algorithms developed for medical imaging systems, such as X-ray digital subtraction angiography (DSA) and magnetic resonance angiography (MRA), generally lack the ability to accurately segment and quantify images of capillary-level microvessels, which are often falsely treated as noise with such algorithms [19, 20]. In addition, the unique multi-scale feature of OR-PAM (imaging of vascular trees up to 6 – 7 orders with vessel diameters ranging from >200 µm to <10 µm) also presents challenges for these algorithms.
In this study, we have developed a dedicated microvascular quantification (MQ) algorithm for OR-PAM to automatically quantify multiple microvascular parameters in parallel, including the vessel diameter distribution, the microvessel density, the vascular tortuosity, and the fractal dimension. The algorithm has been validated with in vivo OR-PAM images of a healthy mouse, demonstrating its ability to accurately segment and quantify the vasculature—including capillary-level microvessels—in OR-PAM images. The results shown in the current study suggest that the developed MQ algorithm may greatly facilitate the application of OR-PAM for quantitative imaging of tumor angiogenesis and many other microcirculation related diseases.
2.1 OR-PAM system and image acquisition
The schematic of our OR-PAM system is illustrated in Fig. 1. For photoacoustic excitation, a pulsed laser beam (pulse width: 1.8 ns) at 532 nm from an Nd:YAG laser source (SPOT-532, Elforlight) was focused to an optical diffraction-limited spot to irradiate the sample. The generated photoacoustic waves were detected by a 75-MHz center frequency ultrasonic transducer via a custom-made optical-acoustic beam combiner. The output electrical signals from the ultrasonic transducer were amplified using a low-noise amplifier (ZFL-500LN-BNC, Mini-Circuits), digitized via a 200-MS/s data acquisition (DAQ) card (CS1422, GaGe), and then stored into a personal computer (PC). For 3D imaging, a 2D motorized translational stage was used to scan the OR-PAM imaging head (Fig. 1). Further details of the OR-PAM system can be found in our previous publication . For in vivo imaging, an anesthetized female BALB/c mouse (around six-week old and weighing ~20 g) was used. OR-PAM images of both the mouse ear and back were acquired. In all in vivo experiments, the laser energy reaching the biological tissue was maintained to be ~100 nJ/pulse, corresponding to an optical fluence of ~19 mJ/cm2 on the tissue surface (when focused ~200 μm below the skin surface), which conforms to the 20 mJ/cm2 ANSI safety standards (the maximum permissible exposure). All experimental animal procedures were carried out in compliance with protocols approved by the Animal Studies Committee of the Shenzhen Institutes of Advanced Technology, the Chinese Academy of Sciences.
2.2 Overall algorithm flow chart
The overall flowchart of our microvascular quantification (MQ) algorithm is shown in Fig. 2.First, blood vessels are extracted (segmented) from the original image, based on a modified Hessian matrix method . Then, the vascular centerlines corresponding to the extracted blood vessels are identified. Finally, using the segmented image and the identified centerlines, the vascular morphological parameters, namely, the vessel diameter distribution, the microvessel density, the vascular tortuosity, and the fractal dimension, are computed. All post-processing of the acquired images in this study were carried out using MATLAB (R2012b, Mathworks) on a PC with an Intel(R) CoreTM 2 Duo CPU E7500@2.93 GHz and a 4-GB RAM.
2.3 Vessel extraction
Step1, feature map calculation. First, before computing the Hessian matrix of the image, high-frequency emphasis (HFE) filtering was used to enhance the contrast of the microvessels in the image. The definition of the HFE filter is as follows:
Second, the Hessian matrix of image I(x,y) under scale s—which essentially extracts the local curvature of the structure—was computed as
Third, the eigenvalues of the Hessian matrix were computed, which represent the curvature along the principal directions corresponding to respective eigenvectors. Since blood vessels in OR-PAM images are shown as bright tubular structures, as discussed in , one of the eigenvalues represents the curvature along the direction of a vessel, which should be close to zero, while the other represents the curvature along the perpendicular direction of the vessel. Mathematically, the eigenvalues of the Hessian matrix representing blood vessels should satisfy the condition of and (let ), where is a small fraction number. More comprehensive descriptions on the relationship between the Hessian matrix eigenvalues and the representation of various structures can be found in .
To further enhance the microvessels and suppress noise, we define the feature map f under scale s as:
Finally, the feature map of image I can be obtained as:
Step2, feature map enhancement by fractional differential. Although most blood vessels can be revealed in the feature map after step 1, capillary regions may still be dim, suggesting that the texture of the feature map needs to be further enhanced. Thus, a fractional differential operator was employed to modulate the horizontal and vertical gradient fields of the feature map, as it could be adjusted appropriately to boost the high-frequency components while retaining the key low-frequency components. Then, using the enhanced gradient fields, a new feature map f’ was obtained through gradient domain reconstruction .
Step3, intensity transformation. To further improve the contrast between the blood vessels and the background, an intensity transformation function was applied to the enhanced feature map f’. The transformed image g is given by:
Step4, region growing. The pixels with the maximum intensity value (which equals to 1) in the ultimately enhanced images were selected as the seed points for region growing, through which the blood vessels were extracted as binary images.
To better illustrate the image processing steps described above, two 2D OR-PAM image slices were selected to show the intermediate results after each major step (Fig. 3). The feature maps of the original images are shown in Figs. 3(a) and 3(d), while the fractional-differential enhanced feature maps are given in Figs. 3(b) and 3(e). Finally, the segmented blood vessel images based on region growing are shown in Figs. 3(c) and 3(f).
2.4 Vascular centerlines
Based on the extracted binary blood vessel images, an augmented fast marching method (AFMM) was employed and adapted to obtain the centerlines of the blood vessels . Briefly, a 3D AFMM method was used to solve the equation describing the evolution of a closed contour surface as a function of time, with a specific speed function in the normal direction at each point on the surface. Specifically, in this study, the centerlines were determined using a two-step AFMM approach. The flowchart of the method is given in Fig. 4, while its major implementation steps are described in detail below.
- Step 1: Separate the vascular networks into individual sub-images, so that each of which includes only one set of interconnected vascular tree; initialize the speed of points inside and outside the vessels as constants 1 and 0, respectively
- Step 2 (First-step AFMM): For one sub-image, the points inside the vessels are used as seeds to expand uniformly at the initial speed, until their contour surfaces reach the closest boundary of the vessel; calculate the distance distribution from every point inside the vessel to its corresponding closest vessel boundary
- Step 3: Find the point with the largest distance as the “global maximum distance point” (usually the center point of a major intersection in a vascular tree); set the distance value of each vessel point as the updated speed of this point for Step 4
- Step 4 (Second-step AFMM): Start the expansion from the global maximum distance point and stops when the expanding contour surfaces reach the entire vessel boundary; similarly, calculate another distance distribution from all vessel points to the global maximum distance point, termed as the flight time map
- Step 5: Set the point of the maximum flight time as the furthest point; calculate the gradient of the flight time map
- Step 6: To extract the centerline of a branch in the vascular tree, set the furthest point as the starting point and then connect along the points in the opposite direction of the fastest gradient descent of the flight time map, termed as the back-tracking method
- Step 7: Repeat Step 6 to extract the centerlines of all branches in one vascular tree to form the skeleton of this vascular tree
- Step 8: Iterate Step 1 to Step 7 until the centerlines of all vascular trees (in all sub-images) are extracted to form the vessel network
To illustrate the concept of vascular centerlines more intuitively, Fig. 5shows a composite image consisting of a segmented 3D OR-PAM vascular image overlaid with the extracted vascular centerlines computed using the algorithm described above. From Fig. 5, it can also be seen that, in general, the computed centerlines are accurately positioned within the vessels.
2.5 Vessel diameters
As the tangential direction at each point on the centerline is perpendicular to the cross-section of the vessel, the diameter of the vessel at a specified point can be calculated as the distance between the two intersection points between the vessel edges (determined in 2.3) and the cross-section line.
2.6 Microvessel density
In this study, the microvessel density (MVD) is defined as the length of the vessel per unit volume, as given below :
For 3D vasculature, there are three widely accepted definitions for tortuosity, namely, the distance metric (DM), the inflection count metric (ICM), and the sum of angles metric (SOAM) [18, 25, 26]. In short, DM is defined as the ratio between the actual path length of a vessel and the linear distance between the end points of this vessel; ICM is defined as the production between the number of a curve’s inflection points and the DM of this curve; SOAM is defined as the sum of the total curvature along a curve normalized by the curve’s path length. In this study, the vascular tortuosity was computed in all three definitions.
2.8 Fractal dimension
The fractal dimension can also be calculated in different methods, including the Hausdorff dimension, the box-counting dimension (or box dimension), and the sand box dimension, among which the box dimension is the most commonly used definition and also the one used in this study [14, 27].
To calculate the box dimension, grids with a size of rrr were first used to cover the extracted vascular image; then, the minimum number of grids occupied by the blood vessels was counted (denoted by N). Upon gradually reducing r (size of the grids), N will increase accordingly. The fractal dimension D was calculated according to Eq. (7) below :
First, we compared the vessel extraction results using our algorithm with those from algorithms developed for conventional medical imaging systems. Because one algorithm (the level set method) was quite time-consuming, only a 2D OR-PAM image slice (acquired from a mouse ear in vivo) consisting of 1600 × 1200 pixels was selected for this comparison. From Fig. 6, it can be seen that, compared with the traditional Hessian matrix algorithms developed for DSA or MRA (Figs. 6(b) and 6(c)), our algorithm has produced better results for microvessel extraction (Fig. 6(d)), including: (1) better extraction of microvessels in the capillary beds regions; (2) better continuity for the extracted microvessels; (3) better intensity uniformity and continuity in individual vessels.
Further, areas 1 and 2 in the dashed boxes in Fig. 6(a) were enlarged (Figs. 7(a) and 7(h)) for further analysis and validation. The processed images of these sub-areas using the level set method (Figs. 7(b) and 7(i)), the traditional Hessian matrix method (Figs. 7(c) and 7(j)), and our method (Figs. 7(d) and 7(k)) were compared side by side. In Fig. 7, from the comparison of the two representative cross-section intensity profiles, it can be clearly seen that, using our method, the majority of the blood vessels from the original images are successfully extracted, and the continuity and intensity uniformity of the extracted vessels are in general superior to the other two methods.
Further, our vessel extraction algorithm was validated using a volumetric OR-PAM image (consisting of 1600 × 1600 × 18 pixels) of a mouse ear in vivo. Figure 8(a) is the original volumetric photoacoustic image with depth encoded in color; Figs. 8(b) and 8(c) are the images with vessels extracted using the simple thresholding algorithm in our previous work  and the algorithm developed in this study, respectively.
Finally, to validate the consistency of the performance of our algorithm, the vasculature in three in vivo OR-PAM images acquired at different regions of an anesthetized mouse were segmented and quantified. Figure 9(a) shows an OR-PAM image acquired from the root region of a mouse ear, which is characterized by relatively large vessels with nicely ordered vascular trees. Figure 9(b) shows an OR-PAM image acquired from a capillary bed region of the same mouse ear. This region, while still shows multiple orders of vascular branches, is essentially dominated by densely packed microvessels. Figure 9(c) shows an OR-PAM image acquired from the superficial dorsal region of the same mouse. The image of this region is characterized by relatively uniform vascular distribution. The calculated vessel diameter distribution using our MQ algorithm is shown in Fig. 9(d), while other parameters are given in Fig. 9(e). It can be seen that, the computed parameters agree well with some visually observed characteristics, for example, Fig. 9(b) has a significantly higher microvessel density due to the densely packed microvessels; meanwhile, they also reveal some information that cannot be easily captured by direct visual observations, for example, Fig. 9(a) has a relatively large tortuosity in all three definitions (DM, ICM, and SOAM). Note that the images in Fig. 9 were acquired by adjusting the focusing depth of the OR-PAM imaging head, in order to image across an extended depth of the various regions of the back of a mouse.
4. Conclusions and discussion
In this study, we developed a Hessian matrix based microvascular quantification algorithm dedicated for optical-resolution photoacoustic microscopy (OR-PAM). The Hessian matrix feature map and subsequent image processing steps were specifically constructed or designed to cater the need for segmenting the intricate microvasculature imaged with OR-PAM. To provide a foundation for computing the vascular morphological parameters, we developed a two-step augmented fast marching method to obtain the vascular centerlines. Using the segmented images and calculated centerlines, four vascular morphological parameters, including the vessel diameter distribution, the microvessel density, the vascular tortuosity, and the fractal dimension, were computed. The developed microvascular quantification (MQ) algorithm has demonstrated great consistency and reliability when tested on multiple in vivo volumetric OR-PAM images acquired from a healthy mouse. The results in this study suggest that, the developed MQ algorithm, in conjunction with OR-PAM or other microvascular imaging technologies such as Doppler OCT, can potentially facilitate the quantitative imaging and study of tumor angiogenesis and many other microcirculation related diseases in vivo.
This work was supported in part by the National Natural Science Foundation of China grant: 61205203; the National Key Basic Research (973) Program of China: 2014CB744503; the Shenzhen Science and Technology Innovation Committee grants: ZDSY20130401165820357, KQCX20120816155844962, CXZZ20120617113635699, and JCYJ20120615125857842.
References and links
5. P. Hajireza, W. Shi, and R. Zemp, “Label-free in vivo GRIN-lens optical resolution photoacoustic micro-endoscopy,” Laser Phys. Lett. 10(5), 055603 (2013). [CrossRef]
7. S. S. Oladipupo, S. Hu, A. C. Santeford, J. Yao, J. R. Kovalski, R. V. Shohet, K. Maslov, L. V. Wang, and J. M. Arbeit, “Conditional HIF-1 induction produces multistage neovascularization with stage-specific sensitivity to VEGFR inhibitors and myeloid cell independence,” Blood 117(15), 4142–4153 (2011). [CrossRef]
8. S. Oladipupo, S. Hu, J. Kovalski, J. J. Yao, A. Santeford, R. E. Sohn, R. Shohet, K. Maslov, L. V. Wang, and J. M. Arbeit, “VEGF is essential for hypoxia-inducible factor-mediated neovascularization but dispensable for endothelial sprouting,” Proc. Natl. Acad. Sci. U. S. A. 108(32), 13264–13269 (2011). [CrossRef]
10. J. Folkman, M. Bach, J. W. Rowe, F. Davidoff, P. Lambert, C. Hirsch, A. Goldberg, H. H. Hiatt, J. Glass, and E. Henshaw, “Tumor angiogenesis: therapeutic implications,” N. Engl. J. Med. 285(21), 1182–1186 (1971). [CrossRef]
13. B. J. Vakoc, R. M. Lanning, J. A. Tyrrell, T. P. Padera, L. A. Bartlett, T. Stylianopoulos, L. L. Munn, G. J. Tearney, D. Fukumura, R. K. Jain, and B. E. Bouma, “Three-dimensional microscopy of the tumor microenvironment in vivo using optical frequency domain imaging,” Nat. Med. 15(10), 1219–1223 (2009). [CrossRef]
14. Y. Gazit, J. W. Baish, N. Safabakhsh, M. Leunig, L. T. Baxter, and R. K. Jain, “Fractal characteristics of tumor vascular architecture during tumor growth and regression,” Microcirculation 4(4), 395–402 (1997). [CrossRef]
15. L. Hlatky, P. Hahnfeldt, and J. Folkman, “Clinical application of antiangiogenic therapy: Microvessel density, what it does and doesn’t tell us,” J. Natl. Cancer Inst. 94(12), 883–893 (2002). [CrossRef]
16. E. Bullitt, D. L. Zeng, G. Gerig, S. Aylward, S. Joshi, J. K. Smith, W. L. Lin, and M. G. Ewend, “Vessel tortuosity and brain tumor malignancy: A blinded study,” Acad. Radiol. 12(10), 1232–1240 (2005). [CrossRef]
17. H. Harada, X. J. Xie, S. Itasaka, L. H. Zeng, Y. X. Zhu, A. Morinibu, K. Shinomiya, and M. Hiraoka, “Diameter of tumor blood vessels is a good parameter to estimate HIF-1-active regions in solid tumors,” Biochem. Biophys. Res. Commun. 373(4), 533–538 (2008). [CrossRef]
18. E. Bullitt, G. Gerig, S. M. Pizer, W. L. Lin, and S. R. Aylward, “Measuring tortuosity of the intracerebral vasculature from MRA images,” IEEE Trans. Med. Imaging 22(9), 1163–1171 (2003). [CrossRef]
19. A. F. Frangi, W. J. Niessen, K. L. Vincken, and M. A. Viergever, “Multiscale vessel enhancement filtering,” Med. Image Comput. Comput. Assisted Intervention 1496, 130–137 (1998).
20. C. Li, R. Huang, Z. Ding, J. Gatenby, D. N. Metaxas, and J. C. Gore, “A level set method for image segmentation in the presence of intensity inhomogeneities with application to MRI,” IEEE Trans. Image Process. 20(7), 2007–2016 (2011). [CrossRef]
22. R. Fattal, D. Lischinski, and M. Werman, “Gradient domain high dynamic range compression,” ACM Trans. Graphics 21, 249–256 (2002).
24. R. K. Jain, N. Safabakhsh, A. Sckell, Y. Chen, P. Jiang, L. Benjamin, F. Yuan, and E. Keshet, “Endothelial cell death, angiogenesis, and microvascular function after castration in an androgen-dependent tumor: role of vascular endothelial growth factor,” Proc. Natl. Acad. Sci. U. S. A. 95(18), 10820–10825 (1998). [CrossRef]
25. E. Bullitt, D. Zeng, G. Gerig, S. Aylward, S. Joshi, J. K. Smith, W. Lin, and M. G. Ewend, “Vessel tortuosity and brain tumor malignancy: a blinded study,” Acad. Radiol. 12(10), 1232–1240 (2005). [CrossRef]
26. A. H. Parikh, J. K. Smith, M. G. Ewend, and E. Bullitt, “Correlation of MR perfusion imaging and vessel tortuosity parameters in assessment of intracranial neoplasms,” Technol. Cancer Res. Treat. 3(6), 585–590 (2004). [PubMed]
27. J. W. Baish and R. K. Jain, “Fractals and cancer,” Cancer Res. 60(14), 3683–3688 (2000). [PubMed]
28. S. P. Lalley and D. Gatzouras, “Hausdorff and box dimensions of certain self-affine fractals,” Indiana Univ. Math. J. 41(2), 533–568 (1992). [CrossRef]