2011년 9월 22일 목요일

Fixing delaunay triangulation

During delaunay triangulation unexpacted line crossing is performed.
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

In this case it is necessary that additional algorithm is to avoid a crossing between lines.

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
....

    // 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,




댓글 없음:

댓글 쓰기