exampleAffineGeometry.cpp

Go to the documentation of this file.

 
#include <iostream>
#include <vector>
#include <random>
#include "DGtal/base/Common.h"
#include "DGtal/geometry/tools/AffineGeometry.h"
#include "ConfigExamples.h"
 
using namespace DGtal;
 
std::random_device rd;
std::mt19937 g(rd());
 
template < typename Point >
std::vector< Point >
makeRandomLatticePointsFromDirVectors( int nb, const std::vector< Point>& V )
{
  std::uniform_int_distribution<int> U(-10, 10);
  std::vector< Point > P;
  int n = V[0].size();
  int m = V.size();
  Point A;
  for ( auto i = 0; i < n; i++ )
    A[ i ] = U( g );
  P.push_back( A );
  for ( auto k = 0; k < nb; k++ )
    {
      Point B = A;
      for ( auto i = 0; i < m; i++ )
        {
          int l = U( g );
          B += l * V[ i ];
        }
      P.push_back( B );
    }
  std::shuffle( P.begin(), P.end(), g );
  return P;
}
 
template <typename T>
std::ostream&
operator<<( std::ostream& out, const std::vector< T >& object )
{
  out << "[";
  for ( auto t : object ) out << " " << t;
  out << " ]";
  return out;
}
 
template < typename Point >
bool checkAffineGeometry( int nb, const std::vector< Point>& V )
{
  auto X = makeRandomLatticePointsFromDirVectors( nb, V );
  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";
  return true;
}
 
 
int main( int argc, char* argv[] )
{
  int dim = argc > 1 ? atoi( argv[ 1 ] ) : 3;  // dimension of the space
  if ( dim > 5 ) dim = 5;
  if ( dim < 2 ) dim = 2;
  int nb  = argc > 2 ? atoi( argv[ 2 ] ) : 10;  // number of points
  if ( dim == 2 )
    {
      typedef PointVector< 2, int64_t> Point;
      std::vector< Point > X = { Point{3,1}  };
      return checkAffineGeometry( nb, X );
    }
  if ( dim == 3 )
    {
      typedef PointVector< 3, int64_t> Point;
      std::vector< Point > X = { Point{3,1,-1}, Point{-1,4,2} };
      return checkAffineGeometry( nb, X );
    }
  if ( dim == 4 )
    {
      typedef PointVector< 4, int64_t> Point;
      std::vector< Point > X = { Point{3,1,-1,2}, Point{-1,4,2,0}, Point{0,5,3,1} };
      return checkAffineGeometry( nb, X );
    }
  if ( dim == 5 )
    {
      typedef PointVector< 5, int64_t> Point;
      std::vector< Point > X = { Point{3,1,-1,2,3}, Point{-1,4,2,0,-2}, Point{0,5,3,1,-1}, Point{-4,-3,2,1,0} };
      return checkAffineGeometry( nb, X );
    }
}