Open Access
23 April 2018 Three-dimensional Hessian matrix-based quantitative vascular imaging of rat iris with optical-resolution photoacoustic microscopy in vivo
Huangxuan Zhao, Guangsong Wang, Riqiang Lin, Xiaojing Gong, Liang Song, Tan Li, Wenjia Wang, Kunya Zhang, Xiuqing Qian, Haixia Zhang, Lin Li, Zhicheng Liu, Chengbo Liu
Author Affiliations +
For the diagnosis and evaluation of ophthalmic diseases, imaging and quantitative characterization of vasculature in the iris are very important. The recently developed photoacoustic imaging, which is ultrasensitive in imaging endogenous hemoglobin molecules, provides a highly efficient label-free method for imaging blood vasculature in the iris. However, the development of advanced vascular quantification algorithms is still needed to enable accurate characterization of the underlying vasculature. We have developed a vascular information quantification algorithm by adopting a three-dimensional (3-D) Hessian matrix and applied for processing iris vasculature images obtained with a custom-built optical-resolution photoacoustic imaging system (OR-PAM). For the first time, we demonstrate in vivo 3-D vascular structures of a rat iris with a the label-free imaging method and also accurately extract quantitative vascular information, such as vessel diameter, vascular density, and vascular tortuosity. Our results indicate that the developed algorithm is capable of quantifying the vasculature in the 3-D photoacoustic images of the iris in-vivo, thus enhancing the diagnostic capability of the OR-PAM system for vascular-related ophthalmic diseases in vivo.



More than 67 million people are diagnosed with glaucoma globally, making it the most commonly seen irreversible ocular disease worldwide.14 Primary angle closure glaucoma (PACG) is a major subtype of the disease and affects up to 16 million patients. Earlier studies have shown that the morphology alteration of the iris, such as rich blood vasculature, vascular rearrangement and hemodynamics changes, occurs with the development of PACG.59 Thus, iris vascular imaging and characterization of vasculature is of great significance to provide valuable information not only for the diagnosis of PACG, but also for understanding the role of blood vasculature alteration in the progression of the disease.

Fluorescein angiography (FA) using exogenously injected contrast agents is currently used in clinics for iris vascular imaging.1012 The primary goal of this method is to reveal the alterations in the morphology of the diseased iris blood vessels and also to assess the hemodynamic changes by evaluating the degree of vascular leakage. Although FA has been widely applied to diagnose PACG and many other eye diseases, such as traumatic iris changes, iris tumor, and iris dysplasia, it is still inadequate for imaging small blood vessels (microvessels) since exogenous contrast agents can hardly reach some of the microvessels due to the leakage of the contrast agents in the artery, the slow blood flow in some of the microvessels, and iris lesions.13,14 In addition, the imaging quality of FA is suboptimal due to the limited resolution and leakage of the contrast agents and is invasive due to the use of exogenous contrast agents, which may pose safety concerns to the patients since the injected contrast agents could potentially cause undesirable complications.14,15 Optical coherence tomography (OCT) is another imaging modality that has been used for ophthalmic vascular imaging,1618 where blood vasculature information and ocular tissue constituents are acquired simultaneously. However, this technology is mostly limited to the posterior segment of the eye, such as imaging retina vessels, and its application to iris imaging is limited due to the slow blood flow in the relatively small blood vessels of the iris and the high background signal interference of the iris pigments.

Imaging of small microvessels is invaluable to diagnose PACG,19 but the existing solutions based upon OCT and FA are suboptimal. The photoacoustic imaging technology, which images small microblood vessels with high fidelity,20,21 is highly suitable for the iris imaging and thus for PACG evaluation. A nanosecond pulsed laser is generally applied in photoacoustic imaging to shed laser into the biological tissues. The absorbed laser energy by the chromophores (such as oxy- and deoxyhemoglobin) in the tissue generates acoustic waves through the process of laser energy absorption, instantaneous temperature rise, transient thermoselastic expansion, and acoustic wave emission. These generated acoustic waves (also termed as photoacoustic signals) are captured by the high sensitive and high resolution ultrasound transducer to form a photoacoustic image. A variety of photoacoustic imaging embodiments have been established and reported so far, including photoacoustic microscopy (PAM), photoacoustic computed tomography, and photoacoustic endoscopy.20,2231 Among the different embodiments, optical-resolution photoacoustic microscopy (OR-PAM) is mostly applicable for iris vascular imaging as it can provide micron or even submicron high resolution images.

In addition to the imaging technology, the development of algorithms to accurately extract quantitative vascular information from the images, such as blood vessel diameter, vascular density (VD), and vascular tortuosity, is also key to directly reflect the diseased state of the underlying imaging tissue. We have previously reported a multiparametric quantitative vascular imaging algorithm32 to accurately extract quantitative vascular information from vasculature images of mouse back and ear obtained with our custom-built OR-PAM system.33 However, the reported quantitative algorithm based upon a two-dimensional (2-D) Hessian matrix is only suitable for vasculature images of objects with flat surfaces, while failing for objects with curved surfaces, such as the iris. In this study, to address the quantification of vascular images of the objects with the curved surfaces, we have developed a vascular quantification algorithm based upon a three-dimensional (3-D) Hessian matrix. We validate the performance of this quantitative vascular imaging method by comparing its performance against the performances of our previously developed 2-D Hessian matrix-based method and also the 3-D Hessian matrix-based method widely used in clinical CT vascular imaging.34

In this study, the in vivo photoacoustic imaging of rat iris was carried out with a custom-built OR-PAM system. The obtained images from this system were processed with our developed algorithm and the quantitative vascular information including vessel diameter distribution, VD, and vascular tortuosity were calculated and analyzed. The results shown in this study indicate that our developed algorithm would greatly facilitate the application of OR-PAM photoacoustic imaging for diagnosing ophthalmic diseases.


Materials and Methods


Photoacoustic Imaging of Rat Iris

A custom-built OR-PAM system was used to acquire all the imaging data in this manuscript. The system is composed of a nanosecond pulsed laser (GKNQL-532, Beijing Guoke Laser Co., Beijing, China), a single-element ultrasound transducer (V2022, Olympus-NDT, Kennewick, Washington), a precision motorized 3-D scanning stage (PLS-85, Micos, Eschbach, Germany), and a high speed data acquisition board (CS1422, Gage Applied Technologies Inc., Lockport). The detailed information of the OR-PAM system can be found in our earlier publication.32,33 One eight weeks healthy female Sprague Dawley (SD) rat (300 g) was selected for photoacoustic imaging, which remained anesthetized throughout the experiment using 1.5% isoflurane gas (Euthanex, Palmer, Pennsylvania) mixed with oxygen. The right eyelid of the rat was flipped inside out and fixed with adhesive tapes to expose the eyeball and 0.4% oxybuprocaine hydrochloride eye drops were used to anesthetize the eyeball during the imaging process. The imaging head of OR-PAM was positioned directly above the eyeball, as shown in Fig. 1(a), and the photoacoustic image of the entire iris was first obtained with the pupil positioned on the top surface of the eyeball, as shown in Fig. 1(b). Next, we manually rotated the eyeball to position the root part of the iris [yellow circle in Fig. 1(b)] on the top surface of the eyeball and acquired another image. This imaging procedure is similar to the real procedure followed in the clinics for the eye examination. In the clinics, patients are typically asked to expose the entire iris with the pupil at the center first and then are asked to rotate their eyeballs so that the region of interest (i.e., the iris root part) comes to the center for close observation. By manually rotating the rat eyeball and moving the iris root part to the top surface, we are not only simulating the real clinical scenario, but we are also obtaining high quality images since both the incident beam for photoacoustic signal excitation as well as the receiving transducer for signal detection are perpendicular to the region of interest [Figs. 1(c) and 1(d)].

Fig. 1

(a) Schematic of the imaging setup. (b) Photo of the rat eyeball during the imaging process. The yellow circle indicates the root part of iris, and the green arrows indicate the illumination laser. (c) and (d) Schematic comparison of the relative position between the incident beam and the sample before and after the rotation of the eyeball.


For the in vivo experiment, the laser frequency was 2 kHz, and the pulse energy was 200 nJ. As the laser focal depth is more than 300  μm below the sample surface to efficiently image iris vessels, and the numerical aperture is 0.1, the surface laser fluence is below 11  mJ/cm2, which is within the American National Standards Institute safety limit.24,35,36 To confirm that our imaging study caused no damage to the iris, the rat was euthanized and the imaged eyeball was removed after the completion of photoacoustic imaging to test for any damages. The 4% paraformaldehyde was used to fix the eyeball, followed by paraffin-cut sectioning and Masson’s trichrome staining to prepare the sample for further optical microscopic observation (Axio Lab.A1, ZEISS, Gottingen, Germany). In this study, all animal handling and experimental procedures conformed to the protocol approved by the Animal Study Committee of Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences.


Vascular Information Quantification Algorithm


Overall algorithm flowchart

The overall flowchart of our proposed algorithm is shown in Fig. 2. The vascular walls from the original 3-D volume data are extracted by background noise degradation, vascular signal identification, vascular signal enhancement, and region growing. On these extracted vessels, centerlines are identified. Finally, the quantitative vascular information, including vessel diameter, VD, and vascular tortuosity, are calculated from the identified centerlines. Details of each module in the flowchart are described in the next few subsections. All the postprocessings of the acquired data were performed using MATLAB (R2012b, Mathworks, Natick, Massachusetts) software on a PC with an Intel(R) Core(TM) i7 CPU at 3.9 GHz and a 64 GB RAM. It takes about 25 s to run a regular maximum amplitude projection (MAP) code and 25 min to run our developed quantitative code.

Fig. 2

The overall flowchart of the vascular information quantification algorithm.



Vascular signal extraction

  • Step 1: Background noise degradation based on open operation.

To remove the nonvascular background noise, such as discrete microvessel wall signals from the raw data, “open imaging” operation was used.37 The open imaging operation is suitable specifically for extracting bright objects in the dark background images and more details about open operation can be found in Ref. 37. We chose a disk-shaped structural element SE(u,v), with a typically used radius of three pixels, to remove the background noise from the input image I(x,y). For 3-D volume, we perform open operation layer-by-layer. On each layer, corrosion, also called erosion, of the image followed by the expansion, also called dilation, of the corroded image is performed using structural element SE(u,v) to get the denoised image. The whole operation can be described using Eq. (1):

Eq. (1)

where corroded image (ISE) with structural element SE is expanded (ISE) to get the final open operated image ISE.

  • Step 2: Vascular signal identification with 3-D Hessian matrix.

Before applying the 3-D Hessian matrix to identify the vascular signals from the 3-D volume, we first enhance the details of the vascular boundary of the open-operated 3-D volume, especially the small blood vessels, by applying a high frequency enhancement (HFE) filter as follows:

Eq. (2)

where Hhp is a high pass filter, a is the offset from the origin, and b is the weighing number of the high frequency signals (a=4, b=2).

The results of image processing on a photoacoustic image of an iris are shown in Fig. 10 in the Appendix for one lateral cross-section (i.e., cross-section perpendicular to the imaging depth direction). When an object with a flat surface is imaged (such as mouse ear in our previous study), most blood vessel segments project longitudinal tubular structures on the lateral cross-section. Our previous 2-D Hessian matrix-based algorithm recognized these longitudinal tubular structures as blood vessels and enhanced their signal. However, blood vessel segments in the iris, due to the curved shape of the eyeball, mostly project elliptical or circular shapes [indicated by yellow arrows in Fig. 10(b) in the Appendix] on the lateral cross-section. If our 2-D Hessian matrix-based method is used to process this image, the algorithm would classify these elliptical or circular shaped blood vessels as noise and will suppress them rather than enhance.

Thus, to overcome the limitations of the 2-D Hessian matrix, we have developed a 3-D Hessian matrix-based algorithm to process iris volume data. Before the 3-D Hessian matrix is computed at each voxel location, we perform 3-D filtering of the input volume data with a Gaussian kernel G(x,y,z,s), where s controls the size of the kernel (1s9), to eliminate high frequency variations in the image generally associated with noise. The 3-D Hessian matrix at each data point in a 3-D volume is a 3×3 matrix consisting of the second-order partial derivative of I(x,y,z,s):

Eq. (3)

where I(x,y,z,s) is the filtered grayscale value at location x, y, z in the input volume.

Three orthogonal vectors (e1,e2,e3) were chosen as the eigenvectors of the Hessian matrix, with one eigenvector corresponding to the axial direction of the blood vessel [or the direction with minimal value descendent of I(x,y,z,s)], and the other two eigen vectors orthogonal to the axial directions. The vascular information is obtained by calculating the three eigenvalues (λ1,λ2,λ3). Since the image feature along the direction of the blood vessels does not change substantially, the eigenvalue λ1 corresponding to this direction should have minimum value and close to 0. The other two eigenvalues, λ2 and λ3, corresponding to the other two orthogonal radial directions should be close to each other and their absolute values should be much larger than that of λ1 (we defined |λ2||λ3|). Thus, the vascular signals within the 3-D volume data can be identified using the following criteria:

Eq. (4)

where ϵ is set empirically to 0.8.

  • Step 3: Vascular signal enhancement.

The vascular signal enhancement, I(s), was carried out using Eq. (5):

Eq. (5)

where RA=|λ2|/|λ3|, RB=|λ1|/|λ2λ3|, and S=λ12+λ22+λ32 (RA, RB, and S are filters to suppress noisy, plate-like, and blob-like structures), and α, β, and γ are thresholds that control the sensitivity of the filter to RA, RB, and S for removing residual nonvascular signals. The previous study has shown that for image with vascular structures, α and β are generally fixed to 0.5 to enable better smoothing effect of the filter and γ is related to the dynamic range of the grayscale image, which is usually set to S/2.

Finally, the extracted feature map of image, i.e., the 3-D map of extracted vascular signals is obtained by repeatedly applying the Gaussian filter with different s to take into account different vessel sizes and then combining them using Eq. (6):

Eq. (6)

where smin=1 and smax=9 are different kernel sizes determined by the distribution of the vessel diameter and I(s) is determined using Eq. (6) for each s. It is to be noted that there will be multiple cases in terms of the selection of the three orthogonal vectors corresponding to each of the vascular voxels. As a result, the yield of Eq. (5) will also change. In practice, we did the calculation corresponding to all the cases in Eq. (5) and obtained multiple yields, and then we selected the maximum among all the yields for vascular signal enhancement, as indicated in Eq. (6).

To further enhance vascular signals and lower background noise, an intensity transformation function was applied to the feature map f to obtain the transformed image g given by

Eq. (7)

where m is an adjustable parameter that represents threshold value to determine background or vessel pixel (grayscale value) of the image and n is an adjustable scaling factor, which was empirically set between 19 and 21 to further scale up or suppress the value. When the intensity of a voxel is lower than m, the voxel will be treated as background and thus will be suppressed. In our experiments, the ideal value of m was obtained automatically by performing the clustering-based image thresholding.

  • Step 4: Region growing.

The voxel with the maximum signal intensity after vascular signal enhancement was selected as the seed point for 26 neighborhood region growing to obtain the binary volume data of the extracted vessels.


Vascular centerlines

We adopted an augmented fast marching method (AFMM) to extract the vascular centerlines.38 However, before applying this method, we separate the binary volume data into individual subdomains so that no connected components (i.e., connected vessels) exist in each subdomain. We apply AFMM to each subdomain to calculate the distance of each voxel present inside the vessel to the vessel boundary at both sides of the vessel to find the global maximum distance point, i.e., the voxel with the largest calculated distance from the boundaries. We apply AFMM again to calculate the distance of each voxel inside the vessel to the global maximum distance point and then update the voxels in the volume dataset with the calculated distance value. The gradient of the updated volume was calculated to get the gradient map of the vessels. Using the global maximum distance point and the gradient map, the vascular centerline in each subdomain was drawn by extending the global maximum distance point along the contour of the fastest descent of the gradient map. Finally, all the subdomains were assembled to obtain the skeleton of the vascular tree for the 3-D volume data.


Vascular parameter extraction

Two 512×512×130 subvolumes at the root part of the iris were selected for extracting quantitative vascular information, and the quantified parameters include vessel diameter, VD, and vascular tortuosity.

  • Parameter 1: Vessel diameter.

Since the vessel diameter is perpendicular to the centerline of the vessel, the vertical distances between the centerline and both sides of the vessel wall can be summed up to get the diameter of the vessel. Thus, by using this technique, we calculate the diameter of all the blood vessels in the selected subvolume.

  • Parameter 2: Vascular density.

The VD was calculated by first taking the 2-D MAP of the selected subvolumes. The VD is calculated from this projected 2-D binary image as follows:

Eq. (8)

VD=The number of pixels with signal intensity value equal to1N,
where N represents the total number of the pixels in the projected 2-D MAP image, which, in our case, is 512×512.

  • Parameter 3: Vascular tortuosity.

Three definitions of tortuosity are widely accepted for 3-D vascular data quantification, namely, the distance metric (DM), the inflection count metric (ICM), and the sum of angles metric (SOAM).39,40 DM is defined as the ratio between the actual path length of a vessel segment in each subdomain and the linear distance between the two ends of the vessel; ICM is defined as the ratio between DM and the number of the vessel’s inflection points; and SOAM is defined as the sum of the curvature at all voxels along the centerline of a vessel normalized by the vessel’s actual path length.


Algorithm Accuracy Evaluation

To evaluate the accuracy of our quantitative vascular parameter extraction using our developed quantitative algorithm, we used Image-J (version 1.51K; National Institutes of Health, Bethesda, Maryland) to manually extract the quantitative vascular parameters and compared results against those obtained with our algorithm. A male colleague blinded from our study manually performed the measurement five times on the same dataset and took the averaged value to exclude subjective influence from the results. The images provided to the colleague for accuracy evaluation are shown in Figs. 6(d) and 6(h).


Results and Discussion

Figure 3 shows the sagittal optical microscopic images of the anterior segment of the SD rat eyeball. The cornea and iris tissues are normally distributed, indicating no lesions or damages were caused due to photoacoustic imaging. Figure 3(b) is the enlarged view of the yellow dashed rectangular region in Fig. 3(a) and boundaries of different layers of the iris are highlighted with different colors. In Fig. 3(b), it can be seen that the iris is divided into three layers near the pupillary end: the stromal layer, the muscle layer, and the pigment epithelium layer, and two layers near the root end: stromal layer and pigment epithelium layer. In addition, we can also see in Fig. 3(b) that the diameters of the stromal vessels (indicated by the white regions in the image) are also different near the pupillary end and the root end of the iris. Thus, deriving quantitative vascular information of the subvolumes of the iris is important rather than that of the entire volume to evaluate the vascularity of the tissue at different regions.

Fig. 3

Sagittal optical microscopic images of the anterior segment of SD rat eyeball (a) using a ×10 objective lens. Black arrows indicate corneas (1), iris (2), ciliary body (3), and lens (4); red arrows indicate the root end (5) and the near pupillary end of iris (6). (b) Using a ×20 objective lens. Enlarged view of the yellow dashed rectangle area in panel (a), with the green arrows indicating the stromal layer (7), the muscle layer (8), and the pigment epithelium layer (9).


Figure 4 shows the depth encoded photoacoustic imaging results of iris vasculature with the pupil at the top of the eyeball. Figure 4(a) is the original image generated with the raw 3-D volume data. Figures 4(b)4(d) are the postprocessed images processed with our previously developed 2-D Hessian matrix-based algorithm [Fig. 4(b)], the conventional 3-D Hessian matrix-based algorithm used in clinical CT imaging [Fig. 4(c)], and the 3-D Hessian matrix-based algorithm developed in this study [Fig. 4(d)], respectively. Figures 4(e)4(j) are the enlarged view of the subareas in Figs. 4(a), 4(c), and 4(d), respectively. The subareas are indicated by the two white rectangles in Fig. 4(a). The smaller rectangular area in Fig. 4(a) corresponds to the area within the depth of focus of OR-PAM, while the larger rectangular area corresponds to the out of focus region. Significant information loss can be seen in Fig. 4(b), confirming that the 2-D Hessian matrix-based method is unsuitable for processing images with large surface curvature, as described in Sec. 2.2.2 of the method part. Compared to Figs. 4(b)4(d), significantly higher quality images are shown, when 3-D Hessian matrix-based methods were used for image processing. Furthermore, by comparing Figs. 4(c), 4(f), and 4(i) with corresponding Figs. 4(d), 4(g), and 4(j), respectively, we find the following advantages for our 3-D Hessian matrix-based algorithm over the conventional 3-D Hessian matrix-based algorithm used in clinical CT imaging: (1) the background noise is suppressed more efficiently, as indicated by the blue arrows in Figs. 4(c) and 4(d) due to the use of open operation and vascular signal enhancement; (2) blood vessels with different diameters are better differentiated, as indicated by the blue and green arrows in Figs. 4(e)4(g) due to the HFE filter, intensity transformation filter, and large Gaussian kernel size; and (3) the extracted blood vessels are more smoothed in Fig. 4(j) compared to Fig. 4(i), due to the larger Gaussian kernel size as well. All these advantages of our algorithm enable better extraction of the vessels from the out of focus region. The depth of focus of the imaging system is 260  μm. Based on our analysis, more than 95% of the blood vessels shown in Fig. 4 are within 1 mm distance from the focal plane of the system. The worst resolution within this range should be around 30  μm, which is good enough to accurately image most of the blood vessels in the iris. Furthermore, Fig. 4(j) shows that the blood vessels can be accurately extracted with our algorithm even for the most defocused region. So, overall, it is believed the performance of our imaging quantification algorithm remains at a high level for the whole iris area. Moreover, to evaluate how the noise affected the performance of our algorithm, we also applied our algorithm to a vascular image with high noises. The result (Fig. 11 in Appendix) shows that the blood vessels can still be accurately extracted even under the noisy condition. To demonstrate the performance of our algorithm in 3-D, we also load the postprocessed volume data (i.e., volume after the vascular signal enhancement step) into VolView (Kitware Inc., Clifton Park, New York) to display the data in 3-D format for easier evaluation (Fig. 5).

Fig. 4

Depth encoded photoacoustic imaging results of iris vasculature for the comparison of different image processing algorithms. (a) The original image generated with the raw 3-D volume data. (b)–(d) The postprocessing images with our previously developed 2-D Hessian matrix-based algorithm, the conventional 3-D Hessian matrix-based algorithm used in clinical CT imaging, and the 3-D Hessian matrix-based algorithm developed in this study. (e)–(g) The enlarged view of the same areas [smaller rectangle area in (a)] in panels (a), (c), and (d). (h)–(j) The enlarged view of the same areas [larger rectangle area in (a)] in panels (a), (c), and (d).


Fig. 5

Representative image frame from Video 1 (Video 1, MP4, 6 KB [URL:]).


Figure 12 in the Appendix shows the MAP image of the entire iris generated using raw 3-D volume data obtained with the root part of the iris on the top surface of the eyeball. Two 512×512×130 subvolumes from the root part of iris were selected for quantitative vascular information extraction and analysis, which correspond to the two green rectangles (labeled as area 1 and 2) in Fig. 12 in the Appendix. As shown in this figure, the two labeled areas are featuring two different vascular characteristics. Area 1 represents the iris part with disorderly distributed vasculature, while area 2 represents the iris part with radially distributed blood vessels. Figure 6 shows the image processing results after the key steps of the vascular information quantification algorithm developed in this study. Figures 6(a) and 6(e) are the original unprocessed MAP images of areas 1 and 2. Figures 6(b)6(d) and 6(f)6(h) are the images obtained after vascular signal enhancement, intensity transformation, and region growing steps of our algorithm. By comparing vascular signals in Figs. 6(a) and 6(e) with the corresponding vascular signals in Figs. 6(b) and 6(f), we can clearly see that the 3-D Hessian matrix-based method enhances the vascular signals throughout in the image. However, microvascular signals are still relatively weak and thus are not diagnostic. But, after intensity transformation, the microvascular signals are further enhanced, as can be clearly seen in Figs. 6(c) and 6(g). The background noises are completely removed and binary images after region growing are shown in Figs. 6(d) and 6(h). To validate the accuracy of our algorithm, the signal intensity profiles along the dashed lines in the original images [Figs. 6(a) and 6(e)] and the binary images [Figs. 6(d) and 6(h)] were plotted and compared with each other in Fig. 7. The blue line in Fig. 7 is from the original image and the red line is from the binary image. The matched peaks of the blue line and the red line indicate that our method did not introduce any distortion in the image and yet resulted in high quality image for accurate vascular information extraction. Furthermore, to validate the accuracy of our algorithm by comparing it to the algorithm used for CT, we also applied the CT algorithm to the images shown in Figs. 6(a) and 6(b) under the same procedure as that in Figs. 6 and 7, and put the final result as Fig. 13 in the Appendix. As can be seen by comparing Fig. 6(d) with Fig. 13(a) in the Appendix and Fig. 6(h) with Fig. 13(b) in the Appendix, the extracted blood vessels look much worse with the algorithm used for CT. Furthermore, by comparing Fig. 7(a) with Fig. 13(c) in the Appendix and Fig. 7(b) with Fig. 13(d) in the Appendix, it can be seen that the extracted images with our algorithm also agree more with the original images compared to the algorithm used for CT.

Fig. 6

The original unprocessed image and postprocessing images after vascular signal enhancement, intensity transformation, and region growing. (a)–(d) Correspond to area 1 in Fig. 12 in the Appendix. (e)–(h) Correspond to area 2 in Fig. 12 in the Appendix.


Fig. 7

Comparison of the intensity profiles along the dashed lines in the original and binary images in Fig. 6.


Finally, the quantified vascular parameters of the selected two subvolumes are shown in Fig. 8. Figure 8(a) shows the distribution of vessels with different diameters and Fig. 8(b) shows the VD and tortuosity information in our input volume. From Fig. 8(a), it can be clearly seen that microvessel (vessel diameter below 20  μm) density in subvolume 1 (area 1) is significantly higher than that in subvolume 2 (area 2), which agrees with the images in Figs. 6(d) and 6(h). In addition, in Fig. 8(b), we can see that the overall vascular densities (VD) in subvolumes 1 and 2 (area 1 and 2) are very close to each other, which is consistent with the findings in the previous reports.8,19 The quantification values shown in these figures also reveal information that is not easily captured by the naked eye, e.g., vascular tortuosity in subvolume 2 (area 2) is more severe than that in subvolume 1 (area 1), as indicated by all three definitions (DM, ICM, and SOAM) in Fig. 8(b).

Fig. 8

Quantified vascular parameters of the two subvolumes, i.e., subvolumes 1 and 2 (areas 1 and 2). (a) Vessel diameter distribution and (b) VD and three different definitions (DM, ICM, and SOAM) of vascular tortuosity.


In Fig. 9, we show the comparison between vessel diameter distributions obtained with our algorithm (series 1) versus manual measurement obtained using Image-J (series 2). In both Figs. 9(a) and 9(b), corresponding to subvolumes 1 and 2, respectively, distribution of vessels for various diameters is very close to each other. The correlation efficiency is very high and is equal to 0.962 and 0.929 (p0.01) for regions 1 and 2, respectively, based on Spearman’s correlation test, demonstrating the high accuracy of our algorithm. The VD obtained with Image-J is 0.3689 and 0.3779 for subvolumes 1 and 2, respectively. These densities are very close to those obtained with our algorithm (0.3766 and 0.3773). The accuracy of vascular tortuosity calculation has been evaluated in many previously reported studies,40,41 and we adopted the same tortuosity calculation method in our algorithm as that in the previous studies. Thus, our method, in addition to being accurate, is automatic with just enough user input thus can be readily deployed clinically. The developed algorithm can be incorporated into the high-end graphics cards with billions of operations to make it a real-time system. This is one of our next goals. Finally, it is to be noted that compared with our algorithm, the Image-J processing takes a much longer time due to the manual measurement.

Fig. 9

Vessel diameter distribution in subvolume 1: (a) and subvolume 2 (b) obtained with our algorithm (series 1) and manual measurement based on Image-J (series 2).




In vivo quantitative vasculature imaging of rat iris was accomplished in this study for the first time with a custom-built OR-PAM system and a developed 3-D Hessian matrix-based image processing algorithm. The quantitative vascular information of rat iris, such as vessel diameter, VD, and vascular tortuosity, was successfully extracted facilitating the IRIS disease diagnosis. The combination of photoacoustic imaging with the developed algorithm has the following advantages: (1) high contrast images of rat iris containing microvascular information can be obtained label free with OR-PAM due to the ultrasensitivity of photoacoustic imaging toward endogenous hemoglobin; (2) compared with the 2-D Hessian matrix-based algorithm, the application of the 3-D Hessian matrix in the algorithm enables significantly improved image quality for further accurate vascular information extraction; and (3) compared to the 3-D Hessian matrix-based algorithm used in the clinical CT imaging, the developed algorithm in this study is capable of suppressing background noise more efficiently, accurately differentiates blood vessels with various diameters, and readily extracts out of focus blood vessels. We confirmed all our findings by using a SD rat for in vivo iris imaging and since densely packed vasculature in rat iris closely relates to the real condition of human eyes, the findings in this study could be extended to human use in the near future. Compared with the previous studies that all used mice,24,35,36 our research has provided promising prospects for further pathological studies as it is more flexible to establish disease models on rat irises. Thus, results in this study indicate that the developed algorithm in conjunction with OR-PAM can potentially facilitate the study of ophthalmology and many other vascular-related diseases in vivo in the near future. In the next step, we will evaluate several diseased iris to establish correlation between quantitative parameters and disease process for easier diagnosis.



Figures in the appendix include the following: iris photoacoustic imaging result for one cross section that is perpendicular to the imaging depth direction (Fig. 10); comparison of a MAP figure with low signal-to-noise before and after processed with our algorithm (Fig. 11); MAP image of the entire iris with the root part of iris at the top surface of the eyeball (Fig. 12); and comparison of the vascular information extraction accuracy for our algorithm and the algorithm used in CT imaging (Fig. 13).

Fig. 10

In the Appendix: (a) iris photoacoustic imaging result for one cross-section that is perpendicular to the imaging depth direction and (b) enlarged view of the rectangle area in panel (a).


Fig. 11

In the Appendix: (a) MAP figure with low signal-to-noise ratio (SNR1.7) and (b) noise degraded figure processed by our algorithm.


Fig. 12

In the Appendix, MAP image of the entire iris with the root part of iris at the top surface of the eyeball.


Fig. 13

In the Appendix, (a) and (b) vessel extracted images (binary images) corresponding to Figs. 6(a) and 6(e) processed by the algorithm used for CT. (c) and (d) comparison of the intensity profiles along the dashed lines in the original images in Figs. 6(a) and 6(e) and binary images in Figs. 13(a), 13(b) in the Appendix.



The authors declare that there are no conflicts of interest related to this article.


The authors gratefully acknowledge the following funds: (1) National Natural Science Foundation of China (NSFC) grants 91739117 and 61405234, National Key Basic Research (973) Program of China grant 2015CB755500, Shenzhen Science and Technology Innovation grants JCYJ20170413153129570, JCYJ20160531175040976, JCYJ20150731154850923 and JCYJ20160422153149834 to C.L.; (2) NSFC grants 81522024 and 81427804, 973 Program of China grant 2014CB744503, Guangdong Natural Science Foundation grant 2014B050505013 and 2014A030312006 to L.S.; (3) Shenzhen Science and Technology Innovation grant JCYJ20150401145529037 to R.L.; (4) NSFC grant 61475182, Shenzhen Science and Technology Innovation grants JCYJ20150521144321005 and JCYJ20160608214524052 to X.G.; and (5) NSFC grant 31570952, Beijing Natural Science Foundation grant 3122010 to Z.L.



H. A. Quigley and A. T. Broman, “The number of people with glaucoma worldwide in 2010 and 2020,” Br. J. Ophthalmol., 90 262 –267 (2006). Google Scholar


Y. C. Tham et al., “Global prevalence of glaucoma and projections of glaucoma burden through 2040: a systematic review and meta-analysis,” Ophthalmology, 121 2081 –2090 (2014). Google Scholar


G. Spaeth, J. Walt and J. Keener, “Evaluation of quality of life for patients with glaucoma,” Am. J. Ophthalmol., 141 3 –14 (2006). Google Scholar


T. Vos et al., “Global, regional, and national incidence, prevalence, and years lived with disability for 310 diseases and injuries, 1990–2015: a systematic analysis for the Global Burden of Disease Study 2015,” Lancet, 388 1545 –1602 (2016). Google Scholar


K. L. Pepple et al., “Quantitative assessment of anterior segment inflammation in a rat model of uveitis using spectral-domain optical coherence tomography,” Invest. Ophthalmol. Visual Sci., 57 3567 –3575 (2016). Google Scholar


H. A. Quigley, “Angle-closure glaucoma–simpler answers to complex mechanisms: LXVI Edward Jackson memorial lecture,” Am. J. Ophthalmol., 148 657 –669.e1 (2009). Google Scholar


F. E. Seager, J. L. Jefferys and H. A. Quigley, “Comparison of dynamic changes in anterior ocular structures examined with anterior segment optical coherence tomography in a cohort of various origins risk factors of primary AC according to geographic origin,” Invest. Ophthalmol. Visual Sci., 55 1672 –1683 (2014). Google Scholar


Y. Zhang et al., “Dynamic iris changes as a risk factor in primary angle closure disease,” Invest. Ophthalmol. Visual Sci., 57 218 –226 (2016). Google Scholar


Y. Zhang et al., “Quantitative analysis of iris changes after physiologic and pharmacologic mydriasis in a rural Chinese population iris changes with mydriasis,” Invest Ophthalmol. Visual Sci., 55 4405 –4412 (2014). Google Scholar


Y. Oya, W. Sugiyama and N. Ando, “Anterior segment fluorescein angiography for evaluating the effect of vitrectomy for neovascular glaucoma,” Nippon Ganka Gakkai Zasshi, 109 741 –747 (2005). Google Scholar


S. Grisanti et al., “Intracameral bevacizumab for iris rubeosis,” Am. J. Ophthalmol., 142 158 –160 (2006). Google Scholar


S. Ishibashi et al., “Angiographic changes in iris and iridocorneal angle neovascularization after intravitreal bevacizumab injection,” Arch. Ophthalmol., 128 1539 –1545 (2010). Google Scholar


S. Dithmar and F. G. Holz, Fluorescence Angiography in Ophthalmology, 1 –26 Springer Science & Business Media, Heidelberg (2008). Google Scholar


A. H. Skalet et al., “Optical coherence tomography angiography characteristics of iris melanocytic tumors,” Ophthalmology, 124 197 –204 (2017). Google Scholar


E. Rosen and D. Lyons, “Microhemangiomas at the pupillary border: demonstrated by fluorescein photography,” Am. J. Ophthalmol., 67 846 –853 (1969). Google Scholar


W. J. Choi, Z. Zhi and R. K. Wang, “In vivo OCT microangiography of rodent iris,” Opt. Lett., 39 2455 –2458 (2014). Google Scholar


P. K. Roberts, D. A. Goldstein and A. A. Fawzi, “Anterior segment optical coherence tomography angiography for identification of iris vasculature and staging of iris neovascularization: a pilot study,” Curr. Eye Res., 42 (8), 1136 –1142 (2017). Google Scholar


R. M. Shtein et al., “Effect of tamsulosin on iris vasculature and morphology,” J. Cataract Refractive Surg., 40 793 –798 (2014). Google Scholar


H. Yang et al., “Quantitative study of the microvasculature and its endothelial cells in the porcine iris,” Exp. Eye Res., 132 249 –258 (2015). Google Scholar


L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science, 335 1458 –1462 (2012). Google Scholar


L. V. Wang and J. Yao, “A practical guide to photoacoustic tomography in the life sciences,” Nat. Methods, 13 627 –638 (2016). Google Scholar


J. Laufer et al., “In vivo preclinical photoacoustic imaging of tumor vasculature development and therapy,” J. Biomed. Opt., 17 056016 (2012). Google Scholar


L. Xi et al., “HER-2/neu targeted delivery of a nanoprobe enables dual photoacoustic and fluorescence tomography of ovarian cancer,” Nanomed. Nanotechnol. Biol. Med., 10 669 –677 (2014). Google Scholar


S. Hu et al., “Label-free photoacoustic ophthalmic angiography,” Opt. Lett., 35 1 –3 (2010). Google Scholar


S. Jeon et al., “In vivo photoacoustic imaging of anterior ocular vasculature: a random sample consensus approach,” Sci. Rep., 7 4318 (2017). Google Scholar


S. W. Lee, H. Kang and T. G. Lee, “Real-time display and in-vivo optical-resolution photoacoustic microscopy for ophthalmic imaging,” in 4th Int. Conf. on Bioimaging, 34 –38 (2017). Google Scholar


S. N. Hennen et al., “Photoacoustic tomography imaging and estimation of oxygen saturation of hemoglobin in ocular tissue of rabbits,” Exp. Eye Res., 138 153 –158 (2015). Google Scholar


W. Liu et al., “In vivo corneal neovascularization imaging by optical-resolution photoacoustic microscopy,” Photoacoustics, 2 81 –86 (2014). Google Scholar


X. Liu et al., “Optical coherence photoacoustic microscopy for in vivo multimodal retinal imaging,” Opt. Lett., 40 1370 –1373 (2015). Google Scholar


S. Y. Nam and S. Y. Emelianov, “Array-based real-time ultrasound and photoacoustic ocular imaging,” J. Opt. Soc. Korea, 18 151 –155 (2014). Google Scholar


C. Tian et al., “Noninvasive chorioretinal imaging in living rabbits using integrated photoacoustic microscopy and optical coherence tomography,” Opt. Express, 25 15947 –15955 (2017). Google Scholar


Z. Yang et al., “Multi-parametric quantitative microvascular imaging with optical-resolution photoacoustic microscopy in vivo,” Opt. Express, 22 1500 –1511 (2014). Google Scholar


J. Chen et al., “Blind-deconvolution optical-resolution photoacoustic microscopy in vivo,” Opt. Express, 21 7316 –7327 (2013). Google Scholar


A. F. Frangi et al., “Multiscale vessel enhancement filtering,” Lect. Notes Comput. Sci., 1496 130 –137 (1998). Google Scholar


N. Wu et al., “High-resolution dual-modality photoacoustic ocular imaging,” Opt. Lett., 39 2451 –2454 (2014). Google Scholar


B. Rao et al., “Hybrid-scanning optical-resolution photoacoustic microscopy for in vivo vasculature imaging,” Opt. Lett., 35 1521 –1523 (2010). Google Scholar


C. R. Gonzalez and R. Woods, Digital Image Processing, 445 –452 Pearson Education Inc., India (2002). Google Scholar


R. Van Uitert and I. Bitter, “Subvoxel precise skeletons of volumetric data based on fast marching methods,” Med. Phys., 34 627 –638 (2007). Google Scholar


E. Bullitt et al., “Measuring tortuosity of the intracerebral vasculature from MRA images,” IEEE Trans. Med. Imaging, 22 1163 –1171 (2003). Google Scholar


A. H. Parikh et al., “Correlation of MR perfusion imaging and vessel tortuosity parameters in assessment of intracranial neoplasms,” Technol. Cancer Res. Treat., 3 585 –590 (2004). Google Scholar


E. Bullitt et al., “Vessel tortuosity and brain tumor malignancy: a blinded study,” Acad. Radiol., 12 1232 –1240 (2005). Google Scholar


Huangxuan Zhao received his bachelor’s degree in mechanical engineering from Wuhan University of Technology in 2013. Currently, he is pursuing his PhD degree in the Department of Biomedical Engineering at Capital Medical University in China. His research interests include utilizing photoacoustic imaging for the detection of diseases and developing new surgical guidance methods for clinical applications.

Xiaojing Gong is an associate professor at Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences. He received his PhD degree in Fine Instruments and Mechanics from the University of Science and Technology of China in 2007. He is currently focusing on the research of biophotonics, including photoacoustic endoscopy, optical coherence tomography/microscopy, terahertz biomedical imaging, etc.

Liang Song is a professor at Shenzhen Institutes of Advanced Technology (SIAT), Chinese Academy of Sciences. He is the founding director of The Research Lab for Biomedical Optics in SIAT, and director of the Shenzhen Key Lab for Molecular Imaging. He studied at Washington University, St. Louis, from 2006 to 2010 and received his PhD degree in biomedical engineering. He has authored more than 50 peer-reviewed journal articles and has been granted more than 10 patents in the field of photoacoustic imaging.

Zhicheng Liu is a professor of biomedical engineering at Capital Medical University (China) and founding director of Beijing Key Laboratory of Fundamental Research on biomechanics in clinical application. His main research interests are the mechanical properties of biological soft tissue and developing new surgical guidance methods for clinical applications.

Chengbo Liu is an associate professor at Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences. He received both his bachelor’s and PhD degrees from Xi’an Jiaotong University, each in 2012 in biophysics and 2007 in biomedical engineering. During 2009 to 2011, he was a visiting scholar at Duke University, working on tissue spectroscopy for early cancer diagnosis. His current research interests include multiscale photoacoustic imaging technology development and photoacoustic translational research.

Biographies for the other authors are not available.

© 2018 Society of Photo-Optical Instrumentation Engineers (SPIE) 1083-3668/2018/$25.00 © 2018 SPIE
Huangxuan Zhao, Guangsong Wang, Riqiang Lin, Xiaojing Gong, Liang Song, Tan Li, Wenjia Wang, Kunya Zhang, Xiuqing Qian, Haixia Zhang, Lin Li, Zhicheng Liu, and Chengbo Liu "Three-dimensional Hessian matrix-based quantitative vascular imaging of rat iris with optical-resolution photoacoustic microscopy in vivo," Journal of Biomedical Optics 23(4), 046006 (23 April 2018).
Received: 16 January 2018; Accepted: 19 March 2018; Published: 23 April 2018 Logo
Cited by 32 scholarly publications.
Get copyright permission  Get copyright permission on Copyright Marketplace
Algorithm development

3D image processing

Iris recognition


Vascular imaging

Blood vessels

Image processing

Back to Top