7 files changed, 202 insertions(+), 58 deletions(-)

M arap-mold.h
M evaluate.cpp
M main.cpp
M mesh.cpp
M mesh.h
M surcut.cpp
M surcut.h

@@ 1,9 1,9 @@
 #pragma once
 
-#include <igl/ARAPEnergyType.h>
-#include <igl/min_quad_with_fixed.h>
 #include <Eigen/Core>
 #include <Eigen/Sparse>
+#include <igl/ARAPEnergyType.h>
+#include <igl/min_quad_with_fixed.h>
 
 #include "mesh.h"
 


@@ 82,7 82,6 @@ struct Data {
   Eigen::Vector3d *mht_plane_n;
   // Map face indices to tet indices. At least one tet, at most two.
   Eigen::MatrixXi mht_tf2t;
-  std::vector<Eigen::Quaternion<double>> mht_tet_rotations;
   // XXX(debug): Register rotation angles for all faces
   std::vector<double> mht_face_mold_angle;
   std::vector<double> mht_face_arap_angle;

@@ 1,48 1,151 @@
-#include <stdio.h>
 #include <Eigen/Core>
+#include <Eigen/Geometry>
+#include <igl/doublearea.h>
+#include <igl/readOBJ.h>
+
+#include <random>
+#include <stdio.h>
 
 #include "mesh.h"
+#include "surcut.h"
 
+using Eigen::AngleAxisd;
+using Eigen::Matrix3d;
 using Eigen::Vector3d;
-using Eigen::Vector4d;
-using Eigen::Matrix4d;
 
 double score_mesh(Mesh &);
+Vector3d face_center(Mesh &, int);
 
-// Give a score to the mesh, to estimate how good it is for molding with stitching 
+Vector3d face_center(Mesh &mesh, int face) {
+  Vector3d x = mesh.vertices.row(mesh.faces(face, 0));
+  Vector3d y = mesh.vertices.row(mesh.faces(face, 1));
+  Vector3d z = mesh.vertices.row(mesh.faces(face, 2));
+  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) {
-  // This is the rough plan:
-  //     d = sample-direction()
-  //     S = orthographic-project(M, d)
-  //     <somehow find cut plane>
-  //     a = sum area of non-moldable faces for d 
-  //     score = area(S) - a
+  const int num_samples = 1000;
+
+  Eigen::VectorXd face_areas;
+  {
+    igl::doublearea(mesh.vertices, mesh.faces, face_areas);
+    face_areas /= 2.0;
+  }
 
-  const int num_samples = 100;
+  /* Even consider adaptive sampling? */
+  std::random_device rd;
+  std::mt19937 gen(rd());
+  std::uniform_real_distribution<double> random(-M_PI, M_PI);
 
   double best_score = 0.0;
-  /* Even consider adaptive sampling? */ 
+  Vector3d best_dir;
+  Vector3d best_plane_point;
+
+#pragma omp parallel for
   for (int i = 0; i < num_samples; i++) {
-    double score;
-    Vector3d dir /* = sample dir */;
-    Matrix4d projection /* = make this matrix */;
+    double score = 0;
 
-    for (int f = 0; f < mesh.faces.rows(); f++) {
-      Vector4d normal = mesh.normals.row(f);
-      Vector4d transformed_normal = projection * normal;
+    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();
+    }
 
-      (void)dir;
-      (void)transformed_normal;
+    Vector3d plane_point;
+    { /* Find the point on which the cut plane is */
+      std::vector<double> projected;
+      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));
+        }
+      }
+      double avg = 0.0;
+      for (double d : projected)
+        avg += d;
+      avg /= projected.size();
+      plane_point = avg * dir;
+    }
+
+    double moldable_area = 0.0;
+    for (int f = 0; f < mesh.faces.rows(); f++) {
+      Vector3d center = face_center(mesh, f);
+      Ray ray;
+      ray.src = center;
+      if ((center - plane_point).dot(dir) >= 0.0) {
+        ray.dir = dir;
+      } else {
+        ray.dir = -dir;
+      }
+      if (!ray_intersect_from_point(ray, mesh, nullptr)) {
+        moldable_area += face_areas(f);
+      }
     }
+    score += moldable_area;
 
-    best_score = std::max(best_score, score);
+    if (score > best_score) {
+      best_score = score;
+      best_dir = dir;
+      best_plane_point = plane_point;
+    }
   }
+  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),
+         best_plane_point(2));
 
   return best_score;
 }
 
-
 int main(int argc, char *argv[]) {
-  printf("Hello, world!\n");
+  if (argc != 2) {
+    fprintf(stderr, "Usage: ./evaluate <path-to-mesh>\n");
+    return 1;
+  }
+
+  TetmeshOut tetmesh_out;
+  {
+    MatrixXd read_verts;
+    MatrixXi read_faces;
+    igl::readOBJ(argv[1], read_verts, read_faces);
+
+    TetmeshIn tetmesh_in;
+    tetmesh_in.vertices = &read_verts;
+    tetmesh_in.faces = &read_faces;
+
+    int error = make_tetmesh(tetmesh_in, tetmesh_out);
+    if (error) {
+
+      fprintf(stderr, "tetmesh returned error: %d\n", error);
+      return 1;
+    }
+  }
+
+  Mesh mesh(tetmesh_out);
+  mesh.boundary = VectorXi(4);
+  mesh.boundary << 0, 1, 2, 3;
+
+  arap_data.faces = &mesh.faces;
+  arap_data.normals = &mesh.normals;
+  arap_data.mht_tf2t = mesh.face2tets;
+  arap_data.mht_plane_p = &plane_p;
+  arap_data.mht_plane_n = &plane_n;
+  arap_data.vert_is_above = &is_point_above;
+  arap_data.moldable = &moldable;
+  arap_data.orientation = &orientation;
+  arap_data.U0 = mesh.vertices;
+  arap_data.with_dynamics = true;
+  // XXX(todo): what's a reasonable value here?
+  arap_data.h = 0.01;
+  arap::precomputation(mesh, arap_data);
+  compute_things(mesh);
+
+  double score = score_mesh(mesh);
+  printf("%lf\n", score);
+
   return 0;
 }

@@ 11,8 11,8 @@
 
 #include "arap-mold.h"
 #include "mesh.h"
-#include "viridis.h"
 #include "surcut.h"
+#include "viridis.h"
 
 using igl::opengl::glfw::Viewer;
 


@@ 132,7 132,7 @@ void draw_menu(Mesh &mesh) {
   // discard very accute triangles, or something.
   // if (ImGui::Button("Remesh")) {
   //   VectorXd field = VectorXd::Zero(mesh.vertices.rows());
-  //   
+  //
   //   MatrixXd new_vertices;
   //   MatrixXi new_faces;
   //   VectorXd new_field;


@@ 296,31 296,19 @@ int main(int argc, char *argv[]) {
     }
   }
 
-  Mesh mesh;
-
-  mesh.vertices = tetmesh_out.vertices;
-  mesh.faces = tetmesh_out.faces;
-  mesh.tets = tetmesh_out.tets;
-  mesh.face2tets = tetmesh_out.face2tets;
-  mesh.tet2faces = tetmesh_out.tet2faces;
-  mesh.neighbors = tetmesh_out.neighbors;
-
-  mesh.rest_state = mesh.vertices;
-  mesh.forces = MatrixXd::Zero(mesh.vertices.rows(), mesh.vertices.cols());
+  Mesh mesh(tetmesh_out);
   mesh.boundary = VectorXi(4);
   mesh.boundary << 0, 1, 2, 3;
 
   mesh_print_sizes(mesh);
 
-  forces_endpoints_ptr = &forces_endpoints; /* this is out of place*/ 
+  forces_endpoints_ptr = &forces_endpoints; /* this is out of place */
 
   arap_data.faces = &mesh.faces;
   arap_data.normals = &mesh.normals;
   arap_data.mht_tf2t = mesh.face2tets;
   arap_data.mht_plane_p = &plane_p;
   arap_data.mht_plane_n = &plane_n;
-  arap_data.mht_tet_rotations = std::vector<Eigen::Quaternion<double>>(
-      mesh.tets.rows(), Eigen::Quaternion<double>(1, 0, 0, 0));
   arap_data.vert_is_above = &is_point_above;
   arap_data.moldable = &moldable;
   arap_data.orientation = &orientation;


@@ 347,7 335,7 @@ int main(int argc, char *argv[]) {
   viewer.callback_pre_draw = [&mesh](Viewer &viewer) {
     draw_plane(viewer, mesh);
     if (run_with_update) {
-      if (step_arap(mesh)) {
+      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) {

@@ 48,7 48,7 @@ int make_tetmesh(TetmeshIn &in, TetmeshOut &out) {
 
   try {
     tetgenbehavior tb;
-    char *switches = (char *)"pq1.414Ynn";
+    char *switches = (char *)"Qpq1.414Ynn";
     bool ok = tb.parse_commandline(switches);
     if (!ok) {
       cerr << "[ERR]: from `tetgen` in " << __FUNCTION__ << ":" << __LINE__


@@ 172,3 172,14 @@ void mesh_print_sizes(Mesh &mesh) {
   Dbg(mesh.boundary);
   Dbg(mesh.forces);
 }
+
+Mesh::Mesh(TetmeshOut &tetmesh_out) {
+  vertices = tetmesh_out.vertices;
+  faces = tetmesh_out.faces;
+  tets = tetmesh_out.tets;
+  face2tets = tetmesh_out.face2tets;
+  tet2faces = tetmesh_out.tet2faces;
+  neighbors = tetmesh_out.neighbors;
+  rest_state = vertices;
+  forces = MatrixXd::Zero(vertices.rows(), vertices.cols());
+}

@@ 51,6 51,8 @@ struct Mesh {
   MatrixXd rest_state;
   VectorXi boundary;
   MatrixXd forces;
+
+  Mesh(TetmeshOut &);
 };
 
 void mesh_print_sizes(Mesh &);

@@ 1,8 1,8 @@
 #include <igl/per_face_normals.h>
 
-#include "surcut.h"
-#include "mesh.h"
 #include "arap-mold.h"
+#include "mesh.h"
+#include "surcut.h"
 
 std::vector<bool> is_point_above;
 std::vector<arap::Moldable::Enum> moldable;


@@ 22,8 22,8 @@ arap::Args arap_args;
 Eigen::Vector3d plane_p(0, 0, 0);
 Eigen::Vector3d plane_n(0, 1, 0);
 
-// Shoot a ray from `p` in direction `v`, and see if it intersects with any of the `faces`. If it
-// does, return the index of the closest vertex to the intersection point. If not, return `-1`.
+// Shoot a `ray` and see if it intersects any of the faces in `mesh`. If it does, return the index
+// of the closest vertex to the intersection point. If not, return `-1`.
 int ray_intersect(Ray ray, int src_vertex, Mesh &mesh) {
   // XXX(todo): Option to intersect with cut plane. Maybe this isn't a good idea - rays from any
   //            non-moldable face will intersect with the plane!


@@ 66,6 66,34 @@ 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)


@@ 99,7 127,8 @@ 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
 
-  *I = ray.src + r * ray.dir; // intersect point of ray and plane
+  if (I)
+    *I = ray.src + r * ray.dir; // intersect point of ray and plane
 
   // is I inside T?
   double uu, uv, vv, wu, wv, D;


@@ 125,7 154,7 @@ int intersect3D_RayTriangle(Ray ray, Vector3d triangle_normal, int triangle, Mat
 
 // Perform one ARAP step. Return `true` if the difference between the previous and current state is
 // negligible.
-bool step_arap(Mesh &mesh, double scale) {
+StepReturn::Enum step_arap(Mesh &mesh, double scale) {
   iteration_count++;
   Eigen::MatrixXd bc(mesh.boundary.rows(), 3);
   for (int i = 0; i < mesh.boundary.rows(); i++)


@@ 136,12 165,14 @@ bool step_arap(Mesh &mesh, double scale) {
   arap::iter(mesh, arap_args, bc, arap_data, mesh.vertices);
   if (mesh.vertices != mesh.vertices) { /* Contains NaN */
     mesh.vertices = prev_vertices;
-    return true;
+    return StepReturn::WouldNaN;
   }
   compute_things(mesh);
 
   double max_diff = (mesh.vertices - prev_vertices).maxCoeff();
-  return max_diff < scale * 5e-4; /* This factor was found experimentally! */
+  if (max_diff < scale * 5e-4)
+    return StepReturn::Converged;
+  return StepReturn::NotDone;
 }
 
 // Compute all the things that we need each time the vertex set changes. Here's what we compute:


@@ 253,7 284,7 @@ void compute_things(Mesh &mesh, double scale, bool only_plane) {
         auto force_direction = other_vertex - vertex;
         if (force_direction.norm() > scale * force_dist_cutoff)
           continue;
-        
+
         // XXX(todo): what should the magnitude of the force be? Should probably take a closer look
         // at the dynamics implementation in `arap-mold.cpp`.
 

@@ 3,12 3,12 @@
 #include <Eigen/Core>
 #include <Eigen/Sparse>
 
-#include "mesh.h"
 #include "arap-mold.h"
+#include "mesh.h"
 
-using Eigen::Vector3d;
 using Eigen::MatrixXd;
 using Eigen::MatrixXi;
+using Eigen::Vector3d;
 
 // Various computed things about the mesh
 // If a point is above the plane.


@@ 32,7 32,8 @@ extern int iteration_count;
 // Whether to stitch together unmoldable geometry
 extern bool stitching;
 
-// Endpoints for all forces calculated. If not `nullptr` then this is written to in `compute_things`.
+// Endpoints for all forces calculated. If not `nullptr` then this is written to in
+// `compute_things`.
 extern MatrixXd *forces_endpoints_ptr;
 
 extern arap::Data arap_data;


@@ 46,7 47,16 @@ struct Ray {
   Vector3d dir;
 };
 
+namespace StepReturn {
+enum Enum {
+  NotDone,
+  Converged,
+  WouldNaN,
+};
+}
+
 int ray_intersect(Ray, int, Mesh &);
-int intersect3D_RayTriangle(Ray , Vector3d , int , MatrixXd &, MatrixXi &, Vector3d *);
-bool step_arap(Mesh &, double = 1.0);
+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);