1.

Introduction

More than 67 million people are diagnosed with glaucoma globally, making it the most commonly seen irreversible ocular disease worldwide.^1–4 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.^5–9 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.^10–12 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,^16–18 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,¹⁹ 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,22–31 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 algorithm³² to accurately extract quantitative vascular information from vasculature images of mouse back and ear obtained with our custom-built OR-PAM system.³³ 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.³⁴

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.

2.

Materials and Methods

2.1.

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\mu \mathrm{m}$ below the sample surface to efficiently image iris vessels, and the numerical aperture is 0.1, the surface laser fluence is below $11\mathrm{mJ}/{\mathrm{cm}}^2$, 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.

2.2.

Vascular Information Quantification Algorithm

2.2.1.

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.

2.2.2.

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.³⁷ 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 $\mathrm{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 $\mathrm{SE}(u,v)$ to get the denoised image. The whole operation can be described using Eq. (1):

Eq. (1)

$$I\circ \mathrm{SE}=(I\ominus \mathrm{SE})\oplus \mathrm{SE},$$

where corroded image ($I\ominus SE$) with structural element SE is expanded ($I\oplus SE$) to get the final open operated image $I\circ SE$.

• 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)

$$H_{\mathrm{hfe}}=a+b{H}_{\mathrm{hp}},$$

where $H_{\mathrm{hp}}$ 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 ($1\le s\le 9$), 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\times 3$ matrix consisting of the second-order partial derivative of $I(x,y,z,s)$:

Eq. (3)

$${▿}^{2}I=\left[\begin{array}{ccc}\frac{{\partial }^{2}I(x,y,z,s)}{\partial x^2}& \frac{{\partial }^{2}I(x,y,z,s)}{\partial x\partial y}& \frac{{\partial }^{2}I(x,y,z,s)}{\partial x\partial z}\\ \frac{{\partial }^{2}I(x,y,z,s)}{\partial y\partial x}& \frac{{\partial }^{2}I(x,y,z,s)}{\partial y^2}& \frac{{\partial }^{2}I(x,y,z,s)}{\partial y\partial z}\\ \frac{{\partial }^{2}I(x,y,z,s)}{\partial z\partial x}& \frac{{\partial }^{2}I(x,y,z,s)}{\partial z\partial y}& \frac{{\partial }^{2}I(x,y,z,s)}{\partial z^2}\end{array}\right],$$

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

Three orthogonal vectors $({\mathrm{e}}_1,{\mathrm{e}}_2,{\mathrm{e}}_3)$ 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 $(\lambda 1,\lambda 2,\lambda 3)$. Since the image feature along the direction of the blood vessels does not change substantially, the eigenvalue $\lambda 1$ corresponding to this direction should have minimum value and close to 0. The other two eigenvalues, $\lambda 2$ and $\lambda 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 $\lambda 1$ (we defined $|\lambda 2|\le |\lambda 3|$). Thus, the vascular signals within the 3-D volume data can be identified using the following criteria:

Eq. (4)

$$|\lambda 1|<ϵ\ll |\lambda 2|;\lambda 2\approx \lambda 3<0,$$

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)

$$I(s)=\{\begin{array}{cc}0,& {\lambda }_2>0\text{or}{\lambda }_3>0\text{or}|{\lambda }_1|\ge 0.8\\ [1-\exp (-\frac{R{A}^2}{2{\alpha }^2}\left)\right]\exp (-\frac{R{B}^2}{2{\beta }^2})\times [1-\exp (-\frac{S^2}{2{\gamma }^2}\left)\right],& \text{else}\end{array},$$

where $RA=|\lambda 2|/|\lambda 3|$, $RB=|\lambda 1|/\sqrt{|\lambda 2\lambda 3|}$, and $S=\sqrt{{\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2}$ (RA, RB, and $S$ are filters to suppress noisy, plate-like, and blob-like structures), and $\alpha$, $\beta$, and $\gamma$ 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, $\alpha$ and $\beta$ are generally fixed to 0.5 to enable better smoothing effect of the filter and $\gamma$ 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)

$$f=\underset{s_{\min }\le s\le s_{\max }}{\max }|I(s)|,$$

where $s_{\min }=1$ and $s_{\max }=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)

$$g=\frac{1}{1+{(m/f)}^n},$$

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.

2.2.3.

Vascular centerlines

We adopted an augmented fast marching method (AFMM) to extract the vascular centerlines.³⁸ 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.

2.2.4.

Vascular parameter extraction

Two $512\times 512\times 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)

$$\mathrm{VD}=\frac{\text{The number of pixels with signal intensity value equal to}1}{N},$$

where $N$ represents the total number of the pixels in the projected 2-D MAP image, which, in our case, is $512\times 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.

3.

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 $\times 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 $\times 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 $\sim 260\mu \mathrm{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\mu \mathrm{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: http://dx.doi.org/10.1117/1.JBO.23.4.046006.1]).

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\times 512\times 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\mu \mathrm{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 ($p\ll 0.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).

4.

Conclusion

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.