cubical-complex-collapse.cpp

Go to the documentation of this file.

 

 

#include <iostream>
#include <map>
#include "DGtal/base/Common.h"
#include "DGtal/helpers/StdDefs.h"
 
#include "DGtal/io/viewers/PolyscopeViewer.h"
#include "DGtal/topology/KhalimskySpaceND.h"
#include "DGtal/topology/CubicalComplex.h"

 
using namespace std;
using namespace DGtal;
using namespace DGtal::Z3i;
 

template <typename CubicalComplex>
struct DiagonalPriority {
  typedef typename CubicalComplex::KSpace KSpace;
  typedef typename CubicalComplex::Point  Point;
  typedef typename CubicalComplex::Cell   Cell;
  typedef typename CubicalComplex::CellMapIterator CellMapIterator;
  
  DiagonalPriority( const CubicalComplex& complex ) : myComplex( complex ) {}
  bool operator()( const CellMapIterator& it1, const CellMapIterator& it2 ) const
  {
    Point k1 = myComplex.space().uKCoords( it1->first );
    Point k2 = myComplex.space().uKCoords( it2->first );
    double d1 = Point::diagonal( 1 ).dot( k1 ) / sqrt( (double) KSpace::dimension );
    double d2 = Point::diagonal( 1 ).dot( k2 ) / sqrt( (double) KSpace::dimension );;
    RealPoint v1( k1[ 0 ] - d1 * k1[ 0 ], k1[ 1 ] - d1 * k1[ 1 ], k1[ 2 ] - d1 * k1[ 2 ] );
    RealPoint v2( k2[ 0 ] - d2 * k2[ 0 ], k2[ 1 ] - d2 * k2[ 1 ], k2[ 2 ] - d2 * k2[ 2 ] );
    double n1 = v1.dot( v1 );
    double n2 = v2.dot( v2 );
    return ( n1 < n2 ) || ( ( n1 == n2 ) && ( it1->first < it2->first ) );
  }
  bool operator()( const CellMapIterator& it1, const CellMapIterator& it2 ) const {…}
  
  const CubicalComplex& myComplex;
};
struct DiagonalPriority {…};

 
int main( int argc, char** argv )
{
  // JOL: unordered_map is approximately twice faster than map for
  // collapsing.
  typedef std::map<Cell, CubicalCellData> Map;
  // typedef boost::unordered_map<Cell, CubicalCellData>   Map;
  typedef CubicalComplex< KSpace, Map >     CC;
 
  trace.beginBlock( "Creating Cubical Complex" );
  KSpace K;
  K.init( Point( 0,0,0 ), Point( 512,512,512 ), true );
  CC complex( K );
  Integer m = 40;
  std::vector<Cell> S;
  for ( Integer x = 0; x <= m; ++x )
    for ( Integer y = 0; y <= m; ++y )
      for ( Integer z = 0; z <= m; ++z )
        {
          Point k1 = Point( x, y, z ); 
          S.push_back( K.uCell( k1 ) );
          double d1 = Point::diagonal( 1 ).dot( k1 ) / (double) KSpace::dimension; // sqrt( (double) KSpace::dimension );
          RealPoint v1( k1[ 0 ], k1[ 1 ], k1[ 2 ] );
          v1 -= d1 * RealPoint::diagonal( 1.0 );
          //RealPoint v1( k1[ 0 ] - d1 * k1[ 0 ], k1[ 1 ] - d1 * k1[ 1 ], k1[ 2 ] - d1 * k1[ 2 ] );
          double n1 = v1.norm();
          bool fixed = ( ( x == 0 ) && ( y == 0 ) && ( z == 0 ) )
            || ( ( x == 0 ) && ( y == m ) && ( z == 0 ) )
            || ( ( x == m ) && ( y == 0 ) && ( z == 0 ) )
            || ( ( x == m ) && ( y == m ) && ( z == 0 ) )
            || ( ( x == m/3 ) && ( y == 2*m/3 ) && ( z == 2*m/3 ) )
            || ( ( x == 0 ) && ( y == 0 ) && ( z == m ) )
            || ( ( x == 0 ) && ( y == m ) && ( z == m ) )
            || ( ( x == m ) && ( y == 0 ) && ( z == m ) )
            || ( ( x == m ) && ( y == m ) && ( z == m ) )
            || ( ( x == 0 ) && ( y == m ) )
            || ( ( x == m ) && ( y == m ) )
            || ( ( z == 0 ) && ( y == m ) )
            || ( ( z == m ) && ( y == m ) );
          complex.insertCell( S.back(), 
                              fixed ? CC::FIXED 
                              : (DGtal::uint32_t) floor(64.0 * n1 ) // This is the priority for collapse 
                              );
        }
  //complex.close();
  trace.info() << "After close: " << complex << std::endl;
  trace.endBlock();
 
  // for 3D display with PolyscopeViewer
  typedef PolyscopeViewer<Space, KSpace> MyViewer;
 
  {
    MyViewer viewer(K);
    typedef CC::CellMapConstIterator CellMapConstIterator;
    for ( Dimension d = 0; d <= 3; ++d )
      for ( CellMapConstIterator it = complex.begin( d ), itE = complex.end( d );
            it != itE; ++it )
        {
          bool fixed = (it->second.data == CC::FIXED);
          if ( fixed ) viewer.drawColor( Color::Red );
          else         viewer.drawColor( Color::White );
          viewer << it->first;
        }
    viewer.show();
  }
  
  trace.beginBlock( "Collapsing complex" );
  CC::DefaultCellMapIteratorPriority P;
  DGtal::uint64_t removed 
    = functions::collapse( complex, S.begin(), S.end(), P, true, true, true );
  trace.info() << "Collapse removed " << removed << " cells." << std::endl;
  trace.info() << "After collapse: " << complex << std::endl;
  trace.endBlock();
 
  {
    MyViewer viewer(K);
    typedef CC::CellMapConstIterator CellMapConstIterator;
    for ( Dimension d = 0; d <= 3; ++d )
      for ( CellMapConstIterator it = complex.begin( d ), itE = complex.end( d );
            it != itE; ++it )
        {
          bool fixed = (it->second.data == CC::FIXED);
          if ( fixed ) viewer.drawColor( Color::Red );
          else         viewer.drawColor( Color::White );
          viewer << it->first;
        }
    viewer.show();
    return 0;
  }
}
int main( int argc, char** argv ) {…}