N=8;

ux(1)=0.561198; uy(1)=0.418403;

vx(1)=0.423177; vy(1)=0.394097;

ux(2)=0.277344; uy(2)=0.428819;

vx(2)=0.174479; vy(2)=0.567708;

ux(3)=0.500000; uy(3)=0.824653;

vx(3)=0.531250; vy(3)=0.781250;

ux(4)=0.680990; uy(4)=0.763889;

vx(4)=0.641927; vy(4)=0.623264;

ux(5)=0.815104; uy(5)=0.343750;

vx(5)=0.549479; vy(5)=0.208333;

ux(6)=0.503906; uy(6)=0.312500;

vx(6)=0.319010; vy(6)=0.331597;

ux(7)=0.161458; uy(7)=0.696181;

vx(7)=0.161458; vy(7)=0.894097;

ux(8)=0.614583; uy(8)=0.137153;

vx(8)=0.376302; vy(8)=0.112847;

for i=1:N

 R(1,i)=ux(i);

 R(2,i)=uy(i);

 R(3,i)=1;

 R(4,i)=vx(i);

 R(5,i)=vy(i);

 R(6,i)=1;

end

A=zeros(N,9);

for i=1:N

 A(i,:)=[ux(i)*vx(i) ux(i)*vy(i) ux(i) uy(i)*vx(i) uy(i)*vy(i) uy(i) vx(i) vy(i) 1];

end;

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

FF=reshape(V(:,end),3,3)';

[S,D,V]=svd(FF);

D(3,3)=0;

F=S*D*V';F=FF;

e1=S(:,3)./S(3,3);

e2=V(:,3)./V(3,3);

%

o1=imread('k1.jpg'); figure(1); 

hold off; imshow(o1); hold on;

os=size(o1); w=os(2); h=os(1);

x1=e1(1); y1=e1(2); x2=1; e1=e1.*[w h 1]'

%axis([x1*w w y1*h h]);

axis([0 w 0 h]);

plot(e1(1),e1(2),'r*');

for i=1:N

 plot(ux(i)*w,uy(i)*h,'go');

 %plot(ux(i)*w,uy(i)*h,'r*');

 l1=F*[vx(i) vy(i) 1]';

 y1=(l1(3)+x1*l1(1))/l1(2);

 y2=(l1(3)+x2*l1(1))/l1(2);

 line([x1*w x2*w]',[-y1*h -y2*h]');

end

%

o2=imread('k2.jpg'); figure(2);

hold off; imshow(o2); hold on;

os=size(o2); w=os(2); h=os(1);

x1=e2(1); y1=e2(2); x2=1; e2=e2.*[w h 1]'

%axis([-5000 w 0 900]);

axis([0 w 0 h]);

plot(e2(1),e2(2),'r*');

for i=1:N

 plot(vx(i)*w,vy(i)*h,'go');

 %plot(vx(i)*w,vy(i)*h,'r*');

 l2=[ux(i) uy(i) 1]*F; 

 y1=(l2(3)+x1*l2(1))/l2(2);

 y2=(l2(3)+x2*l2(1))/l2(2);

 line([x1*w x2*w]',[-y1*h -y2*h]');

end

[S,D,V]=svd(R);

D(5:end,5:end)=0; 

P=S*D;P=P(:,1:4);

X=V';X=X(1:4,:);

X=X./(ones(4,1)*X(4,:))

C1=P(1:3,:);

C2=P(4:6,:);

PP1=diag([0 0 0 1]);

PP2=diag([0 0 0 1]);

P1=diag([1 1 1 1]);

PP1(1:3,:)=P(1:3,:);

PP2(1:3,:)=P(4:6,:);

C1=P(1:3,:);

C2=P(4:6,:);

H=inv(PP1)*P1;

P2=PP2*H;

AA=P2(1:3,1:3);

Tb=-inv(AA)*P2(1:3,4)

S=zeros(3,3);

S(2,1)=Tb(3); S(3,1)=-Tb(2); S(3,2)=Tb(1);

S(1,2)=-Tb(3); S(1,3)=Tb(2); S(2,3)=-Tb(1);

FN=S*inv(AA)

B1=C1*X;

B2=C2*X;

B1=B1./(ones(3,1)*B1(3,:));

B2=B2./(ones(3,1)*B2(3,:));

figure(1); plot(B1(1,:)*w,B1(2,:)*h,'ro');

figure(2); plot(B2(1,:)*w,B2(2,:)*h,'ro');

for i=1:N

 err1(i)=sqrt((w*(ux(i)-B1(1,i)))^2+(h*(uy(i)-B1(2,i)))^2);

end

for i=1:N

 err2(i)=sqrt((w*(vx(i)-B2(1,i)))^2+(h*(vy(i)-B2(2,i)))^2);

end

ser1=sum(err1,2)/N

ser2=sum(err2,2)/N

