blur.m

求解离散病态问题的正则化方法matlab 工具箱

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