[Bf-blender-cvs] [383b4c07274] newboolean: Added floating filters to the initial plane-side tests in tri-tri intersect.
Howard Trickey
noreply at git.blender.org
Sun Jul 19 20:18:15 CEST 2020
Commit: 383b4c07274cba538efac5571e7c519ff45f3e97
Author: Howard Trickey
Date: Sun Jul 19 14:16:39 2020 -0400
Branches: newboolean
https://developer.blender.org/rB383b4c07274cba538efac5571e7c519ff45f3e97
Added floating filters to the initial plane-side tests in tri-tri intersect.
Needed an "abs" function for double3, so added it to all of the
float/double/mpq 2/3 types.
===================================================================
M source/blender/blenlib/BLI_double2.hh
M source/blender/blenlib/BLI_double3.hh
M source/blender/blenlib/BLI_float2.hh
M source/blender/blenlib/BLI_float3.hh
M source/blender/blenlib/BLI_mpq2.hh
M source/blender/blenlib/BLI_mpq3.hh
M source/blender/blenlib/intern/mesh_intersect.cc
M tests/gtests/blenlib/BLI_mesh_intersect_test.cc
===================================================================
diff --git a/source/blender/blenlib/BLI_double2.hh b/source/blender/blenlib/BLI_double2.hh
index c505f860f3c..2685cdc29ab 100644
--- a/source/blender/blenlib/BLI_double2.hh
+++ b/source/blender/blenlib/BLI_double2.hh
@@ -100,6 +100,11 @@ struct double2 {
return a * (1 - t) + b * t;
}
+ static double2 abs(const double2 &a)
+ {
+ return double2(fabsf(a.x), fabsf(a.y));
+ }
+
static double distance(const double2 &a, const double2 &b)
{
return (a - b).length();
diff --git a/source/blender/blenlib/BLI_double3.hh b/source/blender/blenlib/BLI_double3.hh
index 72085df2eb9..5531066ca36 100644
--- a/source/blender/blenlib/BLI_double3.hh
+++ b/source/blender/blenlib/BLI_double3.hh
@@ -217,6 +217,11 @@ struct double3 {
return a * (1 - t) + b * t;
}
+ static double3 abs(const double3 &a)
+ {
+ return double3(fabs(a.x), fabs(a.y), fabs(a.z));
+ }
+
/* orient3d gives the exact result, using multiprecision artihmetic when result
* is close to zero. orient3d_fast just uses double arithmetic, so may be
* wrong if the answer is very close to zero.
diff --git a/source/blender/blenlib/BLI_float2.hh b/source/blender/blenlib/BLI_float2.hh
index b7d07aac23e..8d956955f57 100644
--- a/source/blender/blenlib/BLI_float2.hh
+++ b/source/blender/blenlib/BLI_float2.hh
@@ -123,6 +123,11 @@ struct float2 {
return a * (1 - t) + b * t;
}
+ static float2 abs(const float2 &a)
+ {
+ return float2(fabsf(a.x), fabsf(a.y));
+ }
+
static float distance(const float2 &a, const float2 &b)
{
return (a - b).length();
diff --git a/source/blender/blenlib/BLI_float3.hh b/source/blender/blenlib/BLI_float3.hh
index a36cedad41d..d490108c2b2 100644
--- a/source/blender/blenlib/BLI_float3.hh
+++ b/source/blender/blenlib/BLI_float3.hh
@@ -229,6 +229,11 @@ struct float3 {
{
return a * (1 - t) + b * t;
}
+
+ static float3 abs(const float3 &a)
+ {
+ return float3(fabsf(a.x), fabsf(a.y), fabsf(a.z));
+ }
};
} // namespace blender
diff --git a/source/blender/blenlib/BLI_mpq2.hh b/source/blender/blenlib/BLI_mpq2.hh
index fc4702e4c01..e585c1d5a0a 100644
--- a/source/blender/blenlib/BLI_mpq2.hh
+++ b/source/blender/blenlib/BLI_mpq2.hh
@@ -131,6 +131,13 @@ struct mpq2 {
return a * (1 - t) + b * t;
}
+ static mpq2 abs(const mpq2 &a)
+ {
+ mpq_class abs_x = (a.x >= 0) ? a.x : -a.x;
+ mpq_class abs_y = (a.y >= 0) ? a.y : -a.y;
+ return mpq2(abs_x, abs_y);
+ }
+
static mpq_class distance(const mpq2 &a, const mpq2 &b)
{
return (a - b).length();
diff --git a/source/blender/blenlib/BLI_mpq3.hh b/source/blender/blenlib/BLI_mpq3.hh
index e8f11595eed..41ebe24053d 100644
--- a/source/blender/blenlib/BLI_mpq3.hh
+++ b/source/blender/blenlib/BLI_mpq3.hh
@@ -241,11 +241,19 @@ struct mpq3 {
return mpq3(a.x * s + b.x * t, a.y * s + b.y * t, a.z * s + b.z * t);
}
+ static mpq3 abs(const mpq3 &a)
+ {
+ mpq_class abs_x = (a.x >= 0) ? a.x : -a.x;
+ mpq_class abs_y = (a.y >= 0) ? a.y : -a.y;
+ mpq_class abs_z = (a.z >= 0) ? a.z : -a.z;
+ return mpq3(abs_x, abs_y, abs_z);
+ }
+
static int dominant_axis(const mpq3 &a)
{
- mpq_class x = abs(a[0]);
- mpq_class y = abs(a[1]);
- mpq_class z = abs(a[2]);
+ mpq_class x = (a.x >= 0) ? a.x : -a.x;
+ mpq_class y = (a.y >= 0) ? a.y : -a.y;
+ mpq_class z = (a.z >= 0) ? a.z : -a.z;
return ((x > y) ? ((x > z) ? 0 : 2) : ((y > z) ? 1 : 2));
}
diff --git a/source/blender/blenlib/intern/mesh_intersect.cc b/source/blender/blenlib/intern/mesh_intersect.cc
index 47daa1cbe48..c4befd32f5c 100644
--- a/source/blender/blenlib/intern/mesh_intersect.cc
+++ b/source/blender/blenlib/intern/mesh_intersect.cc
@@ -1181,15 +1181,15 @@ static bool non_trivial_intersect(const ITT_value &itt, Facep tri1, Facep tri2)
static double supremum_cross(const double3 &a, const double3 &b)
{
- double3 abs_a{fabs(a[0]), fabs(a[1]), fabs(a[2])};
- double3 abs_b{fabs(b[0]), fabs(b[1]), fabs(b[2])};
+ double3 abs_a = double3::abs(a);
+ double3 abs_b = double3::abs(b);
double3 c;
/* This is cross(a, b) but using absoluate values for a and b
* and always using + when operation is + or -.
*/
- c[0] = a[1] * b[2] + a[2] * b[1];
- c[1] = a[2] * b[0] + a[0] * b[2];
- c[2] = a[0] * b[1] + a[1] * b[0];
+ c[0] = abs_a[1] * abs_b[2] + abs_a[2] * abs_b[1];
+ c[1] = abs_a[2] * abs_b[0] + abs_a[0] * abs_b[2];
+ c[2] = abs_a[0] * abs_b[1] + abs_a[1] * abs_b[0];
return double3::dot(c, c);
}
@@ -1201,23 +1201,23 @@ constexpr int index_cross = 11;
static double supremum_dot(const double3 &a, const double3 &b)
{
- double3 abs_a{fabs(a[0]), fabs(a[1]), fabs(a[2])};
- double3 abs_b{fabs(b[0]), fabs(b[1]), fabs(b[2])};
+ double3 abs_a = double3::abs(a);
+ double3 abs_b = double3::abs(b);
return double3::dot(abs_a, abs_b);
}
/* This value would be 3 if input values are exact */
-static int index_dot = 5;
+constexpr int index_dot = 5;
static double supremum_orient3d(const double3 &a,
const double3 &b,
const double3 &c,
const double3 &d)
{
- double3 abs_a{fabs(a[0]), fabs(a[1]), fabs(a[2])};
- double3 abs_b{fabs(b[0]), fabs(b[1]), fabs(b[2])};
- double3 abs_c{fabs(c[0]), fabs(c[1]), fabs(c[2])};
- double3 abs_d{fabs(d[0]), fabs(d[1]), fabs(d[2])};
+ 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];
@@ -1242,14 +1242,14 @@ static double supremum_orient3d(const double3 &a,
}
/* This value would be 8 if the input values are exact. */
-static int index_orient3d = 11;
+constexpr int index_orient3d = 11;
/* 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 fliter_orient3d(const double3 &a, const double3 &b, const double3 &c, const double3 &d)
+static int filter_orient3d(const double3 &a, const double3 &b, const double3 &c, const double3 &d)
{
double o3dfast = double3::orient3d_fast(a, b, c, d);
if (o3dfast == 0.0) {
@@ -1269,7 +1269,7 @@ static int fliter_orient3d(const double3 &a, const double3 &b, const double3 &c,
*/
static int filter_tri_plane_vert_orient3d(const Face &tri, Vertp v)
{
- return fliter_orient3d(tri[0]->co, tri[1]->co, tri[2]->co, v->co);
+ return filter_orient3d(tri[0]->co, tri[1]->co, tri[2]->co, v->co);
}
/* Are vectors a and b parallel or nearly parallel?
@@ -1308,6 +1308,34 @@ static bool dot_must_be_positive(const double3 &a, const double3 &b)
return false;
}
+/* 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.
+ * If the answer is 0, we are unsure about which side of the plane (or if it is on the plane).
+ * In exact arithmetic, the answer is just sgn(dot(p - plane_p, plane_no)).
+ */
+
+/* This would be 5 if inputs are exact. */
+constexpr int index_plane_side = 7;
+
+static int filter_plane_side(const double3 &p,
+ const double3 &plane_p,
+ const double3 &plane_no,
+ const double3 &abs_p,
+ const double3 &abs_plane_p,
+ const double3 &abs_plane_no)
+{
+ double d = double3::dot(p - plane_p, plane_no);
+ if (d == 0.0) {
+ return 0;
+ }
+ 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) {
+ return d > 0 ? 1 : -1;
+ }
+ return 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
@@ -1589,7 +1617,7 @@ static ITT_value intersect_tri_tri(const Mesh &tm, uint t1, uint t2)
{
constexpr int dbg_level = 0;
#ifdef PERFDEBUG
- incperfcount(0);
+ incperfcount(2); /* Intersect_tri_tri calls. */
#endif
const Face &tri1 = *tm.face(t1);
const Face &tri2 = *tm.face(t2);
@@ -1609,7 +1637,57 @@ static ITT_value intersect_tri_tri(const Mesh &tm, uint t1, uint t2)
std::cout << " r2 = " << vr2 << "\n";
}
- /* TODO: try doing intersect with double arithmetic first, with error bounds. */
+ /* Get signs of t1's vertices' distances to plane of t2 and vice versa. */
+
+ /* Try first getting signs with double arithmetic, with error bounds.
+ * If the signs calculated in this section are not 0, they are the same
+ * as what they would be using exact arithmetic.
+ */
+ const double3 &d_p1 = vp1->co;
+ const double3 &d_q1 = vq1->co;
+ const double3 &d_r1 = vr1->co;
+ const double3 &d_p2 = vp2->co;
+ const double3 &d_q2 = vq2->co;
+ const double3 &d_r2 = vr2->co;
+ const double3 &d_n2 = tri2.plane.norm;
+
+ const double3 &abs_d_p1 = double3::abs(d_p1);
+ const double3 &abs_d_q1 = double3::abs(d_q1);
+ const double3 &abs_d_r1 = double3::abs(d_r1);
+ const double3 &abs_d_r2 = double3::abs(d_r2);
+ const double3 &abs_d_n2 = double3::abs(d_n2);
+
+ int sp1 = filter_plane_side(d_p1, d_r2, d_n2, abs_d_p1, abs_d_r2, abs_d_n2);
+ int sq1 = filter_plane_side(d_q1, d_r2, d_n2, abs_d_q1, abs_d_r2, abs_d_n2);
+ int sr1 = filter_plane_side(d_r1, d_r2, d_n2, abs_d_r1, abs_d_r2, abs_d_n2);
+ if ((sp1 > 0 && sq1 > 0 && sr1 > 0) || (sp1 < 0 && sq1 < 0 && sr1 < 0)) {
+#ifdef PERFDEBUG
+ incperfcount(3); /* Tri tri intersects decided by filter plane tests. */
+#endif
+ if (dbg_level > 0) {
+ std::cout << "no intersection, all t1's verts above or below t2\n";
+ }
+ return ITT_value(INONE);
+ }
+
+ const double3 &d_n1 = tri1.plane.norm;
+ const double3 &abs_d_p2 = double3::abs(d_p2);
+ const double3 &abs_d_q2 = double3::abs(d_q2);
+ const double3 &abs_d_n1 = double3::abs(d_n1);
+
+ int sp2 = filter_plane_side(d_p2, d_r1
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list