% 3. Homography computation
function H = Homography(p1, p2)
    % Number of points
    num_points = size(p1, 1);
 
    % Build the matrix A for which A * h = 0
    A = zeros(2 * num_points, 9);
    for i = 1:num_points
        x = p1(i, 1);
        y = p1(i, 2);
        xp = p2(i, 1);
        yp = p2(i, 2);
 
        A(2 * i - 1, :) = [-x, -y, -1, 0, 0, 0, x * xp, y * xp, xp];
        A(2 * i, :)     = [0, 0, 0, -x, -y, -1, x * yp, y * yp, yp];
    end
 
    % Perform SVD on A
    [U, D, V] = svd(A);
    disp('The U is:')
    disp(U)

    disp('The D is:')
    disp(D)

    disp('The V is:')
    disp(V)

    % The homography matrix is the last column of V, reshaped to 3x3
    H = reshape(V(:, 9), [3, 3])';

    disp('The H is:')
    disp(H)
end

%Utility function 

function i = geokor (H, f, g)
%        ====================
	[h1 w1 d] = size(f); f = double(f);                     % Image dimensions
	[h2 w2 d] = size(g); g = double(g);

         % Transformation of image corners to determine the resulting size
	cp = norm2(H * [1 1 w1 w1; 1 h1 1 h1; 1 1 1 1]);
	Xpr = min([cp(1, :) 0]) : max([cp(1, :) w2]);              % min x : max x
	Ypr = min([cp(2, :) 0]) : max([cp(2, :) h2]);              % min y : max y
	[Xp, Yp] = ndgrid(Xpr, Ypr);
	[wp hp] = size(Xp);                                           % = size(Yp)
	X = norm2(H \ [Xp(:) Yp(:) ones(wp*hp, 1)]');    % Indirect transformation
	
	xI = reshape(X(1, :), wp, hp)';           % Resampling of intensity values 
	yI = reshape(X(2, :), wp, hp)';
	f2 = zeros(hp, wp); i = f2;
	for k = 1 : d
		f2(1:h1, 1:w1, k) = f(:, :, k);
		i(:, :, k) = interp2(f2(:, :, k), xI, yI);    % bilinear interpolation           
	end
                      % Copy the original image g in the rectified image f
	off = -round([min([cp(1, :) 0]) min([cp(2, :) 0])]);
	Index = find(sum(g, 3) > 9); 
	for k = 1 : d
    		iPart = i(1+off(2) : h2+off(2), 1+off(1) : w2+off(1), k);
    		fChannel = g(:, :, k);
    		iPart(Index) = fChannel(Index);
    		i(1+off(2) : h2+off(2), 1+off(1) : w2+off(1), k) = iPart;
	end
	i(~isfinite(i)) = 0;
	i = uint8(i);

end

function n = norm2(x)
	for i = 1 : 3
		n(i,:) = x(i,:) ./ x(3,:);	
	end

end


% 1. Image acquisition:
%C:\Users\DELL\Desktop\Assignment2

jakobs_a = imread('C:\Users\DELL\Desktop\Assignment2\jakobs_a.jpg');
jakobs_b = imread('C:\Users\DELL\Desktop\Assignment2\jakobs_b.jpg');
jakobs_c = imread('C:\Users\DELL\Desktop\Assignment2\jakobs_c.jpg');

 
%  2. Correspondence analysis: Transfer the three images into the computer (imread) and 
% measure interactively (figure, imshow, ginput) at least four corresponding image points 
% x->x between two neighboring images.

figure;
subplot(1, 3, 1); imshow(jakobs_a); title('jakobs_a.jpg');
subplot(1, 3, 2); imshow(jakobs_b); title('jakobs_b.jpg');
subplot(1, 3, 3); imshow(jakobs_c); title('jakobs_c.jpg');

% Select 4 corresponding points between Image 1 and Image 2
figure; imshow(jakobs_a); title('Select 4 points in jakobs_a');
[x1, y1] = ginput(4);
 
figure; imshow(jakobs_b); title('Select 4 points in jakobs_b');
[x2, y2] = ginput(4);
 
% Store points as homogeneous coordinates
points1 = [x1 y1 ones(4, 1)];
points2 = [x2 y2 ones(4, 1)];



 
%  Compute the homography matrix H12 from first to second image.
H12 = Homography(points1, points2);
 
 

i = geokor(H12, jakobs_a, jakobs_b);     % Transform image f using H and combine with g
% figure; 

imshow(i);title('Select 4 points in attached image 1 and 2');
[x3, y3] = ginput(4);
imshow(img3); title('Select 4 points in attached image 3');
[x4, y4] = ginput(4);
points12 = [x3 y3 ones(4, 1)];
points4 = [x4 y4 ones(4, 1)];
H23 = Homography(points12, points4);

finalResult = geokor(H23, i, jakobs_c); 
imshow(finalResult);