 function A = blur2d(n,z)

 if nargin<2, z=.01; end

 threshold = .1/(pi*z);
 
 K = @(s1,s2,t1,t2) exp(-((s1-t1).^2+(s2-t2).^2)/z)/(pi*z);

 s = ([1:n]-.5)/n;
 t = ([1:n]-.5)/n;

 A = zeros(n*n);

 for k=1:n, for ell=1:n
    svec = [s(k);s(ell)];
    row = (ell-1)*n + k;
    for j=1:n
       Kijkl = K(s(k),s(ell),t,t(j));
       Kijkl = Kijkl.*(Kijkl>=threshold);
       A(row,(j-1)*n+[1:n]) = Kijkl/n^2;
    end
 end, end


