function samplingMask = createSamplingMask(im, rate, type)
%CREATESAMPLINGMASK Funkcja tworzy maskę próbkującą
%
%   im   - obraz z którego za pomoca trnasformaty Fouriera jest wyznaczana
%          funckja gęstości prawdopodobieństwa, dzięki której wyznaczona
%          zostanie adaptacyjna maska próbkujaca, wymiar podanego obrazu
%          jest także wymiarem tworzonej maski
%
%   rate - stosunek pikseli obrazu wejściowego do niezerowych pikseli maski
%
%   type - typ maski wedlug legendy:
%       'a' - losowa o rozkładzie jednostajnym
%       'b' - poziome linie z rozkładu normalnego
%       'c' - pionowe linie z rozkładu normalnego
%       'd' - krzyżujące się linie
%       'e' - koło z rozkładu jednostajnego
%       'f' - koło z rozkładu jednostajnego z kołem w środku
%       'g' - koło z rozkładu normalnego
%       'h' - adaptacyjna - wyznaczona losowo na podstawie transformaty
%             Fouriera


image=imadjust(im2double(im));
%tworzenie maski -  wybrane punkty otrzymuja wartość 1, niewybrane pozostają
%zerami
[M, N]=size(im);
samplingMask=zeros(size(image));
R=numel(samplingMask)/rate;   %liczba probek
counter=0;                    %wskaznik liczby wybranych probek



switch type
    case 'a'
        R=(M*N)/rate;
        counter=0;
        
        while counter<R
            m=randi(M);
            n=randi(N);
            if samplingMask(m,n)==0
                samplingMask(m,n)=1;
                counter=counter+1;
            end
        end
        
    case 'b'
        mu = M/2;
        sigma = 0.17*M;
        pd = makedist('Normal','mu',mu,'sigma',sigma);
        
        x = 1:M;
        y = cdf(pd,x); %dystrybuanta rozkładu normalnego
        
        R = M/rate;
        counter=0;
        while counter<R
            %wybranie losowej wspolrzednej z rozkladu normalnego
            difference=abs(y-rand);
            minM=min(difference);
            Mindex=find(difference==minM);
            
            if samplingMask(Mindex, 1)==0
                samplingMask(Mindex, :)=1;
                counter=counter+1;
            end
        end
        
    case 'c'
        mu = N/2;
        sigma = 0.17*N;
        pd = makedist('Normal','mu',mu,'sigma',sigma);
        
        x = 1:N;
        y = cdf(pd,x); %dystrybuanta rozkładu normalnego
        
        R = N/rate;
        counter=0;
        while counter<R
            %wybranie losowej wspołrzędnej z rozkładu normalnego
            difference=abs(y-rand);
            minN=min(difference);
            Nindex=find(difference==minN);
            
            if samplingMask(1,Nindex)==0
                samplingMask(:,Nindex)=1;
                counter=counter+1;
            end
        end
        
    case 'd'
        mu = [M N]/2;
        sigma = [M N]*0.1;
        pdM = makedist('Normal','mu',mu(1),'sigma',sigma(1));
        pdN = makedist('Normal','mu',mu(2),'sigma',sigma(2));
        
        xM = 1:M;
        yM = cdf(pdM,xM); %dystrybuanta rozkładu normalnego (wiersze)
        
        xN = 1:N;
        yN = cdf(pdN,xN); %dystrybuanta rozkładu normalnego (kolumny)
        
        R = [M N]/(rate*2);
        counterM=0;
        counterN=0;
        while counterM<R(1) || counterN<R(2)
            %wybranie losowej wspolrzednej z rozkladu normalnego (wiersze)
            difference=abs(yM-rand);
            minM=min(difference);
            Mindex=find(difference==minM);
            
            if samplingMask(Mindex, 1)==0
                samplingMask(Mindex, :)=1;
                counterM=counterM+1;
            end
            
            %wybranie losowej wspolrzednej z rozkladu normalnego (kolumny)
            difference=abs(yN-rand);
            minN=min(difference);
            Nindex=find(difference==minN);
            
            if samplingMask(1,Nindex)==0
                samplingMask(:,Nindex)=1;
                counterN=counterN+1;
            end
        end
        
    case 'e'
        S=[M, N]*0.5;
        samplingMask(S(1),S(2)) = 1;
        radius = 0.4*min(M,N);
        
        R=(M*N)/rate;
        counter=0;
        
        while counter<R
            m=randi(M);
            n=randi(N);
            if samplingMask(m,n)==0 && sqrt((S(1)-m)^2+(S(2)-n)^2)<radius
                samplingMask(m,n)=1;
                counter=counter+1;
            end
        end
        
        
    case 'f'
        S=[M, N]*0.5;
        samplingMask(S(1),S(2)) = 1;
        radius1 = 0.4*min(M,N);
        radius2=sqrt(R/pi)/2;
        samplingMask = bwdist(samplingMask) <= radius2;
        
        R=(M*N)/rate;
        counter=sum(sum(samplingMask));
        
        while counter<R
            m=randi(M);
            n=randi(N);
            if samplingMask(m,n)==0 && sqrt((S(1)-m)^2+(S(2)-n)^2)<radius1
                samplingMask(m,n)=1;
                counter=counter+1;
            end
        end
    case 'g'
        S=[M, N]*0.5;
        samplingMask(S(1),S(2)) = 1;
        samplingMask = bwdist(samplingMask) <= 3;
        counter=sum(sum(samplingMask));
        mux = M/2;
        sigmax = 0.15*M;
        pdx = makedist('Normal','mu',mux,'sigma',sigmax);
        
        x = 1:M;
        distx = cdf(pdx,x); %dystrybuanta rozkładu normalnego na osi x
        
        muy = N/2;
        sigmay = 0.15*N;
        pdy = makedist('Normal','mu',muy,'sigma',sigmay);
        
        y=1:N;
        disty = cdf(pdy, y); %dystrybuanta rozkładu normalnego na osi y
        
        while counter<R
            %wybranie losowej wspołrzędnej z rozkładu normalnego na osi x
            differencex=abs(distx-rand);
            minM=min(differencex);
            Mindex=find(differencex==minM);
            
            %wybranie losowej wspołrzędnej z rozkładu normalnego na osi y
            differencey=abs(disty-rand);
            minN=min(differencey);
            Nindex=find(differencey==minN);
            
            if samplingMask(Mindex, Nindex)==0
                samplingMask(Mindex, Nindex)=1;
                counter=counter+1;
            end
        end
        
    case 'h'
        % moduł z transformaty Fouriera wyznacza funkcję gęstości prawdopodobieństwa
        imageFT = abs(fftshift(fft2(image)));
        sumFT=sum(sum(imageFT));
        imagePDF=imageFT./sumFT; 
        % kolumny sa normowane tak, aby suma elementów w każdej kolumnie była równa 1
        imageFTnorm=normalize(imageFT,'norm',1);
        
        % wyznaczany jest rozkład brzegowy zmiennej X a następnie dystrybuanta zmiennej X
        margX=sum(imagePDF,1);
        cumX=cumsum(margX);
        
        while sum(sum(samplingMask))<R
            
            %     wybranie losowo pierwszej współrzędnej punktu maski
            %
            %     losowana jest liczba z przedziału (0,1) ktora odpowiada wartości
            %     dystrybuanty zmiennej X
            %
            %     wartość funkcji odwrotnej w tym punkcie wyznacza pierwszą
            %     współrzędną maski
                          
            difX=abs(cumX-rand);
            minX=min(difX);
            Xind=find(difX==minX);
            
            %     wybranie drugiej współrzędnej  punktu
            %
            %     dla wybranej pierwszej współrzędnej wybierany jest rozkład gęstości
            %     prawdopodbieństwa zmiennej Y wzdłuż wybranej współrzędnej, następnie
            %     druga współrzędna jest wybierana analogicznie jak wyżej
            cumY=cumsum(imageFTnorm(:,Xind));
            difY=abs(cumY-rand);
            minY=min(difY);
            Yind=find(difY==minY);
            if(samplingMask(Yind,Xind)==0)
                samplingMask(Yind,Xind)=1;
                counter=counter+1;
            end
        end
        
end
end