From edb07f1d8e136fba9f06334dc08e0ee2fd47884b Mon Sep 17 00:00:00 2001 From: Yichao Zhou Date: Wed, 7 Feb 2018 23:42:06 -0800 Subject: [PATCH] Compilable postprocessing --- .clang-format | 2 - src/parametrizer.cpp | 357 ++++++++++++------------------------- src/post-solver.cpp | 110 +++++++++--- src/post-solver.hpp | 2 +- src/subdivide.cpp | 416 ++++++++++++++++++++----------------------- 5 files changed, 394 insertions(+), 493 deletions(-) diff --git a/.clang-format b/.clang-format index 6e2176204..3f4066876 100644 --- a/.clang-format +++ b/.clang-format @@ -3,5 +3,3 @@ BasedOnStyle: Google IndentWidth: 4 Standard: Cpp11 ColumnLimit: 99 -BinPackArguments: false -BinPackParameters: false diff --git a/src/parametrizer.cpp b/src/parametrizer.cpp index 85c6dab3c..46ea7f83a 100644 --- a/src/parametrizer.cpp +++ b/src/parametrizer.cpp @@ -16,10 +16,8 @@ #define LOG_OUTPUT //#define PERFORM_TEST -extern void generate_adjacency_matrix_uniform(const MatrixXi& F, - const VectorXi& V2E, - const VectorXi& E2E, - const VectorXi& nonManifold, +extern void generate_adjacency_matrix_uniform(const MatrixXi& F, const VectorXi& V2E, + const VectorXi& E2E, const VectorXi& nonManifold, AdjacentMatrix& adj); inline double fast_acos(double x) { @@ -318,23 +316,9 @@ void Parametrizer::ComputePositionSingularities(int with_scale) { } double inv_scale_x = 1.0 / scale_x, inv_scale_y = 1.0 / scale_y, inv_scale_x_1 = 1.0 / scale_x_1, inv_scale_y_1 = 1.0 / scale_y_1; - std::pair value = compat_position_extrinsic_index_4(v[k], - n[k], - q[k], - o[k], - v[kn], - n[kn], - q[kn], - o[kn], - scale_x, - scale_y, - inv_scale_x, - inv_scale_y, - scale_x_1, - scale_y_1, - inv_scale_x_1, - inv_scale_y_1, - nullptr); + std::pair value = compat_position_extrinsic_index_4( + v[k], n[k], q[k], o[k], v[kn], n[kn], q[kn], o[kn], scale_x, scale_y, inv_scale_x, + inv_scale_y, scale_x_1, scale_y_1, inv_scale_x_1, inv_scale_y_1, nullptr); auto diff = value.first - value.second; index += diff; pos_index(k * 2, f) = diff[0]; @@ -389,75 +373,23 @@ void Parametrizer::EstimateScale() { f = i; len = step; - TravelField(p, - q_xl, - len, - f, - hierarchy.mE2E, - mV, - mF, - Nf, - FQ, - mQ, - mN, - triangle_space, - &tx, - &ty, - &q_yl_unfold); + TravelField(p, q_xl, len, f, hierarchy.mE2E, mV, mF, Nf, FQ, mQ, mN, triangle_space, &tx, + &ty, &q_yl_unfold); f = i; len = step; - TravelField(p, - q_xr, - len, - f, - hierarchy.mE2E, - mV, - mF, - Nf, - FQ, - mQ, - mN, - triangle_space, - &tx, - &ty, - &q_yr_unfold); + TravelField(p, q_xr, len, f, hierarchy.mE2E, mV, mF, Nf, FQ, mQ, mN, triangle_space, &tx, + &ty, &q_yr_unfold); f = i; len = step; - TravelField(p, - q_yl, - len, - f, - hierarchy.mE2E, - mV, - mF, - Nf, - FQ, - mQ, - mN, - triangle_space, - &tx, - &ty, - &q_xl_unfold); + TravelField(p, q_yl, len, f, hierarchy.mE2E, mV, mF, Nf, FQ, mQ, mN, triangle_space, &tx, + &ty, &q_xl_unfold); f = i; len = step; - TravelField(p, - q_yr, - len, - f, - hierarchy.mE2E, - mV, - mF, - Nf, - FQ, - mQ, - mN, - triangle_space, - &tx, - &ty, - &q_xr_unfold); + TravelField(p, q_yr, len, f, hierarchy.mE2E, mV, mF, Nf, FQ, mQ, mN, triangle_space, &tx, + &ty, &q_xr_unfold); double dSx = (q_yr_unfold - q_yl_unfold).dot(q_x) / (2.0f * step); double dSy = (q_xr_unfold - q_xl_unfold).dot(q_y) / (2.0f * step); FS.col(i) = Vector2d(dSx, dSy); @@ -469,8 +401,8 @@ void Parametrizer::EstimateScale() { Vector3d p2 = mV.col(mF(2, i)) - mV.col(mF(0, i)); double area = p1.cross(p2).norm(); for (int j = 0; j < 3; ++j) { - auto index = compat_orientation_extrinsic_index_4( - FQ.col(i), Nf.col(i), mQ.col(mF(j, i)), mN.col(mF(j, i))); + auto index = compat_orientation_extrinsic_index_4(FQ.col(i), Nf.col(i), + mQ.col(mF(j, i)), mN.col(mF(j, i))); double scaleX = FS.col(i).x(), scaleY = FS.col(i).y(); if (index.first != index.second % 2) { std::swap(scaleX, scaleY); @@ -831,36 +763,12 @@ void Parametrizer::ComputeIndexMap(int with_scale) { printf("Fix flip advance...\n"); printf("thanks %ld %ld %zd\n", V.cols(), O.cols(), disajoint_tree.indices.size()); - subdivide_diff(F, - V, - N, - Q, - O, - V2E, - hierarchy.mE2E, - boundary, - nonManifold, - edge_diff, - edge_values, - face_edgeOrients, - face_edgeIds, - singularities); + subdivide_diff(F, V, N, Q, O, V2E, hierarchy.mE2E, boundary, nonManifold, edge_diff, + edge_values, face_edgeOrients, face_edgeIds, singularities); FixFlipAdvance(); - subdivide_diff(F, - V, - N, - Q, - O, - V2E, - hierarchy.mE2E, - boundary, - nonManifold, - edge_diff, - edge_values, - face_edgeOrients, - face_edgeIds, - singularities); + subdivide_diff(F, V, N, Q, O, V2E, hierarchy.mE2E, boundary, nonManifold, edge_diff, + edge_values, face_edgeOrients, face_edgeIds, singularities); for (int i = 0; i < edge_diff.size(); ++i) { if (abs(edge_diff[i][0]) > 1 || abs(edge_diff[i][1]) > 1) { @@ -892,8 +800,8 @@ void Parametrizer::ComputeIndexMap(int with_scale) { if (counter[compact_v] == 0) Q_compact[compact_v] = Q.col(i); else { - auto pairs = compat_orientation_extrinsic_4( - Q_compact[compact_v], N_compact[compact_v], Q.col(i), N.col(i)); + auto pairs = compat_orientation_extrinsic_4(Q_compact[compact_v], N_compact[compact_v], + Q.col(i), N.col(i)); Q_compact[compact_v] = (pairs.first * counter[compact_v] + pairs.second).normalized(); } counter[compact_v] += 1; @@ -1002,8 +910,8 @@ void Parametrizer::ComputeIndexMap(int with_scale) { for (auto& c : quad_cells) { if (c.second.second != Vector3i(-100, -100, -100)) { - F_compact.push_back(Vector4i( - c.second.first[0], c.second.second[2], c.second.first[1], c.second.first[2])); + F_compact.push_back(Vector4i(c.second.first[0], c.second.second[2], c.second.first[1], + c.second.first[2])); } } printf("Fix holes...\n"); @@ -1011,24 +919,12 @@ void Parametrizer::ComputeIndexMap(int with_scale) { // potential bug, not guarantee to have quads at holes! printf("Direct Quad Graph...\n"); - compute_direct_graph_quad( - O_compact, F_compact, V2E_compact, E2E_compact, boundary_compact, nonManifold_compact); + compute_direct_graph_quad(O_compact, F_compact, V2E_compact, E2E_compact, boundary_compact, + nonManifold_compact); printf("Direct Quad Graph finish...\n"); - optimize_quad_positions(O_compact, - N_compact, - Q_compact, - F_compact, - V2E_compact, - E2E_compact, - V, - N, - Q, - O, - F, - V2E, - hierarchy.mE2E, - disajoint_tree); + optimize_quad_positions(O_compact, N_compact, Q_compact, F_compact, V2E_compact, E2E_compact, + V, N, Q, O, F, V2E, hierarchy.mE2E, disajoint_tree, hierarchy.mScale); int count = 0; for (int i = 0; i < F.cols(); ++i) { @@ -1149,11 +1045,11 @@ void Parametrizer::FixHoles() { while (!loop_vertices.empty()) { if (loop_vertices.size() <= 4) { if (loop_vertices.size() == 4) - F_compact.push_back(Vector4i( - loop_vertices[0], loop_vertices[1], loop_vertices[2], loop_vertices[3])); + F_compact.push_back(Vector4i(loop_vertices[0], loop_vertices[1], + loop_vertices[2], loop_vertices[3])); else - F_compact.push_back(Vector4i( - loop_vertices[0], loop_vertices[1], loop_vertices[2], loop_vertices[2])); + F_compact.push_back(Vector4i(loop_vertices[0], loop_vertices[1], + loop_vertices[2], loop_vertices[2])); if (directed_edges.count((long long)loop_vertices[0] * numV + loop_vertices[1])) { std::swap(F_compact.back()[1], F_compact.back()[3]); @@ -1600,12 +1496,8 @@ void Parametrizer::FixFlipAdvance() { auto& l2 = vertices_to_edges[ny][nx]; if (std::find(l1.begin(), l1.end(), i) == l1.end() || std::find(l2.begin(), l2.end(), i) == l2.end()) { - printf("edge %d not indexed in vertices (%d %d) %d %d\n", - i, - nx, - ny, - edge_values[i].x, - edge_values[i].y); + printf("edge %d not indexed in vertices (%d %d) %d %d\n", i, nx, ny, + edge_values[i].x, edge_values[i].y); for (auto& m : l1) { printf("<%d %d> ", edge_values[m].x, edge_values[m].y); } @@ -1668,8 +1560,7 @@ void Parametrizer::FixFlipAdvance() { } if (l.size() != faces_from_edge[f].size()) { - printf("inconsistent edge-face connection! -1 %d %d\n", - (int)l.size(), + printf("inconsistent edge-face connection! -1 %d %d\n", (int)l.size(), (int)faces_from_edge[f].size()); for (auto& p : l) { printf("%d ", p); @@ -1679,10 +1570,7 @@ void Parametrizer::FixFlipAdvance() { printf("%d ", faces_from_edge[f][i]); } printf("\n"); - printf("face %d %d %d %d\n", - f, - tree.Parent(F(0, f)), - tree.Parent(F(1, f)), + printf("face %d %d %d %d\n", f, tree.Parent(F(0, f)), tree.Parent(F(1, f)), tree.Parent(F(2, f))); printf("face origin %d %d %d\n", F(0, f), F(1, f), F(2, f)); // exit(0); @@ -1713,25 +1601,11 @@ void Parametrizer::FixFlipAdvance() { Vector2i total = diff[0] + diff[1] + diff[2]; if (total != Vector2i::Zero()) { printf("zero face constraint violated %d\n", i); - printf("<%d %d> (%d eid %d) <%d %d> (%d eid %d) <%d %d> (%d eid %d)\n", - diff[0][0], - diff[0][1], - orients[0], - pids[0], - diff[1][0], - diff[1][1], - orients[1], - pids[1], - diff[2][0], - diff[2][1], - orients[2], - pids[2]); - printf("f %d (%d %d %d): %d %d %d", - i, - tree.Parent(F(0, i)), - tree.Parent(F(1, i)), - tree.Parent(F(2, i)), - get_parents(parent_edge, face_edgeIds[i][0]), + printf("<%d %d> (%d eid %d) <%d %d> (%d eid %d) <%d %d> (%d eid %d)\n", diff[0][0], + diff[0][1], orients[0], pids[0], diff[1][0], diff[1][1], orients[1], + pids[1], diff[2][0], diff[2][1], orients[2], pids[2]); + printf("f %d (%d %d %d): %d %d %d", i, tree.Parent(F(0, i)), tree.Parent(F(1, i)), + tree.Parent(F(2, i)), get_parents(parent_edge, face_edgeIds[i][0]), get_parents(parent_edge, face_edgeIds[i][1]), get_parents(parent_edge, face_edgeIds[i][2])); // exit(0); @@ -1743,77 +1617,76 @@ void Parametrizer::FixFlipAdvance() { printf("finish...\n"); }; - auto ExtractEdgeSet = - [&](int v1, int v2, int pid, std::vector>& edge_change) { - std::unordered_map edge_set; - edge_change.push_back(std::make_pair(pid, edge_diff[pid])); - edge_set.insert(edge_change.back()); - std::queue faces; - for (auto& f : edge_to_faces[pid]) { - faces.push(f); + auto ExtractEdgeSet = [&](int v1, int v2, int pid, + std::vector>& edge_change) { + std::unordered_map edge_set; + edge_change.push_back(std::make_pair(pid, edge_diff[pid])); + edge_set.insert(edge_change.back()); + std::queue faces; + for (auto& f : edge_to_faces[pid]) { + faces.push(f); + } + std::set modified; + while (!faces.empty()) { + int f = faces.front(); + modified.insert(f); + faces.pop(); + int eids[3]; + int orient[3]; + Vector2i total_diff(0, 0); + for (int i = 0; i < 3; ++i) { + int eid = face_edgeIds[f][i]; + int pid = get_parents(parent_edge, eid); + orient[i] = (get_parents_orient(parent_edge, eid) + face_edgeOrients[f][i]) % 4; + eids[i] = pid; + Vector2i diff = edge_diff[pid]; + if (edge_set.count(pid)) diff -= edge_set[pid]; + total_diff += rshift90(diff, orient[i]); } - std::set modified; - while (!faces.empty()) { - int f = faces.front(); - modified.insert(f); - faces.pop(); - int eids[3]; - int orient[3]; - Vector2i total_diff(0, 0); - for (int i = 0; i < 3; ++i) { - int eid = face_edgeIds[f][i]; - int pid = get_parents(parent_edge, eid); - orient[i] = - (get_parents_orient(parent_edge, eid) + face_edgeOrients[f][i]) % 4; - eids[i] = pid; - Vector2i diff = edge_diff[pid]; - if (edge_set.count(pid)) diff -= edge_set[pid]; - total_diff += rshift90(diff, orient[i]); - } - int count = 0; - for (int i = 0; i < 3; ++i) { - if (!((tree.Parent(edge_values[eids[i]].x) != v1 && - tree.Parent(edge_values[eids[i]].y) != v1) || - edge_set.count(eids[i]))) { - count += 1; - } - } - int next_e = 0; - while ((tree.Parent(edge_values[eids[next_e]].x) != v1 && - tree.Parent(edge_values[eids[next_e]].y) != v1) || - edge_set.count(eids[next_e])) { - next_e += 1; - if (next_e == 3) break; - } - if (total_diff == Vector2i::Zero()) { - continue; - } - if (next_e == 3) { - edge_change.clear(); - return; - } - int e = next_e + 1; - while (e < 3 && eids[next_e] != eids[e]) e += 1; - if (e != 3) { - edge_change.clear(); - return; - } - int change_pid = eids[next_e]; - auto new_diff = rshift90(total_diff, (4 - orient[next_e]) % 4); - if (abs(edge_diff[change_pid][0] - new_diff[0]) > edge_len || - abs(edge_diff[change_pid][1] - new_diff[1]) > edge_len) { - edge_change.clear(); - return; + int count = 0; + for (int i = 0; i < 3; ++i) { + if (!((tree.Parent(edge_values[eids[i]].x) != v1 && + tree.Parent(edge_values[eids[i]].y) != v1) || + edge_set.count(eids[i]))) { + count += 1; } - edge_change.push_back(std::make_pair(change_pid, new_diff)); - edge_set.insert(edge_change.back()); - for (auto& nf : edge_to_faces[change_pid]) { - if (nf != f) { - faces.push(nf); - } + } + int next_e = 0; + while ((tree.Parent(edge_values[eids[next_e]].x) != v1 && + tree.Parent(edge_values[eids[next_e]].y) != v1) || + edge_set.count(eids[next_e])) { + next_e += 1; + if (next_e == 3) break; + } + if (total_diff == Vector2i::Zero()) { + continue; + } + if (next_e == 3) { + edge_change.clear(); + return; + } + int e = next_e + 1; + while (e < 3 && eids[next_e] != eids[e]) e += 1; + if (e != 3) { + edge_change.clear(); + return; + } + int change_pid = eids[next_e]; + auto new_diff = rshift90(total_diff, (4 - orient[next_e]) % 4); + if (abs(edge_diff[change_pid][0] - new_diff[0]) > edge_len || + abs(edge_diff[change_pid][1] - new_diff[1]) > edge_len) { + edge_change.clear(); + return; + } + edge_change.push_back(std::make_pair(change_pid, new_diff)); + edge_set.insert(edge_change.back()); + for (auto& nf : edge_to_faces[change_pid]) { + if (nf != f) { + faces.push(nf); } } - }; + } + }; auto collapse = [&](int v1, int v2) { if (v1 == v2) return; std::set collapsed_faces; @@ -1892,20 +1765,12 @@ void Parametrizer::FixFlipAdvance() { while (orient < 4 && rshift90(diff1, orient) != diff2) orient += 1; if (orient == 4) { - printf("v %d %d %d %d\n", - F(j, f), - F((j + 1) % 3, f), - F(nj, nf), - F((nj + 1) % 3, nf)); + printf("v %d %d %d %d\n", F(j, f), F((j + 1) % 3, f), + F(nj, nf), F((nj + 1) % 3, nf)); printf("edge fail to collapse %d %d %d %d\n", - edge_values[peid].x, - edge_values[peid].y, - edge_values[npeid].x, - edge_values[npeid].y); - printf("%d %d %d %d\n", - diff1[0], - diff1[1], - diff2[0], + edge_values[peid].x, edge_values[peid].y, + edge_values[npeid].x, edge_values[npeid].y); + printf("%d %d %d %d\n", diff1[0], diff1[1], diff2[0], diff2[1]); printf("no orient solution!\n"); exit(0); diff --git a/src/post-solver.cpp b/src/post-solver.cpp index 231f8cd60..3dba10d46 100644 --- a/src/post-solver.cpp +++ b/src/post-solver.cpp @@ -9,60 +9,76 @@ #include "ceres/ceres.h" #include "ceres/rotation.h" +template +T Dot(const T a[3], const T2 b[3]) { + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; +} + +template +T Length2(const T a[3]) { + return Dot(a, a); +} + struct FaceConstraint { - FaceConstraint(double alpha, double beta, Vector3d normal[4], double bias, double length) + FaceConstraint(double alpha, double beta, Vector3d normal[4], double bias[4], double length) : alpha(alpha), beta(beta), - normal0{normal[0], normal[1], normal[2], normal[3]}, - bias0(bias), - length0(length) {} - - template - T dot(const T a[3], const T b[3]) { - return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; - } + length0(length), + bias0{ + bias[0], + bias[1], + bias[2], + bias[3], + }, + normal0{ + normal[0], + normal[1], + normal[2], + normal[3], + } {} template - T length2(const T a[3]) { - return dot(a, a); - } - - template - bool operator()(const T* p0, T* residuals) const { + bool operator()(const T* p0, T* r) const { const T* p[] = {p0, p0 + 3, p0 + 6, p0 + 9}; for (int k = 0; k < 4; ++k) { auto pc = p[k]; auto pa = p[(k + 1) % 4]; - auto pb = p[(k + 2) % 4]; + auto pb = p[(k + 3) % 4]; T a[3]{pa[0] - pc[0], pa[1] - pc[1], pa[2] - pc[2]}; T b[3]{pb[0] - pc[0], pb[1] - pc[1], pb[2] - pc[2]}; + r[3 * k + 0] = alpha * (ceres::sqrt(Length2(a)) - length0); + T normal[3]; ceres::CrossProduct(a, b, normal); - // length2(normal); - // T l2normal = ceres::sqrt(); + T l2normal = ceres::sqrt(Length2(normal)); + if (l2normal == T()) continue; + for (int i = 0; i < 3; ++i) normal[i] /= l2normal; + r[3 * k + 1] = ceres::acos(Dot(normal, &normal0[k][0])); + r[3 * k + 2] = beta * (Dot(pc, normal) - bias0[k]); } return true; } - static ceres::CostFunction* create(double alpha, double beta, Vector3d normal[4], double bias, - double length) { - return new ceres::AutoDiffCostFunction( + static ceres::CostFunction* create(double alpha, double beta, Vector3d normal[4], + double bias[4], double length) { + return new ceres::AutoDiffCostFunction( new FaceConstraint(alpha, beta, normal, bias, length)); } double alpha; double beta; - Vector3d normal0[4]; - double bias0; double length0; + + double bias0[4]; + Vector3d normal0[4]; }; void optimize_quad_positions(std::vector& O_quad, std::vector& N_quad, std::vector& Q_quad, std::vector& F_quad, VectorXi& V2E_quad, VectorXi& E2E_quad, MatrixXd& V, MatrixXd& N, MatrixXd& Q, MatrixXd& O, MatrixXi& F, VectorXi& V2E, VectorXi& E2E, - DisajointTree& disajoint_tree) { + DisajointTree& disajoint_tree, double reference_length) { // Information for the quad mesh printf("Quad mesh info:\n"); printf("Number of vertices with normals and orientations: %d = %d = %d\n", (int)O_quad.size(), @@ -105,4 +121,50 @@ void optimize_quad_positions(std::vector& O_quad, std::vector B_quad(n_quad); + std::vector B_weight(n_quad); + for (int vtrig = 0; vtrig < n_trig; ++vtrig) { + int vquad = disajoint_tree.Index(vtrig); + double b = N.col(vtrig).dot(O.col(vtrig)); + B_quad[vquad] += b; + B_weight[vquad] += 1; + } + for (int vquad = 0; vquad < n_quad; ++vquad) B_quad[vquad] /= B_weight[vquad]; + + ceres::Problem problem; + std::vector solution(n_quad * 3); + for (int vquad = 0; vquad < n_quad; ++vquad) { + solution[3 * vquad + 0] = O_quad[vquad][0]; + solution[3 * vquad + 1] = O_quad[vquad][1]; + solution[3 * vquad + 2] = O_quad[vquad][2]; + } + + for (int fquad = 0; fquad < F_quad.size(); ++fquad) { + auto v = F_quad[fquad]; + + double bias[4]; + Vector3d normal[4]; + std::vector var(12); + for (int k = 0; k < 4; ++k) { + var[3 * k + 0] = &solution[3 * v[k] + 0]; + var[3 * k + 1] = &solution[3 * v[k] + 1]; + var[3 * k + 2] = &solution[3 * v[k] + 2]; + bias[k] = B_quad[k]; + normal[k] = N_quad[k]; + } + ceres::CostFunction* cost_function = + FaceConstraint::create(0.1, 1, normal, bias, reference_length); + problem.AddResidualBlock(cost_function, nullptr, var); + } + + ceres::Solver::Options options; + options.linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY; + options.minimizer_progress_to_stdout = true; + ceres::Solver::Summary summary; + ceres::Solve(options, &problem, &summary); + + std::cout << summary.BriefReport() << std::endl; } diff --git a/src/post-solver.hpp b/src/post-solver.hpp index c8aa5d5e3..c25366a87 100644 --- a/src/post-solver.hpp +++ b/src/post-solver.hpp @@ -54,6 +54,6 @@ void optimize_quad_positions(std::vector& O_quad, std::vector& Q_quad, std::vector& F_quad, VectorXi& V2E_quad, VectorXi& E2E_quad, MatrixXd& V, MatrixXd& N, MatrixXd& Q, MatrixXd& O, MatrixXi& F, VectorXi& V2E, VectorXi& E2E, - DisajointTree& disajoint_tree); + DisajointTree& disajoint_tree, double reference_length); #endif /* post_solver_h */ diff --git a/src/subdivide.cpp b/src/subdivide.cpp index 796f28c6f..094488b88 100644 --- a/src/subdivide.cpp +++ b/src/subdivide.cpp @@ -1,177 +1,160 @@ #include "subdivide.hpp" -#include "dedge.hpp" -#include #include +#include +#include "dedge.hpp" + +void subdivide(MatrixXi &F, MatrixXd &V, VectorXi &V2E, VectorXi &E2E, VectorXi &boundary, + VectorXi &nonmanifold, double maxLength) { + typedef std::pair Edge; + + std::priority_queue queue; + + maxLength *= maxLength; + + for (int i = 0; i < E2E.size(); ++i) { + int v0 = F(i % 3, i / 3), v1 = F((i + 1) % 3, i / 3); + if (nonmanifold[v0] || nonmanifold[v1]) continue; + double length = (V.col(v0) - V.col(v1)).squaredNorm(); + if (length > maxLength) { + int other = E2E[i]; + if (other == -1 || other > i) queue.push(Edge(length, i)); + } + } + + int nV = V.cols(), nF = F.cols(), nSplit = 0; + /* + / v0 \ + v1p 1 | 0 v0p + \ v1 / + + / v0 \ + / 1 | 0 \ + v1p - vn - v0p + \ 2 | 3 / + \ v1 / + + f0: vn, v0p, v0 + f1: vn, v0, v1p + f2: vn, v1p, v1 + f3: vn, v1, v0p + */ + int counter = 0; + while (!queue.empty()) { + counter += 1; + Edge edge = queue.top(); + queue.pop(); + int e0 = edge.second, e1 = E2E[e0]; + bool is_boundary = e1 == -1; + int f0 = e0 / 3, f1 = is_boundary ? -1 : (e1 / 3); + int v0 = F(e0 % 3, f0), v0p = F((e0 + 2) % 3, f0), v1 = F((e0 + 1) % 3, f0); + if ((V.col(v0) - V.col(v1)).squaredNorm() != edge.first) { + continue; + } + int v1p = is_boundary ? -1 : F((e1 + 2) % 3, f1); + int vn = nV++; + nSplit++; + /* Update V */ + if (nV > V.cols()) { + V.conservativeResize(V.rows(), V.cols() * 2); + V2E.conservativeResize(V.cols()); + boundary.conservativeResize(V.cols()); + nonmanifold.conservativeResize(V.cols()); + } + + /* Update V */ + V.col(vn) = (V.col(v0) + V.col(v1)) * 0.5f; + nonmanifold[vn] = false; + boundary[vn] = is_boundary; + + /* Update F and E2E */ + int f2 = is_boundary ? -1 : (nF++); + int f3 = nF++; + if (nF > F.cols()) { + F.conservativeResize(F.rows(), std::max(nF, (int)F.cols() * 2)); + E2E.conservativeResize(F.cols() * 3); + } + + /* Update F */ + F.col(f0) << vn, v0p, v0; + if (!is_boundary) { + F.col(f1) << vn, v0, v1p; + F.col(f2) << vn, v1p, v1; + } + F.col(f3) << vn, v1, v0p; + + /* Update E2E */ + const int e0p = E2E[dedge_prev_3(e0)], e0n = E2E[dedge_next_3(e0)]; -void subdivide(MatrixXi &F, MatrixXd &V, VectorXi &V2E, VectorXi &E2E, - VectorXi &boundary, VectorXi &nonmanifold, double maxLength) { - typedef std::pair Edge; - - std::priority_queue queue; - - maxLength *= maxLength; - - for (int i = 0; i < E2E.size(); ++i) { - int v0 = F(i % 3, i / 3), v1 = F((i + 1) % 3, i / 3); - if (nonmanifold[v0] || nonmanifold[v1]) - continue; - double length = (V.col(v0) - V.col(v1)).squaredNorm(); - if (length > maxLength) { - int other = E2E[i]; - if (other == -1 || other > i) - queue.push(Edge(length, i)); - } - } - - int nV = V.cols(), nF = F.cols(), nSplit = 0; - /* - / v0 \ - v1p 1 | 0 v0p - \ v1 / - - / v0 \ - / 1 | 0 \ - v1p - vn - v0p - \ 2 | 3 / - \ v1 / - - f0: vn, v0p, v0 - f1: vn, v0, v1p - f2: vn, v1p, v1 - f3: vn, v1, v0p - */ - int counter = 0; - while (!queue.empty()) { - counter += 1; - Edge edge = queue.top(); - queue.pop(); - int e0 = edge.second, e1 = E2E[e0]; - bool is_boundary = e1 == -1; - int f0 = e0 / 3, f1 = is_boundary ? -1 : (e1 / 3); - int v0 = F(e0 % 3, f0), v0p = F((e0 + 2) % 3, f0), v1 = F((e0 + 1) % 3, f0); - if ((V.col(v0) - V.col(v1)).squaredNorm() != edge.first) { - continue; - } - int v1p = is_boundary ? -1 : F((e1 + 2) % 3, f1); - int vn = nV++; - nSplit++; - /* Update V */ - if (nV > V.cols()) { - V.conservativeResize(V.rows(), V.cols() * 2); - V2E.conservativeResize(V.cols()); - boundary.conservativeResize(V.cols()); - nonmanifold.conservativeResize(V.cols()); - } - - /* Update V */ - V.col(vn) = (V.col(v0) + V.col(v1)) * 0.5f; - nonmanifold[vn] = false; - boundary[vn] = is_boundary; - - /* Update F and E2E */ - int f2 = is_boundary ? -1 : (nF++); - int f3 = nF++; - if (nF > F.cols()) { - F.conservativeResize(F.rows(), std::max(nF, (int)F.cols() * 2)); - E2E.conservativeResize(F.cols() * 3); - } - - /* Update F */ - F.col(f0) << vn, v0p, v0; - if (!is_boundary) { - F.col(f1) << vn, v0, v1p; - F.col(f2) << vn, v1p, v1; - } - F.col(f3) << vn, v1, v0p; - - /* Update E2E */ - const int e0p = E2E[dedge_prev_3(e0)], - e0n = E2E[dedge_next_3(e0)]; - -#define sE2E(a, b) E2E[a] = b; if (b != -1) E2E[b] = a; - sE2E(3 * f0 + 0, 3 * f3 + 2); - sE2E(3 * f0 + 1, e0p); - sE2E(3 * f3 + 1, e0n); - if (is_boundary) { - sE2E(3 * f0 + 2, -1); - sE2E(3 * f3 + 0, -1); - } - else { - const int e1p = E2E[dedge_prev_3(e1)], - e1n = E2E[dedge_next_3(e1)]; - sE2E(3 * f0 + 2, 3 * f1 + 0); - sE2E(3 * f1 + 1, e1n); - sE2E(3 * f1 + 2, 3 * f2 + 0); - sE2E(3 * f2 + 1, e1p); - sE2E(3 * f2 + 2, 3 * f3 + 0); - } +#define sE2E(a, b) \ + E2E[a] = b; \ + if (b != -1) E2E[b] = a; + sE2E(3 * f0 + 0, 3 * f3 + 2); + sE2E(3 * f0 + 1, e0p); + sE2E(3 * f3 + 1, e0n); + if (is_boundary) { + sE2E(3 * f0 + 2, -1); + sE2E(3 * f3 + 0, -1); + } else { + const int e1p = E2E[dedge_prev_3(e1)], e1n = E2E[dedge_next_3(e1)]; + sE2E(3 * f0 + 2, 3 * f1 + 0); + sE2E(3 * f1 + 1, e1n); + sE2E(3 * f1 + 2, 3 * f2 + 0); + sE2E(3 * f2 + 1, e1p); + sE2E(3 * f2 + 2, 3 * f3 + 0); + } #undef sE2E - /* Update V2E */ - V2E[v0] = 3 * f0 + 2; - V2E[vn] = 3 * f0 + 0; - V2E[v1] = 3 * f3 + 1; - V2E[v0p] = 3 * f0 + 1; - if (!is_boundary) - V2E[v1p] = 3 * f1 + 2; - - auto schedule = [&](int f) { - for (int i = 0; i<3; ++i) { - double length = (V.col(F(i, f)) - V.col(F((i + 1) % 3, f))).squaredNorm(); - if (length > maxLength) - queue.push(Edge(length, f * 3 + i)); - } - }; - - schedule(f0); - if (!is_boundary) { - schedule(f2); - schedule(f1); - }; - schedule(f3); - } - F.conservativeResize(F.rows(), nF); - V.conservativeResize(V.rows(), nV); - V2E.conservativeResize(nV); - boundary.conservativeResize(nV); - nonmanifold.conservativeResize(nV); - E2E.conservativeResize(nF * 3); + /* Update V2E */ + V2E[v0] = 3 * f0 + 2; + V2E[vn] = 3 * f0 + 0; + V2E[v1] = 3 * f3 + 1; + V2E[v0p] = 3 * f0 + 1; + if (!is_boundary) V2E[v1p] = 3 * f1 + 2; + auto schedule = [&](int f) { + for (int i = 0; i < 3; ++i) { + double length = (V.col(F(i, f)) - V.col(F((i + 1) % 3, f))).squaredNorm(); + if (length > maxLength) queue.push(Edge(length, f * 3 + i)); + } + }; + schedule(f0); + if (!is_boundary) { + schedule(f2); + schedule(f1); + }; + schedule(f3); + } + F.conservativeResize(F.rows(), nF); + V.conservativeResize(V.rows(), nV); + V2E.conservativeResize(nV); + boundary.conservativeResize(nV); + nonmanifold.conservativeResize(nV); + E2E.conservativeResize(nF * 3); } -#include "parametrizer.hpp" #include "disajoint-tree.hpp" #include "field-math.hpp" +#include "parametrizer.hpp" -void subdivide_diff(MatrixXi &F, - MatrixXd &V, - MatrixXd &N, - MatrixXd &Q, - MatrixXd &O, - VectorXi &V2E, VectorXi &E2E, - VectorXi &boundary, VectorXi &nonmanifold, - std::vector& edge_diff, - std::vector& edge_values, - std::vector& face_edgeOrients, - std::vector& face_edgeIds, - std::map& singularities - ) { - struct EdgeLink - { +void subdivide_diff(MatrixXi &F, MatrixXd &V, MatrixXd &N, MatrixXd &Q, MatrixXd &O, VectorXi &V2E, + VectorXi &E2E, VectorXi &boundary, VectorXi &nonmanifold, + std::vector &edge_diff, std::vector &edge_values, + std::vector &face_edgeOrients, std::vector &face_edgeIds, + std::map &singularities) { + struct EdgeLink { int id; double length; Vector2i diff; }; - - struct FaceOrient - { + + struct FaceOrient { int orient; Vector3i d; Vector3d q; Vector3d n; }; - + std::vector face_spaces(F.cols()); std::queue queue; std::vector diffs(E2E.size()); @@ -193,10 +176,10 @@ void subdivide_diff(MatrixXi &F, for (int j = 0; j < 3; ++j) { int final_orient = face_edgeOrients[i][j]; int eid = face_edgeIds[i][j]; - auto value = compat_orientation_extrinsic_index_4(Q.col(edge_values[eid].x), N.col(edge_values[eid].x), orient.q, orient.n); + auto value = compat_orientation_extrinsic_index_4( + Q.col(edge_values[eid].x), N.col(edge_values[eid].x), orient.q, orient.n); int target_orient = (value.second - value.first + 4) % 4; - if (F(j, i) == edge_values[eid].y) - target_orient = (target_orient + 2) % 4; + if (F(j, i) == edge_values[eid].y) target_orient = (target_orient + 2) % 4; orient_diff[j] = (final_orient - target_orient + 4) % 4; } if (orient_diff[0] == orient_diff[1]) @@ -205,15 +188,14 @@ void subdivide_diff(MatrixXi &F, orient.orient = orient_diff[2]; else if (orient_diff[1] == orient_diff[2]) orient.orient = orient_diff[1]; - orient.d = Vector3i((orient_diff[0] - orient.orient + 4)%4, - (orient_diff[1] - orient.orient + 4)%4, - (orient_diff[2] - orient.orient + 4)%4); + orient.d = Vector3i((orient_diff[0] - orient.orient + 4) % 4, + (orient_diff[1] - orient.orient + 4) % 4, + (orient_diff[2] - orient.orient + 4) % 4); face_spaces[i] = (orient); } for (int i = 0; i < E2E.size(); ++i) { int v0 = F(i % 3, i / 3), v1 = F((i + 1) % 3, i / 3); - if (nonmanifold[v0] || nonmanifold[v1]) - continue; + if (nonmanifold[v0] || nonmanifold[v1]) continue; double length = (V.col(v0) - V.col(v1)).squaredNorm(); Vector2i diff = diffs[i]; if (abs(diff[0]) > 1 || abs(diff[1]) > 1) { @@ -227,13 +209,13 @@ void subdivide_diff(MatrixXi &F, } } } - auto AnalyzeOrient = [&](int f0, const Vector3i& d) { + auto AnalyzeOrient = [&](int f0, const Vector3i &d) { for (int j = 0; j < 3; ++j) { int orient = face_spaces[f0].orient + d[j]; int v = std::min(F(j, f0), F((j + 1) % 3, f0)); - auto value = compat_orientation_extrinsic_index_4(Q.col(v), N.col(v), face_spaces[f0].q, face_spaces[f0].n); - if (F(j, f0) != v) - orient += 2; + auto value = compat_orientation_extrinsic_index_4( + Q.col(v), N.col(v), face_spaces[f0].q, face_spaces[f0].n); + if (F(j, f0) != v) orient += 2; face_edgeOrients[f0][j] = (orient + value.second - value.first + 4) % 4; } face_spaces[f0].d = d; @@ -249,13 +231,14 @@ void subdivide_diff(MatrixXi &F, auto diff = edge_diff[face_edgeIds[f0][j]]; if (rshift90(diff, face_edgeOrients[f0][j]) != diffs[f0 * 3 + j]) { int orient = 0; - while (orient < 4 && rshift90(diff, orient) != diffs[f0 * 3 + j]) - orient += 1; - face_spaces[f0].d[j] = (face_spaces[f0].d[j] + orient - face_edgeOrients[f0][j]) % 4; + while (orient < 4 && rshift90(diff, orient) != diffs[f0 * 3 + j]) orient += 1; + face_spaces[f0].d[j] = + (face_spaces[f0].d[j] + orient - face_edgeOrients[f0][j]) % 4; face_edgeOrients[f0][j] = orient; } } }; + /* auto Length = [&](int f0) { int l = 0; for (int j = 0; j < 3; ++j) { @@ -267,18 +250,19 @@ void subdivide_diff(MatrixXi &F, printf("\n"); return l; }; + */ int nV = V.cols(), nF = F.cols(), nSplit = 0; /* / v0 \ v1p 1 | 0 v0p \ v1 / - + / v0 \ / 1 | 0 \ v1p - vn - v0p \ 2 | 3 / \ v1 / - + f0: vn, v0p, v0 f1: vn, v0, v1p f2: vn, v1p, v1 @@ -296,8 +280,7 @@ void subdivide_diff(MatrixXi &F, if ((V.col(v0) - V.col(v1)).squaredNorm() != edge.length) { continue; } - if (abs(diffs[e0][0]) < 2 && abs(diffs[e0][1]) < 2) - continue; + if (abs(diffs[e0][0]) < 2 && abs(diffs[e0][1]) < 2) continue; if (f1 != -1) { face_edgeOrients.push_back(Vector3i()); face_edgeIds.push_back(Vector3i()); @@ -314,36 +297,33 @@ void subdivide_diff(MatrixXi &F, boundary.conservativeResize(V.cols()); nonmanifold.conservativeResize(V.cols()); } - - + V.col(vn) = (V.col(v0) + V.col(v1)) * 0.5f; N.col(vn) = N.col(v0); Q.col(vn) = Q.col(v0); O.col(vn) = (O.col(v0) + O.col(v1)) * 0.5; nonmanifold[vn] = false; boundary[vn] = is_boundary; - + int eid = face_edgeIds[f0][e0 % 3]; int eid01 = face_edgeIds[f0][(e0 + 1) % 3]; int eid02 = face_edgeIds[f0][(e0 + 2) % 3]; int eid11 = face_edgeIds[f1][(e1 + 1) % 3]; int eid12 = face_edgeIds[f1][(e1 + 2) % 3]; - + int eid0, eid1, eid0p, eid1p; - eid0 = eid; edge_values[eid0] = DEdge(v0, vn); - + eid1 = edge_values.size(); edge_values.push_back(DEdge(vn, v1)); edge_diff.push_back(Vector2i()); - + eid0p = edge_values.size(); edge_values.push_back(DEdge(vn, v0p)); edge_diff.push_back(Vector2i()); - int f2 = is_boundary ? -1 : (nF++); int f3 = nF++; face_edgeIds.push_back(Vector3i()); @@ -357,24 +337,23 @@ void subdivide_diff(MatrixXi &F, } auto D01 = diffs[e0]; - auto D1p = diffs[e0/3*3+(e0+1)%3]; - auto Dp0 = diffs[e0/3*3+(e0+2)%3]; + auto D1p = diffs[e0 / 3 * 3 + (e0 + 1) % 3]; + auto Dp0 = diffs[e0 / 3 * 3 + (e0 + 2) % 3]; auto D0n = D01 / 2; auto Ds10 = diffs[e1]; - auto Ds0p = diffs[e1/3*3+(e1+1)%3]; - auto Dsp1 = diffs[e1/3*3+(e1+2)%3]; + auto Ds0p = diffs[e1 / 3 * 3 + (e1 + 1) % 3]; + auto Dsp1 = diffs[e1 / 3 * 3 + (e1 + 2) % 3]; int orient = 0; - while (rshift90(D01, orient) != Ds10) - orient += 1; + while (rshift90(D01, orient) != Ds10) orient += 1; auto Dsn0 = rshift90(D0n, orient); - + auto orients1 = face_spaces[f0]; auto orients2 = face_spaces[f1]; F.col(f0) << vn, v0p, v0; face_edgeIds[f0] = Vector3i(eid0p, eid02, eid0); - diffs[f0*3] = D01 + D1p - D0n; - diffs[f0*3+1] = Dp0; - diffs[f0*3+2] = D0n; + diffs[f0 * 3] = D01 + D1p - D0n; + diffs[f0 * 3 + 1] = Dp0; + diffs[f0 * 3 + 2] = D0n; int o1 = e0 % 3, o2 = e1 % 3; AnalyzeOrient(f0, Vector3i(0, orients1.d[(o1 + 2) % 3], orients1.d[o1])); if (!is_boundary) { @@ -383,46 +362,45 @@ void subdivide_diff(MatrixXi &F, edge_values.push_back(DEdge(vn, v1p)); edge_diff.push_back(Vector2i()); face_edgeIds[f1] = (Vector3i(eid0, eid11, eid1p)); - diffs[f1*3] = Dsn0; - diffs[f1*3+1] = Ds0p; - diffs[f1*3+2] = Dsp1 + (Ds10 - Dsn0); + diffs[f1 * 3] = Dsn0; + diffs[f1 * 3 + 1] = Ds0p; + diffs[f1 * 3 + 2] = Dsp1 + (Ds10 - Dsn0); AnalyzeOrient(f1, Vector3i(orients2.d[o2], orients2.d[(o2 + 1) % 3], 0)); face_spaces[f2] = face_spaces[f1]; face_edgeIds[f2] = (Vector3i(eid1p, eid12, eid1)); F.col(f2) << vn, v1p, v1; - diffs[f2*3] = -Dsp1 - (Ds10 - Dsn0); - diffs[f2*3+1] = Dsp1; - diffs[f2*3+2] = Ds10 - Dsn0; + diffs[f2 * 3] = -Dsp1 - (Ds10 - Dsn0); + diffs[f2 * 3 + 1] = Dsp1; + diffs[f2 * 3 + 2] = Ds10 - Dsn0; AnalyzeOrient(f2, Vector3i(0, orients2.d[(o2 + 2) % 3], orients2.d[o2])); } face_spaces[f3] = face_spaces[f0]; face_edgeIds[f3] = (Vector3i(eid1, eid01, eid0p)); F.col(f3) << vn, v1, v0p; - diffs[f3*3] = D01 - D0n; - diffs[f3*3+1] = D1p; - diffs[f3*3+2] = D0n - (D01 + D1p); + diffs[f3 * 3] = D01 - D0n; + diffs[f3 * 3 + 1] = D1p; + diffs[f3 * 3 + 2] = D0n - (D01 + D1p); AnalyzeOrient(f3, Vector3i(orients1.d[o1], orients1.d[(o1 + 1) % 3], 0)); - + FixOrient(f0); FixOrient(f1); FixOrient(f2); FixOrient(f3); - - const int e0p = E2E[dedge_prev_3(e0)], - e0n = E2E[dedge_next_3(e0)]; -#define sE2E(a, b) E2E[a] = b; if (b != -1) E2E[b] = a; + const int e0p = E2E[dedge_prev_3(e0)], e0n = E2E[dedge_next_3(e0)]; + +#define sE2E(a, b) \ + E2E[a] = b; \ + if (b != -1) E2E[b] = a; sE2E(3 * f0 + 0, 3 * f3 + 2); sE2E(3 * f0 + 1, e0p); sE2E(3 * f3 + 1, e0n); if (is_boundary) { sE2E(3 * f0 + 2, -1); sE2E(3 * f3 + 0, -1); - } - else { - const int e1p = E2E[dedge_prev_3(e1)], - e1n = E2E[dedge_next_3(e1)]; + } else { + const int e1p = E2E[dedge_prev_3(e1)], e1n = E2E[dedge_next_3(e1)]; sE2E(3 * f0 + 2, 3 * f1 + 0); sE2E(3 * f1 + 1, e1n); sE2E(3 * f1 + 2, 3 * f2 + 0); @@ -430,27 +408,25 @@ void subdivide_diff(MatrixXi &F, sE2E(3 * f2 + 2, 3 * f3 + 0); } #undef sE2E - - + V2E[v0] = 3 * f0 + 2; V2E[vn] = 3 * f0 + 0; V2E[v1] = 3 * f3 + 1; V2E[v0p] = 3 * f0 + 1; - if (!is_boundary) - V2E[v1p] = 3 * f1 + 2; - + if (!is_boundary) V2E[v1p] = 3 * f1 + 2; + auto schedule = [&](int f) { - for (int i = 0; i<3; ++i) { - if (abs(diffs[f*3+i][0]) > 1 || abs(diffs[f*3+i][1]) > 1) { + for (int i = 0; i < 3; ++i) { + if (abs(diffs[f * 3 + i][0]) > 1 || abs(diffs[f * 3 + i][1]) > 1) { EdgeLink e; e.id = f * 3 + i; - e.length = (V.col(F((i+1)%3, f)) - V.col(F(i, f))).squaredNorm(); - e.diff = diffs[f*3+i]; + e.length = (V.col(F((i + 1) % 3, f)) - V.col(F(i, f))).squaredNorm(); + e.diff = diffs[f * 3 + i]; queue.push(e); } } }; - + schedule(f0); if (!is_boundary) { schedule(f2); @@ -458,7 +434,7 @@ void subdivide_diff(MatrixXi &F, }; schedule(f3); } - + F.conservativeResize(F.rows(), nF); V.conservativeResize(V.rows(), nV); N.conservativeResize(V.rows(), nV); @@ -468,7 +444,7 @@ void subdivide_diff(MatrixXi &F, boundary.conservativeResize(nV); nonmanifold.conservativeResize(nV); E2E.conservativeResize(nF * 3); - + for (int i = 0; i < F.cols(); ++i) { for (int j = 0; j < 3; ++j) { auto diff = edge_diff[face_edgeIds[i][j]];