blur.m 2.31 KB
Newer Older
 Jan-Gerrit Richter committed Jul 29, 2016 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 ``````function [A,b,x] = blur(N,band,sigma) %BLUR Test problem: digital image deblurring. % % function [A,b,x] = blur(N,band,sigma) % % The matrix A is an N*N-by-N*N symmetric, doubly block Toeplitz matrix that % models blurring of an N-by-N image by a Gaussian point spread function. % It is stored in sparse matrix format. % % In each Toeplitz block, only matrix elements within a distance band-1 % from the diagonal are nonzero (i.e., band is the half-bandwidth). % If band is not specified, band = 3 is used. % % The parameter sigma controls the width of the Gaussian point spread % function and thus the amount of smoothing (the larger the sigma, the wider % the function and the more ill posed the problem). If sigma is not % specified, sigma = 0.7 is used. % % The vector x is a columnwise stacked version of a simple test image, while % b holds a columnwise stacked version of the blurrred image; i.e, b = A*x. % Per Christian Hansen, IMM, 11/11/97. % Initialization. if (nargin < 2), band = 3; end band = min(band,N); if (nargin < 3), sigma = 0.7; end % Construct the matrix as a Kronecker product. z = [exp(-((0:band-1).^2)/(2*sigma^2)),zeros(1,N-band)]; A = toeplitz(z); A = sparse(A); A = (1/(2*pi*sigma^2))*kron(A,A); % Generate x and b, if required. if (nargout > 1) % Start with an image of all zeros. x = zeros(N,N); N2 = round(N/2); N3= round(N/3); N6 = round(N/6); N12 = round(N/12); % Add a large ellipse. T = zeros(N6,N3); for i=1:N6 for j=1:N3 if ( (i/N6)^2 + (j/N3)^2 < 1 ), T(i,j) = 1; end end end T = [fliplr(T),T]; T = [flipud(T);T]; x(2+(1:2*N6),N3-1+(1:2*N3)) = T; % Add a smaller ellipse. T = zeros(N6,N3); for i=1:N6 for j=1:N3 if ( (i/N6)^2 + (j/N3)^2 < 0.6 ), T(i,j) = 1; end end end T = [fliplr(T),T]; T = [flipud(T);T]; x(N6+(1:2*N6),N3-1+(1:2*N3)) = x(N6+(1:2*N6),N3-1+(1:2*N3)) + 2*T; % Correct for overlap. f = find(x==3); x(f) = 2*ones(size(f)); % Add a triangle. T = triu(ones(N3,N3)); [mT,nT] = size(T); x(N3+N12+(1:nT),1+(1:mT)) = 3*T; % Add a cross. T = zeros(2*N6+1,2*N6+1); [mT,nT] = size(T); T(N6+1,1:nT) = ones(1,nT); T(1:mT,N6+1) = ones(mT,1); x(N2+N12+(1:mT),N2+(1:nT)) = 4*T; % Make sure x is N-times-N, and stack the columns of x. x = reshape(x(1:N,1:N),N^2,1); % Compute the blurred image. b = A*x; end``````