Robust Predicates and Degeneracy

6.838 March 16 Meeting

Marek Teichmann

1. Robust Predicates

Predicates [Implementations]

See Shewchuk's Robust Predicates page.

IEEE Floating Point Arithmetic.

See: Notes on IEEE 754 Floating Point Arithmetic (postscript)

252 numbers between tick marks. Smallest number > 1 is 1 + 2-52.

"Island" of high precision around 0.

 fl(1+2-53) = 1

Geometric predicates can FAIL if implemented in FP.

Predicates implemented using floating point might give wrong answers:
get unexpected or inconsistent state of a program.

Why? Cannot distinguish the two cases below if a = 1:

[Example of a query]  

| orientation | = (a + e) - a

Example: a = 1, e = 2-53: | orientation | = 0 !

Not invariant w.r.t transformations (a = 0, 1).

Geometric Algorithms.


Geometric algorithm : primitive tests on continuous data --> combinatorial output [Ex. convex hull]

There is a consistency condition between numerical and combinatorial data.

Can geometric algorithms be implemented using floating point arithmetic? At all?


Solution: Exact Arithmetic.

Assumption: input is exact.

Numbers can be represented implicitly.

Evaluate all predicates exactly or up to the required precision (ex. sign = 1bit).

What kind of numbers [subgroups of R] do we need?

Given integer input, what is the smallest class of numbers sufficient for

Number representations.

See survey by Dube and Yap (postscript) for survey and references.

Implementation of exact arithmetic: compute every quantity exactly.

Lazy arithmetic: expression based, precision-driven arithmetic.

Numerical quantities are computed to sufficient precision so that the geometric predicates return correct values. (Underlying combinatorial structure is mathematically exact.)

Expression Arithmetic: some fundamental ideas.

Geometric predicates are expressions over input coordinates.

How to speed up exact arithmetic? Avoid computing all digits.

Must specify desired precision for the expressions: ex. 1 bit = sign.

First idea: evaluate expression in fp, then exactly if needed [Fortune and van Wyk].

Second idea: divide expression into subexpressions, evaluate incrementally.

Expression tree:

Rational Numbers, finitely representable.

Priest & Shewchuk:
non-overlapping, floating point terms.

Idea: allow roundoff error, but account for it after the fact by using more terms.

Example: (3 mantissa bit floating point numbers):

101.0011 = 101 + 11. 2-4
Example: 101 + 0.0011 = 101 + roundoff error terms

Operations: +, -, *, /

Highest order term is a crude approximation.
Evaluate expression using MST first, then add additional terms as needed.

Algebraic Numbers. Real/Expr, LEDA

Example: Need to compare an expression to 0.

Need precision bound for root (ex. 1 bit).

1. Propagate precision bound top down (from root to leafs).

Example: a+b: if precision bound for (a+b) is 0.02, then precision bound for a and b is 0.01.

At each leaf: procedure to extract an approximate value to satisfy precision bound at the leaf.

2. Propagate value up.

Changing precision bound at root may force (incremental) recomputation.

To how many bits do we need to evaluate an expression to get at least 1 significant bit?

Example: Graham's problem:

a - b = 0 ? (a=b to 40 digits!!)

Algebraic numbers are roots of polynomials (characteristic polynomial.)
How close is an expression to 0? --> Root bounds.

For each node in expression tree maintain

Expression tree keeps track of required precision, which is propagated down the tree to know how precisely a given subexpression (or number) has to be computed.

Disadvantages: bit size grows, speed about 4 to 10 times slower than fp.

Packages using exact arithmetic:

Using Finite Precision: Making finite precision predicates robust.


Definition [Fortune]:
Robustness: A geometric algorithm (using f.p.) is robust if

Stability: An algorithm is stable if it is robust and the required perturbation is small.

Applications: line arrangements, 2D Delaunay triangulation (get "almost" delaunay property), convex hull.
Difficulty: careful error analysis of predicates, careful design of algorithm.

Inference of complex predicates from simpler ones.

Predicates : independent and dependent.

Interval Arithmetic.

[Moore, 1960]

Numbers are represented by an interval tex2html_wrap_inline670

Interval arithmetic allows computations to be performed on intervals, such that the result of the computation provides a guaranteed enclosure for the correct mathematical result over the domain defined by the the interval operands.

Basic arithmetic operations tex2html_wrap_inline674 , defined for the the intervals X and Y as all possible results


unless tex2html_wrap_inline680 and tex2html_wrap_inline682 .

Example:    [a,b] + [c,d] = [a+c, b+d]


Implementation rounds each interval endpoint outwards.

Predicates can return YES, NO, MAYBE (it is not possible to tell within precision available.)

Properties: (inclusion isotomy)


Interval geometry

     (thick lines) [qhull]
Fatten geometric objects, use tolerancing zones (conservative estimates.)
Predicates can return "don't know".

Topologicaly consistent distortion

     (polylines replace lines)

Example: Line arrangement construction.
Every intersection between two lines is "snapped" to a grid point.

Rounded geometry

(line coeffs rounded).
Goal: preserve combinatorial structure.
Fact: rounding a collection of disjoint simple polygons while preserving combinatorial structure is NP-complete.



Epsilon balls

Tolerances / epsilons useful in practice, but no proofs of correctness.
Consider an expression = 0 if |expression| < epsilon.
Tweak until your program works :)

Implicit representations

CSG, Grassman Algebra [Stolfi]
Example: represent a line by 2 points, a point by two lines,...
Ultimately comes down to linear algebra.

2. Degeneracy

Partity Algorithm: is a point inside or outside a polygon?

Voronoi Diagram: co-circular points.

Convex Hull: co-planar points. Is a point inside? on a face? on an edge? on a point?

Def.: A degeneracy occurs when arbitrarily small perturbation of the input can result in qualitatively different geometric structures. It is in some sense a discontinuity.

Inherent and algorithm-induced degeneracies.
Ex.: In divide and conquer algorithm for convex hull, colinear points on convex hull are an inherent degeneracy, collinear points inside are an induced degeneracy.

What is general position? (Depends on predicates: no predicates return 0).

Perturbation Schemes.

Algorithm working for generic input --> algorithm working for general input.

General idea:

Let a* be a non-degenerate instance of a problem.
Given a possibly degenerate input instance a, the instance (a + e a*) is non-degenerate for infinitesimal e.

Trivial perturbation: (ex. sorting)

(xj- xi) ? 0. Want >0 or <0

if xi = xj, return "xi < xj" iff i<j

Let e = very small number, with (xi+ei) != xj for any i,j.

( (xj+ ej) - (xi+ei) ) < 0 if xi = xj and i<j,

General comparisons: xij are variables (i=point id, j=coordinate).

Predicates = polynomials in x's and e's.


    linear perturbation    

random perturbation

work with `high' probability

Never compute the values.
Compare x's, then successively higher degrees of e.



Test query point if on edge
Convex hull Merge adjacent faces if coplanar

Remove faces subset of larger faces

Voronoi diagram Merge Voronoi vertices generated by co-circular points
Remove 0 length edges