M=load_corr('six.mat');

n=size(M,2); clear A;

for j=1:n

   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,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);

   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(2,j);

   A(2*j,9)=-M(4,j);

end;

[U,D,V]=svd(A);

h=V(:,end);

H=reshape(h,3,3)';

o1=double(imread('p1.jpg'));

o2=double(imread('p2.jpg'));

o1=o1./255;

o2=o2./255;

h=size(o1,1);

w=size(o1,2);

o=zeros(h,2*w,3);

o(:,w:2*w-1,:)=o1;

for i=1:1:h

   for j=1:1:2*w

      u=[j-w;i;1];

      v=H*u;

      x=v(1)/v(3);

      y=v(2)/v(3);

      if (x>=1) & (x<=w-10) & (y>=1) & (y<=h-10)

         % --- nejblizsi pixel ---

         %

         %ix=round(x); 

         %iy=round(y);

         %o(i,j,:)=o2(iy,ix,:);

         %

	 % --- bilinearni interpolace ---

         ix=floor(x); 

         iy=floor(y);

         fx=x-ix;

         fy=y-iy;

         v1=o2(iy,ix,:);

         v2=o2(iy,ix+1,:);

         v3=o2(iy+1,ix,:);

         v4=o2(iy+1,ix+1,:);

         v=v1+(v2-v1)*fx;

         o(i,j,:)=v+(v3+(v4-v3)*fx-v)*fy;

      end;

   end;

end;

imshow(o)

