function res = rekonstr(corr_file, file1, file2)

% load_corr()
load (corr_file,'-mat');

x              = CORR_EXCHANGE.d.corr(:,:,1)';
y              = CORR_EXCHANGE.d.corr(:,:,2)';
[m n]          = size(x);
M(1:2:2*m,1:n) = x;
M(2:2:2*m,1:n) = y;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

scale = 1000;

M = M/scale;
A = [M(3,:)'.*M(1,:)' M(4,:)'.*M(1,:)' M(1,:)' M(3,:)'.*M(2,:)' M(4,:)'.*M(2,:)' M(2,:)' M(3,:)' M(4,:)' ones(size(M,2),1)];

[U,S,V] = svd(A);
%rank(A)

% fundamentalni matice 
F = reshape(V(:,end), 3, 3)';
[U,S,V] = svd(F);
%rank(F)
S(end,end) = 0;
F = U*S*V';

A = [1/scale, 0, 0; 0, 1/scale, 0; 0, 0, 1];
[U,S,V] = svd(A);
%diag(S)

F = A'*F*A;
[U,S,V] = svd(F);

%rank(F)
% stredy kamer
e1 = U(:,end)./U(end,end)
e2 = V(:,end)./V(end,end)

U1 = [M(1:2,:);ones(1,size(M,2))];
U2 = [M(3:4,:);ones(1,size(M,2))];
sum(abs(diag(U1'*F*U2)));

M = M*scale;

% figure 1
im1 = imread(file1);
figure
imshow(im1);
axis auto;
hold on;

x1 = -1000;
x2 = 1000;
plot(M(1,:),M(2,:),'*m');
plot(e1(1),e1(2),'og');
for i = 1:size(M,2)
	P = [M(3,i) M(4,i) 1];
	line = F*P';
	y1 = -(line(3) + x1*line(1))/line(2);
	y2 = -(line(3) + x2*line(1))/line(2);
	plot([x1,x2],[y1,y2],'r');
end
hold off;

% figure 2
im2 = imread(file2);
figure
imshow(im2);
%axis([-200,900, -200,900]);
axis auto;
hold on;

x1 = -1000;
x2 = 1000;
plot(M(3,:),M(4,:),'*m');
plot(e2(1),e2(2),'og');
%hold on;
for i = 1:size(M,2)
	P = [M(1,i) M(2,i) 1];
	line = P*F;
	y1 = -(line(3) + x1*line(1))/line(2);
	y2 = -(line(3) + x2*line(1))/line(2);
	plot([x1,x2],[y1,y2],'r');
end

R = [M(1:2,:); ones(1,size(M,2)); M(3:4,:); ones(1,size(M,2))];
[U,S,V] = svd(R);
S(5:end,5:end) = 0;
P = U*S;
P1 = P(1:3,1:4)
P2 = P(4:6,1:4)

figure(1);

X = V(:,1:4)';
X = X./(ones(4,1)*X(end,:));
pts = P1*X;
pts = pts./(ones(3,1)*pts(end,:));
hold on;
for i=1:size(pts,2)
   plot(pts(1,i), pts(2,i), 'b*')
end

diff = pts - R(1:3,:);
diff1 = sum(sqrt(sum(diff.^2,1)))/size(diff,2)

figure(2);

X = V(:,1:4)';
X = X./(ones(4,1)*X(end,:));
pts = P2*X;
pts = pts./(ones(3,1)*pts(end,:));
hold on;
for i=1:size(pts,2)
   plot(pts(1,i), pts(2,i), 'b*')
end

diff = pts - R(4:6,:);
diff2 = sum(sqrt(sum(diff.^2,1)))/size(diff,2)
