diff options
Diffstat (limited to 'extern/quadriflow/src/parametrizer-sing.cpp')
| -rw-r--r-- | extern/quadriflow/src/parametrizer-sing.cpp | 142 |
1 files changed, 142 insertions, 0 deletions
diff --git a/extern/quadriflow/src/parametrizer-sing.cpp b/extern/quadriflow/src/parametrizer-sing.cpp new file mode 100644 index 00000000000..2d600e3cda4 --- /dev/null +++ b/extern/quadriflow/src/parametrizer-sing.cpp @@ -0,0 +1,142 @@ +#include "config.hpp" +#include "field-math.hpp" +#include "parametrizer.hpp" + +namespace qflow { + +void Parametrizer::ComputeOrientationSingularities() { + MatrixXd &N = hierarchy.mN[0], &Q = hierarchy.mQ[0]; + const MatrixXi& F = hierarchy.mF; + singularities.clear(); + for (int f = 0; f < F.cols(); ++f) { + int index = 0; + int abs_index = 0; + for (int k = 0; k < 3; ++k) { + int i = F(k, f), j = F(k == 2 ? 0 : (k + 1), f); + auto value = + compat_orientation_extrinsic_index_4(Q.col(i), N.col(i), Q.col(j), N.col(j)); + index += value.second - value.first; + abs_index += std::abs(value.second - value.first); + } + int index_mod = modulo(index, 4); + if (index_mod == 1 || index_mod == 3) { + if (index >= 4 || index < 0) { + Q.col(F(0, f)) = -Q.col(F(0, f)); + } + singularities[f] = index_mod; + } + } +} + +void Parametrizer::ComputePositionSingularities() { + const MatrixXd &V = hierarchy.mV[0], &N = hierarchy.mN[0], &Q = hierarchy.mQ[0], + &O = hierarchy.mO[0]; + const MatrixXi& F = hierarchy.mF; + + pos_sing.clear(); + pos_rank.resize(F.rows(), F.cols()); + pos_index.resize(6, F.cols()); + for (int f = 0; f < F.cols(); ++f) { + Vector2i index = Vector2i::Zero(); + uint32_t i0 = F(0, f), i1 = F(1, f), i2 = F(2, f); + + Vector3d q[3] = {Q.col(i0).normalized(), Q.col(i1).normalized(), Q.col(i2).normalized()}; + Vector3d n[3] = {N.col(i0), N.col(i1), N.col(i2)}; + Vector3d o[3] = {O.col(i0), O.col(i1), O.col(i2)}; + Vector3d v[3] = {V.col(i0), V.col(i1), V.col(i2)}; + + int best[3]; + double best_dp = -std::numeric_limits<double>::infinity(); + for (int i = 0; i < 4; ++i) { + Vector3d v0 = rotate90_by(q[0], n[0], i); + for (int j = 0; j < 4; ++j) { + Vector3d v1 = rotate90_by(q[1], n[1], j); + for (int k = 0; k < 4; ++k) { + Vector3d v2 = rotate90_by(q[2], n[2], k); + double dp = std::min(std::min(v0.dot(v1), v1.dot(v2)), v2.dot(v0)); + if (dp > best_dp) { + best_dp = dp; + best[0] = i; + best[1] = j; + best[2] = k; + } + } + } + } + pos_rank(0, f) = best[0]; + pos_rank(1, f) = best[1]; + pos_rank(2, f) = best[2]; + for (int k = 0; k < 3; ++k) q[k] = rotate90_by(q[k], n[k], best[k]); + + for (int k = 0; k < 3; ++k) { + int kn = k == 2 ? 0 : (k + 1); + double scale_x = hierarchy.mScale, scale_y = hierarchy.mScale, + scale_x_1 = hierarchy.mScale, scale_y_1 = hierarchy.mScale; + if (flag_adaptive_scale) { + scale_x *= hierarchy.mS[0](0, F(k, f)); + scale_y *= hierarchy.mS[0](1, F(k, f)); + scale_x_1 *= hierarchy.mS[0](0, F(kn, f)); + scale_y_1 *= hierarchy.mS[0](1, F(kn, f)); + if (best[k] % 2 != 0) std::swap(scale_x, scale_y); + if (best[kn] % 2 != 0) std::swap(scale_x_1, scale_y_1); + } + 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<Vector2i, Vector2i> 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]; + pos_index(k * 2 + 1, f) = diff[1]; + } + + if (index != Vector2i::Zero()) { + pos_sing[f] = rshift90(index, best[0]); + } + } +} + +void Parametrizer::AnalyzeValence() { + auto& F = hierarchy.mF; + std::map<int, int> sing; + for (auto& f : singularities) { + for (int i = 0; i < 3; ++i) { + sing[F(i, f.first)] = f.second; + } + } + auto& F2E = face_edgeIds; + auto& E2E = hierarchy.mE2E; + auto& FQ = face_edgeOrients; + std::set<int> sing1, sing2; + for (int i = 0; i < F2E.size(); ++i) { + for (int j = 0; j < 3; ++j) { + int deid = i * 3 + j; + int sum_int = 0; + std::vector<int> edges; + std::vector<double> angles; + do { + int deid1 = deid / 3 * 3 + (deid + 2) % 3; + deid = E2E[deid1]; + sum_int += (FQ[deid / 3][deid % 3] + 6 - FQ[deid1 / 3][deid1 % 3]) % 4; + } while (deid != i * 3 + j); + if (sum_int % 4 == 2) { + printf("OMG! valence = 2\n"); + exit(0); + } + if (sum_int % 4 == 1) sing1.insert(F(j, i)); + if (sum_int % 4 == 3) sing2.insert(F(j, i)); + } + } + int count3 = 0, count4 = 0; + for (auto& s : singularities) { + if (s.second == 1) + count3 += 1; + else + count4 += 1; + } + printf("singularity: <%d %d> <%d %d>\n", (int)sing1.size(), (int)sing2.size(), count3, count4); +} + + +} // namespace qflow |