Ensure convexity of tessellation and enable handling of points outside [1,2]x[1,2]

[MiriamHinzen]

Closes #34 and #27.

Non-convexity of the tessellation means we are missing triangles. A non-convex tessellation is obtained, if at least one of the corner points of the initial square lies inside the circumcircle of a missing triangle. This problem is solved by scaling and shifting the points such that the union of the circumcircles of the triangles along the boundary of the tessellation lies inside (1,2)x(1,2). With this modification it is now possible to start off with a point set outside [1,2]x[1,2]. Apply the function scaleShiftPoints to the point set before the tessellation and the function expand to the edge vertices after the tessellation.

Here's an example:

# test point set
points = Array{Point2D,1}(undef,0);
x = [0.5 0.41 0.31 0.2 0.7]; y = [0.2 0.31 0.41 0.5 0.7];
for i in 1:5
    push!( points, Point2D( x[i], y[i] ) )
end
# scale the point set, so that the union of the circumcircles fits in (1,2)x(1,2)
scaledPoints, ranges = scaleShiftPoints( points );
# tessellation
tess = DelaunayTessellation(length(x))
push!( tess, scaledPoints )
# end points of edges, scaled, ...
a = Array{Point2D,1}(undef,0)
b = Array{Point2D,1}(undef,0);
for edge in delaunayedges( tess )
    push!( a, geta(edge) )
    push!( b, getb(edge) )
end
# ... and scaled back to the original scale
a = expand( a, ranges )
b = expand( b, ranges )
# plot
using Gadfly
l = layer( x=getx.(points), y=gety.(points), Theme(default_color=color("orange")), Geom.point )
for i in 1:length(a)
    append!( l, layer( x=[a[i]._x; b[i]._x], y=[a[i]._y; b[i]._y], Geom.line ) )
end
plt = plot( l, Coord.cartesian(fixed=true) )

This produces the following plot:

[MiriamHinzen]

I realised that in a last-minute attempt to simplify the code I introduced a bug in the circumcircleUnion function.

In the previous commit, for each edge of the convex hull the circumcircles of the edge and any of the other points is computed. Then the circumcircle with the largest radius was selected and added to the union of circumcircles. This does not result in the correct union of circumcircles, because the selected circumcircle is often rather in the inner of the convex hull of points and circumcircles that extend the most towards the outside of the convex hull are missed.

The first plot shows the output of the corrected function circumcircleUnion in green. For each edge of the convex hull, the function starts off with the circumcircle of the edge. It then checks all the other points from the point set that happen to lie in the circumcircle and computes the circumcircle of the point and the edge points - if its radius is larger than the radius of the previous circumcircle, then the circumcircle is updated. By restriction to the points inside the circumcircle we can only obtain circumcircles that cover a larger region outside the convex hull than the "edge circumcircle". If we take the union of the computed set and the inner of the convex hull of points, we end up with a superset of the actual union of the circumcircles of the triangles of the Delaunay tessellation, which is sufficient to ensure convexity of the scaled tessellation. The second plot shows the final tessellation.

[MiriamHinzen]

The original tests are now all passed. It's just my own new test that fails and I have no idea why! Can someone help with this? If I run

using VoronoiDelaunay
points = [  Point2D(-1.0563841812533212, -1.4606363138997696) 
            Point2D(0.06312982975128989, -0.48031801152366027)
            Point2D(0.1624918689993189, -0.19919450833195906) 
            Point2D(-1.5293344878962758, -0.7657808444340142) 
            Point2D(0.5319064220493406, 0.6107808132476504)   
            Point2D(-0.3670342825169435, 0.8207427582546951)  
            Point2D(-1.9797019290444364, -0.5066353099040788) 
            Point2D(-1.5811584920674324, 1.0346409888830976)  
            Point2D(1.2185165319349451, 1.4177209167374605)   
            Point2D(-1.5991536318191626, -1.3063986775765466) ];
tess, ranges = DelaunayTessellation( points );
t = locate( tess, Point2D(0.6, 0.6) )
!isexternal( t, ranges )

in the julia REPL, everything is fine. And the test, which is basically the same, fails. Why?

[natschil]

What's the status of this pull-request? I'm encountering some non-convex triangulations in practice, which is causing some code I have to fail.

[natschil]

I think there is still a bug in isexternal, as here ranges is taken to be a ranges::NTuple{4,Float64}, but only its first two values are used. Also, there should be a default for ranges so as to not break the API

[natschil]

It also seems like isexternal doesn't do what I expect it to do, it seems like triangles that are not part of the triangulation don't end up with a sufficiently large coordinate (I saw one that had something like [1.9999998, 1.99998] or so), which means that the isexternal test fails.

[dkarrasch]

What's the status of this pull-request?

It requires a careful review. 😉 Maybe you could try to use the online review interface, so that comments are directly linked to the code they are referring to.

[natschil]

This still doesn't seem to be working, I'm still getting a non-convex triangulation when using the nodes (0.0,0.0),(0.0,1.0),(1.0,0.0),(1.0,1.0),(0.99,0.5).

[MiriamHinzen]

@natschil, did you have a look at the new examples? You probably used the old implementation. Here's how it works:

points = [Point2D(0.0,0.0)
    Point2D(0.0,1.0)
    Point2D(1.0,0.0)
    Point2D(1.0,1.0)
    Point2D(0.99,0.5)]
tess, ranges = DelaunayTessellation( points )
xy = getplotxy( delaunayedges( tess, ranges ) )
l1 = layer( x=getx.(points), y=gety.(points), Theme(default_color=colorant"orange"), Geom.point )
l2 = layer( x=xy[1], y=xy[2], Geom.path )
plot( l1, l2, Coord.cartesian(fixed=true) )

[MiriamHinzen]

I think there is still a bug in isexternal, as here ranges is taken to be a ranges::NTuple{4,Float64}, but only its first two values are used. Also, there should be a default for ranges so as to not break the API

Regarding isexternal only checking the x values: This is actually sufficient. To construct the tessellation, 4 more points are added to the point set that act as a rectangular frame (2 triangles) around the point set, which ensures that a point is always inserted into an already existing triangle. The tessellation for this extended point set now includes triangles with at least one vertex being one of the 4 additional points. These triangles need to be removed again from the tessellation. isexternal checks whether a triangle is to be removed. Checking x values is sufficient, because none of the points from the original point set suffices the condition of being outside the ranges and the 4 "frame" points always suffice the condition for the x values (less than xmin or larger than xmax) as well as the condition for the y values (less than ymin or larger than ymax). So it suffices to check one of these conditions to filter out triangles with one of these "frame" points as vertices. But, as explained in the issue, by removing triangles from the "extended" convex tessellation including the 4 additional points, the tessellation may become non-convex, if the frame wasn't large enough in comparison to the point set.

Regarding default values for ranges: I already set default values in functions that are also used by the old implementation sometimes giving non-convex triangulations, which includes the isexternal function. Setting default values for ranges in the delaunayedges function, however, would only suit catching errors and can not give the right result. I would prefer the user to not even get in touch with ranges at all, by including ranges in the tessellation struct. Currently, it's up to the user to demand receiving ranges as output (and thus using the convexity-corrected implementation) here: tess, ranges = DelaunayTessellation( points ) and also handing ranges to the delaunayedges function like this: xy = getplotxy( delaunayedges( tess, ranges ) ). In a later version one should probably include ranges as a field name in the tessellation struct, so that the user does not get in contact with ranges at all.

[natschil]

My plotting function wasn't specifying ranges to iterate() (i.e had some code that used a for t in tess construct) which was causing the erroneous results shown above. I agree that it would make sense to include ranges as part of the DelaunayTesselation2D struct, this also fixes iterate() so that it works inside for loops. The current version of this pull request also doesn't work with more general AbstractPoint2D objects. I've made some modifications that address these issues and will upload them to github soon.

[MiriamHinzen]

Ah, I get it! I agree that the iterator is not working as intended. Thanks for pointing that out! I am looking forward to your commit.

[natschil]

I'm not sure this is the correct way to go about things on github, but I've now created a pull-request to MiriamHinzen:master with my changes.

[codecov-commenter]

Codecov Report

Merging #50 into master will increase coverage by 5.55%. The diff coverage is 95.23%.

@@            Coverage Diff             @@
##           master      #50      +/-   ##
==========================================
+ Coverage   71.09%   76.65%   +5.55%     
==========================================
  Files           1        2       +1     
  Lines         384      514     +130     
==========================================
+ Hits          273      394     +121     
- Misses        111      120       +9     

Impacted Files	Coverage Δ
src/VoronoiDelaunay.jl	72.16% <88.88%> (+1.07%)	:arrow_up:
src/VoronoiDelaunayExtensions.jl	97.77% <97.77%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update e4adadc...cdb69f8. Read the comment docs.

[natschil]

An issue with this pull-request currently is that it results in some very large ranges if the data has some regularity. Here is an example (inspired by a real-life use case where I want to get a triangulation of a torus)

using VoronoiDelaunay
function test_ranges()
    Random.seed!(8)
    nodes_to_triangulate = VoronoiDelaunay.Point2D[]
    nodes_in = [VoronoiDelaunay.Point2D(rand(), rand()) for x in 1:5]
    #for i in (0, 1, -1), j in (0, 1, -1)
    for i in (0,1,-1), j in (0,1,-1), (index,node) in enumerate(nodes_in)
            new_point = (VoronoiDelaunay.getx(node),VoronoiDelaunay.gety(node)) .+ (i*1.0, j*1.0)
            push!(nodes_to_triangulate,VoronoiDelaunay.Point2D(new_point[1], new_point[2]))
    end
    l = CoherentStructures.VD.scaleShiftPoints(nodes_to_triangulate)
    println(l[2])
end
test_ranges()

I get (and I hope this is reproducible on other machines) ranges (-1.2531755484857024e16, 1.2531755484857024e16, -4.367753289622386, 2.506351096971405e16). This seems too large to be numerically stable.

Edit: here is a picture of the points I'm trying to triangulate:

[natschil]

To comment on my previous comment: perhaps the scaleShiftPoints function can be used to scale (without the cirumcircle checks) the points into a suitable subset of the interval [1,2] x [1,2] after which the approach suggested in https://github.com/JuliaGeometry/VoronoiDelaunay.jl/issues/16 would solve the convexity issue?

I'm not really familiar with the insides of VoronoiDelaunay.jl though, so I cannot comment on the mathematics.

[natschil]

I'm currently working on rewriting the triangulation code to (a) return a convex result and (b) be more readable/documented. I'll push the results to a pull-request soon.

[mkborregaard]

This would be a massive improvement for the usability for this package, but appears to have stalled. What's needed for this to make it into the package?