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:

/* hello_darkness.c */
#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 = -1e30;
   double x3 = 1e-6;
   test_add3(x1,x2,x3);
   test_add3(x1,x3,x2);
   return 0;
}

The first test computes $x_1+x_2+x_3$ and the second one $x_1+x_3+x_2$, mathematically they shoud return the same result, but we got limited precision: in the first expression, everything goes well, we first compute $x_1$ + $x_2$, the result is exactly zero, then add $x_3$, and the final result is $1e-6$. In the second expression, it is another story, when ones computes $x_1+x_3$, that is, $1e30$ + $1e-6$, double precision numbers cannot store all the digits, then the result is truncated (and in this case, it means ignoring $1e-6$. Finally, when we subtract $x_3$, we get exactly zero, which is wrong !

Why should we care ?

If you think of what may happen in the middle of a geometric algorithm, since the result of an operation depends on the order of the operands, 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.

So what can we do ?

The idea is to implement a C++ class expansion_nt that represents numbers exactly. How is it possible in a computer ? First thing, we will implement a very limited set of operations:

  • addition: expansion_nt + expansion_nt $\rightarrow$ expansion_nt
  • subtraction: expansion_nt - expansion_nt $\rightarrow$ expansion_nt
  • product: expansion_nt * expansion_nt $\rightarrow$ expansion_nt
  • sign: expansion_nt $\rightarrow$ {NEGATIVE,ZERO,POSITIVE}

...and that's all folks ! This suffices to implement all the predicates that we need, because all these predicates boil down to computing the signs of polynomials of the point's coordinates.

But wait, we do not have division, how can it work ? You have mesh intersections in geogram right ? How can you do that without divisions ?

As always done in computer graphics and geometry processing, one trick is to use homogeneous coordinates: the point [x/w, y/w, z/w] is represented by a vector [x, y, z, w], each component represented by an expansion_nt. Then, each time one needs to divide, one multiplies $w$ instead. The other trick is to have a rational type, that represents a/b and that stores [a,b], with the four operations implemented. This rational type can be used to store the intermediary results for some operations.

The same example with expansion_nt

Let us now write the same example program, but this time with expansion_nt instead of double:

#include <stdio.h>
#include <geogram/numerics/expansion_nt.h>

void test_add3(double x, double y, double z) {
   GEO::expansion_nt X(x);
   GEO::expansion_nt Y(y);
   GEO::expansion_nt Z(z);
   printf("%f + %f + %f = %f, sign = %d\n",x,y,z,(X+Y+Z).estimate(), (X+Y+Z).sign());
}

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

As can be seen, for both expressions, expansion_nt computes the correct sign. In addition, expansion_nt has an estimate() member function, that gives an (inexact) estimation of the represented number (and in this particular case it gives the correct answer).

I would like to use expansion_nt but geogram is a bit heavy, what could I do ?

You can use the PSM (pluggable software module) geogram.psm.MultiPrecision, it is a .cpp/.h pair of files extracted from geogram. There is also an example program plus instructions to compile it.

The Predicate Constructing Kit (WIP tutorial)

Interval arithmetics (WIP tutorial)

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