Blind Image Deblurring Driven by Nonlinear Processing in the Edge Domain

This work addresses the problem of blind image deblurring, that is, of recovering an original image observed through one or more unknown linear channels and corrupted by additive noise. We resort to an iterative algorithm, belonging to the class of Bussgang algorithms, based on alternating a linear and a nonlinear image estimation stage. In detail, we investigate the design of a novel nonlinear processing acting on the Radon transform of the image edges. This choice is motivated by the fact that the Radon transform of the image edges well describes the structural image features and the e ﬀ ect of blur, thus simplifying the nonlinearity design. The e ﬀ ect of the nonlinear processing is to thin the blurred image edges and to drive the overall blind restoration algorithm to a sharp, focused image. The performance of the algorithm is assessed by experimental results pertaining to restoration of blurred natural images.


INTRODUCTION
Image deblurring has been widely studied in literature because of its theoretical as well as practical importance in fields such as astronomical imaging [1], remote sensing [2], medical imaging [3], to cite only a few. Its goal consists in recovering the original image from a single or multiple blurred observations.
In some application cases, the blur is assumed known, and well-known deconvolution methods, such as Wiener filtering, recursive Kalman filtering, and constrained iterative deconvolution methods, are fruitfully employed for restoration.
However, in many practical situations, the blur is partially known [4] or unknown, because an exact knowledge of the mechanism of the image degradation process is not avail-able. Therefore, the blurring action needs to be characterized on the basis of the available blurred data, and blind image restoration techniques have to be devised for restoration. These techniques aim at the retrieval of the image of interest observed through a nonideal channel whose characteristics are unknown or partially known in the restoration phase. Many blind restoration algorithms have been proposed in the past, and an extended survey can be found in [5,6].
In some applications, the observation system is able to give multiple observations of the original image. In electron microscopy, for example, many differently focused versions of the same image are acquired during a single experiment, due to an intrinsic tradeoff between the bandwidth of the imaging system and the contrast of the resulting image. In other applications, such as telesurveillance, multiple observed images can be acquired in order to better counteract,  in the restoration phase, possible degradation due to motion, defocus, or noise. In remote sensing applications, by employing sensor diversity, different versions of the same scene can be acquired at different times through the atmosphere that can be modeled as a time-variant channel. Different approaches have been proposed in the recent past to face the image deblurring problem. In [7], it is shown that, under some mild assumptions, both the filters and the image can be exactly determined from noise-free observations as well as stably estimated from noisy observations. Both in [7,8], the channel estimation phase precedes the restoration phase. Once the channel has been estimated, image restoration is performed by subspace-based and likelihood-based algorithms [7], or by a bank of finite impulse response (FIR) filters optimized with respect to a deterministic criterion [8].
Different approaches resort to suitable image representation domains. To cite a few, in [9], a wavelet-based edge preserving regularization algorithm is presented, while in [10], the image restoration is accomplished using simulated annealing on a suitably restricted wavelet space. In [11], the authors make use of the Fourier phase for image restoration [12] applying appropriate constraints in the Radon domain.
In [13,14], the authors resort to an iterative algorithm, belonging to the class of Bussgang algorithms, based on alternating a linear and a nonlinear image estimation stage. The nonlinear estimation phase plays a key role in the overall algorithm since it attracts the estimate towards a final restored image possessing some desired structural or statistical characteristics. The design of the nonlinear processing stage is aimed at superimposing the desired characteristics on the restored image. While for class of images with known probabilistic description, such as text images, the nonlinearity design can be conducted on the basis of a Bayesian criterion, for natural images, the image characterization and hence the nonlinearity design is much more difficult. In [14], the authors design the nonlinear processing in a transformed domain that allows a compact representation of the image edges-the edge domain.
In this paper, we investigate the design of the nonlinear processing stage using the Radon Transform (RT) [15] of the image edges. This choice is motivated by the fact that the RT of the image edges well describes the structural image features and the effect of blur, thus simplifying the nonlinearity design.
The herein discussed approach shares some common points with [16] since it exploits a compact multiscale representation of natural images.
The structure of the paper is as follows. In Section 2, the observation model is described. Following the recent literature, a multichannel approach is pursued. Section 3 recalls the basic outline of the Bussgang algorithm, which is described in detail in the appendix. Section 4 is devoted to the description of the image edge extraction as well as to the discussion of the nonlinearity design in the edge domain. Section 5 presents the results of the blind restoration algorithm and Section 6 concludes the paper.

THE OBSERVATION MODEL
The single-input multiple-output (SIMO) observation model of images is represented by M linear observation filters in presence of additive noise. This model, depicted in Figure 1, is given by  processes, that is, Here, E{·} represents the expected value, δ[·] the unit sample, and δ[·, ·] the bidimensional unit sample.

MULTICHANNEL BUSSGANG ALGORITHM
The basic structure of one step of the iterative Bussgang algorithm for blind channel equalization [17,18,19], or blind image restoration [14,20], consists of a linear filtering of the measurements, followed by a nonlinear processing of the filter output, and concluded by updating the filter coefficients using both the measurements and the output of the nonlinear processor. The scheme of the iterative multichannel Bussgang blind deconvolution algorithm, as presented in [14], is depicted in Figure 2. The linear restoration stage is accomplished using a bank of FIR restoration filters f (k) i [m, n], i = 0, . . . , M −1, with finite support of size (2P +1)×(2P +1), namely,

BUSSGANG NONLINEARITY DESIGN IN THE EDGE DOMAIN USING THE RADON TRANSFORM
The quality of the restored image obtained by means of the Bussgang algorithm strictly depends on how the adopted nonlinear processing is able to restore specific characteristics or properties of the original image. If the unknown image is well characterized using a probabilistic description, as for text images, the nonlinearity η(·) can be designed on the basis of a Bayesian criterion, as the "best" estimate of x[m, n] givenx (k) [m, n]. Often, the minimum mean square error (MMSE) criterion is adopted. For natural images, we design the nonlinearity η(·) after having represented the linear estimate 1x [m, n] in a transformed domain in which both the blur effect and the original image structural characteristics are easily understood.
We consider the decomposition of the linear estimatê x[m, n] by means of a filter pair composed of the lowpass filter ψ (0) [m, n] and a bandpass filter ψ (1) [m, n] (see Figure 3) whose impulse responses are Nonlinearity η(·) Locally tuned edge thinningx discrete polar pixel coordinates. These filters belong to the class of the circular harmonic functions (CHFs), whose detailed analysis can be found in [21,22], and possess the interesting characteristic of being invertible by a suitable filter For the values of the form factors σ 0 and σ 1 of interest, the corresponding transfer functions can be well approximated as follows: 2 2 , and γ(ω 1 , ω 2 ) def = arctan ω 2 /ω 1 , being the polar coordinates in the spatial radian frequency domain.
By indicating with (·) the complex conjugate operator, in the experiments, we have chosen to prevent amplification of spurious components occurring at those spatial frequencies, where Ψ (0) (e jω1 , e jω2 ) and Ψ (1) (e jω1 , e jω2 ) are small in magnitude. The optimality of these reconstruction filters is discussed in [23]. The zero-order circular harmonic filter ψ (0) [m, n] extracts a lowpass versionx 0 [m, n] of the input image; the form factor σ 0 is chosen so to retain only very low spatial frequencies, so obtaining a lowpass component exhibiting high spatial correlation. The first-order circular harmonic filter ψ (1) [m, n] is a bandpass filter, with frequency selectivity set by properly choosing the form factor σ 1 . The output of this filter is a complex imagex 1 [m, n], which will be referred to in the following as "edge image," whose magnitude is related to the presence of edges and whose phase is proportional to their orientation.
Coarsely speaking, the edge imagex 1 [m, n] is composed of curves, representing edges occurring in x[m, n], whose width is controlled by the form factor σ 1 , and of low magnitude values representing the interior of uniform or textured regions occurring in x[m, n]. Strong intensity curves inx 1 [m, n] are well analyzed by the local application of the bidimensional RT. This transform maps a straight line into a point in the transformed domain, and therefore it yields a compact and meaningful representation of the image's edges. However, since most image's edges are curves, the analysis must be performed locally by partitioning the image into regions small enough such that in each block, only straight lines may occur. Specifically, after having chosen the region dimensions, the value of the filter parameter σ 1 is set such that the width of the observed curve is a small fraction of its length. In more detail, the evaluation of the edge image is performed by the CH filter of order one ψ (1) [m, n] that can be seen as the cascade of a derivative filter followed by a Gaussian smoothing filter. The response to an abrupt edge of the original image is a line inx 1 [m, n]. The line is centered in correspondence to the edge, whose energy is concentrated in an interval of ±σ −1 1 pixels and that slowly decays to zero in an interval of ±3σ −1 1 pixels. Therefore, by partitioning the image into blocks of 8 × 8 pixels, the choice of σ 1 ≈ 1 yields edge structures that are well readable in the partitions of the edge image.
Then each region is classified as either a "strong edge" region or as a "weak edge" and "textured" region. The proposed enhancement procedures for the different kinds of regions are described in detail in Section 4.2.
It is worth pointing out that our approach shares the local RT as common tool with a family of recently proposed image transforms-the curvelet transforms [16,24,25]that represent a significant alternative to wavelet representation of natural images. In fact, the curvelet transform yields a sparse representation of both smooth image and edges, either straight or curved.

Local Radon transform of the edge image: original image and blur characterization
The edge image is a sparse complex image built by a background of zero or low magnitude areas, in which the objects appearing in the original image domain are sketched by their edges.
We discuss here this representation in more detail. For a continuous image ξ(t 1 , t 2 ), the RT [15] is defined as that represents the summation of ξ(t 1 , t 2 ) along a ray at distance s and angle β.
It is well known [15] that it can be inverted by where Some details about the discrete implementation of the RT follows.
If the image ξ(t 1 , t 2 ) is frequency limited in a circle of diameter D Ω , it can be reconstructed by the samples of its RT taken at spatial sampling interval ∆s ≤ 2π/D Ω , The angular interval ∆β can be chosen so as to assure that the distance between points p ξ β [n] and p ξ β+∆β [n] lying on adjacent diameters remains less than or equal to the chosen spatial sampling interval ∆s, that is, The above condition is satisfied when M ≥ (π/2)·N 1.57· N.
As long as the edge image is concerned, each region is here modeled as obtained by ideal sampling of an original image of the RT p x1 β (s) allow the reconstruction of the image x 1 (t 1 , t 2 ), and hence of any pixel of the selected region.  In Figure 4, some examples of straight edges and their corresponding discrete RT are shown.
We now consider the case of blurred observations. In the edge image, the blur tends to flatten and attenuate the edge peaks, and to smooth the edge contours in directions depending on the blur itself. The effects of blur on the RT of the edge image regions are primarily two. The first effect is that, since the energy of each edge is spread in the spatial domain, the maximum value of the RT is lowered. The second effect is that, since the edge width is typically thickened, it contributes to different tomographic projections, enhancing two triangular regions in the RT. This behavior is illustrated by the example in Figure 5, where a motion blur filter is applied to an original edge.
Stemming from this compact representation of the blur effect, we will devise an effective nonlinearity aimed at restoring the original edge.

Local Radon transform of the edge image: nonlinearity design
The design of the nonlinearity will be conducted after having characterized the blur effect at the output of a first-order CHF bank. By choosing the form factor σ 0 of the zero-order CH filter ψ (0) small enough, in the passband, the blur transfer function is approximately constant, and thus the blur effect on the lowpass component is negligible.
As far as the first-order CH filter's domain is concerned, the blur causes the edges in the spatial domain to be spread along directions depending on the impulse responses of the blurring filters. After having partitioned the edge image into small regions in order to perform a local RT as detailed in Section 4.1, each region has to be classified either as a strong edge area or a weak edge and textured area. Hence, the nonlinearity has to be adapted to the degree of "edgeness" of each region in which the image has been partitioned. The decision rule between the two areas is binary. Specifically, an area characterized as a "strong edge" region has an RT whose coefficients assume significant values only on a subset of directions β m . Therefore, a region is classified as a "strong edge" area by comparing max m n (p ξ βm [n]) 2 with a fixed threshold. If the threshold is exceeded, the area is classified as a strong edge area; otherwise, this implies that either no direction is significant, which corresponds to weak edges, or every direction is equally significant, which corresponds to textured areas.

Strong edges
For significant image edges, characterized by relevant energy concentrated in one direction, the nonlinearity can exploit the spatial memory related to the edge structure. In this case, as above discussed, we use the RT of the edge image. We consider a limited area of the edge imagex 1 [m, n] intersected by an edge, and its RT px 1 βm [m, n], with m, n chosen as discussed in Section 4.1.
The nonlinearity we present aims at focusing the RT both with respect to m and n, and it is given by where max n (px 1 βm [n]) and min n (px 1 βm [n]) represent the maximum and the minimum value, respectively, of the RT for the direction β m under analysis, and max βk,n (px 1 βk [n]) and min βk,n (px 1 βk [n]), with k = 0, . . . , M − 1, the global maximum and the global minimum, respectively, in the Radon domain. Therefore, for each point belonging to the direction β m and having index n, the nonlinearity (13) weights the RT by two gain functions. Specifically, (14) assumes its maximum value (equal to 1) for the direction β Max , where the global maximum occurs and it decreases for the other directions. In other words, (14) assigns a relative weight equal to 1 to the direction β Max whereas attenuates the other directions. Moreover, for a given direction β m , (15) determines the relative weight of the actual displacement n with respect to the others by assigning a weight equal to 1 to the displacement where the maximum occurs and by attenuating the other locations. The factors κ g and κ f in (14) and (15) are defined as κ g = κ 0 σ 2 w (k) and κ f = κ 1 σ 2 w (k) , σ 2 w (k) being the deconvolution noise variance and κ 0 and κ 1 two constants empirically chosen and set for our experiments equal to 2.5 and 0.5, respectively. The deconvolution noise variance σ 2 w (k) depends on both the blur and on the observation noise, and it can be estimated as E{|w (k) [m, n]| 2 } ≈ E{|x (k) [m, n] −x[m, n]| 2 } when the algorithm begins to converge. The deconvolution noise variance gradually decreases at each iteration which guarantees a gradually decreasing action of the nonlinearity as the iterations proceed.
The edge enhancement in the Radon domain is then described by the combined action of (14) and (15), since the first estimates the edge direction and the second performs a thinning operation for that direction.
To depict the effect of the nonlinearity (13) on the edge domain, in Figure 6, the case of a straight edge is illustrated.
The first columns of Figures 7 and 8 show some details extracted from the edge images of blurred versions of the "F16" (Figure 9) and "Cameraman" (Figure 10) images, respectively. For each highlighted block of 8 × 8 pixels, the RT is calculated by considering the block as belonging to a circular region of diameter 8 √ 2 (circumscribed circle). The above discussed nonlinearity is then applied. Then the inverse RT is evaluated for the pixels belonging to the central 8 × 8 pixels square. We observe that, although some pixels belong to two circles, namely, the circle related to the considered block    columns. The edges are clearly enhanced and focused by the processing.

Weak edges and textured regions
If the image is flat or does not exhibit any directional structure able to drive the nonlinearity, we use a spatially zeromemory nonlinearity, acting pointwise on the edge image. Since the edge image is almost zero in every pixel corresponding to the interior of uniform regions, where small values are likely due to noise, the nonlinearity should attenuate low-magnitude values ofx 1 [m, n]. On the other hand, high-magnitude values ofx 1 [m, n], possibly due the presence of structures, should be enhanced. A pointwise nonlinearity performing the said operations is given in the following: The magnitude of (16) is plotted in Figure 11 for different values of the parameter α. The low-gain zone near the origin is controlled by the parameter γ; the parameter α controls the enhancement effect on the edges. Both parameters are set empirically. The nonlinearity (16) has been presented in [14], where the analogy of this nonlinearity with the Bayesian estimator of spiky images in Gaussian observation noise is discussed.
To sum up, the adopted nonlinearity is locally tuned to the image characteristics. When the presence of an edge is detected, an edge thinning in the local RT of the edge image is performed. This operation, which encompasses a spatial   memory in the edge enhancement, is performed directly in the RT domain since the image edges are compactly represented in this domain. When an edge structure is not detected, which may happen for example in textured or uniform regions, the adopted nonlinearity reduces to a pointwise edge enhancement. It is worth pointing out that, as diffusely discussed in [16], the compact representation of an edge in the RT domain is related to the tuning between the size of the local RT transform and the bandpass of the edge extracting filter. After the nonlinear estimatex 1 [m, n] has been computed, the estimatex[m, n] is obtained by reconstructing through the inverse filter bank φ (0) [m, n] and φ (1) [m, n], that is, (see Figure 3) We remark that the nonlinear estimator modifies only the edge image magnitude, leaving the phase restoration to the linear estimation stage, performed by means of the filter bank The action of the nonlinearity on a natural image is presented in Figure 12, where along with the original image, the linear estimationx (1) [m, n], and the nonlinear estimationx (1) [m, n] obtained after the first iteration are shown.

EXPERIMENTAL RESULTS
In Figures 9, 10, and 13, some of the images we have used for our experimentations are reported. The images are blurred In Figure 14, some details belonging to the original image shown in Figure 9 are depicted. The corresponding blurred observations, affected by additive white Gaussian noise at SNR = 20 dB, obtained using the aforementioned blurring filters, are also shown along with the deblurred images. In Figure 15, the same images are reported for blurred images affected by additive white Gaussian noise at SNR = 40 dB. In Figure 16, the MSE, defined as is plotted versus the iterations number at different SNR values for the deblurred image. Similar results are reported in Figure 17 for the image shown in Figure 10 and the corresponding MSE is shown in Figure 18.
Moreover, in Figure 19, along with the restored version of the image depicted in Figure 13 obtained using the proposed method, the corresponding restored images obtained using the method introduced by the authors in [14] are reported. Eventually, in Figure 20, the MSE versus the number of iterations is plotted for both the proposed method and the one presented in [14]. Note that the algorithm converges in few iterations toward a local minimum and then may even diverge. In fact, as well known [17], the global convergence of the Bussgang algorithm cannot be guaranteed. The choice of the desired restoration result can be made by visual inspection of the restored images obtained during the iterations, which is, as also pointed out in [6], typical of blind restoration techniques.
It is worth noting that the deblurring algorithm gives images of improved visual quality, thus significantly reducing the distance, in the mean square sense, from the original unblurred image.

CONCLUSION
This work describes an algorithm for blind image restoration that iteratively performs a linear and a nonlinear processing. The nonlinear stage is based on a novel nonlinear processing technique that is based on a compact representation of the image edges by means of a local Radon transform. In order to focus the blurred image edges, and to drive the overall blind restoration algorithm to a sharp, focused image, a suitable nonlinear processing is designed in the Radon transformed domain. The nonlinear processing is spatially variant, depending on the detection of strong edges versus uniform or textured regions.
The performance of the algorithm is assessed by a number of experiments showing the overall quality of the restored images.

MULTICHANNEL BUSSGANG ALGORITHM
The multichannel Bussgang blind deconvolution algorithm outlined in [14] is here briefly summarized. Figure 19: Daughter image. Left column, first, second, and third row: observations. Fourth row: restored image using the method proposed in [14] (SNR = 20 dB). Fifth row: restored image using the proposed method (SNR = 20 dB). Right column, first, second, and third rows: observations. Fourth row: restored image using the method proposed in [14] (SNR = 40 dB). Fifth row: restored image using the proposed method (SNR = 40 dB). With reference to Figure 2, the kth iteration of the blind deconvolution algorithm is constituted by the following steps.
In general, the nonlinear estimator η(·) exhibits nonzero spatial memory.   Some of the authors have presented a discussion on the convergence of the algorithm in [26], where the analysis is conducted in an error energy reduction sense for the singlechannel case. The condition for convergence given in [26] straightforwardly extends to the multichannel case since the variables involved in the analysis (original image, deconvolved image, and nonlinear estimate) do not change.
Gaetano Scarano was born in Campobasso, Italy. He received the "Laurea" degree in electronic engineering from Università degli Studi di Roma "La Sapienza," Italy, in 1982. In 1982, he joined the Istituto di Acustica, Consiglio Nazionale delle Ricerche, Rome, Italy, as Associate Researcher. Since 1988, he has been teaching digital signal processing at Università degli Studi di Perugia, Perugia, Italy, where in 1991 he became Associate Professor of signal theory. In 1992, he joined the Dipartimento di Scienza e Tecnica dell'Informazione e della Comunicazione (Infocom), Università degli Studi di Roma "La Sapienza," first as an Associate Professor of image processing, then as a Professor of signal theory. His research interests lie in the areas of communications, signal and image processing, and estimation and detection theory, and include channel equalization and estimation, texture synthesis, and image restoration. He has served as an Associate Editor of the IEEE Signal Processing Letters.