function out=blur(m,n,keep)
%
% BLUR(m,n,keep) creates an (m*n) x (m*n) matrix
% that represents -- FOR ILLUSTRATION ONLY -- the
% linear map that replaces each value of an m x n
% matrix by the sum of KEEP times its old value 
% and (1-KEEP)/4 times the values of each neighbor.
%
% Applications -- demo only -- might include the
% blurring/sharpening of images or heat eqn/Laplace 
% eqn. 
% Note that the matrix blur(m,n) is a HUGE, very 
% sparse matrix (at most 5 nonzero entries in each 
% row/col). Thus this implementation is more to
% illustrate THEOETICAL relations between vectors/
% linear transfos and their coordinates/matrix rep's.
%
% Author: Matthias Kawski. July 2000
%         http://math.la.asu.edu/~kawski
%
if nargin < 3  % of no 3rd parameter is provided, the
   keep =0.60; % default is 60% old + 10% each of 
   end         % the neighboring values
smear=(1-keep)/4;  
%
B=keep*eye(m*n); % initialize the diagonal (corners
                     % and edges will be corrected later
%
% interior neighbors
%
for i=1:n-2
   for j=2:m-1
      B(i*m+j,i*m+j-1)=smear;
      B(i*m+j,i*m+j+1)=smear;
      B(i*m+j,(i-1)*m+j)=smear;
      B(i*m+j,(i+1)*m+j)=smear;
   end
end
%
% correct the edges
%
for j=2:m-1
   B(j,j)=keep+smear;
   B(j,m+j)=smear;
   B(j,j-1)=smear; 
   B(j,j+1)=smear;
   B((n-1)*m+j,(n-1)*m+j)=keep+smear;
   B((n-1)*m+j,(n-2)*m+j)=smear;
   B((n-1)*m+j,(n-1)*m+j-1)=smear; 
   B((n-1)*m+j,(n-1)*m+j+1)=smear;
end
for i=1:n-2
   B(i*m+1,i*m+1)=keep+smear;
   B(i*m+1,i*m+2)=smear;
   B(i*m+1,(i-1)*m+1)=smear;
   B(i*m+1,(i+1)*m+1)=smear;
   B((i+1)*m,(i+1)*m)=keep+smear;
   B((i+1)*m,(i+1)*m-1)=smear;
   B((i+1)*m,i*m)=smear;
   B((i+1)*m,(i+2)*m)=smear;
end
%
% correct the corners
%
B(1,1)=keep+2*smear;
B(m,m)=B(1,1);
B((n-1)*m+1,(n-1)*m+1)=B(1,1);
B(m*n,m*n)=B(1,1);
%
B(1,2)=smear;
B(1,m+1)=B(1,2);
B(m,m-1)=B(1,2);
B(m,2*m)=B(1,2);
B((n-1)*m+1,(n-2)*m+1)=B(1,2);
B((n-1)*m+1,(n-1)*m+2)=B(1,2);
B(n*m,n*m-1)=B(1,2);
B(n*m,(n-1)*m)=B(1,2);
%
out=B;