Created 2005-01-03 Modified 2009-04-11
Chelton Evans
Chelton's Convex Mesh Algorithm : Robustness and Numerical Stability
Chelton's Convex Mesh Algorithm: Home
Intro
Hilbert Curve
Mesh Topology
Mesh Minimizer Stability
Travelling salesman problem
My data has been from random data sets that have not been challenging.
See Adding a Point Inside the Convex Mesh for a simple discussion on how the algorithm addresses some of these issues.
I do not know enough, but I acknowledge that this is one of the major components of algorithm design and that this aspect should be looked further into.
The incremental algorithm can be made very robust by randomly inserting points and building the mesh in a random way. This method very good and simple to implement, just randomly sort the points before feeding them to the algorithm.
The orientation or half space tests can fail. For example when the determinant is close to zero or with half spaces the dot product is close to zero.
Jonathan R. Shewchuk [5] uses arbitrary precision arithmetic for performing these tests. While I have not implemented this it is clear that this is a strategy for improving the algorithms correctness. I will however point out that this error may alternatively be managed in having a consistent mesh where the error may exist but the mesh is still consistent and mesh computations can still occur.
The algorithm was designed to maintain a consistent mesh which may give it more stability than other algorithms. See Mesh Consistency.
Essentially the border between simplexes may be non unique. But the mesh must maintain the linked simplex data structure. This ensures that the operators on the mesh will work, even for pathological situations where a line of three points has been made a simplex the mesh should still be able to be transversed so other algorithms can still operate on the mesh.
The choice of the minimizer and how it is implemented effects the meshes stability.
First we can separate the construction from the minimizer. Then we can change the minimizer and use different minimizers on the data set.
I found that adding the DT minimizer to the same data set breaks the tessellator sooner because calculating the DT leads to infinite recursion. For example in 2D if four points lie on the same circle then there are two possible triangle arrangements and DT will go into infinite recursion alternating them.
However this problem is totally avoided with non-recursive minimizers. ie infinite recursion does not happen. So there is a practical application for them in that the degeneracies of DT with respect to infinite recursion can be completely avoided while still producing a good grid with the same circle or sphere criteria.
Extension to DT to eliminate infinite recursion: The DT minimizer could have a counter on the process stack so that each time a simplex is processed the counter is increased. At the start of a new point inserted into the grid let the counter is set to 0. If it goes above some threshold value assume that infinite recursion has occurred and process the stack differently to solve the recursion problem.
So with extra work the minimizers using infinite recursion can also be made as robust as the non-recursive ones.
The order of points inputed into the incremental algorithm can effect the numerical stability of the algorithm. For example if the simplexes generated are degenerate (one dimension less than the space eg line in 2D, triangle in 3D).
I became aware of this when I implemented pre-sorting of the points to improve the complexity of the algorithm. While the algorithm ran faster the algorithm became more unstable and handled fewer points before crashing because of the simplexes being generated. What is worse is that the more points that are used the worse things become because the points are ordered on one axis but increase in the other axes.
A solution to the numerical stability is to order the points as a Hilbert curve and insert them into the algorithm. Unlike an ordering of points on an axis this generates simplexes which are far less likely to be degenerate and also has the property that the previous point is near the next point so constant insertion time in the mesh is ensured.
This idea was suggested in computational geometry notes written by Jack Snoeyink. Unfortunately the notes disappeared from the web. Also this was a suggestion - cache friendly as the property gives us the necessary locality to do fast computations. I have been working the details out myself.
Given n points order along the axes. Assume that the exact number of points to make the Hilbert curve is given, as we can always insert dummy points to make it to the next power of the Hilbert curve, and ignore the dummy points on insertion.
Find a formula to iterate over the curve. This is an integer function in the individual axes. Consider a Hilbert curve in 2D with 16 points.
Looking up the integer sequence at The On-Line Encyclopedia of Integer Sequences in one of the axes gave the following result.
ID Number: A059252
URL: http://www.research.att.com/projects/OEIS?Anum=A059252
Sequence: 0,0,1,1,2,3,3,2,2,3,3,2,1,1,0,0,0,1,1,0,0,0,1,1,2,2,3,3,3,2,
2,3,4,5,5,4,4,4,5,5,6,6,7,7,7,6,6,7,7,7,6,6,5,4,4,5,5,4,4,5,
6,6,7,7
Name: Hilbert's Hamiltonian walk on N X N projected onto x axis: m(3).
Formula: Initially [m(0) = 0, m'(0) = 0]; recursion:
m(2n + 1) = m(2n).w(2n).f(w(2n),2n).c(m(2n),2n + 1);
w(2n + 1) = w(2n).f(m(2n),2n).f(m(2n),2n).mir(w(2n));
m(2n) = m(2n - 1).f(w(2n - 1),2n - 1).f(w(2n - 1),2n - 1).mir(m(2n - 1));
w(2n) = w(2n - 1).m(2n - 1).f(m(2n - 1),2n - 1).c(w(2n - 1),2n);
where
f(m,n) is the alphabetic morphism i := i + 2^n [example: f(0 0 1 1 2 3
3 2 2 3 3 2 1 1 0 0,
2) = 4 4 5 5 6 7 7 6 6 7 7 6 5 5 4 4];
c(m,n) is the
complementation to 2^n - 1
alphabetic morphism [example: c(0 0 1 1 2 3 3 2 2 3 3 2 1 1 0 0,
3) = 7 7 6 6 5 4 4 5 5 4 4
5 6 6 7 7]; and mir(m) is the mirror operator [example:
mir(0 1 1 0 0 0 1 1 2 2 3 3 3 2 2 3) =
3 2 2 3 3 3 2 2 1 1 0 0 0 1 1 0].
Example: [m(1)=0 0 1 1, m'(1)= 0 1 10] [m(2) =0 0 1 1 2 3 3 2 2 3 3 2 1 1 0
0, m'(2)=0 1 1 0 0 0 1 1 2
2 3 3 3 2 2 3]
See also: See also the y-projection, m'(3), A059253.
Sequence in context: A048466 A096838 A096007 this_sequence A030620
A110764 A107901
Adjacent sequences: A059249 A059250 A059251 this_sequence A059253
A059254 A059255
Keywords: nonn
Offset: 0
Author(s): Claude Lenormand (claude.lenormand(AT)free.fr), Jan 23 2001
Let f(k) be the above function, then the Hilbert curve is given by (f(k),f(k+16)) and the first point of the curve is k=0 gives (0,0).
Apart from being a complicated function unfortunately it does not look like the function is computed in constant time because of the recursive calls. Hence it may not be suitable.
As yet I have not implemented the Hilbert curve to improve numerical stability.
Assuming that we have access to Travelling salesman problem algorithms then we can use these to tessellate. For example find a minimal path through all points and feed to the incremental algorithm.
The advantage of a TSP algorithm is the potential to be "random" from the perspective of algorithm stability. For a practical example with tough data sets if the incremental algorithm doesn't break compared with the sorting algorithms then TSP are at advantage when stability is more important than the tessellation speed.