function [obr_fil] = median3D(maska,P,obr)

%% **********************************************************************************
%*  Kamil Sokol		
%*	Mediánová filtrácia 3D s voliteľnou maskou
%*  maska   - vstupná 3D maska pre filtráciu
%*  P       - počet iterácii
%*  obr     - vstupné 3D obrazové dáta pre filtráciu
%*  obr_fil - výstupné 3D obrazové dáta - filtrované
%************************************************************************************

%% Funkcia **************************************************************************
tStart=tic;                 %štart merania času výpočtu 
obr = im2double(obr);
maska(maska==0) = NaN;      %zamena hodnot 0 v maske za Nan 
[m,n,o] = size (obr);       %zistenie rozmerov org. obrazu
[p,r,s] = size (maska);     %zistenie rozmerov masky
obr_fil = zeros(m,n,o);     %vytvorenie prázdnej matice pre výsledny obr.
okno = zeros(p,r,s);
obr_zvacseny = zeros((m+p-1),(n+r-1),(o+s-1));  %vytvorenie prázdnej matice pre zvacseny obr.

fprintf('Medianova filtrácia 3D: START\n');    %Vypisanie informacie
fprintf('Medianova filtrácia 3D hotova z:  00.00 %%');


m_k1 = round((p-1)/2);      %korekcia - o kolko sa ma zvecsit <-m
m_k2 = p - m_k1 - 1;        %korekcia - o kolko sa ma zvecsit m->

n_k1 = round((r-1)/2);      %korekcia - o kolko sa ma zvecsit <-n
n_k2 = r - n_k1 - 1;        %korekcia - o kolko sa ma zvecsit n->

o_k1 = round((s-1)/2);      %korekcia - o kolko sa ma zvecsit <-o
o_k2 = s - o_k1 - 1;        %korekcia - o kolko sa ma zvecsit o->

%************************************************************************************
%% Filtrácia ************************************************************************
for pp = 1:1:P              %opakovanie podla poctu iteracii
    
    obr_zvacseny((1+m_k1):(m+m_k1),(1+n_k1):(n+n_k1),(1+o_k1):(o+o_k1)) = obr(:,:,:);

    for i = 1:1:(m_k1)
        obr_zvacseny(i,:,:) = obr_zvacseny(1+m_k1,:,:);         %zvecsenie <-m
    end

    for i = 1:1:(m_k2)
        obr_zvacseny(m_k1+m+i,:,:) = obr_zvacseny(m_k1+m,:,:);  %zvecsenie m->
    end

    for i = 1:1:(n_k1)
        obr_zvacseny(:,i,:) = obr_zvacseny(:,1+n_k1,:);         %zvecsenie <-n
    end

    for i = 1:1:(n_k2)
        obr_zvacseny(:,n_k1+n+i,:) = obr_zvacseny(:,n_k1+n,:);  %zvecsenie n->
    end

    for i = 1:1:(o_k1)
        obr_zvacseny(:,:,i) = obr_zvacseny(:,:,1+o_k1);         %zvecsenie <-o
    end

    for i = 1:1:(o_k2)
        obr_zvacseny(:,:,o_k1+o+i) = obr_zvacseny(:,:,o_k1+o);  %zvecsenie o->
    end
    
    for i = 1:1:(m)                     %Pohyb okna v 3 osách
        for j = 1:1:(n)
            for k = 1:1:(o)
                okno = obr_zvacseny(i:(i+p-1),j:(j+r-1),k:(k+s-1)); %naplnenie okna hodnotami
                okno2 = okno.*maska;                    %aplikacia masky
                okno3 = (okno2(:))';                    %hodnoty do jednej rady
                obr_fil(i,j,k) = nanmedian (okno3);     %zapisanie medianu ako vysleny bod
            end
        end
        progres = ((pp*i)/(P*m))*100;   %výpočet progresu 
        fprintf('\b\b\b\b\b\b\b%05.2f %%',progres);
    end
    obr(:,:,:) = obr_fil(:,:,:);
end
tElapsed=toc(tStart);                 %koniec merania času výpočtu
fprintf('\nMedianova filtrácia 3D vypočítaná za: %.1f sekúnd \nMedianova filtrácia 3D: HOTOVO \n\n',tElapsed);
%************************************************************************************
