Generic functions for affine geometry

A set of generic functions is provided in header DGtal/geometry/tools/AffineGeometry.h for several affine geometry computations. They can be conveniently used since they deduce automatically the type of their parameters. All these functions support number types for points and vectors that are either integers (int32_t, int64_t, BigInteger) or floating-point numbers (float, double). The computations are exact when types are integral (up to integer overflow), while a tolerance can be provided for floating-point numbers (default tolerance: between -1e-12 and 1e-12 is zero).

Affine dimension and affine subset

functions::computeAffineDimension( const std::vector< TPoint >& X, const double tolerance ) Given a range of points X, returns the affine dimension of its spanned affine subspace.
functions::computeAffineSubset( const std::vector< TPoint >& X, const double tolerance ) Given a range of points X, returns a subset of these points that form an affine basis of X. Equivalently it is a simplex whose affine space spans all the points of X.
functions::computeAffineSubset( const std::vector<TPoint>& X, const TIndexRange& I, const double tolerance ) Given a range of points X and a subset of it given by indices I, returns a subset of these points that form an affine basis of X[I]. Equivalently it is a simplex whose affine space spans all the points of X[I].

Affine basis and independent vectors

The computed bases are in row echelon form (see Wikipedia ), potentially reduced if components are real-value and nnot integers, except when using completeBasis where the last vector may be asked to be orthogonal to all the others.

functions::getAffineBasis( TInputPoint& o, std::vector< TPoint >& basis, const std::vector< TInputPoint >& X, const double tolerance ) Given a range of points X, outputs a point and a range of vectors forming an affine basis containing X.
functions::getAffineBasis( TInputPoint& o, std::vector< TPoint >& basis, const std::vector< TInputPoint >& X, const TIndexRange& I, const double tolerance ) Given a range of points X and the indices I of points in X which form an affine subset of X, returns a point and a range of vectors forming an affine basis containing X.
functions::computeIndependentVector( const std::vector< TPoint >& basis, const double tolerance ) Given a partial basis of vectors, returns a returns a canonic unit vector that is independent.
functions::getCompleteBasis( std::vector< TPoint >& basis, bool normal_vector, const double tolerance ) Complete the vectors basis with independent vectors so as to form a basis of the space. When normal_vector is true, the last added vector is guaranteed to be orthogonal to all the previous vectors, otherwise it is just an independent canonic vector.

Affine basis of codimension 1

functions::getOrthogonalLatticeBasis( std::vector< TPoint >& B, const TInputPoint& N, bool shortened ) Given a lattice vector N, returns a possible basis for its orthogonal d-1 dimensional lattice. The returned basis is guaranteed to span all the lattice points contained in the vector space orthogonal to N. If shortened is false, the basis is in row echelon form, otherwise the vectors are shortened by pairwise linear combinations.

Orthogonal vector to basis of codimension 1

functions::computeOrthogonalVectorToBasis( const std::vector<TPoint>& basis ) Given d-1 independent vectors basis in dD, returns a vector that is orthogonal to each of them.
functions::getOrthogonalVector( TPoint& w, const std::vector< TInputPoint >& X, const TIndexRange& I, const double tolerance ) Given a range of points X, a range of indices I specifying the affine subset of interest, returns a vector that is orthogonal to this affine subset, if it is d-1-dimensional.
functions::getOrthogonalVector( TPoint& w, const std::vector< TInputPoint >& X, const double tolerance ) Given a range of points X, returns a vector that is orthogonal to this affine set, if it is d-1-dimensional.

Utilities

functions::computeSimplifiedVector( const TPoint& v ) Given a vector, returns the aligned vector with its component simplified by the gcd of all components (when the components are integers) or the aligned vector with a maximum oo-norm of 1 (when the components are floating-point numbers).

This example (similar to exampleAffineGeometry.cpp) shows how to process the affine geometry of a set of points that is not full dimensional.

#include "DGtal/base/Common.h"
#include "DGtal/kernel/PointVector.h"
#include "DGtal/geometry/tools/AffineGeometry.h"
 
template <typename T> std::ostream&
operator<<( std::ostream& out, const std::vector< T >& object )
{ // to display vector<T>
  out << "[";
  for ( auto t : object ) out << " " << t;
  out << " ]";
  return out;
}
 
using namespace DGtal;
int main()
{
  typedef PointVector< 3, int64_t> Point;
  std::vector< Point > X = {
    Point{21, -29, -21}, Point{-6, 14, 8}, Point{-19, 1, 9},
    Point{39, 16, -12}, Point{-26, 29, 23}, Point{-4, 45, 19},
    Point{6, 31, 9}, Point{9, 6, -2}, Point{8, 49, 15}
  };
 
  std::cout << "Dimension is " << Point::dimension << "\n";
  std::cout << "X = " << X << "\n";
  std::cout << "Expected dimension of affine set of points X is "
            << (Point::dimension-1) << "\n";
  std::cout << "AffineDim(X) = " << functions::computeAffineDimension( X ) << "\n";
  auto I = functions::computeAffineSubset( X );
  std::cout << "AffineSubset(X) = " << I << "\n";
  Point o;
  std::vector<Point> B;
  functions::getAffineBasis( o, B, X, I );
  std::cout << "AffineBasis(X) =: p+B = " << o << " + " << B << "\n";
  auto e = functions::computeIndependentVector( B );
  std::cout << "Independent(X) =: e = " << e << "\n";
  Point n = functions::computeOrthogonalVectorToBasis( B );
  std::cout << "Orthogonal(X) =: n = " << n << "\n";
  Point ns = functions::computeSimplifiedVector( n );
  std::cout << "Orthogonal(X)/gcd =: ns = " << ns << "\n";  
  for ( auto i = 0; i < B.size(); i++ )
    std::cout << "B[" << i << "] . ns = " << B[i] << " . " << ns
              << " == " << B[i].dot(ns) << " (should be 0)\n";
}

outputs

Dimension is 3
X = [ (-5, -25, 0) (19, 22, 7) (0, -19, 0) (-20, -17, 10) (20, 5, 0) (0, -6, 5) (9, 23, 12) (-14, -41, -2) (22, 10, 1) (34, -38, -23) (12, 11, 6) ]
Expected dimension of affine set of points X is 2
AffineDim(X) = 2
AffineSubset(X) = [ 0 1 2 ]
AffineBasis(X) =: p+B = (-5, -25, 0) + [ (24, 47, 7) (0, -13, -5) ]
Independent(X) =: e = (1, 0, 0)
Orthogonal(X) =: n = (-6, 5, -13)
Orthogonal(X)/gcd =: ns = (-6, 5, -13)
B[0] . ns = (24, 47, 7) . (-6, 5, -13) == 0 (should be 0)
B[1] . ns = (0, -13, -5) . (-6, 5, -13) == 0 (should be 0)

Using class AffineGeometry directly

If you need to finetune how affine geometry computations are performed, you may use the static methods of class AffineGeometry directly. The advantage is that you can specify a type for computations that is more precise than the type of input data.

typedef SpaceND<3, int> Space;
typedef Space::Point    Point;
typedef SpaceND<3, int64_t> ComputationSpace;
typedef ComputationSpace::Point  ComputationPoint;
...
typedef AffineGeometry<ComputationPoint> Affine;
 
std::vector<Point> X = { Point{1, 0, 0}, Point{2, 1, 0}, Point{3, 2, 0},
                         Point{3, 1, 1}, Point{5, 2, 2}, Point{4, 2, 1} };
auto I = Affine::affineSubset( X ); 
auto B = Affine::affineBasis( X );  
auto d = Affine::affineDimension( X ); 
 

The same services are provided as AffineGeometry static methods:

affine dimension and affine subset

AffineGeometry::affineDimension( const std::vector<TInputPoint>& X, const double tolerance ) Given a range of points X, returns the affine dimension of its spanned affine subspace.
AffineGeometry::affineSubset( const std::vector<TInputPoint>& X, const double tolerance ) Given a range of points X, returns a subset of these points that form an affine basis of X. Equivalently it is a simplex whose affine space spans all the points of X.
AffineGeometry::affineSubset( const std::vector<TInputPoint>& X, const TIndexRange& I, const double tolerance ) Given a range of points X and a subset of it given by indices I, returns a subset of these points that form an affine basis of X[I]. Equivalently it is a simplex whose affine space spans all the points of X[I].

affine basis and independent vectors

The computed bases are in row echelon form (see Wikipedia ), even reduced row echelon form if components are real-valued and not integers, except when using completeBasis where the last vector may be asked to be orthogonal to all the others.

AffineGeometry::affineBasis( const std::vector<TInputPoint>& X, const double tolerance ) Given a range of points X, returns a point and a range of vectors forming an affine basis containing X.
AffineGeometry::affineBasis( const std::vector<TInputPoint>& X, const TIndexRange& I, const double tolerance ) Given a range of points X and a range of indices I describing its subset of interest, returns a point and a range of vectors forming an affine basis containing the subset X[I].
AffineGeometry::independentVector Given a partial basis of vectors, returns a new vector that is independent.
AffineGeometry::completeBasis Complete the vectors basis with independent vectors so as to form a basis of the space. When normal_vector is true, the last added vector is guaranteed to be orthogonal to all the previous vectors.

orthogonal vector to basis of codimension 1

AffineGeometry::orthogonalVector Given d-1 independent vectors in dD, returns a vector that is orthogonal to each of them.

utilities

Reducing a vector onto a basis is a Gaussian elimination but slightly modified to handle integers (in case of lattice points and vectors).

AffineGeometry::reductionOnBasis Reduces the vector v on the (partial or not) basis of vectors basis, and returns the part of v that cannot be expressed as a linear combination of vectors of basis, and hence a null vector if v is a linear combination of the vectors of the basis.
AffineGeometry::reduceVector Reduces vector w by the vector b, and returns the coefficients \( (\alpha,\beta) \) for reduction such that \( \alpha w - \beta b \) is the returned reduced vector, \( \alpha \ge 0 \).
AffineGeometry::simplifiedVector Given a vector, returns the aligned vector with its component simplified by the gcd of all components (when the components are integers) or the aligned vector with unit L2-norm (when the components are floating-point numbers).

Projecting points onto an affine basis

AffineBasis is a convenience class that represents an affine basis (integral or not) in a form chosen by the user: either echelon reduced (AffineBasis::Type::ECHELON_REDUCED) or LLL reduced (AffineBasis::Type::ECHELON_REDUCED). It is mostly useful for projecting points and vectors onto it, and providing a system of coordinates on the affine space.

Warning: For now, you cannot project onto the basis when it is LLL-reduced (AffineBasis::Type::LLL_REDUCED). Use AffineBasis::Type::ECHELON_REDUCED instead.

Constructors

You can construct an affine basis from a set of points or from an origin and arbitrary vectors, while choosing its form in AffineBasis::Type.

AffineBasis::AffineBasis( const double tolerance ) Constructs the canonic basis
AffineBasis::AffineBasis( const std::vector< TInputPoint >& points, AffineBasis::Type type, const double delta, const double tolerance ) Constructs an affine basis spanning the set of points.
AffineBasis::AffineBasis( const TInputPoint& origin, const std::vector<TInputPoint>& basis, AffineBasis::Type type, bool is_reduced, const double delta, const double tolerance ) Constructs an affine basis spanning the given vectors and going through the given origin.
AffineBasis::AffineBasis( const TInputPoint& origin, const TInputPoint& normal, AffineBasis::Type type, const double delta, const double tolerance ) Creates an affine basis of codimension 1, going through origin and orthogonal to the lattice vector normal.

Standard services

AffineBasis::dimension() returns the dimension of the spanned affine subspace.
AffineBasis::origin() returns a point on the spanned affine subspace.
AffineBasis::basis() returns the vectors spanning the affine subspace.

Geometric services

AffineBasis::isParallel returns 'true' iff the vectors w lies in this affine space.
AffineBasis::isOnAffineSpace returns 'true' iff the point p lies in this affine space.

Projection services

AffineBasis::decompose Decomposes a point p into a tuple \( (d,L,R) \), such that: \( d (p - o) = L \cdot B + R \), where R is independent from this basis \( B=(b_0, \ldots, b_i) \), with origin o. d is a scalar, \( L/d \) are the rational coordinates of \( (p-o) \) in the basis, R is the remainer vector. Note that R is not null whenever p does not belong to the affine space.
AffineBasis::decomposeVector Same as above but for a vector (ignores the origin).
AffineBasis::recompose and AffineBasis::recomposeVector are the inverse functions of the two previous ones.
AffineBasis::projectPoints Projects a range of points onto the affine basis and outputs it. It also returns the maximum dilation of projected rational coordinates which was used to create lattice coordinates.

The example below shows how to use AffineBasis to project 4D points onto a 2D affine subspace.

typedef SpaceND< 4, int64_t >            Space4;
typedef Space4::Point                    Point4;
typedef SpaceND< 2, int64_t >            Space2;
typedef Space2::Point                    Point2;
typedef AffineBasis< Point4 >            Basis;
...
std::vector< Point4 > V = { Point{ 3, 4, 0, 2 }, Point{ -2, -1, 5, -7 } };
Point4 o( 1,4,2,-1 );
std::vector< Point4 > X = { o, o + V[0], o + V[1], o + 3*V[0] - V[1], o - V[0] + V[1] };
Basis B( X, Basis::Type::SHORTEST_ECHELON_REDUCED ); // affine basis of dimension 2
std::vector< Point2 >     proj_X;
auto d = B.projectPoints( proj_X, X ); // proj_X contains the projection of X onto B
// d is the dilatation factor, hence proj_X/d corresponds to rational coordinates

See also: testAffineBasis.cpp

Notes

Several functions that works with lattice vectors use functions specialized to integer computations (defined in math/linalg/IntegerMatrixFunctions.h:

Table of Contents

Generic functions for affine geometry

Using class AffineGeometry directly

Projecting points onto an affine basis

Notes