Triangulation of regular point grids produces odd patterns

[Kevin-Mattheus-Moerman]

I'm working on a method to create unstructured quad grids like on the right here:

I start with the Delaunay triangulation on the left and join adjacent triangles when they would make nice quads when grouped. Next I split the remaining triangles into 3 quads and split each quad.

Anyway, usually I use an nearly equilateral grid of triangles as you can see on the left. To make the quads more regular I wanted to use a normal grid. However now I'm noticing there are interesting patterns in the resulting triangulation. Can these be avoided somehow? E.g. could you just do a "left slash" if the left and right slash would be equivalent numerically (both satisfying Delaunay? As you can see in the below image the connectivity is not always great. For instance there are a lot of "4-connected" nodes. Usually I filter these out actually, producing the second figure below, but in this case that makes it even weirder. It would be nice if a regular triangulation was achieved for a regular grid. Is this something that is possible?

[DanielVandH]

I have thought about something like this too in the past... it's a bit annoying.

E.g. could you just do a "left slash" if the left and right slash would be equivalent numerically

I think there is a way to always make a "consistent" choice when it comes to choosing these edges. It is described in Chapter 2.9 of this book https://people.eecs.berkeley.edu/~jrs/meshbook.html (Chapter 2 is freely available if you are curious) using some idea about "symbolic weight perturbations". I suspect implementing this is what is needed to really get this right (other than some post-processing step of the triangulation which I think wouldn't really fit in this package) and would probably choose consistent "slashes".

Unfortunately I've never really understood how to take the algorithm described in this chapter and actually implement it in practice, nor have I found a good reference that also discusses this idea (what little I could find is either an exercise in this very book, or a similar exercise in Shewchuk's (one of the authors..) lecture notes)... Trying to understand this method was actually my main motivation to even implement weighted triangulations in this package in the first place!

So, it probably is possible but for whatever reason it's out of my depth anytime I've looked at it.

[DanielVandH]

For your example at least - does it matter if you don't randomise the insertion order? It might not depending on how your points are sorted to start with, but worth looking at.

[Kevin-Mattheus-Moerman]

@DanielVandH thanks for the quick reply. I understand, and yes for my algorithm relating to the quads this issue does not matter, the groupings to quads still work. I do think there is an opportunity to add some behaviour to have more regular meshes for gridded points, or to switch the orientation if it improves connectivity (e.g. if it can avoid four connectedness in favour of homogeneous 6-connectedness, provided the triangulation is a Delaunay one both ways for instance).

FYI one "ducktape" solution is to give the points a tiny tiny shear. The shear causes your algorithm to use a single "slash" direction consistently throughout a gridded point region, see below (ignore the boundary effects). After the Delaunay triangulation is formed one can perform the inverse shear to get the proper coordinates back. I suppose the shear should be small enough to not alter the Delaunay property and large enough to be of influence to the algorithm in terms of favouring a "slash direction". If this combination is at all possible 😄

So, for the moment, I've added something like this when a user requests a Cartesian grid infill:

γ = 1e-6 # Tiny Shear 
f = [1.0 γ 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0] # Deformation gradient tensor
Vn = [f*v for v in Vn] # Deform grid

[DanielVandH]

orientation if it improves connectivity (e.g. if it can avoid four connectedness in favour of homogeneous 6-connectedness, provided the triangulation is a Delaunay one both ways for instance).

This is an interesting idea. I could try and take a look at whether that's feasible to implement (I'd have to try and remember where in the code it should be checked, since nothing at the moment in the code directly addresses point sets where the triangulation is not unique).

FYI one "ducktape" solution is to give the points a tiny tiny shear.

Good idea. You might have to do that for now unless I can actually come up with some solution. FWIW, I think this solution is like a numerical implementation of the symbolic approach I mentioned in the book (thus coming with the limitation, as you already mention, that the shear factor might make things non-Delaunay). As long as you roughly know the length scales you're working with though I imagine coming up with a reasonable shear factor is not too bad.

It might be a while before I can come up with anything for your idea to favor small-connectivity, but I'll let you know if I can find a solution.

[Kevin-Mattheus-Moerman]

It might be a while before I can come up with anything for your idea to favor small-connectivity, but I'll let you know if I can find a solution.

Sounds good. Just to be clear I meant optimised connectivity. Which might mean to increase connectivity. Very low connectivity, like 3 and 4 connected nodes should be avoided if possible. Nodes that are "over connected" e.g. 12 edges meeting at a node, are also bad so a middle ground e.g. having mostly 6-connected nodes is good.
Perhaps this can be done by doing the normal Delaunay as you do now but then do edge flips that keep it Delaunay but improve connectivity. Anyway, I think I'm out of the woods with my "tiny shear" fix so feel free to close this one.

[DanielVandH]

Ah, optimised. I see what you mean now. That would probably be hard to do without post-processing which, as mentioned, maybe doesn't fit in this package (although might be an interesting tutorial example?). I can probably write an example loop of how to do such post-processing... could be fun. Can you send me code that reproduces this figure?

[DanielVandH]

I had a loop like

for i in each_solid_vertex(tri)
    is_boundary_node(tri, i)[1] && continue # only do interior nodes
    if num_neighbours(tri, i) ≠ 6 # issue is this probably
        S = get_surrounding_polygon!(tri, i) # not each_neighbour since we would be modifying the iterator while iterating over it
        for j in S
            k, ℓ = get_adjacent(tri, i, j), get_adjacent(tri, j, i)
            cert = point_position_relative_to_circumcircle(tri, i, j, k, ℓ)
            if is_on(cert)
                flip_edge!(tri, i, j)
            end
        end
    end
end

that I threw together in a few minutes but, for reasons I'd probably see if I looked a bit longer, doesn't really fix e.g. a triangulation given by retriangulate(triangulate_rectangle(0, 10, 0, 10, 11, 11, single_boundary=true)). Maybe gives you some ideas though if you want to experiment. Feel free to close though if you are satisfied enough with the shear.

[DanielVandH]

Will close, but feel free to open if you want to discuss more.