Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2885291

Formats

Article sections

- Abstract
- I. INTRODUCTION
- II. DATA ACQUISITION AND MODEL
- III. ML AND MAP ESTIMATION METHODS
- IV. RESULTS OF PROCESSING REAL AND SIMULATED DATA
- V. CONCLUSION AND DISCUSSION
- REFERENCES

Authors

Related links

IEEE Trans Biomed Eng. Author manuscript; available in PMC 2010 June 14.

Published in final edited form as:

Published online 2009 April 21. doi: 10.1109/TBME.2009.2020505

PMCID: PMC2885291

NIHMSID: NIHMS206125

Isa Yildirim, IEEE, Student Member,^{} Rashid Ansari, IEEE, Fellow, Justin Wanek, Imam Samil Yetik, IEEE, Senior Member, and Mahnaz Shahidi

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

The publisher's final edited version of this article is available at IEEE Trans Biomed Eng

See other articles in PMC that cite the published article.

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 [1]–[3]. A novel system for optical section phosphorescence imaging [1] was devised in our earlier work for acquiring image data to estimate phosphorescence lifetimes and oxygen tension (*p*O_{2}).

The estimation methods for phosphorescence lifetime imaging in the retinal tissue [1]–[3] rely on a linear model that was developed earlier by Lakowicz *et al.* for the purpose of estimating fluorescence lifetime values [4].

In the studies carried out in [1]–[3], 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 *p*O_{2} 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 [5]–[12]. Rudin *et al.* [5] 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 [6] proposed new approaches for choosing the regularization parameter and estimating noise variance in the images. Kang and Katsaggelos [7] 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.* [8] 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 [9] 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 [10] 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 [11] 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.* [12] proposed an image denoising method based on regularization. They used total-variation regularization for denoising as in [5], but modeled noise as Poisson and showed the superiority of their method when compared with the method in [5] 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 [13].

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 [1], [14] 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 [14]. 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 *p*O_{2}, the femoral artery was cannulated and arterial blood (200 µL) was drawn through the catheter without air exposure. Systemic arterial *p*O_{2} 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:

$$I(t)={I}_{0}\phantom{\rule{thinmathspace}{0ex}}\text{exp}\phantom{\rule{thinmathspace}{0ex}}\left(\frac{-t}{\tau}\right)$$

(1)

where *I*_{0} 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 *p*O_{2} in this study, is given by the Stern–Volmer expression

$$\frac{{\tau}_{0}}{\tau}=1+{K}_{}$$

(2)

where *p*O_{2} (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 *p*O_{2}, given constants τ_{0} and *K*_{ϕ}. To estimate *p*O_{2} values, we need to determine τ in (2). As previously described [3], [4], τ 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

$$\text{tan}\phantom{\rule{thinmathspace}{0ex}}\theta =\omega \tau $$

(3)

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} [4]

$$I({\theta}_{n})=k[\text{Pd}]\phantom{\rule{thinmathspace}{0ex}}\left(1+\frac{1}{2}{m}_{n}\phantom{\rule{thinmathspace}{0ex}}m\phantom{\rule{thinmathspace}{0ex}}\text{cos}\{\theta -{\theta}_{n}\}\right)\phantom{\rule{thinmathspace}{0ex}},\text{}n=1,2,\dots ,{n}_{p}$$

(4)

where *m _{n}* is the modulation profile of the image intensifier, θ

Using the trigonometric identity,

$$\text{cos}(\theta -{\theta}_{n})=\text{cos}(\theta )\phantom{\rule{thinmathspace}{0ex}}\text{cos}({\theta}_{n})+\text{sin}(\theta )\phantom{\rule{thinmathspace}{0ex}}\text{sin}({\theta}_{n})$$

(5)

the intensity values in (4) can be written as [4]

$$I({\theta}_{n})=k[\text{Pd}]+\frac{1}{2}k[\text{Pd}]{m}_{n}\phantom{\rule{thinmathspace}{0ex}}m(\text{cos}(\theta )\phantom{\rule{thinmathspace}{0ex}}\text{cos}({\theta}_{n})+\text{sin}(\theta )\phantom{\rule{thinmathspace}{0ex}}\text{sin}({\theta}_{n})).$$

(6)

Defining *a*_{0} = *k*[Pd], *a*_{1} = (1/2)*k*[Pd]*m _{n} m* cos(θ), and

$$I({\theta}_{n})={a}_{0}+{a}_{1}\phantom{\rule{thinmathspace}{0ex}}\text{cos}({\theta}_{n})+{b}_{1}\phantom{\rule{thinmathspace}{0ex}}\text{sin}({\theta}_{n}).$$

(7)

Noting that *b*_{1}/*a*_{1} = tan(θ), the phase angle θ in (6) is obtained as

$$\theta ={\text{tan}}^{-1}\phantom{\rule{thinmathspace}{0ex}}\left(\frac{{b}_{1}}{{a}_{1}}\right).$$

(8)

The *p*O_{2} values at each observation location can, therefore, be obtained from *a*_{1} and *b*_{1} using (2), (3), and (8).

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

$${[\begin{array}{ccc}{I}_{1}^{i}\hfill & {I}_{2}^{i}\hfill & {I}_{{n}_{p}}^{i}\hfill \hfill \\ ]\end{array}T}^{=}$$

(9)

In order to determine the three unknowns *a*_{0} , *a*_{1} , and *b*_{1} in (9), we need at least three intensity values obtained at three different gain modulation phases for each pixel. To obtain more reliable *p*O_{2} maps in our experiments, ten intensity values (*n _{p}* = 10) were obtained at each pixel for the instrument phase delays θ

$${[\begin{array}{ccc}{a}_{0}^{i}\hfill & {a}_{1}^{i}\hfill & {b}_{1}^{i}\hfill \end{array}]}^{T}=Q{\underset{\xaf}{I}}^{i}$$

(10)

where *Q* is the pseudoinverse of the matrix in (9). The estimation of parameters *a*_{1} and *b*_{1} 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

$$\underset{\xaf}{Y}=S\underset{\xaf}{X}+\underset{\xaf}{N}$$

(11)

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

ML estimation yields the value of *Y̲* that is most likely to have been produced by *X̲* [15]. 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

$${\underset{\xaf}{X}}_{\text{ML}}=\underset{\underset{\xaf}{X}}{\text{arg max}}\phantom{\rule{thinmathspace}{0ex}}p\phantom{\rule{thinmathspace}{0ex}}\left(\frac{\underset{\xaf}{Y}}{\underset{\xaf}{X}}\right)\phantom{\rule{thinmathspace}{0ex}}.$$

(12)

If the additive noise consists of independent identically distributed (i.i.d.) Gaussian samples with zero mean, then the ML estimate of *X̲* becomes

$${\underset{\xaf}{X}}_{\text{ML}}=\underset{\underset{\xaf}{X}}{\text{arg max}}\phantom{\rule{thinmathspace}{0ex}}l(\underset{\xaf}{X})=\underset{\underset{\xaf}{X}}{\text{arg min}}\phantom{\rule{thinmathspace}{0ex}}{\Vert \underset{\xaf}{Y}-S\underset{\xaf}{X}\Vert}_{2}^{2}$$

(13)

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 [16].

The MAP estimate is obtained by maximizing the *a posteriori* probability density function, the conditional distribution of *X̲* given *Y̲*

$${\underset{\xaf}{X}}_{\text{MAP}}=\underset{\underset{\xaf}{X}}{\text{arg max}}\frac{p(\underset{\xaf}{Y}/\underset{\xaf}{X})p(\underset{\xaf}{X})}{p(\underset{\xaf}{Y})}.$$

(14)

Taking the logarithm and omitting the terms that do not depend on *X̲*, we obtain

$${\underset{\xaf}{X}}_{\text{MAP}}=\underset{\underset{\xaf}{X}}{\text{arg max}}\phantom{\rule{thinmathspace}{0ex}}\text{log}\phantom{\rule{thinmathspace}{0ex}}\left(p\phantom{\rule{thinmathspace}{0ex}}\left(\frac{\underset{\xaf}{Y}}{\underset{\xaf}{X}}\right)\right)+\text{log}(p(\underset{\xaf}{X})).$$

(15)

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

$$\underset{\underset{\xaf}{X}}{\text{arg min}}\phantom{\rule{thinmathspace}{0ex}}\{{(\underset{\xaf}{Y}-S\underset{\xaf}{X})}^{T}\phantom{\rule{thinmathspace}{0ex}}(\underset{\xaf}{Y}-S\underset{\xaf}{X})+\beta {(\underset{\xaf}{X}-\underset{\xaf}{{\mu}_{X}})}^{T}\phantom{\rule{thinmathspace}{0ex}}(\underset{\xaf}{X}-\underset{\xaf}{{\mu}_{X}})\}$$

(16)

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 [17]; therefore, there exists a bias-variance tradeoff, which is controlled by β in (16).

We are interested only in parameters *a*_{1} and *b*_{1} , since only the ratio of values of *a*_{1} and *b*_{1} is needed to calculate the *p*O_{2} values. Parameters *a*_{1} and *b*_{1} will be estimated by defining a cost function as follows:

$${f}_{{a}_{1}}^{i}={({a}_{1}^{i}-{\underset{\xaf}{AI}}^{i})}^{2}+\beta \phantom{\rule{thinmathspace}{0ex}}{\left({a}_{1}^{i}-\overline{{a}_{1}^{i}}\right)}^{2},\text{}i=1,2,\dots ,J$$

(17)

where *I̲ ^{i}* is the observed intensity vector, $\overline{{a}_{1}^{i}}$ is the mean value of the parameter to be estimated for the

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 *p*O_{2} 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 *p*O_{2} values. On the other hand, if β is chosen too large, then the local variation of interest in the *p*O_{2} 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

$${F}_{{a}_{1}}={\displaystyle \sum _{i=1}^{J}{f}_{{a}_{1}}^{i}}$$

(18)

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:

$${a}_{1}\phantom{\rule{thinmathspace}{0ex}}(k+1)={a}_{1}\phantom{\rule{thinmathspace}{0ex}}(k)-\alpha (k)F({a}_{1}(k))$$

(19)

where α(*k*) is step size at the *k*th 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 *p*O_{2} 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 *p*O_{2} 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 *p*O_{2} maps were computed for 100 repetitive simulations, and the bias and variance of these 100 estimates were calculated.

Fig. 1 shows oxygen tension (*p*O_{2}) 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.

Estimated *p*O_{2} values in retinal artery and veins of experimental data using the (left) LS and (right) RLS methods. The color bar represents estimated *p*O_{2} values in millimeters of mercury.

The *p*O_{2} measurements estimated by the LS and RLS methods in the retinal vasculatures, and the systemic arterial *p*O_{2} for each rat are listed in Table I (21% oxygen breathing condition) and Table II (10% oxygen breathing condition). The mean *p*O_{2} 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 µm^{2} ) were chosen and the sample standard deviation of *p*O_{2} values for each region over the arteries and veins was computed. The average of the sample standard deviation of *p*O_{2} 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, *p*O_{2} 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, *p*O_{2} variance in retinal blood vessels was significantly lower (*p* < 0.0001) with the RLS method.

Mean of Standard Deviation (SD) of *p*O_{2} Values in Retinal Vasculatures at 21% Oxygen Breathing Condition

Mean of Standard Deviation (SD) of *p*O_{2} Values in Retinal Vasculatures at 10% Oxygen Breathing Condition

To evaluate the accuracy of the RLS method, a relationship between systemic and retinal arterial *p*O_{2} was determined using linear regression analysis. The systemic and retinal arterial *p*O_{2} 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 *p*O_{2} was lower than the systemic arterial *p*O_{2}, which was in agreement with previously reported results [18]. This observed *p*O_{2} 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 *p*O_{2} map was generated with certain assumed values of the parameters *a*_{1} and *b*_{1} . Based on the mathematical model, which calculates these parameters from the intensity values, intensity images of the simulated *p*O_{2} map were generated from *a*_{0} , *a*_{1} , and *b*_{1} parameters associated with the simulated *p*O_{2} map. In the noise-free case, the LS method can estimate the *p*O_{2} map from the intensity images perfectly as expected, but as noise at increasing levels is added to the image, the difference between the original *p*O_{2} values and those estimated by the LS method increases.

The simulated *p*O_{2} 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 *p*O_{2} 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 *p*O_{2} values are expected to be higher. The RLS method produces estimates very close to the original values, whereas the *p*O_{2} map of the LS method shows high variability with a marked difference from the original one, and it produces a larger error.

(Left) Original simulated *p*O_{2} map (60 × 50 pixels), and its estimates in the presence of noise with 20 dB SNR using the (middle) LS and (right) RLS methods, respectively. The color bar represents *p*O_{2} values.

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.

Standard deviation comparison of the (left) LS and (right) RLS methods. The color bar represents standard deviation values.

For different SNR values, MAE of the *p*O_{2} 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 *p*O_{2} values of retinal vessels with a phosphorescence lifetime imaging model. The method decreased the variance of the classical regression analysis in the generation of *p*O_{2} 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 *p*O_{2} estimates with significantly reduced variance while preserving the mean *p*O_{2} 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 [1]. 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 *p*O_{2} 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).

1. Shahidi M, Shakoor A, Shonat R, Mori M, Blair NP. A method for measurement of chorio-retinal oxygen tension. Curr. Eye Res. 2006;vol. 31:357–366. [PMC free article] [PubMed]

2. Shahidi M, Blair NP, Mori M, Zelkha R. Feasibility of noninvasive imaging of chorio-retinal oxygenation. Ophthalmic Surg., Lasers Imag. 2004;vol. 35(no. 5):415–422. [PubMed]

3. Shonat RD, Kight AC. Oxygen tension imaging in the mouse retina. Ann. Biomed. Eng. 2003;vol. 31:1084–1096. [PubMed]

4. Lakowicz JR, Szmacinski H, Nowaczyk K, Berndt KW, Johnson M. Fluorescence lifetime imaging. Anal. Biochem. 1992;vol. 202:316–330. [PubMed]

5. Rudin L, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. Phys. D. 1992;vol. 60:259–268.

6. Galatsanos N, Katsaggelos AK. Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation. IEEE Trans. Image Process. 1992 Jul.vol. 1(no. 3):332–336. [PubMed]

7. Kang M, Katsaggelos AK. General choice of the regularization functional in regularized image restoration. IEEE Trans. Image Process. 1995 May;vol. 4(no. 5):594–602. [PubMed]

8. Belge M, Kilmer ME, Miller EL. Wavelet domain image restoration with adaptive edge-preserving regularization. IEEE Trans. Image Process. 2000 Apr.vol. 9(no. 4):597–608. [PubMed]

9. Cetin M, Karl WC. Superresolution and edge-preserving reconstruction of complex-valued synthetic aperture radar images. Proc. 2000 IEEE Int. Conf. Image Process.; Sep. 10–13; Vancouver, BC, Canada. pp. 701–704.

10. Yu DF, Fessler JA. Edge-preserving tomographic reconstruction with nonlocal regularization. IEEE Trans. Med. Imag. 2002 Feb.vol. 21(no. 2):159–173. [PubMed]

11. Mignotte M. A segmentation-based regularization term for image deconvolution. IEEE Trans. Image Process. 2006 Jul.vol. 15(no. 7):1973–1984. [PubMed]

12. Le T, Chartrand R, Asaki TJ. A variational approach to reconstructing images corrupted by Poisson noise. J.Math. Imag. Vis. 2007;vol. 27(no. 3):257–263.

13. Yildirim I, Ansari R, Wanek J, Yetik IS, Shahidi M. Retinal oxygen tension estimation in phosphorescence lifetime images using regularized least squares. Proc. IEEE Int. Conf. Electr./Inf. Technol.; May 18–20, 2008; Ames, IA. pp. 465–469.

14. Shahidi M, Wanek J, Blair NP, Mori M. Three-dimensional mapping of chorioretinal vascular oxygen tension in the rat. Invest. Ophthalmol. Vis. Sci. 2008;vol. 50:820–825. [PMC free article] [PubMed]

15. Kamen EW, Su JK. Introduction to Optimal Estimation. London, U.K: Springer-Verlag; 1999.

16. Poor HV. An Introduction to Signal Detection and Estimation. New York: Springer-Verlag; 1994.

17. Cucker F, Smale S. Best choices for regularization parameters in learning theory: On the bias-variance problem. Found. Comput. Math. 2002;vol. 2(no. 4):413–428.

18. Yu D, Cringle SJ, Alder V, Su EN. Intraretinal oxygen distribution in the rat with graded systemic hyperoxia and hypercapnia. Invest. Ophthalmol. Vis. Sci. 1999;vol. 40:2082–2087. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |