|Home | About | Journals | Submit | Contact Us | Français|
The level of retinal oxygenation is potentially an important cue to the onset or presence of some common retinal diseases. An improved method for assessing oxygen tension in retinal blood vessels from phosphorescence lifetime imaging data is reported in this paper. The optimum estimate for phosphorescence lifetime and oxygen tension is obtained by regularizing the least-squares (LS) method. The estimation method is implemented with an iterative algorithm to minimize a regularized LS cost function. The effectiveness of the proposed method is demonstrated by applying it to simulated data as well as image data acquired from rat retinas. The method is shown to yield estimates that are robust to noise and whose variance is lower than that obtained with the classical LS method.
Abnormalities in retinal oxygenation are thought to play a significant role in the development of common eye diseases, such as diabetic retinopathy, glaucoma, and age-related macular degeneration (AMD). However, our knowledge of oxygen’s function in many retinal diseases is incomplete. Accurate assessment of oxygen tension in the retinal vasculature is therefore important to establish the role of oxygen in the development of retinal disease and aid in the onset or presence of hypoxia. Oxygen tension can be estimated using the technique of phosphorescence lifetime imaging –. A novel system for optical section phosphorescence imaging  was devised in our earlier work for acquiring image data to estimate phosphorescence lifetimes and oxygen tension (pO2).
The estimation methods for phosphorescence lifetime imaging in the retinal tissue – rely on a linear model that was developed earlier by Lakowicz et al. for the purpose of estimating fluorescence lifetime values .
In the studies carried out in –, the classical least-squares (LS) method was used to calculate model parameters. Although the method is computationally efficient, it suffers from high variance in the estimation of values, and produces spikes in pO2 maps, which can be attributed to the lack of an appropriate prior model for the data. To overcome these shortcomings and utilize knowledge of the prior distribution of the parameters, we propose the use of regularization of the LS method with appropriate signal model assumptions. Regularization has been extensively used in inverse problems in different fields such as image denoising, astronomical imaging, reconstruction, and synthetic aperture radar (SAR) imaging –. Rudin et al.  proposed an image denoising algorithm based on total-variation regularization. The image was denoised by minimizing total variation norm of the estimated solution. The method is suitable for images corrupted with additive Gaussian noise. In performing regularization, it is critical to properly choose the regularization parameter according to the noise variance. Galatsanos and Katsaggelos  proposed new approaches for choosing the regularization parameter and estimating noise variance in the images. Kang and Katsaggelos  proposed an iterative image restoration method by using an adaptive regularization parameter. Their method updates the regularization parameter at each step of the iteration.
Regularization was also applied in the transform domain. Belge et al.  proposed an edge-preserving regularization method performed after applying a wavelet transform to the image. They modeled the wavelet coefficients of the desired image as generalized Gaussian and applied regularization in the wavelet domain. Regularization applied to real-valued images was adapted to complex-valued SAR image data. Cetin and Karl  modeled SAR image reconstruction as an optimization problem by imposing nonquadratic constraints on the cost function. Their model for prior information constraint consisted of the sum of two nonquadratic functions of the absolute value and the derivative of the desired signal that was used to control resolution and smoothness simultaneously. Yu and Fessler  proposed an edge-preserving regularization method for tomographic image reconstruction. Their method was driven by the need to include nonlocal boundary information into the regularization method, unlike other methods that use only local boundary information to avoid smoothing of edges. Mignotte  used regularization for restoration of images degraded by space-invariant blur and additive Gaussian noise. Regularization has also been applied in an iterative restoration process with an adaptive regularization parameter based on the prior segmentation information. Le et al.  proposed an image denoising method based on regularization. They used total-variation regularization for denoising as in , but modeled noise as Poisson and showed the superiority of their method when compared with the method in  for signal-dependent noise with Poisson distribution. The approach we propose here is also one of regularizing the LS estimation in spatial domain.
This extensive use of regularization methods for improved estimation in different image analysis applications provided a strong motivation for us to examine the usefulness of such methods in retinal oxygen tension estimation, for reducing the effect of noise obtained with the classical LS method. In our regularized LS (RLS) method, the optimum solution is obtained with the use of an iterative procedure to minimize a suitably defined RLS cost function. The effectiveness of the RLS method is demonstrated by applying it to simulated data as well as real image data acquired from rat retinas. The method is shown to be robust to the presence of noise, and it yields estimates that are more in the physiologically expected range and whose variance is lower than that obtained with classical regression analysis. Some preliminary results of our study were reported in .
The rest of the paper is organized as follows. In Section II, we describe the method of our data acquisition and the model used in the estimation. The maximum likelihood (ML) estimation method, which reduces to the LS method under certain noise model assumption, is briefly explained in Section III. The MAP estimation and its application to our problem are also explained in Section III. In Section IV, results obtained using the LS and RLS methods are compared and discussed. Results show improved performance of the RLS method in estimating the oxygen tension maps. Variance is significantly reduced without altering the mean when compared with the classical approach. Improvement is demonstrated for both experimental data, acquired from rat retinas, and simulated data. We conclude by discussing the results and describing some issues for future investigation in Section V.
A novel system for optical section phosphorescence imaging ,  developed in our earlier work is briefly described. A diode laser at a wavelength of 532 nm was expanded to a line with a cylindrical lens and focused at an oblique angle on the retina following intravenous injection of an oxygen-sensitive molecular probe. The imaging optics of a slit-lamp biomicroscope, fitted with an infrared filter with transmission overlapping the phosphorescence emission of the oxygen probe, was utilized to form an image of the retina on a digital intensified charge-coupled device (ICCD) camera (Roper Scientific, Trenton, NJ). Since the incident laser beam was not coaxial with the imaging axis of the instrument, retinal features at different depths were laterally displaced on the phosphorescence optical section retinal image according to their depth location. A galvanometer scanner was integrated into the system for scanning the laser line horizontally on the retina. A series of optical section phosphorescence images from spatially adjacent locations were acquired. The optical section phosphorescence images were then combined to generate en face phosphorescence images of the retinal vasculatures in the retinal plane using an image reconstruction technique described previously . During image acquisition, the excitation laser beam was modulated with an optical chopper (Photon Technology, London, ON, Canada) positioned in the beam path and operated at a frequency of 1600 Hz. The digital output from the chopper triggered the timing controller of the ICCD camera. At each location, ten phase-delayed optical section phosphorescence images were acquired with delay increments of 74 µs.
A total of five male Long Evans pigmented rats (450–650 g) were used for this study. The animals were treated in compliance with the Association for Research in Vision and Ophthalmology (ARVO) Statement for the Use of Animals in Ophthalmic and Vision Research. An intraperitoneal (IP) infusion of ketamine (85 mg/kg IP) and xylazine (3.5 mg/kg IP) at the rate of 0.5 and 0.02 mg/(kg·min), respectively, was used to anesthetize the rats. The percentage of inspired oxygen was controlled by means of a high-flow mask system. Gas mixtures containing 21% oxygen (room air, normoxia) and 10% oxygen (hypoxia) were administered to rats for 5 min before and during imaging. For measurements of systemic arterial pO2, the femoral artery was cannulated and arterial blood (200 µL) was drawn through the catheter without air exposure. Systemic arterial pO2 was measured with use of a blood gas analyzer (Radiometer, Westlake, OH).
To image the retinal vasculatures, pupils were dilated with 2.5% phenylephrine and 1% tropicamide. In order to eliminate the cornea refractive power and prevent dehydration, a glass cover slip was placed on the cornea after application of 1% hydroxypropyl methylcellulose. An oxygen-sensitive molecular probe, Pd-porphine (Frontier Scientific, Logan, UT), was dissolved (12 mg/mL) in bovine serum albumin solution (60 mg/mL), physiological saline buffered to a pH of 7, and injected intravenously (20 mg/kg). Before imaging, the rat was positioned in front of the imaging instrument and the laser power was adjusted to 100 µW, which is safe for viewing according to the American National Standard Institute for Safety Standards. Imaging was performed at areas within two-disk diameters (600 µm) from the edge of the optic nerve head.
The variation of the intensity of the emitted phosphorescence, denoted by I(t), can be expressed as a function of time as follows:
where I0 is the maximum intensity at time zero and τ is the lifetime. The relationship between the phosphorescence lifetime and the quenching agent, which is the oxygen tension pO2 in this study, is given by the Stern–Volmer expression
where pO2 (in millimeters of mercury) is the oxygen tension, Kϕ (in inverse of millimeters of mercury·microseconds) is the quenching constant, and τ0 is the lifetime in zero oxygen environments. Based on (2), our goal is to estimate pO2, given constants τ0 and Kϕ. To estimate pO2 values, we need to determine τ in (2). As previously described , , τ can be calculated from the distribution of intensity values in different modulation phase images, since the intensity depends on the phase angle between the modulated excitation laser light and the emitted phosphorescence, denoted by θ. The relation between θ and τ is
where ω is the modulation frequency.
The intensity I(t) depends on the concentration [Pd] of the probe, k, θ, and the modulation m, which are unknown. With some abuse of notation, we represent intensity as a function of θn 
where mn is the modulation profile of the image intensifier, θn denotes the phase of the gain modulation, and np is the number of measurements at each pixel location for different phase values θn .
Defining a0 = k[Pd], a1 = (1/2)k[Pd]mn m cos(θ), and b1 = (1/2)k[Pd]mn m sin(θ), (6) can be expressed as
Noting that b1/a1 = tan(θ), the phase angle θ in (6) is obtained as
The intensity values in (7) are shown in the absence of noise. In the presence of additive noise, the observed intensity vector I̲i for the ith pixel can be written in a matrix form as follows:
In order to determine the three unknowns a0 , a1 , and b1 in (9), we need at least three intensity values obtained at three different gain modulation phases for each pixel. To obtain more reliable pO2 maps in our experiments, ten intensity values (np = 10) were obtained at each pixel for the instrument phase delays θn equal to (n − 1)(0.24π). The optimum solution for the ith pixel in the LS sense is obtained from
where Q is the pseudoinverse of the matrix in (9). The estimation of parameters a1 and b1 is considered next.
We first describe the mathematical model for estimation of the parameters. In the presence of additive noise, the observed intensity vector can be expressed as
where X̲ ¡M denotes the vector of M parameters of interest, Y̲ ¡L is the vector of measurements, N̲ ¡L is the vector of the additive noise, and S̲ ¡L×M is the system matrix defined according to the imaging model in (7). Note that in our analysis, vector Y̲ is the concatenation of all pixel-wise observation vectors I̲i defined in (9) at each observation location and X̲ is the concatenation of all pixel-wise parameter vectors shown in (9). If J is the total number of pixels, then M = 3J and L = np J in our analysis. We wish to estimate the parameter vector X̲ from the observation vector Y̲. We briefly review the ML and MAP estimation methods and describe the RLS method in the subsequent sections.
ML estimation yields the value of Y̲ that is most likely to have been produced by X̲ . The ML estimate of X̲ is obtained by maximizing the conditional probability distribution function of Y̲ given X̲, which is also called the likelihood function
If the additive noise consists of independent identically distributed (i.i.d.) Gaussian samples with zero mean, then the ML estimate of X̲ becomes
where l(X̲) is the log-likelihood function of X̲ and X̲ML is equivalent to the LS solution. If the additive noise samples are independent and Gaussian, but not identically distributed, then X̲ML becomes equivalent to the weighted LS estimation, with the weights equal to noise variances. The ML method provides an estimate assuming a uniform prior distribution for the desired parameters. In practice, it is often not an appropriate assumption; therefore, the ML estimate is more vulnerable to noise than the MAP estimate, which utilizes the prior knowledge of the distribution of the signal to diminish the effect of noise on the estimates.
If prior knowledge of the distribution of the desired parameters is available, then MAP estimation improves upon the performance of ML. These two estimators become equivalent when the prior distribution is uniform .
The MAP estimate is obtained by maximizing the a posteriori probability density function, the conditional distribution of X̲ given Y̲
Taking the logarithm and omitting the terms that do not depend on X̲, we obtain
In our problem, we model the distribution of X̲ as multivariate Gaussian with uncorrelated random variables that have the same variance, and noise that is i.i.d. with multivariate Gaussian distribution. Under these assumptions, the MAP estimation becomes
where β is equal to the ratio of noise and signal variance. The first term in the cost function in (16) is called the data fidelity term and the second term is the regularization term. The weights of these two terms are controlled by the regularization parameter β that can theoretically vary over (0, ∞). As noise variance increases, the weight of the data fidelity term ideally decreases, and the estimation gives more weight to the regularization term, which is sometimes called the smoothness term. The regularization term leads to some bias in the estimate while decreasing the variance in the estimate ; therefore, there exists a bias-variance tradeoff, which is controlled by β in (16).
We are interested only in parameters a1 and b1 , since only the ratio of values of a1 and b1 is needed to calculate the pO2 values. Parameters a1 and b1 will be estimated by defining a cost function as follows:
where I̲i is the observed intensity vector, is the mean value of the parameter to be estimated for the ith pixel, A̲ is the second row of the matrix Q in (10), β is the regularization parameter, and J is the number of pixels in the image. A similar cost function can be defined for the parameter b1 when A̲ would be modified to be the third row of the matrix Q in (10).
In our algorithm, the mean in (17) is replaced with the sample mean computed in a neighborhood, which is chosen in our implementation as a 3 × 3 window. The choice of the size of the neighborhood was inspired by the fact that the variation among the pO2 values in a neighborhood of this size within the retinal arteries and veins should physiologically be minimal. We also used this fact when choosing the regularization parameter. If β is not chosen sufficiently large, then the RLS method results in high variance in the estimated pO2 values. On the other hand, if β is chosen too large, then the local variation of interest in the pO2 maps is lost due to excessive smoothing.
In (17), the cost function at each pixel location is dependent on the values of its neighboring pixels due to our definition of the mean value of the parameter. Therefore, there is no pixel-wise isolated solution, and a joint optimization is performed for the overall cost. The total cost is equal to the sum of the individual costs of each pixel in (17), and it is expressed as
where, as before, J denotes the total number of pixels.
The minimum of this cost function does not have a closed-form solution. Therefore, a gradient-based iterative approach is applied to find the minimum. The iterative algorithm starts with an initial choice of the parameter values, which, in our implementation, is the classical LS estimate. Next, the gradient of the cost equation is computed and used to determine the search direction for the unknown parameters. The updated value is then obtained as follows:
where α(k) is step size at the kth step that is found with a line search and F is the estimate of the gradient of the cost function.
In order to illustrate improvement in the performance due to the RLS method, we used both simulated and experimental datasets. Experimental data were acquired in animals under hypoxic (10%) as well as normal (21%) oxygen breathing conditions.
Simulated data were generated to mimic features of the experimental data. The pO2 maps and intensity images of simulated data were created based on the real image model, and performance can be better assessed due the availability of the true synthesized pO2 maps. Gaussian noise at different SNR values was added to synthesized intensity images. We then compared the mean absolute error (MAE) computed for the LS and RLS methods. To compare bias and variance performances of the two methods, we used Monte Carlo simulations. The pO2 maps were computed for 100 repetitive simulations, and the bias and variance of these 100 estimates were calculated.
Fig. 1 shows oxygen tension (pO2) maps in a rat under 21% oxygen breathing condition that is generated using the LS and RLS methods, respectively. In this map, three vessels are shown. The outer two vessels are retinal veins and the middle vessel is an artery. The map generated by the LS method produces spikes in the map due to the presence of noise and results in high variance, whereas the RLS method suppresses the spikes while not producing a significant bias, which allows reliable analysis of local variation in retinal vascular oxygenation.
The pO2 measurements estimated by the LS and RLS methods in the retinal vasculatures, and the systemic arterial pO2 for each rat are listed in Table I (21% oxygen breathing condition) and Table II (10% oxygen breathing condition). The mean pO2 values estimated in retinal arteries and veins under different breathing conditions by both methods do not differ significantly (p > 0.4; N = 10), which shows that the RLS method does not produce a noticeable bias.
Small pixel regions with the size of 2 × 2 pixels (10 × 10 µm2 ) were chosen and the sample standard deviation of pO2 values for each region over the arteries and veins was computed. The average of the sample standard deviation of pO2 values in retinal vessels calculated using both methods for 21% and 10% oxygen breathing conditions are given in Tables III and andIV,IV, respectively. In animals, pO2 variance in retinal arteries and veins generated by the LS method was (8 ± 4 mmHg; N = 10), and significantly (p < 0.0001) higher than the variance measured by the RLS method (2 ± 1 mmHg; N = 10) under 21% breathing condition. Similarly, under 10% breathing condition, pO2 variance in retinal blood vessels was significantly lower (p < 0.0001) with the RLS method.
To evaluate the accuracy of the RLS method, a relationship between systemic and retinal arterial pO2 was determined using linear regression analysis. The systemic and retinal arterial pO2 measurements were highly correlated (R = 0.9; p = 0.003; N = 10). The slope of the regression line was 0.6, indicating that retinal arterial pO2 was lower than the systemic arterial pO2, which was in agreement with previously reported results . This observed pO2 reduction is likely due to tissue oxygen consumption that reduces the oxygen content of the blood.
In order to show robustness of the RLS method in the presence of different levels of noise, a synthetic pO2 map was generated with certain assumed values of the parameters a1 and b1 . Based on the mathematical model, which calculates these parameters from the intensity values, intensity images of the simulated pO2 map were generated from a0 , a1 , and b1 parameters associated with the simulated pO2 map. In the noise-free case, the LS method can estimate the pO2 map from the intensity images perfectly as expected, but as noise at increasing levels is added to the image, the difference between the original pO2 values and those estimated by the LS method increases.
The simulated pO2 map and its estimation by the LS and RLS methods in the presence of noise at 20 dB SNR are shown in Fig. 2. The simulated pO2 map was generated in a manner similar to the procedure described earlier by generating three large vessels: two retinal veins and one retinal artery placed between the veins in which pO2 values are expected to be higher. The RLS method produces estimates very close to the original values, whereas the pO2 map of the LS method shows high variability with a marked difference from the original one, and it produces a larger error.
After simulating 100 repeated independent trials, the performance of the two methods were compared, and the results for the bias and variance are shown in Figs. 3 and and4,4, respectively. The RLS method yields a variance that is much smaller than that of the LS method, whereas there is no significant difference in the bias performance of the two methods.
For different SNR values, MAE of the pO2 maps estimated by using the RLS and LS methods are shown in Fig. 5. Instead of mean square error, we preferred MAE since it is more robust to outliers. The SNR range in horizontal axis was from 15 to 45 dB. As SNR increases, the LS method starts to give improved performance, and the difference in the MAE of the two methods decreases as expected.
We proposed an RLS method to estimate pO2 values of retinal vessels with a phosphorescence lifetime imaging model. The method decreased the variance of the classical regression analysis in the generation of pO2 maps of retinal vessels while preserving mean values. This was illustrated using experimental and simulated data under different levels of noise.
From the results given in Figs. 1–5 and Tables I–IV, it is observed that the RLS method produces more reliable pO2 estimates with significantly reduced variance while preserving the mean pO2 values of the retinal vessels.
In this study, only one dataset per animal was collected; therefore, reproducibility or intrasubject variability of measurements is not reported. However, we have previously published the reproducibility of oxygen tension measurements to be 2–3 mmHg in the artery and vein, under 21% and 10% breathing conditions . Based on the data in the current study, the intersubject variability in artery and vein oxygen tension measurements were comparable using the LS and RLS methods, ranging between 3 and 9 mmHg. Since oxygen tension is calculated based on phosphorescence lifetime, and not on intensity, changes in phosphorescence intensity due to pigmentation or blood vessel diameter are not expected to influence the measurements.
Abnormalities in retinal oxygen delivery and consumption are thought to play a significant role in retinal diseases, such as diabetic retinopathy and vascular occlusions. Application of this method for improved oxygen tension measurements in the retinal vasculatures will provide information on fundamental mechanisms that implicate hypoxia in the development of retinal pathologies. This information is greatly needed to broaden knowledge of disease pathophysiology, and thereby advance diagnostic and therapeutic procedures.
In future work, we plan to modify our method to estimate pO2 values in microvascular structures of retina. Since microvascular data contain very small features separated by sharp edges, the method must be made data-adaptive and appropriate preprocessing must be performed.
This work was supported by the National Institutes of Health (NIH) under Grant EY01792 (Department Core). The work of M. Shahidi was supported in part by the NIH under Grant EY017918. The work of I. Yildirim was supported by Istanbul Technical University.
Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.
Isa Yildirim, Department of Electrical and Computer Engineering, University of Illinois, Chicago, IL 60607 USA, and also with the Informatics Institute, Istanbul Technical University, 34469 Istanbul, Turkey (Email: ude.ciu.ece@iridliyi).
Rashid Ansari, Department of Electrical and Computer Engineering, University of Illinois, Chicago, IL 60607 USA (Email: ude.ciu@irasnar).
Justin Wanek, Department of Ophthalmology and Visual Sciences, University of Illinois, Chicago, IL 60607 USA (Email: ude.ciu@kenawmj).
Imam Samil Yetik, Medical Imaging Research Center, Illinois Institute of Technology, Chicago, IL 60616 USA (Email: ude.tii@kitey).
Mahnaz Shahidi, Department of Ophthalmology and Visual Sciences, University of Illinois, Chicago, IL 60607 USA (Email: ude.ciu@hahsnham).