function mat = myMap(type,N,sig,x,y)

% creates a predefined map

switch type
    case 'phantom'
        mat = phantom('Modified Shepp-Logan',N);
    case 'noise'
        mat = rand(N);
    case 'mickey'
        temp = myMapElem('disk',N,0.24,0,0.2)- ...
            myMapElem('disk',N,0.28,0,0);
        temp = temp.*double(temp>=0);
        mat = myMapElem('disk',N,0.2,-0.3,-0.2)+ ...
            myMapElem('disk',N,0.2,0.3,-0.2)+ ...
            myMapElem('disk',N,0.4,0,0.2)-temp;
        mat = mat-myMapElem('disk',N,0.02,-0.16,0.06)- ...
            myMapElem('disk',N,0.02,0.16,0.06);
    case 'nicedisks' % x plays #disks, y = distance from the center
        mat = myMapElem('zero',N,0.12,0.55,0.6);
        pts = [0., y*exp(2*pi*i*(1/2/x:1/x:1-1/2/x))];
        for j=1:x+1,
            mat = mat + myMapElem('disk',N,sig,real(pts(j)),imag(pts(j)));
        end
    case 'nicebumps' % x plays #bumps, y = distance from the center
        mat = zeros(N);
        pts = [0.,y*exp(2*pi*i*(1/2/x:1/x:1-1/2/x))];
        for j=1:x+1,
            mat = mat + myMapElem('gaussian',N,sig,real(pts(j)),imag(pts(j)));
        end
    case 'disks'
        mat = zeros(N);
        for j=1:size(x,2),
            mat = mat + myMapElem('disk',N,sig,x(j),y(j));
        end
    case 'gaussians'
        mat = zeros(N);
        for j=1:size(x,2),
            mat = mat + myMapElem('gaussian',N,sig,x(j),y(j));
        end
    case 'stripes'
        mat = myMapElem('line',N,sig,0,0);
        mat = mat + myMapElem('line',N,sig,0,0.5);
        mat = mat + myMapElem('line',N,sig,0,-0.5);
    case 'potatoes'
        mat = myMapElem('potatoe',N,sig,x(1),y(1));
        for j=2:size(x,2),
            mat = mat + myMapElem('potatoe',N,sig,x(j),y(j));
        end
    case 'spiral'
        mat = myMapElem('disk',N,sig,0.,0.);
        for j=1:8,
            mat = mat + myMapElem('disk',N,sig,j/10*cos(j*pi/6),j/10*sin(j*pi/6));
            mat = mat + myMapElem('disk',N,sig,-j/10*cos(j*pi/6),-j/10*sin(j*pi/6));
        end
        
    otherwise
        mat = myMapElem(type,N,sig,x,y);
end

end


%--------------------------------------------------------------------------
function mat = myMapElem(type,N,sig,x,y)

% creates a predefined map
% x,y: coordinates of a center in the square [0,1]x[0,1]
% N: size of the resulting image
% sig: radius, width, theta for the cross...

xs = -1+1/N:2/N:1.-1/N;
[X,Y] = meshgrid(xs,xs);

switch type
    case 'one'
        mat = ones(N);
    case 'gaussian'
        mat = exp(-((X-x).^2 +(Y-y).^2)/2/sig^2);
    case 'disk'
        mat = double(((X-x).^2+(Y-y).^2)<=sig^2);
    case 'square'
        mat = double(abs(X)<=0.5).*double(abs(Y)<=0.5);
    case 'sines'
        mat = 1. + sin(2*pi*X).*sin(2*pi*Y);
    case 'cross' % sig stands for theta in this case
        mat = double(abs(cos(sig)*X +sin(sig)*Y)<0.01)+ ...
            double(abs(-sin(sig)*X+cos(sig)*Y)<0.01);
        mat = mat.*double((X.^2+Y.^2)<0.04);
        mat = mat + double((X.^2+Y.^2)<1);
    case 'y'
        mat = Y;
    case 'line'
        mat = double(abs(X-x + Y-y)<=sig);
    case 'potatoe'
        % (x,y) play the center position
        % sig = [r0, r1, .., rn; 0, t1, .., tn]
        ncoeffs = size(sig,2)-1;
        Z = X-x + i*(Y-y);
        argZ = angle(temp)/2/pi;
        rs = sig(1,1)*ones(N);
        for j=1:ncoeffs,
            rs = rs + sig(1,j+1)*cos(2*pi*j*(argZ-sig(2,j+1)));
        end
        mat = double(abs(Z)<=rs);
    otherwise
        mat = zeros(N);
end

end