Voronoi tessellation fails when four cells meet at a point.

[mbostock]

I’m not sure the subject is an accurate description of when the error occurs, but this bug is a mirror of d3/d3#1578. A small example that produces the error is:

162.625000000000000,  97.42857142857143
-12.022138010542562, 426.24997828655780
-208.632321624577860, 239.89177832355410
-164.515150982886550,  97.42857142857147

The expected result is:

However, it crashes when trying to close the cell because none of the four branches are taken to initialize the other end of the border edge.

[gorhill]

Looking into this using your coordinates. Ok, using [-400,+400;-300,+300] for the bounding box I can reproduce the infinite loop.

[mbostock]

The same example, modified slightly:

http://bl.ocks.org/mbostock/1402461a872def87504c/fd2052b3094994b789506d0f52aded4d225958ee (crashes)
http://bl.ocks.org/mbostock/1402461a872def87504c/b3a1d88e8d8c8066e770a32c824328aee1218650 (nudged to fix)

To simplify the test I moved all the points to positive x and y:

var vertices = [
  [488.000, 426.24997828655780],
  [292.000, 239.89177832355410],
  [335.500,  97.42857142857147],
  [662.625,  97.42857142857143] /* y = 97.42857142857147 fixes crash */
];

And the bounding box is [[40, 40], [920, 460]].

[gorhill]

Hmm it appears the beachsections making up the beachline for some reasons become unordered, which breaks everything. Now I have to find out why this happens.

[mbostock]

I’m just speculating here, but there are a few places that depend on exact match rather than using epsilon, such as within leftBreakPoint (e.g., rhill-voronoi-core.js:702). I tried using epsilon instead and it caused a different error, so my observation might not be helpful. :\

[gorhill]

Problem is likely here: https://github.com/gorhill/Javascript-Voronoi/blob/master/rhill-voronoi-core.js#L720-L722

Using your original set of sites, at some point I get:

 pby2 = -328.82140685798635
plby2 = -328.8214068579863

Which means:

aby2 = 4.336808689942018e-19

This is the problem with parabola and finite precision arithmetic. The farther the directrix, the more severe the loss of precision when calculating cut points (edit: this sentence was nonsense), which are used to order the beachsections on the beachline, and when this fails this break the whole diagram computation.

This case is likely to happen when using sites which uses all the available arithmetic precision -- hence your nudge worked. Possible solutions:

1. Quantize input sites to a maximum precision in order to prevent users to use up all the available arithmetic precision, like say, all digits below epsilon are zeroed. I like this solution because it allows algorithm to stay efficient, and it is easy to implement. But this means changing the values passed by the user -- not a big deal I suppose as long as the behavior is documented.

2. Rewrite leftBreakPoint() to increase arithmetic precision. I would need to read more on this (ex.: http://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format#Double-double_arithmetic), and of course algorithm become less efficient.

Edit: Regarding not using epsilon in leftBreakPoint(), this is on purpose, because we want the most accurate result possible in this case -- approximation is undesirable (hence the bug here).

[gorhill]

I don't know if it is acceptable for you while I figure the best fix... In Voronoi.compute(), replacing:

// Initialize site event queue
var siteEvents = sites.slice(0);

with:

// Initialize site event queue
var ε = this.EPSILON;
var siteEvents = sites.map(function(site) {
    site.x = Math.floor(site.x / ε) * ε;
    site.y = Math.floor(site.y / ε) * ε;
    return site;
    });

solved the issue.

[gorhill]

Here is a nice name to refer to this bug: "catastrophic cancellation" => http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html#700

[mbostock]

Quantizing the inputs to ε sounds reasonable to me! Especially since it can be done before detecting coincident points, and thus points within some ε will also be considered coincident.

[oderby]

In case you're looking for another test case exercising this precision bug, I've stumbled across another failure case which appears to be due to the same problem.

Running the following with the latest version of the library throws the error Voronoi.closeCells() > this makes no sense!. Note that these 3 points are collinear.

var sites = [
        {"x":168.77278407280934,"y":243.94568810816415},
        {"x":243.77278407280934,"y":168.94568810816415},
        {"x":311.2727840728093, "y":101.44568810816413}
    ],
    bbox = {
        "xl":86, "xr":384,
        "yt":26,"yb":324};
var voronoi = new Voronoi();
voronoi.compute(sites, bbox);

If I apply the quantization patch from above, it still fails, but now with the error Voronoi.closeCells() > Could not close the Voronoi cell! \n(See https://github.com/gorhill/Javascript-Voronoi/issues/15)

[gorhill]

Thanks for reporting, I will look into this tomorrow.

[mbostock]

FWIW, this error does not occur with the D3 port, as shown in bl.ocks.org/1402461a872def87504c:

[oderby]

var sites = [
        {"x":457.5781,"y":49.0705},
        {"x":311.669,"y":194.9796},
        {"x":386.669,"y":119.9796}
    ],
    bbox = {"xl":232.5781,"xr":529.563,"yt":-25.9295,"yb":271.0554};
var voronoi = new Voronoi();
voronoi.quantizeSites(sites);
voronoi.compute(sites, bbox);

Is another case that fails for me (whether or not I call quantizeSites).

[gorhill]

Ok, the "problem" here is the vertex is so far away that the arithmetic in Liang-Barsky breaks when it tries to clip:

va:
  x =  426594690820137500
  y =  426594690820137200

vb:
  x =  0
  y = -191.68939999999998

Which causes Liang-Barsky's t0 and t1 to end up as 0.9999999999999982 and 0.9999999999999996, which means important digits have been lost and the resulting clipped line unexpectedly ends up not being connected to the bounding box.

[gorhill]

Best fix I can think of so far is to constraint "world" domain to 1/ε, when a Vertex object is created:

Voronoi.prototype.Vertex = function(x, y) {
    this.x = Math.max(Math.min(x, 1/ε), -1/ε);
    this.y = Math.max(Math.min(y, 1/ε), -1/ε);
    };

[oderby]

95bfbd5 doesn't fix it for me :/ I can confirm that your constraint to 1/epsilon is occurring, but I still get an error of Voronoi.closeCells() > this makes no sense! with both cases I provided above.

[gorhill]

Ugh, I made silly changes after testing.

But then, fixing the silliness doesn't fix the issue if I use exactly your bbox (originally I just used the default one since I could reproduce the bug with it).

[gorhill]

Even when constraining vertices to within [1/ε,-1/ε] domain on both axes, the Liang-Barsky arithmetic can loose enough precision to fail to connect to one the edges of the bbox.

[gorhill]

I can get your last case to work using (1/ε)/100 but this is getting way too arbitrary. It makes sense to constraint the domain to 1/ε, but for a sensible further fix I need to look at the clipping code: if a line was intersecting with the bbox, then it is a fact that the line must connect to the bbox where it was intersecting as a result.

[jesta88]

Ported the library to C# and getting the same error here. Seems to happen no matter what domain (large or small) I use.

[s4l4x]

Is the C# port publicly available? Thanks!
.xX

On Feb 20, 2014, at 8:12 PM, jesta88 notifications@github.com wrote:

Ported the library to C# and getting the same error here. Seems to happen no matter what domain (large or small) I use.

—
Reply to this email directly or view it on GitHub.

[jesta88]

It will be pretty soon. I'm testing it with Unity first and as soon as this bug is fixed, I'll make it available.

[gorhill]

Ported the library to C# and getting the same error here

Err, more... details? Which set of values are causing the errors? Are these values causing the same errors using this js version? I can't investigate this code here without any specifics.

[jesta88]

Sorry, I guess I was expecting you to magically fix all my problems. Please excuse my stupidity.

I actually took the time to really understand the clipping method and turns out I only had to put Epsilon to a higher value (0.1 for example) to fit my particular site coordinates.

If anyone's interested, the Unity port is here (wait for commit #4 though): https://github.com/jesta88/Unity-Voronoi

Thanks for your great work.

[gorhill]

I only had to put Epsilon to a higher value (0.1 for example) to fit my particular site coordinates

I still would like to know which set of points gave you problem, to find out where the algorithm fails.