M CMakeLists.txt => CMakeLists.txt +10 -4
@@ 14,10 14,16 @@ find_package(LIBIGL REQUIRED QUIET)
igl_download_stb()
find_package(OpenMP REQUIRED)
-file(GLOB SRCFILES *.cpp)
-add_executable(${PROJECT_NAME} ${SRCFILES} mesh.cpp)
-target_link_libraries(${PROJECT_NAME} igl::core igl::opengl igl::opengl_glfw igl::opengl_glfw_imgui
+set(SRCFILES surcut.cpp arap-mold.cpp mesh.cpp)
+
+add_executable(surcut ${SRCFILES} main.cpp)
+target_link_libraries(surcut igl::core igl::opengl igl::opengl_glfw igl::opengl_glfw_imgui
+ igl::tetgen OpenMP::OpenMP_CXX igl_stb_image)
+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_compile_options(${PROJECT_NAME} PRIVATE -g -O2 -Wall -Werror -Wpedantic -DDEBUG)
+target_compile_options(evaluate PRIVATE -g -O2 -Wall -Werror -Wpedantic -DDEBUG)
igl_download_tutorial_data()
M README.md => README.md +6 -1
@@ 1,7 1,12 @@
+# To be done
+
+ 1. Find/make models in which ARAP performs well, and in which stitching performs well. Probably
+ want to combine the two.
+ 2. Find some way of making good looking renderings of these results
+ 3. Actually write the report.
# Future Work
- Handle intersections when having forces
- - Visualize forces as springs: color code magnitude
M arap-mold.cpp => arap-mold.cpp +14 -12
@@ 1,15 1,15 @@
-#include "igl/arap_rhs.h"
-#include "igl/columnize.h"
-#include "igl/cotmatrix.h"
-#include "igl/covariance_scatter_matrix.h"
-#include "igl/fit_rotations.h"
-#include "igl/group_sum_matrix.h"
-#include "igl/massmatrix.h"
-#include "igl/mode.h"
-#include "igl/project_isometrically_to_plane.h"
-#include "igl/repdiag.h"
-#include "igl/slice.h"
-#include "igl/speye.h"
+#include <igl/arap_rhs.h>
+#include <igl/columnize.h>
+#include <igl/cotmatrix.h>
+#include <igl/covariance_scatter_matrix.h>
+#include <igl/fit_rotations.h>
+#include <igl/group_sum_matrix.h>
+#include <igl/massmatrix.h>
+#include <igl/mode.h>
+#include <igl/project_isometrically_to_plane.h>
+#include <igl/repdiag.h>
+#include <igl/slice.h>
+#include <igl/speye.h>
#include <cassert>
#include <iostream>
@@ 67,6 67,7 @@ bool arap::precomputation(const Mesh &mesh, arap::Data &data) {
data.n = n;
assert((b.size() == 0 || b.maxCoeff() < n) && "b out of bounds");
assert((b.size() == 0 || b.minCoeff() >= 0) && "b out of bounds");
+ assert(V == V && "arap-mold::pre V contains NaN");
data.b = b;
@@ 268,6 269,7 @@ void arap::iter(const Mesh &mesh, const arap::Args &args, const Eigen::MatrixXd
Eigen::MatrixXd arap_rotation(Rotations.block(0, 3 * tet_i, 3, 3));
face_normal = arap_rotation * face_normal;
}
+ assert(face_normal == face_normal && "arap-mod::iter face_normal contains NaN");
auto m = data.moldable->at(face_i);
switch (m) {
M arap-mold.h => arap-mold.h +2 -2
@@ 1,7 1,7 @@
#pragma once
-#include "igl/ARAPEnergyType.h"
-#include "igl/min_quad_with_fixed.h"
+#include <igl/ARAPEnergyType.h>
+#include <igl/min_quad_with_fixed.h>
#include <Eigen/Core>
#include <Eigen/Sparse>
A evaluate.cpp => evaluate.cpp +48 -0
@@ 0,0 1,48 @@
+#include <stdio.h>
+#include <Eigen/Core>
+
+#include "mesh.h"
+
+using Eigen::Vector3d;
+using Eigen::Vector4d;
+using Eigen::Matrix4d;
+
+double score_mesh(Mesh &);
+
+// 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 = 100;
+
+ double best_score = 0.0;
+ /* Even consider adaptive sampling? */
+ for (int i = 0; i < num_samples; i++) {
+ double score;
+ Vector3d dir /* = sample dir */;
+ Matrix4d projection /* = make this matrix */;
+
+ for (int f = 0; f < mesh.faces.rows(); f++) {
+ Vector4d normal = mesh.normals.row(f);
+ Vector4d transformed_normal = projection * normal;
+
+ (void)dir;
+ (void)transformed_normal;
+ }
+
+ best_score = std::max(best_score, score);
+ }
+
+ return best_score;
+}
+
+
+int main(int argc, char *argv[]) {
+ printf("Hello, world!\n");
+ return 0;
+}
M libigl => libigl +1 -1
@@ 1,1 1,1 @@
-Subproject commit f6b406427400ed7ddb56cfc2577b6af571827c8c
+Subproject commit 84c886b48294e80fed802a9f68ad64fb76ebb48c
M main.cpp => main.cpp +63 -303
@@ 1,9 1,9 @@
#include <igl/opengl/glfw/Viewer.h>
#include <igl/opengl/glfw/imgui/ImGuiHelpers.h>
#include <igl/opengl/glfw/imgui/ImGuiMenu.h>
-#include <igl/per_face_normals.h>
#include <igl/png/writePNG.h>
#include <igl/readOBJ.h>
+#include <igl/remesh_along_isoline.h>
#include <imgui/imgui.h>
#include <algorithm>
@@ 12,44 12,22 @@
#include "arap-mold.h"
#include "mesh.h"
#include "viridis.h"
-
-#define BOUNDARY_SIZE 4
+#include "surcut.h"
using igl::opengl::glfw::Viewer;
+using Eigen::Matrix4d;
using Eigen::MatrixXd;
using Eigen::MatrixXi;
using Eigen::Vector3d;
+using Eigen::Vector4d;
+using Eigen::VectorXd;
using Eigen::VectorXi;
-struct Ray {
- Vector3d src;
- Vector3d dir;
-};
-
-bool step_arap(Mesh &);
void draw_plane(Viewer &, Mesh &);
void draw_menu(Mesh &);
void update_vectors();
-int ray_intersect(Ray, int, Mesh &);
-void compute_things(Mesh &, bool = false);
-int intersect3D_RayTriangle(Ray, Vector3d, int, MatrixXd &, MatrixXi &, Vector3d *);
-
-// 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;
// Run arap, and when it converges update the rest state.
@@ 58,23 36,10 @@ bool run_with_update = 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;
-// Scale factor for all stitching forces
-float force_scale = 1.0;
-// Multiply forces with the face area
-bool area_force = false;
-// Force distance cutoff
-float force_dist_cutoff = 1.0;
-// The number of ARAP iterations performed
-int iteration_count = 0;
// Cut plane and point with `float`, for ImGui.
Eigen::Vector3f fplane_p;
Eigen::Vector3f fplane_n;
-// The real deal
-Eigen::Vector3d plane_p(0, 0, 0);
-Eigen::Vector3d plane_n(0, 1, 0);
// Model size properties
float scale; /* Minimal span for any of the three dimensions of the model */
@@ 92,6 57,23 @@ void update_vectors() {
plane_n = fplane_n.cast<double>().normalized();
}
+void update_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);
+ double ff = force.norm() / scale;
+ if (ff > 1.0)
+ ff = 1.0;
+ double *c = plasma(ff);
+ for (int i = 0; i < 3; i++)
+ forces_colors(v, i) = c[i];
+ }
+
+ viewer.data().set_mesh(mesh.vertices, mesh.faces);
+ viewer.data().set_normals(mesh.normals);
+ viewer.data().set_colors(mesh.face_colors);
+}
+
// Draw the custom menu
void draw_menu(Mesh &mesh) {
bool recompute = false;
@@ 114,12 96,18 @@ void draw_menu(Mesh &mesh) {
if (ImGui::Button("Run ARAP"))
run_arap = true;
} else {
- if (ImGui::Button("Stop ARAP"))
+ if (ImGui::Button("Stop [Run ARAP]"))
run_arap = false;
}
- if (ImGui::Button("Run and Update")) {
- run_with_update = true;
+ if (!run_with_update) {
+ if (ImGui::Button("Run and Update")) {
+ run_with_update = true;
+ }
+ } else {
+ if (ImGui::Button("Stop [Run and Update]")) {
+ run_with_update = false;
+ }
}
ImGui::Text("iter count: %d", iteration_count);
@@ 137,8 125,29 @@ 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
+ // 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;
+ // VectorXi birth_triangles;
+ // Eigen::SparseMatrix<double> bary_coords;
+ // VectorXi face_above_field;
+
+ // igl::remesh_along_isoline(mesh.vertices, mesh.faces, field, 0.0,
+ // new_vertices, new_faces, new_field, birth_triangles, bary_coords, face_above_field);
+
+ // mesh.vertices = new_vertices;
+ // mesh.faces = new_faces;
+ // recompute |= true;
+ // }
+
ImGui::Checkbox("Blend rotations", &arap_args.should_blend);
if (arap_args.should_blend) {
ImGui::SliderFloat("blend", &arap_args.blend, 0.0, 1.0, "%.2f");
@@ 198,160 207,7 @@ void draw_menu(Mesh &mesh) {
if (recompute) {
compute_things(mesh);
- }
-}
-
-// Perform one ARAP step. Return `true` if the difference between the previous and current state is
-// negligible.
-bool step_arap(Mesh &mesh) {
- iteration_count++;
- Eigen::MatrixXd bc(BOUNDARY_SIZE, 3);
- for (int i = 0; i < BOUNDARY_SIZE; i++)
- bc.row(i) = mesh.rest_state.row(mesh.boundary(i));
-
- arap_args.skip_mold_rotations = regular_arap;
- Eigen::MatrixXd prev_vertices = mesh.vertices;
- arap::iter(mesh, arap_args, bc, arap_data, mesh.vertices);
- compute_things(mesh);
-
- double max_diff = (mesh.vertices - prev_vertices).maxCoeff();
- return max_diff < scale * 5e-4; /* This factor was found experimentally! */
-}
-
-// Compute all the things that we need each time the vertex set changes. Here's
-// what we compute:
-// - `TN` (normals)
-// - `is_point_above`
-// - `colors`
-void compute_things(Mesh &mesh, bool only_plane) {
-
- if (!only_plane) { /* `normals`: depends on vertices */
- igl::per_face_normals(mesh.vertices, mesh.faces, mesh.normals);
- viewer.data().set_vertices(mesh.vertices);
- viewer.data().set_normals(mesh.normals);
- }
-
- { /* `is_point_above`: depends on vertices and the plane */
- is_point_above.resize(mesh.vertices.rows());
- for (int v_i = 0; v_i < mesh.vertices.rows(); v_i++) {
- Vector3d v(mesh.vertices.row(v_i));
- is_point_above[v_i] = (v - plane_p).dot(plane_n) >= 0.0;
- }
- }
-
- { /* `moldable`, `orientation`, and `dots`: depends on vertices and the plane
- */
- moldable.resize(mesh.faces.rows());
- orientation.resize(mesh.faces.rows());
- dots.resize(mesh.faces.rows());
- for (int f = 0; f < mesh.faces.rows(); f++) {
- double dot;
- Eigen::Vector3d removal_dir = plane_n;
- moldable[f] =
- arap::is_moldable(removal_dir, is_point_above, mesh.faces, f, mesh.normals.row(f), &dot);
- if (moldable[f] == arap::Moldable::Crossing) {
- orientation[f] = arap::Orientation::Crossing;
- } else {
- if (removal_dir == plane_n) {
- orientation[f] = arap::Orientation::Above;
- } else {
- orientation[f] = arap::Orientation::Below;
- }
- }
- dots[f] = dot;
- }
- }
-
- { /* `colors`: depends on dot */
- mesh.face_colors.resize(mesh.faces.rows(), 3);
- for (int f = 0; f < mesh.faces.rows(); f++) {
- switch (moldable[f]) {
- case arap::Moldable::Yes: {
- mesh.face_colors.row(f) = Eigen::Array3d(1.0, 1.0, 1.0);
- } break;
- case arap::Moldable::Crossing: {
- mesh.face_colors.row(f) = Eigen::Array3d(0.0, 0.0, 0.0);
- } break;
- default: {
- double r = -dots[f];
- r = sqrt(r);
- mesh.face_colors.row(f) = Eigen::Array3d(1.0, (1.0 - r), (1.0 - r));
- } break;
- }
- }
- viewer.data().set_colors(mesh.face_colors);
- }
-
- if (stitching) {
- mesh.forces = MatrixXd::Zero(mesh.vertices.rows(), 3);
- forces_colors = MatrixXd::Zero(mesh.vertices.rows(), 3);
- forces_endpoints = mesh.vertices;
- std::vector<bool> seen(mesh.vertices.rows());
- for (int f = 0; f < mesh.faces.rows(); f++) {
- auto m = moldable[f];
- // XXX(todo): what do we do about `Crossing`? For now, skip as well.
- if (m == arap::Moldable::Yes)
- continue;
- if (m == arap::Moldable::Crossing) {
- for (int i = 0; i < 3; i++) {
- int v = mesh.faces(f, i);
- seen[v] = true;
- mesh.forces.row(v) = Vector3d::Zero();
- }
- continue;
- }
- Eigen::Vector3d ray_dir =
- (orientation[f] == arap::Orientation::Below) ? plane_n : Vector3d(-plane_n);
-
- double area = 1.0;
- if (area_force) {
- Vector3d a = mesh.vertices.row(mesh.faces(f, 1)) - mesh.vertices.row(mesh.faces(f, 0));
- Vector3d b = mesh.vertices.row(mesh.faces(f, 2)) - mesh.vertices.row(mesh.faces(f, 0));
- area = (a.cross(b).norm()) / 2.0;
- }
-
- for (int i = 0; i < 3; i++) {
- int v = mesh.faces(f, i);
- // Only do a vertex one time; this is okay since the direction is +- plane_n, and not the
- // triangle normal.
- if (seen[v])
- continue;
- seen[v] = true;
-
- Vector3d vertex = mesh.vertices.row(v);
- Ray ray;
- ray.src = vertex;
- ray.dir = ray_dir;
-
- int other_i = ray_intersect(ray, v, mesh);
- if (other_i == -1)
- continue;
-
- Vector3d other_vertex = mesh.vertices.row(other_i);
- 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`.
-
- // Symmetric forces - think springs
- Vector3d force =
- force_direction * area * (double)force_scale * force_direction.squaredNorm() / scale;
-
- mesh.forces.row(v) += force;
- forces_endpoints.row(v) = other_vertex;
- double ff = force.norm() / scale;
- if (ff > 1.0)
- ff = 1.0;
- double *c = plasma(ff);
- for (int i = 0; i < 3; i++)
- forces_colors(v, i) = c[i];
-
- mesh.forces.row(other_i) -= force; // Should this depend on area of `other_i`?
- }
- }
-
- arap_data.f_ext = mesh.forces;
+ update_viewer(mesh);
}
}
@@ 414,107 270,6 @@ void draw_plane(Viewer &viewer, Mesh &mesh) {
}
}
-// 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(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!
-
- Vector3d collision_point;
- double best_dist = std::numeric_limits<double>::max();
- int best_vertex = -1;
-
- for (int f = 0; f < mesh.faces.rows(); f++) {
- int res = intersect3D_RayTriangle(ray, mesh.normals.row(f), f, mesh.vertices, mesh.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 (mesh.faces(f, i) == src_vertex)
- is_own_face = true;
- if (is_own_face)
- continue;
-
- for (int i = 0; i < 3; i++) {
- int v = mesh.faces(f, i);
- double norm = (collision_point - Vector3d(mesh.vertices.row(v))).norm();
- if (norm < dist_to_collision) {
- dist_to_collision = norm;
- closest_vertex = v;
- }
- }
-
- double dist_to_ray_pos = (Vector3d(mesh.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 best_vertex;
-}
-
-// 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(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))),
- v2(verts.row(faces(triangle, 2)));
- Vector3d u = v1 - v0;
- Vector3d v = v2 - v0;
-
- Vector3d w0 = ray.src - v0;
- double a = -triangle_normal.dot(w0);
- 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;
- 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.src + r * ray.dir; // 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
-}
-
int main(int argc, char *argv[]) {
int error;
using namespace std;
@@ 557,6 312,8 @@ int main(int argc, char *argv[]) {
mesh_print_sizes(mesh);
+ 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;
@@ 593,18 350,20 @@ int main(int argc, char *argv[]) {
if (step_arap(mesh)) {
static int last_iter_update = -99;
// If we have converged with rest state updates, stop.
- if (iteration_count - last_iter_update < 3) { /* This number was also found experimentally */
- run_with_update = false;
+ if (iteration_count - last_iter_update == 1) {
+ // run_with_update = false;
}
fprintf(stderr, "Rest state update (diff=%d)\n", iteration_count - last_iter_update);
mesh.rest_state = mesh.vertices;
arap::precomputation(mesh, arap_data);
compute_things(mesh);
+ update_viewer(mesh);
last_iter_update = iteration_count;
}
} else {
if (run_arap) {
step_arap(mesh);
+ update_viewer(mesh);
}
}
return false;
@@ 615,10 374,12 @@ int main(int argc, char *argv[]) {
run_arap = !run_arap;
} else if (key == '.') {
step_arap(mesh);
+ update_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);
} 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);
@@ 631,11 392,10 @@ int main(int argc, char *argv[]) {
};
viewer.data().clear();
- viewer.data().set_mesh(mesh.vertices, mesh.faces);
viewer.data().set_face_based(true);
viewer.core().is_animating = true;
-
compute_things(mesh);
+ update_viewer(mesh);
viewer.launch();
}
A surcut.cpp => surcut.cpp +273 -0
@@ 0,0 1,273 @@
+#include <igl/per_face_normals.h>
+
+#include "surcut.h"
+#include "mesh.h"
+#include "arap-mold.h"
+
+std::vector<bool> is_point_above;
+std::vector<arap::Moldable::Enum> moldable;
+std::vector<arap::Orientation::Enum> orientation;
+std::vector<double> dots;
+bool regular_arap = false;
+float force_scale = 1.0;
+bool area_force = false;
+float force_dist_cutoff = 1.0;
+int iteration_count = 0;
+bool stitching = false;
+MatrixXd *forces_endpoints_ptr = nullptr;
+
+arap::Data arap_data;
+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`.
+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!
+
+ Vector3d collision_point;
+ double best_dist = std::numeric_limits<double>::max();
+ int best_vertex = -1;
+
+ for (int f = 0; f < mesh.faces.rows(); f++) {
+ int res = intersect3D_RayTriangle(ray, mesh.normals.row(f), f, mesh.vertices, mesh.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 (mesh.faces(f, i) == src_vertex)
+ is_own_face = true;
+ if (is_own_face)
+ continue;
+
+ for (int i = 0; i < 3; i++) {
+ int v = mesh.faces(f, i);
+ double norm = (collision_point - Vector3d(mesh.vertices.row(v))).norm();
+ if (norm < dist_to_collision) {
+ dist_to_collision = norm;
+ closest_vertex = v;
+ }
+ }
+
+ double dist_to_ray_pos = (Vector3d(mesh.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 best_vertex;
+}
+
+// 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(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))),
+ v2(verts.row(faces(triangle, 2)));
+ Vector3d u = v1 - v0;
+ Vector3d v = v2 - v0;
+
+ Vector3d w0 = ray.src - v0;
+ double a = -triangle_normal.dot(w0);
+ 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;
+ 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.src + r * ray.dir; // 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
+}
+
+// Perform one ARAP step. Return `true` if the difference between the previous and current state is
+// negligible.
+bool step_arap(Mesh &mesh, double scale) {
+ iteration_count++;
+ Eigen::MatrixXd bc(mesh.boundary.rows(), 3);
+ for (int i = 0; i < mesh.boundary.rows(); i++)
+ bc.row(i) = mesh.rest_state.row(mesh.boundary(i));
+
+ arap_args.skip_mold_rotations = regular_arap;
+ Eigen::MatrixXd prev_vertices = mesh.vertices;
+ arap::iter(mesh, arap_args, bc, arap_data, mesh.vertices);
+ if (mesh.vertices != mesh.vertices) { /* Contains NaN */
+ mesh.vertices = prev_vertices;
+ return true;
+ }
+ compute_things(mesh);
+
+ double max_diff = (mesh.vertices - prev_vertices).maxCoeff();
+ return max_diff < scale * 5e-4; /* This factor was found experimentally! */
+}
+
+// Compute all the things that we need each time the vertex set changes. Here's what we compute:
+// - `TN` (normals)
+// - `is_point_above`
+// - `colors`
+void compute_things(Mesh &mesh, double scale, bool only_plane) {
+ if (!only_plane) { /* `normals`: depends on vertices */
+ igl::per_face_normals(mesh.vertices, mesh.faces, mesh.normals);
+ }
+
+ { /* `is_point_above`: depends on vertices and the plane */
+ is_point_above.resize(mesh.vertices.rows());
+ for (int v_i = 0; v_i < mesh.vertices.rows(); v_i++) {
+ Vector3d v(mesh.vertices.row(v_i));
+ is_point_above[v_i] = (v - plane_p).dot(plane_n) >= 0.0;
+ }
+ }
+
+ { /* `moldable`, `orientation`, and `dots`: depends on vertices and the plane
+ */
+ moldable.resize(mesh.faces.rows());
+ orientation.resize(mesh.faces.rows());
+ dots.resize(mesh.faces.rows());
+ for (int f = 0; f < mesh.faces.rows(); f++) {
+ double dot;
+ Eigen::Vector3d removal_dir = plane_n;
+ moldable[f] =
+ arap::is_moldable(removal_dir, is_point_above, mesh.faces, f, mesh.normals.row(f), &dot);
+ if (moldable[f] == arap::Moldable::Crossing) {
+ orientation[f] = arap::Orientation::Crossing;
+ } else {
+ if (removal_dir == plane_n) {
+ orientation[f] = arap::Orientation::Above;
+ } else {
+ orientation[f] = arap::Orientation::Below;
+ }
+ }
+ dots[f] = dot;
+ }
+ }
+
+ { /* `colors`: depends on dot */
+ mesh.face_colors.resize(mesh.faces.rows(), 3);
+ for (int f = 0; f < mesh.faces.rows(); f++) {
+ switch (moldable[f]) {
+ case arap::Moldable::Yes: {
+ mesh.face_colors.row(f) = Eigen::Array3d(1.0, 1.0, 1.0);
+ } break;
+ case arap::Moldable::Crossing: {
+ mesh.face_colors.row(f) = Eigen::Array3d(0.0, 0.0, 0.0);
+ } break;
+ default: {
+ double r = -dots[f];
+ r = sqrt(r);
+ mesh.face_colors.row(f) = Eigen::Array3d(1.0, (1.0 - r), (1.0 - r));
+ } break;
+ }
+ }
+ }
+
+ if (forces_endpoints_ptr) {
+ *forces_endpoints_ptr = mesh.vertices;
+ }
+ if (stitching) {
+ mesh.forces = MatrixXd::Zero(mesh.vertices.rows(), 3);
+ std::vector<bool> seen(mesh.vertices.rows());
+ for (int f = 0; f < mesh.faces.rows(); f++) {
+ auto m = moldable[f];
+ // XXX(todo): what do we do about `Crossing`? For now, skip as well.
+ if (m == arap::Moldable::Yes)
+ continue;
+ if (m == arap::Moldable::Crossing) {
+ for (int i = 0; i < 3; i++) {
+ int v = mesh.faces(f, i);
+ seen[v] = true;
+ mesh.forces.row(v) = Vector3d::Zero();
+ }
+ continue;
+ }
+ Eigen::Vector3d ray_dir =
+ (orientation[f] == arap::Orientation::Below) ? plane_n : Vector3d(-plane_n);
+
+ double area = 1.0;
+ if (area_force) {
+ Vector3d a = mesh.vertices.row(mesh.faces(f, 1)) - mesh.vertices.row(mesh.faces(f, 0));
+ Vector3d b = mesh.vertices.row(mesh.faces(f, 2)) - mesh.vertices.row(mesh.faces(f, 0));
+ area = (a.cross(b).norm()) / 2.0;
+ }
+
+ for (int i = 0; i < 3; i++) {
+ int v = mesh.faces(f, i);
+ // Only do a vertex one time; this is okay since the direction is +- plane_n, and not the
+ // triangle normal.
+ if (seen[v])
+ continue;
+ seen[v] = true;
+
+ Vector3d vertex = mesh.vertices.row(v);
+ Ray ray;
+ ray.src = vertex;
+ ray.dir = ray_dir;
+
+ int other_i = ray_intersect(ray, v, mesh);
+ if (other_i == -1)
+ continue;
+
+ Vector3d other_vertex = mesh.vertices.row(other_i);
+ 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`.
+
+ // Symmetric forces - think springs
+ Vector3d force =
+ force_direction * area * (double)force_scale * force_direction.squaredNorm() / scale;
+
+ mesh.forces.row(v) += force;
+ 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`?
+ }
+ }
+
+ arap_data.f_ext = mesh.forces;
+ }
+}
A surcut.h => surcut.h +52 -0
@@ 0,0 1,52 @@
+#pragma once
+
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+#include "mesh.h"
+#include "arap-mold.h"
+
+using Eigen::Vector3d;
+using Eigen::MatrixXd;
+using Eigen::MatrixXi;
+
+// Various computed things about the mesh
+// If a point is above the plane.
+extern std::vector<bool> is_point_above;
+// `Moldable` values for all faces
+extern std::vector<arap::Moldable::Enum> moldable;
+// `Orientation` values for all faces
+extern std::vector<arap::Orientation::Enum> orientation;
+// dot products with the removal direction for each face
+extern std::vector<double> dots;
+// Don't use the ARAP enforcing moldability. (only makes sense to use wiht `stitching`).
+extern bool regular_arap;
+// Scale factor for all stitching forces
+extern float force_scale;
+// Multiply forces with the face area
+extern bool area_force;
+// Force distance cutoff
+extern float force_dist_cutoff;
+// The number of ARAP iterations performed
+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`.
+extern MatrixXd *forces_endpoints_ptr;
+
+extern arap::Data arap_data;
+extern arap::Args arap_args;
+
+extern Eigen::Vector3d plane_p;
+extern Eigen::Vector3d plane_n;
+
+struct Ray {
+ Vector3d src;
+ Vector3d dir;
+};
+
+int ray_intersect(Ray, int, Mesh &);
+int intersect3D_RayTriangle(Ray , Vector3d , int , MatrixXd &, MatrixXi &, Vector3d *);
+bool step_arap(Mesh &, double = 1.0);
+void compute_things(Mesh &mesh, double = 1.0, bool only_plane = false);