topology/homotopicThinning3D.cpp

An homotopic thinning is an iterative removal of simple points from a given digital object.

See also

Simple points

Resulting 3d thinning with the 18_6 object

Resulting 3d thinning with the 6_26 object

 

 

#include <iostream>
#include <queue>
#include "DGtal/base/Common.h"
#include "DGtal/io/viewers/PolyscopeViewer.h"
#include "DGtal/io/Color.h"
#include "DGtal/shapes/Shapes.h"
#include "DGtal/helpers/StdDefs.h"

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

 
 
 
int main( int argc, char** argv )
{
 
  trace.beginBlock ( "Example simple example of 3DViewer" );
  
  PolyscopeViewer<> viewer;
  
  // Domain cretation from two bounding points.
  Point c( 0, 0, 0 );
  Point p1( -50, -50, -50 );
  Point p2( 50, 50, 50 );
  Domain domain( p1, p2 );
  
  trace.warning() << "Constructing a ring DigitalSet  ... ";
  DigitalSet shape_set( domain );
  for (Domain::ConstIterator it = domain.begin(); it != domain.end(); ++it )
    {
      if ( ((*it - c ).norm() <= 25) && ((*it - c ).norm() >= 18)
     && ( (((*it)[0] <= 3)&& ((*it)[0] >= -3))|| (((*it)[1] <= 3)&& ((*it)[1] >= -3)))){
  shape_set.insertNew( *it );
      }
    }
  trace.warning() << "  [Done]" << std::endl;
  
  trace.beginBlock ( "Thinning" );
  Object18_6 shape( dt18_6,  shape_set );
  int nb_simple=0; 
  DigitalSet::Iterator it, itE;
  do 
    {
      DigitalSet & S = shape.pointSet();
      std::queue<DigitalSet::Iterator> Q;
      it = S.begin(); 
      itE = S.end();
#ifdef DGTAL_WITH_OPENMP
      std::vector<DigitalSet::Iterator> v( S.size() );
      std::vector<uint8_t> b( v.size() );
      for ( size_t i = 0; it != itE; ++it, ++i )
        v[ i ] = it;
#pragma omp parallel for schedule(dynamic)
      for ( size_t i = 0; i < v.size(); ++i )
        b[ i ] = shape.isSimple( *(v[ i ]) );
 
      for ( size_t i = 0; i < v.size(); ++i )
        if ( b[ i ] ) Q.push( v[ i ] ); 
#else
      for ( ; it != itE; ++it )
        if ( shape.isSimple( *it ) )
          Q.push( it );
#endif
      nb_simple = 0;
      while ( ! Q.empty() )
        {
          DigitalSet::Iterator itt = Q.front();
          Q.pop();
          if ( shape.isSimple( *itt ) )
            {
              S.erase( *itt );
              ++nb_simple;
            }
        }
    }
  while ( nb_simple != 0 );
  DigitalSet & S = shape.pointSet();
  trace.endBlock();
 
  // Display by using two different list to manage OpenGL transparency.
 
  viewer << Color(25,25,255, 255);
  viewer << S ; 
 
  viewer << Color(250, 0,0, 25);
  viewer << shape_set;
 
  
  trace.endBlock();
  viewer.show();  
  return 0;
 
}
//                                                                           //