Skip to content

Exact

Jump to bottom
Bruno edited this page Jun 17, 2024 · 12 revisions

Exact numbers and exact predicates

What is the problem ?

The geometric algorithms in Geogram, such as 2D triangulations, 3D triangulations, mesh intersection and mesh CSG depend on on more "elementary questions", such as whether a point is above or below a plane (in a certain sense).

Such "elementary questions", called predicates, are functions that take as an argument a (small) number of points (or simple geometric objects) and that returns a set of discrete values.

For instance, consider four points ${\bf p}_1,{\bf p}_2,{\bf p}_3,{\bf p}_4$. One may want to know the position of ${\bf p}_4$ relative to the supporting plane of ${\bf p}_1,{\bf p}_2,{\bf p}_3$, that can be one of ABOVE, BELOW, ON_PLANE.

These predicates are the "nevralgic" point of mesh intersection methods: if at one moment the algorithm "thinks" that ${\bf p}_4$ is above the supporting plane of ${\bf p}_1,{\bf p}_2,{\bf p}_3$, it is important that at another moment the predicate does not say that ${\bf p}_4$ is below the same plane. How could this happen ?

In fact, these predicates can be expressed as the sign of a polynomial in the coordinates of the points, and due to the limited precision of floating point numbers, the output of the predicate can be different from the exact mathematical result, especially around zero, and it can depend on the order of the points. Try this: write the following little C program:

#include <stdio.h>

void test_add3(double x, double y, double z) {
   printf("%f+%f+%f=%f\n",x,y,z,x+y+z);
}

int main() {
   double x1 = 1e30;
   double x2 = 1e-6;
   double x3 = -1e30;
   test_add3(x1,x2,x3);
   test_add3(x1,x3,x2);
   return 0;
}

in floating point arithmetic, imagine you compute $(x_1 + x_2) + x_3$, where $x_1 = 1e30$, $x_2 = 1e-6$ and $x_3 = -1e30$ (the result should be $1e-6$). When the computer first evaluates $x_1 + x_2$, it gets $1e30$ (because $x_2$ is too small relative to $x_1$), and in the end you get $0$. Now if you compute $(x_1 + x_3) + x_2$, you will get a different result ($1e-6$). Because of that, you may obtain a different result when you ask for the position of ${\bf p}_4$ relative to the supporting plane of ${\bf p}_1, {\bf p}_2, {\bf p}_3$ or when you ask for the position of ${\bf p}_4$ relative to ${\bf p}_3, {\bf p}_2, {\bf p}_1$ ! It can have catastrophic consequences, such as generating an incorrect mesh.

See also...

The optional geogramplus expansion package

To gain more speed and more robustness in the extreme cases, read about GeogramPlus.

Clone this wiki locally