Vertices of triangles as indexes of original points

[dfdx]

Is it possible, given instance of DelaunayTriangle, to retrieve indexes of original points, making its vertices? For example, if we added 20 points, and first triangle is made by points (3), (5) and (11) (in order of pushing them), is it possible to get these indexes?

[dfdx]

I ended up creating my own point type with built-in indexing:

type IndexedPoint2D <: AbstractPoint2D
    _x::Float64
    _y::Float64
    _idx::Int64
    IndexedPoint2D(x, y, idx) = new(x, y, idx)
    IndexedPoint2D(x, y) = new(x, y, 0)
end

getx(p::IndexedPoint2D) = p._x
gety(p::IndexedPoint2D) = p._y
getidx(p::IndexedPoint2D) = p._idx 

This, however, is a kind of a hack and doesn't work well in all circumstances. In particular, with this approach tesselation normally includes one fake point with zero index (created by default 2 parameter constructor).

Anyway, currently I'm using following code to get list of trignales as a matrix of size (num_triangles, 3), where each row represents 3 indexes of its vertices:

# convert matrix of shape (npoints, 2) to an array of IndexedPoints-s, 
# scaling everything to allowed interval and returning scaling factor
function to_points(shape::Matrix{Float64})
    min_value, max_value = (minimum(shape), maximum(shape))
    sc = (max_value - min_value) / (max_coord - min_coord)
    sc_shape = (shape .- min_value) / sc + min_coord
    points = IndexedPoint2D[IndexedPoint2D(sc_shape[i, 1], sc_shape[i, 2], i)
                            for i in 1:size(shape, 1)]
    return (points, sc)
end

function gettriangles(shape::Shape)
    npoints = size(shape, 1)
    tess = DelaunayTessellation2D{IndexedPoint2D}(npoints)
    points, _ = to_points(shape)
    push!(tess, points)
    return collect(tess)[1:end-1]   # excluding extra point       
end

Since I don't have any more time to spend on this issue and such (hacky) fix works for me, I'm leaving it as it is, hoping that some brave mind will be able to make it stable in all circumstances.

[skariel]

This is a bug, I'll fix it- your solution should work just fine, And there shouldn't be a fake point with an index of zero.

So I'll fix it, remember to change your code when you upgrade (so you don't miss a real point)

[dfdx]

Great, thank you!

[skariel]

The problem is that the user tessellation is built inside a preexisting tessellation of 4 points at (1,1),(1,2),(2,1) and (2,2). These are "virtual" points and need to be removed and currently this library does not remove these correctly.

I'm fixing this but it will take a bit more time since its not a trivial fix -- there should be a post-processing pass.

On the good side, I've improved performance by a bout 15% by removing a few type instabilities at GeometricalPredicates. And more performance is possible when using in-place BigInt operations. I'm working on a PR for that...

[dfdx]

Any news on this issue? I'm planning to release library that uses Delaunay triangulation, so it will be nice to know status of the bug and either fix it or warn users about some "strange" behavior.

Note, it's not critical for that library, so no hurry here. It's just for keeping users informed.

[skariel]

No proper fix yet, but maybe you could just filter out those triangles containing the virtual point (those with index 0) can this solve the problem? If not, then I guess best thing as of now would be to warn the users...

[skariel]

Actually filtering these triangles would still sometimes give weird results. Another approach for a "fast" fix would be to make a dense "frame" of points and let the user add points only inside. Again, not a real solution good for all circumstances but maybe fits you need.

I do intend to fix this issue...

[dfdx]

I already have kind of a fast fix for this, so I will go with it for now and just warn users until this issue is resolved normally. Thanks for your time!

[cortner]

Can I just confirm something: when I create a tessellation and then plot it, the reason what I see is not convex are the virtual points at the corners?

[skariel]

@cortner can you attach an example plot? I think Delaunay tessellations should not be convex, and Voronoy ones just fill the whole space... so I don't know what you mean. Maybe in your case the corners so affect the shape that you expect, it's possible. Just work further away from the corners where they don't affect the result

[cortner]

@skariel i am pretty sure the Union of triangles should be the convex hull of the vertices. I can attach a figure later.

[skariel]

I know a theorem that says that a Delaunay tessellation is the convex hull of the n+1 dimensional elevation for e.g. a tessellation in 2D has points (x,y) and 3D points (x,y, x*x+y*y) in this case the projection of the convex hull of the 3D configuration is the 2D tessellation. But the 2D tessellation does not need to be convex. Imagine using just 3 points... the Delaunay tessellation would be a triangle which is of course not convex...

[cortner]

@skariel

I can only imagine that we are having a clash of terminology, you are probably referring to the set of edges, whereas I am referring to the union of the triangles (incl interior), I guess I should have said triangulation instead of tesselation?

Anyhow, checking on Wikipedia, the Delauynay triangulation triangulates the convex hull of the set of vertices.

[cortner]

@skariel

In trying to debug my example to post it here I ran into another problem. If you think this should go somewhere else, I can delete the post and re-post as a separate issue, (but possibly it is just due to a misunderstanding of mine).

using VoronoiDelaunay
import PyPlot
X = 1.0 + rand(2, 10)
tess = DelaunayTessellation()
array2points(X) = [ Point2D(X[1,n], X[2,n]) for n = 1:size(X, 2) ]
push!(tess, array2points(1.0 + rand(2, 10)))
for edge in delaunayedges(tess)
    x1 = geta(edge)
    x2 = getb(edge)
    PyPlot.plot( [x1._x; x2._x], [x1._y; x2._y], "b-" )
end
PyPlot.plot(X[1,:][:], X[2,:][:], "r.", markersize=20)

Here is the output:

The red vertices don't co-incide with the triangulation vertices ?! where is my mistake?!

[skariel]

@cortner you are pushing random points into the tessellation instead of points from X:

push!(tess, array2points(1.0 + rand(2, 10)))

other than that, you are right about Delaunay tessellations being convex. Indeed the function delaunayedges does not iterate over edges going to the corner "external" points. This was the case for the voronoiedges method too but in master it's changed to iterate over everything

[skariel]

@cortner here is what it looks like with my code:

and the code:

[cortner]

sorry - that is embarrassing, and thank you for catching it. I was just about to post the correction.

In your picture you see exactly my issue. The second point from top in the top left corner is a re-entrant corner. This is possible for constrained Delaunay triangulations, but not for actual Delaunay triangulations. I ran another example; you can see the result below.

VoronoiDelaunay

Qhull (or Triangle gives same result)

[cortner]

@skariel I am pretty sure what happened is that you removed all triangles that have the corners of the cube as one of their vertices. I am also pretty sure that the result is therefore not a Delaunay triangulation.

[skariel]

@cortner it's just the iteration function delaunayedges that does not iterate over the corners, and originally it was on purpose. I already changed this for Voronoi edges, I'll change it now for delaunay too.

The tessellation includes the corners too, I didn;y implement yet a method to remove points hence I'm not removing the corners

[skariel]

@cortner here's my result when commenting out the single line that is responsible for skipping the corners. As you can see this is a perfectly good Delaunay tessellation:

(the line is isexternal(tr) && continue in the delaunayedges function:

[cortner]

@skariel Unfortunately, this is not useful for me - since I need to know that I am getting the convex hull of my vertices. So if I understand correctly then you can fix this by implementing a method to remove points from the triangulation. open an issue?

[skariel]

@cortner yes please open an issue that's an important missing feature!

[skariel]

@cortner actually removing the corners is much easier than removing a general point, so there are 2 mising features here:

1) remove a general point -- this is a bit hard 2) remove only the corners -- this is easy

[skariel]

@cortner I opened an issue

[skariel]

@cortner the most simple solution to this is just scale the points to an inner square, like in the image I attach below. Probably I'll fix this whole issue by just changing mincoord and maxcoord values:

[skariel]

@cortner actually, this doesn't work perfectly. I'll just have to implement point removal the hard way...

please don't use this method. Sorry for that

[tomasaschan]

@cortner If you only needed to get the complex hull of a set of points (and not the triangulation of the interior), you could take a look at my new package (it will be available by Pkg.add("QuickHull") once this PR merges, but for now you'll have to use Pkg.clone()).

[cortner]

I am actually after the triangulation. I just want to make sure I get the Delaunay triangulation and not something that I am not sure how it's a defined.

I now actually think that I can make VoronoiDelaunay produce exactly what I want - I jut haven't had time to make the switch.

[juliohm]

This package is great, but it definitely misses the option of indexed points. Very often after the tesselation one needs to post-process the result, discard points and edges for example.

I am gonna try the hack proposed by @dfdx , but it would be great to have this as a main feature of the package.

[robertdj]

Admittedly, I haven't read the whole thread, but consider the https://github.com/JuliaGeometry/VoronoiCells.jl package for working with indexed points.

[juliohm]

Hi @robertdj, I have a solution already with the proposed IndexedPoint2D from @dfdx. My comment was just to emphasize that very often we need to post-process the tessellation using the indexes. It is strange that the package doesn't support it by default. Actually it would be much more simple in my opinion to use raw Julia Arrays for storing vertices and edges instead of this AbstractPoint2D wrappers, they don't help.

[robertdj]

I completely agree with

we need to post-process the tessellation using the indexes

This is one reason why I made VoronoiCells :-)

I don't know if this was @skariel 's motivation, but raw Julia Arrays are/were not the most lightweight data type in Julia.

[gabor1]

Is there an easy way to get indexable Delaunay triangles? what I need to do is to do a tessellation, then pick a new point, and find out which three points (as indices into the original generators) form the triangle that my new point is in.

[cortner]

that should be standard - but if VoronoiDelaunay doesn't do it, then as a fallback scipy.spatial.Delaunay for sure has this option via tri.find_simplex.