343 namespace NekMeshUtils
381 double dxoa, dxda, dxod;
382 double dyoa, dyda, dyod;
383 double oalen, dalen, odlen;
386 dxoa = triorg[0] - triapex[0];
387 dyoa = triorg[1] - triapex[1];
388 dxda = tridest[0] - triapex[0];
389 dyda = tridest[1] - triapex[1];
390 dxod = triorg[0] - tridest[0];
391 dyod = triorg[1] - tridest[1];
393 oalen = dxoa * dxoa + dyoa * dyoa;
394 dalen = dxda * dxda + dyda * dyda;
395 odlen = dxod * dxod + dyod * dyod;
397 maxlen = (dalen > oalen) ? dalen : oalen;
398 maxlen = (odlen > maxlen) ? odlen : maxlen;
400 if (maxlen > 0.05 * (triorg[0] * triorg[0] + triorg[1] * triorg[1]) + 0.02)
429 memptr = (
void *)malloc((
unsigned int)size);
430 if (memptr == (
void *)NULL)
432 printf(
"Error: Out of memory.\n");
459 printf(
" Please report this bug to jrs@cs.berkeley.edu\n");
460 printf(
" Include the message above, your input data set, and the exact\n");
461 printf(
" command line you used to run Triangle.\n");
475 char workstring[2048];
483 for (i = 0; i < argc; i++)
485 for (j = 0; argv[i][j] !=
'\0'; j++)
487 if (argv[i][j] ==
'p')
491 if (argv[i][j] ==
'q')
494 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
495 (argv[i][j + 1] ==
'.'))
499 ((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
500 (argv[i][j + 1] ==
'.'))
503 workstring[k] = argv[i][j];
506 workstring[k] =
'\0';
507 b->
minangle = (double)strtod(workstring, (
char **)NULL);
514 if (argv[i][j] ==
'u')
519 if (argv[i][j] ==
'w')
523 if (argv[i][j] ==
'W')
527 if (argv[i][j] ==
'j')
531 if (argv[i][j] ==
'Y')
602 unsigned long alignptr;
610 alignptr = (
unsigned long)(pool->
nowblock + 1);
613 (alignptr % (
unsigned long)pool->
alignbytes));
649 if (alignment >
sizeof(
void *))
660 if (firstitemcount == 0)
706 unsigned long alignptr;
721 if (*(pool->
nowblock) == (
void *)NULL)
728 *(pool->
nowblock) = (
void *)newblock;
730 *newblock = (
void *)NULL;
737 alignptr = (
unsigned long)(pool->
nowblock + 1);
740 (
void *)(alignptr + (
unsigned long)pool->
alignbytes -
741 (alignptr % (
unsigned long)pool->
alignbytes));
783 unsigned long alignptr;
788 alignptr = (
unsigned long)(pool->
pathblock + 1);
791 (alignptr % (
unsigned long)pool->
alignbytes));
813 unsigned long alignptr;
828 alignptr = (
unsigned long)(pool->
pathblock + 1);
831 (alignptr % (
unsigned long)pool->
alignbytes));
876 unsigned long alignptr;
1023 8 *
sizeof(
triangle) +
sizeof(
int),
1065 if (newtriangle == (
triangle *)NULL)
1069 }
while (
deadtri(newtriangle));
1100 if (newsubseg == (
subseg *)NULL)
1135 if (newvertex == (
vertex)NULL)
1196 unsigned long alignptr;
1205 getblock = (
void **)*getblock;
1209 getblock = (
void **)*getblock;
1215 alignptr = (
unsigned long)(getblock + 1);
1283 for (i = 0; i < m->
eextras; i++)
1339 #define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
1355 #define Fast_Two_Sum_Tail(a, b, x, y) \
1359 #define Fast_Two_Sum(a, b, x, y) \
1360 x = (double)(a + b); \
1361 Fast_Two_Sum_Tail(a, b, x, y)
1363 #define Two_Sum_Tail(a, b, x, y) \
1364 bvirt = (double)(x - a); \
1365 avirt = x - bvirt; \
1366 bround = b - bvirt; \
1367 around = a - avirt; \
1370 #define Two_Sum(a, b, x, y) \
1371 x = (double)(a + b); \
1372 Two_Sum_Tail(a, b, x, y)
1374 #define Two_Diff_Tail(a, b, x, y) \
1375 bvirt = (double)(a - x); \
1376 avirt = x + bvirt; \
1377 bround = bvirt - b; \
1378 around = a - avirt; \
1381 #define Two_Diff(a, b, x, y) \
1382 x = (double)(a - b); \
1383 Two_Diff_Tail(a, b, x, y)
1385 #define Split(a, ahi, alo) \
1386 c = (double)(splitter * a); \
1387 abig = (double)(c - a); \
1391 #define Two_Product_Tail(a, b, x, y) \
1392 Split(a, ahi, alo); \
1393 Split(b, bhi, blo); \
1394 err1 = x - (ahi * bhi); \
1395 err2 = err1 - (alo * bhi); \
1396 err3 = err2 - (ahi * blo); \
1397 y = (alo * blo) - err3
1399 #define Two_Product(a, b, x, y) \
1400 x = (double)(a * b); \
1401 Two_Product_Tail(a, b, x, y)
1406 #define Two_Product_Presplit(a, b, bhi, blo, x, y) \
1407 x = (double)(a * b); \
1408 Split(a, ahi, alo); \
1409 err1 = x - (ahi * bhi); \
1410 err2 = err1 - (alo * bhi); \
1411 err3 = err2 - (ahi * blo); \
1412 y = (alo * blo) - err3
1416 #define Square_Tail(a, x, y) \
1417 Split(a, ahi, alo); \
1418 err1 = x - (ahi * ahi); \
1419 err3 = err1 - ((ahi + ahi) * alo); \
1420 y = (alo * alo) - err3
1422 #define Square(a, x, y) \
1423 x = (double)(a * a); \
1424 Square_Tail(a, x, y)
1429 #define Two_One_Sum(a1, a0, b, x2, x1, x0) \
1430 Two_Sum(a0, b, _i, x0); \
1431 Two_Sum(a1, _i, x2, x1)
1433 #define Two_One_Diff(a1, a0, b, x2, x1, x0) \
1434 Two_Diff(a0, b, _i, x0); \
1435 Two_Sum(a1, _i, x2, x1)
1437 #define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \
1438 Two_One_Sum(a1, a0, b0, _j, _0, x0); \
1439 Two_One_Sum(_j, _0, b1, x3, x2, x1)
1441 #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
1442 Two_One_Diff(a1, a0, b0, _j, _0, x0); \
1443 Two_One_Diff(_j, _0, b1, x3, x2, x1)
1447 #define Two_One_Product(a1, a0, b, x3, x2, x1, x0) \
1448 Split(b, bhi, blo); \
1449 Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
1450 Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
1451 Two_Sum(_i, _0, _k, x1); \
1452 Fast_Two_Sum(_j, _k, x3, x2)
1476 double check, lastcheck;
1496 every_other = !every_other;
1498 }
while ((check != 1.0) && (check != lastcheck));
1528 int elen,
double *e,
int flen,
double *f,
double *h)
1534 double avirt, bround, around;
1535 int eindex, findex, hindex;
1540 eindex = findex = 0;
1541 if ((fnow > enow) == (fnow > -enow))
1552 if ((eindex < elen) && (findex < flen))
1554 if ((fnow > enow) == (fnow > -enow))
1569 while ((eindex < elen) && (findex < flen))
1571 if ((fnow > enow) == (fnow > -enow))
1588 while (eindex < elen)
1598 while (findex < flen)
1608 if ((Q != 0.0) || (hindex == 0))
1639 double avirt, bround, around;
1642 double ahi, alo, bhi, blo;
1643 double err1, err2, err3;
1652 for (eindex = 1; eindex < elen; eindex++)
1656 Two_Sum(Q, product0, sum, hh);
1667 if ((Q != 0.0) || (hindex == 0))
1688 for (eindex = 1; eindex < elen; eindex++)
1717 double acx, acy, bcx, bcy;
1718 double acxtail, acytail, bcxtail, bcytail;
1719 double detleft, detright;
1720 double detlefttail, detrighttail;
1721 double det, errbound;
1722 double B[4], C1[8], C2[12], D[16];
1724 int C1length, C2length, Dlength;
1731 double avirt, bround, around;
1734 double ahi, alo, bhi, blo;
1735 double err1, err2, err3;
1739 acx = (double)(pa[0] - pc[0]);
1740 bcx = (double)(pb[0] - pc[0]);
1741 acy = (double)(pa[1] - pc[1]);
1742 bcy = (double)(pb[1] - pc[1]);
1748 detleft, detlefttail, detright, detrighttail, B3, B[2], B[1], B[0]);
1753 if ((det >= errbound) || (-det >= errbound))
1763 if ((acxtail == 0.0) && (acytail == 0.0) && (bcxtail == 0.0) &&
1770 det += (acx * bcytail + bcy * acxtail) - (acy * bcxtail + bcx * acytail);
1771 if ((det >= errbound) || (-det >= errbound))
1794 return (D[Dlength - 1]);
1800 double detleft, detright, det;
1801 double detsum, errbound;
1805 detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);
1806 detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);
1807 det = detleft - detright;
1811 if (detright <= 0.0)
1817 detsum = detleft + detright;
1820 else if (detleft < 0.0)
1822 if (detright >= 0.0)
1828 detsum = -detleft - detright;
1837 if ((det >= errbound) || (-det >= errbound))
1867 double adx, bdx, cdx, ady, bdy, cdy;
1868 double det, errbound;
1870 double bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
1871 double bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
1872 double bc[4], ca[4], ab[4];
1873 double bc3, ca3, ab3;
1874 double axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32];
1875 int axbclen, axxbclen, aybclen, ayybclen, alen;
1876 double bxca[8], bxxca[16], byca[8], byyca[16], bdet[32];
1877 int bxcalen, bxxcalen, bycalen, byycalen, blen;
1878 double cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32];
1879 int cxablen, cxxablen, cyablen, cyyablen, clen;
1882 double fin1[1152], fin2[1152];
1883 double *finnow, *finother, *finswap;
1886 double adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
1887 double adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
1888 double adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
1889 double aa[4], bb[4], cc[4];
1890 double aa3, bb3, cc3;
1895 double temp8[8], temp16a[16], temp16b[16], temp16c[16];
1896 double temp32a[32], temp32b[32], temp48[48], temp64[64];
1897 int temp8len, temp16alen, temp16blen, temp16clen;
1898 int temp32alen, temp32blen, temp48len, temp64len;
1899 double axtbb[8], axtcc[8], aytbb[8], aytcc[8];
1900 int axtbblen, axtcclen, aytbblen, aytcclen;
1901 double bxtaa[8], bxtcc[8], bytaa[8], bytcc[8];
1902 int bxtaalen, bxtcclen, bytaalen, bytcclen;
1903 double cxtaa[8], cxtbb[8], cytaa[8], cytbb[8];
1904 int cxtaalen, cxtbblen, cytaalen, cytbblen;
1905 double axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8];
1906 int axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen;
1907 double axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16],
1909 int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
1910 double axtbctt[8], aytbctt[8], bxtcatt[8];
1911 double bytcatt[8], cxtabtt[8], cytabtt[8];
1912 int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
1913 double abt[8], bct[8], cat[8];
1914 int abtlen, bctlen, catlen;
1915 double abtt[4], bctt[4], catt[4];
1916 int abttlen, bcttlen, cattlen;
1917 double abtt3, bctt3, catt3;
1921 double avirt, bround, around;
1924 double ahi, alo, bhi, blo;
1925 double err1, err2, err3;
1929 adx = (double)(pa[0] - pd[0]);
1930 bdx = (double)(pb[0] - pd[0]);
1931 cdx = (double)(pc[0] - pd[0]);
1932 ady = (double)(pa[1] - pd[1]);
1933 bdy = (double)(pb[1] - pd[1]);
1934 cdy = (double)(pc[1] - pd[1]);
1938 Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
1948 Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
1958 Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
1971 if ((det >= errbound) || (-det >= errbound))
1982 if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) &&
1983 (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0))
1989 det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail) -
1990 (bdy * cdxtail + cdx * bdytail)) +
1991 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx)) +
1992 ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail) -
1993 (cdy * adxtail + adx * cdytail)) +
1994 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx)) +
1995 ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail) -
1996 (ady * bdxtail + bdx * adytail)) +
1997 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
1998 if ((det >= errbound) || (-det >= errbound))
2006 if ((bdxtail != 0.0) || (bdytail != 0.0) || (cdxtail != 0.0) ||
2009 Square(adx, adxadx1, adxadx0);
2010 Square(ady, adyady1, adyady0);
2012 adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]);
2015 if ((cdxtail != 0.0) || (cdytail != 0.0) || (adxtail != 0.0) ||
2018 Square(bdx, bdxbdx1, bdxbdx0);
2019 Square(bdy, bdybdy1, bdybdy0);
2021 bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]);
2024 if ((adxtail != 0.0) || (adytail != 0.0) || (bdxtail != 0.0) ||
2027 Square(cdx, cdxcdx1, cdxcdx0);
2028 Square(cdy, cdycdy1, cdycdy0);
2030 cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]);
2047 temp16alen, temp16a, temp16blen, temp16b, temp32a);
2049 temp16clen, temp16c, temp32alen, temp32a, temp48);
2051 finlength, finnow, temp48len, temp48, finother);
2069 temp16alen, temp16a, temp16blen, temp16b, temp32a);
2071 temp16clen, temp16c, temp32alen, temp32a, temp48);
2073 finlength, finnow, temp48len, temp48, finother);
2091 temp16alen, temp16a, temp16blen, temp16b, temp32a);
2093 temp16clen, temp16c, temp32alen, temp32a, temp48);
2095 finlength, finnow, temp48len, temp48, finother);
2113 temp16alen, temp16a, temp16blen, temp16b, temp32a);
2115 temp16clen, temp16c, temp32alen, temp32a, temp48);
2117 finlength, finnow, temp48len, temp48, finother);
2135 temp16alen, temp16a, temp16blen, temp16b, temp32a);
2137 temp16clen, temp16c, temp32alen, temp32a, temp48);
2139 finlength, finnow, temp48len, temp48, finother);
2157 temp16alen, temp16a, temp16blen, temp16b, temp32a);
2159 temp16clen, temp16c, temp32alen, temp32a, temp48);
2161 finlength, finnow, temp48len, temp48, finother);
2167 if ((adxtail != 0.0) || (adytail != 0.0))
2169 if ((bdxtail != 0.0) || (bdytail != 0.0) || (cdxtail != 0.0) ||
2174 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
2180 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
2186 Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]);
2206 temp16alen, temp16a, temp32alen, temp32a, temp48);
2208 finlength, finnow, temp48len, temp48, finother);
2218 finlength, finnow, temp16alen, temp16a, finother);
2229 finlength, finnow, temp16alen, temp16a, finother);
2240 axtbcttlen, axtbctt, 2.0 * adx, temp16a);
2244 temp16alen, temp16a, temp16blen, temp16b, temp32b);
2246 temp32alen, temp32a, temp32blen, temp32b, temp64);
2248 finlength, finnow, temp64len, temp64, finother);
2261 temp16alen, temp16a, temp32alen, temp32a, temp48);
2263 finlength, finnow, temp48len, temp48, finother);
2273 aytbcttlen, aytbctt, 2.0 * ady, temp16a);
2277 temp16alen, temp16a, temp16blen, temp16b, temp32b);
2279 temp32alen, temp32a, temp32blen, temp32b, temp64);
2281 finlength, finnow, temp64len, temp64, finother);
2287 if ((bdxtail != 0.0) || (bdytail != 0.0))
2289 if ((cdxtail != 0.0) || (cdytail != 0.0) || (adxtail != 0.0) ||
2294 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
2300 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
2306 Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]);
2326 temp16alen, temp16a, temp32alen, temp32a, temp48);
2328 finlength, finnow, temp48len, temp48, finother);
2338 finlength, finnow, temp16alen, temp16a, finother);
2349 finlength, finnow, temp16alen, temp16a, finother);
2360 bxtcattlen, bxtcatt, 2.0 * bdx, temp16a);
2364 temp16alen, temp16a, temp16blen, temp16b, temp32b);
2366 temp32alen, temp32a, temp32blen, temp32b, temp64);
2368 finlength, finnow, temp64len, temp64, finother);
2381 temp16alen, temp16a, temp32alen, temp32a, temp48);
2383 finlength, finnow, temp48len, temp48, finother);
2393 bytcattlen, bytcatt, 2.0 * bdy, temp16a);
2397 temp16alen, temp16a, temp16blen, temp16b, temp32b);
2399 temp32alen, temp32a, temp32blen, temp32b, temp64);
2401 finlength, finnow, temp64len, temp64, finother);
2407 if ((cdxtail != 0.0) || (cdytail != 0.0))
2409 if ((adxtail != 0.0) || (adytail != 0.0) || (bdxtail != 0.0) ||
2414 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
2420 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
2426 Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]);
2446 temp16alen, temp16a, temp32alen, temp32a, temp48);
2448 finlength, finnow, temp48len, temp48, finother);
2458 finlength, finnow, temp16alen, temp16a, finother);
2469 finlength, finnow, temp16alen, temp16a, finother);
2480 cxtabttlen, cxtabtt, 2.0 * cdx, temp16a);
2484 temp16alen, temp16a, temp16blen, temp16b, temp32b);
2486 temp32alen, temp32a, temp32blen, temp32b, temp64);
2488 finlength, finnow, temp64len, temp64, finother);
2501 temp16alen, temp16a, temp32alen, temp32a, temp48);
2503 finlength, finnow, temp48len, temp48, finother);
2513 cytabttlen, cytabtt, 2.0 * cdy, temp16a);
2517 temp16alen, temp16a, temp16blen, temp16b, temp32b);
2519 temp32alen, temp32a, temp32blen, temp32b, temp64);
2521 finlength, finnow, temp64len, temp64, finother);
2528 return finnow[finlength - 1];
2538 double adx, bdx, cdx, ady, bdy, cdy;
2539 double bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
2540 double alift, blift, clift;
2542 double permanent, errbound;
2546 adx = pa[0] - pd[0];
2547 bdx = pb[0] - pd[0];
2548 cdx = pc[0] - pd[0];
2549 ady = pa[1] - pd[1];
2550 bdy = pb[1] - pd[1];
2551 cdy = pc[1] - pd[1];
2555 alift = adx * adx + ady * ady;
2559 blift = bdx * bdx + bdy * bdy;
2563 clift = cdx * cdx + cdy * cdy;
2565 det = alift * (bdxcdy - cdxbdy) + blift * (cdxady - adxcdy) +
2566 clift * (adxbdy - bdxady);
2572 if ((det > errbound) || (-det > errbound))
2612 double adx, bdx, cdx, ady, bdy, cdy, adheight, bdheight, cdheight;
2613 double det, errbound;
2615 double bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
2616 double bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
2617 double bc[4], ca[4], ab[4];
2618 double bc3, ca3, ab3;
2619 double adet[8], bdet[8], cdet[8];
2620 int alen, blen, clen;
2623 double *finnow, *finother, *finswap;
2624 double fin1[192], fin2[192];
2627 double adxtail, bdxtail, cdxtail;
2628 double adytail, bdytail, cdytail;
2629 double adheighttail, bdheighttail, cdheighttail;
2630 double at_blarge, at_clarge;
2631 double bt_clarge, bt_alarge;
2632 double ct_alarge, ct_blarge;
2633 double at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4];
2634 int at_blen, at_clen, bt_clen, bt_alen, ct_alen, ct_blen;
2635 double bdxt_cdy1, cdxt_bdy1, cdxt_ady1;
2636 double adxt_cdy1, adxt_bdy1, bdxt_ady1;
2637 double bdxt_cdy0, cdxt_bdy0, cdxt_ady0;
2638 double adxt_cdy0, adxt_bdy0, bdxt_ady0;
2639 double bdyt_cdx1, cdyt_bdx1, cdyt_adx1;
2640 double adyt_cdx1, adyt_bdx1, bdyt_adx1;
2641 double bdyt_cdx0, cdyt_bdx0, cdyt_adx0;
2642 double adyt_cdx0, adyt_bdx0, bdyt_adx0;
2643 double bct[8], cat[8], abt[8];
2644 int bctlen, catlen, abtlen;
2645 double bdxt_cdyt1, cdxt_bdyt1, cdxt_adyt1;
2646 double adxt_cdyt1, adxt_bdyt1, bdxt_adyt1;
2647 double bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0;
2648 double adxt_cdyt0, adxt_bdyt0, bdxt_adyt0;
2649 double u[4], v[12], w[16];
2651 int vlength, wlength;
2655 double avirt, bround, around;
2658 double ahi, alo, bhi, blo;
2659 double err1, err2, err3;
2663 adx = (double)(pa[0] - pd[0]);
2664 bdx = (double)(pb[0] - pd[0]);
2665 cdx = (double)(pc[0] - pd[0]);
2666 ady = (double)(pa[1] - pd[1]);
2667 bdy = (double)(pb[1] - pd[1]);
2668 cdy = (double)(pc[1] - pd[1]);
2669 adheight = (double)(aheight - dheight);
2670 bdheight = (double)(bheight - dheight);
2671 cdheight = (double)(cheight - dheight);
2675 Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
2681 Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
2687 Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
2696 if ((det >= errbound) || (-det >= errbound))
2711 if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) &&
2712 (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0) &&
2713 (adheighttail == 0.0) && (bdheighttail == 0.0) && (cdheighttail == 0.0))
2719 det += (adheight * ((bdx * cdytail + cdy * bdxtail) -
2720 (bdy * cdxtail + cdx * bdytail)) +
2721 adheighttail * (bdx * cdy - bdy * cdx)) +
2722 (bdheight * ((cdx * adytail + ady * cdxtail) -
2723 (cdy * adxtail + adx * cdytail)) +
2724 bdheighttail * (cdx * ady - cdy * adx)) +
2725 (cdheight * ((adx * bdytail + bdy * adxtail) -
2726 (ady * bdxtail + bdx * adytail)) +
2727 cdheighttail * (adx * bdy - ady * bdx));
2728 if ((det >= errbound) || (-det >= errbound))
2749 at_b[1] = at_blarge;
2752 at_c[1] = at_clarge;
2761 at_b[1] = at_blarge;
2765 at_c[1] = at_clarge;
2780 at_b[3] = at_blarge;
2792 at_c[3] = at_clarge;
2809 bt_c[1] = bt_clarge;
2812 bt_a[1] = bt_alarge;
2821 bt_c[1] = bt_clarge;
2825 bt_a[1] = bt_alarge;
2840 bt_c[3] = bt_clarge;
2852 bt_a[3] = bt_alarge;
2869 ct_a[1] = ct_alarge;
2872 ct_b[1] = ct_blarge;
2881 ct_a[1] = ct_alarge;
2885 ct_b[1] = ct_blarge;
2900 ct_a[3] = ct_alarge;
2912 ct_b[3] = ct_blarge;
2941 if (adheighttail != 0.0)
2945 finlength, finnow, vlength, v, finother);
2950 if (bdheighttail != 0.0)
2954 finlength, finnow, vlength, v, finother);
2959 if (cdheighttail != 0.0)
2963 finlength, finnow, vlength, v, finother);
2973 Two_Product(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0);
2975 adxt_bdyt1, adxt_bdyt0, cdheight, u3, u[2], u[1], u[0]);
2982 if (cdheighttail != 0.0)
2985 adxt_bdyt1, adxt_bdyt0, cdheighttail, u3, u[2], u[1], u[0]);
2988 finlength, finnow, 4, u, finother);
2997 Two_Product(negate, cdytail, adxt_cdyt1, adxt_cdyt0);
2999 adxt_cdyt1, adxt_cdyt0, bdheight, u3, u[2], u[1], u[0]);
3006 if (bdheighttail != 0.0)
3009 adxt_cdyt1, adxt_cdyt0, bdheighttail, u3, u[2], u[1], u[0]);
3012 finlength, finnow, 4, u, finother);
3023 Two_Product(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0);
3025 bdxt_cdyt1, bdxt_cdyt0, adheight, u3, u[2], u[1], u[0]);
3032 if (adheighttail != 0.0)
3035 bdxt_cdyt1, bdxt_cdyt0, adheighttail, u3, u[2], u[1], u[0]);
3038 finlength, finnow, 4, u, finother);
3047 Two_Product(negate, adytail, bdxt_adyt1, bdxt_adyt0);
3049 bdxt_adyt1, bdxt_adyt0, cdheight, u3, u[2], u[1], u[0]);
3056 if (cdheighttail != 0.0)
3059 bdxt_adyt1, bdxt_adyt0, cdheighttail, u3, u[2], u[1], u[0]);
3062 finlength, finnow, 4, u, finother);
3073 Two_Product(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0);
3075 cdxt_adyt1, cdxt_adyt0, bdheight, u3, u[2], u[1], u[0]);
3082 if (bdheighttail != 0.0)
3085 cdxt_adyt1, cdxt_adyt0, bdheighttail, u3, u[2], u[1], u[0]);
3088 finlength, finnow, 4, u, finother);
3097 Two_Product(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0);
3099 cdxt_bdyt1, cdxt_bdyt0, adheight, u3, u[2], u[1], u[0]);
3106 if (adheighttail != 0.0)
3109 cdxt_bdyt1, cdxt_bdyt0, adheighttail, u3, u[2], u[1], u[0]);
3112 finlength, finnow, 4, u, finother);
3120 if (adheighttail != 0.0)
3124 finlength, finnow, wlength, w, finother);
3129 if (bdheighttail != 0.0)
3133 finlength, finnow, wlength, w, finother);
3138 if (cdheighttail != 0.0)
3142 finlength, finnow, wlength, w, finother);
3148 return finnow[finlength - 1];
3162 double adx, bdx, cdx, ady, bdy, cdy, adheight, bdheight, cdheight;
3163 double bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
3165 double permanent, errbound;
3169 adx = pa[0] - pd[0];
3170 bdx = pb[0] - pd[0];
3171 cdx = pc[0] - pd[0];
3172 ady = pa[1] - pd[1];
3173 bdy = pb[1] - pd[1];
3174 cdy = pc[1] - pd[1];
3175 adheight = aheight - dheight;
3176 bdheight = bheight - dheight;
3177 cdheight = cheight - dheight;
3188 det = adheight * (bdxcdy - cdxbdy) + bdheight * (cdxady - adxcdy) +
3189 cdheight * (adxbdy - bdxady);
3195 if ((det > errbound) || (-det > errbound))
3201 pa, pb, pc, pd, aheight, bheight, cheight, dheight, permanent);
3231 return incircle(m, b, pa, pb, pc, pd);
3241 pa[0] * pa[0] + pa[1] * pa[1] - pa[2],
3242 pb[0] * pb[0] + pb[1] * pb[1] - pb[2],
3243 pc[0] * pc[0] + pc[1] * pc[1] - pc[2],
3244 pd[0] * pd[0] + pd[1] * pd[1] - pd[2]);
3248 return orient3d(m, b, pa, pb, pc, pd, pa[2], pb[2], pc[2], pd[2]);
3276 double xdo, ydo, xao, yao;
3277 double dodist, aodist, dadist;
3279 double dx, dy, dxoff, dyoff;
3284 xdo = tdest[0] - torg[0];
3285 ydo = tdest[1] - torg[1];
3286 xao = tapex[0] - torg[0];
3287 yao = tapex[1] - torg[1];
3288 dodist = xdo * xdo + ydo * ydo;
3289 aodist = xao * xao + yao * yao;
3290 dadist = (tdest[0] - tapex[0]) * (tdest[0] - tapex[0]) +
3291 (tdest[1] - tapex[1]) * (tdest[1] - tapex[1]);
3300 dx = (yao * dodist - ydo * aodist) * denominator;
3301 dy = (xdo * aodist - xao * dodist) * denominator;
3308 if ((dodist < aodist) && (dodist < dadist))
3318 if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy)
3325 else if (aodist < dadist)
3333 if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy)
3344 dxoff = 0.5 * (tapex[0] - tdest[0]) -
3346 dyoff = 0.5 * (tapex[1] - tdest[1]) +
3350 if (dxoff * dxoff + dyoff * dyoff <
3351 (dx - xdo) * (dx - xdo) + (dy - ydo) * (dy - ydo))
3359 circumcenter[0] = torg[0] + dx;
3360 circumcenter[1] = torg[1] + dy;
3367 *xi = (yao * dx - xao * dy) * (2.0 * denominator);
3368 *eta = (xdo * dy - ydo * dx) * (2.0 * denominator);
3435 double length, multiplier;
3436 int exponent, expincrement;
3443 if (badtri->
key >= 1.0)
3445 length = badtri->
key;
3452 length = 1.0 / badtri->
key;
3458 while (length > 2.0)
3463 while (length * multiplier * multiplier > 1.0)
3466 multiplier *= multiplier;
3469 exponent += expincrement;
3470 length *= multiplier;
3479 queuenumber = 2047 - exponent;
3483 queuenumber = 2048 + exponent;
3501 i = queuenumber + 1;
3536 struct otri *enqtri,
3547 newbad->
key = minedge;
3607 struct osub *testsubseg)
3609 struct otri neighbortri;
3610 struct osub testsym;
3615 vertex eorg, edest, eapex;
3621 sorg(*testsubseg, eorg);
3622 sdest(*testsubseg, edest);
3624 stpivot(*testsubseg, neighbortri);
3630 apex(neighbortri, eapex);
3636 dotproduct = (eorg[0] - eapex[0]) * (edest[0] - eapex[0]) +
3637 (eorg[1] - eapex[1]) * (edest[1] - eapex[1]);
3638 if (dotproduct < 0.0)
3640 if (dotproduct * dotproduct >=
3642 ((eorg[0] - eapex[0]) * (eorg[0] - eapex[0]) +
3643 (eorg[1] - eapex[1]) * (eorg[1] - eapex[1])) *
3644 ((edest[0] - eapex[0]) * (edest[0] - eapex[0]) +
3645 (edest[1] - eapex[1]) * (edest[1] - eapex[1])))
3652 ssym(*testsubseg, testsym);
3653 stpivot(testsym, neighbortri);
3659 apex(neighbortri, eapex);
3662 dotproduct = (eorg[0] - eapex[0]) * (edest[0] - eapex[0]) +
3663 (eorg[1] - eapex[1]) * (edest[1] - eapex[1]);
3664 if (dotproduct < 0.0)
3666 if (dotproduct * dotproduct >=
3668 ((eorg[0] - eapex[0]) * (eorg[0] - eapex[0]) +
3669 (eorg[1] - eapex[1]) * (eorg[1] - eapex[1])) *
3670 ((edest[0] - eapex[0]) * (edest[0] - eapex[0]) +
3671 (edest[1] - eapex[1]) * (edest[1] - eapex[1])))
3683 if (encroached == 1)
3712 struct otri tri1, tri2;
3713 struct osub testsub;
3714 vertex torg, tdest, tapex;
3716 vertex org1, dest1, org2, dest2;
3718 double dxod, dyod, dxda, dyda, dxao, dyao;
3719 double dxod2, dyod2, dxda2, dyda2, dxao2, dyao2;
3720 double apexlen, orglen, destlen, minedge;
3723 double dist1, dist2;
3727 org(*testtri, torg);
3728 dest(*testtri, tdest);
3729 apex(*testtri, tapex);
3730 dxod = torg[0] - tdest[0];
3731 dyod = torg[1] - tdest[1];
3732 dxda = tdest[0] - tapex[0];
3733 dyda = tdest[1] - tapex[1];
3734 dxao = tapex[0] - torg[0];
3735 dyao = tapex[1] - torg[1];
3736 dxod2 = dxod * dxod;
3737 dyod2 = dyod * dyod;
3738 dxda2 = dxda * dxda;
3739 dyda2 = dyda * dyda;
3740 dxao2 = dxao * dxao;
3741 dyao2 = dyao * dyao;
3743 apexlen = dxod2 + dyod2;
3744 orglen = dxda2 + dyda2;
3745 destlen = dxao2 + dyao2;
3747 if ((apexlen < orglen) && (apexlen < destlen))
3752 angle = dxda * dxao + dyda * dyao;
3753 angle = angle * angle / (orglen * destlen);
3758 else if (orglen < destlen)
3763 angle = dxod * dxao + dyod * dyao;
3764 angle = angle * angle / (apexlen * destlen);
3767 lnext(*testtri, tri1);
3774 angle = dxod * dxda + dyod * dyda;
3775 angle = angle * angle / (apexlen * orglen);
3778 lprev(*testtri, tri1);
3837 joinvertex = (
vertex)NULL;
3838 if ((dest1[0] == org2[0]) && (dest1[1] == org2[1]))
3842 else if ((org1[0] == dest2[0]) && (org1[1] == dest2[1]))
3846 if (joinvertex != (
vertex)NULL)
3852 dist1 = ((base1[0] - joinvertex[0]) *
3853 (base1[0] - joinvertex[0]) +
3854 (base1[1] - joinvertex[1]) *
3855 (base1[1] - joinvertex[1]));
3856 dist2 = ((base2[0] - joinvertex[0]) *
3857 (base2[0] - joinvertex[0]) +
3858 (base2[1] - joinvertex[1]) *
3859 (base2[1] - joinvertex[1]));
3862 if ((dist1 < 1.001 * dist2) && (dist1 > 0.999 * dist2))
3900 struct otri triangleloop;
3911 org(triangleloop, triorg);
3988 struct otri *searchtri,
3989 int stopatsubsegment)
3991 struct otri backtracktri;
3992 struct osub checkedge;
3993 vertex forg, fdest, fapex;
3994 double orgorient, destorient;
4000 org(*searchtri, forg);
4001 dest(*searchtri, fdest);
4002 apex(*searchtri, fapex);
4006 if ((fapex[0] == searchpoint[0]) && (fapex[1] == searchpoint[1]))
4017 if (destorient > 0.0)
4019 if (orgorient > 0.0)
4030 (fapex[0] - searchpoint[0]) * (fdest[0] - forg[0]) +
4031 (fapex[1] - searchpoint[1]) * (fdest[1] - forg[1]) >
4041 if (orgorient > 0.0)
4050 if (destorient == 0.0)
4055 if (orgorient == 0.0)
4069 lprev(*searchtri, backtracktri);
4074 lnext(*searchtri, backtracktri);
4077 sym(backtracktri, *searchtri);
4082 tspivot(backtracktri, checkedge);
4086 otricopy(backtracktri, *searchtri);
4094 otricopy(backtracktri, *searchtri);
4098 apex(*searchtri, fapex);
4141 struct otri *searchtri)
4145 struct otri sampletri;
4147 unsigned long alignptr;
4148 double searchdist, dist;
4150 long samplesperblock, totalsamplesleft, samplesleft;
4151 long population, totalpopulation;
4156 org(*searchtri, torg);
4157 searchdist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0]) +
4158 (searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
4167 if ((torg[0] == searchpoint[0]) && (torg[1] == searchpoint[1]))
4172 dist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0]) +
4173 (searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
4174 if (dist < searchdist)
4204 totalsamplesleft = m->
samples;
4209 while (totalsamplesleft > 0)
4212 if (population > totalpopulation)
4214 population = totalpopulation;
4217 alignptr = (
unsigned long)(sampleblock + 1);
4231 org(sampletri, torg);
4232 dist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0]) +
4233 (searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
4234 if (dist < searchdist)
4243 }
while ((samplesleft > 0) && (totalsamplesleft > 0));
4245 if (totalsamplesleft > 0)
4247 sampleblock = (
void **)*sampleblock;
4248 samplesleft = samplesperblock;
4249 totalpopulation -= population;
4255 org(*searchtri, torg);
4256 dest(*searchtri, tdest);
4258 if ((torg[0] == searchpoint[0]) && (torg[1] == searchpoint[1]))
4262 if ((tdest[0] == searchpoint[0]) && (tdest[1] == searchpoint[1]))
4276 else if (ahead == 0.0)
4279 if (((torg[0] < searchpoint[0]) == (searchpoint[0] < tdest[0])) &&
4280 ((torg[1] < searchpoint[1]) == (searchpoint[1] < tdest[1])))
4312 struct otri oppotri;
4313 struct osub newsubseg;
4319 dest(*tri, tridest);
4346 tsbond(oppotri, newsubseg);
4347 setmark(newsubseg, subsegmark);
4351 if (
mark(newsubseg) == 0)
4353 setmark(newsubseg, subsegmark);
4408 struct otri botleft, botright;
4409 struct otri topleft, topright;
4411 struct otri botlcasing, botrcasing;
4412 struct otri toplcasing, toprcasing;
4413 struct osub botlsubseg, botrsubseg;
4414 struct osub toplsubseg, toprsubseg;
4415 vertex leftvertex, rightvertex, botvertex;
4421 org(*flipedge, rightvertex);
4422 dest(*flipedge, leftvertex);
4423 apex(*flipedge, botvertex);
4424 sym(*flipedge, top);
4425 apex(top, farvertex);
4428 lprev(top, topleft);
4429 sym(topleft, toplcasing);
4430 lnext(top, topright);
4431 sym(topright, toprcasing);
4432 lnext(*flipedge, botleft);
4433 sym(botleft, botlcasing);
4434 lprev(*flipedge, botright);
4435 sym(botright, botrcasing);
4437 bond(topleft, botlcasing);
4438 bond(botleft, botrcasing);
4439 bond(botright, toprcasing);
4440 bond(topright, toplcasing);
4447 tspivot(botright, botrsubseg);
4448 tspivot(topright, toprsubseg);
4455 tsbond(topright, toplsubseg);
4463 tsbond(topleft, botlsubseg);
4471 tsbond(botleft, botrsubseg);
4479 tsbond(botright, toprsubseg);
4484 setorg(*flipedge, farvertex);
4485 setdest(*flipedge, botvertex);
4486 setapex(*flipedge, rightvertex);
4527 struct otri botleft, botright;
4528 struct otri topleft, topright;
4530 struct otri botlcasing, botrcasing;
4531 struct otri toplcasing, toprcasing;
4532 struct osub botlsubseg, botrsubseg;
4533 struct osub toplsubseg, toprsubseg;
4534 vertex leftvertex, rightvertex, botvertex;
4540 org(*flipedge, rightvertex);
4541 dest(*flipedge, leftvertex);
4542 apex(*flipedge, botvertex);
4543 sym(*flipedge, top);
4544 apex(top, farvertex);
4547 lprev(top, topleft);
4548 sym(topleft, toplcasing);
4549 lnext(top, topright);
4550 sym(topright, toprcasing);
4551 lnext(*flipedge, botleft);
4552 sym(botleft, botlcasing);
4553 lprev(*flipedge, botright);
4554 sym(botright, botrcasing);
4556 bond(topleft, toprcasing);
4557 bond(botleft, toplcasing);
4558 bond(botright, botlcasing);
4559 bond(topright, botrcasing);
4566 tspivot(botright, botrsubseg);
4567 tspivot(topright, toprsubseg);
4574 tsbond(botleft, toplsubseg);
4582 tsbond(botright, botlsubseg);
4590 tsbond(topright, botrsubseg);
4598 tsbond(topleft, toprsubseg);
4603 setorg(*flipedge, botvertex);
4604 setdest(*flipedge, farvertex);
4605 setapex(*flipedge, leftvertex);
4661 struct otri *searchtri,
4662 struct osub *splitseg,
4668 struct otri botleft, botright;
4669 struct otri topleft, topright;
4670 struct otri newbotleft, newbotright;
4671 struct otri newtopright;
4672 struct otri botlcasing, botrcasing;
4673 struct otri toplcasing, toprcasing;
4674 struct otri testtri;
4675 struct osub botlsubseg, botrsubseg;
4676 struct osub toplsubseg, toprsubseg;
4677 struct osub brokensubseg;
4678 struct osub checksubseg;
4679 struct osub rightsubseg;
4680 struct osub newsubseg;
4684 vertex leftvertex, rightvertex, botvertex, topvertex, farvertex;
4685 vertex segmentorg, segmentdest;
4696 if (splitseg == (
struct osub *)NULL)
4707 intersect =
locate(m, b, newvertex, &horiz);
4750 sym(horiz, testtri);
4775 lprev(horiz, botright);
4776 sym(botright, botrcasing);
4777 sym(horiz, topright);
4784 sym(topright, toprcasing);
4796 org(horiz, rightvertex);
4797 dest(horiz, leftvertex);
4798 apex(horiz, botvertex);
4799 setorg(newbotright, botvertex);
4800 setdest(newbotright, rightvertex);
4801 setapex(newbotright, newvertex);
4802 setorg(horiz, newvertex);
4803 for (i = 0; i < m->
eextras; i++)
4811 dest(topright, topvertex);
4812 setorg(newtopright, rightvertex);
4813 setdest(newtopright, topvertex);
4814 setapex(newtopright, newvertex);
4815 setorg(topright, newvertex);
4816 for (i = 0; i < m->
eextras; i++)
4827 tspivot(botright, botrsubseg);
4831 tsbond(newbotright, botrsubseg);
4835 tspivot(topright, toprsubseg);
4839 tsbond(newtopright, toprsubseg);
4845 bond(newbotright, botrcasing);
4847 bond(newbotright, botright);
4851 bond(newtopright, toprcasing);
4853 bond(newtopright, topright);
4855 bond(newtopright, newbotright);
4858 if (splitseg != (
struct osub *)NULL)
4862 segorg(*splitseg, segmentorg);
4863 segdest(*splitseg, segmentdest);
4865 spivot(*splitseg, rightsubseg);
4867 tspivot(newbotright, newsubseg);
4870 sbond(*splitseg, newsubseg);
4872 sbond(newsubseg, rightsubseg);
4897 lnext(horiz, botleft);
4898 lprev(horiz, botright);
4899 sym(botleft, botlcasing);
4900 sym(botright, botrcasing);
4905 org(horiz, rightvertex);
4906 dest(horiz, leftvertex);
4907 apex(horiz, botvertex);
4908 setorg(newbotleft, leftvertex);
4909 setdest(newbotleft, botvertex);
4910 setapex(newbotleft, newvertex);
4911 setorg(newbotright, botvertex);
4912 setdest(newbotright, rightvertex);
4913 setapex(newbotright, newvertex);
4915 for (i = 0; i < m->
eextras; i++)
4931 tsbond(newbotleft, botlsubseg);
4933 tspivot(botright, botrsubseg);
4937 tsbond(newbotright, botrsubseg);
4942 bond(newbotleft, botlcasing);
4943 bond(newbotright, botrcasing);
4946 bond(newbotleft, newbotright);
4948 bond(botleft, newbotleft);
4950 bond(botright, newbotright);
4970 rightvertex = first;
4971 dest(horiz, leftvertex);
5009 apex(top, farvertex);
5030 m, b, newvertex, rightvertex, farvertex) > 0.0;
5044 m, b, farvertex, leftvertex, newvertex) > 0.0;
5071 lprev(top, topleft);
5072 sym(topleft, toplcasing);
5073 lnext(top, topright);
5074 sym(topright, toprcasing);
5075 lnext(horiz, botleft);
5076 sym(botleft, botlcasing);
5077 lprev(horiz, botright);
5078 sym(botright, botrcasing);
5081 bond(topleft, botlcasing);
5082 bond(botleft, botrcasing);
5083 bond(botright, toprcasing);
5084 bond(topright, toplcasing);
5091 tspivot(botright, botrsubseg);
5092 tspivot(topright, toprsubseg);
5099 tsbond(topright, toplsubseg);
5107 tsbond(topleft, botlsubseg);
5115 tsbond(botleft, botrsubseg);
5123 tsbond(botright, toprsubseg);
5127 setorg(horiz, farvertex);
5133 for (i = 0; i < m->
eextras; i++)
5157 leftvertex = farvertex;
5171 sym(horiz, testtri);
5177 if ((leftvertex == first) || (testtri.
tri == m->
dummytri))
5181 lnext(horiz, *searchtri);
5186 lnext(testtri, horiz);
5187 rightvertex = leftvertex;
5188 dest(horiz, leftvertex);
5259 struct otri *firstedge,
5260 struct otri *lastedge,
5265 struct otri testtri;
5266 struct otri besttri;
5267 struct otri tempedge;
5268 vertex leftbasevertex, rightbasevertex;
5276 apex(*lastedge, leftbasevertex);
5277 dest(*firstedge, rightbasevertex);
5280 onext(*firstedge, besttri);
5281 dest(besttri, bestvertex);
5284 for (i = 2; i <= edgecount - 2; i++)
5287 dest(testtri, testvertex);
5290 m, b, leftbasevertex, rightbasevertex, bestvertex, testvertex) >
5294 bestvertex = testvertex;
5302 oprev(besttri, tempedge);
5304 m, b, firstedge, &tempedge, bestnumber + 1, 1, triflaws);
5306 if (bestnumber < edgecount - 2)
5309 sym(besttri, tempedge);
5311 m, b, &besttri, lastedge, edgecount - bestnumber, 1, triflaws);
5313 sym(tempedge, besttri);
5318 flip(m, b, &besttri);
5322 sym(besttri, testtri);
5346 struct otri countingtri;
5347 struct otri firstedge, lastedge;
5348 struct otri deltriright;
5349 struct otri lefttri, righttri;
5350 struct otri leftcasing, rightcasing;
5351 struct osub leftsubseg, rightsubseg;
5358 org(*deltri, delvertex);
5363 onext(*deltri, countingtri);
5365 while (!
otriequal(*deltri, countingtri))
5376 onext(*deltri, firstedge);
5377 oprev(*deltri, lastedge);
5379 m, b, &firstedge, &lastedge, edgecount, 0, !b->
nobisect);
5382 lprev(*deltri, deltriright);
5383 dnext(*deltri, lefttri);
5384 sym(lefttri, leftcasing);
5385 oprev(deltriright, righttri);
5386 sym(righttri, rightcasing);
5387 bond(*deltri, leftcasing);
5388 bond(deltriright, rightcasing);
5392 tsbond(*deltri, leftsubseg);
5394 tspivot(righttri, rightsubseg);
5397 tsbond(deltriright, rightsubseg);
5401 org(lefttri, neworg);
5426 struct otri fliptri;
5427 struct otri botleft, botright, topright;
5428 struct otri botlcasing, botrcasing, toprcasing;
5429 struct otri gluetri;
5430 struct osub botlsubseg, botrsubseg, toprsubseg;
5431 vertex botvertex, rightvertex;
5453 dprev(fliptri, botleft);
5455 onext(fliptri, botright);
5457 sym(botleft, botlcasing);
5458 sym(botright, botrcasing);
5459 dest(botleft, botvertex);
5463 bond(fliptri, botlcasing);
5465 tsbond(fliptri, botlsubseg);
5467 bond(fliptri, botrcasing);
5468 tspivot(botright, botrsubseg);
5469 tsbond(fliptri, botrsubseg);
5479 lprev(fliptri, gluetri);
5480 sym(gluetri, botright);
5482 sym(botright, botrcasing);
5483 dest(botright, rightvertex);
5485 setorg(fliptri, rightvertex);
5486 bond(gluetri, botrcasing);
5487 tspivot(botright, botrsubseg);
5488 tsbond(gluetri, botrsubseg);
5493 sym(fliptri, gluetri);
5497 dnext(gluetri, topright);
5498 sym(topright, toprcasing);
5500 setorg(gluetri, rightvertex);
5501 bond(gluetri, toprcasing);
5502 tspivot(topright, toprsubseg);
5503 tsbond(gluetri, toprsubseg);
5576 double pivotx, pivoty;
5582 if ((sortarray[0][0] > sortarray[1][0]) ||
5583 ((sortarray[0][0] == sortarray[1][0]) &&
5584 (sortarray[0][1] > sortarray[1][1])))
5586 temp = sortarray[1];
5587 sortarray[1] = sortarray[0];
5588 sortarray[0] = temp;
5594 pivotx = sortarray[pivot][0];
5595 pivoty = sortarray[pivot][1];
5599 while (left < right)
5605 }
while ((left <= right) && ((sortarray[left][0] < pivotx) ||
5606 ((sortarray[left][0] == pivotx) &&
5607 (sortarray[left][1] < pivoty))));
5612 }
while ((left <= right) && ((sortarray[right][0] > pivotx) ||
5613 ((sortarray[right][0] == pivotx) &&
5614 (sortarray[right][1] > pivoty))));
5618 temp = sortarray[left];
5619 sortarray[left] = sortarray[right];
5620 sortarray[right] = temp;
5628 if (right < arraysize - 2)
5631 vertexsort(&sortarray[right + 1], arraysize - right - 1);
5651 double pivot1, pivot2;
5657 if ((sortarray[0][axis] > sortarray[1][axis]) ||
5658 ((sortarray[0][axis] == sortarray[1][axis]) &&
5659 (sortarray[0][1 - axis] > sortarray[1][1 - axis])))
5661 temp = sortarray[1];
5662 sortarray[1] = sortarray[0];
5663 sortarray[0] = temp;
5669 pivot1 = sortarray[pivot][axis];
5670 pivot2 = sortarray[pivot][1 - axis];
5674 while (left < right)
5680 }
while ((left <= right) && ((sortarray[left][axis] < pivot1) ||
5681 ((sortarray[left][axis] == pivot1) &&
5682 (sortarray[left][1 - axis] < pivot2))));
5687 }
while ((left <= right) && ((sortarray[right][axis] > pivot1) ||
5688 ((sortarray[right][axis] == pivot1) &&
5689 (sortarray[right][1 - axis] > pivot2))));
5693 temp = sortarray[left];
5694 sortarray[left] = sortarray[right];
5695 sortarray[right] = temp;
5705 if (right < median - 1)
5709 arraysize - right - 1,
5730 divider = arraysize >> 1;
5740 if (arraysize - divider >= 2)
5746 alternateaxes(&sortarray[divider], arraysize - divider, 1 - axis);
5787 struct otri *farleft,
5788 struct otri *innerleft,
5789 struct otri *innerright,
5790 struct otri *farright,
5793 struct otri leftcand, rightcand;
5794 struct otri baseedge;
5795 struct otri nextedge;
5796 struct otri sidecasing, topcasing, outercasing;
5797 struct otri checkedge;
5800 vertex innerleftapex, innerrightapex;
5801 vertex farleftpt, farrightpt;
5802 vertex farleftapex, farrightapex;
5803 vertex lowerleft, lowerright;
5804 vertex upperleft, upperright;
5809 int leftfinished, rightfinished;
5812 dest(*innerleft, innerleftdest);
5813 apex(*innerleft, innerleftapex);
5814 org(*innerright, innerrightorg);
5815 apex(*innerright, innerrightapex);
5819 org(*farleft, farleftpt);
5820 apex(*farleft, farleftapex);
5821 dest(*farright, farrightpt);
5822 apex(*farright, farrightapex);
5826 while (farleftapex[1] < farleftpt[1])
5830 farleftpt = farleftapex;
5831 apex(*farleft, farleftapex);
5833 sym(*innerleft, checkedge);
5834 apex(checkedge, checkvertex);
5835 while (checkvertex[1] > innerleftdest[1])
5837 lnext(checkedge, *innerleft);
5838 innerleftapex = innerleftdest;
5839 innerleftdest = checkvertex;
5840 sym(*innerleft, checkedge);
5841 apex(checkedge, checkvertex);
5843 while (innerrightapex[1] < innerrightorg[1])
5847 innerrightorg = innerrightapex;
5848 apex(*innerright, innerrightapex);
5850 sym(*farright, checkedge);
5851 apex(checkedge, checkvertex);
5852 while (checkvertex[1] > farrightpt[1])
5854 lnext(checkedge, *farright);
5855 farrightapex = farrightpt;
5856 farrightpt = checkvertex;
5857 sym(*farright, checkedge);
5858 apex(checkedge, checkvertex);
5867 m, b, innerleftdest, innerleftapex, innerrightorg) > 0.0)
5871 innerleftdest = innerleftapex;
5872 apex(*innerleft, innerleftapex);
5877 m, b, innerrightapex, innerrightorg, innerleftdest) > 0.0)
5881 innerrightorg = innerrightapex;
5882 apex(*innerright, innerrightapex);
5885 }
while (changemade);
5887 sym(*innerleft, leftcand);
5888 sym(*innerright, rightcand);
5892 bond(baseedge, *innerleft);
5894 bond(baseedge, *innerright);
5896 setorg(baseedge, innerrightorg);
5897 setdest(baseedge, innerleftdest);
5901 org(*farleft, farleftpt);
5902 if (innerleftdest == farleftpt)
5904 lnext(baseedge, *farleft);
5906 dest(*farright, farrightpt);
5907 if (innerrightorg == farrightpt)
5909 lprev(baseedge, *farright);
5912 lowerleft = innerleftdest;
5913 lowerright = innerrightorg;
5915 apex(leftcand, upperleft);
5916 apex(rightcand, upperright);
5931 if (leftfinished && rightfinished)
5935 setorg(nextedge, lowerleft);
5936 setdest(nextedge, lowerright);
5939 bond(nextedge, baseedge);
5941 bond(nextedge, rightcand);
5943 bond(nextedge, leftcand);
5948 org(*farleft, farleftpt);
5949 apex(*farleft, farleftapex);
5950 dest(*farright, farrightpt);
5951 apex(*farright, farrightapex);
5952 sym(*farleft, checkedge);
5953 apex(checkedge, checkvertex);
5957 while (checkvertex[0] < farleftpt[0])
5959 lprev(checkedge, *farleft);
5960 farleftapex = farleftpt;
5961 farleftpt = checkvertex;
5962 sym(*farleft, checkedge);
5963 apex(checkedge, checkvertex);
5965 while (farrightapex[0] > farrightpt[0])
5969 farrightpt = farrightapex;
5970 apex(*farright, farrightapex);
5979 lprev(leftcand, nextedge);
5981 apex(nextedge, nextapex);
5984 if (nextapex != (
vertex)NULL)
5988 incircle(m, b, lowerleft, lowerright, upperleft, nextapex) >
5997 sym(nextedge, topcasing);
5999 sym(nextedge, sidecasing);
6000 bond(nextedge, topcasing);
6001 bond(leftcand, sidecasing);
6003 sym(leftcand, outercasing);
6005 bond(nextedge, outercasing);
6007 setorg(leftcand, lowerleft);
6014 upperleft = nextapex;
6018 apex(nextedge, nextapex);
6019 if (nextapex != (
vertex)NULL)
6041 lnext(rightcand, nextedge);
6043 apex(nextedge, nextapex);
6046 if (nextapex != (
vertex)NULL)
6051 m, b, lowerleft, lowerright, upperright, nextapex) >
6060 sym(nextedge, topcasing);
6062 sym(nextedge, sidecasing);
6063 bond(nextedge, topcasing);
6064 bond(rightcand, sidecasing);
6066 sym(rightcand, outercasing);
6068 bond(nextedge, outercasing);
6071 setdest(rightcand, lowerright);
6073 setorg(nextedge, upperright);
6077 upperright = nextapex;
6081 apex(nextedge, nextapex);
6082 if (nextapex != (
vertex)NULL)
6102 (
incircle(m, b, upperleft, lowerleft, lowerright, upperright) >
6107 bond(baseedge, rightcand);
6108 lprev(rightcand, baseedge);
6110 lowerright = upperright;
6111 sym(baseedge, rightcand);
6112 apex(rightcand, upperright);
6118 bond(baseedge, leftcand);
6119 lnext(leftcand, baseedge);
6120 setorg(baseedge, lowerright);
6121 lowerleft = upperleft;
6122 sym(baseedge, leftcand);
6123 apex(leftcand, upperleft);
6150 struct otri *farleft,
6151 struct otri *farright)
6153 struct otri midtri, tri1, tri2, tri3;
6154 struct otri innerleft, innerright;
6163 setorg(*farleft, sortarray[0]);
6164 setdest(*farleft, sortarray[1]);
6167 setorg(*farright, sortarray[1]);
6168 setdest(*farright, sortarray[0]);
6170 bond(*farleft, *farright);
6173 bond(*farleft, *farright);
6176 bond(*farleft, *farright);
6179 lprev(*farright, *farleft);
6182 else if (vertices == 3)
6195 setorg(midtri, sortarray[0]);
6196 setdest(midtri, sortarray[1]);
6197 setorg(tri1, sortarray[1]);
6199 setorg(tri2, sortarray[2]);
6201 setorg(tri3, sortarray[1]);
6227 setorg(midtri, sortarray[0]);
6229 setorg(tri3, sortarray[0]);
6234 setdest(midtri, sortarray[1]);
6235 setorg(tri1, sortarray[1]);
6237 setapex(midtri, sortarray[2]);
6238 setorg(tri2, sortarray[2]);
6244 setdest(midtri, sortarray[2]);
6245 setorg(tri1, sortarray[2]);
6247 setapex(midtri, sortarray[1]);
6248 setorg(tri2, sortarray[1]);
6275 lnext(*farleft, *farright);
6284 divider = vertices >> 1;
6286 divconqrecurse(m, b, sortarray, divider, 1 - axis, farleft, &innerleft);
6289 &sortarray[divider],
6296 mergehulls(m, b, farleft, &innerleft, &innerright, farright, axis);
6302 struct otri searchedge;
6303 struct otri dissolveedge;
6304 struct otri deadtriangle;
6311 lprev(*startghost, searchedge);
6315 otricopy(*startghost, dissolveedge);
6320 lnext(dissolveedge, deadtriangle);
6332 org(dissolveedge, markorg);
6342 sym(deadtriangle, dissolveedge);
6345 }
while (!
otriequal(dissolveedge, *startghost));
6362 struct otri hullleft, hullright;
6379 if ((sortarray[i][0] == sortarray[j][0]) &&
6380 (sortarray[i][1] == sortarray[j][1]))
6388 sortarray[i] = sortarray[j];
6395 if (i - divider >= 2)
6472 struct otri *searchtri,
6475 struct otri checktri;
6477 vertex leftvertex, rightvertex;
6478 double leftccw, rightccw;
6479 int leftflag, rightflag;
6482 org(*searchtri, startvertex);
6483 dest(*searchtri, rightvertex);
6484 apex(*searchtri, leftvertex);
6487 leftflag = leftccw > 0.0;
6490 rightflag = rightccw > 0.0;
6491 if (leftflag && rightflag)
6497 onext(*searchtri, checktri);
6513 printf(
"Internal error in finddirection(): Unable to find a\n");
6514 printf(
" triangle leading from (%.12g, %.12g) to",
6517 printf(
" (%.12g, %.12g).\n", searchpoint[0], searchpoint[1]);
6520 apex(*searchtri, leftvertex);
6523 leftflag = leftccw > 0.0;
6531 printf(
"Internal error in finddirection(): Unable to find a\n");
6532 printf(
" triangle leading from (%.12g, %.12g) to",
6535 printf(
" (%.12g, %.12g).\n", searchpoint[0], searchpoint[1]);
6538 dest(*searchtri, rightvertex);
6542 rightflag = rightccw > 0.0;
6548 else if (rightccw == 0.0)
6577 struct otri *splittri,
6578 struct osub *splitsubseg,
6581 struct osub opposubseg;
6584 vertex leftvertex, rightvertex;
6590 double split, denom;
6596 apex(*splittri, endpoint1);
6597 org(*splittri, torg);
6598 dest(*splittri, tdest);
6600 tx = tdest[0] - torg[0];
6601 ty = tdest[1] - torg[1];
6602 ex = endpoint2[0] - endpoint1[0];
6603 ey = endpoint2[1] - endpoint1[1];
6604 etx = torg[0] - endpoint2[0];
6605 ety = torg[1] - endpoint2[1];
6606 denom = ty * ex - tx * ey;
6609 printf(
"Internal error in segmentintersection():");
6610 printf(
" Attempt to find intersection of parallel segments.\n");
6613 split = (ey * etx - ex * ety) / denom;
6617 for (i = 0; i < 2 + m->
nextras; i++)
6619 newvertex[i] = torg[i] + split * (tdest[i] - torg[i]);
6625 success =
insertvertex(m, b, newvertex, splittri, splitsubseg, 0, 0);
6628 printf(
"Internal error in segmentintersection():\n");
6629 printf(
" Failure to split a segment.\n");
6641 spivot(*splitsubseg, opposubseg);
6659 dest(*splittri, rightvertex);
6660 apex(*splittri, leftvertex);
6661 if ((leftvertex[0] == endpoint1[0]) && (leftvertex[1] == endpoint1[1]))
6665 else if ((rightvertex[0] != endpoint1[0]) ||
6666 (rightvertex[1] != endpoint1[1]))
6668 printf(
"Internal error in segmentintersection():\n");
6669 printf(
" Topological inconsistency after splitting a segment.\n");
6701 struct otri *searchtri,
6705 struct otri crosstri;
6706 struct osub crosssubseg;
6707 vertex leftvertex, rightvertex;
6712 dest(*searchtri, rightvertex);
6713 apex(*searchtri, leftvertex);
6714 if (((leftvertex[0] == endpoint2[0]) && (leftvertex[1] == endpoint2[1])) ||
6715 ((rightvertex[0] == endpoint2[0]) && (rightvertex[1] == endpoint2[1])))
6718 if ((leftvertex[0] == endpoint2[0]) && (leftvertex[1] == endpoint2[1]))
6733 return scoutsegment(m, b, searchtri, endpoint2, newmark);
6742 return scoutsegment(m, b, searchtri, endpoint2, newmark);
6746 lnext(*searchtri, crosstri);
6747 tspivot(crosstri, crosssubseg);
6760 return scoutsegment(m, b, searchtri, endpoint2, newmark);
6790 struct otri searchtri1, searchtri2;
6791 struct osub brokensubseg;
6793 vertex midvertex1, midvertex2;
6801 for (i = 0; i < 2 + m->
nextras; i++)
6803 newvertex[i] = 0.5 * (endpoint1[i] + endpoint2[i]);
6817 org(searchtri1, newvertex);
6825 tspivot(searchtri1, brokensubseg);
6827 insertvertex(m, b, newvertex, &searchtri1, &brokensubseg, 0, 0);
6830 printf(
"Internal error in conformingedge():\n");
6831 printf(
" Failure to split a segment.\n");
6848 if (!
scoutsegment(m, b, &searchtri1, endpoint1, newmark))
6852 org(searchtri1, midvertex1);
6855 if (!
scoutsegment(m, b, &searchtri2, endpoint2, newmark))
6859 org(searchtri2, midvertex2);
6904 struct otri *fixuptri,
6907 struct otri neartri;
6909 struct osub faredge;
6910 vertex nearvertex, leftvertex, rightvertex, farvertex;
6914 lnext(*fixuptri, neartri);
6915 sym(neartri, fartri);
6927 apex(neartri, nearvertex);
6928 org(neartri, leftvertex);
6929 dest(neartri, rightvertex);
6930 apex(fartri, farvertex);
6956 if (
incircle(m, b, leftvertex, farvertex, rightvertex, nearvertex) <=
6963 flip(m, b, &neartri);
7026 struct otri *starttri,
7030 struct otri fixuptri, fixuptri2;
7031 struct osub crosssubseg;
7040 org(*starttri, endpoint1);
7041 lnext(*starttri, fixuptri);
7042 flip(m, b, &fixuptri);
7049 org(fixuptri, farvertex);
7052 if ((farvertex[0] == endpoint2[0]) && (farvertex[1] == endpoint2[1]))
7054 oprev(fixuptri, fixuptri2);
7071 oprev(fixuptri, fixuptri2);
7081 oprev(fixuptri, fixuptri2);
7104 tspivot(fixuptri, crosssubseg);
7118 m, b, &fixuptri, &crosssubseg, endpoint2);
7131 if (!
scoutsegment(m, b, &fixuptri, endpoint2, newmark))
7150 struct otri searchtri1, searchtri2;
7158 checkvertex = (
vertex)NULL;
7162 decode(encodedtri, searchtri1);
7163 org(searchtri1, checkvertex);
7165 if (checkvertex != endpoint1)
7174 printf(
"Internal error in insertsegment(): Unable to locate PSLG "
7176 printf(
" (%.12g, %.12g) in triangulation.\n",
7186 if (
scoutsegment(m, b, &searchtri1, endpoint2, newmark))
7193 org(searchtri1, endpoint1);
7196 checkvertex = (
vertex)NULL;
7200 decode(encodedtri, searchtri2);
7201 org(searchtri2, checkvertex);
7203 if (checkvertex != endpoint2)
7212 printf(
"Internal error in insertsegment(): Unable to locate PSLG "
7214 printf(
" (%.12g, %.12g) in triangulation.\n",
7224 if (
scoutsegment(m, b, &searchtri2, endpoint1, newmark))
7232 org(searchtri2, endpoint2);
7247 struct otri hulltri;
7248 struct otri nexttri;
7249 struct otri starttri;
7265 oprev(hulltri, nexttri);
7269 oprev(hulltri, nexttri);
7271 }
while (!
otriequal(hulltri, starttri));
7287 int *segmentmarkerlist,
7288 int numberofsegments)
7290 char polyfilename[6];
7292 vertex endpoint1, endpoint2;
7300 strcpy(polyfilename,
"input");
7302 segmentmarkers = segmentmarkerlist != (
int *)NULL;
7323 end1 = segmentlist[index++];
7324 end2 = segmentlist[index++];
7327 boundmarker = segmentmarkerlist[i];
7334 else if ((end2 < 0) ||
7344 if ((endpoint1[0] == endpoint2[0]) &&
7345 (endpoint1[1] == endpoint2[1]))
7386 struct otri hulltri;
7387 struct otri nexttri;
7388 struct otri starttri;
7389 struct osub hullsubseg;
7416 *deadtriangle = hulltri.
tri;
7423 if (
mark(hullsubseg) == 0)
7427 dest(hulltri, hdest);
7441 oprev(hulltri, nexttri);
7445 oprev(hulltri, nexttri);
7447 }
while (!
otriequal(hulltri, starttri));
7469 struct otri testtri;
7470 struct otri neighbor;
7473 struct osub neighborsubseg;
7484 while (virusloop != (
triangle **)NULL)
7486 testtri.
tri = *virusloop;
7499 sym(testtri, neighbor);
7501 tspivot(testtri, neighborsubseg);
7533 *deadtriangle = neighbor.
tri;
7541 if (
mark(neighborsubseg) == 0)
7545 org(neighbor, norg);
7546 dest(neighbor, ndest);
7567 while (virusloop != (
triangle **)NULL)
7569 testtri.
tri = *virusloop;
7576 org(testtri, testvertex);
7578 if (testvertex != (
vertex)NULL)
7584 onext(testtri, neighbor);
7607 oprev(testtri, neighbor);
7639 sym(testtri, neighbor);
7686 struct otri testtri;
7687 struct otri neighbor;
7690 struct osub neighborsubseg;
7700 while (virusloop != (
triangle **)NULL)
7702 testtri.
tri = *virusloop;
7715 sym(testtri, neighbor);
7717 tspivot(testtri, neighborsubseg);
7728 *regiontri = neighbor.
tri;
7741 while (virusloop != (
triangle **)NULL)
7743 testtri.
tri = *virusloop;
7770 struct otri searchtri;
7771 struct otri *regiontris;
7774 vertex searchorg, searchdest;
7787 regiontris = (
struct otri *)NULL;
7802 for (i = 0; i < 2 * holes; i += 2)
7805 if ((holelist[i] >= m->
xmin) && (holelist[i] <= m->
xmax) &&
7806 (holelist[i + 1] >= m->
ymin) && (holelist[i + 1] <= m->
ymax))
7815 org(searchtri, searchorg);
7816 dest(searchtri, searchdest);
7818 m, b, searchorg, searchdest, &holelist[i]) > 0.0)
7821 intersect =
locate(m, b, &holelist[i], &searchtri);
7830 *holetri = searchtri.
tri;
7847 for (i = 0; i < regions; i++)
7852 if ((regionlist[4 * i] >= m->
xmin) &&
7853 (regionlist[4 * i] <= m->
xmax) &&
7854 (regionlist[4 * i + 1] >= m->
ymin) &&
7855 (regionlist[4 * i + 1] <= m->
ymax))
7865 org(searchtri, searchorg);
7866 dest(searchtri, searchdest);
7868 m, b, searchorg, searchdest, ®ionlist[4 * i]) > 0.0)
7871 intersect =
locate(m, b, ®ionlist[4 * i], &searchtri);
7876 otricopy(searchtri, regiontris[i]);
7892 for (i = 0; i < regions; i++)
7894 if (regiontris[i].tri != m->
dummytri)
7898 if (!
deadtri(regiontris[i].tri))
7903 *regiontri = regiontris[i].
tri;
7906 m, b, regionlist[4 * i + 2], regionlist[4 * i + 3]);
7939 struct osub subsegloop;
7944 while (subsegloop.
ss != (
subseg *)NULL)
7960 printf(
"Try increasing the area criterion and/or reducing the minimum\n");
7961 printf(
" allowable angle so that tiny triangles are not created.\n");
7963 printf(
"Alternatively, try recompiling me with double precision\n");
7964 printf(
" arithmetic (by removing \"#define SINGLE\" from the\n");
7965 printf(
" source file or \"-DSINGLE\" from the makefile).\n");
7986 struct otri testtri;
7988 struct osub currentenc;
7990 vertex eorg, edest, eapex;
7993 double segmentlength, nearestpoweroftwo;
7995 double multiplier, divisor;
7996 int acuteorg, acuteorg2, acutedest, acutedest2;
8010 sorg(currentenc, eorg);
8011 sdest(currentenc, edest);
8038 lnext(enctri, testtri);
8049 if (!acuteorg && !acutedest)
8051 apex(enctri, eapex);
8053 ((eorg[0] - eapex[0]) * (edest[0] - eapex[0]) +
8054 (eorg[1] - eapex[1]) * (edest[1] - eapex[1]) <
8059 apex(enctri, eapex);
8060 lprev(enctri, testtri);
8067 sym(enctri, testtri);
8074 acutedest = acutedest || acutedest2;
8079 acuteorg = acuteorg || acuteorg2;
8083 if (!acuteorg2 && !acutedest2)
8085 org(testtri, eapex);
8088 ((eorg[0] - eapex[0]) * (edest[0] - eapex[0]) +
8089 (eorg[1] - eapex[1]) * (edest[1] - eapex[1]) <
8093 sym(enctri, testtri);
8094 apex(testtri, eapex);
8103 if (acuteorg || acutedest)
8106 sqrt((edest[0] - eorg[0]) * (edest[0] - eorg[0]) +
8107 (edest[1] - eorg[1]) * (edest[1] - eorg[1]));
8112 nearestpoweroftwo = 1.0;
8113 while (segmentlength > 3.0 * nearestpoweroftwo)
8115 nearestpoweroftwo *= 2.0;
8117 while (segmentlength < 1.5 * nearestpoweroftwo)
8119 nearestpoweroftwo *= 0.5;
8122 split = nearestpoweroftwo / segmentlength;
8125 split = 1.0 - split;
8138 for (i = 0; i < 2 + m->
nextras; i++)
8140 newvertex[i] = eorg[i] + split * (edest[i] - eorg[i]);
8150 divisor = ((eorg[0] - edest[0]) * (eorg[0] - edest[0]) +
8151 (eorg[1] - edest[1]) * (eorg[1] - edest[1]));
8152 if ((multiplier != 0.0) && (divisor != 0.0))
8154 multiplier = multiplier / divisor;
8156 if (multiplier == multiplier)
8158 newvertex[0] += multiplier * (edest[1] - eorg[1]);
8159 newvertex[1] += multiplier * (eorg[0] - edest[0]);
8167 if (((newvertex[0] == eorg[0]) && (newvertex[1] == eorg[1])) ||
8168 ((newvertex[0] == edest[0]) && (newvertex[1] == edest[1])))
8170 printf(
"Error: Ran out of precision at (%.12g, %.12g).\n",
8173 printf(
"I attempted to split a segment to a smaller size "
8176 " can be accommodated by the finite precision of\n");
8177 printf(
" floating point arithmetic.\n");
8183 m, b, newvertex, &enctri, ¤tenc, 1, triflaws);
8187 printf(
"Internal error in splitencsegs():\n");
8188 printf(
" Failure to split a segment.\n");
8216 struct otri triangleloop;
8240 struct otri badotri;
8241 vertex borg, bdest, bapex;
8250 dest(badotri, bdest);
8251 apex(badotri, bapex);
8266 if (((newvertex[0] == borg[0]) && (newvertex[1] == borg[1])) ||
8267 ((newvertex[0] == bdest[0]) && (newvertex[1] == bdest[1])) ||
8268 ((newvertex[0] == bapex[0]) && (newvertex[1] == bapex[1])))
8274 for (i = 2; i < 2 + m->
nextras; i++)
8277 newvertex[i] = borg[i] + xi * (bdest[i] - borg[i]) +
8278 eta * (bapex[i] - borg[i]);
8303 m, b, newvertex, &badotri, (
struct osub *)NULL, 1, 1);
8337 "This probably means that I am trying to refine triangles\n");
8339 " to a smaller size than can be accommodated by the finite\n");
8340 printf(
" precision of floating point arithmetic. (You can be\n");
8341 printf(
" sure of this if I fail to terminate.)\n");
8384 for (i = 0; i < 4096; i++)
8444 double *pointattriblist,
8445 int *pointmarkerlist,
8447 int numberofpointattribs)
8457 m->
nextras = numberofpointattribs;
8461 printf(
"Error: Input must have at least three input vertices.\n");
8478 x = vertexloop[0] = pointlist[coordindex++];
8479 y = vertexloop[1] = pointlist[coordindex++];
8481 for (j = 0; j < numberofpointattribs; j++)
8483 vertexloop[2 + j] = pointattriblist[attribindex++];
8485 if (pointmarkerlist != (
int *)NULL)
8528 double **pointattriblist,
8529 int **pointmarkerlist)
8551 if (*pointlist == (
double *)NULL)
8554 (
double *)
trimalloc((
int)(outvertices * 2 *
sizeof(double)));
8557 if ((m->
nextras > 0) && (*pointattriblist == (
double *)NULL))
8560 (
int)(outvertices * m->
nextras *
sizeof(double)));
8563 if (*pointmarkerlist == (
int *)NULL)
8565 *pointmarkerlist = (
int *)
trimalloc((
int)(outvertices *
sizeof(int)));
8568 palist = *pointattriblist;
8569 pmlist = *pointmarkerlist;
8576 while (vertexloop != (
vertex)NULL)
8581 plist[coordindex++] = vertexloop[0];
8582 plist[coordindex++] = vertexloop[1];
8584 for (i = 0; i < m->
nextras; i++)
8586 palist[attribindex++] = vertexloop[2 + i];
8590 pmlist[vertexnumber] =
vertexmark(vertexloop);
8617 while (vertexloop != (
vertex)NULL)
8637 double **triangleattriblist)
8643 struct otri triangleloop;
8649 if (*trianglelist == (
int *)NULL)
8656 if ((m->
eextras > 0) && (*triangleattriblist == (
double *)NULL))
8658 *triangleattriblist = (
double *)
trimalloc(
8661 tlist = *trianglelist;
8662 talist = *triangleattriblist;
8672 org(triangleloop, p1);
8673 dest(triangleloop, p2);
8674 apex(triangleloop, p3);
8682 for (i = 0; i < m->
eextras; i++)
8701 int **segmentmarkerlist)
8706 struct osub subsegloop;
8707 vertex endpoint1, endpoint2;
8711 if (*segmentlist == (
int *)NULL)
8717 if (*segmentmarkerlist == (
int *)NULL)
8719 *segmentmarkerlist =
8722 slist = *segmentlist;
8723 smlist = *segmentmarkerlist;
8730 while (subsegloop.
ss != (
subseg *)NULL)
8732 sorg(subsegloop, endpoint1);
8733 sdest(subsegloop, endpoint2);
8739 smlist[subsegnumber] =
mark(subsegloop);
8801 in.pointattributelist,
8804 in.numberofpointattributes);
8822 in.segmentmarkerlist,
8823 in.numberofsegments);
8828 holearray = in.holelist;
8829 m.
holes = in.numberofholes;
8830 regionarray = in.regionlist;
8831 m.
regions = in.numberofregions;
#define setdest(otri, vertexptr)
void triangledeinit(struct mesh *m, struct behavior *b)
void badsubsegdealloc(struct mesh *m, struct badsubseg *dyingseg)
struct badtriang * dequeuebadtriang(struct mesh *m)
struct flipstacker * lastflip
void alternateaxes(vertex *sortarray, int arraysize, int axis)
void vertexmedian(vertex *sortarray, int arraysize, int median, int axis)
#define segdest(osub, vertexptr)
void writepoly(struct mesh *m, struct behavior *b, int **segmentlist, int **segmentmarkerlist)
void writenodes(struct mesh *m, struct behavior *b, double **pointlist, double **pointattriblist, int **pointmarkerlist)
void tallyencs(struct mesh *m, struct behavior *b)
void formskeleton(struct mesh *m, struct behavior *b, int *segmentlist, int *segmentmarkerlist, int numberofsegments)
#define setmark(osub, value)
#define bond(otri1, otri2)
void vertexsort(vertex *sortarray, unsigned int arraysize)
int numberoftriangleattributes
void writeelements(struct mesh *m, struct behavior *b, int **trianglelist, double **triangleattriblist)
#define setsegdest(osub, vertexptr)
void subsegdealloc(struct mesh *m, subseg *dyingsubseg)
double incircle(struct mesh *m, struct behavior *b, vertex pa, vertex pb, vertex pc, vertex pd)
#define setapex(otri, vertexptr)
void pooldeinit(struct memorypool *pool)
unsigned long randomnation(unsigned int choices)
void * traverse(struct memorypool *pool)
#define sbond(osub1, osub2)
int triunsuitable(vertex triorg, vertex tridest, vertex triapex, double area)
#define stpivot(osub, otri)
struct badsubseg * badsubsegtraverse(struct mesh *m)
void * poolalloc(struct memorypool *pool)
#define Two_One_Product(a1, a0, b, x3, x2, x1, x0)
#define Two_Product_Presplit(a, b, bhi, blo, x, y)
struct memorypool subsegs
struct memorypool triangles
void * trimalloc(int size)
#define Fast_Two_Sum(a, b, x, y)
struct badtriang * queuetail[4096]
void maketriangle(struct mesh *m, struct behavior *b, struct otri *newotri)
struct triangulateio in out
void makevertexmap(struct mesh *m, struct behavior *b)
#define setsegorg(osub, vertexptr)
void poolinit(struct memorypool *pool, int bytecount, int itemcount, int firstitemcount, int alignment)
void poolzero(struct memorypool *pool)
#define decode(ptr, otri)
void carveholes(struct mesh *m, struct behavior *b, double *holelist, int holes, double *regionlist, int regions)
void triangulatepolygon(struct mesh *m, struct behavior *b, struct otri *firstedge, struct otri *lastedge, int edgecount, int doflip, int triflaws)
enum insertvertexresult insertvertex(struct mesh *m, struct behavior *b, vertex newvertex, struct otri *searchtri, struct osub *splitseg, int segmentflaws, int triflaws)
#define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0)
#define segorg(osub, vertexptr)
#define setvertextype(vx, value)
void enqueuebadtriang(struct mesh *m, struct behavior *b, struct badtriang *badtri)
#define Two_Diff_Tail(a, b, x, y)
double orient3dadapt(vertex pa, vertex pb, vertex pc, vertex pd, double aheight, double bheight, double cheight, double dheight, double permanent)
void traversalinit(struct memorypool *pool)
#define otriequal(otri1, otri2)
void triangledealloc(struct mesh *m, triangle *dyingtriangle)
void initializevertexpool(struct mesh *m, struct behavior *b)
void insertsegment(struct mesh *m, struct behavior *b, vertex endpoint1, vertex endpoint2, int newmark)
long delaunay(struct mesh *m, struct behavior *b)
#define setorg(otri, vertexptr)
struct memorypool badsubsegs
#define sym(otri1, otri2)
void testtriangle(struct mesh *m, struct behavior *b, struct otri *testtri)
struct badtriang * queuefront[4096]
vertex getvertex(struct mesh *m, struct behavior *b, int number)
void vertexdealloc(struct mesh *m, vertex dyingvertex)
#define tsbond(otri, osub)
void enqueuebadtri(struct mesh *m, struct behavior *b, struct otri *enqtri, double minedge, vertex enqapex, vertex enqorg, vertex enqdest)
#define lprev(otri1, otri2)
void unflip(struct mesh *m, struct behavior *b, struct otri *flipedge)
void segmentintersection(struct mesh *m, struct behavior *b, struct otri *splittri, struct osub *splitsubseg, vertex endpoint2)
#define setsdest(osub, vertexptr)
void mergehulls(struct mesh *m, struct behavior *b, struct otri *farleft, struct otri *innerleft, struct otri *innerright, struct otri *farright, int axis)
void trifree(void *memptr)
double estimateTRI(int elen, double *e)
void numbernodes(struct mesh *m, struct behavior *b)
#define setvertex2tri(vx, value)
void flip(struct mesh *m, struct behavior *b, struct otri *flipedge)
#define dnext(otri1, otri2)
long removeghosts(struct mesh *m, struct behavior *b, struct otri *startghost)
void markhull(struct mesh *m, struct behavior *b)
void insertsubseg(struct mesh *m, struct behavior *b, struct otri *tri, int subsegmark)
void divconqrecurse(struct mesh *m, struct behavior *b, vertex *sortarray, int vertices, int axis, struct otri *farleft, struct otri *farright)
double counterclockwiseadapt(vertex pa, vertex pb, vertex pc, double detsum)
struct memorypool splaynodes
int numberofpointattributes
struct memorypool vertices
void delaunayfixup(struct mesh *m, struct behavior *b, struct otri *fixuptri, int leftside)
void negate(NekPoint< DataType > &rhs)
enum locateresult preciselocate(struct mesh *m, struct behavior *b, vertex searchpoint, struct otri *searchtri, int stopatsubsegment)
#define sdest(osub, vertexptr)
long divconqdelaunay(struct mesh *m, struct behavior *b)
double counterclockwise(struct mesh *m, struct behavior *b, vertex pa, vertex pb, vertex pc)
#define setelemattribute(otri, attnum, value)
void triangleinit(struct mesh *m)
int fast_expansion_sum_zeroelimTRI(int elen, double *e, int flen, double *f, double *h)
#define FLIPSTACKERPERBLOCK
#define onext(otri1, otri2)
double orient3d(struct mesh *m, struct behavior *b, vertex pa, vertex pb, vertex pc, vertex pd, double aheight, double bheight, double cheight, double dheight)
double * pointattributelist
#define apex(otri, vertexptr)
void plague(struct mesh *m, struct behavior *b)
double * triangleattributelist
void poolrestart(struct memorypool *pool)
#define sorg(osub, vertexptr)
int scale_expansion_zeroelimTRI(int elen, double *e, double b, double *h)
void transfernodes(struct mesh *m, struct behavior *b, double *pointlist, double *pointattriblist, int *pointmarkerlist, int numberofpoints, int numberofpointattribs)
double nonregular(struct mesh *m, struct behavior *b, vertex pa, vertex pb, vertex pc, vertex pd)
void splittriangle(struct mesh *m, struct behavior *b, struct badtriang *badtri)
void undovertex(struct mesh *m, struct behavior *b)
#define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0)
triangle * triangletraverse(struct mesh *m)
void deletevertex(struct mesh *m, struct behavior *b, struct otri *deltri)
enum finddirectionresult finddirection(struct mesh *m, struct behavior *b, struct otri *searchtri, vertex searchpoint)
#define setvertexmark(vx, value)
#define elemattribute(otri, attnum)
#define oprev(otri1, otri2)
void splitencsegs(struct mesh *m, struct behavior *b, int triflaws)
void constrainededge(struct mesh *m, struct behavior *b, struct otri *starttri, vertex endpoint2, int newmark)
void infecthull(struct mesh *m, struct behavior *b)
void enforcequality(struct mesh *m, struct behavior *b)
void dummyinit(struct mesh *m, struct behavior *b, int trianglebytes, int subsegbytes)
#define otricopy(otri1, otri2)
#define tspivot(otri, osub)
#define lnext(otri1, otri2)
#define dest(otri, vertexptr)
#define Split(a, ahi, alo)
struct memorypool flipstackers
#define Two_Product(a, b, x, y)
void makesubseg(struct mesh *m, struct osub *newsubseg)
struct badtriang * nexttriang
int checkseg4encroach(struct mesh *m, struct behavior *b, struct osub *testsubseg)
vertex vertextraverse(struct mesh *m)
void initializetrisubpools(struct mesh *m, struct behavior *b)
struct memorypool badtriangles
void tallyfaces(struct mesh *m, struct behavior *b)
struct flipstacker * prevflip
#define BADSUBSEGPERBLOCK
subseg * subsegtraverse(struct mesh *m)
#define spivot(osub1, osub2)
#define setsorg(osub, vertexptr)
double incircleadaptTRI(vertex pa, vertex pb, vertex pc, vertex pd, double permanent)
void findcircumcenter(struct mesh *m, struct behavior *b, vertex torg, vertex tdest, vertex tapex, vertex circumcenter, double *xi, double *eta, int offcenter)
enum locateresult locate(struct mesh *m, struct behavior *b, vertex searchpoint, struct otri *searchtri)
int scoutsegment(struct mesh *m, struct behavior *b, struct otri *searchtri, vertex endpoint2, int newmark)
void pooldealloc(struct memorypool *pool, void *dyingitem)
#define org(otri, vertexptr)
void conformingedge(struct mesh *m, struct behavior *b, vertex endpoint1, vertex endpoint2, int newmark)
#define Two_Sum(a, b, x, y)
#define ssym(osub1, osub2)
#define dprev(otri1, otri2)
#define sdecode(sptr, osub)
void regionplague(struct mesh *m, struct behavior *b, double attribute, double area)
void parsecommandline(int argc, char **argv, struct behavior *b)