function P=psfkernel(ny,nx,g,par1,par2);

% build a sparse matrix K for the psf-kernel
% g   - string which contains the blur method
%       'gauss' : out of focus blur simulated by Gaussian
%               par1  -  sigma of Gaussian 
%               par2  -  r for cutoff of square of sidelength 2*r [r=3*sigma]
%       'line': spread along line, needs two more parameters
%               par1  - length of spread
%               par2  - unit vector for the direction of the spread

h=1/max([nx,ny]);
[xx,yy]=meshgrid(linspace(0,(nx-1)*h,nx),linspace(0,(ny-1)*h,ny));
switch g
    case 'gauss'
        if exist('par1')
            sigma=par1;
        else
            sigma=h;
        end
        if exist('par2')
            r=par2;
        else
            r=3*sigma;
        end
        na=floor(r/h);
        [xx,yy]=meshgrid(h*[-na:na]);
        Pcell=h^2/(2*pi*sigma^2)*exp(-(xx.^2+yy.^2)/(2*sigma^2));
        Pcell=Pcell/(sum(sum(Pcell)));  % renormalization due to cutoff
        nx1=-na; % first index of k in x-direction
        nx2=na;  % last  index of k in x-direction
        ny1=-na; % first index of k in y-direction
        ny2=na;  % last  index of k in y-direction
    case 'line'
        if exist('par1')
            L=par1;
        else
            L=h;
        end
        if exist('par2') % must then be a unit vector 
            v=par2;
            v=v/norm(v);
        else
            v=[1,0];  % horizontal direction of the blur
        end
        alpha=atan2(v(2),v(1));  % angle of the direction v
        na=floor(L/2*cos(alpha)/h);
        nb=floor(L/2*sin(alpha)/h);
        [xx,yy]=meshgrid(h*[-na:na],h*[-nb:nb]);
        Pcell=(abs(-v(2)*xx+v(1)*yy) <h/2).*(h/2-abs(-v(2)*xx+v(1)*yy)) ;
        Pcell=Pcell/(sum(sum(Pcell)));  % renormalization due to cutoff
        nx1=-na; % first index of k in x-direction
        nx2=na;  % last  index of k in x-direction
        ny1=-nb; % first index of k in y-direction
        ny2=nb;  % last  index of k in y-direction
    otherwise
        error('illegal string argument in psfkernel.m')
    end

num=0;
i=zeros(1,nx*ny*prod(size(Pcell)));
j=zeros(1,nx*ny*prod(size(Pcell)));
s=zeros(1,nx*ny*prod(size(Pcell)));
for k=1:nx
    kk=(k-1)*ny;
    nxlow=max(nx1,1-k);
    nxhigh=min(nx2,nx-k);
    tmpx=Pcell(:,(nxlow-nx1+1):(nxhigh-nx1+1));
    diffx=[nxlow:nxhigh];
    lenx=length(diffx);

    for l=1:ny
        kkl=kk+l; % row number of psf-kernel P
        nylow=max(ny1,1-l);
        nyhigh=min(ny2,ny-l);
        tmp=tmpx((nylow-ny1+1):(nyhigh-ny1+1),:);
        tmp=tmp(:);
        tmp=tmp/sum(tmp);  % renormalization near boundary
        diffy=[nylow:nyhigh];
        leny=length(diffy);        
        diffkkl=repmat(ny*diffx,leny,1)+repmat(diffy',1,lenx);
        diffkkl=diffkkl(:);
        len=lenx*leny;                    % number of entries of psf-kernel P in row kkl
        i(num+1:num+len)=kkl*ones(1,len); % put row number in next len entries
        j(num+1:num+len)=kkl+diffkkl';    % put column number in next len entries
        s(num+1:num+len)=tmp';            % store the values of next len entries
        num=num+len;
    end
end
numend=num;
i=i(1:numend);
j=j(1:numend);
s=s(1:numend);
P=sparse(i,j,s,nx*ny,nx*ny,numend);