~mht/surcut

0dbc35cd48d7238ffa3ae1b1cd3effe59f237622 — Martin Hafskjold Thoresen 1 year, 2 months ago 16669cc
Only optionally do the force cutoff thing

The scaling is still very off - bunny worked fine with the previous
cutoff, but fertility did not.
7 files changed, 64 insertions(+), 28 deletions(-)

M CMakeLists.txt
M README.md
M evaluate.cpp
M find-good-models.sh
M main.cpp
M surcut.cpp
M surcut.h
M CMakeLists.txt => CMakeLists.txt +4 -0
@@ 29,4 29,8 @@ 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)

add_executable(2obj 2obj.cpp)
target_link_libraries(2obj igl::core)
target_compile_options(2obj PRIVATE -g -O2 -Wall -Werror -Wpedantic -DDEBUG)

igl_download_tutorial_data()

M README.md => README.md +20 -0
@@ 10,3 10,23 @@

 - Handle intersections when having forces


## Scripts

Convert all files from one format to `obj`:

```bash
OUT_DIR="./meshes"
FMT="ply"
DIR="./libigl/tutorial/data"
MODELS=$(find "$DIR" -type f -iname "*.$FMT")
for model in $MODELS; do
  f=${model##*/}
  build/2obj $model "$OUT_DIR"/${f/$FMT/obj}
done
```

Get scores from the output of `find-good-models.sh`:
```
awk '/file/ {printf "%45s", $0} /score/{print $0}'
```

M evaluate.cpp => evaluate.cpp +11 -9
@@ 73,20 73,22 @@ double score_mesh(Mesh &mesh, int iterations, int range) {

    Vector3d plane_point;
    { /* Find the point on which the cut plane is */
      std::vector<Vector3d> projecteds;
      std::vector<double> dots;
      Vector3d moldables(0, 0, 0), non_moldables(0, 0, 0);
      int n_moldables = 0, n_non_moldables = 0;
      for (int f = 0; f < mesh.faces.rows(); f++) { /* O(F) */
        Vector3d normal = mesh.normals.row(f);
        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);
          moldables += face_center(mesh, f);
          n_moldables++;
        } else {
          non_moldables += face_center(mesh, f);
          n_non_moldables++;
        }
      }
      double avg_dot = avg(dots, 0.0);
      plane_point = avg(projecteds, Vector3d(0.0, 0.0, 0.0)) + avg_dot * dir;
      moldables /= n_moldables;
      non_moldables /= n_non_moldables;

      plane_point = (moldables + non_moldables) / 2.0;
    }

    double moldable_area = 0.0;

M find-good-models.sh => find-good-models.sh +7 -2
@@ 1,8 1,13 @@
#!/bin/bash

MODEL_DIR="./libigl/tutorial/data"
MODEL_DIR="./meshes"
MODELS=$(find "$MODEL_DIR" -type f -iname "*.obj")

for model in $MODELS; do
  build/size -t $model && build/evaluate $model -i1
  build/size -t $model 2>/dev/null
  if [[ $? == 0 ]]; then
    build/evaluate $model -i100 -r8
  else
    printf ">> File '%s' failed\n" "$model"
  fi
done

M main.cpp => main.cpp +14 -14
@@ 57,7 57,7 @@ void update_vectors() {
  plane_n = fplane_n.cast<double>().normalized();
}

void update_viewer(Mesh &mesh) {
void update_viewer(Viewer &viewer, Mesh &mesh) {
  forces_colors = MatrixXd::Zero(mesh.vertices.rows(), 3);
  for (int v = 0; v < mesh.vertices.rows(); v++) {
    auto force = mesh.forces.row(v);


@@ 69,9 69,14 @@ void update_viewer(Mesh &mesh) {
      forces_colors(v, i) = c[i];
  }

  draw_plane(viewer, mesh);
  viewer.data().set_mesh(mesh.vertices, mesh.faces);
  viewer.data().set_normals(mesh.normals);
  viewer.data().set_colors(mesh.face_colors);

  if (show_forces) {
    viewer.data().add_edges(mesh.vertices, forces_endpoints, forces_colors);
  }
}

// Draw the custom menu


@@ 125,7 130,6 @@ void draw_menu(Mesh &mesh) {

  if (ImGui::Button("Step ARAP")) {
    step_arap(mesh);
    update_viewer(mesh);
  }

  // Attempt at remeshing: Doesn't really work, since we probably need to tell it that we should


@@ 158,6 162,7 @@ void draw_menu(Mesh &mesh) {
      mesh.forces = MatrixXd::Zero(mesh.vertices.rows(), mesh.vertices.cols());
    recompute |= true;
  }

  if (stitching) {
    ImGui::Indent();
    ImGui::Checkbox("Show Forces", &show_forces);


@@ 165,7 170,8 @@ void draw_menu(Mesh &mesh) {
      recompute |= true;
    if (ImGui::Checkbox("Force ~ Area", &area_force))
      recompute |= true;
    if (ImGui::SliderFloat("Distance Cutoff", &force_dist_cutoff, 0.1, 10.0, "%.2f"))
    ImGui::Checkbox("Distance Cutoff", &force_cutoff);
    if (force_cutoff && ImGui::SliderFloat("Distance Cutoff", &force_dist_cutoff, 0.1, 10.0, "%.2f"))
      recompute |= true;
    ImGui::Unindent();
  }


@@ 207,7 213,6 @@ void draw_menu(Mesh &mesh) {

  if (recompute) {
    compute_things(mesh);
    update_viewer(mesh);
  }
}



@@ 264,10 269,6 @@ void draw_plane(Viewer &viewer, Mesh &mesh) {

  Eigen::RowVector3d c(1.0, 1.0, 1.0);
  viewer.data().add_edges(planeEA, planeEB, c);

  if (show_forces) {
    viewer.data().add_edges(mesh.vertices, forces_endpoints, forces_colors);
  }
}

int main(int argc, char *argv[]) {


@@ 333,7 334,6 @@ int main(int argc, char *argv[]) {
  menu.callback_draw_viewer_menu = [&menu]() { menu.draw_viewer_menu(); };

  viewer.callback_pre_draw = [&mesh](Viewer &viewer) {
    draw_plane(viewer, mesh);
    if (run_with_update) {
      if (step_arap(mesh) == StepReturn::Converged) {
        mesh.rest_state = mesh.vertices;


@@ 345,7 345,7 @@ int main(int argc, char *argv[]) {
        step_arap(mesh);
      }
    }
    update_viewer(mesh);
    update_viewer(viewer, mesh);
    return false;
  };



@@ 354,12 354,12 @@ int main(int argc, char *argv[]) {
      run_arap = !run_arap;
    } else if (key == '.') {
      step_arap(mesh);
      update_viewer(mesh);
      update_viewer(viewer, mesh);
    } else if (key == 'R') {
      mesh.vertices = mesh.rest_state = tetmesh_out.vertices;
      arap::precomputation(mesh, arap_data);
      compute_things(mesh);
      update_viewer(mesh);
      update_viewer(viewer, mesh);
    } else if (key == 'P') {
      Eigen::Matrix<unsigned char, Eigen::Dynamic, Eigen::Dynamic> R(1280, 800);
      Eigen::Matrix<unsigned char, Eigen::Dynamic, Eigen::Dynamic> G(1280, 800);


@@ 372,10 372,10 @@ int main(int argc, char *argv[]) {
  };

  viewer.data().clear();
  viewer.data().set_face_based(true);
  // viewer.data().set_face_based(true);
  viewer.data().set_mesh(mesh.vertices, mesh.faces);
  viewer.core().is_animating = true;
  compute_things(mesh);
  update_viewer(mesh);

  viewer.launch();
}

M surcut.cpp => surcut.cpp +6 -3
@@ 11,6 11,7 @@ std::vector<double> dots;
bool regular_arap = false;
float force_scale = 1.0;
bool area_force = false;
bool force_cutoff = false;
float force_dist_cutoff = 1.0;
int iteration_count = 0;
bool stitching = false;


@@ 253,8 254,9 @@ void compute_things(Mesh &mesh, double scale, bool only_plane) {

        Vector3d other_vertex = mesh.vertices.row(other_i);
        auto force_direction = other_vertex - vertex;
        if (force_direction.norm() > scale * force_dist_cutoff)
          continue;
        if (force_cutoff)
          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`.


@@ 264,8 266,9 @@ void compute_things(Mesh &mesh, double scale, bool only_plane) {
            force_direction * area * (double)force_scale * force_direction.squaredNorm() / scale;

        mesh.forces.row(v) += force;
        if (forces_endpoints_ptr)
        if (forces_endpoints_ptr) {
          forces_endpoints_ptr->row(v) = other_vertex;
        }
        mesh.forces.row(other_i) -= force; // Should this depend on area of `other_i`?
      }
    }

M surcut.h => surcut.h +2 -0
@@ 25,6 25,8 @@ extern bool regular_arap;
extern float force_scale;
// Multiply forces with the face area
extern bool area_force;
// Perform force distance cutoff
extern bool force_cutoff;
// Force distance cutoff
extern float force_dist_cutoff;
// The number of ARAP iterations performed