DGtal  1.5.beta
Plane recognition and plane width computation
Author(s) of this documentation:
Jacques-Olivier Lachaud

Part of the Geometry package.

This part of the manual describes classes and functions related to the planarity of a set of points in 3D. For instance, it provides solutions to the following problems:

  • planarity decision: checks whether or not a given set of (digital or not) points belongs to a parallel strip of specified width. Strips may be specified with an axis width (like naive planes) or with a diagonal width (like standard planes).
  • plane recognition: decide if a given set of (digital or not) points is planar and provides an explicit solution (i.e. gives the normal and shift of some parallel strip containing all the points).
  • incremental or additive plane recognition: same as above, except that planarity of points can be checked progressively by adding points one at a time (incremental) or by groups (additive).
  • width computation: computes the minimal width of a parallel strip that containes a given set of (digital or not) points.

This module thus extends to the 3D case the recognition of digital straight segments and of blurred segments.

For now, these problems are solved thanks to two algorithms: the COBA algorithm by Charrier and Buzer [22] and the Chord algorithm by Gerard, Debled-Rennesson and Zimmermann [61] .

More precisely, these modules tackle naive digital planes (with or without a specified axis) and digital planes of arbitrary rational axis width.

Note
A current limitation is that these algorithms do not compute the minimal characteristics of the digital plane containing the input points, but only a valid parallel strip of approximated normal.

Related examples are viewer3D-7-planes.cpp, viewer3D-7bis-planes.cpp, viewer3D-7-stdplane.cpp, greedy-plane-segmentation.cpp, greedy-plane-segmentation-ex2.cpp.

Planarity as a parallel strip

Definition of parallel strip

A parallel strip is the intersection of two parallel half-planes such that each half-plane includes the other. Half-planes may be closed or open. Hence, parallel strips take one of the following form:

Here are a few useful definitions:

  • The normal to the strip is the vector N, which is a unit vector in the class ParallelStrip.
  • The width of the strip is the scalar \(\epsilon\).
  • The axis width of the strip is the quantity \(|\epsilon / N_i|\) where \(|N_i|\) is the greatest component. The main axis is the i-th axis.
  • The diagonal width is the quantity \( \frac{\epsilon \sqrt{n}}{\|N\|_1} \), where n is the dimension of the space. The main diagonal is the vector \(( \pm 1, \ldots, \pm 1) \), where the \( \pm 1 \) are the signs of the components of N.

Digital planes are specific cases of parallel strips. We give below their 3D definitions, but their n-dimensional definition is straightforward.

Naive digital planes

A naive digital plane is a set of digital points \((x,y,z) \in Z^3\) such that \( d \le ax+by+cz < d + \omega \), with a, b, c integer numbers and \( \omega = \max(|a|,|b|,|c|) \).

It is easily seen that naive planes are parallel strips of axis width strictly smaller than 1, by a simple division of both inequalities by \(\omega\).

Standard digital planes

A standard digital plane is a set of digital points \((x,y,z) \in Z^3\) such that \( d \le ax+by+cz < d + \omega \), with a, b, c integer numbers and \( \omega = |a|+|b|+|c| \).

It is easily seen that standard planes are parallel strips of diagonal width strictly smaller than \( \sqrt{3} \), by a simple division of both inequalities by \(\sqrt{a^2+b^2+c^2}\).

Naive digital plane recognition (and naive planes with rational width)

Since plane recognition solves also planarity decision, we focus on plane recognition. Note that both COBA algorithm and Chord algorithm takes the same time for a planarity decision and for a plane recognition. In this section, the user specifies the (strict) maximal axis width as a rational p/q. Naive digital planes are thus specified by 1/1.

Naive plane recognition (known axis) with COBA algorithm

The user should instantiate a COBANaivePlaneComputer with the appropriate digital space. The axis and the maximal axis width (as a rational number p/q) are specified by method COBANaivePlaneComputer::init. The diameter of the set of points S must also be specified (i.e. \( \max_{p,q \in S} \|p-q\|_\infty \)). Then, method COBANaivePlaneComputer::extend is called with the appropriate range of points. The method COBANaivePlaneComputer::primitive returns the solution (if any) as parallel strip of axis width smaller than p/q.

#include "DGtal/geometry/surfaces/COBANaivePlaneComputer.h"
...
std::vector<Z3i::Point> pts; // any container of 3D digital points
... // fill pts with your points
COBANaivePlaneComputer<Z3i::Space, int64_t> plane; // int64_t is enough for diameter 100
plane.init( 2, 100, 1, 1 ); // axis is z (2), diameter is 100, epsilon=1/1
bool isPlane = plane.extend( pts.begin(), pts.end() );
if ( isPlane )
std::cout << "This is a naive plane, well approximated by " << plane.primitive() << std::endl;
else
std::cout << "This is a not naive plane." << std::endl;
Space::Point Point
Definition: StdDefs.h:168
Note
The second template parameter for class COBANaivePlaneComputer gives the internal integer type for computation. Detailed explanations can be found in Detailed explanation of COBA plane recognition algorithm. The type should be able to hold integers of order \((2D^3)^2\) if D is the diameter of the set of digital points. In practice, diameter is limited to 20 for int32_t, diameter is approximately 500 for int64_t, and whatever with BigInteger / GMP integers. For huge diameters, the slow-down is polylogarithmic with respect to the diameter.

Naive plane recognition (known axis) with Chord algorithm

The user should instantiate a ChordNaivePlaneComputer with the appropriate digital space. The axis and the maximal axis width (as a rational number p/q) are specified by method ChordNaivePlaneComputer::init. Then, method ChordNaivePlaneComputer::extend is called with the appropriate range of points. The method ChordNaivePlaneComputer::primitive returns the solution (if any) as parallel strip of axis width smaller than p/q.

#include "DGtal/geometry/surfaces/ChordNaivePlaneComputer.h"
...
std::vector<Z3i::Point> pts; // any container of 3D digital points
... // fill pts with your points
ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, int64_t> plane; // int64_t is enough for diameter 440000000 !
plane.init( 2, 1, 1 ); // axis is z (2), epsilon=1/1
bool isPlane = plane.extend( pts.begin(), pts.end() );
if ( isPlane )
std::cout << "This is a naive plane, well approximated by " << plane.primitive() << std::endl;
else
std::cout << "This is a not naive plane." << std::endl;
Note
The second template parameter of the class ChordNaivePlaneComputer gives the type of input points. When recognizing digital planes, you shoud set it to Space::Point, where Space is your digital space. Note that the Chord algorithm works also for real-valued points: in this case, you can specify another kind of Point like RealPoint.
The third template parameter of the class ChordNaivePlaneComputer gives the internal scalar (or integer) type for computation. If the input points are digital points, then this type should be integral and hold integers of order \((2D)^2\) if D is the diameter of the set of digital points. In practice, diameter is limited to 14000 for int32_t, diameter is approximately 440000000 for int64_t, and whatever with DGtal::BigInteger / GMP integers. For huge diameters, the slow-down is polylogarithmic with respect to the diameter. If the input points are not digital, you should provide a scalar type for intermediate computations (e.g. double).

Naive plane recognition (unknown axis) with COBA algorithm

The user should instantiate a COBAGenericNaivePlaneComputer with the appropriate digital space. The maximal axis width (as a rational number p/q) is specified by method COBAGenericNaivePlaneComputer::init. The diameter of the set of points S must also be specified (i.e. \( \max_{p,q \in S} \|p-q\|_\infty \)). Then, method COBAGenericNaivePlaneComputer::extend is called with the appropriate range of points. The method COBAGenericNaivePlaneComputer::primitive returns the solution (if any) as parallel strip of axis width smaller than p/q.

#include "DGtal/geometry/surfaces/COBAGenericNaivePlaneComputer.h"
...
std::vector<Z3i::Point> pts; // any container of 3D digital points
... // fill pts with your points
COBAGenericNaivePlaneComputer<Z3i::Space, int64_t> plane; // int64_t is enough for diameter 100
plane.init( 100, 1, 1 ); // any axis, diameter is 100, epsilon=1/1
bool isPlane = plane.extend( pts.begin(), pts.end() );
if ( isPlane )
std::cout << "This is a naive plane, well approximated by " << plane.primitive() << std::endl;
else
std::cout << "This is a not naive plane." << std::endl;

The advantage is that the object detects progressively what is the correct main axis. You may know what is a correct main axis by calling COBAGenericNaivePlaneComputer::active().

Note
The principle of COBAGenericNaivePlaneComputer is to have three instances of COBANaivePlaneComputer at the beginning, one per possible axis. When extending the object by adding points, all active instances are extended. If any one of the active instance fails, it is removed from the active instances. It is guaranteed that if there is a naive plane (of specified width) that contains the given set of points, then the object COBAGenericNaivePlaneComputer will have at least one active COBANaivePlaneComputer instance. COBAGenericNaivePlaneComputer is thus a correct recognizer of arbitrary pieces of naive planes.
Template parameters have the same role as for class COBANaivePlaneComputer (see above).

Naive plane recognition (unknown axis) with Chord algorithm

The user should instantiate a ChordGenericNaivePlaneComputer with the appropriate digital space. The maximal axis width (as a rational number p/q) is specified by method ChordGenericNaivePlaneComputer::init. Then, method ChordGenericNaivePlaneComputer::extend is called with the appropriate range of points. The method ChordGenericNaivePlaneComputer::primitive returns the solution (if any) as parallel strip of axis width smaller than p/q.

#include "DGtal/geometry/surfaces/ChordGenericNaivePlaneComputer.h"
...
std::vector<Z3i::Point> pts; // any container of 3D digital points
... // fill pts with your points
ChordGenericNaivePlaneComputer<Z3i::Space, Z3i::Point, int64_t> plane; // int64_t is enough for diameter 440000000 !
plane.init( 1, 1 ); // epsilon=1/1
bool isPlane = plane.extend( pts.begin(), pts.end() );
if ( isPlane )
std::cout << "This is a naive plane, well approximated by " << plane.primitive() << std::endl;
else
std::cout << "This is a not naive plane." << std::endl;

The advantage is that the object detects progressively what is the correct main axis. You may know what is a correct main axis by calling ChordGenericNaivePlaneComputer::active().

Note
Similarly to COBAGenericNaivePlaneComputer, the principle of ChordGenericNaivePlaneComputer is to have three instances of ChordNaivePlaneComputer at the beginning, one per possible axis.

Standard plane recognition with COBA algorithm

The user should instantiate a COBAGenericStandardPlaneComputer with the appropriate digital space. The maximal diagonal width (as a rational number \(p/q \times \sqrt{3}\) ) is specified by method COBAGenericStandardPlaneComputer::init. The diameter of the set of points S must also be specified (i.e. \( \max_{p,q \in S} \|p-q\|_\infty \)). Then, method COBAGenericStandardPlaneComputer::extend is called with the appropriate range of points. The method COBAGenericStandardPlaneComputer::primitive returns the solution (if any) as parallel strip of axis width smaller than p/q.

#include "DGtal/geometry/surfaces/COBAGenericStandardPlaneComputer.h"
...
std::vector<Z3i::Point> pts; // any container of 3D digital points
... // fill pts with your points
COBAGenericStandardPlaneComputer<Z3i::Space, int64_t> plane; // int64_t is enough for diameter 100
plane.init( 100, 1, 1 ); // any axis, diameter is 100, epsilon=1/1 x sqrt(3)
bool isPlane = plane.extend( pts.begin(), pts.end() );
if ( isPlane )
std::cout << "This is a standard plane, well approximated by " << plane.primitive() << std::endl;
else
std::cout << "This is a not standard plane." << std::endl;

The advantage is that the object detects progressively what is the correct main diagonal. You may know what is a correct main diagonal by calling COBAGenericStandardPlaneComputer::active().

Note
The principle of COBAGenericNaivePlaneComputer is to have four instances of COBANaivePlaneComputer at the beginning, all of z axis. When extending the object by adding points, the points are transformed such as: \((x,y,z) \mapsto (x \pm z, y \pm z, z)\). We need four instances for the four possible sign combinations. It is guaranteed that if there is a naive plane (of specified width) that contains the given set of points, then it corresponds to a standard plane before the transformation. Then the object COBAGenericStandardPlaneComputer will have at least one active COBANaivePlaneComputer instance. COBAGenericStandardPlaneComputer is thus a correct recognizer of arbitrary pieces of standard planes.
Template parameters have the same role as for class COBANaivePlaneComputer (see above).

Standard plane recognition with Chord algorithm

The user should instantiate a ChordGenericStandardPlaneComputer with the appropriate digital space. The maximal diagonal width (as a rational number \(p/q \times \sqrt{3}\) ) is specified by method ChordGenericStandardPlaneComputer::init. Then, method ChordGenericStandardPlaneComputer::extend is called with the appropriate range of points. The method ChordGenericStandardPlaneComputer::primitive returns the solution (if any) as parallel strip of axis width smaller than p/q.

#include "DGtal/geometry/surfaces/ChordGenericStandardPlaneComputer.h"
...
std::vector<Z3i::Point> pts; // any container of 3D digital points
... // fill pts with your points
ChordGenericStandardPlaneComputer<Z3i::Space, Z3i::Point, int64_t> plane;
plane.init( 1, 1 ); // any axis, epsilon=1/1 x sqrt(3)
bool isPlane = plane.extend( pts.begin(), pts.end() );
if ( isPlane )
std::cout << "This is a standard plane, well approximated by " << plane.primitive() << std::endl;
else
std::cout << "This is a not standard plane." << std::endl;

The advantage is that the object detects progressively what is the correct main diagonal. You may know what is a correct main diagonal by calling ChordGenericStandardPlaneComputer::active().

Note
The principle of ChordGenericNaivePlaneComputer is to have four instances of ChordNaivePlaneComputer at the beginning, all of z axis. When extending the object by adding points, the points are transformed such as: \((x,y,z) \mapsto (x \pm z, y \pm z, z)\). We need four instances for the four possible sign combinations. It is guaranteed that if there is a naive plane (of specified width) that contains the given set of points, then it corresponds to a standard plane before the transformation. Then the object ChordGenericStandardPlaneComputer will have at least one active ChordNaivePlaneComputer instance. ChordGenericStandardPlaneComputer is thus a correct recognizer of arbitrary pieces of standard planes.
Template parameters have the same role as for class ChordNaivePlaneComputer (see above).

Incremental or additive plane recognition

Incremental and additive plane recognition are done in exactly the same way for all six classes COBANaivePlaneComputer, COBAGenericNaivePlaneComputer, COBAGenericStandardPlaneComputer, ChordNaivePlaneComputer, ChordGenericNaivePlaneComputer, ChordGenericStandardPlaneComputer. They are illustrated in examples viewer3D-7-planes.cpp, viewer3D-7bis-planes.cpp, viewer3D-7-stdplane.cpp, and viewer3D-7bis-stdplane.cpp.

Incremental plane recognition

These six classes are models of concepts::CIncrementalPrimitiveComputer. They all have an internal state that store a current parallel strip that contains all points previously inserted. Furthermore, they provide methods extend( p ) and isExtendable( p ) to insert a new point p and to check whether there is a some parallel strip that contains both the previously inserted points and the new point p.

  • extend( Point p ) : bool. Tries to add point p to the current plane computer. If it is possible, update potentially the parallel strip and returns true, otherwise returns false and the object is left unchanged.
  • isExtendable( Point p ) : bool. Tries to add point p to the current plane computer. Returns true iff it is possible. The object is always left unchanged.
typedef ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, int32_t> PlaneComputer;
typedef PlaneComputer::Primitive Primitive;
PlaneComputer plane;
plane.init( 2, 1, 1 ); // Z-axis
Point pt0( 0, 0, 0 );
bool pt0_inside = plane.extend( pt0 ); // true
Point pt1( 8, 1, 3 );
bool pt1_inside = plane.extend( pt1 ); // true
Point pt2( 2, 7, 1 );
bool pt2_inside = plane.extend( pt2 ); // true
Point pt3( 0, 5, 12 );
bool pt3_inside = plane.extend( pt3 ); // false
Point pt4( -5, -5, 10 );
bool pt4_inside = plane.extend( pt4 ); // false
Point pt5 = pt0 + pt1 + pt2 + Point( 0, 0, 1 );
bool pt5_inside = plane.extend( pt5 ); // true
Point pt6 = Point( 1, 0, 1 );
bool pt6_inside = plane.extend( pt6 ); // true
Primitive strip = plane.primitive();
trace.info() << "strip=" << strip
<< " axis=" << strip.mainAxis()
<< " axiswidth=" << strip.axisWidth()
<< " diag=" << strip.mainDiagonal()
<< " diagwidth=" << strip.diagonalWidth()
<< std::endl;
std::ostream & info()
Trace trace
Definition: Common.h:153
MyPointD Point
Definition: testClone2.cpp:383
Note
Calling isExtendable before extend does not speed up the execution of extend. Hence the user should always prefer to call extend directly whenever possible.

Additive plane recognition

These six classes are models of concepts::CAdditivePrimitiveComputer. They all have an internal state that store a current parallel strip that contains all points previously inserted. Furthermore, they provide templated methods extend( it, itE ) and isExtendable( it, itE ) to insert a range [it, itE) of new points and to check whether there is a some parallel strip that contains both the previously inserted points and all the new points.

  • extend( PointIterator it, PointIterator itE ) : bool. Tries to add the range [it, itE) of points to the current plane computer. If it is possible, update potentially the parallel strip and returns true, otherwise returns false and the object is left unchanged.
  • isExtendable( PointIterator it, PointIterator itE ) : bool. Tries to add the range [it, itE) of points to the current plane computer. Returns true iff it is possible. The object is always left unchanged.
Note
Calling isExtendable before extend does not speed up the execution of extend. Hence the user should always prefer to call extend directly whenever possible.

Width of a set of points

Computing the axis width of a set of points is easily done with the static method ChordNaivePlaneComputer::computeAxisWidth. The axis width is returned as a pair (p,q), such that the width is p/q.

#include "DGtal/geometry/surfaces/ChordNaivePlaneComputer.h"
...
std::set<Z3i::Point> pts; // any container of 3D digital points
... // fill pts with your points
typedef ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, int64_t> PlaneComputer;
std::pair<Z3i::Integer,Z3i::Integer> width =
PlaneComputer::computeAxisWidth( 2, pts.begin(), pts.end() ); // z-axis

Comparative evaluation of COBA and Chord algorithm

For now, there are only two computers for computing a ParallelStrip primitive: the ChordNaivePlaneComputer (CHORD) and the COBANaivePlaneComputer (COBA).

  1. Complexity: (COBA) has a better worst time complexity than (CHORD), but neither (CHORD) nor (COBA) has an easy bound on the number of global recomputation (traversal of all input points to recompute a valid direction). From experiments, asymptotic behavior seems quasi-linear.
  2. Big integers: (CHORD) requires (significantly) smaller integers than (COBA). For instance, int64_t are required for diameter greater than 25 for (COBA) instead of 14000 for (CHORD).
  3. Practical speed: both algorithms are very comparable. (CHORD) seems slightly faster than (COBA) on average (but this was not tested on many architecture). According to the graph below, obtained with programs testCOBANaivePlaneComputer-benchmark.cpp and testChordNaivePlaneComputer-benchmark.cpp, (CHORD) is faster for a number of points below 200 while (COBA) is best for a big number of points. However, if the diameter is not too big, ChordNaivePlaneComputer with int64_t seems always faster.
  4. Exactness: Both algorithms do not return the smallest possible arithmetic parameters for the plane, but only a rational approximation.
  5. Services: Algorithm (CHORD) can be transformed to find the exact axis width of a given set of points, (COBA) is not suited for that task.
Evaluation of computation times of COBA and Chord algorithms according to the number of points. Times are in ms. We did not put the graph of COBANaivePlaneComputer with int64_t since we cannot exceed a diameter of 500. In this benchmark, the diameter was always 10 times the number of points.