## 1.

## Introduction

The ${\mathrm{SO}}_{2}$ in retinal vessels can provide doctors with important clinical information about the metabolic state of the retina.^{1} 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.^{2} Monitoring ${\mathrm{SO}}_{2}$ changes in the retina could evaluate relationships between oxygen consumption, blood sugar levels, and vascular function in diabetic retinopathy.^{3} 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.^{4} 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.^{5} 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.

However, the aforementioned measurement method may lead to dramatic misinterpretations of the measured data.^{6} 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.^{6} 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^{7} 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.^{8} 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^{9} 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.^{10} 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^{11} 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.

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.^{12} 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.^{17} 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.^{18} 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)^{19} and block-matching and three-dimensional filtering (BM3D).^{20} It is surprisingly simple to implement and has been successfully employed in applications such as synthetic aperture radar image^{21} and 3D video denoising.^{22} Several works such as nonlocal dual denoising^{23} and data adaptive dual-domain denoising^{24} 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.1.

### Mixed Poisson–Gaussian Noise Model Estimation

The MPG noise model for dual-wavelength retinal images is described as

where $x$ and $y$ are the original and noised signals, respectively. ${n}_{P}$ is the signal-dependent Poisson noise component satisfying where $a>0$ is the parameter related to the noise level of the Poisson component. According to the elementary properties of the Poisson distribution, we can obtain the mean and variance of ${n}_{p}(x)$.^{15}It follows that $E\{{n}_{p}(x)\}=0$ and $\mathrm{var}\{{n}_{p}(x)\}=ax$. ${n}_{G}$ is the additive zero-mean Gaussian noise component, for which $b$ is the signal-independent variance. Thus, the total noise variance can be derived as

^{15}The parameters $a$ and $b$ for the MPG model can be estimated from the homogeneous regions in the observed image.

^{15}At first, the image is transformed to the wavelet domain and then segmented into nonoverlapping level sets, where each level set stands for a homogeneous region. The signal and noise components are approximately extracted in each homogeneous region. Then the mean value of signal ${x}_{i}$ and the variance of noise ${\sigma}^{2}({x}_{i})$ in each region are estimated. Finally, a curve fitting method is used to estimate the parameters $a$ and $b$. Following the method proposed in Ref. 15, wavelet transform is used for signal and noise extraction. The 2D discrete wavelet transform (2D-DWT) generates the approximation coefficient matrix WA and three detail coefficient matrices WH, WV, and WD for horizontal, vertical, and diagonal directions, respectively. For a homogeneous region, WA is considered as the coefficients of the signal component, while WD is regarded as the coefficients of the noise component. The whole estimation procedure can be summarized as Algorithm 1. For dual-wavelength retinal images, $a$ and $b$ are independently estimated for each image.

## Algorithm 1

MPG noise model estimation.

Input: Observed image $I$. |

Output: Estimated parameters $\hat{a}$ and $\hat{b}$ of MPG noise model. |

1 Perform 2D-DWT on the observed image $I$. |

2 Remove the pixels with large gradient in WA domain so that the smooth regions remain. |

3 Segment the smooth region into $N$ level sets ${S}_{1},{S}_{2},\dots ,{S}_{N}$, where each level set ${S}_{j}$ stands for a homogeneous region. |

4 For each homogeneous region ${S}_{j}$, estimate ${\tilde{x}}_{j}$ from WA and estimate $\tilde{{\sigma}_{j}}({x}_{j})$ from WD. |

5 Estimate $a$ and $b$ by curve fitting method according to Eq. (3). |

## 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:^{18}

^{18}

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.^{23} 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:^{25}

## (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}.$$## 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 $\{\tilde{{x}_{i}},\tilde{{\sigma}_{i}}({x}_{i})\}$ are fitted into a line determined by the parameters $a$ and $b$ via maximum-likelihood curve fitting method.

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.^{26}

## 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.^{25} The two denoised results for Fig. 1 can be seen in Figs. 6 and 7, respectively.

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

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_{2}

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{\hspace{0.17em}\hspace{0.17em}}\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%.

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},$$## 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.

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.

## Acknowledgments

This work was supported by National Scientific Instruments and Equipment Development Special Foundation of China (No. 2013YQ49085903) and Science Foundation of Tianjin Municipal Bureau of Health in China (No. 2013KR19).

## References

## Biography

**Yong-Li Xian** received her master’s degree from the University of Science and Technology of China in 2011. She is now a doctoral candidate in the School of Optoelectronic Information at the University of Electronic Science and Technology of China. Her current research interests include biomedical optics and image processing.

**Yun Dai** received his PhD from the Institute of Optics and Electronics, Chinese Academy of Sciences in 2005. He is now a researcher in the Institute of Optics and Electronics, Chinese Academy of Sciences, Chengdu. His main research interests include adaptive optics, visual optics, and biomedical optics and instrumentation.

**Chun-Ming Gao** received his PhD from Nan Jing University in 2004. He is now a professor in the School of Optoelectronic Information at the University of Electronic Science and Technology of China. His main research interests include optical, acoustic, and photoacoustic technology.

**Rui Du** received his BSc and PhD degrees from Huazhong University of Science and Technology in 2006 and 2011, respectively. He is now an associated researcher in the Institute of Optics and Electronics, Chinese Academy of Sciences, Chengdu. His main research interests include biomedical optics and instrumentation.