testAffineGeometry.cpp

Go to the documentation of this file.

#include <iostream>
#include <vector>
#include <algorithm>
#include <random>
#include "DGtal/base/Common.h"
#include "DGtal/kernel/SpaceND.h"
#include "DGtal/geometry/tools/AffineGeometry.h"
#include "DGtal/geometry/tools/AffineBasis.h"
#include "DGtalCatch.h"
 
using namespace std;
using namespace DGtal;
 
std::random_device rd;
std::mt19937 g(rd());
std::uniform_real_distribution<double> uniform(-1.0, 1.0);
 
template <typename RealPoint>
void perturb( RealPoint& x, double perturbation )
{
  for ( auto& c : x ) c += uniform( g ) * perturbation;
}
 
template <typename RealPoint>
void perturb( std::vector< RealPoint >& X, double perturbation )
{
  for ( auto& x : X ) perturb( x, perturbation );
}
 
template < typename Point >
std::vector< Point >
makeRandomVectors( int nb, int amplitude )
{
  std::uniform_int_distribution<int> U(-amplitude, amplitude);
  std::vector< Point > P;
  for ( auto n = 0; n < nb; ++n )
    {
      Point A;
      for ( auto i = 0; i < Point::dimension; i++ )
        A[ i ] = U( g );
      P.push_back( A );
    }
  return P;
}
 
template < typename Point >
std::vector< Point >
makeRandomLatticePointsFromDirVectors( int nb, const vector< Point>& V )
{
  std::uniform_int_distribution<int> U(-10, 10);
  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 Point >
std::vector< Point >
makeRandomRealPointsFromDirVectors( int nb, const vector< Point>& V )
{
  std::uniform_real_distribution<double> U(-1., 1.);
  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++ )
        {
          double l = 5.* U( g );
          B += l * V[ i ];
        }
      P.push_back( B );
    }
  std::shuffle( P.begin(), P.end(), g );
  return P;
}

// Functions for testing class AffineGeometry in 2D.
 
SCENARIO( "AffineGeometry< Point2i > unit tests", "[affine_subset][2i]" )
{
  typedef SpaceND< 2, int >                Space;
  typedef Space::Point                     Point;
  typedef AffineGeometry< Point >            Affine;
  GIVEN( "Given X = { (0,0), (-4,-1), (16,4), (-3,5), (7,3), (5, -2) } of affine dimension 2" ) {
    std::vector<Point> X
      = { Point(0,0), Point(-4,-1), Point(16,4), Point(-3,5), Point(7,3), Point(5, -2) };
    auto I = Affine::affineSubset( X );
    THEN( "It has an affine basis of 3 points" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 3 );
      REQUIRE( I[ 2 ] == 3 );
    }
  }
  GIVEN( "Given X = { (0,0), (-4,-1), (-8,-2), (8,2), (16,4), (200,50) } of affine dimension 1" ) {
    std::vector<Point> X
      = { Point(0,0), Point(-4,-1), Point(-8,-2), Point(8,2), Point(16,4), Point(200,50) };
    auto I = Affine::affineSubset( X );
    THEN( "It has an affine basis of 2 points" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 2 );
    }
  }
}
 
SCENARIO( "AffineGeometry< Point2d > unit tests", "[affine_subset][2d]" )
{
  typedef SpaceND< 2, int >                Space;
  typedef Space::RealPoint                 Point;
  typedef AffineGeometry< Point >            Affine;
  GIVEN( "Given X = { (0,0), (-4,-1), (16,4), (-3,5), (7,3), (5, -2) } of affine dimension 2" ) {
    std::vector<Point> X
      = { Point(0,0), Point(-4,-1), Point(16,4), Point(-3,5), Point(7,3), Point(5, -2) };
    auto I = Affine::affineSubset( X );
    THEN( "It has an affine basis of 3 points [0,1,3]" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 3 );
      REQUIRE( I[ 2 ] == 3 );
    }
  }
  GIVEN( "Given X = { (0,0), (-4,-1), (-8,-2), (8,2), (16,4), (200,50) } of affine dimension 1" ) {
    std::vector<Point> X
      = { Point(0,0), Point(-4,-1), Point(-8,-2), Point(8,2), Point(16,4), Point(200,50) };
    auto I = Affine::affineSubset( X );
    THEN( "It has an affine basis of 2 points [0,1]" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 2 );
    }
  }
  GIVEN( "Given a perturbated X = { (0,0), (-4,-1), (16,4), (-3,5), (7,3), (5, -2) } of affine dimension 2 by U[-1e-4,1e-4]" ) {
    std::vector<Point> X
      = { Point(0,0), Point(-4,-1), Point(16,4), Point(-3,5), Point(7,3), Point(5, -2) };
    perturb( X, 1e-4 );
    auto I = Affine::affineSubset( X, 1e-12 );
    THEN( "It has an affine basis of 3 points [0,1,x]" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 3 );
    }
  }
  GIVEN( "Given a perturbated X = { (0,0), (-4,-1), (-8,-2), (8,2), (16,4), (200,50) } of affine dimension 1 by U[-1e-4,1e-4]" ) {
    std::vector<Point> X
      = { Point(0,0), Point(-4,-1), Point(-8,-2), Point(8,2), Point(16,4), Point(200,50) };
    perturb( X, 1e-4 );
    auto I = Affine::affineSubset( X, 1e-12 );
    THEN( "It has an affine basis of 3 points [0,1,x]" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 3 );
    }
  }
  GIVEN( "Given a perturbated X = { (0,0), (-4,-1), (-8,-2), (8,2), (16,4), (200,50) } of affine dimension 1 by U[-1e-11,1e-11]" ) {
    std::vector<Point> X
      = { Point(0,0), Point(-4,-1), Point(-8,-2), Point(8,2), Point(16,4), Point(200,50) };
    perturb( X, 1e-11 );
    auto I = Affine::affineSubset( X, 1e-7 );
    THEN( "It has an affine basis of 2 points [0,1] if tolerance is 1e-7" ) {
      CAPTURE( X );
      CAPTURE( I );
      REQUIRE( I.size() == 2 );
    }
  }
}
 
SCENARIO( "AffineGeometry< Point2i > orthogonal tests", "[orthogonal_vector][2i]" )
{
  typedef SpaceND< 2, int >                Space;
  typedef Space::Point                     Point;
  typedef AffineGeometry< Point >          Affine;
  GIVEN( "Given basis B = { (7,3) } " ) {
    std::vector<Point> B = { Point(7,3) };
    Affine::completeBasis( B, true );
    THEN( "The complete basis has dimension 2" ) {
      CAPTURE( B );
      REQUIRE( B.size() == 2 );
    }
    THEN( "The last vector is non null and orthogonal to all the others" ) {
      CAPTURE( B.back() );
      REQUIRE( B.back().normInfinity() > 0 );
      REQUIRE( B.back().dot( B[ 0 ] ) == 0 );
    }
  }
}

// Functions for testing class AffineGeometry in 3D.
 
SCENARIO( "AffineGeometry< Point3i > unit tests", "[affine_subset][3i]" )
{
  typedef SpaceND< 3, int >                Space;
  typedef Space::Point                     Point;
  typedef AffineGeometry< Point >            Affine;
  GIVEN( "Given X = { (1, 0, 0), (2, 1, 0), (3, 1, 1), (3, 2, 0), (5, 2, 2), (4, 2, 1)} of affine dimension 2" ) {
    std::vector<Point> X
      = { Point{1, 0, 0}, Point{2, 1, 0}, Point{3, 1, 1}, Point{3, 2, 0}, Point{5, 2, 2}, Point{4, 2, 1} };
 
    auto I = Affine::affineSubset( X );
    THEN( "It has an affine basis of 3 points" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 3 );
      REQUIRE( I[ 2 ] == 2 );
    }
  }
  GIVEN( "Given X = { (1, 0, 0), (2, 1, 0), (3, 1, 1), (3, 2, 0), (5, 2, 2), (4, 2, 1), (7, 3, 2)} of affine dimension 3" ) {
    std::vector<Point> X
      = { Point{1, 0, 0}, Point{2, 1, 0}, Point{3, 1, 1}, Point{3, 2, 0}, Point{5, 2, 2}, Point{4, 2, 1}, Point{7, 3, 2} };
 
    auto I = Affine::affineSubset( X );
    THEN( "It has an affine basis of 4 points" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 4 );
      REQUIRE( I[ 2 ] == 2 );
      REQUIRE( I[ 3 ] == 6 );
    }
  }
  GIVEN( "Given X a set of randomly generated points by adding linear combinations of 1 lattice vectors" ) {
    std::vector< Point > V = { Point{ 3, 1, 0 } };
    auto X = makeRandomLatticePointsFromDirVectors( 20, V );
    auto I = Affine::affineSubset( X );
    THEN( "It has an affine basis of 2 points" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 2 );
    }
  }
  GIVEN( "Given X a set of randomly generated points by adding linear combinations of 2 lattice vectors" ) {
    std::vector< Point > V = { Point{ 3, 1, 0 }, Point{ -2, -1, 2 } };
    auto X = makeRandomLatticePointsFromDirVectors( 20, V );
    auto I = Affine::affineSubset( X );
    THEN( "It has an affine basis of 3 points" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 3 );
    }
  }
  GIVEN( "Given X a set of randomly generated points by adding linear combinations of 3 lattice vectors" ) {
    std::vector< Point > V = { Point{ 3, 1, 0 }, Point{ -2, -1, 2 }, Point{ -1, 4, 3 } };
    auto X = makeRandomLatticePointsFromDirVectors( 20, V );
    auto I = Affine::affineSubset( X );
    THEN( "It has an affine basis of 4 points" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 4 );
    }
  }
}
 
SCENARIO( "AffineGeometry< Point3i > orthogonal tests", "[orthogonal_vector][3i]" )
{
  typedef SpaceND< 3, int >                Space;
  typedef Space::Point                     Point;
  typedef AffineGeometry< Point >          Affine;
  GIVEN( "Given basis B = { (7,3,1), (2,-5,1) } " ) {
    std::vector<Point> B = { Point(7,3,1), Point( 2,-5,1) };
    Point n = B[ 0 ].crossProduct( B[ 1 ] );
    Affine::completeBasis( B, true );
    THEN( "The complete basis has dimension 3" ) {
      CAPTURE( B );
      REQUIRE( B.size() == 3 );
    }
    THEN( "The last vector is non null and orthogonal to all the others" ) {
      CAPTURE( B.back() );
      REQUIRE( B.back().normInfinity() > 0 );
      REQUIRE( B.back().dot( B[ 0 ] ) == 0 );
      REQUIRE( B.back().dot( B[ 1 ] ) == 0 );
    }
    THEN( "The last vector is the cross product of the two others" ) {
      CAPTURE( B.back() );
      REQUIRE( B.back() == n );
    }
  }
  GIVEN( "Given basis B = { (7,3,1) } " ) {
    std::vector<Point> B = { Point(7,3,1) };
    Affine::completeBasis( B, true );
    THEN( "The complete basis has dimension 3" ) {
      CAPTURE( B );
      REQUIRE( B.size() == 3 );
    }
    THEN( "The mid vector is non null and trivial" ) {
      CAPTURE( B[ 1 ] );
      REQUIRE( B[ 1 ].norm1() == 1 );
    }
    THEN( "The last vector is non null and orthogonal to all the others" ) {
      CAPTURE( B.back() );
      REQUIRE( B.back().normInfinity() > 0 );
      REQUIRE( B.back().dot( B[ 0 ] ) == 0 );
      REQUIRE( B.back().dot( B[ 1 ] ) == 0 );
    }
    THEN( "The last vector is the cross product of the two others" ) {
      Point n = B[ 0 ].crossProduct( B[ 1 ] );
      CAPTURE( B.back() );
      REQUIRE( B.back() == n );
    }
  }
}

// Functions for testing class AffineGeometry in 4D.
 
SCENARIO( "AffineGeometry< Point4i > unit tests", "[affine_subset][4i]" )
{
  typedef SpaceND< 4, int >                Space;
  typedef Space::Point                     Point;
  GIVEN( "Given X a set of randomly generated points by adding linear combinations of 1 lattice vectors" ) {
    std::vector< Point > V = { Point{ 3, 1, 0, 2 } };
    auto X = makeRandomLatticePointsFromDirVectors( 20, V );
    auto I = functions::computeAffineSubset( X );
    Point o;
    std::vector< Point > B;
    functions::getAffineBasis( o, B, X );
    THEN( "It has an affine basis of 2 points" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 2 );
    }
    THEN( "It has an affine basis of 1 vector" ) {
      REQUIRE( B.size() == 1 );
    }
 
  }
  GIVEN( "Given X a set of randomly generated points by adding linear combinations of 2 lattice vectors" ) {
    std::vector< Point > V = { Point{ 3, 1, 0, 2 }, Point{ -2, -1, 2, 7 } };
    auto X = makeRandomLatticePointsFromDirVectors( 20, V );
    auto I = functions::computeAffineSubset( X );
    Point o;
    std::vector< Point > B;
    functions::getAffineBasis( o, B, X );
    THEN( "It has an affine basis of 3 points" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 3 );
    }
    THEN( "It has an affine basis of 2 vectors" ) {
      REQUIRE( B.size() == 2 );
    }
  }
  GIVEN( "Given X a set of randomly generated points by adding linear combinations of 3 lattice vectors" ) {
    std::vector< Point > V = { Point{ 3, 1, 0, 2 }, Point{ -2, -1, 2, 7 }, Point{ -1, 4, 3, -1  } };
    auto X = makeRandomLatticePointsFromDirVectors( 20, V );
    auto I = functions::computeAffineSubset( X );
    Point o;
    std::vector< Point > B;
    functions::getAffineBasis( o, B, X );
    THEN( "It has an affine basis of 4 points" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 4 );
    }
    THEN( "It has an affine basis of 3 vectors" ) {
      REQUIRE( B.size() == 3 );
    }
  }
  GIVEN( "Given X a set of randomly generated points by adding linear combinations of 4 lattice vectors" ) {
    std::vector< Point > V = { Point{ 3, 1, 0, 2 }, Point{ -2, -1, 2, 7 }, Point{ -1, 4, 3, -1 }, Point{ 2, 1, -3, -4 } };
    auto X = makeRandomLatticePointsFromDirVectors( 20, V );
    auto I = functions::computeAffineSubset( X );
    Point o;
    std::vector< Point > B;
    functions::getAffineBasis( o, B, X, 1e-10 );
    THEN( "It has an affine basis of 5 points" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 5 );
    }
    THEN( "It has an affine basis of 4 vectors" ) {
      REQUIRE( B.size() == 4 );
    }
  }
}
 
SCENARIO( "AffineGeometry< Point4d > unit tests", "[affine_subset][4d]" )
{
  // NB: 1e-10 in tolerance is ok for these examples.
  // max norm of rejected vectors are 2e-12.
  typedef SpaceND< 4, int >                Space;
  typedef Space::RealPoint                 Point;
  GIVEN( "Given X a set of randomly generated points by adding linear combinations of 1 lattice vectors" ) {
    std::vector< Point > V = { Point{ 3, 1, 0, 2 } };
    auto X = makeRandomRealPointsFromDirVectors( 20, V );
    auto I = functions::computeAffineSubset( X, 1e-10 );
    Point o;
    std::vector< Point > B;
    functions::getAffineBasis( o, B, X, 1e-10 );
    THEN( "It has an affine subset of 2 points" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 2 );
    }
    THEN( "It has an affine basis of 1 vector" ) {
      REQUIRE( B.size() == 1 );
    }
  }
  GIVEN( "Given X a set of randomly generated points by adding linear combinations of 2 lattice vectors" ) {
    std::vector< Point > V = { Point{ 3, 1, 0, 2 }, Point{ -2, -1, 2, 7 } };
    auto X = makeRandomRealPointsFromDirVectors( 20, V );
    auto I = functions::computeAffineSubset( X, 1e-10 );
    Point o;
    std::vector< Point > B;
    functions::getAffineBasis( o, B, X, 1e-10 );
    THEN( "It has an affine subset of 3 points" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 3 );
    }
    THEN( "It has an affine basis of 2 vectors" ) {
      REQUIRE( B.size() == 2 );
    }
  }
  GIVEN( "Given X a set of randomly generated points by adding linear combinations of 3 lattice vectors" ) {
    std::vector< Point > V = { Point{ 3, 1, 0, 2 }, Point{ -2, -1, 2, 7 }, Point{ -1, 4, 3, -1  } };
    auto X = makeRandomRealPointsFromDirVectors( 20, V );
    auto I = functions::computeAffineSubset( X, 1e-10 );
    Point o;
    std::vector< Point > B;
    functions::getAffineBasis( o, B, X, 1e-10 );
    THEN( "It has an affine subset of 4 points" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 4 );
    }
    THEN( "It has an affine basis of 3 vectors" ) {
      REQUIRE( B.size() == 3 );
    }
  }
  GIVEN( "Given X a set of randomly generated points by adding linear combinations of 4 lattice vectors" ) {
    std::vector< Point > V = { Point{ 3, 1, 0, 2 }, Point{ -2, -1, 2, 7 }, Point{ -1, 4, 3, -1 }, Point{ 2, 1, -3, -4 } };
    auto X = makeRandomRealPointsFromDirVectors( 20, V );
    auto I = functions::computeAffineSubset( X, 1e-10 );
    Point o;
    std::vector< Point > B;
    functions::getAffineBasis( o, B, X, 1e-10 );
    THEN( "It has an affine subset of 5 points" ) {
      CAPTURE( I );
      REQUIRE( I.size() == 5 );
    }
    THEN( "It has an affine basis of 4 vectors" ) {
      REQUIRE( B.size() == 4 );
    }
  }
}
 
SCENARIO( "AffineGeometry< Point4i > orthogonal tests", "[orthogonal_vector][4i]" )
{
  typedef SpaceND< 4, int >                Space;
  typedef Space::Point                     Point;
  typedef AffineGeometry< Point >          Affine;
  GIVEN( "Given basis B = { (7,3,1,0), (2,-5,1,2), (-1,2,2,-3) } " ) {
    std::vector<Point> B = { Point(7,3,1,0), Point(2,-5,1,2), Point(-1,2,2,-3) };
    Affine::completeBasis( B, true );
    THEN( "The complete basis has dimension 4" ) {
      CAPTURE( B );
      REQUIRE( B.size() == 4 );
    }
    THEN( "The last vector is non null and orthogonal to all the others" ) {
      CAPTURE( B.back() );
      REQUIRE( B.back().normInfinity() > 0 );
      REQUIRE( B.back().dot( B[ 0 ] ) == 0 );
      REQUIRE( B.back().dot( B[ 1 ] ) == 0 );
      REQUIRE( B.back().dot( B[ 2 ] ) == 0 );
    }
  }
  GIVEN( "Given basis B = { (7,3,1,0), (2,-5,1,2) } " ) {
    std::vector<Point> B = { Point(7,3,1,0), Point(2,-5,1,2) };
    Affine::completeBasis( B, true );
    THEN( "The complete basis has dimension 4" ) {
      CAPTURE( B );
      REQUIRE( B.size() == 4 );
    }
    THEN( "The third vector is non null and trivial" ) {
      CAPTURE( B[ 2 ] );
      REQUIRE( B[ 2 ].norm1() == 1 );
    }
    THEN( "The last vector is non null and orthogonal to all the others" ) {
      CAPTURE( B.back() );
      REQUIRE( B.back().normInfinity() > 0 );
      REQUIRE( B.back().dot( B[ 0 ] ) == 0 );
      REQUIRE( B.back().dot( B[ 1 ] ) == 0 );
      REQUIRE( B.back().dot( B[ 2 ] ) == 0 );
    }
  }
}
 
 
SCENARIO( "AffineGeometry< Z3 > bug", "[affine_geom][3d]" )
{
  typedef SpaceND<3,int>          Space;
  typedef Space::Point            Point;
  std::vector< Point > X = { {-46, 38, -43}, {27, -89, 20}, {53, 26, -57} };
  Point o;
  std::vector< Point > B;
  functions::getAffineBasis ( o, B, X );
  auto  C         = ( X[1]-X[0] ).crossProduct( X[2]-X[0] );
  auto  sC        = functions::computeSimplifiedVector( C );
  WHEN( "Computing orthogonal vector" ) {
    auto N = functions::computeOrthogonalVectorToBasis( B );
    THEN( "It is non null" ) {
      CAPTURE( N );
      REQUIRE( N != Point::zero );
    }
    THEN( "It corresponds to the reduced cross product" ) {
      CAPTURE( C );
      CAPTURE( sC );
      REQUIRE( N == sC );
    }
  }
}
 
SCENARIO( "AffineGeometry< Z3 > orthogonality", "[affine_geom][3d]" )
{
  typedef SpaceND<3,int64_t>          Space;
  typedef Space::Point            Point;
 
  WHEN( "Computing orthogonal vector to multiple random points" ) {
    std::vector< Point > X = makeRandomVectors<Point>( 100, 50 );
    std::size_t nb        = 300;
    std::size_t nb_equal  = 0;
    std::size_t nb_C_zero = 0;
    std::size_t nb_N_zero = 0;
    for ( auto i = 0; i < nb; i++ )
      {
        std::vector< Point > Y = { X[ rand() % 100 ], X[ rand() % 100 ], X[ rand() % 100 ] };
        auto  C         = ( Y[1]-Y[0] ).crossProduct( Y[2]-Y[0] );
        auto  sC        = functions::computeSimplifiedVector( C );
        Point N;
        functions::getOrthogonalVector( N, Y );
        nb_equal += (sC == N) ? 1 : 0;
        if ( sC != N )
          std::cout << "Y = " << Y[0] << " " << Y[1] << " " << Y[ 2 ]
                    << " sC = " << sC << " N = " << N << "\n";
        nb_C_zero += ( C == Point::zero ) ? 1 : 0;
        nb_N_zero += ( N == Point::zero ) ? 1 : 0;
      }
    THEN( "They corresponds to the reduced cross product" ) {
      REQUIRE( nb_equal == nb );
      REQUIRE( nb_C_zero == nb_N_zero );
    }
  }
}