% Calcul d'une mosaique d'image à partir de 2 images : I1 et I2
% connaissant l'homographie H entre les 2.
% L'image resultat est stockee dans Res.
% On choisit de projeter I2 dans I1 pour construire la mosaique.
% Attention !!!
% On suppose un axe de rotation parallèle aux colonnes.
% C'est la raison pour laquelle on inverse les lignes et les colonnes
% dans la reconstruction de la mosaique.

function [Imos] = mosaiqueter(homos, imgs)

    % On recupere la taille des deux images.
    Iref = im2double(imgs{1});
    [nblIref, nbcIref, nc] = size(Iref);

    xmin = Inf;
    ymin = Inf;
    xmax = -Inf;
    ymax = -Inf;

    for k=1:length(imgs)
        [nblIk, nbcIk, ~] = size(imgs{k});
        
        % On calcule l'homographie inverse, normalisee,
        % pour nous permettre d'effectuer la transformation de I2 vers I1.
        Hinv = inv(homos{k});
        Hinv = Hinv ./ Hinv(3, 3);

        % On calcule les coordonnees des 4 coins de I2 dans I1.
        xy_coinsIk_Rk = [1 1; nbcIk 1; nbcIk nblIk; 1 nblIk];
        
        % Application de l'homographie Hinv sur ces coins.
        % Calcul des images des coins dans I1.
        xy_coinsIk_Rref = appliquerHomographie(Hinv, xy_coinsIk_Rk);

        % Determination des dimensions de l'image mosaique,
        % les xmin ymin xmax ymax, ou :
        %  - xmin represente la plus petite abscisse parmi les abscisses des images
        %    des coins de I2 projetes dans I1 et les coins dans I1,
        %  - etc
        % Lignes et colonnes sont inversees.
        xmin = min([xy_coinsIk_Rref(:, 1)' xmin]);
        ymin = min([xy_coinsIk_Rref(:, 2)' ymin]);
        xmax = max([xy_coinsIk_Rref(:, 1)' xmax]);
        ymax = max([xy_coinsIk_Rref(:, 2)' ymax]); 
    end

    % On arrondit de maniere a etre certain d'avoir les coordonnees initiales
    % bien comprises dans l'image.
    xmin = floor(xmin);
    ymin = floor(ymin);
    xmax = ceil(xmax);
    ymax = ceil(ymax);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CONTRUCTION DE LA MOSAIQUE %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Calcul de la taille de la mosaique.
    % Lignes et colonnes sont inversees.
    nblImos = ymax - ymin + 1;
    nbcImos = xmax - xmin + 1;

    Imos = zeros(nblImos, nbcImos, nc);

    % Calcul de l'origine de l'image I1 dans le repere de la mosaique Imos.
    O1x_Rmos = 1 - (xmin - 1);
    O1y_Rmos = 1 - (ymin - 1);

    % Copie de l'image I1.
    % Lignes et colonnes sont inversees.
    Imos(O1y_Rmos:O1y_Rmos + nblIref - 1, O1x_Rmos:O1x_Rmos + nbcIref - 1, :) = Iref;

    % Copie de l'image I2 transformee par l'homographie H.
    for x = 1:nbcImos
        fprintf('%d/%d\r', x, nbcImos);
        for y = 1:nblImos
            % Calcul des coordonnees dans I1 connaissant les coordonnees du point origine de I1 dans Imos.
            y_Rref = y - O1y_Rmos;
            x_Rref = x - O1x_Rmos;
            
            distances = zeros(1, length(imgs));
        
            dref = min([nbcIref - x_Rref, nblIref - y_Rref, x_Rref, y_Rref]);
            if dref <= 0 || all(Iref(y_Rref, x_Rref, :) == 0)
                dref = 0;
            end
            distances(1) = dref;

            % Dans le repere attache a l'image I1,
            % nous estimons les coordonnees du point image de (y_R1,x_R1)
            % par l'homographie H : (xy_R2).
            xy_Rs = zeros(length(imgs), 2);
            xy_Rs(1,:) = [x_Rref y_Rref];
            for k = 1:length(imgs)
                
                xy_Rk = appliquerHomographie(homos{k}, [x_Rref y_Rref]);
                xy_Rs(k,:) = round(xy_Rk);
                % Il existe plusieurs strategies, mais, ici,
                % pour estimer les coordonnees (entieres) ,
                % on choisit : sans interpolation, le plus proche voisin.
                x_Rk = round(xy_Rk(1));
                y_Rk = round(xy_Rk(2));
                [nblIk, nbcIk, ~] = size(imgs{k});
                dk = min([nbcIk - x_Rk, nblIk - y_Rk, x_Rk, y_Rk]);

                if dk <= 0 || all(imgs{k}(y_Rk, x_Rk, :) == 0)
                    dk = 0;
                end

                distances(k) = dk;
            end
            
            poids = zeros(1, length(imgs));
            if sum(distances) ~= 0
                poids = distances / sum(distances);
                poids = max(poids, 0);
            end

            % On verifie que xy_R2 appartient bien a l'image I2
            % avant d'affecter cette valeur a Imos
            % Lignes et colonnes sont inversees.
            ts = zeros(length(imgs), 3);
            for k = 1:length(imgs)
                if poids(k) ~= 0
                    ts(k, :) = poids(k) * double(imgs{k}(xy_Rs(k,2), xy_Rs(k,1), :)) / 255;
                else
                    ts(k, :) = [0 0 0];
                end
            end
            Imos(y, x, :) = sum(ts, 1);
            
        end
        
    end
    fprintf('\n');

end