M arap-mold.cpp => arap-mold.cpp +69 -59
@@ 7,7 7,6 @@
#include "igl/group_sum_matrix.h"
#include "igl/massmatrix.h"
#include "igl/mode.h"
-#include "igl/per_face_normals.h"
#include "igl/project_isometrically_to_plane.h"
#include "igl/repdiag.h"
#include "igl/slice.h"
@@ 24,18 23,19 @@ typedef Eigen::MatrixXi Derivedb;
typedef Eigen::MatrixXd Derivedbc;
typedef Eigen::MatrixXd DerivedU;
-Eigen::Quaternion<double> quat_angle_multiply(Eigen::Quaternion<double> qt, double f) {
+Eigen::Quaternion<double> quat_angle_multiply(Eigen::Quaternion<double> qt,
+ double f) {
Eigen::AngleAxis<double> aa(qt);
aa.angle() *= f;
return Eigen::Quaternion<double>(aa);
}
-int color_from_dot(double dot) { return 1 + (int(100 * -dot)); }
-
// Check whether `face` is moldable, given the plane
-arap::Moldable arap::is_moldable(Eigen::Vector3d &remove_dir, std::vector<bool> &is_point_above,
- Eigen::MatrixXi &faces, int face, Eigen::Vector3d normal,
- double *out_dot) {
+arap::Moldable::Enum arap::is_moldable(Eigen::Vector3d &remove_dir,
+ std::vector<bool> &is_point_above,
+ Eigen::MatrixXi &faces, int face,
+ Eigen::Vector3d normal,
+ double *out_dot) {
bool x_above = is_point_above[faces(face, 0)];
bool y_above = is_point_above[faces(face, 1)];
@@ 46,19 46,20 @@ arap::Moldable arap::is_moldable(Eigen::Vector3d &remove_dir, std::vector<bool>
if (all_below)
remove_dir = -remove_dir;
else if (!all_above)
- return CROSSING;
+ return Moldable::Crossing;
auto dot = normal.dot(remove_dir);
if (out_dot)
*out_dot = dot;
if (dot < 0.0)
- return NO;
- return YES;
+ return Moldable::No;
+ return Moldable::Yes;
}
-bool arap::precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const int dim,
- const Eigen::MatrixXi &b, arap::Data &data) {
+bool arap::precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F,
+ const int dim, const Eigen::MatrixXi &b,
+ arap::Data &data) {
using namespace std;
using namespace Eigen;
typedef typename DerivedV::Scalar Scalar;
@@ 170,13 171,16 @@ bool arap::precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, co
data.vel = MatrixXd::Zero(n, data.dim);
}
- return min_quad_with_fixed_precompute(Q, b, SparseMatrix<double>(), true, data.solver_data);
+ return min_quad_with_fixed_precompute(Q, b, SparseMatrix<double>(), true,
+ data.solver_data);
}
-// # Note: if ran with `use_dynamics` you must set `ARAPData::U0` to be your initial `U` before
-// calling this. Otherwise it'll be 0, an you will have a bad time.
-void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &data, int num_tets,
- Eigen::MatrixXi &TNeigh, Eigen::MatrixXd &U) {
+// # Note: if ran with `use_dynamics` you must set `ARAPData::U0` to be your
+// initial `U` before calling this. Otherwise it'll be 0, an you will have a bad
+// time.
+void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc,
+ arap::Data &data, int num_tets, Eigen::MatrixXi &TNeigh,
+ Eigen::MatrixXd &U) {
using namespace Eigen;
using namespace std;
assert(data.b.size() == bc.rows());
@@ 187,7 191,8 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &d
if (U.size() == 0) {
// terrible initial guess.. should at least copy input mesh
#ifndef NDEBUG
- cerr << "arap_solve: Using terrible initial guess for U. Try U = V." << endl;
+ cerr << "arap_solve: Using terrible initial guess for U. Try U = V."
+ << endl;
#endif
U = MatrixXd::Zero(data.n, data.dim);
} else {
@@ 215,8 220,8 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &d
StackCovMats /= StackCovMats.array().abs().maxCoeff();
const int Rdim = data.dim;
- // This is a big matrix, which contains stacked 3x3 rotation matrices for each *tet* in the
- // mesh.
+ // This is a big matrix, which contains stacked 3x3 rotation matrices for each
+ // *tet* in the mesh.
MatrixXd Rotations(Rdim, data.CSM.rows());
if (Rdim == 2) {
@@ 238,11 243,8 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &d
for (int i = 0; i < U.rows(); i++)
is_point_above[i] = above_plane(U.row(i));
- // XXX(todo): compute this somewhere else
- Eigen::MatrixXd normals;
- igl::per_face_normals(U, data.mht_triangles, normals);
-
- // If we have multiple outward faces for a single tet we must combine rotations.
+ // If we have multiple outward faces for a single tet we must combine
+ // rotations.
std::vector<std::vector<Eigen::Quaternion<double>>> tet_rotations(num_tets);
for (int i = 0; i < num_tets; i++) {
tet_rotations.push_back(std::vector<Eigen::Quaternion<double>>());
@@ 251,10 253,11 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &d
data.mht_face_mold_angle.clear();
data.mht_face_arap_angle.clear();
// #1: All faces returned by tetgen are at the boundary.
- for (int face_i = 0; face_i < data.mht_triangles.rows(); face_i++) {
+ for (int face_i = 0; face_i < data.faces->rows(); face_i++) {
- // We need to check if multiple boundary faces belong to the same tet. Since we're on the
- // boundary, one of these tets will be `-1`, and the other is the one we want.
+ // We need to check if multiple boundary faces belong to the same tet.
+ // Since we're on the boundary, one of these tets will be `-1`, and the
+ // other is the one we want.
auto tets = data.mht_tf2t.row(face_i);
auto tet_i = tets(0) == -1 ? tets(1) : tets(0);
assert(tet_i < num_tets && tet_i >= 0);
@@ 266,35 269,33 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &d
data.mht_face_arap_angle.push_back(arap_aa.angle() * 180 / M_PI);
// Get the normal of this face, possible rotated if we want that.
- Eigen::Vector3d face_normal(normals.row(face_i));
+ Eigen::Vector3d face_normal(data.normals->row(face_i));
if (args.rotate_before) {
Eigen::MatrixXd arap_rotation(Rotations.block(0, 3 * tet_i, 3, 3));
face_normal = arap_rotation * face_normal;
}
- double dot;
- // XXX(todo): this is awkward, but we need to know if we've negated d
- Eigen::Vector3d remove_dir = *data.mht_plane_n;
- auto moldable =
- is_moldable(remove_dir, is_point_above, data.mht_triangles, face_i, face_normal, &dot);
-
- switch (moldable) {
- case YES: {
- data.mht_face_classes[face_i] = 0;
+ auto m = data.moldable->at(face_i);
+ switch (m) {
+ case Moldable::Yes: {
data.mht_face_mold_angle.push_back(0.0);
} break;
- case CROSSING: {
- data.mht_face_classes[face_i] = -1;
+ case Moldable::Crossing: {
data.mht_face_mold_angle.push_back(0.0);
continue;
} break;
- case NO: {
+ case Moldable::No: {
// Must rotate the tet before finding how much further we must rotate.
// Note that we check the `dot` with the unrotated normal!
- // Eigen::Vector3d face_normal_rot = arap_rotation * normals.row(face_i).transpose();
- data.mht_face_classes[face_i] = color_from_dot(dot);
- // XXX(note): this is oposite from the paper, but I'm pretty sure this is the way to go.
- Eigen::Vector3d rotation_axis = face_normal.cross(remove_dir).normalized();
+ // Eigen::Vector3d face_normal_rot = arap_rotation *
+ // data.normals->row(face_i).transpose();
+ // XXX(note): this is oposite from the paper, but I'm pretty sure this
+ // is the way to go.
+ auto remove_dir = *data.mht_plane_n;
+ if (data.orientation->at(face_i) == Orientation::Below)
+ remove_dir = -remove_dir;
+ Eigen::Vector3d rotation_axis =
+ face_normal.cross(remove_dir).normalized();
auto x = remove_dir.dot(face_normal);
auto angle_between = atan2(1.0, x);
@@ 305,10 306,12 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &d
angle_between -= M_PI / 2.0;
data.mht_face_mold_angle.push_back(angle_between * 180 / M_PI);
- // #4: In the case where a tet has multiple boundary faces, we'll get multiple rotations.
- // XXX(note): We're assuming this uses right-hand rotation, but the Eigen docs doens't say
- // anything about it(!!).
- Eigen::Quaternion<double> qt(Eigen::AngleAxis<double>(angle_between, rotation_axis));
+ // #4: In the case where a tet has multiple boundary faces, we'll get
+ // multiple rotations.
+ // XXX(note): We're assuming this uses right-hand rotation, but the
+ // Eigen docs doens't say anything about it(!!).
+ Eigen::Quaternion<double> qt(
+ Eigen::AngleAxis<double>(angle_between, rotation_axis));
tet_rotations[tet_i].push_back(qt);
} break;
}
@@ 335,8 338,10 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &d
}
} break;
default: {
- // XXX(note): This probably doesn't make much sense, but I'm not sure what to do.
- fprintf(stderr, "[WARN]: we're interpolating quats here (tet with >2 boundary faces)\n");
+ // XXX(note): This probably doesn't make much sense, but I'm not sure
+ // what to do.
+ fprintf(stderr, "[WARN]: we're interpolating quats here (tet with >2 "
+ "boundary faces)\n");
Eigen::Quaternion<double> q(1, 0, 0, 0);
double d = 1.0 / (double)tet_rotations[i].size();
for (size_t ii = 0; ii < tet_rotations[i].size(); ii++) {
@@ 364,8 369,10 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &d
for (int n_i = 0; n_i < 4; n_i++) {
int neigh = TNeigh(tet_i, n_i);
if (neigh >= 0) {
- Eigen::Quaternion<double> neigh_qt = merged_rotations[TNeigh(tet_i, n_i)];
- neigh_qt = quat_angle_multiply(neigh_qt, args.blend / (double)num_neighbors);
+ Eigen::Quaternion<double> neigh_qt =
+ merged_rotations[TNeigh(tet_i, n_i)];
+ neigh_qt = quat_angle_multiply(neigh_qt,
+ args.blend / (double)num_neighbors);
qt = neigh_qt * qt;
}
}
@@ 392,14 399,16 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &d
} else {
eff_R.resize(Rdim, num_rots * Rdim);
for (int r = 0; r < num_rots; r++) {
- eff_R.block(0, Rdim * r, Rdim, Rdim) = Rotations.block(0, Rdim * data.G(r), Rdim, Rdim);
+ eff_R.block(0, Rdim * r, Rdim, Rdim) =
+ Rotations.block(0, Rdim * data.G(r), Rdim, Rdim);
}
}
MatrixXd Dl;
if (data.with_dynamics) {
- assert(data.M.rows() == n &&
- "No mass matrix. Call arap_precomputation if changing with_dynamics");
+ assert(
+ data.M.rows() == n &&
+ "No mass matrix. Call arap_precomputation if changing with_dynamics");
const double h = data.h;
assert(h != 0);
// Dl = 1./(h*h*h)*M*(-2.*V0 + Vm1) - fext;
@@ 414,8 423,8 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &d
VectorXd Rcol;
igl::columnize(eff_R, num_rots, 2, Rcol);
VectorXd Bcol = -data.K * Rcol;
- // `Bcol` is just the sum of edges in a tet, rotated and weighted by cot and length.
- // `Bc` becomes the `c` component.
+ // `Bcol` is just the sum of edges in a tet, rotated and weighted by cot and
+ // length. `Bc` becomes the `c` component.
assert(Bcol.size() == data.n * data.dim);
for (int c = 0; c < data.dim; c++) {
VectorXd Uc, Bc, bcc, Beq;
@@ 426,7 435,8 @@ void arap::iter(const arap::Args &args, const Eigen::MatrixXd &bc, arap::Data &d
if (bc.size() > 0) {
bcc = bc.col(c);
}
- bool ok = igl::min_quad_with_fixed_solve(data.solver_data, Bc, bcc, Beq, Uc);
+ bool ok =
+ igl::min_quad_with_fixed_solve(data.solver_data, Bc, bcc, Beq, Uc);
if (!ok) {
fprintf(stderr, "`igl::min_quad_with_fixed_solve` returned err!!\n");
}
M arap-mold.h => arap-mold.h +46 -22
@@ 13,15 13,34 @@ typedef Eigen::MatrixXi Derivedb;
typedef Eigen::MatrixXd Derivedbc;
typedef Eigen::MatrixXd DerivedU;
+namespace Moldable {
+enum Enum {
+ Yes,
+ No,
+ Crossing,
+};
+}
+
+namespace Orientation {
+enum Enum {
+ // This is in the direction of the normal
+ Above,
+ // This is the other direction.
+ Below,
+ Crossing
+};
+}
+
struct Data {
// n #V
int n;
- // G #V list of group indices (1 to k) for each vertex, such that vertex i is assigned to group
- // G(i)
+ // G #V list of group indices (1 to k) for each vertex, such that vertex i is
+ // assigned to group G(i)
Eigen::VectorXi G;
// energy type of energy to use
igl::ARAPEnergyType energy;
- // with_dynamics whether using dynamics (need to call arap_precomputation after changing)
+ // with_dynamics whether using dynamics (need to call arap_precomputation
+ // after changing)
bool with_dynamics;
// f_ext #V by dim list of external forces
Eigen::MatrixXd f_ext;
@@ 45,8 64,17 @@ struct Data {
// dim dimension being used for solving
int dim;
- // All triangle faces
- Eigen::MatrixXi mht_triangles;
+ // Boundary faces
+ Eigen::MatrixXi *faces;
+ // Face normals
+ Eigen::MatrixXd *normals;
+ // Is this face moldable?
+ std::vector<Moldable::Enum> *moldable;
+ // Whether the face is above or below the plane
+ std::vector<Orientation::Enum> *orientation;
+ // Is this point above or below the plane?
+ std::vector<bool> *vert_is_above;
+
// Is the ith face facing outwards?
Eigen::Vector3d *mht_plane_p;
Eigen::Vector3d *mht_plane_n;
@@ 56,14 84,12 @@ struct Data {
// XXX(debug): Register rotation angles for all faces
std::vector<double> mht_face_mold_angle;
std::vector<double> mht_face_arap_angle;
- // XXX(todo): remove this, and rather have `main` make color from the `dot`.
- std::vector<int> mht_face_classes;
// if `with_dynamics` is on, we need the initial `U0` across ARAP steps.
Eigen::MatrixXd U0;
Data()
- : n(0), G(), energy(igl::ARAP_ENERGY_TYPE_DEFAULT), with_dynamics(false), f_ext(), h(1),
- ym(1), max_iter(10), K(), CSM(), solver_data(), b(),
+ : n(0), G(), energy(igl::ARAP_ENERGY_TYPE_DEFAULT), with_dynamics(false),
+ f_ext(), h(1), ym(1), max_iter(10), K(), CSM(), solver_data(), b(),
dim(-1) // force this to be set by _precomputation
{};
};
@@ 74,7 100,9 @@ struct Args {
bool slerp_2_rots; /* SLERP rotations to join 2. Otherwise multiply by 0.5 */
bool rotate_before; /* Rotate normal before checking moldability */
- Args() : should_blend(false), blend(0.0), slerp_2_rots(false), rotate_before(false){};
+ Args()
+ : should_blend(false), blend(0.0), slerp_2_rots(false),
+ rotate_before(false){};
};
// Compute necessary information to start using an ARAP deformation
@@ 87,23 115,19 @@ struct Args {
// b #b list of "boundary" fixed vertex indices into V
// Outputs:
// data struct containing necessary precomputation
-bool precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F, const int dim,
- const Eigen::MatrixXi &b, Data &data);
+bool precomputation(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F,
+ const int dim, const Eigen::MatrixXi &b, Data &data);
// Inputs:
// bc #b by dim list of boundary conditions
// data struct containing necessary precomputation and parameters
// U #V by dim initial guess
-void iter(const Args &, const Eigen::MatrixXd &bc, Data &data, int num_tets, Eigen::MatrixXi &,
- Eigen::MatrixXd &U);
-
-enum Moldable {
- YES,
- NO,
- CROSSING,
-};
+void iter(const Args &, const Eigen::MatrixXd &bc, Data &data, int num_tets,
+ Eigen::MatrixXi &, Eigen::MatrixXd &U);
-Moldable is_moldable(Eigen::Vector3d &remove_dir, std::vector<bool> &is_point_above,
- Eigen::MatrixXi &faces, int face, Eigen::Vector3d normal, double *out_dot);
+Moldable::Enum is_moldable(Eigen::Vector3d &remove_dir,
+ std::vector<bool> &is_point_above,
+ Eigen::MatrixXi &faces, int face,
+ Eigen::Vector3d normal, double *out_dot);
} // namespace arap
M main.cpp => main.cpp +164 -98
@@ 16,6 16,11 @@
using igl::opengl::glfw::Viewer;
+using Eigen::MatrixXd;
+using Eigen::MatrixXi;
+using Eigen::Vector3d;
+using Eigen::VectorXi;
+
void step_arap();
bool pre_draw(Viewer &);
void draw_plane(Viewer &);
@@ 24,28 29,37 @@ void draw_menu();
void update_vectors();
void ensure_moldability();
void remesh();
-int ray_intersect(Eigen::Vector3d, Eigen::Vector3d, Eigen::MatrixXd &, Eigen::MatrixXi &);
-void compute_forces(Eigen::MatrixXd &, Eigen::MatrixXi &);
+int ray_intersect(Vector3d, Vector3d, MatrixXd &, MatrixXi &);
+void compute_forces(MatrixXd &, MatrixXi &);
+void compute_things(MatrixXd &, MatrixXi &, bool = false);
// Model data
-Eigen::MatrixXd V;
-Eigen::MatrixXi F;
-Eigen::MatrixXd N;
-Eigen::MatrixXd colors;
+MatrixXd V;
+MatrixXi F;
+MatrixXd colors;
// Tetrahedralized data
-Eigen::MatrixXd TV; // Vertices
-Eigen::MatrixXi TT; // Tets
-Eigen::MatrixXi TF; // Triangles
-Eigen::MatrixXi TF2T; // Face to tets mapping
-Eigen::MatrixXi TT2F; // Tet to faces mapping
-Eigen::MatrixXi TNeigh; // Neighbor tets for each tet
+MatrixXd TV; // Vertices
+MatrixXi TT; // Tets
+MatrixXi TF; // Triangles
+MatrixXd TN; // Normals
+MatrixXi TF2T; // Face to tets mapping
+MatrixXi TT2F; // Tet to faces mapping
+MatrixXi TNeigh; // Neighbor tets for each tet
// Indices of the fixed vertices in TV.
-Eigen::VectorXi boundary(BOUNDARY_SIZE);
+VectorXi boundary(BOUNDARY_SIZE);
// ARAP'd vertices
-Eigen::MatrixXd U; /* Deformed vertices */
-Eigen::MatrixXd U0; /* Original U (this is really just TV) */
+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 */
arap::Data arap_data;
arap::Args arap_args;
@@ 77,9 91,12 @@ void update_vectors() {
// 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");
@@ 90,6 107,7 @@ void draw_menu() {
plane_n = fplane_n.cast<double>().normalized();
update_vectors();
draw_plane(viewer);
+ compute_things(U, TF, true);
}
if (!run_arap) {
@@ 100,14 118,18 @@ void draw_menu() {
run_arap = false;
}
- if (ImGui::Button("Reset"))
+ if (ImGui::Button("Reset")) {
U = TV;
+ compute_things(U, TF);
+ }
- 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("Ensure Moldability"))
+ // ensure_moldability();
if (ImGui::Button("Remesh"))
remesh();
@@ 142,7 164,8 @@ 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);
@@ 157,22 180,20 @@ 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);
}
return false;
}
// Perform one ARAP step.
void step_arap() {
- static int count = 0;
- count++;
- printf("arap #%d\r", count);
- fflush(stdout);
Eigen::MatrixXd bc(BOUNDARY_SIZE, 3);
for (int i = 0; i < BOUNDARY_SIZE; i++)
bc.row(i) = TV.row(boundary(i));
- fill(arap_data.mht_face_classes.begin(), arap_data.mht_face_classes.end(), 0);
+
// 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++) {
@@ 180,36 201,84 @@ void step_arap() {
// // forces.row(i) = (U.row(0) - U.row(i)) * 0.1;
// }
// arap_data.f_ext = forces;
- Eigen::MatrixXd U_prev = U;
+
arap::iter(arap_args, bc, arap_data, TT.rows(), TNeigh, U);
}
+// 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(Eigen::MatrixXd &vertices, Eigen::MatrixXi &faces,
+ bool only_plane) {
+
+ if (!only_plane) { /* `normals`: depends on vertices */
+ igl::per_face_normals(vertices, faces, TN);
+ viewer.data().set_vertices(vertices);
+ viewer.data().set_normals(TN);
+ }
+
+ { /* `is_point_above`: depends on vertices and the plane */
+ is_point_above.resize(vertices.rows());
+ for (int v_i = 0; v_i < vertices.rows(); v_i++) {
+ Vector3d v(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(faces.rows());
+ orientation.resize(faces.rows());
+ dots.resize(faces.rows());
+ 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);
+ 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 */
+ colors.resize(faces.rows(), 3);
+ for (int f = 0; f < faces.rows(); f++) {
+ switch (moldable[f]) {
+ case arap::Moldable::Yes: {
+ colors.row(f) = Eigen::Array3d(1.0, 1.0, 1.0);
+ } break;
+ case arap::Moldable::Crossing: {
+ colors.row(f) = Eigen::Array3d(0.0, 0.0, 0.0);
+ } break;
+ default: {
+ double r = -dots[f];
+ r = sqrt(r);
+ colors.row(f) = Eigen::Array3d(1.0, (1.0 - r), (1.0 - r));
+ } break;
+ }
+ }
+ viewer.data().set_colors(colors);
+ }
+}
+
// 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
// the rest is pretty constant.
bool pre_draw(Viewer &viewer) {
if (run_arap) {
step_arap();
+ compute_things(U, TF);
}
- viewer.data().set_vertices(U);
- viewer.data().compute_normals();
- for (ssize_t i = 0; i < TF.rows(); i++) {
- int c = arap_data.mht_face_classes[i];
- switch (c) {
- case 0: {
- colors.row(i) = Eigen::Array3d(1.0, 1.0, 1.0);
- } break;
- case -1: {
- colors.row(i) = Eigen::Array3d(0.0, 0.0, 0.0);
- } break;
- default: {
- double r = (double)(c - 1) / 100.0;
- r = sqrt(r);
- colors.row(i) = Eigen::Array3d(1.0, (1.0 - r), (1.0 - r));
- } break;
- }
- }
- viewer.data().set_colors(colors);
return false;
}
@@ 273,10 342,11 @@ void draw_plane(Viewer &viewer) {
// Deform the object so that moldability is ensured.
void ensure_moldability() {
- // We'll do this by finding all faces that are not moldable, and change their positions in `U`
- // such that their faces would be moldable. This is weird, because altering the position of a
- // vertex which is adjacent to a non-moldable face may make another face adjacent to that vertex,
- // which was previsously moldable, non-moldable. It'll be something like this:
+ // We'll do this by finding all faces that are not moldable, and change their
+ // positions in `U` such that their faces would be moldable. This is weird,
+ // because altering the position of a vertex which is adjacent to a
+ // non-moldable face may make another face adjacent to that vertex, which was
+ // previsously moldable, non-moldable. It'll be something like this:
//
// q = make_queue()>
// for face in faces {
@@ 294,12 364,14 @@ void ensure_moldability() {
// }
//
Eigen::MatrixXd normals;
- // XXX(todo): move this to after all places in which U is modified. We pretty much
- // always need to have this up to date.
+ // XXX(todo): move this to after all places in which U is modified. We pretty
+ // much always need to have this up to date.
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));
}
@@ 363,23 435,25 @@ 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);
- if (moldable == arap::NO) {
+ auto moldable =
+ arap::is_moldable(removal_dir, is_point_above, TF, i, normal, NULL);
+ if (moldable == arap::Moldable::No) {
indices.push_back(i);
}
}
while (indices.size() > 0) {
- // We have to drag all points towards the point that's highest along `plane_n` in order to make
- // it moldable.
+ // We have to drag all points towards the point that's highest along
+ // `plane_n` in order to make it moldable.
int face = indices.back();
indices.pop_back();
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::NO)
+ if (moldable != arap::Moldable::No)
continue;
auto drag_dir = (normal - plane_n.dot(normal) * plane_n).normalized();
{ /* XXX(debug): NaN check */
@@ 397,54 471,44 @@ void ensure_moldability() {
// auto top = highest_vertex(face);
auto far = furthest_vertex(drag_dir, face);
- (void) far;
+ (void)far;
- // XXX(todo): check neighboring faces of these vertices and add to `indices` as needed.
+ // XXX(todo): check neighboring faces of these vertices and add to `indices`
+ // as needed.
}
-
- igl::per_face_normals(TV, TF, N);
}
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.
- std::vector<bool> is_point_above(vertices.rows());
- { /* refactorable */
- auto above_plane = [&](Eigen::Vector3d v) {
- return (v - plane_p).dot(plane_n) >= 0.0;
- };
- for (int i = 0; i < vertices.rows(); i++)
- is_point_above[i] = above_plane(vertices.row(i));
- }
-
for (int f = 0; f < TF.rows(); f++) {
- auto normal = N.row(f);
- double dot;
- // XXX(refactor): this is done quite a few times all over the place now.
- auto moldable = arap::is_moldable(remove_dir, is_point_above, TF, f, normal, &dot);
- // XXX(todo): what do we do about `CROSSING`? For now, skip as well.
- if (moldable != arap::NO) continue;
+ auto m = moldable[f];
+ // XXX(todo): what do we do about `Crossing`? For now, skip as well.
+ if (m != arap::Moldable::No)
+ continue;
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.
+ // 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);
- if (other_i == -1) continue;
+ if (other_i == -1)
+ continue;
auto 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`.
+ // 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 (?)
}
}
}
-// 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) {
+// 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
assert(false);
return -1;
@@ 500,9 564,7 @@ int main(int argc, char *argv[]) {
TF2T = tetmesh_out.face2tets;
TT2F = tetmesh_out.tet2faces;
TNeigh = tetmesh_out.neighbors;
-
- igl::per_face_normals(TV, TF, N);
- colors = Eigen::MatrixXd(N.rows(), 3);
+ U = TV;
boundary << 0, 1, 2, 3; // 19, 22, 121, 122, 221, 224, 321, 326;
@@ 515,19 577,21 @@ int main(int argc, char *argv[]) {
// Dbg(TF2T);
// Dbg(TNeigh);
- arap_data.mht_triangles = TF;
+ arap_data.faces = &TF;
+ arap_data.normals = &TN;
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.mht_face_classes = std::vector<int>(TF.size(), 0);
+ 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;
- U = TV;
+ arap_data.vert_is_above = &is_point_above;
+ arap_data.moldable = &moldable;
+ arap_data.orientation = &orientation;
maxs = V.colwise().maxCoeff().cast<float>();
mins = V.colwise().minCoeff().cast<float>();
@@ 536,7 600,6 @@ int main(int argc, char *argv[]) {
fplane_p = ((maxs + mins) / 2.0);
fplane_n = Eigen::Vector3f(0.0, 0.0, 1.0);
update_vectors();
- draw_plane(viewer);
igl::opengl::glfw::imgui::ImGuiMenu menu;
viewer.plugins.push_back(&menu);
@@ 550,6 613,9 @@ int main(int argc, char *argv[]) {
viewer.data().set_mesh(TV, TF);
viewer.data().set_face_based(true);
viewer.core.is_animating = true;
+
draw_plane(viewer);
+ compute_things(TV, TF);
+
viewer.launch();
}
M mesh.cpp => mesh.cpp +5 -4
@@ 57,7 57,8 @@ int make_tetmesh(TetmeshIn &in, TetmeshOut &out) {
}
::tetrahedralize(&tb, &tetgen_in, &tetgen_out);
} catch (int e) {
- cerr << "[ERR]: from `tetgen` in " << __FUNCTION__ << ":" << __LINE__ << ": " << e << endl;
+ cerr << "[ERR]: from `tetgen` in " << __FUNCTION__ << ":" << __LINE__
+ << ": " << e << endl;
return 1;
}
@@ 137,9 138,9 @@ int make_tetmesh(TetmeshIn &in, TetmeshOut &out) {
out.face2tets(i, 1) = tetgen_out.face2tetlist[2 * i + 1];
}
- // XXX(note): we don't really need this, and this requires the `f` flag to tetgen, which again
- // outputs all faces, and not just the boundary faces (which really are the only we care about,
- // for visualization).
+ // XXX(note): we don't really need this, and this requires the `f` flag to
+ // tetgen, which again outputs all faces, and not just the boundary faces
+ // (which really are the only we care about, for visualization).
return 0;
// T2F
if (tetgen_out.tet2facelist == NULL) {