diff options
author | Howard Trickey <howard.trickey@gmail.com> | 2020-11-21 19:55:14 +0300 |
---|---|---|
committer | Howard Trickey <howard.trickey@gmail.com> | 2020-11-21 19:55:14 +0300 |
commit | df8cc5662b9a03fb0cbec05bb9f9bad103b8870b (patch) | |
tree | 9e4b543ec5e10d629c8a354239636495e4fa3616 | |
parent | 38fe962d9542296e94b5881f45043ae5afe8e20e (diff) |
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.
-rw-r--r-- | source/blender/blenlib/BLI_double2.hh | 1 | ||||
-rw-r--r-- | source/blender/blenlib/BLI_mpq2.hh | 1 | ||||
-rw-r--r-- | source/blender/blenlib/intern/delaunay_2d.cc | 350 | ||||
-rw-r--r-- | source/blender/blenlib/intern/math_vec.cc | 16 | ||||
-rw-r--r-- | source/blender/blenlib/tests/BLI_delaunay_2d_test.cc | 3 |
5 files changed, 299 insertions, 72 deletions
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; } @@ -521,30 +590,187 @@ 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, const vec2<T> &b, const vec2<T> &c) +template<typename T> static bool in_line(const FatCo<T> &a, const FatCo<T> &b, const FatCo<T> &c); + +#ifdef WITH_GMP +template<> +bool in_line<mpq_class>(const FatCo<mpq_class> &a, + const FatCo<mpq_class> &b, + const FatCo<mpq_class> &c) +{ + vec2<double> ab = b.approx - a.approx; + vec2<double> bc = c.approx - b.approx; + vec2<double> ac = c.approx - a.approx; + vec2<double> supremum_ab = b.abs_approx + a.abs_approx; + vec2<double> supremum_bc = c.abs_approx + b.abs_approx; + vec2<double> supremum_ac = c.abs_approx + a.abs_approx; + double dot_ab_ac = ab.x * ac.x + ab.y * ac.y; + double supremum_dot_ab_ac = supremum_ab.x * supremum_ac.x + supremum_ab.y * supremum_ac.y; + constexpr double index = 6; + double err_bound = supremum_dot_ab_ac * index * DBL_EPSILON; + if (dot_ab_ac < -err_bound) { + return false; + } + double dot_bc_ac = bc.x * ac.x + bc.y * ac.y; + double supremum_dot_bc_ac = supremum_bc.x * supremum_ac.x + supremum_bc.y * supremum_ac.y; + err_bound = supremum_dot_bc_ac * index * DBL_EPSILON; + if (dot_bc_ac < -err_bound) { + return false; + } + vec2<mpq_class> exact_ab = b.exact - a.exact; + vec2<mpq_class> exact_ac = c.exact - a.exact; + if (vec2<mpq_class>::dot(exact_ab, exact_ac) < 0) { + return false; + } + vec2<mpq_class> exact_bc = c.exact - b.exact; + return vec2<mpq_class>::dot(exact_bc, exact_ac) >= 0; +} +#endif + +template<> +bool in_line<double>(const FatCo<double> &a, const FatCo<double> &b, const FatCo<double> &c) { - vec2<T> ab = b - a; - vec2<T> bc = c - b; - vec2<T> ac = c - a; - if (vec2<T>::dot(ab, ac) < 0) { + vec2<double> ab = b.approx - a.approx; + vec2<double> ac = c.approx - a.approx; + if (vec2<double>::dot(ab, ac) < 0) { return false; } - return vec2<T>::dot(bc, ac) >= 0; + vec2<double> bc = c.approx - b.approx; + return vec2<double>::dot(bc, ac) >= 0; +} + +template<> CDTVert<double>::CDTVert(const vec2<double> &pt) +{ + this->co.exact = pt; + this->co.approx = pt; + this->co.abs_approx = pt; /* Not used, so does't matter. */ + this->input_ids = nullptr; + this->symedge = nullptr; + this->index = -1; + this->merge_to_index = -1; + this->visit_index = 0; } -template<typename T> CDTVert<T>::CDTVert(const vec2<T> &pt) +#ifdef WITH_GMP +template<> CDTVert<mpq_class>::CDTVert(const vec2<mpq_class> &pt) { - this->co = pt; + this->co.exact = pt; + this->co.approx = double2(pt.x.get_d(), pt.y.get_d()); + this->co.abs_approx = double2(fabs(this->co.approx.x), fabs(this->co.approx.y)); this->input_ids = nullptr; this->symedge = nullptr; this->index = -1; this->merge_to_index = -1; this->visit_index = 0; } +#endif template<typename T> CDTVert<T> *CDTArrangement<T>::add_vert(const vec2<T> &pt) { @@ -805,8 +1031,8 @@ CDTEdge<T> *CDTArrangement<T>::connect_separate_parts(SymEdge<T> *se1, SymEdge<T template<typename T> CDTEdge<T> *CDTArrangement<T>::split_edge(SymEdge<T> *se, T lambda) { /* Split e at lambda. */ - const vec2<T> *a = &se->vert->co; - const vec2<T> *b = &se->next->vert->co; + const vec2<T> *a = &se->vert->co.exact; + const vec2<T> *b = &se->next->vert->co.exact; SymEdge<T> *sesym = sym(se); SymEdge<T> *sesymprev = prev(sesym); SymEdge<T> *sesymprevsym = sym(sesymprev); @@ -918,8 +1144,8 @@ template<typename T> class SiteInfo { */ template<typename T> bool site_lexicographic_sort(const SiteInfo<T> &a, const SiteInfo<T> &b) { - const vec2<T> &co_a = a.v->co; - const vec2<T> &co_b = b.v->co; + const vec2<T> &co_a = a.v->co.exact; + const vec2<T> &co_b = b.v->co.exact; if (co_a[0] < co_b[0]) { return true; } @@ -945,7 +1171,7 @@ template<typename T> void find_site_merges(Array<SiteInfo<T>> &sites) int n = sites.size(); for (int i = 0; i < n - 1; ++i) { int j = i + 1; - while (j < n && sites[j].v->co == sites[i].v->co) { + while (j < n && sites[j].v->co.exact == sites[i].v->co.exact) { sites[j].v->merge_to_index = sites[i].orig_index; ++j; } @@ -957,19 +1183,19 @@ template<typename T> void find_site_merges(Array<SiteInfo<T>> &sites) template<typename T> inline bool vert_left_of_symedge(CDTVert<T> *v, SymEdge<T> *se) { - return orient2d(v->co, se->vert->co, se->next->vert->co) > 0; + return filtered_orient2d(v->co, se->vert->co, se->next->vert->co) > 0; } template<typename T> inline bool vert_right_of_symedge(CDTVert<T> *v, SymEdge<T> *se) { - return orient2d(v->co, se->next->vert->co, se->vert->co) > 0; + return filtered_orient2d(v->co, se->next->vert->co, se->vert->co) > 0; } /* Is se above basel? */ template<typename T> inline bool dc_tri_valid(SymEdge<T> *se, SymEdge<T> *basel, SymEdge<T> *basel_sym) { - return orient2d(se->next->vert->co, basel_sym->vert->co, basel->vert->co) > 0; + return filtered_orient2d(se->next->vert->co, basel_sym->vert->co, basel->vert->co) > 0; } /** @@ -1013,7 +1239,7 @@ void dc_tri(CDTArrangement<T> *cdt, } CDTVert<T> *v3 = sites[start + 2].v; CDTEdge<T> *eb = cdt->add_vert_to_symedge_edge(v3, &ea->symedges[1]); - int orient = orient2d(v1->co, v2->co, v3->co); + int orient = filtered_orient2d(v1->co, v2->co, v3->co); if (orient > 0) { cdt->add_diagonal(&eb->symedges[0], &ea->symedges[0]); *r_le = &ea->symedges[0]; @@ -1101,10 +1327,10 @@ void dc_tri(CDTArrangement<T> *cdt, std::cout << "found valid lcand\n"; std::cout << " lcand" << lcand << "\n"; } - while (incircle(basel_sym->vert->co, - basel->vert->co, - lcand->next->vert->co, - lcand->rot->next->vert->co) > 0.0) { + while (filtered_incircle(basel_sym->vert->co, + basel->vert->co, + lcand->next->vert->co, + lcand->rot->next->vert->co) > 0.0) { if (dbg_level > 1) { std::cout << "incircle says to remove lcand\n"; std::cout << " lcand" << lcand << "\n"; @@ -1120,10 +1346,10 @@ void dc_tri(CDTArrangement<T> *cdt, std::cout << "found valid rcand\n"; std::cout << " rcand" << rcand << "\n"; } - while (incircle(basel_sym->vert->co, - basel->vert->co, - rcand->next->vert->co, - sym(rcand)->next->next->vert->co) > 0.0) { + while (filtered_incircle(basel_sym->vert->co, + basel->vert->co, + rcand->next->vert->co, + sym(rcand)->next->next->vert->co) > 0.0) { if (dbg_level > 0) { std::cout << "incircle says to remove rcand\n"; std::cout << " rcand" << rcand << "\n"; @@ -1147,10 +1373,10 @@ void dc_tri(CDTArrangement<T> *cdt, } /* The next cross edge to be connected is to either `lcand->next->vert` or `rcand->next->vert`; * if both are valid, choose the appropriate one using the #incircle test. */ - if (!valid_lcand || - (valid_rcand && - incircle(lcand->next->vert->co, lcand->vert->co, rcand->vert->co, rcand->next->vert->co) > - 0)) { + if (!valid_lcand || (valid_rcand && filtered_incircle(lcand->next->vert->co, + lcand->vert->co, + rcand->vert->co, + rcand->next->vert->co) > 0)) { if (dbg_level > 0) { std::cout << "connecting rcand\n"; std::cout << " se1=basel_sym" << basel_sym << "\n"; @@ -1265,7 +1491,7 @@ template<typename T> static void re_delaunay_triangulate(CDTArrangement<T> *cdt, SymEdge<T> *cse = first; for (SymEdge<T> *ss = first->next; ss != se; ss = ss->next) { CDTVert<T> *v = ss->vert; - if (incircle(a->co, b->co, c->co, v->co) > 0) { + if (filtered_incircle(a->co, b->co, c->co, v->co) > 0) { c = v; cse = ss; } @@ -1290,7 +1516,7 @@ template<typename T> static void re_delaunay_triangulate(CDTArrangement<T> *cdt, template<typename T> inline int tri_orient(const SymEdge<T> *t) { - return orient2d(t->vert->co, t->next->vert->co, t->next->next->vert->co); + return filtered_orient2d(t->vert->co, t->next->vert->co, t->next->next->vert->co); } /** @@ -1419,7 +1645,7 @@ void fill_crossdata_for_through_vert(CDTVert<T> *v, * case. We need to fill in cd's 'out' edge if it was a vert case. */ template<typename T> -void fill_crossdata_for_intersect(const vec2<T> &curco, +void fill_crossdata_for_intersect(const FatCo<T> &curco, const CDTVert<T> *v2, SymEdge<T> *t, CrossData<T> *cd, @@ -1434,7 +1660,7 @@ void fill_crossdata_for_intersect(const vec2<T> &curco, BLI_assert(se_vcva->vert == vc && se_vcva->next->vert == va); BLI_assert(se_vcvb->vert == vc && se_vcvb->next->vert == vb); UNUSED_VARS_NDEBUG(vc); - auto isect = vec2<T>::isect_seg_seg(va->co, vb->co, curco, v2->co); + auto isect = vec2<T>::isect_seg_seg(va->co.exact, vb->co.exact, curco.exact, v2->co.exact); T &lambda = isect.lambda; switch (isect.kind) { case vec2<T>::isect_result::LINE_LINE_CROSS: { @@ -1443,7 +1669,7 @@ void fill_crossdata_for_intersect(const vec2<T> &curco, #else if (true) { #endif - T len_ab = vec2<T>::distance(va->co, vb->co); + double len_ab = vec2<double>::distance(va->co.approx, vb->co.approx); if (lambda * len_ab <= epsilon) { fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next); } @@ -1497,7 +1723,8 @@ void fill_crossdata_for_intersect(const vec2<T> &curco, break; } case vec2<T>::isect_result::LINE_LINE_COLINEAR: { - if (vec2<T>::distance_squared(va->co, v2->co) <= vec2<T>::distance_squared(vb->co, v2->co)) { + if (vec2<double>::distance_squared(va->co.approx, v2->co.approx) <= + vec2<double>::distance_squared(vb->co.approx, v2->co.approx)) { fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next); } else { @@ -1537,14 +1764,14 @@ bool get_next_crossing_from_vert(CDT_state<T> *cdt_state, } CDTVert<T> *va = t->next->vert; CDTVert<T> *vb = t->next->next->vert; - int orient1 = orient2d(t->vert->co, va->co, v2->co); + int orient1 = filtered_orient2d(t->vert->co, va->co, v2->co); if (orient1 == 0 && in_line<T>(vcur->co, va->co, v2->co)) { fill_crossdata_for_through_vert(va, t, cd, cd_next); ok = true; break; } if (t->face != cdt_state->cdt.outer_face) { - int orient2 = orient2d(vcur->co, vb->co, v2->co); + int orient2 = filtered_orient2d(vcur->co, vb->co, v2->co); /* Don't handle orient2 == 0 case here: next rotation will get it. */ if (orient1 > 0 && orient2 < 0) { /* Segment intersection. */ @@ -1574,15 +1801,16 @@ void get_next_crossing_from_edge(CrossData<T> *cd, { CDTVert<T> *va = cd->in->vert; CDTVert<T> *vb = cd->in->next->vert; - vec2<T> curco = vec2<T>::interpolate(va->co, vb->co, cd->lambda); + vec2<T> curco = vec2<T>::interpolate(va->co.exact, vb->co.exact, cd->lambda); + FatCo<T> fat_curco(curco); SymEdge<T> *se_ac = sym(cd->in)->next; CDTVert<T> *vc = se_ac->next->vert; - int orient = orient2d(curco, v2->co, vc->co); + int orient = filtered_orient2d(fat_curco, v2->co, vc->co); if (orient < 0) { - fill_crossdata_for_intersect<T>(curco, v2, se_ac->next, cd, cd_next, epsilon); + fill_crossdata_for_intersect<T>(fat_curco, v2, se_ac->next, cd, cd_next, epsilon); } else if (orient > 0.0) { - fill_crossdata_for_intersect(curco, v2, se_ac, cd, cd_next, epsilon); + fill_crossdata_for_intersect(fat_curco, v2, se_ac, cd, cd_next, epsilon); } else { *cd_next = CrossData<T>{0.0, vc, se_ac->next, nullptr}; @@ -1682,7 +1910,7 @@ void add_edge_constraint( ok = get_next_crossing_from_vert(cdt_state, cd, cd_next, v2); } else { - get_next_crossing_from_edge(cd, cd_next, v2, cdt_state->epsilon); + get_next_crossing_from_edge<T>(cd, cd_next, v2, cdt_state->epsilon); ok = true; } constexpr int unreasonably_large_crossings = 100000; @@ -2057,7 +2285,7 @@ template<typename T> void remove_non_constraint_edges(CDT_state<T> *cdt_state) * For sorting edges by decreasing length (squared). */ template<typename T> struct EdgeToSort { - T len_squared = T(0); + double len_squared = 0.0; CDTEdge<T> *e{nullptr}; EdgeToSort() = default; @@ -2098,9 +2326,9 @@ template<typename T> void remove_non_constraint_edges_leave_valid_bmesh(CDT_stat if (!is_deleted_edge(e) && !is_constrained_edge(e)) { dissolvable_edges.append(EdgeToSort<T>()); dissolvable_edges[i].e = e; - const vec2<T> &co1 = e->symedges[0].vert->co; - const vec2<T> &co2 = e->symedges[1].vert->co; - dissolvable_edges[i].len_squared = vec2<T>::distance_squared(co1, co2); + const vec2<double> &co1 = e->symedges[0].vert->co.approx; + const vec2<double> &co2 = e->symedges[1].vert->co.approx; + dissolvable_edges[i].len_squared = vec2<double>::distance_squared(co1, co2); i++; } } @@ -2268,7 +2496,7 @@ CDT_result<T> get_cdt_output(CDT_state<T> *cdt_state, for (int i = 0; i < verts_size; ++i) { CDTVert<T> *v = cdt->verts[i]; if (v->merge_to_index == -1) { - result.vert[i_out] = v->co; + result.vert[i_out] = v->co.exact; if (i < cdt_state->input_vert_tot) { result.vert_orig[i_out].append(i); } diff --git a/source/blender/blenlib/intern/math_vec.cc b/source/blender/blenlib/intern/math_vec.cc index 54926f84eb9..84fa6c69d17 100644 --- a/source/blender/blenlib/intern/math_vec.cc +++ b/source/blender/blenlib/intern/math_vec.cc @@ -70,14 +70,13 @@ double2::isect_result double2::isect_seg_seg(const double2 &v1, double div = (v2[0] - v1[0]) * (v4[1] - v3[1]) - (v2[1] - v1[1]) * (v4[0] - v3[0]); if (div == 0.0) { ans.lambda = 0.0; - ans.mu = 0.0; ans.kind = double2::isect_result::LINE_LINE_COLINEAR; } else { ans.lambda = ((v1[1] - v3[1]) * (v4[0] - v3[0]) - (v1[0] - v3[0]) * (v4[1] - v3[1])) / div; - ans.mu = ((v1[1] - v3[1]) * (v2[0] - v1[0]) - (v1[0] - v3[0]) * (v2[1] - v1[1])) / div; - if (ans.lambda >= 0.0 && ans.lambda <= 1.0 && ans.mu >= 0.0 && ans.mu <= 1.0) { - if (ans.lambda == 0.0 || ans.lambda == 1.0 || ans.mu == 0.0 || ans.mu == 1.0) { + double mu = ((v1[1] - v3[1]) * (v2[0] - v1[0]) - (v1[0] - v3[0]) * (v2[1] - v1[1])) / div; + if (ans.lambda >= 0.0 && ans.lambda <= 1.0 && mu >= 0.0 && mu <= 1.0) { + if (ans.lambda == 0.0 || ans.lambda == 1.0 || mu == 0.0 || mu == 1.0) { ans.kind = double2::isect_result::LINE_LINE_EXACT; } else { @@ -101,14 +100,15 @@ mpq2::isect_result mpq2::isect_seg_seg(const mpq2 &v1, mpq_class div = (v2[0] - v1[0]) * (v4[1] - v3[1]) - (v2[1] - v1[1]) * (v4[0] - v3[0]); if (div == 0.0) { ans.lambda = 0.0; - ans.mu = 0.0; ans.kind = mpq2::isect_result::LINE_LINE_COLINEAR; } else { ans.lambda = ((v1[1] - v3[1]) * (v4[0] - v3[0]) - (v1[0] - v3[0]) * (v4[1] - v3[1])) / div; - ans.mu = ((v1[1] - v3[1]) * (v2[0] - v1[0]) - (v1[0] - v3[0]) * (v2[1] - v1[1])) / div; - if (ans.lambda >= 0 && ans.lambda <= 1 && ans.mu >= 0 && ans.mu <= 1) { - if (ans.lambda == 0 || ans.lambda == 1 || ans.mu == 0 || ans.mu == 1) { + /* Avoid dividing mu by div: it is expensive in multiprecision. */ + mpq_class mudiv = ((v1[1] - v3[1]) * (v2[0] - v1[0]) - (v1[0] - v3[0]) * (v2[1] - v1[1])); + if (ans.lambda >= 0 && ans.lambda <= 1 && + ((div > 0 && mudiv >= 0 && mudiv <= div) || (div < 0 && mudiv <= 0 && mudiv >= div))) { + if (ans.lambda == 0 || ans.lambda == 1 || mudiv == 0 || mudiv == div) { ans.kind = mpq2::isect_result::LINE_LINE_EXACT; } else { diff --git a/source/blender/blenlib/tests/BLI_delaunay_2d_test.cc b/source/blender/blenlib/tests/BLI_delaunay_2d_test.cc index 338d6f96bef..487afc095f9 100644 --- a/source/blender/blenlib/tests/BLI_delaunay_2d_test.cc +++ b/source/blender/blenlib/tests/BLI_delaunay_2d_test.cc @@ -21,6 +21,7 @@ extern "C" { #include "BLI_array.hh" #include "BLI_double2.hh" +#include "BLI_math_boolean.hh" #include "BLI_math_mpq.hh" #include "BLI_mpq2.hh" #include "BLI_vector.hh" @@ -1696,7 +1697,7 @@ void rand_delaunay_test(int test_kind, in.vert[ic][1] = T((param * sin(angle3))); /* Put the coordinates in ccw order. */ in.face[i].append(ia); - int orient = vec2<T>::orient2d(in.vert[ia], in.vert[ib], in.vert[ic]); + int orient = orient2d(in.vert[ia], in.vert[ib], in.vert[ic]); if (orient >= 0) { in.face[i].append(ib); in.face[i].append(ic); |