clear;
close all;

nb_images = 36; % Nombre d'images

% chargement des images
for i = 1:nb_images

    if i <= 10
        nom = sprintf('images/viff.00%d.ppm', i - 1);
    else
        nom = sprintf('images/viff.0%d.ppm', i - 1);
    end

    % L'ensemble des images de taille : nb_lignes x nb_colonnes x nb_canaux
    % x nb_images
    im(:, :, :, i) = imread(nom);
end

% chargement des points 2D suivis
% pts de taille nb_points x (2 x nb_images)
% sur chaque ligne de pts
% tous les appariements possibles pour un point 3D donne
% on affiche les coordonnees (xi,yi) de Pi dans les colonnes 2i-1 et 2i
% tout le reste vaut -1
pts = load('viff.xy');

% Chargement des matrices de projection
% Chaque P{i} contient la matrice de projection associee a l'image i
% RAPPEL : P{i} est de taille 3 x 4
load dino_Ps;

% chargement des masques (pour l'elimination des fonds bleus)
% de taille nb_lignes x nb_colonnes x nb_images
% A COMPLETER quand vous aurez termine la premiere partie permettant de
% binariser les images
load("mask.mat");

% cleaning des masks
im_mask = im_mask(:, :, :) == 0;
im_mask(1:8, :, :) = 0;
im_mask(end - 5:end, :, :) = 0;
im_mask(:, 1:5, :) = 0;
im_mask(:, end - 30:end, :) = 0;

%% Affichage de l'image que l'on souhaite segmenter

K = 100;
m = 0.1;
n = 3;
seuil_E = 10;
q_max = 5;
ind_img = 1;

figure(1);
imshow(im(:, :, :, ind_img)); title("Image " + num2str(ind_img));

%% Calculs des superpixels

figure(2);
imshow(im(:, :, :, ind_img)); title("Image " + num2str(ind_img));
hold on;
[germes, image_labelise, E] = super_pixel(im(:, :, :, ind_img), K, m, n, seuil_E, q_max);

%% Binarisation de l'image à partir des superpixels

figure(3);
R = germes(:, 3);
B = germes(:, 5);
scatter(R, B);
hold on

% définition des coefficients de la droite pour le seuillage
a = (0.7 - 0.1) / (0.6 - 0.2);
b = 0.1 - a * 0.2;

% affichage de la droite (pour vérification visuelle du seuillage)
plot([0 1], [b a + b]);

W = a * R + b - B > 0;
germes = [germes W];

%% Affichage et amélioration de la binarisation

figure(4);
tiledlayout(2, 2, 'Padding', 'none', 'TileSpacing', 'compact');

bw_img = reshape(germes(image_labelise, 6), size(im, 1), []);

% affichage de l'image binarisée
nexttile;
imshow(bw_img);

% calcul de la plus grande zone convexe blanche
CC = bwconncomp(bw_img, 4);
bw_img = zeros(size(bw_img));
bw_img(CC.PixelIdxList{1}) = 1;

% affichage de la plus grande zone convexe blanche
nexttile;
imshow(bw_img);

% calcul des zones convexes noires
bw_img = ~bw_img;

CC = bwconncomp(bw_img, 4);
bw_img = zeros(size(bw_img));
bw_img(CC.PixelIdxList{1}) = 1;

bw_img = ~bw_img;

% suppressions des zones convexes noires
nexttile;
imshow(bw_img);

% affichage du masque binaire fourni sur Moodle
bw_img = im_mask(:, :, ind_img);
nexttile;
imshow(bw_img);

%% Squeletisation de l'image

figure(5);
tiledlayout(2, 2, 'Padding', 'none', 'TileSpacing', 'compact');

% affichage de l'image
nexttile;
imshow(bw_img);
hold on

% affichage du contour de la binarisation
pixel_b = find(bw_img == 1);
[r, c] = ind2sub(size(bw_img), pixel_b(1));
contour = bwtraceboundary(bw_img, [r c], 'W', 8);
plot(contour(:, 2), contour(:, 1), 'g', 'LineWidth', 3);

% affichage du diagramme de voronoi du contour
T = 10; % échantillonage du contour matlab
[vx, vy] = voronoi(contour(1:T:end, 1), contour(1:T:end, 2));
plot(vy, vx, 'b');

% Selection des segments qui ont leurs extrémités dans l'image
ok = vx(1, :) > 1 & vx(1, :) < size(bw_img, 1) & ...
    vx(2, :) > 1 & vx(2, :) < size(bw_img, 1) & ...
    vy(1, :) > 1 & vy(1, :) < size(bw_img, 2) & ...
    vy(2, :) > 1 & vy(2, :) < size(bw_img, 2);
vx = floor(vx(:, ok));
vy = floor(vy(:, ok));

% affichage des points de voronoi visibles
nexttile;
imshow(bw_img);
hold on
plot(vy, vx, 'b');

% selection des segments avec les extremités dans la forme
ind1 = sub2ind(size(bw_img), vx(1, :), vy(1, :));
ok1 = bw_img(ind1) > 0;

ind2 = sub2ind(size(bw_img), vx(2, :), vy(2, :));
ok2 = bw_img(ind2) > 0;

ok = ok1 & ok2;
vx = vx(:, ok);
vy = vy(:, ok);

% affichage des points de voronoi uniquement dans la forme binaire
nexttile;
imshow(bw_img);
hold on

% mise en forme de vx et vy pour les prochains calculs
vx_ = vx';
vx_ = [vx_(:, 1); vx_(:, 2)];

vy_ = vy';
vy_ = [vy_(:, 1); vy_(:, 2)];

V_ = [vx_ vy_];

% calcul des rayons des points de voronoi au contour
contour_ = contour';
R = complex(V_(:, 1), V_(:, 2)) - complex(contour_(1, :), contour_(2, :));
R = abs(R);
R = min(R, [], 2);
R = R';

% affichage des cercles
vx_vy = [vy(1, :) vy(2, :); vx(1, :) vx(2, :)]';
T = 1;
viscircles(vx_vy(1:T:end, :), R(1:T:end), 'Color',[0 1 0 0.3]);

% affichage du squellette
plot(vy, vx, 'b');

% Scale Axis Transform
R_scaled = 1.05 * R;

dist = abs(complex(V_(:, 1), V_(:, 2)) - transpose(complex(V_(:, 1), V_(:, 2))));

R_vertical = ones(length(R_scaled), 1) * R_scaled;
R_horizontal = R_scaled' * ones(1, length(R_scaled));

[~, c] = ind2sub(size(dist), find(dist + R_vertical < R_horizontal));

vx(:, mod(c - 1, length(vx)) + 1) = [];
vy(:, mod(c - 1, length(vy)) + 1) = [];

R = [R(1:length(R) / 2); R(length(R) / 2 + 1:end)];
R(:, mod(c - 1, length(R)) + 1) = [];
R = [R(1, :) R(2, :)];

% affichage du squelette filtré par SAT
nexttile;
imshow(bw_img);
hold on
vx_vy = [vy(1, :) vy(2, :); vx(1, :) vx(2, :)]';
viscircles(vx_vy(1:T:end, :), R(1:T:end), 'Color',[0 1 0 0.3]); 
plot(vy, vx, 'b');

%% Reconstruction des points 3D

X = []; % Contient les coordonnees des points en 3D
color = []; % Contient la couleur associee
% Pour chaque coupple de points apparies
for i = 1:size(pts, 1)
    % Recuperation des ensembles de points apparies
    l = find(pts(i, 1:2:end) ~= -1);
    % Verification qu'il existe bien des points apparies dans cette image
    if size(l, 2) > 1 & max(l) - min(l) > 1 & max(l) - min(l) < 36
        A = [];
        R = 0;
        G = 0;
        B = 0;
        % Pour chaque point recupere, calcul des coordonnees en 3D
        for j = l
            A = [A; P{j}(1, :) - pts(i, (j - 1) * 2 + 1) * P{j}(3, :);
                P{j}(2, :) - pts(i, (j - 1) * 2 + 2) * P{j}(3, :)];
            R = R + double(im(int16(pts(i, (j - 1) * 2 + 1)), int16(pts(i, (j - 1) * 2 + 2)), 1, j));
            G = G + double(im(int16(pts(i, (j - 1) * 2 + 1)), int16(pts(i, (j - 1) * 2 + 2)), 2, j));
            B = B + double(im(int16(pts(i, (j - 1) * 2 + 1)), int16(pts(i, (j - 1) * 2 + 2)), 3, j));
        end;

        [U, S, V] = svd(A);
        X = [X V(:, end) / V(end, end)];
        color = [color [R / size(l, 2); G / size(l, 2); B / size(l, 2)]];
    end;

end;

fprintf('Calcul des points 3D termine : %d points trouves. \n', size(X, 2));

figure(6);
hold on;

for i = 1:size(X, 2)
    plot3(X(1, i), X(2, i), X(3, i), '.', 'col', color(:, i) / 255);
end;

axis equal;
view(80, -10);

%% Tetraedrisation de Delaunay

T = delaunayTriangulation(X(1, :)', X(2, :)', X(3, :)');
fprintf('Tetraedrisation terminee : %d tetraedres trouves. \n',size(T,1));

% Affichage de la tetraedrisation de Delaunay
% figure(7);
% tetramesh(T);

%% Filtrage de la tetraedrisation de Delaunay

% Calcul des barycentres de chacun des tetraedres
poids = [1 1 1 1] / 4;
nb_barycentres = size(T.ConnectivityList, 1);

for i = 1:size(T, 1)
    % Calcul des barycentres differents en fonction des poids differents
    % En commencant par le barycentre avec poids uniformes
    C_g(1:3, i) = poids * T.Points(T.ConnectivityList(i, :), :);
    C_g(4, i) = 1;
end

% Visualisation pour vérifier le bon calcul des barycentres
% for i = 1:nb_images
%    for k = 1:nb_barycentres
%        o = P{i}*C_g(:,k);
%        o = o./repmat(o(3,:),3,1);
%        imshow(im_mask(:,:,i));
%        hold on;
%        plot(o(2,:),o(1,:),'rx');
%        pause;
%        close;
%    end
% end

% Copie de la triangulation pour pouvoir supprimer des tetraedres
tri = T.ConnectivityList;
% Retrait des tetraedres dont au moins un des barycentres
% ne se trouvent pas dans au moins un des masques des images de travail
% Pour chaque barycentre
to_save = [];

for k = 1:nb_barycentres
    valide = 0;

    for i = 1:nb_images
        o = P{i} * C_g(:, k);
        o = o / o(3);
        x = floor(o(1));
        y = floor(o(2));

        if im_mask(x, y, i) == 0
            valide = 1;
            break
        end

    end

    if valide
        continue
    end

    to_save = [to_save k];
end

triBis = tri(to_save, :);
nb_barycentres = length(to_save);

fprintf('Retrait des tetraedres exterieurs a la forme 3D termine : %d tetraedres restants. \n', size(T, 1));

figure(8);
tetramesh(triBis, T.Points);
view(80, -10);

% Sauvegarde des donnees
save donnees;