Optimal wavelengths for subdiffuse scanning laser oximetry of the human retina

Abstract. Retinal blood vessel oxygenation is considered to be an important marker for numerous eye diseases. Oxygenation is typically assessed by imaging the retinal vessels at different wavelengths using multispectral imaging techniques, where the choice of wavelengths will affect the achievable measurement accuracy. Here, we present a detailed analysis of the error propagation of measurement noise in retinal oximetry, to identify optimal wavelengths that will yield the lowest uncertainty in saturation estimation for a given measurement noise level. In our analysis, we also investigate the effect of hemoglobin packing in discrete blood vessels (pigment packaging), which may result in a nonnegligible bias in saturation estimation if unaccounted for under specific geometrical conditions, such as subdiffuse sampling of smaller blood vessels located deeper within the retina. Our analyses show that using 470, 506, and 592 nm, a fairly accurate estimation of the whole oxygen saturation regime [0 1] can be realized, even in the presence of the pigment packing effect. To validate the analysis, we developed a scanning laser ophthalmoscope to produce high contrast images with a maximum pixel rate of 60 kHz and a maximum 30-deg imaging field of view. Confocal reflectance measurements were then conducted on a tissue-mimicking scattering phantom with optical properties similar to retinal tissue including narrow channels filled with absorbing dyes to mimic blood vessels. By imaging at three optimal wavelengths, the saturation of the dye combination was calculated. The experimental values show good agreement with our theoretical derivations.


Introduction
Hemoglobin oxygen saturation in the blood flowing through the retinal vasculature is an important physiological parameter involved in the pathophysiology of numerous retinal diseases including diabetic retinopathy, 1,2 glaucoma, 3-5 retinitis pigmentosa, [6][7][8] age-related macular degeneration, 9 and retinal vessel occlusions. 10 Several retinal diseases result in reduced oxygen circulation in the retina causing hypoxia, 11 which is one of the key drivers of angiogenesis. As such, measuring the retinal blood oxygen saturation will enable monitoring of the development of hypoxia and thus will allow preventing the loss of retinal tissue due to vasoproliferation through timely therapeutic interventions. Furthermore, although damage occurring to the retina in the advanced stages of many retinal diseases associated with hypoxia has been established, oxygenation measurements could also aid in understanding the onset or presence of systemic diseases. For example, a recent study reported that elevated oxygen saturation in retinal blood vessels is an indicator of Alzheimer's disease. 12,13 The earliest attempts to demonstrate retinal vessel oximetry date back several decades. 14,15 Various imaging modalities, such as fundus camera, [16][17][18][19][20][21][22][23] hyperspectral imaging, [24][25][26] scanning laser ophthalmoscope (SLO), 27,28 adaptive optics SLO, 29 and optical coherence tomography (OCT) 30, 31 have since then been employed to estimate oxygen saturation in retinal blood vessels. These noninvasive techniques all involve recording light reflected from the fundus and measuring the amount of light absorbed by the blood flowing through the retinal vessels at multiple wavelengths. Retinal oxygen saturation is typically calculated from reflection intensities recorded by a point or array detector, and the error on the measured intensities propagates to a saturation error. Using a detailed error analysis, we aim to identify optimum wavelengths that will yield minimum uncertainty on the saturation estimation for a given measurement error in the intensities.
One factor that has not been considered in retinal oximetry to date is the fact that hemoglobin is not distributed homogeneously in the tissue but confined to discrete blood vessels. This effect is particularly relevant for smaller blood vessels that occupy only a small portion of the total probed volume in a quasiconfocal/subdiffuse scheme. As a result, a fraction of the incident photons are reflected from the tissue without ever encountering a blood vessel, causing the apparent absorption coefficient of blood packed in small vessels to be different compared with when the blood would have been homogeneously distributed through the probed tissue. Such an effect is called the pigment packaging effect, and several research studies have experimentally verified this effect. [32][33][34] We have included this effect in retinal oximetry algorithms to increase the accuracy of the saturation estimation, in particular for small blood vessels located deeper in the retina, probed using a (sub)diffuse measurement technique. Large and superficially located blood vessels, whose diameters are comparable with the detection aperture, are not expected to suffer from this effect as the probability that the detected photons have not traveled through the blood vessel is negligible. In this paper, we present modified retinal oximetry equations that include the pigment packaging effect and evaluate the effect of pigment packaging on saturation estimation. We then identify optimum wavelengths for retinal oximetry using a detailed error analysis of the oximetry equations including the pigment packaging effect and identify the wavelengths that give minimum saturation error for both discussed geometries [Figs. 1(a) and 1(b)] and minimize the saturation offset in the presence of the pigment packing effect.
To validate our theoretical predictions, we have fabricated a thin, multilayer retina-mimicking phantom with narrow channels representing the blood vessels in the retina. The reduced scattering coefficient ðμ 0 s Þ of the scattering medium in the layers was chosen to be similar to the expected values for the retina. The retina-mimicking phantom was placed in a model eye. An SLO using a supercontinuum laser as the light source was used to scan the artificial retina. Two dyes of suitable concentrations were tailored to have absorption coefficients comparable with oxy-and deoxyhemoglobin with an artificial isosbestic point at 548 nm. The procedure for selecting optimum wavelengths to estimate the saturation within this artificial vessel was applied to this phantom and three optimal wavelengths (488, 548, and 612 nm) were selected based on this analysis. Using the resulting wavelengths, the experiments were conducted on the phantom. There was a good agreement between the predicted and real values of saturation of the different dye combinations using these three optimal wavelengths.

Theory of Retinal Oximetry and Identifying
Optimum Wavelengths

Theory of Retinal Oximetry
Figure 2(a) shows the absorption spectrum of the oxy-and deoxyhemoglobin having various isosbestic points, i.e., wavelengths where the absorption of oxy and deoxyhemoglobin is equal. Figure 2(b) shows a typical fundus image of a healthy human eye obtained with 570-nm illumination. A general method for measuring the absorbance of blood at a particular blood vessel in the retina is shown in the inset of Fig. 2 where the intensity profile of the blood vessel [dotted line in Fig. 2(b)] is plotted as a function of the location on the retina. The scattering properties of blood have been reported in Bosschaart et al., 35 where it was shown that although blood has a high scattering coefficient, it also has a high scattering anisotropy, resulting in a reduced scattering coefficient of ≈1 to 2 mm −1 . Further, the reduced scattering coefficient of blood is close to the reduced scattering coefficient of the surrounding tissue 36 and about an order of magnitude smaller than the absorption coefficient of blood in the 450-to 650-nm wavelength range. In a retinal imaging experiment, the measured intensity Iðx o ; λÞ (see Fig. 1) from a spatial location x o in the retina for a given illumination wavelength λ depends on the incident illumination intensity distribution I i ðx i ; λÞ at the point x i , and the backscattered reflectance function of the retina Rðx o ; x i ; λÞ. We consider a quasiconfocal or subdiffuse scheme, where a narrow, focused beam illuminates a point x i in the retina and the detection pinhole are collected from a larger volume around a point x o . Letx o be the location of an averaged volume of tissue surrounding x o that contributes to the collected intensity, Iðx o ; λÞ. The extent of volume averaging for any location concerning the spatial distribution of the optical properties at that location is determined by the incident illumination profile I i ðx i ; λÞ and the detection pinhole. Under these conditions, we write the measured intensity I as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 1 4 9 Iðx o ; λÞ ¼ The function Rðx o ; λÞ describes the reflectivity of the tissue averaged over a small tissue volume at the location x o due to scattering and absorption. Let x b and x t denote the recorded Fig. 1 Two different vessel geometries in subdiffuse retinal oximetry. (a) A large blood vessel with a diameter comparable to the collection area present superficially in the retina. In this case, the probability that a detected photon has not passed through the blood vessel is negligible, and hence, the pigment packing effect is not observed. (b) A small blood vessel embedded in the scattering medium (retinal tissue) is illuminated with a narrow beam, whereas the collection aperture is much larger than the illumination aperture and the vessel diameter. In this case, the blood vessel occupies only a small portion of the total probed volume and hence the pigment packing effect should be considered.
locations at the center of a blood vessel and in the adjacent tissue, respectively. Then, the light collected from a tissue location x t is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 6 3 ; 3 8 8 Iðx t ;λÞ ∝ I i ðx t ;λÞ · Rðx t ;λÞ ¼ I i ðx t ;λÞ · R½μ s ðx t ;λÞ;pðx t ;θ;λÞ; (2) where we have stated that Rðx t ; λÞ is a function of the average scattering coefficient ½μ s ðx t ; λÞ of the volume of tissue around x t (denoted byx t ) at the wavelength λ and of the scattering phase function [pðx t ; θ; λÞ] of the same tissue volume. We have assumed that there is no blood present in tissue volumē x t and that any other absorbing molecules within the tissue volume x t in the neural retina have a negligible contribution to the reflected intensities. Similarly, the light collected from the center of a blood vessel location x b is given as where we have stated that Rðx b ; λÞ is a function of the average scattering coefficient μ s ðx b ; λÞ of the volume of tissue around x b (denoted byx b ) at the wavelength λ, of the scattering phase function pðx b ; θ; λÞ of the same tissue volume, and the effective absorption coefficient of the tissue volume containing the blood vessel μ a ðx b ; λÞ. We can now write the relative optical density OD λ of the tissue at the blood vessel location compared with the surrounding tissue at a particular wavelength as the logarithm of the ratio of tissue and blood vessel reflectivity at the same wavelength, i.e., Here, we have assumed that the incident light intensity was the same at both locations and that the attenuation of the reflected light due to the blood withinx b is governed by modified Beer-Lamberts law, 37,38 with hL eff ðλÞi, the effective path length of photons traveling through the probed volume before reaching the detector. The effective path length hL eff ðλÞi depends on the absorption and scattering properties of the probed volume.
GðλÞ is a factor that accounts for any apparent increase or decrease in the ODs purely due to scattering differences within the volumes x b and x t i.e., due to the difference in scattering properties of blood within the blood vessel compared with the surrounding tissue. As in the wavelength range 450 to 650 nm, the reduced scattering coefficient of blood is similar to that of the surrounding tissue, we do not expect G to be a significant factor, in particular for the case shown in Fig. 1(b), where the blood vessel occupies only a small fraction of the interrogated tissue volume. Equation (4) assumes that all additional reflections occurring internally in the system due to optics and the stray reflections from the cornea are accounted for, e.g., by subtracting a reference measurement. The absorption coefficient μ a ðx b ; λÞ of the blood is a function of saturation (S) of the blood (i.e., the fraction of the oxygenated hemoglobin concentration to the total concentration of hemoglobin) and is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 9 2 μ a ðx b ; λÞ ¼ ½S · μ HbO 2 a ðλÞ þ ð1 − SÞ · μ Hb a ðλÞ: The parameters μ HbO 2 a ðλÞ and μ Hb a ðλÞ are the absorption coefficients of oxy-and deoxyhemoglobin, respectively (assuming a concentration of 150 mg of hemoglobin in 1 mL of blood) and are well-known functions of wavelength. Combining Eqs. (4) and (5) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 6 3 ; 6 9 7 OD λ ¼ hL eff ðλÞi · ν · ½S · μ HbO 2 a ðλÞ þ ð1 − SÞ · μ Hb a ðλÞ þ GðλÞ; where ν is the fraction of interrogated tissue occupied by the blood vessel. 33,34 We define the ratio ρ λ 1 jλ 2 of ODs at two wavelengths λ 1 and λ 2 as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 6 3 ; 6 1 8 If the scattering and absorption coefficients at the two wavelengths are close to each other, L eff ðλ 1 Þ ≃ L eff ðλ 2 Þ and Gðλ 1 Þ ≃ Gðλ 2 Þ, we will later discuss the validity of this assumption. Under this condition, we can write S as a function of ρ λ 1 jλ 2 as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 6 3 ; 4 6  ) is an offset, if any, to the saturation due to scattering. As discussed earlier, for visible wavelenghts we expect GðλÞ → 0, and in the remainder of the analyses, we only focus on the first term in the right-hand side of Eq. (8).
Although not mathematically necessary, it is advantageous to choose one of the wavelengths as an isosbestic wavelength to simplify the calculations, especially for error propagation purposes. If λ 2 is taken to be the isosbestic wavelength, In that case, Eq. (8) can be written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 6 3 ; 2 9 2 S λ a jλ i ¼ where we have denoted the saturation sensitive "anisosbestic" wavelength as λ a . We will now first consider the geometry as shown in Fig. 1(a), where the blood vessel especially a large, superficial blood vessel occupies the majority of the probed volume.

Optimum wavelengths without pigment packing
If there is a large blood vessel with a diameter close to or greater than the size of the collection aperture directly below the illumination point, as shown in Fig. 1(a), the assumption that the blood vessel occupies only a portion of the probed volume does not hold. In such geometries, the blood vessel occupies majority of the probed volume. Monte Carlo simulations 39,40 have suggested that in this type of geometry with a quasiconfocal aperture, the photons reaching the detector have predominantly been backscattered from within the blood vessel. In this case, hL eff ðλÞi in Eq. (6) becomes the average backscatter path length hL bs ðλÞi from within the blood vessel (and blood volume fraction ν becomes 1) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 3 2 6 ; 7 1 9 OD λ ¼ hL bs ðλÞi · ½S · μ HbO 2 a ðλÞ þ ð1 − SÞ · μ Hb a ðλÞ þ GðλÞ: Under the condition that scattering coefficient and absorption coefficient of blood is approximately equal at both wavelengths λ a and λ i , hL bs ðλ a Þi ≈ hL bs ðλ i Þi. The saturation in this case is given by Eq. (9). Turning now to the error on the saturation, this error is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 3 2 6 ; 6 0 7 where Δρ λ a jλ i is given by propagating the error on OD a and OD i in Eq. (7) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 3 2 6 ; 5 4 2 The error in the ODs is given by the underlying error of the reflected intensities Iðx t ; λÞ and Iðx b ; λÞ as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 3 2 6 ; 4 6 8 Due to very low absorption at wavelengths >650 nm, there is a poor contrast between the recorded intensity at the tissue location and the blood vessel location that results in OD → 0. This results in large errors in ρ [following Eq. (12)] and therefore, in saturation estimation. Absorption is very high in the Soret band (430 nm) resulting in most of the light being absorbed by the blood resulting in a very poor signal-to-noise ratio (SNR) at the blood vessel locations. Most importantly, the blue light hazard for the wavelengths <450 nm substantially limits the maximum permissible exposure of the retina for light 41 of these wavelengths. Thus, we exclude the wavelengths <450 nm and <650 nm for the remainder of the analysis. From Eqs. (11)- (13), it is evident that the standard error of the mean recorded intensities by the detector (ΔI∕I or standard deviation/mean) propagates to an error in the ODs and then to an error in ρ, and thus contributes to the saturation error ΔS. For example, a 1.0% error in the intensities ðΔI∕I ¼ 0.01Þ propagates to an error in the ODs ðΔOD ¼ 0.014Þ. If OD λ 1 ¼ 0.5 and OD λ 2 ¼ 0.7, this gives 3.4% error in ρ according to Eq. (12), this error in ρ propagates to saturation error in a complex manner ΔS through Eq. (11). Figure 3 shows the calculated saturation error for different oxygen saturation levels in the retina when 506 nm was used as the isosbestic wavelength and the anisosbestic wavelength varied over the 450-to 650-nm ranges, assuming a 1% measurement error in the intensity. For low saturation values [0.0 and 0.5], choosing the anisosbestic wavelength in the ranges 458 to 474 nm gives a saturation error <5%, while for high saturation values [0.5 and 1.0], choosing the anisosbestic wavelength in the ranges 591 to 599 nm gives a saturation error <5%. A similar trend was found for other isosbestic points shown in Fig. 2(a). Note, however, that the condition that the absorption coefficients are similar at μ a ðλ a Þ ≈ μ a ðλ i Þ, so that hL bs ðλ a Þi ≈ hL bs ðλ i Þi is met best by choosing 506 nm as the isosbestic wavelength [ Fig. 2(a)]. Figure 1(b) shows a geometry, where a small blood vessel is embedded in the scattering retina. In this case, the blood vessel occupies only a fraction of the total probed volume by the imager. As a result, Eq. (6) (which assumes a homogeneous distribution of hemoglobin throughout the probed tissue volume) does not hold as the absorbing hemoglobin molecules are not homogeneously distributed but are contained in a discrete package (the blood vessel) that occupies only part of the probed tissue volume. This packaging of absorbing molecules effectively flattens the apparent absorption spectrum for high absorption coefficients, and is called the pigment packaging effect. It has previously been shown that this effect can be modeled through the incorporation of a correction factor to the absorption spectrum. [32][33][34] Equation (5) in this case is modified to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 6 3 ; 2 9 7 μ a ðx b ; λÞ ¼ ν · ½S · μ HbO 2 a ðλÞ þ ð1 − SÞ · μ Hb a ðλÞ · C λ ; (14) where the correction factor C λ 34,42 is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 6 3 ; 2 5 3 C λ ðS; λ; D ves Þ

Optimum wavelengths including pigment packing
where D ves is the blood vessel diameter. The OD can now be written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 6 3 ; 1 5 1 The saturation S can be written as [ignoring the term α · G (α ¼ ) as was done for Eq. (9)] E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 3 2 6 ; 7 5 2 where S is now also a function of the blood vessel diameter D ves , due to the addition of the diameter-dependent correction factor. A detailed error analysis of Eq. (13) is necessary to find the optimal wavelengths. The error propagation equations gives the error in saturation ΔS as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 3 2 6 ; 6 5 9 The error ΔS is a function of S, D ves , and λ a for a chosen λ i . Now we proceed to analyze each term in Eq. (18) and its contribution to the saturation error ΔS. For an image-based estimation of the blood vessel diameter D ves (explained in Sec. 3.3) using an SLO technique, we chose a constant error on the estimation of the diameter ðΔD ves Þ as 12 μm based on spatial resolution of the system (see Table 1). The error on ρ λ a jλ i is given by error propagation in Eq. (12). Figures 4 and 5 show the color-coded graph of the calculated errors in saturation as a function of D ves and λ a considering a 1% error ðΔI∕IÞ in the intensities when 506 nm (lowest μ i a in the 450-to 650-nm interval) and 548 nm (highest μ i a in the 450-to 650-nm interval) were used as isosbestic points, respectively. For simplicity, all the saturation error values ðΔSÞ > 0.10 were set to 0.10 to emphasize the combinations with accuracy better than 0.1. From Fig. 4, we see that when 506 nm is used as an isosbestic wavelength, there are two spectral bands optimal for oximetry. For low saturation values (0.00 to 0.25), either of the bands, 460 to 480 nm or 590 to 600 nm could be used. But for higher saturation levels (0.75 and 1.00), only the 589 to 600 nm gives ΔS < 0.1. When 548 nm is used as an isosbestic wavelength (Fig. 5), 589 to 600 nm is not ideal for detecting the low saturation values and 460 to 480 nm is not suitable for high saturation values. Figure 6 shows the saturation error for a 50-μm [Figs. 6(a)-6(c)] and a 100-μm [Figs. 6(d)-6(f)] blood vessel for different saturation levels. From these graphs, we see that that all the five isosbestic wavelengths perform similarly for estimating fully oxygenated blood (S ¼ 1.00) when the 590-to 600-nm wavelength range is used as an anisosbestic wavelength. However, for estimating low saturation values, 522, 586, and 506 nm perform marginally better than 548 and 569 nm when the 460-to 480-nm wavelength range is used as an anisosbestic wavelength. Figure 7 shows ΔS values as a function of λ a ðλ i ¼ 506 nmÞ for a fully oxygenated 50-and 100-μm blood vessel for different 0.00 0.25 0.50 0.75 1.00 S Fig. 3 Saturation error ΔS as a function of λ a when estimating different oxygenation levels from 0% to 100%, with 1% error on the recorded intensities. A wavelength of 506 nm was chosen as the isosbestic wavelength. This graph holds for geometry described in Fig. 1(a), where the photons are mostly backscattered from a blood vessel.  It can be seen in this figure that to reach saturation accuracy levels better than 5% for 100% oxygenation when measuring near 594 nm, the required accuracy level of the mean intensity was 1% for both the 50-and 100-μm vessels.
We note that the optimal wavelength bands calculated in the presence, as well as in the absence of pigment packaging, are very similar. In fact, it is observed that if the ratio of correction factors C λ i ∕C λ a in Eq. (17)    [590 and 600] bands, where C 506 ∕C λ a ≈ 1 would serve the purpose of minimizing the saturation error independent of the effect of pigment packing. For such wavelength combinations, Eq. (9) can be used even in the presence of pigment packing without applying any correction. The dark blue region in the color map in Figs. 8(a) and 8(d) denotes the situation when C 506 ∕C λ a ≈ 1, i.e., when pigment packing has no effect. Unfortunately, this condition is satisfied for different wavelengths in the [460 and 480] and [590 and 600] wavelength bands for different saturation levels. Figures 8(b) and 8(e) show the saturation error [Eq. (11)] in the [460 and 480] and [590 and 600] wavelength bands, respectively, assuming a 1% measurement error of intensities and when 506 nm was used as an isosbestic wavelength. Figures 8(c) and 8(f) show the combined effect of saturation error from a 1% measurement error and ignoring pigment packing, indicating that using the 470-to 506-nm pair, the combined error is minimized for saturation in the range [0.00 and 0.55], and using the 592-and 506-nm pair, the combined error is minimized for saturation in the ranges [0.40 and 1.00]. Thus, the triad of wavelengths 470, 506, and 592 nm appear to be optimal for retinal oximetry. At these wavelengths, (i) the blood absorption coefficients are approximately equal, and as a consequence, C 506 ∕C λ a ≈ 1, which minimizes the effect of pigment packaging on saturation estimation. Similar absorption and scattering coefficients at these wavelengths satisfy the assumption that hL bs ðλ a Þi ≃ hL bs ðλ i Þi and hL eff ðλ a Þi ≃ hL eff ðλ i Þi, which are conditions for our algorithms to hold, (ii) the blood absorption coefficients are sufficiently high to yield a good vascular contrast in the images and, finally, (iii) the saturation error is low because the error in measured intensities has a minimized propagation of the error to OD and ρ.

Experimental Validation
To validate our theoretical predictions, we have fabricated a thin, multilayer retina-mimicking phantom with narrow channels representing blood vessels in the retina. The retina-mimicking phantom was placed in a model eye. A scanning laser ophthalmoscope (SLO) built in-house using a supercontinuum laser as the light source was used to scan the artificial retina at different wavelengths, to experimentally validate our theory.

Scanning Laser Ophthalmoscope-Description of the System
A multiwavelength SLO [43][44][45][46] was developed for the continuous acquisition of en face images of a phantom retina as shown in Fig. 9. The pixel rate of the system for imaging the phantoms was 30 kHz, although the maximum pixel rate of the system was 60 kHz. The size of a single image frame was 512 × 512 pixels. Thus, each image took 8.7 s to record. The essentially static nature of the phantoms precluded motion artifacts in the recorded image. The key system parameters are shown in Table 1. Using a relatively large core diameter of 100 μm in the detection compared with a single-mode core diameter of 2.5 μm in illumination, eliminated the speckle noise in the final image at the expense of spatial resolution due to reduced confocality.

Measurements in a Model Eye Using a Retina Mimicking Phantom
Uniformly scattering phantoms were manufactured using a method previously described by de Bruin et al. 47 Using mixtures with different wt. % of TiO 2 , layered phantoms were made in a Petri dish, each layer cured in a thermal oven at 80°C for 60 min. Thin metallic wires of different diameter thicknesses (90 and 140 μm) were embedded in the topmost scattering layer of the phantom while still in the viscous state and then cured. Then, a transparent layer of silicone without any TiO 2 particles was poured finally to provide mechanical stability to the phantom. After solidifying completely, the metallic wires were removed from the phantom. For our experiments, we used three-layer phantoms [ Fig. 10(a)] with a transparent top layer, a thin 200-μm layer with μ 0 s ¼ 1 mm −1 with channels in the middle, and a bottom layer with μ 0 s ¼ 5 mm −1 , respectively. Figure 10(b) shows an OCT B-scan (cross section) of the phantom used for the measurements in the model eye. Two water-soluble dyes-Evans Blue (Sigma Aldrich B.V., The Netherlands) and red food colorant (Lebenmittelfarbe, Birkmanns Voedelkleurstof, Germany) were used to mimic the deoxyhemoglobin and oxyhemoglobin absorptions, respectively. Using a suitable concentration of these dyes, an artificial isosbestic point was created at 548 nm as shown in Fig. 10(c). Using a syringe, different calibrated proportions of the dyes were injected into the channels. If [EB] is the concentration of Evans Blue and [FC] is the concentration of the red food colorant, we define a parameter ς as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 3 2 6 ; 7 1 9 ς ¼ ½FC ½FC þ ½EB : (19) An analysis similar to the analysis in Sec. 2 was performed to find the wavelengths resulting in low errors when estimating the parameter ς. The optimum wavelengths for the dye combination were found to be 500, 548 (isosbestic), and 610 nm, respectively. Due to the choice of the dichroic mirrors (DM1 and DM2 in Fig. 9) in the detection channel, 488 and 612 nm were chosen as the two anisosbestic wavelengths to measure the saturation. Using a suitable ratio of [FC] and [EB], solutions with five different ς (0.00, 0.25, 0.50, 0.75, and 1.00) were made and injected into the channels using a syringe. The phantoms were imaged simultaneously using the three different wavelengths, 488-nm [bandwidth, full width at half maximum (FWHM): 2 nm], 612-nm (bandwidth, FWHM: 3 nm), and 548-nm (isosbestic) (bandwidth, FWHM: 2 nm), and en face face images were created for each wavelength. The number of reflected photons from TiO 2 and the dye-filled channel was estimated from the image and the parameter ς was then calculated after estimating the ODs, ρ and D C (channel diameter) values from the images, respectively.  Figure 11 shows the typical intensity profile of a phantom blood vessel. The phantom blood vessel diameter was determined by the distance (in pixels) between the points with maximum slope magnitude along the blood vessel profiles v 1 and v 2 , respectively. The intensity at the blood vessel location was directly estimated by the minima of the vessel profile.

Estimating Vessel Diameter from the Images
The intensity at the tissue location t was estimated using the points t 1 and t 2 by a linear approximation. The points t 1 and t 2 represent the locations where the light has not interacted with the blood and are at a fixed number of pixels from v 1 and v 2 , respectively, such that the locations t and b are affected by similar factors except the absorption by hemoglobin. In our analysis, we took the points t 1 and t 2 , at a distance of 2 × diameter from on each side of v 1 and v 2 , respectively. To increase the accuracy of the estimated intensities (points t and b), the intensities were averaged along multiple locations in the tissue and in the middle of the blood vessel as shown in Fig. 11.
The intensity Iðx b ; λÞ was obtained by averaging the intensities of the point b along the valley (along with the y-dimension that is denoted by dotted lines). The point t and thereby Iðx t ; λÞ is obtained by averaging the points in an area as shown in Fig. 11. Figures 12(a)-12(c) show the reflectance measurements for 100% food coloring (ς ¼ 1.00) showing high absorption at 488 nm [ Fig. 12(a)]. Figures 12(d)-12(f) show the reflectance measurements for 0% food coloring (ς ¼ 0.00). Measurement points were chosen along the center of the channel and in the tissue (TiO 2 layer in this case) to estimate the reflected light intensity (I b and I t ). To bring the standard error on the mean intensity in both the tissue and along the blood vessels to acceptable values of ≤1%, up to 400 points were averaged in the tissue location and 40 points were averaged in the blood vessel location. The intensity profiles along the lines indicated in Figs. 12(d)-12(f) are shown in Fig. 12(g). The dip in the intensity profile due to the absorption of blood is obvious in Fig. 12(g). The channels in the phantom are slightly curved due to the process of manufacturing the phantoms as described in Sec. 3.2. The phantom was slightly angled during measurements to reduce specular reflections. This resulted in a slight change in the nonuniformity of the collected intensities within the imaging field. The bright spots in Figs. 12(a) and 12(d) are due to specular reflections from the surface of the sample and could be caused by dust on the surface.

Experimental Results with Retinal Phantoms
The comparison between the theoretical (with and without pigment packing correction) and experimentally obtained average ρ values in the phantom channel is shown in Fig. 13, for both the 90-and 140-μm channels used in the phantom. Our phantom geometry resembles the geometry in Fig. 1(a) here the pigment packing effect is not expected to play a large role. Indeed, Fig. 13 clearly shows that the estimated values closely follow the theoretical values without correction for pigment packaging. The estimated diameter of the channels from the images for 140 and 90 μm were 144 AE 6 μm and 92 AE 7 μm, respectively.
In Fig. 13, both the fitted curves with and without pigment packing describe the measurements reasonably well as the difference between the uncorrected and corrected calculated ρ values is small. This illustrates that by picking optimum wavelengths that fulfill the conditions stated in Sec. 2, the deviation between the ρ values with and without pigment packing could be minimized. For nonoptimal wavelengths, the difference in ρ values with and without correction is much higher (data not shown) and hence, using the simplified model presented without pigment packing would result in large offsets to the saturation estimation for smaller blood vessels embedded deeper in tissue. Figure 14 shows the pseudocolor map of the saturation overlaid on the corresponding reflectance image from the SLO for the 140-μm blood vessel for five different calibrated values. The estimated values show the overall mean and standard deviation of the saturation in the channel. For the estimation of ς, the 488-to 548-nm pair was used for ς ¼ 0.00; 0.25, and 0.50, whereas the 548-to 612-nm pair was used for ς ¼ 0.75 and 1.00, respectively.

Discussion
SLO-based retinal oximetry has advantages over conventional fundus camera-based method regarding resolution and contrast. A critical factor for clinically relevant retinal oximetry is able to determine blood oxygen saturation in small retinal vesselscapillaries, venules, and arterioles. It is in these microvessels that the oxygen saturation is expected to decrease in response to increased metabolic demand or decreased oxygen delivery capacity. The larger retinal vessels (>100 μm) are expected to be less sensitive to changes in tissue metabolic demand or microvascular dysfunction and are therefore not ideal as early hypoxia markers. The clinical value of fundus camera-based oximeters is most likely limited due to the low spatial resolution, even if the oxygen saturation estimates would be accurate and robust taking into account our optimized wavelengths and our proposed algorithms that include pigment packaging. SLObased oximeters have been investigated before, but their performance was limited due to suboptimal choices for the wavelengths used, combined with oximetry algorithms that did not take pigment packaging into account, which is likely a nonnegligible factor for smaller blood vessels located deeper in the retina. Journal of Biomedical Optics 086003-10 August 2018 • Vol. 23 (8) In this paper, for the first time, we have addressed factors that affect the accuracy of SLO-based oximeters, i.e., measurement noise and pigment packaging when measuring small, embedded blood vessels. Although we have used relatively large diameter blood vessels in our validation experiments, this was a limitation of the phantom fabrication technique yielding unstable phantoms at low channel diameters and is not a limitation of the imaging technique as a current SLO technology can image blood vessels <50 μm.
In our SLO design, we used a single-mode fiber illumination and multimode fiber detection. Hence the imaging is not truly confocal, but quasiconfocal or subdiffuse. We believe that quasiconfocal, subdiffuse detection reduces the central light reflex, a bright specular reflection from the center of a blood vessel that is common to many fundus imaging modalities that do not use cross polarizers. When the central light reflex is present in an image, the intensity in the center of the vessel could be extracted by reconstructing the blood vessel profile using analytical models, such as a fourth-order polynomial. 39,48 This adds additional error to the blood vessel OD estimation and thereby, the saturation estimation. Longitudinal chromatic aberration (LCA) is an another potential source of error in oximetry due to the dispersion in the human eye. 49 LCA causes the spatial location of each pixel in the multispectral image to be in a slightly different axial plane for different wavelengths. This effect can lead up to a 250-μm longitudinal focal shift between the 470-and 592-nm wavelengths; however, this effect can be corrected as has been done previously by others 50,51 using an achromatizing lens in the beam path. The optimum wavelength analysis was made with a 1-nm bandwidth around each central wavelength. The filter bandwidth choice can be incorporated into the oximetry analysis by multiplying the absorption spectrum of hemoglobin and the filter transmission spectrum and obtaining an average absorption coefficient for that spectral band, and then using this averaged absorption coefficients in Eqs. (9) and (17) to obtain a saturation estimation, and therefore the saturation error.
As explained previously, the noise in the images propagates into the saturation error. To achieve acceptable levels of saturation error (we aim for <5%, but it depends on the clinical application), a large number of points along the blood vessel and the tissue location have to be averaged. Although this is relatively easy in our experiments due to the homogeneous nature of the phantom, an in-vivo image contains structural information and choosing many points might lead to an incorrect estimation of the intensity in the absence of a blood vessel. A smaller number of points around the blood vessel should be averaged to overcome this problem, and this requires multiple images to be produced within a short amount of time. Fundus cameras are based on snapshot imaging, and it might not always be possible to acquire multiple images continuously without causing discomfort to the subject or exceeding the safety limit for radiant exposure. Although our current SLO images at relatively slow speeds, significant improvements in speed can be achieved using a resonant scanner, with a reduced SNR. Inhomogeneities in absorption within a vessel due to the distribution of red blood cells is less relevant in the final saturation estimation due to the required averaging and resulting image acquisition speeds to achieve desired SNR. Current commercial SLO technology has an imaging speed of about 30 Hz, which facilitates the averaging of multiple consecutive images to accomplish noise reduction.
The effective path length ðhL eff iÞ that photons travel through a tissue volume before being collected depends on the illumination and collection geometry and is also a function of the optical properties of the tissue and is, therefore, wavelength dependent. We have assumed that the effective path length for the two wavelengths used to calculate the factor ρ [Eq. (12)] is almost the same. In the wavelength range 450 to 600 nm, absorption remains in the same order and scattering varies slowly with wavelength (roughly as 1∕λ). However, if there is a slight mismatch in the path lengths, this results in an offset in the saturation estimation. We have calculated that a 20% mismatch in path lengths ðhL λ p i∕hL λ i i ¼ 1.2Þ could lead up to a 30% offset in the saturation estimates. For the phantoms that we fabricated, scattering was a slowly varying function of wavelength, similar to tissue, and the absorption coefficients of the dyes in the channels were of the same order of magnitude as the absorption coefficients of oxy-and deoxyhemoglobin. For our phantoms, we have shown the ability to accurately recover the saturation in the channels, suggesting that there were no significant changes in the path length due to scattering and absorption in a tissue-mimicking phantom with relevant scattering and absorbing properties. A limiting factor of our phantom was that the dyes that mimic the oxy-and deoxyhemoglobin were nonscattering, which is different from real blood.
Additionally, in a subdiffuse approach, pigments that are present in the retina, especially in the retinal pigment epithelium (RPE) and absorption background due to choroid blood, can affect the recorded intensity. The effect of the melanin on backscattered light varies with concentration, 52 and an a priori knowledge of the pigment concentration in RPE melanin and its contribution to the backscattered light can aid in removing the influence of RPE in the backscattered light. OCT-based oximetry methods might have an advantage here that layers of the retina free from the influence of RPE can be used for extracting the intensities. But, OCT being a coherent detection method suffers from speckles, and a significant amount of averaging has to be performed to improve the image quality and reduce the error on intensity down to 1%. We must emphasize that the insufficient removal of the system reflections from the recorded intensities will result in a wrong estimation of the ODs and hence the saturation. Therefore, in our experiments, all the internal system reflections were subtracted by taking a reference measurement without any sample.
There have been several attempts to choose the optimal wavelengths for retinal oximetry, but none of them have accounted for the hemoglobin packaging effect. Most of these attempts were aimed at two wavelength oximetries 15,53,54 using one isosbestic wavelength and another nonisosbestic wavelength, where jμ HbO 2 a ðλÞ − μ Hb a ðλÞj is maximum, the popular choice being 570 (isosbestic) and 600 nm, respectively. Although this wavelength combination might give sufficiently accurate saturation estimation in arteries, our analysis shows that they are not reliable for low saturation values of S < 60%. Schweitzer et al. 55 reported a physiologically relevant retinal oxygenation range of 57% to 92%. However, lower saturations may be expected due to various pathological conditions. Smith et al. 56 chose 803 (isosbestic) and 670 nm to estimate the saturation. The infrared wavelengths have weak blood absorption and therefore give a poor estimate of the blood absorption leading to a large error in the saturation estimations according to our calculations. In 1988, Delori 14 used three wavelengths-558, 570 (isosbestic), and 548 nm (isosbestic)-to scan a small area of the retina with blood vessel tracking to reduce eye motion artefacts. This wavelength choice, similar to the 570to 600-nm wavelength pair, results in poor accuracy for low saturation values, S <60%. Our assessment of the modified oximetry equations suggests that the 470, 506 (isosbestic), and 592 nm triad is the best combination for oximetry. Although this choice spans a relatively broad wavelength range ≈130 nm, the absorption coefficients and scattering coefficients in this wavelength range are similar to those we have used in our phantoms, suggesting that the effective path length differences within this wavelength range are of minor importance.
The accuracy of saturation estimation can be potentially improved by including more wavelengths. To this end, hyperspectral imaging 17,[23][24][25]57 approaches toward in vivo retinal oximetry have also been undertaken by various groups. However, the challenge of acquiring images with sufficient SNR from multiple wavelengths at high spatial resolution, in a large field of view, with high speed and with low power levels remains.

Conclusions
We have presented a method to select optimal wavelengths for oximetry measurements in the retina through an error analysis on a modified oximetry equation that includes the effect of pigment packaging. The result of the analysis yielded 470, 506 (isosbestic), and 592 nm as three suitable wavelengths for accurate estimation of oxygenation in the retinal blood vessels. We have developed a quasiconfocal, multispectral SLO capable of simultaneously imaging the retina at all these wavelengths and demonstrated the validity of our approach in tissue-mimicking phantoms. In vivo human experiments are pending for ethical committee approval.