[Bf-blender-cvs] [df8cc5662b9] master: Improve speed of Constrained Delaunay Triangulation with exact arith.
Howard Trickey
noreply at git.blender.org
Sat Nov 21 17:55:36 CET 2020
Commit: df8cc5662b9a03fb0cbec05bb9f9bad103b8870b
Author: Howard Trickey
Date: Sat Nov 21 11:55:14 2020 -0500
Branches: master
https://developer.blender.org/rBdf8cc5662b9a03fb0cbec05bb9f9bad103b8870b
Improve speed of Constrained Delaunay Triangulation with exact arith.
By using floating point filters, the speed improves by a factor of 2 to 10.
This will help speed up some cases of the Exact Boolean modifier.
Changed the interface of mpq2::isect_seg_seg to not return mu, as it was
not needed and not calculating it saved 15% time.
===================================================================
M source/blender/blenlib/BLI_double2.hh
M source/blender/blenlib/BLI_mpq2.hh
M source/blender/blenlib/intern/delaunay_2d.cc
M source/blender/blenlib/intern/math_vec.cc
M source/blender/blenlib/tests/BLI_delaunay_2d_test.cc
===================================================================
diff --git a/source/blender/blenlib/BLI_double2.hh b/source/blender/blenlib/BLI_double2.hh
index 313afd06d1b..621ac4d01fc 100644
--- a/source/blender/blenlib/BLI_double2.hh
+++ b/source/blender/blenlib/BLI_double2.hh
@@ -131,7 +131,6 @@ struct double2 {
LINE_LINE_CROSS = 2,
} kind;
double lambda;
- double mu;
};
static isect_result isect_seg_seg(const double2 &v1,
diff --git a/source/blender/blenlib/BLI_mpq2.hh b/source/blender/blenlib/BLI_mpq2.hh
index 6261b50466b..40e227019ce 100644
--- a/source/blender/blenlib/BLI_mpq2.hh
+++ b/source/blender/blenlib/BLI_mpq2.hh
@@ -167,7 +167,6 @@ struct mpq2 {
LINE_LINE_CROSS = 2,
} kind;
mpq_class lambda;
- mpq_class mu;
};
static isect_result isect_seg_seg(const mpq2 &v1,
diff --git a/source/blender/blenlib/intern/delaunay_2d.cc b/source/blender/blenlib/intern/delaunay_2d.cc
index daa006fd5cc..ae161fc2d1c 100644
--- a/source/blender/blenlib/intern/delaunay_2d.cc
+++ b/source/blender/blenlib/intern/delaunay_2d.cc
@@ -117,11 +117,80 @@ template<typename T> inline SymEdge<T> *prev(const SymEdge<T> *se)
return se->rot->next->rot;
}
-template<typename Arith_t> struct CDTVert {
+/** A coordinate class with extra information for fast filtered orient tests. */
+
+template<typename T> struct FatCo {
+ vec2<T> exact;
+ vec2<double> approx;
+ vec2<double> abs_approx;
+
+ FatCo();
+ FatCo(const vec2<mpq_class> &v);
+ FatCo(const vec2<double> &v);
+};
+
+#ifdef WITH_GMP
+template<> struct FatCo<mpq_class> {
+ vec2<mpq_class> exact;
+ vec2<double> approx;
+ vec2<double> abs_approx;
+
+ FatCo()
+ : exact(vec2<mpq_class>(0, 0)), approx(vec2<double>(0, 0)), abs_approx(vec2<double>(0, 0))
+ {
+ }
+
+ FatCo(const vec2<mpq_class> &v)
+ {
+ exact = v;
+ approx = vec2<double>(v.x.get_d(), v.y.get_d());
+ abs_approx = vec2<double>(fabs(approx.x), fabs(approx.y));
+ }
+
+ FatCo(const vec2<double> &v)
+ {
+ exact = vec2<mpq_class>(v.x, v.y);
+ approx = v;
+ abs_approx = vec2<double>(fabs(approx.x), fabs(approx.y));
+ }
+};
+#endif
+
+template<> struct FatCo<double> {
+ vec2<double> exact;
+ vec2<double> approx;
+ vec2<double> abs_approx;
+
+ FatCo() : exact(vec2<double>(0, 0)), approx(vec2<double>(0, 0)), abs_approx(vec2<double>(0, 0))
+ {
+ }
+
+ FatCo(const vec2<mpq_class> &v)
+ {
+ exact = vec2<double>(v.x.get_d(), v.y.get_d());
+ approx = exact;
+ abs_approx = vec2<double>(fabs(approx.x), fabs(approx.y));
+ }
+
+ FatCo(const vec2<double> &v)
+ {
+ exact = v;
+ approx = v;
+ abs_approx = vec2<double>(fabs(approx.x), fabs(approx.y));
+ }
+};
+
+template<typename T> std::ostream &operator<<(std::ostream &stream, const FatCo<T> &co)
+{
+ stream << "(" << co.approx.x << ", " << co.approx.y << ")";
+ return stream;
+}
+
+template<typename T> struct CDTVert {
/** Coordinate. */
- vec2<Arith_t> co;
+ FatCo<T> co;
/** Some edge attached to it. */
- SymEdge<Arith_t> *symedge{nullptr};
+ SymEdge<T> *symedge{nullptr};
/** List of corresponding vertex input ids. */
LinkNode *input_ids{nullptr};
/** Index into array that #CDTArrangement keeps. */
@@ -132,7 +201,7 @@ template<typename Arith_t> struct CDTVert {
int visit_index{0};
CDTVert() = default;
- explicit CDTVert(const vec2<Arith_t> &pt);
+ explicit CDTVert(const vec2<T> &pt);
};
template<typename Arith_t> struct CDTEdge {
@@ -431,7 +500,7 @@ template<typename T> void cdt_draw(const std::string &label, const CDTArrangemen
vec2<double> vmax(-DBL_MAX, -DBL_MAX);
for (const CDTVert<T> *v : cdt.verts) {
for (int i = 0; i < 2; ++i) {
- double dvi = math_to_double(v->co[i]);
+ double dvi = v->co.approx[i];
if (dvi < vmin[i]) {
vmin[i] = dvi;
}
@@ -457,8 +526,8 @@ template<typename T> void cdt_draw(const std::string &label, const CDTArrangemen
}
double scale = view_width / width;
-# define SX(x) ((math_to_double(x) - minx) * scale)
-# define SY(y) ((maxy - math_to_double(y)) * scale)
+# define SX(x) (((x)-minx) * scale)
+# define SY(y) ((maxy - (y)) * scale)
std::ofstream f;
if (append) {
@@ -485,8 +554,8 @@ template<typename T> void cdt_draw(const std::string &label, const CDTArrangemen
}
const CDTVert<T> *u = e->symedges[0].vert;
const CDTVert<T> *v = e->symedges[1].vert;
- const vec2<T> &uco = u->co;
- const vec2<T> &vco = v->co;
+ const vec2<double> &uco = u->co.approx;
+ const vec2<double> &vco = v->co.approx;
int strokew = e->input_ids == nullptr ? thin_line : thick_line;
f << "<line fill=\"none\" stroke=\"black\" stroke-width=\"" << strokew << "\" x1=\""
<< SX(uco[0]) << "\" y1=\"" << SY(uco[1]) << "\" x2=\"" << SX(vco[0]) << "\" y2=\""
@@ -503,13 +572,13 @@ template<typename T> void cdt_draw(const std::string &label, const CDTArrangemen
int i = 0;
for (const CDTVert<T> *v : cdt.verts) {
- f << "<circle fill=\"black\" cx=\"" << SX(v->co[0]) << "\" cy=\"" << SY(v->co[1]) << "\" r=\""
- << vert_radius << "\">\n";
- f << " <title>[" << i << "]" << v->co << "</title>\n";
+ f << "<circle fill=\"black\" cx=\"" << SX(v->co.approx[0]) << "\" cy=\"" << SY(v->co.approx[1])
+ << "\" r=\"" << vert_radius << "\">\n";
+ f << " <title>[" << i << "]" << v->co.approx << "</title>\n";
f << "</circle>\n";
if (draw_vert_labels) {
- f << "<text x=\"" << SX(v->co[0]) + vert_radius << "\" y=\"" << SY(v->co[1]) - vert_radius
- << "\" font-size=\"small\">[" << i << "]</text>\n";
+ f << "<text x=\"" << SX(v->co.approx[0]) + vert_radius << "\" y=\""
+ << SY(v->co.approx[1]) - vert_radius << "\" font-size=\"small\">[" << i << "]</text>\n";
}
++i;
}
@@ -520,31 +589,188 @@ template<typename T> void cdt_draw(const std::string &label, const CDTArrangemen
}
#endif
+/**
+ * A filtered version of orient2d, which will usually be much faster when using exact arithmetic.
+ * See EXACT GEOMETRIC COMPUTATION USING CASCADING, by Burnikel, Funke, and Seel.
+ */
+template<typename T>
+static int filtered_orient2d(const FatCo<T> &a, const FatCo<T> &b, const FatCo<T> &c);
+
+#ifdef WITH_GMP
+template<>
+int filtered_orient2d<mpq_class>(const FatCo<mpq_class> &a,
+ const FatCo<mpq_class> &b,
+ const FatCo<mpq_class> &c)
+{
+ double det = (a.approx[0] - c.approx[0]) * (b.approx[1] - c.approx[1]) -
+ (a.approx[1] - c.approx[1]) * (b.approx[0] - c.approx[0]);
+ double supremum = (a.abs_approx[0] + c.abs_approx[0]) * (b.abs_approx[1] + c.abs_approx[1]) +
+ (a.abs_approx[1] + c.abs_approx[1]) * (b.abs_approx[0] + c.abs_approx[0]);
+ constexpr double index_orient2d = 6;
+ double err_bound = supremum * index_orient2d * DBL_EPSILON;
+ if (fabs(det) > err_bound) {
+ return det > 0 ? 1 : -1;
+ }
+ return orient2d(a.exact, b.exact, c.exact);
+}
+#endif
+
+template<>
+int filtered_orient2d<double>(const FatCo<double> &a,
+ const FatCo<double> &b,
+ const FatCo<double> &c)
+{
+ return orient2d(a.approx, b.approx, c.approx);
+}
+
+/**
+ * A filtered version of incircle.
+ */
+template<typename T>
+static int filtered_incircle(const FatCo<T> &a,
+ const FatCo<T> &b,
+ const FatCo<T> &c,
+ const FatCo<T> &d);
+
+#ifdef WITH_GMP
+template<>
+int filtered_incircle<mpq_class>(const FatCo<mpq_class> &a,
+ const FatCo<mpq_class> &b,
+ const FatCo<mpq_class> &c,
+ const FatCo<mpq_class> &d)
+{
+ double adx = a.approx[0] - d.approx[0];
+ double bdx = b.approx[0] - d.approx[0];
+ double cdx = c.approx[0] - d.approx[0];
+ double ady = a.approx[1] - d.approx[1];
+ double bdy = b.approx[1] - d.approx[1];
+ double cdy = c.approx[1] - d.approx[1];
+ double bdxcdy = bdx * cdy;
+ double cdxbdy = cdx * bdy;
+ double alift = adx * adx + ady * ady;
+ double cdxady = cdx * ady;
+ double adxcdy = adx * cdy;
+ double blift = bdx * bdx + bdy * bdy;
+ double adxbdy = adx * bdy;
+ double bdxady = bdx * ady;
+ double clift = cdx * cdx + cdy * cdy;
+ double det = alift * (bdxcdy - cdxbdy) + blift * (cdxady - adxcdy) + clift * (adxbdy - bdxady);
+
+ double sup_adx = a.abs_approx[0] + d.abs_approx[0]; /* index 2. */
+ double sup_bdx = b.abs_approx[0] + d.abs_approx[0];
+ double sup_cdx = c.abs_approx[0] + d.abs_approx[0];
+ double sup_ady = a.abs_approx[1] + d.abs_approx[1];
+ double sup_bdy = b.abs_approx[1] + d.abs_approx[1];
+ double sup_cdy = c.abs_approx[1] + d.abs_approx[1];
+ double sup_bdxcdy = sup_bdx * sup_cdy; /* index 5. */
+ double sup_cdxbdy = sup_cdx * sup_bdy;
+ double sup_alift = sup_adx * sup_adx + sup_ady * sup_ady; /* index 6. */
+ double sup_cdxady = sup_cdx * sup_ady;
+ double sup_adxcdy = sup_adx * sup_cdy;
+ double sup_blift = sup_bdx * sup_bdx + sup_bdy * sup_bdy;
+ double sup_adxbdy = sup_adx * sup_bdy;
+ double sup_bdxady = sup_bdx * sup_ady;
+ double sup_clift = sup_cdx * sup_cdx + sup_cdy * sup_cdy;
+ double sup_det = sup_alift * (sup_bdxcdy + sup_cdxbdy) + sup_blift * (sup_cdxady + sup_adxcdy) +
+ sup_clift * (sup_adxbdy + sup_bdxady);
+ int index = 14;
+ double err_bound = sup_det * index * DBL_EPSILON;
+ if (fabs(det) > err_bound) {
+ return det < 0.0 ? -1 : (det > 0.0 ? 1 : 0);
+ }
+ return incircle(a.exact, b.exact, c.exact, d.exact);
+}
+#endif
+
+template<>
+int filtered_incircle<double>(const FatCo<double> &a,
+ const FatCo<double> &b,
+ const FatCo<double> &c,
+ const FatCo<double> &d)
+{
+ return incircle(a.approx, b.approx, c.approx, d.approx);
+}
+
/**
* Return true if `a -- b -- c` are in that order, assuming they are on a straight line according
* to #orient2d and we know the order is either `abc` or `bac`.
* This means `ab . ac` and `bc . ac` must both be non-negative.
+ * Use filtering to speed this up when using exact arithmetic.
*/
-template<typename T> bool in_line(const vec2<T> &a
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list