W Randolph Franklin home page
... (old version)
AdvancedComputerGraphicsSpring2008/ home page Login


Theoretical Real numbers form a field and satisfy 11 axioms, such as:

  a+0=0+a
a+b = b+a
a(b+c) = ab+ac

Most of these are usually false for computer floating point because of roundoff. This causes many algorithms that are theoretically fine sometimes to fail in practice. The bigger the problem, the greater the chance of failure.

E.g., the point where two lines intersect is probably not exactly on either line. A polyhedron union algorithm in a CAD package that implicitly assumes this might crash or give the wrong answer.

I've personally observed this when writing a routine to test if the line segment AB intersected the line segment CD. Because of roundoff, my program sometimes reported that both points A and B were above the line CD and simultaneously both points C and D were above the line AB. That is not ever possible for any four points.

That's what motivates the work in robustness. Many people have spent much effort, and it's still not completely solved.

One solution for straight lines and straight-edges polygons is to compute with a new class of rational numbers represented as an integer numerator divided by an integer denominator:

  class {
    public:
    int numerator;
    int denominator;
  }

E.g. {$\frac{2}{5}+\frac{1}{4}=\frac{13}{20}$}.

You have to write your own routines for all the usual operations, such as long division.

Rationals are slower and the number of digits doubles whenever two rationals are added. Therefore the numerator and denominator have to be stored as arrays of digits. Any operations get very slow.

  class {
    public:
    vector<int> numerator;
    vector<int> denominator;
  }

However, when you want to intersect circles and lines exactly, you need to represent square roots exactly. In abstract algebra terms, you want to compute in the field of the rationals extended with a square root, {$(Q,\sqrt{\cdot})$}. Operationally, a number is represented as a polynomial that has the number as one of its roots, plus some way of identifying that root. All the arithmetic operations turn into operations on polynomials. Even telling if a number is zero is hard. Is the following number positive, negative, or zero?

{$ \sqrt{3+2\sqrt{2}} - \sqrt{3-2\sqrt{2}} - 2\sqrt{2}$}

One solution is to combine this with interval arithmetic. Here you represent each number as a range that contains that number. Instead of 3, use [2.9,3.1], or [2.9999,3.0001] or whatever. You can add etc:

[2.9,3.1] + [4.9,5.1] = [7.8,8.2]

The ranges get wider quickly.

Often you need to know only the sign of the result of a computation. This might determines whether two polygons intersect. If the interval doesn't contain zero, you're good. That strategy is this:

  1. Compute a number with interval arithmetic.
  2. If the interval doesn't contain zero, you have the answer.
  3. Otherwise, use a more powerful and slower technique, such as rationals, polynomials or whatever, perhaps evaluated to some precision.
  4. Keep increasing the precision until the interval doesn't contain zero.
  5. If the interval is small enough and still contains zero, then eventually you can use a theorem that says that the number is identically zero.
    It's called a zero separation theorem. Given a polynomial of degree n whose integer coefficients have d digits, then any root that isn't zero, has at least some minimum absolute value.

Computing in {$(Q,\sqrt{\cdot})$} still doesn't handle parametric splines. They need to be at least cubic to match the radius of curvature at the joints, i.e.

{$$ x = \sum_{i=0}^{i=3} \sum_{j=0}^{j=3} a_{ij} u^i v^j, \ \ y = \sum_{i=0}^{i=3} \sum_{j=0}^{j=3} b_{ij} u^i v^j, \ \ z = \sum_{i=0}^{i=3} \sum_{j=0}^{j=3} c_{ij} u^i v^j $$}

Mathmatically, this surface is a piece of an implicit function f(x,y,z)=0. There's a nice process using resultants for eliminating variables from sets of polynomial equations just as you can do for linear equations. However, unlike for linear equations, the degree of the resulting polynomial is something like the product of the degrees of the input polynomials. (I may have omitted a 2 or whatever.) IIRC f has degree 18.

Two patches intersect at a curve of degree 324 or so. Three patch intersect at a point that's a root of a polynomial of degree of about 5000. In a practical sense, working with polynomials like this is completely impossible. E.g., very small changes in the coefficients cause very large changes in the roots because the terms almost cancel out. So, approximations are still required.

Many years ago I wrote a paper Problems with raster graphics algorithms on the difference between ints and floats in computer graphics. Counterintuitively, ints were messier.