function tv(M, W, C);
  printf("Trucco & Verri Method\n");
  r=1:3;
  % scale the matrix
  gamma = sqrt( dot(M(3,r),M(3,r)) );
  M = M ./ gamma;
  M = sign(M(3,4))*M; % force Tz > 0
  ox = dot(M(1,r),M(3,r));
  oy = dot(M(2,r),M(3,r));
  printf("  (Ox,Oy): %f %f\n",ox,oy);
  fx = sqrt( dot(M(1,r),M(1,r)) - ox^2 );
  fy = sqrt( dot(M(2,r),M(2,r)) - oy^2 );
  T = zeros(3,1);
  T(3) = -M(3,4);
  T(1) = (-M(1,4) - ox*T(3))/fx;
  T(2) = (-M(2,4) - oy*T(3))/fy;
  R = zeros(3);
  R(1,:) = (M(1,r) - ox*M(3,r)) ./fx;
  R(2,:) = (M(2,r) - oy*M(3,r)) ./fy;
  R(3,:) = M(3,r);
  R
  xy = fx/fy;
  %radial([fy T(3) 0], W, C, R, ox, oy, T(2) );
  %radial([fy T(3) -3.9842e-07], W, C, R, ox, oy, T(2) );
  k = fmins( 'radial', [fy T(3) 0], [], 0, W, C, R, ox, oy, T(2) );
  fx = xy*k(1);
  fy = k(1);
  T(3) = k(2);
  T
  printf("  (fx,fy): %f %f\n",fx,fy);
  printf("  (k1): %f\n\n",k(3));
%  fmins( 'radial', [fy T(3) 0], [], 0, W, C, R, ox, oy, T(2) )