The International Conference on Image Processing, Paris, France, October 27-30 2014Directly applying circular de-convolution to real-world blurred images usually results in boundary artifacts. Classic boundary extension techniques fail to provide likely results, in terms of a circular boundary-condition observation model. Boundary reflection gives raise to non-smooth features, especially when oblique oriented features encounter the image boundaries. Tapering the boundaries of the image support, or similar strategies (like constrained diffusion), provides smoothness on the toroidal support; however this does not guarantee consistency with the spectral properties of the blur (in particular, to its zeros). Here we propose a simple, yet effective, model-derived method for extending real-world blurred images, so that they become likely in terms of a Gaussian circular boundary-condition observation model. We achieve artifact-free results, even under highly unfavorable conditions, when other methods fail.Peer Reviewe

Portilla, Javier

Digital.CSIC

MAXIMUM LIKELIHOOD EXTENSION FOR NON-CIRCULANT DECONVOLUTIONJavier PortillaInstituto de O´pticaConsejo Superior de Investigaciones Cientı´ficas (CSIC), MadridABSTRACTDirectly applying circular de-convolution to real-world blurredimages usually results in boundary artifacts. Classic boundaryextension techniques fail to provide likely results, in terms ofa circular boundary-condition observation model. Boundaryreflection gives raise to non-smooth features, especially whenoblique oriented features encounter the image boundaries.Tapering the boundaries of the image support, or similarstrategies (like constrained diffusion), provides smoothnesson the toroidal support; however this does not guarantee con-sistency with the spectral properties of the blur (in particular,to its zeros). Here we propose a simple, yet effective, model-derived method for extending real-world blurred images,so that they become likely in terms of a Gaussian circularboundary-condition observation model. We achieve artifact-free results, even under highly unfavorable conditions, whenother methods fail.Index Terms— non-circulant deconvolution, boundaryartifacts, image restoration, maximum likelihood extension1. INTRODUCTIONFilters used in image restoration are typically expressed asa quotient in the discrete Fourier domain, and they are usu-ally applied in that domain. This applies both to linear andnon-linear (possibly state-of-the-art) techniques, as these alsotypically apply iterative linear filtering, combined with somenon-linear operations. By filtering through the DFT (DiscreteFourier Transform) we establish circular boundary-conditionsto the observation model. In real-world degraded images thisunnatural constraint, of course, does not hold. Mismatchingof circular and real boundary condition typically plagues therestoration result with boundary artifacts.In recent years, a growing interest has emerged in theboundary condition problem in image restoration. Reeves [1]considered the specific case of deconvolution being doneusing a regularized inverse linear operator with its classicalstructure (as the ones used in Richardson-Lucy, Tikhonov,Wiener filtering, etc.) to separate the (non-circulant) lineardeconvolution problem as a sequence of two steps: first, toThanks to Spanish Government funding grant TIN2012-38102-C03-01.compute a set of boundary pixels extending the observa-tion, fulfilling some constraints (a non-translation-invariantproblem, but with relatively few unknowns); and, second, toperform the filtering of the extended image, assuming circularboundary conditions. Recently Sorel [2] adapted that methodto an efficient non-linear restoration technique [3]. For iter-ative non-linear restoration methods the modification of theoriginal circulant deconvolution matrices also needs to bedone iteratively over partial results. Other authors (e.g., [4])have followed this same path of modifying previous restora-tion algorithms assuming circulant convolutions to deal withthe boundary condition problem in effective, but ad-hoc (i.e.,specific to the restoration method) ways. Here we follow thealternative approach of obtaining a single appropriated exten-sion of the observed image, independently of the restorationalgorithms (assumed DFT-based) to be applied.Liu and Jia [5] investigated the extension, through mir-ror reflection, of the image to another, four times bigger insurface, providing continuity on a toroidal support (a tilingimage, in computer graphics terminology). Then they solvedfor a variational problem in order to further smooth bound-ary transitions. Note that classical partial mirror extensiondoes not preserve continuity on a toroidal support. In gen-eral, whenever reflection is considered, strength of artifactswill critically depend on (1) the existence of oblique fea-tures near the image boundaries (as they generate new spatialfrequencies after being reflected); and (2) the lack of verti-cal/horizontal symmetry of the blurring kernel, which pro-vokes that the mirrored extensions can no longer be describedas filtered by the same kernel.Another simple universal strategy to avoid boundary arti-facts consists of using a spatially variant circulant smoothingoperator acting only on the vicinity of the image support’sboundaries, thus imposing smoothness on a toroidal support.This is implemented by the popular rMatlab function ed-getaper.m. However, methods like this or the referred [5] dis-regard the fact that smoothness on the toroidal support, onits own, does not fully comply with the circulant observationmodel. It is also necessary that the extended image respectsthe spectral valleys of the blurring kernel. This latter obser-vation is a central motivation of this work.978-1-4799-5751-4/14/$31.00 ©2014 IEEE ICIP 201442762. METHODS2.1. Realistic vs. circulant observation modelsA useful observation model for noisy convolved images isy = Hx+w, (1)where the H matrix represents an M2 ×N2 (with M < N ),non-circulant convolution, and w is the added noise term.Whereas using a non-circulant matrix is the correct choicefor real-world image convolution, most restoration algorithmassume, for convenience, a circular boundary condition ob-servation model:z = Hcx′ +w′, (2)where Hc is a square (N ′)2 × (N ′)2 block-circulant convo-lution matrix. We force the mean and autocovariance of theextended images to be identical to those of the original im-age (the latter computed/estimated with non-circular bound-ary conditions). A way of making compatible Eqs. (1) and (2)is to force Sz = y, where S is a M2 × (N ′)2 matrix, withN ′ ≥ N , selecting the observed M2 pixels of y from its ex-tended version. Hc is a square circulant extension of H. Therestoration problem of non-circulant deconvolution has beenthus split into three sequential steps: (1) find a suitable exten-sion zˆ(y) of y; (2) perform a circulant deconvolution of zˆ(y);(3) discard the added pixels.2.2. Problem assessmentThere are (N ′)2 −M2 degrees of freedom for extending yto z. Here we follow the criterion of choosing the exten-sion which maximizes the likelihood of an adapted, Gaussiancircular boundary-condition observation model. We have as-sumed that both the added noise w (considered white) andthe original image x vector terms follow an homogenous(translation-invariant) distribution. Then the Discrete FourierTransform diagonalizes the circulant covariance matrices oftheir extensions. In their diagonals appear the correspondingPower Spectral Densities (PSDs):PY (fi) = E|Z(fi)|2= |H(fi)|2E|X ′(fi)|2 + E|W ′(fi)|2= |H(fi)|2PX(fi) + σ2W , (3)for each of the (N ′)2 discrete frequencies fi, where X ′(f)represents here the Fourier transform of x′ (analogous with zand w′). Considering, for simplicity, Gaussian distributionsfor all the implied vectors, the observation’s minus log likeli-hood is (ignoring constant terms) ∝ ∑N ′i=1 |Z(fi)|2/PY (fi),and our maximum likelihood (ML) extension problem is:zˆ(y) = argminz||D−1/2PY F∗z||2, s.t. Sz = y, (4)where F∗ is a matrix performing a 2D-DFT, and DPY is adiagonal matrix whose entries correspond to PY (u, v), bothlexico-graphically ordered.2.3. Estimating PYThere are two alternative paths for estimating PY . First one isto compute it as the DFT of the sample auto-covariance of theobservation y. This possibility is appealing, as it allows to at-tack Eq.(4) without knowing the blurring kernel h or the noisevariance σ2w. However, estimating sample statistics poses atrade-off between variance and bias (noisy vs. blurred esti-mation) whose optimal balance point may be difficult to knowin advance. Second alternative consists of estimating/settingPX first, and then applying Eq.(3). The advantage of the lat-ter approach is that, even if our PX estimate is rough (e.g., ageneric signal model), degradation will nevertheless be reli-ably accounted for (provided we know h and σ2w with someaccuracy). In particular, it will suffice to identify the set of fre-quencies of the extended observation which are more heavilypenalized by the cost function (low PY (fi)). Here we havefollowed this second strategy.We have chosen a Gaussian Markov Random Field modelfor the original image x, a separable AR-1 (auto-regressive,one tap) model, whose PSD is:PX(u, v) =4σ2x log(ρ)2(log(ρ)2 + 4pi2u2)(log(ρ)2 + 4pi2v2). (5)We have hand-optimized its parameters, obtaining σx = 30and ρ = 0.65, for the images and experiments described inSection 3. We note that the behavior of the method was robustover change in these parameters, over a wide range.2.4. Unconstrained reformulationFirst, we split zP = Pz = [zTe zTi ]T into two non-overlappingvectors: zi, corresponding to theM2 original observed pixels,and ze, corresponding to the (N ′)2 −M2 extended pixels. Pis a permutation matrix. Analogously, the Fourier transformcan be re-ordered such that FP∗ = F∗P−1 = [F∗e F∗i ], corre-sponding to the two component of the Fourier transform ap-plied only to the external and the internal pixels, respectively.F∗z = FP∗zP = F∗eze + F∗i zi, and the extension constraintimplies zi = y. Now Eq.(4) simplifies tozˆ(y) = argminze||D−1/2PY (Fe∗ze + Fi∗y)||2, (6)whose solution is:zˆ(y) = −(FeD−1PY F∗e)−1FeD−1PY F∗iy. (7)In terms of practical implementation, it is useful to re-expressF∗e = F∗ET , and F∗i = F∗ST, whereET ((N ′)2× ((N ′)2−M2))) is the extension of external pixels to a vector contain-ing the whole image, the center being filled with zeros, ST(((N ′)2×M2) performs analogously, but preserving the cen-tral region. Their transpose do the reverse (selection) opera-tions. Thus, the difficulty in Eq. (7) lies only in the involvedmatrix inversion.978-1-4799-5751-4/14/$31.00 ©2014 IEEE ICIP 201442772.5. Recursive computationLet us call A = FeD−1PY F∗e and b(y) = FeD−1PYF∗iy. Firstnote thatA always exists when σ2w > 0. In this case existenceof A implies A−1 also exists. Thus, we can use Neumann’sseries [6] A−1 =∑∞n=0 (I−A)n to construct the followingrecursive method for computing zˆ(y) = A−1b(y):z(n+1) = z(n) + λ(b(y)−Az(n)),with, in principle, λ = 1. However, we accelerate conver-gence by using λ = 2/||A)|| (||A)|| is the largest eigenvalueof A in absolute value), without affecting the fixed-point so-lution. We have used a close approximation: λ = 2σ2w. As astarting point we use a mirror-extended image, edge-taperedwith a Gaussian kernel (see details in next section). The ob-servation’s sample mean was subtracted, and added back afterthe extension.3. EXPERIMENTS AND RESULTSWe have tested our method with three 8-bit gray-level 2562pixel images: two typical images (Cameraman and House),and an image from Brodatz photo texture album, (Straw),having a lot of energy along diagonal and oblique orienta-tions. We have used 8 degradations, with 4 blurring ker-nels, each adding two levels (low and medium) of noise.We have used the following point spread functions (PSFs),all of them with normalized sum: PSF1 is hi,j = (1 +i2 + j2)−1, for i, j = −7...7. PSF2 is a 9 × 9 uniformkernel. PSF3 simulates a uniform camera vertical move-ment 9 × 1. PSF4 is the oblique blurring kernel given by[0000111; 0012321; 0134310; 1232100; 1110000].We obtained y by cropping the M ×M valid central re-gion after performing circular convolution and noise addi-tion to the original. Note that M = N − 2L, being 2L +1 the width/height of the kernel support. We have appliedour extension to two very different restoration methods us-ing circular boundary conditions: the classic Wiener filter(by using the same spectral model of Eqs. (5) and (3)), andthe L2-relaxed L0 method ConDy, optimized for 10 itera-tions [7]. We have compared our results, measured as incre-ment in signal-to-noise ratio (ISNR, in decibels), to the resultsof: (1) Mirror extension, using Le = 8 (N ′ = N + 2Le =M + 2L + 2Le), being Le an extra extension of the origi-nal support (which we have tested it helps to reduce the ar-tifacts); (2) Full mirror extension, without cropping the mir-rored copies ((N ′ = 2M ); (3) Same boundary reflection as(1), but now followed by edge tapering using rMatlab ed-getaper.m, with a Gaussian kernel of σ = 6 (hand-optimizedvalue); (4) Boundary artifact-free oracle: obtained by restor-ing the whole N × N circulant convolved image, and thencropping the M ×M central region.To understand the influence of the degradation on the re-sults we have plotted the average mean square improvementratio (expressed as ISNR) in Figure 1, for every simulateddegradation. The oracle result is depicted as a white box indi-cating a practical upper bound performance. First, we observethat mirror extension-based methods are not robust. Lim-ited reflection (labeled as MIRROR) only performs accept-ably with PSF1, which is the only amongst the 4 used kernelsthat has no zeros in Fourier. Full reflection extension (labeledMIRROR-T, ”T” for tiling) is far more robust, but also failsfor PSF4. Lacking this kernel vertical/horizontal symmetry,the reflected observation is highly unlikely in terms of the ob-servation model. In contrast, edgetaper rarely provokes verystrong artifacts. However, its results are highly dependent onthe kernel and observed image, and do not get close, on aver-age, to the oracle performance, especially when the PSF haszeros in Fourier and the noise level is low. In such cases, theonly method amongst the compared ones which behaves closeto the reference oracle is our maximum likelihood extension.Table 1 shows, for reference purposes, the complete re-sults of our MLE method, as compared to the boundaryartifact-free oracle and the edgetaper rMatlab function.There we can see, among other things, how the oblique spa-tial frequencies, dominant in Straw, affect very negatively tothe edgetaper method, as implemented here (using mirror ex-tension prior to the edgetapering). The advantage of MLE isespecially prominent for the degradation experiments havinglow noise and PSFs having zeros in Fourier, and for those nonHV-symmetric (degradation 3, 5, 7 and 8). Figure 3 shows avisual comparison of the edgetaper extended result vs. MLEfor Straw, using Wiener restoration.Finally, Figure 2 shows how the quality gap between theedgetaper technique and the artifact-free oracle performanceis effectively saved by the MLE method. We see how 10 to50 iterations of the algorithm suffice to obtain good results inpractice. ML extension using 25 iterations in our implementa-tion (rMatlab over double quadrIntelrXeon architecture)takes around 0.5s for images this size.4. CONCLUSIONBoundary artifacts in image restoration are due to the mis-match between non-circulant convolved real-life observa-tions and the circulant observation model usually assumed byrestoration methods. To overcome this mismatch, we extendthe observation to make it consistent with a Gaussian circu-lant observation model, while keeping its auto-covariance.By maximizing likelihood on this model, we have obtained,in all studied cases, restoration results very close to those,artifact-free, obtained using circulant convolutions for sim-ulating the observation. We have tested the method againsta variety of PSFs, noise levels, original images, and evenagainst two different restoration algorithms.978-1-4799-5751-4/14/$31.00 ©2014 IEEE ICIP 201442781 2 3 4 5 6 7 80123456789#degradationaverage performance (ISNR, in dB)  MIRRORMIRROR?TEDGETAPERMLEFig. 1. Average performance, for the compared extensionmethods, and the 8 degradation (blur + noise) experiments.Restor→ Wiener ConDy10Method→ E-T MLE Ora E-T MLE OraPSF σ2w HOUSE1 .25 7.41 7.63 7.82 8.55 8.86 9.151 2 5.01 5.12 5.16 7.51 7.80 7.902 .31 4.95 6.56 6.69 9.14 10.25 10.342 4 4.09 4.11 4.07 7.08 7.29 7.423 1 5.18 6.20 6.27 7.68 9.14 9.273 4 4.10 4.35 4.34 6.38 7.24 7.404 .25 2.73 3.10 3.01 4.93 5.59 5.614 4 1.97 2.11 2.09 4.13 4.55 4.56PSF σ2w CAMERAMAN1 .25 8.67 8.70 8.71 10.56 10.58 10.621 2 5.56 5.55 5.56 8.06 8.11 8.122 .31 5.73 5.95 6.24 8.80 9.12 9.252 4 3.89 3.87 3.90 5.40 5.48 5.503 1 5.76 6.40 6.60 8.51 9.66 9.993 4 4.29 4.42 4.49 6.87 7.32 7.434 .25 3.77 3.90 3.90 6.58 6.85 6.914 4 2.32 2.35 2.36 3.66 3.76 3.78PSF σ2w STRAW1 .25 8.97 9.26 9.39 7.63 7.72 7.821 2 5.85 5.97 6.01 5.33 5.41 5.462 .31 3.06 4.56 4.84 3.37 4.29 4.492 4 2.31 2.53 2.60 1.92 2.08 2.133 1 5.73 7.93 8.46 5.90 8.03 8.563 4 5.00 5.82 6.03 4.67 5.53 5.734 .25 2.37 4.91 4.96 -0.67 3.17 3.194 4 2.62 3.15 3.19 1.79 2.19 2.20Table 1. Performance comparison, measured in Increment ofSignal-to Noise Ratio, in decibels. ”E-T” refers to ”edgeta-per”, ”Ora” to the oracle or upper bound reference used, and”MLE” to our ML extension method.100 101 102 10333.544.555.566.5#iterationsaverage performance (ISNR, in dB)  MLEORACLEEDGETAPERMIRROR?TFig. 2. Average performance, for ConDy10 restoration, as afunction of the number of iterations.Fig. 3. Upper left corner of Straw. From left to right: blurred(degradation 7), E-T + Wiener, MLE + Wiener.5. REFERENCES[1] S.J. Reeves, “Fast image restoration without boundary artifacts,”Image Processing, IEEE Transactions on, vol. 14, no. 10, pp.1448–1453, 2005.[2] M. Sorel, “Removing boundary artifacts for real-time iteratedshrinkage deconvolution,” Image Processing, IEEE Transac-tions on, vol. 21, no. 4, pp. 2329–2334, 2012.[3] J. Portilla, “Image restoration though `0 analysis-based sparseoptimization in tight frames,” in Image Processing (ICIP), 2009IEEE International Conference on, 2009, pp. 3909–3912.[4] A. Matakos, S. Ramani, and J.A. Fessler, “Image restoration us-ing non-circulant shift-invariant system models,” in Image Pro-cessing (ICIP), 2012 IEEE International Conference on, 2012,pp. 3061–3064.[5] Renting Liu and Jiaya Jia, “Reducing boundary artifacts in im-age deconvolution,” in Image Processing, (ICIP), 2008 IEEEInternational Conference on, 2008, pp. 505–508.[6] G.W. Stewart, Matrix Algorithms: Volume 1: Basic Decompo-sitions, Society for Industrial and Applied Mathematics, 1998.[7] J. Portilla, E. Gil-Rodrigo, D. Miraut, and R. Suarez-Mesa,“ConDy: Ultra-fast high performance restoration using multi-frame L2-relaxed-L0 sparsity and constrained dynamic heuris-tics,” in Image Processing, (ICIP), 2011 IEEE InternationalConference on, 2011, pp. 1837–1840.978-1-4799-5751-4/14/$31.00 ©2014 IEEE ICIP 20144279

Maximum likelihood extension for non-circulant deconvolution

Abstract

Similar works

Full text

Available Versions

Digital.CSIC