diff --git a/modules/core/src/main/java/org/locationtech/jts/triangulate/quadedge/TrianglePredicate.java b/modules/core/src/main/java/org/locationtech/jts/triangulate/quadedge/TrianglePredicate.java index 9397408103..5b45a02cd7 100644 --- a/modules/core/src/main/java/org/locationtech/jts/triangulate/quadedge/TrianglePredicate.java +++ b/modules/core/src/main/java/org/locationtech/jts/triangulate/quadedge/TrianglePredicate.java @@ -111,7 +111,8 @@ private static double triArea(Coordinate a, Coordinate b, Coordinate c) { /** * Tests if a point is inside the circle defined by * the triangle with vertices a, b, c (oriented counter-clockwise). - * This method uses more robust computation. + * This method uses a fast errorbound filter and falls back to a + * precise computation if it cannot guarantee correctness. * * @param a a vertex of the triangle * @param b a vertex of the triangle @@ -123,9 +124,47 @@ public static boolean isInCircleRobust( Coordinate a, Coordinate b, Coordinate c, Coordinate p) { - //checkRobustInCircle(a, b, c, p); -// return isInCircleNonRobust(a, b, c, p); - return isInCircleNormalized(a, b, c, p); + double adx, bdx, cdx, ady, bdy, cdy; + double bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady; + double alift, blift, clift; + double det; + double errbound, permanent; + + adx = a.x - p.x; + bdx = b.x - p.x; + cdx = c.x - p.x; + ady = a.y - p.y; + bdy = b.y - p.y; + cdy = c.y - p.y; + + bdxcdy = bdx * cdy; + cdxbdy = cdx * bdy; + alift = adx * adx + ady * ady; + + cdxady = cdx * ady; + adxcdy = adx * cdy; + blift = bdx * bdx + bdy * bdy; + + adxbdy = adx * bdy; + bdxady = bdx * ady; + clift = cdx * cdx + cdy * cdy; + + det = alift * (bdxcdy - cdxbdy) + + blift * (cdxady - adxcdy) + + clift * (adxbdy - bdxady); + + permanent = (Math.abs(bdxcdy) + Math.abs(cdxbdy)) * alift + + (Math.abs(cdxady) + Math.abs(adxcdy)) * blift + + (Math.abs(adxbdy) + Math.abs(bdxady)) * clift; + errbound = iccerrboundA * permanent; + if (Math.abs(det) >= errbound) { + return det > 0; + } + double[] pa = new double[] {a.x, a.y}; + double[] pb = new double[] {b.x, b.y}; + double[] pc = new double[] {c.x, c.y}; + double[] pd = new double[] {p.x, p.y}; + return incircleadapt(pa, pb, pc, pd, permanent) > 0; } /** @@ -312,5 +351,874 @@ private static void checkRobustInCircle(Coordinate a, Coordinate b, Coordinate c } */ + /** + * Expansion arithmetic methods from Shewchuk for exact IsInCircle Predicate. + * + * This is ported directly from the public domain implementation at + * https://www.cs.cmu.edu/afs/cs/project/quake/public/code/predicates.c . + * + * Code is intentionally kept as close to the original as possible to allow + * reasonable debugging by comparison to the original C code. No documentation + * of these private methods, expect what is copied from the original. + */ + private final static double epsilon = Math.ulp(0.5); + private final static double resulterrbound = (3.0 + 8.0 * epsilon) * epsilon; + private final static double iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon; + private final static double iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon; + private final static double iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon; + + private static double[] Fast_Two_Sum(double a, double b) + { + double x = a + b; + double bvirt = x - a; + double y = b - bvirt; + return new double[] {x, y}; + } + + private static double[] Two_Sum(double a, double b) + { + double x = a + b; + double bvirt = x - a; + double avirt = x - bvirt; + double bround = b - bvirt; + double around = a - avirt; + double y = around + bround; + return new double[] {x, y}; + } + + private static double Two_Diff_Tail(double a, double b, double x) + { + double bvirt = a - x; + double avirt = x + bvirt; + double bround = bvirt - b; + double around = a - avirt; + return around + bround; + } + + private static double[] Two_Diff(double a, double b) + { + double x = a - b; + double y = Two_Diff_Tail(a, b, x); + return new double[] {x, y}; + } + + private static double[] Two_Product(double a, double b) { + double x = a * b; + double y = Math.fma(a, b, -x); + return new double[] {x, y}; + } + + private static double[] Two_One_Sum(double a1, double a0, double b) + { + double[] _1x0 = Two_Sum(a0, b); + double[] x2x1 = Two_Sum(a1, _1x0[0]); + return new double[] {x2x1[0], x2x1[1], _1x0[1]}; + } + + private static double[] Two_One_Diff(double a1, double a0, double b) + { + double[] _1x0 = Two_Diff(a0, b); + double[] x2x1 = Two_Sum( a1, _1x0[0]); + return new double[] {x2x1[0], x2x1[1], _1x0[1]}; + } + + private static double[] Two_Two_Sum(double a1, double a0, double b1, double b0) + { + double[] _j_0x0 = Two_One_Sum(a1, a0, b0); + double[] x3x2x1 = Two_One_Sum(_j_0x0[0], _j_0x0[1], b1); + return new double[] {x3x2x1[0], x3x2x1[1], x3x2x1[2], _j_0x0[2]}; + } + + private static double[] Two_Two_Diff(double a1, double a0, double b1, double b0) + { + double[] _j_0x0 = Two_One_Diff(a1, a0, b0); + double[] x3x2x1 = Two_One_Diff(_j_0x0[0], _j_0x0[1], b1); + return new double[] {x3x2x1[0], x3x2x1[1], x3x2x1[2], _j_0x0[2]}; + } + + /*****************************************************************************/ + /* */ + /* fast_expansion_sum_zeroelim() Sum two expansions, eliminating zero */ + /* components from the output expansion. */ + /* */ + /* Sets h = e + f. See the long version of my paper for details. */ + /* */ + /* If round-to-even is used (as with IEEE 754), maintains the strongly */ + /* nonoverlapping property. (That is, if e is strongly nonoverlapping, h */ + /* will be also.) Does NOT maintain the nonoverlapping or nonadjacent */ + /* properties. */ + /* */ + /*****************************************************************************/ + + private static double[] fast_expansion_sum_zeroelim(double[] e, double[] f) /* h cannot be e or f. */ + { + double Q; + double Qnew; + double hh; + int eindex, findex, hindex; + double enow, fnow; + + double[] h = new double[e.length + f.length]; + + enow = e[0]; + fnow = f[0]; + eindex = findex = 0; + if ((fnow > enow) == (fnow > -enow)) { + Q = enow; + if(++eindex < e.length) + { + enow = e[eindex]; + } + } else { + Q = fnow; + if(++findex < f.length) + { + fnow = f[findex]; + } + } + hindex = 0; + if ((eindex < e.length) && (findex < f.length)) { + if ((fnow > enow) == (fnow > -enow)) { + double[] Qnhh = Fast_Two_Sum(enow, Q); + Qnew = Qnhh[0]; + hh = Qnhh[1]; + if(++eindex < e.length) + { + enow = e[eindex]; + } + } else { + double[] Qnhh = Fast_Two_Sum(fnow, Q); + Qnew = Qnhh[0]; + hh = Qnhh[1]; + if(++findex < f.length) + { + fnow = f[findex]; + } + } + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + while ((eindex < e.length) && (findex < f.length)) { + if ((fnow > enow) == (fnow > -enow)) { + double[] Qnhh = Two_Sum(Q, enow); + Qnew = Qnhh[0]; + hh = Qnhh[1]; + if(++eindex < e.length) + { + enow = e[eindex]; + } + } else { + double[] Qnhh = Two_Sum(Q, fnow); + Qnew = Qnhh[0]; + hh = Qnhh[1]; + if(++findex < f.length) + { + fnow = f[findex]; + } + } + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + } + } + while (eindex < e.length) { + double[] Qnhh = Two_Sum(Q, enow); + Qnew = Qnhh[0]; + hh = Qnhh[1]; + if(++eindex < e.length) + { + enow = e[eindex]; + } + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + } + while (findex < f.length) { + double[] Qnhh = Two_Sum(Q, fnow); + Qnew = Qnhh[0]; + hh = Qnhh[1]; + if(++findex < f.length) + { + fnow = f[findex]; + } + Q = Qnew; + if (hh != 0.0) { + h[hindex++] = hh; + } + } + if ((Q != 0.0) || (hindex == 0)) { + h[hindex++] = Q; + } + + double[] h_out = new double[hindex]; + System.arraycopy(h, 0, h_out, 0, hindex); + return h_out; + } + + /*****************************************************************************/ + /* */ + /* scale_expansion_zeroelim() Multiply an expansion by a scalar, */ + /* eliminating zero components from the */ + /* output expansion. */ + /* */ + /* Sets h = be. See either version of my paper for details. */ + /* */ + /* Maintains the nonoverlapping property. If round-to-even is used (as */ + /* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */ + /* properties as well. (That is, if e has one of these properties, so */ + /* will h.) */ + /* */ + /*****************************************************************************/ + + private static double[] scale_expansion_zeroelim(double[] e, double b) + { + double Q, sum; + double hh; + double product1; + double product0; + int eindex, hindex; + double enow; + + double[] h = new double[e.length * 2]; + + double[] Qhh = Two_Product(e[0], b); + Q = Qhh[0]; + hh = Qhh[1]; + hindex = 0; + if (hh != 0) { + h[hindex++] = hh; + } + for (eindex = 1; eindex < e.length; eindex++) { + enow = e[eindex]; + double[] p1p0 = Two_Product(enow, b); + product1 = p1p0[0]; + product0 = p1p0[1]; + double[] sumhh = Two_Sum(Q, product0); + sum = sumhh[0]; + hh = sumhh[1]; + if (hh != 0) { + h[hindex++] = hh; + } + Qhh = Fast_Two_Sum(product1, sum); + Q = Qhh[0]; + hh = Qhh[1]; + if (hh != 0) { + h[hindex++] = hh; + } + } + if ((Q != 0.0) || (hindex == 0)) { + h[hindex++] = Q; + } + double[] h_out = new double[hindex]; + System.arraycopy(h, 0, h_out, 0, hindex); + return h_out; + } + + /*****************************************************************************/ + /* */ + /* estimate() Produce a one-word estimate of an expansion's value. */ + /* */ + /* See either version of my paper for details. */ + /* */ + /*****************************************************************************/ + + private static double estimate(double[] e) + { + double Q; + int eindex; + + Q = e[0]; + for (eindex = 1; eindex < e.length; eindex++) { + Q += e[eindex]; + } + return Q; + } + + /*****************************************************************************/ + /* */ + /* incircleadapt() Adaptive exact 2D incircle test. Robust. */ + /* */ + /* Return a positive value if the point pd lies inside the */ + /* circle passing through pa, pb, and pc; a negative value if */ + /* it lies outside; and zero if the four points are cocircular.*/ + /* The points pa, pb, and pc must be in counterclockwise */ + /* order, or the sign of the result will be reversed. */ + /* */ + /* Only the first and last routine should be used; the middle two are for */ + /* timings. */ + /*****************************************************************************/ + + private static double incircleadapt(double[] pa, double[] pb, double[] pc, double[] pd, double permanent) + { + double adx, bdx, cdx, ady, bdy, cdy; + double det; + double errbound; + + adx = pa[0] - pd[0]; + bdx = pb[0] - pd[0]; + cdx = pc[0] - pd[0]; + ady = pa[1] - pd[1]; + bdy = pb[1] - pd[1]; + cdy = pc[1] - pd[1]; + + double bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1; + double bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0; + double[] bc = new double[4]; + double[] ca = new double[4]; + double[] ab = new double[4]; + double bc3, ca3, ab3; + double[] axbc, axxbc, aybc, ayybc, adet; + double[] bxca, bxxca, byca, byyca, bdet; + double[] cxab, cxxab, cyab, cyyab, cdet; + double[] abdet; + double[] fin1, finnow, finother, finswap; + + double adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail; + double adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1; + double adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0; + double[] aa = new double[4]; + double[] bb = new double[4]; + double[] cc = new double[4]; + double aa3, bb3, cc3; + double ti1, tj1; + double ti0, tj0; + double[] u = new double[4]; + double[] v = new double[4]; + double u3, v3; + double[] temp8, temp16a, temp16b, temp16c, temp32a, temp32b, temp48, temp64; + double[] axtbb, axtcc, aytbb, aytcc; + double[] bxtaa, bxtcc, bytaa, bytcc; + double[] cxtaa, cxtbb, cytaa, cytbb; + double[] axtbc = new double[8]; + double[] aytbc = new double[8]; + double[] bxtca = new double[8]; + double[] bytca = new double[8]; + double[] cxtab = new double[8]; + double[] cytab = new double[8]; + double[] axtbct, aytbct, bxtcat, bytcat, cxtabt, cytabt; + double[] axtbctt, aytbctt, bxtcatt, bytcatt, cxtabtt, cytabtt; + double[] abt, bct, cat; + double[] abtt = new double[4]; + double[] bctt = new double[4]; + double[] catt = new double[4]; + double bctt3, catt3; + double negate; + + double tmp[] = Two_Product(bdx, cdy); + bdxcdy1 = tmp[0]; + bdxcdy0 = tmp[1]; + tmp = Two_Product(cdx, bdy); + cdxbdy1 = tmp[0]; + cdxbdy0 = tmp[1]; + tmp = Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0); + bc3 = tmp[0]; + bc[2] = tmp[1]; + bc[1] = tmp[2]; + bc[0] = tmp[3]; + bc[3] = bc3; + axbc = scale_expansion_zeroelim(bc, adx); + axxbc = scale_expansion_zeroelim(axbc, adx); + aybc = scale_expansion_zeroelim(bc, ady); + ayybc = scale_expansion_zeroelim(aybc, ady); + adet = fast_expansion_sum_zeroelim(axxbc, ayybc); + + tmp = Two_Product(cdx, ady); + cdxady1 = tmp[0]; + cdxady0 = tmp[1]; + tmp = Two_Product(adx, cdy); + adxcdy1 = tmp[0]; + adxcdy0 = tmp[1]; + tmp = Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0); + ca3 = tmp[0]; + ca[2] = tmp[1]; + ca[1] = tmp[2]; + ca[0] = tmp[3]; + ca[3] = ca3; + bxca = scale_expansion_zeroelim(ca, bdx); + bxxca = scale_expansion_zeroelim(bxca, bdx); + byca = scale_expansion_zeroelim(ca, bdy); + byyca = scale_expansion_zeroelim(byca, bdy); + bdet = fast_expansion_sum_zeroelim(bxxca, byyca); + + tmp = Two_Product(adx, bdy); + adxbdy1 = tmp[0]; + adxbdy0 = tmp[1]; + tmp = Two_Product(bdx, ady); + bdxady1 = tmp[0]; + bdxady0 = tmp[1]; + tmp = Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0); + ab3 = tmp[0]; + ab[2] = tmp[1]; + ab[1] = tmp[2]; + ab[0] = tmp[3]; + ab[3] = ab3; + cxab = scale_expansion_zeroelim(ab, cdx); + cxxab = scale_expansion_zeroelim(cxab, cdx); + cyab = scale_expansion_zeroelim(ab, cdy); + cyyab = scale_expansion_zeroelim(cyab, cdy); + cdet = fast_expansion_sum_zeroelim(cxxab, cyyab); + + abdet = fast_expansion_sum_zeroelim(adet, bdet); + fin1 = fast_expansion_sum_zeroelim(abdet, cdet); + + det = estimate(fin1); + errbound = iccerrboundB * permanent; + if (Math.abs(det) >= errbound) { + return det; + } + + adxtail = Two_Diff_Tail(pa[0], pd[0], adx); + adytail = Two_Diff_Tail(pa[1], pd[1], ady); + bdxtail = Two_Diff_Tail(pb[0], pd[0], bdx); + bdytail = Two_Diff_Tail(pb[1], pd[1], bdy); + cdxtail = Two_Diff_Tail(pc[0], pd[0], cdx); + cdytail = Two_Diff_Tail(pc[1], pd[1], cdy); + if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) + && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)) { + return det; + } + + errbound = iccerrboundC * permanent + resulterrbound * Math.abs(det); + det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail) + - (bdy * cdxtail + cdx * bdytail)) + + 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx)) + + ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail) + - (cdy * adxtail + adx * cdytail)) + + 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx)) + + ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail) + - (ady * bdxtail + bdx * adytail)) + + 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx)); + if (Math.abs(det) >= errbound) { + return det; + } + + finnow = fin1; + + if ((bdxtail != 0.0) || (bdytail != 0.0) + || (cdxtail != 0.0) || (cdytail != 0.0)) { + tmp = Two_Product(adx, adx); + adxadx1 = tmp[0]; + adxadx0 = tmp[1]; + tmp = Two_Product(ady, ady); + adyady1 = tmp[0]; + adyady0 = tmp[1]; + tmp = Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0); + aa3 = tmp[0]; + aa[2] = tmp[1]; + aa[1] = tmp[2]; + aa[0] = tmp[3]; + aa[3] = aa3; + } + if ((cdxtail != 0.0) || (cdytail != 0.0) + || (adxtail != 0.0) || (adytail != 0.0)) { + tmp = Two_Product(bdx, bdx); + bdxbdx1 = tmp[0]; + bdxbdx0 = tmp[1]; + tmp = Two_Product(bdy, bdy); + bdybdy1 = tmp[0]; + bdybdy0 = tmp[1]; + tmp = Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0); + bb3 = tmp[0]; + bb[2] = tmp[1]; + bb[1] = tmp[2]; + bb[0] = tmp[3]; + bb[3] = bb3; + } + if ((adxtail != 0.0) || (adytail != 0.0) + || (bdxtail != 0.0) || (bdytail != 0.0)) { + tmp = Two_Product(cdx, cdx); + cdxcdx1 = tmp[0]; + cdxcdx0 = tmp[1]; + tmp = Two_Product(cdy, cdy); + cdycdy1 = tmp[0]; + cdycdy0 = tmp[1]; + tmp = Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0); + cc3 = tmp[0]; + cc[2] = tmp[1]; + cc[1] = tmp[2]; + cc[0] = tmp[3]; + cc[3] = cc3; + } + + if (adxtail != 0.0) { + axtbc = scale_expansion_zeroelim(bc, adxtail); + temp16a = scale_expansion_zeroelim(axtbc, 2.0 * adx); + + axtcc = scale_expansion_zeroelim(cc, adxtail); + temp16b = scale_expansion_zeroelim(axtcc, bdy); + + axtbb = scale_expansion_zeroelim(bb, adxtail); + temp16c = scale_expansion_zeroelim(axtbb, -cdy); + + temp32a = fast_expansion_sum_zeroelim(temp16a, temp16b); + temp48 = fast_expansion_sum_zeroelim(temp16c, temp32a); + finother = fast_expansion_sum_zeroelim(finnow, temp48); + finswap = finnow; finnow = finother; finother = finswap; + } + if (adytail != 0.0) { + aytbc = scale_expansion_zeroelim(bc, adytail); + temp16a = scale_expansion_zeroelim(aytbc, 2.0 * ady); + + aytbb = scale_expansion_zeroelim(bb, adytail); + temp16b = scale_expansion_zeroelim(aytbb, cdx); + + aytcc = scale_expansion_zeroelim(cc, adytail); + temp16c = scale_expansion_zeroelim(aytcc, -bdx); + + temp32a = fast_expansion_sum_zeroelim(temp16a, temp16b); + temp48 = fast_expansion_sum_zeroelim(temp16c, temp32a); + finother = fast_expansion_sum_zeroelim(finnow, temp48); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdxtail != 0.0) { + bxtca = scale_expansion_zeroelim(ca, bdxtail); + temp16a = scale_expansion_zeroelim(bxtca, 2.0 * bdx); + + bxtaa = scale_expansion_zeroelim(aa, bdxtail); + temp16b = scale_expansion_zeroelim(bxtaa, cdy); + + bxtcc = scale_expansion_zeroelim(cc, bdxtail); + temp16c = scale_expansion_zeroelim(bxtcc, -ady); + + temp32a = fast_expansion_sum_zeroelim(temp16a, temp16b); + temp48 = fast_expansion_sum_zeroelim(temp16c, temp32a); + finother = fast_expansion_sum_zeroelim(finnow, temp48); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdytail != 0.0) { + bytca = scale_expansion_zeroelim(ca, bdytail); + temp16a = scale_expansion_zeroelim(bytca, 2.0 * bdy); + bytcc = scale_expansion_zeroelim(cc, bdytail); + temp16b = scale_expansion_zeroelim(bytcc, adx); + + bytaa = scale_expansion_zeroelim(aa, bdytail); + temp16c = scale_expansion_zeroelim(bytaa, -cdx); + + temp32a = fast_expansion_sum_zeroelim(temp16a, temp16b); + temp48 = fast_expansion_sum_zeroelim(temp16c, temp32a); + finother = fast_expansion_sum_zeroelim(finnow, temp48); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdxtail != 0.0) { + cxtab = scale_expansion_zeroelim(ab, cdxtail); + temp16a = scale_expansion_zeroelim(cxtab, 2.0 * cdx); + + cxtbb = scale_expansion_zeroelim(bb, cdxtail); + temp16b = scale_expansion_zeroelim(cxtbb, ady); + + cxtaa = scale_expansion_zeroelim(aa, cdxtail); + temp16c = scale_expansion_zeroelim(cxtaa, -bdy); + + temp32a = fast_expansion_sum_zeroelim(temp16a, temp16b); + temp48 = fast_expansion_sum_zeroelim(temp16c, temp32a); + finother = fast_expansion_sum_zeroelim(finnow, temp48); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdytail != 0.0) { + cytab = scale_expansion_zeroelim(ab, cdytail); + temp16a = scale_expansion_zeroelim(cytab, 2.0 * cdy); + + cytaa = scale_expansion_zeroelim(aa, cdytail); + temp16b = scale_expansion_zeroelim(cytaa, bdx); + + cytbb = scale_expansion_zeroelim(bb, cdytail); + temp16c = scale_expansion_zeroelim(cytbb, -adx); + + temp32a = fast_expansion_sum_zeroelim(temp16a, temp16b); + temp48 = fast_expansion_sum_zeroelim(temp16c, temp32a); + finother = fast_expansion_sum_zeroelim(finnow, temp48); + finswap = finnow; finnow = finother; finother = finswap; + } + + if ((adxtail != 0.0) || (adytail != 0.0)) { + if ((bdxtail != 0.0) || (bdytail != 0.0) + || (cdxtail != 0.0) || (cdytail != 0.0)) { + tmp = Two_Product(bdxtail, cdy); + ti1 = tmp[0]; + ti0 = tmp[1]; + tmp = Two_Product(bdx, cdytail); + tj1 = tmp[0]; + tj0 = tmp[1]; + tmp = Two_Two_Sum(ti1, ti0, tj1, tj0); + u3 = tmp[0]; + u[2] = tmp[1]; + u[1] = tmp[2]; + u[0] = tmp[3]; + u[3] = u3; + negate = -bdy; + tmp = Two_Product(cdxtail, negate); + ti1 = tmp[0]; + ti0 = tmp[1]; + negate = -bdytail; + tmp = Two_Product(cdx, negate); + tj1 = tmp[0]; + tj0 = tmp[1]; + tmp = Two_Two_Sum(ti1, ti0, tj1, tj0); + v3 = tmp[0]; + v[2] = tmp[1]; + v[1] = tmp[2]; + v[0] = tmp[3]; + v[3] = v3; + bct = fast_expansion_sum_zeroelim(u, v); + + tmp = Two_Product(bdxtail, cdytail); + ti1 = tmp[0]; + ti0 = tmp[1]; + tmp = Two_Product(cdxtail, bdytail); + tj1 = tmp[0]; + tj0 = tmp[1]; + tmp = Two_Two_Diff(ti1, ti0, tj1, tj0); + bctt3 = tmp[0]; + bctt[2] = tmp[1]; + bctt[1] = tmp[2]; + bctt[0] = tmp[3]; + bctt[3] = bctt3; + } else { + bct = new double[] {0.0}; + bctt = new double[] {0.0}; + } + + if (adxtail != 0.0) { + temp16a = scale_expansion_zeroelim(axtbc, adxtail); + axtbct = scale_expansion_zeroelim(bct, adxtail); + temp32a = scale_expansion_zeroelim(axtbct, 2.0 * adx); + temp48 = fast_expansion_sum_zeroelim(temp16a, temp32a); + finother = fast_expansion_sum_zeroelim(finnow, temp48); + finswap = finnow; finnow = finother; finother = finswap; + if (bdytail != 0.0) { + temp8 = scale_expansion_zeroelim(cc, adxtail); + temp16a = scale_expansion_zeroelim(temp8, bdytail); + finother = fast_expansion_sum_zeroelim(finnow, temp16a); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdytail != 0.0) { + temp8 = scale_expansion_zeroelim(bb, -adxtail); + temp16a = scale_expansion_zeroelim(temp8, cdytail); + finother = fast_expansion_sum_zeroelim(finnow, temp16a); + finswap = finnow; finnow = finother; finother = finswap; + } + + temp32a = scale_expansion_zeroelim(axtbct, adxtail); + axtbctt = scale_expansion_zeroelim(bctt, adxtail); + temp16a = scale_expansion_zeroelim(axtbctt, 2.0 * adx); + temp16b = scale_expansion_zeroelim(axtbctt, adxtail); + temp32b = fast_expansion_sum_zeroelim(temp16a, temp16b); + temp64 = fast_expansion_sum_zeroelim(temp32a, temp32b); + finother = fast_expansion_sum_zeroelim(finnow, temp64); + finswap = finnow; finnow = finother; finother = finswap; + } + if (adytail != 0.0) { + temp16a = scale_expansion_zeroelim(aytbc, adytail); + aytbct = scale_expansion_zeroelim(bct, adytail); + temp32a = scale_expansion_zeroelim(aytbct, 2.0 * ady); + temp48 = fast_expansion_sum_zeroelim(temp16a, temp32a); + finother = fast_expansion_sum_zeroelim(finnow, temp48); + finswap = finnow; finnow = finother; finother = finswap; + + + temp32a = scale_expansion_zeroelim(aytbct, adytail); + aytbctt = scale_expansion_zeroelim(bctt, adytail); + temp16a = scale_expansion_zeroelim(aytbctt, 2.0 * ady); + temp16b = scale_expansion_zeroelim(aytbctt, adytail); + temp32b = fast_expansion_sum_zeroelim(temp16a, temp16b); + temp64 = fast_expansion_sum_zeroelim(temp32a, temp32b); + finother = fast_expansion_sum_zeroelim(finnow, temp64); + finswap = finnow; finnow = finother; finother = finswap; + } + } + if ((bdxtail != 0.0) || (bdytail != 0.0)) { + if ((cdxtail != 0.0) || (cdytail != 0.0) + || (adxtail != 0.0) || (adytail != 0.0)) { + tmp = Two_Product(cdxtail, ady); + ti1 = tmp[0]; + ti0 = tmp[1]; + tmp = Two_Product(cdx, adytail); + tj1 = tmp[0]; + tj0 = tmp[1]; + tmp = Two_Two_Sum(ti1, ti0, tj1, tj0); + u3 = tmp[0]; + u[2] = tmp[1]; + u[1] = tmp[2]; + u[0] = tmp[3]; + u[3] = u3; + negate = -cdy; + tmp = Two_Product(adxtail, negate); + ti1 = tmp[0]; + ti0 = tmp[1]; + negate = -cdytail; + tmp = Two_Product(adx, negate); + tj1 = tmp[0]; + tj0 = tmp[1]; + tmp = Two_Two_Sum(ti1, ti0, tj1, tj0); + v3 = tmp[0]; + v[2] = tmp[1]; + v[1] = tmp[2]; + v[0] = tmp[3]; + v[3] = v3; + cat = fast_expansion_sum_zeroelim(u, v); + + tmp = Two_Product(cdxtail, adytail); + ti1 = tmp[0]; + ti0 = tmp[1]; + tmp = Two_Product(adxtail, cdytail); + tj1 = tmp[0]; + tj0 = tmp[1]; + tmp = Two_Two_Diff(ti1, ti0, tj1, tj0); + catt3 = tmp[0]; + catt[2] = tmp[1]; + catt[1] = tmp[2]; + catt[0] = tmp[3]; + catt[3] = catt3; + } else { + cat = new double[] {0.}; + catt = new double[] {0.}; + } + + if (bdxtail != 0.0) { + temp16a = scale_expansion_zeroelim(bxtca, bdxtail); + bxtcat = scale_expansion_zeroelim(cat, bdxtail); + temp32a = scale_expansion_zeroelim(bxtcat, 2.0 * bdx); + temp48 = fast_expansion_sum_zeroelim(temp16a, temp32a); + finother = fast_expansion_sum_zeroelim(finnow, temp48); + finswap = finnow; finnow = finother; finother = finswap; + if (cdytail != 0.0) { + temp8 = scale_expansion_zeroelim(aa, bdxtail); + temp16a = scale_expansion_zeroelim(temp8, cdytail); + finother = fast_expansion_sum_zeroelim(finnow, temp16a); + finswap = finnow; finnow = finother; finother = finswap; + } + if (adytail != 0.0) { + temp8 = scale_expansion_zeroelim(cc, -bdxtail); + temp16a = scale_expansion_zeroelim(temp8, adytail); + finother = fast_expansion_sum_zeroelim(finnow, temp16a); + finswap = finnow; finnow = finother; finother = finswap; + } + + temp32a = scale_expansion_zeroelim(bxtcat, bdxtail); + bxtcatt = scale_expansion_zeroelim(catt, bdxtail); + temp16a = scale_expansion_zeroelim(bxtcatt, 2.0 * bdx); + temp16b = scale_expansion_zeroelim(bxtcatt, bdxtail); + temp32b = fast_expansion_sum_zeroelim(temp16a, temp16b); + temp64 = fast_expansion_sum_zeroelim(temp32a, temp32b); + finother = fast_expansion_sum_zeroelim(finnow, temp64); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdytail != 0.0) { + temp16a = scale_expansion_zeroelim(bytca, bdytail); + bytcat = scale_expansion_zeroelim(cat, bdytail); + temp32a = scale_expansion_zeroelim(bytcat, 2.0 * bdy); + temp48 = fast_expansion_sum_zeroelim(temp16a, temp32a); + finother = fast_expansion_sum_zeroelim(finnow, temp48); + finswap = finnow; finnow = finother; finother = finswap; + + + temp32a = scale_expansion_zeroelim(bytcat, bdytail); + bytcatt = scale_expansion_zeroelim(catt, bdytail); + temp16a = scale_expansion_zeroelim(bytcatt, 2.0 * bdy); + temp16b = scale_expansion_zeroelim(bytcatt, bdytail); + temp32b = fast_expansion_sum_zeroelim(temp16a, temp16b); + temp64 = fast_expansion_sum_zeroelim(temp32a, temp32b); + finother = fast_expansion_sum_zeroelim(finnow, temp64); + finswap = finnow; finnow = finother; finother = finswap; + } + } + if ((cdxtail != 0.0) || (cdytail != 0.0)) { + if ((adxtail != 0.0) || (adytail != 0.0) + || (bdxtail != 0.0) || (bdytail != 0.0)) { + tmp = Two_Product(adxtail, bdy); + ti1 = tmp[0]; + ti0 = tmp[1]; + tmp = Two_Product(adx, bdytail); + tj1 = tmp[0]; + tj0 = tmp[1]; + tmp = Two_Two_Sum(ti1, ti0, tj1, tj0); + u3 = tmp[0]; + u[2] = tmp[1]; + u[1] = tmp[2]; + u[0] = tmp[3]; + u[3] = u3; + negate = -ady; + tmp = Two_Product(bdxtail, negate); + ti1 = tmp[0]; + ti0 = tmp[1]; + negate = -adytail; + tmp = Two_Product(bdx, negate); + tj1 = tmp[0]; + tj0 = tmp[1]; + tmp = Two_Two_Sum(ti1, ti0, tj1, tj0); + v3 = tmp[0]; + v[2] = tmp[1]; + v[1] = tmp[2]; + v[0] = tmp[3]; + v[3] = v3; + abt = fast_expansion_sum_zeroelim(u, v); + + tmp = Two_Product(adxtail, bdytail); + ti1 = tmp[0]; + ti0 = tmp[1]; + tmp = Two_Product(bdxtail, adytail); + tj1 = tmp[0]; + tj0 = tmp[1]; + tmp = Two_Two_Diff(ti1, ti0, tj1, tj0); + abtt = new double[] {tmp[3], tmp[2], tmp[1], tmp[0]}; + } else { + abt = new double[] {0.0}; + abtt = new double[] {0.0}; + } + + if (cdxtail != 0.0) { + temp16a = scale_expansion_zeroelim(cxtab, cdxtail); + cxtabt = scale_expansion_zeroelim(abt, cdxtail); + temp32a = scale_expansion_zeroelim(cxtabt, 2.0 * cdx); + temp48 = fast_expansion_sum_zeroelim(temp16a, temp32a); + finother = fast_expansion_sum_zeroelim(finnow, temp48); + finswap = finnow; finnow = finother; finother = finswap; + if (adytail != 0.0) { + temp8 = scale_expansion_zeroelim(bb, cdxtail); + temp16a = scale_expansion_zeroelim(temp8, adytail); + finother = fast_expansion_sum_zeroelim(finnow, temp16a); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdytail != 0.0) { + temp8 = scale_expansion_zeroelim(aa, -cdxtail); + temp16a = scale_expansion_zeroelim(temp8, bdytail); + finother = fast_expansion_sum_zeroelim(finnow, temp16a); + finswap = finnow; finnow = finother; finother = finswap; + } + + temp32a = scale_expansion_zeroelim(cxtabt, cdxtail); + cxtabtt = scale_expansion_zeroelim(abtt, cdxtail); + temp16a = scale_expansion_zeroelim(cxtabtt, 2.0 * cdx); + temp16b = scale_expansion_zeroelim(cxtabtt, cdxtail); + temp32b = fast_expansion_sum_zeroelim(temp16a, temp16b); + temp64 = fast_expansion_sum_zeroelim(temp32a, temp32b); + finother = fast_expansion_sum_zeroelim(finnow, temp64); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdytail != 0.0) { + temp16a = scale_expansion_zeroelim(cytab, cdytail); + cytabt = scale_expansion_zeroelim(abt, cdytail); + temp32a = scale_expansion_zeroelim(cytabt, 2.0 * cdy); + temp48 = fast_expansion_sum_zeroelim(temp16a, temp32a); + finother = fast_expansion_sum_zeroelim(finnow, temp48); + finswap = finnow; finnow = finother; finother = finswap; + + + temp32a = scale_expansion_zeroelim(cytabt, cdytail); + cytabtt = scale_expansion_zeroelim(abtt, cdytail); + temp16a = scale_expansion_zeroelim(cytabtt, 2.0 * cdy); + temp16b = scale_expansion_zeroelim(cytabtt, cdytail); + temp32b = fast_expansion_sum_zeroelim(temp16a, temp16b); + temp64 = fast_expansion_sum_zeroelim(temp32a, temp32b); + finother = fast_expansion_sum_zeroelim(finnow, temp64); + finswap = finnow; finnow = finother; finother = finswap; + } + } + + return finnow[finnow.length - 1]; + } } diff --git a/modules/core/src/test/java/org/locationtech/jts/triangulate/DelaunayTest.java b/modules/core/src/test/java/org/locationtech/jts/triangulate/DelaunayTest.java index be0b90ead6..fbee700c43 100644 --- a/modules/core/src/test/java/org/locationtech/jts/triangulate/DelaunayTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/triangulate/DelaunayTest.java @@ -63,7 +63,7 @@ public void testGrid() public void testCircle() { String wkt = "POLYGON ((42 30, 41.96 29.61, 41.85 29.23, 41.66 28.89, 41.41 28.59, 41.11 28.34, 40.77 28.15, 40.39 28.04, 40 28, 39.61 28.04, 39.23 28.15, 38.89 28.34, 38.59 28.59, 38.34 28.89, 38.15 29.23, 38.04 29.61, 38 30, 38.04 30.39, 38.15 30.77, 38.34 31.11, 38.59 31.41, 38.89 31.66, 39.23 31.85, 39.61 31.96, 40 32, 40.39 31.96, 40.77 31.85, 41.11 31.66, 41.41 31.41, 41.66 31.11, 41.85 30.77, 41.96 30.39, 42 30))"; - String expected = "MULTILINESTRING ((41.66 31.11, 41.85 30.77), (41.41 31.41, 41.66 31.11), (41.11 31.66, 41.41 31.41), (40.77 31.85, 41.11 31.66), (40.39 31.96, 40.77 31.85), (40 32, 40.39 31.96), (39.61 31.96, 40 32), (39.23 31.85, 39.61 31.96), (38.89 31.66, 39.23 31.85), (38.59 31.41, 38.89 31.66), (38.34 31.11, 38.59 31.41), (38.15 30.77, 38.34 31.11), (38.04 30.39, 38.15 30.77), (38 30, 38.04 30.39), (38 30, 38.04 29.61), (38.04 29.61, 38.15 29.23), (38.15 29.23, 38.34 28.89), (38.34 28.89, 38.59 28.59), (38.59 28.59, 38.89 28.34), (38.89 28.34, 39.23 28.15), (39.23 28.15, 39.61 28.04), (39.61 28.04, 40 28), (40 28, 40.39 28.04), (40.39 28.04, 40.77 28.15), (40.77 28.15, 41.11 28.34), (41.11 28.34, 41.41 28.59), (41.41 28.59, 41.66 28.89), (41.66 28.89, 41.85 29.23), (41.85 29.23, 41.96 29.61), (41.96 29.61, 42 30), (41.96 30.39, 42 30), (41.85 30.77, 41.96 30.39), (41.66 31.11, 41.96 30.39), (41.41 31.41, 41.96 30.39), (41.41 28.59, 41.96 30.39), (41.41 28.59, 41.41 31.41), (38.59 28.59, 41.41 28.59), (38.59 28.59, 41.41 31.41), (38.59 28.59, 38.59 31.41), (38.59 31.41, 41.41 31.41), (38.59 31.41, 39.61 31.96), (39.61 31.96, 41.41 31.41), (39.61 31.96, 40.39 31.96), (40.39 31.96, 41.41 31.41), (40.39 31.96, 41.11 31.66), (38.04 30.39, 38.59 28.59), (38.04 30.39, 38.59 31.41), (38.04 30.39, 38.34 31.11), (38.04 29.61, 38.59 28.59), (38.04 29.61, 38.04 30.39), (39.61 28.04, 41.41 28.59), (38.59 28.59, 39.61 28.04), (38.89 28.34, 39.61 28.04), (40.39 28.04, 41.41 28.59), (39.61 28.04, 40.39 28.04), (41.96 29.61, 41.96 30.39), (41.41 28.59, 41.96 29.61), (41.66 28.89, 41.96 29.61), (40.39 28.04, 41.11 28.34), (38.04 29.61, 38.34 28.89), (38.89 31.66, 39.61 31.96))"; + String expected = "MULTILINESTRING ((38 30, 38.04 29.61), (38 30, 38.04 30.39), (38.04 29.61, 38.04 30.39), (38.04 29.61, 38.15 29.23), (38.04 29.61, 38.34 28.89), (38.04 29.61, 38.59 28.59), (38.04 30.39, 38.15 30.77), (38.04 30.39, 38.34 31.11), (38.04 30.39, 38.59 28.59), (38.04 30.39, 38.59 31.41), (38.15 29.23, 38.34 28.89), (38.15 30.77, 38.34 31.11), (38.34 28.89, 38.59 28.59), (38.34 31.11, 38.59 31.41), (38.59 28.59, 38.59 31.41), (38.59 28.59, 38.89 28.34), (38.59 28.59, 39.61 28.04), (38.59 28.59, 40.39 28.04), (38.59 28.59, 41.41 28.59), (38.59 31.41, 38.89 31.66), (38.59 31.41, 39.61 31.96), (38.59 31.41, 40.39 31.96), (38.59 31.41, 41.41 28.59), (38.59 31.41, 41.41 31.41), (38.89 28.34, 39.23 28.15), (38.89 28.34, 39.61 28.04), (38.89 31.66, 39.23 31.85), (38.89 31.66, 39.61 31.96), (39.23 28.15, 39.61 28.04), (39.23 31.85, 39.61 31.96), (39.61 28.04, 40 28), (39.61 28.04, 40.39 28.04), (39.61 31.96, 40 32), (39.61 31.96, 40.39 31.96), (40 28, 40.39 28.04), (40 32, 40.39 31.96), (40.39 28.04, 40.77 28.15), (40.39 28.04, 41.11 28.34), (40.39 28.04, 41.41 28.59), (40.39 31.96, 40.77 31.85), (40.39 31.96, 41.11 31.66), (40.39 31.96, 41.41 31.41), (40.77 28.15, 41.11 28.34), (40.77 31.85, 41.11 31.66), (41.11 28.34, 41.41 28.59), (41.11 31.66, 41.41 31.41), (41.41 28.59, 41.41 31.41), (41.41 28.59, 41.66 28.89), (41.41 28.59, 41.96 29.61), (41.41 31.41, 41.66 31.11), (41.41 31.41, 41.96 29.61), (41.41 31.41, 41.96 30.39), (41.66 28.89, 41.85 29.23), (41.66 28.89, 41.96 29.61), (41.66 31.11, 41.85 30.77), (41.66 31.11, 41.96 30.39), (41.85 29.23, 41.96 29.61), (41.85 30.77, 41.96 30.39), (41.96 29.61, 41.96 30.39), (41.96 29.61, 42 30), (41.96 30.39, 42 30))"; checkDelaunayEdges(wkt, expected); }