tutorial-examples/polyhedralizer.cpp

Example of tutorial 2: making a polyhedron from a digital object

See also

Tutorial: making a polyhedron from a digital object

The polyhedral surface approaching Al Capone digital object, for width=3/1.

 
#include <iostream>
#include <vector>
#include <set>
#include <map>
#include <queue>
 
#include "DGtal/base/Common.h"
#include "DGtal/helpers/StdDefs.h"
#include "ConfigExamples.h"
 
#include "DGtal/io/readers/VolReader.h"
#include "DGtal/images/ImageSelector.h"
#include "DGtal/images/ImageContainerBySTLVector.h"
#include "DGtal/images/SimpleThresholdForegroundPredicate.h"
 
#include "DGtal/io/Display3D.h"
#include "DGtal/io/viewers/Viewer3D.h"
#include "DGtal/io/DrawWithDisplay3DModifier.h"
#include "DGtal/images/imagesSetsUtils/SetFromImage.h"
 
#include "DGtal/topology/DigitalSurface.h"
#include "DGtal/topology/helpers/Surfaces.h"
#include "DGtal/topology/ImplicitDigitalSurface.h"
 
#include "DGtal/graph/BreadthFirstVisitor.h"
#include "DGtal/geometry/surfaces/COBANaivePlaneComputer.h"
#include "DGtal/geometry/surfaces/ChordNaivePlaneComputer.h"
#include "DGtal/geometry/surfaces/ChordGenericNaivePlaneComputer.h"
 
#include "DGtal/math/linalg/SimpleMatrix.h"
#include "DGtal/math/linalg/EigenDecomposition.h"
 
 
using namespace std;
using namespace DGtal;
using namespace Z3i;
 
template <typename T1, typename T2>
struct PairSorted2nd
{
  typedef PairSorted2nd<T1,T2> Self;
  inline PairSorted2nd( const T1& t1, const T2& t2 ) : first( t1 ), second( t2 ) {}
  bool operator<( const Self& other ) const
  {
    return second < other.second;
  }
  T1 first;
  T2 second;
};
 
template <typename T1, typename T2, typename T3>
struct Triple
{
  T1 first;
  T2 second;
  T3 third;
  Triple( T1 t1 = T1(), T2 t2 = T2(), T3 t3 = T3() )
    : first( t1 ), second( t2 ), third( t3 )
  {}
};
 
// Least-Square Fit of a plane N.x=mu on points [itB,itE). Returns mu.
template <typename RealVector,
          typename ConstIterator>
double LSF( RealVector& N, ConstIterator itB, ConstIterator itE )
{
  typedef typename RealVector::Component Component;
  typedef SimpleMatrix<Component,3,3> Matrix;
  Matrix A; A.clear();
  unsigned int nb = 0;
  RealVector G = RealVector::zero; // center of gravity
  for ( ConstIterator it = itB; it != itE; ++it )
    {
      G += RealVector( (*it)[ 0 ], (*it)[ 1 ], (*it)[ 2 ] );
      ++nb;
    }
  G /= nb;
  for ( ConstIterator it = itB; it != itE; ++it )
    {
      RealVector p( (*it)[ 0 ], (*it)[ 1 ], (*it)[ 2 ] );
      p -= G;
      for ( Dimension i = 0; i < 3; ++i )
        for ( Dimension j = 0; j < 3; ++j )
          A.setComponent( i, j, A( i, j ) + p[ i ] * p[ j ] );
    }
  // A is Gram matrix. We look for V such that V^t A V / |V|^2 is
  // minimal. It is thus the first eigenvalue.
  Matrix V;
  RealVector values;
  EigenDecomposition<3,Component>::getEigenDecomposition( A, V, values );
  N = V.column( 0 ); // first eigenvector;
  double mu = 0.0;
  for ( ConstIterator it = itB; it != itE; ++it )
    mu += N.dot( *it );
  return mu/(double)nb;
}
 
 
int main( int argc, char** argv )
{
  QApplication application(argc,argv);
  string inputFilename = argc > 1 ? argv[ 1 ] : examplesPath+"/samples/Al.100.vol";
  int threshold = argc > 2 ? atoi( argv[ 2 ] ) : 0;
  int widthNum = argc > 3 ? atoi( argv[ 3 ] ) : 2;
  int widthDen = argc > 4 ? atoi( argv[ 4 ] ) : 1;
 
  trace.beginBlock( "Reading vol file into an image." );
  typedef ImageContainerBySTLVector< Domain, int> Image;
  Image image = VolReader<Image>::importVol(inputFilename);
  typedef functors::SimpleThresholdForegroundPredicate<Image> DigitalObject;
  DigitalObject digitalObject( image, threshold );
  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."<<endl;
      return 2;
    }
  trace.endBlock();
 
  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
  MySurfelAdjacency surfAdj( false ); // exterior in all directions.
 
  trace.beginBlock( "Extracting boundary by tracking the surface. " );
  typedef KSpace::Surfel Surfel;
  Surfel start_surfel = Surfaces<KSpace>::findABel( ks, digitalObject, 100000 );
  typedef ImplicitDigitalSurface< KSpace, DigitalObject > MyContainer;
  typedef DigitalSurface< MyContainer > MyDigitalSurface;
  typedef MyDigitalSurface::ConstIterator ConstIterator;
  MyContainer container( ks, digitalObject, surfAdj, start_surfel );
  MyDigitalSurface digSurf( container );
  trace.info() << "Digital surface has " << digSurf.size() << " surfels."
               << endl;
  trace.endBlock();
 
  // First pass to find biggest planes.
  trace.beginBlock( "Decomposition first pass. Computes all planes so as to sort vertices by the plane size." );
  typedef BreadthFirstVisitor<MyDigitalSurface> Visitor;
  typedef ChordGenericNaivePlaneComputer<Z3,Z3::Point, DGtal::int64_t> NaivePlaneComputer;
  map<Surfel,unsigned int> v2size;
  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
    v2size[ *it ] = 0;
  int j = 0;
  int nb = digSurf.size();
  NaivePlaneComputer planeComputer;
  vector<Point> layer;
  vector<Surfel> layer_surfel;
  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
    {
      if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
      Surfel v = *it;
      planeComputer.init( widthNum, widthDen );
      // The visitor takes care of all the breadth-first traversal.
      Visitor visitor( digSurf, v );
      layer.clear();
      layer_surfel.clear();
      Visitor::Size currentSize = visitor.current().second;
      while ( ! visitor.finished() )
        {
          Visitor::Node node = visitor.current();
          v = node.first;
          int axis = ks.sOrthDir( v );
          Point p = ks.sCoords( ks.sDirectIncident( v, axis ) );
          if ( node.second != currentSize )
            {
              bool isExtended = planeComputer.extend( layer.begin(), layer.end() );
              if ( isExtended )
                {
                  for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
                          it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
                    {
                      ++v2size[ *it_layer ];
                    }
                  layer_surfel.clear();
                  layer.clear();
                  currentSize = node.second;
                }
              else
                break;
            }
          layer_surfel.push_back( v );
          layer.push_back( p );
          visitor.expand();
        }
    }
  // Prepare queue
  typedef PairSorted2nd<Surfel,int> SurfelWeight;
  priority_queue<SurfelWeight> Q;
  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
    Q.push( SurfelWeight( *it, v2size[ *it ] ) );
  trace.endBlock();
 
  // Segmentation into planes
  trace.beginBlock( "Decomposition second pass. Visits vertices from the one with biggest plane to the one with smallest plane." );
  typedef Triple<NaivePlaneComputer, Color, pair<RealVector,double> > RoundPlane;
  set<Surfel> processedVertices;
  vector<RoundPlane*> roundPlanes;
  map<Surfel,RoundPlane*> v2plane;
  j = 0;
  while ( ! Q.empty() )
    {
      if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
      Surfel v = Q.top().first;
      Q.pop();
      if ( processedVertices.find( v ) != processedVertices.end() ) // already in set
        continue; // process to next vertex
 
      RoundPlane* ptrRoundPlane = new RoundPlane;
      roundPlanes.push_back( ptrRoundPlane ); // to delete them afterwards.
      v2plane[ v ] = ptrRoundPlane;
      ptrRoundPlane->first.init( widthNum, widthDen );
      ptrRoundPlane->third = make_pair( RealVector::zero, 0.0 );
      // The visitor takes care of all the breadth-first traversal.
      Visitor visitor( digSurf, v );
      layer.clear();
      layer_surfel.clear();
      Visitor::Size currentSize = visitor.current().second;
      while ( ! visitor.finished() )
        {
          Visitor::Node node = visitor.current();
          v = node.first;
          Dimension axis = ks.sOrthDir( v );
          Point p = ks.sCoords( ks.sDirectIncident( v, axis ) );
          if ( node.second != currentSize )
            {
              bool isExtended = ptrRoundPlane->first.extend( layer.begin(), layer.end() );
              if ( isExtended )
                {
                  for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
                          it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
                    {
                      Surfel s = *it_layer;
                      processedVertices.insert( s );
                      if ( v2plane.find( s ) == v2plane.end() )
                        v2plane[ s ] = ptrRoundPlane;
                    }
                  layer.clear();
                  layer_surfel.clear();
                  currentSize = node.second;
                }
              else break;
            }
          layer_surfel.push_back( v );
          layer.push_back( p );
          if ( processedVertices.find( v ) != processedVertices.end() )
            // surfel is already in some plane.
            visitor.ignore();
          else
            visitor.expand();
        }
      if ( visitor.finished() )
        {
          for ( vector<Surfel>::const_iterator it_layer = layer_surfel.begin(),
                  it_layer_end = layer_surfel.end(); it_layer != it_layer_end; ++it_layer )
            {
              Surfel s = *it_layer;
              processedVertices.insert( s );
              if ( v2plane.find( s ) == v2plane.end() )
                v2plane[ s ] = ptrRoundPlane;
            }
        }
      // Assign random color for each plane.
      ptrRoundPlane->second = Color( rand() % 192 + 64, rand() % 192 + 64, rand() % 192 + 64, 255 );
    }
  trace.endBlock();
 
  for ( vector<RoundPlane*>::iterator
          it = roundPlanes.begin(), itE = roundPlanes.end();
        it != itE; ++it )
    {
      NaivePlaneComputer& computer = (*it)->first;
      RealVector normal;
      double mu = LSF( normal, computer.begin(), computer.end() );
      (*it)->third = make_pair( normal, mu );
    }
 
  map<Surfel, RealPoint> coordinates;
  for ( map<Surfel,RoundPlane*>::const_iterator
          it = v2plane.begin(), itE = v2plane.end();
        it != itE; ++it )
    {
      Surfel v = it->first;
      RoundPlane* rplane = it->second;
      Point p = ks.sKCoords( v );
      RealPoint rp( (double)p[ 0 ]/2.0, (double)p[ 1 ]/2.0, (double)p[ 2 ]/2.0 );
      double mu = rplane->third.second;
      RealVector normal = rplane->third.first;
      double lambda = mu - rp.dot( normal );
      coordinates[ v ] = rp + lambda*normal;
    }
  typedef vector<Surfel> SurfelRange;
  map<Surfel, RealPoint> new_coordinates;
  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
    {
      Surfel s = *it;
      SurfelRange neighbors;
      back_insert_iterator<SurfelRange> writeIt = back_inserter( neighbors );
      digSurf.writeNeighbors( writeIt, *it );
      RealPoint x = RealPoint::zero;
      for ( SurfelRange::const_iterator its = neighbors.begin(), itsE = neighbors.end();
            its != itsE; ++its )
        x += coordinates[ *its ];
      new_coordinates[ s ] = x / neighbors.size();
    }
 
  typedef unsigned int Number;
  typedef Mesh<RealPoint> MyMesh;
  typedef MyMesh::MeshFace MeshFace;
  typedef MyDigitalSurface::FaceSet FaceSet;
  typedef MyDigitalSurface::VertexRange VertexRange;
  map<Surfel, Number> index;   // Numbers all vertices.
  Number nbv = 0;
  MyMesh polyhedron( true );
  // Insert all projected surfels as vertices of the polyhedral surface.
  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
    {
      polyhedron.addVertex( new_coordinates[ *it ] );
      index[ *it ] = nbv++;
    }
  // Define faces of the mesh. Outputs closed faces.
  FaceSet faces = digSurf.allClosedFaces();
  for ( typename FaceSet::const_iterator itf = faces.begin(), itf_end = faces.end();
        itf != itf_end; ++itf )
    {
      MeshFace mface( itf->nbVertices );
      VertexRange vtcs = digSurf.verticesAroundFace( *itf );
      int i = 0;
      for ( typename VertexRange::const_iterator itv = vtcs.begin(), itv_end = vtcs.end();
            itv != itv_end; ++itv )
        {
          mface[ i++ ] = index[ *itv ];
        }
      polyhedron.addFace( mface, Color( 255, 243, 150, 255 ) );
    }
 
  typedef Viewer3D<Space,KSpace> MyViewer3D;
  MyViewer3D viewer( ks );
  viewer.show();
  bool isOK = polyhedron >> "test.off";
  bool isOK2 = polyhedron >> "test.obj";
  viewer << polyhedron;
  viewer << MyViewer3D::updateDisplay;
  application.exec();
 
  for ( vector<RoundPlane*>::iterator
          it = roundPlanes.begin(), itE = roundPlanes.end();
        it != itE; ++it )
    delete *it;
 
  if (isOK && isOK2)
    return 0;
  else
    return 1;
}