201 #define FILENAMESIZE 2048
206 #define INPUTLINESIZE 1024
212 #define TRIPERBLOCK 4092
213 #define SUBSEGPERBLOCK 508
214 #define VERTEXPERBLOCK 4092
215 #define VIRUSPERBLOCK 1020
216 #define BADSUBSEGPERBLOCK 252
217 #define BADTRIPERBLOCK 4092
218 #define FLIPSTACKERPERBLOCK 252
219 #define SPLAYNODEPERBLOCK 508
225 #define INPUTVERTEX 0
226 #define SEGMENTVERTEX 1
228 #define DEADVERTEX -32768
229 #define UNDEADVERTEX -32767
237 #define SAMPLEFACTOR 11
243 #define SAMPLERATE 10
247 #define PI 3.141592653589793238462643383279502884197169399375105820974944592308
251 #define SQUAREROOTTWO 1.4142135623730950488016887242096980785696718753769480732
255 #define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333
262 #include "triangle.h"
268 enum locateresult {INTRIANGLE, ONEDGE, ONVERTEX, OUTSIDE};
276 enum insertvertexresult {SUCCESSFULVERTEX, ENCROACHINGVERTEX, VIOLATINGVERTEX,
284 enum finddirectionresult {WITHIN, LEFTCOLLINEAR, RIGHTCOLLINEAR};
403 typedef float **triangle;
420 typedef float **subseg;
437 typedef float *vertex;
444 vertex subsegorg, subsegdest;
454 vertex triangorg, triangdest, triangapex;
455 struct badtriang *nexttriang;
464 struct flipstacker *prevflip;
497 struct splaynode *lchild, *rchild;
524 int **firstblock, **nowblock;
533 long items, maxitems;
534 int unallocateditems;
543 float resulterrbound;
544 float ccwerrboundA, ccwerrboundB, ccwerrboundC;
545 float iccerrboundA, iccerrboundB, iccerrboundC;
546 float o3derrboundA, o3derrboundB, o3derrboundC;
550 unsigned long randomseed;
562 struct memorypool triangles;
563 struct memorypool subsegs;
564 struct memorypool vertices;
565 struct memorypool viri;
566 struct memorypool badsubsegs;
567 struct memorypool badtriangles;
568 struct memorypool flipstackers;
569 struct memorypool splaynodes;
574 struct badtriang *queuefront[4096];
575 struct badtriang *queuetail[4096];
576 int nextnonemptyq[4096];
581 struct flipstacker *lastflip;
585 float xmin, xmax, ymin, ymax;
610 long counterclockcount;
613 long circumcentercount;
618 vertex infvertex1, infvertex2, infvertex3;
623 triangle *dummytribase;
630 subseg *dummysubbase;
635 struct otri recenttri;
678 int poly, refine, quality, vararea, fixedarea, usertest;
679 int regionattrib, convex, weighted, jettison;
681 int edgesout, voronoi, neighbors, geomview;
682 int nobound, nopolywritten, nonodewritten, noelewritten, noiterationnum;
683 int noholes, noexact, conformdel;
684 int incremental, sweepline, dwyer;
692 float minangle, goodangle, offconstant;
812 int plus1mod3[3] = {1, 2, 0};
813 int minus1mod3[3] = {2, 0, 1};
822 #define decode(ptr, otri) \
823 (otri).orient = (int) ((unsigned long) (ptr) & (unsigned long) 3l); \
824 (otri).tri = (triangle *) \
825 ((unsigned long) (ptr) ^ (unsigned long) (otri).orient)
831 #define encode(otri) \
832 (triangle) ((unsigned long) (otri).tri | (unsigned long) (otri).orient)
842 #define sym(otri1, otri2) \
843 ptr = (otri1).tri[(otri1).orient]; \
846 #define symself(otri) \
847 ptr = (otri).tri[(otri).orient]; \
852 #define lnext(otri1, otri2) \
853 (otri2).tri = (otri1).tri; \
854 (otri2).orient = plus1mod3[(otri1).orient]
856 #define lnextself(otri) \
857 (otri).orient = plus1mod3[(otri).orient]
861 #define lprev(otri1, otri2) \
862 (otri2).tri = (otri1).tri; \
863 (otri2).orient = minus1mod3[(otri1).orient]
865 #define lprevself(otri) \
866 (otri).orient = minus1mod3[(otri).orient]
872 #define onext(otri1, otri2) \
873 lprev(otri1, otri2); \
876 #define onextself(otri) \
884 #define oprev(otri1, otri2) \
888 #define oprevself(otri) \
896 #define dnext(otri1, otri2) \
900 #define dnextself(otri) \
908 #define dprev(otri1, otri2) \
909 lnext(otri1, otri2); \
912 #define dprevself(otri) \
920 #define rnext(otri1, otri2) \
925 #define rnextself(otri) \
934 #define rprev(otri1, otri2) \
939 #define rprevself(otri) \
947 #define org(otri, vertexptr) \
948 vertexptr = (vertex) (otri).tri[plus1mod3[(otri).orient] + 3]
950 #define dest(otri, vertexptr) \
951 vertexptr = (vertex) (otri).tri[minus1mod3[(otri).orient] + 3]
953 #define apex(otri, vertexptr) \
954 vertexptr = (vertex) (otri).tri[(otri).orient + 3]
956 #define setorg(otri, vertexptr) \
957 (otri).tri[plus1mod3[(otri).orient] + 3] = (triangle) vertexptr
959 #define setdest(otri, vertexptr) \
960 (otri).tri[minus1mod3[(otri).orient] + 3] = (triangle) vertexptr
962 #define setapex(otri, vertexptr) \
963 (otri).tri[(otri).orient + 3] = (triangle) vertexptr
967 #define bond(otri1, otri2) \
968 (otri1).tri[(otri1).orient] = encode(otri2); \
969 (otri2).tri[(otri2).orient] = encode(otri1)
976 #define dissolve(otri) \
977 (otri).tri[(otri).orient] = (triangle) m->dummytri
981 #define otricopy(otri1, otri2) \
982 (otri2).tri = (otri1).tri; \
983 (otri2).orient = (otri1).orient
987 #define otriequal(otri1, otri2) \
988 (((otri1).tri == (otri2).tri) && \
989 ((otri1).orient == (otri2).orient))
994 #define infect(otri) \
995 (otri).tri[6] = (triangle) \
996 ((unsigned long) (otri).tri[6] | (unsigned long) 2l)
998 #define uninfect(otri) \
999 (otri).tri[6] = (triangle) \
1000 ((unsigned long) (otri).tri[6] & ~ (unsigned long) 2l)
1004 #define infected(otri) \
1005 (((unsigned long) (otri).tri[6] & (unsigned long) 2l) != 0l)
1009 #define elemattribute(otri, attnum) \
1010 ((float *) (otri).tri)[m->elemattribindex + (attnum)]
1012 #define setelemattribute(otri, attnum, value) \
1013 ((float *) (otri).tri)[m->elemattribindex + (attnum)] = value
1017 #define areabound(otri) ((float *) (otri).tri)[m->areaboundindex]
1019 #define setareabound(otri, value) \
1020 ((float *) (otri).tri)[m->areaboundindex] = value
1027 #define deadtri(tria) ((tria)[1] == (triangle) NULL)
1029 #define killtri(tria) \
1030 (tria)[1] = (triangle) NULL; \
1031 (tria)[3] = (triangle) NULL
1042 #define sdecode(sptr, osub) \
1043 (osub).ssorient = (int) ((unsigned long) (sptr) & (unsigned long) 1l); \
1044 (osub).ss = (subseg *) \
1045 ((unsigned long) (sptr) & ~ (unsigned long) 3l)
1051 #define sencode(osub) \
1052 (subseg) ((unsigned long) (osub).ss | (unsigned long) (osub).ssorient)
1056 #define ssym(osub1, osub2) \
1057 (osub2).ss = (osub1).ss; \
1058 (osub2).ssorient = 1 - (osub1).ssorient
1060 #define ssymself(osub) \
1061 (osub).ssorient = 1 - (osub).ssorient
1066 #define spivot(osub1, osub2) \
1067 sptr = (osub1).ss[(osub1).ssorient]; \
1068 sdecode(sptr, osub2)
1070 #define spivotself(osub) \
1071 sptr = (osub).ss[(osub).ssorient]; \
1077 #define snext(osub1, osub2) \
1078 sptr = (osub1).ss[1 - (osub1).ssorient]; \
1079 sdecode(sptr, osub2)
1081 #define snextself(osub) \
1082 sptr = (osub).ss[1 - (osub).ssorient]; \
1088 #define sorg(osub, vertexptr) \
1089 vertexptr = (vertex) (osub).ss[2 + (osub).ssorient]
1091 #define sdest(osub, vertexptr) \
1092 vertexptr = (vertex) (osub).ss[3 - (osub).ssorient]
1094 #define setsorg(osub, vertexptr) \
1095 (osub).ss[2 + (osub).ssorient] = (subseg) vertexptr
1097 #define setsdest(osub, vertexptr) \
1098 (osub).ss[3 - (osub).ssorient] = (subseg) vertexptr
1100 #define segorg(osub, vertexptr) \
1101 vertexptr = (vertex) (osub).ss[4 + (osub).ssorient]
1103 #define segdest(osub, vertexptr) \
1104 vertexptr = (vertex) (osub).ss[5 - (osub).ssorient]
1106 #define setsegorg(osub, vertexptr) \
1107 (osub).ss[4 + (osub).ssorient] = (subseg) vertexptr
1109 #define setsegdest(osub, vertexptr) \
1110 (osub).ss[5 - (osub).ssorient] = (subseg) vertexptr
1116 #define mark(osub) (* (int *) ((osub).ss + 8))
1118 #define setmark(osub, value) \
1119 * (int *) ((osub).ss + 8) = value
1123 #define sbond(osub1, osub2) \
1124 (osub1).ss[(osub1).ssorient] = sencode(osub2); \
1125 (osub2).ss[(osub2).ssorient] = sencode(osub1)
1130 #define sdissolve(osub) \
1131 (osub).ss[(osub).ssorient] = (subseg) m->dummysub
1135 #define subsegcopy(osub1, osub2) \
1136 (osub2).ss = (osub1).ss; \
1137 (osub2).ssorient = (osub1).ssorient
1141 #define subsegequal(osub1, osub2) \
1142 (((osub1).ss == (osub2).ss) && \
1143 ((osub1).ssorient == (osub2).ssorient))
1150 #define deadsubseg(sub) ((sub)[1] == (subseg) NULL)
1152 #define killsubseg(sub) \
1153 (sub)[1] = (subseg) NULL; \
1154 (sub)[2] = (subseg) NULL
1162 #define tspivot(otri, osub) \
1163 sptr = (subseg) (otri).tri[6 + (otri).orient]; \
1169 #define stpivot(osub, otri) \
1170 ptr = (triangle) (osub).ss[6 + (osub).ssorient]; \
1175 #define tsbond(otri, osub) \
1176 (otri).tri[6 + (otri).orient] = (triangle) sencode(osub); \
1177 (osub).ss[6 + (osub).ssorient] = (subseg) encode(otri)
1181 #define tsdissolve(otri) \
1182 (otri).tri[6 + (otri).orient] = (triangle) m->dummysub
1186 #define stdissolve(osub) \
1187 (osub).ss[6 + (osub).ssorient] = (subseg) m->dummytri
1193 #define vertexmark(vx) ((int *) (vx))[m->vertexmarkindex]
1195 #define setvertexmark(vx, value) \
1196 ((int *) (vx))[m->vertexmarkindex] = value
1198 #define vertextype(vx) ((int *) (vx))[m->vertexmarkindex + 1]
1200 #define setvertextype(vx, value) \
1201 ((int *) (vx))[m->vertexmarkindex + 1] = value
1203 #define vertex2tri(vx) ((triangle *) (vx))[m->vertex2triindex]
1205 #define setvertex2tri(vx, value) \
1206 ((triangle *) (vx))[m->vertex2triindex] = value
1216 void triexit(
int status)
1221 int *trimalloc(
int size)
1225 memptr = (
int *) malloc((
unsigned int) size);
1226 if (memptr == (
int *) NULL) {
1227 printf(
"Error: Out of memory.\n");
1233 void trifree(
int *memptr)
1248 void internalerror()
1250 printf(
" Please report this bug to jrs@cs.berkeley.edu\n");
1251 printf(
" Include the message above, your input data set, and the exact\n");
1252 printf(
" command line you used to run Triangle.\n");
1263 void parsecommandline(
int argc,
char **argv,
struct behavior *b) {
1265 char workstring[FILENAMESIZE];
1267 b->poly = b->refine = b->quality = 0;
1268 b->vararea = b->fixedarea = b->usertest = 0;
1269 b->regionattrib = b->convex = b->weighted = b->jettison = 0;
1271 b->edgesout = b->voronoi = b->neighbors = b->geomview = 0;
1272 b->nobound = b->nopolywritten = b->nonodewritten = b->noelewritten = 0;
1273 b->noiterationnum = 0;
1274 b->noholes = b->noexact = 0;
1275 b->incremental = b->sweepline = 0;
1285 b->quiet = b->verbose = 0;
1287 for (i = 0; i < argc; i++) {
1288 for (j = 0; argv[i][j] !=
'\0'; j++) {
1289 if (argv[i][j] ==
'p') {
1292 if (argv[i][j] ==
'A') {
1293 b->regionattrib = 1;
1295 if (argv[i][j] ==
'c') {
1298 if (argv[i][j] ==
'w') {
1301 if (argv[i][j] ==
'W') {
1304 if (argv[i][j] ==
'j') {
1307 if (argv[i][j] ==
'z') {
1310 if (argv[i][j] ==
'e') {
1313 if (argv[i][j] ==
'v') {
1316 if (argv[i][j] ==
'n') {
1319 if (argv[i][j] ==
'g') {
1322 if (argv[i][j] ==
'B') {
1325 if (argv[i][j] ==
'P') {
1326 b->nopolywritten = 1;
1328 if (argv[i][j] ==
'N') {
1329 b->nonodewritten = 1;
1331 if (argv[i][j] ==
'E') {
1332 b->noelewritten = 1;
1334 if (argv[i][j] ==
'O') {
1337 if (argv[i][j] ==
'X') {
1340 if (argv[i][j] ==
'o') {
1341 if (argv[i][j + 1] ==
'2') {
1346 if (argv[i][j] ==
'l') {
1349 if (argv[i][j] ==
'Q') {
1352 if (argv[i][j] ==
'V') {
1357 b->usesegments = b->poly || b->refine || b->quality || b->convex;
1358 b->goodangle = cos(b->minangle * PI / 180.0);
1359 if (b->goodangle == 1.0) {
1360 b->offconstant = 0.0;
1362 b->offconstant = 0.475 * sqrt((1.0 + b->goodangle) / (1.0 - b->goodangle));
1364 b->goodangle *= b->goodangle;
1365 if (b->refine && b->noiterationnum) {
1367 "Error: You cannot use the -I switch when refining a triangulation.\n");
1372 if (!b->refine && !b->poly) {
1377 if (b->refine || !b->poly) {
1378 b->regionattrib = 0;
1382 if (b->weighted && (b->poly || b->quality)) {
1385 printf(
"Warning: weighted triangulations (-w, -W) are incompatible\n");
1386 printf(
" with PSLGs (-p) and meshing (-q, -a, -u). Weights ignored.\n"
1390 if (b->jettison && b->nonodewritten && !b->quiet) {
1391 printf(
"Warning: -j and -N switches are somewhat incompatible.\n");
1392 printf(
" If any vertices are jettisoned, you will need the output\n");
1393 printf(
" .node file to reconstruct the new node indices.");
1416 void printtriangle(
struct mesh *m,
struct behavior *b,
struct otri *t)
1418 struct otri printtri;
1419 struct osub printsh;
1422 printf(
"triangle x%lx with orientation %d:\n", (
unsigned long) t->tri,
1424 decode(t->tri[0], printtri);
1425 if (printtri.tri == m->dummytri) {
1426 printf(
" [0] = Outer space\n");
1428 printf(
" [0] = x%lx %d\n", (
unsigned long) printtri.tri,
1431 decode(t->tri[1], printtri);
1432 if (printtri.tri == m->dummytri) {
1433 printf(
" [1] = Outer space\n");
1435 printf(
" [1] = x%lx %d\n", (
unsigned long) printtri.tri,
1438 decode(t->tri[2], printtri);
1439 if (printtri.tri == m->dummytri) {
1440 printf(
" [2] = Outer space\n");
1442 printf(
" [2] = x%lx %d\n", (
unsigned long) printtri.tri,
1446 org(*t, printvertex);
1447 if (printvertex == (vertex) NULL)
1448 printf(
" Origin[%d] = NULL\n", (t->orient + 1) % 3 + 3);
1450 printf(
" Origin[%d] = x%lx (%.12g, %.12g)\n",
1451 (t->orient + 1) % 3 + 3, (
unsigned long) printvertex,
1452 printvertex[0], printvertex[1]);
1453 dest(*t, printvertex);
1454 if (printvertex == (vertex) NULL)
1455 printf(
" Dest [%d] = NULL\n", (t->orient + 2) % 3 + 3);
1457 printf(
" Dest [%d] = x%lx (%.12g, %.12g)\n",
1458 (t->orient + 2) % 3 + 3, (
unsigned long) printvertex,
1459 printvertex[0], printvertex[1]);
1460 apex(*t, printvertex);
1461 if (printvertex == (vertex) NULL)
1462 printf(
" Apex [%d] = NULL\n", t->orient + 3);
1464 printf(
" Apex [%d] = x%lx (%.12g, %.12g)\n",
1465 t->orient + 3, (
unsigned long) printvertex,
1466 printvertex[0], printvertex[1]);
1468 if (b->usesegments) {
1469 sdecode(t->tri[6], printsh);
1470 if (printsh.ss != m->dummysub) {
1471 printf(
" [6] = x%lx %d\n", (
unsigned long) printsh.ss,
1474 sdecode(t->tri[7], printsh);
1475 if (printsh.ss != m->dummysub) {
1476 printf(
" [7] = x%lx %d\n", (
unsigned long) printsh.ss,
1479 sdecode(t->tri[8], printsh);
1480 if (printsh.ss != m->dummysub) {
1481 printf(
" [8] = x%lx %d\n", (
unsigned long) printsh.ss,
1487 printf(
" Area constraint: %.4g\n", areabound(*t));
1502 void printsubseg(
struct mesh *m,
struct behavior *b,
struct osub *s)
1504 struct osub printsh;
1505 struct otri printtri;
1508 printf(
"subsegment x%lx with orientation %d and mark %d:\n",
1509 (
unsigned long) s->ss, s->ssorient, mark(*s));
1510 sdecode(s->ss[0], printsh);
1511 if (printsh.ss == m->dummysub) {
1512 printf(
" [0] = No subsegment\n");
1514 printf(
" [0] = x%lx %d\n", (
unsigned long) printsh.ss,
1517 sdecode(s->ss[1], printsh);
1518 if (printsh.ss == m->dummysub) {
1519 printf(
" [1] = No subsegment\n");
1521 printf(
" [1] = x%lx %d\n", (
unsigned long) printsh.ss,
1525 sorg(*s, printvertex);
1526 if (printvertex == (vertex) NULL)
1527 printf(
" Origin[%d] = NULL\n", 2 + s->ssorient);
1529 printf(
" Origin[%d] = x%lx (%.12g, %.12g)\n",
1530 2 + s->ssorient, (
unsigned long) printvertex,
1531 printvertex[0], printvertex[1]);
1532 sdest(*s, printvertex);
1533 if (printvertex == (vertex) NULL)
1534 printf(
" Dest [%d] = NULL\n", 3 - s->ssorient);
1536 printf(
" Dest [%d] = x%lx (%.12g, %.12g)\n",
1537 3 - s->ssorient, (
unsigned long) printvertex,
1538 printvertex[0], printvertex[1]);
1540 decode(s->ss[6], printtri);
1541 if (printtri.tri == m->dummytri) {
1542 printf(
" [6] = Outer space\n");
1544 printf(
" [6] = x%lx %d\n", (
unsigned long) printtri.tri,
1547 decode(s->ss[7], printtri);
1548 if (printtri.tri == m->dummytri) {
1549 printf(
" [7] = Outer space\n");
1551 printf(
" [7] = x%lx %d\n", (
unsigned long) printtri.tri,
1555 segorg(*s, printvertex);
1556 if (printvertex == (vertex) NULL)
1557 printf(
" Segment origin[%d] = NULL\n", 4 + s->ssorient);
1559 printf(
" Segment origin[%d] = x%lx (%.12g, %.12g)\n",
1560 4 + s->ssorient, (
unsigned long) printvertex,
1561 printvertex[0], printvertex[1]);
1562 segdest(*s, printvertex);
1563 if (printvertex == (vertex) NULL)
1564 printf(
" Segment dest [%d] = NULL\n", 5 - s->ssorient);
1566 printf(
" Segment dest [%d] = x%lx (%.12g, %.12g)\n",
1567 5 - s->ssorient, (
unsigned long) printvertex,
1568 printvertex[0], printvertex[1]);
1588 void poolzero(
struct memorypool *pool)
1590 pool->firstblock = (
int **) NULL;
1591 pool->nowblock = (
int **) NULL;
1592 pool->nextitem = (
int *) NULL;
1593 pool->deaditemstack = (
int *) NULL;
1594 pool->pathblock = (
int **) NULL;
1595 pool->pathitem = (
int *) NULL;
1596 pool->alignbytes = 0;
1597 pool->itembytes = 0;
1598 pool->itemsperblock = 0;
1599 pool->itemsfirstblock = 0;
1602 pool->unallocateditems = 0;
1603 pool->pathitemsleft = 0;
1616 void poolrestart(
struct memorypool *pool)
1618 unsigned long alignptr;
1624 pool->nowblock = pool->firstblock;
1626 alignptr = (
unsigned long) (pool->nowblock + 1);
1628 pool->nextitem = (
int *)
1629 (alignptr + (
unsigned long) pool->alignbytes -
1630 (alignptr % (
unsigned long) pool->alignbytes));
1632 pool->unallocateditems = pool->itemsfirstblock;
1634 pool->deaditemstack = (
int *) NULL;
1656 void poolinit(
struct memorypool *pool,
int bytecount,
int itemcount,
1657 int firstitemcount,
int alignment)
1663 if (alignment >
sizeof(
int *)) {
1664 pool->alignbytes = alignment;
1666 pool->alignbytes =
sizeof(
int *);
1668 pool->itembytes = ((bytecount - 1) / pool->alignbytes + 1) *
1670 pool->itemsperblock = itemcount;
1671 if (firstitemcount == 0) {
1672 pool->itemsfirstblock = itemcount;
1674 pool->itemsfirstblock = firstitemcount;
1680 pool->firstblock = (
int **)
1681 trimalloc(pool->itemsfirstblock * pool->itembytes + (
int)
sizeof(
int *) +
1684 *(pool->firstblock) = (
int *) NULL;
1694 void pooldeinit(
struct memorypool *pool)
1696 while (pool->firstblock != (
int **) NULL) {
1697 pool->nowblock = (
int **) *(pool->firstblock);
1698 trifree((
int *) pool->firstblock);
1699 pool->firstblock = pool->nowblock;
1709 int *poolalloc(
struct memorypool *pool)
1713 unsigned long alignptr;
1717 if (pool->deaditemstack != (
int *) NULL) {
1718 newitem = pool->deaditemstack;
1719 pool->deaditemstack = * (
int **) pool->deaditemstack;
1722 if (pool->unallocateditems == 0) {
1724 if (*(pool->nowblock) == (
int *) NULL) {
1726 newblock = (
int **) trimalloc(pool->itemsperblock * pool->itembytes +
1727 (
int)
sizeof(
int *) +
1729 *(pool->nowblock) = (
int *) newblock;
1731 *newblock = (
int *) NULL;
1735 pool->nowblock = (
int **) *(pool->nowblock);
1738 alignptr = (
unsigned long) (pool->nowblock + 1);
1740 pool->nextitem = (
int *)
1741 (alignptr + (
unsigned long) pool->alignbytes -
1742 (alignptr % (
unsigned long) pool->alignbytes));
1744 pool->unallocateditems = pool->itemsperblock;
1748 newitem = pool->nextitem;
1750 pool->nextitem = (
int *) ((
char *) pool->nextitem + pool->itembytes);
1751 pool->unallocateditems--;
1766 void pooldealloc(
struct memorypool *pool,
int *dyingitem)
1769 *((
int **) dyingitem) = pool->deaditemstack;
1770 pool->deaditemstack = dyingitem;
1782 void traversalinit(
struct memorypool *pool)
1784 unsigned long alignptr;
1787 pool->pathblock = pool->firstblock;
1789 alignptr = (
unsigned long) (pool->pathblock + 1);
1791 pool->pathitem = (
int *)
1792 (alignptr + (
unsigned long) pool->alignbytes -
1793 (alignptr % (
unsigned long) pool->alignbytes));
1795 pool->pathitemsleft = pool->itemsfirstblock;
1812 int *traverse(
struct memorypool *pool)
1815 unsigned long alignptr;
1818 if (pool->pathitem == pool->nextitem) {
1819 return (
int *) NULL;
1823 if (pool->pathitemsleft == 0) {
1825 pool->pathblock = (
int **) *(pool->pathblock);
1827 alignptr = (
unsigned long) (pool->pathblock + 1);
1829 pool->pathitem = (
int *)
1830 (alignptr + (
unsigned long) pool->alignbytes -
1831 (alignptr % (
unsigned long) pool->alignbytes));
1833 pool->pathitemsleft = pool->itemsperblock;
1836 newitem = pool->pathitem;
1838 pool->pathitem = (
int *) ((
char *) pool->pathitem + pool->itembytes);
1839 pool->pathitemsleft--;
1871 void dummyinit(
struct mesh *m,
struct behavior *b,
int trianglebytes,
1874 unsigned long alignptr;
1877 m->dummytribase = (triangle *) trimalloc(trianglebytes +
1878 m->triangles.alignbytes);
1880 alignptr = (
unsigned long) m->dummytribase;
1881 m->dummytri = (triangle *)
1882 (alignptr + (
unsigned long) m->triangles.alignbytes -
1883 (alignptr % (
unsigned long) m->triangles.alignbytes));
1888 m->dummytri[0] = (triangle) m->dummytri;
1889 m->dummytri[1] = (triangle) m->dummytri;
1890 m->dummytri[2] = (triangle) m->dummytri;
1892 m->dummytri[3] = (triangle) NULL;
1893 m->dummytri[4] = (triangle) NULL;
1894 m->dummytri[5] = (triangle) NULL;
1896 if (b->usesegments) {
1900 m->dummysubbase = (subseg *) trimalloc(subsegbytes +
1901 m->subsegs.alignbytes);
1903 alignptr = (
unsigned long) m->dummysubbase;
1904 m->dummysub = (subseg *)
1905 (alignptr + (
unsigned long) m->subsegs.alignbytes -
1906 (alignptr % (
unsigned long) m->subsegs.alignbytes));
1911 m->dummysub[0] = (subseg) m->dummysub;
1912 m->dummysub[1] = (subseg) m->dummysub;
1914 m->dummysub[2] = (subseg) NULL;
1915 m->dummysub[3] = (subseg) NULL;
1916 m->dummysub[4] = (subseg) NULL;
1917 m->dummysub[5] = (subseg) NULL;
1919 m->dummysub[6] = (subseg) m->dummytri;
1920 m->dummysub[7] = (subseg) m->dummytri;
1922 * (
int *) (m->dummysub + 8) = 0;
1926 m->dummytri[6] = (triangle) m->dummysub;
1927 m->dummytri[7] = (triangle) m->dummysub;
1928 m->dummytri[8] = (triangle) m->dummysub;
1942 void initializevertexpool(
struct mesh *m,
struct behavior *b)
1949 m->vertexmarkindex = ((m->mesh_dim + m->nextras) *
sizeof(
float) +
1952 vertexsize = (m->vertexmarkindex + 2) *
sizeof(
int);
1956 m->vertex2triindex = (vertexsize +
sizeof(triangle) - 1) /
1958 vertexsize = (m->vertex2triindex + 1) *
sizeof(triangle);
1962 poolinit(&m->vertices, vertexsize, VERTEXPERBLOCK,
1963 m->invertices > VERTEXPERBLOCK ? m->invertices : VERTEXPERBLOCK,
1978 void initializetrisubpools(
struct mesh *m,
struct behavior *b)
1986 m->highorderindex = 6 + (b->usesegments * 3);
1988 trisize = ((b->order + 1) * (b->order + 2) / 2 + (m->highorderindex - 3)) *
1992 m->elemattribindex = (trisize +
sizeof(float) - 1) /
sizeof(float);
1996 m->areaboundindex = m->elemattribindex + m->eextras + b->regionattrib;
2000 trisize = (m->areaboundindex + 1) *
sizeof(
float);
2001 }
else if (m->eextras + b->regionattrib > 0) {
2002 trisize = m->areaboundindex *
sizeof(float);
2008 if ((b->voronoi || b->neighbors) &&
2009 (trisize < 6 *
sizeof(triangle) +
sizeof(
int))) {
2010 trisize = 6 *
sizeof(triangle) +
sizeof(
int);
2014 poolinit(&m->triangles, trisize, TRIPERBLOCK,
2015 (2 * m->invertices - 2) > TRIPERBLOCK ? (2 * m->invertices - 2) :
2018 if (b->usesegments) {
2021 poolinit(&m->subsegs, 8 *
sizeof(triangle) +
sizeof(
int),
2022 SUBSEGPERBLOCK, SUBSEGPERBLOCK, 4);
2025 dummyinit(m, b, m->triangles.itembytes, m->subsegs.itembytes);
2028 dummyinit(m, b, m->triangles.itembytes, 0);
2038 void triangledealloc(
struct mesh *m, triangle *dyingtriangle)
2042 killtri(dyingtriangle);
2043 pooldealloc(&m->triangles, (
int *) dyingtriangle);
2052 triangle *triangletraverse(
struct mesh *m)
2054 triangle *newtriangle;
2057 newtriangle = (triangle *) traverse(&m->triangles);
2058 if (newtriangle == (triangle *) NULL) {
2059 return (triangle *) NULL;
2061 }
while (deadtri(newtriangle));
2071 void subsegdealloc(
struct mesh *m, subseg *dyingsubseg)
2075 killsubseg(dyingsubseg);
2076 pooldealloc(&m->subsegs, (
int *) dyingsubseg);
2085 subseg *subsegtraverse(
struct mesh *m)
2090 newsubseg = (subseg *) traverse(&m->subsegs);
2091 if (newsubseg == (subseg *) NULL) {
2092 return (subseg *) NULL;
2094 }
while (deadsubseg(newsubseg));
2104 void vertexdealloc(
struct mesh *m, vertex dyingvertex)
2108 setvertextype(dyingvertex, DEADVERTEX);
2109 pooldealloc(&m->vertices, (
int *) dyingvertex);
2118 vertex vertextraverse(
struct mesh *m)
2123 newvertex = (vertex) traverse(&m->vertices);
2124 if (newvertex == (vertex) NULL) {
2125 return (vertex) NULL;
2127 }
while (vertextype(newvertex) == DEADVERTEX);
2143 vertex getvertex(
struct mesh *m,
struct behavior *b,
int number)
2147 unsigned long alignptr;
2150 getblock = m->vertices.firstblock;
2151 current = b->firstnumber;
2154 if (current + m->vertices.itemsfirstblock <= number) {
2155 getblock = (
int **) *getblock;
2156 current += m->vertices.itemsfirstblock;
2157 while (current + m->vertices.itemsperblock <= number) {
2158 getblock = (
int **) *getblock;
2159 current += m->vertices.itemsperblock;
2164 alignptr = (
unsigned long) (getblock + 1);
2165 foundvertex = (
char *) (alignptr + (
unsigned long) m->vertices.alignbytes -
2166 (alignptr % (
unsigned long) m->vertices.alignbytes));
2167 return (vertex) (foundvertex + m->vertices.itembytes * (number - current));
2176 void triangledeinit(
struct mesh *m,
struct behavior *b)
2178 pooldeinit(&m->triangles);
2179 trifree((
int *) m->dummytribase);
2180 if (b->usesegments) {
2181 pooldeinit(&m->subsegs);
2182 trifree((
int *) m->dummysubbase);
2184 pooldeinit(&m->vertices);
2201 void maketriangle(
struct mesh *m,
struct behavior *b,
struct otri *newotri)
2205 newotri->tri = (triangle *) poolalloc(&m->triangles);
2207 newotri->tri[0] = (triangle) m->dummytri;
2208 newotri->tri[1] = (triangle) m->dummytri;
2209 newotri->tri[2] = (triangle) m->dummytri;
2211 newotri->tri[3] = (triangle) NULL;
2212 newotri->tri[4] = (triangle) NULL;
2213 newotri->tri[5] = (triangle) NULL;
2214 if (b->usesegments) {
2217 newotri->tri[6] = (triangle) m->dummysub;
2218 newotri->tri[7] = (triangle) m->dummysub;
2219 newotri->tri[8] = (triangle) m->dummysub;
2221 for (i = 0; i < m->eextras; i++) {
2222 setelemattribute(*newotri, i, 0.0);
2225 setareabound(*newotri, -1.0);
2228 newotri->orient = 0;
2237 void makesubseg(
struct mesh *m,
struct osub *newsubseg)
2239 newsubseg->ss = (subseg *) poolalloc(&m->subsegs);
2242 newsubseg->ss[0] = (subseg) m->dummysub;
2243 newsubseg->ss[1] = (subseg) m->dummysub;
2245 newsubseg->ss[2] = (subseg) NULL;
2246 newsubseg->ss[3] = (subseg) NULL;
2247 newsubseg->ss[4] = (subseg) NULL;
2248 newsubseg->ss[5] = (subseg) NULL;
2250 newsubseg->ss[6] = (subseg) m->dummytri;
2251 newsubseg->ss[7] = (subseg) m->dummytri;
2253 setmark(*newsubseg, 0);
2255 newsubseg->ssorient = 0;
2279 #define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
2295 #define Fast_Two_Sum_Tail(a, b, x, y) \
2299 #define Fast_Two_Sum(a, b, x, y) \
2300 x = (float) (a + b); \
2301 Fast_Two_Sum_Tail(a, b, x, y)
2303 #define Two_Sum_Tail(a, b, x, y) \
2304 bvirt = (float) (x - a); \
2305 avirt = x - bvirt; \
2306 bround = b - bvirt; \
2307 around = a - avirt; \
2310 #define Two_Sum(a, b, x, y) \
2311 x = (float) (a + b); \
2312 Two_Sum_Tail(a, b, x, y)
2314 #define Two_Diff_Tail(a, b, x, y) \
2315 bvirt = (float) (a - x); \
2316 avirt = x + bvirt; \
2317 bround = bvirt - b; \
2318 around = a - avirt; \
2321 #define Two_Diff(a, b, x, y) \
2322 x = (float) (a - b); \
2323 Two_Diff_Tail(a, b, x, y)
2325 #define Split(a, ahi, alo) \
2326 c = (float) (splitter * a); \
2327 abig = (float) (c - a); \
2331 #define Two_Product_Tail(a, b, x, y) \
2332 Split(a, ahi, alo); \
2333 Split(b, bhi, blo); \
2334 err1 = x - (ahi * bhi); \
2335 err2 = err1 - (alo * bhi); \
2336 err3 = err2 - (ahi * blo); \
2337 y = (alo * blo) - err3
2339 #define Two_Product(a, b, x, y) \
2340 x = (float) (a * b); \
2341 Two_Product_Tail(a, b, x, y)
2346 #define Two_Product_Presplit(a, b, bhi, blo, x, y) \
2347 x = (float) (a * b); \
2348 Split(a, ahi, alo); \
2349 err1 = x - (ahi * bhi); \
2350 err2 = err1 - (alo * bhi); \
2351 err3 = err2 - (ahi * blo); \
2352 y = (alo * blo) - err3
2356 #define Square_Tail(a, x, y) \
2357 Split(a, ahi, alo); \
2358 err1 = x - (ahi * ahi); \
2359 err3 = err1 - ((ahi + ahi) * alo); \
2360 y = (alo * alo) - err3
2362 #define Square(a, x, y) \
2363 x = (float) (a * a); \
2364 Square_Tail(a, x, y)
2369 #define Two_One_Sum(a1, a0, b, x2, x1, x0) \
2370 Two_Sum(a0, b , _i, x0); \
2371 Two_Sum(a1, _i, x2, x1)
2373 #define Two_One_Diff(a1, a0, b, x2, x1, x0) \
2374 Two_Diff(a0, b , _i, x0); \
2375 Two_Sum( a1, _i, x2, x1)
2377 #define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \
2378 Two_One_Sum(a1, a0, b0, _j, _0, x0); \
2379 Two_One_Sum(_j, _0, b1, x3, x2, x1)
2381 #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
2382 Two_One_Diff(a1, a0, b0, _j, _0, x0); \
2383 Two_One_Diff(_j, _0, b1, x3, x2, x1)
2387 #define Two_One_Product(a1, a0, b, x3, x2, x1, x0) \
2388 Split(b, bhi, blo); \
2389 Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
2390 Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
2391 Two_Sum(_i, _0, _k, x1); \
2392 Fast_Two_Sum(_j, _k, x3, x2)
2416 float check, lastcheck;
2433 every_other = !every_other;
2434 check = 1.0 + epsilon;
2435 }
while ((check != 1.0) && (check != lastcheck));
2438 resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
2439 ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;
2440 ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;
2441 ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;
2442 iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon;
2443 iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon;
2444 iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon;
2445 o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;
2446 o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;
2447 o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
2464 int fast_expansion_sum_zeroelim(
int elen,
float *e,
int flen,
float *f,
float *h)
2470 float avirt, bround, around;
2471 int eindex, findex, hindex;
2476 eindex = findex = 0;
2477 if ((fnow > enow) == (fnow > -enow)) {
2485 if ((eindex < elen) && (findex < flen)) {
2486 if ((fnow > enow) == (fnow > -enow)) {
2487 Fast_Two_Sum(enow, Q, Qnew, hh);
2490 Fast_Two_Sum(fnow, Q, Qnew, hh);
2497 while ((eindex < elen) && (findex < flen)) {
2498 if ((fnow > enow) == (fnow > -enow)) {
2499 Two_Sum(Q, enow, Qnew, hh);
2502 Two_Sum(Q, fnow, Qnew, hh);
2511 while (eindex < elen) {
2512 Two_Sum(Q, enow, Qnew, hh);
2519 while (findex < flen) {
2520 Two_Sum(Q, fnow, Qnew, hh);
2527 if ((Q != 0.0) || (hindex == 0)) {
2548 int scale_expansion_zeroelim(
int elen,
float *e,
float b,
float *h)
2557 float avirt, bround, around;
2560 float ahi, alo, bhi, blo;
2561 float err1, err2, err3;
2564 Two_Product_Presplit(e[0], b, bhi, blo, Q, hh);
2569 for (eindex = 1; eindex < elen; eindex++) {
2571 Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
2572 Two_Sum(Q, product0, sum, hh);
2576 Fast_Two_Sum(product1, sum, Q, hh);
2581 if ((Q != 0.0) || (hindex == 0)) {
2595 float estimate(
int elen,
float *e)
2600 for (eindex = 1; eindex < elen; eindex++) {
2626 float counterclockwiseadapt(vertex pa, vertex pb, vertex pc,
float detsum)
2628 float acx, acy, bcx, bcy;
2629 float acxtail, acytail, bcxtail, bcytail;
2630 float detleft, detright;
2631 float detlefttail, detrighttail;
2632 float det, errbound;
2633 float B[4], C1[8], C2[12], D[16];
2635 int C1length, C2length, Dlength;
2642 float avirt, bround, around;
2645 float ahi, alo, bhi, blo;
2646 float err1, err2, err3;
2650 acx = (float) (pa[0] - pc[0]);
2651 bcx = (float) (pb[0] - pc[0]);
2652 acy = (float) (pa[1] - pc[1]);
2653 bcy = (float) (pb[1] - pc[1]);
2655 Two_Product(acx, bcy, detleft, detlefttail);
2656 Two_Product(acy, bcx, detright, detrighttail);
2658 Two_Two_Diff(detleft, detlefttail, detright, detrighttail,
2659 B3, B[2], B[1], B[0]);
2662 det = estimate(4, B);
2663 errbound = ccwerrboundB * detsum;
2664 if ((det >= errbound) || (-det >= errbound)) {
2668 Two_Diff_Tail(pa[0], pc[0], acx, acxtail);
2669 Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail);
2670 Two_Diff_Tail(pa[1], pc[1], acy, acytail);
2671 Two_Diff_Tail(pb[1], pc[1], bcy, bcytail);
2673 if ((acxtail == 0.0) && (acytail == 0.0)
2674 && (bcxtail == 0.0) && (bcytail == 0.0)) {
2678 errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det);
2679 det += (acx * bcytail + bcy * acxtail)
2680 - (acy * bcxtail + bcx * acytail);
2681 if ((det >= errbound) || (-det >= errbound)) {
2685 Two_Product(acxtail, bcy, s1, s0);
2686 Two_Product(acytail, bcx, t1, t0);
2687 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
2689 C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1);
2691 Two_Product(acx, bcytail, s1, s0);
2692 Two_Product(acy, bcxtail, t1, t0);
2693 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
2695 C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2);
2697 Two_Product(acxtail, bcytail, s1, s0);
2698 Two_Product(acytail, bcxtail, t1, t0);
2699 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
2701 Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D);
2703 return(D[Dlength - 1]);
2706 float counterclockwise(
struct mesh *m,
struct behavior *b,
2707 vertex pa, vertex pb, vertex pc)
2709 float detleft, detright, det;
2710 float detsum, errbound;
2712 m->counterclockcount++;
2714 detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);
2715 detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);
2716 det = detleft - detright;
2722 if (detleft > 0.0) {
2723 if (detright <= 0.0) {
2726 detsum = detleft + detright;
2728 }
else if (detleft < 0.0) {
2729 if (detright >= 0.0) {
2732 detsum = -detleft - detright;
2738 errbound = ccwerrboundA * detsum;
2739 if ((det >= errbound) || (-det >= errbound)) {
2743 return counterclockwiseadapt(pa, pb, pc, detsum);
2765 float incircleadapt(vertex pa, vertex pb, vertex pc, vertex pd,
float permanent)
2767 float adx, bdx, cdx, ady, bdy, cdy;
2768 float det, errbound;
2770 float bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
2771 float bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
2772 float bc[4], ca[4], ab[4];
2773 float bc3, ca3, ab3;
2774 float axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32];
2775 int axbclen, axxbclen, aybclen, ayybclen, alen;
2776 float bxca[8], bxxca[16], byca[8], byyca[16], bdet[32];
2777 int bxcalen, bxxcalen, bycalen, byycalen, blen;
2778 float cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32];
2779 int cxablen, cxxablen, cyablen, cyyablen, clen;
2782 float fin1[1152], fin2[1152];
2783 float *finnow, *finother, *finswap;
2786 float adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
2787 float adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
2788 float adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
2789 float aa[4], bb[4], cc[4];
2790 float aa3, bb3, cc3;
2795 float temp8[8], temp16a[16], temp16b[16], temp16c[16];
2796 float temp32a[32], temp32b[32], temp48[48], temp64[64];
2797 int temp8len, temp16alen, temp16blen, temp16clen;
2798 int temp32alen, temp32blen, temp48len, temp64len;
2799 float axtbb[8], axtcc[8], aytbb[8], aytcc[8];
2800 int axtbblen, axtcclen, aytbblen, aytcclen;
2801 float bxtaa[8], bxtcc[8], bytaa[8], bytcc[8];
2802 int bxtaalen, bxtcclen, bytaalen, bytcclen;
2803 float cxtaa[8], cxtbb[8], cytaa[8], cytbb[8];
2804 int cxtaalen, cxtbblen, cytaalen, cytbblen;
2805 float axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8];
2806 int axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen;
2807 float axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16];
2808 int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
2809 float axtbctt[8], aytbctt[8], bxtcatt[8];
2810 float bytcatt[8], cxtabtt[8], cytabtt[8];
2811 int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
2812 float abt[8], bct[8], cat[8];
2813 int abtlen, bctlen, catlen;
2814 float abtt[4], bctt[4], catt[4];
2815 int abttlen, bcttlen, cattlen;
2816 float abtt3, bctt3, catt3;
2820 float avirt, bround, around;
2823 float ahi, alo, bhi, blo;
2824 float err1, err2, err3;
2828 adx = (float) (pa[0] - pd[0]);
2829 bdx = (float) (pb[0] - pd[0]);
2830 cdx = (float) (pc[0] - pd[0]);
2831 ady = (float) (pa[1] - pd[1]);
2832 bdy = (float) (pb[1] - pd[1]);
2833 cdy = (float) (pc[1] - pd[1]);
2835 Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
2836 Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
2837 Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
2839 axbclen = scale_expansion_zeroelim(4, bc, adx, axbc);
2840 axxbclen = scale_expansion_zeroelim(axbclen, axbc, adx, axxbc);
2841 aybclen = scale_expansion_zeroelim(4, bc, ady, aybc);
2842 ayybclen = scale_expansion_zeroelim(aybclen, aybc, ady, ayybc);
2843 alen = fast_expansion_sum_zeroelim(axxbclen, axxbc, ayybclen, ayybc, adet);
2845 Two_Product(cdx, ady, cdxady1, cdxady0);
2846 Two_Product(adx, cdy, adxcdy1, adxcdy0);
2847 Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
2849 bxcalen = scale_expansion_zeroelim(4, ca, bdx, bxca);
2850 bxxcalen = scale_expansion_zeroelim(bxcalen, bxca, bdx, bxxca);
2851 bycalen = scale_expansion_zeroelim(4, ca, bdy, byca);
2852 byycalen = scale_expansion_zeroelim(bycalen, byca, bdy, byyca);
2853 blen = fast_expansion_sum_zeroelim(bxxcalen, bxxca, byycalen, byyca, bdet);
2855 Two_Product(adx, bdy, adxbdy1, adxbdy0);
2856 Two_Product(bdx, ady, bdxady1, bdxady0);
2857 Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
2859 cxablen = scale_expansion_zeroelim(4, ab, cdx, cxab);
2860 cxxablen = scale_expansion_zeroelim(cxablen, cxab, cdx, cxxab);
2861 cyablen = scale_expansion_zeroelim(4, ab, cdy, cyab);
2862 cyyablen = scale_expansion_zeroelim(cyablen, cyab, cdy, cyyab);
2863 clen = fast_expansion_sum_zeroelim(cxxablen, cxxab, cyyablen, cyyab, cdet);
2865 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
2866 finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
2868 det = estimate(finlength, fin1);
2869 errbound = iccerrboundB * permanent;
2870 if ((det >= errbound) || (-det >= errbound)) {
2874 Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
2875 Two_Diff_Tail(pa[1], pd[1], ady, adytail);
2876 Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
2877 Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
2878 Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
2879 Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
2880 if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
2881 && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)) {
2885 errbound = iccerrboundC * permanent + resulterrbound * Absolute(det);
2886 det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail)
2887 - (bdy * cdxtail + cdx * bdytail))
2888 + 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx))
2889 + ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail)
2890 - (cdy * adxtail + adx * cdytail))
2891 + 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx))
2892 + ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail)
2893 - (ady * bdxtail + bdx * adytail))
2894 + 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
2895 if ((det >= errbound) || (-det >= errbound)) {
2902 if ((bdxtail != 0.0) || (bdytail != 0.0)
2903 || (cdxtail != 0.0) || (cdytail != 0.0)) {
2904 Square(adx, adxadx1, adxadx0);
2905 Square(ady, adyady1, adyady0);
2906 Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]);
2909 if ((cdxtail != 0.0) || (cdytail != 0.0)
2910 || (adxtail != 0.0) || (adytail != 0.0)) {
2911 Square(bdx, bdxbdx1, bdxbdx0);
2912 Square(bdy, bdybdy1, bdybdy0);
2913 Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]);
2916 if ((adxtail != 0.0) || (adytail != 0.0)
2917 || (bdxtail != 0.0) || (bdytail != 0.0)) {
2918 Square(cdx, cdxcdx1, cdxcdx0);
2919 Square(cdy, cdycdy1, cdycdy0);
2920 Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]);
2924 if (adxtail != 0.0) {
2925 axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc);
2926 temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx,
2929 axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc);
2930 temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b);
2932 axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb);
2933 temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c);
2935 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2936 temp16blen, temp16b, temp32a);
2937 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2938 temp32alen, temp32a, temp48);
2939 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2941 finswap = finnow; finnow = finother; finother = finswap;
2943 if (adytail != 0.0) {
2944 aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc);
2945 temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady,
2948 aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb);
2949 temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b);
2951 aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc);
2952 temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c);
2954 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2955 temp16blen, temp16b, temp32a);
2956 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2957 temp32alen, temp32a, temp48);
2958 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2960 finswap = finnow; finnow = finother; finother = finswap;
2962 if (bdxtail != 0.0) {
2963 bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca);
2964 temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx,
2967 bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa);
2968 temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b);
2970 bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc);
2971 temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c);
2973 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2974 temp16blen, temp16b, temp32a);
2975 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2976 temp32alen, temp32a, temp48);
2977 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2979 finswap = finnow; finnow = finother; finother = finswap;
2981 if (bdytail != 0.0) {
2982 bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca);
2983 temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy,
2986 bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc);
2987 temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b);
2989 bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa);
2990 temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c);
2992 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2993 temp16blen, temp16b, temp32a);
2994 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2995 temp32alen, temp32a, temp48);
2996 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2998 finswap = finnow; finnow = finother; finother = finswap;
3000 if (cdxtail != 0.0) {
3001 cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab);
3002 temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx,
3005 cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb);
3006 temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b);
3008 cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa);
3009 temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c);
3011 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3012 temp16blen, temp16b, temp32a);
3013 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
3014 temp32alen, temp32a, temp48);
3015 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
3017 finswap = finnow; finnow = finother; finother = finswap;
3019 if (cdytail != 0.0) {
3020 cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab);
3021 temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy,
3024 cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa);
3025 temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b);
3027 cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb);
3028 temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c);
3030 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3031 temp16blen, temp16b, temp32a);
3032 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
3033 temp32alen, temp32a, temp48);
3034 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
3036 finswap = finnow; finnow = finother; finother = finswap;
3039 if ((adxtail != 0.0) || (adytail != 0.0)) {
3040 if ((bdxtail != 0.0) || (bdytail != 0.0)
3041 || (cdxtail != 0.0) || (cdytail != 0.0)) {
3042 Two_Product(bdxtail, cdy, ti1, ti0);
3043 Two_Product(bdx, cdytail, tj1, tj0);
3044 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
3047 Two_Product(cdxtail, negate, ti1, ti0);
3049 Two_Product(cdx, negate, tj1, tj0);
3050 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
3052 bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct);
3054 Two_Product(bdxtail, cdytail, ti1, ti0);
3055 Two_Product(cdxtail, bdytail, tj1, tj0);
3056 Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]);
3066 if (adxtail != 0.0) {
3067 temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, adxtail, temp16a);
3068 axtbctlen = scale_expansion_zeroelim(bctlen, bct, adxtail, axtbct);
3069 temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, 2.0 * adx,
3071 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3072 temp32alen, temp32a, temp48);
3073 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
3075 finswap = finnow; finnow = finother; finother = finswap;
3076 if (bdytail != 0.0) {
3077 temp8len = scale_expansion_zeroelim(4, cc, adxtail, temp8);
3078 temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
3080 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
3082 finswap = finnow; finnow = finother; finother = finswap;
3084 if (cdytail != 0.0) {
3085 temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8);
3086 temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
3088 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
3090 finswap = finnow; finnow = finother; finother = finswap;
3093 temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail,
3095 axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt);
3096 temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx,
3098 temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail,
3100 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3101 temp16blen, temp16b, temp32b);
3102 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3103 temp32blen, temp32b, temp64);
3104 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
3106 finswap = finnow; finnow = finother; finother = finswap;
3108 if (adytail != 0.0) {
3109 temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, adytail, temp16a);
3110 aytbctlen = scale_expansion_zeroelim(bctlen, bct, adytail, aytbct);
3111 temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, 2.0 * ady,
3113 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3114 temp32alen, temp32a, temp48);
3115 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
3117 finswap = finnow; finnow = finother; finother = finswap;
3120 temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail,
3122 aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt);
3123 temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady,
3125 temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail,
3127 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3128 temp16blen, temp16b, temp32b);
3129 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3130 temp32blen, temp32b, temp64);
3131 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
3133 finswap = finnow; finnow = finother; finother = finswap;
3136 if ((bdxtail != 0.0) || (bdytail != 0.0)) {
3137 if ((cdxtail != 0.0) || (cdytail != 0.0)
3138 || (adxtail != 0.0) || (adytail != 0.0)) {
3139 Two_Product(cdxtail, ady, ti1, ti0);
3140 Two_Product(cdx, adytail, tj1, tj0);
3141 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
3144 Two_Product(adxtail, negate, ti1, ti0);
3146 Two_Product(adx, negate, tj1, tj0);
3147 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
3149 catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat);
3151 Two_Product(cdxtail, adytail, ti1, ti0);
3152 Two_Product(adxtail, cdytail, tj1, tj0);
3153 Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]);
3163 if (bdxtail != 0.0) {
3164 temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, bdxtail, temp16a);
3165 bxtcatlen = scale_expansion_zeroelim(catlen, cat, bdxtail, bxtcat);
3166 temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, 2.0 * bdx,
3168 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3169 temp32alen, temp32a, temp48);
3170 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
3172 finswap = finnow; finnow = finother; finother = finswap;
3173 if (cdytail != 0.0) {
3174 temp8len = scale_expansion_zeroelim(4, aa, bdxtail, temp8);
3175 temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
3177 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
3179 finswap = finnow; finnow = finother; finother = finswap;
3181 if (adytail != 0.0) {
3182 temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8);
3183 temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
3185 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
3187 finswap = finnow; finnow = finother; finother = finswap;
3190 temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail,
3192 bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt);
3193 temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx,
3195 temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail,
3197 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3198 temp16blen, temp16b, temp32b);
3199 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3200 temp32blen, temp32b, temp64);
3201 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
3203 finswap = finnow; finnow = finother; finother = finswap;
3205 if (bdytail != 0.0) {
3206 temp16alen = scale_expansion_zeroelim(bytcalen, bytca, bdytail, temp16a);
3207 bytcatlen = scale_expansion_zeroelim(catlen, cat, bdytail, bytcat);
3208 temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, 2.0 * bdy,
3210 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3211 temp32alen, temp32a, temp48);
3212 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
3214 finswap = finnow; finnow = finother; finother = finswap;
3217 temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail,
3219 bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt);
3220 temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy,
3222 temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail,
3224 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3225 temp16blen, temp16b, temp32b);
3226 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3227 temp32blen, temp32b, temp64);
3228 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
3230 finswap = finnow; finnow = finother; finother = finswap;
3233 if ((cdxtail != 0.0) || (cdytail != 0.0)) {
3234 if ((adxtail != 0.0) || (adytail != 0.0)
3235 || (bdxtail != 0.0) || (bdytail != 0.0)) {
3236 Two_Product(adxtail, bdy, ti1, ti0);
3237 Two_Product(adx, bdytail, tj1, tj0);
3238 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
3241 Two_Product(bdxtail, negate, ti1, ti0);
3243 Two_Product(bdx, negate, tj1, tj0);
3244 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
3246 abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt);
3248 Two_Product(adxtail, bdytail, ti1, ti0);
3249 Two_Product(bdxtail, adytail, tj1, tj0);
3250 Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]);
3260 if (cdxtail != 0.0) {
3261 temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, cdxtail, temp16a);
3262 cxtabtlen = scale_expansion_zeroelim(abtlen, abt, cdxtail, cxtabt);
3263 temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, 2.0 * cdx,
3265 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3266 temp32alen, temp32a, temp48);
3267 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
3269 finswap = finnow; finnow = finother; finother = finswap;
3270 if (adytail != 0.0) {
3271 temp8len = scale_expansion_zeroelim(4, bb, cdxtail, temp8);
3272 temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
3274 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
3276 finswap = finnow; finnow = finother; finother = finswap;
3278 if (bdytail != 0.0) {
3279 temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8);
3280 temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
3282 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
3284 finswap = finnow; finnow = finother; finother = finswap;
3287 temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail,
3289 cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt);
3290 temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx,
3292 temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail,
3294 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3295 temp16blen, temp16b, temp32b);
3296 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3297 temp32blen, temp32b, temp64);
3298 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
3300 finswap = finnow; finnow = finother; finother = finswap;
3302 if (cdytail != 0.0) {
3303 temp16alen = scale_expansion_zeroelim(cytablen, cytab, cdytail, temp16a);
3304 cytabtlen = scale_expansion_zeroelim(abtlen, abt, cdytail, cytabt);
3305 temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, 2.0 * cdy,
3307 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3308 temp32alen, temp32a, temp48);
3309 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
3311 finswap = finnow; finnow = finother; finother = finswap;
3314 temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail,
3316 cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt);
3317 temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy,
3319 temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail,
3321 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3322 temp16blen, temp16b, temp32b);
3323 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3324 temp32blen, temp32b, temp64);
3325 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
3327 finswap = finnow; finnow = finother; finother = finswap;
3331 return finnow[finlength - 1];
3334 float incircle(
struct mesh *m,
struct behavior *b,
3335 vertex pa, vertex pb, vertex pc, vertex pd)
3337 float adx, bdx, cdx, ady, bdy, cdy;
3338 float bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
3339 float alift, blift, clift;
3341 float permanent, errbound;
3345 adx = pa[0] - pd[0];
3346 bdx = pb[0] - pd[0];
3347 cdx = pc[0] - pd[0];
3348 ady = pa[1] - pd[1];
3349 bdy = pb[1] - pd[1];
3350 cdy = pc[1] - pd[1];
3354 alift = adx * adx + ady * ady;
3358 blift = bdx * bdx + bdy * bdy;
3362 clift = cdx * cdx + cdy * cdy;
3364 det = alift * (bdxcdy - cdxbdy)
3365 + blift * (cdxady - adxcdy)
3366 + clift * (adxbdy - bdxady);
3372 permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift
3373 + (Absolute(cdxady) + Absolute(adxcdy)) * blift
3374 + (Absolute(adxbdy) + Absolute(bdxady)) * clift;
3375 errbound = iccerrboundA * permanent;
3376 if ((det > errbound) || (-det > errbound)) {
3380 return incircleadapt(pa, pb, pc, pd, permanent);
3405 float orient3dadapt(vertex pa, vertex pb, vertex pc, vertex pd,
3406 float aheight,
float bheight,
float cheight,
float dheight,
3409 float adx, bdx, cdx, ady, bdy, cdy, adheight, bdheight, cdheight;
3410 float det, errbound;
3412 float bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
3413 float bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
3414 float bc[4], ca[4], ab[4];
3415 float bc3, ca3, ab3;
3416 float adet[8], bdet[8], cdet[8];
3417 int alen, blen, clen;
3420 float *finnow, *finother, *finswap;
3421 float fin1[192], fin2[192];
3424 float adxtail, bdxtail, cdxtail;
3425 float adytail, bdytail, cdytail;
3426 float adheighttail, bdheighttail, cdheighttail;
3427 float at_blarge, at_clarge;
3428 float bt_clarge, bt_alarge;
3429 float ct_alarge, ct_blarge;
3430 float at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4];
3431 int at_blen, at_clen, bt_clen, bt_alen, ct_alen, ct_blen;
3432 float bdxt_cdy1, cdxt_bdy1, cdxt_ady1;
3433 float adxt_cdy1, adxt_bdy1, bdxt_ady1;
3434 float bdxt_cdy0, cdxt_bdy0, cdxt_ady0;
3435 float adxt_cdy0, adxt_bdy0, bdxt_ady0;
3436 float bdyt_cdx1, cdyt_bdx1, cdyt_adx1;
3437 float adyt_cdx1, adyt_bdx1, bdyt_adx1;
3438 float bdyt_cdx0, cdyt_bdx0, cdyt_adx0;
3439 float adyt_cdx0, adyt_bdx0, bdyt_adx0;
3440 float bct[8], cat[8], abt[8];
3441 int bctlen, catlen, abtlen;
3442 float bdxt_cdyt1, cdxt_bdyt1, cdxt_adyt1;
3443 float adxt_cdyt1, adxt_bdyt1, bdxt_adyt1;
3444 float bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0;
3445 float adxt_cdyt0, adxt_bdyt0, bdxt_adyt0;
3446 float u[4], v[12], w[16];
3448 int vlength, wlength;
3452 float avirt, bround, around;
3455 float ahi, alo, bhi, blo;
3456 float err1, err2, err3;
3460 adx = (float) (pa[0] - pd[0]);
3461 bdx = (float) (pb[0] - pd[0]);
3462 cdx = (float) (pc[0] - pd[0]);
3463 ady = (float) (pa[1] - pd[1]);
3464 bdy = (float) (pb[1] - pd[1]);
3465 cdy = (float) (pc[1] - pd[1]);
3466 adheight = (float) (aheight - dheight);
3467 bdheight = (float) (bheight - dheight);
3468 cdheight = (float) (cheight - dheight);
3470 Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
3471 Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
3472 Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
3474 alen = scale_expansion_zeroelim(4, bc, adheight, adet);
3476 Two_Product(cdx, ady, cdxady1, cdxady0);
3477 Two_Product(adx, cdy, adxcdy1, adxcdy0);
3478 Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
3480 blen = scale_expansion_zeroelim(4, ca, bdheight, bdet);
3482 Two_Product(adx, bdy, adxbdy1, adxbdy0);
3483 Two_Product(bdx, ady, bdxady1, bdxady0);
3484 Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
3486 clen = scale_expansion_zeroelim(4, ab, cdheight, cdet);
3488 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
3489 finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
3491 det = estimate(finlength, fin1);
3492 errbound = o3derrboundB * permanent;
3493 if ((det >= errbound) || (-det >= errbound)) {
3497 Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
3498 Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
3499 Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
3500 Two_Diff_Tail(pa[1], pd[1], ady, adytail);
3501 Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
3502 Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
3503 Two_Diff_Tail(aheight, dheight, adheight, adheighttail);
3504 Two_Diff_Tail(bheight, dheight, bdheight, bdheighttail);
3505 Two_Diff_Tail(cheight, dheight, cdheight, cdheighttail);
3507 if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) &&
3508 (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0) &&
3509 (adheighttail == 0.0) &&
3510 (bdheighttail == 0.0) &&
3511 (cdheighttail == 0.0)) {
3515 errbound = o3derrboundC * permanent + resulterrbound * Absolute(det);
3516 det += (adheight * ((bdx * cdytail + cdy * bdxtail) -
3517 (bdy * cdxtail + cdx * bdytail)) +
3518 adheighttail * (bdx * cdy - bdy * cdx)) +
3519 (bdheight * ((cdx * adytail + ady * cdxtail) -
3520 (cdy * adxtail + adx * cdytail)) +
3521 bdheighttail * (cdx * ady - cdy * adx)) +
3522 (cdheight * ((adx * bdytail + bdy * adxtail) -
3523 (ady * bdxtail + bdx * adytail)) +
3524 cdheighttail * (adx * bdy - ady * bdx));
3525 if ((det >= errbound) || (-det >= errbound)) {
3532 if (adxtail == 0.0) {
3533 if (adytail == 0.0) {
3540 Two_Product(negate, bdx, at_blarge, at_b[0]);
3541 at_b[1] = at_blarge;
3543 Two_Product(adytail, cdx, at_clarge, at_c[0]);
3544 at_c[1] = at_clarge;
3548 if (adytail == 0.0) {
3549 Two_Product(adxtail, bdy, at_blarge, at_b[0]);
3550 at_b[1] = at_blarge;
3553 Two_Product(negate, cdy, at_clarge, at_c[0]);
3554 at_c[1] = at_clarge;
3557 Two_Product(adxtail, bdy, adxt_bdy1, adxt_bdy0);
3558 Two_Product(adytail, bdx, adyt_bdx1, adyt_bdx0);
3559 Two_Two_Diff(adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0,
3560 at_blarge, at_b[2], at_b[1], at_b[0]);
3561 at_b[3] = at_blarge;
3563 Two_Product(adytail, cdx, adyt_cdx1, adyt_cdx0);
3564 Two_Product(adxtail, cdy, adxt_cdy1, adxt_cdy0);
3565 Two_Two_Diff(adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0,
3566 at_clarge, at_c[2], at_c[1], at_c[0]);
3567 at_c[3] = at_clarge;
3571 if (bdxtail == 0.0) {
3572 if (bdytail == 0.0) {
3579 Two_Product(negate, cdx, bt_clarge, bt_c[0]);
3580 bt_c[1] = bt_clarge;
3582 Two_Product(bdytail, adx, bt_alarge, bt_a[0]);
3583 bt_a[1] = bt_alarge;
3587 if (bdytail == 0.0) {
3588 Two_Product(bdxtail, cdy, bt_clarge, bt_c[0]);
3589 bt_c[1] = bt_clarge;
3592 Two_Product(negate, ady, bt_alarge, bt_a[0]);
3593 bt_a[1] = bt_alarge;
3596 Two_Product(bdxtail, cdy, bdxt_cdy1, bdxt_cdy0);
3597 Two_Product(bdytail, cdx, bdyt_cdx1, bdyt_cdx0);
3598 Two_Two_Diff(bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0,
3599 bt_clarge, bt_c[2], bt_c[1], bt_c[0]);
3600 bt_c[3] = bt_clarge;
3602 Two_Product(bdytail, adx, bdyt_adx1, bdyt_adx0);
3603 Two_Product(bdxtail, ady, bdxt_ady1, bdxt_ady0);
3604 Two_Two_Diff(bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0,
3605 bt_alarge, bt_a[2], bt_a[1], bt_a[0]);
3606 bt_a[3] = bt_alarge;
3610 if (cdxtail == 0.0) {
3611 if (cdytail == 0.0) {
3618 Two_Product(negate, adx, ct_alarge, ct_a[0]);
3619 ct_a[1] = ct_alarge;
3621 Two_Product(cdytail, bdx, ct_blarge, ct_b[0]);
3622 ct_b[1] = ct_blarge;
3626 if (cdytail == 0.0) {
3627 Two_Product(cdxtail, ady, ct_alarge, ct_a[0]);
3628 ct_a[1] = ct_alarge;
3631 Two_Product(negate, bdy, ct_blarge, ct_b[0]);
3632 ct_b[1] = ct_blarge;
3635 Two_Product(cdxtail, ady, cdxt_ady1, cdxt_ady0);
3636 Two_Product(cdytail, adx, cdyt_adx1, cdyt_adx0);
3637 Two_Two_Diff(cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0,
3638 ct_alarge, ct_a[2], ct_a[1], ct_a[0]);
3639 ct_a[3] = ct_alarge;
3641 Two_Product(cdytail, bdx, cdyt_bdx1, cdyt_bdx0);
3642 Two_Product(cdxtail, bdy, cdxt_bdy1, cdxt_bdy0);
3643 Two_Two_Diff(cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0,
3644 ct_blarge, ct_b[2], ct_b[1], ct_b[0]);
3645 ct_b[3] = ct_blarge;
3650 bctlen = fast_expansion_sum_zeroelim(bt_clen, bt_c, ct_blen, ct_b, bct);
3651 wlength = scale_expansion_zeroelim(bctlen, bct, adheight, w);
3652 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
3654 finswap = finnow; finnow = finother; finother = finswap;
3656 catlen = fast_expansion_sum_zeroelim(ct_alen, ct_a, at_clen, at_c, cat);
3657 wlength = scale_expansion_zeroelim(catlen, cat, bdheight, w);
3658 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
3660 finswap = finnow; finnow = finother; finother = finswap;
3662 abtlen = fast_expansion_sum_zeroelim(at_blen, at_b, bt_alen, bt_a, abt);
3663 wlength = scale_expansion_zeroelim(abtlen, abt, cdheight, w);
3664 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
3666 finswap = finnow; finnow = finother; finother = finswap;
3668 if (adheighttail != 0.0) {
3669 vlength = scale_expansion_zeroelim(4, bc, adheighttail, v);
3670 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
3672 finswap = finnow; finnow = finother; finother = finswap;
3674 if (bdheighttail != 0.0) {
3675 vlength = scale_expansion_zeroelim(4, ca, bdheighttail, v);
3676 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
3678 finswap = finnow; finnow = finother; finother = finswap;
3680 if (cdheighttail != 0.0) {
3681 vlength = scale_expansion_zeroelim(4, ab, cdheighttail, v);
3682 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
3684 finswap = finnow; finnow = finother; finother = finswap;
3687 if (adxtail != 0.0) {
3688 if (bdytail != 0.0) {
3689 Two_Product(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0);
3690 Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdheight, u3, u[2], u[1], u[0]);
3692 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
3694 finswap = finnow; finnow = finother; finother = finswap;
3695 if (cdheighttail != 0.0) {
3696 Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdheighttail,
3697 u3, u[2], u[1], u[0]);
3699 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
3701 finswap = finnow; finnow = finother; finother = finswap;
3704 if (cdytail != 0.0) {
3706 Two_Product(negate, cdytail, adxt_cdyt1, adxt_cdyt0);
3707 Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdheight, u3, u[2], u[1], u[0]);
3709 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
3711 finswap = finnow; finnow = finother; finother = finswap;
3712 if (bdheighttail != 0.0) {
3713 Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdheighttail,
3714 u3, u[2], u[1], u[0]);
3716 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
3718 finswap = finnow; finnow = finother; finother = finswap;
3722 if (bdxtail != 0.0) {
3723 if (cdytail != 0.0) {
3724 Two_Product(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0);
3725 Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adheight, u3, u[2], u[1], u[0]);
3727 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
3729 finswap = finnow; finnow = finother; finother = finswap;
3730 if (adheighttail != 0.0) {
3731 Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adheighttail,
3732 u3, u[2], u[1], u[0]);
3734 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
3736 finswap = finnow; finnow = finother; finother = finswap;
3739 if (adytail != 0.0) {
3741 Two_Product(negate, adytail, bdxt_adyt1, bdxt_adyt0);
3742 Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdheight, u3, u[2], u[1], u[0]);
3744 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
3746 finswap = finnow; finnow = finother; finother = finswap;
3747 if (cdheighttail != 0.0) {
3748 Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdheighttail,
3749 u3, u[2], u[1], u[0]);
3751 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
3753 finswap = finnow; finnow = finother; finother = finswap;
3757 if (cdxtail != 0.0) {
3758 if (adytail != 0.0) {
3759 Two_Product(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0);
3760 Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdheight, u3, u[2], u[1], u[0]);
3762 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
3764 finswap = finnow; finnow = finother; finother = finswap;
3765 if (bdheighttail != 0.0) {
3766 Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdheighttail,
3767 u3, u[2], u[1], u[0]);
3769 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
3771 finswap = finnow; finnow = finother; finother = finswap;
3774 if (bdytail != 0.0) {
3776 Two_Product(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0);
3777 Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adheight, u3, u[2], u[1], u[0]);
3779 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
3781 finswap = finnow; finnow = finother; finother = finswap;
3782 if (adheighttail != 0.0) {
3783 Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adheighttail,
3784 u3, u[2], u[1], u[0]);
3786 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
3788 finswap = finnow; finnow = finother; finother = finswap;
3793 if (adheighttail != 0.0) {
3794 wlength = scale_expansion_zeroelim(bctlen, bct, adheighttail, w);
3795 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
3797 finswap = finnow; finnow = finother; finother = finswap;
3799 if (bdheighttail != 0.0) {
3800 wlength = scale_expansion_zeroelim(catlen, cat, bdheighttail, w);
3801 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
3803 finswap = finnow; finnow = finother; finother = finswap;
3805 if (cdheighttail != 0.0) {
3806 wlength = scale_expansion_zeroelim(abtlen, abt, cdheighttail, w);
3807 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
3809 finswap = finnow; finnow = finother; finother = finswap;
3812 return finnow[finlength - 1];
3815 float orient3d(
struct mesh *m,
struct behavior *b,
3816 vertex pa, vertex pb, vertex pc, vertex pd,
3817 float aheight,
float bheight,
float cheight,
float dheight)
3819 float adx, bdx, cdx, ady, bdy, cdy, adheight, bdheight, cdheight;
3820 float bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
3822 float permanent, errbound;
3826 adx = pa[0] - pd[0];
3827 bdx = pb[0] - pd[0];
3828 cdx = pc[0] - pd[0];
3829 ady = pa[1] - pd[1];
3830 bdy = pb[1] - pd[1];
3831 cdy = pc[1] - pd[1];
3832 adheight = aheight - dheight;
3833 bdheight = bheight - dheight;
3834 cdheight = cheight - dheight;
3845 det = adheight * (bdxcdy - cdxbdy)
3846 + bdheight * (cdxady - adxcdy)
3847 + cdheight * (adxbdy - bdxady);
3853 permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * Absolute(adheight)
3854 + (Absolute(cdxady) + Absolute(adxcdy)) * Absolute(bdheight)
3855 + (Absolute(adxbdy) + Absolute(bdxady)) * Absolute(cdheight);
3856 errbound = o3derrboundA * permanent;
3857 if ((det > errbound) || (-det > errbound)) {
3861 return orient3dadapt(pa, pb, pc, pd, aheight, bheight, cheight, dheight,
3883 float nonregular(
struct mesh *m,
struct behavior *b,
3884 vertex pa, vertex pb, vertex pc, vertex pd)
3886 if (b->weighted == 0) {
3887 return incircle(m, b, pa, pb, pc, pd);
3888 }
else if (b->weighted == 1) {
3889 return orient3d(m, b, pa, pb, pc, pd,
3890 pa[0] * pa[0] + pa[1] * pa[1] - pa[2],
3891 pb[0] * pb[0] + pb[1] * pb[1] - pb[2],
3892 pc[0] * pc[0] + pc[1] * pc[1] - pc[2],
3893 pd[0] * pd[0] + pd[1] * pd[1] - pd[2]);
3895 return orient3d(m, b, pa, pb, pc, pd, pa[2], pb[2], pc[2], pd[2]);
3913 void findcircumcenter(
struct mesh *m,
struct behavior *b,
3914 vertex torg, vertex tdest, vertex tapex,
3915 vertex circumcenter,
float *xi,
float *eta,
int offcenter)
3917 float xdo, ydo, xao, yao;
3918 float dodist, aodist, dadist;
3920 float dx, dy, dxoff, dyoff;
3922 m->circumcentercount++;
3925 xdo = tdest[0] - torg[0];
3926 ydo = tdest[1] - torg[1];
3927 xao = tapex[0] - torg[0];
3928 yao = tapex[1] - torg[1];
3929 dodist = xdo * xdo + ydo * ydo;
3930 aodist = xao * xao + yao * yao;
3931 dadist = (tdest[0] - tapex[0]) * (tdest[0] - tapex[0]) +
3932 (tdest[1] - tapex[1]) * (tdest[1] - tapex[1]);
3934 denominator = 0.5 / (xdo * yao - xao * ydo);
3939 denominator = 0.5 / counterclockwise(m, b, tdest, tapex, torg);
3941 m->counterclockcount--;
3943 dx = (yao * dodist - ydo * aodist) * denominator;
3944 dy = (xdo * aodist - xao * dodist) * denominator;
3951 if ((dodist < aodist) && (dodist < dadist)) {
3952 if (offcenter && (b->offconstant > 0.0)) {
3954 dxoff = 0.5 * xdo - b->offconstant * ydo;
3955 dyoff = 0.5 * ydo + b->offconstant * xdo;
3958 if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy) {
3963 }
else if (aodist < dadist) {
3964 if (offcenter && (b->offconstant > 0.0)) {
3965 dxoff = 0.5 * xao + b->offconstant * yao;
3966 dyoff = 0.5 * yao - b->offconstant * xao;
3969 if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy) {
3975 if (offcenter && (b->offconstant > 0.0)) {
3976 dxoff = 0.5 * (tapex[0] - tdest[0]) -
3977 b->offconstant * (tapex[1] - tdest[1]);
3978 dyoff = 0.5 * (tapex[1] - tdest[1]) +
3979 b->offconstant * (tapex[0] - tdest[0]);
3982 if (dxoff * dxoff + dyoff * dyoff <
3983 (dx - xdo) * (dx - xdo) + (dy - ydo) * (dy - ydo)) {
3990 circumcenter[0] = torg[0] + dx;
3991 circumcenter[1] = torg[1] + dy;
3998 *xi = (yao * dx - xao * dy) * (2.0 * denominator);
3999 *eta = (xdo * dy - ydo * dx) * (2.0 * denominator);
4012 void triangleinit(
struct mesh *m)
4014 poolzero(&m->vertices);
4015 poolzero(&m->triangles);
4016 poolzero(&m->subsegs);
4018 poolzero(&m->badsubsegs);
4019 poolzero(&m->badtriangles);
4020 poolzero(&m->flipstackers);
4021 poolzero(&m->splaynodes);
4023 m->recenttri.tri = (triangle *) NULL;
4026 m->checksegments = 0;
4027 m->checkquality = 0;
4028 m->incirclecount = m->counterclockcount = m->orient3dcount = 0;
4029 m->hyperbolacount = m->circletopcount = m->circumcentercount = 0;
4045 unsigned long randomnation(
unsigned int choices)
4047 randomseed = (randomseed * 1366l + 150889l) % 714025l;
4048 return randomseed / (714025l / choices + 1);
4069 void makevertexmap(
struct mesh *m,
struct behavior *b)
4071 struct otri triangleloop;
4075 printf(
" Constructing mapping from vertices to triangles.\n");
4077 traversalinit(&m->triangles);
4078 triangleloop.tri = triangletraverse(m);
4079 while (triangleloop.tri != (triangle *) NULL) {
4081 for (triangleloop.orient = 0; triangleloop.orient < 3;
4082 triangleloop.orient++) {
4083 org(triangleloop, triorg);
4084 setvertex2tri(triorg, encode(triangleloop));
4086 triangleloop.tri = triangletraverse(m);
4157 enum locateresult preciselocate(
struct mesh *m,
struct behavior *b,
4158 vertex searchpoint,
struct otri *searchtri,
4159 int stopatsubsegment)
4161 struct otri backtracktri;
4162 struct osub checkedge;
4163 vertex forg, fdest, fapex;
4164 float orgorient, destorient;
4169 if (b->verbose > 2) {
4170 printf(
" Searching for point (%.12g, %.12g).\n",
4171 searchpoint[0], searchpoint[1]);
4174 org(*searchtri, forg);
4175 dest(*searchtri, fdest);
4176 apex(*searchtri, fapex);
4178 if (b->verbose > 2) {
4179 printf(
" At (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
4180 forg[0], forg[1], fdest[0], fdest[1], fapex[0], fapex[1]);
4183 if ((fapex[0] == searchpoint[0]) && (fapex[1] == searchpoint[1])) {
4184 lprevself(*searchtri);
4189 destorient = counterclockwise(m, b, forg, fapex, searchpoint);
4192 orgorient = counterclockwise(m, b, fapex, fdest, searchpoint);
4193 if (destorient > 0.0) {
4194 if (orgorient > 0.0) {
4200 moveleft = (fapex[0] - searchpoint[0]) * (fdest[0] - forg[0]) +
4201 (fapex[1] - searchpoint[1]) * (fdest[1] - forg[1]) > 0.0;
4206 if (orgorient > 0.0) {
4211 if (destorient == 0.0) {
4212 lprevself(*searchtri);
4215 if (orgorient == 0.0) {
4216 lnextself(*searchtri);
4227 lprev(*searchtri, backtracktri);
4230 lnext(*searchtri, backtracktri);
4233 sym(backtracktri, *searchtri);
4235 if (m->checksegments && stopatsubsegment) {
4237 tspivot(backtracktri, checkedge);
4238 if (checkedge.ss != m->dummysub) {
4240 otricopy(backtracktri, *searchtri);
4245 if (searchtri->tri == m->dummytri) {
4247 otricopy(backtracktri, *searchtri);
4251 apex(*searchtri, fapex);
4291 enum locateresult locate(
struct mesh *m,
struct behavior *b,
4292 vertex searchpoint,
struct otri *searchtri)
4296 struct otri sampletri;
4298 unsigned long alignptr;
4299 float searchdist, dist;
4301 long samplesperblock, totalsamplesleft, samplesleft;
4302 long population, totalpopulation;
4305 if (b->verbose > 2) {
4306 printf(
" Randomly sampling for a triangle near point (%.12g, %.12g).\n",
4307 searchpoint[0], searchpoint[1]);
4311 org(*searchtri, torg);
4312 searchdist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0]) +
4313 (searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
4314 if (b->verbose > 2) {
4315 printf(
" Boundary triangle has origin (%.12g, %.12g).\n",
4321 if (m->recenttri.tri != (triangle *) NULL) {
4322 if (!deadtri(m->recenttri.tri)) {
4323 org(m->recenttri, torg);
4324 if ((torg[0] == searchpoint[0]) && (torg[1] == searchpoint[1])) {
4325 otricopy(m->recenttri, *searchtri);
4328 dist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0]) +
4329 (searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
4330 if (dist < searchdist) {
4331 otricopy(m->recenttri, *searchtri);
4333 if (b->verbose > 2) {
4334 printf(
" Choosing recent triangle with origin (%.12g, %.12g).\n",
4345 while (SAMPLEFACTOR * m->samples * m->samples * m->samples <
4346 m->triangles.items) {
4354 samplesperblock = (m->samples * TRIPERBLOCK - 1) / m->triangles.maxitems + 1;
4357 samplesleft = (m->samples * m->triangles.itemsfirstblock - 1) /
4358 m->triangles.maxitems + 1;
4359 totalsamplesleft = m->samples;
4360 population = m->triangles.itemsfirstblock;
4361 totalpopulation = m->triangles.maxitems;
4362 sampleblock = m->triangles.firstblock;
4363 sampletri.orient = 0;
4364 while (totalsamplesleft > 0) {
4366 if (population > totalpopulation) {
4367 population = totalpopulation;
4370 alignptr = (
unsigned long) (sampleblock + 1);
4371 firsttri = (
char *) (alignptr +
4372 (
unsigned long) m->triangles.alignbytes -
4374 (
unsigned long) m->triangles.alignbytes));
4378 sampletri.tri = (triangle *) (firsttri +
4379 (randomnation((
unsigned int) population) *
4380 m->triangles.itembytes));
4381 if (!deadtri(sampletri.tri)) {
4382 org(sampletri, torg);
4383 dist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0]) +
4384 (searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
4385 if (dist < searchdist) {
4386 otricopy(sampletri, *searchtri);
4388 if (b->verbose > 2) {
4389 printf(
" Choosing triangle with origin (%.12g, %.12g).\n",
4397 }
while ((samplesleft > 0) && (totalsamplesleft > 0));
4399 if (totalsamplesleft > 0) {
4400 sampleblock = (
int **) *sampleblock;
4401 samplesleft = samplesperblock;
4402 totalpopulation -= population;
4403 population = TRIPERBLOCK;
4408 org(*searchtri, torg);
4409 dest(*searchtri, tdest);
4411 if ((torg[0] == searchpoint[0]) && (torg[1] == searchpoint[1])) {
4414 if ((tdest[0] == searchpoint[0]) && (tdest[1] == searchpoint[1])) {
4415 lnextself(*searchtri);
4419 ahead = counterclockwise(m, b, torg, tdest, searchpoint);
4423 symself(*searchtri);
4424 }
else if (ahead == 0.0) {
4426 if (((torg[0] < searchpoint[0]) == (searchpoint[0] < tdest[0])) &&
4427 ((torg[1] < searchpoint[1]) == (searchpoint[1] < tdest[1]))) {
4431 return preciselocate(m, b, searchpoint, searchtri, 0);
4453 void insertsubseg(
struct mesh *m,
struct behavior *b,
struct otri *tri,
4456 struct otri oppotri;
4457 struct osub newsubseg;
4458 vertex triorg, tridest;
4463 dest(*tri, tridest);
4465 if (vertexmark(triorg) == 0) {
4466 setvertexmark(triorg, subsegmark);
4468 if (vertexmark(tridest) == 0) {
4469 setvertexmark(tridest, subsegmark);
4472 tspivot(*tri, newsubseg);
4473 if (newsubseg.ss == m->dummysub) {
4475 makesubseg(m, &newsubseg);
4476 setsorg(newsubseg, tridest);
4477 setsdest(newsubseg, triorg);
4478 setsegorg(newsubseg, tridest);
4479 setsegdest(newsubseg, triorg);
4484 tsbond(*tri, newsubseg);
4486 ssymself(newsubseg);
4487 tsbond(oppotri, newsubseg);
4488 setmark(newsubseg, subsegmark);
4489 if (b->verbose > 2) {
4490 printf(
" Inserting new ");
4491 printsubseg(m, b, &newsubseg);
4494 if (mark(newsubseg) == 0) {
4495 setmark(newsubseg, subsegmark);
4548 void flip(
struct mesh *m,
struct behavior *b,
struct otri *flipedge)
4550 struct otri botleft, botright;
4551 struct otri topleft, topright;
4553 struct otri botlcasing, botrcasing;
4554 struct otri toplcasing, toprcasing;
4555 struct osub botlsubseg, botrsubseg;
4556 struct osub toplsubseg, toprsubseg;
4557 vertex leftvertex, rightvertex, botvertex;
4563 org(*flipedge, rightvertex);
4564 dest(*flipedge, leftvertex);
4565 apex(*flipedge, botvertex);
4566 sym(*flipedge, top);
4567 apex(top, farvertex);
4570 lprev(top, topleft);
4571 sym(topleft, toplcasing);
4572 lnext(top, topright);
4573 sym(topright, toprcasing);
4574 lnext(*flipedge, botleft);
4575 sym(botleft, botlcasing);
4576 lprev(*flipedge, botright);
4577 sym(botright, botrcasing);
4579 bond(topleft, botlcasing);
4580 bond(botleft, botrcasing);
4581 bond(botright, toprcasing);
4582 bond(topright, toplcasing);
4584 if (m->checksegments) {
4586 tspivot(topleft, toplsubseg);
4587 tspivot(botleft, botlsubseg);
4588 tspivot(botright, botrsubseg);
4589 tspivot(topright, toprsubseg);
4590 if (toplsubseg.ss == m->dummysub) {
4591 tsdissolve(topright);
4593 tsbond(topright, toplsubseg);
4595 if (botlsubseg.ss == m->dummysub) {
4596 tsdissolve(topleft);
4598 tsbond(topleft, botlsubseg);
4600 if (botrsubseg.ss == m->dummysub) {
4601 tsdissolve(botleft);
4603 tsbond(botleft, botrsubseg);
4605 if (toprsubseg.ss == m->dummysub) {
4606 tsdissolve(botright);
4608 tsbond(botright, toprsubseg);
4613 setorg(*flipedge, farvertex);
4614 setdest(*flipedge, botvertex);
4615 setapex(*flipedge, rightvertex);
4616 setorg(top, botvertex);
4617 setdest(top, farvertex);
4618 setapex(top, leftvertex);
4619 if (b->verbose > 2) {
4620 printf(
" Edge flip results in left ");
4621 printtriangle(m, b, &top);
4622 printf(
" and right ");
4623 printtriangle(m, b, flipedge);
4660 void unflip(
struct mesh *m,
struct behavior *b,
struct otri *flipedge)
4662 struct otri botleft, botright;
4663 struct otri topleft, topright;
4665 struct otri botlcasing, botrcasing;
4666 struct otri toplcasing, toprcasing;
4667 struct osub botlsubseg, botrsubseg;
4668 struct osub toplsubseg, toprsubseg;
4669 vertex leftvertex, rightvertex, botvertex;
4675 org(*flipedge, rightvertex);
4676 dest(*flipedge, leftvertex);
4677 apex(*flipedge, botvertex);
4678 sym(*flipedge, top);
4679 apex(top, farvertex);
4682 lprev(top, topleft);
4683 sym(topleft, toplcasing);
4684 lnext(top, topright);
4685 sym(topright, toprcasing);
4686 lnext(*flipedge, botleft);
4687 sym(botleft, botlcasing);
4688 lprev(*flipedge, botright);
4689 sym(botright, botrcasing);
4691 bond(topleft, toprcasing);
4692 bond(botleft, toplcasing);
4693 bond(botright, botlcasing);
4694 bond(topright, botrcasing);
4696 if (m->checksegments) {
4698 tspivot(topleft, toplsubseg);
4699 tspivot(botleft, botlsubseg);
4700 tspivot(botright, botrsubseg);
4701 tspivot(topright, toprsubseg);
4702 if (toplsubseg.ss == m->dummysub) {
4703 tsdissolve(botleft);
4705 tsbond(botleft, toplsubseg);
4707 if (botlsubseg.ss == m->dummysub) {
4708 tsdissolve(botright);
4710 tsbond(botright, botlsubseg);
4712 if (botrsubseg.ss == m->dummysub) {
4713 tsdissolve(topright);
4715 tsbond(topright, botrsubseg);
4717 if (toprsubseg.ss == m->dummysub) {
4718 tsdissolve(topleft);
4720 tsbond(topleft, toprsubseg);
4725 setorg(*flipedge, botvertex);
4726 setdest(*flipedge, farvertex);
4727 setapex(*flipedge, leftvertex);
4728 setorg(top, farvertex);
4729 setdest(top, botvertex);
4730 setapex(top, rightvertex);
4731 if (b->verbose > 2) {
4732 printf(
" Edge unflip results in left ");
4733 printtriangle(m, b, flipedge);
4734 printf(
" and right ");
4735 printtriangle(m, b, &top);
4786 enum insertvertexresult insertvertex(
struct mesh *m,
struct behavior *b,
4787 vertex newvertex,
struct otri *searchtri,
4788 struct osub *splitseg,
4789 int segmentflaws,
int triflaws)
4793 struct otri botleft, botright;
4794 struct otri topleft, topright;
4795 struct otri newbotleft, newbotright;
4796 struct otri newtopright;
4797 struct otri botlcasing, botrcasing;
4798 struct otri toplcasing, toprcasing;
4799 struct otri testtri;
4800 struct osub botlsubseg, botrsubseg;
4801 struct osub toplsubseg, toprsubseg;
4802 struct osub brokensubseg;
4803 struct osub checksubseg;
4804 struct osub rightsubseg;
4805 struct osub newsubseg;
4806 struct badsubseg *encroached;
4807 struct flipstacker *newflip;
4809 vertex leftvertex, rightvertex, botvertex, topvertex, farvertex;
4810 vertex segmentorg, segmentdest;
4813 enum insertvertexresult success;
4814 enum locateresult intersect;
4822 if (b->verbose > 1) {
4823 printf(
" Inserting (%.12g, %.12g).\n", newvertex[0], newvertex[1]);
4826 if (splitseg == (
struct osub *) NULL) {
4829 if (searchtri->tri == m->dummytri) {
4831 horiz.tri = m->dummytri;
4835 intersect = locate(m, b, newvertex, &horiz);
4838 otricopy(*searchtri, horiz);
4839 intersect = preciselocate(m, b, newvertex, &horiz, 1);
4844 otricopy(*searchtri, horiz);
4848 if (intersect == ONVERTEX) {
4851 otricopy(horiz, *searchtri);
4852 otricopy(horiz, m->recenttri);
4853 return DUPLICATEVERTEX;
4855 if ((intersect == ONEDGE) || (intersect == OUTSIDE)) {
4857 if (m->checksegments && (splitseg == (
struct osub *) NULL)) {
4859 tspivot(horiz, brokensubseg);
4860 if (brokensubseg.ss != m->dummysub) {
4863 enq = b->nobisect != 2;
4864 if (enq && (b->nobisect == 1)) {
4867 sym(horiz, testtri);
4868 enq = testtri.tri != m->dummytri;
4872 encroached = (
struct badsubseg *) poolalloc(&m->badsubsegs);
4873 encroached->encsubseg = sencode(brokensubseg);
4874 sorg(brokensubseg, encroached->subsegorg);
4875 sdest(brokensubseg, encroached->subsegdest);
4876 if (b->verbose > 2) {
4878 " Queueing encroached subsegment (%.12g, %.12g) (%.12g, %.12g).\n",
4879 encroached->subsegorg[0], encroached->subsegorg[1],
4880 encroached->subsegdest[0], encroached->subsegdest[1]);
4886 otricopy(horiz, *searchtri);
4887 otricopy(horiz, m->recenttri);
4888 return VIOLATINGVERTEX;
4894 lprev(horiz, botright);
4895 sym(botright, botrcasing);
4896 sym(horiz, topright);
4898 mirrorflag = topright.tri != m->dummytri;
4900 lnextself(topright);
4901 sym(topright, toprcasing);
4902 maketriangle(m, b, &newtopright);
4907 maketriangle(m, b, &newbotright);
4910 org(horiz, rightvertex);
4911 dest(horiz, leftvertex);
4912 apex(horiz, botvertex);
4913 setorg(newbotright, botvertex);
4914 setdest(newbotright, rightvertex);
4915 setapex(newbotright, newvertex);
4916 setorg(horiz, newvertex);
4917 for (i = 0; i < m->eextras; i++) {
4919 setelemattribute(newbotright, i, elemattribute(botright, i));
4923 setareabound(newbotright, areabound(botright));
4926 dest(topright, topvertex);
4927 setorg(newtopright, rightvertex);
4928 setdest(newtopright, topvertex);
4929 setapex(newtopright, newvertex);
4930 setorg(topright, newvertex);
4931 for (i = 0; i < m->eextras; i++) {
4933 setelemattribute(newtopright, i, elemattribute(topright, i));
4937 setareabound(newtopright, areabound(topright));
4943 if (m->checksegments) {
4944 tspivot(botright, botrsubseg);
4945 if (botrsubseg.ss != m->dummysub) {
4946 tsdissolve(botright);
4947 tsbond(newbotright, botrsubseg);
4950 tspivot(topright, toprsubseg);
4951 if (toprsubseg.ss != m->dummysub) {
4952 tsdissolve(topright);
4953 tsbond(newtopright, toprsubseg);
4959 bond(newbotright, botrcasing);
4960 lprevself(newbotright);
4961 bond(newbotright, botright);
4962 lprevself(newbotright);
4964 bond(newtopright, toprcasing);
4965 lnextself(newtopright);
4966 bond(newtopright, topright);
4967 lnextself(newtopright);
4968 bond(newtopright, newbotright);
4971 if (splitseg != (
struct osub *) NULL) {
4973 setsdest(*splitseg, newvertex);
4974 segorg(*splitseg, segmentorg);
4975 segdest(*splitseg, segmentdest);
4976 ssymself(*splitseg);
4977 spivot(*splitseg, rightsubseg);
4978 insertsubseg(m, b, &newbotright, mark(*splitseg));
4979 tspivot(newbotright, newsubseg);
4980 setsegorg(newsubseg, segmentorg);
4981 setsegdest(newsubseg, segmentdest);
4982 sbond(*splitseg, newsubseg);
4983 ssymself(newsubseg);
4984 sbond(newsubseg, rightsubseg);
4985 ssymself(*splitseg);
4988 if (vertexmark(newvertex) == 0) {
4989 setvertexmark(newvertex, mark(*splitseg));
4993 if (m->checkquality) {
4994 poolrestart(&m->flipstackers);
4995 m->lastflip = (
struct flipstacker *) poolalloc(&m->flipstackers);
4996 m->lastflip->flippedtri = encode(horiz);
4997 m->lastflip->prevflip = (
struct flipstacker *) &insertvertex;
4999 if (b->verbose > 2) {
5000 printf(
" Updating bottom left ");
5001 printtriangle(m, b, &botright);
5003 printf(
" Updating top left ");
5004 printtriangle(m, b, &topright);
5005 printf(
" Creating top right ");
5006 printtriangle(m, b, &newtopright);
5008 printf(
" Creating bottom right ");
5009 printtriangle(m, b, &newbotright);
5017 lnext(horiz, botleft);
5018 lprev(horiz, botright);
5019 sym(botleft, botlcasing);
5020 sym(botright, botrcasing);
5021 maketriangle(m, b, &newbotleft);
5022 maketriangle(m, b, &newbotright);
5025 org(horiz, rightvertex);
5026 dest(horiz, leftvertex);
5027 apex(horiz, botvertex);
5028 setorg(newbotleft, leftvertex);
5029 setdest(newbotleft, botvertex);
5030 setapex(newbotleft, newvertex);
5031 setorg(newbotright, botvertex);
5032 setdest(newbotright, rightvertex);
5033 setapex(newbotright, newvertex);
5034 setapex(horiz, newvertex);
5035 for (i = 0; i < m->eextras; i++) {
5037 attrib = elemattribute(horiz, i);
5038 setelemattribute(newbotleft, i, attrib);
5039 setelemattribute(newbotright, i, attrib);
5043 area = areabound(horiz);
5044 setareabound(newbotleft, area);
5045 setareabound(newbotright, area);
5050 if (m->checksegments) {
5051 tspivot(botleft, botlsubseg);
5052 if (botlsubseg.ss != m->dummysub) {
5053 tsdissolve(botleft);
5054 tsbond(newbotleft, botlsubseg);
5056 tspivot(botright, botrsubseg);
5057 if (botrsubseg.ss != m->dummysub) {
5058 tsdissolve(botright);
5059 tsbond(newbotright, botrsubseg);
5064 bond(newbotleft, botlcasing);
5065 bond(newbotright, botrcasing);
5066 lnextself(newbotleft);
5067 lprevself(newbotright);
5068 bond(newbotleft, newbotright);
5069 lnextself(newbotleft);
5070 bond(botleft, newbotleft);
5071 lprevself(newbotright);
5072 bond(botright, newbotright);
5074 if (m->checkquality) {
5075 poolrestart(&m->flipstackers);
5076 m->lastflip = (
struct flipstacker *) poolalloc(&m->flipstackers);
5077 m->lastflip->flippedtri = encode(horiz);
5078 m->lastflip->prevflip = (
struct flipstacker *) NULL;
5080 if (b->verbose > 2) {
5081 printf(
" Updating top ");
5082 printtriangle(m, b, &horiz);
5083 printf(
" Creating left ");
5084 printtriangle(m, b, &newbotleft);
5085 printf(
" Creating right ");
5086 printtriangle(m, b, &newbotright);
5092 success = SUCCESSFULVERTEX;
5098 rightvertex = first;
5099 dest(horiz, leftvertex);
5105 if (m->checksegments) {
5107 tspivot(horiz, checksubseg);
5108 if (checksubseg.ss != m->dummysub) {
5117 if (top.tri == m->dummytri) {
5122 apex(top, farvertex);
5128 if ((leftvertex == m->infvertex1) || (leftvertex == m->infvertex2) ||
5129 (leftvertex == m->infvertex3)) {
5134 doflip = counterclockwise(m, b, newvertex, rightvertex, farvertex)
5136 }
else if ((rightvertex == m->infvertex1) ||
5137 (rightvertex == m->infvertex2) ||
5138 (rightvertex == m->infvertex3)) {
5143 doflip = counterclockwise(m, b, farvertex, leftvertex, newvertex)
5145 }
else if ((farvertex == m->infvertex1) ||
5146 (farvertex == m->infvertex2) ||
5147 (farvertex == m->infvertex3)) {
5153 doflip = incircle(m, b, leftvertex, newvertex, rightvertex,
5160 lprev(top, topleft);
5161 sym(topleft, toplcasing);
5162 lnext(top, topright);
5163 sym(topright, toprcasing);
5164 lnext(horiz, botleft);
5165 sym(botleft, botlcasing);
5166 lprev(horiz, botright);
5167 sym(botright, botrcasing);
5169 bond(topleft, botlcasing);
5170 bond(botleft, botrcasing);
5171 bond(botright, toprcasing);
5172 bond(topright, toplcasing);
5173 if (m->checksegments) {
5175 tspivot(topleft, toplsubseg);
5176 tspivot(botleft, botlsubseg);
5177 tspivot(botright, botrsubseg);
5178 tspivot(topright, toprsubseg);
5179 if (toplsubseg.ss == m->dummysub) {
5180 tsdissolve(topright);
5182 tsbond(topright, toplsubseg);
5184 if (botlsubseg.ss == m->dummysub) {
5185 tsdissolve(topleft);
5187 tsbond(topleft, botlsubseg);
5189 if (botrsubseg.ss == m->dummysub) {
5190 tsdissolve(botleft);
5192 tsbond(botleft, botrsubseg);
5194 if (toprsubseg.ss == m->dummysub) {
5195 tsdissolve(botright);
5197 tsbond(botright, toprsubseg);
5201 setorg(horiz, farvertex);
5202 setdest(horiz, newvertex);
5203 setapex(horiz, rightvertex);
5204 setorg(top, newvertex);
5205 setdest(top, farvertex);
5206 setapex(top, leftvertex);
5207 for (i = 0; i < m->eextras; i++) {
5209 attrib = 0.5 * (elemattribute(top, i) + elemattribute(horiz, i));
5210 setelemattribute(top, i, attrib);
5211 setelemattribute(horiz, i, attrib);
5214 if ((areabound(top) <= 0.0) || (areabound(horiz) <= 0.0)) {
5220 area = 0.5 * (areabound(top) + areabound(horiz));
5222 setareabound(top, area);
5223 setareabound(horiz, area);
5226 if (m->checkquality) {
5227 newflip = (
struct flipstacker *) poolalloc(&m->flipstackers);
5228 newflip->flippedtri = encode(horiz);
5229 newflip->prevflip = m->lastflip;
5230 m->lastflip = newflip;
5232 if (b->verbose > 2) {
5233 printf(
" Edge flip results in left ");
5235 printtriangle(m, b, &topleft);
5236 printf(
" and right ");
5237 printtriangle(m, b, &horiz);
5243 leftvertex = farvertex;
5251 sym(horiz, testtri);
5255 if ((leftvertex == first) || (testtri.tri == m->dummytri)) {
5257 lnext(horiz, *searchtri);
5258 lnext(horiz, m->recenttri);
5262 lnext(testtri, horiz);
5263 rightvertex = leftvertex;
5264 dest(horiz, leftvertex);
5333 void triangulatepolygon(
struct mesh *m,
struct behavior *b,
5334 struct otri *firstedge,
struct otri *lastedge,
5335 int edgecount,
int doflip,
int triflaws)
5337 struct otri testtri;
5338 struct otri besttri;
5339 struct otri tempedge;
5340 vertex leftbasevertex, rightbasevertex;
5348 apex(*lastedge, leftbasevertex);
5349 dest(*firstedge, rightbasevertex);
5350 if (b->verbose > 2) {
5351 printf(
" Triangulating interior polygon at edge\n");
5352 printf(
" (%.12g, %.12g) (%.12g, %.12g)\n", leftbasevertex[0],
5353 leftbasevertex[1], rightbasevertex[0], rightbasevertex[1]);
5356 onext(*firstedge, besttri);
5357 dest(besttri, bestvertex);
5358 otricopy(besttri, testtri);
5360 for (i = 2; i <= edgecount - 2; i++) {
5362 dest(testtri, testvertex);
5364 if (incircle(m, b, leftbasevertex, rightbasevertex, bestvertex,
5365 testvertex) > 0.0) {
5366 otricopy(testtri, besttri);
5367 bestvertex = testvertex;
5371 if (b->verbose > 2) {
5372 printf(
" Connecting edge to (%.12g, %.12g)\n", bestvertex[0],
5375 if (bestnumber > 1) {
5377 oprev(besttri, tempedge);
5378 triangulatepolygon(m, b, firstedge, &tempedge, bestnumber + 1, 1,
5381 if (bestnumber < edgecount - 2) {
5383 sym(besttri, tempedge);
5384 triangulatepolygon(m, b, &besttri, lastedge, edgecount - bestnumber, 1,
5387 sym(tempedge, besttri);
5391 flip(m, b, &besttri);
5394 otricopy(besttri, *lastedge);
5446 void vertexsort(vertex *sortarray,
int arraysize)
5450 float pivotx, pivoty;
5453 if (arraysize == 2) {
5455 if ((sortarray[0][0] > sortarray[1][0]) ||
5456 ((sortarray[0][0] == sortarray[1][0]) &&
5457 (sortarray[0][1] > sortarray[1][1]))) {
5458 temp = sortarray[1];
5459 sortarray[1] = sortarray[0];
5460 sortarray[0] = temp;
5465 pivot = (int) randomnation((
unsigned int) arraysize);
5466 pivotx = sortarray[pivot][0];
5467 pivoty = sortarray[pivot][1];
5471 while (left < right) {
5475 }
while ((left <= right) && ((sortarray[left][0] < pivotx) ||
5476 ((sortarray[left][0] == pivotx) &&
5477 (sortarray[left][1] < pivoty))));
5481 }
while ((left <= right) && ((sortarray[right][0] > pivotx) ||
5482 ((sortarray[right][0] == pivotx) &&
5483 (sortarray[right][1] > pivoty))));
5486 temp = sortarray[left];
5487 sortarray[left] = sortarray[right];
5488 sortarray[right] = temp;
5493 vertexsort(sortarray, left);
5495 if (right < arraysize - 2) {
5497 vertexsort(&sortarray[right + 1], arraysize - right - 1);
5513 void vertexmedian(vertex *sortarray,
int arraysize,
int median,
int axis)
5517 float pivot1, pivot2;
5520 if (arraysize == 2) {
5522 if ((sortarray[0][axis] > sortarray[1][axis]) ||
5523 ((sortarray[0][axis] == sortarray[1][axis]) &&
5524 (sortarray[0][1 - axis] > sortarray[1][1 - axis]))) {
5525 temp = sortarray[1];
5526 sortarray[1] = sortarray[0];
5527 sortarray[0] = temp;
5532 pivot = (int) randomnation((
unsigned int) arraysize);
5533 pivot1 = sortarray[pivot][axis];
5534 pivot2 = sortarray[pivot][1 - axis];
5538 while (left < right) {
5542 }
while ((left <= right) && ((sortarray[left][axis] < pivot1) ||
5543 ((sortarray[left][axis] == pivot1) &&
5544 (sortarray[left][1 - axis] < pivot2))));
5548 }
while ((left <= right) && ((sortarray[right][axis] > pivot1) ||
5549 ((sortarray[right][axis] == pivot1) &&
5550 (sortarray[right][1 - axis] > pivot2))));
5553 temp = sortarray[left];
5554 sortarray[left] = sortarray[right];
5555 sortarray[right] = temp;
5560 if (left > median) {
5562 vertexmedian(sortarray, left, median, axis);
5564 if (right < median - 1) {
5566 vertexmedian(&sortarray[right + 1], arraysize - right - 1,
5567 median - right - 1, axis);
5582 void alternateaxes(vertex *sortarray,
int arraysize,
int axis)
5586 divider = arraysize >> 1;
5587 if (arraysize <= 3) {
5593 vertexmedian(sortarray, arraysize, divider, axis);
5595 if (arraysize - divider >= 2) {
5597 alternateaxes(sortarray, divider, 1 - axis);
5599 alternateaxes(&sortarray[divider], arraysize - divider, 1 - axis);
5638 void mergehulls(
struct mesh *m,
struct behavior *b,
struct otri *farleft,
5639 struct otri *innerleft,
struct otri *innerright,
5640 struct otri *farright,
int axis)
5642 struct otri leftcand, rightcand;
5643 struct otri baseedge;
5644 struct otri nextedge;
5645 struct otri sidecasing, topcasing, outercasing;
5646 struct otri checkedge;
5647 vertex innerleftdest;
5648 vertex innerrightorg;
5649 vertex innerleftapex, innerrightapex;
5650 vertex farleftpt, farrightpt;
5651 vertex farleftapex, farrightapex;
5652 vertex lowerleft, lowerright;
5653 vertex upperleft, upperright;
5658 int leftfinished, rightfinished;
5661 dest(*innerleft, innerleftdest);
5662 apex(*innerleft, innerleftapex);
5663 org(*innerright, innerrightorg);
5664 apex(*innerright, innerrightapex);
5666 if (b->dwyer && (axis == 1)) {
5667 org(*farleft, farleftpt);
5668 apex(*farleft, farleftapex);
5669 dest(*farright, farrightpt);
5670 apex(*farright, farrightapex);
5674 while (farleftapex[1] < farleftpt[1]) {
5675 lnextself(*farleft);
5677 farleftpt = farleftapex;
5678 apex(*farleft, farleftapex);
5680 sym(*innerleft, checkedge);
5681 apex(checkedge, checkvertex);
5682 while (checkvertex[1] > innerleftdest[1]) {
5683 lnext(checkedge, *innerleft);
5684 innerleftapex = innerleftdest;
5685 innerleftdest = checkvertex;
5686 sym(*innerleft, checkedge);
5687 apex(checkedge, checkvertex);
5689 while (innerrightapex[1] < innerrightorg[1]) {
5690 lnextself(*innerright);
5691 symself(*innerright);
5692 innerrightorg = innerrightapex;
5693 apex(*innerright, innerrightapex);
5695 sym(*farright, checkedge);
5696 apex(checkedge, checkvertex);
5697 while (checkvertex[1] > farrightpt[1]) {
5698 lnext(checkedge, *farright);
5699 farrightapex = farrightpt;
5700 farrightpt = checkvertex;
5701 sym(*farright, checkedge);
5702 apex(checkedge, checkvertex);
5709 if (counterclockwise(m, b, innerleftdest, innerleftapex, innerrightorg) >
5711 lprevself(*innerleft);
5712 symself(*innerleft);
5713 innerleftdest = innerleftapex;
5714 apex(*innerleft, innerleftapex);
5718 if (counterclockwise(m, b, innerrightapex, innerrightorg, innerleftdest) >
5720 lnextself(*innerright);
5721 symself(*innerright);
5722 innerrightorg = innerrightapex;
5723 apex(*innerright, innerrightapex);
5726 }
while (changemade);
5728 sym(*innerleft, leftcand);
5729 sym(*innerright, rightcand);
5731 maketriangle(m, b, &baseedge);
5733 bond(baseedge, *innerleft);
5734 lnextself(baseedge);
5735 bond(baseedge, *innerright);
5736 lnextself(baseedge);
5737 setorg(baseedge, innerrightorg);
5738 setdest(baseedge, innerleftdest);
5740 if (b->verbose > 2) {
5741 printf(
" Creating base bounding ");
5742 printtriangle(m, b, &baseedge);
5745 org(*farleft, farleftpt);
5746 if (innerleftdest == farleftpt) {
5747 lnext(baseedge, *farleft);
5749 dest(*farright, farrightpt);
5750 if (innerrightorg == farrightpt) {
5751 lprev(baseedge, *farright);
5754 lowerleft = innerleftdest;
5755 lowerright = innerrightorg;
5757 apex(leftcand, upperleft);
5758 apex(rightcand, upperright);
5765 leftfinished = counterclockwise(m, b, upperleft, lowerleft, lowerright) <=
5767 rightfinished = counterclockwise(m, b, upperright, lowerleft, lowerright)
5769 if (leftfinished && rightfinished) {
5771 maketriangle(m, b, &nextedge);
5772 setorg(nextedge, lowerleft);
5773 setdest(nextedge, lowerright);
5776 bond(nextedge, baseedge);
5777 lnextself(nextedge);
5778 bond(nextedge, rightcand);
5779 lnextself(nextedge);
5780 bond(nextedge, leftcand);
5781 if (b->verbose > 2) {
5782 printf(
" Creating top bounding ");
5783 printtriangle(m, b, &nextedge);
5786 if (b->dwyer && (axis == 1)) {
5787 org(*farleft, farleftpt);
5788 apex(*farleft, farleftapex);
5789 dest(*farright, farrightpt);
5790 apex(*farright, farrightapex);
5791 sym(*farleft, checkedge);
5792 apex(checkedge, checkvertex);
5796 while (checkvertex[0] < farleftpt[0]) {
5797 lprev(checkedge, *farleft);
5798 farleftapex = farleftpt;
5799 farleftpt = checkvertex;
5800 sym(*farleft, checkedge);
5801 apex(checkedge, checkvertex);
5803 while (farrightapex[0] > farrightpt[0]) {
5804 lprevself(*farright);
5806 farrightpt = farrightapex;
5807 apex(*farright, farrightapex);
5813 if (!leftfinished) {
5815 lprev(leftcand, nextedge);
5817 apex(nextedge, nextapex);
5820 if (nextapex != (vertex) NULL) {
5822 badedge = incircle(m, b, lowerleft, lowerright, upperleft, nextapex) >
5827 lnextself(nextedge);
5828 sym(nextedge, topcasing);
5829 lnextself(nextedge);
5830 sym(nextedge, sidecasing);
5831 bond(nextedge, topcasing);
5832 bond(leftcand, sidecasing);
5833 lnextself(leftcand);
5834 sym(leftcand, outercasing);
5835 lprevself(nextedge);
5836 bond(nextedge, outercasing);
5838 setorg(leftcand, lowerleft);
5839 setdest(leftcand, NULL);
5840 setapex(leftcand, nextapex);
5841 setorg(nextedge, NULL);
5842 setdest(nextedge, upperleft);
5843 setapex(nextedge, nextapex);
5845 upperleft = nextapex;
5847 otricopy(sidecasing, nextedge);
5848 apex(nextedge, nextapex);
5849 if (nextapex != (vertex) NULL) {
5851 badedge = incircle(m, b, lowerleft, lowerright, upperleft,
5861 if (!rightfinished) {
5863 lnext(rightcand, nextedge);
5865 apex(nextedge, nextapex);
5868 if (nextapex != (vertex) NULL) {
5870 badedge = incircle(m, b, lowerleft, lowerright, upperright, nextapex) >
5875 lprevself(nextedge);
5876 sym(nextedge, topcasing);
5877 lprevself(nextedge);
5878 sym(nextedge, sidecasing);
5879 bond(nextedge, topcasing);
5880 bond(rightcand, sidecasing);
5881 lprevself(rightcand);
5882 sym(rightcand, outercasing);
5883 lnextself(nextedge);
5884 bond(nextedge, outercasing);
5886 setorg(rightcand, NULL);
5887 setdest(rightcand, lowerright);
5888 setapex(rightcand, nextapex);
5889 setorg(nextedge, upperright);
5890 setdest(nextedge, NULL);
5891 setapex(nextedge, nextapex);
5893 upperright = nextapex;
5895 otricopy(sidecasing, nextedge);
5896 apex(nextedge, nextapex);
5897 if (nextapex != (vertex) NULL) {
5899 badedge = incircle(m, b, lowerleft, lowerright, upperright,
5908 if (leftfinished || (!rightfinished &&
5909 (incircle(m, b, upperleft, lowerleft, lowerright, upperright) >
5913 bond(baseedge, rightcand);
5914 lprev(rightcand, baseedge);
5915 setdest(baseedge, lowerleft);
5916 lowerright = upperright;
5917 sym(baseedge, rightcand);
5918 apex(rightcand, upperright);
5922 bond(baseedge, leftcand);
5923 lnext(leftcand, baseedge);
5924 setorg(baseedge, lowerright);
5925 lowerleft = upperleft;
5926 sym(baseedge, leftcand);
5927 apex(leftcand, upperleft);
5929 if (b->verbose > 2) {
5930 printf(
" Connecting ");
5931 printtriangle(m, b, &baseedge);
5953 void divconqrecurse(
struct mesh *m,
struct behavior *b, vertex *sortarray,
5954 int vertices,
int axis,
5955 struct otri *farleft,
struct otri *farright)
5957 struct otri midtri, tri1, tri2, tri3;
5958 struct otri innerleft, innerright;
5962 if (b->verbose > 2) {
5963 printf(
" Triangulating %d vertices.\n", vertices);
5965 if (vertices == 2) {
5968 maketriangle(m, b, farleft);
5969 setorg(*farleft, sortarray[0]);
5970 setdest(*farleft, sortarray[1]);
5972 maketriangle(m, b, farright);
5973 setorg(*farright, sortarray[1]);
5974 setdest(*farright, sortarray[0]);
5976 bond(*farleft, *farright);
5977 lprevself(*farleft);
5978 lnextself(*farright);
5979 bond(*farleft, *farright);
5980 lprevself(*farleft);
5981 lnextself(*farright);
5982 bond(*farleft, *farright);
5983 if (b->verbose > 2) {
5984 printf(
" Creating ");
5985 printtriangle(m, b, farleft);
5986 printf(
" Creating ");
5987 printtriangle(m, b, farright);
5990 lprev(*farright, *farleft);
5992 }
else if (vertices == 3) {
5996 maketriangle(m, b, &midtri);
5997 maketriangle(m, b, &tri1);
5998 maketriangle(m, b, &tri2);
5999 maketriangle(m, b, &tri3);
6000 area = counterclockwise(m, b, sortarray[0], sortarray[1], sortarray[2]);
6003 setorg(midtri, sortarray[0]);
6004 setdest(midtri, sortarray[1]);
6005 setorg(tri1, sortarray[1]);
6006 setdest(tri1, sortarray[0]);
6007 setorg(tri2, sortarray[2]);
6008 setdest(tri2, sortarray[1]);
6009 setorg(tri3, sortarray[1]);
6010 setdest(tri3, sortarray[2]);
6027 otricopy(tri1, *farleft);
6029 otricopy(tri2, *farright);
6033 setorg(midtri, sortarray[0]);
6034 setdest(tri1, sortarray[0]);
6035 setorg(tri3, sortarray[0]);
6039 setdest(midtri, sortarray[1]);
6040 setorg(tri1, sortarray[1]);
6041 setdest(tri2, sortarray[1]);
6042 setapex(midtri, sortarray[2]);
6043 setorg(tri2, sortarray[2]);
6044 setdest(tri3, sortarray[2]);
6047 setdest(midtri, sortarray[2]);
6048 setorg(tri1, sortarray[2]);
6049 setdest(tri2, sortarray[2]);
6050 setapex(midtri, sortarray[1]);
6051 setorg(tri2, sortarray[1]);
6052 setdest(tri3, sortarray[1]);
6070 otricopy(tri1, *farleft);
6073 otricopy(tri2, *farright);
6075 lnext(*farleft, *farright);
6078 if (b->verbose > 2) {
6079 printf(
" Creating ");
6080 printtriangle(m, b, &midtri);
6081 printf(
" Creating ");
6082 printtriangle(m, b, &tri1);
6083 printf(
" Creating ");
6084 printtriangle(m, b, &tri2);
6085 printf(
" Creating ");
6086 printtriangle(m, b, &tri3);
6091 divider = vertices >> 1;
6093 divconqrecurse(m, b, sortarray, divider, 1 - axis, farleft, &innerleft);
6094 divconqrecurse(m, b, &sortarray[divider], vertices - divider, 1 - axis,
6095 &innerright, farright);
6096 if (b->verbose > 1) {
6097 printf(
" Joining triangulations with %d and %d vertices.\n", divider,
6098 vertices - divider);
6101 mergehulls(m, b, farleft, &innerleft, &innerright, farright, axis);
6105 long removeghosts(
struct mesh *m,
struct behavior *b,
struct otri *startghost)
6107 struct otri searchedge;
6108 struct otri dissolveedge;
6109 struct otri deadtriangle;
6115 printf(
" Removing ghost triangles.\n");
6118 lprev(*startghost, searchedge);
6119 symself(searchedge);
6120 m->dummytri[0] = encode(searchedge);
6122 otricopy(*startghost, dissolveedge);
6126 lnext(dissolveedge, deadtriangle);
6127 lprevself(dissolveedge);
6128 symself(dissolveedge);
6133 if (dissolveedge.tri != m->dummytri) {
6134 org(dissolveedge, markorg);
6135 if (vertexmark(markorg) == 0) {
6136 setvertexmark(markorg, 1);
6141 dissolve(dissolveedge);
6143 sym(deadtriangle, dissolveedge);
6145 triangledealloc(m, deadtriangle.tri);
6146 }
while (!otriequal(dissolveedge, *startghost));
6160 long divconqdelaunay(
struct mesh *m,
struct behavior *b)
6163 struct otri hullleft, hullright;
6168 printf(
" Sorting vertices.\n");
6172 sortarray = (vertex *) trimalloc(m->invertices * (
int)
sizeof(vertex));
6173 traversalinit(&m->vertices);
6174 for (i = 0; i < m->invertices; i++) {
6175 sortarray[i] = vertextraverse(m);
6178 vertexsort(sortarray, m->invertices);
6181 for (j = 1; j < m->invertices; j++) {
6182 if ((sortarray[i][0] == sortarray[j][0])
6183 && (sortarray[i][1] == sortarray[j][1])) {
6186 "Warning: A duplicate vertex at (%.12g, %.12g) appeared and was ignored.\n",
6187 sortarray[j][0], sortarray[j][1]);
6189 setvertextype(sortarray[j], UNDEADVERTEX);
6193 sortarray[i] = sortarray[j];
6200 if (i - divider >= 2) {
6202 alternateaxes(sortarray, divider, 1);
6204 alternateaxes(&sortarray[divider], i - divider, 1);
6209 printf(
" Forming triangulation.\n");
6213 divconqrecurse(m, b, sortarray, i, 0, &hullleft, &hullright);
6214 trifree((
int *) sortarray);
6216 return removeghosts(m, b, &hullleft);
6233 long delaunay(
struct mesh *m,
struct behavior *b)
6238 initializetrisubpools(m, b);
6242 "Constructing Delaunay triangulation by divide-and-conquer method.\n");
6244 hulledges = divconqdelaunay(m, b);
6246 if (m->triangles.items == 0) {
6279 enum finddirectionresult finddirection(
struct mesh *m,
struct behavior *b,
6280 struct otri *searchtri,
6283 struct otri checktri;
6285 vertex leftvertex, rightvertex;
6286 float leftccw, rightccw;
6287 int leftflag, rightflag;
6290 org(*searchtri, startvertex);
6291 dest(*searchtri, rightvertex);
6292 apex(*searchtri, leftvertex);
6294 leftccw = counterclockwise(m, b, searchpoint, startvertex, leftvertex);
6295 leftflag = leftccw > 0.0;
6297 rightccw = counterclockwise(m, b, startvertex, searchpoint, rightvertex);
6298 rightflag = rightccw > 0.0;
6299 if (leftflag && rightflag) {
6302 onext(*searchtri, checktri);
6303 if (checktri.tri == m->dummytri) {
6311 onextself(*searchtri);
6312 if (searchtri->tri == m->dummytri) {
6313 printf(
"Internal error in finddirection(): Unable to find a\n");
6314 printf(
" triangle leading from (%.12g, %.12g) to", startvertex[0],
6316 printf(
" (%.12g, %.12g).\n", searchpoint[0], searchpoint[1]);
6319 apex(*searchtri, leftvertex);
6321 leftccw = counterclockwise(m, b, searchpoint, startvertex, leftvertex);
6322 leftflag = leftccw > 0.0;
6326 oprevself(*searchtri);
6327 if (searchtri->tri == m->dummytri) {
6328 printf(
"Internal error in finddirection(): Unable to find a\n");
6329 printf(
" triangle leading from (%.12g, %.12g) to", startvertex[0],
6331 printf(
" (%.12g, %.12g).\n", searchpoint[0], searchpoint[1]);
6334 dest(*searchtri, rightvertex);
6336 rightccw = counterclockwise(m, b, startvertex, searchpoint, rightvertex);
6337 rightflag = rightccw > 0.0;
6339 if (leftccw == 0.0) {
6340 return LEFTCOLLINEAR;
6341 }
else if (rightccw == 0.0) {
6342 return RIGHTCOLLINEAR;
6365 void segmentintersection(
struct mesh *m,
struct behavior *b,
6366 struct otri *splittri,
struct osub *splitsubseg,
6369 struct osub opposubseg;
6372 vertex leftvertex, rightvertex;
6374 enum insertvertexresult success;
6375 enum finddirectionresult collinear;
6385 apex(*splittri, endpoint1);
6386 org(*splittri, torg);
6387 dest(*splittri, tdest);
6389 tx = tdest[0] - torg[0];
6390 ty = tdest[1] - torg[1];
6391 ex = endpoint2[0] - endpoint1[0];
6392 ey = endpoint2[1] - endpoint1[1];
6393 etx = torg[0] - endpoint2[0];
6394 ety = torg[1] - endpoint2[1];
6395 denom = ty * ex - tx * ey;
6397 printf(
"Internal error in segmentintersection():");
6398 printf(
" Attempt to find intersection of parallel segments.\n");
6401 split = (ey * etx - ex * ety) / denom;
6403 newvertex = (vertex) poolalloc(&m->vertices);
6405 for (i = 0; i < 2 + m->nextras; i++) {
6406 newvertex[i] = torg[i] + split * (tdest[i] - torg[i]);
6408 setvertexmark(newvertex, mark(*splitsubseg));
6409 setvertextype(newvertex, INPUTVERTEX);
6410 if (b->verbose > 1) {
6412 " Splitting subsegment (%.12g, %.12g) (%.12g, %.12g) at (%.12g, %.12g).\n",
6413 torg[0], torg[1], tdest[0], tdest[1], newvertex[0], newvertex[1]);
6416 success = insertvertex(m, b, newvertex, splittri, splitsubseg, 0, 0);
6417 if (success != SUCCESSFULVERTEX) {
6418 printf(
"Internal error in segmentintersection():\n");
6419 printf(
" Failure to split a segment.\n");
6423 setvertex2tri(newvertex, encode(*splittri));
6424 if (m->steinerleft > 0) {
6429 ssymself(*splitsubseg);
6430 spivot(*splitsubseg, opposubseg);
6431 sdissolve(*splitsubseg);
6432 sdissolve(opposubseg);
6434 setsegorg(*splitsubseg, newvertex);
6435 snextself(*splitsubseg);
6436 }
while (splitsubseg->ss != m->dummysub);
6438 setsegorg(opposubseg, newvertex);
6439 snextself(opposubseg);
6440 }
while (opposubseg.ss != m->dummysub);
6444 collinear = finddirection(m, b, splittri, endpoint1);
6445 dest(*splittri, rightvertex);
6446 apex(*splittri, leftvertex);
6447 if ((leftvertex[0] == endpoint1[0]) && (leftvertex[1] == endpoint1[1])) {
6448 onextself(*splittri);
6449 }
else if ((rightvertex[0] != endpoint1[0]) ||
6450 (rightvertex[1] != endpoint1[1])) {
6451 printf(
"Internal error in segmentintersection():\n");
6452 printf(
" Topological inconsistency after splitting a segment.\n");
6482 int scoutsegment(
struct mesh *m,
struct behavior *b,
struct otri *searchtri,
6483 vertex endpoint2,
int newmark)
6485 struct otri crosstri;
6486 struct osub crosssubseg;
6487 vertex leftvertex, rightvertex;
6488 enum finddirectionresult collinear;
6491 collinear = finddirection(m, b, searchtri, endpoint2);
6492 dest(*searchtri, rightvertex);
6493 apex(*searchtri, leftvertex);
6494 if (((leftvertex[0] == endpoint2[0]) && (leftvertex[1] == endpoint2[1])) ||
6495 ((rightvertex[0] == endpoint2[0]) && (rightvertex[1] == endpoint2[1]))) {
6497 if ((leftvertex[0] == endpoint2[0]) && (leftvertex[1] == endpoint2[1])) {
6498 lprevself(*searchtri);
6501 insertsubseg(m, b, searchtri, newmark);
6503 }
else if (collinear == LEFTCOLLINEAR) {
6506 lprevself(*searchtri);
6507 insertsubseg(m, b, searchtri, newmark);
6509 return scoutsegment(m, b, searchtri, endpoint2, newmark);
6510 }
else if (collinear == RIGHTCOLLINEAR) {
6512 insertsubseg(m, b, searchtri, newmark);
6514 lnextself(*searchtri);
6516 return scoutsegment(m, b, searchtri, endpoint2, newmark);
6518 lnext(*searchtri, crosstri);
6519 tspivot(crosstri, crosssubseg);
6521 if (crosssubseg.ss == m->dummysub) {
6525 segmentintersection(m, b, &crosstri, &crosssubseg, endpoint2);
6526 otricopy(crosstri, *searchtri);
6527 insertsubseg(m, b, searchtri, newmark);
6529 return scoutsegment(m, b, searchtri, endpoint2, newmark);
6572 void delaunayfixup(
struct mesh *m,
struct behavior *b,
6573 struct otri *fixuptri,
int leftside)
6575 struct otri neartri;
6577 struct osub faredge;
6578 vertex nearvertex, leftvertex, rightvertex, farvertex;
6582 lnext(*fixuptri, neartri);
6583 sym(neartri, fartri);
6585 if (fartri.tri == m->dummytri) {
6588 tspivot(neartri, faredge);
6589 if (faredge.ss != m->dummysub) {
6593 apex(neartri, nearvertex);
6594 org(neartri, leftvertex);
6595 dest(neartri, rightvertex);
6596 apex(fartri, farvertex);
6599 if (counterclockwise(m, b, nearvertex, leftvertex, farvertex) <= 0.0) {
6605 if (counterclockwise(m, b, farvertex, rightvertex, nearvertex) <= 0.0) {
6611 if (counterclockwise(m, b, rightvertex, leftvertex, farvertex) > 0.0) {
6616 if (incircle(m, b, leftvertex, farvertex, rightvertex, nearvertex) <=
6622 flip(m, b, &neartri);
6623 lprevself(*fixuptri);
6625 delaunayfixup(m, b, fixuptri, leftside);
6626 delaunayfixup(m, b, &fartri, leftside);
6683 void constrainededge(
struct mesh *m,
struct behavior *b,
6684 struct otri *starttri, vertex endpoint2,
int newmark)
6686 struct otri fixuptri, fixuptri2;
6687 struct osub crosssubseg;
6696 org(*starttri, endpoint1);
6697 lnext(*starttri, fixuptri);
6698 flip(m, b, &fixuptri);
6704 org(fixuptri, farvertex);
6707 if ((farvertex[0] == endpoint2[0]) && (farvertex[1] == endpoint2[1])) {
6708 oprev(fixuptri, fixuptri2);
6710 delaunayfixup(m, b, &fixuptri, 0);
6711 delaunayfixup(m, b, &fixuptri2, 1);
6717 area = counterclockwise(m, b, endpoint1, endpoint2, farvertex);
6721 oprev(fixuptri, fixuptri2);
6723 delaunayfixup(m, b, &fixuptri, 0);
6724 delaunayfixup(m, b, &fixuptri2, 1);
6728 oprev(fixuptri, fixuptri2);
6731 delaunayfixup(m, b, &fixuptri2, 1);
6735 lprevself(fixuptri);
6737 delaunayfixup(m, b, &fixuptri, 0);
6741 oprevself(fixuptri);
6744 tspivot(fixuptri, crosssubseg);
6745 if (crosssubseg.ss == m->dummysub) {
6746 flip(m, b, &fixuptri);
6751 segmentintersection(m, b, &fixuptri, &crosssubseg, endpoint2);
6758 insertsubseg(m, b, &fixuptri, newmark);
6763 if (!scoutsegment(m, b, &fixuptri, endpoint2, newmark)) {
6764 constrainededge(m, b, &fixuptri, endpoint2, newmark);
6775 void insertsegment(
struct mesh *m,
struct behavior *b,
6776 vertex endpoint1, vertex endpoint2,
int newmark)
6778 struct otri searchtri1, searchtri2;
6779 triangle encodedtri;
6783 if (b->verbose > 1) {
6784 printf(
" Connecting (%.12g, %.12g) to (%.12g, %.12g).\n",
6785 endpoint1[0], endpoint1[1], endpoint2[0], endpoint2[1]);
6789 checkvertex = (vertex) NULL;
6790 encodedtri = vertex2tri(endpoint1);
6791 if (encodedtri != (triangle) NULL) {
6792 decode(encodedtri, searchtri1);
6793 org(searchtri1, checkvertex);
6795 if (checkvertex != endpoint1) {
6797 searchtri1.tri = m->dummytri;
6798 searchtri1.orient = 0;
6799 symself(searchtri1);
6801 if (locate(m, b, endpoint1, &searchtri1) != ONVERTEX) {
6803 "Internal error in insertsegment(): Unable to locate PSLG vertex\n");
6804 printf(
" (%.12g, %.12g) in triangulation.\n",
6805 endpoint1[0], endpoint1[1]);
6810 otricopy(searchtri1, m->recenttri);
6813 if (scoutsegment(m, b, &searchtri1, endpoint2, newmark)) {
6819 org(searchtri1, endpoint1);
6822 checkvertex = (vertex) NULL;
6823 encodedtri = vertex2tri(endpoint2);
6824 if (encodedtri != (triangle) NULL) {
6825 decode(encodedtri, searchtri2);
6826 org(searchtri2, checkvertex);
6828 if (checkvertex != endpoint2) {
6830 searchtri2.tri = m->dummytri;
6831 searchtri2.orient = 0;
6832 symself(searchtri2);
6834 if (locate(m, b, endpoint2, &searchtri2) != ONVERTEX) {
6836 "Internal error in insertsegment(): Unable to locate PSLG vertex\n");
6837 printf(
" (%.12g, %.12g) in triangulation.\n",
6838 endpoint2[0], endpoint2[1]);
6843 otricopy(searchtri2, m->recenttri);
6846 if (scoutsegment(m, b, &searchtri2, endpoint1, newmark)) {
6852 org(searchtri2, endpoint2);
6855 constrainededge(m, b, &searchtri1, endpoint2, newmark);
6864 void markhull(
struct mesh *m,
struct behavior *b)
6866 struct otri hulltri;
6867 struct otri nexttri;
6868 struct otri starttri;
6872 hulltri.tri = m->dummytri;
6876 otricopy(hulltri, starttri);
6880 insertsubseg(m, b, &hulltri, 1);
6883 oprev(hulltri, nexttri);
6884 while (nexttri.tri != m->dummytri) {
6885 otricopy(nexttri, hulltri);
6886 oprev(hulltri, nexttri);
6888 }
while (!otriequal(hulltri, starttri));
6901 void formskeleton(
struct mesh *m,
struct behavior *b,
int *segmentlist,
6902 int *segmentmarkerlist,
int numberofsegments)
6904 char polyfilename[6];
6906 vertex endpoint1, endpoint2;
6914 printf(
"Recovering segments in Delaunay triangulation.\n");
6916 strcpy(polyfilename,
"input");
6917 m->insegments = numberofsegments;
6918 segmentmarkers = segmentmarkerlist != (
int *) NULL;
6922 if (m->triangles.items == 0) {
6928 if (m->insegments > 0) {
6929 makevertexmap(m, b);
6931 printf(
" Recovering PSLG segments.\n");
6937 for (i = 0; i < m->insegments; i++) {
6938 end1 = segmentlist[index++];
6939 end2 = segmentlist[index++];
6940 if (segmentmarkers) {
6941 boundmarker = segmentmarkerlist[i];
6943 if ((end1 < b->firstnumber) ||
6944 (end1 >= b->firstnumber + m->invertices)) {
6946 printf(
"Warning: Invalid first endpoint of segment %d in %s.\n",
6947 b->firstnumber + i, polyfilename);
6949 }
else if ((end2 < b->firstnumber) ||
6950 (end2 >= b->firstnumber + m->invertices)) {
6952 printf(
"Warning: Invalid second endpoint of segment %d in %s.\n",
6953 b->firstnumber + i, polyfilename);
6957 endpoint1 = getvertex(m, b, end1);
6958 endpoint2 = getvertex(m, b, end2);
6959 if ((endpoint1[0] == endpoint2[0]) && (endpoint1[1] == endpoint2[1])) {
6961 printf(
"Warning: Endpoints of segment %d are coincident in %s.\n",
6962 b->firstnumber + i, polyfilename);
6965 insertsegment(m, b, endpoint1, endpoint2, boundmarker);
6972 if (b->convex || !b->poly) {
6975 printf(
" Enclosing convex hull with segments.\n");
6997 void infecthull(
struct mesh *m,
struct behavior *b)
6999 struct otri hulltri;
7000 struct otri nexttri;
7001 struct otri starttri;
7002 struct osub hullsubseg;
7003 triangle **deadtriangle;
7009 printf(
" Marking concavities (external triangles) for elimination.\n");
7012 hulltri.tri = m->dummytri;
7016 otricopy(hulltri, starttri);
7020 if (!infected(hulltri)) {
7022 tspivot(hulltri, hullsubseg);
7023 if (hullsubseg.ss == m->dummysub) {
7025 if (!infected(hulltri)) {
7027 deadtriangle = (triangle **) poolalloc(&m->viri);
7028 *deadtriangle = hulltri.tri;
7032 if (mark(hullsubseg) == 0) {
7033 setmark(hullsubseg, 1);
7035 dest(hulltri, hdest);
7036 if (vertexmark(horg) == 0) {
7037 setvertexmark(horg, 1);
7039 if (vertexmark(hdest) == 0) {
7040 setvertexmark(hdest, 1);
7047 oprev(hulltri, nexttri);
7048 while (nexttri.tri != m->dummytri) {
7049 otricopy(nexttri, hulltri);
7050 oprev(hulltri, nexttri);
7052 }
while (!otriequal(hulltri, starttri));
7072 void plague(
struct mesh *m,
struct behavior *b)
7074 struct otri testtri;
7075 struct otri neighbor;
7076 triangle **virusloop;
7077 triangle **deadtriangle;
7078 struct osub neighborsubseg;
7081 vertex deadorg, deaddest, deadapex;
7087 printf(
" Marking neighbors of marked triangles.\n");
7091 traversalinit(&m->viri);
7092 virusloop = (triangle **) traverse(&m->viri);
7093 while (virusloop != (triangle **) NULL) {
7094 testtri.tri = *virusloop;
7100 if (b->verbose > 2) {
7104 org(testtri, deadorg);
7105 dest(testtri, deaddest);
7106 apex(testtri, deadapex);
7107 printf(
" Checking (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
7108 deadorg[0], deadorg[1], deaddest[0], deaddest[1],
7109 deadapex[0], deadapex[1]);
7112 for (testtri.orient = 0; testtri.orient < 3; testtri.orient++) {
7114 sym(testtri, neighbor);
7116 tspivot(testtri, neighborsubseg);
7118 if ((neighbor.tri == m->dummytri) || infected(neighbor)) {
7119 if (neighborsubseg.ss != m->dummysub) {
7123 subsegdealloc(m, neighborsubseg.ss);
7124 if (neighbor.tri != m->dummytri) {
7128 tsdissolve(neighbor);
7133 if (neighborsubseg.ss == m->dummysub) {
7136 if (b->verbose > 2) {
7137 org(neighbor, deadorg);
7138 dest(neighbor, deaddest);
7139 apex(neighbor, deadapex);
7141 " Marking (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
7142 deadorg[0], deadorg[1], deaddest[0], deaddest[1],
7143 deadapex[0], deadapex[1]);
7147 deadtriangle = (triangle **) poolalloc(&m->viri);
7148 *deadtriangle = neighbor.tri;
7151 stdissolve(neighborsubseg);
7153 if (mark(neighborsubseg) == 0) {
7154 setmark(neighborsubseg, 1);
7156 org(neighbor, norg);
7157 dest(neighbor, ndest);
7158 if (vertexmark(norg) == 0) {
7159 setvertexmark(norg, 1);
7161 if (vertexmark(ndest) == 0) {
7162 setvertexmark(ndest, 1);
7170 virusloop = (triangle **) traverse(&m->viri);
7174 printf(
" Deleting marked triangles.\n");
7177 traversalinit(&m->viri);
7178 virusloop = (triangle **) traverse(&m->viri);
7179 while (virusloop != (triangle **) NULL) {
7180 testtri.tri = *virusloop;
7185 for (testtri.orient = 0; testtri.orient < 3; testtri.orient++) {
7186 org(testtri, testvertex);
7188 if (testvertex != (vertex) NULL) {
7191 setorg(testtri, NULL);
7193 onext(testtri, neighbor);
7195 while ((neighbor.tri != m->dummytri) &&
7196 (!otriequal(neighbor, testtri))) {
7197 if (infected(neighbor)) {
7199 setorg(neighbor, NULL);
7205 onextself(neighbor);
7208 if (neighbor.tri == m->dummytri) {
7210 oprev(testtri, neighbor);
7212 while (neighbor.tri != m->dummytri) {
7213 if (infected(neighbor)) {
7215 setorg(neighbor, NULL);
7221 oprevself(neighbor);
7225 if (b->verbose > 1) {
7226 printf(
" Deleting vertex (%.12g, %.12g)\n",
7227 testvertex[0], testvertex[1]);
7229 setvertextype(testvertex, UNDEADVERTEX);
7237 for (testtri.orient = 0; testtri.orient < 3; testtri.orient++) {
7238 sym(testtri, neighbor);
7239 if (neighbor.tri == m->dummytri) {
7253 triangledealloc(m, testtri.tri);
7254 virusloop = (triangle **) traverse(&m->viri);
7257 poolrestart(&m->viri);
7275 void regionplague(
struct mesh *m,
struct behavior *b,
7276 float attribute,
float area)
7278 struct otri testtri;
7279 struct otri neighbor;
7280 triangle **virusloop;
7281 triangle **regiontri;
7282 struct osub neighborsubseg;
7283 vertex regionorg, regiondest, regionapex;
7287 if (b->verbose > 1) {
7288 printf(
" Marking neighbors of marked triangles.\n");
7293 traversalinit(&m->viri);
7294 virusloop = (triangle **) traverse(&m->viri);
7295 while (virusloop != (triangle **) NULL) {
7296 testtri.tri = *virusloop;
7302 if (b->regionattrib) {
7304 setelemattribute(testtri, m->eextras, attribute);
7308 setareabound(testtri, area);
7310 if (b->verbose > 2) {
7314 org(testtri, regionorg);
7315 dest(testtri, regiondest);
7316 apex(testtri, regionapex);
7317 printf(
" Checking (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
7318 regionorg[0], regionorg[1], regiondest[0], regiondest[1],
7319 regionapex[0], regionapex[1]);
7322 for (testtri.orient = 0; testtri.orient < 3; testtri.orient++) {
7324 sym(testtri, neighbor);
7326 tspivot(testtri, neighborsubseg);
7329 if ((neighbor.tri != m->dummytri) && !infected(neighbor)
7330 && (neighborsubseg.ss == m->dummysub)) {
7331 if (b->verbose > 2) {
7332 org(neighbor, regionorg);
7333 dest(neighbor, regiondest);
7334 apex(neighbor, regionapex);
7335 printf(
" Marking (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
7336 regionorg[0], regionorg[1], regiondest[0], regiondest[1],
7337 regionapex[0], regionapex[1]);
7342 regiontri = (triangle **) poolalloc(&m->viri);
7343 *regiontri = neighbor.tri;
7349 virusloop = (triangle **) traverse(&m->viri);
7353 if (b->verbose > 1) {
7354 printf(
" Unmarking marked triangles.\n");
7356 traversalinit(&m->viri);
7357 virusloop = (triangle **) traverse(&m->viri);
7358 while (virusloop != (triangle **) NULL) {
7359 testtri.tri = *virusloop;
7361 virusloop = (triangle **) traverse(&m->viri);
7364 poolrestart(&m->viri);
7379 void carveholes(
struct mesh *m,
struct behavior *b,
float *holelist,
int holes,
7380 float *regionlist,
int regions)
7382 struct otri searchtri;
7383 struct otri triangleloop;
7384 struct otri *regiontris;
7386 triangle **regiontri;
7387 vertex searchorg, searchdest;
7388 enum locateresult intersect;
7392 if (!(b->quiet || (b->noholes && b->convex))) {
7393 printf(
"Removing unwanted triangles.\n");
7394 if (b->verbose && (holes > 0)) {
7395 printf(
" Marking holes for elimination.\n");
7401 regiontris = (
struct otri *) trimalloc(regions *
7402 (
int)
sizeof(
struct otri));
7404 regiontris = (
struct otri *) NULL;
7407 if (((holes > 0) && !b->noholes) || !b->convex || (regions > 0)) {
7410 poolinit(&m->viri,
sizeof(triangle *), VIRUSPERBLOCK, VIRUSPERBLOCK, 0);
7419 if ((holes > 0) && !b->noholes) {
7421 for (i = 0; i < 2 * holes; i += 2) {
7423 if ((holelist[i] >= m->xmin) && (holelist[i] <= m->xmax)
7424 && (holelist[i + 1] >= m->ymin) && (holelist[i + 1] <= m->ymax)) {
7426 searchtri.tri = m->dummytri;
7427 searchtri.orient = 0;
7432 org(searchtri, searchorg);
7433 dest(searchtri, searchdest);
7434 if (counterclockwise(m, b, searchorg, searchdest, &holelist[i]) >
7437 intersect = locate(m, b, &holelist[i], &searchtri);
7438 if ((intersect != OUTSIDE) && (!infected(searchtri))) {
7442 holetri = (triangle **) poolalloc(&m->viri);
7443 *holetri = searchtri.tri;
7458 for (i = 0; i < regions; i++) {
7459 regiontris[i].tri = m->dummytri;
7461 if ((regionlist[4 * i] >= m->xmin) && (regionlist[4 * i] <= m->xmax) &&
7462 (regionlist[4 * i + 1] >= m->ymin) &&
7463 (regionlist[4 * i + 1] <= m->ymax)) {
7465 searchtri.tri = m->dummytri;
7466 searchtri.orient = 0;
7471 org(searchtri, searchorg);
7472 dest(searchtri, searchdest);
7473 if (counterclockwise(m, b, searchorg, searchdest, ®ionlist[4 * i]) >
7476 intersect = locate(m, b, ®ionlist[4 * i], &searchtri);
7477 if ((intersect != OUTSIDE) && (!infected(searchtri))) {
7480 otricopy(searchtri, regiontris[i]);
7487 if (m->viri.items > 0) {
7495 if (b->regionattrib) {
7497 printf(
"Spreading regional attributes and area constraints.\n");
7499 printf(
"Spreading regional attributes.\n");
7502 printf(
"Spreading regional area constraints.\n");
7505 if (b->regionattrib && !b->refine) {
7507 traversalinit(&m->triangles);
7508 triangleloop.orient = 0;
7509 triangleloop.tri = triangletraverse(m);
7510 while (triangleloop.tri != (triangle *) NULL) {
7511 setelemattribute(triangleloop, m->eextras, 0.0);
7512 triangleloop.tri = triangletraverse(m);
7515 for (i = 0; i < regions; i++) {
7516 if (regiontris[i].tri != m->dummytri) {
7519 if (!deadtri(regiontris[i].tri)) {
7521 infect(regiontris[i]);
7522 regiontri = (triangle **) poolalloc(&m->viri);
7523 *regiontri = regiontris[i].tri;
7525 regionplague(m, b, regionlist[4 * i + 2], regionlist[4 * i + 3]);
7530 if (b->regionattrib && !b->refine) {
7537 if (((holes > 0) && !b->noholes) || !b->convex || (regions > 0)) {
7538 pooldeinit(&m->viri);
7541 trifree((
int *) regiontris);
7555 void highorder(
struct mesh *m,
struct behavior *b)
7557 struct otri triangleloop, trisym;
7558 struct osub checkmark;
7566 printf(
"Adding vertices for second-order triangles.\n");
7573 m->vertices.deaditemstack = (
int *) NULL;
7575 traversalinit(&m->triangles);
7576 triangleloop.tri = triangletraverse(m);
7583 while (triangleloop.tri != (triangle *) NULL) {
7584 for (triangleloop.orient = 0; triangleloop.orient < 3;
7585 triangleloop.orient++) {
7586 sym(triangleloop, trisym);
7587 if ((triangleloop.tri < trisym.tri) || (trisym.tri == m->dummytri)) {
7588 org(triangleloop, torg);
7589 dest(triangleloop, tdest);
7592 newvertex = (vertex) poolalloc(&m->vertices);
7593 for (i = 0; i < 2 + m->nextras; i++) {
7594 newvertex[i] = 0.5 * (torg[i] + tdest[i]);
7598 setvertexmark(newvertex, trisym.tri == m->dummytri);
7599 setvertextype(newvertex,
7600 trisym.tri == m->dummytri ? FREEVERTEX : SEGMENTVERTEX);
7601 if (b->usesegments) {
7602 tspivot(triangleloop, checkmark);
7604 if (checkmark.ss != m->dummysub) {
7605 setvertexmark(newvertex, mark(checkmark));
7606 setvertextype(newvertex, SEGMENTVERTEX);
7609 if (b->verbose > 1) {
7610 printf(
" Creating (%.12g, %.12g).\n", newvertex[0], newvertex[1]);
7613 triangleloop.tri[m->highorderindex + triangleloop.orient] =
7614 (triangle) newvertex;
7615 if (trisym.tri != m->dummytri) {
7616 trisym.tri[m->highorderindex + trisym.orient] = (triangle) newvertex;
7620 triangleloop.tri = triangletraverse(m);
7634 void transfernodes(
struct mesh *m,
struct behavior *b,
float *pointlist,
7635 float *pointattriblist,
int *pointmarkerlist,
7636 int numberofpoints,
int numberofpointattribs)
7644 m->invertices = numberofpoints;
7646 m->nextras = numberofpointattribs;
7647 m->readnodefile = 0;
7648 if (m->invertices < 3) {
7649 printf(
"Error: Input must have at least three input vertices.\n");
7652 if (m->nextras == 0) {
7656 initializevertexpool(m, b);
7661 for (i = 0; i < m->invertices; i++) {
7662 vertexloop = (vertex) poolalloc(&m->vertices);
7664 x = vertexloop[0] = pointlist[coordindex++];
7665 y = vertexloop[1] = pointlist[coordindex++];
7667 for (j = 0; j < numberofpointattribs; j++) {
7668 vertexloop[2 + j] = pointattriblist[attribindex++];
7670 if (pointmarkerlist != (
int *) NULL) {
7672 setvertexmark(vertexloop, pointmarkerlist[i]);
7675 setvertexmark(vertexloop, 0);
7677 setvertextype(vertexloop, INPUTVERTEX);
7680 m->xmin = m->xmax = x;
7681 m->ymin = m->ymax = y;
7683 m->xmin = (x < m->xmin) ? x : m->xmin;
7684 m->xmax = (x > m->xmax) ? x : m->xmax;
7685 m->ymin = (y < m->ymin) ? y : m->ymin;
7686 m->ymax = (y > m->ymax) ? y : m->ymax;
7692 m->xminextreme = 10 * m->xmin - 9 * m->xmax;
7704 void writenodes(
struct mesh *m,
struct behavior *b,
float **pointlist,
7705 float **pointattriblist,
int **pointmarkerlist)
7718 outvertices = m->vertices.items - m->undeads;
7720 outvertices = m->vertices.items;
7724 printf(
"Writing vertices.\n");
7727 if (*pointlist == (
float *) NULL) {
7728 *pointlist = (
float *) trimalloc((
int) (outvertices * 2 *
sizeof(float)));
7731 if ((m->nextras > 0) && (*pointattriblist == (
float *) NULL)) {
7732 *pointattriblist = (
float *) trimalloc((
int) (outvertices * m->nextras *
7736 if (!b->nobound && (*pointmarkerlist == (
int *) NULL)) {
7737 *pointmarkerlist = (
int *) trimalloc((
int) (outvertices *
sizeof(int)));
7740 palist = *pointattriblist;
7741 pmlist = *pointmarkerlist;
7744 traversalinit(&m->vertices);
7745 vertexnumber = b->firstnumber;
7746 vertexloop = vertextraverse(m);
7747 while (vertexloop != (vertex) NULL) {
7748 if (!b->jettison || (vertextype(vertexloop) != UNDEADVERTEX)) {
7750 plist[coordindex++] = vertexloop[0];
7751 plist[coordindex++] = vertexloop[1];
7753 for (i = 0; i < m->nextras; i++) {
7754 palist[attribindex++] = vertexloop[2 + i];
7758 pmlist[vertexnumber - b->firstnumber] = vertexmark(vertexloop);
7760 setvertexmark(vertexloop, vertexnumber);
7763 vertexloop = vertextraverse(m);
7777 void numbernodes(
struct mesh *m,
struct behavior *b)
7782 traversalinit(&m->vertices);
7783 vertexnumber = b->firstnumber;
7784 vertexloop = vertextraverse(m);
7785 while (vertexloop != (vertex) NULL) {
7786 setvertexmark(vertexloop, vertexnumber);
7787 if (!b->jettison || (vertextype(vertexloop) != UNDEADVERTEX)) {
7790 vertexloop = vertextraverse(m);
7800 void writeelements(
struct mesh *m,
struct behavior *b,
7801 int **trianglelist,
float **triangleattriblist)
7807 struct otri triangleloop;
7809 vertex mid1, mid2, mid3;
7814 printf(
"Writing triangles.\n");
7817 if (*trianglelist == (
int *) NULL) {
7818 *trianglelist = (
int *) trimalloc((
int) (m->triangles.items *
7819 ((b->order + 1) * (b->order + 2) /
7823 if ((m->eextras > 0) && (*triangleattriblist == (
float *) NULL)) {
7824 *triangleattriblist = (
float *) trimalloc((
int) (m->triangles.items *
7828 tlist = *trianglelist;
7829 talist = *triangleattriblist;
7832 traversalinit(&m->triangles);
7833 triangleloop.tri = triangletraverse(m);
7834 triangleloop.orient = 0;
7835 elementnumber = b->firstnumber;
7836 while (triangleloop.tri != (triangle *) NULL) {
7837 org(triangleloop, p1);
7838 dest(triangleloop, p2);
7839 apex(triangleloop, p3);
7840 if (b->order == 1) {
7841 tlist[vertexindex++] = vertexmark(p1);
7842 tlist[vertexindex++] = vertexmark(p2);
7843 tlist[vertexindex++] = vertexmark(p3);
7845 mid1 = (vertex) triangleloop.tri[m->highorderindex + 1];
7846 mid2 = (vertex) triangleloop.tri[m->highorderindex + 2];
7847 mid3 = (vertex) triangleloop.tri[m->highorderindex];
7848 tlist[vertexindex++] = vertexmark(p1);
7849 tlist[vertexindex++] = vertexmark(p2);
7850 tlist[vertexindex++] = vertexmark(p3);
7851 tlist[vertexindex++] = vertexmark(mid1);
7852 tlist[vertexindex++] = vertexmark(mid2);
7853 tlist[vertexindex++] = vertexmark(mid3);
7856 for (i = 0; i < m->eextras; i++) {
7857 talist[attribindex++] = elemattribute(triangleloop, i);
7859 triangleloop.tri = triangletraverse(m);
7870 void writepoly(
struct mesh *m,
struct behavior *b,
7871 int **segmentlist,
int **segmentmarkerlist)
7876 struct osub subsegloop;
7877 vertex endpoint1, endpoint2;
7881 printf(
"Writing segments.\n");
7884 if (*segmentlist == (
int *) NULL) {
7885 *segmentlist = (
int *) trimalloc((
int) (m->subsegs.items * 2 *
7889 if (!b->nobound && (*segmentmarkerlist == (
int *) NULL)) {
7890 *segmentmarkerlist = (
int *) trimalloc((
int) (m->subsegs.items *
7893 slist = *segmentlist;
7894 smlist = *segmentmarkerlist;
7897 traversalinit(&m->subsegs);
7898 subsegloop.ss = subsegtraverse(m);
7899 subsegloop.ssorient = 0;
7900 subsegnumber = b->firstnumber;
7901 while (subsegloop.ss != (subseg *) NULL) {
7902 sorg(subsegloop, endpoint1);
7903 sdest(subsegloop, endpoint2);
7905 slist[index++] = vertexmark(endpoint1);
7906 slist[index++] = vertexmark(endpoint2);
7909 smlist[subsegnumber - b->firstnumber] = mark(subsegloop);
7911 subsegloop.ss = subsegtraverse(m);
7922 void writeedges(
struct mesh *m,
struct behavior *b,
7923 int **edgelist,
int **edgemarkerlist)
7928 struct otri triangleloop, trisym;
7929 struct osub checkmark;
7936 printf(
"Writing edges.\n");
7939 if (*edgelist == (
int *) NULL) {
7940 *edgelist = (
int *) trimalloc((
int) (m->edges * 2 *
sizeof(int)));
7943 if (!b->nobound && (*edgemarkerlist == (
int *) NULL)) {
7944 *edgemarkerlist = (
int *) trimalloc((
int) (m->edges *
sizeof(int)));
7947 emlist = *edgemarkerlist;
7950 traversalinit(&m->triangles);
7951 triangleloop.tri = triangletraverse(m);
7952 edgenumber = b->firstnumber;
7959 while (triangleloop.tri != (triangle *) NULL) {
7960 for (triangleloop.orient = 0; triangleloop.orient < 3;
7961 triangleloop.orient++) {
7962 sym(triangleloop, trisym);
7963 if ((triangleloop.tri < trisym.tri) || (trisym.tri == m->dummytri)) {
7964 org(triangleloop, p1);
7965 dest(triangleloop, p2);
7966 elist[index++] = vertexmark(p1);
7967 elist[index++] = vertexmark(p2);
7972 if (b->usesegments) {
7973 tspivot(triangleloop, checkmark);
7974 if (checkmark.ss == m->dummysub) {
7975 emlist[edgenumber - b->firstnumber] = 0;
7977 emlist[edgenumber - b->firstnumber] = mark(checkmark);
7980 emlist[edgenumber - b->firstnumber] = trisym.tri == m->dummytri;
7986 triangleloop.tri = triangletraverse(m);
8006 void writevoronoi(
struct mesh *m,
struct behavior *b,
float **vpointlist,
8007 float **vpointattriblist,
int **vpointmarkerlist,
8008 int **vedgelist,
int **vedgemarkerlist,
float **vnormlist)
8016 struct otri triangleloop, trisym;
8017 vertex torg, tdest, tapex;
8018 float circumcenter[2];
8020 long vnodenumber, vedgenumber;
8026 printf(
"Writing Voronoi vertices.\n");
8029 if (*vpointlist == (
float *) NULL) {
8030 *vpointlist = (
float *) trimalloc((
int) (m->triangles.items * 2 *
8034 if (*vpointattriblist == (
float *) NULL) {
8035 *vpointattriblist = (
float *) trimalloc((
int) (m->triangles.items *
8036 m->nextras *
sizeof(float)));
8038 *vpointmarkerlist = (
int *) NULL;
8039 plist = *vpointlist;
8040 palist = *vpointattriblist;
8044 traversalinit(&m->triangles);
8045 triangleloop.tri = triangletraverse(m);
8046 triangleloop.orient = 0;
8047 vnodenumber = b->firstnumber;
8048 while (triangleloop.tri != (triangle *) NULL) {
8049 org(triangleloop, torg);
8050 dest(triangleloop, tdest);
8051 apex(triangleloop, tapex);
8052 findcircumcenter(m, b, torg, tdest, tapex, circumcenter, &xi, &eta, 0);
8055 plist[coordindex++] = circumcenter[0];
8056 plist[coordindex++] = circumcenter[1];
8057 for (i = 2; i < 2 + m->nextras; i++) {
8059 palist[attribindex++] = torg[i] + xi * (tdest[i] - torg[i])
8060 + eta * (tapex[i] - torg[i]);
8063 * (
int *) (triangleloop.tri + 6) = (int) vnodenumber;
8064 triangleloop.tri = triangletraverse(m);
8069 printf(
"Writing Voronoi edges.\n");
8072 if (*vedgelist == (
int *) NULL) {
8073 *vedgelist = (
int *) trimalloc((
int) (m->edges * 2 *
sizeof(int)));
8075 *vedgemarkerlist = (
int *) NULL;
8077 if (*vnormlist == (
float *) NULL) {
8078 *vnormlist = (
float *) trimalloc((
int) (m->edges * 2 *
sizeof(float)));
8081 normlist = *vnormlist;
8084 traversalinit(&m->triangles);
8085 triangleloop.tri = triangletraverse(m);
8086 vedgenumber = b->firstnumber;
8093 while (triangleloop.tri != (triangle *) NULL) {
8094 for (triangleloop.orient = 0; triangleloop.orient < 3;
8095 triangleloop.orient++) {
8096 sym(triangleloop, trisym);
8097 if ((triangleloop.tri < trisym.tri) || (trisym.tri == m->dummytri)) {
8099 p1 = * (
int *) (triangleloop.tri + 6);
8100 if (trisym.tri == m->dummytri) {
8101 org(triangleloop, torg);
8102 dest(triangleloop, tdest);
8104 elist[coordindex] = p1;
8105 normlist[coordindex++] = tdest[1] - torg[1];
8106 elist[coordindex] = -1;
8107 normlist[coordindex++] = torg[0] - tdest[0];
8110 p2 = * (
int *) (trisym.tri + 6);
8112 elist[coordindex] = p1;
8113 normlist[coordindex++] = 0.0;
8114 elist[coordindex] = p2;
8115 normlist[coordindex++] = 0.0;
8120 triangleloop.tri = triangletraverse(m);
8125 void writeneighbors(
struct mesh *m,
struct behavior *b,
int **neighborlist)
8129 struct otri triangleloop, trisym;
8131 int neighbor1, neighbor2, neighbor3;
8135 printf(
"Writing neighbors.\n");
8138 if (*neighborlist == (
int *) NULL) {
8139 *neighborlist = (
int *) trimalloc((
int) (m->triangles.items * 3 *
8142 nlist = *neighborlist;
8145 traversalinit(&m->triangles);
8146 triangleloop.tri = triangletraverse(m);
8147 triangleloop.orient = 0;
8148 elementnumber = b->firstnumber;
8149 while (triangleloop.tri != (triangle *) NULL) {
8150 * (
int *) (triangleloop.tri + 6) = (int) elementnumber;
8151 triangleloop.tri = triangletraverse(m);
8154 * (
int *) (m->dummytri + 6) = -1;
8156 traversalinit(&m->triangles);
8157 triangleloop.tri = triangletraverse(m);
8158 elementnumber = b->firstnumber;
8159 while (triangleloop.tri != (triangle *) NULL) {
8160 triangleloop.orient = 1;
8161 sym(triangleloop, trisym);
8162 neighbor1 = * (
int *) (trisym.tri + 6);
8163 triangleloop.orient = 2;
8164 sym(triangleloop, trisym);
8165 neighbor2 = * (
int *) (trisym.tri + 6);
8166 triangleloop.orient = 0;
8167 sym(triangleloop, trisym);
8168 neighbor3 = * (
int *) (trisym.tri + 6);
8169 nlist[index++] = neighbor1;
8170 nlist[index++] = neighbor2;
8171 nlist[index++] = neighbor3;
8173 triangleloop.tri = triangletraverse(m);
8188 void quality_statistics(
struct mesh *m,
struct behavior *b)
8190 struct otri triangleloop;
8192 float cossquaretable[8];
8193 float ratiotable[16];
8195 float edgelength[3];
8199 float shortest, longest;
8201 float smallestarea, biggestarea;
8202 float triminaltitude2;
8206 float smallestangle, biggestangle;
8207 float radconst, degconst;
8209 int aspecttable[16];
8215 printf(
"Mesh quality statistics:\n\n");
8216 radconst = PI / 18.0;
8217 degconst = 180.0 / PI;
8218 for (i = 0; i < 8; i++) {
8219 cossquaretable[i] = cos(radconst * (
float) (i + 1));
8220 cossquaretable[i] = cossquaretable[i] * cossquaretable[i];
8222 for (i = 0; i < 18; i++) {
8226 ratiotable[0] = 1.5; ratiotable[1] = 2.0;
8227 ratiotable[2] = 2.5; ratiotable[3] = 3.0;
8228 ratiotable[4] = 4.0; ratiotable[5] = 6.0;
8229 ratiotable[6] = 10.0; ratiotable[7] = 15.0;
8230 ratiotable[8] = 25.0; ratiotable[9] = 50.0;
8231 ratiotable[10] = 100.0; ratiotable[11] = 300.0;
8232 ratiotable[12] = 1000.0; ratiotable[13] = 10000.0;
8233 ratiotable[14] = 100000.0; ratiotable[15] = 0.0;
8234 for (i = 0; i < 16; i++) {
8239 minaltitude = m->xmax - m->xmin + m->ymax - m->ymin;
8240 minaltitude = minaltitude * minaltitude;
8241 shortest = minaltitude;
8243 smallestarea = minaltitude;
8246 smallestangle = 0.0;
8250 traversalinit(&m->triangles);
8251 triangleloop.tri = triangletraverse(m);
8252 triangleloop.orient = 0;
8253 while (triangleloop.tri != (triangle *) NULL) {
8254 org(triangleloop, p[0]);
8255 dest(triangleloop, p[1]);
8256 apex(triangleloop, p[2]);
8259 for (i = 0; i < 3; i++) {
8262 dx[i] = p[j][0] - p[k][0];
8263 dy[i] = p[j][1] - p[k][1];
8264 edgelength[i] = dx[i] * dx[i] + dy[i] * dy[i];
8265 if (edgelength[i] > trilongest2) {
8266 trilongest2 = edgelength[i];
8268 if (edgelength[i] > longest) {
8269 longest = edgelength[i];
8271 if (edgelength[i] < shortest) {
8272 shortest = edgelength[i];
8276 triarea = counterclockwise(m, b, p[0], p[1], p[2]);
8277 if (triarea < smallestarea) {
8278 smallestarea = triarea;
8280 if (triarea > biggestarea) {
8281 biggestarea = triarea;
8283 triminaltitude2 = triarea * triarea / trilongest2;
8284 if (triminaltitude2 < minaltitude) {
8285 minaltitude = triminaltitude2;
8287 triaspect2 = trilongest2 / triminaltitude2;
8288 if (triaspect2 > worstaspect) {
8289 worstaspect = triaspect2;
8292 while ((triaspect2 > ratiotable[aspectindex] * ratiotable[aspectindex])
8293 && (aspectindex < 15)) {
8296 aspecttable[aspectindex]++;
8298 for (i = 0; i < 3; i++) {
8301 dotproduct = dx[j] * dx[k] + dy[j] * dy[k];
8302 cossquare = dotproduct * dotproduct / (edgelength[j] * edgelength[k]);
8304 for (ii = 7; ii >= 0; ii--) {
8305 if (cossquare > cossquaretable[ii]) {
8309 if (dotproduct <= 0.0) {
8310 angletable[tendegree]++;
8311 if (cossquare > smallestangle) {
8312 smallestangle = cossquare;
8314 if (acutebiggest && (cossquare < biggestangle)) {
8315 biggestangle = cossquare;
8318 angletable[17 - tendegree]++;
8319 if (acutebiggest || (cossquare > biggestangle)) {
8320 biggestangle = cossquare;
8325 triangleloop.tri = triangletraverse(m);
8328 shortest = sqrt(shortest);
8329 longest = sqrt(longest);
8330 minaltitude = sqrt(minaltitude);
8331 worstaspect = sqrt(worstaspect);
8332 smallestarea *= 0.5;
8334 if (smallestangle >= 1.0) {
8335 smallestangle = 0.0;
8337 smallestangle = degconst * acos(sqrt(smallestangle));
8339 if (biggestangle >= 1.0) {
8340 biggestangle = 180.0;
8343 biggestangle = degconst * acos(sqrt(biggestangle));
8345 biggestangle = 180.0 - degconst * acos(sqrt(biggestangle));
8349 printf(
" Smallest area: %16.5g | Largest area: %16.5g\n",
8350 smallestarea, biggestarea);
8351 printf(
" Shortest edge: %16.5g | Longest edge: %16.5g\n",
8353 printf(
" Shortest altitude: %12.5g | Largest aspect ratio: %8.5g\n\n",
8354 minaltitude, worstaspect);
8356 printf(
" Triangle aspect ratio histogram:\n");
8357 printf(
" 1.1547 - %-6.6g : %8d | %6.6g - %-6.6g : %8d\n",
8358 ratiotable[0], aspecttable[0], ratiotable[7], ratiotable[8],
8360 for (i = 1; i < 7; i++) {
8361 printf(
" %6.6g - %-6.6g : %8d | %6.6g - %-6.6g : %8d\n",
8362 ratiotable[i - 1], ratiotable[i], aspecttable[i],
8363 ratiotable[i + 7], ratiotable[i + 8], aspecttable[i + 8]);
8365 printf(
" %6.6g - %-6.6g : %8d | %6.6g - : %8d\n",
8366 ratiotable[6], ratiotable[7], aspecttable[7], ratiotable[14],
8368 printf(
" (Aspect ratio is longest edge divided by shortest altitude)\n\n");
8370 printf(
" Smallest angle: %15.5g | Largest angle: %15.5g\n\n",
8371 smallestangle, biggestangle);
8373 printf(
" Angle histogram:\n");
8374 for (i = 0; i < 9; i++) {
8375 printf(
" %3d - %3d degrees: %8d | %3d - %3d degrees: %8d\n",
8376 i * 10, i * 10 + 10, angletable[i],
8377 i * 10 + 90, i * 10 + 100, angletable[i + 9]);
8388 void statistics(
struct mesh *m,
struct behavior *b)
8390 printf(
"\nStatistics:\n\n");
8391 printf(
" Input vertices: %d\n", m->invertices);
8393 printf(
" Input triangles: %d\n", m->inelements);
8396 printf(
" Input segments: %d\n", m->insegments);
8398 printf(
" Input holes: %d\n", m->holes);
8402 printf(
"\n Mesh vertices: %ld\n", m->vertices.items - m->undeads);
8403 printf(
" Mesh triangles: %ld\n", m->triangles.items);
8404 printf(
" Mesh edges: %ld\n", m->edges);
8405 printf(
" Mesh exterior boundary edges: %ld\n", m->hullsize);
8406 if (b->poly || b->refine) {
8407 printf(
" Mesh interior boundary edges: %ld\n",
8408 m->subsegs.items - m->hullsize);
8409 printf(
" Mesh subsegments (constrained edges): %ld\n",
8415 quality_statistics(m, b);
8416 printf(
"Memory allocation statistics:\n\n");
8417 printf(
" Maximum number of vertices: %ld\n", m->vertices.maxitems);
8418 printf(
" Maximum number of triangles: %ld\n", m->triangles.maxitems);
8419 if (m->subsegs.maxitems > 0) {
8420 printf(
" Maximum number of subsegments: %ld\n", m->subsegs.maxitems);
8422 if (m->viri.maxitems > 0) {
8423 printf(
" Maximum number of viri: %ld\n", m->viri.maxitems);
8425 if (m->badsubsegs.maxitems > 0) {
8426 printf(
" Maximum number of encroached subsegments: %ld\n",
8427 m->badsubsegs.maxitems);
8429 if (m->badtriangles.maxitems > 0) {
8430 printf(
" Maximum number of bad triangles: %ld\n",
8431 m->badtriangles.maxitems);
8433 if (m->flipstackers.maxitems > 0) {
8434 printf(
" Maximum number of stacked triangle flips: %ld\n",
8435 m->flipstackers.maxitems);
8437 if (m->splaynodes.maxitems > 0) {
8438 printf(
" Maximum number of splay tree nodes: %ld\n",
8439 m->splaynodes.maxitems);
8441 printf(
" Approximate heap memory use (bytes): %ld\n\n",
8442 m->vertices.maxitems * m->vertices.itembytes +
8443 m->triangles.maxitems * m->triangles.itembytes +
8444 m->subsegs.maxitems * m->subsegs.itembytes +
8445 m->viri.maxitems * m->viri.itembytes +
8446 m->badsubsegs.maxitems * m->badsubsegs.itembytes +
8447 m->badtriangles.maxitems * m->badtriangles.itembytes +
8448 m->flipstackers.maxitems * m->flipstackers.itembytes +
8449 m->splaynodes.maxitems * m->splaynodes.itembytes);
8451 printf(
"Algorithmic statistics:\n\n");
8453 printf(
" Number of incircle tests: %ld\n", m->incirclecount);
8455 printf(
" Number of 3D orientation tests: %ld\n", m->orient3dcount);
8457 printf(
" Number of 2D orientation tests: %ld\n", m->counterclockcount);
8458 if (m->hyperbolacount > 0) {
8459 printf(
" Number of right-of-hyperbola tests: %ld\n",
8462 if (m->circletopcount > 0) {
8463 printf(
" Number of circle top computations: %ld\n",
8466 if (m->circumcentercount > 0) {
8467 printf(
" Number of triangle circumcenter computations: %ld\n",
8468 m->circumcentercount);
8499 void triangulate(
char *triswitches,
struct triangulateio *in,
8500 struct triangulateio *out,
struct triangulateio *vorout)
8508 parsecommandline(1, &triswitches, &b);
8509 m.steinerleft = b.steiner;
8511 transfernodes(&m, &b, in->pointlist, in->pointattributelist,
8512 in->pointmarkerlist, in->numberofpoints,
8513 in->numberofpointattributes);
8515 m.hullsize = delaunay(&m, &b);
8518 m.infvertex1 = (vertex) NULL;
8519 m.infvertex2 = (vertex) NULL;
8520 m.infvertex3 = (vertex) NULL;
8522 if (b.usesegments) {
8523 m.checksegments = 1;
8526 formskeleton(&m, &b, in->segmentlist,
8527 in->segmentmarkerlist, in->numberofsegments);
8531 if (b.poly && (m.triangles.items > 0)) {
8532 holearray = in->holelist;
8533 m.holes = in->numberofholes;
8534 regionarray = in->regionlist;
8535 m.regions = in->numberofregions;
8538 carveholes(&m, &b, holearray, m.holes, regionarray, m.regions);
8549 m.edges = (3l * m.triangles.items + m.hullsize) / 2l;
8559 out->numberofpoints = m.vertices.items - m.undeads;
8561 out->numberofpoints = m.vertices.items;
8563 out->numberofpointattributes = m.nextras;
8564 out->numberoftriangles = m.triangles.items;
8565 out->numberofcorners = (b.order + 1) * (b.order + 2) / 2;
8566 out->numberoftriangleattributes = m.eextras;
8567 out->numberofedges = m.edges;
8568 if (b.usesegments) {
8569 out->numberofsegments = m.subsegs.items;
8571 out->numberofsegments = m.hullsize;
8573 if (vorout != (
struct triangulateio *) NULL) {
8574 vorout->numberofpoints = m.triangles.items;
8575 vorout->numberofpointattributes = m.nextras;
8576 vorout->numberofedges = m.edges;
8580 if (b.nonodewritten || (b.noiterationnum && m.readnodefile)) {
8582 printf(
"NOT writing vertices.\n");
8584 numbernodes(&m, &b);
8587 writenodes(&m, &b, &out->pointlist, &out->pointattributelist,
8588 &out->pointmarkerlist);
8590 if (b.noelewritten) {
8592 printf(
"NOT writing triangles.\n");
8595 writeelements(&m, &b, &out->trianglelist, &out->triangleattributelist);
8599 if (b.poly || b.convex) {
8601 if (b.nopolywritten || b.noiterationnum) {
8603 printf(
"NOT writing segments.\n");
8606 writepoly(&m, &b, &out->segmentlist, &out->segmentmarkerlist);
8607 out->numberofholes = m.holes;
8608 out->numberofregions = m.regions;
8610 out->holelist = in->holelist;
8611 out->regionlist = in->regionlist;
8613 out->holelist = (
float *) NULL;
8614 out->regionlist = (
float *) NULL;
8619 writeedges(&m, &b, &out->edgelist, &out->edgemarkerlist);
8622 writevoronoi(&m, &b, &vorout->pointlist, &vorout->pointattributelist,
8623 &vorout->pointmarkerlist, &vorout->edgelist,
8624 &vorout->edgemarkerlist, &vorout->normlist);
8627 writeneighbors(&m, &b, &out->neighborlist);
8634 triangledeinit(&m, &b);