"Introduction to Computer and Human Vision Exercise Set 2"
Introduction to Computer and Human Vision Exercise Set 2 Deblurring using Natural Image Priors Submission: individual or in pairs Due date: Thursday, December 31st , 2009 In this assignment we will investigate the solution of inverse problems using natural image priors. In particular we will focus on the deblurring problem and explore solutions using Gaussian and sparse priors. Brief overview: In class we discussed the problem of deblurring an m × n image y. We assume that the blurred image y is generated by y =x k+n where k is a blur kernel and n ∼ N (0, η 2 ) is zero mean Gaussian noise. Introducing a prior on natural images 1 y x |gij (x)|α +|gij (x)|α p(x) ∝ e− σ2 i,j y x where σ 2 is a parameter of the prior and gij (x) and gij (x) are the responses at coordinate (i, j) of the image to the ﬁlters (−1, 1) and (−1, 1)T respectively. We then formulated the problem of deblurring using a prior on the distribution of the edges in the image as ﬁnding a minimizer of 1 1 y L(x) = 2 y − Ax 2 + 22 x |gij (x)|α + |gij (x)|α 2η 2σ 1≤i≤m,1≤j≤n Where A is the mn × mn blur matrix corresponding to k, y is the column-stacked blurred image and η 2 is the variance of the noise in the blurred image. Setting α = 2 we get the Gaussian prior and using α = 0.8 we get the sparse prior. The “iteratively reweighted least squares” algorithm is an iterative procedure to ﬁnd a minimizer of a function of the form n E(x) = ρi (|bi x − zi |) i=1 Given an estimate of a minimizer xt deﬁne an iteration by n xt+1 = arg min q it (|bi x − zi |) i=1 1 Where q it (c) is a quadratic (scaler) function, which coincides with ρi at ct = bi xt − zi and which is strictly greater than ρi at all other points. From this we can derive three conditions on q it 1. q it (ct ) = ρi (ct ) 2. q it (ct ) = ρi (ct ) 3. q it is symmetric around the y-axis (thus q it (c) = ac2 + d for some constants a, d). These conditions determine q it . Since all we care about is ﬁnding the minimizer of q it we can neglect the additive term d and use ρ (ct ) 2 q it (c) = c 2ct Instructions: This is a Matlab programming assignment. • It is very important that you thoroughly document your code. In particular each function should begin with a short description of what it does and what are the input and output. • Remember that Matlab is eﬃcient when the code is vectorized, so try to avoid using loops. Speciﬁcally linear systems should be formulated in matrix notation and solved without using loops (looping over IRLS iterations is ﬁne). • You will most probably want to use the following Matlab functions: help, imread, imwrite, imshow, ﬁgure, conv2, reshape, spdiags... • In this exercise you will have to use Matlab’s sparse matrix representation. If you are not familiar with sparse matrices please consult the Matlab premier (you can ﬁnd a link on the “links” page of the course website). The important thing to note here is that if you don’t keep all your matrices sparse you will run out of memory. • You can use the function getConvMat.m which is posted on the website to generate the convolution matrices. • Save each matlab function in a separate ﬁle. 1. Write a function function x_star = solveGaussianPrior(A,y,eta,sigma) which computes x_star the minimizer of 1 2 1 y L(x) = y − Ax 2 + x |gij (x)|2 + |gij (x)|2 2η 2 2σ 2 1≤i≤m,1≤i≤n Where A is an mn × mn matrix (which could correspond to deblurring super resolution or any linear operator), y is an m × n image and eta and sigma are scalars as deﬁned above. The output is of size m × n. The optimization should be done in the spatial domain, solving a linear system of equations using Matlab’s backslash. 2. Write a function 2 function x_hat = solveSparsePrior(A,y,eta,sigma,x0,n_iteration) which computes x_hat an (approximate) minimizer of 1 2 1 y L(x) = y − Ax 2 + x |gij (x)|0.8 + |gij (x)|0.8 2η 2 2σ 2 1≤i≤m,1≤j≤n by running n_iteration iterations of “iteratively reweighted least squares” with an initial guess of x0 (also of size m × n). Thus the iteration should be 1 1t 1 y xt+1 = arg min 2 x q (|yij − aij x|) + 2 q 2t (|gij (x)|) + q 2t (|gij (x)|) 2η 2σ 1≤i≤m,1≤j≤n where q 1t (c) = c2 , q 2t (c) = |ct |−1.2 c2 and aij is the ij’th row of A. Or equivalently 1 2 1 y y xt+1 = arg min y − Ax 2 + x x wij |gij (x)|2 + wij |gij (x)|2 2η 2 2σ 2 1≤i≤m,1≤j≤n y y x x where wij = |gij (xt )|−1.2 and wij = |gij (xt )|−1.2 . To avoid division by zero we modify the weights to x x wij = max(|gij (xt )|, )−1.2 y y wij = max(|gij (xt )|, )−1.2 In this exercise we use = 1e − 5 and for the initial estimate x0 we use the solution of the Gaussian prior minimization. Technical note: assume that the input to the function is σ 2 = 2obsvar (which is the correct estimate for the Gaussian case). Use this value to compute x0 . Then to compensate for the inaccuracy of this estimation for the non- Gaussian case multiply σ by 5 and proceed with the IRLS iterations. In all the questions in this exercise we will perform 3 iterations of IRLS for the sparse deblurring. 3. Download the ﬁle ex2_q3.mat, this ﬁle contains three blur kernels ker_1,ker_2 and ker_3. It also contains 4 versions of the same image. im_1 the original image and y1 y2 y3 image im_1 blurred with the corresponding kernels plus noise. Use the clean image to estimate σ 2 . A good estimate is 2obsvar where obsvar is the mean of the squared y x ﬁlter responses (gij (xt ))2 and (gij (xt ))2 . You can also use the clean image to estimate the noise η 2 . (a) Estimate the noise for each one of the noisy images and the prior parameter σ 2 . Include the estimated values in your report along with a brief explanation of how you obtained the noise estimators. (b) Deblur the blurry images using a Gaussian image prior and then a sparse image prior. For each blur kernel j = 1, 2, 3 generate two images im_1_deblur_ker_j_Gauss and im_1_deblur_ker_j_sparse where j is and index which should be substituted with the value it assumes. Also generate a (3 × 3) image array of the following form. Row 1 consists of the three blurred images. Row 2 shows the results of the Gaussian deblurring and row 3 the results of sparse deblurring. Save this image array as a .png image called ex2_q3_b.png. 3 (c) Investigate the eﬀect of under and over estimating the noise parameter on debluring. Use y2 to generate a (1 × 3) image array of the following form. The left most image shows an example of deblurring when the noise is underestimated (multiply η by 0.2). The central image is deblurring with the correct noise estimation and the rightmost should show a result of deblurring with overestimation (multiply η by 10). Save this image array as a .png image called ex2_q3_c.png. Also include a brief discussion which explains your results and why you think they turn out as they do. (d) Investigate the eﬀect of using diﬀerent blur kernels on debluring. Generate a (2 × 2) image array of the following form. The left most image shows deblurring of y1 with ker_1. The next image shows the results of deblurring y1 with ker_3. The next row is deblurring of y3 with ker_1 and y3 with ker_3 Save this image array as a .png image called ex2_q3_d.png. Also include a brief discussion which explains your results and why you think they turn out as they do. 4. Download the ﬁle ex2_q4.mat, this ﬁle contains a blur kernel and a blurry-noisy image called ker and y respectively. Deblur using a Gaussian prior and a sparse prior. Use η = 0.01 and σ = 0.09. Save the results as im_2_deblur_sparse and im_2_deblur_Gauss Also save a (1×3) image array (blurry, Gaussian, sparse) as a .png image called ex2_q4.png. In this exercise we perform the debluring in the spatial domain, thus the computations are pretty demanding. In particular you will need more memory than is usually available on a home PC. Thus in this question you will have to perform the deblurring in overlapping blocks which should be stitched together to form the complete image. There should be an overlap so as to avoid boundary problems within the image (of course you will still get boundary eﬀects on the boundary of the stitched image). Another thing to note is that it is important to give the function solveSparsePrior an initial guess which comes from the ﬁnal stitched image of the Gaussian deblurring part. Submission instructions: Submit one .zip archive called CVFall09_ex2_Name1_FamilyName1_Name2_FamilyName2.zip by sending an email to firstname.lastname@example.org with the subject “CVfall2009 ex2”. This zip archive should contain 1. A discussion of the results in pdf, ps or rtf format. The ﬁrst page should include the names of the submitters and your ID numbers. 2. A .mat ﬁle called ex2.mat which includes the deblurred images from question 3 (im_1_deblur_ker_j_Gauss and im_1_deblur_ker_j_sparse for j = 1, 2, 3) and from question 4 (im_2_deblur_sparse and im_2_deblur_Gauss). 3. all the .png images namely: ex2_q3_b.png ex2_q3_c.png ex2_q3_d.png ex2_q4.png. 4. Two .m ﬁles implementing the functions from questions 1 and 2. and the scripts / functions which you used to answer questions 3 and 4. 5. Please submit a hard copy of your code and report to the “computer vision” mailbox. 4