- Open Access
A fast regularized least-squares method for retinal vascular oxygen tension estimation using a phosphorescence lifetime imaging model
BioMedical Engineering OnLinevolume 12, Article number: 106 (2013)
Monitoring retinal oxygenation is of primary importance in detecting the presence of some common eye diseases. To improve the estimation of oxygen tension in retinal vessels, regularized least-squares (RLS) method was shown to be very effective compared with the conventional least-squares (LS) estimation. In this study, we propose an accelerated RLS estimation method for the problem of assessing the oxygenation of retinal vessels from phosphorescence lifetime images.
In the previous work, gradient descent algorithms were used to find the minimum of the RLS cost function. This approach is computationally expensive, especially when the oxygen tension map is large. In this study, using a closed-form solution of the RLS estimation and some inherent properties of the problem at hand, the RLS process is reduced to the weighted averaging of the LS estimates. This decreases the computational complexity of the RLS estimation considerably without sacrificing its performance.
Performance analyses are conducted using both real and simulated data sets. In terms of computational complexity, the proposed RLS estimation method is significantly better than RLS methods that use gradient descent algorithms to find the minimum of the cost function. Additionally, there is no significant difference between the estimates acquired by the proposed and conventional RLS estimation methods.
The proposed RLS estimation method for computing the retinal oxygen tension is computationally efficient, and produces estimates with negligible difference from those obtained by iterative RLS methods. Further, the results of this study can be applied to other lifetime imaging problems that have similar properties.
Retinal tissue requires regular oxygenation to prevent some devastating eye diseases, such as diabetic retinopathy, glaucoma, and age-related macular degeneration [1, 2]. Oxygen tension in retinal vessels should therefore be monitored for use in the early diagnosis of these devastating eye diseases.
Oxygen sensitive microelectrodes can be used to obtain retinal oxygenation . However, this is an invasive method, as the electrodes damage the microenvironment of the retina. Magnetic resonance imaging, a noninvasive imaging technique, has been used to measure retinal oxygen tension . This approach is limited by its low resolution compared with optical imaging methods.
The oxygen tension of retinal vessels can also be measured using a phosphorescence lifetime imaging model (PLIM) [5–7]. PLIM is a noninvasive method, though it requires intravenous injection, and does not damage the microenvironment of the retina. The least-squares (LS) estimation method has been used to estimate the oxygenation of retinal vessels based on PLIM. However, the LS method is not robust against noise contamination, producing high variance and artificial peaks in the estimates, as well as values outside of the physiological range, which can be attributed to the no use of a suitable prior model for the data.
To overcome these shortcomings of the LS method, a regularized LS (RLS) estimation method was proposed . The RLS method was developed by utilizing knowledge of the prior distribution of the model parameters. By applying the RLS method to simulated data as well as real image data acquired from rat retinas, the method was shown to be robust to the presence of noise, and generated estimates that were more in the physiologically expected range and whose variance was lower than that obtained with the LS method which does not use a suitable prior model for the data.
Successful applications of regularization in image processing , biomedical engineering and imaging [10–13], and astronomical imaging  motivated the development of an RLS estimation method for oxygen tension in retinal blood vessels using PLIM . In , the regularization window was chosen as a 3×3 window, equivalent to 15x15 μm2. The choice of this size of the window was inspired by the fact that variation among oxygen tension values in a neighborhood of this size of window should physiologically be negligible . The regularization term in the RLS cost function was developed using this physiological information and assuming that the mean of a pixel in an oxygen tension map of retinal blood vessels is equal to the sample mean of oxygen tension values in neighboring pixels. The RLS estimation method was shown to be superior to the conventional LS estimation method as mentioned before. However, the gradient-based iterative methods used to find the minimum of the RLS cost function are computationally costly, especially when considering large oxygen tension maps.
Existing applications have a typical oxygen tension map frame size of around 512 × 160 pixels. In the early diagnosis of some eye diseases, analyses of small capillaries are of great importance and require imaging systems with higher resolutions and larger frame sizes. However, as the map frame size becomes larger, the computational cost of implementing the RLS estimation method increases drastically. Therefore, to realize RLS estimation for large oxygen tension maps, faster estimation methods are needed.
In the literature, fast regularized estimation methods have been proposed. For instance, Toh and Yun proposed an accelerated proximal gradient method to minimize a non-smooth convex regularized cost function , and Lampe and Voss proposed a fast algorithm for Tikhonov-based regularized total LS estimation problems . Mastronardi, Lemmerling and Van Huffel provided a fast regularized total LS algorithm for solving the basic deconvolution problem .
In this study, we propose a computationally efficient RLS method for estimating retinal oxygen tension using a closed-form solution. The proposed method is shown to be much faster than the iterative solution, and generates almost identical results. Additionally, the approach derived in this study can be applied to other imaging problems using PLIM or FLIM with some minor changes.
Past RLS estimation method
LS estimation of retinal oxygen tension
The Stern–Volmer equation, giving the relationship between the oxygen tension pO2 and the phosphorescence lifetime τ, is as follows:
where pO2 (in millimeters of mercury) denotes oxygen tension, K ϕ is the quenching constant, and τ0 denotes the lifetime in a zero oxygen environment. Based on Eq. (1), we need to estimate τ to find the value of pO2. The phosphorescence lifetime τ can be calculated from the distribution of phosphorescence intensity values in different modulation phase images, because the intensity depends on the phase angle, denoted by θ, between the modulated excitation laser light and the emitted phosphorescence [5, 6]. The relation between θ and τ is
where ω is the modulation frequency. The intensity can be represented as a function of θ n ,
where mn, θ n , and n p denote the modulation profile of the image intensifier, the phase of the gain modulation, and the number of phosphorescence intensity observations at each pixel location for different phase values of θ n , respectively. The concentration of the probe [Pd], k, and the modulation m are unavailable.
If we define a0 = k[Pd], a1 = (1/2)k[Pd]mm n cos(θ), and b1 = (1/2)k[Pd]mm n sin(θ), we can express Eq. (3) as:
The phase angle θ in Eq. (3) is obtained as:
Therefore, using Equations (1), (2), and (5), we are able to acquire values for pO2 at each observation location. Here it should be stressed that pO2 in the Equation 1 can also be found using phosphorescence quenching measurement. However, for multiple and closely spaced lifetime applications, the PLIM method satisfies needs of the application more .
In the presence of additive noise and for multiple observations, Eq. (4) can be written in matrix–vector form as:
where i denotes the i-th pixel and n denotes the number of phosphorescence intensity observations for each pixel. We can rewrite Eq. (6) as follows:
where y, A, x, and n are the phosphorescence lifetime intensity observations, system matrix, PLIM parameters, and additive noise, respectively. Clearly, to find the model parameters a 0 , a 1 , and b 1 , at least three phosphorescence intensity observations must be obtained at three different gain modulation phases for each pixel, though in practice we need more than three observations because of noise contamination. The LS estimate of the model parameters is obtained as follows:
where Q is the pseudo-inverse of the system matrix in Eq. (7).
RLS estimation of retinal oxygen tension
In , the RLS cost function for the parameter a1 was defined as:
where Q(2,:)yi is the LS estimate of the parameter Q(2,:) is the second row of the pseudo-inverse of the system matrix A in Eq. (7), yi is the phosphorescence intensity observation vector of the i-th pixel, and β is the regularization parameter. denotes the mean value of the parameter to be estimated for the i-th pixel. The global cost function was defined as:
where M denotes the total number of pixels in the image. A gradient-based iterative approach was then used to find the minimum of the RLS cost function.
Derivation of the proposed solution for the RLS estimation
The RLS cost function for the model parameters a 0 , a 1 , and b 1 for the i-th pixel is defined as:
where yi is the noisy observation vector, A is the system matrix, and γ is a regularization coefficient. In Eq. (11), xi is the parameter vector and is its sample mean, which are defined as:
where a0, a1, and b1 are vectors of the phosphorescence lifetime imaging model parameters and K is a weighting matrix defining the relations between these parameters. K is defined as follows:
where l, p, and q denote the weights of a pixel with respect to itself, to directly adjacent pixels, and to cross-adjacent pixels, respectively. These are chosen such that l+4p+4q=1. K (j,k) denotes the weight coefficient of the k-th pixel on the j-th pixel in the oxygen tension map, and R denotes the number of rows in the map. Note that K is a sparse positive-definite symmetric Toeplitz matrix.
From Equations (11) and (12), we see that the cost function for a pixel is dependent on the values of its neighboring pixels. Therefore, there is no pixel-wise solution, and the problem must be dealt with globally. The global cost function is written as:
where M is the number of pixels in the oxygen tension map.
For any ℜIxJ, ǁX ǁF is called the Frobenius norm in ℜIxJ space, and
where tr denotes the trace of a matrix. The operation tr(XTX) is called the Frobenius inner product . The global cost function using the Frobenius inner product can now be written as:
X and Y in Eq. (16) are defined as:
where and S are the s-th phosphorescence intensity observation of the i-th pixel and the number of observations per pixel, respectively. Given the definition of X, Eqs. (12) and (16) can be rewritten as follows:
For the model parameters, the regularization term in the cost function (19) can be expanded as:
The least-squares term in (19) can be rewritten as:
Because ATA is equal to
From the definition of the Frobenius inner product, tr(YTAX) can be rewritten as:
Considering the above equations, the global cost function can now be written as:
We need the model parameters a1 and b1 to estimate the oxygen tension, because this is nonlinearly dependent on the ratio of b1 and a1. Taking the gradient of the cost function with respect to the model parameters and equating to zero, we obtain the RLS estimates of the model parameters.
where and are the LS estimates of these parameters. By defining a new matrix L as:
where I is the identity matrix, the RLS estimates of a1 and b1 can be rewritten in a simpler form as:
The matrix K in Eq. (32) is a symmetric Toeplitz and positive-definite sparse matrix. Therefore, L is a symmetric positive-definite matrix. However, it is not in Toeplitz form because the matrix KTK is not in Toeplitz form. The elements of KTK can be written as:
Because K(j,k) = 0 if |j − k| > R + 1, all non-zero elements are within the range 2R+2, where R is the number of rows in the oxygen tension map. Thus, (KTK) i j = 0 if |i −j| > 2R + 2.
If we ignore the first and last 2R+2 rows of the matrix KTK, it is a positive-definite symmetric Toeplitz matrix. Additionally, as all matrices in L are in Toeplitz form except for KTK, the matrix L has the same properties as KTK. That is, except for the first and last 2R+2 rows, it is a purely positive-definite symmetric Toeplitz matrix. The inverse of L shares the same properties as L. Therefore, if we choose a row between the first and last 2R+2 rows of the inverse of L, we can reconstruct it with negligible errors in the first and last 2R+2 rows.
Further, even if we change the size of L, there occurs only a negligible difference in the most significant elements of its inverse. This makes the matrix L and its inverse very useful when large pO2 images are being used. For a 640 × 480 pO2 map, L is a 307200 × 307200 matrix whose inverse cannot be calculated using a regular computer. In such cases, iterative approaches must be employed to find the optimum points of the RLS cost function, and these are also computationally expensive. Because the most significant elements in L−1 do not vary considerably when its dimension is reduced, we can calculate a small and acquire its most significant units, then extend the solution to larger pO2 maps with negligible errors.
The proposed method is visualized in Figure 1 as a flow chart: (A) Using the phosphorescence lifetime images, (B-C) the LS estimates of the model parameters a1 and b1 of the original pO2 map are found in matrix form. (D), whose dimension is R s 2× R s 2, is found for a small R s × R s pO2 map, and (E) its middle row is chosen. R s is chosen as an odd number to ensure symmetry with respect to the middle unit in the chosen row. Because is in Toeplitz form except for the first and last two rows, all other rows can be recovered exactly using this row with appropriate shifts. Additionally, elements in the middle row that determine the weights to obtain approximate RLS estimates of the model parameters (F) can be reshaped to give an R s × R s weighted averaging window. (G) Finally, this weighted averaging window is applied to the LS estimates of the a1 and b1 parameter maps to find approximate RLS estimates of the model parameters.
Simulated and experimental data results
In the RLS method, the standard Newton–Raphson method is employed to find the minimum of the RLS cost function. The LS estimates are used as the initial values of the model parameters to be estimated. In the proposed RLS method, we first determine as a 169 × 169 matrix for a 13 × 13 pO2 map. We then follow the procedure given in the previous section. The computer system used for these simulations is a standard notebook with an Intel Core Duo 2.27 GHz processor and 4 GB memory. Unless otherwise stated, the pre-set values of β, window size, and signal-to-noise ratio (SNR) are 5, 13 × 13, and 20 dB, respectively.
We use both simulated and real oxygen tension maps (Figures 2 and 3, respectively) in our experiments. The real retinal oxygen tension map used in this study was acquired from a Long Evans pigmented rat (~500 g) using a novel system described in . The animal was treated according to the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research. In order to anesthetize the rat, an intra-peritoneal infusion of Ketamine (85 mg/kg IP) and Xylazine (3.5 mg/kg IP) at the respective rates of 0.5 mg/kg/min and 0.02 mg/kg/min was used. Before and during the imaging, gas mixture containing 21% oxygen (room air, normoxia) was administered to rats for 5 minutes. Pd-porphine (Frontier Scientific, Logan, Utah), an oxygen-sensitive molecular probe, was dissolved (12 mg/ml) in bovine serum albumin solution (60 mg/ml) and physiological saline buffered to a pH of 7, and injected intravenously (20 mg/kg). Imaging was conducted at areas within two-disk diameters (600 microns) from the edge of the optic nerve head. For each pixel location, 10 phase-delayed optical section phosphorescence lifetime intensity images were acquired with a corresponding phase delay increments of 74 μs, and then these intensity images were brought together to form en-face phosphorescence images of the retinal vasculatures in the retinal plane. To generate the simulated data we used in this work, physiological features and topology of real retinal vessels were followed based on the previous studies [6–8, 15]. Near the optic disc, which has a diameter nearly 300 microns , venous and arterial oxygen tensions of rat retina vary around 35 mm-Hg and 60 mm-Hg, respectively in the normoxia condition. As getting far from the optic disc, arterial oxygen tension decreases almost linearly while venous oxygen tension remains almost the same . Since the iterative method requires a considerable amount of time as size of the pO2 map gets larger, the simulated pO2 map was generated to be 485x600 pixels. But it should be noted that this size is not a limitation to the proposed method. In the simulations, we add i.i.d. white Gaussian noise with 15, 20, and 25 dB SNR to the phosphorescence lifetime images.
We examined the MAE performance of the LS, iterative RLS, and proposed RLS methods for different regularization coefficient values. As shown in Figure 4, there is a negligible difference between the MAE of the iterative and proposed RLS estimation methods. On the other hand, there is a significant difference in computation time for the iterative and proposed RLS methods (see Figure 5). As mentioned previously, the Newton–Raphson method is used to minimize the RLS cost functions, and its step size, for our problem, is as follows:
Thus, the step size decreases for higher values of β, and this results in a proportional increase in computation time.
The time difference between the two methods would be much more notable if the pO2 map to be estimated were of higher dimensions. Therefore, it is clear that the proposed method is preferable to the iterative one in the sense of MAE performance versus computational complexity.
From Figure 6, we can see that there is a small difference between the proposed and iterative approaches in terms of MAE performance for different sizes of the small pO2 maps used in the estimation of . As their size increases, the MAE of the proposed method converges to that of the standard iterative RLS method.
In Figure 7, we present the computation times of the iterative and proposed RLS estimation methods for a real oxygen tension map which was extended from original form 120x85 pixels to 600x425 pixels. As seen from Figures 4 and 6, the results for the real and simulated oxygen tension maps are proportional in terms of computation time.
Figure 8 compares the visual results of the proposed and LS estimation methods for the frame shown in the rectangle in Figure 1. This illustrates the artifacts of the LS estimation more clearly. As mentioned in , the RLS estimation generates smoother pO2 maps whose values fall within the physiologically expected range. The acquisition of the real oxygen tension map of a rat’s retina, shown in Figure 3, is described in . In Figure 3, we compare the visual results of the proposed RLS, iterative RLS, and LS estimation methods for the real data. When compared with the LS result, the iterative and proposed RLS methods generate smoother oxygen tension maps.
Discussions and conclusion
Presently, by use of PLIM method, available oxygen tension maps have frame sizes around 512x160 pixels. As the imaging instruments become enhanced, this size will rise considerably. The RLS method implemented in an iterative way can handle oxygen tension images having this size of frames but as shown in the Simulated and Experimental Data Results Section, it becomes slower as size of the images gets larger. The proposed method addresses this problem and can be used even for much larger oxygen tension maps. On the other hand, the proposed method is limited by the assumption made when developing the regularization window that variation among oxygen tension values of retinal vessels within this size of window should be minimal. For the problems in which this assumption is not valid, the proposed method can generate unsatisfactory results. In this study, by exploiting intrinsic properties of the problem, we have developed a fast RLS estimation method that is applicable to large oxygen tension maps. By comparing the results of the available and proposed approaches, it was shown that, although there is a minor difference in estimation performance, the proposed method is significantly faster. Additionally, the proposed method is not restricted to the problem of oxygen tension estimation, and is applicable to other problems where a similarly strong relationship exists between neighboring pixels.
Regularized least squares
Phosphorescence lifetime imaging
Fluorescence lifetime imaging
Alder VA, Su EN, Yu DY, Cringle SJ, Yu PK: Diabetic retinopathy: early functional changes. Clin Exp Pharmacol Physiol 1997, 24: 785–788. 10.1111/j.1440-1681.1997.tb02133.x
Berkowitz BA, Kowluru RA, Frank RN, Kern TS, Hohman TC, Prakash M: Subnormal retinal oxygenation response precedes diabetic-like retinopathy. Invest Ophthalmol Vis Sci 1999, 40: 2100–2105.
Wangsa-Wirawan ND, Linsenmeier RA: Retinal oxygen: fundamental and clinical aspects. Arch Ophthalmology 2003, 121: 547–557. 10.1001/archopht.121.4.547
Ito Y, Berkowitz BA: MR studies of retinal oxygenation. Vision Res 2001, 41: 1307–1311. 10.1016/S0042-6989(00)00258-3
Shonat RD, Kight AC: Oxygen tension imaging in the mouse retina. Ann Biomed Eng 2003, 31: 1084–1096.
Shahidi M, Shakoor A, Shonat R, Mori M, Blair NP: A method for measurement of chorio-retinal oxygen tension. Curr Eye Res 2006, 31: 357–366. 10.1080/02713680600599446
Shahidi M, Wanek J, Blair NP, Mori M: Three-dimensional mapping of chorioretinal vascular oxygen tension in the rat. Invest Ophthalmol Vis Sci 2009, 50: 820–825.
Yildirim I, Ansari R, Wanek J, Yetik IS, Shahidi M: Regularized estimation of retinal vascular oxygen tension from phosphorescence images. IEEE Trans Biomed Eng 2009, 56(8):1989–1995.
Elmoataz A, Lezoray O, Bougleux S: Nonlocal discrete regularization on weighted graphs: a framework for image and manifold processing. IEEE Trans Image Process 2008, 17: 1047–1060.
Funai AK, Fessler JA, Yeo DTB, Olafsson VT, Noll DC: Regularized field map. Estimation in MRI. IEEE Trans Med Imaging 2008, 27: 1484–1494.
Zhu W, Wang Y, Deng Y, Yao Y, Barbour RL: A wavelet-based multiresolution regularized least squares reconstruction approach for optical tomography. IEEE Trans Med Imaging 1997, 16: 210–217. 10.1109/42.563666
Tarvainen MP, Ranta-aho PO, Karjalainen PA: An advanced detrending method with application to HRV analysis. IEEE Trans Biomed Eng 2002, 49: 172–175. 10.1109/10.979357
Greensite F, Huiskamp G: An improved method for estimating epicardial potentials from the body surface. IEEE Trans Biomed Eng 1998, 45: 98–104. 10.1109/10.650360
Frazin RA: Tomography of the solar corona, a robust, regularized, positive estimation method. Astrophysical J 2008, 530: 1026–1035.
Shakoor A, Blair NP, Mori M, Shahidi M: Chorioretinal vascular oxygen tension changes in response to light flicker. Invest Ophthalmol Vis Sci 2006, 47(11):4962–4965. 10.1167/iovs.06-0291
Toh KC, Yun S: An accelerated proximal gradient algorithm for nuclear norm regularized linear least squares problems. Pac J Optimization 2008, 6: 615–640.
Lampe J, Voss H: A fast algorithm for solving regularized total least squares problems. Electron Trans Numerical Anal 2008, 31: 12–24.
Mastronardi N, Lemmerling P, Huffel SV: Fast regularized structured total least squares algorithm for solving the basic deconvolution problem. Numerical Linear Algebra Appl 2004, 12: 201–209.
Shonat RD, Wachman ES, Niu WH, Koretsky AP, Farkas DL: Near-simultaneous hemoglobin saturation and oxygen tension maps in mouse brain using an AOTF microscope. Biophys J 1997, 73: 1223–1231. 10.1016/S0006-3495(97)78155-4
Horn RA, Johnson CR: Matrix Analysis. Cambridge, UK: Cambridge University Press; 1985.
Cohan BE, Pearch AC, Jokelainen PT, Bohr DF: Optic disc imaging in conscious rats and mice. Invest Ophthalmol Vis Sci 2003, 44: 160–163. 10.1167/iovs.02-0105
We would like to thank Dr. Mahnaz Shahidi for providing us the experimental data which was obtained from rat retinas in her lab.
The authors declare that they have no competing interests.
Each Author has contributed substantially to the preparation of the paper and approves of its submission to the Journal. Both authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.