I’ve been playing with some triangulation code, and got sidetracked by a
FIXME: Vectorize this for loop... block at the end of
delaunayn.m (line 140). The point of the block is to identify and remove trivially small (or even zero volume) elements created in the space by calculating the volumes of all created figures (area for 2D triangle, volume for 3D, and ‘volume’ equivalent for >3D) then:
Any simplex with a relative volume less than some arbitrary criteria is rejected. The criteria we use is the volume of the simplex corresponding to an orthogonal simplex is equal edge length all equal to the edge length of the original simplex. If the relative volume is 1e3*eps then the simplex is rejected.
this makes sense. however, the code does not calculate a relative volume. the 2D case is already broken out and vectorized for fast computation, and it’s easy to see that after the areas of the triangles are calculated, those areas are separately divided by the two lengths of each side used for area calculation, and the two results are compared to the tolerance. This is not a relative area, it’s two separate area/length comparisons
for the >2D cases, the same thing happens. the “volume” is calculated, then that “volume” divided by each side length is compared with the same tolerance. so for 3D it’s a volume/length, for 4D its (4Dspace/length). In no case is a relative ‘volume’ being compared against the tolerance.
I would think each would be divided by some product of the side-lengths to give a dimensionally appropriate ‘relative volume’. otherwise it just doesn’t seem to make sense.
Thoughts on this? not sure if anyone here was involved in writing the delaunay code or if there’s a reference that provides this as a justified approach.