indefinite hang on push!

[mkborregaard]

I have a weird bug where push! hangs indefinitely on a certain data set. I pushed one by one and found that any point added after the first 66 points hangs. Here are the points

x = [1.03186, 1.05258, 1.0001, 1.02761, 1.5616, 1.55889, 1.38433, 1.27905, 1.34451, 1.37068, 1.26181, 1.5115, 1.33991, 1.26601, 1.50506, 1.34784, 1.56958, 1.34287, 1.38122, 1.20729, 1.47628, 1.21273, 1.37316, 1.3211, 1.47897, 1.52702, 1.57301, 1.38779, 1.46979, 1.36821, 1.34085, 1.52872, 1.48673, 1.18732, 1.18903, 1.26833, 1.32804, 1.4833, 1.27449, 1.31384, 1.19336, 1.23589, 1.55747, 1.33041, 1.26193, 1.24199, 1.51688, 1.50591, 1.17781, 1.32349, 1.2099, 1.59301, 1.86091, 1.65146, 1.67671, 1.82353, 1.60677, 1.93322, 1.88519, 1.80052, 1.78107, 1.61865, 1.61219, 1.66273, 1.52872, 1.74341, 1.79085, 1.9901, 1.60177, 1.57926, 1.19508, 1.19577, 1.17234, 1.16354, 1.17483, 1.16712, 1.15652, 1.1951, 1.19052, 1.31302, 1.27429]
y = [1.41492, 1.40206, 1.40788, 1.40209, 1.3553, 1.24067, 1.35118, 1.3278, 1.17381, 1.31145, 1.36195, 1.25964, 1.07412, 1.05174, 1.43915,1.14804, 1.19189, 1.10733, 1.0001, 1.23957, 1.42062, 1.19433, 1.27763, 1.40607, 1.24241, 1.42738, 1.29633, 1.22027, 1.11073, 1.26872, 1.38729, 1.28611, 1.13272, 1.27253, 1.32449, 1.39246, 1.22622, 1.1121, 1.05529, 1.1763, 1.23462, 1.32394, 1.26295, 1.15142, 1.05591, 1.24609, 1.3416, 1.44114, 1.2875, 1.08118, 1.32153, 1.42507, 1.04641, 1.23174, 1.16237, 1.23365, 1.26502, 1.06616, 1.29185, 1.0097, 1.28653, 1.21709, 1.3947, 1.29859, 1.28611, 1.30817, 1.14595, 1.10608, 1.2296, 1.42348, 1.18309, 1.09683, 1.09064, 1.28974, 1.15095, 1.16269, 1.16553, 1.1083, 1.25221, 1.08263, 1.13242]

vd = DelaunayTessellation(length(x))
for i in eachindex(x)
    println("$i of $(length(x))")
    push!(vd, Point2D(x[i], y[i])) # I know I can add all the points at once, but this is better for testing
end

After easily adding all the black points, adding any of the red points will cause the session to hang/crash

[mkborregaard]

Ah. I think it happens when two points are within eps() distance from each other.

[mkborregaard]

Should there be a check that either errors or removes duplicate points with a warning, to prevent the indefinite hang?

[robertdj]

Well spotted! Unfortunately, I don't have time to look into this at the moment.
But what would be the preferred behavior? Error or warning?

[mkborregaard]

I think the preferred behaviour would be an error but with an option to thin. In my case for instance the identity of the points are important (they belong to different genetic lineages) so I'd prefer to get an error then do the thinning manually to be in control, but for most geometric uses I guess automatic thinning is preferable.

[jondeuce]

Is there any progress on this issue? I, too, have run into this problem. As @mkborregaard suggested, an error with an option to thin would be great. I would be willing to help, of course, though I'm not sure how the thinning would be implemented.

[dlfivefifty]

Here seems to be a minimal example:

julia> tess = DelaunayTessellation();

julia> push!(tess, Point2D(0.9510565162951536, 0.1)) # hangs

[robertdj]

@dlfivefifty : This example probably fails because both coordinates have to be between 1 and 2.

[dlfivefifty]

Ah ok, in that case an ArgumentError should be thrown.

[robertdj]

@mkborregaard : I'm trying out your data and also see that the computations never seem to finish. I have narrowed it down to the function findindex with a never ending loop (while true) that has to break inside the loop.

One question: How did you reach the conclusion about eps() distance? All your points seem further apart than this.

[robertdj]

@dlfivefifty : Good point. This input restriction has been discussed quite a few times :-) I'm actually not sure what the status is.

[robertdj]

@skariel : It seems that mkborregaard's point number 67 cannot be assigned to a triangle in findindex, as GeometricalPredicates.intriangle alternates between the non-acceptable outputs -2 and -3.
Can you shed lights on why we don't eventually exit?

[mkborregaard]

Ah ok, in that case an ArgumentError should be thrown.

I think the method that only takes (1,2)-interval coordinates should be an internal low-level function for connoisseurs and there should be a convenience method automatically rescaling.

[mkborregaard]

@robertdj I don't remember? But afair I had some identical points, these created errors.

[robertdj]

@mkborregaard : My moment of not having time also lasted a while :-) Regarding this bug I would need to spend quite some time to dig into the internals of both VoronoiDelaunay and GeometricalPredicates, so I cannot make any promises about a resolution.

BTW, I agree that it would be a good idea to automatically scale points. I do this in VoronoiCells, but it would make sense to move it to VoronoiDelaunay.

[natschil]

There is a pull-request that addresses the rescaling issue: #50

[jstarczewski]

@robertdj Hello, what the status of this?

[robertdj]

@natschil I had not seen your PR -- sorry! But I'm so out of sync with this package that I don't think it makes sense for me to make decisions about a PR

[bchaber]

I'm investigating the root cause of the error, with no luck so far. Looking at @mkborregaard's point set I can see that after adding the first 66 triangles, the resulting tesselation contains triangles composed of the same vertices (however in a different order). This confuses findindex as it is stuck in an infinite loop. When I manually set vd._last_trig_index before pushinig 67th point it works correctly. I've tried rescaling y-coordinates to cover almost the whole (1, 2) range but it still fails.

[natschil]

@bchaber not directly related, but I've re-written the package to deal with the issue addressed in #50 (non-convex tesselations), that kind of does rescaling (but I think poorly, there is still a bug somewhere) and also have some fixes caused by co-linear points. I would like to improve the code and document it some more, but in case you are interested, it can be found in https://github.com/natschil/VoronoiDelaunay.jl/blob/tmpbranch/src/tesselation.jl . I've cleaned up some of the old code, so it should be slightly easier to see what is going on.

The fix for co-linear points requires my branch of GeometricalPredicates.jl which can be found at https://github.com/natschil/GeometricalPredicates.jl

I have yet to make these into pull-requests because I'd like to clean up the code and document it more. However, as it has been sitting on my laptop for the past few months without any progress on my part, I thought it might make sense to share it here.

[bchaber]

Thank you @natschil, I've checked out your forks. They of course work for the data points in the first post.

[natschil]

Well, that's nice to hear I guess :) Maybe it's an issue with colinear points? I had some problems previously with triangles being added that were degenerate, and GeometricalPredicates.jl didn't return the correct return value. I convinced myself that my "fix" actually fixes the issue, but it is somewhat hacky so it might still fail sometimes.

Note that the space of "allowed" points is smaller than [1,2]x[1,2], so I do some rescaling. Basically the point is that instead of square initially, there is a triangle that point should be in. Whenever a point is added I check whether it is in the domain, and then rescale the domain appropriately. I am still not 100% happy with this step, as (a) the rescaling is much more extreme than it needs to be, probably some basic linear algebra bug somewhere and (b) the edges to the three points of the big triangle (treated as "points at infinity") may need to be updated, as I think that rescaling can result in a non-triangle being formed as the points at infinity aren't rescaled like the rest of the mesh.

What also needs to be fixed is treatment of duplicate points. This might be straight-forward.

[robertdj]

I have another 5 cents to add: The Deldir package handles the problematic point set. It seems to be because it first sorts the points in a suitable way.
Using the same sorting also allows VoronoiDelaunay to handle the points:

using Deldir
da = Deldir.DeldirArguments(x, y, [0.0; 2.0; 0.0; 2.0], 1e-9)
pts = [VoronoiDelaunay.Point2D(x[n], y[n]) for n in 1:length(x)]
pts2 = pts[da.indices]
vd = DelaunayTessellation(length(pts));
push!(vd, pts2)

The sorting is performed in the wrapped Fortran code, so I don't know how it actually works.
(Note that the DeldirArguments stuff is not exported from Deldir)