AffineBasis.h

 
#pragma once
 
#if defined(AffineBasis_RECURSES)
#error Recursive header files inclusion detected in AffineBasis.h
#else // defined(AffineBasis_RECURSES)
#define AffineBasis_RECURSES
 
#if !defined AffineBasis_h
#define AffineBasis_h
 
// Inclusions
#include <type_traits>
#include <vector>
#include "DGtal/base/Common.h"
#include "DGtal/geometry/tools/AffineGeometry.h"
 
namespace DGtal
{
  // template class AffineBasis
 
  template < typename TPoint >
  struct AffineBasis
  {
    typedef AffineBasis<TPoint>        Self;
    typedef TPoint                     Point;
    typedef typename Point::Coordinate Scalar;
    typedef std::vector< Point >       Points;
    typedef AffineGeometry< Point >    Affine;
 
    enum struct Type {
      INVALID = 0,     
      ECHELON_REDUCED, 
      SHORTEST_ECHELON_REDUCED, 
      LLL_REDUCED      
    };
    
    // ----------------------- standard services --------------------------
  public:
 
    AffineBasis( const double tolerance = 1e-12)
      : epsilon( tolerance )
    {
      first  = Point::zero;
      second.resize( Point::dimension );
      for ( auto k = 0; k < Point::dimension; k++ )
        second[ k ] = Point::base( k );
      _type = Type::SHORTEST_ECHELON_REDUCED;
    }
 
    template <typename TInputPoint>
    AffineBasis( const std::vector< TInputPoint >& points,
                 AffineBasis::Type type,
                 const double delta = 0.99,
                 const double tolerance = 1e-12 )
      : epsilon( tolerance )
    {
      if ( points.size() == 0 ) return;
      first  = Affine::transform( points[ 0 ] );
      std::vector< TInputPoint > basis( points.size() - 1 );
      for ( std::size_t i = 0; i < basis.size(); i++ )
        basis[ i ] = ( points[ i+1 ] - first );
      initBasis( basis );
      reduce( type, delta );
    }
    
    template <typename TInputPoint>
    AffineBasis( const TInputPoint& origin,
                 const std::vector<TInputPoint>& basis,
                 AffineBasis::Type type,
                 bool is_reduced = false,
                 const double delta = 0.99,
                 const double tolerance = 1e-12 )
      : epsilon( tolerance )
    {
      first = Affine::transform( origin );
      initBasis( basis );
      if ( ! is_reduced ) reduce( type, delta );
      else _type = type;
    }
 
    template <typename TInputPoint>
    AffineBasis( const TInputPoint& origin,
                 const TInputPoint& normal,
                 AffineBasis::Type type = Type::ECHELON_REDUCED,
                 const double tolerance = 1e-12 )
      : epsilon( tolerance )
    {
      first = Affine::transform( origin );
      // basis is shortened is type is LLL.
      second = Affine::orthogonalLatticeBasis( normal, type == Type::LLL_REDUCED );
      _type = ( type == Type::LLL_REDUCED )
        ? Type::LLL_REDUCED : Type::ECHELON_REDUCED;
    }
    
    void reduce( AffineBasis::Type type, double delta )
    {
      if ( type == Type::SHORTEST_ECHELON_REDUCED )
        sortBasis();
      if ( type == Type::ECHELON_REDUCED || type == Type::SHORTEST_ECHELON_REDUCED )
        reduceAsEchelon( type );
      else if ( type == Type::LLL_REDUCED )
        {
          std::vector< Point > X( second.size()+1 );
          X[ 0 ] = first;
          for ( auto i = 0; i < second.size(); i++ )
            X[ i+1 ] = second[ i ] + first;
          std::tie( first, second ) = AffineGeometry<Point>::affineBasis( X, epsilon );
          reduceAsLLL( delta, (Scalar) 0 );
        }
    }
    
    Dimension dimension() const
    {
      return second.size();
    }
 
    const Point& origin() const
    {
      return first;
    }
    
    const Points& basis() const
    {
      return second;
    }
 
    // ----------------------- geometry services --------------------------
  public:
 
    bool isParallel( const Self& other ) const
    {
      if ( dimension() != other.dimension() ) return false;
      if ( ( _type != Type::ECHELON_REDUCED )
           && ( _type != Type::SHORTEST_ECHELON_REDUCED ) )
        trace.error() << "[AffineBasis::isParallel] Requires type=*_ECHELON_REDUCED\n"
                      << " type=" << reductionTypeName() << "\n" ;
      for ( const auto& b : other.second )
        if ( ! isParallel( b ) ) return false;
      return true;
    }
 
    std::pair< Scalar, Point > rationalCoordinates( const Point& p ) const
    {
      const auto [d, lambda, remainder] = decompose( p );
      return std::make_pair( d, lambda );
    }
 
    bool isOnAffineSpace( const Point& p ) const
    {
      const auto [d, lambda, r] = decompose( p );
      return ! Affine::ScalarOps::isNonZero( r.normInfinity(), epsilon );
    }
 
    bool isParallel( const Point& w ) const
    {
      const auto [d, lambda, r] = decomposeVector( w );
      return ! Affine::ScalarOps::isNonZero( r.normInfinity(), epsilon );
    }
 
    Point recompose( Scalar d, const Point& lambda,
                     const Point& r = Point::zero ) const
    {
      return first + recomposeVector( d, lambda, r );
    }
 
    Point recomposeVector( Scalar d, const Point& lambda,
                           const Point& r = Point::zero ) const
    {
      Point w = r;
      for ( std::size_t i = 0; i < second.size(); i++ )
        w += lambda[ i ] * second[ i ];
      return w / d;
    }
    
    std::tuple< Scalar, Point, Point > decompose( const Point& p ) const
    {
      return decomposeVector( p - first );
    }
 
    std::tuple< Scalar, Point, Point > decomposeVector( Point w ) const
    {
      Point r;
      Scalar alphas = 1;
      for ( auto i = 0; i < second.size(); i++ )
        {
          std::pair< Scalar, Scalar > c
            = Affine::reduceVector( w, second[ i ], i, epsilon );
          for ( auto j = 0; j < i; j++ )
            r[ j ] *= c.first;
          r[ i ]  = c.second;
          alphas *= c.first;
        }
      return alphas >= 0
        ? std::make_tuple(  alphas,  r,  w )
        : std::make_tuple( -alphas, -r, -w );
    }
 
    template <typename ProjectedPoint>
    Scalar projectPoints( std::vector< ProjectedPoint >& result,
                          const Points& input )
    {
      Scalar lcm = 1;
      std::vector< Scalar > denoms ( input.size() );
      std::vector< Point >  lambdas( input.size() );
      Point  r;
      // get points rational coordinates.
      for ( std::size_t i = 0; i < input.size(); i++ )
        std::tie( denoms[ i ], lambdas[ i ], r ) = decompose( input[ i ] );
      // compute ppcm
      for ( std::size_t i = 0; i < denoms.size(); i++ )
        {
          const Scalar d = denoms[ i ];
          if ( d == 1 ) continue;
          lcm = Affine::ScalarOps::lcmPositive( lcm, d );
        }
      // project points
      result.resize( input.size() );
      for ( std::size_t i = 0; i < input.size(); i++ )
        dilatedTransform( result[ i ], lambdas[ i ], lcm / denoms[ i ] );
      return lcm;
    }
 
    template <typename OtherPoint>
    static
    void transform( OtherPoint& pp, const Point& p )
    {
      BOOST_STATIC_ASSERT( OtherPoint::dimension <= Point::dimension );
      typedef typename OtherPoint::Coordinate Scalar;
      for ( std::size_t i = 0; i < pp.dimension; ++i )
        pp[ i ] = Scalar( p[ i ] );
    }
      
    template <typename OtherPoint>
    static
    void dilatedTransform( OtherPoint& pp, const Point& p, Scalar m )
    {
      BOOST_STATIC_ASSERT( OtherPoint::dimension <= Point::dimension );
      for ( std::size_t i = 0; i < pp.dimension; ++i )
        pp[ i ] = m * Scalar( p[ i ] );
    }
    
 
    // ----------------------- debug and I/O services --------------------------
  public:
 
    void selfDisplay( std::ostream& out ) const
    {
      out << "[ AffineBasis o=" << first
          << " type=" << reductionTypeName()
          << " B=";
      for ( auto b : second ) std::cout << "\n  " << b;
      out << " ]";
    }
 
    bool isValid() const
    {
      return _type != Type::INVALID;
    }  
 
    std::string reductionTypeName() const
    {
      if ( _type == Type::INVALID )              return "INVALID";
      else if ( _type == Type::ECHELON_REDUCED ) return "ECHELON_REDUCED";
      else if ( _type == Type::SHORTEST_ECHELON_REDUCED ) return "SHORTEST_ECHELON_REDUCED";
      else if ( _type == Type::LLL_REDUCED )     return "LLL_REDUCED";
      else return "";
    }
 
    // ----------------------- public data --------------------------
  public:
 
    Point  first;  
    Points second; 
    double epsilon {1e-12};
    AffineBasis::Type _type; 
 
    // ----------------------- protected services --------------------------
  protected:
 
    void normalize()
    {
      for ( auto& v : second )
        v = Affine::simplifiedVector( v );
    }
    
    void reduceAsEchelon( Type type )
    {
      std::vector< bool > is_independent( second.size(), false );
      std::vector< std::vector< Scalar > > U( second.size() );
      Dimension k = 0; // the current column to put in echelon form.
      for ( std::size_t i = 0; i < second.size(); i++ )
        {
          std::size_t row = findIndexWithSmallestNonNullComponent( k, i, second );
          if ( row != i && row != second.size() )
            std::swap( second[ i ], second[ row ] );            
          Point& w            = second[ i ];
          // check if this vector is independent from the previous ones
          is_independent[ i ] = true;
          for ( std::size_t j = 0; j < i; j++ )
            if ( is_independent[ j ] )
              Affine::reduceVector( w, second[ j ], epsilon );
          if ( ! Affine::ScalarOps::isNonZero( w.normInfinity(), epsilon ) )
            is_independent[ i ] = false; // not independent, forget it
          else
            { // independent, make sure it is reduced.
              w = Affine::simplifiedVector( w );
              k++;
            }
        }
      Points new_basis;
      for ( std::size_t i = 0; i < second.size(); i++ )
        if ( is_independent[ i ] )
          new_basis.push_back( second[ i ] );
      std::swap( second, new_basis );
      orderEchelonBasis();
      _type = type;
    }
 
    void orderEchelonBasis()
    {
      auto compare = [this] ( const Point& v, const Point& w ) -> bool
      {
        // Note: curiously std::sort sometimes test v against itself
        // and must return false in this case.
        for ( auto k = 0; k < Point::dimension; ++k )
          {
            bool v_non_null = Affine::ScalarOps::isNonZero( v[ k ], epsilon );
            bool w_non_null = Affine::ScalarOps::isNonZero( w[ k ], epsilon );
            if ( v_non_null && ! w_non_null )      return true;
            else if ( ! v_non_null && w_non_null ) return false;
            else if ( v_non_null && w_non_null )   return false; // ==
          }
        return false; // 0 == 0
      };
      std::sort( second.begin(), second.end(), compare );
    }
    
    void reduceAsLLL( double delta, Scalar )
    {
      if constexpr( std::is_floating_point< Scalar >::value == true )
        {
          trace.error() << "[AffineBasis::reduceAsLLL]"
                        << " It has no meaning to use LLL algorithm on matrix with double coefficients\n";
          _type = Type::INVALID;
          return;
        }
      for ( auto& v : second )
        v = Affine::simplifiedVector( v );
      std::vector< std::vector< Scalar > > B( second.size() );
      for ( auto i = 0; i < second.size(); i++ )
        {
          B[ i ] = std::vector< Scalar >( Point::dimension );
          for ( auto j = 0; j < Point::dimension; j++ )
            B[ i ][ j ] = second[ i ][ j ];
        }
      functions::reduceBasisWithLLL( B, delta );
      second.clear();
      for ( std::size_t i = 0; i < B.size(); i++ )
        {
          Point b;
          for ( auto j = 0; j < Point::dimension; j++ )
            b[ j ] = B[ i ][ j ];
          if ( b != Point::zero )
            second.push_back( Affine::simplifiedVector( b ) );
        }
      _type = Type::LLL_REDUCED;
    }
 
    template <typename TInputPoint>
    void initBasis( const std::vector<TInputPoint>& basis )
    {
      second.reserve( basis.size() );
      for ( auto i = 0; i < basis.size(); i++ )
        {
          Point b = Affine::transform( basis[ i ] );
          if ( b != Point::zero ) second.push_back( b );
        }
    }
 
    void sortBasis()
    {
      // Reduces all vectors
      normalize();
      // Purge duplicates
      std::sort( second.begin(), second.end() );
      second.erase( std::unique( second.begin(), second.end() ), second.end() );
      // Sort according to size of components.
      auto compare = []( const Point& u, const Point& v ) -> bool
      {
        const auto n1_u = u.norm1();
        const auto n1_v = v.norm1();
        if ( n1_u < n1_v ) return true;
        else if ( n1_v < n1_u ) return false;
        const auto noo_u = u.normInfinity();
        const auto noo_v = v.normInfinity();
        if ( noo_u < noo_v ) return true;
        else if ( noo_v < noo_u ) return false;
        return u < v;
      };
      std::sort( second.begin(), second.end(), compare );
    }
 
    std::size_t
    findIndexWithSmallestNonNullComponent( Dimension k,
                                           std::size_t i,
                                           const std::vector< Point >& basis )
    {
      ASSERT( ! basis.empty() );
      ASSERT( k < Point::dimension );
      std::size_t index = i;
      Scalar      v     = 0;
      for ( ; index < basis.size(); index++ )
        {
          v = abs( basis[ index ][ k ] );
          if ( Affine::ScalarOps::isNonZero( v, epsilon ) )
            break;
        }
      for ( auto j = index + 1; j < basis.size(); j++ )
        {
          Scalar vj = abs( basis[ j ][ k ] );
          if ( vj != 0 && vj < v )
            {
              index = j;
              v     = vj;
            }
        }
      return index;
    }
    
  }; // struct AffineBasis
    
} // namespace DGtal
 
// Includes inline functions.
//                                                                           //
 
template <typename TPoint>
std::ostream&
operator<<( std::ostream& out, const DGtal::AffineBasis<TPoint>& B )
{
  B.selfDisplay( out );
  return out;
}
 
 
#endif // !defined AffineBasis_h
 
#undef AffineBasis_RECURSES
#endif // else defined(AffineBasis_RECURSES)