diff options
author | Erick Abrahammson <ecke101@gmail.com> | 2021-06-01 00:03:48 +0300 |
---|---|---|
committer | Howard Trickey <howard.trickey@gmail.com> | 2021-06-01 00:03:48 +0300 |
commit | 8180d478e1aefbbe538bd54b42dda388b482abf5 (patch) | |
tree | 4574fb0daaaa720a04ed37b9cff119932589f84a /source/blender | |
parent | 73967e2047b7d46168355a7158cd1c40d20afc4e (diff) |
Speedup exact boolean by avoiding some mallocs and frees.
This is from patch D11432 from Erik Abrahamsson. He found that
in some mpq3 functions called frequently from loops, passing in
buffers for termporary mpq3 values can save substantial time.
On my machine, his example in that patch went from 9.48s to 7.50s
for the boolean part of the calculation. On his machine, a running
time went from 17s to 10.3s.
Diffstat (limited to 'source/blender')
-rw-r--r-- | source/blender/blenlib/BLI_mpq3.hh | 16 | ||||
-rw-r--r-- | source/blender/blenlib/intern/mesh_boolean.cc | 94 | ||||
-rw-r--r-- | source/blender/blenlib/intern/mesh_intersect.cc | 91 |
3 files changed, 146 insertions, 55 deletions
diff --git a/source/blender/blenlib/BLI_mpq3.hh b/source/blender/blenlib/BLI_mpq3.hh index fb5e2b61cdb..b9eda2ad7e1 100644 --- a/source/blender/blenlib/BLI_mpq3.hh +++ b/source/blender/blenlib/BLI_mpq3.hh @@ -218,6 +218,15 @@ struct mpq3 { return a.x * b.x + a.y * b.y + a.z * b.z; } + static mpq_class dot_with_buffer(const mpq3 &a, const mpq3 &b, mpq3 &buffer) + { + buffer = a; + buffer *= b; + buffer.x += buffer.y; + buffer.x += buffer.z; + return buffer.x; + } + static mpq3 cross(const mpq3 &a, const mpq3 &b) { return mpq3(a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]); @@ -246,6 +255,13 @@ struct mpq3 { return mpq3::dot(diff, diff); } + static mpq_class distance_squared_with_buffer(const mpq3 &a, const mpq3 &b, mpq3 &buffer) + { + buffer = a; + buffer -= b; + return mpq3::dot(buffer, buffer); + } + static mpq3 interpolate(const mpq3 &a, const mpq3 &b, mpq_class t) { mpq_class s = 1 - t; diff --git a/source/blender/blenlib/intern/mesh_boolean.cc b/source/blender/blenlib/intern/mesh_boolean.cc index 431720761bc..cc27e9238c3 100644 --- a/source/blender/blenlib/intern/mesh_boolean.cc +++ b/source/blender/blenlib/intern/mesh_boolean.cc @@ -1695,9 +1695,24 @@ static int find_containing_cell(const Vert *v, * If the closest point is on an edge, return 0, 1, or 2 * for edges ab, bc, or ca in *r_edge; else -1. * (Adapted from #closest_on_tri_to_point_v3()). + * The arguments ab, ac, ..., r are used as temporaries + * in this routine. Passing them in from the caller can + * avoid many allocs and frees of temporary mpq3 values + * and the mpq_class values within them. */ -static mpq_class closest_on_tri_to_point( - const mpq3 &p, const mpq3 &a, const mpq3 &b, const mpq3 &c, int *r_edge, int *r_vert) +static mpq_class closest_on_tri_to_point(const mpq3 &p, + const mpq3 &a, + const mpq3 &b, + const mpq3 &c, + mpq3 &ab, + mpq3 &ac, + mpq3 &ap, + mpq3 &bp, + mpq3 &cp, + mpq3 &m, + mpq3 &r, + int *r_edge, + int *r_vert) { constexpr int dbg_level = 0; if (dbg_level > 0) { @@ -1705,11 +1720,15 @@ static mpq_class closest_on_tri_to_point( std::cout << " a = " << a << ", b = " << b << ", c = " << c << "\n"; } /* Check if p in vertex region outside a. */ - mpq3 ab = b - a; - mpq3 ac = c - a; - mpq3 ap = p - a; - mpq_class d1 = mpq3::dot(ab, ap); - mpq_class d2 = mpq3::dot(ac, ap); + ab = b; + ab -= a; + ac = c; + ac -= a; + ap = p; + ap -= a; + + mpq_class d1 = mpq3::dot_with_buffer(ab, ap, m); + mpq_class d2 = mpq3::dot_with_buffer(ac, ap, m); if (d1 <= 0 && d2 <= 0) { /* Barycentric coordinates (1,0,0). */ *r_edge = -1; @@ -1717,12 +1736,13 @@ static mpq_class closest_on_tri_to_point( if (dbg_level > 0) { std::cout << " answer = a\n"; } - return mpq3::distance_squared(p, a); + return mpq3::distance_squared_with_buffer(p, a, m); } /* Check if p in vertex region outside b. */ - mpq3 bp = p - b; - mpq_class d3 = mpq3::dot(ab, bp); - mpq_class d4 = mpq3::dot(ac, bp); + bp = p; + bp -= b; + mpq_class d3 = mpq3::dot_with_buffer(ab, bp, m); + mpq_class d4 = mpq3::dot_with_buffer(ac, bp, m); if (d3 >= 0 && d4 <= d3) { /* Barycentric coordinates (0,1,0). */ *r_edge = -1; @@ -1730,25 +1750,28 @@ static mpq_class closest_on_tri_to_point( if (dbg_level > 0) { std::cout << " answer = b\n"; } - return mpq3::distance_squared(p, b); + return mpq3::distance_squared_with_buffer(p, b, m); } /* Check if p in region of ab. */ mpq_class vc = d1 * d4 - d3 * d2; if (vc <= 0 && d1 >= 0 && d3 <= 0) { mpq_class v = d1 / (d1 - d3); /* Barycentric coordinates (1-v,v,0). */ - mpq3 r = a + v * ab; + r = ab; + r *= v; + r += a; *r_vert = -1; *r_edge = 0; if (dbg_level > 0) { std::cout << " answer = on ab at " << r << "\n"; } - return mpq3::distance_squared(p, r); + return mpq3::distance_squared_with_buffer(p, r, m); } /* Check if p in vertex region outside c. */ - mpq3 cp = p - c; - mpq_class d5 = mpq3::dot(ab, cp); - mpq_class d6 = mpq3::dot(ac, cp); + cp = p; + cp -= c; + mpq_class d5 = mpq3::dot_with_buffer(ab, cp, m); + mpq_class d6 = mpq3::dot_with_buffer(ac, cp, m); if (d6 >= 0 && d5 <= d6) { /* Barycentric coordinates (0,0,1). */ *r_edge = -1; @@ -1756,49 +1779,54 @@ static mpq_class closest_on_tri_to_point( if (dbg_level > 0) { std::cout << " answer = c\n"; } - return mpq3::distance_squared(p, c); + return mpq3::distance_squared_with_buffer(p, c, m); } /* Check if p in edge region of ac. */ mpq_class vb = d5 * d2 - d1 * d6; if (vb <= 0 && d2 >= 0 && d6 <= 0) { mpq_class w = d2 / (d2 - d6); /* Barycentric coordinates (1-w,0,w). */ - mpq3 r = a + w * ac; + r = ac; + r *= w; + r += a; *r_vert = -1; *r_edge = 2; if (dbg_level > 0) { std::cout << " answer = on ac at " << r << "\n"; } - return mpq3::distance_squared(p, r); + return mpq3::distance_squared_with_buffer(p, r, m); } /* Check if p in edge region of bc. */ mpq_class va = d3 * d6 - d5 * d4; if (va <= 0 && (d4 - d3) >= 0 && (d5 - d6) >= 0) { mpq_class w = (d4 - d3) / ((d4 - d3) + (d5 - d6)); /* Barycentric coordinates (0,1-w,w). */ - mpq3 r = c - b; - r = w * r; - r = r + b; + r = c; + r -= b; + r *= w; + r += b; *r_vert = -1; *r_edge = 1; if (dbg_level > 0) { std::cout << " answer = on bc at " << r << "\n"; } - return mpq3::distance_squared(p, r); + return mpq3::distance_squared_with_buffer(p, r, m); } /* p inside face region. Compute barycentric coordinates (u,v,w). */ mpq_class denom = 1 / (va + vb + vc); mpq_class v = vb * denom; mpq_class w = vc * denom; - ac = w * ac; - mpq3 r = a + v * ab; - r = r + ac; + ac *= w; + r = ab; + r *= v; + r += a; + r += ac; *r_vert = -1; *r_edge = -1; if (dbg_level > 0) { std::cout << " answer = inside at " << r << "\n"; } - return mpq3::distance_squared(p, r); + return mpq3::distance_squared_with_buffer(p, r, m); } struct ComponentContainer { @@ -1837,6 +1865,9 @@ static Vector<ComponentContainer> find_component_containers(int comp, if (dbg_level > 0) { std::cout << "test vertex in comp: " << test_v << "\n"; } + + mpq3 buf[7]; + for (int comp_other : components.index_range()) { if (comp == comp_other) { continue; @@ -1861,6 +1892,13 @@ static Vector<ComponentContainer> find_component_containers(int comp, tri[0]->co_exact, tri[1]->co_exact, tri[2]->co_exact, + buf[0], + buf[1], + buf[2], + buf[3], + buf[4], + buf[5], + buf[6], &close_edge, &close_vert); if (dbg_level > 1) { diff --git a/source/blender/blenlib/intern/mesh_intersect.cc b/source/blender/blenlib/intern/mesh_intersect.cc index 636209883c3..97f856476c5 100644 --- a/source/blender/blenlib/intern/mesh_intersect.cc +++ b/source/blender/blenlib/intern/mesh_intersect.cc @@ -1165,13 +1165,19 @@ static int filter_plane_side(const double3 &p, * Assumes ab is not perpendicular to n. * This works because the ratio of the projections of ab and ac onto n is the same as * the ratio along the line ab of the intersection point to the whole of ab. + * The ab, ac, and dotbuf arguments are used as a temporaries; declaring them + * in the caller can avoid many allocs and frees of mpq3 and mpq_class structures. */ -static inline mpq3 tti_interp(const mpq3 &a, const mpq3 &b, const mpq3 &c, const mpq3 &n) -{ - mpq3 ab = a - b; - mpq_class den = mpq3::dot(ab, n); +static inline mpq3 tti_interp( + const mpq3 &a, const mpq3 &b, const mpq3 &c, const mpq3 &n, mpq3 &ab, mpq3 &ac, mpq3 &dotbuf) +{ + ab = a; + ab -= b; + ac = a; + ac -= c; + mpq_class den = mpq3::dot_with_buffer(ab, n, dotbuf); BLI_assert(den != 0); - mpq_class alpha = mpq3::dot(a - c, n) / den; + mpq_class alpha = mpq3::dot_with_buffer(ac, n, dotbuf) / den; return a - alpha * ab; } @@ -1179,11 +1185,28 @@ static inline mpq3 tti_interp(const mpq3 &a, const mpq3 &b, const mpq3 &c, const * Return +1, 0, -1 as a + ad is above, on, or below the oriented plane containing a, b, c in CCW * order. This is the same as -oriented(a, b, c, a + ad), but uses fewer arithmetic operations. * TODO: change arguments to `const Vert *` and use floating filters. + * The ba, ca, n, and dotbuf arguments are used as temporaries; declaring them + * in the caller can avoid many allocs and frees of mpq3 and mpq_class structures. */ -static inline int tti_above(const mpq3 &a, const mpq3 &b, const mpq3 &c, const mpq3 &ad) +static inline int tti_above(const mpq3 &a, + const mpq3 &b, + const mpq3 &c, + const mpq3 &ad, + mpq3 &ba, + mpq3 &ca, + mpq3 &n, + mpq3 &dotbuf) { - mpq3 n = mpq3::cross(b - a, c - a); - return sgn(mpq3::dot(ad, n)); + ba = b; + ba -= a; + ca = c; + ca -= a; + + n.x = ba.y * ca.z - ba.z * ca.y; + n.y = ba.z * ca.x - ba.x * ca.z; + n.z = ba.x * ca.y - ba.y * ca.x; + + return sgn(mpq3::dot_with_buffer(ad, n, dotbuf)); } /** @@ -1227,20 +1250,21 @@ static ITT_value itt_canon2(const mpq3 &p1, mpq3 p1p2 = p2 - p1; mpq3 intersect_1; mpq3 intersect_2; + mpq3 buf[4]; bool no_overlap = false; /* Top test in classification tree. */ - if (tti_above(p1, q1, r2, p1p2) > 0) { + if (tti_above(p1, q1, r2, p1p2, buf[0], buf[1], buf[2], buf[3]) > 0) { /* Middle right test in classification tree. */ - if (tti_above(p1, r1, r2, p1p2) <= 0) { + if (tti_above(p1, r1, r2, p1p2, buf[0], buf[1], buf[2], buf[3]) <= 0) { /* Bottom right test in classification tree. */ - if (tti_above(p1, r1, q2, p1p2) > 0) { + if (tti_above(p1, r1, q2, p1p2, buf[0], buf[1], buf[2], buf[3]) > 0) { /* Overlap is [k [i l] j]. */ if (dbg_level > 0) { std::cout << "overlap [k [i l] j]\n"; } /* i is intersect with p1r1. l is intersect with p2r2. */ - intersect_1 = tti_interp(p1, r1, p2, n2); - intersect_2 = tti_interp(p2, r2, p1, n1); + intersect_1 = tti_interp(p1, r1, p2, n2, buf[0], buf[1], buf[2]); + intersect_2 = tti_interp(p2, r2, p1, n1, buf[0], buf[1], buf[2]); } else { /* Overlap is [i [k l] j]. */ @@ -1248,8 +1272,8 @@ static ITT_value itt_canon2(const mpq3 &p1, std::cout << "overlap [i [k l] j]\n"; } /* k is intersect with p2q2. l is intersect is p2r2. */ - intersect_1 = tti_interp(p2, q2, p1, n1); - intersect_2 = tti_interp(p2, r2, p1, n1); + intersect_1 = tti_interp(p2, q2, p1, n1, buf[0], buf[1], buf[2]); + intersect_2 = tti_interp(p2, r2, p1, n1, buf[0], buf[1], buf[2]); } } else { @@ -1262,7 +1286,7 @@ static ITT_value itt_canon2(const mpq3 &p1, } else { /* Middle left test in classification tree. */ - if (tti_above(p1, q1, q2, p1p2) < 0) { + if (tti_above(p1, q1, q2, p1p2, buf[0], buf[1], buf[2], buf[3]) < 0) { /* No overlap: [i j] [k l]. */ if (dbg_level > 0) { std::cout << "no overlap: [i j] [k l]\n"; @@ -1271,14 +1295,14 @@ static ITT_value itt_canon2(const mpq3 &p1, } else { /* Bottom left test in classification tree. */ - if (tti_above(p1, r1, q2, p1p2) >= 0) { + if (tti_above(p1, r1, q2, p1p2, buf[0], buf[1], buf[2], buf[3]) >= 0) { /* Overlap is [k [i j] l]. */ if (dbg_level > 0) { std::cout << "overlap [k [i j] l]\n"; } /* i is intersect with p1r1. j is intersect with p1q1. */ - intersect_1 = tti_interp(p1, r1, p2, n2); - intersect_2 = tti_interp(p1, q1, p2, n2); + intersect_1 = tti_interp(p1, r1, p2, n2, buf[0], buf[1], buf[2]); + intersect_2 = tti_interp(p1, q1, p2, n2, buf[0], buf[1], buf[2]); } else { /* Overlap is [i [k j] l]. */ @@ -1286,8 +1310,8 @@ static ITT_value itt_canon2(const mpq3 &p1, std::cout << "overlap [i [k j] l]\n"; } /* k is intersect with p2q2. j is intersect with p1q1. */ - intersect_1 = tti_interp(p2, q2, p1, n1); - intersect_2 = tti_interp(p1, q1, p2, n2); + intersect_1 = tti_interp(p2, q2, p1, n1, buf[0], buf[1], buf[2]); + intersect_2 = tti_interp(p1, q1, p2, n2, buf[0], buf[1], buf[2]); } } } @@ -1438,6 +1462,7 @@ static ITT_value intersect_tri_tri(const IMesh &tm, int t1, int t2) return ITT_value(INONE); } + mpq3 buf[2]; const mpq3 &p1 = vp1->co_exact; const mpq3 &q1 = vq1->co_exact; const mpq3 &r1 = vr1->co_exact; @@ -1447,13 +1472,19 @@ static ITT_value intersect_tri_tri(const IMesh &tm, int t1, int t2) const mpq3 &n2 = tri2.plane->norm_exact; if (sp1 == 0) { - sp1 = sgn(mpq3::dot(p1 - r2, n2)); + buf[0] = p1; + buf[0] -= r2; + sp1 = sgn(mpq3::dot_with_buffer(buf[0], n2, buf[1])); } if (sq1 == 0) { - sq1 = sgn(mpq3::dot(q1 - r2, n2)); + buf[0] = q1; + buf[0] -= r2; + sq1 = sgn(mpq3::dot_with_buffer(buf[0], n2, buf[1])); } if (sr1 == 0) { - sr1 = sgn(mpq3::dot(r1 - r2, n2)); + buf[0] = r1; + buf[0] -= r2; + sr1 = sgn(mpq3::dot_with_buffer(buf[0], n2, buf[1])); } if (dbg_level > 1) { @@ -1473,13 +1504,19 @@ static ITT_value intersect_tri_tri(const IMesh &tm, int t1, int t2) /* Repeat for signs of t2's vertices with respect to plane of t1. */ const mpq3 &n1 = tri1.plane->norm_exact; if (sp2 == 0) { - sp2 = sgn(mpq3::dot(p2 - r1, n1)); + buf[0] = p2; + buf[0] -= r1; + sp2 = sgn(mpq3::dot_with_buffer(buf[0], n1, buf[1])); } if (sq2 == 0) { - sq2 = sgn(mpq3::dot(q2 - r1, n1)); + buf[0] = q2; + buf[0] -= r1; + sq2 = sgn(mpq3::dot_with_buffer(buf[0], n1, buf[1])); } if (sr2 == 0) { - sr2 = sgn(mpq3::dot(r2 - r1, n1)); + buf[0] = r2; + buf[0] -= r1; + sr2 = sgn(mpq3::dot_with_buffer(buf[0], n1, buf[1])); } if (dbg_level > 1) { |