[OpendTect_Developers] Avoid the "On the Edge" Case in Geometrical Computing

[Jaap Glas]

Dear fellow developers,


How to avoid the "On the Edge" Case in Geometrical Computing
============================================================

Software packages like OpendTect internally have to do quite some geometrical
computing to manipulate their data. A recurring question is what to do in
case something is exactly "on the edge", for example when checking whether
a point is inside or outside an area or volume. What you often see is that
programmers start treating "on the edge" as a special case, and that usually
is where the shit hits the fan. Because of numerical instability, there is
little consistency in which points are exactly on the edge and which are not.
It is tempting to "solve" this by introducing some kind of "epsilon"-margin,
but in fact this makes things even worse. The complexity of the algorithm
keeps increasing, and it is not unlikely that finally there are still
configurations in which the algorithm fails, especially for regularly
spaced input data.

I plea to get rid of "on the edge" cases at all. If something is not in,
it is out. This requires a careful design of basic geometrical functions,
with extra attention to the effects of numerical instability, but it will
yield algorithms with less complexity and more stability. I will try to
show this with an example from 2D Delaunay triangulation. The whole code
can be found in src/Algo/delaunay.cc.

The following geometrical issue has to be solved again and again in
2D Delaunay triangulation. Given a current set of adjacent triangles
and one particular point, in which triangle is it located? If this point
happens to be on an edge, I do not mind to which triangle it is assigned,
as long as it is always assigned to the same one. It must be prevented
that numerical instability causes any point to be assigned to either
multiple triangles or no triangle at all.

The first requirement is to order the corner points of any triangle UVW
in the same mathematical sense, either clockwise or counter-clockwise.
If this is not already the case from the start, the sense can be checked
by the sign of the determinant det(V-U,W-U). The second requirement is
that corner points shared between triangles originate from the same source.
This guarantees that their numerical representation is exactly the same.

Life is made easy by first checking whether a point exactly coincides
with one of the corner points of any triangle in the set. By definition,
it will be assigned to the first triangle for which this is the case.

A point will be inside a counter-clockwise triangle UVW if and only if
the point is left of directional line UV, left of directional line VW, and
left of directional line WU. So the basic geometrical function to implement
can be defined as the boolean isPointLeftOfLine( point P, from A, to B ).
Now the third requirement is the numerical stability of this function.
It must be guaranteed that for any triplet (P, A, B) with P!=A, P!=B,
and A!=B, the boolean result negates if A and B are exchanged.

isPointLeftOfLine( P, A, B ) != isPointLeftOfLine( P, B, A )

This function will be called twice for every edge AB in the triangle set.
Because all triangle points are ordered in the same mathematical sense,
once with A and B in order, and once with A and B exchanged.

The implementation of isPointLeftOfLine( P, A, B ) is listed below. The
SIGN of any of the two calculated determinants det(P-A,B-A) and det(A-B,P-B)
is already sufficient to decide on which side of the directed line the point
is located. The essential observation to meet the third requirement is that
neither of their MAGNITUDES is invariant to exchanging A and B. This becomes
a source of numerical instability if point P is almost on the line. However,
the SUM of both determinants happens to be invariant to exchanging A and B.
The reader may verify this himself/herself by symbolically evaluating the
two determinants.

===============================================================================
static bool isPointLeftOfLine( const Coord& p, const Coord& a, const Coord& b )
{
    const double det_ap_ab = (p.x-a.x)*(b.y-a.y) - (p.y-a.y)*(b.x-a.x);
    const double det_ba_bp = (p.y-b.y)*(a.x-b.x) - (p.x-b.x)*(a.y-b.y);

    // |sum| is numerically invariant to swap of a and b, only sign toggles.
    const double sum = det_ap_ab + det_ba_bp;
    if ( sum != 0.0 )
        return sum > 0.0;

    // exactly on the line: make sure that swap of a and b negates result.
    return ( a.x>b.x ) || ( a.x==b.x && a.y>b.y );
}
===============================================================================

Of course this is only one example of avoiding the "on the edge" case in
geometrical computing. It may not always be this easy to guarantee numerical
stability, especially not in 3D geometrical computing. However, it is my
opinion that treating "on the edge" as a special case in geometrical
algorithms will only make things worse.


Best regards,

Jaap Glas

-- 
-- dr. Jaap C. Glas
-- Software Engineer
-- dGB Earth Sciences
-- Nijverheidstraat 11-2
-- 7511 JM Enschede, The Netherlands
-- jaap.glas at dgbes.com
-- http://www.dgbes.com
-- Tel: +31 534315155, Fax: +31 534315104