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