clear all; close all; clc
tic
% load('obr_na_segmentaciu_br.mat')
load('obr_na_segmentaciu_br2.mat')
% load('obr_na_segmentaciu_br3.mat')
% load('obr_na_segmentaciu_br4.mat')

a = obraz;

% imshow(a,[])
[s1, s2] = size(a);

%% LBP
metoda_lbp = 1; % naprogramovana funkcia R = 1, N = 8
% metoda_lbp = 2; % naprogramovana funkcia R = 2, N = 12
% metoda_lbp = 3; % naprogramovana funkcia R = 3, N = 16
% metoda_lbp = 4; % prevzata funkcia LBP - potrebne definovat R, N
% metoda_lbp = 5; % prevzata funkcia LBP u2 - potrebne definovat R, N
% metoda_lbp = 6; % prevzata funkcia LBP ri - potrebne definovat R, N
% metoda_lbp = 7; % prevzata funkcia LBP riu2 - potrebne definovat R, N
R = 1; N = 8;

%% kontrast
metoda_kontrast = 1; % naprogramovana funkcia R = 1, N = 8
% metoda_kontrast = 2; % naprogramovana funkcia R = 2, N = 12
% metoda_kontrast = 3; % naprogramovana funkcia R = 3, N = 16
% metoda_kontrast = 4; % prevzata funkcia - potrebne definovat R, N

%% počet kvantizačných úrovní kontrastu
pq = 4; % pre 1,1

%% spôsob kvantizácie
% q_k = 1; % kvantizácia celého rozsahu od min po max
% q_k = 2; % kvantizácia rozsahu <median-1.5*IQR, median+1.5*IQR>
q_k = 3; % kvantizácia podľa rozdelenia celych dat - percentily
% q_k = 4; % kvantizácia rozsahu <kvantil1-1.5*IQR, kvantil3+1.5*IQR>

%% konštanta
K = 100;

%% prah pre delenie oblastí
X = 1.4; % prah pre delenie oblastí na kvadranty

%% prah pre spájanie oblastí
Yn = 2.5;

%% veľkosť oblasti pri klasifikácii hraničných pixelov
r_lbp = 25; 

%% veľkosť najmenších elementov po hierarchickom delení
s_min = 20; 
s_max = s_min*4;

%% zmena veľkosti obrazu v závislosti na zvolenom s_min
% orezanie
rd = (s1-floor(s1/s_max)*s_max)/2; % rozdiel riadky
cd = (s2-floor(s2/s_max)*s_max)/2; % rozdiel stĺpce
a = a(1+rd : end - rd, 1+cd : end - cd);
a = a(1:floor(s1/s_max)*s_max, 1:floor(s2/s_max)*s_max);

% % doplnenie
% rd = s1/s_max; % rozdiel riadky
% cd = s2/s_max; % rozdiel stĺpce
% 
% if rd - floor(rd)== 0
%     rd = rd-1;
%     aar = ((rd+1)*s_max-s1)/2;
% else
%     rd = floor(rd);
%     aar = ((rd+1)*s_max-s1)/2;
% end
% 
% if cd - floor(cd) == 0
%     cd = cd-1;
%     aac = ((cd+1)*s_max-s2)/2;
% else
%     cd = floor(cd);
%     aac = ((cd+1)*s_max-s2)/2;
% end
% aa = zeros((rd+1)*s_max, (cd+1)*s_max);
% 
% aa(aar+1:end-aar, aac+1:end-aac) = a;
% aa(1:aar, aac+1:end-aac) = flipud(a(1:aar, 1:end));
% aa(end-aar:end, aac+1:end-aac) = flipud(a(end-aar:end, 1:end)); 
% aa(1:end, 1:aac) = fliplr(aa(1:end, aac+1:2*aac));
% aa(1:end, end-aac:end) = fliplr(aa(1:end, end-2*aac:end-aac));
% 
% a = aa;

[s1, s2] = size(a);

%% analyza distribucie kontrastu/VAR v celom obraze

switch(metoda_kontrast)
    case 1
        C_z = kontrast_8_1(a);
    case 2
        C_z = kontrast_12_2(a);
    case 3
        C_z = kontrast_16_3(a);
    case 4
        C_z = cont(a,R,N);
end

%% vytvorenie intervalov pre kvantovanie podľa spôsobu kvantovania
C_z = C_z'; C_z = C_z(:)'; % prevod matice priznakov na vektor

switch(q_k)
    case 2
        iqr_C = iqr(C_z);
    	C_median = median(C_z);
        int_Cm1 = C_median-1.5*iqr_C;
        int_Cm2 = C_median+1.5*iqr_C;
        krok_m = (int_Cm2-int_Cm1)/pq;
        C_med = int_Cm1:krok_m:int_Cm2;
        C_med(1) = min(min(C_z));
        C_med(end) = max(max(C_z));
    case 3
        krok_p = 100/(pq);
        C_perc = prctile(C_z,krok_p:krok_p:100-krok_p);
        C_perc = [min(min(C_z)), C_perc, max(max(C_z))];
    case 4
        iqr_C = iqr(C_z);
        k1 = quantile(C_z,0.25);
        k3 = quantile(C_z,0.75);
        int_Ck1 = k1-1.5*iqr_C;
        int_Ck2 = k3+1.5*iqr_C;
        krok_k = (int_Ck2-int_Ck1)/pq;
        C_k = int_Ck1:krok_k:int_Ck2;
        C_k(1) = min(min(C_z));
        C_k(end) = max(max(C_z));
end

%% hierarchické delenie obrazu
% X = 1.3; % prah pre delenie oblastí na kvadranty
% imshow(a,[]);

b = zeros(s1,s2);

s11 = s1/s_min;
s21 = s2/s_min;

s12 = s1/s_max;
s22 = s2/s_max;
po = s11*s22; % pocet moznych oblasti;

krok = (s_max/s_min)^2;
pozicie = [4,1];
p_pozicie = 1;
k = 0;

% rozdelenie obrazu na bloky počiatočnej veľkosti s_max
for i = 1:s12 
    for j = 1:s22
        A = 1+(i-1)*s_max;
        B = 1+(j-1)*s_max;
        b(A: A+s_max - 1, B: B+s_max - 1) = k*krok;
        k = k+1;
    end
end

b_max = max(max(b));
zoznam = (0:krok:b_max);
zoznam2 = zeros(1,1);
tr = s_max;
while tr > 2*s_min
    pp = 1;
    for m = zoznam % prechadzanie stvorcov 64 x 64
        [r, c] = find(b == m);
        r(2:end-1) = [];
        c(2:end-1) = [];

        % delenie stvorcov az do splnenia podmienky alebo dosiahnutia s_min
        
        t = a(r(1):r(2), c(1):c(2));
        [tr, tc] = size(t);
        [t1, t2, t3, t4] = stvrtenie(t, tr, tc);

        switch(metoda_lbp)
            case 1
                lbp_m1 = lbp_8_1(t1);
                lbp_m2 = lbp_8_1(t2);
                lbp_m3 = lbp_8_1(t3);
                lbp_m4 = lbp_8_1(t4);
                lbp_q = 256;
            case 2
                lbp_m1 = lbp_12_2(t1);
                lbp_m2 = lbp_12_2(t2);
                lbp_m3 = lbp_12_2(t3);
                lbp_m4 = lbp_12_2(t4);
                lbp_q = 4096;
            case 3
                lbp_m1 = lbp_16_3(t1);
                lbp_m2 = lbp_16_3(t2);
                lbp_m3 = lbp_16_3(t3);
                lbp_m4 = lbp_16_3(t4);

                lbp_m1 = kvantizacia_lbp(lbp_m1, 4096);
                lbp_m2 = kvantizacia_lbp(lbp_m2, 4096);
                lbp_m3 = kvantizacia_lbp(lbp_m3, 4096);
                lbp_m4 = kvantizacia_lbp(lbp_m4, 4096);
                lbp_q = 4096;
            case 4
                lbp_m1 = lbp(t1,R,N,0,0);
                lbp_m2 = lbp(t2,R,N,0,0);
                lbp_m3 = lbp(t3,R,N,0,0);
                lbp_m4 = lbp(t4,R,N,0,0);
                lbp_q = 2^N;
            case 5
                mapping = getmapping(N,'u2');
                lbp_m1 = lbp(t1,R,N,mapping,0);
                lbp_m2 = lbp(t2,R,N,mapping,0);
                lbp_m3 = lbp(t3,R,N,mapping,0);
                lbp_m4 = lbp(t4,R,N,mapping,0);
                lbp_q = mapping.num;
            case 6
                mapping = getmapping(N,'ri');
                lbp_m1 = lbp(t1,R,N,mapping,0);
                lbp_m2 = lbp(t2,R,N,mapping,0);
                lbp_m3 = lbp(t3,R,N,mapping,0);
                lbp_m4 = lbp(t4,R,N,mapping,0);
                lbp_q = mapping.num;
            case 7
                mapping = getmapping(N,'riu2');
                lbp_m1 = lbp(t1,R,N,mapping,0);
                lbp_m2 = lbp(t2,R,N,mapping,0);
                lbp_m3 = lbp(t3,R,N,mapping,0);
                lbp_m4 = lbp(t4,R,N,mapping,0);
                lbp_q = mapping.num;
        end
        
        switch(metoda_kontrast)
            case 1
                c_m1 = kontrast_8_1(t1);
                c_m2 = kontrast_8_1(t2);
                c_m3 = kontrast_8_1(t3);
                c_m4 = kontrast_8_1(t4);
            case 2
                c_m1 = kontrast_12_2(t1);
                c_m2 = kontrast_12_2(t2);
                c_m3 = kontrast_12_2(t3);
                c_m4 = kontrast_12_2(t4);
            case 3
                c_m1 = kontrast_16_3(t1);
                c_m2 = kontrast_16_3(t2);
                c_m3 = kontrast_16_3(t3);
                c_m4 = kontrast_16_3(t4);
            case 4
                c_m1 = cont(t1,R,N);
                c_m2 = cont(t2,R,N);
                c_m3 = cont(t3,R,N);
                c_m4 = cont(t4,R,N);
        end
        
        switch(q_k)
            case 1
            c_m1 = kvantizacia_C(c_m1, pq);
            c_m2 = kvantizacia_C(c_m2, pq);
            c_m3 = kvantizacia_C(c_m3, pq);
            c_m4 = kvantizacia_C(c_m4, pq);  
            case 2
            c_m1 = kvantizacia_C_median_iqr(c_m1, C_med);
            c_m2 = kvantizacia_C_median_iqr(c_m2, C_med);
            c_m3 = kvantizacia_C_median_iqr(c_m3, C_med);
            c_m4 = kvantizacia_C_median_iqr(c_m4, C_med);  
            case 3
            c_m1 = kvantizacia_C_percentily(c_m1, C_perc);
            c_m2 = kvantizacia_C_percentily(c_m2, C_perc);
            c_m3 = kvantizacia_C_percentily(c_m3, C_perc);
            c_m4 = kvantizacia_C_percentily(c_m4, C_perc);  
            case 4
            c_m1 = kvantizacia_C_kvantily_iqr(c_m1, C_k);
            c_m2 = kvantizacia_C_kvantily_iqr(c_m2, C_k);
            c_m3 = kvantizacia_C_kvantily_iqr(c_m3, C_k);
            c_m4 = kvantizacia_C_kvantily_iqr(c_m4, C_k);  
        end
        
        h1 = (joint_histogram_d(c_m1, lbp_m1, pq, lbp_q))*K; % kontrast vertikalna os, lbp horizontalna
        h2 = (joint_histogram_d(c_m2, lbp_m2, pq, lbp_q))*K; % kontrast vertikalna os, lbp horizontalna
        h3 = (joint_histogram_d(c_m3, lbp_m3, pq, lbp_q))*K; % kontrast vertikalna os, lbp horizontalna
        h4 = (joint_histogram_d(c_m4, lbp_m4, pq, lbp_q))*K; % kontrast vertikalna os, lbp horizontalna
        
        h1(h1<1)=1;
        h2(h2<1)=1;
        h3(h3<1)=1;
        h4(h4<1)=1;
        
        G = zeros(1,6);
        G(1,1) = g_stat(h1, h2);
        G(1,2) = g_stat(h1, h3);
        G(1,3) = g_stat(h1, h4);
        G(1,4) = g_stat(h2, h3);
        G(1,5) = g_stat(h2, h4);
        G(1,6) = g_stat(h3, h4);

        G_max = max(G);
        G_min = min(G);

        p1 = m;
        R = G_max/G_min;

        if R > X
           b(r(1):r(1)+tr/2-1, c(1):c(1)+tr/2-1) = p1+0*pozicie(1,p_pozicie);
           b(r(1):r(1)+tr/2-1, c(1)+tr/2:c(1)+tr-1)= p1+1*pozicie(1,p_pozicie);
           b(r(1)+tr/2:r(1)+tr-1, c(1):c(1)+tr/2-1)= p1+2*pozicie(1,p_pozicie);
           b(r(1)+tr/2:r(1)+tr-1, c(1)+tr/2:c(1)+tr-1)= p1+3*pozicie(1,p_pozicie);
           
           zoznam2(1,pp) = p1+0*pozicie(1,p_pozicie);
           zoznam2(1,pp+1) = p1+1*pozicie(1,p_pozicie);
           zoznam2(1,pp+2) = p1+2*pozicie(1,p_pozicie);
           zoznam2(1,pp+3) = p1+3*pozicie(1,p_pozicie);
           pp = pp+4;
        end
    end
    zoznam = zoznam2;
    p_pozicie = p_pozicie+1;
    if zoznam == 0
        break
    end
end

%% vykreslenie výsledkov hierarchického delenia
% s = b;
% mx = max(s(:));
% cmap = hsv2rgb([(0:mx-1)'./max(mx,1) ones(mx,2)]);
% 
% imshow(s,cmap)
% colorbar
% datacursormode on
% % 
BW1 = edge(b,'canny',0.000001);
% 
% figure
% imshow(BW1);

%% spájanie oblastí
list = (unique(b))'; %zoznam všetkých hodnôt (segmentov) v matici
lb = length(list);
% algoritmus prehlada susedov smermi vpravo a dole - posledny stlpec vpravo
% a posledny riadok sú prepísané na NaN, nie na 0, pretože 0 môže byť index
% segmentu
b_n = b;
% b_n(1,:) = NaN;
b_n(end,:) = NaN;
% b_n(:,1) = NaN;
b_n(:,end) = NaN;

% zoznam_s = zeros(s11*s21*4,4);
zoznam_s = zeros(lb,8);
p_z = 1; % pozicia v zozname
for m = list % prechádzanie všetkých oblasí
    [r_n, c_n] = find(b_n == m);
    r_n(2:end-1) = [];
    c_n(2:end-1) = [];
    r_n(2) = r_n(2)+1; % suradnice riadkov okolia
    c_n(2) = c_n(2)+1; % suradnice stlpcov okolia
    okolie = b(r_n(1):r_n(2), c_n(1):c_n(2));
    h = unique(okolie); % zoznam hodnôt vyskytujúcich sa v matici = hodnoty
    % segmentu + hodnoty susedov
    % vynechanie susedov po diagonalach
    hs = histc(okolie(:),h); % počet výskytov daných hodnôt v matici
    ind = find(hs == 1);
    h(ind) = [];
    h(h == m) = []; % vynechanie hodnoty testovaneho segmentu
    h = h';
    ukoncenie = isempty(h);
    if ukoncenie == 1;
       break
    end
    [r, c] = find(b == m);
    r(2:end-1) = [];
    c(2:end-1) = [];
    
    o1 = a(r(1):r(2), c(1):c(2));
    o1_s = size(o1);
    
    switch(metoda_lbp)
        case 1
            o1_lbp = lbp_8_1(o1);
            lbp_q = 256;
        case 2
            o1_lbp = lbp_12_2(o1);
            lbp_q = 4096;
        case 3
            o1_lbp = lbp_16_3(o1);
            o1_lbp = kvantizacia_lbp(o1_lbp, 4096);
            lbp_q = 4096;
        case 4
            o1_lbp = lbp(o1,R,N,0,0);
            lbp_q = 2^N;
        case 5
            mapping = getmapping(N,'u2');
            o1_lbp = lbp(o1,R,N,mapping,0);
            lbp_q = mapping.num;
        case 6
            mapping = getmapping(N,'ri');
            o1_lbp = lbp(o1,R,N,mapping,0);
            lbp_q = mapping.num;
        case 7
            mapping = getmapping(N,'riu2');
            o1_lbp = lbp(o1,R,N,mapping,0);
            lbp_q = mapping.num;
    end

    switch(metoda_kontrast)
        case 1
            o1_c = kontrast_8_1(o1);
        case 2
            o1_c = kontrast_12_2(o1);
        case 3
            o1_c = kontrast_16_3(o1);
        case 4
            o1_c = cont(o1,R,N);
    end

    switch(q_k)
        case 1
        o1_c = kvantizacia_C(o1_c, pq);
        case 2
        o1_c = kvantizacia_C_median_iqr(o1_c, C_med);
        case 3
        o1_c = kvantizacia_C_percentily(o1_c, C_perc);
        case 4
        o1_c = kvantizacia_C_kvantily_iqr(o1_c, C_k);
    end
    
    o_h1 = (joint_histogram_d(o1_c, o1_lbp, pq, lbp_q))*K;  % lbp histogram testovaneho segmentu
    
    for n = h
        [rt, ct] = find(b == n);
        rt(2:end-1) = [];
        ct(2:end-1) = [];
        
        o = o1;
        o_h = o_h1; % lbp histogram testovaneho segmentu
        o_s = o1_s;
        
        sused = a(rt(1):rt(2), ct(1):ct(2));
        sused_s = size(sused);
        
        switch(metoda_lbp)
            case 1
                sused_lbp = lbp_8_1(sused);
                lbp_q = 256;
            case 2
                sused_lbp = lbp_12_2(sused);
                lbp_q = 4096;
            case 3
                sused_lbp = lbp_16_3(sused);
                sused_lbp = kvantizacia_lbp(sused_lbp, 4096);
                lbp_q = 4096;
            case 4
                sused_lbp = lbp(sused,R,N,0,0);
                lbp_q = 2^N;
            case 5
                mapping = getmapping(N,'u2');
                sused_lbp = lbp(sused,R,N,mapping,0);
                lbp_q = mapping.num;
            case 6
                mapping = getmapping(N,'ri');
                sused_lbp = lbp(sused,R,N,mapping,0);
                lbp_q = mapping.num;
            case 7
                mapping = getmapping(N,'riu2');
                sused_lbp = lbp(sused,R,N,mapping,0);
                lbp_q = mapping.num;
        end

        switch(metoda_kontrast)
            case 1
                sused_c = kontrast_8_1(sused);
            case 2
                sused_c = kontrast_12_2(sused);
            case 3
                sused_c = kontrast_16_3(sused);
            case 4
                sused_c = cont(sused,R,N);
        end

        switch(q_k)
            case 1
            sused_c = kvantizacia_C(sused_c, pq);
            case 2
            sused_c = kvantizacia_C_median_iqr(sused_c, C_med);
            case 3
            sused_c = kvantizacia_C_percentily(sused_c, C_perc);
            case 4
            sused_c = kvantizacia_C_kvantily_iqr(sused_c, C_k);
        end

        sused_h = (joint_histogram_d(sused_c, sused_lbp, pq, lbp_q))*K; % histogram suseda

        pps = o1_s(1)^2; % pocet pixelov segment
        ppn = sused_s(1)^2; % pocet pixelov sused
        
        o_h(o_h<1) = 1;
        sused_h(sused_h<1) = 1;
        G_S = g_stat(o_h, sused_h); % miera rozdielnosti medzi segm. a susedom
%         min_s = (min([o_s(1,1), sused_s(1,1)]))^2;
        min_s = min(pps, ppn);
        zoznam_s(p_z,1) = G_S;
        zoznam_s(p_z,2) = m;
        zoznam_s(p_z,3) = n;
        zoznam_s(p_z,4) = pps/s_min^2;
        zoznam_s(p_z,5) = ppn/s_min^2;
        % miera rozdielnosti * počet pix. menšej z oblastí
        zoznam_s(p_z,6) = G_S*min_s;
        
        zoznam_s(p_z,7) = pps; 
        zoznam_s(p_z,8) = ppn; 
        p_z = p_z+1;
    end
end

% nájdenie MI min - aktuálna dvojica, ktorá sa zlúči do jedného segmentu
min_z_s = zoznam_s(:,6);
min_z = min(min_z_s);
mi_p = 1;
MI = zeros(1,1000);
MI(1) = min_z;
Y = 10;
% Yn = 2;
MIR = min_z/max(MI);
MIR_z = zeros(1,1000);
p10 = round(0.1*length(min_z_s));

while MIR <= Y
    min_z_s = zoznam_s(:,6);
    min_z = min(min_z_s);
    MIR = min_z/max(MI);
    if MIR >= Y
       MIR_z(1,mi_p) = MIR;
        break
    end
    if mi_p > p10
        Y = Yn;
    end
    MI(1,mi_p) = min_z;
    MIR_z(1,mi_p) = MIR; % zoznam hodnôt MIR
    mi_p = mi_p+1;
    [rs, cs] = find(min_z_s == min_z);

    ms = zoznam_s(rs,2:3);
    ms = ms(1,:); % uprava v pripade ze v zoznam_s je viac riadkov s rovnakym min_z

    des_c1 = ms(1)-floor(ms(1));
    des_c2 = ms(2)-floor(ms(2));
    
    if des_c1 == 0 && des_c2 == 0
       [ri, ci] = find(b == ms(1));
       ri(2:end-1) = [];
       ci(2:end-1) = [];
       b(ri(1):ri(2), ci(1):ci(2))= ms(1)+0.001;
       [ri, ci] = find(b == ms(2));
       ri(2:end-1) = [];
       ci(2:end-1) = [];
       b(ri(1):ri(2), ci(1):ci(2))= ms(1)+0.002;
       v = [ms(1)+0.001, ms(1)+0.002];
       % vymazanie spojenych segmentov zo zoznam_s
       for i = ms
           [rz, cz] = find(zoznam_s == i);
           for j = rz
               zoznam_s(rz,:) = [];
           end
       end
       
    elseif des_c1 ~= 0 && des_c2 == 0
        v1 = unique(b);
        v1 = v1(v1>=floor(ms(1)) & v1<=ceil(ms(1))-0.001);

        [ri, ci] = find(b == ms(2));
        ri(2:end-1) = [];
        ci(2:end-1) = [];
        b(ri(1):ri(2), ci(1):ci(2))= max(v1)+0.001;

        % zoznam hodnôt ktorými je označený nový segment
        v = unique(b);
        v = v(v>=floor(ms(1)) & v<=ceil(ms(1))-0.001);
        v = v';
        % alternativa                                                otestovat!
    %     v = [v1, max(v1)+0.001];
        v1 = v1';
        vv = unique([v1, ms]);
        for i = vv
            [rz, cz] = find(zoznam_s == i);
            for j = rz
                zoznam_s(rz,:) = [];
            end
        end
        
    elseif des_c2 ~= 0 && des_c1 == 0
        v1 = unique(b);
        v1 = v1(v1>=floor(ms(2)) & v1<=ceil(ms(2))-0.001);

        [ri, ci] = find(b == ms(1));
        ri(2:end-1) = [];
        ci(2:end-1) = [];
        b(ri(1):ri(2), ci(1):ci(2))= max(v1)+0.001;

        % zoznam hodnôt ktorými je označený nový segment
        v = unique(b);
        v = v(v>=floor(ms(2)) & v<=ceil(ms(2))-0.001);
        v = v';
        % alternativa                                                otestovat!
    %     v = [v1, max(v1)+0.001];
        % vymazanie spojených segmentov zo zoznam_s
        v1 = v1';
        vv = unique([v1, ms]);
        for i = vv
            [rz, cz] = find(zoznam_s == i);
            for j = rz
                zoznam_s(rz,:) = [];
            end
        end
        
    elseif des_c1 ~= 0 && des_c2 ~= 0
        v1 = unique(b);
        v2 = v1;
        v1 = v1(v1>=floor(ms(1)) & v1<=ceil(ms(1))-0.001);
        v2 = v2(v2>=floor(ms(2)) & v2<=ceil(ms(2))-0.001);
        ii = 1;
        v2 = v2';
        v1 = v1';
        % prepis označenia druhého z pôvodnych segmentov
        for i = v2
            [ri, ci] = find(b == i);
             ri(2:end-1) = [];
             ci(2:end-1) = [];
             b(ri(1):ri(2), ci(1):ci(2))= max(v1)+ii*0.001;
             ii = ii+1;
        end
        % zoznam hodnôt ktorými je označený nový segment
        v = unique(b);
        v = v(v>=floor(ms(1)) & v<=ceil(ms(1))-0.001);
        v = v';
        
         % vymazanie spojenych segmentov zo zoznam_s
        vv = [v1, v2];
        for i = vv
            [rz, cz] = find(zoznam_s == i);
            for j = rz
                zoznam_s(rz,:) = [];
            end
        end
    end
    
    switch(q_k)
        case 2
            kvant = C_med;
        case 3
            kvant = C_perc;
        case 4
            kvant = C_k;
    end
    
    [susedia,h_s, h_s_s] = n_susedia_s(a,b,v, pq, K, q_k, metoda_lbp, metoda_kontrast, kvant, R, N, lbp_q);
    % susedia - zoznam susedov
    % h_s - histogram nového segmentu
    % h_s_s - počet pixelov nového segmentu

    [p_z, p_z2] = size(zoznam_s);
    p_z = p_z+1;
    for i = susedia
        des_c = i-floor(i);
        if des_c ~= 0
           sus = unique(b);
           sus = sus(sus>=floor(i) & sus<=ceil(i)-0.001);
           sus = sus';
           n_sused_s = 0;
           n_sused_h = zeros(pq,lbp_q);
           for j = sus
                [r_sus, c_sus] = find(b == j);
                r_sus(2:end-1) = [];
                c_sus(2:end-1) = [];
                n_sus = a(r_sus(1):r_sus(2), c_sus(1):c_sus(2));
                               
                switch(metoda_lbp)
                    case 1
                        n_sus_lbp = lbp_8_1(n_sus);
                        lbp_q = 256;
                    case 2
                        n_sus_lbp = lbp_12_2(n_sus);
                        lbp_q = 4096;
                    case 3
                        n_sus_lbp = lbp_16_3(n_sus);
                        n_sus_lbp = kvantizacia_lbp(n_sus_lbp, 4096);
                        lbp_q = 4096;
                    case 4
                        n_sus_lbp = lbp(n_sus,R,N,0,0);
                        lbp_q = 2^N;
                    case 5
                        mapping = getmapping(N,'u2');
                        n_sus_lbp = lbp(n_sus,R,N,mapping,0);
                        lbp_q = mapping.num;
                    case 6
                        mapping = getmapping(N,'ri');
                        n_sus_lbp = lbp(n_sus,R,N,mapping,0);
                        lbp_q = mapping.num;
                    case 7
                        mapping = getmapping(N,'riu2');
                        n_sus_lbp = lbp(n_sus,R,N,mapping,0);
                        lbp_q = mapping.num;
                end

                switch(metoda_kontrast)
                    case 1
                        n_sus_c = kontrast_8_1(n_sus);
                    case 2
                        n_sus_c = kontrast_12_2(n_sus);
                    case 3
                        n_sus_c = kontrast_16_3(n_sus);
                    case 4
                        n_sus_c = cont(n_sus,R,N);
                end

                switch(q_k)
                    case 1
                    n_sus_c = kvantizacia_C(n_sus_c, pq);
                    case 2
                    n_sus_c = kvantizacia_C_median_iqr(n_sus_c, C_med);
                    case 3
                    n_sus_c = kvantizacia_C_percentily(n_sus_c, C_perc);
                    case 4
                    n_sus_c = kvantizacia_C_kvantily_iqr(n_sus_c, C_k);
                end   
                
                n_sus_h = (joint_histogram_d(n_sus_c, n_sus_lbp, pq, lbp_q))*K; % histogram suseda
                n_sused_h = n_sused_h + n_sus_h;
                s = size(n_sus);
                n_sused_s = n_sused_s+(s(1))^2;
           end
        else
          [r_su, c_su] = find(b == i);
          r_su(2:end-1) = [];
          c_su(2:end-1) = [];
          n_sused = a(r_su(1):r_su(2), c_su(1):c_su(2));
          
            switch(metoda_lbp)
                case 1
                    n_sused_lbp = lbp_8_1(n_sused);
                    lbp_q = 256;
                case 2
                    n_sused_lbp = lbp_12_2(n_sused);
                    lbp_q = 4096;
                case 3
                    n_sused_lbp = lbp_16_3(n_sused);
                    n_sused_lbp = kvantizacia_lbp(n_sused_lbp, 4096);
                    lbp_q = 4096;
                case 4
                    n_sused_lbp = lbp(n_sused,R,N,0,0);
                    lbp_q = 2^N;
                case 5
                    mapping = getmapping(N,'u2');
                    n_sused_lbp = lbp(n_sused,R,N,mapping,0);
                    lbp_q = mapping.num;
                case 6
                    mapping = getmapping(N,'ri');
                    n_sused_lbp = lbp(n_sused,R,N,mapping,0);
                    lbp_q = mapping.num;
                case 7
                    mapping = getmapping(N,'riu2');
                    n_sused_lbp = lbp(n_sused,R,N,mapping,0);
                    lbp_q = mapping.num;
            end

            switch(metoda_kontrast)
                case 1
                    n_sused_c = kontrast_8_1(n_sused);
                case 2
                    n_sused_c = kontrast_12_2(n_sused);
                case 3
                    n_sused_c = kontrast_16_3(n_sused);
                case 4
                    n_sused_c = cont(n_sused,R,N);
            end

            switch(q_k)
                case 1
                n_sused_c = kvantizacia_C(n_sused_c, pq);
                case 2
                n_sused_c = kvantizacia_C_median_iqr(n_sused_c, C_med);
                case 3
                n_sused_c = kvantizacia_C_percentily(n_sused_c, C_perc);
                case 4
                n_sused_c = kvantizacia_C_kvantily_iqr(n_sused_c, C_k);
            end

          n_sused_h = (joint_histogram_d(n_sused_c, n_sused_lbp, pq, lbp_q))*K;  % lbp histogram testovaneho segmentu
          n_sused_s = size(n_sused);
          n_sused_s = (n_sused_s(1))^2;
        end
        
        pps = h_s_s;
        ppn = n_sused_s;
        h_s1 = h_s;
        
        h_s1(h_s1<1) = 1;
        n_sused_h(n_sused_h<1) = 1;
        G_S_n = g_stat(h_s1, n_sused_h);
        min_s_n = min([pps, ppn]);
        zoznam_s(p_z,1) = G_S_n;
        zoznam_s(p_z,2) = v(1);
        zoznam_s(p_z,3) = i;
        
        zoznam_s(p_z,4) = pps/s_min^2;
        zoznam_s(p_z,5) = ppn/s_min^2;
        zoznam_s(p_z,6) = G_S_n*min_s_n;
        
        zoznam_s(p_z,7) = pps; 
        zoznam_s(p_z,8) = ppn; 
        p_z = p_z+1;
    end
%     b_round = round(b);
%     b_ru = unique(b_round);
%     b_rul = length(b_ru);
end
toc

bf = floor(b);
BW2 = edge(bf,'canny',0.000001);
imshow(BW2);

%% upresnenie hraníc segmentov - pixelwise classification
tic

% a = a(aar+1:end-aar, aac+1:end-aac);
% b = b(aar+1:end-aar, aac+1:end-aac);
% bf = bf(aar+1:end-aar, aac+1:end-aac);
% r_lbp = 16;


% r_lbp - veľkosť oblasti pre výpočet LBP/C histogramu

f_s = p_klasifikacia(a, b, bf, s_min, pq, r_lbp, K, q_k, metoda_lbp, metoda_kontrast, kvant, R, N, lbp_q);
toc

%% orezanie výsledkov na pôvodnú veľkosť
% f_s = f_s(aar+1:end-aar, aac+1:end-aac);
% BW1 = BW1(aar+1:end-aar, aac+1:end-aac);
% BW2 = BW2(aar+1:end-aar, aac+1:end-aac);
%% vykreslenie výsledkov
figure
subplot(1,2,1)
imshow(f_s,[])
subplot(1,2,2)
imshow(a,[])

% figure
% subplot(1,2,1)
% imshow(f_s,[])
% subplot(1,2,2)
% imshow(segmenty,[])

figure
subplot(1,2,1)
imshow(BW1)
subplot(1,2,2)
imshow(BW2)

MIR_z(MIR_z==0) = [];
length(MIR_z);
xx = 1:length(MIR_z);
figure;
plot(xx, MIR_z);
axis([0 255 0 10]);
axis 'auto x'
