remove points on corners

[skariel]

see #6

[cortner]

just leaving a comment so I get an update. (FYI, This is one of the issues why I don't use VoronoiDelaunay yet.)

[tomasaschan]

I have a need for a meshing algorithm like this one, and started investigating if this package would fill my needs. It seems wonderful, except for this issue :)

I'm almost entirely unfamiliar with the mathematics behind the algorithm here, but as far as I can understand from reading through #6, the problem arises because of these two facts:

• A Delaunay tesselation requires at least three points (otherwise there's no way to form a triangle)
• This package works around this by introducing "virtual" points at the corners of the valid interval.

The workaround has a couple of drawbacks, mainly that the introduction of these virtual points results in the tesselation being a different one than if one had generated one for all inner points without the corners. (From what I understand, the only difference is that a couple of boundary edges are missing, but still - it's a different result.)

Instead of just enabling removal of corners, I think a better approach would be to handle the initial steps of generating the tesselation differently. I suggest the following approach:

1. The virtual corners are not included. (IIUC this might also yield the valid domain for all user-inserted points to [1,2]^2 rather than [1+eps(), 2-2eps()]^2...?)
2. For the first point, no edge-adding logic is run at all, since there are (obviously) no other points to add edges to.
3. For the second point, an edge is added between the first and the second.
4. For the third point, edges are added to the both previous, thus forming the first triangle.
5. For the fourth point onwards, the normal algorithm applies.

This would give the package a more intuitive behavior; until I have four or more points, which edges should exist is pretty straight-forward, as there aren't any real choices (even if those aren't a strict Delaunay tessalation - again, I'm not familiar with the maths).

The main drawback with this approach is that it would require logic for adding a point outside of the current convex hull, which I don't know if it exists yet. I can certainly help with implementing the steps above, but I'd need help to understand how to add new points outside (if it requires any changes).

[skariel]

@tlycken see how the VoronoiCells package removes corners in this document.

Basically, you can use an inner square region 1.25 .. 1.75 (instead of 1 .. 2) then all connections between your nodes are real, and you can just ignore the four outside nodes when iterating.

I'm thinking maybe this is the easiest way to fix this issue, just changing the allowed region to use by changing the values of maxcoord and mincoord. It will also be the fastest of course. Code that uses these values could go on essentially unchaged, but it depends on how they currently treat iteration. Anyway, the fixes should be easy

Thoughts on this?

[tomasaschan]

Well, that might work for Voronoi cells, but IIUC it won't for Delaunay (which is what my use case needs). Maybe I'm wrong, but won't the Delaunay tessalation still miss the outer edges?

[tomasaschan]

Btw, @cortner: You can click the "subscribe" button to the low right and get notifications without having to comment :)

[skariel]

@tlycken if Voronoi is good Delaunay has to be too, and vice versa. They can be built from each other, this is how this library works: the internal representation is Delaunay triangles. So this should work, you can try and test this

[skariel]

and of course keep me updated ho it works for you :)

[tomasaschan]

I think my use case could be enlightening on why the current implementation doesn't really work the way I think it should. I'm working on implementing linear interpolation of an irregular 2D dataset, and my approach was going to be the following:

1. Take the (x,y) coordinates of the data and construct a Delaunay tesselation from them.
2. For each corner in the tesselation, store the corresponding z-value in a Dict{Tuple{Float,Float}, Float}, i.e. (x,y) => z. (Side note: due to the domain of this library being constrained to [1,2]^2, I rescaled (x,y) and used the rescaled values as keys.)
3. To interpolate, use locate to find the corners of the triangle which contains the given (x,y), look up the corresponding z values, and use some vector identities to calculate the equation of the plane through the three corners. Then use that to get the sought z-value.

This approach fails with the tesselation produced by this library, because the triangles close to the corners include the corner point, for which I have no z value. This might be easier to illustrate with a figure:

The blue lines are the tesselation edges returned by getplotxy(delaunayedges(tess)) with ten randomly picked points in the domain. When looking for the triangle containing the red point (1.225,1.6), locate returns the green triangle - as you see, including the corner point. What I would have wanted, would be the magenta edge (and, for correctness, one of the black-dashed ones) to be included in the tesselation, so that the triangle containing the red dot would consist of three corners from the original set of points. (I don't seem to be the only one expecting this.)

As you see, I'm interested in the tesselation mainly as a means of systematically selecting which three points to use for the interpolation; voronoi tesselations are of no value in this (since their edges don't hit the original points), and the approach taken in this library to ensure enough points is effectively destroying its usefulness for me.

I'm not saying the tesselations produced by this library are incorrect for the set of points it considers. I'm saying the library considers the wrong set of points. I'm also suggesting a means of instead considering the right set of points :)

[skariel]

@tlycken in the example above you didn't follow the suggested solution. Your points are still outside of the square 1.25 .. 1.75. If you put the points inside and ignore the corners all connections are what you need and expect. Try this, do the same plot but scale the points a bit further, you'll see it works

[skariel]

actually, let me draw this, just one minute

[skariel]

here, see image below, it works just perfect:

[skariel]

ahhh wait a sec... Im wrong.

[skariel]

OK solved it -- 1.25 .. 1.75 is just not enough. Try 1.45 .. 1.55 like in the image below. Now everything is perfect:

[tomasaschan]

But is this always good enough to get a convex hull? Is it provably so? Or is it just good enough to hide from plain sight the fact that it isn't a convex hull anymore? :)

(I'm asking out of honest curiosity, not out of spite...)

[skariel]

I think it's always good to create a convex hull... I guess this can be proven mathematically, and probably I'll look into it but for now I run several configurations and all are convex

[skariel]

also -- I'll mathematically find the largest region in which points are safe to being inserted in, as a prerequisite for really solving this issue

[skariel]

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]

Well, we might not be able to remove points; it might be enough to just not start doing the actual edge flipping until the user adds the fourth point. Then you have to handle adding external points instead, but I don't know which is hardest to implement.

[cortner]

@tlycken if you need a temporary replacement, have a look at fem.jl

My memory is this is just a Qhull wrapper. I am hoping to replace this eventually with VoronoiDelaunay

[tomasaschan]

@skariel I've spent some time reading up on the maths, and it seems it's valid to just remove the corners (and all edges to them) if you start with a triangle rather than a rectangle. In other words, a simple fix might be to do two things:

• Remove the top corners from the starting state, and instead add a single point at (1.5, 2), so that the initial state before the user starts adding vertices is a triangle with corners at (1,1), (1.5, 2), (2,1).
• Change min_coord and max_coord to be the largest region which fits inside this triangle, for example min_coord, max_coord = 4/3, 5/3 (those are the mathematically largest available; to avoid issues with float rounding, we might want to go make the required box just a little smaller).

Now, since the boundary is just a triangle and not a square, the diagram with the "virtual" vertices removed will still be a Delaunay diagram.

[skariel]

That's interesting... I did't know this. This might just be the easiest way to fix this issue. I'll implement it when I have some time

Thanks!

[mschauer]

I just ran into this, just affirming that there is still interest in this.