%nahrejeme data korespodujicich bodu
M = load_corr ('panoram2.mat');

%vytvorime matici pro reseni soustavy
n=size(M,2);
for j=1:n

   %liche
   A(2*j-1,1)=M(1,j);
   A(2*j-1,2)=M(2,j);
   A(2*j-1,3)=1;
   A(2*j-1,4)=0;
   A(2*j-1,5)=0;
   A(2*j-1,6)=0;
   A(2*j-1,7)=-M(3,j)*M(1,j);
   A(2*j-1,8)=-M(3,j)*M(2,j);
   A(2*j-1,9)=-M(3,j);

   %sude	
   A(2*j,1)=0;
   A(2*j,2)=0;
   A(2*j,3)=0;
   A(2*j,4)=M(1,j);
   A(2*j,5)=M(2,j);
   A(2*j,6)=1;
   A(2*j,7)=-M(4,j)*M(1,j);
   A(2*j,8)=-M(4,j)*M(1,j);
   A(2*j,9)=-M(4,j);
end

%resime soustavu pomoci SVD
[U,S,V]=svd(A);
h=V(:,9);

%vysledek jeste upravime do matice 3x3
H=reshape(h,3,3)';

%a jeste matici homogenity zinvertujeme
Hi=inv(H);

%nacteme oba zdrojove obrazky
im1=imread('panorama1.jpg');
im2=imread('panorama2.jpg');

%zjistime velikosti
vel1=size(im1)
vel2=size(im2);

%pripravime vysledny obrazek
im=zeros(2*vel2(1),2*vel2(2),3);

%obrazek 2 nakopirujeme do vysledku
im(size(im2,1)/2:size(im2,1)*3/2-1, size(im2,2):2*size(im2,2)-1,:)=im2(:,:,:);

%vypocitame hranice pro transformovani obrazku
%musime je vydelit treti souradnici
edge1=H*[1,1,1]';
edge1=round(edge1./edge1(3))
edge2=H*[1,vel1(1),1]';
edge2=round(edge2./edge2(3))
edge3=H*[vel1(2),vel1(1),1]';
edge3=round(edge3./edge3(3))
edge4=H*[vel1(2),1,1]';
edge4=round(edge4./edge4(3))

for i=edge1(1):359
	for j=edge1(2):287
      
     	buff=Hi*[i;j;1];
		%zde provadime zaokrouhleni vysledku
		%inverzni transformace, lepsi by byla interpolace
		buffx = round(buff(1)/buff(3));
		buffy = round(buff(2)/buff(3));
		%kontrola jestli jsme v obraze
   	if ((buffx>=1) & (buffx<=size(im1,1)) & (buffy>=1) & (buffy<=size(im1,2)))
			im(j+144,i+360,:)=im1(buffy, buffx,:);
		end
	end
end

%zobrazime si krasny vysledek
im=im./255;
imshow(im,[]);
hold on;