topology/volTrackBoundary.cpp

In many circumstances, it is better to use the presented graph structure of digital surfaces. For instance it may be used to find the surface just by searching it by adjacencies. This process is called tracking. This is done for you by static method Surfaces::trackBoundary.

See also

Constructing digital surfaces by tracking

On the lobser.vol volume, volTrackBoundary.cpp extracts 148364 surfels in 351ms.

* # Commands
* $ ./examples/topology/volTrackBoundary ../examples/samples/lobster.vol 50 255
* 

Digital surface that is the boundary of a (6,18)-connected component in image lobst* er.vol, extracted by tracking from an initial surfel in 351ms.

 

#include <iostream>
#include <queue>
#include "DGtal/base/Common.h"
#include "DGtal/io/viewers/PolyscopeViewer.h"
#include "DGtal/io/readers/VolReader.h"
#include "DGtal/io/Color.h"
#include "DGtal/images/ImageSelector.h"
#include "DGtal/images/imagesSetsUtils/SetFromImage.h"
#include "DGtal/shapes/Shapes.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/topology/helpers/Surfaces.h"

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

 
void usage( int /*argc*/, char** argv )
{
  std::cerr << "Usage: " << argv[ 0 ] << " <fileName.vol> <minT> <maxT>" << std::endl;
  std::cerr << "\t - displays the boundary of the shape stored in vol file <fileName.vol>." << std::endl;
  std::cerr << "\t - voxel v belongs to the shape iff its value I(v) follows minT <= I(v) <= maxT." << std::endl;
}
 
int main( int argc, char** argv )
{
  if ( argc < 4 )
    {
      usage( argc, argv );
      return 1;
    }
  std::string inputFilename = argv[ 1 ];
  unsigned int minThreshold = atoi( argv[ 2 ] );
  unsigned int maxThreshold = atoi( argv[ 3 ] );

  trace.beginBlock( "Reading vol file into an image." );
  typedef ImageSelector < Domain, int>::Type Image;
  Image image = VolReader<Image>::importVol(inputFilename);
  DigitalSet set3d (image.domain());
  SetFromImage<DigitalSet>::append<Image>(set3d, image, 
                                          minThreshold, maxThreshold);
  trace.endBlock();
  
  
  trace.beginBlock( "Construct the Khalimsky space from the image domain." );
  KSpace ks;
  bool space_ok = ks.init( image.domain().lowerBound(), 
                           image.domain().upperBound(), true );
  if (!space_ok)
    {
      trace.error() << "Error in the Khamisky space construction."<<std::endl;
      return 2;
    }
  trace.endBlock();

  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
  MySurfelAdjacency surfAdj( true ); // interior in all directions.

  trace.beginBlock( "Extracting boundary by tracking from an initial bel." );
  KSpace::SCellSet boundary;
  SCell bel = Surfaces<KSpace>::findABel( ks, set3d, 100000 );
  Surfaces<KSpace>::trackBoundary( boundary, ks, 
                                   surfAdj,
                                   set3d, bel );
  trace.endBlock();

  trace.beginBlock( "Displaying surface in PolyscopeViewer." );
  PolyscopeViewer<> viewer( ks );
  viewer << Color( 128, 128, 128 );
  unsigned long nbSurfels = 0;
  for ( KSpace::SCellSet::const_iterator it = boundary.begin(),
          it_end = boundary.end(); it != it_end; ++it, ++nbSurfels )
    viewer << *it;
  trace.info() << "nb surfels = " << nbSurfels << std::endl;
  trace.endBlock();
  viewer.show(); 
  return 0;
}