@@ 28,10 28,12 @@ bool key_down(Viewer &, unsigned char, int);
void draw_menu();
void update_vectors();
void ensure_moldability();
+void update_rest_state();
void remesh();
int ray_intersect(Vector3d, Vector3d, MatrixXd &, MatrixXi &);
void compute_forces(MatrixXd &, MatrixXi &);
void compute_things(MatrixXd &, MatrixXi &, bool = false);
+int intersect3D_RayTriangle(Vector3d, Vector3d, Vector3d, int, MatrixXd &, MatrixXi &, Vector3d *);
// Model data
MatrixXd V;
@@ 52,18 54,24 @@ VectorXi boundary(BOUNDARY_SIZE);
// ARAP'd vertices
MatrixXd U; /* Deformed vertices */
MatrixXd U0; /* Original U (this is really just TV) */
-
-std::vector<bool>
- is_point_above; /* whether this point is above or below the plane */
-std::vector<arap::Moldable::Enum>
- moldable; /* `Moldable` values for all faces */
-std::vector<arap::Orientation::Enum>
- orientation; /* `Moldable` values for all faces */
-std::vector<double> dots; /* dot products for each face */
+MatrixXd Forces;
+
+// Various computed things about the mesh
+// If a point is above the plane.
+std::vector<bool> is_point_above;
+// `Moldable` values for all faces
+std::vector<arap::Moldable::Enum> moldable;
+// `Orientation` values for all faces
+std::vector<arap::Orientation::Enum> orientation;
+// dot products with the removal direction for each face
+std::vector<double> dots;
arap::Data arap_data;
arap::Args arap_args;
+
+// Whether to stitch together unmoldable geometry
+bool stitching = false;
// Toggle running of ARAP
bool run_arap = false;
// Show rotations for all *faces*
@@ 89,14 97,15 @@ void update_vectors() {
plane_n = fplane_n.cast<double>().normalized();
}
+void update_rest_state() {
+ arap::precomputation(U, TT, TV.cols(), boundary, arap_data);
+}
+
// Draw the custom menu
void draw_menu() {
- bool px =
- ImGui::SliderFloat("plane_p.x", &fplane_p[0], mins.x(), maxs.x(), "%.4f");
- bool py =
- ImGui::SliderFloat("plane_p.y", &fplane_p[1], mins.y(), maxs.y(), "%.4f");
- bool pz =
- ImGui::SliderFloat("plane_p.z", &fplane_p[2], mins.z(), maxs.z(), "%.4f");
+ bool px = ImGui::SliderFloat("plane_p.x", &fplane_p[0], mins.x(), maxs.x(), "%.4f");
+ bool py = ImGui::SliderFloat("plane_p.y", &fplane_p[1], mins.y(), maxs.y(), "%.4f");
+ bool pz = ImGui::SliderFloat("plane_p.z", &fplane_p[2], mins.z(), maxs.z(), "%.4f");
bool nx = ImGui::SliderFloat("plane_n.x", &fplane_n[0], -1.0, 1.0, "%.2f");
bool ny = ImGui::SliderFloat("plane_n.y", &fplane_n[1], -1.0, 1.0, "%.2f");
@@ 124,25 133,24 @@ void draw_menu() {
}
if (ImGui::Button("Update Rest State")) {
- arap::precomputation(U, TT, TV.cols(), boundary, arap_data);
+ update_rest_state();
}
- if (ImGui::Button("Step ARAP"))
+ if (ImGui::Button("Step ARAP")) {
step_arap();
- compute_things(U, TF);
}
- // if (ImGui::Button("Ensure Moldability"))
- // ensure_moldability();
-
- // if (ImGui::Button("Remesh"))
- // remesh();
+ if (ImGui::Button("Compute Forces")) {
+ compute_forces(U, TF);
+ }
ImGui::Checkbox("Blend rotations", &arap_args.should_blend);
if (arap_args.should_blend) {
ImGui::SliderFloat("blend", &arap_args.blend, 0.0, 1.0, "%.2f");
}
+
+ ImGui::Checkbox("Stiching", &stitching);
ImGui::Checkbox("SLERP 2 rotations", &arap_args.slerp_2_rots);
ImGui::Checkbox("Rotate before mold check", &arap_args.rotate_before);
@@ 168,8 176,7 @@ void draw_menu() {
float aa = arap_data.mht_face_arap_angle[i];
int offset = sprintf(buffer, "%zu: ", i);
float limit = 0.01;
- if ((either && ((ma > limit) || (aa > limit))) ||
- (both && ((ma > limit) && (aa > limit)))) {
+ if ((either && ((ma > limit) || (aa > limit))) || (both && ((ma > limit) && (aa > limit)))) {
sprintf(buffer + offset, "%3.2f / %3.2f", ma, aa);
}
ImGui::Text("%s", buffer);
@@ 184,7 191,6 @@ bool key_down(Viewer &viewer, unsigned char key, int mod) {
run_arap = !run_arap;
} else if (key == '.') {
step_arap();
- compute_things(U, TF);
} else if (key == 'R') {
U = TV;
compute_things(U, TF);
@@ 198,15 204,13 @@ void step_arap() {
for (int i = 0; i < BOUNDARY_SIZE; i++)
bc.row(i) = TV.row(boundary(i));
- // XXX(note): this is for adding forces
- // Eigen::MatrixXd forces = Eigen::MatrixXd::Zero(TV.rows(), TV.cols());
- // for (int i = 0; i < TV.rows(); i++) {
- // // Go to the one fixed point
- // // forces.row(i) = (U.row(0) - U.row(i)) * 0.1;
- // }
- // arap_data.f_ext = forces;
-
arap::iter(arap_args, bc, arap_data, TT.rows(), TNeigh, U);
+ compute_things(U, TF);
+ if (stitching) {
+ compute_forces(U, TF);
+ } else {
+ Forces = MatrixXd::Zero(U.rows(), 3);
+ }
}
// Compute all the things that we need each time the vertex set changes. Here's
@@ 214,8 218,7 @@ void step_arap() {
// - `TN` (normals)
// - `is_point_above`
// - `colors`
-void compute_things(Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces,
- bool only_plane) {
+void compute_things(Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces, bool only_plane) {
if (!only_plane) { /* `normals`: depends on vertices */
igl::per_face_normals(vertices, faces, TN);
@@ 239,8 242,7 @@ void compute_things(Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces,
for (int f = 0; f < faces.rows(); f++) {
double dot;
Eigen::Vector3d removal_dir = plane_n;
- moldable[f] = arap::is_moldable(removal_dir, is_point_above, faces, f,
- TN.row(f), &dot);
+ moldable[f] = arap::is_moldable(removal_dir, is_point_above, faces, f, TN.row(f), &dot);
if (moldable[f] == arap::Moldable::Crossing) {
orientation[f] = arap::Orientation::Crossing;
} else {
@@ 277,9 279,7 @@ void compute_things(Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces,
// Here we look at the faces that are not moldable, and try to match them up with nearby faces; by
// joining up faces like this we hope to fill in cavities, such that the final object is moldalbe.
-void merge_step() {
-
-}
+void merge_step() {}
// This is called before each draw call. Compute normals and set right colors.
// XXX(note): much of this should probably just be done after `step_arap`, since
@@ 379,9 379,7 @@ void ensure_moldability() {
igl::per_face_normals(U, TF, normals);
std::vector<bool> is_point_above(U.rows());
- auto above_plane = [&](Eigen::Vector3d v) {
- return (v - plane_p).dot(plane_n) >= 0.0;
- };
+ auto above_plane = [&](Eigen::Vector3d v) { return (v - plane_p).dot(plane_n) >= 0.0; };
for (int i = 0; i < U.rows(); i++) {
is_point_above[i] = above_plane(U.row(i));
}
@@ 445,8 443,7 @@ void ensure_moldability() {
for (int i = 0; i < TF.rows(); i++) {
auto normal = normals.row(i);
Eigen::Vector3d removal_dir(plane_n);
- auto moldable =
- arap::is_moldable(removal_dir, is_point_above, TF, i, normal, NULL);
+ auto moldable = arap::is_moldable(removal_dir, is_point_above, TF, i, normal, NULL);
if (moldable == arap::Moldable::No) {
indices.push_back(i);
}
@@ 460,8 457,7 @@ void ensure_moldability() {
Eigen::Vector3d normal = normals.row(face);
Eigen::Vector3d removal_dir(plane_n);
- auto moldable =
- arap::is_moldable(removal_dir, is_point_above, TF, face, normal, NULL);
+ auto moldable = arap::is_moldable(removal_dir, is_point_above, TF, face, normal, NULL);
// XXX(todo): how to handle crossing faces?
if (moldable != arap::Moldable::No)
continue;
@@ 489,41 485,111 @@ void ensure_moldability() {
}
void compute_forces(Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces) {
- Eigen::Vector3d remove_dir = plane_n;
- // XXX(refactor): this is computed *a couple* of other places. Refactor.
+ Forces.resize(vertices.rows(), 3);
for (int f = 0; f < TF.rows(); f++) {
auto m = moldable[f];
// XXX(todo): what do we do about `Crossing`? For now, skip as well.
if (m != arap::Moldable::No)
continue;
+ Eigen::Vector3d ray_dir = (orientation[f] == arap::Orientation::Below) ? plane_n : Vector3d(-plane_n);
- for (int v_i = 0; v_i < 3; v_i++) {
- // XXX(todo): this will check each vertex once per non-moldable adjacent
- // face, when we only need to check once.
- auto vertex = vertices.row(faces(f, v_i));
- auto other_i = ray_intersect(vertex, remove_dir, vertices, faces);
+ for (int i = 0; i < 3; i++) {
+ int v_i = faces(f, i);
+ // XXX(todo): this will check each vertex once per non-moldable adjacent face, when we only
+ // need to check once.
+ Vector3d vertex = vertices.row(v_i);
+
+ int other_i = ray_intersect(ray_dir, vertex, vertices, faces);
if (other_i == -1)
continue;
- auto other_vertex = vertices.row(other_i);
+
+ Vector3d other_vertex = vertices.row(other_i);
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 */
// XXX(todo): these forces can be visualized with libigl (?)
+ Forces.row(v_i) = force_direction;
}
}
+ arap_data.f_ext = Forces;
}
// 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) {
+int ray_intersect(Eigen::Vector3d v, Eigen::Vector3d p, Eigen::MatrixXd &vertices,
+ Eigen::MatrixXi &faces) {
// XXX(todo): implement this
- assert(false);
+ Vector3d collision_point;
+ for (int f = 0; f < faces.rows(); f++) {
+
+ 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);
+ }
+ }
return -1;
}
+// 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)
+// Return: -1 = triangle is degenerate (a segment or point)
+// 0 = disjoint (no intersect)
+// 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) {
+#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))),
+ v2(verts.row(faces(triangle, 2)));
+ Vector3d u = v1 - v0;
+ Vector3d v = v2 - v0;
+
+ Vector3d w0 = ray_p - v0;
+ double a = -triangle_normal.dot(w0);
+ double b = triangle_normal.dot(ray_d);
+ if (abs(b) < SMALL_NUM) { // ray is parallel to triangle plane
+ if (abs(a) < SMALL_NUM) // ray lies in triangle plane
+ return 2;
+ else
+ return 0; // ray disjoint from plane
+ }
+
+ // get intersect point of ray with triangle plane
+ double r = a / b;
+ if (r < 0.0) // ray goes away from triangle
+ 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
+
+ // is I inside T?
+ double uu, uv, vv, wu, wv, D;
+ uu = u.dot(u);
+ uv = u.dot(v);
+ vv = v.dot(v);
+ Vector3d w = *I - v0;
+ wu = w.dot(u);
+ wv = w.dot(v);
+ D = uv * uv - uu * vv;
+
+ // get and test parametric coords
+ double s, t;
+ s = (uv * wv - vv * wu) / D;
+ if (s < 0.0 || s > 1.0) // I is outside T
+ return 0;
+ t = (uv * wu - uu * wv) / D;
+ if (t < 0.0 || (s + t) > 1.0) // I is outside T
+ return 0;
+
+ return 1; // I is in T
+}
+
// // Remesh the model based on the deformed version.
// // WARNING: THIS WILL DESTROY THE ORIGINAL MODEL, SO `reset` DOES NOTHING!
// void remesh() { NOTE: This is not remeshing, but triangulating..
@@ 531,13 597,13 @@ int ray_intersect(Eigen::Vector3d v, Eigen::Vector3d p,
// TetmeshOut tetmesh_out;
// tetmesh_in.vertices = &U;
// tetmesh_in.faces = &TF;
-//
+//
// int error = make_tetmesh(tetmesh_in, tetmesh_out);
// if (error) {
// std::cerr << "[ERR] tetesh returned " << error << std::endl;
// exit(1);
// }
-//
+//
// TV = tetmesh_out.vertices;
// U = TV;
// TT = tetmesh_out.tets;
@@ 592,16 658,17 @@ int main(int argc, char *argv[]) {
arap_data.mht_tf2t = TF2T;
arap_data.mht_plane_p = &plane_p;
arap_data.mht_plane_n = &plane_n;
- arap_data.mht_tet_rotations = std::vector<Eigen::Quaternion<double>>(
- TT.size(), Eigen::Quaternion<double>(1, 0, 0, 0));
- // arap_data.with_dynamics = true;
- // arap_data.h = 0.01;
- arap::precomputation(TV, TT, TV.cols(), boundary, arap_data);
- // arap_data.U0 = TV;
-
+ arap_data.mht_tet_rotations =
+ std::vector<Eigen::Quaternion<double>>(TT.size(), Eigen::Quaternion<double>(1, 0, 0, 0));
arap_data.vert_is_above = &is_point_above;
arap_data.moldable = &moldable;
arap_data.orientation = &orientation;
+ arap_data.U0 = TV;
+ arap_data.with_dynamics = true;
+ // XXX(todo): what's a reasonable value here?
+ arap_data.h = 0.01;
+ arap::precomputation(TV, TT, TV.cols(), boundary, arap_data);
+
maxs = V.colwise().maxCoeff().cast<float>();
mins = V.colwise().minCoeff().cast<float>();