 Research
 Open Access
 Published:
A fast regularized leastsquares method for retinal vascular oxygen tension estimation using a phosphorescence lifetime imaging model
BioMedical Engineering OnLinevolume 12, Article number: 106 (2013)
Abstract
Background
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 leastsquares (RLS) method was shown to be very effective compared with the conventional leastsquares (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.
Methods
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 closedform 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.
Results
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.
Conclusion
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.
Background
Retinal tissue requires regular oxygenation to prevent some devastating eye diseases, such as diabetic retinopathy, glaucoma, and agerelated 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 [3]. 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 [4]. 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 leastsquares (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 [8]. 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 [9], biomedical engineering and imaging [10–13], and astronomical imaging [14] motivated the development of an RLS estimation method for oxygen tension in retinal blood vessels using PLIM [8]. In [8], the regularization window was chosen as a 3×3 window, equivalent to 15x15 μm^{2}. 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 [15]. 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 gradientbased 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 nonsmooth convex regularized cost function [16], and Lampe and Voss proposed a fast algorithm for Tikhonovbased regularized total LS estimation problems [17]. Mastronardi, Lemmerling and Van Huffel provided a fast regularized total LS algorithm for solving the basic deconvolution problem [18].
In this study, we propose a computationally efficient RLS method for estimating retinal oxygen tension using a closedform 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 pO_{2} and the phosphorescence lifetime τ, is as follows:
where pO_{2} (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 pO_{2}. 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 m_{n}, θ_{ 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 a_{0} = k[Pd], a_{1} = (1/2)k[Pd]mm_{ n }cos(θ), and b_{1} = (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 pO_{2} at each observation location. Here it should be stressed that pO_{2} 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 [19].
In the presence of additive noise and for multiple observations, Eq. (4) can be written in matrix–vector form as:
where i denotes the ith 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 pseudoinverse of the system matrix in Eq. (7).
RLS estimation of retinal oxygen tension
In [14], the RLS cost function for the parameter a_{1} was defined as:
where Q(2,:)y^{i} is the LS estimate of the parameter ${a}_{1}^{i},$ Q(2,:) is the second row of the pseudoinverse of the system matrix A in Eq. (7), y^{i} is the phosphorescence intensity observation vector of the ith pixel, and β is the regularization parameter. ${\overline{a}}_{1}^{i}$ denotes the mean value of the parameter to be estimated for the ith pixel. The global cost function was defined as:
where M denotes the total number of pixels in the image. A gradientbased 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 ith pixel is defined as:
where y^{i} is the noisy observation vector, A is the system matrix, and γ is a regularization coefficient. In Eq. (11), x^{i} is the parameter vector and ${\overline{x}}^{i}$ is its sample mean, which are defined as:
where a_{0}, a_{1}, and b_{1} 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 crossadjacent pixels, respectively. These are chosen such that l+4p+4q=1. K (j,k) denotes the weight coefficient of the kth pixel on the jth pixel in the oxygen tension map, and R denotes the number of rows in the map. Note that K is a sparse positivedefinite 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 pixelwise 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(X^{T}X) is called the Frobenius inner product [20]. The global cost function using the Frobenius inner product can now be written as:
X and Y in Eq. (16) are defined as:
where ${I}_{s}^{i}$ and S are the sth phosphorescence intensity observation of the ith 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 leastsquares term in (19) can be rewritten as:
Because A^{T}A is equal to
tr(X^{T}A^{T}AX) becomes
From the definition of the Frobenius inner product, tr(Y^{T}AX) can be rewritten as:
Considering the above equations, the global cost function can now be written as:
We need the model parameters a_{1} and b_{1} to estimate the oxygen tension, because this is nonlinearly dependent on the ratio of b_{1} and a_{1}. 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.
In Eqs. (27) and (29), $\left({\scriptscriptstyle \frac{2}{S}}{Y}^{T}A\left(:,2\right)\right)\phantom{\rule{0.25em}{0ex}}\mathrm{and}\phantom{\rule{0.25em}{0ex}}\left({\scriptscriptstyle \frac{2}{S}}{Y}^{T}A\left(:,3\right)\right)$ are the LS estimates of a_{1} and b_{1}, respectively. Therefore, the RLS estimates of a_{1} and b_{1} can be written as:
where ${\widehat{a}}_{1}$ and ${\widehat{b}}_{1}$ are the LS estimates of these parameters. By defining a new matrix L as:
where I is the identity matrix, the RLS estimates of a_{1} and b_{1} can be rewritten in a simpler form as:
The matrix K in Eq. (32) is a symmetric Toeplitz and positivedefinite sparse matrix. Therefore, L is a symmetric positivedefinite matrix. However, it is not in Toeplitz form because the matrix K^{T}K is not in Toeplitz form. The elements of K^{T}K can be written as:
Because K(j,k) = 0 if j − k > R + 1, all nonzero elements are within the range 2R+2, where R is the number of rows in the oxygen tension map. Thus, (K^{T}K)_{ i j } = 0 if i −j > 2R + 2.
If we ignore the first and last 2R+2 rows of the matrix K^{T}K, it is a positivedefinite symmetric Toeplitz matrix. Additionally, as all matrices in L are in Toeplitz form except for K^{T}K, the matrix L has the same properties as K^{T}K. That is, except for the first and last 2R+2 rows, it is a purely positivedefinite 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 pO_{2} images are being used. For a 640 × 480 pO_{2} 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 ${L}_{s}^{1}$ and acquire its most significant units, then extend the solution to larger pO_{2} maps with negligible errors.
The proposed method is visualized in Figure 1 as a flow chart: (A) Using the phosphorescence lifetime images, (BC) the LS estimates of the model parameters a_{1} and b_{1} of the original pO_{2} map are found in matrix form. (D)${L}_{s}^{1}$, whose dimension is R_{ s }^{2}× R_{ s }^{2}, is found for a small R_{ s }× R_{ s }pO_{2} 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 ${L}_{s}^{1}$ 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 a_{1} and b_{1} 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 ${L}_{s}^{1}$ as a 169 × 169 matrix for a 13 × 13 pO_{2} 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 preset values of β, window size, and signaltonoise 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 [6]. 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 intraperitoneal 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. Pdporphine (Frontier Scientific, Logan, Utah), an oxygensensitive 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 twodisk diameters (600 microns) from the edge of the optic nerve head. For each pixel location, 10 phasedelayed 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 enface 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 [21], venous and arterial oxygen tensions of rat retina vary around 35 mmHg and 60 mmHg, 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 [15]. Since the iterative method requires a considerable amount of time as size of the pO_{2} map gets larger, the simulated pO_{2} 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 pO_{2} 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 pO_{2} maps used in the estimation of ${L}_{s}^{1}$. 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 [14], the RLS estimation generates smoother pO_{2} 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 [6]. 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.
Abbreviations
 LS:

Least squares
 RLS:

Regularized least squares
 PLIM:

Phosphorescence lifetime imaging
 FLIM:

Fluorescence lifetime imaging
 pO2:

Oxygen tension.
References
 1.
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.14401681.1997.tb02133.x
 2.
Berkowitz BA, Kowluru RA, Frank RN, Kern TS, Hohman TC, Prakash M: Subnormal retinal oxygenation response precedes diabeticlike retinopathy. Invest Ophthalmol Vis Sci 1999, 40: 2100–2105.
 3.
WangsaWirawan ND, Linsenmeier RA: Retinal oxygen: fundamental and clinical aspects. Arch Ophthalmology 2003, 121: 547–557. 10.1001/archopht.121.4.547
 4.
Ito Y, Berkowitz BA: MR studies of retinal oxygenation. Vision Res 2001, 41: 1307–1311. 10.1016/S00426989(00)002583
 5.
Shonat RD, Kight AC: Oxygen tension imaging in the mouse retina. Ann Biomed Eng 2003, 31: 1084–1096.
 6.
Shahidi M, Shakoor A, Shonat R, Mori M, Blair NP: A method for measurement of chorioretinal oxygen tension. Curr Eye Res 2006, 31: 357–366. 10.1080/02713680600599446
 7.
Shahidi M, Wanek J, Blair NP, Mori M: Threedimensional mapping of chorioretinal vascular oxygen tension in the rat. Invest Ophthalmol Vis Sci 2009, 50: 820–825.
 8.
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.
 9.
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.
 10.
Funai AK, Fessler JA, Yeo DTB, Olafsson VT, Noll DC: Regularized field map. Estimation in MRI. IEEE Trans Med Imaging 2008, 27: 1484–1494.
 11.
Zhu W, Wang Y, Deng Y, Yao Y, Barbour RL: A waveletbased multiresolution regularized least squares reconstruction approach for optical tomography. IEEE Trans Med Imaging 1997, 16: 210–217. 10.1109/42.563666
 12.
Tarvainen MP, Rantaaho PO, Karjalainen PA: An advanced detrending method with application to HRV analysis. IEEE Trans Biomed Eng 2002, 49: 172–175. 10.1109/10.979357
 13.
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
 14.
Frazin RA: Tomography of the solar corona, a robust, regularized, positive estimation method. Astrophysical J 2008, 530: 1026–1035.
 15.
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.060291
 16.
Toh KC, Yun S: An accelerated proximal gradient algorithm for nuclear norm regularized linear least squares problems. Pac J Optimization 2008, 6: 615–640.
 17.
Lampe J, Voss H: A fast algorithm for solving regularized total least squares problems. Electron Trans Numerical Anal 2008, 31: 12–24.
 18.
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.
 19.
Shonat RD, Wachman ES, Niu WH, Koretsky AP, Farkas DL: Nearsimultaneous hemoglobin saturation and oxygen tension maps in mouse brain using an AOTF microscope. Biophys J 1997, 73: 1223–1231. 10.1016/S00063495(97)781554
 20.
Horn RA, Johnson CR: Matrix Analysis. Cambridge, UK: Cambridge University Press; 1985.
 21.
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.020105
Acknowledgement
We would like to thank Dr. Mahnaz Shahidi for providing us the experimental data which was obtained from rat retinas in her lab.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
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.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Received
Accepted
Published
DOI
Keywords
 Accelerated estimation
 Closedform solution
 Regularized estimation
 Retinal oxygen tension
 Phosphorescence lifetime imaging