[Trifocal+P2Pt] normalization of all line params not just normals

rfabbri · Oct 29, 2023 · 4f9d56c · 4f9d56c
1 parent 7a5d40f
commit 4f9d56c
Show file tree

Hide file tree

Showing 2 changed files with 15 additions and 14 deletions.
diff --git a/minus/chicago14a.hxx b/minus/chicago14a.hxx
@@ -7060,18 +7060,17 @@ lines2params(const F plines[pp::nvislines][io::ncoords2d_h], C<F> * __restrict p
     const F *l0 = plines[triple_intersections[l][0]];
     const F *l1 = plines[triple_intersections[l][1]];
     const F *l2 = plines[triple_intersections[l][2]];
-    double l0l0 = vec::dot(l0,l0), l0l1 = vec::dot(l0,l1), l1l1 = vec::dot(l1,l1),
-    l2l0 = vec::dot(l2,l0), l2l1 = vec::dot(l2,l1);
+    double l0l1 = vec::dot(l0,l1), l2l0 = vec::dot(l2,l0), l2l1 = vec::dot(l2,l1);
     // cross([l0l0 l1l0 l2l0], [l0l1 l1l1 l2l1], l2_l0l1);
     double l2_l0l1[3]; 
     {
       F v1[3], v2[3];
-      v1[0] = l0l0; v1[1] = l0l1; v1[2] = l2l0;
-      v2[0] = l0l1; v2[1] = l1l1; v2[2] = l2l1;
+      v1[0] = 1.; v1[1] = l0l1; v1[2] = l2l0;
+      v2[0] = l0l1; v2[1] = 1.; v2[2] = l2l1;
       vec::cross(v1, v2, l2_l0l1);
     }
     params_lines[l][0] = l2_l0l1[0]/l2_l0l1[2]; // divide by the last coord (see cross prod formula, plug direct)
-    params_lines[l][1] = l2_l0l1[1]/l2_l0l1[2];
+    params_lines[l][1] = l2_l0l1[1]/l2_l0l1[2]; // TODO: normalize this to unit/ instead
   }
   //        
   //    pChart: just unit rands 17x1
@@ -7143,7 +7142,7 @@ gammify(C<F> * __restrict params /*[ chicago: M::nparams]*/)
   static unsigned constexpr triple_intersect[6][2] = {{0,3},{0+1,3+1},{0+2,3+2},{0,6},{0+1,6+1},{0+2,6+2}};
 
   // diag1 --> pTriple in params -----------------------------------------------
-  unsigned i = 9*3;
+  unsigned i = 9*3; // TODO: move to unsigned char
   for (unsigned tl=0; tl < 6; ++tl) {  // for each tangent line
     params[i++] *= std::conj(gammas[triple_intersect[tl][0]]);
     params[i++] *= std::conj(gammas[triple_intersect[tl][1]]);
@@ -7210,20 +7209,20 @@ point_tangents2lines(const F p[pp::nviews][pp::npoints][io::ncoords2d], const F
   assert (i0 < i1 && i1 < 3);
   unsigned i2 = (i0 == 0) ? ((i1 == 1) ? 2 : 1) : 0;
 
-  static constexpr double eps = 1e-4;
+  static constexpr double eps = 1e-5;
 
   if (v::area2(p[0][i0],p[0][i1],p[0][i2])  < eps || 
       v::area2(p[1][i0],p[1][i1],p[1][i2])  < eps || 
       v::area2(p[2][i0],p[2][i1],p[2][i2])  < eps) {
-//    std::cerr << "MINUS: area error ----XXXXXXXXXXXXX-------\n";
+//    std::cerr << "MINUS: area error ------------------------\n";
 //    std::cerr << "Areas: " << 
 //      v::area2(p[0][i0],p[0][i1],p[0][i2]) << " "  << 
 //      v::area2(p[1][i0],p[1][i1],p[1][i2]) << " " << 
 //      v::area2(p[2][i0],p[2][i1],p[2][i2]) << std::endl;
     return false;
   }
-
 
+  // TODO: require inpunt points be on the view-sphere (unit-normalized), not h-normalized
   vec::cross2(p[0][i0], p[0][i1], plines[0]);
   vec::cross2(p[1][i0], p[1][i1], plines[1]);
   vec::cross2(p[2][i0], p[2][i1], plines[2]);
@@ -7245,7 +7244,6 @@ point_tangents2lines(const F p[pp::nviews][pp::npoints][io::ncoords2d], const F
   minus_3d<F>::point_tangent2line(p[0][i1], t[0][i1], plines[12]);
   minus_3d<F>::point_tangent2line(p[1][i1], t[1][i1], plines[13]);
   minus_3d<F>::point_tangent2line(p[2][i1], t[2][i1], plines[14]);
-  // TODO: test normalize to unit vectors for better numerics
 
   if (v::abs_angle_between_lines(plines[0], plines[9])  < eps || 
       v::abs_angle_between_lines(plines[1], plines[10]) < eps ||
@@ -7262,7 +7260,7 @@ point_tangents2lines(const F p[pp::nviews][pp::npoints][io::ncoords2d], const F
       v::abs_angle_between_lines(plines[0], plines[12]) < eps ||
       v::abs_angle_between_lines(plines[1], plines[13]) < eps ||
       v::abs_angle_between_lines(plines[2], plines[14]) < eps) {
-//    std::cerr << "MINUS: angle error ----AAAAAAAAAAAAA-------\n";
+//    std::cerr << "MINUS: angle error ------------------------\n";
 //    std::cerr << "Angles: " << 
 //          v::abs_angle_between_lines(plines[0], plines[9])  << " " << 
 //          v::abs_angle_between_lines(plines[1], plines[10]) << " " <<

diff --git a/minus/minus.h b/minus/minus.h
@@ -214,8 +214,11 @@ struct minus_io_common {
   // Functions -----------------------------------------------------------------
   static void invert_intrinsics(const F K[/*3 or 2 ignoring last line*/][ncoords2d_h], const double pix_coords[][ncoords2d], double normalized_coords[][ncoords2d], unsigned npts);
   static void invert_intrinsics_tgt(const F K[/*3 or 2 ignoring last line*/][ncoords2d_h], const double pix_tgt_coords[][ncoords2d], double normalized_tgt_coords[][ncoords2d], unsigned npts);
-  static void normalize_line(F l[ncoords2d_h]) {
-    const F nrm = std::hypot(l[0], l[1]);
+  static void normalize_line(F l[ncoords2d_h]) { // we were making it unit
+                                                 // normal, but the scale of the
+                                                 // whole line is more important
+                                                 // for certain numerics
+    const F nrm = l[0]*l[0] + l[1]*l[1] + l[2]*l[2];
     l[0] /= nrm; l[1] /= nrm; l[2] /= nrm;
   }
   static void normalize_lines(F lines[][ncoords2d_h], unsigned nlines);
@@ -247,7 +250,7 @@ struct minus_io_14a : public minus_io_common<F> {
     typedef minus_util<F> u;
     // camera 0 (2nd camera relative to 1st)
     u::quat2rotm(rs, (F *) cameras[0]);
-    F n = sqrt(rs[8]*rs[8] + rs[9]*rs[9] + rs[10]*rs[10]) + sqrt(rs[11]*rs[11] + rs[12]*rs[12] + rs[13]*rs[13]);
+    const F n = sqrt(rs[8]*rs[8] + rs[9]*rs[9] + rs[10]*rs[10]) + sqrt(rs[11]*rs[11] + rs[12]*rs[12] + rs[13]*rs[13]);
     cameras[0][3][0] = rs[8]/n;
     cameras[0][3][1] = rs[9]/n;
     cameras[0][3][2] = rs[10]/n;