A fast regularized least-squares method for retinal vascular oxygen tension estimation using a phosphorescence lifetime imaging model

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

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

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 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 [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 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 [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². 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 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 [16], and Lampe and Voss proposed a fast algorithm for Tikhonov-based 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 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 pO₂ and the phosphorescence lifetime τ, is as follows:

$$\frac{{\tau }_0}{\tau }=1+K_{\phi }{\tau }_0\left(p{O}_2\right),$$

(1)

where pO₂ (in millimeters of mercury) denotes oxygen tension, K _ϕ is the quenching constant, and τ₀ denotes the lifetime in a zero oxygen environment. Based on Eq. (1), we need to estimate τ to find the value of pO₂. 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

$$\mathit{tan}\theta =\mathit{\omega \tau }\mathit{,}$$

(2)

where ω is the modulation frequency. The intensity can be represented as a function of θ _n ,

$$\begin{array}{ll}I\left(\mathit{\theta n}\right)=& k\left[\mathrm{Pd}\right]\\ +\frac{1}{2}k\left[\mathrm{Pd}\right]m{m}_n\left(cos\left(\theta \right)cos\left(\mathit{\theta n}\right)+sin\left(\theta \right)sin\left(\mathit{\theta n}\right)\right)\cdots \cdots n\mathit{=}1,2,\dots ,n_p\end{array}$$

(3)

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₀ = k[Pd], a₁ = (1/2)k[Pd]mm _n cos(θ), and b₁ = (1/2)k[Pd]mm _n sin(θ), we can express Eq. (3) as:

$$I\left(\mathit{\theta n}\right)=a_0+a_{1}cos\left(\mathit{\theta n}\right)+b_{1}sin\left(\mathit{\theta n}\right).$$

(4)

The phase angle θ in Eq. (3) is obtained as:

$$\theta =\mathit{ta}{n}^{-1}\left(b_1/a_1\right).$$

(5)

Therefore, using Equations (1), (2), and (5), we are able to acquire values for pO₂ at each observation location. Here it should be stressed that pO₂ 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:

$$\left[\begin{array}{c}\hfill I_1^i\hfill \\ \hfill I_2^i\hfill \\ \hfill ⋮\hfill \\ \hfill ⋮\hfill \\ \hfill I_n^i\hfill \end{array}\right]=\left[\begin{array}{c}\hfill 1\cos \left({\theta }_1\right)\sin \left({\theta }_1\right)\hfill \\ \hfill 1\cos \left({\theta }_2\right)\sin \left({\theta }_2\right)\hfill \\ \hfill ⋮\hfill \\ \hfill ⋮\hfill \\ \hfill 1\cos \left({\theta }_n\right)\sin \left({\theta }_n\right)\hfill \end{array},\right]\left[\begin{array}{l}{a}_0^i\\ a_1^i\\ b_1^i\end{array}\right]+\left[\begin{array}{c}\hfill N_1^i\hfill \\ \hfill N_2^i\hfill \\ \hfill ⋮\hfill \\ \hfill ⋮\hfill \\ \hfill N_n^i\hfill \end{array}\right]$$

(6)

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:

$$y\mathit{=}\mathit{Ax}\mathit{+}n$$

(7)

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 ₀ , a ₁ , and b ₁ , 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:

$${\left[{\widehat{a}}_0{\widehat{a}}_1{\widehat{b}}_1\right]}^T=Q{\underset{¯}{I}}^{\mathit{i}},$$

(8)

where Q is the pseudo-inverse of the system matrix in Eq. (7).

RLS estimation of retinal oxygen tension

In [14], the RLS cost function for the parameter a₁ was defined as:

$$f_{a_1}^i={\left(a_1^i-Q\left(2,:\right)y^i\right)}^2+\beta {\left(a_1^i-{\overline{a}}_1^i\right)}^2$$

(9)

where Q(2,:)yⁱ is the LS estimate of the parameter $a_1^i,$ Q(2,:) is the second row of the pseudo-inverse of the system matrix A in Eq. (7), yⁱ is the phosphorescence intensity observation vector of the i-th pixel, and β is the regularization parameter. ${\overline{a}}_1^i$ denotes the mean value of the parameter to be estimated for the i-th pixel. The global cost function was defined as:

$$F_{a_1}={\displaystyle \sum _1^M{f}_{a_1}^i,}$$

(10)

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 ₀ , a ₁ , and b ₁ for the i-th pixel is defined as:

$$C_{\mathit{RLS}}^i=∥y^i+{A{x}^i∥}_2^2+\gamma ∥x^i-{{\overline{x}}^i∥}_2^2,$$

(11)

where yⁱ is the noisy observation vector, A is the system matrix, and γ is a regularization coefficient. In Eq. (11), xⁱ is the parameter vector and ${\overline{x}}^i$ is its sample mean, which are defined as:

$$x^i={\left[a_0^i{a}_1^i{b}_1^i\right]}^T\mathrm{and}{\overline{x}}^i={\left[K\left(i\mathit{,}:\right)a_{0}K\left(i\mathit{,}:\right)a_{1}K\left(i\mathit{,}:\right)b_1\right]}^T,$$

(12)

where a₀, a₁, and b₁ 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:

$$K\left(j,k\right)=\left\{\begin{array}{l}l\mathit{if}\left|j-\left.k\right|=0\right.\\ p\mathit{if}\left|j-\left.k\right|=1\mathit{or}R\right.\\ q\mathit{if}\left|j-\left.k\right|=R\pm 1\right.\\ 0\mathit{otherwise}\end{array}\right.$$

(13)

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:

$$C_{\mathit{RLS}}={\displaystyle \sum _{i\mathit{=}1}^M{C}_{\mathit{RLS}}^i,}$$

(14)

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

$$\parallel X{\parallel }_F=\sqrt{\mathit{tr}\left(X^{T}X\right)}=\sqrt{\mathit{tr}\left(X{X}^T\right)}=\sqrt{{\displaystyle \sum _{i\mathit{=}1}^I{\displaystyle \sum _{j=1}^J{X}_{\mathit{ij}}{X}_{\mathit{ij}}},}}$$

(15)

where tr denotes the trace of a matrix. The operation tr(X^TX) is called the Frobenius inner product [20]. The global cost function using the Frobenius inner product can now be written as:

$$C_{\mathit{RLS}}=\parallel Y-\mathit{AX}{\parallel }_F^2+\gamma \parallel X-\overline{X}{\parallel }_F^2.$$

(16)

X and Y in Eq. (16) are defined as:

$$x=\left|\begin{array}{c}\hfill a_0^1\cdots a_0^i\cdots a_0^M\hfill \\ \hfill a_1^1\cdots a_1^i\cdots a_1^M\hfill \\ \hfill b_1^1\cdots b_1^i\cdots b_1^M\hfill \end{array}\right|=\left[x^1\dots x^i\dots x^M\right]\mathrm{and}Y=\left|\begin{array}{c}\hfill I_1^1\cdots I_1^i\cdots I_1^M\hfill \\ \hfill ⋮\cdots ⋮\cdots ⋮\hfill \\ \hfill I_S^1\cdots I_S^i\cdots I_S^M\hfill \end{array}\right|,$$

(17)

where $I_s^i$ 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:

$${\overline{x}}^i={\left(K\left(i\mathit{,}:\right)X^T\right)}^T,$$

(18)

$$C_{\mathit{RLS}}=\parallel Y-\mathit{AX}{\parallel }_F^2+\gamma \parallel X-{\left(K{X}^T\right)}^T{\parallel }_F^2.$$

(19)

For the model parameters, the regularization term in the cost function (19) can be expanded as:

$$\parallel X-{\left(K{X}^T\right)}^T{\parallel }_F^2=\parallel a_0-K{a}_0{\parallel }_2^2+\parallel a_1-K{a}_1{\parallel }_2^2+\parallel b_1-K{b}_1{\parallel }_2^2.$$

(20)

The least-squares term in (19) can be rewritten as:

$$\parallel Y-\mathit{AX}{\parallel }_F^2=\mathit{tr}\left(Y^{T}Y\right)-2\mathit{tr}\left(Y^T\mathit{AX}\right)+\mathit{tr}\left(X^T{A}^T\mathit{AX}\right).$$

(21)

Because A^TA is equal to

$$A^{T}A=\left|\begin{array}{c}\hfill S00\hfill \\ \hfill 0S\mathit{/}20\hfill \\ \hfill 00S\mathit{/}2\hfill \end{array}\right|,$$

(22)

tr(X^TA^TAX) becomes

$$\mathit{tr}\left(X^T{A}^T\mathit{AX}\right)=S{a}_0^T{a}_0+\frac{S}{2}{a}_1^T{a}_1+\frac{S}{2}{b}_1^{\mathit{T}}{b}_1.$$

(23)

From the definition of the Frobenius inner product, tr(Y^TAX) can be rewritten as:

$$\mathit{tr}\left(Y^T\mathit{AX}\right)=\mathit{tr}\left(X{Y}^{T}A\right)=a_0^T{Y}^{T}A\left(:,1\right)+a_1^T{Y}^{T}A\left(:,2\right)+b_1^T{Y}^{T}A\left(:,3\right).$$

(24)

Considering the above equations, the global cost function can now be written as:

$$\begin{array}{l}{C}_{\mathit{RLS}}=\mathit{tr}\left(Y^{T}Y\right)+a_0^T\left(S{a}_0-2Y^{T}A\left(:,1\right)\right)+\hfill \\ a_1^T\left(\frac{S}{2}{a}_1-2Y^{T}A\left(:,2\right)\right)+b_1^T\left(\frac{S}{2}{b}_1-2Y^{T}A\left(:,3\right)\right)+\hfill \\ \gamma \left(\parallel a_0-K{a}_0{\parallel }_2^2+\parallel a_1-K{a}_1{\parallel }_2^2+\parallel b_1-K{b}_1{\parallel }_2^2\right).\hfill \end{array}$$

(25)

We need the model parameters a₁ and b₁ to estimate the oxygen tension, because this is nonlinearly dependent on the ratio of b₁ and a₁. 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.

$${\nabla }_{a_1}{C}_{\mathit{RLS}}=\left(S{a}_1-2Y^{T}A\left(:,2\right)\right)+\mathit{S\beta }\left(I-2K+K^2\right)a_1=0,$$

(26)

$${\widehat{a}}_{1-\mathit{RLS}}={\left(I+\beta \left(I-2K+K^2\right)\right)}^{-1}\left(\frac{2}{S}{Y}^{T}A\left(:,2\right)\right),$$

(27)

$${\nabla }_{b_1}{C}_{\mathit{RLS}}=\left(S{b}_1-2Y^{T}A\left(:,3\right)\right)+\mathit{S\beta }\left(I-2K+K^2\right)b_1=0,$$

(28)

$${\widehat{b}}_{1-\mathit{RLS}}={\left(I+\beta \left(I-2K+K^2\right)\right)}^{-1}\left(\frac{2}{S}{Y}^{T}A\left(:,3\right)\right).$$

(29)

In Eqs. (27) and (29), $\left(\frac{2}{S}{Y}^{T}A\left(:,2\right)\right)\mathrm{and}\left(\frac{2}{S}{Y}^{T}A\left(:,3\right)\right)$ are the LS estimates of a₁ and b₁, respectively. Therefore, the RLS estimates of a₁ and b₁ can be written as:

$${\widehat{a}}_{1-\mathit{RLS}}=(\mathrm{I}+\beta {\left(I+K^{T}K-K-K^T\right)}^{-1}{\widehat{a}}_1,$$

(30)

$${\widehat{b}}_{1-\mathit{RLS}}=(\mathrm{I}+\beta {\left(I+K^{T}K-K-K^T\right)}^{-1}{\widehat{b}}_1,$$

(31)

where ${\widehat{a}}_1$ and ${\widehat{b}}_1$ are the LS estimates of these parameters. By defining a new matrix L as:

$$L_{\mathit{MxM}}=I+\beta \left(I+K^{T}K-K-K^T\right),$$

(32)

where I is the identity matrix, the RLS estimates of a₁ and b₁ can be rewritten in a simpler form as:

$${\widehat{a}}_{1-\mathit{RLS}}=L^{-1}{\widehat{a}}_1,$$

(33)

$${\widehat{b}}_{1-\mathit{RLS}}=L^{-1}{\widehat{b}}_1.$$

(34)

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 K^TK is not in Toeplitz form. The elements of K^TK can be written as:

$${\left(K^{T}K\right)}_{\mathit{ij}}={\displaystyle \sum _{k\mathit{=}1}^M{K}_{ik}{K}_{kj}}={\displaystyle \sum _{k\mathit{=}1}^M{K}_{ik}{K}_{jk}.}$$

(35)

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, (K^TK) _{i j} = 0 if |i −j| > 2R + 2.

If we ignore the first and last 2R+2 rows of the matrix K^TK, it is a positive-definite symmetric Toeplitz matrix. Additionally, as all matrices in L are in Toeplitz form except for K^TK, the matrix L has the same properties as K^TK. 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 pO₂ images are being used. For a 640 × 480 pO₂ 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⁻¹ 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₂ 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 a₁ and b₁ of the original pO₂ map are found in matrix form. (D)$L_s^{-1}$, whose dimension is R _s ²× R _s ², is found for a small R _s × R _s pO₂ 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₁ and b₁ parameter maps to find approximate RLS estimates of the model parameters.

Figure 1

Flowchart of the proposed method.

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₂ 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 [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 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 [21], 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 [15]. Since the iterative method requires a considerable amount of time as size of the pO₂ map gets larger, the simulated pO₂ 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.

Figure 2

Simulated pO ₂ map (485 × 600 pixels). Color bar represents pO₂ values in millimeters of mercury.

Figure 3

Oxygen tension map (600 × 425 pixels) of a rat generated using the LS (1), iterative RLS (2), and the proposed RLS (3) methods. The color bar shows oxygen tension in millimeters of mercury.

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:

$$\delta ={\left(2+2\beta \left(l\mathit{+}4p\mathit{+}4q\right)\right)}^{-1}$$

(36)

Figure 4

MAE of the LS (er-ls), iterative RLS (er-rls-it), and proposed RLS estimation (er-rls-pr) methods.

Figure 5

Computation times of the iterative RLS (t-it) and proposed RLS (t-pr) estimation methods for different values of β.

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

Figure 6

MAE of the iterative RLS (er-it), and proposed RLS (er-pr) methods for different sizes of pO ₂ maps used in the estimation of ${\mathit{L}}_{\mathit{s}}^{-\mathbf{1}}$ .

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 7

Computation times of iterative RLS and proposed RLS estimation methods for different values of β for a real oxygen tension map.

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

Figure 8

The rectangular frame in Figure 1(1) and its estimates in the presence of noise with 20 dB SNR using the LS (2) and proposed RLS (3) methods. The color bar shows oxygen tension in millimeters of mercury.

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.