polyhedralizer.cpp File Reference

#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/viewers/PolyscopeViewer.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"

Include dependency graph for polyhedralizer.cpp:

Go to the source code of this file.

Data Structures
struct	PairSorted2nd< T1, T2 >
	[polyhedralizer-typedefs] More...
struct	Triple< T1, T2, T3 >

Functions
template<typename RealVector, typename ConstIterator>
double	LSF (RealVector &N, ConstIterator itB, ConstIterator itE)
int	main (int argc, char **argv)

Detailed Description

This program is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.

Author

Jacques-Olivier Lachaud (jacqu.nosp@m.es-o.nosp@m.livie.nosp@m.r.la.nosp@m.chaud.nosp@m.@uni.nosp@m.v-sav.nosp@m.oie..nosp@m.fr ) Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France

Date

2014/06/21

An example file named polyhedralizer.

This file is part of the DGtal library.

Definition in file polyhedralizer.cpp.

Function Documentation

LSF()

template<typename RealVector, typename ConstIterator>

double LSF	(	RealVector &	N,
		ConstIterator	itB,
		ConstIterator	itE )

Examples

examples/tutorial-examples/polyhedralizer.cpp, and tutorial-examples/polyhedralizer.cpp.

Definition at line 123 of file polyhedralizer.cpp.

{
  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;
}

References DGtal::PointVector< dim, TEuclideanRing, TContainer >::dot(), DGtal::EigenDecomposition< TN, TComponent, TMatrix >::getEigenDecomposition(), and DGtal::PointVector< dim, double >::zero.

Referenced by main().

main()

int main	(	int	argc,
		char **	argv )

[polyhedralizer-readVol]

[polyhedralizer-readVol]

[polyhedralizer-KSpace]

[polyhedralizer-KSpace]

[polyhedralizer-SurfelAdjacency]

[polyhedralizer-SurfelAdjacency]

[polyhedralizer-ExtractingSurface]

[polyhedralizer-ExtractingSurface]

[polyhedralizer-ComputingPlaneSize]

[polyhedralizer-ComputingPlaneSize]

[polyhedralizer-segment]

[polyhedralizer-segment]

[polyhedralizer-lsf]

[polyhedralizer-lsf]

[polyhedralizer-projection]

[polyhedralizer-projection]

[polyhedralizer-MakeMesh]

[polyhedralizer-MakeMesh]

[polyhedralizer-visualization]

[polyhedralizer-visualization]

[polyhedralizer-freeMemory]

[polyhedralizer-freeMemory]

Definition at line 157 of file polyhedralizer.cpp.

{
  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 PolyscopeViewer<Space,KSpace> MyViewer3D;
  MyViewer3D viewer( ks );
  // bool isOK = polyhedron >> "test.off";
  // bool isOK2 = polyhedron >> "test.obj";
  viewer << polyhedron;
  viewer.show();

  for ( vector<RoundPlane*>::iterator
          it = roundPlanes.begin(), itE = roundPlanes.end();
        it != itE; ++it )
    delete *it;
  return true;
  }

References DGtal::DigitalSurface< TDigitalSurfaceContainer >::allClosedFaces(), DGtal::COBANaivePlaneComputer< TSpace, TInternalInteger >::begin(), DGtal::DigitalSurface< TDigitalSurfaceContainer >::begin(), DGtal::BreadthFirstVisitor< TGraph, TMarkSet >::current(), DGtal::PointVector< dim, TEuclideanRing, TContainer >::dot(), DGtal::COBANaivePlaneComputer< TSpace, TInternalInteger >::end(), DGtal::DigitalSurface< TDigitalSurfaceContainer >::end(), DGtal::BreadthFirstVisitor< TGraph, TMarkSet >::expand(), DGtal::COBANaivePlaneComputer< TSpace, TInternalInteger >::extend(), DGtal::Surfaces< TKSpace >::findABel(), DGtal::BreadthFirstVisitor< TGraph, TMarkSet >::finished(), DGtal::BreadthFirstVisitor< TGraph, TMarkSet >::ignore(), image(), DGtal::VolReader< TImageContainer, TFunctor >::importVol(), index(), DGtal::COBANaivePlaneComputer< TSpace, TInternalInteger >::init(), DGtal::KhalimskySpaceND< dim, TInteger >::init(), LSF(), DGtal::KhalimskySpaceND< dim, TInteger >::sCoords(), DGtal::KhalimskySpaceND< dim, TInteger >::sDirectIncident(), DGtal::DigitalSurface< TDigitalSurfaceContainer >::size(), DGtal::PointVector< dim, TEuclideanRing, TContainer >::size(), DGtal::KhalimskySpaceND< dim, TInteger >::sKCoords(), DGtal::KhalimskySpaceND< dim, TInteger >::sOrthDir(), DGtal::trace, DGtal::DigitalSurface< TDigitalSurfaceContainer >::verticesAroundFace(), DGtal::DigitalSurface< TDigitalSurfaceContainer >::writeNeighbors(), DGtal::PointVector< 3, double >::zero, and DGtal::PointVector< dim, double >::zero.