1.

Introduction

The ${\mathrm{SO}}_2$ in retinal vessels can provide doctors with important clinical information about the metabolic state of the retina.¹ For example, assessment of ${\mathrm{SO}}_2$ in the optic nerve head may help with early diagnosis of glaucoma and is crucial for timely effective treatment.² Monitoring ${\mathrm{SO}}_2$ changes in the retina could evaluate relationships between oxygen consumption, blood sugar levels, and vascular function in diabetic retinopathy.³ Therefore, it is important to obtain retinal oxygen saturation, which could be used to monitor retinal disorders and improve the efficiency and accuracy of diagnosis.

Among numerous ${\mathrm{SO}}_2$ measurement techniques in retinal vessels, dual-wavelength retinal oximetry based on fundus camera has been proven effective.⁴ It is based on the different optical properties of oxygenated hemoglobin (${\mathrm{HbO}}_2$) and hemoglobin (Hb). Several versions of the oximeters, which are based on the same principle, have been reviewed in literature.⁵ The system simultaneously collects images at two wavelengths: one is sensitive to changes of ${\mathrm{HbO}}_2$ in the blood and the other is an isosbestic wavelength for ${\mathrm{HbO}}_2$ and Hb. The ${\mathrm{SO}}_2$ values in the retinal vessels are finally measured after vessel segmentation, image registration, and calculation of the optical density ratio of two images. Figure 1 shows two retinal images of human eye with 45-deg view field acquired simultaneously at 570 and 600 nm. The 570-nm image, which is insensitive to ${\mathrm{HbO}}_2$ change, has a higher contrast between the blood vessels and the background than does the image at 600 nm.

Fig. 1

Dual-wavelength fundus images: (a) retinal image at 570 nm and (b) retinal image at 600 nm.

However, the aforementioned measurement method may lead to dramatic misinterpretations of the measured data.⁶ They are caused by many factors such as image noise, the imaging errors of the detection system, vessel reflex, fundus pigmentation extinction, and so on. In particular, noise can affect the subsequent feature extraction and calculation accuracy of ${\mathrm{SO}}_2$. Thamm et al.⁶ proposed a data processing method for the improvement of the primary information leading to a reduction of the uncertainties of the derived ${\mathrm{SO}}_2$. The approach may be the most accurate one because they took vessel reflex, imaging errors of the detection system, and the noise of the signal into account and developed a sophisticated model to simultaneously measure over a broad wavelength range using a spectrometer. But the major shortcoming is the limitation on measurements at a small region on the fundus including a retinal vessel and its surroundings. They could not give a complete, two-dimensional (2D) mapping of the ${\mathrm{SO}}_2$ in the retinal vascular tree, which is needed for clinical diagnostics. Hammer and Schweitzer⁷ used polarized light to enhance the reflection signal and improve the accuracy of oximetry. Actually, differences in choroidal pigmentation influence the measurement of ${\mathrm{SO}}_2$ and should be taken into account. Beach et al.⁸ have modified the calculation formula for ${\mathrm{SO}}_2$ to reduce the influence of choroidal pigmentation and have significantly reduced variation in the arterial ${\mathrm{SO}}_2$ measurements among subjects. In order to increase the measurement accuracy of ${\mathrm{SO}}_2$, Smith⁹ developed a series of oximetric equations that explicitly consider the effects of multiple light paths, including the effects of backscattering by red blood cells and lateral diffusion of light in the ocular fundus. Hammer et al.¹⁰ also took the influence of vessel diameter as well as fundus pigmentation on the ${\mathrm{SO}}_2$ into account and compensated for that by introducing linear compensation terms into the oximetry equation. Recently, Hardarson and Jonsson¹¹ studied whether and how image quality affects measurements of ${\mathrm{SO}}_2$ in retinal vessels. They employed a newly developed software tool to automatically grade the images on a scale of 0 to 1 according to the quality of the images. The quality grade was composed of the assessment of focus and contrast. They concluded that the poor image quality could lead to lower measured venous ${\mathrm{SO}}_2$ and could also affect measurements of arteries in more extreme cases. But they did not fully analyze the reasons for the decrease of image quality, nor did they address the noise problem of spectrophotometric images. It is necessary and important to study a denoising algorithm for the dual-wavelength images, since the image noise can affect the subsequent image processing and lead to substantial uncertainties in the calculated ${\mathrm{SO}}_2$.

In this study, we proposed a hybrid denoising method for dual-wavelength retinal images. First, noise model and parameter estimation were based on a mixed Poisson–Gaussian (MPG) noise model, then an MPG denoising algorithm was proposed based on variance stabilizing transform (VST) and improved dual-domain filter, which was called VST + dual-domain image denoising (DDID). The experimental results and quantitative analysis show that the proposed algorithm could more effectively remove the noise while preserving image details. The following simulation and analysis indicate that MPG noise in the retinal images could lead to erroneously low measured ${\mathrm{SO}}_2$ values. The denoised images could provide more accurate grayscale values for retinal oximetry.

This paper is organized as follows. Section 2 introduces the spectrophotometric retinal image. Section 3 describes the proposed denoising algorithm (VST + DDID). Section 4 presents the corresponding denoised results and evaluates the VST + DDID algorithm from two aspects, including image quality and the effect of noise on ${\mathrm{SO}}_2$ calculation. Finally, conclusions are drawn in Sec. 5.

2.

Spectrophotometric Retinal Image

The dual-wavelength image acquisition system developed by our lab consists of a commercial fundus camera, relay imaging optics (L1 and L2/L3), a beam splitter, two interference filters, and CCD cameras. The retinal image at the exit pupil of the fundus camera is relayed and separated by the beam splitter into two optical paths and simultaneously reimaged on each CCD camera at the wavelengths of 570 and 600 nm through narrow bandpass (bandwidth 10 nm) interference filters. The optical path of the dual-wavelength retinal image acquisition system can be seen in Fig. 2(a). The outputs of the CCDs are sent to a computer where the dual-wavelength images are stored. The images are postprocessed by ${\mathrm{SO}}_2$ computation software in the computer, and ${\mathrm{SO}}_2$ values are finally obtained after vessel segmentation, image registration, and calculation of optical density ratio of the two images. Figure 2(b) is a picture of the dual-wavelength retinal image acquisition system.

Fig. 2

(a) Schematic of the dual-wavelength retinal image acquisition subsystem. (b) Picture of dual-wavelength retinal image acquisition subsystem.

From the schematic of the dual-wavelength retinal imaging system, there exists a relatively large loss of light energy because the output retinal image of the fundus camera is separated by a $5:5$ beam splitter (50% transmission and 50% reflection) into two parts. Meanwhile, considering safety for the human eye, the maximum permissible exposure must be limited. All these factors make the dual-wavelength retinal images be captured under the condition of weak light. Dual-wavelength low-light retinal images are inevitably degraded by photon random noise, which shows graininess in the CCD image and follows Poisson noise distribution.¹² Due to cost considerations, industrial scan cameras are used; these cameras are plagued by readout noise which further degrades retinal image quality. In our retinal oximetry, the readout noise of CCD satisfies additive Gaussian distribution as determined by a camera performance evaluation. Because the calculation of ${\mathrm{SO}}_2$ is based on retinal image grayscale values, the MPG noise may dramatically affect subsequent image processing and lead to serious measurement error for ${\mathrm{SO}}_2$. A hybrid denoising method is developed for dual-wavelength retinal images denoising.

3.

Hybrid Denoising Algorithm

Although denoising algorithms designed for MPG noise have been proposed (e.g., Refs. 1314.15.–16), the removal of MPG noise is often performed by employing a VST such as the Anscombe transform.¹⁷ By VST, the MPG noise can be stabilized to an additive Gaussian noise. VST provides a simple but effective way to remove MPG noise. The denoising algorithm can be executed by four steps: (1) MPG noise model estimation, (2) apply VST to standardize the image noise, (3) denoise the image with an additive white Gaussian noise (AWGN) filter, and (4) return the image to its original range via inverse transformation.

In this paper, DDID is used in step 3, which was proposed by Knaus and Zwicker.¹⁸ DDID is a new class of image denoising methods which produce high-quality results and has been proven to be highly competitive for natural and grayscale images with state-of-the-art algorithms such as nonlocal bayes (NL-bayes)¹⁹ and block-matching and three-dimensional filtering (BM3D).²⁰ It is surprisingly simple to implement and has been successfully employed in applications such as synthetic aperture radar image²¹ and 3D video denoising.²² Several works such as nonlocal dual denoising²³ and data adaptive dual-domain denoising²⁴ have been proposed to improve the DDID. The proposed algorithm utilizes the advantages of both VST and DDID and is expected to remove MPG noise and contribute to a better preservation of details and edges in the dual-wavelength retinal images. To our knowledge, such studies in hybridization of VST and DDID for MPG noise reduction in digital image have not yet been reported.

3.2.

Dual-Domain Image Denoising

DDID is originally designed for AWGN. Its aim is to estimate the original image $x$ from a noise-contaminated image $y=x+\eta$ with a stationary variance ${\sigma }^2=\mathrm{Var}[\eta ]$. DDID splits the image into two layers and denoises them separately, then they are handled with a spatial method (bilateral filter) and a frequency-domain method [short-time Fourier transform (STFT)], respectively. The original image can be approximated as the sum of the two denoised layers:¹⁸

(4)

$$\tilde{x}=\tilde{s}+\tilde{S},$$

where $\tilde{s}$ and $\tilde{S}$ are the denoised high-contrast and low-contrast images. A bilateral filter is used in high-contrast image denoising, STFT is used in the low-contrast image denoising, and the bilateral filter protects large amplitudes of the signal, while wavelet shrinkage discards small amplitudes. Therefore, all amplitudes from large to small can be denoised by iterating. Figure 3 describes the block diagram of DDID. Starting with the noisy input $y$, the denoised result of an iteration is used to guide the subsequent iteration. It is worth mentioning that the noise image is used as the guide image in the first iteration. A more complete introduction for DDID is available in Knaus’s paper.¹⁸

Fig. 3

The block diagram for DDID.

The main drawback for DDID is that it produces typical frequency-domain artifacts, since the method itself was developed to avoid the artifacts of frequency-based methods.²³ The reason is that the guide image used in the first iteration is noisy, and the kernel in the bilateral filter is computed from it; “parasite” information is retained and propagated in the following iterations. Therefore, it is necessary to improve DDID in order to reduce the artifacts. In our study, we employed NL-bayes denoising algorithm to provide a clean guide for DDID.

3.3.

Variance Stabilizing Transform and Anscombe Transform

Once the parameters $a$ and $b$ in the MPG noise model are estimated, it is expected that a constant noise variance for every pixel is obtained by VST. One of the most popular VSTs is the generalized Anscombe transformation:²⁵

(5)

$$f(x)=\{\begin{array}{ll}\frac{2}{a}\sqrt{ax+\frac{3}{8}{a}^2+b},& x>-\frac{3}{8}a-\frac{b}{a}\\ 0,& x\le -\frac{3}{8}a-\frac{b}{a}\end{array}.$$

Applying Eq. (5) to MPG distributed data gives a signal whose noise is asymptotically additive standard normal. The original noise variance is stabilized to a unity. Now, the noise can be seen as AWGN with mean 0 and constant variance. Thus, the denoising methods based on the AWGN model can be applied, including DDID. The improved DDID was applied to the transformed image. The final denoised image can be obtained by inverse VST [$f_{\mathrm{VST}}^{-1}(x)$]. The whole procedure of our proposed method for MPG noise removal can be seen in Fig. 4.

Fig. 4

The block diagram of the proposed method (VST + DDID).

4.

Results and Evaluation

4.1.

Parameter Estimation Results for Dual-Wavelength Retinal Images

First, the parameters $a$ and $b$ in the MPG noise model are estimated for the dual-wavelength retinal images. Figure 5 shows the parameter estimation results for Fig. 1, in which the estimated sample points $\{\widetilde{x_i},\widetilde{{\sigma }_i}(x_i)\}$ are fitted into a line determined by the parameters $a$ and $b$ via maximum-likelihood curve fitting method.

Fig. 5

(a) Parameter estimation for Fig. 1(a). (b) Parameter estimation for Fig. 1(b).

According to Fig. 5, $a$ is 10 times greater than $b$ in each image, thus we can conclude that the influence of Poisson noise on retinal image is dramatically greater than Gaussian noise. We can also see that parameters $a$ and $b$ may vary at different wavelengths. This is because various wavelengths always cause different photon noises.²⁶

4.2.

Denoising Results for Dual-Wavelength Images

To evaluate the proposed scheme, we applied two approaches to dual-wavelength retinal images: VST + DDID (the version of our proposed algorithm) and VST + BM3D.²⁵ The two denoised results for Fig. 1 can be seen in Figs. 6 and 7, respectively.

Fig. 6

(a) Denoising result by VST + DDID for Fig. 1(a). (b) Denoising result by VST + DDID for Fig. 1(b).

Fig. 7

(a) Denoising result by VST + BM3D for Fig. 1(a). (b) Denoising result by VST + BM3D for Fig. 1(b).

4.3.

Evaluation

The proposed denoising algorithm is evaluated in two aspects, including image quality evaluation and evaluation of the effect of noise on calculation of ${\mathrm{SO}}_2$.

4.3.1.

Image quality evaluation

In order to evaluate the denoising algorithm, the power spectra in Figs. 1, 6, and 7 are calculated. The corresponding horizontal cross section of the power spectrum can be seen in Figs. 8(a)–8(f). The power spectra of images in Fig. 1 (original retinal images) are compared with the power spectra in Figs. 6 (VST + DDID) and 7 (VST + BM3D). It can be seen that at high frequencies, which probably correspond to noise, the power spectra of Figs. 8(c)–8(f) are decreased. The power spectra indicate that both VST + DDID and VST + BM3D can remove the Poisson noise and enhance the image contrast.

Fig. 8

(a) Horizontal cross section of the power spectrum of original retinal image at 570 nm; (b) horizontal cross section of the power spectrum of original retinal image at 600 nm; (c) horizontal cross section of the power spectrum of retinal image at 570 nm denoised by VST + DDID; (d) horizontal cross section of the power spectrum of retinal image at 600 nm denoised by VST + DDID; (e) horizontal cross section of the power spectrum of retinal image at 570 nm denoised by VST + BM3D; and (f) horizontal cross section of the power spectrum of retinal image at 600 nm denoised by VST + BM3D.

It is worth noting that there is edge blurring at low-contrast details in Fig. 7 (VST + BM3D) and no edge blurring in Fig. 6 (VST + DDID), as can be directly seen in the four amplified rectangular areas. Both VST + DDID and VST + BM3D can reduce noise and preserve features like edges; however, VST + BM3D has difficulty in preserving low-contrast details. In this aspect, VST + DDID is better than VST + BM3D, which is a state-of-the-art MPG denoising algorithm.

4.3.2.

Effect of noise on calculation of SO₂

The focus of our study is on the impact of MPG noise on the calculation of ${\mathrm{SO}}_2$. We assume that the results of blood vessel segmentation and image registration are perfect and have no effect on ${\mathrm{SO}}_2$. In order to objectively evaluate the influence, we first simulate artificial retinal images at 570 and 600 nm, with known blood vessel position, gray information, by which the ${\mathrm{SO}}_2$ can be directly obtained. The original artificial retinal image with a size of $512\times 512\text{pixels}$ and in a gray scale ranging from 0 to 255 can be seen in Fig. 9(a). Two veins and two arteries are the observable vessels (as large as 16 pixels in diameter) in the horizontal and vertical directions, respectively. The gray values in the radial direction of vessels are subjected to Gaussian distribution with a mean of 8.5 and variance of 8. The gray values in the axial direction of vessels are set to meet increment for ${\mathrm{SO}}_2$. Here, the axial direction is the direction of blood flow in a blood vessel, and the radial direction is substantially orthogonal to blood flow. The MPG noise with estimated parameters $a$ and $b$ was added in the original artificial retinal images, as shown in Fig. 9(b). We can see that the images contaminated by MPG noise have graininess which leads to detail loss and contrast reduction. The proposed algorithm was used for the artificial images contaminated with MPG noise. Figure 9(c) displays the denoised artificial images. The corresponding pseudo-color ${\mathrm{SO}}_2$ maps are shown in Fig. 10. Colors indicate ${\mathrm{SO}}_2$ in retinal vessels (scale to the right of the image). Veins can vary from blue to yellow but normally are light green, indicating ${\mathrm{SO}}_2$ approximately 50% to 60%. Arteries generally are orange to red, indicating ${\mathrm{SO}}_2$ approximately 80% to 100%.

Fig. 9

(a) The artificial retinal vessel images, (b) corresponding images contaminated by MPG noise, and (c) denoised images by the proposed method at 570 and 600 nm, respectively.

Fig. 10

Pseudo-color ${\mathrm{SO}}_2$ map (a) for the original artificial images, (b) for the images contaminated with MPG, and (c) for the images denoised by the proposed method, respectively.

Caused by MPG noise, ${\mathrm{SO}}_2$ values in both veins and arteries gravely deviate from the true values and are even smaller than the true values, as can be seen in Fig. 11. ${\mathrm{SO}}_2$ values in the denoised images are much closer to the true ${\mathrm{SO}}_2$ values. In addition, we employed mean square error (MSE) to further evaluate the results, as can be seen in Table 1:

(6)

$$\mathrm{MSE}=\frac{1}{N}\sum _{i=1}^N{[{\mathrm{S}}_{\text{original}}(i)-{\mathrm{S}}_{\text{measured}}(i)]}^2,$$

where $S_{\text{original}}$ represents ${\mathrm{SO}}_2$ values in the original image, and $S_{\text{measured}}$ represents ${\mathrm{SO}}_2$ values in the noisy image or denoised image. $N$ is the arterious or venous section numbers.

Fig. 11

(a) The ${\mathrm{SO}}_2$ values in all venous sections for Figs. 10(a)–10(c). (b) The ${\mathrm{SO}}_2$ values for all arterial sections for Figs. 10(a)–10(c).

Table 1

MSEs of SO2 for Figs. 10(a)/10(b) and 10(a)/10(c).

	Artery	Vein
${\mathrm{MSE}}_1$	0.0455	0.0074
${\mathrm{MSE}}_2$	0.0010	$1.2918\times {10}^{-4}$

${\mathrm{MSE}}_1$ is the MSE of ${\mathrm{SO}}_2$ between the original image and noisy image, while ${\mathrm{MSE}}_2$ is the MSE of ${\mathrm{SO}}_2$ between the original image and denoised image. We can find that ${\mathrm{MSE}}_2$ is significantly smaller than ${\mathrm{MSE}}_1$ in both artery and vein. From the above, we can conclude that compared with ${\mathrm{SO}}_2$ in artificial images corrupted by MPG noise, the proposed denoising algorithm can significantly improve the accuracy of ${\mathrm{SO}}_2$.

We give two pseudo-color fundus ${\mathrm{SO}}_2$ maps generated automatically by our ${\mathrm{SO}}_2$ calculation software, as shown in Fig. 12. They correspond to the original retinal images (Fig. 1) and the retinal images denoised by VST + DDID (Fig. 6), respectively. Colors indicate ${\mathrm{SO}}_2$ in retinal vessels (scale to the right of the image). Based on the above analysis, the values of ${\mathrm{SO}}_2$ are more accurate after denoising the dual-wavelength retinal images.

Fig. 12

Pseudo-color fundus ${\mathrm{SO}}_2$ map for the original retinal images and the retinal images denoised by VST + DDID, respectively.

A study in healthy volunteers, including 10 [6 males, 4 females, mean age: 22.9 years, standard deviation (SD): 1.7 years] subjects, was performed for evaluating the effect of noise on the reproducibility of the retinal oximetry. All procedures were in accordance with the tenets of the Declaration of Helsinki. Before their participation, signed informed consent was obtained from the subjects. Five fundus images of each subject were recorded. For each of the five images, two different radii of circles centered on optic disc center were created in order to reduce the influence of optic disc on the calculated ${\mathrm{SO}}_2$. The ${\mathrm{SO}}_2$ was determined for corresponding arterial and venous sections in the annular areas. The annular area for blood oxygenation analysis can be determined manually, and our ${\mathrm{SO}}_2$ calculation software can automatically calculate the ${\mathrm{SO}}_2$ values for each of the arterial and venous sections. In order to facilitate the observation of the results, average ${\mathrm{SO}}_2$ values for one arterial and one venous section in the annular areas are shown in Fig. 12. Table 2 shows the mean and SD of the five repeated measurements of the same vessel section. Table 3 shows the mean and SD of the five repeated measurements of the same vessel section after denoising. According to Tables 2 and 3, we can obtain a higher ${\mathrm{SO}}_2$ value using the denoised retinal images. This phenomenon agrees well with the simulation results in Sec. 4.3.2.

Table 2

The mean and SD of the five original images of the same individual.

	Artery	Vein
Mean (%)	90.32	52.68
	81.49 to 100.26	46.75 to 63.48
SD (%)	3.24	4.05
	1.57 to 5.89	1.81 to 6.48

Table 3

The mean and SD of the five denoised images of the same individual.

	Artery	Vein
Mean (%)	93.82	56.73
	84.76 to 102.69	50.28 to 68.85
SD (%)	2.57	3.58
	1.36 to 5.38	1.25 to 5.93

For 10 individuals, one artery and one vein were measured in each. Mean and SD were calculated respectively for each individual (from five images) as a measure of the reproducibility. Tables 2 and 3 show the means and ranges for these individual means and SDs.

5.

Conclusions

MPG noise exists in the dual-wavelength retinal images because of the maximal permissible exposure and light energy loss in the retinal oximetry. It affects further processing such as vessel segmentation and ${\mathrm{SO}}_2$ calculation. Inspired by the methodology of Refs. 15 and 18, we have proposed a hybrid method called VST + DDID for denoising MPG noise. VST + DDID can denoise the dual-wavelength retinal images and reduce the error in calculated ${\mathrm{SO}}_2$ for retinal vessels. To the best of our knowledge, this is the first time the image quality in retinal oximetry has been improved through noise estimation and a specific denoising method. The experiments with synthetic artificial images demonstrate that the proposed algorithm is effective in improving the ${\mathrm{SO}}_2$ calculation accuracy. We also demonstrate systematically higher saturation levels measured from a series of 10 healthy subjects using VST + DDID.

Our work was based on well-defined assumption that only MPG noise could affect the calculation of ${\mathrm{SO}}_2$. However, ongoing work is required to further improve the deviation of the calculated ${\mathrm{SO}}_2$. The deviation may be caused by nonuniformity of the CCD response, calibration, insufficient adjustment of the rotation and the lateral deviation between the dual-wavelength images, and so on.