Fix two bugs in delaunay blenlib function.

Bugs were: (1) needed an epsilon test in CCW test in order to
handle new costraint edge that intersects an existing point
but only within epsilon; (2) the "valid bmesh" output mode
sometimes left a face that included outside frame point.
This commit is contained in:
Howard Trickey
2019-09-07 20:28:03 +05:30
parent 0b2d1badec
commit b380a98887
2 changed files with 56 additions and 17 deletions

View File

@@ -115,8 +115,9 @@ static void validate_face_centroid(SymEdge *se);
static void validate_cdt(CDT_state *cdt, bool check_all_tris); static void validate_cdt(CDT_state *cdt, bool check_all_tris);
#endif #endif
/** return 1 if a,b,c forms CCW angle, -1 if a CW angle, 0 if straight */ /** return 1 if a,b,c forms CCW angle, -1 if a CW angle, 0 if straight.
static int CCW_test(const double a[2], const double b[2], const double c[2]) * For straight test, allow b to be withing eps of line. */
static int CCW_test(const double a[2], const double b[2], const double c[2], const double eps)
{ {
double det; double det;
double ab; double ab;
@@ -124,14 +125,14 @@ static int CCW_test(const double a[2], const double b[2], const double c[2])
/* This is twice the signed area of triangle abc. */ /* This is twice the signed area of triangle abc. */
det = (b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1]); det = (b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1]);
ab = len_v2v2_db(a, b); ab = len_v2v2_db(a, b);
if (ab < DBL_EPSILON) { if (ab <= eps) {
return 0; return 0;
} }
det /= ab; det /= ab;
if (det > DBL_EPSILON) { if (det > eps) {
return 1; return 1;
} }
else if (det < -DBL_EPSILON) { else if (det < -eps) {
return -1; return -1;
} }
return 0; return 0;
@@ -662,7 +663,7 @@ static bool locate_point_final(const double p[2],
} }
else { else {
dist_inside[i] = len_close_p; dist_inside[i] = len_close_p;
dist_inside[i] = CCW_test(a, b, p) >= 0 ? len_close_p : -len_close_p; dist_inside[i] = CCW_test(a, b, p, epsilon) >= 0 ? len_close_p : -len_close_p;
} }
i++; i++;
se = se->next; se = se->next;
@@ -813,7 +814,8 @@ static LocateResult locate_point(CDT_state *cdt, const double p[2])
a = cur_se->vert->co; a = cur_se->vert->co;
b = cur_se->next->vert->co; b = cur_se->next->vert->co;
c = cur_se->next->next->vert->co; c = cur_se->next->next->vert->co;
if (CCW_test(a, b, p) >= 0 && CCW_test(b, c, p) >= 0 && CCW_test(c, a, p) >= 0) { if (CCW_test(a, b, p, epsilon) >= 0 && CCW_test(b, c, p, epsilon) >= 0 &&
CCW_test(c, a, p, epsilon) >= 0) {
#ifdef DEBUG_CDT #ifdef DEBUG_CDT
if (dbglevel > 1) { if (dbglevel > 1) {
fprintf(stderr, "p in current triangle\n"); fprintf(stderr, "p in current triangle\n");
@@ -837,7 +839,7 @@ static LocateResult locate_point(CDT_state *cdt, const double p[2])
} }
#endif #endif
next_se_sym = sym(next_se); next_se_sym = sym(next_se);
if (CCW_test(a, b, p) <= 0 && next_se->face != cdt->outer_face) { if (CCW_test(a, b, p, epsilon) <= 0 && next_se->face != cdt->outer_face) {
#ifdef DEBUG_CDT #ifdef DEBUG_CDT
if (dbglevel > 1) { if (dbglevel > 1) {
fprintf(stderr, "CCW_test(a, b, p) <= 0\n"); fprintf(stderr, "CCW_test(a, b, p) <= 0\n");
@@ -1431,6 +1433,7 @@ static void add_edge_constraint(
int ccw1, ccw2, isect; int ccw1, ccw2, isect;
int i, search_count; int i, search_count;
double lambda; double lambda;
const double epsilon = cdt->epsilon;
bool done, state_through_vert; bool done, state_through_vert;
LinkNodePair edge_list = {NULL, NULL}; LinkNodePair edge_list = {NULL, NULL};
typedef struct CrossData { typedef struct CrossData {
@@ -1541,8 +1544,8 @@ static void add_edge_constraint(
do { do {
va = t->next->vert; va = t->next->vert;
vb = t->next->next->vert; vb = t->next->next->vert;
ccw1 = CCW_test(t->vert->co, va->co, v2->co); ccw1 = CCW_test(t->vert->co, va->co, v2->co, epsilon);
ccw2 = CCW_test(t->vert->co, vb->co, v2->co); ccw2 = CCW_test(t->vert->co, vb->co, v2->co, epsilon);
#ifdef DEBUG_CDT #ifdef DEBUG_CDT
if (dbg_level > 1) { if (dbg_level > 1) {
fprintf(stderr, "non-final through vert case\n"); fprintf(stderr, "non-final through vert case\n");
@@ -1591,7 +1594,7 @@ static void add_edge_constraint(
} }
#endif #endif
} while (t != tstart); } while (t != tstart);
BLI_assert(tout != NULL); /* TODO: something sensivle for "this can't happen" */ BLI_assert(tout != NULL); /* TODO: something sensible for "this can't happen" */
crossings[BLI_array_len(crossings) - 1].out = tout; crossings[BLI_array_len(crossings) - 1].out = tout;
} }
} }
@@ -1634,7 +1637,7 @@ static void add_edge_constraint(
/* 'tout' is 'symedge' from 'vb' to third vertex, 'vc'. */ /* 'tout' is 'symedge' from 'vb' to third vertex, 'vc'. */
BLI_assert(tout->vert == va); BLI_assert(tout->vert == va);
vc = tout->next->vert; vc = tout->next->vert;
ccw1 = CCW_test(v1->co, v2->co, vc->co); ccw1 = CCW_test(v1->co, v2->co, vc->co, epsilon);
#ifdef DEBUG_CDT #ifdef DEBUG_CDT
if (dbg_level > 1) { if (dbg_level > 1) {
fprintf(stderr, "now searching with third vertex "); fprintf(stderr, "now searching with third vertex ");
@@ -1911,7 +1914,7 @@ static void remove_non_constraint_edges(CDT_state *cdt, const bool valid_bmesh)
dissolve = !is_deleted_edge(e) && !is_constrained_edge(e); dissolve = !is_deleted_edge(e) && !is_constrained_edge(e);
if (dissolve) { if (dissolve) {
se = &e->symedges[0]; se = &e->symedges[0];
if (valid_bmesh) { if (valid_bmesh && !edge_touches_frame(e)) {
fleft = se->face; fleft = se->face;
fright = sym(se)->face; fright = sym(se)->face;
if (fleft != cdt->outer_face && fright != cdt->outer_face && if (fleft != cdt->outer_face && fright != cdt->outer_face &&
@@ -2227,7 +2230,7 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty
CDTEdge *face_edge; CDTEdge *face_edge;
SymEdge *face_symedge; SymEdge *face_symedge;
#ifdef DEBUG_CDT #ifdef DEBUG_CDT
int dbg_level = 1; int dbg_level = 0;
#endif #endif
if ((nv > 0 && input->vert_coords == NULL) || (ne > 0 && input->edges == NULL) || if ((nv > 0 && input->vert_coords == NULL) || (ne > 0 && input->edges == NULL) ||
@@ -2285,6 +2288,13 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty
} }
add_edge_constraint(cdt, verts[v1], verts[v2], i, NULL); add_edge_constraint(cdt, verts[v1], verts[v2], i, NULL);
} }
#ifdef DEBUG_CDT
if (dbg_level > 2) {
cdt_draw(cdt, "after edge constraints");
dump_cdt(cdt, "after edge constraints");
validate_cdt(cdt, true);
}
#endif
cdt->face_edge_offset = ne; cdt->face_edge_offset = ne;
for (f = 0; f < nf; f++) { for (f = 0; f < nf; f++) {
int flen = input->faces_len_table[f]; int flen = input->faces_len_table[f];
@@ -2317,6 +2327,11 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty
F2(cdt_e->symedges[1].vert->co)); F2(cdt_e->symedges[1].vert->co));
} }
} }
if (dbg_level > 2) {
cdt_draw(cdt, "after a face edge");
dump_cdt(cdt, "after a face edge");
validate_cdt(cdt, true);
}
#endif #endif
if (i == 0) { if (i == 0) {
face_edge = (CDTEdge *)edge_list->link; face_edge = (CDTEdge *)edge_list->link;

View File

@@ -15,8 +15,6 @@ extern "C" {
#include <fstream> #include <fstream>
#include <sstream> #include <sstream>
#define DLNY_EPSILON 1e-8
static void fill_input_verts(CDT_input *r_input, float (*vcos)[2], int nverts) static void fill_input_verts(CDT_input *r_input, float (*vcos)[2], int nverts)
{ {
r_input->verts_len = nverts; r_input->verts_len = nverts;
@@ -27,7 +25,7 @@ static void fill_input_verts(CDT_input *r_input, float (*vcos)[2], int nverts)
r_input->faces = NULL; r_input->faces = NULL;
r_input->faces_start_table = NULL; r_input->faces_start_table = NULL;
r_input->faces_len_table = NULL; r_input->faces_len_table = NULL;
r_input->epsilon = 1e-6f; r_input->epsilon = 1e-5f;
} }
static void add_input_edges(CDT_input *r_input, int (*edges)[2], int nedges) static void add_input_edges(CDT_input *r_input, int (*edges)[2], int nedges)
@@ -643,6 +641,32 @@ TEST(delaunay, TwoSquaresOverlap)
BLI_delaunay_2d_cdt_free(out); BLI_delaunay_2d_cdt_free(out);
} }
TEST(delaunay, TriCutoff)
{
CDT_input in;
CDT_result *out;
float p[][2] = {
{-3.53009f, 1.29403f},
{-4.11844f, -1.08375f},
{1.56893f, 1.29403f},
{0.621034f, 0.897734f},
{0.549125f, 1.29403f},
};
int f[] = {0, 2, 1};
int fstart[] = {0};
int flen[] = {3};
int e[][2] = {{3, 4}};
fill_input_verts(&in, p, 5);
add_input_faces(&in, f, fstart, flen, 1);
add_input_edges(&in, e, 1);
out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS_VALID_BMESH);
EXPECT_EQ(out->verts_len, 5);
EXPECT_EQ(out->edges_len, 6);
EXPECT_EQ(out->faces_len, 2);
BLI_delaunay_2d_cdt_free(out);
}
enum { enum {
RANDOM_PTS, RANDOM_PTS,
RANDOM_SEGS, RANDOM_SEGS,