function M=mask(siz,sigma,X,Y,type,sigma2);
% siz can be even or odd
% X,Y denotes the subpixel distance of the centre of the
% gaussian from the centre of the mask
% M=mask(siz,sigma,X,Y,type,sigma2);

M=zeros(siz);

[x,y] = meshgrid(-(siz-1)/2:(siz-1)/2,-(siz-1)/2:(siz-1)/2);
x=x+X;
y=y+Y;

if type=='gau',
    h = (1/(2*pi*sigma2^2))*exp(-(x.*x + y.*y)/(2*sigma2*sigma2));
    %h = h/sum(sum(h)); % can keep this line if want to be accurate
    
elseif type=='dog',
    h1 = (1/(2*pi*sigma^2))*exp(-(x.*x + y.*y)/(2*sigma*sigma));
    %h1 = h1/sum(sum(h1)); % can keep this line if want to be accurate
    h2 = (1/(2*pi*sigma2^2))*exp(-(x.*x + y.*y)/(2*sigma2*sigma2));
    %h2 = h2/sum(sum(h2)); % can keep this line if want to be accurate
   
    h=h1-h2;
    
elseif type=='cir',
    h1 = exp(-(x.*x + y.*y)/(2*sigma*sigma));
    h1 = h1/sum(sum(h1));
    h2 = exp(-(x.*x + y.*y)/(2*sigma2*sigma2));
    h2 = h2/sum(sum(h2));
    
    h=h1-h2;
    h=-(h.*(h>0));
    
elseif type=='lap',
    std = sigma*sigma;
    h = (x.*x + y.*y - 2*std)./((pi^2)*std);
    h = h - sum(h(:))/prod(size(h));   
    
elseif type=='gab',
    %[C,S] = gabormask(size,sigma,[],sigma2);
    %h=C;
        
    orient = sigma2;
    period = [];
    Size = siz;
    
    sy=Size-1; sx=Size-1;
    period = sigma*2*sqrt(2); 
    
    % Basic grid
    hy = sy/2; hx = sx/2;
    [x, y] = meshgrid(-hx:sx-hx, -hy:sy-hy);
    x=x+X;
    y=y+Y;
    
    % Parameters
    omega = 2*pi/period;
    cs = omega * cos(orient);
    sn = omega * sin(orient);
    k = -1/(2*sigma*sigma);
    
    % Main computations
    g = exp(k * (x.*x + y.*y));     % Gaussian mask
    xp = x * cs + y * sn;           % Rotated x coords, phase
    cx = cos(xp);                   % cos grating
    cmask = g .* cx;                % modulated cos grating
    sx = sin(xp);                   % sin grating
    smask = g .* sx;                % modulated sin grating
    
    % Normalise so that convolution of mask with a harmonic curve of the
    % matching frequency gives unity peaks
    cnorm = sum(sum(cmask.*cx));
    cmask = cmask/cnorm;
    snorm = sum(sum(smask.*sx));
    smask = smask/snorm;
    
    C=cmask; S=smask;
    h{1}=C;
    h{2}=S;
    
    
end

%  [sy,sx]=size(h);
%  element=h(round(sy/2),1);
%  M=h.*(h>=element);
% 
%  M=M/sum(sum(M));

M=h;