DGtal
1.2.beta

Several algorithms of the Geometry package rely on geometric predicates. For instance, any convex hull algorithm relies on an orientation test (see Convex hull and alphashape in the plane for algorithms computing the convex hull of a finite planar set of points).
If the test is not exact or if it is not robust to possible overflows, the algorithm can return an incorrect result or can even infinitely loop. We provide below several solutions to cope with such problems.
In our framework, we consider points in a space of dimension \( n \). The orientation of \( k+1 \) points \( (p_1, \ldots, p_{k+1}) \) is given by the sign of the algebraic distance between \( p_{k+1} \) and an algebraic variety, chosen such that it is uniquely defined by the first \( k \) points \( (p_1, \ldots, p_{k}) \).
An algebraic variety is the set of points where a given \(n\)variate polynomial \( \mathcal{P} \) is evaluated to zero, ie. \( \{ x \in \mathbb{R}^n  \mathcal{P}(x) = 0 \} \). We consider only \( n\)variate polynomials that are sums of exactly \( k+1 \) monomials, so that they are uniquely defined by exactly \( k \) points. The algebraic distance of \( p_{k+1} \) to the algebraic variety of polynomial \( \mathcal{P} \) is
equal to \( \mathcal{P}(p_{k+1}) \).
Examples of algebraic varieties are:
The main concept of this module is COrientationFunctor, but there are also refinements, which are specific to small polynomials: COrientationFunctor2 and COrientationFunctor3.
Models of COrientationFunctor provide a method init() taking a static array of \( k \) points \( (p_1, \ldots, p_{k}) \) as argument. This initialization step is geometrically equivalent to defining the unique shape \( \mathcal{S} \) passing by these \( k \) points. Models of COrientationFunctor are also refinements of CPointFunctor. As such, they provide a parenthesis operator, which takes an extra point \( p_{k+1} \) as argument and returns a signed value that is:
PredicateFromOrientationFunctor2 is an example of such adapter.
In order to determine whether 3 given points are aligned, clockwise oriented or counterclockwise oriented, we provide several predicates. They are actually adaptations of some models of COrientationFunctor2, devoted to the computation of the distance of a point to a line.
Instances of such classes must be initialized from two points \( P, Q \) and provide a parenthesis operator taking a third point \( R \) as argument and returning a signed value.
The resulting value may be interpreted as follows:
 equal to 0 if \( P, Q, R \) belong to the same line
 strictly positive if \( R \) belongs to the open halfplane lying on the left of the oriented line \( (PQ) \) (ie. \( P, Q, R \) are counterclockwise oriented)
 strictly negative if \( R \) belongs to the open halfplane lying on the right of the oriented line \( (PQ) \) (ie. \( P, Q, R \) are clockwise oriented)
Let us assume that we are working with the following domain:
To determine the orientation of the three following points...
...you must construct an orientation functor as follows:
Then, you can adapt this functor in order to get a predicate:
The default behavior of PredicateFromOrientationFunctor2 is to return 'true' for strictly positive functor values.
The test can be done in one or two separate steps as follows:
Two useful classes have been implemented:
\( \begin{vmatrix} P_x & Q_x & R_x \\ P_y & Q_y & R_y \\ 1 & 1 & 1 \end{vmatrix} \)
\( \begin{vmatrix} Q_x  P_x & R_x  P_x \\ Q_y  P_y & R_y  P_y \end{vmatrix} \)
InHalfPlaneBy2x2DetComputer delegates the computation of this determinant to a model of C2x2DetComputer.
 Simple2x2DetComputer, which merely computes \( (Q_x  P_x)(R_y  P_y)  (Q_y  P_y)(R_x  P_x) \).
 AvnaimEtAl2x2DetSignComputer, an implementation of [Avnaim et.al., 1997 : [8]] that returns the sign of the determinant without increasing the size of the matrix entries.
Most classes are template classes parametrized by a type for the points (or its coordinates) and an integral type for the computations. All these implementations return an exact value (or sign), provided that the integral type used for the computations is well chosen with respect to the coordinates of the points.
Let \( x \) and \( x' \) be respectively \( b \)bits and \( b' \)bits integers. The sum \( x+x' \) may require \( \max(b,b')+1 \) bits and the product \( xx' \) may require \( b+b' \) bits. Consequently, we can determine the number of bits required for the different computations.
If the coordinates of the points are \( b \)bits integers, both InHalfPlaneBySimple3x3Matrix and InHalfPlaneBy2x2DetComputer with Simple2x2DetComputer may return determinant values of \( 2b + 3 \) bits. However, InHalfPlaneBy2x2DetComputer with AvnaimEtAl2x2DetSignComputer only require integers of \( b+1 \) bits.
For coordinates of \( 30 \) bits, lying within the range \( ]2^{30}; 2^{30}[ \), we recommand to use InHalfPlaneBySimple3x3Matrix with DGtal::int64_t as result type.
For coordinates of \( 52 \) bits, lying within the range \( ]2^{52}; 2^{52}[ \), we recommand to use InHalfPlaneBy2x2DetComputer with a lazy implementation of [Avnaim et.al., 1997 : [8]] using the \( 53 \) bits of the mantissa of the doubleprecision floatingpoint data type.
For coordinates of \( 62 \) bits, lying within the range \( ]2^{62}; 2^{62}[ \), we recommand to use InHalfPlaneBy2x2DetComputer with an implementation of [Avnaim et.al., 1997 : [8]] using DGtal::int64_t as a working type.
For greater coordinates, we recommand to use InHalfPlaneBy2x2DetComputer together with Simple2x2DetComputer using DGtal::BigInteger as integral types (WITH_GMP ON).
Experimental tests justify the above recommendations. In testInHalfPlanebenchmark.cpp, we compare several methods:
3x3
: InHalfPlaneBySimple3x3Matrix2x2
: InHalfPlaneBy2x2DetComputer with Simple2x2DetComputer2x2avnaim
: InHalfPlaneBy2x2DetComputer with AvnaimEtAl2x2DetSignComputer2x2avnaim++
: a combination of InHalfPlaneBy2x2DetComputer, Filtered2x2DetComputer and AvnaimEtAl2x2DetSignComputerInput and output types follow the name of the method. For instance, the 3x3int32int64
method runs on coordinates of type DGtal::int32_t and return a determinant of type DGtal::int64_t.
We ran all the above methods on a laptop (Intel core i5 2.50GHz, 4GB RAM) with Ubuntu 12.04. testInHalfPlanebenchmark.cpp was compiled in Release mode with gcc 4.6.3.
We performed 1 million orientation tests on five different kinds of inputs:
1. Three points \( P, Q, R \) of random coordinates.
method vs input  1.  2.  3.  4.  5. 

3x3int32int64  0.28  0.18  0.18  0.13  0.18 
3x3int32BigInt  2.54  2.3  2.4  2.34  2.36 
2x2int32int64  0.28  0.18  0.19  0.12  0.2 
2x2int32BigInt  0.61  0.5  0.46  0.45  0.52 
2x2avnaimint32int32  0.32  0.2  0.18  0.3  0.36 
2x2avnaimint32double  0.33  0.21  0.19  0.38  0.46 
2x2avnaim++int32double  0.28  0.21  0.2  0.38  0.19 
Best methods are 3x3int32int64
and 2x2int32int64
. We therefore recommend to use 3x3int32int64
, whose type is simpler to define.
2x2avnaim++int32double
performs well, excepted for null determinants (columns 2, 3 and 4), because more steps are required in such cases to take a decision. 2x2int32BigInt
outperforms 3x3int32BigInt
because the allocation/desallocation of BigIntegers, which represent a substantial part of running time, is minimized in Simple2x2DetComputer.For coordinates of \( 52 \) bits, we obtain the following results:
method vs input  1.  2.  3.  4.  5. 

3x3doubleBigInt  3.32  3.02  3.01  2.7  2.75 
2x2doubleBigInt  1.08  0.91  0.72  0.76  0.84 
2x2avnaimint64int64  0.55  0.36  0.34  0.87  0.94 
2x2avnaimdoubleint64  0.56  0.37  0.35  0.87  0.94 
2x2avnaimint64double  0.52  0.34  0.33  0.6  0.65 
2x2avnaimdoubledouble  0.5  0.33  0.31  0.58  0.68 
2x2avnaim++int64double  0.5  0.35  0.33  0.62  0.24 
2x2avnaim++doubledouble  0.46  0.34  0.32  0.67  0.26 
The best method is the lazy implementation of [Avnaim et.al., 1997 : [8]] using the \( 53 \) bits of the mantissa of the double data type, ie. 2x2avnaim++doubledouble
(last line).
For coordinates of \( 62 \) bits, we obtain the following results:
method vs input  1.  2.  3.  4.  5. 

3x3BigIntBigInt  7.4  6.62  4.89  2.4  2.54 
2x2BigIntBigInt  4.76  3.54  3.49  1.28  1.44 
2x2avnaimint64int64  0.65  0.42  0.4  0.98  1.06 
2x2avnaim++int64longdouble  0.6  0.42  0.4  0.78  0.26 
Best methods are 2x2avnaimint64int64
and 2x2avnaim++int64longdouble
, where long double is implemented as the 80bit extended precision type. If this implementation is available on your system, you should use
2x2avnaim++int64longdouble
in order to perform faster computations in the case of (quasi)collinear (but not confunded) points (columns 4 and 5).