Given points are like this.
And Its result is
Its code is
void delaunay_list_create_2d_network(delaunay_list_t *list, point_list_t *points, bitmap_t *roi)
{
dlink_t *x, *y, *z;
point_t *o, *a, *b, *c;
point_t *v1, *v2, *v3, *org;
delaunay_net_t *neta, *netb, *netc, *net;
//dlist_t *triangles;
triangle_list_t *triangles;
triangle_t *tri;
//dlist_t *lines;
line_list_t *lines;
line_t *line, *l;
real_t da, db, dc, s;
real_t val, vmin;
real_t radius, area;
assert(list);
assert(points);
// allocate delaunay network list
//list = delaunay_list_new();
for (x = points->tail->next; x != points->head; x = x->next) {
net = delaunay_net_new((point_t *)x->object);
delaunay_list_insert(net, list);
}
// starting point
neta = delaunay_list_get_net(0, list);
a = delaunay_net_get_coordinate(neta);
// find minimum distance between points
for (netb = NULL, x = list->tail->next->next; x != list->head; x = x->next) {
net = (delaunay_net_t *)x->object;
o = delaunay_net_get_coordinate(net);
val = distance_between_points(a, o);
if (val < REAL_EPSILON) continue;
if (netb == NULL || val < vmin) { netb = net; vmin = val; }
}
b = delaunay_net_get_coordinate(netb);
// find minimum radius of circumscribed circle for triangle
for (netc = NULL, x = list->tail->next->next; x != list->head; x = x->next) {
net = (delaunay_net_t *)x->object;
if (net == netb) continue;
o = delaunay_net_get_coordinate(net);
// Heron's formula to find the area formed by three points
da = distance_between_points(a, b);
db = distance_between_points(b, o);
dc = distance_between_points(o, a);
s = (da + db + dc) / 2;
area = sqrt(s * (s - da) * (s - db) * (s - dc));
if (area < REAL_EPSILON) continue;
// calculate radius of circumscribed circle for triangle
radius = (da * db * dc) / (4 * area);
if (netc == NULL || radius < vmin) { netc = net; vmin = radius; }
}
c = delaunay_net_get_coordinate(netc);
v1 = point_new();
v2 = point_new();
v3 = point_new();
org = point_new();
// re-ordering
point_subtract(v1, b, a);
point_subtract(v2, c, a);
point_xproduct(v3, v1, v2);
if (v3->z > 0) { // counterclockwise following a to b to c ?
net = netb, netb = netc, netc = net;
o = b, b = c, c = o;
}
// connect delaunay network
delaunay_net_connect(neta, netb);
delaunay_net_connect(netb, netc);
delaunay_net_connect(netc, neta);
// Synthesize a triangle and insert it.
//triangles = dlist_new();
triangles = triangle_list_new();
tri = triangle_new_and_assign(a, b, c);
triangle_list_insert(tri, triangles);
//triangle_dump(tri);
// Synthesize three line and insert it.
lines = line_list_new();
line = line_new_and_assign(a, b);
line_list_insert(line, lines);
line = line_new_and_assign(b, c);
line_list_insert(line, lines);
line = line_new_and_assign(c, a);
line_list_insert(line, lines);
// main routine
for (x = lines->tail->next; x != lines->head; x = x->next) {
line = (line_t *)x->object;
//line_dump(line);
for (y = list->tail->next; y != list->head; y = y->next) {
netc = (delaunay_net_t *)y->object;
c = delaunay_net_get_coordinate(netc);
if (point_cmp(c, line->a) == 0 || point_cmp(c, line->b) == 0)
continue;
// counter-clockwise rotation of the triangle formed by a, b, and c
point_subtract(v1, line->b, line->a);
point_subtract(v2, c, line->a);
point_xproduct(v3, v1, v2);
if (v3->z < REAL_EPSILON) continue;
// Does the triangle exist in the pool of delaunay triangles before ?
tri = triangle_new_and_assign(line->a, c, line->b);
//if (is_triangle_in_list(tri, triangles)) {
if (triangle_list_query(tri, triangles) >= 0) {
triangle_destroy(tri);
tri = NULL;
continue;
}
// Does the unique triangle form the delaunay triangle ?
radius = triangle_circumcircle(org, tri);
for (z = list->tail->next; z != list->head; z = z->next) {
net = (delaunay_net_t *)z->object;
o = delaunay_net_get_coordinate(net);
if (point_cmp(o, line->a) == 0 ||
point_cmp(o, line->b) == 0 ||
point_cmp(o, c) == 0)
continue;
if (distance_between_points(o, org) < radius) {
triangle_destroy(tri);
tri = NULL;
break;
}
}
if (tri) {
/* put the trinangle into triangle's list,
* establish connection between delaunay networks, and
* allocate and put two lines into line's list */
triangle_list_insert(tri, triangles);
l = line_new_and_assign(line->a, c);
if (line_list_query(l, lines) < 0) {
line_list_insert(l, lines);
for (z = list->tail->next; z != list->head; z = z->next) {
net = (delaunay_net_t *)z->object;
o = delaunay_net_get_coordinate(net);
if (point_cmp(o, line->a) == 0) {
delaunay_net_connect(netc, net);
break;
}
}
} else line_destroy(l);
l = line_new_and_assign(c, line->b);
if (line_list_query(l, lines) < 0) {
line_list_insert(l, lines);
for (z = list->tail->next; z != list->head; z = z->next) {
net = (delaunay_net_t *)z->object;
o = delaunay_net_get_coordinate(net);
if (point_cmp(o, line->b) == 0) {
delaunay_net_connect(netc, net);
break;
}
}
} else line_destroy(l);
break;
}
}
}
// de-allocation lines and triangles
/*
while (triangles->count) {
tri = triangle_list_pop(triangles);
//triangle_dump(tri);
triangle_destroy(tri);
}
triangle_list_destroy(triangles);
*/
list->head->object = (void *)triangles;
while (lines->count) {
line = line_list_pop(lines);
//line_dump(line);
line_destroy(line);
}
line_list_destroy(lines);
point_destroy(v1);
point_destroy(v2);
point_destroy(v3);
point_destroy(org);
// return list;
}
Its reason is because of four points on the same circle. Delaunay Triangulation is determined by whether there is a point in circumference circle, but for this four points it is impossible to construct the unique delaunay triangulation.
The triangle with red color has properties like this;
three vertex: (333, 311), (333, 183), (308, 315)
center of circumference circle: (309.94, 247)
radius of circumference circle: 68.027688
The triangle with blue color has properties like this;
three vertex: (308, 315), (333, 311), (308, 179)
center of circumference circle: (309.94, 247)
radius of circumference circle: 68.027688
The subroutine to detect whether two lines has intersecting point or not. It is referenced at http://paulbourke.net/geometry/lineline2d/
int intersect_point_between_lines(point_t *p, line_t *l, line_t *m)
{
real_t denom, ua, ub;
assert(p);
assert(l);
assert(m);
// Pl = l->a + Ua * (l->b - l->a)
// Pm = m->a + Ub * (m->b - m->a)
// for Pl = Pm
denom = (m->b->y - m->a->y) * (l->b->x - l->a->x) - (m->b->x - m->a->x) * (l->b->y - l->a->y);
ua = (m->b->x - m->a->x) * (l->a->y - m->a->y) - (m->b->y - m->a->y) * (l->a->x - m->a->x);
ub = (l->b->x - l->a->x) * (l->a->y - m->a->y) - (l->b->y - l->a->y) * (l->a->x - m->a->x);
if (abs(denom) < REAL_EPSILON) {
if (abs(ua) < REAL_EPSILON && abs(ub) < REAL_EPSILON)
return 0; // COINCIDENT
return -1; // PARALLEL
}
ua = ua / denom;
ub = ub / denom;
if (ua < 0 || ua > 1 || ub < 0 || ub > 1)
return -2; // Not intersecting
p->x = l->a->x + ua * (l->b->x - l->a->x);
p->y = l->a->y + ua * (l->b->y - l->a->y);
if ((abs(ua) < REAL_EPSILON || abs(ua - 1.0) < REAL_EPSILON) ||
(abs(ub) < REAL_EPSILON || abs(ub - 1.0) < REAL_EPSILON))
return 2; // ENDPOINT
return 1; // INTERSECTING
}
And the code of modified delaunay triangulation is here.
void delaunay_list_create_2d_network(delaunay_list_t *list, point_list_t *points, bitmap_t *roi)
{
dlink_t *x, *y, *z;
point_t *o, *a, *b, *c;
point_t *v1, *v2, *v3, *org;
delaunay_net_t *neta, *netb, *netc, *net;
//dlist_t *triangles;
triangle_list_t *triangles;
triangle_t *tri;
//dlist_t *lines;
line_list_t *lines;
line_t *line, *l, *m;
real_t da, db, dc, s;
real_t val, vmin;
real_t radius, area;
assert(list);
assert(points);
// allocate delaunay network list
//list = delaunay_list_new();
for (x = points->tail->next; x != points->head; x = x->next) {
net = delaunay_net_new((point_t *)x->object);
delaunay_list_insert(net, list);
}
// starting point
neta = delaunay_list_get_net(0, list);
a = delaunay_net_get_coordinate(neta);
// find minimum distance between points
for (netb = NULL, x = list->tail->next->next; x != list->head; x = x->next) {
net = (delaunay_net_t *)x->object;
o = delaunay_net_get_coordinate(net);
val = distance_between_points(a, o);
if (val < REAL_EPSILON) continue;
if (netb == NULL || val < vmin) { netb = net; vmin = val; }
}
b = delaunay_net_get_coordinate(netb);
// find minimum radius of circumscribed circle for triangle
for (netc = NULL, x = list->tail->next->next; x != list->head; x = x->next) {
net = (delaunay_net_t *)x->object;
if (net == netb) continue;
o = delaunay_net_get_coordinate(net);
// Heron's formula to find the area formed by three points
da = distance_between_points(a, b);
db = distance_between_points(b, o);
dc = distance_between_points(o, a);
s = (da + db + dc) / 2;
area = sqrt(s * (s - da) * (s - db) * (s - dc));
if (area < REAL_EPSILON) continue;
// calculate radius of circumscribed circle for triangle
radius = (da * db * dc) / (4 * area);
if (netc == NULL || radius < vmin) { netc = net; vmin = radius; }
}
c = delaunay_net_get_coordinate(netc);
v1 = point_new();
v2 = point_new();
v3 = point_new();
org = point_new();
// re-ordering
point_subtract(v1, b, a);
point_subtract(v2, c, a);
point_xproduct(v3, v1, v2);
if (v3->z > 0) { // counterclockwise following a to b to c ?
net = netb, netb = netc, netc = net;
o = b, b = c, c = o;
}
// connect delaunay network
delaunay_net_connect(neta, netb);
delaunay_net_connect(netb, netc);
delaunay_net_connect(netc, neta);
// Synthesize a triangle and insert it.
//triangles = dlist_new();
triangles = triangle_list_new();
tri = triangle_new_and_assign(a, b, c);
triangle_list_insert(tri, triangles);
//triangle_dump(tri);
// Synthesize three line and insert it.
lines = line_list_new();
line = line_new_and_assign(a, b);
line_list_insert(line, lines);
line = line_new_and_assign(b, c);
line_list_insert(line, lines);
line = line_new_and_assign(c, a);
line_list_insert(line, lines);
// main routine
for (x = lines->tail->next; x != lines->head; x = x->next) {
line = (line_t *)x->object;
//line_dump(line);
for (y = list->tail->next; y != list->head; y = y->next) {
netc = (delaunay_net_t *)y->object;
c = delaunay_net_get_coordinate(netc);
if (point_cmp(c, line->a) == 0 || point_cmp(c, line->b) == 0)
continue;
// counter-clockwise rotation of the triangle formed by a, b, and c
point_subtract(v1, line->b, line->a);
point_subtract(v2, c, line->a);
point_xproduct(v3, v1, v2);
if (v3->z < REAL_EPSILON) continue;
// Does any line crossing with line in the line list ?
l = line_new_and_assign(line->a, c);
m = line_new_and_assign(c, line->b);
for (z = lines->tail->next; z != lines->head; z = z->next) {
// Is there an intersecting ?
if ((intersect_point_between_lines(v3, l, (line_t *)z->object) == 1) ||
(intersect_point_between_lines(v3, m, (line_t *)z->object) == 1)) {
line_destroy(l); l = NULL;
line_destroy(m); m = NULL;
break;
}
}
if (l == NULL && m == NULL) continue;
// Does the triangle exist in the pool of delaunay triangles before ?
tri = triangle_new_and_assign(line->a, c, line->b);
//if (is_triangle_in_list(tri, triangles)) {
if (triangle_list_query(tri, triangles) >= 0) {
triangle_destroy(tri);
tri = NULL;
continue;
}
// Does the unique triangle form the delaunay triangle ?
radius = triangle_circumcircle(org, tri);
for (z = list->tail->next; z != list->head; z = z->next) {
net = (delaunay_net_t *)z->object;
o = delaunay_net_get_coordinate(net);
if (point_cmp(o, line->a) == 0 ||
point_cmp(o, line->b) == 0 ||
point_cmp(o, c) == 0)
continue;
if (distance_between_points(o, org) < radius) {
triangle_destroy(tri);
tri = NULL;
break;
}
}
if (tri) {
/* put the trinangle into triangle's list,
* establish connection between delaunay networks, and
* allocate and put two lines into line's list */
triangle_list_insert(tri, triangles);
//triangle_dump(tri);
//point_dump(org);
//printf("radius: %lf\n", radius);
//l = line_new_and_assign(line->a, c);
if (line_list_query(l, lines) < 0) {
line_list_insert(l, lines);
//line_dump(l);
for (z = list->tail->next; z != list->head; z = z->next) {
net = (delaunay_net_t *)z->object;
o = delaunay_net_get_coordinate(net);
if (point_cmp(o, line->a) == 0) {
delaunay_net_connect(netc, net);
break;
}
}
} else line_destroy(l);
//m = line_new_and_assign(c, line->b);
if (line_list_query(m, lines) < 0) {
line_list_insert(m, lines);
//line_dump(m);
for (z = list->tail->next; z != list->head; z = z->next) {
net = (delaunay_net_t *)z->object;
o = delaunay_net_get_coordinate(net);
if (point_cmp(o, line->b) == 0) {
delaunay_net_connect(netc, net);
break;
}
}
} else line_destroy(m);
//keyhit();
break;
}
}
}
// de-allocation lines and triangles
/*
while (triangles->count) {
tri = triangle_list_pop(triangles);
//triangle_dump(tri);
triangle_destroy(tri);
}
triangle_list_destroy(triangles);
*/
list->head->object = (void *)triangles;
while (lines->count) {
line = line_list_pop(lines);
//line_dump(line);
line_destroy(line);
}
line_list_destroy(lines);
point_destroy(v1);
point_destroy(v2);
point_destroy(v3);
point_destroy(org);
// return list;
}
At result the line crossing is disappeared.
Unfortunately, two delaunay triangles are deleted in the figure. In the figure dot is the mean center of delaunay triangle. As you see, you can find two triangles without mean center. It is remaining side effect of four points on the same circle.
So, it needs to change comparison between variables expressed in floating points.
Change
....
// Does the unique triangle form the delaunay triangle ?
radius = triangle_circumcircle(org, tri);
for (z = list->tail->next; z != list->head; z = z->next) {
net = (delaunay_net_t *)z->object;
o = delaunay_net_get_coordinate(net);
if (point_cmp(o, line->a) == 0 ||
point_cmp(o, line->b) == 0 ||
point_cmp(o, c) == 0)
continue;
if (distance_between_points(o, org) < radius) {
triangle_destroy(tri);
tri = NULL;
break;
}
}
....
to
....
....
to
....
// Does the unique triangle form the delaunay triangle ?
radius = triangle_circumcircle(org, tri);
for (z = list->tail->next; z != list->head; z = z->next) {
net = (delaunay_net_t *)z->object;
o = delaunay_net_get_coordinate(net);
if (point_cmp(o, line->a) == 0 ||
point_cmp(o, line->b) == 0 ||
point_cmp(o, c) == 0)
continue;
val = distance_between_points(o, org);
if ((val < radius) && (round(100*val) < round(100*radius)) {
triangle_destroy(tri);
tri = NULL;
break;
}
}
....
Then,
....
Then,
댓글 없음:
댓글 쓰기