Greatly improve speed of Delaunay when have a lot of holes.

Using part of a patch from Erik Abrahamsson, this replaces the
use of linked lists for original id tracking by Sets.
I had thought that the lists were unlikely to grow to more than
a few elements, but when the mesh has a lot of holes (whose
original ids go *outside* the hole, and therefore, most of the
mesh), this assumption can be very wrong.
On a Text regression test, the time went from 11.67s to 0.16s
with this fix. I also tested to make sure that Boolean didn't
slow down with this, and found it actually had a very slight speedup.

Using Sets exposed a dependency on the ordering of the items
in the id lists, luckily caught by a mesh intersect regression test,
so fixed that.
This commit is contained in:
Howard Trickey 2021-07-18 10:48:57 -04:00
parent e82c5c6607
commit 3c8d261557
3 changed files with 206 additions and 77 deletions

View File

@ -29,6 +29,7 @@
#include "BLI_math_boolean.hh"
#include "BLI_math_mpq.hh"
#include "BLI_mpq2.hh"
#include "BLI_set.hh"
#include "BLI_vector.hh"
#include "BLI_delaunay_2d.h"
@ -77,8 +78,8 @@ template<> double math_to_double<double>(const double v)
* Define a templated 2D arrangement of vertices, edges, and faces.
* The #SymEdge data structure is the basis for a structure that allows
* easy traversal to neighboring (by topology) geometric elements.
* Each of #CDTVert, #CDTEdge, and #CDTFace have an input_id linked list,
* whose nodes contain integers that keep track of which input verts, edges,
* Each of #CDTVert, #CDTEdge, and #CDTFace have an input_id set,
* which contain integers that keep track of which input verts, edges,
* and faces, respectively, that the element was derived from.
*
* While this could be cleaned up some, it is usable by other routines in Blender
@ -195,8 +196,8 @@ template<typename T> struct CDTVert {
FatCo<T> co;
/** Some edge attached to it. */
SymEdge<T> *symedge{nullptr};
/** List of corresponding vertex input ids. */
LinkNode *input_ids{nullptr};
/** Set of corresponding vertex input ids. */
blender::Set<int> input_ids;
/** Index into array that #CDTArrangement keeps. */
int index{-1};
/** Index of a CDTVert that this has merged to. -1 if no merge. */
@ -209,8 +210,8 @@ template<typename T> struct CDTVert {
};
template<typename Arith_t> struct CDTEdge {
/** List of input edge ids that this is part of. */
LinkNode *input_ids{nullptr};
/** Set of input edge ids that this is part of. */
blender::Set<int> input_ids;
/** The directed edges for this edge. */
SymEdge<Arith_t> symedges[2]{SymEdge<Arith_t>(), SymEdge<Arith_t>()};
@ -220,8 +221,8 @@ template<typename Arith_t> struct CDTEdge {
template<typename Arith_t> struct CDTFace {
/** A symedge in face; only used during output, so only valid then. */
SymEdge<Arith_t> *symedge{nullptr};
/** List of input face ids that this is part of. */
LinkNode *input_ids{nullptr};
/** Set of input face ids that this is part of. */
blender::Set<int> input_ids;
/** Used by algorithms operating on CDT structures. */
int visit_index{0};
/** Marks this face no longer used. */
@ -342,19 +343,19 @@ template<typename T> CDTArrangement<T>::~CDTArrangement()
{
for (int i : this->verts.index_range()) {
CDTVert<T> *v = this->verts[i];
BLI_linklist_free(v->input_ids, nullptr);
v->input_ids.clear();
delete v;
this->verts[i] = nullptr;
}
for (int i : this->edges.index_range()) {
CDTEdge<T> *e = this->edges[i];
BLI_linklist_free(e->input_ids, nullptr);
e->input_ids.clear();
delete e;
this->edges[i] = nullptr;
}
for (int i : this->faces.index_range()) {
CDTFace<T> *f = this->faces[i];
BLI_linklist_free(f->input_ids, nullptr);
f->input_ids.clear();
delete f;
this->faces[i] = nullptr;
}
@ -495,6 +496,7 @@ template<typename T> void cdt_draw(const std::string &label, const CDTArrangemen
constexpr int vert_radius = 3;
constexpr bool draw_vert_labels = true;
constexpr bool draw_edge_labels = false;
constexpr bool draw_face_labels = false;
if (cdt.verts.size() == 0) {
return;
@ -559,7 +561,7 @@ template<typename T> void cdt_draw(const std::string &label, const CDTArrangemen
const CDTVert<T> *v = e->symedges[1].vert;
const vec2<double> &uco = u->co.approx;
const vec2<double> &vco = v->co.approx;
int strokew = e->input_ids == nullptr ? thin_line : thick_line;
int strokew = e->input_ids.size() == 0 ? thin_line : thick_line;
f << R"(<line fill="none" stroke="black" stroke-width=")" << strokew << "\" x1=\""
<< SX(uco[0]) << "\" y1=\"" << SY(uco[1]) << "\" x2=\"" << SX(vco[0]) << "\" y2=\""
<< SY(vco[1]) << "\">\n";
@ -586,6 +588,54 @@ template<typename T> void cdt_draw(const std::string &label, const CDTArrangemen
++i;
}
if (draw_face_labels) {
for (const CDTFace<T> *face : cdt.faces) {
if (!face->deleted) {
/* Since may not have prepared output yet, need a slow find of a SymEdge for this face. */
const SymEdge<T> *se_face_start = nullptr;
for (const CDTEdge<T> *e : cdt.edges) {
if (e->symedges[0].face == face) {
se_face_start = &e->symedges[0];
break;
}
if (e->symedges[1].face == face) {
se_face_start = &e->symedges[1];
}
}
if (se_face_start != nullptr) {
/* Find center of face. */
int face_nverts = 0;
vec2<double> cen(0.0, 0.0);
if (face == cdt.outer_face) {
cen.x = minx;
cen.y = miny;
}
else {
const SymEdge<T> *se = se_face_start;
do {
if (se->face == face) {
cen = cen + se->vert->co.approx;
face_nverts++;
}
} while ((se = se->next) != se_face_start);
if (face_nverts > 0) {
cen = cen / double(face_nverts);
}
}
f << "<text x=\"" << SX(cen[0]) << "\" y=\"" << SY(cen[1]) << "\""
<< " font-size=\"small\">[";
f << trunc_ptr(face);
if (face->input_ids.size() > 0) {
for (int id : face->input_ids) {
f << " " << id;
}
}
f << "]</text>\n";
}
}
}
}
append = true;
# undef SX
# undef SY
@ -754,7 +804,6 @@ template<> CDTVert<double>::CDTVert(const vec2<double> &pt)
this->co.exact = pt;
this->co.approx = pt;
this->co.abs_approx = pt; /* Not used, so doesn't matter. */
this->input_ids = nullptr;
this->symedge = nullptr;
this->index = -1;
this->merge_to_index = -1;
@ -767,7 +816,6 @@ 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;
@ -833,26 +881,10 @@ CDT_state<T>::CDT_state(int num_input_verts, int num_input_edges, int num_input_
this->visit_count = 0;
}
static bool id_in_list(const LinkNode *id_list, int id)
{
const LinkNode *ln;
for (ln = id_list; ln != nullptr; ln = ln->next) {
if (POINTER_AS_INT(ln->link) == id) {
return true;
}
}
return false;
}
/* Is any id in (range_start, range_start+1, ... , range_end) in id_list? */
static bool id_range_in_list(const LinkNode *id_list, int range_start, int range_end)
static bool id_range_in_list(const blender::Set<int> &id_list, int range_start, int range_end)
{
const LinkNode *ln;
int id;
for (ln = id_list; ln != nullptr; ln = ln->next) {
id = POINTER_AS_INT(ln->link);
for (int id : id_list) {
if (id >= range_start && id <= range_end) {
return true;
}
@ -860,19 +892,15 @@ static bool id_range_in_list(const LinkNode *id_list, int range_start, int range
return false;
}
static void add_to_input_ids(LinkNode **dst, int input_id)
static void add_to_input_ids(blender::Set<int> &dst, int input_id)
{
if (!id_in_list(*dst, input_id)) {
BLI_linklist_prepend(dst, POINTER_FROM_INT(input_id));
}
dst.add(input_id);
}
static void add_list_to_input_ids(LinkNode **dst, const LinkNode *src)
static void add_list_to_input_ids(blender::Set<int> &dst, const blender::Set<int> &src)
{
const LinkNode *ln;
for (ln = src; ln != nullptr; ln = ln->next) {
add_to_input_ids(dst, POINTER_AS_INT(ln->link));
for (int value : src) {
dst.add(value);
}
}
@ -883,7 +911,7 @@ template<typename T> inline bool is_border_edge(const CDTEdge<T> *e, const CDT_s
template<typename T> inline bool is_constrained_edge(const CDTEdge<T> *e)
{
return e->input_ids != nullptr;
return e->input_ids.size() > 0;
}
template<typename T> inline bool is_deleted_edge(const CDTEdge<T> *e)
@ -979,7 +1007,7 @@ template<typename T> CDTEdge<T> *CDTArrangement<T>::add_diagonal(SymEdge<T> *s1,
for (SymEdge<T> *se = s2; se != sdiag; se = se->next) {
se->face = fnew;
}
add_list_to_input_ids(&fnew->input_ids, fold->input_ids);
add_list_to_input_ids(fnew->input_ids, fold->input_ids);
return ediag;
}
@ -1058,7 +1086,7 @@ template<typename T> CDTEdge<T> *CDTArrangement<T>::split_edge(SymEdge<T> *se, T
if (newsesym->vert->symedge == sesym) {
newsesym->vert->symedge = newsesym;
}
add_list_to_input_ids(&e->input_ids, se->edge->input_ids);
add_list_to_input_ids(e->input_ids, se->edge->input_ids);
return e;
}
@ -1880,7 +1908,7 @@ void add_edge_constraint(
SymEdge<T> *t = find_symedge_between_verts(v1, v2);
if (t != nullptr) {
/* Segment already there. */
add_to_input_ids(&t->edge->input_ids, input_id);
add_to_input_ids(t->edge->input_ids, input_id);
if (r_edges != nullptr) {
BLI_linklist_append(&edge_list, t->edge);
*r_edges = edge_list.list;
@ -2041,7 +2069,7 @@ void add_edge_constraint(
BLI_assert(cd_prev->lambda == 0.0);
BLI_assert(cd_prev->out->next->vert == cd->vert);
edge = cd_prev->out->edge;
add_to_input_ids(&edge->input_ids, input_id);
add_to_input_ids(edge->input_ids, input_id);
if (r_edges != nullptr) {
BLI_linklist_append(&edge_list, edge);
}
@ -2054,7 +2082,7 @@ void add_edge_constraint(
else {
edge = cdt_state->cdt.add_diagonal(tstart, t);
}
add_to_input_ids(&edge->input_ids, input_id);
add_to_input_ids(edge->input_ids, input_id);
if (r_edges != nullptr) {
BLI_linklist_append(&edge_list, edge);
}
@ -2132,7 +2160,7 @@ void add_face_ids(
continue;
}
face->visit_index = visit;
add_to_input_ids(&face->input_ids, face_id);
add_to_input_ids(face->input_ids, face_id);
SymEdge<T> *se_start = se;
for (se = se->next; se != se_start; se = se->next) {
if (!id_range_in_list(se->edge->input_ids, fedge_start, fedge_end)) {
@ -2349,7 +2377,7 @@ template<typename T> void remove_non_constraint_edges_leave_valid_bmesh(CDT_stat
CDTFace<T> *fleft = se->face;
CDTFace<T> *fright = sym(se)->face;
if (fleft != cdt->outer_face && fright != cdt->outer_face &&
(fleft->input_ids != nullptr || fright->input_ids != nullptr)) {
(fleft->input_ids.size() > 0 || fright->input_ids.size() > 0)) {
/* Is there another #SymEdge with same left and right faces?
* Or is there a vertex not part of e touching the same left and right faces? */
for (SymEdge<T> *se2 = se->next; dissolve && se2 != se; se2 = se2->next) {
@ -2632,7 +2660,7 @@ CDT_result<T> get_cdt_output(CDT_state<T> *cdt_state,
CDTVert<T> *v = cdt->verts[i];
if (v->merge_to_index != -1) {
if (i < cdt_state->input_vert_tot) {
add_to_input_ids(&cdt->verts[v->merge_to_index]->input_ids, i);
add_to_input_ids(cdt->verts[v->merge_to_index]->input_ids, i);
}
vert_to_output_map[i] = vert_to_output_map[v->merge_to_index];
}
@ -2648,8 +2676,8 @@ CDT_result<T> get_cdt_output(CDT_state<T> *cdt_state,
if (i < cdt_state->input_vert_tot) {
result.vert_orig[i_out].append(i);
}
for (LinkNode *ln = v->input_ids; ln; ln = ln->next) {
result.vert_orig[i_out].append(POINTER_AS_INT(ln->link));
for (int vert : v->input_ids) {
result.vert_orig[i_out].append(vert);
}
++i_out;
}
@ -2667,8 +2695,8 @@ CDT_result<T> get_cdt_output(CDT_state<T> *cdt_state,
int vo1 = vert_to_output_map[e->symedges[0].vert->index];
int vo2 = vert_to_output_map[e->symedges[1].vert->index];
result.edge[e_out] = std::pair<int, int>(vo1, vo2);
for (LinkNode *ln = e->input_ids; ln; ln = ln->next) {
result.edge_orig[e_out].append(POINTER_AS_INT(ln->link));
for (int edge : e->input_ids) {
result.edge_orig[e_out].append(edge);
}
++e_out;
}
@ -2690,8 +2718,8 @@ CDT_result<T> get_cdt_output(CDT_state<T> *cdt_state,
result.face[f_out].append(vert_to_output_map[se->vert->index]);
se = se->next;
} while (se != se_start);
for (LinkNode *ln = f->input_ids; ln; ln = ln->next) {
result.face_orig[f_out].append(POINTER_AS_INT(ln->link));
for (int face : f->input_ids) {
result.face_orig[f_out].append(face);
}
++f_out;
}

View File

@ -1875,6 +1875,16 @@ static void do_cdt(CDT_data &cd)
}
}
/* Find an original edge index that goes with the edge that the CDT output edge
* that goes between verts i0 and i1 (in the CDT output vert indexing scheme).
* There may be more than one: if so, prefer one that was originally a face edge.
* The input to CDT for a triangle with some intersecting segments from other triangles
* will have both edges and a face constraint for the main triangle (this is redundant
* but allows us to discover which face edge goes with which output edges).
* If there is any face edge, return one of those as the original.
* If there is no face edge but there is another edge in the input problem, then that
* edge must have come from intersection with another triangle, so set *r_is_intersect
* to true in that case. */
static int get_cdt_edge_orig(
int i0, int i1, const CDT_data &cd, const IMesh &in_tm, bool *r_is_intersect)
{
@ -1900,35 +1910,50 @@ static int get_cdt_edge_orig(
}
/* Pick an arbitrary orig, but not one equal to NO_INDEX, if we can help it. */
/* TODO: if edge has origs from more than on part of the nary input,
/* TODO: if edge has origs from more than one part of the nary input,
* then want to set *r_is_intersect to true. */
int face_eorig = NO_INDEX;
bool have_non_face_eorig = false;
for (int orig_index : cd.cdt_out.edge_orig[e]) {
/* orig_index encodes the triangle and pos within the triangle of the input edge. */
if (orig_index >= foff) {
int in_face_index = (orig_index / foff) - 1;
int pos = orig_index % foff;
/* We need to retrieve the edge orig field from the Face used to populate the
* in_face_index'th face of the CDT, at the pos'th position of the face. */
int in_tm_face_index = cd.input_face[in_face_index];
BLI_assert(in_tm_face_index < in_tm.face_size());
const Face *facep = in_tm.face(in_tm_face_index);
BLI_assert(pos < facep->size());
bool is_rev = cd.is_reversed[in_face_index];
int eorig = is_rev ? facep->edge_orig[2 - pos] : facep->edge_orig[pos];
if (eorig != NO_INDEX) {
return eorig;
if (face_eorig == NO_INDEX) {
int in_face_index = (orig_index / foff) - 1;
int pos = orig_index % foff;
/* We need to retrieve the edge orig field from the Face used to populate the
* in_face_index'th face of the CDT, at the pos'th position of the face. */
int in_tm_face_index = cd.input_face[in_face_index];
BLI_assert(in_tm_face_index < in_tm.face_size());
const Face *facep = in_tm.face(in_tm_face_index);
BLI_assert(pos < facep->size());
bool is_rev = cd.is_reversed[in_face_index];
int eorig = is_rev ? facep->edge_orig[2 - pos] : facep->edge_orig[pos];
if (eorig != NO_INDEX) {
face_eorig = eorig;
}
}
}
else {
/* This edge came from an edge input to the CDT problem,
* so it is an intersect edge. */
*r_is_intersect = true;
/* TODO: maybe there is an orig index:
* This happens if an input edge was formed by an input face having
* an edge that is co-planar with the cluster, while the face as a whole is not. */
return NO_INDEX;
if (!have_non_face_eorig) {
have_non_face_eorig = true;
}
if (face_eorig != NO_INDEX && have_non_face_eorig) {
/* Only need at most one orig for each type. */
break;
}
}
}
if (face_eorig != NO_INDEX) {
return face_eorig;
}
else if (have_non_face_eorig) {
/* This must have been an input to the CDT problem that was an intersection edge. */
/* TODO: maybe there is an orig index:
* This happens if an input edge was formed by an input face having
* an edge that is co-planar with the cluster, while the face as a whole is not. */
*r_is_intersect = true;
return NO_INDEX;
}
return NO_INDEX;
}

View File

@ -241,7 +241,7 @@ template<typename T> std::ostream &operator<<(std::ostream &os, const CDT_result
for (int i : r.edge.index_range()) {
os << "e" << i << " = (" << r.edge[i].first << ", " << r.edge[i].second << ")\n";
os << " orig: ";
for (int j : r.edge_orig[i].size()) {
for (int j : r.edge_orig[i].index_range()) {
os << r.edge_orig[i][j] << " ";
}
os << "\n";
@ -797,6 +797,55 @@ template<typename T> void crosssegs_test()
}
}
template<typename T> void cutacrosstri_test()
{
/* Right triangle with horizontal segment exactly crossing in the middle. */
const char *spec = R"(5 1 1
0.0 0.0
1.0 0.0
0.0 1.0
0.0 0.5
0.5 0.5
3 4
0 1 2
)";
CDT_input<T> in = fill_input_from_string<T>(spec);
CDT_result<T> out = delaunay_2d_calc(in, CDT_FULL);
EXPECT_EQ(out.vert.size(), 5);
EXPECT_EQ(out.edge.size(), 7);
EXPECT_EQ(out.face.size(), 3);
int v0_out = get_orig_index(out.vert_orig, 0);
int v1_out = get_orig_index(out.vert_orig, 1);
int v2_out = get_orig_index(out.vert_orig, 2);
int v3_out = get_orig_index(out.vert_orig, 3);
int v4_out = get_orig_index(out.vert_orig, 4);
EXPECT_TRUE(v0_out != -1 && v1_out != -1 && v2_out != -1 && v3_out != -1 && v4_out != -1);
if (out.face.size() == 3) {
int e0_out = get_orig_index(out.edge_orig, 0);
EXPECT_NE(e0_out, -1);
int fe0_out = get_output_edge_index(out, v0_out, v1_out);
EXPECT_NE(fe0_out, -1);
int fe1a_out = get_output_edge_index(out, v1_out, v4_out);
EXPECT_NE(fe1a_out, -1);
int fe1b_out = get_output_edge_index(out, v4_out, v2_out);
EXPECT_NE(fe1b_out, -1);
if (fe1a_out != 0 && fe1b_out != 0) {
EXPECT_EQ(e0_out, get_orig_index(out.edge_orig, 0));
EXPECT_TRUE(out.edge_orig[fe1a_out].size() == 1 && out.edge_orig[fe1a_out][0] == 11);
EXPECT_TRUE(out.edge_orig[fe1b_out].size() == 1 && out.edge_orig[fe1b_out][0] == 11);
}
int e_diag = get_output_edge_index(out, v0_out, v4_out);
EXPECT_NE(e_diag, -1);
if (e_diag != -1) {
EXPECT_EQ(out.edge_orig[e_diag].size(), 0);
}
}
if (DO_DRAW) {
graph_draw<T>("CutAcrossTri", out.vert, out.edge, out.face);
}
}
template<typename T> void diamondcross_test()
{
/* Diamond with constraint edge from top to bottom. Some dup verts. */
@ -1470,6 +1519,11 @@ TEST(delaunay_d, CrossSegs)
crosssegs_test<double>();
}
TEST(delaunay_d, CutAcrossTri)
{
cutacrosstri_test<double>();
}
TEST(delaunay_d, DiamondCross)
{
diamondcross_test<double>();
@ -1610,6 +1664,11 @@ TEST(delaunay_m, CrossSegs)
crosssegs_test<mpq_class>();
}
TEST(delaunay_m, CutAcrossTri)
{
cutacrosstri_test<mpq_class>();
}
TEST(delaunay_m, DiamondCross)
{
diamondcross_test<mpq_class>();
@ -1876,6 +1935,23 @@ TEST(delaunay_d, TextB10_10_10)
text_test<double>(10, 10, 10, CDT_INSIDE_WITH_HOLES);
}
# ifdef WITH_GMP
TEST(delaunay_m, TextB10)
{
text_test<mpq_class>(10, 1, 1, CDT_INSIDE_WITH_HOLES);
}
TEST(delaunay_m, TextB200)
{
text_test<mpq_class>(200, 1, 1, CDT_INSIDE_WITH_HOLES);
}
TEST(delaunay_m, TextB10_10_10)
{
text_test<mpq_class>(10, 10, 10, CDT_INSIDE_WITH_HOLES);
}
# endif
#endif
#if DO_RANDOM_TESTS