 # Delaunayn - trivial triangle removal criteria

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.

Ok. Next attempt in trying to understand the algorithm.

If I understand the description of the algorithm correctly, there are two volumes that are important:

1. The volume of the original simplex.
2. The volume of an orthogonal simplex where the first `n` side are equal to the original simplex (with `n` being the number of dimensions).

The first volume is trivial to calculate and a reference to the source is given in the comment:
`V = abs (det (X));`

It doesn’t look like the second volume is actually calculated.
Instead, the following system of linear equations is solved:
`Y = V / L`
Like Nick wrote, the denominator `L = sqrt (sumsq (X, 2))` is a vector with the lengths of the first `n` sides of the original simplex.
So, the question is: What is `Y` in that equation?
For that, it might be a good idea to multiply both sides of that equation by `L`:
`V = Y * L`

IIUC (at least for the 2d case), `Y` is the vector that is needed to form a rectangle with `L` that has the same volume as the original simplex.
If any of the components of `Y` is small, that means (at least) one of the original side lengths was very small compared to the others. That seems to be the condition that is used for discarding it.

There is probably more to it. But it seems to make some sense for me now…

I hope this is finally somewhat correct.

Yes, I think you mostly have it. It took me a moment to figure out what the linear equation was doing, but comparing it to the 2D case, it’s just doing a parallel/vectorized version of the check. The area or volume is divided by the length of the side, then a comparison operation is performed that only returns true if all, not just any, of the comparisons are smaller than `tol`. the
`if (abs (det (X)) / sqrt (sumsq (X, 2)) < tol)` is just vectorizing the same check, as does the same division it only returns true if the vector is all ones.

so I get what it’s doing, but it’s definitely not doing a relative volume check, and it means that the value of tol does not have consistent meaning across different dimensions like they seem to assume it should.

I’m not sure if I misunderstood your comment. But the condition uses `mrdivide`, not `rdivide`.

yes, it is a matrix divide. and now that I’m playing with different values I’m realizing my chosen trivial test values made it look almost equivalent to rdivide. error on my part. But this just emphasizes that the broken out 2D version isn’t doing the same thing.

The 2D code takes the area, then does a `./` of the volume by each vector independently

The >2D code does the det(volume)/[vector lengths] operation as:
`[1x1] / [n x 1] = [1 x n]`

and the elements of the `[1xn]` output, say `V / [X1;X2;X3] = [a b c]` satisfies
`a*X1+b*X2+c*X3 = V`

this means the abc aren’t independently determined by the vectors X1,2,3.

a simple 2D test case:

``````x = [0 4 0 5];
y = [0 0 3 4];

Tri =
3 2 1
3 2 4

Volume =
12
19
``````

inside the 2D code block:

``````abs (notdet(:,3) ./ sqrt (sumsq (p12, 2))) < tol &
abs (notdet(:,3) ./ sqrt (sumsq (p23, 2))) < tol
-->
[ [2.4; 3.8] < tol &
[3.0; 4.6082] < tol]
``````

forcing the 2D code to run through the nD block:

``````X(1) =
-4 3
4 0

X(2) =
-4  3
-1 -4

abs (det (X)) / sqrt (sumsq (X, 2)) < tol
-->
X1:  [1.4634   1.1707] < tol

X2:  [2.2619   1.8652] < tol

``````

so the split out 2D and nD cases aren’t running the same test, and I don’t think either test equates to calculating a relative volume.

I was going to play around with ways to vectorize at least the 3D and maybe 4D cases since it’s fairly trivial to do vector based volume calculations in parallel like the 2D case, and maybe adjust the for loop to iterate over each dimension instead of each triangle, but I got stuck on the apparent inconsistency between the current code paths, neither of which seem to match the stated intent.

The 2d path was added here for this bug:
GNU Octave - Bugs: bug #53689, vectorize delaunayn for the… [Savannah]

I’m still trying to wrap my head around the comment explaining the algorithm:

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.

That sentence doesn’t sound like correct English to me. (That might be because I’m not a native speaker.) Could you please describe in other words what that sentence means?

doesn’t sound like correct english to me either i think what they’re trying to get across:

they use a volume constructed as a simplex with all side lengths equal to the test volume’s edge length.

i had assumed when reading it this meant they were going to generate an ‘equilateral’ triangle and (nD equivalent) from each side, and see if ‘test volume’ / ‘equilateral volume’ < tol, calling that ‘too small’ and throwing it away.

looks like that bug fix was an earlier attempt to do the same thing I’m looking at doing for >2d (as it says, the speedup is significant for large meshes) but it doesn’t provide much clarity on what it should be doing. it didn’t actually break out a compatible test, and the test doesn’t match what I think is the intent of the text.

I was sort of hoping someone might have a simple explanation of how that linear equation solve is an acceptable substitute for an actual volume/volume comparison.

I don’t know about the 2-d case. But the n-d case might be correct.

Like I currently understand it (for simplicity in 2d):
A simple sketch - not to scale:

The original simplex is the triangle with the dark blue, light blue, and brown sides.
The algorithm selects two of the sides (in the example the dark and light blue ones).
It forms a new vector (what I called `L` above) with those two sides but orthogonal to each other. That is the green dashed line in the plot.
It then finds a vector that creates a rectangle that has the same area as the original triangle. That is what I called `Y` above. In reality that vector is colinear to `L`. But I’ve drawn it perpendicular (in black) to illustrate the area.

The last part might be the most tricky one: If all of the components of `Y` are small, the simplex is discarded.
IIUC, the components of `Y` are small only if the original edges have been almost colinear (or an n-d simplex is very “thin”). But that is just a geometric argument from varying that sketch in my head…
It might be that the used algorithm is just a clever way of probing that.

Edited to reflect how `if` actually works with vector arguments…

ok, I think I follow.

so rather than calculating a volume and dividing, the output Y=V/L is the nDim components of the vector that would be needed to form that comparison volume. And the <tol is checking for any of those components to be arbitrarily small. so in all cases, you’re comparing tol with a linear measurement, not something that changes arbitrarily with dimension.

I don’t think the 2d case does the same. But it might be that it comes to the same conclusion with a different algorithm…

ok. playing through it with my test values, your interpretation of L and Y appear correct -

``````x = [0 4 0 5];
y = [0 0 3 4];

Tri =
3 2 1
3 2 4

Volume =
12
19
``````

just looking at the first triangle,

``````tri(1) = [3 2 1],
pts(tri(1) =
0 3
4 0
0 0

X1 = [0 3; 4 0] - [4 0; 0 0] = [-4 3; 4 0]   ## two vectors defining triangle (already perpendicular in this test case)

sqrt(sumsq(X1)) = [5;4]  ##  that is your `L`

V = 12

12 / [5;4] = [1.4634   1.1707]  ## your `Y`

(verifying [1.4634   1.1707] * [5;4] = 12)

sqrt(sumsq(Y)) = 1.8741

sqrt(sumsq(L)) = 6.4031

1.8741 * 6.4031 = 12
``````

soooo…

``````if (abs (det (X)) / sqrt (sumsq (X, 2)) < tol)
--> if ( L < tol)
--> if ([1.4634 1.1707] < tol)
``````

is checking whether all of the components of Y are less than tol. (makes sense for `all` since a vector or some components could parallel an axis and legitimately be zero)

now the 2D case is testing whether:

``````(Volume / side_length1) & (Volume / side_length2) < tol

just looking at triangle 1 again (ignoring the extra z-dimension zeros):

side1 = p12 = p1 - p2 = [-4 3]
side2 = p23 = p2 - p3 = [4 0]    ## (X1 in nD is just [p12; p23])

sidelength1 = sqrt(sumsq(p12)) = 5
sidelength2 = sqrt(sumsq(p23)) = 4

V / sidelenth1 = 12/5 = 2.4
V / sidelength2 = 12/4 = 3.0
``````

so for just the first triangle it’s checking:

`````` (2.4 < tol & 3.0 < tol)
``````

whether both of those are less than tol.

OK. so they’re both arbitrary definitions of a length related to the volume and simplex vectors, but they aren’t quite the same. the nD case could be made equivalent by changing

`if (abs (det (X)) / sqrt (sumsq (X, 2)) < tol)`

to
`if (abs (det (X)) ./ sqrt (sumsq (X, 2)) < tol)`

or changing the 2D case to something like

`sqrt (sumsq (abs (det(:,3) / p12(1,1:2)'))) & sqrt (sumsq (abs (det(:,3) / p23(1,1:2)')))`

(the vectorized layout is definitely not right for more than one triangle, and is a bit tricky for `mrdivide`)

my recommendation would be that we make them consistent. then that will make it easier to vectorize the nD codepaths.

converting to the current 2D approach is much easier, and easier to vectorize. But when taking it to nD you’re now comparing volume/length < tol, and the dimensionality of volume is increasing, making this a bit nonsensical. keeping the nD case with the component comparison is much more consistent, as tol is always being compared to 1D component lengths, but as a matrix operation it will be much harder to vectorize. (hmmm… maybe implementing a vectorized `pagedmtimes` will be faster…)

also stepped back through delaunay.m changes and that nD algorithm comparing components to tol is part of the original function add by D. Batemen in 2007.

https://hg.savannah.gnu.org/hgweb/octave/rev/9fddcc586065

now that we’ve concluded there is a difference between the 2D and nD code paths, and the GNU Octave - Bugs: bug #53689, vectorize delaunayn for the... [Savannah] didn’t actually faithfully break out the 2D case, I think I’ll open a bug report on this to capture some elements of the discussion and make a placeholder for any useful patches.

see bug #60818

I wrote above that `L` and `Y` are colinear. And it looks like in Octave they (currently) are.
But as a solution of the underdetermined system of linear equations, `Y` doesn’t need to be colinear with `L`.
And it looks like Matlab gives a different solution.
In Octave:

``````>> 6 / [12; 2]
ans =

0.486486   0.081081

>> (6 / [12; 2]) ./ [12 2]
ans =

0.040541   0.040541
``````

In Matlab:

``````>> 6 / [12; 2]

ans =

0.5000         0
``````

Which is obviously not colinear with `[12 2]`.

Both are valid solutions for the system of equations.

If Octave should ever change which solution it returns for `mrdivide` for underdetermined systems, I don’t think the algorithm will still work.
But that is probably a big “if”. (On the other hand: Some users might complain about Matlab compatibility…)

those two vectors do have the same length, though. so i suspect any rotation of the 0.5 vector will probably satisfy the underdetermined linear equation, and I think would produce a valid answer for this application.

that said, I was trying to think through how to make sense of the nD cases, and if there’s a simple way of making a consistent version of the 2D check to apply to them.

(also looking up a few delaunay papers to see what their sliver removal processes might be. I obviously haven’t looked at the code, but according to a StackExchange post Matlab only does ‘zero volume’ removals. not sure if it tries to define ‘relative volume’ or just looks for vol < eps^n)

I don’t think that is correct. The example might have been a bit misleading in that with it that is almost true in that particular case.

But take for example the following linear equation:
`X * [1e20; 1e-20] = 1`
Since it is underdetermined (and solvable), it has infinitely many distinct solutions. Two trivial solutions among those are `X_1 = [0 1e20]` and `X_2 = [1e-20 0]`. Those two solutions obviously don’t have the same 2-norm. Comparing the components (or the 2-norm) of those two valid solutions indiscriminately to a given value doesn’t make much sense here.

Another argument for why the chosen solution `Y` must be colinear to the vector `L` is the following identity for the scalar product:
`a*b = |a| * |b| * cos(theta)` where `|x|` is the 2-norm of `x` and `theta` is the angle between `a` and `b`.
Here, `L*Y` is used. The components of `Y` are meaningful in the used algorithm if and only if `abs(L*Y) = |L| * |Y|`. I.e., `abs(cos(theta))` must be 1. I.e., `Y` and `L` must be colinear. Otherwise, the used condition doesn’t make much sense. (See the example just given.)

Since a pre-condition for the algorithm is that `L` and `Y` are colinear, it might be possible to reduce everything to a 1-d problem which is maybe more easily to implement in a vectorized form.

It might be possible to rewrite the condition `abs (det (X)) / sqrt (sumsq (X, 2)) < tol` as `abs (det (X)) / sqrt (sum (sumsq (X, 2))) < sqrt (N) * tol` where `N` is the simplex dimension.
That is just a scalar division. So that part would be easy to vectorize (using `./`).
But I don’t know if there is an efficient way to vectorize `det`. Or if there is another way to calculate the volume of the simplices that is more easily vectorizable.

Edit: It might make sense to make that change to the current implementation in any case so that it no longer relies on (arbitrarily) selecting a colinear solution…

We could probably use the Rule of Sarrus - Wikipedia to calculate the determinant for 3-d simplices.
I’d guess that this could be vectorized quite easily.
But we’d need to check if that actually leads to a better performance than the currently used loop…

That is a special case of the Leibniz formula for determinants - Wikipedia that would work for n-d simplices. But that will most probably be less performant than calling `det` in a loop for larger dimensions.

Before realizing the mismatch for the small-volume check I had written the following for a 3D case. the vector calculation for volume breaks out very easily just like for triangles using the scalar scalar triple product:

``````else
## 3D
p4 = pts(T(:,4), :);
p34 = p4 - p3;
if (nd == 4)
tet_vol = abs (dot (p12, cross (p23, p34, 2), 2));
idx = (tet_vol ./ sqrt (sumsq (p12, 2))) < tol & ...
(tet_vol ./ sqrt (sumsq (p23, 2))) < tol & ...
(tet_vol ./ sqrt (sumsq (p23, 2))) < tol;
``````

I didn’t do a time test to the for loop because that’s when I realized the `tol` comparisons didn’t match.