Author(s) of this documentation:: Tristan Roussillon

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 alpha-shape 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.

Introduction

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:

in dimension two
- line ( \( k = 2 \))
- circle ( \( k = 3 \))
- circle of given radius and orientation ( \( k = 2 \))
- circle passing by a given point ( \( k = 2 \))
in dimension three
- plane ( \( k = 3 \))
- sphere ( \( k = 4 \))

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:

0 if \( p_{k+1} \) is lying on \( \mathcal{S} \).
strictly negative if \( p_{k+1} \) belongs to the interior of \( \mathcal{S} \).
strictly positive if \( p_{k+1} \) belongs to the exterior of \( \mathcal{S} \).

Note: Interior and exterior are not geometric concepts, but algebraic ones: they depend on the lexicographic order of the points. For example, for a circle passing by three points in the plane, the geometric interior coincides with the algebraic interior iff the three points are counter-clockwise oriented.

Geometric predicates, which return a boolean, usually adapt models of COrientationFunctor, which return an algebraic distance, so that either

strictly positive values,
non-negative values,
strictly negative values,
or non-positive values are accepted.

PredicateFromOrientationFunctor2 is an example of such adapter.

Orientation of three points in the plane

In order to determine whether 3 given points are aligned, clockwise oriented or counter-clockwise 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 half-plane lying on the left of the oriented line \( (PQ) \) (ie. \( P, Q, R \) are counter-clockwise oriented)
strictly negative if \( R \) belongs to the open half-plane lying on the right of the oriented line \( (PQ) \) (ie. \( P, Q, R \) are clockwise oriented)

Basic usage

Let us assume that we are working with the following domain:

  //domain definition
  typedef SpaceND< 2, DGtal::int16_t > DigitalSpace;
  typedef HyperRectDomain< DigitalSpace > Domain; 
  typedef Domain::Point Point; 
  Domain domain( Point(-32767,-32767), Point(32767,32767) ); 

To determine the orientation of the three following points...

  Point P, Q, R; 
  P = Point(0,0); 
  Q = Point(5,2); 
  R = Point(2,1); 
  //problem: are P, Q, R counter-clockwise oriented ?

...you must construct an orientation functor as follows:

  //orientation functor
  typedef InHalfPlaneBySimple3x3Matrix<Point, DGtal::int32_t> OrientationFunctor; 
  OrientationFunctor orientationFunctor;

Then, you can adapt this functor in order to get a predicate:

  //geometric predicate
  PredicateFromOrientationFunctor2<OrientationFunctor> 
    pointPredicate( orientationFunctor ); 

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:

  //initialization
  pointPredicate.init( P, Q ); 
  //which is equivalent to 
  //orientationFunctor.init( P, Q );
  //because the predicate stores a pointer to the functor

  //decision
  isCCW = pointPredicate( R ); 
  //which is equivalent to the following shortcut:
  //isCCW = pointPredicate( P, Q, R ); 

Note: This small example, which may be found in exampleInHalfPlane.cpp, requires the following headers:
#include "DGtal/geometry/tools/determinant/PredicateFromOrientationFunctor2.h"

#include "DGtal/geometry/tools/determinant/InHalfPlaneBySimple3x3Matrix.h"

List of available functors

Two useful classes have been implemented:

InHalfPlaneBySimple3x3Matrix, an orientation functor based on the determinant of a 3x3 matrix:

\( \begin{vmatrix} P_x & Q_x & R_x \\ P_y & Q_y & R_y \\ 1 & 1 & 1 \end{vmatrix} \)

InHalfPlaneBy2x2DetComputer, an orientation functor that reduces the computation to the determinant of a 2x2 matrix:

\( \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.
Filtered2x2DetComputer, which is a lazy adapter of any other determinant computer: the adaptee is only used for determinants close to zero.

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.

How to avoid overflows ?

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 30 (not zero) bits
  typedef PointVector<2, DGtal::int32_t> Point; 
  typedef InHalfPlaneBySimple3x3Matrix<Point, DGtal::int64_t> Functor; 

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 double-precision floating-point data type.

  //for coordinates of 52 (not zero) bits
  typedef PointVector<2, DGtal::int64_t> Point;  
  typedef AvnaimEtAl2x2DetSignComputer<double> DetComputer; 
  typedef Filtered2x2DetComputer<DetComputer> FDetComputer; 
  typedef InHalfPlaneBy2x2DetComputer<Point, FDetComputer> Functor; 

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 coordinates of 62 (not zero) bits
  typedef PointVector<2, DGtal::int64_t> Point; 
  typedef AvnaimEtAl2x2DetSignComputer<DGtal::int64_t> DetComputer; 
  typedef InHalfPlaneBy2x2DetComputer<Point, DetComputer> Functor; 

For greater coordinates, we recommand to use InHalfPlaneBy2x2DetComputer together with Simple2x2DetComputer using DGtal::BigInteger as integral types (WITH_GMP ON).

  //for arbitrary coordinates
  typedef PointVector<2, DGtal::BigInteger> Point; 
  typedef Simple2x2DetComputer<DGtal::BigInteger> DetComputer; 
  typedef InHalfPlaneBy2x2DetComputer<Point, DetComputer> Functor; 
  Functor *a= new Functor();
  BOOST_VERIFY( a);

Benchmark

Experimental tests justify the above recommendations. In testInHalfPlane-benchmark.cpp, we compare several methods:

3x3 : InHalfPlaneBySimple3x3Matrix
2x2 : InHalfPlaneBy2x2DetComputer with Simple2x2DetComputer
2x2-avnaim : InHalfPlaneBy2x2DetComputer with AvnaimEtAl2x2DetSignComputer
2x2-avnaim++ : a combination of InHalfPlaneBy2x2DetComputer, Filtered2x2DetComputer and AvnaimEtAl2x2DetSignComputer

Input and output types follow the name of the method. For instance, the 3x3-int32-int64-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. testInHalfPlane-benchmark.cpp was compiled in Release mode with gcc 4.6.3.

We performed 1 million orientation tests on five different kinds of inputs:

Three points \( P, Q, R \) of random coordinates.
Random \( P \) and \( Q \) and \( R = Q \) (null determinants).
Random \( P \) and \( Q \) and \( R = P \) (null determinants).
\( P = (0,0) \) and \( Q-P = k\vec{u}, R-P = l\vec{u} \) with random \( k, l, \vec{u}_x \vec{u}_y \) (null determinants).
\( P = (0,0) \) and \( Q-P = k\vec{u} + \epsilon, R-P = l\vec{u} + \epsilon \) with random \( k, l, \vec{u}_x \vec{u}_y \) and \( \epsilon \) (determinants close to zero).

For coordinates of \( 30 \) bits, we obtain the following results (running times are in seconds):

method vs input	1.	2.	3.	4.	5.
3x3-int32-int64	0.28	0.18	0.18	0.13	0.18
3x3-int32-BigInt	2.54	2.3	2.4	2.34	2.36
2x2-int32-int64	0.28	0.18	0.19	0.12	0.2
2x2-int32-BigInt	0.61	0.5	0.46	0.45	0.52
2x2-avnaim-int32-int32	0.32	0.2	0.18	0.3	0.36
2x2-avnaim-int32-double	0.33	0.21	0.19	0.38	0.46
2x2-avnaim++-int32-double	0.28	0.21	0.2	0.38	0.19

Best methods are 3x3-int32-int64 and 2x2-int32-int64. We therefore recommend to use 3x3-int32-int64, whose type is simpler to define.

Note

As expected 2x2-avnaim++-int32-double performs well, excepted for null determinants (columns 2, 3 and 4), because more steps are required in such cases to take a decision.
2x2-int32-BigInt outperforms 3x3-int32-BigInt 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.
3x3-double-BigInt	3.32	3.02	3.01	2.7	2.75
2x2-double-BigInt	1.08	0.91	0.72	0.76	0.84
2x2-avnaim-int64-int64	0.55	0.36	0.34	0.87	0.94
2x2-avnaim-double-int64	0.56	0.37	0.35	0.87	0.94
2x2-avnaim-int64-double	0.52	0.34	0.33	0.6	0.65
2x2-avnaim-double-double	0.5	0.33	0.31	0.58	0.68
2x2-avnaim++-int64-double	0.5	0.35	0.33	0.62	0.24
2x2-avnaim++-double-double	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. 2x2-avnaim++-double-double (last line).

For coordinates of \( 62 \) bits, we obtain the following results:

method vs input	1.	2.	3.	4.	5.
3x3-BigInt-BigInt	7.4	6.62	4.89	2.4	2.54
2x2-BigInt-BigInt	4.76	3.54	3.49	1.28	1.44
2x2-avnaim-int64-int64	0.65	0.42	0.4	0.98	1.06
2x2-avnaim++-int64-longdouble	0.6	0.42	0.4	0.78	0.26

Best methods are 2x2-avnaim-int64-int64 and 2x2-avnaim++-int64-longdouble, where long double is implemented as the 80-bit extended precision type. If this implementation is available on your system, you should use
2x2-avnaim++-int64-longdouble in order to perform faster computations in the case of (quasi)-collinear (but not confunded) points (columns 4 and 5).