diff options
Diffstat (limited to 'extern/quadriflow/src/parametrizer.cpp')
| -rw-r--r-- | extern/quadriflow/src/parametrizer.cpp | 247 |
1 files changed, 247 insertions, 0 deletions
diff --git a/extern/quadriflow/src/parametrizer.cpp b/extern/quadriflow/src/parametrizer.cpp new file mode 100644 index 00000000000..b85383566c9 --- /dev/null +++ b/extern/quadriflow/src/parametrizer.cpp @@ -0,0 +1,247 @@ +#include "parametrizer.hpp" +#include "config.hpp" +#include "dedge.hpp" +#include "field-math.hpp" +#include "flow.hpp" +#include "localsat.hpp" +#include "optimizer.hpp" +#include "subdivide.hpp" + +#include "dset.hpp" + +#include <Eigen/Sparse> +#include <fstream> +#include <list> +#include <map> +#include <queue> +#include <set> + +namespace qflow { + +void Parametrizer::ComputeIndexMap(int with_scale) { + // build edge info + auto& V = hierarchy.mV[0]; + auto& F = hierarchy.mF; + auto& Q = hierarchy.mQ[0]; + auto& N = hierarchy.mN[0]; + auto& O = hierarchy.mO[0]; + auto& S = hierarchy.mS[0]; + // ComputeOrientationSingularities(); + + BuildEdgeInfo(); + + if (flag_preserve_sharp) { + // ComputeSharpO(); + } + for (int i = 0; i < sharp_edges.size(); ++i) { + if (sharp_edges[i]) { + int e = face_edgeIds[i / 3][i % 3]; + if (edge_diff[e][0] * edge_diff[e][1] != 0) { + Vector3d d = O.col(edge_values[e].y) - O.col(edge_values[e].x); + Vector3d q = Q.col(edge_values[e].x); + Vector3d n = N.col(edge_values[e].x); + Vector3d qy = n.cross(q); + if (abs(q.dot(d)) > qy.dot(d)) + edge_diff[e][1] = 0; + else + edge_diff[e][0] = 0; + } + } + } + std::map<int, std::pair<Vector3d, Vector3d>> sharp_constraints; + std::set<int> sharpvert; + for (int i = 0; i < sharp_edges.size(); ++i) { + if (sharp_edges[i]) { + sharpvert.insert(F(i % 3, i / 3)); + sharpvert.insert(F((i + 1) % 3, i / 3)); + } + } + + allow_changes.resize(edge_diff.size() * 2, 1); + for (int i = 0; i < sharp_edges.size(); ++i) { + int e = face_edgeIds[i / 3][i % 3]; + if (sharpvert.count(edge_values[e].x) && sharpvert.count(edge_values[e].y)) { + if (sharp_edges[i] != 0) { + for (int k = 0; k < 2; ++k) { + if (edge_diff[e][k] == 0) { + allow_changes[e * 2 + k] = 0; + } + } + } + } + } +#ifdef LOG_OUTPUT + printf("Build Integer Constraints...\n"); +#endif + BuildIntegerConstraints(); + + ComputeMaxFlow(); + // potential bug +#ifdef LOG_OUTPUT + printf("subdivide...\n"); +#endif + subdivide_edgeDiff(F, V, N, Q, O, &hierarchy.mS[0], V2E, hierarchy.mE2E, boundary, nonManifold, + edge_diff, edge_values, face_edgeOrients, face_edgeIds, sharp_edges, + singularities, 1); + + allow_changes.clear(); + allow_changes.resize(edge_diff.size() * 2, 1); + for (int i = 0; i < sharp_edges.size(); ++i) { + if (sharp_edges[i] == 0) continue; + int e = face_edgeIds[i / 3][i % 3]; + for (int k = 0; k < 2; ++k) { + if (edge_diff[e][k] == 0) allow_changes[e * 2 + k] = 0; + } + } + +#ifdef LOG_OUTPUT + printf("Fix flip advance...\n"); + int t1 = GetCurrentTime64(); +#endif + FixFlipHierarchy(); + subdivide_edgeDiff(F, V, N, Q, O, &hierarchy.mS[0], V2E, hierarchy.mE2E, boundary, nonManifold, + edge_diff, edge_values, face_edgeOrients, face_edgeIds, sharp_edges, + singularities, 1); + FixFlipSat(); + +#ifdef LOG_OUTPUT + int t2 = GetCurrentTime64(); + printf("Flip use %lf\n", (t2 - t1) * 1e-3); + printf("Post Linear Solver...\n"); +#endif + std::set<int> sharp_vertices; + for (int i = 0; i < sharp_edges.size(); ++i) { + if (sharp_edges[i] == 1) { + sharp_vertices.insert(F(i % 3, i / 3)); + sharp_vertices.insert(F((i + 1) % 3, i / 3)); + } + } + + Optimizer::optimize_positions_sharp(hierarchy, edge_values, edge_diff, sharp_edges, + sharp_vertices, sharp_constraints, with_scale); + + Optimizer::optimize_positions_fixed(hierarchy, edge_values, edge_diff, sharp_vertices, + sharp_constraints, flag_adaptive_scale); + + AdvancedExtractQuad(); + + FixValence(); + + std::vector<int> sharp_o(O_compact.size(), 0); + std::map<int, std::pair<Vector3d, Vector3d>> compact_sharp_constraints; + for (int i = 0; i < Vset.size(); ++i) { + int sharpv = -1; + for (auto& p : Vset[i]) { + if (sharp_constraints.count(p)) { + sharpv = p; + sharp_o[i] = 1; + if (compact_sharp_constraints.count(i) == 0 || + compact_sharp_constraints[i].second != Vector3d::Zero()) { + compact_sharp_constraints[i] = sharp_constraints[sharpv]; + O_compact[i] = O.col(sharpv); + compact_sharp_constraints[i].first = O_compact[i]; + } + } + } + } + + std::map<std::pair<int, int>, int> o2e; + for (int i = 0; i < F_compact.size(); ++i) { + for (int j = 0; j < 4; ++j) { + int v1 = F_compact[i][j]; + int v2 = F_compact[i][(j + 1) % 4]; + o2e[std::make_pair(v1, v2)] = i * 4 + j; + } + } + std::vector<std::vector<int>> v2o(V.cols()); + for (int i = 0; i < Vset.size(); ++i) { + for (auto v : Vset[i]) { + v2o[v].push_back(i); + } + } + std::vector<Vector3d> diffs(F_compact.size() * 4, Vector3d(0, 0, 0)); + std::vector<int> diff_count(F_compact.size() * 4, 0); + for (int i = 0; i < F.cols(); ++i) { + for (int j = 0; j < 3; ++j) { + int v1 = F(j, i); + int v2 = F((j + 1) % 3, i); + if (v1 != edge_values[face_edgeIds[i][j]].x) continue; + if (edge_diff[face_edgeIds[i][j]].array().abs().sum() != 1) continue; + if (v2o[v1].size() > 1 || v2o[v2].size() > 1) continue; + for (auto o1 : v2o[v1]) { + for (auto o2 : v2o[v2]) { + auto key = std::make_pair(o1, o2); + if (o2e.count(key)) { + int dedge = o2e[key]; + Vector3d q_1 = Q.col(v1); + Vector3d q_2 = Q.col(v2); + Vector3d n_1 = N.col(v1); + Vector3d n_2 = N.col(v2); + Vector3d q_1_y = n_1.cross(q_1); + Vector3d q_2_y = n_2.cross(q_2); + auto index = compat_orientation_extrinsic_index_4(q_1, n_1, q_2, n_2); + double s_x1 = S(0, v1), s_y1 = S(1, v1); + double s_x2 = S(0, v2), s_y2 = S(1, v2); + int rank_diff = (index.second + 4 - index.first) % 4; + if (rank_diff % 2 == 1) std::swap(s_x2, s_y2); + Vector3d qd_x = 0.5 * (rotate90_by(q_2, n_2, rank_diff) + q_1); + Vector3d qd_y = 0.5 * (rotate90_by(q_2_y, n_2, rank_diff) + q_1_y); + double scale_x = (with_scale ? 0.5 * (s_x1 + s_x2) : 1) * hierarchy.mScale; + double scale_y = (with_scale ? 0.5 * (s_y1 + s_y2) : 1) * hierarchy.mScale; + Vector2i diff = edge_diff[face_edgeIds[i][j]]; + Vector3d C = diff[0] * scale_x * qd_x + diff[1] * scale_y * qd_y; + + diff_count[dedge] += 1; + diffs[dedge] += C; + auto key = std::make_pair(o2, o1); + if (o2e.count(key)) { + int dedge = o2e[key]; + diff_count[dedge] += 1; + diffs[dedge] -= C; + } + } + } + } + } + } + + for (int i = 0; i < F.cols(); ++i) { + Vector2i d1 = rshift90(edge_diff[face_edgeIds[i][0]], face_edgeOrients[i][0]); + Vector2i d2 = rshift90(edge_diff[face_edgeIds[i][1]], face_edgeOrients[i][1]); + if (d1[0] * d2[1] - d1[1] * d2[0] < 0) { + for (int j = 0; j < 3; ++j) { + int v1 = F(j, i); + int v2 = F((j + 1) % 3, i); + for (auto o1 : v2o[v1]) { + for (auto o2 : v2o[v2]) { + auto key = std::make_pair(o1, o2); + if (o2e.count(key)) { + int dedge = o2e[key]; + diff_count[dedge] = 0; + diffs[dedge] = Vector3d(0, 0, 0); + } + } + } + } + } + } + + for (int i = 0; i < diff_count.size(); ++i) { + if (diff_count[i] != 0) { + diffs[i] /= diff_count[i]; + diff_count[i] = 1; + } + } + + Optimizer::optimize_positions_dynamic(F, V, N, Q, Vset, O_compact, F_compact, V2E_compact, + E2E_compact, sqrt(surface_area / F_compact.size()), + diffs, diff_count, o2e, sharp_o, + compact_sharp_constraints, flag_adaptive_scale); + + // 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, false); +} + +} // namespace qflow |