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