1.1.

Contribution

In this work, we propose a method for removing blur from retinal images. We consider images degraded with SV blur, which may be due to factors like aberrations in the eye or relative camera-eye motion. Because restoring a single blurred image is an ill-posed problem, we make use of two blurred retinal images from the same eye fundus to accurately estimate the SV PSF. Before the PSF estimation and restoration stages take place, we preprocess the images to accurately register them and compensate for illumination variations not caused by blur, but by the lighting system of the fundus camera. This is depicted in the block diagram shown in Fig. 2. The individual stages of the method are explained in Sec. 3.

Fig. 2

Block diagram illustrating the proposed method. $z$ is the degraded image, $g$ is an auxiliary image of the same eye fundus used for the point-spread function (PSF) estimation, and $u$ is the restored image. The other variables are intermediate outputs of every stage; their meaning is given in the text.

We assume that in small image patches, the SV blur can be approximated by a spatially invariant PSF. In other words, that in a small region, the wavefront aberrations remain relatively constant; the so-called isoplanatic patch.⁶ An important aspect of our approach is that instead of deblurring each patch with its corresponding space-invariant PSF—and later stitching together the results—we sew the individual PSFs by interpolation and restore the image globally. This is intended to reduce some artifacts that otherwise would likely appear at the seams of the restored patches. The estimation of the local space-invariant PSFs may fail in patches with hardly any structural information (e.g., such as blood vessels). These poorly estimated or nonvalid PSFs introduce artifacts in the restored image. Detecting such artifacts and inferring the nonvalid PSFs is a difficult problem. Recently, Tallón et al.¹⁴ developed a strategy for detecting these patches in an SV deconvolution and denoising algorithm from a pair of images acquired with different exposures: a sharp noisy image with a short exposure and a blurry image with a long exposure. Because they had two distinct input images that were able to: (i) Identify patches where the blur estimates were poor based on a comparison (via a thresholding operation) of the deconvolved patches with the sharp noisy patches. (ii) In those patches, instead of correcting the local PSFs and deconvolving the patches again, they performed denoising in the noisy sharp image patch. The end result is a patchwork-image of deconvolved patches stitched together with denoised patches. Their method is mainly oriented at motion blur, this is the reason for a dual exposure strategy. This is not readily implementable in the retinal imaging scenario where the SV blur is generally caused by factors like aberrations, including those belonging to the patient’s eye optical system, the eye fundus shape, and the retinal camera. In this paper, we address the question “how to identify PSF estimation failure to improve the SV deconvolution of retinal images?” Retinal imaging provides a constrained imaging scenario from which we can formulate a restoration approach that incorporates prior knowledge of blur through the optics of the eye. The novelty in our approach is in the strategy based on eye-domain knowledge for identifying the nonvalid local PSFs and replacing them with appropriate ones. Even though methods for processing retinal images in a space-dependent way (like locally adaptive filtering techniques^15,16) have been proposed in the literature; to the best of our knowledge, this is the first time a method for SV deblurring of retinal images is proposed.

2.1.

Representation of SV PSF

An obvious problem of spatially varying blur is that the PSF is now a function of four variables. Except trivial cases, it is hard to express it by an explicit formula. Even if the PSF is known, we must solve the problem of a computationally efficient representation.

In practice, we work with a discrete representation, where the same notation can be used but with the following differences: the PSF $h$ is defined on a discrete set of coordinates, the integral sign in Eq. (1) becomes a sum, operator $H$ corresponds to a sparse matrix and $u$ to a vector obtained by stacking the columns of the image into one long vector. For example in the case of standard convolution, $H$ is a block-Toeplitz matrix with Toeplitz blocks and each column of $H$ corresponds to the same kernel $h(s,t)$.¹⁷ In the SV case that we address here, as each column of $H$ corresponds to a different position $(x,y)$, it may contain a different kernel $h(s,t,x,y)$.

In retinal imaging, all typical causes of blur change in a continuous gradual way,¹⁸ which is why we assume the blur to be locally constant. Therefore, we can make the approximation that locally the PSFs are space-invariant. By taking advantage of this property, we do not have to estimate local PSFs for every pixel. Instead, we divide the image into rectangular windows and estimate only a small set of local PSFs [see Fig. 3(a)] following the method described in Ref. 4 and outlined in Sec. 3. The estimated PSFs are assigned to the centers of the windows from where they were computed. In the rest of the image, the PSF $h$ is approximated by bilinear interpolation from the four adjacent local PSFs. This procedure is explained in further detail in the following section.

Fig. 3

(a) Retinal image degraded with artificial SV blur given by (b) grid of PSFs. The grid and the image patches shown in (a) are used for local PSF estimation.

3.1.

Preprocessing

Because we use a multichannel scheme for the estimation of the local PSFs, the images are preprocessed so that they meet the requirements imposed by the space-invariant convolutional model given by Eq. (2). This consists in registering the images and adjusting their illumination distribution following Ref. 4. By carrying out this procedure, the remaining radiometric differences between the images are assumed to be caused by blur and noise. Unlike the case considered in Ref. 4, where a relatively long lapse of time between the two image acquisitions may possibly involve a potential structural change, in this study, both the images $z$ and $g$ are originated from the same image or, in practice, they are acquired one shortly after the other. Therefore, no structural change is expected. Since image $g$ is registered and its illumination matched to $z$, we denote this transformed auxiliary image as $\widehat{g}$.

3.2.

Estimation of the Local PSFs

In Sec. 2, we described the model for a spatially varying blur in which we assume the PSF $h$ to vary gradually, which means that within small regions the blur can be locally approximated by convolution with a space-invariant PSF. For this reason, we approximate the global function $h$ from Eq. (1) by interpolating local PSFs estimated on a set of discrete positions. The main advantage of this approach is that the global PSF needs not be computed on a perpixel basis which is inherently time-consuming.

The procedure for estimating the local PSFs is the following. We divide the images $z$ and $\widehat{g}$ with a grid of $m\times m$ patches [Fig. 3(a)]. In each patch $p$, we assume a convolutional blurring model where an ideal sharp patch $u_p$ originates from two degraded patches $z_p$ and ${\widehat{g}}_p$ (for $p=1,\dots ,m\times m$). The local blurring model is

From this model, we can estimate the local PSFs with an alternating minimization procedure as described in Ref. 4 but applied locally. The general guideline is that the patch size should be large enough to include retinal structures and much larger than the size of the local PSF. In Sec. 4, we show further analysis on the robustness of the method to these parameters. Every local PSF is computed on each patch $p$ by minimizing the functional

3.3.

Identifying and Correcting Nonvalid PSFs

3.3.1.

Strategy based on eye-domain knowledge

The local PSF estimation procedure does not always succeeds. Consequently, such nonvalid PSFs must be identified, removed, and replaced. In our case, we replace them by an average of adjacent valid kernels. The main reason why the PSF estimation may fail is due to the existence of textureless or nearly homogenous regions bereft of structures with edges (e.g., blood vessels) to provide sufficient information.¹⁴ To identify these nonvalid PSFs, we devised an eye-domain knowledge strategy. The incorporation of proper a priori assumptions and domain knowledge about the blur into the method provides an effective mechanism for a successful identification of poorly estimated PSFs.

The optics of the eye is part of the imaging system, therefore it is reasonable to assume that the PSF of the imaging system is determined by the PSF of the eye. The retinal camera can indeed be close to diffraction limited with a very narrow PSF, but the optics of the eye is governed by optical aberrations that change across the visual field¹⁸ that lead to an SV PSF. The typical PSFs of the human eye, as reported in the literature,^18,20 display distinct shapes in many cases displaying long tails’ evidence of the inherent optical aberrations. What is common to all PSFs of the human eye is that the energy is spread from a central lobe and decreases outwardly. In this paper, we assume that PSFs that do not follow this general pattern are nonvalid PSFs which correspond to patches where the estimation failed.

In order to prove this, we designed an artificial experiment where we compare the estimated PSFs with the ground-truth PSFs by using a kernel similarity measure $S$ proposed by Hu and Yang.²¹ The measure is based on the peak of the normalized cross correlation between the two PSFs. The authors showed in the paper that this measure is more accurate than the root mean square error, especially because it is shift invariant. The measure is defined as the blur kernel similarity $S(h,\tilde{h})$ of two kernels, $h$ and $\tilde{h}$,

Eq. (4)

$$S(h,\tilde{h})=\underset{\gamma }{\max }\rho (h,\tilde{h},\gamma ),$$

where $\rho (.)$ is the normalized cross-correlation function and $\gamma$ is the possible shift between the two kernels. Let $\tau$ represent element coordinates, $\rho (.)$ is given by

Eq. (5)

$$\rho (h_p,{\tilde{h}}_p,\gamma )=\frac{\sum _{\tau }{h}_p(\tau )·{\tilde{h}}_p(\tau +\gamma )}{\parallel h_p\parallel ·\parallel {\tilde{h}}_p\parallel },$$

where $\parallel ·\parallel$ is the ${\ell }_2$-norm. Larger similarity values reflect more accurate PSF estimation, thus better image restoration. The graphical representation of the similarity measure $S$ is shown in Fig. 4(a). The darkest squares correspond to PSFs with the lowest similarity score.

Fig. 4

(a) PSF similarity measure. (b) Characterization of estimated local PSFs by energy distribution. Histograms plotted in white bars in (b) have been labeled as valid PSFs based on the histogram probability descriptor ($P$) described in the text.

The way we determine a nonvalid PSF is by characterizing the energy distribution along the local PSF space. We add the PSF values along concentric squares of radius $r$ to build an energy distribution histogram $f(r)$, which is normalized to sum to 1. In Fig. 4(b), we show the histograms for the energy distribution characterization of the estimated PSFs from the artificially degraded retinal image [Fig. 3(a)]. To identify nonvalid PSFs, i.e., PSFs that do not follow the pattern we have previously described, we compute a shape descriptor for $f(r)$ defined as the probability $P(r\le r_m)$,

Eq. (6)

$$P(r\le r_{\mathrm{m}})=\sum _0^{r_{\mathrm{m}}}f(r)·r,$$

where $r_{\mathrm{m}}$ is the mode of $f(r)$. This descriptor gives an indication of how much a histogram is spread in relation to the peak ($r_{\mathrm{m}}$) of the histogram. Particularly, it yields low values when the mode is located in the first few bins of the histogram and large values otherwise. This descriptor can be correlated with the PSF similarity measure to determine how it discriminates valid and nonvalid PSFs.

In Fig. 5, we plot the similarity measure $S$ against $P(r\le r_{\mathrm{m}})$ for the estimated PSFs from two different retinal images which have been artificially degraded. The correlation coefficient for the two variables is $-0.64$. In addition to the correlation, from the plot we note that most PSFs with high $P(r\le r_{\mathrm{m}})$ have low kernel similarity ($S$) values, which is an indication that these are nonvalid or poorly estimated PSFs. A machine learning algorithm or clustering technique that automatically classifies nonvalid PSFs is out of the scope of this paper. Instead, our aim is to show that the energy characterization approach is sufficient to identify nonvalid PSFs, even at the expense of a few valid ones. As we show in Sec. 4.1, this is not critical because the PSF changes smoothly throughout the FOV.

Fig. 5

PSF similarity measure ($S$) versus histogram probability ($P$). An $S$ value of 1 represents an accurately estimated kernel. Low $P$ values are associated with correctly estimated kernels.

The histograms plotted in white bars in Fig. 4(b) correspond to PSFs with $P(r\le r_{\mathrm{m}})<0.3$, which we have labeled as valid PSFs. This means that we are favoring histograms skewed toward the left side. This is based on our assumption of the pattern that valid PSFs should follow. This correlates well with the similarity measure. Note that most PSFs with low similarity measure [darkest squares in Fig. 4(a)] have been correctly identified as nonvalid PSFs [black histograms in Fig. 4(b) and PSFs labeled with boxes in Fig. 6(a)].

Fig. 6

(a) Estimated set of $5\times 5$ local PSFs. (b) New set of local PSFs with nonvalid PSFs replaced [compared with ground-truth set Fig. 3(b)]. Nonvalid PSFs have been labeled with a red square.

The procedure for correcting the nonvalid local PSFs consists of replacing them with the average of adjacent valid kernels. Without this correction, the reconstruction develops ringing artifacts [see for example Fig. 11(b)]. The new set of valid local PSFs after replacing the nonvalid ones for the artificially degraded image is shown in Fig. 6(b).

3.4.

PSF Interpolation

The computation of the SV PSF $h$ is carried out by interpolating the local PSFs estimated on the regular grid of positions. The PSF values at intermediate positions are computed by bilinear interpolation of four adjacent known PSFs,²² as shown in Fig. 7. Indexing any four adjacent grid points as $p=1,\dots ,4$ (starting from the top-left corner and continuing clockwise), the SV PSF in the position between them is defined as

Fig. 7

Because the blur changes gradually, we can estimate convolution kernels on a grid of positions and approximate the PSF in the rest of the image (bottom kernel) by interpolation from four adjacent kernels.

In light of the definition for an SV PSF in Eq. (7), we can compute the convolution of Eq. (1) as a sum of four convolutions of the image weighted by coefficients ${\alpha }_p(x,y)$

In the same fashion, the operator adjoint to $H$ (SV counterpart of correlation) denoted by $H^{*}$ can also be defined in terms of the sums of four convolutions weighted by the ${\alpha }_p$ coefficients. These two operators are needed in all first-order minimization algorithms as the one used in the restoration stage (see Ref. 23 for further details).

3.5.

Restoration

Having estimated a reliable SV PSF, we proceed to deblur the image. Image restoration is typically formulated within the Bayesian paradigm, in which the restored image is sought as the most probable solution to an optimization problem. The restoration can be described as the minimization of the functional

The functional of Eq. (15) bounds the original function in Eq. (14) and has the same value and gradient in the current $u_i$, which guarantees convergence to the global minimum. To solve Eq. (15), we used the conjugate gradient method.¹⁷ The initial value of $u_i$ for $i=0$ is set to be equal to $z$. In order to avoid numerical instability for areas with small gradient ($|\nabla u_i|$ approaching zero), we use a relaxed $ϵ$-form of the minimized functional in Eq. (14), which implies that $|\nabla u_i|<ϵ$ is equal to $ϵ$.

As regards the restoration of color RGB retinal images, we consider the following. The most suitable channel for PSF estimation is the green because it provides the best contrast.²⁶ This is mainly due to the spectral absorption of the blood in this band, which yields the dark and well contrasted blood vessels.²⁷ Conversely, the blue channel encompasses the wavelengths most scattered and absorbed by the optical media of the eye,²⁸ therefore the image in this band has very low energy and a relatively high level of noise. In the spectral zone of wavelengths larger than 590 nm, the light scattering on the red blood cells and the light reflection from the eye structures behind the vessel are dominant.²⁹ This produces the red band to be saturated and of poor contrast. As a result, we estimate the SV PSF from the green channel of the RGB color image, and later deconvolve every R, G, and B channels with the estimated SV PSF to obtain a restored RGB color image.

4.1.

Artificially Degraded Images

For the artificial experiment, we take a pair of images and degrade them with a $5\times 5$ grid of realistic PSFs plus Gaussian noise ($\sigma ={10}^{-6}$). The grid of PSFs was built upon realistic PSFs estimated from real degraded retinal images following the approach of Ref. 4. From the two input images, we restore one, and the other is used exclusively for the purpose of PSF estimation. We estimate the local PSFs by dividing the image into overlapping patches on a $5\times 5$ grid [as shown in Fig. 3(a)]. The estimated PSFs are shown in Fig. 6(a). Because the PSF estimation may fail, we identify the nonvalid PSFs as described in Sec. 3.3. We replace them with with the average of adjacent valid kernels.

In Fig. 8(a), we show the restored artificial image with the directly estimated PSFs. The effect of nonvalid PSFs is evident in the poor quality of the restoration and the ringing artifacts. In Fig. 8(b), we show the restoration with the proposed method, where the nonvalid PSFs have been identified, removed, and replaced by the average of adjacent PSFs. To evaluate the restoration, we use the cumulative error histogram on a patch basis. The error⁵ is the difference between a recovered image $I_r$ with the estimated kernels and the known ground-truth sharp image $I_g$ over the difference between the deblurred image $I_{kg}$ with the ground-truth kernels. The error is given by $\parallel I_r-I_g\parallel /\parallel I_{kg}-I_g\parallel$. In Fig. 8(c), we show the cumulative error histogram for three restorations. H1 is the restoration with the directly estimated PSFs. It is important to note that shifted local PSFs warp the image which introduce additional artifacts and is the reason for such a low performance with approximately 40% of patches with an error lower than 2.5. After shifting the centroid of the PSFs to the geometrical center [restoration H2 in Fig. 8(c)], the reconstruction error is reduced significantly, about 60% of patches have an error lower than 1.5. Finally, the restoration (H3) with the removal of nonvalid PSFs increases significantly with all patches now displaying an error lower than 1.5. This means that after the nonvalid PSFs have been replaced the restoration quality is significantly increased.

Fig. 8

(a) Restoration with the set of directly estimated PSFs shown in Fig. 6(a) (notice the artifacts due to nonvalid PSFs) and (b) restoration with the new set of PSFs with nonvalid PSFs replaced shown in Fig. 6. (c) Error histogram for evaluating the reconstruction using: H1—directly estimated PSFs, H2—PSFs shifted toward the geometrical center, and H3—the new set of valid PSFs.

To determine the limitations and robustness of the proposed method, we carried out several tests. First, we need to determine an optimal patch size for accurately estimating the local PSF. A patch that is too small compared to the kernel size may not have enough information and is likely to favor the trivial solution (convolution with a delta). Conversely, a patch that is bigger than necessary may hinder the SV approach in addition to increasing the computational burden. In Fig. 9(a), we show the mean similarity measure versus the image patch size. It is important to note that a patch size of roughly four times the kernel size may be sufficient for an accurate PSF estimation. Increasing too much the patch size hinders the PSF estimation.

Fig. 9

(a) Mean similarity measure versus patch size. (b) Mean reconstruction error versus local PSF size. (c) Reconstruction error versus PSF grid size.

The second aspect to consider is under- and over-estimation of the local PSF size. In Fig. 9(b), we show the mean patch reconstruction error versus PSF size. As expected the error is minimum when the size is similar to the true size, yet under-estimating the PSF size is worse than over-estimating.

In relation to the SV characterization of the blur, we performed the PSF estimation with a different number of kernels. Initially with a coarse $2\times 2$ grid increasing up to a $10\times 10$ fine grid (the $5\times 5$ grid is the ground-truth). A similar behavior is observed. Under-estimating the proper PSF grid size has a negative effect in that the SV nature of the blur is hardly identified which yields an error above 6 for the whole image. In this case, the error is not computed on a patch-basis because of the variable grid size.

4.2.

Naturally Degraded Images

All of the naturally degraded images used in the experiments were acquired in pairs, typically with a time span between acquisitions of several minutes. Initially, to show the limits of the space-invariant approach we restored the blurred retinal image from Fig. 10(a) with a single global PSF with the space-invariant method we proposed in Ref. 4. This image corresponds to the eye fundus of a patient with strong astigmatism, which induces an SV blur as depicted by the image details shown in Figs. 10(b)–10(e). The restoration is shown in Fig. 11(a) and we can clearly observe various artifacts despite an increase in sharpness in a small number of areas. In view of this, it is evident that the space-invariant assumption does not hold in such cases. In the following, we move to the SV approach.

Fig. 10

(a) Retinal image degraded with real SV blur given by strong astigmatism. (b), (c), (d), and (e) zoomed regions to show the SV nature of the blur.

Fig. 11

(a) Space-invariant restoration, (b) SV restoration with directly estimated PSFs, and (c) SV restoration with the new set of PSFs. The reader is strongly encouraged to view these images in full resolution at http://www.goapi.upc.edu/usr/andre/sv-restoration/index.html.

To carry out the SV restoration, we estimated the local PSFs on a $5\times 5$ grid of image patches. From the estimated PSFs shown in Fig. 12(a), we notice a clear variation in shape mainly from the top-right corner where they are quite narrow, to the bottom left corner where they are more spread and wide. This variation is consistent with the spatial variation of the blur observed in the retinal image of Fig. 10(a). We restored the image with these local PSFs that were estimated directly without any adjustment. The restored image is shown in Fig. 11(b). One immediately obvious feature is that in several areas the restoration is rather poor, displaying ringing artifacts, whereas in others it is to some extent satisfactory. The local poor-restoration is linked to areas where the PSF estimation failed. By removing and correcting these nonvalid local PSFs, we obtained a noteworthy restoration shown in Fig. 11(c). Notice the overall improvement in sharpness and resolution with small blood vessels properly defined as shown by the image-details in the third column of Fig. 13. It could be said that without the replacement of the nonvalid PSFs the image quality after restoration is certainly worse than the original degraded image (see second column of Fig. 13).

Fig. 12

(a) Estimated set of local PSFs from naturally degraded retinal image shown in Fig. 10(a). (b) Set of local PSFs after replacing nonvalid PSFs. Nonvalid PSFs have been labeled with a red square.

Fig. 13

Details of restoration. From left to right: the original degraded image, the SV restoration without correction of PSFs and the SV restoration with the correction.

To further demonstrate the capabilities of our method, additional restoration results on real cases from the clinical practice are shown in the following figures. As we mentioned in Sec. 1, a typical source of retinal image degradation comes from patients with corneal defects in which the cornea has an irregular structure [Fig. 1(a)]. This induces optical aberrations, which are mainly responsible for the SV blur observed in the retinal image. The image details shown in Fig. 1(b) reveal a significant improvement in which the retinal structures are much sharper and enhanced. In Fig. 14(a), a full color retinal image is shown, in which three small hemorrhages are more easily discernible in the restored image, along with small blood vessels. Another retinal image, shown in Fig. 14(b), reveals a clear improvement in resolution with a much finer definition of blood vessels.

Fig. 14

Other retinal images restored with the proposed method. (a) First row: original and restored full-size retinal images. (b) Second and third rows: image details.

In addition, we processed retinal angiography images to test our method against a different imaging modality. Ocular angiography is a diagnostic test that documents, by means of photographs, the dynamic flow of dye in the blood vessels of the eye.³⁰ The ophthalmologists use these photographs both for diagnosis and as a guide to patient treatment. Ocular angiography differs from fundus photography in that it requires an exciter–barrier filter set (for further details see Ref. 30). The retinal angiography shown in Fig. 15 is degraded with a mild SV blur that hinders the resolution of small—yet important—details. The restoration serves to overcome this impediment; this can be observed from the zoomed-detail of the restored image. The image enhancement may be useful for the improvement of recent analysis techniques for automated flow dynamics and identification of clinical relevant anatomy in angiographies.³¹

Fig. 15

Restoration of a retinal angiography. First row: original and restored full retinal images. Second row: image details.

Finally, another way to demonstrate the added value of deblurring the retinal images is to extract important features, in this case detection of blood vessels. Such a procedure is commonly used in many automated disease detection algorithms. The improvement in resolution paves the way for a better segmentation of structures with edges. This is in great part due to the effect of the total variation regularization because it preserves the edge information in the image. To carry out the detection of the retinal vasculature, we used Kirsch’s method.³² It is a matched filter algorithm that computes the gradient by convolution with the image and eight templates to account for all possible directions. This algorithm has been widely used for detecting the blood vessels in retinal images.³³ In Fig. 16, we show the detection of the blood vessels from a real image of poor quality image and its restored version using our proposed method. A significant improvement in blood vessel detection is achieved. Smaller blood vessels are detected in the restored image, whereas the detection from the original image barely covers the main branch of the vasculature.

Fig. 16

First row (from left to right): original and restored retinal images. Second row: detection of blood vessels. Notice how the small blood vessels are better detected in the restored image.