A= ReadList["greg.dat", {Number,Number,Number,Number}];

c1=1.0;                       (* distance O C1  in meters          *)
c2=1.0;                       (* distance O C2  in meters          *)
r1=0.5/500;                   (* relation picture width/width pixel*)
r2=0.5/500;                   (* relation picture width/width pixel*)
f1=60*Pi/180;                 (* angle of first camera             *)
f2=100*Pi/180;                (* angle of second camera            *)
v1={Cos[f1],Sin[f1],0};       (* vector from camera1 to object     *)
v2={Cos[f2],Sin[f2],0};       (* vector in direction of object     *)
v3={0,0,1};                   (* vector orthogonal (vcross(v1,v2)) *)
w1={-Sin[f1],Cos[f1],0};      (* vector in film direction          *)
w2={-Sin[f2],Cos[f2],0};      (* vector in film direction          *)
C1=-c1*v1;                    (* position of first camera          *)
C2=-c2*v2;                    (* position of second camera         *)
a=1;                          (* half of film width                *)
b=1;                          (* half of film height               *)

D1 =C1+a*w1-b*{0,0,1};        (* corner of first film              *)
DD1=C1+a*w1+b*{0,0,1};        (* corner of first film              *)
E1 =C1-a*w1-b*{0,0,1};        (* corner of first film              *)
EE1=C1-a*w1+b*{0,0,1};        (* corner of first film              *)

D2 =C2+a*w2-b*{0,0,1};        (* corner of second film             *)
DD2=C2+a*w2+b*{0,0,1};        (* corner of second film             *)
E2 =C2-a*w2-b*{0,0,1};        (* corner of second film             *)
EE2=C2-a*w2+b*{0,0,1};        (* corner of second film             *)


f[{xx1_,yy1_,xx2_,yy2_}]:=Module[{}, 
  x1=xx1*r1; y1=yy1*r2;
  x2=xx2*r1; y2=yy2*r2;
  A1=C1+x1*w1+y1*{0,0,1};       (* point on film1                    *)
  A2=C2+x2*w2+y2*{0,0,1};       (* point on film2                    *)
  B1=A1+2*c1*v1;                (* on light ray                      *)
  B2=A2+2*c2*v2;                (* on light ray                      *)
  A1x=A1[[1]]; A1y=A1[[2]]; A1z=A1[[3]];
  A2x=A2[[1]]; A2y=A2[[2]]; A2z=A2[[3]];
  v1x=v1[[1]]; v1y=v1[[2]]; v1z=v1[[3]];
  v2x=v2[[1]]; v2y=v2[[2]]; v2z=v2[[3]];
  v3x=v3[[1]]; v3y=v3[[2]]; v3z=v3[[3]];
  t1=(A1z*v2y*v3x-A2z*v2y*v3x-A1y*v2z*v3x+A2y*v2z*v3x-A1z*v2x*v3y+A2z*v2x*v3y+A1x*v2z*v3y-A2x*v2z*v3y+A1y*v2x*v3z -
      A2y*v2x*v3z-A1x*v2y*v3z+A2x*v2y*v3z)/(-(v1z*v2y*v3x)+v1y*v2z*v3x+v1z*v2x*v3y-v1x*v2z*v3y-v1y*v2x*v3z+v1x*v2y*v3z);
  t2=(A1z*v1y*v3x-A2z*v1y*v3x-A1y*v1z*v3x+A2y*v1z*v3x-A1z*v1x*v3y+A2z*v1x*v3y+A1x*v1z*v3y-A2x*v1z*v3y+A1y*v1x*v3z -
      A2y*v1x*v3z-A1x*v1y*v3z+A2x*v1y*v3z)/(-(v1z*v2y*v3x)+v1y*v2z*v3x+v1z*v2x*v3y-v1x*v2z*v3y-v1y*v2x*v3z+v1x*v2y*v3z);
  t3=(A1z*v1y*v2x-A2z*v1y*v2x-A1y*v1z*v2x+A2y*v1z*v2x-A1z*v1x*v2y+A2z*v1x*v2y+A1x*v1z*v2y-A2x*v1z*v2y+A1y*v1x*v2z -
      A2y*v1x*v2z-A1x*v1y*v2z+A2x*v1y*v2z)/(-(v1z*v2y*v3x)+v1y*v2z*v3x+v1z*v2x*v3y-v1x*v2z*v3y-v1y*v2x*v3z+v1x*v2y*v3z);
  PP={A1x+t1*v1x+t3*v3x/2,A1y+t1*v1y+t3*v3y/2,A1z+t1*v1z+t3*v3z/2}; PP
];

Do[ P=f[A[[k]]]; WriteString["3d.dat",P[[1]]," ",P[[2]]," ",P[[3]],"\n"],{k,Length[A]}]

