Welcome to mirror list, hosted at ThFree Co, Russian Federation.

git.blender.org/blender.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
path: root/source/blender/blenlib/intern/mesh_intersect.cc
diff options
context:
space:
mode:
Diffstat (limited to 'source/blender/blenlib/intern/mesh_intersect.cc')
-rw-r--r--source/blender/blenlib/intern/mesh_intersect.cc458
1 files changed, 9 insertions, 449 deletions
diff --git a/source/blender/blenlib/intern/mesh_intersect.cc b/source/blender/blenlib/intern/mesh_intersect.cc
index c36dfa80be7..1eea58e4eb2 100644
--- a/source/blender/blenlib/intern/mesh_intersect.cc
+++ b/source/blender/blenlib/intern/mesh_intersect.cc
@@ -24,6 +24,7 @@
# include <algorithm>
# include <fstream>
# include <iostream>
+# include <memory>
# include "BLI_allocator.hh"
# include "BLI_array.hh"
@@ -519,7 +520,7 @@ class IMeshArena::IMeshArenaImpl : NonCopyable, NonMovable {
IMeshArena::IMeshArena()
{
- pimpl_ = std::unique_ptr<IMeshArenaImpl>(new IMeshArenaImpl());
+ pimpl_ = std::make_unique<IMeshArenaImpl>();
}
IMeshArena::~IMeshArena()
@@ -1107,254 +1108,6 @@ static mpq2 project_3d_to_2d(const mpq3 &p3d, int proj_axis)
}
/**
- Is a point in the interior of a 2d triangle or on one of its
- * edges but not either endpoint of the edge?
- * orient[pi][i] is the orientation test of the point pi against
- * the side of the triangle starting at index i.
- * Assume the triangle is non-degenerate and CCW-oriented.
- * Then answer is true if p is left of or on all three of triangle a's edges,
- * and strictly left of at least on of them.
- */
-static bool non_trivially_2d_point_in_tri(const int orients[3][3], int pi)
-{
- int p_left_01 = orients[pi][0];
- int p_left_12 = orients[pi][1];
- int p_left_20 = orients[pi][2];
- return (p_left_01 >= 0 && p_left_12 >= 0 && p_left_20 >= 0 &&
- (p_left_01 + p_left_12 + p_left_20) >= 2);
-}
-
-/**
- * Given orients as defined in non_trivially_2d_intersect, do the triangles
- * overlap in a "hex" pattern? That is, the overlap region is a hexagon, which
- * one gets by having, each point of one triangle being strictly right-of one
- * edge of the other and strictly left of the other two edges; and vice versa.
- * In addition, it must not be the case that all of the points of one triangle
- * are totally on the outside of one edge of the other triangle, and vice versa.
- */
-static bool non_trivially_2d_hex_overlap(int orients[2][3][3])
-{
- for (int ab = 0; ab < 2; ++ab) {
- for (int i = 0; i < 3; ++i) {
- bool ok = orients[ab][i][0] + orients[ab][i][1] + orients[ab][i][2] == 1 &&
- orients[ab][i][0] != 0 && orients[ab][i][1] != 0 && orients[i][2] != 0;
- if (!ok) {
- return false;
- }
- int s = orients[ab][0][i] + orients[ab][1][i] + orients[ab][2][i];
- if (s == -3) {
- return false;
- }
- }
- }
- return true;
-}
-
-/**
- * Given orients as defined in non_trivially_2d_intersect, do the triangles
- * have one shared edge in a "folded-over" configuration?
- * As well as a shared edge, the third vertex of one triangle needs to be
- * right-of one and left-of the other two edges of the other triangle.
- */
-static bool non_trivially_2d_shared_edge_overlap(int orients[2][3][3],
- const mpq2 *a[3],
- const mpq2 *b[3])
-{
- for (int i = 0; i < 3; ++i) {
- int in = (i + 1) % 3;
- int inn = (i + 2) % 3;
- for (int j = 0; j < 3; ++j) {
- int jn = (j + 1) % 3;
- int jnn = (j + 2) % 3;
- if (*a[i] == *b[j] && *a[in] == *b[jn]) {
- /* Edge from a[i] is shared with edge from b[j]. */
- /* See if a[inn] is right-of or on one of the other edges of b.
- * If it is on, then it has to be right-of or left-of the shared edge,
- * depending on which edge it is. */
- if (orients[0][inn][jn] < 0 || orients[0][inn][jnn] < 0) {
- return true;
- }
- if (orients[0][inn][jn] == 0 && orients[0][inn][j] == 1) {
- return true;
- }
- if (orients[0][inn][jnn] == 0 && orients[0][inn][j] == -1) {
- return true;
- }
- /* Similarly for `b[jnn]`. */
- if (orients[1][jnn][in] < 0 || orients[1][jnn][inn] < 0) {
- return true;
- }
- if (orients[1][jnn][in] == 0 && orients[1][jnn][i] == 1) {
- return true;
- }
- if (orients[1][jnn][inn] == 0 && orients[1][jnn][i] == -1) {
- return true;
- }
- }
- }
- }
- return false;
-}
-
-/**
- * Are the triangles the same, perhaps with some permutation of vertices?
- */
-static bool same_triangles(const mpq2 *a[3], const mpq2 *b[3])
-{
- for (int i = 0; i < 3; ++i) {
- if (a[0] == b[i] && a[1] == b[(i + 1) % 3] && a[2] == b[(i + 2) % 3]) {
- return true;
- }
- }
- return false;
-}
-
-/**
- * Do 2d triangles (a[0], a[1], a[2]) and (b[0], b[1], b2[2]) intersect at more than just shared
- * vertices or a shared edge? This is true if any point of one triangle is non-trivially inside the
- * other. NO: that isn't quite sufficient: there is also the case where the verts are all mutually
- * outside the other's triangle, but there is a hexagonal overlap region where they overlap.
- */
-static bool non_trivially_2d_intersect(const mpq2 *a[3], const mpq2 *b[3])
-{
- /* TODO: Could experiment with trying bounding box tests before these.
- * TODO: Find a less expensive way than 18 orient tests to do this. */
-
- /* `orients[0][ai][bi]` is orient of point `a[ai]` compared to segment starting at `b[bi]`.
- * `orients[1][bi][ai]` is orient of point `b[bi]` compared to segment starting at `a[ai]`. */
- int orients[2][3][3];
- for (int ab = 0; ab < 2; ++ab) {
- for (int ai = 0; ai < 3; ++ai) {
- for (int bi = 0; bi < 3; ++bi) {
- if (ab == 0) {
- orients[0][ai][bi] = orient2d(*b[bi], *b[(bi + 1) % 3], *a[ai]);
- }
- else {
- orients[1][bi][ai] = orient2d(*a[ai], *a[(ai + 1) % 3], *b[bi]);
- }
- }
- }
- }
- return non_trivially_2d_point_in_tri(orients[0], 0) ||
- non_trivially_2d_point_in_tri(orients[0], 1) ||
- non_trivially_2d_point_in_tri(orients[0], 2) ||
- non_trivially_2d_point_in_tri(orients[1], 0) ||
- non_trivially_2d_point_in_tri(orients[1], 1) ||
- non_trivially_2d_point_in_tri(orients[1], 2) || non_trivially_2d_hex_overlap(orients) ||
- non_trivially_2d_shared_edge_overlap(orients, a, b) || same_triangles(a, b);
- return true;
-}
-
-/**
- * Does triangle t in tm non-trivially non-co-planar intersect any triangle
- * in `CoplanarCluster cl`? Assume t is known to be in the same plane as all
- * the triangles in cl, and that proj_axis is a good axis to project down
- * to solve this problem in 2d.
- */
-static bool non_trivially_coplanar_intersects(const IMesh &tm,
- int t,
- const CoplanarCluster &cl,
- int proj_axis,
- const Map<std::pair<int, int>, ITT_value> &itt_map)
-{
- const Face &tri = *tm.face(t);
- mpq2 v0 = project_3d_to_2d(tri[0]->co_exact, proj_axis);
- mpq2 v1 = project_3d_to_2d(tri[1]->co_exact, proj_axis);
- mpq2 v2 = project_3d_to_2d(tri[2]->co_exact, proj_axis);
- if (orient2d(v0, v1, v2) != 1) {
- mpq2 tmp = v1;
- v1 = v2;
- v2 = tmp;
- }
- for (const int cl_t : cl) {
- if (!itt_map.contains(std::pair<int, int>(t, cl_t)) &&
- !itt_map.contains(std::pair<int, int>(cl_t, t))) {
- continue;
- }
- const Face &cl_tri = *tm.face(cl_t);
- mpq2 ctv0 = project_3d_to_2d(cl_tri[0]->co_exact, proj_axis);
- mpq2 ctv1 = project_3d_to_2d(cl_tri[1]->co_exact, proj_axis);
- mpq2 ctv2 = project_3d_to_2d(cl_tri[2]->co_exact, proj_axis);
- if (orient2d(ctv0, ctv1, ctv2) != 1) {
- mpq2 tmp = ctv1;
- ctv1 = ctv2;
- ctv2 = tmp;
- }
- const mpq2 *v[] = {&v0, &v1, &v2};
- const mpq2 *ctv[] = {&ctv0, &ctv1, &ctv2};
- if (non_trivially_2d_intersect(v, ctv)) {
- return true;
- }
- }
- return false;
-}
-
-/* Keeping this code for a while, but for now, almost all
- * trivial intersects are found before calling intersect_tri_tri now.
- */
-# if 0
-/**
- * Do tri1 and tri2 intersect at all, and if so, is the intersection
- * something other than a common vertex or a common edge?
- * The \a itt value is the result of calling intersect_tri_tri on tri1, tri2.
- */
-static bool non_trivial_intersect(const ITT_value &itt, const Face * tri1, const Face * tri2)
-{
- if (itt.kind == INONE) {
- return false;
- }
- const Face * tris[2] = {tri1, tri2};
- if (itt.kind == IPOINT) {
- bool has_p_as_vert[2] {false, false};
- for (int i = 0; i < 2; ++i) {
- for (const Vert * v : *tris[i]) {
- if (itt.p1 == v->co_exact) {
- has_p_as_vert[i] = true;
- break;
- }
- }
- }
- return !(has_p_as_vert[0] && has_p_as_vert[1]);
- }
- if (itt.kind == ISEGMENT) {
- bool has_seg_as_edge[2] = {false, false};
- for (int i = 0; i < 2; ++i) {
- const Face &t = *tris[i];
- for (int pos : t.index_range()) {
- int nextpos = t.next_pos(pos);
- if ((itt.p1 == t[pos]->co_exact && itt.p2 == t[nextpos]->co_exact) ||
- (itt.p2 == t[pos]->co_exact && itt.p1 == t[nextpos]->co_exact)) {
- has_seg_as_edge[i] = true;
- break;
- }
- }
- }
- return !(has_seg_as_edge[0] && has_seg_as_edge[1]);
- }
- BLI_assert(itt.kind == ICOPLANAR);
- /* TODO: refactor this common code with code above. */
- int proj_axis = mpq3::dominant_axis(tri1->plane.norm_exact);
- mpq2 tri_2d[2][3];
- for (int i = 0; i < 2; ++i) {
- mpq2 v0 = project_3d_to_2d((*tris[i])[0]->co_exact, proj_axis);
- mpq2 v1 = project_3d_to_2d((*tris[i])[1]->co_exact, proj_axis);
- mpq2 v2 = project_3d_to_2d((*tris[i])[2]->co_exact, proj_axis);
- if (mpq2::orient2d(v0, v1, v2) != 1) {
- mpq2 tmp = v1;
- v1 = v2;
- v2 = tmp;
- }
- tri_2d[i][0] = v0;
- tri_2d[i][1] = v1;
- tri_2d[i][2] = v2;
- }
- const mpq2 *va[] = {&tri_2d[0][0], &tri_2d[0][1], &tri_2d[0][2]};
- const mpq2 *vb[] = {&tri_2d[1][0], &tri_2d[1][1], &tri_2d[1][2]};
- return non_trivially_2d_intersect(va, vb);
-}
-# endif
-
-/**
* The sup and index functions are defined in the paper:
* EXACT GEOMETRIC COMPUTATION USING CASCADING, by
* Burnikel, Funke, and Seel. They are used to find absolute
@@ -1414,119 +1167,6 @@ constexpr int index_dot_plane_coords = 15;
*/
constexpr int index_dot_cross = 11;
-/* Not using this at the moment. Leaving it for a bit in case we want it again. */
-# if 0
-static double supremum_dot(const double3 &a, const double3 &b)
-{
- double3 abs_a = double3::abs(a);
- double3 abs_b = double3::abs(b);
- return double3::dot(abs_a, abs_b);
-}
-
-static double supremum_orient3d(const double3 &a,
- const double3 &b,
- const double3 &c,
- const double3 &d)
-{
- double3 abs_a = double3::abs(a);
- double3 abs_b = double3::abs(b);
- double3 abs_c = double3::abs(c);
- double3 abs_d = double3::abs(d);
- double adx = abs_a[0] + abs_d[0];
- double bdx = abs_b[0] + abs_d[0];
- double cdx = abs_c[0] + abs_d[0];
- double ady = abs_a[1] + abs_d[1];
- double bdy = abs_b[1] + abs_d[1];
- double cdy = abs_c[1] + abs_d[1];
- double adz = abs_a[2] + abs_d[2];
- double bdz = abs_b[2] + abs_d[2];
- double cdz = abs_c[2] + abs_d[2];
-
- double bdxcdy = bdx * cdy;
- double cdxbdy = cdx * bdy;
-
- double cdxady = cdx * ady;
- double adxcdy = adx * cdy;
-
- double adxbdy = adx * bdy;
- double bdxady = bdx * ady;
-
- double det = adz * (bdxcdy + cdxbdy) + bdz * (cdxady + adxcdy) + cdz * (adxbdy + bdxady);
- return det;
-}
-
-/** Actually index_orient3d = 10 + 4 * (max degree of input coordinates) */
-constexpr int index_orient3d = 14;
-
-/**
- * Return the approximate orient3d of the four double3's, with
- * the guarantee that if the value is -1 or 1 then the underlying
- * mpq3 test would also have returned that value.
- * When the return value is 0, we are not sure of the sign.
- */
-static int filter_orient3d(const double3 &a, const double3 &b, const double3 &c, const double3 &d)
-{
- double o3dfast = orient3d_fast(a, b, c, d);
- if (o3dfast == 0.0) {
- return 0;
- }
- double err_bound = supremum_orient3d(a, b, c, d) * index_orient3d * DBL_EPSILON;
- if (fabs(o3dfast) > err_bound) {
- return o3dfast > 0.0 ? 1 : -1;
- }
- return 0;
-}
-
-/**
- * Return the approximate orient3d of the triangle plane points and v, with
- * the guarantee that if the value is -1 or 1 then the underlying
- * mpq3 test would also have returned that value.
- * When the return value is 0, we are not sure of the sign.
- */
-static int filter_tri_plane_vert_orient3d(const Face &tri, const Vert *v)
-{
- return filter_orient3d(tri[0]->co, tri[1]->co, tri[2]->co, v->co);
-}
-
-/**
- * Are vectors a and b parallel or nearly parallel?
- * This routine should only return false if we are certain
- * that they are not parallel, taking into account the
- * possible numeric errors and input value approximation.
- */
-static bool near_parallel_vecs(const double3 &a, const double3 &b)
-{
- double3 cr = double3::cross_high_precision(a, b);
- double cr_len_sq = cr.length_squared();
- if (cr_len_sq == 0.0) {
- return true;
- }
- double err_bound = supremum_dot_cross(a, b) * index_dot_cross * DBL_EPSILON;
- if (cr_len_sq > err_bound) {
- return false;
- }
- return true;
-}
-
-/**
- * Return true if we are sure that dot(a,b) > 0, taking into
- * account the error bounds due to numeric errors and input value
- * approximation.
- */
-static bool dot_must_be_positive(const double3 &a, const double3 &b)
-{
- double d = double3::dot(a, b);
- if (d <= 0.0) {
- return false;
- }
- double err_bound = supremum_dot(a, b) * index_dot_plane_coords * DBL_EPSILON;
- if (d > err_bound) {
- return true;
- }
- return false;
-}
-# endif
-
/**
* Return the approximate side of point p on a plane with normal plane_no and point plane_p.
* The answer will be 1 if p is definitely above the plane, -1 if it is definitely below.
@@ -1550,89 +1190,12 @@ static int filter_plane_side(const double3 &p,
}
double supremum = double3::dot(abs_p + abs_plane_p, abs_plane_no);
double err_bound = supremum * index_plane_side * DBL_EPSILON;
- if (d > err_bound) {
+ if (fabs(d) > err_bound) {
return d > 0 ? 1 : -1;
}
return 0;
}
-/* Not using this at the moment. Leave it here for a while in case we want it again. */
-# if 0
-/**
- * A fast, non-exhaustive test for non_trivial intersection.
- * If this returns false then we are sure that tri1 and tri2
- * do not intersect. If it returns true, they may or may not
- * non-trivially intersect.
- * We assume that bounding box overlap tests have already been
- * done, so don't repeat those here. This routine is checking
- * for the very common cases (when doing mesh self-intersect)
- * where triangles share an edge or a vertex, but don't
- * otherwise intersect.
- */
-static bool may_non_trivially_intersect(Face *t1, Face *t2)
-{
- Face &tri1 = *t1;
- Face &tri2 = *t2;
- BLI_assert(t1->plane_populated() && t2->plane_populated());
- Face::FacePos share1_pos[3];
- Face::FacePos share2_pos[3];
- int n_shared = 0;
- for (Face::FacePos p1 = 0; p1 < 3; ++p1) {
- const Vert *v1 = tri1[p1];
- for (Face::FacePos p2 = 0; p2 < 3; ++p2) {
- const Vert *v2 = tri2[p2];
- if (v1 == v2) {
- share1_pos[n_shared] = p1;
- share2_pos[n_shared] = p2;
- ++n_shared;
- }
- }
- }
- if (n_shared == 2) {
- /* t1 and t2 share an entire edge.
- * If their normals are not parallel, they cannot non-trivially intersect. */
- if (!near_parallel_vecs(tri1.plane->norm, tri2.plane->norm)) {
- return false;
- }
- /* The normals are parallel or nearly parallel.
- * If the normals are in the same direction and the edges have opposite
- * directions in the two triangles, they cannot non-trivially intersect. */
- bool erev1 = tri1.prev_pos(share1_pos[0]) == share1_pos[1];
- bool erev2 = tri2.prev_pos(share2_pos[0]) == share2_pos[1];
- if (erev1 != erev2) {
- if (dot_must_be_positive(tri1.plane->norm, tri2.plane->norm)) {
- return false;
- }
- }
- }
- else if (n_shared == 1) {
- /* t1 and t2 share a vertex, but not an entire edge.
- * If the two non-shared verts of t2 are both on the same
- * side of tri1's plane, then they cannot non-trivially intersect.
- * (There are some other cases that could be caught here but
- * they are more expensive to check). */
- Face::FacePos p = share2_pos[0];
- const Vert *v2a = p == 0 ? tri2[1] : tri2[0];
- const Vert *v2b = (p == 0 || p == 1) ? tri2[2] : tri2[1];
- int o1 = filter_tri_plane_vert_orient3d(tri1, v2a);
- int o2 = filter_tri_plane_vert_orient3d(tri1, v2b);
- if (o1 == o2 && o1 != 0) {
- return false;
- }
- p = share1_pos[0];
- const Vert *v1a = p == 0 ? tri1[1] : tri1[0];
- const Vert *v1b = (p == 0 || p == 1) ? tri1[2] : tri1[1];
- o1 = filter_tri_plane_vert_orient3d(tri2, v1a);
- o2 = filter_tri_plane_vert_orient3d(tri2, v1b);
- if (o1 == o2 && o1 != 0) {
- return false;
- }
- }
- /* We weren't able to prove that any intersection is trivial. */
- return true;
-}
-# endif
-
/*
* interesect_tri_tri and helper functions.
* This code uses the algorithm of Guigue and Devillers, as described
@@ -2460,12 +2023,12 @@ class TriOverlaps {
if (two_trees_no_self) {
BLI_bvhtree_balance(tree_b_);
/* Don't expect a lot of trivial intersects in this case. */
- overlap_ = BLI_bvhtree_overlap(tree_, tree_b_, &overlap_tot_, NULL, NULL);
+ overlap_ = BLI_bvhtree_overlap(tree_, tree_b_, &overlap_tot_, nullptr, nullptr);
}
else {
CBData cbdata{tm, shape_fn, nshapes, use_self};
if (nshapes == 1) {
- overlap_ = BLI_bvhtree_overlap(tree_, tree_, &overlap_tot_, NULL, NULL);
+ overlap_ = BLI_bvhtree_overlap(tree_, tree_, &overlap_tot_, nullptr, nullptr);
}
else {
overlap_ = BLI_bvhtree_overlap(
@@ -2637,8 +2200,7 @@ struct SubdivideTrisData {
tm(tm),
itt_map(itt_map),
overlap(overlap),
- arena(arena),
- overlap_tri_range{}
+ arena(arena)
{
}
};
@@ -2771,7 +2333,7 @@ static CDT_data calc_cluster_subdivided(const CoplanarClusterInfo &clinfo,
std::pair<int, int> key = canon_int_pair(t, t_other);
if (itt_map.contains(key)) {
ITT_value itt = itt_map.lookup(key);
- if (itt.kind != INONE && itt.kind != ICOPLANAR) {
+ if (!ELEM(itt.kind, INONE, ICOPLANAR)) {
itts.append(itt);
if (dbg_level > 0) {
std::cout << " itt = " << itt << "\n";
@@ -2852,7 +2414,6 @@ static CoplanarClusterInfo find_clusters(const IMesh &tm,
if (dbg_level > 0) {
std::cout << "already has " << curcls.size() << " clusters\n";
}
- int proj_axis = mpq3::dominant_axis(tplane.norm_exact);
/* Partition `curcls` into those that intersect t non-trivially, and those that don't. */
Vector<CoplanarCluster *> int_cls;
Vector<CoplanarCluster *> no_int_cls;
@@ -2860,8 +2421,7 @@ static CoplanarClusterInfo find_clusters(const IMesh &tm,
if (dbg_level > 1) {
std::cout << "consider intersecting with cluster " << cl << "\n";
}
- if (bbs_might_intersect(tri_bb[t], cl.bounding_box()) &&
- non_trivially_coplanar_intersects(tm, t, cl, proj_axis, itt_map)) {
+ if (bbs_might_intersect(tri_bb[t], cl.bounding_box())) {
if (dbg_level > 1) {
std::cout << "append to int_cls\n";
}
@@ -3317,6 +2877,6 @@ static void dump_perfdata(void)
}
# endif
-}; // namespace blender::meshintersect
+} // namespace blender::meshintersect
#endif // WITH_GMP
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-18 21:35:20 +0300