regular grid errors

Dec 9, 2014 at 12:53 PM
I am having an issue obtaining a 3d triangulation tetrahedra set for a 3d regular grid. While the suggestion is to impart some random noise to the system, and on some level works, this is not really a suitable resolution to the problem. Delaunay triangulation functionality should include regular grid triangulation functionality. I understand this is a limitation of the library, but has anyone come to circumvent this problem as of yet, or alter the code in some way as to provide a suitable workaround?

Nick
Developer
Dec 13, 2014 at 7:48 PM
Hi Nick,

since you know in advance you are computing triangulation on regular data, do you really need to use a library that implements a general convex hull algorithm for n-dimensional data?

When the grid is regular, 8 points that define one grid cell can be easily split into 6 tetrahedra (http://paulbourke.net/geometry/polygonise/polytetra1.gif). Implementing this is rather straightforward. If you also need tetrahedra adjacency, filling in the neighbourhood isn't particularly difficult to implement either.

Moreover, using this custom approach, your algorithm will run in linear time compared to O(n log n) of MIConvexHull and you won't have to resort to the "perturbation trick".

David
Dec 13, 2014 at 8:08 PM
I very much appreciate your response.

Yes, it is the case that I could perform triangulation using a regular grid algorithm, however, that algorithm implementation would be used only in a small subset of all cases. In other words, at runtime, there is no telling whether the grid is 100% regular or only partially regular. I agree that I could work around this with some sort of conditional evaluating the regularity of the grid before triangulation, but my question was more in line with assessing if anyone had come up with a solution that makes this already robust convex hull and triangulation library more robust by solving the issue it has with regular grids.

Best regards,
nick
Developer
Dec 14, 2014 at 8:24 PM
Well, you can always use this trick:
  • Initialize Random with specific seed
  • Shift all your positions using this random instance
  • Compute the triangulation
  • Initialize another Random with the same seed and shift the positions back
Alternatively, when you implement the IVertex interface, you keep the position twice:
  • one for the triangulation which is slightly shifted in random direction (I recommend using "deterministic random" by supplying a fixed seed to the Random class) which you supply to IVertex.Position
  • and the original position
Lastly, you can try http://is.muni.cz/th/374149/fi_b that was done by one of my students. It is a C# implementation of the "Delaunay hierarchy" algorithm using "Hilbert curve sorting" for 3D triangulations. The text is in Czech, but you should not have trouble using the code (found in Projekt.zip file). It's about 4 times faster than MIConvexHull for the 3D triangulations (but only does that).

Maybe it will handle your use regular cases better. However, I haven't tested that because in my use cases the degenerate ones do not occur. Either way, for 3D triangulations it is better than MIConvexHull and you should be able to use the same "randomization trick". Let me know if you decide to use the code, I have a slightly modified version in my private repository which fixes some minor issues with his code.

David
Developer
Dec 15, 2014 at 12:40 AM
I was thinking about a bit more and there are issues with the shift trick that I previously didn't understand well enough. I.e. when you shift the points back, some of the tetrahedra might end up having zero volume. This is obviously a problem.

However, there is also good news, because I think I've found a way to correct the problem:
  • Do the "random" shift (or some deterministic one -- it can just be a function that provides points on a unit sphere that are "irregular enough").
  • Compute the triangulation.
  • Shift the points back and detect the zero volume simplexes.
  • Remove these zero volume simplexes: On the boundary can just throw them away, inside need to add additional ones.
  • The removal of an internal simplex might cause a violation of the Dalaunay property, so the triangulation needs to be corrected using "flips" (http://en.wikipedia.org/wiki/Delaunay_triangulation#Flip_algorithms). In most cases this should only need O(1) steps.
This obviously introduces a non-trivial performance cost so I guess this should be done only if the user wishes to. I could also add an interface that would control the way the points are shifted and provide a few implementations to play with.

I am also thinking if a similar approach could be used for the convex hulls. Like, each point can be tested for planarity with it's neighbours and removed.

Now just the small detail of actually implementing it :) Maybe I will get to it within a few weeks. It about time for 3.0 version anyway. Might even do it in F# just to be fancy :) [But no promises, I've got a lot of other work to do nowadays]
Developer
Dec 20, 2014 at 3:14 AM
So a new version is available. Maybe it will help fix your problems.
Dec 22, 2014 at 5:39 PM
Dave,

I apologize for not getting back to you earlier as i have been quite busy, but have just gotten around to testing the new algorithm. Your modification works like a charm given the degenerate condition. Thank you so much for modifying it! I would love to show you some images of my application of this algorithm if you're interested. please send me your email if you would like to see some snapshots.

Best regards,
Nick