[Bf-blender-cvs] [6f63417b500] master: Fix T84493 et al: New Boolean on Suzanne.
Howard Trickey
noreply at git.blender.org
Sun Feb 7 17:33:07 CET 2021
Commit: 6f63417b500d0893e89fef1ecddb9ff345322e96
Author: Howard Trickey
Date: Sun Feb 7 11:25:07 2021 -0500
Branches: master
https://developer.blender.org/rB6f63417b500d0893e89fef1ecddb9ff345322e96
Fix T84493 et al: New Boolean on Suzanne.
While Boolean is not guaranteed to work if the operands are not
volume-enclosing (technically: PWN - piecewise constant winding number),
it needs to do something in those cases. This change makes
more cases meet user expectations in T84493, T64544, T83403,
T82642 (though very slow on that one).
The original new boolean code used "generalized winding number"
for this fallback; replaced this with code that uses raycasting.
Raycasting would have been faster, but for unfortunately also
switchd to per-triangle tests rather than per-patch tests since
it is possible (e.g., with Suzanne) to have patches that are
both inside and outside the other shape. That can make it much
slower in some cases, sadly.
===================================================================
M source/blender/blenlib/intern/mesh_boolean.cc
M source/blender/bmesh/tools/bmesh_bevel.c
===================================================================
diff --git a/source/blender/blenlib/intern/mesh_boolean.cc b/source/blender/blenlib/intern/mesh_boolean.cc
index 88d90a7816f..85c64c3c855 100644
--- a/source/blender/blenlib/intern/mesh_boolean.cc
+++ b/source/blender/blenlib/intern/mesh_boolean.cc
@@ -27,10 +27,14 @@
# include "BLI_array.hh"
# include "BLI_assert.h"
# include "BLI_delaunay_2d.h"
+# include "BLI_double3.hh"
+# include "BLI_float3.hh"
# include "BLI_hash.hh"
+# include "BLI_kdopbvh.h"
# include "BLI_map.hh"
# include "BLI_math.h"
# include "BLI_math_boolean.hh"
+# include "BLI_math_geom.h"
# include "BLI_math_mpq.hh"
# include "BLI_mesh_intersect.hh"
# include "BLI_mpq3.hh"
@@ -45,6 +49,7 @@
# include "BLI_mesh_boolean.hh"
// # define PERFDEBUG
+
namespace blender::meshintersect {
/**
@@ -2328,159 +2333,216 @@ static const char *bool_optype_name(BoolOpType op)
}
}
-static mpq3 calc_point_inside_tri(const Face &tri)
+static double3 calc_point_inside_tri_db(const Face &tri)
{
const Vert *v0 = tri.vert[0];
const Vert *v1 = tri.vert[1];
const Vert *v2 = tri.vert[2];
- mpq3 ans = v0->co_exact / 3 + v1->co_exact / 3 + v2->co_exact / 3;
+ double3 ans = v0->co / 3 + v1->co / 3 + v2->co / 3;
return ans;
}
+class InsideShapeTestData {
+ public:
+ const IMesh &tm;
+ std::function<int(int)> shape_fn;
+ int nshapes;
+ /* A per-shape vector of parity of hits of that shape. */
+ Array<int> hit_parity;
+
+ InsideShapeTestData(const IMesh &tm, std::function<int(int)> shape_fn, int nshapes)
+ : tm(tm), shape_fn(shape_fn), nshapes(nshapes)
+ {
+ }
+};
-/**
- * Return the Generalized Winding Number of point \a testp with respect to the
- * volume implied by the faces for which shape_fn returns the value shape.
- * See "Robust Inside-Outside Segmentation using Generalized Winding Numbers"
- * by Jacobson, Kavan, and Sorkine-Hornung.
- * This is like a winding number in that if it is positive, the point
- * is inside the volume. But it is tolerant of not-completely-watertight
- * volumes, still doing a passable job of classifying inside/outside
- * as we intuitively understand that to mean.
- *
- * TOOD: speed up this calculation using the hierarchical algorithm in that paper.
- */
-static double generalized_winding_number(const IMesh &tm,
- std::function<int(int)> shape_fn,
- const double3 &testp,
- int shape)
+static void inside_shape_callback(void *userdata,
+ int index,
+ const BVHTreeRay *ray,
+ BVHTreeRayHit *UNUSED(hit))
{
- constexpr int dbg_level = 0;
+ const int dbg_level = 0;
if (dbg_level > 0) {
- std::cout << "GENERALIZED_WINDING_NUMBER testp = " << testp << ", shape = " << shape << "\n";
+ std::cout << "inside_shape_callback, index = " << index << "\n";
}
- double gwn = 0;
- for (int t : tm.face_index_range()) {
- const Face *f = tm.face(t);
- const Face &tri = *f;
- if (shape_fn(tri.orig) == shape) {
- if (dbg_level > 0) {
- std::cout << "accumulate for tri t = " << t << " = " << f << "\n";
- }
- const Vert *v0 = tri.vert[0];
- const Vert *v1 = tri.vert[1];
- const Vert *v2 = tri.vert[2];
- double3 a = v0->co - testp;
- double3 b = v1->co - testp;
- double3 c = v2->co - testp;
- /* Calculate the solid angle of abc relative to origin.
- * See "The Solid Angle of a Plane Triangle" by Oosterom and Strackee
- * for the derivation of the formula. */
- double alen = a.length();
- double blen = b.length();
- double clen = c.length();
- double3 bxc = double3::cross_high_precision(b, c);
- double num = double3::dot(a, bxc);
- double denom = alen * blen * clen + double3::dot(a, b) * clen + double3::dot(a, c) * blen +
- double3::dot(b, c) * alen;
- if (denom == 0.0) {
- if (dbg_level > 0) {
- std::cout << "denom == 0, skipping this tri\n";
- }
- continue;
- }
- double x = atan2(num, denom);
- double fgwn = 2.0 * x;
- if (dbg_level > 0) {
- std::cout << "tri contributes " << fgwn << "\n";
- }
- gwn += fgwn;
- }
+ InsideShapeTestData *data = static_cast<InsideShapeTestData *>(userdata);
+ const Face &tri = *data->tm.face(index);
+ int shape = data->shape_fn(tri.orig);
+ if (shape == -1) {
+ return;
+ }
+ float dist;
+ float fv0[3];
+ float fv1[3];
+ float fv2[3];
+ for (int i = 0; i < 3; ++i) {
+ fv0[i] = float(tri.vert[0]->co[i]);
+ fv1[i] = float(tri.vert[1]->co[i]);
+ fv2[i] = float(tri.vert[2]->co[i]);
}
- gwn = gwn / (M_PI * 4.0);
if (dbg_level > 0) {
- std::cout << "final gwn = " << gwn << "\n";
+ std::cout << " fv0=(" << fv0[0] << "," << fv0[1] << "," << fv0[2] << ")\n";
+ std::cout << " fv1=(" << fv1[0] << "," << fv1[1] << "," << fv1[2] << ")\n";
+ std::cout << " fv2=(" << fv2[0] << "," << fv2[1] << "," << fv2[2] << ")\n";
+ }
+ if (isect_ray_tri_epsilon_v3(
+ ray->origin, ray->direction, fv0, fv1, fv2, &dist, NULL, FLT_EPSILON)) {
+ /* Count parity as +1 if ray is in the same direction as tri's normal,
+ * and -1 if the directions are opposite. */
+ double3 o_db{double(ray->origin[0]), double(ray->origin[1]), double(ray->origin[2])};
+ int parity = orient3d(tri.vert[0]->co, tri.vert[1]->co, tri.vert[2]->co, o_db);
+ if (dbg_level > 0) {
+ std::cout << "origin at " << o_db << ", parity = " << parity << "\n";
+ }
+ data->hit_parity[shape] += parity;
}
- return gwn;
}
/**
- * Return true if point \a testp is inside the volume implied by the
- * faces for which the shape_fn returns the value shape.
- * If \a high_confidence is true then we want a higher degree
- * of "insideness" than if it is false.
+ * Test the triangle with index \a t_index to see which shapes it is inside,
+ * and fill in \a in_shape with a confidence value between 0 and 1 that says
+ * how likely we think it is that it is inside.
+ * This is done by casting some rays from just on the positive side of a test
+ * face in various directions and summing the parity of crossing faces of each face.
+ * The BVHtree \a tree contains all the triangles of \a tm and can be used for
+ * fast raycasting.
*/
-static bool point_is_inside_shape(const IMesh &tm,
- std::function<int(int)> shape_fn,
- const double3 &testp,
- int shape,
- bool high_confidence)
+static void test_tri_inside_shapes(const IMesh &tm,
+ std::function<int(int)> shape_fn,
+ int nshapes,
+ int test_t_index,
+ BVHTree *tree,
+ Array<float> &in_shape)
{
- double gwn = generalized_winding_number(tm, shape_fn, testp, shape);
- /* Due to floating point error, an outside point should get a value
- * of zero for gwn, but may have a very slightly positive value instead.
- * It is not important to get this epsilon very small, because practical
- * cases of interest will have gwn at least 0.2 if it is not zero. */
- if (high_confidence) {
- return (gwn > 0.9);
+ const int dbg_level = 0;
+ if (dbg_level > 0) {
+ std::cout << "test_point_inside_shapes, t_index = " << test_t_index << "\n";
+ }
+ Face &tri_test = *tm.face(test_t_index);
+ int shape = shape_fn(tri_test.orig);
+ if (shape == -1) {
+ in_shape.fill(0.0f);
+ return;
+ }
+ double3 test_point = calc_point_inside_tri_db(tri_test);
+ /* Offset the test point a tiny bit in the tri_test normal direcction. */
+ tri_test.populate_plane(false);
+ double3 norm = tri_test.plane->norm.normalized();
+ const double offset_amount = 1e-5;
+ double3 offset_test_point = test_point + offset_amount * norm;
+ if (dbg_level > 0) {
+ std::cout << "test tri is in shape " << shape << "\n";
+ std::cout << "test point = " << test_point << "\n";
+ std::cout << "offset_test_point = " << offset_test_point << "\n";
+ }
+ /* Try six test rays almost along orthogonal axes.
+ * Perturb their directions slightly to make it less likely to hit a seam.
+ * Raycast assumes they have unit length, so use r1 near 1 and
+ * ra near 0.5, and rb near .01, but normalized so sqrt(r1^2 + ra^2 + rb^2) == 1. */
+ constexpr int num_rays = 6;
+ constexpr float r1 = 0.9987025295199663f;
+ constexpr float ra = 0.04993512647599832f;
+ constexpr float rb = 0.009987025295199663f;
+ const float test_rays[num_rays][3] = {
+ {r1, ra, rb}, {-r1, -ra, -rb}, {rb, r1, ra}, {-rb, -r1, -ra}, {ra, rb, r1}, {-ra, -rb, -r1}};
+ InsideShapeTestData data(tm, shape_fn, nshapes);
+ data.hit_parity = Array<int>(nshapes, 0);
+ Array<int> count_insides(nshapes, 0);
+ const float co[3] = {
+ float(offset_test_point[0]), float(offset_test_point[1]), float(offset_test_point[2])};
+ for (int i = 0; i < num_rays; ++i) {
+ if (dbg_level > 0) {
+ std::cout << "shoot ray " << i << "(" << test_rays[i][0] << "," << test_rays[i][1] << ","
+ << test_rays[i][2] << ")\n";
+ }
+ BLI_bvhtree_ray_cast_all(tree, co, test_rays[i], 0.0f, FLT_MAX, inside_shape_callback, &data);
+ if (dbg_level > 0) {
+ std::cout << "ray " << i << " result:";
+ for (int j = 0; j < nshapes; ++j) {
+ std::cout << " " << data.hit_parity[j];
+ }
+ std::cout << "\n";
+ }
+ for (int j = 0; j < nshapes; ++j) {
+ if (j != shape && data.hit_parity[j] > 0) {
+ ++count_insides[j];
+ }
+ }
+ data.hit_parity.fill(0);
+ }
+ for (int j = 0; j < nshapes; ++j) {
+ if (j == shape) {
+ in_shape[j] = 1.0f; /* Let's say a shape is always inside itself. */
+ }
+ else {
+ in_shape[j] = float(count_insides[j]) / float(num_rays);
+ }
+ if (dbg_level > 0) {
+ std::cout << "shape " << j << " inside = " << in_shape[j] << "\n";
+ }
}
-
- return (gwn > 0.01);
}
/**
- * Use the Gen
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list