5 files changed, 146 insertions(+), 32 deletions(-)

M CMakeLists.txt
A eval.sh
M evaluate.cpp
A find-good-models.sh
A size.cpp

@@ 22,8 22,11 @@ target_link_libraries(surcut igl::core igl::opengl igl::opengl_glfw igl::opengl_
 target_compile_options(surcut PRIVATE -g -O2 -Wall -Werror -Wpedantic -DDEBUG)
 
 add_executable(evaluate ${SRCFILES} evaluate.cpp)
-target_link_libraries(evaluate igl::core igl::opengl igl::opengl_glfw igl::opengl_glfw_imgui
-  igl::tetgen OpenMP::OpenMP_CXX igl_stb_image)
+target_link_libraries(evaluate igl::core igl::tetgen OpenMP::OpenMP_CXX igl_stb_image)
 target_compile_options(evaluate PRIVATE -g -O2 -Wall -Werror -Wpedantic -DDEBUG)
 
+add_executable(size mesh.cpp size.cpp)
+target_link_libraries(size igl::core igl::tetgen)
+target_compile_options(size PRIVATE -g -O2 -Wall -Werror -Wpedantic -DDEBUG)
+
 igl_download_tutorial_data()

@@ 0,0 1,7 @@
+#!/bin/bash
+
+mkdir eval-outputs 2>/dev/null
+
+for mesh in $(find meshes/small -type f); do
+  build/evaluate $mesh -i32768 | xargs -I'{}' printf "%s: %s\n" $mesh '{}' > eval-outputs/${mesh##*/}.output
+done

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


@@ 23,8 23,7 @@ Vector3d face_center(Mesh &mesh, int face) {
   return (x + y + z) / 3.0;
 }
 
-template<typename T>
-T avg(std::vector<T> v, T t) {
+template <typename T> T avg(std::vector<T> v, T t) {
   for (T e : v) {
     t += e;
   }


@@ 32,9 31,8 @@ T avg(std::vector<T> v, T t) {
   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) {
+double score_mesh(Mesh &mesh, int iterations, int range) {
   Eigen::VectorXd face_areas;
   {
     igl::doublearea(mesh.vertices, mesh.faces, face_areas);


@@ 46,17 44,17 @@ double score_mesh(Mesh &mesh, int iterations) {
 
   /* Even consider adaptive sampling? */
   std::mt19937 gen(0);
-  std::uniform_real_distribution<double> rand_zero_2pi(0, M_PI);
-  std::uniform_real_distribution<double> rand_unit(-1, 1);
+  std::uniform_int_distribution<int> rand_positive(0, range);
+  std::uniform_int_distribution<int> rand_centered(-range / 2, range / 2);
 
   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);
+    double z = ((double) rand_centered(gen)) * 2.0 / (double) range;
+    double theta = ((double) rand_positive(gen)) / (double) range * M_PI;
+    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;


@@ 68,48 66,64 @@ double score_mesh(Mesh &mesh, int iterations) {
   Vector3d best_plane_point;
 
 #pragma omp parallel for
-  for (int i = 0; i < iterations; i++) {
+  for (int i = 0; i < iterations; i++) { /* O(I) */
     double score = 0;
 
     Vector3d dir = random_directions.row(i);
+
     Vector3d plane_point;
     { /* Find the point on which the cut plane is */
-      std::vector<double> dots;
       std::vector<Vector3d> projecteds;
-      for (int f = 0; f < mesh.faces.rows(); f++) {
+      std::vector<double> dots;
+      for (int f = 0; f < mesh.faces.rows(); f++) { /* O(F) */
         Vector3d normal = mesh.normals.row(f);
-        if (normal.dot(dir) > 0.0) {
+        if (normal.dot(dir) >= 0.0) {
           Vector3d center = face_center(mesh, f);
           double dot = center.dot(dir);
           Vector3d projected = dir * dot;
+          dots.push_back(dot);
           projecteds.push_back(center - projected);
         }
       }
-      plane_point = avg(projecteds, Vector3d(0.0, 0.0, 0.0));
+      double avg_dot = avg(dots, 0.0);
+      plane_point = avg(projecteds, Vector3d(0.0, 0.0, 0.0)) + avg_dot * dir;
     }
 
     double moldable_area = 0.0;
-    for (int f = 0; f < mesh.faces.rows(); f++) {
+    for (int f = 0; f < mesh.faces.rows(); f++) { /* O(F^2) */
+      Vector3d normal = mesh.normals.row(f);
       Vector3d center = face_center(mesh, f);
+
       Ray ray;
-      ray.src = center;
+      // Move 1% away from the face to avoid intersecting with triangle boundaries
+      ray.src = center + normal * 0.01 * sqrt(face_areas(f));
       if ((center - plane_point).dot(dir) >= 0.0) {
         ray.dir = dir;
       } else {
         ray.dir = -dir;
       }
 
-      // 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;
+      double dot = normal.dot(ray.dir);
+      if (dot >= -0.0001) {
+        // It can be molded. But is it obstructed? We approximate this by shooting a ray from the
+        // triangle center.
+        for (int other_f = 0; other_f < mesh.faces.rows(); other_f++) { /* O(F) */
+          if (other_f == f) continue;
+          Vector3d other_normal = mesh.normals.row(other_f);
+          // If we're going along its face, don't check ;) 
+          if (abs(ray.dir.dot(other_normal)) < 0.1) {
+            continue;
+          }
+          Vector3d _intersection_point;
+          if (intersect3D_RayTriangle(ray, other_normal, other_f, mesh.vertices, mesh.faces, 
+                &_intersection_point) == 1) {
+            goto obstructed;
+          }
         }
-      }
-      if (!hit) 
         moldable_area += face_areas(f);
+obstructed:
+        (void)0;
+      }
     }
     score += moldable_area / total_area;
 


@@ 118,18 132,21 @@ double score_mesh(Mesh &mesh, int iterations) {
       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("  %4.2lf  ", score);
+    printf("[%7.2lf, %7.2lf, %7.2lf]  ", plane_point(0), plane_point(1), plane_point(2));
+    printf("[%7.2lf, %7.2lf, %7.2lf]\n", dir(0), dir(1), dir(2));
   }
   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),
+  printf("plane point: [%7.2lf, %7.2lf, %7.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));
+  printf("  direction: [%7.2lf, %7.2lf, %7.2lf]\n", best_dir(0), best_dir(1), best_dir(2));
 
   return best_score;
 }
 
 int main(int argc, char *argv[]) {
   int iterations = 100;
+  int range = 100;
   char *filename = nullptr;
   for (int i = 1; i < argc; i++) {
     char *arg = argv[i];


@@ 154,6 171,25 @@ int main(int argc, char *argv[]) {
           }
         }
       } break;
+      case 'r': {
+        if (arg[2] != 0) {
+          int res = sscanf(arg + 2, "%d", &range);
+          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", &range);
+          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;


@@ 201,7 237,7 @@ int main(int argc, char *argv[]) {
   arap::precomputation(mesh, arap_data);
   compute_things(mesh);
 
-  score_mesh(mesh, iterations);
+  score_mesh(mesh, iterations, range);
 
   return 0;
 }

@@ 0,0 1,8 @@
+#!/bin/bash
+
+MODEL_DIR="./libigl/tutorial/data"
+MODELS=$(find "$MODEL_DIR" -type f -iname "*.obj")
+
+for model in $MODELS; do
+  build/size -t $model && build/evaluate $model -i1
+done

@@ 0,0 1,60 @@
+#include <igl/readOBJ.h>
+#include <Eigen/Core>
+
+#include "mesh.h"
+
+int main(int argc, char *argv[]) {
+  bool mesh = false;
+  char *filename = nullptr;
+
+  for (int i = 1; i < argc; i++) {
+    char *arg = argv[i];
+    if (arg[0] == '-') {
+      switch (arg[1]) {
+      case 't': {
+        mesh = true;
+      } break;
+      default: {
+        fprintf(stderr, "Unknown argument '-%c'\n", arg[1]);
+        return 1;
+      }
+      }
+    } else {
+      filename = arg;
+    }
+  }
+
+  if (filename == nullptr) {
+    fprintf(stderr, "Missing input file!\n");
+    return 1;
+  }
+
+  if (mesh) {
+    TetmeshOut tetmesh_out;
+    Eigen::MatrixXd read_verts;
+    Eigen::MatrixXi read_faces;
+    igl::readOBJ(filename, 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;
+    }
+    printf("file: %s\n", filename);
+    printf("vertices: %6zu\n", tetmesh_out.vertices.rows());
+    printf("   faces: %6zu\n", tetmesh_out.faces.rows());
+    printf("    tets: %6zu\n", tetmesh_out.tets.rows());
+  } else {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    igl::readOBJ(filename, V, F);
+    printf("file: %s\n", filename);
+    printf("vertices: %6zu\n", V.rows());
+    printf("   faces: %6zu\n", F.rows());
+  }
+  return 0;
+}