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.
This commit is contained in:
parent
38fe962d95
commit
df8cc5662b
Notes:
blender-bot
2023-02-14 06:00:51 +01:00
Referenced by issue #82914, Manipulators don't render properly
|
|
@ -131,7 +131,6 @@ struct double2 {
|
|||
LINE_LINE_CROSS = 2,
|
||||
} kind;
|
||||
double lambda;
|
||||
double mu;
|
||||
};
|
||||
|
||||
static isect_result isect_seg_seg(const double2 &v1,
|
||||
|
|
|
|||
|
|
@ -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,
|
||||
|
|
|
|||
|
|
@ -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,25 +589,168 @@ 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<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> 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;
|
||||
}
|
||||
return vec2<T>::dot(bc, ac) >= 0;
|
||||
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<double> ab = b.approx - a.approx;
|
||||
vec2<double> ac = c.approx - a.approx;
|
||||
if (vec2<double>::dot(ab, ac) < 0) {
|
||||
return false;
|
||||
}
|
||||
vec2<double> bc = c.approx - b.approx;
|
||||
return vec2<double>::dot(bc, ac) >= 0;
|
||||
}
|
||||
|
||||
template<typename T> CDTVert<T>::CDTVert(const vec2<T> &pt)
|
||||
template<> CDTVert<double>::CDTVert(const vec2<double> &pt)
|
||||
{
|
||||
this->co = 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;
|
||||
|
|
@ -546,6 +758,20 @@ template<typename T> CDTVert<T>::CDTVert(const vec2<T> &pt)
|
|||
this->visit_index = 0;
|
||||
}
|
||||
|
||||
#ifdef WITH_GMP
|
||||
template<> CDTVert<mpq_class>::CDTVert(const vec2<mpq_class> &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)
|
||||
{
|
||||
CDTVert<T> *v = new CDTVert<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);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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 {
|
||||
|
|
|
|||
|
|
@ -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);
|
||||
|
|
|
|||
Loading…
Reference in New Issue