VoronoiDiagramBuilder gives incorrect results in JTS 1.20.0

[james-willis]

Voronoi results incorrect in jts 1.20.0

The Voronoi diagram results are incorrect in JTS 1.20.0. It seems the results are not correctly clipped to the clip envelope.

The following code snippet demonstrates the issue:

jshell> import org.locationtech.jts.geom.Envelope;
   ...> import org.locationtech.jts.geom.GeometryFactory;
   ...> import org.locationtech.jts.geom.Coordinate;
   ...> import org.locationtech.jts.geom.MultiPoint;
   ...> import org.locationtech.jts.triangulate.VoronoiDiagramBuilder;
   ...>
   ...> GeometryFactory g = new GeometryFactory();
   ...> MultiPoint geom = g.createMultiPointFromCoords(new Coordinate[] {new Coordinate(0, 0), new Coordinate(2, 2)});
   ...>
   ...>
   ...> VoronoiDiagramBuilder builder = new VoronoiDiagramBuilder();
   ...> builder.setSites(geom);
   ...> builder.setTolerance(0);
   ...> builder.setClipEnvelope(new Envelope(-2.0, 4.0, -2.0, 4.0));
   ...> builder.getDiagram(geom.getFactory());
   ...>
   ...>
g ==> org.locationtech.jts.geom.GeometryFactory@1f17ae12
geom ==> MULTIPOINT ((0 0), (2 2))
builder ==> org.locationtech.jts.triangulate.VoronoiDiagramBuilder@255316f2
$12 ==> GEOMETRYCOLLECTION (POLYGON ((-2 -2, -2 3.9999999999999996, 4 -1.9999999999999996, 4 -2, -2 -2)), POLYGON ((-2 3.9999999999999996, -2 4, 4 4, 4 -1.9999999999999996, -2 3.9999999999999996)))

JTS version 1.19.0 gives the correct result:

$12 ==> GEOMETRYCOLLECTION (POLYGON ((-2 -2, -2 4, 4 -2, -2 -2)), POLYGON ((-2 4, 4 4, 4 -2, -2 4)))

For convenience, the code I fed into jshell:

import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.MultiPoint;
import org.locationtech.jts.triangulate.VoronoiDiagramBuilder;

GeometryFactory g = new GeometryFactory();
MultiPoint geom = g.createMultiPointFromCoords(new Coordinate[] {new Coordinate(0, 0), new Coordinate(2, 2)});

VoronoiDiagramBuilder builder = new VoronoiDiagramBuilder();
builder.setSites(geom);
builder.setTolerance(0);
builder.setClipEnvelope(new Envelope(-2.0, 4.0, -2.0, 4.0));
builder.getDiagram(geom.getFactory());

[dr-jts]

The Voronoi result has slight discrepancies from the result theoretically expected:

GEOMETRYCOLLECTION (
     POLYGON ((-2 3.9999999999999996, 4 -1.9999999999999996, 4 -2, -2 -2, -2 3.9999999999999996)), 
     POLYGON ((-2 4, 4 4, 4 -1.9999999999999996, -2 3.9999999999999996, -2 4)))

The non-integral coordinate POINT (-2 3.9999999999999996) is computed from the intersection of the following line segments (the first is the edge of the clipping rectangle, and the second is computed as the bisector between the two input points):

LINESTRING (-2 -2, -2 4)
LINESTRING (-30.484126984126984 32.48412698412698, 32.492063492063494 -30.49206349206349)

By comparing the ordinates of the bisector line, it is can be seen to be slightly off the expected angle of -45, due to the slight discrepancy in low-order decimal places:

-30.484126984126984   -30.49206349206349
 32.48412698412698     32.492063492063494 

This means that it cannot intersect the clipping square POLYGON ((-2 -2, -2 4, 4 4, 4 -2, -2 -2)) at corner points. This is visible when the topology is magnified:

This is likely caused by a change to make line intersection computation more robust by using DoubleDouble arithmetic (in #989). Ironically, in this case computing a more exact result is less accurate than the rounded result produced by standard floating-point arithmetic.

As with all floating-point computation, it's not possible to provide absolutely exact results, even in apparently simple cases where the theoretically correct is obvious by inspection. The most that can be expected is that results are within a small tolerance distance of the exact value. This is the case here. This is why VoronoiTest uses a distance tolerance.