1.1.

Overview of Proposed Approach

In this paper, we propose a new strategy for retinal image deblurring where we consider the most general image degradation case: blurred retinal images acquired in different moments in time, ranging from minutes to months; hence, disease progression is also considered. The main reason for this general image degradation case that considers long time lapses comes from the potential need to restore a degraded image acquired in the past being the only one available at that stage of the disease. This problem arises quite often in clinical practice. A correct assessment of a patient's state evolution requires sharp images from all moments in time; the method proposed here enables such opportunity. Disease progression characterization is embedded in the algorithm with the identification of areas of structural change (see Sec. 3.3).

Our restoration method is based on a technique called blind deconvolution (BD).^{16, 17} The goal of BD is to recover the original scene from a single image or a set of blurred images in the presence of a poorly determined or unknown point-spread function (PSF). The main assumption is that the blur can be described by a convolution of a sharp image with the unknown PSF. Restoration by deconvolution improves contrast and resolution of digital images (i.e., it is easier to resolve and distinguish features in the restored image). To avoid confusion with super-resolution, we briefly describe what we mean by resolution improvement. Digital deconvolution can be described as any scheme that sharpens up the PSF, while the spatial frequency bandwidth remains unchanged. This means that the spatial frequency response and the two-point resolution is improved, but the cutoff frequency is unchanged;¹⁸ in the super-resolution context, the goal is to increase the cutoff frequency.

BD algorithms can be of single input [single-image blind deconvolution (SBD)] or of multiple images [multichannel blind deconvolution (MBD)]. Despite the fact that SBD is one of the most ill-posed problems, there are several reliable SBD algorithms,¹⁹ although most of them require that the blurred image be governed by relatively strong edges, which is not case here. In Sec. 4.1 we compare our approach to a recent state-of-the-art SBD method.²⁰ The computational overhead from MBD (all of the preprocessing to adjust the time-sequence of images) in comparison to SBD is practically negligible, and the robustness of MBD is far superior and worth applying because SBD fails to produce a valid restoration. By the same token, the additional processing enables the identification of structural changes in the retina over time—a central task in medical practice. As a result, we have chosen a multichannel approach for the restoration of blurred retinal images.

An overview of the proposed approach is described in Fig. 1. We consider as input two-color retinal images acquired with a conventional fundus camera within a time lapse that can span from several minutes to months given by routine patient checkups. The images correspond to the same retina but can differ with respect to illumination distribution, blur, and local structural changes given by pathological developments. These differences cannot solely be accounted for by the convolutional model described in Sec. 2. For that reason, the images must be preprocessed before the blind deconvolution stage can take place. We register the images and compensate for interimage illumination variation and structural changes. In fact, this preprocessing work becomes a great opportunity to meet one of the main concerns of ophthalmologists when they visually compare fundus images of the same retina over time: To identify true structural or morphological changes pertaining to possible pathological damage and, consequently, disregarding other changes merely caused by variation of illumination or blur. Ours is a two-stage blind deconvolution strategy. The first stage consists in the estimation of the PSFs following a multichannel scheme, and the second stage is the image deconvolution, where we restore every image with its corresponding PSF, independently. This has several advantages that will be explained in detail Sec. 3.5. The multichannel scheme is based on the method described in Ref. 21, which has proved to work well in practice with sufficient experimental data. It is an alternating minimization scheme based on a maximum a posteriori (MAP) estimation, with a priori distribution of blurs derived from the multichannel framework and a priori distribution of the ideal sharp image defined by regularization with the total variation of the image.²² MAP is formulated as an optimization problem, where regularization terms are directly related to priors. Regularization involves the introduction of additional information in order to solve an ill-posed problem in the form of a penalty or restriction in the minimization routine (see Sec. 3.4). This provides good quality of restoration—significantly better than, for example, Lucy–Richardson algorithm,²³ still widely used in biomedical applications. We have modified the algorithm in Ref. 21 to leave out regions where the eye fundus has structurally changed (it only takes into account one image in these regions) with the use of a masking operator, similarly to the solution proposed in Ref. 24 within the super-resolution context. This enabled us to restore both degraded input images.

Fig. 1

Block diagram illustrating the proposed method. [TeX:] \documentclass[12pt]{minimal}\begin{document}$\breve{z}_i$\end{document} ${ž}_i$ are the unregistered degraded input images, and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{u}_i$\end{document} ${û}_i$ are their restored versions. The other variables are intermediate outputs of every stage; their meaning is given in the text.

In this work, our novel contributions to the retinal image processing task are twofold. First, we propose a degradation model for time-series retinal images, which captures the underlying distortions resulting from instrument limitations and changes between patient visits; we are also able to identify and highlight such changes. Second, we propose a restoration strategy based on blind deconvolution that is able to obtain image enhancement and resolution improvement using inexpensive digital methods applied to images acquired with a conventional fundus camera.

3.1.

Image Registration

Image registration is a procedure that consists of spatial alignment of two or more images. General and application-specific image registration, such as in retinal imaging, has been investigated from the beginning of image-processing research. The interested reader is referred to the image registration review of Zitová and Flusser²⁵ and the recent work by Lee ²⁶ for objective validation of several retinal image registration algorithms. Image-registration techniques are usually divided into two groups: intensity-based and feature-based methods. Intensity-based methods have the drawback of poor performance under varying illumination conditions. Feature-based methods are robust to such effects but rely on accurate and repeatable extraction of the features. The retinal vasculature is known to provide a stable set of features for registration.

For registering the images, we use the robust dual-bootstrap iterative closest-point algorithm. We briefly describe it here; for a full description, of the method the reader is referred to Ref. 27. The vasculature from each image is automatically traced; starting from initial seed points extracted from a 1-D edge detection and, later, recursively tracking the vessels using directional templates. The vessel branching and crossover points are used as landmarks to register the images to subpixel accuracy. The registration algorithm starts from initial low-order estimates that are accurate only in small image regions called bootstrap regions. The transformation is then refined using constraints in the region, and the bootstrap region is expanded iteratively. The algorithm stops when the bootstrap region expands to cover the overlap between the images, and uses 12-dimensional quadratic mapping. This transformation model includes rotation, scale, translation, a shearing term, and a quadratic term that describes the spherical shape of the retina. We refer the interested reader to Ref. 28 for details on the model derivation. This registration algorithm is very robust to local changes and low overlap between images as demonstrated by its high success rate on test images with at least one common landmark point and overlaps even as low as 35%.²⁷ Even though the reported accuracy in Ref. 27 is of subpixel accuracy, in our case of degraded images this can be slightly worse without compromising the outcome. Minor local misregistration errors may occur when landmark points do not match precisely, but they will not be taken into account in the restoration because they will be masked out before the PSF estimation and image deconvolution stages (see Sec. 3.3).

To confirm the registration outcome, the pair of images before and after registration are shown in Fig. 3 in checkerboard representation, where the images are merged together in a chesslike pattern, where each square alternates information from one image to the other. Note how after registration the images have been correctly aligned, especially the blood vessel distribution.

Fig. 3

Registration of images from Fig. 2 in checkerboard representation. (a) Before and (b) after registration.

3.2.

Compensation of Uneven Illumination

Despite controlled conditions in retinal image acquisition, such as optical stops to prevent glares and provide a diffuse illumination, there are many patient-dependent aspects that are difficult to control and mainly affect the illumination component with gradual nonuniform spatial variations. Some of the contributing factors are (i) the curved surface of the retina (as a consequence, all regions cannot be illuminated uniformly); (ii) imaging requires either a naturally or an artificially dilated pupil (The degree of dilation is highly variable across patients); (iii) unexpected movements of the patient's eye; and (iv) presence of diseases. This nonuniform illumination across the image results in shading artifacts and vignetting. This effect hinders both quantitative image analysis and the reliable operation of subsequent global operators.

In our model, described by Eq. 1, the relative changes in intensity between the two fundus images cannot be described exclusively by convolution with different PSFs and must be compensated by k. A number of general-purpose techniques have been investigated to attenuate the variation of illumination. However, most techniques are oriented toward single-image compensation,² for instance, using the red channel to estimate background illumination.²⁹ Therefore, no consistency between two images is guaranteed. For our case, this uneven illumination can be compensated by properly adjusting the intensity values on one image to approximately match that of the other while satisfying a predetermined illumination model. This can be carried out if the blurring is not too large and the illumination changes smoothly, which is usually the case for fundus images. This assumption can be expressed mathematically as

Fig. 4

Illumination compensation function k(x, y).

3.3.

Segmentation of Areas with Structural Changes

The pathological region is actually a structural change and cannot be taken as a variation of illumination. Image change analysis is of interest in various fields and many algorithms have been developed for change detection.^{32, 33} A survey of change detection methods can be found in Ref. 34. An initial step in order to identify these changes comes from computing the difference from the two registered images including the illumination compensation as

Fig. 5

Intermediate outputs from the algorithm: (a) image difference Δz(x, y) in absolute value, (b) image difference in pseudocolor on top of gray-scale fundus image, and (c) mask m for avoiding areas with structural changes.

3.4.

Point-Spread Function Estimation

In this section, we describe the basic principles of the blind deconvolution method used for the estimation of the PSFs. For this purpose, we have chosen one of the best working MBD methods.²¹ MATLAB implementation of this method is available on the web of the authors.³⁶ The algorithm can be viewed as a Bayesian MAP estimation of the most probable sharp image and blur kernels. For our purposes, we used a modification of the original method that ignores regions affected by structural changes, which improves stability and precision of the computation. Without this modification, represented by the mask m in Eq. 5, the algorithm does not work reliably. The algorithm can be described as a minimization of the functional

The functional is alternately minimized in the subspaces corresponding to the image and the PSFs. The advantage of this scheme lies in its simplicity, this alternating minimization approach is actually a variation of the steepest-descent algorithm. The minimization in the PSF subspace is equivalent to the solution of a system of linear equations in the least-squares sense with the non-negativity constraint, in our implementation solved by the MATLAB fmincon function. The nonblind deconvolution realized by the minimization in the image subspace, is solved by half-quadratic iterative scheme,³⁷ replacing the total variation by [TeX:] \documentclass[12pt]{minimal}\begin{document}$\int \sqrt{|\nabla u|^2+\epsilon ^2}$\end{document} $\int \sqrt{{\left|\nabla u\right|}^2+{\epsilon }^2}$ , where ε is an auxiliary variable in the range 0 < ε ≪ 1. It is a small relaxation parameter that makes total variation differentiable around zero. A typical value for ε is 10⁻¹.

The main difference with respect to the original method²¹ is the introduction of the masking function m, which is computed in the beginning of the algorithm as described in Sec. 3.3. During the minimization, the multiplication by m is included in the operator corresponding to the convolution with u (in the PSF minimization step) and in the operator corresponding to the convolution with h ₂ (in the image minimization step). Because of the simplicity of this masking operation, the speed is practically the same as the speed of the original algorithm. In addition, even though we work with a complicated set of pixels, we can use the standard operation of convolution, which can eventually be speeded up using Fast Fourier transform (FFT).

3.5.

Image Restoration

The aim of our algorithm is to restore both images as much as possible. Note that from Eq. 5 the restored version of z ₁ ( [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{u}_1$\end{document} ${û}_1$ ) is obtained because z ₂ is masked; [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{u}_2$\end{document} ${û}_2$ could be obtained by minimizing Eq. 5 again with fixed PSFs and masking z ₁. This procedure has the disadvantage that both images are restored only within the common object field. Therefore, an appropriate solution is to restore each image z _i via single-channel deconvolution with their corresponding PSF h _i (estimated from the previous step) by the minimization of the functional

Finally, it should also be noted that the whole process of PSF estimation plus deconvolution can be computed for every channel of the RGB fundus image. However, in spite of the increase in computational burden, tests showed no real advantage to estimate the PSF for each channel. Moreover, the most suitable channel for PSF estimation is the green because it provides the best contrast. Whereas the blue channel encompasses the wavelengths most scattered and absorbed by the optical media of the eye; hence, the image has very low energy and a relatively high level of noise. As a result, the RGB deconvolved fundus image was computed by deconvolving every R, G, and B channel from the green channel PSF.

4.1.

Synthetic Images

In this section, we use synthetically degraded retinal images to test the performance of the proposed method. We use blurred signal-to-noise ratio (BSNR) to measure the noise contained in the degraded image, and improvement in signal-to-noise ratio (ISNR) to measure the quality of restored images.³⁸ They are defined as follows:

The first example is shown in Fig. 6, where the degraded images are synthesized from a sharp real image and the kernels shown in Fig. 6c and 6d plus Gaussian noise with zero mean and variance σ² = 10⁻⁶ (BSNR=40 dB). The recovered image and PSFs are shown in Fig. 7. The restoration provides an ISNR=4.45 dB. In this case, for synthetically degraded images the masking operation of Sec. 3.3 was not applied. Visual inspection of the details shown in Fig. 8 clearly reveal the accuracy of the method. Under these circumstances, the algorithm is able to produce a significant restoration of fine details like small blood vessels around the optic disc.

Fig. 6

(a, b) Degraded images (BSNR = 40 dB) and (c, d) PSFs.

Fig. 7

(a) Restored image (ISNR = 4.45 dB) and (b) Estimated PSFs.

Fig. 8

Details from (a) degraded image, (b) restored image, and (c) original image.

To further test our approach under a more realistic degradation, we produced an initial geometrical distortion, via a quadratic model^{26, 28} as the one used for registration (Fig. 9). After the geometric distortion, the degradation (blur plus noise) is produced on both images (BSNR=40 dB). They are then registered, and the restored image is recovered via MBD. The restored image and the estimated PSFs are shown in Fig. 10. The ISNR is slightly less (4.11 dB) than in the previous case, but still sufficient to produce a significant restoration. To corroborate our assumption that MBD methods seem better suited for this type of images, we tried to restore the image with a recent SBD method proposed in Ref. 20. The result is shown in Fig. 10e and visually reveals that it does not follow the true nature of the blurring with artifacts around the blood vessels, thus being prone to produce a poor restoration evidenced by an ISNR=−0.72 dB.

Fig. 9

(a) Original image and (b) geometrically distorted image.

Fig. 10

Image restoration from degraded and geometrically distorted images: (a) restored image by the proposed method (ISNR = 4.11 dB); (b) estimated PSFs; and (c) image detail, restored image by the method in Ref. 20 (ISNR = −0.72 dB); and (d) image detail.

Concerning parameter setting, in Fig. 11 we show the sensitivity of the two parameters λ_u and λ_h for the minimization of Eq. 5 in ISNR of the restored images. In Fig. 11a, we fix the value of λ_h to 10 and check the ISNR of the restored images for different initial values of λ_u = {10⁰, 10¹, 10², 10³, 10⁴, 10⁵}. The best restoration is obtained with λ_u = 10³; thus, in Fig. 11b we carried out the same procedure by fixing the value of λ_u to 10³ and checking the ISNR of the restored image for different values of λ_h = {1, 10, 20, 30, 40, 50}. The best restoration was obtained with an initial value of λ_h = 30. For this type of image, when scaled to the interval 〈0, 1〉, we find 20 < λ_h < 40 to be a suitable range to produce an optimal restoration.

Fig. 11

Test on parameter setting (BSNR = 40 dB). Average ISNR with respect to different initial values of (a) λ_u and (b) λ_h.

4.2.

Real Images

The experiments shown in this section aim to demonstrate the applicability of the proposed method for retinal image deblurring in real scenarios. Three different cases are shown in Fig. 12, including the retinal images that were used to illustrate the method (Fig. 2). The estimated PSFs are shown at the bottom of the restored images. All images contain some pathological damage and have been acquired within considerable lapses of time (several months). In all examples, the resolution improvement can be visually assessed by the clear distinction of details, such as small blood vessels or the increase in sharpness of edges, especially in the pathological areas. We emphasize the fact that these images correspond to real routine patient follow-up and were not intentionally degraded. From a clinical viewpoint, the enhancement can be used for a more precise assessment of a patient's state. Likewise, the images are more suitable for subsequent processing such as for the detection of retinal pathology.^{29, 39}

Fig. 12

Original and restored color retinal images; (a–c) indicate three separate cases arranged from left to right following our notation for degraded (z _i) and restored ( [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{u}_i$\end{document} ${û}_i$ ) images. The images are cropped to represent the region of interest given by the pathological area. The estimated PSF is shown at the bottom of the restored image. Video files are also included for change detection in cases (a) and (b). (Video 1, Quicktime, 0.5 MB); (Video 2, Quicktime, 0.4 MB)

In Fig. 13, the same images are shown but in gray scale to highlight the areas of structural change in pseudocolor. As mentioned earlier, this is an important result for its potential impact in medical practice. Subtle changes can be identified by this approach, such as the ones in Fig. 13b and the hemorrhage in the region of the optic disk in Fig. 13c. Another technique to rapidly identify changes from the two images is by alternating both restored images in a video sequence. Videos 1 and 2 (Fig. 12) correspond to the first two real cases.

Fig. 13

Visualization of structural changes in pseudo-color for the images of Fig. 12.