AffineGeometry.h

 
#pragma once
 
#if defined(AffineGeometry_RECURSES)
#error Recursive header files inclusion detected in AffineGeometry.h
#else // defined(AffineGeometry_RECURSES)
#define AffineGeometry_RECURSES
 
#if !defined AffineGeometry_h
#define AffineGeometry_h
 
// Inclusions
#include <vector>
#include "DGtal/base/Common.h"
#include "DGtal/kernel/IntegerConverter.h"
#include "DGtal/kernel/PointVector.h"
#include "DGtal/arithmetic/IntegerComputer.h"
#include "DGtal/math/linalg/SimpleMatrix.h"
#include "DGtal/math/linalg/IntegerMatrixFunctions.h"
 
namespace DGtal
{
  // Forward declaration of AffineGeometry (needed for friend declaration).
  template < typename T > struct AffineGeometry;
  // Forward declaration of AffineBasis (needed for friend declaration).
  template < typename T > struct AffineBasis;
 
  namespace detail {
 
    // template class AffineGeometryPointOperations
 
    template < DGtal::Dimension dim,
               typename TEuclideanRing,
               typename TContainer=std::array<TEuclideanRing,dim> >
    struct AffineGeometryPointOperations
    {
      template <typename T> friend struct DGtal::AffineGeometry;
      template <typename T> friend struct DGtal::AffineBasis;
      typedef PointVector< dim, TEuclideanRing, TContainer > Point;
      typedef typename Point::Coordinate Scalar;
      typedef Scalar Integer; 
 
      // ----------------------- internal services --------------------------
    private:
      
      template <typename TInteger >
      static void normalizeVector( Point& w, TInteger )
      {
        DGtal::Dimension i = 0;
        while ( i < Point::dimension && w[ i ] == 0 ) i++;
        if ( i == Point::dimension ) return;
        TInteger g = abs( w[ i ] );
        for ( DGtal::Dimension k = i+1; k < Point::dimension; k++ )
          g = DGtal::IntegerComputer< TInteger >::staticGcd( g, abs( w[ k ] ) );
        for ( DGtal::Dimension k = i; k < Point::dimension; k++ )
          w[ k ] /= g;
      }
      
      template <typename TOtherPoint>
      static
      Point cast( const TOtherPoint& other )
      {
        return IntegerConverter< Point::dimension, Scalar >::cast( other );
      }
    }; // end of class AffineGeometryPointOperations
    
    
    template < DGtal::Dimension dim,
               typename TContainer >
    struct AffineGeometryPointOperations< dim, double, TContainer >
    {
      template <typename T> friend struct DGtal::AffineGeometry;
      template <typename T> friend struct DGtal::AffineBasis;
      typedef PointVector< dim, double, TContainer > Point;
      typedef typename Point::Coordinate Scalar;
 
      // ----------------------- internal services --------------------------
    private:
 
      static
      void normalizeVector( Point& w, double x )
      {
        w /= sqrt( x );
      }
 
      template <typename TOtherPoint>
      static
      Point cast( const TOtherPoint& other )
      {
        typedef typename TOtherPoint::Coordinate OtherScalar;
        Point result;
        for ( std::size_t i = 0; i < Point::dimension; i++ )
          result[ i ] = NumberTraits<OtherScalar>::castToDouble( other[ i ] );
        return result;
      }
    }; // end of class AffineGeometryPointOperations
 
    template < DGtal::Dimension dim,
               typename TContainer >
    struct AffineGeometryPointOperations< dim, float, TContainer >
    {
      template <typename T> friend struct DGtal::AffineGeometry;
      template <typename T> friend struct DGtal::AffineBasis;
      typedef PointVector< dim, float, TContainer > Point;
      typedef typename Point::Coordinate Scalar;
 
      // ----------------------- internal services --------------------------
    private:
 
      static
      void normalizeVector( Point& w, float x )
      {
        w /= sqrt( x );
      }
 
      template <typename TOtherPoint>
      static
      Point cast( const TOtherPoint& other )
      {
        Point result;
        for ( std::size_t i = 0; i < Point::dimension; i++ )
          result[ i ] = float( other[ i ] );
        return result;
      }
    }; // end of class AffineGeometryPointOperations
    
    
 
    // template class AffineGeometryScalarOperations
 
    template < typename TScalar >
    struct AffineGeometryScalarOperations
    {
      template <typename T> friend struct DGtal::AffineGeometry;
      template <typename T> friend struct DGtal::AffineBasis;
 
      typedef TScalar Scalar;
      typedef Scalar Integer; 
      
      // ----------------------- internal services --------------------------
    private:
 
      static
      std::pair< Integer, Integer > getMultipliers( Integer a, Integer b )
      {
        const Integer g = gcd( a, b );
        // return std::make_pair( b / g, a / g );
        // left multiplier should be positive.
        return (b >= 0)
          ? std::make_pair( b / g, a / g )
          : std::make_pair( -b / g, -a / g );
      }
 
      static
      Integer gcd( Integer a, Integer b )
      {
        return IntegerComputer< Integer >::staticGcd( abs( a ), abs( b ) );
      }
 
      static
      Integer lcmPositive( Integer a, Integer b )
      {
        const Integer g = gcd( a, b );
        return (a / g) * b;
      }
      
      static
      bool isNonZero( Integer x, double )
      {
        return x != Integer( 0 );
      }
 
    }; // end of class AffineGeometryScalarOperations
 
    template <>
    struct AffineGeometryScalarOperations< double >
    {
      template <typename T> friend struct DGtal::AffineGeometry;
      template <typename T> friend struct DGtal::AffineBasis;
    
      static
      std::pair< double, double > getMultipliers( double a, double b )
      {
        return (b >= 0 )
          ? std::make_pair( b, a )
          : std::make_pair( -b, -a );
      }
 
      static
      double gcd( double, double )
      {
        return 1.0;
      }
 
      static
      double lcmPositive( double, double )
      {
        return 1.0;
      }
 
      
      static
      bool isNonZero( double x, double tol )
      {
        return ( x > tol ) || ( x < -tol );
      }
 
    }; // end of class AffineGeometryScalarOperations< double >
 
    template <>
    struct AffineGeometryScalarOperations< float >
    {
      template <typename T> friend struct DGtal::AffineGeometry;
      template <typename T> friend struct DGtal::AffineBasis;
    
      static
      std::pair< float, float > getMultipliers( float a, float b )
      {
        return (b >= 0 )
          ? std::make_pair( b, a )
          : std::make_pair( -b, -a );
      }
 
      static
      float gcd( float, float )
      {
        return 1.0f;
      }
 
      static
      float lcmPositive( float, float )
      {
        return 1.0f;
      }
 
      
      static
      bool isNonZero( float x, double tol )
      {
        const double dx = double(x);
        return (dx > tol) || ( dx < -tol );
      }
      
    }; // end of class AffineGeometryScalarOperations< float >
 
    template < typename TScalar, bool Robust >
    struct AffineGeometryInternalNumber {
      typedef TScalar type; 
    };
    template <>
    struct AffineGeometryInternalNumber<int32_t, false> {
      typedef int64_t type; 
    };
    template <>
    struct AffineGeometryInternalNumber<int32_t, true> {
      typedef BigInteger type; 
    };
    template <>
    struct AffineGeometryInternalNumber<int64_t, false> {
      typedef int64_t type; 
    };
    template <>
    struct AffineGeometryInternalNumber<int64_t, true> {
      typedef BigInteger type; 
    };
    template <>
    struct AffineGeometryInternalNumber<float, false> {
      typedef double type; 
    };
    template <>
    struct AffineGeometryInternalNumber<float, true> {
      typedef double type; 
    };
    template <>
    struct AffineGeometryInternalNumber<double, false> {
      typedef double type; 
    };
    template <>
    struct AffineGeometryInternalNumber<double, true> {
      typedef double type; 
    };
    
  } // namespace detail
 
} // namespace DGtal
 
 
namespace DGtal
{
  // template class AffineGeometry
 
  template < typename TPoint >
  struct AffineGeometry
  {
    typedef TPoint               Point;
    typedef typename Point::Coordinate Scalar;
    typedef typename Point::Container  Container;
    typedef std::size_t                Size;
    typedef std::vector< Point > Points; 
    typedef std::vector< Size >        Sizes;  
    static const Size  dimension  = Point::dimension;
    typedef DGtal::detail::AffineGeometryPointOperations
    < dimension, Scalar, Container >   PointOps;
    typedef DGtal::detail::AffineGeometryScalarOperations< Scalar > ScalarOps;
    
    // ----------------------- standard services --------------------------
  public:
 
    template <typename TInputPoint>
    static
    DGtal::int64_t
    affineDimension( const std::vector<TInputPoint>& X, const double tolerance = 1e-12 )
    {
      return DGtal::int64_t( affineSubset( X, tolerance ).size() ) - 1;
    }
    
    template <typename TInputPoint>
    static
    std::vector< Size >
    affineSubset( const std::vector<TInputPoint>& X, const double tolerance = 1e-12 )
    {
      Size m = X.size();
      // Process trivial cases.
      if ( m == 0 ) return { };
      if ( m == 1 ) return { Size( 0 ) };
      // Process general case.
      Points basis;  //< direction vectors
      Sizes  chosen; //< selected points
      chosen.push_back( 0 ); //< reference point (first one, as it may be any one)
      for ( Size i = 1; i < m; i++ )
        {
          Point v = transform( X[ i ] - X[ 0 ] );
          if ( addIfIndependent( basis, v, tolerance ) )
            chosen.push_back( i );
          if ( chosen.size() > dimension ) break;
        }
      return chosen;
    }
 
    template <typename TInputPoint, typename TIndexRange>
    static
    std::vector< Size >
    affineSubset( const std::vector<TInputPoint>& X,
                  const TIndexRange& I,
                  const double tolerance = 1e-12 )
    {
      Size m = I.size();
      // Process trivial cases.
      if ( m == 0 ) return { };
      if ( m == 1 ) return { I[ 0 ] };
      // Process general case.
      Points basis;  //< direction vectors
      Sizes  chosen; //< selected points
      chosen.push_back( I[ 0 ] ); //< reference point (first one, as it may be any one)
      for ( Size i = 1; i < m; i++ )
        {
          Point v = transform( X[ I[ i ] ] - X[ I[ 0 ] ] );
          if ( addIfIndependent( basis, v, tolerance ) )
            chosen.push_back( I[ i ] );
          if ( chosen.size() > dimension ) break;
        }
      return chosen;
    }
    
    template <typename TInputPoint>
    static
    std::pair< Point, Points >
    affineBasis( const std::vector<TInputPoint>& X, const double tolerance = 1e-12 )
    {
      Points basis;  //< direction vectors
      Size m = X.size();
      // Process trivial cases.
      if ( m == 0 ) return std::make_pair( Point::zero, basis );
      const Point o = transform( X[ 0 ] );
      if ( m == 1 ) return std::make_pair( o, basis );
      // Process general case.
      basis.reserve( dimension );
      for ( Size i = 1; i < m; i++ )
        {
          Point v = transform( X[ i ] ) - o;
          if ( addIfIndependent( basis, v, tolerance ) )
            if ( basis.size() >= dimension ) break;
        }
      return std::make_pair( o, basis );
    }
 
    template < typename TInputPoint, typename TIndexRange >
    static
    std::pair< TInputPoint, Points >
    affineBasis( const std::vector<TInputPoint>& X,
                 const TIndexRange& I,
                 const double tolerance = 1e-12 )
    {
      Points basis;  //< direction vectors
      Size m = I.size();
      // Process trivial cases.
      if ( m == 0 ) return std::make_pair( TInputPoint::zero, basis );
      if ( m == 1 ) return std::make_pair( X[ I[ 0 ] ], basis );
      // Process general case.
      const Point  o = transform( X[ I[ 0 ] ] );
      basis.reserve( dimension );
      for ( Size i = 1; i < m; i++ )
        {
          Point v = transform( X[ I[ i ] ] ) - o;
          if ( addIfIndependent( basis, v, tolerance ) )
            if ( basis.size() > dimension ) break;
        }
      return std::make_pair( X[ I[ 0 ] ], basis );
    }
    
    static
    Point
    independentVector( const Points& basis, const double tolerance = 1e-12 )
    {
      // If basis has already d independent vectors, then there is no
      // other independent vector.
      if ( basis.size() >= dimension ) return Point::zero;
      // At least one trivial canonic vector should be independant.
      Dimension k = 0;
      for ( ; k < dimension; k++ )
        {
          Point e_k = Point::base( k );
          Point w   = reductionOnBasis( e_k, basis, tolerance );
          Scalar sql2;
          functions::getSquaredNormL2( sql2, w );
          if ( ScalarOps::isNonZero( sql2, tolerance ) )
            return e_k;
        }
      trace.error() << "[AffineGeometry::independentVector]"
                    << " Unable to find independent vector." << std::endl;
      return Point::zero;
    }
 
    template < typename TOtherPoint >
    static
    TOtherPoint
    independentVector( const Points& basis, const double tolerance = 1e-12 )
    {
      // If basis has already d independent vectors, then there is no
      // other independent vector.
      if ( basis.size() >= dimension ) return TOtherPoint::zero;
      // At least one trivial canonic vector should be independant.
      Dimension k = 0;
      for ( ; k < dimension; k++ )
        {
          Point e_k = Point::base( k );
          Point w   = reductionOnBasis( e_k, basis, tolerance );
          Scalar sql2;
          functions::getSquaredNormL2( sql2, w );
          if ( ScalarOps::isNonZero( sql2, tolerance ) )
            return TOtherPoint::base( k ); 
        }
      trace.error() << "[AffineGeometry::independentVector]"
                    << " Unable to find independent vector." << std::endl;
      return TOtherPoint::zero;
    }
 
    static
    void
    completeBasis( Points& basis,
                   bool normal_vector,
                   const double tolerance = 1e-12 )
    {
      if ( basis.size() >= dimension ) return;
      while ( ( basis.size() + 1 ) < dimension )
        { // add an independent vector
          const Point u = independentVector( basis, tolerance );
          basis.push_back( u );
        }
      if ( normal_vector )
        basis.push_back( orthogonalVector( basis ) );
      else
        basis.push_back( independentVector( basis, tolerance ) );
    }
 
    template < typename TInputPoint >
    static
    std::vector< Point >
    orthogonalLatticeBasis( const TInputPoint& N, bool shortened = false )
    {
      std::vector< Scalar > Np( N.begin(), N.end() );
      Np.resize( Point::dimension );
      auto B = functions::computeOrthogonalLattice( Np );
      if ( shortened )
        functions::shortenBasis( B );
      std::vector< Point > R( B.size() );
      for ( std::size_t i = 0; i < B.size(); i++ )
        for ( std::size_t j = 0; j < Point::dimension; j++ )
          R[ i ][ j ] = B[ i ][ j ];
      return R;
    }
      
 
    // ----------------------- specific services --------------------------
  public:
 
    static
    Point
    reductionOnBasis( const Point& v,
                      const Points& basis,
                      const double tolerance )
    {
      Point w( v );
      for ( const auto& b : basis ) reduceVector( w, b, tolerance );
      return w;
    }
 
    static
    bool addIfIndependent( Points& basis,
                           const Point& v,
                           const double tolerance )
    {
      Point  w = reductionOnBasis( v, basis, tolerance );
      Scalar x;
      functions::getSquaredNormL2( x, w );
      if ( ScalarOps::isNonZero( x, tolerance ) )
        {
          // Useful to reduce the norm of vectors for lattice vectors
          // and necessary for real vectors so that `tolerance` keeps
          // the same meaning.
          PointOps::normalizeVector( w, x ); 
          basis.push_back( w );
          return true;
        }
      return false;
    }
    
    static
    std::pair< Scalar, Scalar >
    reduceVector( Point& w, const Point& b, const double tolerance )
    {
      Scalar mul_w = 0;
      Scalar mul_b = 0;
      Size n = w.size();
      // Find index of first non null pivot in b.
      Size lead = n;
      for ( Size j = 0; j < n; j++)
        if ( ScalarOps::isNonZero( b[j], tolerance ) ) { lead = j; break; }
      if ( lead == n ) return std::make_pair( mul_w, mul_b ); // b is null vector
 
      std::tie( mul_w, mul_b ) = ScalarOps::getMultipliers( w[ lead ], b[ lead ] );
      
      for (Size j = 0; j < n; j++) 
        w[j] = mul_w * w[j] - mul_b * b[j];
      return std::make_pair( mul_w, mul_b );
    }
 
    static
    std::pair< Scalar, Scalar >
    reduceVector( Point& w, const Point& b,
                  Dimension start, const double tolerance )
    {
      Scalar mul_w = 0;
      Scalar mul_b = 0;
      Size n = w.size();
      // Find index of first non null pivot in b.
      Size lead = n;
      for ( Size j = start; j < n; j++)
        if ( ScalarOps::isNonZero( b[j], tolerance ) ) { lead = j; break; }
      if ( lead == n ) return std::make_pair( mul_w, mul_b ); // b is null vector
 
      std::tie( mul_w, mul_b ) = ScalarOps::getMultipliers( w[ lead ], b[ lead ] );
      
      for (Size j = 0; j < n; j++) 
        w[j] = mul_w * w[j] - mul_b * b[j];
      return std::make_pair( mul_w, mul_b );
    }
 
    
    static
    const Point&
    transform( const Point& w )
    {
      return w;
    }
 
    template <typename TInputPoint>
    static
    Point
    transform( const TInputPoint& w )
    {
      return PointOps::cast( w );
    }
    
    static
    Point
    orthogonalVector( const Points& basis )
    {
      Point w;
      const std::size_t n = dimension;
      if ( ( basis.size() + 1 ) != dimension ) return w;
      const std::size_t m = basis.size();
      SimpleMatrix< Scalar, dimension-1, dimension> A;
      for ( std::size_t i = 0; i < m; ++i )
        for ( std::size_t j = 0; j < n; ++j )
          A( i, j ) = basis[ i ][ j ];
      for ( std::size_t col = 0; col < n; ++col)
        { // construct sub-matrix removing column col
          SimpleMatrix< Scalar, dimension-1, dimension-1> M;
          for ( std::size_t i = 0; i < m; ++i )
            {
              std::size_t c = 0;
              for ( std::size_t j = 0; j < n; ++j)
                if ( j != col ) M( i, c++ ) = A( i, j );
            }
          functions::getDeterminantBareiss( w[ col ], M );
          if ( (col+dimension) % 2 == 0 ) w[ col ] = -w[ col ];
        }
      return simplifiedVector( w );
    }
 
    template <typename TInternalNumber>
    static
    Point
    orthogonalVector( const Points& basis )
    {
      Point w;
      typedef  PointVector< dimension, TInternalNumber > InternalPoint;
      InternalPoint W;
      const std::size_t n = dimension;
      if ( ( basis.size() + 1 ) != dimension ) return w;
      const std::size_t m = basis.size();
      SimpleMatrix< Scalar, dimension-1, dimension> A;
      for ( std::size_t i = 0; i < m; ++i )
        for ( std::size_t j = 0; j < n; ++j )
          A( i, j ) = basis[ i ][ j ];
      for ( std::size_t col = 0; col < n; ++col)
        { // construct sub-matrix removing column col
          SimpleMatrix< Scalar, dimension-1, dimension-1> M;
          for ( std::size_t i = 0; i < m; ++i )
            {
              std::size_t c = 0;
              for ( std::size_t j = 0; j < n; ++j)
                if ( j != col ) M( i, c++ ) = A( i, j );
            }
          functions::getDeterminantBareiss( W[ col ], M );
          if ( (col+dimension) % 2 == 0 ) W[ col ] = -W[ col ];
        }
      W = AffineGeometry<InternalPoint>::simplifiedVector( W );
      return transform( W );
    }
    
    static 
    Point simplifiedVector( Point v )
    {
      Scalar x;
      functions::getSquaredNormL2( x, v );
      PointOps::normalizeVector( v, x );
      return v;
    }
    
    
  };
 
  namespace functions {
    
    template <typename TPoint>
    DGtal::int64_t
    computeAffineDimension( const std::vector< TPoint >& X, const double tolerance = 1e-12 )
    {
      return AffineGeometry<TPoint>::affineDimension( X, tolerance );
    }
 
    template <typename TPoint>
    std::vector< std::size_t >
    computeAffineSubset( const std::vector< TPoint >& X, const double tolerance = 1e-12 )
    {
      return AffineGeometry<TPoint>::affineSubset( X, tolerance );
    }
 
    template <typename TPoint, typename TIndexRange>
    static
    std::vector< std::size_t >
    computeAffineSubset( const std::vector<TPoint>& X,
                         const TIndexRange& I,
                         const double tolerance = 1e-12 )
    {
      return AffineGeometry<TPoint>::affineSubset( X, I, tolerance );
    }
 
    
    template <typename TPoint, typename TInputPoint>
    void
    getAffineBasis( TInputPoint& o,
                    std::vector< TPoint >& basis,
                    const std::vector< TInputPoint >& X,
                    const double tolerance = 1e-12 )
    {
      std::tie( o, basis ) = AffineGeometry<TPoint>::affineBasis( X, tolerance );
    }
 
    template < typename TPoint, typename TInputPoint, typename TIndexRange >
    static
    void
    getAffineBasis( TInputPoint& o,
                    std::vector< TPoint >& basis,
                    const std::vector< TInputPoint >& X,
                    const TIndexRange& I,
                    const double tolerance = 1e-12 )
    {
      std::tie( o, basis ) = AffineGeometry<TPoint>::affineBasis( X, I, tolerance );
    }
    
    template <typename TPoint>
    TPoint
    computeIndependentVector( const std::vector< TPoint >& basis,
                              const double tolerance = 1e-12 )
    {
      return AffineGeometry<TPoint>::independentVector( basis, tolerance );
    }
 
    template <typename TPoint>
    static
    void
    getCompleteBasis( std::vector< TPoint >& basis,
                      bool normal_vector,
                      const double tolerance = 1e-12 )
    {
      AffineGeometry<TPoint>::completeBasis( basis, normal_vector, tolerance );
    }
 
    template < typename TPoint >
    static
    TPoint
    computeOrthogonalVectorToBasis( const std::vector<TPoint>& basis )
    {
      if ( ( basis.size() + 1 ) != TPoint::dimension ) return TPoint::zero;
      return AffineGeometry<TPoint>::orthogonalVector( basis );
    }
 
    template < typename TPoint, typename TInputPoint, typename TIndexRange >
    static
    void
    getOrthogonalVector( TPoint& w,
                         const std::vector< TInputPoint >& X,
                         const TIndexRange& I,
                         const double tolerance = 1e-12 )
    {
      typedef AffineGeometry<TPoint> Affine;
      w = TPoint::zero;
      TInputPoint o;
      auto subset = Affine::affineSubset( X, I, tolerance );
      if ( ( subset.size() ) != TPoint::dimension ) return;
      std::vector<TPoint> basis( TPoint::dimension - 1);
      for ( auto i = 0; i < basis.size(); i++ )
        basis[ i ] = Affine::transform( X[ subset[ i+1 ] ] - X[ subset[ 0 ] ] );
      w = Affine::orthogonalVector( basis );
    }
 
    template < typename TPoint, typename TInputPoint>
    static
    void
    getOrthogonalVector( TPoint& w,
                         const std::vector< TInputPoint >& X,
                         const double tolerance = 1e-12 )
    { 
      typedef AffineGeometry<TPoint> Affine;
      w = TPoint::zero;
      TInputPoint o;
      auto subset = AffineGeometry<TPoint>::affineSubset( X, tolerance );
      if ( ( subset.size() ) != TPoint::dimension ) return;
      std::vector<TPoint> basis( TPoint::dimension - 1);
      for ( auto i = 0; i < basis.size(); i++ )
        basis[ i ] = Affine::transform( X[ subset[ i+1 ] ] - X[ subset[ 0 ] ] );
      w = Affine::orthogonalVector( basis );
    }
 
    template < typename TPoint, typename TInputPoint>
    static
    void
    getOrthogonalLatticeBasis( std::vector< TPoint >& B,
                               const TInputPoint& N, bool shortened = false )
    {
      B = AffineGeometry<TPoint>::orthogonalLatticeBasis( N, shortened );
    }
 
    
    template <typename TPoint>
    TPoint
    computeSimplifiedVector( const TPoint& v )
    {
      return AffineGeometry<TPoint>::simplifiedVector( v );
    }
    
  } // namespace functions
} // namespace DGtal
 
// Includes inline functions.
//                                                                           //
 
#endif // !defined AffineGeometry_h
 
#undef AffineGeometry_RECURSES
#endif // else defined(AffineGeometry_RECURSES)