function A = GaussBlur(N, band, sigma)
%GAUSSBLUR: GIDE digital image deblurring.
%
% function A = GaussBlur(N, band, sigma)
%
% This function creates an N^2 x N^2 blurring matrix A 
% corresponding to a Gaussian point-spread function (PSF) 
% assuming zero boundary conditions.
%
% Input: 
%   N     - The image is N x N.
%   band  - 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.
%   sigma - The parameter sigma controls the width of 
%           the Gaussian point spread function and thus 
%           the amount of smoothing. Larger sigma means
%           broader blur and more ill-posedness.
%           If sigma is not specified, sigma = 0.7 is used.
%
% The matrix A is a symmetric, doubly block Toeplitz matrix 
% stored in sparse matrix format. It models blurring of an 
% N x N image by a Gaussian point spread function.
%
%
% Code based on that of Per Christian Hansen, IMM, 11/11/97,
% Brianna Cash and Dianne O'Leary 11/2014
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Defaults:

if (nargin < 2), 
   band = 3; 
end

if (nargin < 3), 
   sigma = 0.7; 
end

% Construct the matrix A as a Kronecker product.

band = min(band,N);

z = [exp(-((0:band-1).^2)/(2*sigma^2)),zeros(1,N-band)];

A = toeplitz(z);

A = sparse(A);

A = kron(A,A);

A = A / max(sum(A,2)); % Normalize

