4 files changed, 94 insertions(+), 76 deletions(-)

M evaluate.cpp
M main.cpp
M surcut.cpp
M surcut.h

@@ 13,7 13,7 @@ using Eigen::AngleAxisd;
 using Eigen::Matrix3d;
 using Eigen::Vector3d;
 
-double score_mesh(Mesh &);
+double score_mesh(Mesh &, int);
 Vector3d face_center(Mesh &, int);
 
 Vector3d face_center(Mesh &mesh, int face) {


@@ 23,52 23,69 @@ Vector3d face_center(Mesh &mesh, int face) {
   return (x + y + z) / 3.0;
 }
 
-// Give a score to the mesh, to estimate how good it is for molding with stitching
-double score_mesh(Mesh &mesh) {
-  const int num_samples = 1000;
+template<typename T>
+T avg(std::vector<T> v, T t) {
+  for (T e : v) {
+    t += e;
+  }
+  t /= v.size();
+  return t;
+}
 
+
+// Give a score to the mesh, to estimate how good it is for molding with stitching
+double score_mesh(Mesh &mesh, int iterations) {
   Eigen::VectorXd face_areas;
   {
     igl::doublearea(mesh.vertices, mesh.faces, face_areas);
     face_areas /= 2.0;
   }
+  double total_area = 0.0;
+  for (int i = 0; i < face_areas.rows(); i++)
+    total_area += face_areas(i);
 
   /* Even consider adaptive sampling? */
-  std::random_device rd;
-  std::mt19937 gen(rd());
-  std::uniform_real_distribution<double> random(-M_PI, M_PI);
+  std::mt19937 gen(0);
+  std::uniform_real_distribution<double> rand_zero_2pi(0, M_PI);
+  std::uniform_real_distribution<double> rand_unit(-1, 1);
+
+  MatrixXd random_directions(iterations, 3);
+  // Sample directions deterministically (since we're using OpenMP in the next `for` loop)
+  for (int i = 0; i < iterations; i++) {
+    // See:
+    // https://math.stackexchange.com/questions/44689/how-to-find-a-random-axis-or-unit-vector-in-3d
+    double theta = rand_zero_2pi(gen);
+    double z = rand_unit(gen);
+    double s = sqrt(1.0 - z*z);
+    random_directions(i, 0) = s * cos(theta);
+    random_directions(i, 1) = s * sin(theta);
+    random_directions(i, 2) = z;
+    random_directions.row(i) = random_directions.row(i).normalized();
+  }
 
   double best_score = 0.0;
   Vector3d best_dir;
   Vector3d best_plane_point;
 
 #pragma omp parallel for
-  for (int i = 0; i < num_samples; i++) {
+  for (int i = 0; i < iterations; i++) {
     double score = 0;
 
-    Vector3d dir;
-    { /* sample direction */
-      double y_theta = random(gen);
-      double z_theta = random(gen);
-      Matrix3d rotation;
-      rotation = AngleAxisd(y_theta, Vector3d::UnitY()) * AngleAxisd(z_theta, Vector3d::UnitZ());
-      dir = rotation * Vector3d::UnitX();
-    }
-
+    Vector3d dir = random_directions.row(i);
     Vector3d plane_point;
     { /* Find the point on which the cut plane is */
-      std::vector<double> projected;
+      std::vector<double> dots;
+      std::vector<Vector3d> projecteds;
       for (int f = 0; f < mesh.faces.rows(); f++) {
         Vector3d normal = mesh.normals.row(f);
         if (normal.dot(dir) > 0.0) {
-          projected.push_back(face_center(mesh, f).dot(dir));
+          Vector3d center = face_center(mesh, f);
+          double dot = center.dot(dir);
+          Vector3d projected = dir * dot;
+          projecteds.push_back(center - projected);
         }
       }
-      double avg = 0.0;
-      for (double d : projected)
-        avg += d;
-      avg /= projected.size();
-      plane_point = avg * dir;
+      plane_point = avg(projecteds, Vector3d(0.0, 0.0, 0.0));
     }
 
     double moldable_area = 0.0;


@@ 81,37 98,77 @@ double score_mesh(Mesh &mesh) {
       } else {
         ray.dir = -dir;
       }
-      if (!ray_intersect_from_point(ray, mesh, nullptr)) {
-        moldable_area += face_areas(f);
+
+      // No intersections roughly means it's moldable in this direction
+      bool hit = false;
+      for (int f = 0; f < mesh.faces.rows(); f++) {
+        Vector3d _intersection_point;
+        if (intersect3D_RayTriangle(ray, mesh.normals.row(f), f, mesh.vertices, mesh.faces, &_intersection_point) == 1) {
+          hit = true;
+          break;
+        }
       }
+      if (!hit) 
+        moldable_area += face_areas(f);
     }
-    score += moldable_area;
+    score += moldable_area / total_area;
 
     if (score > best_score) {
       best_score = score;
       best_dir = dir;
       best_plane_point = plane_point;
     }
+    printf("  %4.2lf [%5.2lf, %5.2lf, %5.2lf]\n", score, dir(0), dir(1), dir(2));
   }
-  printf("best score: %lf\n", best_score);
-  printf("best direction: [%4.2lf, %4.2lf, %4.2lf]\n", best_dir(0), best_dir(1), best_dir(2));
-  printf("best plane point: [%4.2lf, %4.2lf, %4.2lf]\n", best_plane_point(0), best_plane_point(1),
+  printf("      score:  %4.2lf\n", best_score);
+  printf("plane point: [%4.2lf, %4.2lf, %4.2lf]\n", best_plane_point(0), best_plane_point(1),
          best_plane_point(2));
+  printf("  direction: [%4.2lf, %4.2lf, %4.2lf]\n", best_dir(0), best_dir(1), best_dir(2));
 
   return best_score;
 }
 
 int main(int argc, char *argv[]) {
-  if (argc != 2) {
-    fprintf(stderr, "Usage: ./evaluate <path-to-mesh>\n");
-    return 1;
+  int iterations = 100;
+  char *filename = nullptr;
+  for (int i = 1; i < argc; i++) {
+    char *arg = argv[i];
+    if (arg[0] == '-') {
+      switch (arg[1]) {
+      case 'i': {
+        if (arg[2] != 0) {
+          int res = sscanf(arg + 2, "%d", &iterations);
+          if (res < 1) {
+            fprintf(stderr, "failed to read number '%s'\n", arg + 2);
+            return 1;
+          }
+        } else {
+          if (i + 1 == argc) {
+            fprintf(stderr, "Ran out of arguments! Missing number to -i\n");
+          }
+          char *arg = argv[++i];
+          int res = sscanf(arg, "%d", &iterations);
+          if (res < 1) {
+            fprintf(stderr, "failed to read number '%s'\n", arg);
+            return 1;
+          }
+        }
+      } break;
+      default: {
+        fprintf(stderr, "Unknown argument '-%c'\n", arg[1]);
+        return 1;
+      }
+      }
+    } else {
+      filename = arg;
+    }
   }
 
   TetmeshOut tetmesh_out;
   {
     MatrixXd read_verts;
     MatrixXi read_faces;
-    igl::readOBJ(argv[1], read_verts, read_faces);
+    igl::readOBJ(filename, read_verts, read_faces);
 
     TetmeshIn tetmesh_in;
     tetmesh_in.vertices = &read_verts;


@@ 144,8 201,7 @@ int main(int argc, char *argv[]) {
   arap::precomputation(mesh, arap_data);
   compute_things(mesh);
 
-  double score = score_mesh(mesh);
-  printf("%lf\n", score);
+  score_mesh(mesh, iterations);
 
   return 0;
 }

@@ 336,24 336,16 @@ int main(int argc, char *argv[]) {
     draw_plane(viewer, mesh);
     if (run_with_update) {
       if (step_arap(mesh) == StepReturn::Converged) {
-        static int last_iter_update = -99;
-        // If we have converged with rest state updates, stop.
-        if (iteration_count - last_iter_update == 1) {
-          // run_with_update = false;
-        }
-        fprintf(stderr, "Rest state update (diff=%d)\n", iteration_count - last_iter_update);
         mesh.rest_state = mesh.vertices;
         arap::precomputation(mesh, arap_data);
         compute_things(mesh);
-        update_viewer(mesh);
-        last_iter_update = iteration_count;
       }
     } else {
       if (run_arap) {
         step_arap(mesh);
-        update_viewer(mesh);
       }
     }
+    update_viewer(mesh);
     return false;
   };
 

@@ 66,34 66,6 @@ int ray_intersect(Ray ray, int src_vertex, Mesh &mesh) {
   return best_vertex;
 }
 
-// Checks whether the Ray intersects any faces, and stores the intersection point in `I` if it
-// does (and if `I != nullptr`)
-bool ray_intersect_from_point(Ray ray, Mesh &mesh, Vector3d *I) {
-  Vector3d closest_ip;
-  double closest_dist2 = std::numeric_limits<double>::max();
-  bool did_hit = false;
-  for (int f = 0; f < mesh.faces.rows(); f++) {
-    Vector3d point;
-    int res =
-        intersect3D_RayTriangle(ray, mesh.normals.row(f), f, mesh.vertices, mesh.faces, &point);
-    if (res == 1) {
-      double dist = (ray.src - point).squaredNorm();
-      if (dist > 0.0 && dist < closest_dist2) {
-        closest_dist2 = dist;
-        closest_ip = point;
-        did_hit = true;
-      }
-    }
-  }
-
-  if (!did_hit)
-    return false;
-
-  if (I)
-    *I = closest_ip;
-  return true;
-}
-
 // intersect3D_RayTriangle(): find the 3D intersection of a ray with a triangle
 //    Input:  a ray R, and a triangle T
 //    Output: *I = intersection point (when it exists)


@@ 127,8 99,7 @@ int intersect3D_RayTriangle(Ray ray, Vector3d triangle_normal, int triangle, Mat
     return 0;  // => no intersect
   // for a segment, also test if (r > 1.0) => no intersect
 
-  if (I)
-    *I = ray.src + r * ray.dir; // intersect point of ray and plane
+  *I = ray.src + r * ray.dir; // intersect point of ray and plane
 
   // is I inside T?
   double uu, uv, vv, wu, wv, D;

@@ 56,7 56,6 @@ enum Enum {
 }
 
 int ray_intersect(Ray, int, Mesh &);
-bool ray_intersect_from_point(Ray, Mesh &, Vector3d *);
 int intersect3D_RayTriangle(Ray, Vector3d, int, MatrixXd &, MatrixXi &, Vector3d *);
 StepReturn::Enum step_arap(Mesh &, double = 1.0);
 void compute_things(Mesh &mesh, double = 1.0, bool only_plane = false);