Overlapping Delaunay triangles

[dbaston]

The following multipoint (WKB) produces overlapping Delaunay triangles:

01040000004F0000000101000000244AD52CA286504058F4A4D3AB255140010100000062359967C8556240C56A32CF906B5E40010100000011DA2800CE1B60405A908C4AC1505D400101000000264AD52CA286504059F4A4D3AB255140010100000079E1CE20AB316540E21BBCE547446240010100000040496AF647A85F4031A081A090A36640010100000010DA2800CE9B52402C4846A560A8534001010000006040034121475E4060DBCA04DC2B634001010000009CBFFCBEDEB858409F2435FB23546740010100000079D2FDCC90544A4036F0ED5FDA964F40010100000079D2FDCC9054554035F0ED5FDA96554001010000003E533FDCB333604040102E99DBBC65400101000000EE19A0276FE15640602D7D25154B64400101000000D8F176DFB5276140AFE3EDBE6B0F5C400101000000C0B69509B8575740CF5F7E5F6FDC634001010000000000000000E065400000000000C06240010100000055E103C018095F40ADF001608C0467400101000000E4CAD299530D58403F98951A548557400101000000A3BFFCBEDEB858405EDBCA04DC2B63400101000000B8BB7C33D97E5D4051E8E48F47625B400101000000AF4E6833720C5D403AC4413D17115B4001010000008C4B663885294B406D4AF2FB35274F400101000000A6E153330BE33E4012504F75F3DC474001010000004EC3A76616C65A4048403DD5CD7359400101000000D86BB8A33572434090F27AC2234C4A40010100000062359967C8556240C36A32CF906B5E400101000000AA1EFC3FE7F65740530FFE9F737B634001010000008B4B663885294B406C4AF2FB35274F400101000000D8F176DFB5276140B0E3EDBE6B0F5C400101000000B04E6833F23064403AC4413D178E614001010000008559814798985640C0EFD16624C364400101000000A1DB5D6F204E5C40DE30542AEF6D57400101000000856E77BD817853407AC350A9BCB75240010100000000000000000049400000000000004E4001010000007AD2FDCC9054554036F0ED5FDA9655400101000000A6E153330BE3444025A09EEAE6B94B400101000000E492194E616A56409C92FC7ECD495440010100000042B7BBDE405C5940BD61A854DEDB5540010100000059E103C018095F40540FFE9F737B63400101000000E3CAD299538D65403E98951A54856240010100000045C3A76616C64F4021A09EEAE6B9514001010000007BD2FDCC90D4624036F0ED5FDA9660400101000000FBFFFFFFFF7F5B400000000000C06740010100000009F32F6C480F60409FD282DAEA34664001010000000000000000805640FEFFFFFFFF3F654001010000000AF32F6C480F6040622D7D25154B64400101000000EC78BBEFDA836340EC78BBEFDA636040010100000084598147989856403C102E99DBBC65400101000000A41EFC3FE7F65740AAF001608C0467400101000000846E77BD8178534079C350A9BCB752400101000000465693662F786140633834053B3F5F400101000000000000000040604000000000004065400101000000BF5AFA4A2A96594008F32F6C488F67400101000000A6E153330BE3444024A09EEAE6B94B400101000000EB19A0276FE156409BD282DAEA346640010100000057AFE18ED6C83B401DE5F58447984640010100000042496AF647A85F40D15F7E5F6FDC634001010000007EDFA3CD48865A403D533FDCB3B36740010100000076D2FDCC90544A4032F0ED5FDA964F4001010000000EDA2800CE9B52402B4846A560A85340010100000086DFA3CD48865A40C2ACC0234CCC6240010100000037A505B5D5695D400AF32F6C488F674001010000000000000000405F4000000000000059400101000000CA6BB8A33572434087F27AC2234C4A4001010000007AD2FDCC90544A4035F0ED5FDA964F4001010000003E533FDCB3336040C2EFD16624C3644001010000000000000000003440000000000000444001010000003DA505B5D5695D40F70CD093B7F062400101000000AF4E6833F23064403AC4413D178E614001010000004DC3A76616C64F4024A09EEAE6B9514001010000009B5CA9AE46F35F409B5CA9AE46B359400101000000E292194E616A56409B92FC7ECD4954400101000000C75AFA4A2A965940F60CD093B7F062400101000000BCB69509B85757402DA081A090A36640010100000078205C32B7795C403E533FDCB3B3674001010000005C40034121475E40A22435FB2354674001010000007E205C32B7795C40C2ACC0234CCC624001010000000200000000805B400000000000C06240010100000076BCDD77EDB1644076BCDD77ED916140

[dr-jts]

Yes, because there are a lot of nearly coinicdent points, well below the robustness tolerance threshold of the Delaunay code.

For instance,
POINT (63.54755862418691 70.90471902361652)
and
POINT (63.54755862418697 70.90471902361656).

If you use the Delaunay with a distance tolerance (which can be quite small - e.g. 0.00000001), I think this avoids the issue.

[Komzpa]

Same case minimized to 5 points:

010400000005000000010100000045C3A76616C64F4021A09EEAE6B9514001010000004DC3A76616C64F4024A09EEAE6B951400101000000244AD52CA286504058F4A4D3AB2551400101000000856E77BD817853407AC350A9BCB75240010100000009F32F6C480F60409FD282DAEA346640

[Komzpa]

Minimized to other 5 points:
01040000000500000001010000008C4B663885294B406D4AF2FB35274F40010100000045C3A76616C64F4021A09EEAE6B9514001010000004DC3A76616C64F4024A09EEAE6B9514001010000008559814798985640C0EFD16624C3644001010000000000000000805640FEFFFFFFFF3F6540

[Komzpa]

Mapbox Delaunator is able to triangulate all above cases:
https://twitter.com/komzpa/status/1035480307665907713

[dr-jts]

Do you have the actual output so it can be verified that the Mapbox code worked correctly? All I see is an image, which doesn't prove correctness.

Not to cast too much doubt on the correct result. Different algorithms, or even different implementations of different algorithms, can have different robustness characteristics.

I note that there is a comment that it is not known how robust Delaunator actually is. Robustness is hard to prove, of course. I would expect to see the code doing some standard things to improve robustness, however (such as using a robust approach to computing the inCircle predicate).

[dr-jts]

I see that Delaunator uses standard FP computation for the inCircle predicate. This is known to be non-robust (see Shewchuck's extensive analysis of this). It would be interesting to inspect the actual computation on this case in detail, to see how (whether) it handles it.

I think that Delaunator uses a sweep-line approach, rather than the incremental approach of the JTS algorithm. This is highly likely to have different robustness characteristics (and may be more robust in general - if it is possible for robustness to be a percentage rather a boolean).

[Komzpa]

@dr-jts thank you for the hint that issue can be in inCircle predicate.

Replacing double -> long double fixed the case for me in GEOS.

libgeos/geos@6abe491

[Komzpa]

@dr-jts
I read through Shewchuck's code. I like their adaptive approach, even though it's broken on modern CPUs :)

I see a way this issue can be addressed in JTS, even though there are no 80-bit floats in Java.

In this publuc domain snippet:
https://www.cs.cmu.edu/afs/cs/project/quake/public/code/predicates.c
there is incircle() implementation that measures error margin along the calculation and retreats to longer implementation in case if problems:

REAL incircle(pa, pb, pc, pd)
REAL *pa;
REAL *pb;
REAL *pc;
REAL *pd;
{
  REAL adx, bdx, cdx, ady, bdy, cdy;
  REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
  REAL alift, blift, clift;
  REAL det;
  REAL permanent, 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];

  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 = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift
            + (Absolute(cdxady) + Absolute(adxcdy)) * blift
            + (Absolute(adxbdy) + Absolute(bdxady)) * clift;
  errbound = iccerrboundA * permanent;
  if ((det > errbound) || (-det > errbound)) {
    return det;
  }

  return incircleadapt(pa, pb, pc, pd, permanent);
}

While porting the whole adapt system can take ages, I see that you already have doubledouble implementation. Is there a reason you can't retreat to it in the same fashion?

[dr-jts]

@Komzpa Have to admit that I thought the JTS Delaunay code was already using the inCircleDD algorithm. I know I tried it at one point, but perhaps didn't make if the default because I didn't have a good failure case. And maybe also because it is slower than the DP version.

So it would be good to make it an option, at least. Or even the default, so that there will not be unexpected failures.

[dr-jts]

@dbaston Do you have code that validates a Delaunay Triangulation (and/or set of Delaunay Edges) ? It would be useful to have this in JTS.

[dr-jts]

In partial answer to the question of how to validate Delaunay output, one way to compute partial validation is to use Geometry.isSimple() to test whether the set of Delaunay edges has any self-intersections at non-node points.

Note that this test is necessary but not sufficient. A computed Delaunay triangulation might have no edge crossings but still contain triangles that overlap. However the simple test suffices to detect the failure in the case in this ticket.

[dr-jts]

@Komzpa I tried the TrianglePredicate.inCircleFastDD, and it allows computing the Delaunay correctly.

It would still be interesting to understand how the Mapbox code works, even though the FP inCircle presumably fails on some of the points.

[Komzpa]

@dr-jts Delaunator only attempts building new triangles from edges of current convex hull. Supposedly with this approach you will have some triangulation even with random() instead of incircle().

What can fail in this approach is initial sorting by floating point distance, though.

[dr-jts]

Okay, so Delaunator will always produce a valid triangulation, but it may not be Delaunay? (Not that it matters for many applications).

Offhand I can't see why the JTS Incremental Delaunay doesn't also always produce a valid triangulation. Would be interesting to find out. It might be possible to detect when an invalid edge was commited, and then fail at that point. The inCircleDD could be provided as an option for when the triangulation fails (similar to the way overlay works with snapping).

[dbaston]

@dr-jts sorry, I don't have anything for validating a triangulation. I just posted this ticket after seeing a bug report in PostGIS.

[dr-jts]

For future reference here's a roadmap for addressing this issue:

1. Add option to IncrementalDelaunayTriangulator to allow using a robust inCircle predicate (i.e. inCircleDDFast)
2. Add unit tests covering this failure case (for DT and inCircle predicate). This will require implementing a DT validator
3. [optional science project] determine if possible to detect triangulation failures when they occur and fail fast, so that user can choose to fall back to using robust inCircle
4. Do performance testing to see performance impact of using inCircleDDFast. If not too bad then make it the default (and keep option to use FP inCircle if performance is needed).

[Komzpa]

@dr-jts I don't quite see why do you want to make it an option when Shewchuk provides a method of quickly choosing between implementation, most of time staying on fast one.

I hope this sketch will make the idea more clear:
#311