diff options
| -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) { |