function[kmask,dmask]=segmentacia(path,varargin)
%Funkcia vykonava segmentaciu spinalneho kanalu a intervertebralnych diskov
%v T2-vahovanych objemovych MRI snimkach chrbtice. Data musia byt ulozene
%vo formate NIFTI, prvym vstupom funkcie je cesta k datovemu suboru (retazec
%znakov 'path'), volitelnym vstupnym parametrom je suradnica (kladne cele cislo) 
%oznacujuca polohu sagitalneho rezu, ktory ma byt pouzity pre vyber inicializacnych 
%semienok pre segmentaciu. Pokial tento parameter nie je uvedeny, automaticky sa
%zvoli prostredny sagitalny rez objemovych dat. V tomto reze sa semienka urcuju 
%interaktivne pocitacovou mysou (kliknutim sa zvoli 1 bod, potiahnutim a vyznacenim 
%linie viac bodov). Vystupom funkcie su binarne masky pre spinalny kanal ('kmask') 
%a IVD ('dmask') oznacujuce polohu vymedzenych objektov v povodnom obraze (1=objekt, 
%0=pozadie). Taktiez je zobrazeny vysledok segmentacie osobitne pre kanal a IVD. 

%Vytvorene v programovom prostredi MATLAB R2015b.

%kontrola vstupov

if ~ischar(path)
    error('Zadajte prosím vstup(y) v správnom formáte')
elseif length(varargin)>1
       error('Zadajte prosím vstup(y) v správnom formáte')
end

%parametre

beta=90;
P=0.6;
r=8;
Irw=5;
Ils=10;
smoothfactor=0.5;
contractionbias=-0.6;
Ier=5;
IKdil=10;
IDdil=5;

%nacitanie dat

A1=load_nii(path);
o=A1.img;
[y,x,z]=size(o);
oo=zeros(y,x,z);
for i=1:z
    oo(:,:,i)=o(:,:,z+1-i);
end
im=zeros(z,x,y);
for j=1:y
a=oo(j,:,:);
aa=zeros(z,x);
    for i=1:z
        aa(i,:)=a(1,:,i);
    end
im(:,:,j)=aa;
end
obr=oo;

%normalizacia

oo=(oo-min(oo(:)))/(max(oo(:))-min(oo(:)));
obr=(obr-min(obr(:)))/(max(obr(:))-min(obr(:)));
im=(im-min(im(:)))/(max(im(:))-min(im(:)));

%--------------------------------------------------------------------------
%------------------------------SPINALNY KANAL------------------------------
%--------------------------------------------------------------------------

%---------------------POCIATOCNE SEMIENKA A SEGMENTACIA--------------------

%pociatocne semienka

if isempty(varargin)
   V=round(y/2);
elseif ~isnumeric(varargin{1})
       error('Zadajte prosím vstup(y) v správnom formáte')
elseif varargin{1}<1
       error('Zadajte prosím vstup(y) v správnom formáte')
elseif varargin{1}-round(varargin{1})~=0
       error('Zadajte prosím vstup(y) v správnom formáte')
elseif varargin{1}>y
       error('Zvolená súradnica presahuje rozmery obrazu')
else V=varargin{1};
end
im1=im(:,:,V);
figure('Name','Semienka')
imshow(im1,[])
title('Zadajte počiatočný bod (body) pre segmentáciu')
h=imfreehand('Closed',0);
pos=round(getPosition(h));
if pos(1,2)>pos(end,2)
   pos=flipud(pos);
end
%osetrenie duplicity
[~,IA,~]=unique(pos(:,2),'stable');
pos=pos(IA,:);
xx=pos(:,1)';
yy=pos(:,2)';
if numel(pos)==2
   X=xx;
   Y=yy;
else N=pos(1,2):1:pos(end,2);
     pp=spline(yy,xx,N);
     pp=round(pp);
     X=pp;
     Y=N;
end

%random walker

fmask=zeros(y,x,z);
LX=length(X);
for i=1:LX
    img1=obr(:,:,Y(i));
    seed2x=zeros(50);
    seed2y=zeros(50);
    %semienka objektu
    seed1x=X(i);
    seed1y=V;
    seed1=sub2ind([y,x],seed1y,seed1x);
    %semienka pozadia
    
    %kruznica
    for j=0:2*pi/50:2*pi-2*pi/50
        k=round(j*50/(2*pi)+1);
        seed2x(k)=round(X(i)+r*cos(j));
        seed2y(k)=round(V+r*sin(j));
    end
    %celistva kruznica
    pom1=seed2y-1;
    pom2=seed2y+1;
    pom3=seed2x+1;
    pom4=seed2x-1;
    seed2x=[seed2x pom3   pom4   seed2x seed2x];
    seed2y=[seed2y seed2y seed2y pom1   pom2  ];
    %osetrenie rozsahu
    p1=(find(seed2x>x))';
    p2=(find(seed2y>y))';
    p3=(find(seed2x<1))';
    p4=(find(seed2y<1))';
    seed2x([p1,p2,p3,p4])=[];
    seed2y([p1,p2,p3,p4])=[];
    %doplnok kruznice
    seed2=sub2ind([y,x],seed2y,seed2x);
    pom=zeros(y,x);
    pom(seed2)=1;
    pom=imfill(pom,'holes');
    seed2=(find(pom==0))';

    %algoritmus rw
    [mask,~] = random_walker1(img1,[seed1,seed2],...
    [ones(1,length(seed1)),zeros(1,length(seed2))],beta,P);
    %vytvorenie masky
    fmask(:,:,Y(i))=mask;
end

%--------------------------ITERACIE RW + LEVELSET--------------------------

for q=1:Irw
    for i=1:LX
        
        %random walker
        
        img1=obr(:,:,Y(i));
        seed2x=zeros(50);
        seed2y=zeros(50);
        %semienka objektu
        seed1=(find(fmask(:,:,Y(i))))';
        %centroid
        [seed1y,seed1x]=ind2sub([y x],seed1);
        Yc=mean(seed1y);
        Xc=mean(seed1x);
        %semienka pozadia
        
        %kruznica
        for j=0:2*pi/50:2*pi-2*pi/50
            k=round(j*50/(2*pi)+1);
            seed2x(k)=round(Xc+r*cos(j));
            seed2y(k)=round(Yc+r*sin(j));
        end
        %celistva kruznica
        pom1=seed2y-1;
        pom2=seed2y+1;
        pom3=seed2x+1;
        pom4=seed2x-1;
        seed2x=[seed2x pom3   pom4   seed2x seed2x];
        seed2y=[seed2y seed2y seed2y pom1   pom2  ];
        %osetrenie rozsahu
        p1=(find(seed2x>x))';
        p2=(find(seed2y>y))';
        p3=(find(seed2x<1))';
        p4=(find(seed2y<1))';
        seed2x([p1,p2,p3,p4])=[];
        seed2y([p1,p2,p3,p4])=[];
        %doplnok kruznice
        seed2=sub2ind([y,x],seed2y,seed2x);
        pom=zeros(y,x);
        pom(seed2)=1;
        pom=imfill(pom,'holes');
        seed2=(find(pom==0))';

        %osetrenie duplicity (objekt vs pozadie)
        for m=1:length(seed1)
            for n=1:length(seed2)
                if seed1(m)==seed2(n)
                   seed1(m)=0;
                end
            end
        end
        seed1(seed1==0)=[];
        %algoritmus rw
        [mask,~] = random_walker1(img1,[seed1,seed2],...
        [ones(1,length(seed1)),zeros(1,length(seed2))],beta,P);
        %vytvorenie masky
        fmask(:,:,Y(i))=mask;
    end
end
fmask1=zeros(y,x,z);
for i=1:LX  
    
    %levelset
    
    img1=obr(:,:,Y(i));
    %inicializacna maska
    mask=fmask(:,:,Y(i));
    maski=(find(mask))';
    %centroid
    [masky,maskx]=ind2sub([y x],maski);
    Yc=mean(masky);
    Xc=mean(maskx);
    %algoritmus ls 
    mask1=activecontour(img1,mask,Ils,'Chan-Vese','SmoothFactor',smoothfactor,...
           'ContractionBias',contractionbias);
    mask1=double(mask1);
    %obmedzenie
         
    %median euklidovskych vzdialenosti pixelov povodnej masky od centroidu
    h=median(sqrt((maskx-Xc).^2+(masky-Yc).^2));
    %euklidovske vzdialenosti pixelov novej masky od centroidu
    mask1i=find(mask1);
    [mask1y,mask1x]=ind2sub([y x],mask1i);
    v=sqrt((mask1x-Xc).^2+(mask1y-Yc).^2);
    %vyradenie odlahlych pixelov
    mask1(mask1i(v>h+3))=0;

    %vytvorenie masky
    fmask1(:,:,Y(i))=mask1;
end

%------------------------PROPAGACIA NADOL (RW + LS)------------------------

mask=fmask(:,:,Y(LX));
for i=Y(LX)+1:z
    
    %random walker
    
    img1=obr(:,:,i);
    seed2x=zeros(50);
    seed2y=zeros(50);
    %semienka objektu
    [fx,fy]=gradient(mask);
    xc=find(fx);
    yc=find(fy);
    mask([xc',yc'])=0;
    [fx,fy]=gradient(mask);
    xc=find(fx);
    yc=find(fy);
    mask([xc',yc'])=0;
    seed1=(find(mask))';
    if sum(sum(mask))==0
        break
    end
    %centroid
    [seed1y,seed1x]=ind2sub([y x],seed1);
    Yc=mean(seed1y);
    Xc=mean(seed1x);
    %semienka pozadia
    
    %kruznica  
    for j=0:2*pi/50:2*pi-2*pi/50
        k=round(j*50/(2*pi)+1);
        seed2x(k)=round(Xc+r*cos(j));
        seed2y(k)=round(Yc+r*sin(j));
    end
    %celistva kruznica
    pom1=seed2y-1;
    pom2=seed2y+1;
    pom3=seed2x+1;
    pom4=seed2x-1;
    seed2x=[seed2x pom3   pom4   seed2x seed2x];
    seed2y=[seed2y seed2y seed2y pom1   pom2  ];
    %osetrenie rozsahu
    p1=(find(seed2x>x))';
    p2=(find(seed2y>y))';
    p3=(find(seed2x<1))';
    p4=(find(seed2y<1))';
    seed2x([p1,p2,p3,p4])=[];
    seed2y([p1,p2,p3,p4])=[];
    %doplnok kruznice
    seed2=sub2ind([y,x],seed2y,seed2x);
    pom=zeros(y,x);
    pom(seed2)=1;
    pom=imfill(pom,'holes');
    seed2=(find(pom==0))';

    %osetrenie duplicity (objekt vs pozadie)
    for m=1:length(seed1)
        for n=1:length(seed2)
            if seed1(m)==seed2(n)
               seed1(m)=0;
            end
        end
    end
    seed1(seed1==0)=[];
    %algoritmus rw
    [mask,~] = random_walker1(img1,[seed1,seed2],...
    [ones(1,length(seed1)),zeros(1,length(seed2))],beta,P);
    if sum(sum(mask))==0
        break
    end
    %obmedzenie
    
    %median euklidovskych vzdialenosti semienok od centroidu
    h=median(sqrt((seed1x-Xc).^2+(seed1y-Yc).^2));
    %euklidovske vzdialenosti pixelov novej masky od centroidu
    maski=find(mask);
    [masky,maskx]=ind2sub([y x],maski);
    v=sqrt((maskx-Xc).^2+(masky-Yc).^2);
    %vyradenie odlahlych pixelov
    mask(maski(v>h+3))=0;

    %levelset
    
    mask1=activecontour(img1,mask,Ils,'Chan-Vese','SmoothFactor',smoothfactor,...
        'ContractionBias',contractionbias);
    mask1=double(mask1);
    %obmedzenie
    
    %euklidovske vzdialenosti pixelov novej masky od centroidu
    maski=find(mask1);
    [masky,maskx]=ind2sub([y x],maski);
    v=sqrt((maskx-Xc).^2+(masky-Yc).^2);
    %vyradenie odlahlych pixelov
    mask1(maski(v>h+3))=0;
    
    %vytvorenie masky
    fmask1(:,:,i)=mask1;
end

%-----------------------PROPAGACIA NAHOR (RW + LS)-------------------------

mask=fmask(:,:,Y(1));
for i=Y(1)-1:-1:1
    
    %random walker
    
    img1=obr(:,:,i);
    seed2x=zeros(50);
    seed2y=zeros(50);
    %semienka objektu
    [fx,fy]=gradient(mask);
    xc=find(fx);
    yc=find(fy);
    mask([xc',yc'])=0;
    seed1=(find(mask))';
    if sum(sum(mask))==0
        break
    end
    %centroid
    [seed1y,seed1x]=ind2sub([y x],seed1);
    Yc=mean(seed1y);
    Xc=mean(seed1x);
    %semienka pozadia
    
    %kruznica  
    for j=0:2*pi/50:2*pi-2*pi/50
        k=round(j*50/(2*pi)+1);
        seed2x(k)=round(Xc+r*cos(j));
        seed2y(k)=round(Yc+r*sin(j));
    end
    %celistva kruznica
    pom1=seed2y-1;
    pom2=seed2y+1;
    pom3=seed2x+1;
    pom4=seed2x-1;
    seed2x=[seed2x pom3   pom4   seed2x    seed2x];
    seed2y=[seed2y seed2y seed2y pom1      pom2  ];
    %osetrenie rozsahu
    p1=(find(seed2x>x))';
    p2=(find(seed2y>y))';
    p3=(find(seed2x<1))';
    p4=(find(seed2y<1))';
    seed2x([p1,p2,p3,p4])=[];
    seed2y([p1,p2,p3,p4])=[];
    %doplnok kruznice
    seed2=sub2ind([y,x],seed2y,seed2x);
    pom=zeros(y,x);
    pom(seed2)=1;
    pom=imfill(pom,'holes');
    seed2=(find(pom==0))';

    %osetrenie duplicity (objekt vs pozadie)
    for m=1:length(seed1)
        for n=1:length(seed2)
            if seed1(m)==seed2(n)
               seed1(m)=0;
            end
        end
    end
    seed1(seed1==0)=[];
    %algoritmus rw
    [mask,~] = random_walker1(img1,[seed1,seed2],...
    [ones(1,length(seed1)),zeros(1,length(seed2))],beta,P);
    if sum(sum(mask))==0
        break
    end
    %obmedzenie
    
    %median euklidovskych vzdialenosti semienok od centroidu
    h=median(sqrt((seed1x-Xc).^2+(seed1y-Yc).^2));
    %euklidovske vzdialenosti pixelov novej masky od centroidu
    maski=find(mask);
    [masky,maskx]=ind2sub([y x],maski);
    v=sqrt((maskx-Xc).^2+(masky-Yc).^2);
    %vyradenie odlahlych pixelov
    mask(maski(v>h+3))=0;
    
    %levelset
    
    mask1=activecontour(img1,mask,Ils,'Chan-Vese','SmoothFactor',smoothfactor,...
        'ContractionBias',contractionbias);
    mask1=double(mask1);
    %obmedzenie
    
    %euklidovske vzdialenosti pixelov novej masky od centroidu
    maski=find(mask1);
    [masky,maskx]=ind2sub([y x],maski);
    v=sqrt((maskx-Xc).^2+(masky-Yc).^2);
    %vyradenie odlahlych pixelov
    mask1(maski(v>h+3))=0;
    
    %vytvorenie masky
    fmask1(:,:,i)=mask1;
end

%-------------------------VYHLADENIE (1. iteracia)-------------------------

[px,py,pz] = ndgrid(-3:3);
se=strel(sqrt(px.^2+py.^2+pz.^2)<=1);
kmask=fmask1;
dprah=mean(oo(fmask1==1))-std(oo(fmask1==1));
hprah=mean(oo(fmask1==1))+std(oo(fmask1==1));

%erozia

for i=1:Ier
    er=imerode(kmask,se);
    fmask1(er==1)=0;
    pos=find(fmask1==1);
    kor=ones(y,x,z);
    kor(pos)=oo(pos);
    kor(kor<dprah)=0;
    kmask(kor==0)=0;
end

%dilatacia

for i=1:IKdil
    dil=imdilate(kmask,se);
    dil(kmask==1)=0;
    pos=find(dil==1);
    kor=zeros(y,x,z);
    kor(pos)=oo(pos);
    kor(kor>hprah)=1;
    kmask(kor==1)=1;
    %obmedzenie + vyplnenie medzier
    
    for j=1:z
        mask=kmask(:,:,j);
        maski=find(mask);
        if isempty(maski)
            continue
        end
        [masky,maskx]=ind2sub([y x],maski);
        %centroid
        Yc=mean(masky);
        Xc=mean(maskx);
        %median euklidovskych vzdialenosti pixelov masky od centroidu
        v=sqrt((maskx-Xc).^2+(masky-Yc).^2);
        h=median(v);
        %vyradenie odlahlych pixelov
        mask(maski(v>h+3))=0;
        %vyplnenie medzier 
        mask=imfill(mask,'holes');
        for k=round(Yc)-10:1:round(Yc)+10
            if k<2 || k>y-1
               continue
            end
            for l=round(Xc)-10:1:round(Xc)+10
                if l<2 || l>x-1
                   continue
                end
                okolie=[mask(k-1,l-1),mask(k-1,l),mask(k-1,l+1),mask(k,l-1),...
                       mask(k,l+1),mask(k+1,l-1),mask(k+1,l),mask(k+1,l+1)];
                if sum(okolie)>3
                   mask(k,l)=1;
                end
            end
        end
        kmask(:,:,j)=mask;
    end
end

%-------------------------LEVELSET (2. iteracia)---------------------------

for i=1:z
    img1=obr(:,:,i);
    mask=kmask(:,:,i);
    maski=find(mask);
    [masky,maskx]=ind2sub([y x],maski);
    
    %centroid
    
    Yc=mean(masky);
    Xc=mean(maskx);
    
    %median euklidovskych vzdialenosti pixelov povodnej masky od centroidu
    
    h=median(sqrt((maskx-Xc).^2+(masky-Yc).^2));
    
    %inicializacna maska
    
    [fx,fy]=gradient(mask);
    xc=find(fx);
    yc=find(fy);
    mask([xc',yc'])=0;
    [fx,fy]=gradient(mask);
    xc=find(fx);
    yc=find(fy);
    mask([xc',yc'])=0;
    if sum(sum(mask))==0
        continue
    end

    %levelset
    
    mask=activecontour(img1,mask,Ils,'Chan-Vese','SmoothFactor',smoothfactor,...
          'ContractionBias',contractionbias);
    mask=double(mask);
    %obmedzenie
    
    %euklidovske vzdialenosti pixelov novej masky od centroidu
    maski=find(mask);
    [masky,maskx]=ind2sub([y x],maski);
    v=sqrt((maskx-Xc).^2+(masky-Yc).^2);
    %vyradenie odlahlych pixelov
    mask(maski(v>h+3))=0;
    
    %vytvorenie masky
    kmask(:,:,i)=mask;
end

%-------------------------VYHLADENIE (2. iteracia)-------------------------

fmask2=kmask;
dprah=mean(oo(kmask==1))-std(oo(kmask==1));
hprah=mean(oo(kmask==1))+std(oo(kmask==1));

%erozia

for i=1:Ier
    er=imerode(kmask,se);
    fmask2(er==1)=0;
    pos=find(fmask2==1);
    kor=ones(y,x,z);
    kor(pos)=oo(pos);
    kor(kor<dprah)=0;
    kmask(kor==0)=0;
end

%dilatacia

for i=1:IKdil
    dil=imdilate(kmask,se);
    dil(kmask==1)=0;
    pos=find(dil==1);
    kor=zeros(y,x,z);
    kor(pos)=oo(pos);
    kor(kor>hprah)=1;
    kmask(kor==1)=1;
    %obmedzenie + vyplnenie medzier
    
    for j=1:z
        mask=kmask(:,:,j);
        maski=find(mask);
        if isempty(maski)
            continue
        end
        [masky,maskx]=ind2sub([y x],maski);
        %centroid
        Yc=mean(masky);
        Xc=mean(maskx);
        %median euklidovskych vzdialenosti masky od centroidu
        v=sqrt((maskx-Xc).^2+(masky-Yc).^2);
        h=median(v);
        %vyradenie odlahlych pixelov
        mask(maski(v>h+3))=0;
        %vyplnenie medzier 
        mask=imfill(mask,'holes');
        for k=round(Yc)-10:1:round(Yc)+10
            if k<2 || k>y-1
               continue
            end
            for l=round(Xc)-10:1:round(Xc)+10
                if l<2 || l>x-1
                   continue
                end
                okolie=[mask(k-1,l-1),mask(k-1,l),mask(k-1,l+1),mask(k,l-1),...
                       mask(k,l+1),mask(k+1,l-1),mask(k+1,l),mask(k+1,l+1)];
                if sum(okolie)>3
                   mask(k,l)=1;
                end
            end
        end
        kmask(:,:,j)=mask;
    end
end
    
%zobrazenie 

%axialna rovina
va=oo;
va(kmask==1)=max(max(max(oo)));
figure('Name','Spinálny kanál (axiálna rovina)')
imshow3D(va)
%sagitalna rovina
vs=zeros(z,x,y);
for j=1:y
    a=va(j,:,:);
    aa=zeros(z,x);
    for i=1:z
        aa(i,:)=a(1,:,i);
    end
    vs(:,:,j)=aa;
end
figure('Name','Spinálny kanál (sagitálna rovina)')
imshow3D(vs)

%--------------------------------------------------------------------------
%--------------------------INTERVERTEBRALNE DISKY--------------------------
%--------------------------------------------------------------------------

%vymedzenie prednej linie spinalneho kanalu

mss=kmask(V,:,:);
linia=[];
for i=1:z
    pom=mss(1,:,i);
    poz=find(pom,1,'first');
    if isempty(poz)
    else linia=[linia;i poz];
    end
end
xx=linia(:,2)';
yy=linia(:,1)';
N=yy(1):yy(end);
pp=spline(yy,xx,N);
pp=round(pp);
liniax=pp;
liniay=N;
kk=zeros(z,x);
for i=1:length(liniay)
    kk(liniay(i),liniax(i)-20:liniax(i)-12)=im1(liniay(i),liniax(i)-20:liniax(i)-12);
end

%prahovanie linie

prah=multithresh(im1(kk~=0),2);
kk(kk<prah(1))=0;
kk(kk~=0)=1;
kk=imfill(kk,'holes');

%distancna mapa linie (jej doplnku)

D=bwdist(~kk);

%watershed

D=-D;
D(~kk)=-Inf;
L=watershed(D);

%odstranenie hranic oblasti

L(1,:)=1;
L(end,:)=1;
L(:,1)=1;
L(:,end)=1;
Q=L;
for i=1:z
    for j=1:x
        if L(i,j)==0
           okolie=[L(i-1,j-1),L(i-1,j),L(i-1,j+1),L(i,j-1),L(i,j+1),L(i+1,j-1),L(i+1,j),L(i+1,j+1)];
           okolie(okolie==1)=[];
           okolie(okolie==0)=[];
           Q(i,j)=mode(okolie); 
        end
    end
end
L=Q;

%zlucenie dotykajucich sa oblasti

for q=1:3
    for i=1:z
        for j=1:x
            if L(i,j)~=1
               okolie=[L(i-1,j),L(i,j-1)];
               okolie(okolie==1)=[];
               okolie(okolie==L(i,j))=[];
               if isempty(okolie)
               else L(i,j)=okolie(1); 
               end
            end
        end
    end
    for i=z:-1:1
        for j=x:-1:1
            if L(i,j)~=1
               okolie=[L(i,j+1),L(i+1,j)];
               okolie(okolie==1)=[];
               okolie(okolie==L(i,j))=[];
               if isempty(okolie)
               else L(i,j)=okolie(1); 
               end
            end
        end
    end
end

%uprava indexovania oblasti

L(L==1)=0;
ind=unique(L);
ind(1)=[];
M=ind(end);
pom=zeros(1,length(ind));
LL=L';
for i=1:length(ind)
    pom(i)=find(LL==ind(i),1);
end
[~,I]=sort(pom);
for i=1:length(ind)
    L(L==ind(I(i)))=i+M;
end
L=L-M;
L(L==-M)=0;
L=double(L);

%vyplnenie chybajucich IVD

%stanovenie vzdialenosti medzi oblastami
D=[];
for q=1:max(max(L))-1
    f=find(L==q);
    [I,J]=ind2sub([z x],f);
    m=max(I);
    if (I==m-1)==0
        pom=find(I==m);
    elseif (I==m-2)==0
            pom=find(I==m|I==m-1);
    else pom=find(I==m|I==m-1|I==m-2);
    end
    I=I(pom);
    J=J(pom);
    ff=find(L==q+1);
    [II,JJ]=ind2sub([z x],ff);
    d=[];
    for i=1:length(I)
        for ii=1:length(II)
            dd=sqrt((I(i)-II(ii))^2+(J(i)-JJ(ii))^2);
            d=[d dd];
        end
    end
    DD=min(d);
    D=[D DD];
end
%vyplnenie IVD
M=mode(round(D));
v=find(D>M+2);
for i=1:length(v)
    poz=v(i);
    p1=find(L==poz);
    [p1r,p1c]=ind2sub([z x],p1);
    r1=max(p1r);
    p11=[p1c(p1r==r1);p1c(p1r==r1-1);p1c(p1r==r1-2)];
    c1=min(p11);
    c2=max(p11);
    p2=find(L==poz+1);
    [p2r,~]=ind2sub([z x],p2);
    r2=min(p2r);
    if ~isempty(r1+2:r2-2)
    L(r1+2:r2-2,c1:c2)=poz+0.5;
    else v(i)=[];
    end
end
%uprava indexovania
for i=1:length(v)
    pom=find(L>=v(i)+1);
    L(pom)=L(pom)+1;
    v=v+1;
end
L=round(L);
 
%rozdelenie zlucenych oblasti (prve 3 objekty)

PP1=find(L==1);
PP2=find(L==2);
PP3=find(L==3);
[PP1r,~]=ind2sub([z x],PP1);
[PP2r,~]=ind2sub([z x],PP2);
[PP3r,~]=ind2sub([z x],PP3);
r1min=min(PP1r);
r1max=max(PP1r);
r2min=min(PP2r);
r2max=max(PP2r);
r3min=min(PP3r);
r3max=max(PP3r);
PPP1=zeros(1,length(r1min:r1max));
PPP2=zeros(1,length(r2min:r2max));
PPP3=zeros(1,length(r3min:r3max));
%medialna linia-1
        
%1.objekt
DIV1=[];
for i=r1min:r1max
    PPP1(i-r1min+1)=L(i,liniax(i-liniay(1)+1)-17);
end
i1=find(PPP1,1,'first');
i2=find(PPP1,1,'last');
PPP1=PPP1(i1:i2);
POZ=find(PPP1==0);
if ~isempty(POZ)
    dpoz=diff(POZ);
    if isempty(dpoz)
       DIV1=POZ+r1min-1;
    else pom=find(dpoz~=1);
         if isempty(pom)
            POZ=POZ(round(length(POZ)/2));
            DIV1=POZ+r1min-1;
        else pom=[0 pom length(POZ)];
             DIV1=zeros(1,length(pom)-1);
             for i=1:length(pom)-1
                 ppom=(pom(i)+1):pom(i+1);
                 pppom=round(length(ppom)/2);
                 DIV1(i)=POZ(ppom(pppom));
             end
             DIV1=DIV1+r1min-1+i1-1;
        end
    end
end
%2.objekt
DIV2=[];
for i=r2min:r2max
    PPP2(i-r2min+1)=L(i,liniax(i-liniay(1)+1)-17);
end
i1=find(PPP2,1,'first');
i2=find(PPP2,1,'last');
PPP2=PPP2(i1:i2);
POZ=find(PPP2==0);
if ~isempty(POZ)
    dpoz=diff(POZ);
    if isempty(dpoz)
       DIV2=POZ+r2min-1;
    else pom=find(dpoz~=1);
         if isempty(pom)
            POZ=POZ(round(length(POZ)/2));
            DIV2=POZ+r2min-1;
        else pom=[0 pom length(POZ)];
             DIV2=zeros(1,length(pom)-1);
             for i=1:length(pom)-1
                 ppom=(pom(i)+1):pom(i+1);
                 pppom=round(length(ppom)/2);
                 DIV2(i)=POZ(ppom(pppom));
             end
             DIV2=DIV2+r2min-1+i1-1;
        end
    end
end
%3.objekt
DIV3=[];
for i=r3min:r3max
    PPP3(i-r3min+1)=L(i,liniax(i-liniay(1)+1)-17);
end
i1=find(PPP3,1,'first');
i2=find(PPP3,1,'last');
PPP3=PPP3(i1:i2);
POZ=find(PPP3==0);
if ~isempty(POZ)
    dpoz=diff(POZ);
    if isempty(dpoz)
       DIV3=POZ+r3min-1;
    else pom=find(dpoz~=1);
         if isempty(pom)
            POZ=POZ(round(length(POZ)/2));
            DIV3=POZ+r3min-1;
        else pom=[0 pom length(POZ)];
             DIV3=zeros(1,length(pom)-1);
             for i=1:length(pom)-1
                 ppom=(pom(i)+1):pom(i+1);
                 pppom=round(length(ppom)/2);
                 DIV3(i)=POZ(ppom(pppom));
             end
             DIV3=DIV3+r3min-1+i1-1;
        end
    end
end
%medialna linia+1

%1.objekt
DIV4=[];
for i=r1min:r1max
    PPP1(i-r1min+1)=L(i,liniax(i-liniay(1)+1)-15);
end
i1=find(PPP1,1,'first');
i2=find(PPP1,1,'last');
PPP1=PPP1(i1:i2);
POZ=find(PPP1==0);
if ~isempty(POZ)
    dpoz=diff(POZ);
    if isempty(dpoz)
       DIV4=POZ+r1min-1;
    else pom=find(dpoz~=1);
         if isempty(pom)
            POZ=POZ(round(length(POZ)/2));
            DIV4=POZ+r1min-1;
        else pom=[0 pom length(POZ)];
             DIV4=zeros(1,length(pom)-1);
             for i=1:length(pom)-1
                 ppom=(pom(i)+1):pom(i+1);
                 pppom=round(length(ppom)/2);
                 DIV4(i)=POZ(ppom(pppom));
             end
             DIV4=DIV4+r1min-1+i1-1;
        end
    end
end
%2.objekt
DIV5=[];
for i=r2min:r2max
    PPP2(i-r2min+1)=L(i,liniax(i-liniay(1)+1)-15);
end
i1=find(PPP2,1,'first');
i2=find(PPP2,1,'last');
PPP2=PPP2(i1:i2);
POZ=find(PPP2==0);
if ~isempty(POZ)
    dpoz=diff(POZ);
    if isempty(dpoz)
       DIV5=POZ+r2min-1;
    else pom=find(dpoz~=1);
         if isempty(pom)
            POZ=POZ(round(length(POZ)/2));
            DIV5=POZ+r2min-1;
        else pom=[0 pom length(POZ)];
             DIV5=zeros(1,length(pom)-1);
             for i=1:length(pom)-1
                 ppom=(pom(i)+1):pom(i+1);
                 pppom=round(length(ppom)/2);
                 DIV5(i)=POZ(ppom(pppom));
             end
             DIV5=DIV5+r2min-1+i1-1;
        end
    end
end
%3.objekt
DIV6=[];
for i=r3min:r3max
    PPP3(i-r3min+1)=L(i,liniax(i-liniay(1)+1)-15);
end
i1=find(PPP3,1,'first');
i2=find(PPP3,1,'last');
PPP3=PPP3(i1:i2);
POZ=find(PPP3==0);
if ~isempty(POZ)
    dpoz=diff(POZ);
    if isempty(dpoz)
       DIV6=POZ+r3min-1;
    else pom=find(dpoz~=1);
         if isempty(pom)
            POZ=POZ(round(length(POZ)/2));
            DIV6=POZ+r3min-1;
        else pom=[0 pom length(POZ)];
             DIV6=zeros(1,length(pom)-1);
             for i=1:length(pom)-1
                 ppom=(pom(i)+1):pom(i+1);
                 pppom=round(length(ppom)/2);
                 DIV6(i)=POZ(ppom(pppom));
             end
             DIV6=DIV6+r3min-1+i1-1;
        end
    end
end
%rozdelenie a uprava indexovania
DIV=[DIV1 DIV2 DIV3 DIV4 DIV5 DIV6];
DIV=unique(DIV);
pom=diff(DIV);
DIV(find(pom<3)+1)=[];
for i=1:length(DIV)
    L(DIV(i),:)=0;
    Lpom=L;
    Lpom(1:DIV(i),:)=0;
    poz=find(Lpom~=0);
    L(poz)=L(poz)+1;
end
ind=unique(L);
ind(1)=[];

%rozdelenie zlucenych oblasti (ostatne objekty)

for q=1:2
    %vyhladanie IVD
    plochy=zeros(1,length(ind)-3);
    plochy(1)=1000000000;
    for i=2:length(ind)-3
        plochy(i)=length(find(L==i));
    end
    [~,ivd]=min(plochy);
    %identifikacia oblasti pre rozdelenie
    por=zeros(1,length(ind)-3);
    i=ivd;
    while i>2
        if plochy(i)>plochy(i-1)
           por(i)=1;
        end
        i=i-2;
    end
    i=ivd-1;  
    while i>2
        if plochy(i)<plochy(i-1)
           por(i)=1;
        end
        i=i-2;
    end
    i=ivd;
    while i<length(ind)-3
        if plochy(i)>plochy(i+1)
           por(i)=1;
        end
        i=i+2;
    end
    i=ivd+1;
    while i<length(ind)-3
        if plochy(i)<plochy(i+1)
           por(i)=1;
        end
        i=i+2;
    end
    P1=find(por(1:ivd),1,'last');
    P2=find(por(ivd:end),1,'first')+ivd-1;
    %rozdelenie oblasti (smer nahor)
    if ~isempty(P1)
        S=[plochy(P1-1),plochy(P1),plochy(P1+1)];
        [~,index]=min(S);
        if index==1
           PP1=find(L==P1);
           PP2=find(L==P1+1);
        else PP1=find(L==P1-1);
             PP2=find(L==P1);
        end
        [PP1r,~]=ind2sub([z x],PP1);
        [PP2r,~]=ind2sub([z x],PP2);
        r1min=min(PP1r);
        r1max=max(PP1r);
        r2min=min(PP2r);
        r2max=max(PP2r);
        PPP1=zeros(1,length(r1min:r1max));
        PPP2=zeros(1,length(r2min:r2max));
        %medialna linia-2
        
        %1.objekt
        DIV1=[];
        for i=r1min:r1max
            PPP1(i-r1min+1)=L(i,liniax(i-liniay(1)+1)-18);
        end
        i1=find(PPP1,1,'first');
        i2=find(PPP1,1,'last');
        PPP1=PPP1(i1:i2);
        POZ=find(PPP1==0);
        if ~isempty(POZ)
            dpoz=diff(POZ);
            if isempty(dpoz)
               DIV1=POZ+r1min-1;
            else pom=find(dpoz~=1);
                 if isempty(pom)
                    POZ=POZ(round(length(POZ)/2));
                    DIV1=POZ+r1min-1;
                 else pom=[0 pom length(POZ)];
                      DIV1=zeros(1,length(pom)-1);
                      for i=1:length(pom)-1
                          ppom=(pom(i)+1):pom(i+1);
                          pppom=round(length(ppom)/2);   
                          DIV1(i)=POZ(ppom(pppom));
                      end
                      DIV1=DIV1+r1min-1+i1-1;
                 end
            end
        end 
        %2.objekt
        DIV2=[];
        for i=r2min:r2max
            PPP2(i-r2min+1)=L(i,liniax(i-liniay(1)+1)-18);
        end
        i1=find(PPP2,1,'first');
        i2=find(PPP2,1,'last');
        PPP2=PPP2(i1:i2);
        POZ=find(PPP2==0);
        if ~isempty(POZ)
            dpoz=diff(POZ);
            if isempty(dpoz)
               DIV2=POZ+r2min-1;
            else pom=find(dpoz~=1);
                 if isempty(pom)
                    POZ=POZ(round(length(POZ)/2));
                    DIV2=POZ+r2min-1;
                 else pom=[0 pom length(POZ)];
                      DIV2=zeros(1,length(pom)-1);
                      for i=1:length(pom)-1
                          ppom=(pom(i)+1):pom(i+1);
                          pppom=round(length(ppom)/2);   
                          DIV2(i)=POZ(ppom(pppom));
                      end
                      DIV2=DIV2+r2min-1+i1-1;
                 end
            end
        end    
        %medialna linia+2
        
        %1.objekt
        DIV3=[];
        for i=r1min:r1max
            PPP1(i-r1min+1)=L(i,liniax(i-liniay(1)+1)-14);
        end
        i1=find(PPP1,1,'first');
        i2=find(PPP1,1,'last');
        PPP1=PPP1(i1:i2);
        POZ=find(PPP1==0);
        if ~isempty(POZ)
            dpoz=diff(POZ);
            if isempty(dpoz)
               DIV3=POZ+r1min-1;
            else pom=find(dpoz~=1);
                 if isempty(pom)
                    POZ=POZ(round(length(POZ)/2));
                    DIV3=POZ+r1min-1;
                 else pom=[0 pom length(POZ)];
                      DIV3=zeros(1,length(pom)-1);
                      for i=1:length(pom)-1
                          ppom=(pom(i)+1):pom(i+1);
                          pppom=round(length(ppom)/2);   
                          DIV3(i)=POZ(ppom(pppom));
                      end
                      DIV3=DIV3+r1min-1+i1-1;
                 end
            end
        end 
        %2.objekt
        DIV4=[];
        for i=r2min:r2max
            PPP2(i-r2min+1)=L(i,liniax(i-liniay(1)+1)-14);
        end
        i1=find(PPP2,1,'first');
        i2=find(PPP2,1,'last');
        PPP2=PPP2(i1:i2);
        POZ=find(PPP2==0);
        if ~isempty(POZ)
            dpoz=diff(POZ);
            if isempty(dpoz)
               DIV4=POZ+r2min-1;
            else pom=find(dpoz~=1);
                 if isempty(pom)
                    POZ=POZ(round(length(POZ)/2));
                    DIV4=POZ+r2min-1;
                 else pom=[0 pom length(POZ)];
                      DIV4=zeros(1,length(pom)-1);
                      for i=1:length(pom)-1
                          ppom=(pom(i)+1):pom(i+1);
                          pppom=round(length(ppom)/2);   
                          DIV4(i)=POZ(ppom(pppom));
                      end
                      DIV4=DIV4+r2min-1+i1-1;
                 end
            end
        end    
        %rozdelenie a uprava indexovania
        DIV=[DIV1 DIV2 DIV3 DIV4];
        DIV=unique(DIV);
        pom=diff(DIV);
        DIV(find(pom<3)+1)=[];
        for i=1:length(DIV)
            L(DIV(i),:)=0;
            Lpom=L;
            Lpom(1:DIV(i),:)=0;
            poz=find(Lpom~=0);
            L(poz)=L(poz)+1;
        end
    end
    %rozdelenie oblasti(smer nadol)
    if ~isempty(P2)
        S=[plochy(P2-1),plochy(P2),plochy(P2+1)];
        [~,index]=min(S);
        if index==1
           PP1=find(L==P2);
           PP2=find(L==P2+1);
        else PP1=find(L==P2-1);
             PP2=find(L==P2);
        end
        [PP1r,~]=ind2sub([z x],PP1);
        [PP2r,~]=ind2sub([z x],PP2);
        r1min=min(PP1r);
        r1max=max(PP1r);
        r2min=min(PP2r);
        r2max=max(PP2r);
        PPP1=zeros(1,length(r1min:r1max));
        PPP2=zeros(1,length(r2min:r2max));
        %medialna linia-2
        
        %1.objekt
        DIV1=[];
        for i=r1min:r1max
            PPP1(i-r1min+1)=L(i,liniax(i-liniay(1)+1)-18);
        end
        i1=find(PPP1,1,'first');
        i2=find(PPP1,1,'last');
        PPP1=PPP1(i1:i2);
        POZ=find(PPP1==0);
        if ~isempty(POZ)
            dpoz=diff(POZ);
            if isempty(dpoz)
               DIV1=POZ+r1min-1;
            else pom=find(dpoz~=1);
                 if isempty(pom)
                    POZ=POZ(round(length(POZ)/2));
                    DIV1=POZ+r1min-1;
                 else pom=[0 pom length(POZ)];
                      DIV1=zeros(1,length(pom)-1);
                      for i=1:length(pom)-1
                          ppom=(pom(i)+1):pom(i+1);
                          pppom=round(length(ppom)/2);   
                          DIV1(i)=POZ(ppom(pppom));
                      end
                      DIV1=DIV1+r1min-1+i1-1;
                 end
            end
        end 
        %2.objekt
        DIV2=[];
        for i=r2min:r2max
            PPP2(i-r2min+1)=L(i,liniax(i-liniay(1)+1)-18);
        end
        i1=find(PPP2,1,'first');
        i2=find(PPP2,1,'last');
        PPP2=PPP2(i1:i2);
        POZ=find(PPP2==0);
        if ~isempty(POZ)
            dpoz=diff(POZ);
            if isempty(dpoz)
               DIV2=POZ+r2min-1;
            else pom=find(dpoz~=1);
                 if isempty(pom)
                    POZ=POZ(round(length(POZ)/2));
                    DIV2=POZ+r2min-1;
                 else pom=[0 pom length(POZ)];
                      DIV2=zeros(1,length(pom)-1);
                      for i=1:length(pom)-1
                          ppom=(pom(i)+1):pom(i+1);
                          pppom=round(length(ppom)/2);   
                          DIV2(i)=POZ(ppom(pppom));
                      end
                      DIV2=DIV2+r2min-1+i1-1;
                 end
            end
        end    
        %medialna linia+2
        
        %1.objekt
        DIV3=[];
        for i=r1min:r1max
            PPP1(i-r1min+1)=L(i,liniax(i-liniay(1)+1)-14);
        end
        i1=find(PPP1,1,'first');
        i2=find(PPP1,1,'last');
        PPP1=PPP1(i1:i2);
        POZ=find(PPP1==0);
        if ~isempty(POZ)
            dpoz=diff(POZ);
            if isempty(dpoz)
               DIV3=POZ+r1min-1;
            else pom=find(dpoz~=1);
                 if isempty(pom)
                    POZ=POZ(round(length(POZ)/2));
                    DIV3=POZ+r1min-1;
                 else pom=[0 pom length(POZ)];
                      DIV3=zeros(1,length(pom)-1);
                      for i=1:length(pom)-1
                          ppom=(pom(i)+1):pom(i+1);
                          pppom=round(length(ppom)/2);   
                          DIV3(i)=POZ(ppom(pppom));
                      end
                      DIV3=DIV3+r1min-1+i1-1;
                 end
            end
        end 
        %2.objekt
        DIV4=[];
        for i=r2min:r2max
            PPP2(i-r2min+1)=L(i,liniax(i-liniay(1)+1)-14);
        end
        i1=find(PPP2,1,'first');
        i2=find(PPP2,1,'last');
        PPP2=PPP2(i1:i2);
        POZ=find(PPP2==0);
        if ~isempty(POZ)
            dpoz=diff(POZ);
            if isempty(dpoz)
               DIV4=POZ+r2min-1;
            else pom=find(dpoz~=1);
                 if isempty(pom)
                    POZ=POZ(round(length(POZ)/2));
                    DIV4=POZ+r2min-1;
                 else pom=[0 pom length(POZ)];
                      DIV4=zeros(1,length(pom)-1);
                      for i=1:length(pom)-1
                          ppom=(pom(i)+1):pom(i+1);
                          pppom=round(length(ppom)/2);   
                          DIV4(i)=POZ(ppom(pppom));
                      end
                      DIV4=DIV4+r2min-1+i1-1;
                 end
            end
        end    
        %rozdelenie a uprava indexovania
        DIV=[DIV1 DIV2 DIV3 DIV4];
        DIV=unique(DIV);
        pom=diff(DIV);
        DIV(find(pom<3)+1)=[];
        for i=1:length(DIV)
            L(DIV(i),:)=0;
            Lpom=L;
            Lpom(1:DIV(i),:)=0;
            poz=find(Lpom~=0);
            L(poz)=L(poz)+1;
        end
    end
    ind=unique(L);
    ind(1)=[];
end

%opakovane vyhladanie IVD

plochy=zeros(1,length(ind)-3);
plochy(1)=1000000000;
for i=2:length(ind)-3
    plochy(i)=length(find(L==i));
end
[~,ivd]=min(plochy);

%pociatocne semienka 

if mod(ivd,2)==0
   f=2;
   b=1;
else f=1;
     b=2;
end
%objekt
seed1=[];
for i=f:2:length(ind)
    seed1=[seed1,(find(L==i))'];
end
%pozadie
seed2=[];
for i=b:2:length(ind)
    seed2=[seed2,(find(L==i))'];
end

%random walker 3D

%semienka objektu
[seed1z,seed1x]=ind2sub([z,x],seed1);
seed1=[];
for i=-5:5
seed1=[seed1,sub2ind([y,x,z],ones(1,length(seed1x))*(V+i),seed1x,seed1z)];
end
%semienka pozadia
[seed2z,seed2x]=ind2sub([z,x],seed2);
seed2=[];
for i=-5:5
seed2=[seed2,sub2ind([y,x,z],ones(1,length(seed2x))*(V+i),seed2x,seed2z)];
end
%algoritmus rw 3D
[mask,~] = myRandomWalker1(obr,[seed1,seed2],...
[ones(1,length(seed1)),zeros(1,length(seed2))],beta,P);

%odstranenie artefaktov

dmask=imfill(~mask,1);
dmask=imfill(dmask,y*(x-1)+1);
dmask=imfill(dmask,y*x*(z-1)+y);
dmask=imfill(dmask,y*x*z);
dmask=~dmask;
for i=1:length(N)
    rez=dmask(:,:,N(i));
    [sy,sx]=find(rez);
    sy(sx<liniax(i))=[];
    sx(sx<liniax(i))=[];
    rez(sub2ind([y,x],sy,sx))=0;
    dmask(:,:,N(i))=rez;
end
for i=N(end)+1:z
    dmask(:,:,i)=0;
end
for i=1:z
    rez=dmask(:,:,i);
    [sy,sx]=find(rez);
    sy(sx<max(liniax))=[];
    sx(sx<max(liniax))=[];
    rez(sub2ind([y,x],sy,sx))=0;
    dmask(:,:,i)=rez;
end
dmask=bwareaopen(dmask,50);

%-------------------------------VYHLADENIE---------------------------------

for i=1:IDdil
    
    %dilatacia
    
    dil=imdilate(dmask,se);
    dil(dmask==1)=0;
    pos=find(dil==1);
    kor=ones(y,x,z);
    kor(pos)=oo(pos);
    kor(kor<=prah(1))=0;
    dmask(kor==0)=1;
    
    %vyplnenie medzier
    
    for j=1:z
        mask=dmask(:,:,j);
        maski=find(mask);
        if isempty(maski)
            continue
        end
        [masky,maskx]=ind2sub([y x],maski);
        %centroid
        Yc=mean(masky);
        Xc=mean(maskx);
        %vyplnenie medzier - axialne rezy
        mask=imfill(mask,'holes');
        for k=round(Yc)-15:1:round(Yc)+15
            if k<2 || k>y-1
               continue
            end
            for l=round(Xc)-15:1:round(Xc)+15
                if l<2 || l>x-1
                   continue
                end
                okolie=[mask(k-1,l-1),mask(k-1,l),mask(k-1,l+1),mask(k,l-1),...
                       mask(k,l+1),mask(k+1,l-1),mask(k+1,l),mask(k+1,l+1)];
                if sum(okolie)>4
                   mask(k,l)=1;
                end
            end
        end
        dmask(:,:,j)=mask;
    end
    for j=1:y
        mask=dmask(j,:,:);
        [mask,~]=shiftdim(mask);
        maski=find(mask,1);
        if isempty(maski)
            continue
        end
        %vyplnenie medzier - sagitalne rezy
        mask=imfill(mask,'holes');
        for k=2:z-1
            for l=2:x-1
                okolie=[mask(k-1,l-1),mask(k-1,l),mask(k-1,l+1),mask(k,l-1),...
                       mask(k,l+1),mask(k+1,l-1),mask(k+1,l),mask(k+1,l+1)];
                if sum(okolie)>4
                   mask(k,l)=1;
                end
            end
        end
        dmask(j,:,:)=shiftdim(mask,-1);
    end
end

%dodatocne vyhladenie - sagitalne rezy

for i=1:10
    for j=1:y
        mask=dmask(j,:,:);
        [mask,~]=shiftdim(mask);
        maski=find(mask,1);
        if isempty(maski)
            continue
        end
        %vyplnenie medzier 
        mask=imfill(mask,'holes');
        for k=2:z-1
            for l=2:x-1
                okolie=[mask(k-1,l-1),mask(k-1,l),mask(k-1,l+1),mask(k,l-1),...
                       mask(k,l+1),mask(k+1,l-1),mask(k+1,l),mask(k+1,l+1)];
                if sum(okolie)>4
                   mask(k,l)=1;
                end
            end
        end
        dmask(j,:,:)=shiftdim(mask,-1);
    end
end
dmask=double(dmask);

%zobrazenie

%axialna rovina
va=oo;
va(dmask==1)=max(max(max(oo)));
figure('Name','Intervertebrálne disky (axiálna rovina)')
imshow3D(va)
%sagitalna rovina
vs=zeros(z,x,y);
for j=1:y
    a=va(j,:,:);
    aa=zeros(z,x);
    for i=1:z
        aa(i,:)=a(1,:,i);
    end
    vs(:,:,j)=aa;
end
figure('Name','Intervertebrálne disky (sagitálna rovina)')
imshow3D(vs)
end