## ~mht/surcut

a4ac05d43ad83df81b7ee6fb255f34ce429e8197 — Martin Hafskjold Thoresen 5 years ago
```Intersect with the first triangle, and set force to a vertex

Not sure the vertex is really what we want, but here it is..
```
```4 files changed, 71 insertions(+), 20 deletions(-)

A .clang-format
M arap-mold.cpp
M arap-mold.h
M main.cpp
```
`A .clang-format => .clang-format +1 -0`
```@@ 0,0 1,1 @@
+ColumnLimit: 99

```
`M arap-mold.cpp => arap-mold.cpp +1 -1`
```@@ 223,7 223,7 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &d
igl::fit_rotations(StackCovMats, true, Rotations);
}

-  { /* ~~ Moldability ~~ */
+  if (!args.skip_mold_rotations) { /* ~~ Moldability ~~ */
// 1. Go through all the boundary faces of the tet mesh
// 2. Skip if the face satisfies direction constraint
// 3. Compute minimal rotation to satisfy constraint

```
`M arap-mold.h => arap-mold.h +1 -0`
```@@ 99,6 99,7 @@ struct Args {
float blend;        /* blend factor [0, 1] */
bool slerp_2_rots;  /* SLERP rotations to join 2. Otherwise multiply by 0.5 */
bool rotate_before; /* Rotate normal before checking moldability */
+  bool skip_mold_rotations; /* Don't do the rotation stuff to it make moldable */

Args() : should_blend(false), blend(0.0), slerp_2_rots(false), rotate_before(false){};
};

```
`M main.cpp => main.cpp +68 -19`
```@@ 21,6 21,11 @@ using Eigen::MatrixXi;
using Eigen::Vector3d;
using Eigen::VectorXi;

+struct Ray {
+  Vector3d src;
+  Vector3d dir;
+};
+
void step_arap();
bool pre_draw(Viewer &);
void draw_plane(Viewer &);

@@ 30,9 35,9 @@ void update_vectors();
void ensure_moldability();
void update_rest_state();
void remesh();
-int ray_intersect(Vector3d, Vector3d, MatrixXd &, MatrixXi &);
+int ray_intersect(Ray, int, MatrixXd &, MatrixXi &);
void compute_things(MatrixXd &, MatrixXi &, bool = false);
-int intersect3D_RayTriangle(Vector3d, Vector3d, Vector3d, int, MatrixXd &, MatrixXi &, Vector3d *);
+int intersect3D_RayTriangle(Ray, Vector3d, int, MatrixXd &, MatrixXi &, Vector3d *);

// Model data
MatrixXd V;

@@ 77,6 82,8 @@ bool run_arap = false;
bool show_list = false;
// Visualize the forces on each vertex when stitching
bool show_forces = false;
+// Don't use the ARAP enforcing moldability.  (only makes sense to use wiht `stitching`).
+bool regular_arap = false;

// Cut plane and point with `float`, for ImGui.
Eigen::Vector3f fplane_p;

@@ 155,6 162,7 @@ void draw_menu() {
Forces = MatrixXd::Zero(U.rows(), 3);
recompute |= true;
}
+  ImGui::Checkbox("Skip Moldability Rotations", &regular_arap);

ImGui::Checkbox("SLERP 2 rotations", &arap_args.slerp_2_rots);
ImGui::Checkbox("Rotate before mold check", &arap_args.rotate_before);

@@ 212,6 220,7 @@ void step_arap() {
for (int i = 0; i < BOUNDARY_SIZE; i++)
bc.row(i) = RestState.row(boundary(i));

+  arap_args.skip_mold_rotations = regular_arap;
arap::iter(arap_args, bc, arap_data, TT.rows(), TNeigh, U);
compute_things(U, TF);
}

@@ 285,8 294,16 @@ void compute_things(Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces, bool only
for (int f = 0; f < faces.rows(); f++) {
auto m = moldable[f];
// XXX(todo): what do we do about `Crossing`? For now, skip as well.
-      if (m != arap::Moldable::No)
+      if (m == arap::Moldable::Yes)
+        continue;
+      if (m == arap::Moldable::Crossing) {
+        for (int i = 0; i < 3; i++) {
+          int v = faces(f, i);
+          seen[v] = true;
+          Forces.row(v) = Vector3d::Zero();
+        }
continue;
+      }
Eigen::Vector3d ray_dir =
(orientation[f] == arap::Orientation::Below) ? plane_n : Vector3d(-plane_n);

@@ 299,8 316,11 @@ void compute_things(Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces, bool only
seen[v] = true;

Vector3d vertex = vertices.row(v);
+        Ray ray;
+        ray.src = vertex;
+        ray.dir = ray_dir;

-        int other_i = ray_intersect(ray_dir, vertex, vertices, faces);
+        int other_i = ray_intersect(ray, v, vertices, faces);
if (other_i == -1)
continue;

@@ 308,10 328,13 @@ void compute_things(Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces, bool only
auto force_direction = other_vertex - vertex;
// XXX(todo): what should the magnitude of the force be? Should probably take a closer look
// at the dynamics implementation in `arap-mold.cpp`.
-        (void)force_direction; /* unused warning */
-        Forces.row(v) = force_direction;
+
+        // Symmetric forces - think springs
+        Forces.row(v) += force_direction * 0.48;
+        Forces.row(other_i) -= force_direction * 0.48;
}
}
+
arap_data.f_ext = Forces;
}
}

@@ 391,19 414,45 @@ void draw_plane(Viewer &viewer) {
// 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`.
-int ray_intersect(Eigen::Vector3d v, Eigen::Vector3d p, Eigen::MatrixXd &vertices,
-                  Eigen::MatrixXi &faces) {
-  // XXX(todo): implement this
+int ray_intersect(Ray ray, int src_vertex, Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces) {
+  // 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!
+
Vector3d collision_point;
+  double best_dist = std::numeric_limits<double>::max();
+  int best_vertex = -1;
+
for (int f = 0; f < faces.rows(); f++) {
+    int res = intersect3D_RayTriangle(ray, TN.row(f), f, vertices, faces, &collision_point);
+    if (res != 1)
+      continue;
+
+    double dist_to_collision = std::numeric_limits<double>::max();
+    int closest_vertex = -1;
+
+    bool is_own_face = false;
+    for (int i = 0; i < 3; i++)
+      if (faces(f, i) == src_vertex)
+        is_own_face = true;
+    if (is_own_face)
+      continue;
+
+    for (int i = 0; i < 3; i++) {
+      int v = faces(f, i);
+      double norm = (collision_point - Vector3d(vertices.row(v))).norm();
+      if (norm < dist_to_collision) {
+        dist_to_collision = norm;
+        closest_vertex = v;
+      }
+    }

-    int res = intersect3D_RayTriangle(p, v, TN.row(f), f, vertices, faces, &collision_point);
-    if (res == 1) {
-      // XXX(todo): return closest to `collision_point`
-      return faces(f, 0);
+    double dist_to_ray_pos = (Vector3d(vertices.row(closest_vertex)) - ray.src).norm();
+    if (dist_to_ray_pos < best_dist) {
+      best_dist = dist_to_ray_pos;
+      best_vertex = closest_vertex;
}
}
-  return -1;
+  return best_vertex;
}

// intersect3D_RayTriangle(): find the 3D intersection of a ray with a triangle

@@ 414,8 463,8 @@ int ray_intersect(Eigen::Vector3d v, Eigen::Vector3d p, Eigen::MatrixXd &vertice
//             1 =  intersect in unique point I1
//             2 =  are in the same plane
// Found at http://geomalgorithms.com/a06-_intersect-2.html#intersect3D_RayTriangle
-int intersect3D_RayTriangle(Vector3d ray_p, Vector3d ray_d, Vector3d triangle_normal, int triangle,
-                            MatrixXd &verts, MatrixXi &faces, Vector3d *I) {
+int intersect3D_RayTriangle(Ray ray, Vector3d triangle_normal, int triangle, MatrixXd &verts,
+                            MatrixXi &faces, Vector3d *I) {
#define SMALL_NUM 0.00001
// get triangle edge vectors and plane normal
Vector3d v0(verts.row(faces(triangle, 0))), v1(verts.row(faces(triangle, 1))),

@@ 423,9 472,9 @@ int intersect3D_RayTriangle(Vector3d ray_p, Vector3d ray_d, Vector3d triangle_no
Vector3d u = v1 - v0;
Vector3d v = v2 - v0;

-  Vector3d w0 = ray_p - v0;
+  Vector3d w0 = ray.src - v0;
double a = -triangle_normal.dot(w0);
-  double b = triangle_normal.dot(ray_d);
+  double b = triangle_normal.dot(ray.dir);
if (abs(b) < SMALL_NUM) { // ray is  parallel to triangle plane
if (abs(a) < SMALL_NUM) // ray lies in triangle plane
return 2;

@@ 439,7 488,7 @@ int intersect3D_RayTriangle(Vector3d ray_p, Vector3d ray_d, Vector3d triangle_no
return 0;  // => no intersect
// for a segment, also test if (r > 1.0) => no intersect

-  *I = ray_p + r * ray_d; // 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;

```