1.

Introduction

Oxygen utilization is altered in many disorders of the retina.^{1, 2, 3, 4, 5, 6} The present work is motivated by the need for developing automatic and noninvasive ways of measuring the percent blood ${\mathrm{O}}_2$ saturation from retinal images. Evaluations of blood oxygenation along vascular segments of individual artery and vein networks can show, respectively, the distributions of oxygen supply and extraction from local regions. Here we describe automated image analysis techniques based on dual-wavelength oximetry⁷ and present results from normal subjects of indices from retinal vessels that are proportional to ${\mathrm{O}}_2$ saturation. Additional information concerning the vessel diameter, which has previously been correlated with blood flow changes,⁸ is forthcoming in the analysis. Future applications for an automatic system could include large-scale screening for assessing risk for progression of the stages of diabetic retinopathy, postoperative monitoring or image-guided intervention with patients receiving laser treatments for proliferative retinal vessel disease, and application to change detection studies of the retinal blood oxygen saturation.

Existing methods of vessel oximetry acquire images at multiple wavelengths, where measurements are limited to a small region in the image.^{7, 9, 10, 11, 12} In the present method, digital images of continuous vessel segments of several millimeters length, containing artery and vein pairs, are acquired. The automated image analysis identifies serial segments of the vessels in two images, one sensitive to changes in the percent of $\mathrm{Hb}{\mathrm{O}}_2$ present in the blood, and the other a reference image that is not sensitive.⁷ Figure 1 shows two images acquired simultaneously at $570\mathrm{nm}$ and $600\mathrm{nm}$ . The $570\text{-}\mathrm{nm}$ image, which is insensitive to $\mathrm{Hb}{\mathrm{O}}_2$ change, has higher contrast between the blood vessels and the background than does the image at $600\mathrm{nm}$ , where differential contrast between the artery and vein comes from the higher and lower $\mathrm{Hb}{\mathrm{O}}_2$ content in each vessel. The apparent OD of the blood inside vessels is measured from reflectance values taken from inside and outside of the vessels. Algorithms that trace the blood vessels, reconstruct the vascular hierarchy, and measure the apparent OD of the blood contained in the retinal vessels are described below.

Fig. 1

Sample images from subject 3 at oxygen-insensitive and oxygen-sensitive wavelengths. (a) Insensitive image acquired at $570\mathrm{nm}$ . (b) Sensitive image acquired at $600\mathrm{nm}$ and registered to the $570\text{-}\mathrm{nm}$ image. The $570\text{-}\mathrm{nm}$ image shows higher contrast between the blood vessels and the background, and reflected light from the vein (upper vessel) and artery (lower vessel) are comparable at this wavelength. At $600\mathrm{nm}$ , the vessel contrast is low and the artery appears lighter than the vein.

Many advances and refinements in technique for retinal oximetry have been made over the past four decades.^{7, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20} Hickam ^{9, 10} were the first to report studies of retinal oxygenation in man by noninvasive measurements, using photographic methods to record vessel OD at oxygen-sensitive and insensitive wavelengths, and ratio analysis to obtain relative measures of oxygen saturation. Later, Laing ¹¹ extended this method by showing that blood saturation varied linearly with the ratio over physiological and hypoxic ranges of saturation. Delori developed a three-wavelength oximetry method to obtain absolute oxygen saturation values from single vessel segments.¹² Delori’s oximeter was the first to demonstrate altered blood saturation in optic atrophy.¹³ A method to perform retinal vessel oximetry by spectroscopic recording was introduced by Schweitzer ¹⁴ Relative measures of oxygen saturation along retinal vessels using a ratiometric measurement were reported by Beach, ⁷ using a digital imaging technique to simultaneously record oxygen-sensitive and insensitive images. Jensen¹⁵ and Critten ¹⁶ have employed the ratio technique successfully to measure oxygen saturation changes in retinal vessels. Recently Khoobehi ¹⁷ employed hyperspectral imaging to monitor the relative spatial changes in retinal oxygen saturation of the retina and optic disc tissue in the monkey. Uzumcu ¹⁸ obtained measurements at six wavelengths and estimated the parameters related to the oxygen saturation using an iterative algorithm. A scanning laser ophthalmoscope that images at four wavelengths was used by Smith ^{19, 20} They suggested that the choice of wavelengths is critical and used an optimized wavelength set. The analysis was done on a small selected region from the images. Studies of the relationship between oxygen utilization and retinal disorders are also extensively discussed in the literature.^{1, 2, 3, 4, 5, 6}

The motivation for our image analyses was to enable higher levels of automation and precision using the dual-wavelength approach with digital imaging. We have also added the capability to automatically distinguish the oxygen saturation values from the different levels of the vessel hierarchy. Our vessel hierarchy algorithm automatically reconstructs the vascular structure and hence enables the comparison at the corresponding levels in the vascular hierarchy. This ability should be especially useful in venous networks for localization of areas of relative hypoxia.

Several methods have been proposed in the literature for detection of vascular features from retinal images.^{21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46} These methods can be broadly classified into four categories.The first approach is based on adaptive filtering or segmentation.^{21, 22, 23, 24, 25, 26, 27, 28, 29, 30} In the second method, referred to as vessel tracking or tracing, local image properties are used to trace the vasculature recursively starting from an initial point.^{31, 32, 33, 43, 44, 45, 46} Some of these algorithms exploit a priori knowledge on the position of one or more reference points belonging to the vessels. In some of the methods commonly used in quantitative coronary analysis, the initial and end points of the vessel are entered manually.³² The disadvantage of these methods is that they are not fully automated. The third approach uses mathematical morphology to detect long tubular objects representing blood vessels.^{38, 39, 40, 41} The fourth approach uses classification methods to segment vessels.⁴² Chaudhuri assume a Gaussian cross-sectional profile for the vessels and detect them by convolving with a matched filter.²¹ Hoover use the matched filter response with an adaptive threshold probing.²² Gang detect vessels using amplitude modified second-order Gaussian filter.²³ Other methods based on matched filters are described in Refs. 24, 25. Directional edge detectors were used by Li, ²⁸ while local image characteristics were exploited for detecting the blood vessels in window-based methods described in Refs. 29, 30. A parameterized variation of the exploratory approach known as “snakes” is discussed in Refs. 34, 35. Aylward and Bullitt³⁶ describe a method based on intensity ridge traversal to extract centerlines of tubular objects. Jiang and Mojon³⁷ implement adaptive local threshold using multiple threshold verification and a curvilinear structure model based on angle, width, and contrast. Zana and Klein use mathematical morphology to detect vessel-like structures.³⁸ Morphological processing is also described in Refs. 39, 40, 41 for finding blood vessels from retinal images. An unsupervised classification algorithm was used by Fontaine for segmenting the vasculature.⁴² In this paper we use the exploratory tracing algorithm developed by Can. ⁴³ and further developed in Refs. 44, 45, 46. The algorithm will be described in detail in the next section.

2.

Fully-Automatic Retinal Vessel Tracing Algorithm

Previous research in this group has developed algorithms for rapid automated tracing and feature extraction from retinal images.^{43, 44, 45, 46} These algorithms were originally inspired by the need to construct computer vision tools for assisting ophthalmic procedures such as laser retinal surgery and perimetry. The tracing method is based on adaptive exploratory processing of the image, directly at the image intensity level, avoiding image-wide pixel-processing operations such as gray-scale filters, edge operators, morphological filters, etc. The algorithm is designed such that the computations are restricted to the most relevant and promising locality of pixels. In the present work we extend these algorithms to develop an automated system to measure the oxygen content in the blood from images of the retina obtained at two different wavelengths.

The tracing algorithm is based on a parallel edge model. Over short distances, the vessel segments can be assumed to be linear and having parallel edges oriented along a particular direction. These parallel edges can be detected by using a set of two-dimensional correlation kernels. The exploratory algorithm starts from seed points that are selected automatically using a two-step detection and verification stage from a uniformly spaced grid. Each seed is assigned a direction during the verification stage and the tracing algorithm proceeds by tracing in these initial directions from each valid seed point. The algorithm takes a step in this direction and then determines the direction to proceed at that point. Given a point ${\vec{p}}^k$ with initial direction $s^k$ , the directional templates are used to estimate the next point ${\vec{p}}^{k+1}$ and direction $s^{k+1}$ in a recursive fashion. Figure 2(a) shows how the algorithm takes a step in the current direction and then estimates an orientation at the new point. The exploratory algorithm terminates if it meets certain stopping conditions. Interested readers are referred to Refs. 43, 44 for details of the tracing algorithm. At each point in the vessel centerline, the algorithm estimates the width, boundary points and also the direction. From the centerline network, the branching and crossover points are detected as the intersection points of the centerline trace segments. These points are referred to as landmarks and help in determining the hierarchy of the vascular tree as described below. Each landmark has two or more vessel segments associated with it and also stores the angles that these segments make with the landmark. This is illustrated in Fig. 2(b).

Fig. 2

Recursive tracing algorithm and landmarks. (a) The vessel edges can be modeled to be parallel over small distances. The templates for two different orientations are shown in panel (a). Separate templates are used for the left and right edges. The direction of tracing at point $p^k$ is $s^k$ . From that point the algorithm takes a step in the direction $s^k$ . The new direction $s^{k+1}$ is found at $p^{k+1}$ by the application of the templates. (b) Landmarks obtained from the vessel centerlines. Each landmark has a set of traces associated with it and also the angle the vessel extends at it. This information is used by the vessel linking algorithm to reconstruct the vascular hierarchy.

3.

Reconstructing the Vascular Hierarchy

The detected trace segments from the exploratory tracing algorithm do not give any information about the vascular hierarchy. Also the same vessel might be traced in parts as several fragments. If we are able to automatically identify fragments that belong to the same vessel and also its branches, we would be able to make much more accurate measurements from a single connected vessel rather than making the measurements from the individual fragments. The information about the vascular hierarchy is also desirable for picking appropriate vessels to compare the oxygen saturation for detecting changes. The traces are therefore processed using a recursive algorithm that reconstructs the vascular structure.

The traces are first sorted according to the decreasing order of the product of their lengths and average width and put in a priority list. This is done to make sure that the vessel linking algorithm starts from segments that are significantly wide and long. Starting from the first vessel in the sorted list, the algorithm checks whether this segment has any landmarks associated with it. If there are no landmarks associated with this parent segment, it is given a unique name and removed from the list. On the other hand, if the parent segment has one or more landmarks associated with its end points, the segments that are associated with these landmarks are processed and named one at a time as follows. While selecting from multiple segments, associated through the same landmark, the segments that are closer in direction to the current segment are processed first followed by segments that make larger angles with the current segment. Also if there are multiple approximately collinear segments associated with the current segment, it is considered probable that the other segments are part of another vessel and are therefore, not named. When naming the segments, the following rules are followed.

In order to avoid naming a single segment more than once, the named segments are removed from the list and eliminated from further consideration. Once all the connected segments under one parent segment have been processed, the algorithm checks the list and process the next segment. The algorithm proceeds until there are no more elements in the list to be processed. A pseudo code description of the above algorithm is provided in Fig. 3 , and Fig. 4 shows the original image, traced image, initial trace segments and their numbers and the reconstructed vessel structure. The vessel linking algorithm is limited by the accuracy of the initial tracing and landmark extraction and cannot recover from errors in the tracing. Another source of error is significant changes in the width or intensity of parts of the same vessel. For example, segment 10 in Fig. 4(c) is really a daughter of the vein marked as V. But in this case, the intensity of segment 10 was significantly higher from that of V, and hence the algorithm mistook it to be part of another vessel. Subsequent segments (64 and 15) are also named based on the name given for segment 10. Hence naming errors at the higher levels of the hierarchy are propagated to the lower levels.

Fig. 3

Pseudo code for the recursive vessel segment naming algorithm. Starting with the longest traced segment, find the segments connected to that segment and name them in a recursive fashion. This results in a hierarchical representation of the vascular network.

Fig. 4

(a) The original oxygen-insensitive image from subject 1 to be processed using the recursive vessel linking algorithm. Vein (top vessel) and artery (small vessel below the vein). (b) Traces and landmarks for the image shown in (a). The green points are the vessel centerlines and the red lines indicate the detected landmarks (bifurcations and crossovers). (c) Initial trace fragments and trace numbers for the traces shown in (b). Each fragment is colored differently for clarity. (d) Illustrating the hierarchical linking of vessel traces to reconstruct the vascular tree. The hierarchical names given to the fragments shown in (c) are shown in red along the original fragments. The starting number for the vessel names have been replaced with A and V for clarity in differentiating arteries and veins.

4.

Optical Density Measurements

Images of the eye are obtained at $570\mathrm{nm}$ and $600\mathrm{nm}$ . The tracing algorithm is applied to the image obtained at $570\mathrm{nm}$ . This is because at this wavelength, there is better contrast between the vessels and the background resulting in better detection of the vasculature. The $600\text{-}\mathrm{nm}$ image is registered to the $570\text{-}\mathrm{nm}$ image using a similarity transformation. This is a rigid transformation that allows for translation, rotation, and scaling between two images. The registration parameters for our images were determined using the image registration feature in IPLab (Scanalytics).

The apparent OD can be calculated on the vessel segments by finding the minimum intensity inside the vessel and the average outside intensity. For each point ${\vec{p}}^k$ in a vessel centerline, let ${\vec{r}}^k$ and ${\vec{l}}^k$ denote the right and left boundary points. In order to find the pixel with the minimum reflection at this particular cross section of the vessel, the pixels along the straight line joining ${\vec{r}}^k$ and ${\vec{l}}^k$ are examined. The minimum intensity inside the vessel at ${\vec{p}}^k$ is given by:

Extravascular light reflection at point ${\vec{p}}^k$ is measured from points lying at one vessel diameter distance on either side of the vessel boundary. It is given by:

Fig. 5

Illustration of the ODR calculation method. The blue dots on the vessel indicate the centerline point ${\vec{p}}^k$ , the left boundary point, ${\vec{l}}^k$ and the right boundary point, ${\vec{r}}^k$ . The unit vector in the direction perpendicular to the vessel orientation is denoted by $u_{\perp }^k$ and $D$ is the diameter of the vessel at ${\vec{p}}^k$ . The two red dots on the outside of the vessel indicate the points at which the extra vascular reflection is measured and the black boxes represent the local neighborhood over which the average is calculated. The white dot represents the minimum intensity point inside the vessel.

Fig. 6

Illustration of the OD calculation method. (a) The centerlines in white and the pixels across each cross section used for finding the minimum intensity in blue (b) The minimum intensity pixel found across each cross section in white and the points that are used for calculating the extravascular reflectance in red.

Case 1. In this method, the apparent OD is calculated by averaging the inside intensities and the outside intensities and finding one OD for the segment. Let $I_{\mathrm{in},\mathrm{seg}}$ and $I_{\mathrm{out},\mathrm{seg}}$ be the average inside and outside values for trace segment $i$ with $M$ points:

Case 2. This variation of the method is aimed at dealing with any illumination gradient across the image. In this method, we calculate the OD for each centerline point in the trace segment and then average the ODs to obtain the overall OD for the segment. The OD at centerline point $p^k$ is given by:

The $600\text{-}\mathrm{nm}$ image that was registered to the $570\text{-}\mathrm{nm}$ image is then processed as described above to get ${\mathit{OD}}_{600}$ and the ${\mathit{SD}}_{600}$ . The centerlines and boundary points used for this are the same as the ones obtained from tracing the $570\text{-}\mathrm{nm}$ image. The ratio of the OD at $600\mathrm{nm}$ to that at $570\mathrm{nm}$ is computed as the ODR:

5.

Description of the Retinal Image Analysis

Experimental data were obtained from seven volunteer subjects and the recordings confirmed to an approved institutional review board protocol. Subjects’ eyes were dilated with 1% topical tropicamide and 2.5% phenylephrine. Recordings were made simultaneously from retinal artery vein pairs using a dual-wavelength camera that is described in Ref. 7. The present instrument uses an inexpensive $16\text{-}\text{bit}$ digital camera (Apogee U-2) and upgraded optics to improve wavelength selection and the pixel resolution is $9\times 9\mu \mathrm{m}$ . Images containing larger retinal vessels and side-branch vessels were recorded during room air breathing and again while the subject breathed pure oxygen $(10\mathrm{L}∕\min )$ from a mask that was attached over the nose and mouth. The same retinal areas were selected for low and high oxygen. Arterial saturation was read with a pulse oximeter using a finger probe to assess changes in arterial saturation. The retina was illuminated with the xenon flash from the fundus camera (Topcon TRC).

Recorded dual-wavelength retinal images were processed first by separating them into $570\text{-}\mathrm{nm}$ and $600\text{-}\mathrm{nm}$ images as shown in Fig. 1, and these separate images were registered. The centerline traces of vessels in the $570\text{-}\mathrm{nm}$ image were found using the tracing algorithm. Results from the tracing were used as inputs to the vessel hierarchy algorithm. Figures 4(a) and 4(b) show, respectively, a $570\text{-}\mathrm{nm}$ vessel image with individually identified vessel segments overlaid by centerline tracings and connecting landmarks that are used to link segments based on their hierarchical structure. The full complement of detected segments in the image was then resolved into serialized segments along the vascular tree as shown in Fig. 4(d). Vessel segment hierarchy is defined by us using a naming scheme where the large parent vessel is identified with a single number, daughter vessels having the same number and the designation D which follows, and multiple daughters or side branches that occur along the vessel moving distally being assigned higher numbers. This naming scheme is marked on the image of vessels to confirm validity.

The algorithm for vessel hierarchy has been validated by an ophthalmologist (one of the authors). For each image, the vessel names were superimposed on the segments and presented for validation. Major vessels and its branches were selected and the segments that were given a name inconsistent with the convention were marked out. The success rate of the vessel linking algorithm is defined as the number of trace segments correctly named by the algorithm to the total number of trace segments. For the images considered, the success rate was found to be 87%. It has to be noted that the success of the vessel linking algorithm is limited by accuracy of the vessel tracing algorithm. Specifically, out of the 13% of the errors made by the algorithm, 7% were caused due to the errors in finding the bifurcation and crossing points by the tracing algorithm. The other errors were caused due to the big differences in width and intensity of the branches of the same vessel, in which case the algorithm mistook the branch to be part of a different vessel.

Along centerline traces, lines perpendicular to the vessel orientation were followed to sample minimum-valued pixels that estimate the reflected light from the blood column as shown in Fig. 6(a). Along the same line, approximately one vessel diameter from the outer wall, the extravascular reflection was averaged from neighboring pixels. OD was determined either by averaging values from inside and outside the vessel and then applying Eq. 7 for Case 1, or by finding the OD for each pair of values along the length of the segment and then averaging these as in Eq. 9 for Case 2, in either case ending up with a particular OD for each segment. Case 2 should reduce the influence of light gradients on OD values since these are determined at points along vessels rather than from area-averaged values. Centerline traces from $570\text{-}\mathrm{nm}$ images were applied to the registered $600\text{-}\mathrm{nm}$ and ODs were obtained for $600\text{-}\mathrm{nm}$ images. ODRs were then calculated as described above.

6.

Determination of Oxygen Change Sensitivity

Data from six subjects were used to determine the oximeter response resulting from changing the inspired gas from room air to pure oxygen. Table 1 compares ODRs from normal and high oxygen conditions in ten artery and eleven vein segments of the superior temporal vessel arcades. Systemic arterial oxygen saturation averaged 97% in the subjects during room air breathing and 100% during pure oxygen. ODs were determined from averaged intensities along the vessel segments (Case 1). These vessels began one disk diameter from the rim of the optic disk. The ODR of the arteries and veins decreases when breathing gas was switched from room air to pure oxygen. The mean change for the artery was 0.0226 ODR units. This change is very close to the expected 0.0236 ODR units for a blood saturation increase of three percentages, assuming an instrument sensitivity of 0.00787 ODR/change in percent saturation that was established with a similar methodology and vessel analysis using Case 1.⁷ The high degree of agreement between the two methods confirms that our automated analysis gives the same result as the previous vessel tracking technique which required user intervention to pick vessels and determine the size of each scanned segment before performing the analysis.

Table 1

Comparison of ODR (Case 1) at different SO2 levels for the selected arteries and veins. SO2 was read from a pulse oximeter with a finger probe. The ODR decreases as the systemic SO2 increases. The vessel naming convention is that A (V) represents a first level vessel, and A-D-* (V-D-*) represents a daughter vessel branching off from the main vessel. The first number in the name given by the vessel linking algorithm has been replaced with A and V for clarity in distinguishing arteries (A*) and veins (V*) . SO2 in room air=97% , in pure O2=100% .

Oxygen breathing caused a greater change in the ODRs of veins; the decrease observed for the veins was 3.31 times that of the artery, which corresponds to an increase in venous saturation of approximately 10%. The venous saturation changes should significantly exceed those in arteries since pure oxygen causes the tissue oxygen source to switch from hemoglobin to dissolved oxygen in plasma, with the result that there is less desaturation of the venous blood. Table 2 gives results when the vessel ODs were found by the alternate method of first determining OD at each location along the segments and then averaging these (Case 2). Case 2 had little effect on the outcome; ODRs of arteries and veins (room air) were 2.3% higher and 1.2%lower, respectively, than the corresponding values from Case 1. In arteries, there was no significant difference between the two methods in the response to breathing oxygen. In veins, Case 2 resulted in a 9% greater decrease in the ODR. Figure 7 shows vessels from two subjects for room air (left panels) and pure ${\mathrm{O}}_2$ (right panels), with the ODR calculated using Case 2 and vessel diameters marked on selected veins and arteries. Both the datasets show the reduction in the ODR, which corresponds to increased saturation, and in the reduced diameter, which resulted from oxygen-induced vasoconstriction. In the lower panels, the responses to ${\mathrm{O}}_2$ were greater; reduced vessel diameters are easily discerned from the images and both parameters changed by larger degrees. Table 3 gives the intermediate steps for the ODR calculation for the two subjects shown in Fig. 7. Although the table has many entries, all of the values are computed by the automatic algorithm, and the user only needs to analyze the results.

Fig. 7

Optical density ratios (ODR) and the measured vessel widths (pixels) in parentheses for selected arteries (A) and veins (V) from subject 5 (a,b) and subject 6 (c,d). (a,c) Room air breathing, (b,d) pure ${\mathrm{O}}_2$ breathing.

Table 2

Comparison of ODR (Case 2) at different SO2 levels for the selected arteries (A*) and veins (V*) . The vessel naming convention is the same as in Table 1. SO2 in room air=97% , in pure O2=100% .

Table 3

Intermediate ODR calculation (Case 2) for the segments shown in Fig. 7. The table shows the name of the segment, the index of the original trace segments that were linked together by the algorithm, the individual ODs at 570nm and 600nm , and the ODR. The standard deviation (SD) is calculated over the points in each segment. The vessel naming convention is the same as in Table 1. SO2 in room air=97% , in pure O2=100% .

Repeatability of this method should be assessed from artery data since the saturation change is more consistent than in veins. For results from arteries in Table 1, the coefficient of variation $(\mathrm{COV}=\mathrm{SE}∕\mathrm{MEAN})$ for saturation changes equals 0.21 $(N=10)$ . This measure of variability in part reflects differences in the actual rise in retinal arterial saturation among the subjects that resulted from breathing pure oxygen.

The automated method also allows studying the changes in vessel diameters with oxygen saturation. The vessel diameters are estimated by the tracing algorithm at each trace point and the mean width for each segment is calculated. Table 4 shows the change in diameter of the blood vessels under room air and pure oxygen breathing for the six subjects. The mean decrease was 7.5% for the arteries and 9.93% for the veins for the increase of 3% in oxygen saturation.

Table 4

Comparison of the vessel diameters at different O2 levels. The percentage decrease in the diameter is calculated as (D1−D2)∕D1100* , where D1 and D2 are the measured diameters of the vessel at room air and pure O2 breathing, respectively, in pixels. The vessel naming convention is the same as in Table 1. SO2 in room air=97% , in pure O2=100% .

Plotted in Fig. 8 are values of the ODR versus distance along a vessel segment using Case 2, for air and pure ${\mathrm{O}}_2$ breathing. The expected reduction in the ODR with increased blood oxygen saturation is evident, while a gradient in saturation along the vessel is also observed for both cases. These data show how it would be possible to measure oxygen gradients along retinal vessels, or even discontinuous changes resulting from vessel occlusion.

Fig. 8

ODR plotted against position along a vein segment during room air (diamond symbols) and pure ${\mathrm{O}}_2$ (circles) breathing for subject 4. In this case, pure ${\mathrm{O}}_2$ increased the vein saturation (reduced ODR). The vessel showed a continuous positive gradient in saturation left to right (negative ODR gradient).

Vessel traces have been applied to larger images of the retina that were taken with the full frame of the CCD camera using the $55\text{-}\mathrm{deg}$ field-of-view setting on the fundus camera. Figure 9 shows images of temporal vessels with trace points included. In Fig. 9(a), the $570\text{-}\mathrm{nm}$ recording wavelength produces high contrast between both artery and vein against the fundus background. These traces followed the larger vessels and the numerous smaller side branches faithfully. In Fig. 9(b) the $600\text{-}\mathrm{nm}$ recording, which produces lower vessel contrast, is also faithfully traced. In this, the only case of a thinned retina that we have worked with, our method correctly identified retinal vessels in both the oxygen-insensitive and oxygen-sensitive images. Most of the small traces seen in Fig. 9(b) between the large artery and the vein (darker vessel) do connect with one of the large vessels, hence these are likely to be smaller branching arteries. Even though the myopic eye in Fig. 9 was traced faithfully, there is a possibility that the algorithm might mistake choroidal pigment bands in myopic subjects as vessels. Hence, more tests with this kind of subjects are necessary to know whether the algorithm will always succeed where there are choroidal pigment bands.

Fig. 9

(a) Illustration over larger areas of the retina with traces of vessel interior (white points) and exterior (red points). Oxygen-insensitive image showing vessels arch over the macula (dark region at bottom) from subject 7. The nerve fiber layer is visible near vessels. (b) Oxygen-sensitive image $\left(600\mathrm{nm}\right)$ from subject 7, where vessels arch below the macula. Retinal thinning as a condition of myopia produces intensity variations (light areas). The traces were obtained from the corresponding $570\text{-}\mathrm{nm}$ image and are superimposed on the registered $600\text{-}\mathrm{nm}$ image. Based on these traces, the minimum intensity points inside the vessels (indicated in white) were then determined by looking at the intensities in the $600\text{-}\mathrm{nm}$ image.

7.

Conclusions

The main contribution of this work is the adaptation of automated tracing algorithms for oxygen measurement from dual wavelength images. Automation for identification and analysis of retinal vessels is important for development of clinically useful technology for diagnosis. The steps taken thus far show that oxygen changes are forthcoming from segments of prominent vessels, including arteries and veins in dual-wavelength images. The vessel tracing and linking algorithms make it possible to obtain OD and oxygen change measurements from larger areas than were previously described in the literature. The vessel linking algorithm starts from a main vessel segment and names its side branches and bifurcations hierarchically. The algorithm was found to perform well in all the subjects under study, and we believe that this is the right start for solving complex problems associated with vessel hierarchies. Also our vessel tracing algorithm is able to find the vessel diameters with sufficient precision to detect oxygen vasoconstriction. This is significant since the vessel diameter is an important parameter in the analysis of retinal blood flow. Automated OD measurements have also been successful in smaller $(<50\mu \mathrm{m})$ veins that branch from the main vessels. Although there is room for improvement, the result show it is feasible for automatic image analysis to indicate areas of relative hypoxia based on distributions of venous blood saturation. One problem that remains is to successfully determine ${\mathrm{O}}_2$ saturation values in the smaller artery segments with automation. The very low vessel OD, strong central reflection and variation in the extravascular reflection at $600\mathrm{nm}$ make this determination difficult.

The choice of measurement wavelengths in the described method are such that the retinal vessels can be easily traced by measuring density against the fundus background, and sensitivity to saturation change is high. Retinal vessels become more visible against the fundus background at wavelengths where hemoglobin light absorption is strongest, between 540 and $580\mathrm{nm}$ where our reference wavelength is set. The sensing wavelength, which is longer than $580\mathrm{nm}$ , provides less vessel contrast but gives the highest ratio of light absorption between saturated and desaturated blood. At wavelengths longer than $610\mathrm{nm}$ , the retinal arteries become difficult to detect. Turbid media causes significant signal loss by scattering as measuring wavelength decreases toward blue wavelengths. Thus the present choice of wavelengths should perform reasonably well in patients with moderately turbid media, although high turbidity limits our method as it does color fundus imaging.

Three wavelength methods have been reported that provide absolute calibration of blood oxygen saturation from retinal vessels.¹² By comparison, the two wavelength method described here provides relative saturation indices using the optical density ratio (ODR). The absolute sensitivity of the ODR to saturation changes was previously determined.⁷ Thus, the method can be used to determine saturation changes in single vessels, or between two vessels. The artery-vein saturation difference (A-V difference) is an important parameter for oxygen utilization. [Venous saturation can be estimated from the A-V difference and the sensitivity of the ODR if the arterial saturation is known, such as during the breathing of pure oxygen.] Venous saturation during normoxia can then be assessed from the decrease in the ODR relative to hyperoxia. Simultaneous imagery of the fundus at two wavelengths can be easily acquired using image splitting techniques; this is significantly more difficult to achieve for higher numbers of wavelength.

Results from the automated analysis were essentially similar to earlier findings with dual-wavelength imaging that employed an automated analysis that had to be initiated manually.⁷ An area for improvement is in the registration of the $570\text{-}\mathrm{nm}$ and $600\text{-}\mathrm{nm}$ images. Highly accurate feature based registration algorithms have been developed by our group for retinal images.^{47, 48, 49} These algorithms depend on reliable feature extraction from the images. In the case of dual wavelength images, the central reflex for arteries in the $600\text{-}\mathrm{nm}$ image is so strong that the two sides of the arteries get traced separately as two segments. Although minimum values we seek are along either of these traces, it would be better if the centerline of the vessel was traced as one segment. Since in the $570\text{-}\mathrm{nm}$ image, this same vessel gets traced as a single segment through the centerline, the registration algorithm would tend to associate one or the other of the traces on each side of center in the $600\text{-}\mathrm{nm}$ image to the centerline from the $570\text{-}\mathrm{nm}$ image. This could result in wrong correspondences during registration, causing convergence to a wrong solution. Further research is planned to make the feature extraction better in retinal images with “hollow” vessels, which would allow the registration to be done automatically. Our results indicate that any light gradients present with the current method of fundus camera imaging should not significantly influence the calculation of vessel ODRs, which are the measure here for relative oxygen saturation.

Future work will also concentrate on refining the methods to enable analysis of smaller vessels, particularly the more distal networks of veins that drain blood from local areas of the retina. These vessels carry blood away from the retinal areas that could be affected by altered metabolic activity. This information would allow identification of local hypoxia that is associated with disorders of the retina, and in conjunction with automated techniques, allow for screening of disorders in clinical rather than research environments. The ODR values, while not absolute measures of saturation, are highly sensitive to saturation differences that would be present between areas of normal and low oxygenation. Thus, a diagnostic image that automatically reveals areas of relative hypoxia on the fundus image is one of the important goals of this work. For the other goal of investigating oxygen utilization in the normal and diseased retina, the present stage of development can lead to an image-based measurement of saturations and vessel diameters across the retina. We have demonstrated the feasibility of the methods on normal retina in this paper. Future measurements in patients with retinal disorders will be performed to assess clinical potential for this method.