function h = function_GetPSF(type,par1,par2)

% Create special 2-D filters
%
% Syntax
% 
% h = function_GetPSF(type)
% h = function_GetPSF(type,par1)
% h = function_GetPSF(type,par1,par2)
%
% Description
% 
% h = function_GetPSF(type) creates a two-dimensional filter h of the 
% specified type. 
% 
% the following types are allowed
%  
% h = function_GetPSF('ggh_45') returns an averaging filter of size 5x5
% which is rotated by 45 degree.
% 
% h = function_GetPSF('inv_quad',par1) returns an inverse-quadratic filter
% h(x1,x2) = (1+x1^2+x2^2)^-1 with support par1 x par1
% 
% h = function_GetPSF('random',par1) returns a liner filter with support of
% size par1 x par1 where coefficients a generated randomly. The sum of
% filter coefficients is equal to 1.
% 
% h = function_GetPSF('test1',par1) returns a circular averaging filter 
% (pillbox) with with the raduis equal to par1 (see help fspecial). The
% sector from 0 till pi/2 is equal to 0.
%  
% other options replicate MATLAB function 'fspecial' of
% image processing toolbox. For more details see help fspecial

type = lower(type);
if strcmp(type,'ggh_45'),
    load ggh_45
    v = ggh_45; 
    h = v/sum(v(:));
    
elseif strcmp(type,'inv_quad'),
    s1=0;
    for a1=-par1:par1;
        s1=s1+1;
        s2=0;
        for a2=-par1:par1;
            s2=s2+1; v(s1,s2)=1/(a1^2+a2^2+1);
        end,
    end;
    h = v./sum(v(:));
    
elseif strcmp(type,'sroubek1'),
    a = [0 0 1 0 0; 
         0 1 1 1 0;
         1 1 1 1 1;
         0 1 1 1 0;
         0 0 1 0 0];
    h = a./sum(a(:));

elseif strcmp(type,'sroubek2'),
    a = [0 0 1 0 0; 
         0 0 1 0 0;
         0 0 1 1 1;
         1 1 1 1 1;
         1 1 1 1 1];
    h = a./sum(a(:));

elseif strcmp(type,'test1'),
    a = fspecial('disk',par1);
    [yN,xN] = size(a);
    yN2 = ceil(yN/2); xN2 = ceil(xN/2);
    a(1:yN2-1,xN2+1:end) = 0;
    h = a./sum(a(:));

elseif strcmp(type,'test2'),
    a = fspecial('average',par1);
    [yN,xN] = size(a);
    yN4 = ceil(yN/4); xN3 = ceil(xN/3);
    a(1:yN4,xN3:end) = 0;
    a((end-yN4+1):end,1:(end-xN3+1)) = 0;
    h = a./sum(a(:));
    
elseif strcmp(type,'random'),
    a = randn([par1 par1]);
    h = a./sum(a(:));
     
else
    if nargin<3,
        h = fspecial(type,par1);
    else
        h = fspecial(type,par1,par2);
    end;
end;