% 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] = mosaiqueter2(homos, imgs)

    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;

    % 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 I2 transformee par l'homographie H.
    Imos = [];
    n_bloc = 4;
    bloc_size = floor(nblImos*nbcImos/n_bloc);
    for i = 1:n_bloc
       fprintf('%d/%d\r', i, n_bloc);
        [y, x] = ind2sub([nblImos, nbcImos], i*bloc_size+1:bloc_size*(i+1));
        y = y';
        x = x';
        
        distances = zeros(bloc_size, length(imgs));

        % 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;

        % 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(bloc_size, 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], [], 2);
            dk(dk < 0) = 0;

            distances(:,k) = dk;
        end

        poids = distances ./ sum(distances, 2);
        poids(isnan(poids)) = 0;
        poids(:) = max(poids(:), 0);

        % 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(bloc_size, length(imgs), 3);
        for k = 1:length(imgs)
            imgk = im2double(imgs{k});
            [nblIk, nbcIk, ~] = size(imgk);
            imgk = reshape(imgk, nblIk*nbcIk, 3);
            
            valid_pixels = (xy_Rs(:,k,1) >= 1) & (xy_Rs(:,k,1) <= nbcIk) & (xy_Rs(:,k,2) >= 1) & (xy_Rs(:,k,2) <= nblIk);
            ts(valid_pixels, k, :) = poids(valid_pixels,k) .* imgk(sub2ind([nblIk, nbcIk], xy_Rs(valid_pixels,k,2), xy_Rs(valid_pixels,k,1)), :);
        end
        Imos = cat(1, Imos, sum(ts, 2));
    end
    Imos = cat(1, Imos, zeros(nblImos*nbcImos - size(Imos, 1), 1, 3));
    Imos = reshape(Imos, nblImos, nbcImos, 3);
    fprintf('\n');
end