% Skript pro demonstraci implementovaných tvarových deskriptorů na kontuře
% z reálných CT dat.
%
%                                                Vypracoval: Tomáš Janovič

%% Inicializace
close all;
clear all;
clc;

%% Načteni, extrakce kontury, sledování a výpočet těžiště
load img_seg.mat;
thresh = graythresh(img_seg); %otsu práh
img_bw = im2bw(img_seg,thresh); %prahováni
[L,num] = bwlabel(img_bw,4); %označení regionů
img_r = bwmorph(L==6,'remove'); %odstranění vnitřních pixelů
[row column] = find(img_r); %souřadnice pixelů kontury
[~,poz] = max(row); %pozice pro startovní pixel
start = [row(poz) column(poz)]; %souřadnice startovního pixelu
contour = bwtraceboundary(img_r,start,'E',8,Inf,'clockwise'); %sledování kontury
contour = contour'; %transponování matice
centroid = regionprops(img_r,'centroid'); %výpočet těžiště
centroid = centroid.Centroid; %zisk těžště ze struktury

figure(1)
subplot(221)
imshow(img_seg)
title('Originálni obraz');
subplot(222)
imshow(img_r)
title('Extrakce kontury');
subplot(223)
plot(contour(2,:),contour(1,:))
hold on
plot(centroid(1),centroid(2),'-green*','LineWidth',5);
title('Sledování kontury');
set(gca,'Ydir', 'reverse')
xlabel('x');
ylabel('y');

%% Funkce vzdálenosti od těžiště
xosa = 0:length(contour)-1; %x osa
dist = sqrt(((contour(2,:)-centroid(1))').^2 + ((contour(1,:)-centroid(2))').^2); %Euklidovská vzdálenost mezi body a těžištěm

figure(2)
plot(xosa,dist)
title('Funkce vzdálenosti od těžiště');
xlabel('Body kontury podél obvodu');
ylabel('r(n) [px]');

%% Kumulativní úhlová funkce
X = contour(2,:); %x souřadnice
Y = contour(1,:); %y souřadnice
M = size(X,2);

okno = 10; %okno pro výpočet diferencí
A = zeros(1,M);
for i=1:M
    x1 = 0;
    x2 = 0;
    y1 = 0;
    y2 = 0;
    for j=1:okno
        pa = i-j;
        pb = i+j;
        if pa<1
            pa = M+pa;
        end
        if pb>M
            pb = pb-M;
        end
        x1 = x1+X(pa);
        y1 = y1+Y(pa);
        x2 = x2+X(pb);
        y2 = y2+Y(pb);
    end
    x1 = x1/okno;
    x2 = x2/okno;
    y1 = y1/okno;
    y2 = y2/okno;
    dx = x2-x1;
    dy = y2-y1;
    if dx==0
        dx = 1e-5;
    end
    if dx>0 && dy>0
        A(i) = atan(dy/dx);
    elseif dx>0 && dy<0
        A(i) = atan(dy/dx)+2*pi;
    else
        A(i) = atan(dy/dx)+pi;
    end
end

G = zeros(1,M);
G(1) = 0;
for i=2:M
    d = min(abs(A(i)-A(i-1)), abs(abs(A(i)-A(i-1))-2*pi));
    if d>0.5
        G(i) = G(i-1);
    elseif (A(i)-A(i-1))<-pi
        G(i) = G(i-1)-(A(i)-A(i-1)+2*pi);
    elseif (A(i)-A(i-1))>pi
        G(i) = G(i-1)-(A(i)-A(i-1)-2*pi);
    else
        G(i) = G(i-1)-(A(i)-A(i-1));
    end
end

figure(3)
subplot(121)
plot(xosa,A)
title('Úhlová funkce');
xlabel('Body kontury podél obvodu');
ylabel('\theta [rad]');
subplot(122)
plot(xosa,G)
title('Kumulatní úhlová funkce');
xlabel('Body kontury podél obvodu');
ylabel('$$\varphi$$ [rad]','Interpreter','latex')

%% Funkce popisující směrnici tečny
phi = atan2(contour(1,:),contour(2,:)); %výpočet tangens úhlu
phi = unwrap(phi); %rozbalení fáze

figure(4)
plot(xosa,phi)
xlabel('Body kontury podél obvodu');
ylabel('\theta [rad]');
title('Funkce popisující směrnici tečny');

%% Fourierovy deskriptory (FDs)
len_G = length(G); %délka kumulativní úhlové funkce
pocet_FDs = [8 16 32];
figure(5)
for i=1:length(pocet_FDs)
    FD_count = pocet_FDs(i); %počet FDs (doporučeno 16-64)
    spektrum = fft(G,len_G); %spektrum funkce
    FDs = spektrum(1:FD_count); %FDs
    rep_spectrum = [FDs zeros(1,len_G-2*FD_count+1) conj(fliplr(FDs(2:end)))]; %spektrum kumulativní úhlové funkce tvaru dle FDs
    rep = real(ifft(rep_spectrum)); %rekonstruovaná kumulativní úhlová funkce dle FDs
    subplot(1,3,i)
    plot(xosa,G)
    hold on
    plot(xosa,rep,'red')
    xlabel('Body kontury podél obvodu');
    ylabel('$$\varphi$$ [rad]','Interpreter','latex');
    title(['Rekonstruovaná funkce dle ' num2str(pocet_FDs(i)) ' FDs']);
    legend('Originálni funkce','Rekonstruovaná funkce dle FDs');
end

%% Vlnkové deskriptory (WDs)
delka = round(length(G)/20); %délka prodlouženych částí
poc = G(1)*ones(1,delka); %opakování první hodnoty funkce
kon = G(end)*ones(1,delka); %opakování poslední hodnoty funkce
fun = [poc G kon]; %prodloužení funkce kvůli okrajovým jevům

CWT = cwt(fun,15,'bior1.5'); %spojitá vlnková transformace
CWT = CWT(delka+1:end-delka); %odstranění prodloužených částí

figure(6)
subplot(121)
plot(xosa,G)
title('Kumulativní úhlová funkce');
xlabel('Body kontury podél obvodu');
ylabel('$$\varphi$$ [rad]','Interpreter','latex')
subplot(122)
plot(xosa,CWT)
title('Koeficienty CWT');
xlabel('Body kontury podél obvodu');

