delaunay.c (20598B)
1 /* 2 ** delaunay.c : compute 2D delaunay triangulation in the plane. 3 ** Copyright (C) 2005 Wael El Oraiby <wael.eloraiby@gmail.com> 4 ** 5 ** 6 ** This program is free software: you can redistribute it and/or modify 7 ** it under the terms of the GNU Affero General Public License as 8 ** published by the Free Software Foundation, either version 3 of the 9 ** License, or (at your option) any later version. 10 ** 11 ** This program is distributed in the hope that it will be useful, 12 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 13 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 ** GNU Affero General Public License for more details. 15 ** 16 ** You should have received a copy of the GNU Affero General Public License 17 ** along with this program. If not, see <http://www.gnu.org/licenses/>. 18 */ 19 20 #include <stdio.h> 21 #include <stdlib.h> 22 #include <math.h> 23 #include <string.h> 24 #include <assert.h> 25 26 #include "delaunay.h" 27 28 #define ON_RIGHT 1 29 #define ON_SEG 0 30 #define ON_LEFT -1 31 32 #define OUTSIDE -1 33 #define ON_CIRCLE 0 34 #define INSIDE 1 35 36 struct point2d_s; 37 struct face_s; 38 struct halfedge_s; 39 struct delaunay_s; 40 41 42 43 #define REAL_ZERO 0.0l 44 #define REAL_ONE 1.0l 45 #define REAL_TWO 2.0l 46 #define REAL_FOUR 4.0l 47 48 49 typedef struct point2d_s point2d_t; 50 typedef struct face_s face_t; 51 typedef struct halfedge_s halfedge_t; 52 typedef struct delaunay_s delaunay_t; 53 typedef struct working_set_s working_set_t; 54 55 typedef long double lreal; 56 typedef lreal mat3_t[3][3]; 57 58 struct point2d_s { 59 real x, y; /* point coordinates */ 60 halfedge_t* he; /* point halfedge */ 61 unsigned int idx; /* point index in input buffer */ 62 }; 63 64 struct face_s { 65 halfedge_t* he; /* a pointing half edge */ 66 unsigned int num_verts; /* number of vertices on this face */ 67 }; 68 69 struct halfedge_s { 70 point2d_t* vertex; /* vertex */ 71 halfedge_t* pair; /* pair */ 72 halfedge_t* next; /* next */ 73 halfedge_t* prev; /* next^-1 */ 74 face_t* face; /* halfedge face */ 75 }; 76 77 struct delaunay_s { 78 halfedge_t* rightmost_he; /* right most halfedge */ 79 halfedge_t* leftmost_he; /* left most halfedge */ 80 point2d_t* points; /* pointer to points */ 81 face_t* faces; /* faces of delaunay */ 82 unsigned int num_faces; /* face count */ 83 unsigned int start_point; /* start point index */ 84 unsigned int end_point; /* end point index */ 85 }; 86 87 struct working_set_s { 88 halfedge_t* edges; /* all the edges (allocated in one shot) */ 89 face_t* faces; /* all the faces (allocated in one shot) */ 90 91 unsigned int max_edge; /* maximum edge count: 2 * 3 * n where n is point count */ 92 unsigned int max_face; /* maximum face count: 2 * n where n is point count */ 93 94 unsigned int num_edges; /* number of allocated edges */ 95 unsigned int num_faces; /* number of allocated faces */ 96 97 halfedge_t* free_edge; /* pointer to the first free edge */ 98 face_t* free_face; /* pointer to the first free face */ 99 }; 100 101 /* 102 * 3x3 matrix determinant 103 */ 104 static lreal det3(mat3_t m) 105 { 106 lreal res = m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) 107 - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) 108 + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]); 109 110 return res; 111 } 112 113 /* 114 * allocate a halfedge 115 */ 116 static halfedge_t* halfedge_alloc() 117 { 118 halfedge_t* d; 119 120 d = (halfedge_t*)malloc(sizeof(halfedge_t)); 121 assert( NULL != d ); 122 memset(d, 0, sizeof(halfedge_t)); 123 124 return d; 125 } 126 127 /* 128 * free a halfedge 129 */ 130 static void halfedge_free( halfedge_t* d ) 131 { 132 assert( d != NULL ); 133 memset(d, 0, sizeof(halfedge_t)); 134 free(d); 135 } 136 137 /* 138 * free all delaunay halfedges 139 */ 140 void del_free_halfedges( delaunay_t *del ) 141 { 142 unsigned int i; 143 halfedge_t *d, *sig; 144 145 /* if there is nothing to do */ 146 if( del->points == NULL ) 147 return; 148 149 for( i = 0; i <= (del->end_point - del->start_point); i++ ) 150 { 151 /* free all the halfedges around the point */ 152 d = del->points[i].he; 153 if( d != NULL ) 154 { 155 do { 156 sig = d->next; 157 halfedge_free( d ); 158 d = sig; 159 } while( d != del->points[i].he ); 160 del->points[i].he = NULL; 161 } 162 } 163 } 164 165 /* 166 * compare 2 points when sorting 167 */ 168 static int cmp_points( const void *_pt0, const void *_pt1 ) 169 { 170 point2d_t *pt0, *pt1; 171 172 pt0 = (point2d_t*)(_pt0); 173 pt1 = (point2d_t*)(_pt1); 174 175 if( pt0->x < pt1->x ) 176 return -1; 177 else if( pt0->x > pt1->x ) 178 return 1; 179 else if( pt0->y < pt1->y ) 180 return -1; 181 else if( pt0->y > pt1->y ) 182 return 1; 183 printf("same: pt0 %f %f pt1 %f %f\n", 184 pt0->x, pt0->y, pt1->x, pt1->y); 185 assert(0 && "2 or more points share the same exact coordinate"); 186 return 0; /* Should not be given! */ 187 } 188 189 /* 190 * classify a point relative to a segment 191 */ 192 static int classify_point_seg( point2d_t *s, point2d_t *e, point2d_t *pt ) 193 { 194 lreal se_x, se_y, spt_x, spt_y; 195 lreal res; 196 197 se_x = e->x - s->x; 198 se_y = e->y - s->y; 199 200 spt_x = pt->x - s->x; 201 spt_y = pt->y - s->y; 202 203 res = (( se_x * spt_y ) - ( se_y * spt_x )); 204 if( res < REAL_ZERO ) 205 return ON_RIGHT; 206 else if( res > REAL_ZERO ) 207 return ON_LEFT; 208 209 return ON_SEG; 210 } 211 212 /* 213 * classify a point relative to a halfedge, -1 is left, 0 is on, 1 is right 214 */ 215 static int del_classify_point( halfedge_t *d, point2d_t *pt ) 216 { 217 point2d_t *s, *e; 218 219 s = d->vertex; 220 e = d->pair->vertex; 221 222 return classify_point_seg(s, e, pt); 223 } 224 225 /* 226 * test if a point is inside a circle given by 3 points, 1 if inside, 0 if outside 227 */ 228 static int in_circle( point2d_t *pt0, point2d_t *pt1, point2d_t *pt2, point2d_t *p ) 229 { 230 // reduce the computational complexity by substracting the last row of the matrix 231 // ref: https://www.cs.cmu.edu/~quake/robust.html 232 lreal p0p_x, p0p_y, p1p_x, p1p_y, p2p_x, p2p_y, p0p, p1p, p2p, res; 233 mat3_t m; 234 235 p0p_x = pt0->x - p->x; 236 p0p_y = pt0->y - p->y; 237 238 p1p_x = pt1->x - p->x; 239 p1p_y = pt1->y - p->y; 240 241 p2p_x = pt2->x - p->x; 242 p2p_y = pt2->y - p->y; 243 244 p0p = p0p_x * p0p_x + p0p_y * p0p_y; 245 p1p = p1p_x * p1p_x + p1p_y * p1p_y; 246 p2p = p2p_x * p2p_x + p2p_y * p2p_y; 247 248 m[0][0] = p0p_x; 249 m[0][1] = p0p_y; 250 m[0][2] = p0p; 251 252 m[1][0] = p1p_x; 253 m[1][1] = p1p_y; 254 m[1][2] = p1p; 255 256 m[2][0] = p2p_x; 257 m[2][1] = p2p_y; 258 m[2][2] = p2p; 259 260 res = -det3(m); 261 262 if( res < REAL_ZERO ) 263 return INSIDE; 264 else if( res > REAL_ZERO ) 265 return OUTSIDE; 266 267 return ON_CIRCLE; 268 } 269 270 /* 271 * initialize delaunay segment 272 */ 273 static int del_init_seg( delaunay_t *del, int start ) 274 { 275 halfedge_t *d0, *d1; 276 point2d_t *pt0, *pt1; 277 278 /* init delaunay */ 279 del->start_point = start; 280 del->end_point = start + 1; 281 282 /* setup pt0 and pt1 */ 283 pt0 = &(del->points[start]); 284 pt1 = &(del->points[start + 1]); 285 286 /* allocate the halfedges and setup them */ 287 d0 = halfedge_alloc(); 288 d1 = halfedge_alloc(); 289 290 d0->vertex = pt0; 291 d1->vertex = pt1; 292 293 d0->next = d0->prev = d0; 294 d1->next = d1->prev = d1; 295 296 d0->pair = d1; 297 d1->pair = d0; 298 299 pt0->he = d0; 300 pt1->he = d1; 301 302 del->rightmost_he = d1; 303 del->leftmost_he = d0; 304 305 306 return 0; 307 } 308 309 /* 310 * initialize delaunay triangle 311 */ 312 static int del_init_tri( delaunay_t *del, int start ) 313 { 314 halfedge_t *d0, *d1, *d2, *d3, *d4, *d5; 315 point2d_t *pt0, *pt1, *pt2; 316 317 /* initiate delaunay */ 318 del->start_point = start; 319 del->end_point = start + 2; 320 321 /* setup the points */ 322 pt0 = &(del->points[start]); 323 pt1 = &(del->points[start + 1]); 324 pt2 = &(del->points[start + 2]); 325 326 /* allocate the 6 halfedges */ 327 d0 = halfedge_alloc(); 328 d1 = halfedge_alloc(); 329 d2 = halfedge_alloc(); 330 d3 = halfedge_alloc(); 331 d4 = halfedge_alloc(); 332 d5 = halfedge_alloc(); 333 334 if( classify_point_seg(pt0, pt2, pt1) == ON_LEFT ) /* first case */ 335 { 336 /* set halfedges points */ 337 d0->vertex = pt0; 338 d1->vertex = pt2; 339 d2->vertex = pt1; 340 341 d3->vertex = pt2; 342 d4->vertex = pt1; 343 d5->vertex = pt0; 344 345 /* set points halfedges */ 346 pt0->he = d0; 347 pt1->he = d2; 348 pt2->he = d1; 349 350 /* next and next -1 setup */ 351 d0->next = d5; 352 d0->prev = d5; 353 354 d1->next = d3; 355 d1->prev = d3; 356 357 d2->next = d4; 358 d2->prev = d4; 359 360 d3->next = d1; 361 d3->prev = d1; 362 363 d4->next = d2; 364 d4->prev = d2; 365 366 d5->next = d0; 367 d5->prev = d0; 368 369 /* set halfedges pair */ 370 d0->pair = d3; 371 d3->pair = d0; 372 373 d1->pair = d4; 374 d4->pair = d1; 375 376 d2->pair = d5; 377 d5->pair = d2; 378 379 del->rightmost_he = d1; 380 del->leftmost_he = d0; 381 382 } else /* 2nd case */ 383 { 384 /* set halfedges points */ 385 d0->vertex = pt0; 386 d1->vertex = pt1; 387 d2->vertex = pt2; 388 389 d3->vertex = pt1; 390 d4->vertex = pt2; 391 d5->vertex = pt0; 392 393 /* set points halfedges */ 394 pt0->he = d0; 395 pt1->he = d1; 396 pt2->he = d2; 397 398 /* next and next -1 setup */ 399 d0->next = d5; 400 d0->prev = d5; 401 402 d1->next = d3; 403 d1->prev = d3; 404 405 d2->next = d4; 406 d2->prev = d4; 407 408 d3->next = d1; 409 d3->prev = d1; 410 411 d4->next = d2; 412 d4->prev = d2; 413 414 d5->next = d0; 415 d5->prev = d0; 416 417 /* set halfedges pair */ 418 d0->pair = d3; 419 d3->pair = d0; 420 421 d1->pair = d4; 422 d4->pair = d1; 423 424 d2->pair = d5; 425 d5->pair = d2; 426 427 del->rightmost_he = d2; 428 del->leftmost_he = d0; 429 } 430 431 return 0; 432 } 433 434 /* 435 * remove an edge given a halfedge 436 */ 437 static void del_remove_edge( halfedge_t *d ) 438 { 439 halfedge_t *next, *prev, *pair, *orig_pair; 440 441 orig_pair = d->pair; 442 443 next = d->next; 444 prev = d->prev; 445 pair = d->pair; 446 447 assert(next != NULL); 448 assert(prev != NULL); 449 450 next->prev = prev; 451 prev->next = next; 452 453 454 /* check to see if we have already removed pair */ 455 if( pair ) 456 pair->pair = NULL; 457 458 /* check to see if the vertex points to this halfedge */ 459 if( d->vertex->he == d ) 460 d->vertex->he = next; 461 462 d->vertex = NULL; 463 d->next = NULL; 464 d->prev = NULL; 465 d->pair = NULL; 466 467 next = orig_pair->next; 468 prev = orig_pair->prev; 469 pair = orig_pair->pair; 470 471 assert(next != NULL); 472 assert(prev != NULL); 473 474 next->prev = prev; 475 prev->next = next; 476 477 478 /* check to see if we have already removed pair */ 479 if( pair ) 480 pair->pair = NULL; 481 482 /* check to see if the vertex points to this halfedge */ 483 if( orig_pair->vertex->he == orig_pair ) 484 orig_pair->vertex->he = next; 485 486 orig_pair->vertex = NULL; 487 orig_pair->next = NULL; 488 orig_pair->prev = NULL; 489 orig_pair->pair = NULL; 490 491 492 /* finally free the halfedges */ 493 halfedge_free(d); 494 halfedge_free(orig_pair); 495 } 496 497 /* 498 * pass through all the halfedges on the left side and validate them 499 */ 500 static halfedge_t* del_valid_left( halfedge_t* b ) 501 { 502 point2d_t *g, *d, *u, *v; 503 halfedge_t *c, *du, *dg; 504 505 g = b->vertex; /* base halfedge point */ 506 dg = b; 507 508 d = b->pair->vertex; /* pair(halfedge) point */ 509 b = b->next; 510 511 u = b->pair->vertex; /* next(pair(halfedge)) point */ 512 du = b->pair; 513 514 v = b->next->pair->vertex; /* pair(next(next(halfedge)) point */ 515 516 if( classify_point_seg(g, d, u) == ON_LEFT ) 517 { 518 /* 3 points aren't colinear */ 519 /* as long as the 4 points belong to the same circle, do the cleaning */ 520 assert( v != u && "1: floating point precision error"); 521 while( v != d && v != g && in_circle(g, d, u, v) == INSIDE ) 522 { 523 c = b->next; 524 du = b->next->pair; 525 del_remove_edge(b); 526 b = c; 527 u = du->vertex; 528 v = b->next->pair->vertex; 529 } 530 531 assert( v != u && "2: floating point precision error"); 532 if( v != d && v != g && in_circle(g, d, u, v) == ON_CIRCLE ) 533 { 534 du = du->prev; 535 del_remove_edge(b); 536 } 537 } else /* treat the case where the 3 points are colinear */ 538 du = dg; 539 540 assert(du->pair); 541 return du; 542 } 543 544 /* 545 * pass through all the halfedges on the right side and validate them 546 */ 547 static halfedge_t* del_valid_right( halfedge_t *b ) 548 { 549 point2d_t *rv, *lv, *u, *v; 550 halfedge_t *c, *dd, *du; 551 552 b = b->pair; 553 rv = b->vertex; 554 dd = b; 555 lv = b->pair->vertex; 556 b = b->prev; 557 u = b->pair->vertex; 558 du = b->pair; 559 560 v = b->prev->pair->vertex; 561 562 if( classify_point_seg(lv, rv, u) == ON_LEFT ) 563 { 564 assert( v != u && "1: floating point precision error"); 565 while( v != lv && v != rv && in_circle(lv, rv, u, v) == INSIDE ) 566 { 567 c = b->prev; 568 du = c->pair; 569 del_remove_edge(b); 570 b = c; 571 u = du->vertex; 572 v = b->prev->pair->vertex; 573 } 574 575 assert( v != u && "1: floating point precision error"); 576 if( v != lv && v != rv && in_circle(lv, rv, u, v) == ON_CIRCLE ) 577 { 578 du = du->next; 579 del_remove_edge(b); 580 } 581 } else 582 du = dd; 583 584 assert(du->pair); 585 return du; 586 } 587 588 589 /* 590 * validate a link 591 */ 592 static halfedge_t* del_valid_link( halfedge_t *b ) 593 { 594 point2d_t *g, *g_p, *d, *d_p; 595 halfedge_t *gd, *dd, *new_gd, *new_dd; 596 int a; 597 598 g = b->vertex; 599 gd = del_valid_left(b); 600 g_p = gd->vertex; 601 602 assert(b->pair); 603 d = b->pair->vertex; 604 dd = del_valid_right(b); 605 d_p = dd->vertex; 606 assert(b->pair); 607 608 if( g != g_p && d != d_p ) { 609 a = in_circle(g, d, g_p, d_p); 610 611 if( a != ON_CIRCLE ) { 612 if( a == INSIDE ) { 613 g_p = g; 614 gd = b; 615 } else { 616 d_p = d; 617 dd = b->pair; 618 } 619 } 620 } 621 622 /* create the 2 halfedges */ 623 new_gd = halfedge_alloc(); 624 new_dd = halfedge_alloc(); 625 626 /* setup new_gd and new_dd */ 627 628 new_gd->vertex = gd->vertex; 629 new_gd->pair = new_dd; 630 new_gd->prev = gd; 631 new_gd->next = gd->next; 632 gd->next->prev = new_gd; 633 gd->next = new_gd; 634 635 new_dd->vertex = dd->vertex; 636 new_dd->pair = new_gd; 637 new_dd->prev = dd->prev; 638 dd->prev->next = new_dd; 639 new_dd->next = dd; 640 dd->prev = new_dd; 641 642 return new_gd; 643 } 644 645 /* 646 * find the lower tangent between the two delaunay, going from left to right (returns the left half edge) 647 */ 648 static halfedge_t* del_get_lower_tangent( delaunay_t *left, delaunay_t *right ) 649 { 650 point2d_t *pl, *pr; 651 halfedge_t *right_d, *left_d, *new_ld, *new_rd; 652 int sl, sr; 653 654 left_d = left->rightmost_he; 655 right_d = right->leftmost_he; 656 657 do { 658 pl = left_d->prev->pair->vertex; 659 pr = right_d->pair->vertex; 660 661 if( (sl = classify_point_seg(left_d->vertex, right_d->vertex, pl)) == ON_RIGHT ) { 662 left_d = left_d->prev->pair; 663 } 664 665 if( (sr = classify_point_seg(left_d->vertex, right_d->vertex, pr)) == ON_RIGHT ) { 666 right_d = right_d->pair->next; 667 } 668 669 } while( sl == ON_RIGHT || sr == ON_RIGHT ); 670 671 /* create the 2 halfedges */ 672 new_ld = halfedge_alloc(); 673 new_rd = halfedge_alloc(); 674 675 /* setup new_gd and new_dd */ 676 new_ld->vertex = left_d->vertex; 677 new_ld->pair = new_rd; 678 new_ld->prev = left_d->prev; 679 left_d->prev->next = new_ld; 680 new_ld->next = left_d; 681 left_d->prev = new_ld; 682 683 new_rd->vertex = right_d->vertex; 684 new_rd->pair = new_ld; 685 new_rd->prev = right_d->prev; 686 right_d->prev->next = new_rd; 687 new_rd->next = right_d; 688 right_d->prev = new_rd; 689 690 return new_ld; 691 } 692 693 /* 694 * link the 2 delaunay together 695 */ 696 static void del_link( delaunay_t *result, delaunay_t *left, delaunay_t *right ) 697 { 698 point2d_t *u, *v, *ml, *mr; 699 halfedge_t *base; 700 701 assert( left->points == right->points ); 702 703 /* save the most right point and the most left point */ 704 ml = left->leftmost_he->vertex; 705 mr = right->rightmost_he->vertex; 706 707 base = del_get_lower_tangent(left, right); 708 709 u = base->next->pair->vertex; 710 v = base->pair->prev->pair->vertex; 711 712 while( del_classify_point(base, u) == ON_LEFT || 713 del_classify_point(base, v) == ON_LEFT ) 714 { 715 base = del_valid_link(base); 716 u = base->next->pair->vertex; 717 v = base->pair->prev->pair->vertex; 718 } 719 720 right->rightmost_he = mr->he; 721 left->leftmost_he = ml->he; 722 723 /* TODO: this part is not needed, and can be optimized */ 724 while( del_classify_point( right->rightmost_he, right->rightmost_he->prev->pair->vertex ) == ON_RIGHT ) 725 right->rightmost_he = right->rightmost_he->prev; 726 727 while( del_classify_point( left->leftmost_he, left->leftmost_he->prev->pair->vertex ) == ON_RIGHT ) 728 left->leftmost_he = left->leftmost_he->prev; 729 730 result->leftmost_he = left->leftmost_he; 731 result->rightmost_he = right->rightmost_he; 732 result->points = left->points; 733 result->start_point = left->start_point; 734 result->end_point = right->end_point; 735 } 736 737 /* 738 * divide and conquer delaunay 739 */ 740 void del_divide_and_conquer( delaunay_t *del, int start, int end ) 741 { 742 delaunay_t left, right; 743 int i, n; 744 745 n = (end - start + 1); 746 747 if( n > 3 ) { 748 i = (n / 2) + (n & 1); 749 left.points = del->points; 750 right.points = del->points; 751 del_divide_and_conquer( &left, start, start + i - 1 ); 752 del_divide_and_conquer( &right, start + i, end ); 753 del_link( del, &left, &right ); 754 } else { 755 if( n == 3 ) { 756 del_init_tri( del, start ); 757 } else { 758 if( n == 2 ) { 759 del_init_seg( del, start ); 760 } 761 } 762 } 763 } 764 765 static void build_halfedge_face( delaunay_t *del, halfedge_t *d ) 766 { 767 halfedge_t *curr; 768 769 /* test if the halfedge has already a pointing face */ 770 if( d->face != NULL ) 771 return; 772 773 /* TODO: optimize this */ 774 del->faces = (face_t*)realloc(del->faces, (del->num_faces + 1) * sizeof(face_t)); 775 assert( NULL != del->faces ); 776 777 face_t *f = &(del->faces[del->num_faces]); 778 curr = d; 779 f->he = d; 780 f->num_verts = 0; 781 do { 782 curr->face = f; 783 (f->num_verts)++; 784 curr = curr->pair->prev; 785 } while( curr != d ); 786 787 (del->num_faces)++; 788 } 789 790 /* 791 * build the faces for all the halfedge 792 */ 793 void del_build_faces( delaunay_t *del ) 794 { 795 unsigned int i; 796 halfedge_t *curr; 797 798 del->num_faces = 0; 799 del->faces = NULL; 800 801 /* build external face first */ 802 build_halfedge_face(del, del->rightmost_he->pair); 803 804 for( i = del->start_point; i <= del->end_point; i++ ) 805 { 806 curr = del->points[i].he; 807 808 do { 809 build_halfedge_face( del, curr ); 810 curr = curr->next; 811 } while( curr != del->points[i].he ); 812 } 813 } 814 815 /* 816 */ 817 delaunay2d_t* delaunay2d_from(del_point2d_t *points, unsigned int num_points) { 818 delaunay2d_t* res = NULL; 819 delaunay_t del; 820 unsigned int i, j, fbuff_size = 0; 821 unsigned int* faces = NULL; 822 823 /* allocate the points */ 824 del.points = (point2d_t*)malloc(num_points * sizeof(point2d_t)); 825 assert( NULL != del.points ); 826 memset(del.points, 0, num_points * sizeof(point2d_t)); 827 828 /* copy the points */ 829 for( i = 0; i < num_points; i++ ) 830 { 831 del.points[i].idx = i; 832 del.points[i].x = points[i].x; 833 del.points[i].y = points[i].y; 834 } 835 836 qsort(del.points, num_points, sizeof(point2d_t), cmp_points); 837 838 if( num_points >= 3 ) { 839 del_divide_and_conquer( &del, 0, num_points - 1 ); 840 841 del_build_faces( &del ); 842 843 fbuff_size = 0; 844 for( i = 0; i < del.num_faces; i++ ) 845 fbuff_size += del.faces[i].num_verts + 1; 846 847 faces = (unsigned int*)malloc(sizeof(unsigned int) * fbuff_size); 848 assert( NULL != faces ); 849 850 j = 0; 851 for( i = 0; i < del.num_faces; i++ ) 852 { 853 halfedge_t *curr; 854 855 faces[j] = del.faces[i].num_verts; 856 j++; 857 858 curr = del.faces[i].he; 859 do { 860 faces[j] = curr->vertex->idx; 861 j++; 862 curr = curr->pair->prev; 863 } while( curr != del.faces[i].he ); 864 } 865 866 del_free_halfedges( &del ); 867 868 free(del.faces); 869 free(del.points); 870 } 871 872 res = (delaunay2d_t*)malloc(sizeof(delaunay2d_t)); 873 assert( NULL != res ); 874 res->num_points = num_points; 875 res->points = (del_point2d_t*)malloc(sizeof(del_point2d_t) * num_points); 876 assert( NULL != res->points ); 877 memcpy(res->points, points, sizeof(del_point2d_t) * num_points); 878 res->num_faces = del.num_faces; 879 res->faces = faces; 880 881 return res; 882 } 883 884 void delaunay2d_release(delaunay2d_t *del) { 885 free(del->faces); 886 free(del->points); 887 free(del); 888 } 889 890 891 tri_delaunay2d_t* tri_delaunay2d_from(delaunay2d_t* del) { 892 unsigned int v_offset = del->faces[0] + 1; /* ignore external face */ 893 unsigned int dst_offset = 0; 894 unsigned int i; 895 896 tri_delaunay2d_t* tdel = (tri_delaunay2d_t*)malloc(sizeof(tri_delaunay2d_t)); 897 assert( NULL != tdel ); 898 tdel->num_triangles = 0; 899 900 /* count the number of triangles */ 901 if( 1 == del->num_faces ) { /* degenerate case: only external face exists */ 902 unsigned int nv = del->faces[0]; 903 tdel->num_triangles += nv - 2; 904 } else { 905 for( i = 1; i < del->num_faces; ++i ) { 906 unsigned int nv = del->faces[v_offset]; 907 tdel->num_triangles += nv - 2; 908 v_offset += nv + 1; 909 } 910 } 911 912 /* copy points */ 913 tdel->num_points = del->num_points; 914 tdel->points = (del_point2d_t*)malloc(sizeof(del_point2d_t) * del->num_points); 915 assert( NULL != tdel->points ); 916 memcpy(tdel->points, del->points, sizeof(del_point2d_t) * del->num_points); 917 918 /* build the triangles */ 919 tdel->tris = (unsigned int*)malloc(sizeof(unsigned int) * 3 * tdel->num_triangles); 920 assert( NULL != tdel->tris ); 921 922 v_offset = del->faces[0] + 1; /* ignore external face */ 923 924 if( 1 == del->num_faces ) { 925 /* handle the degenerated case where only the external face exists */ 926 unsigned int nv = del->faces[0]; 927 unsigned int j = 0; 928 v_offset = 1; 929 for( ; j < nv - 2; ++j ) { 930 tdel->tris[dst_offset] = del->faces[v_offset + j]; 931 tdel->tris[dst_offset + 1] = del->faces[(v_offset + j + 1) % nv]; 932 tdel->tris[dst_offset + 2] = del->faces[v_offset + j]; 933 dst_offset += 3; 934 } 935 } else { 936 for( i = 1; i < del->num_faces; ++i ) { 937 unsigned int nv = del->faces[v_offset]; 938 unsigned int j = 0; 939 unsigned int first = del->faces[v_offset + 1]; 940 941 942 for( ; j < nv - 2; ++j ) { 943 tdel->tris[dst_offset] = first; 944 tdel->tris[dst_offset + 1] = del->faces[v_offset + j + 2]; 945 tdel->tris[dst_offset + 2] = del->faces[v_offset + j + 3]; 946 dst_offset += 3; 947 } 948 949 v_offset += nv + 1; 950 } 951 } 952 953 return tdel; 954 } 955 956 957 void tri_delaunay2d_release(tri_delaunay2d_t* tdel) { 958 free(tdel->tris); 959 free(tdel->points); 960 free(tdel); 961 }