DGtal  0.9.3beta
greedy-plane-segmentation-ex2.cpp
1 
30 #include <iostream>
32 #include <vector>
33 #include <set>
34 #include <map>
35 #include <queue>
36 #include "DGtal/base/Common.h"
37 #include "DGtal/io/readers/VolReader.h"
38 
39 #include "DGtal/io/Display3D.h"
40 #include "DGtal/io/viewers/Viewer3D.h"
41 #include "DGtal/io/DrawWithDisplay3DModifier.h"
42 #include "DGtal/images/ImageSelector.h"
43 #include "DGtal/images/imagesSetsUtils/SetFromImage.h"
44 #include "DGtal/topology/DigitalSurface.h"
45 #include "DGtal/topology/DigitalSetBoundary.h"
46 #include "DGtal/graph/BreadthFirstVisitor.h"
47 #include "DGtal/geometry/surfaces/COBANaivePlaneComputer.h"
48 #include "DGtal/helpers/StdDefs.h"
49 #include "ConfigExamples.h"
50 
52 
53 using namespace std;
54 using namespace DGtal;
55 
57 using namespace Z3i;
58 typedef DGtal::int64_t InternalInteger;
60 // We choose the DigitalSetBoundary surface container in order to
61 // segment connected or unconnected surfaces.
64 typedef MyDigitalSurface::ConstIterator ConstIterator;
65 typedef MyDigitalSurface::Vertex Vertex;
66 typedef MyDigitalSurface::SurfelSet SurfelSet;
67 typedef SurfelSet::iterator SurfelSetIterator;
73 struct SegmentedPlane {
74  NaivePlaneComputer plane;
75  Color color;
76 };
77 struct VertexSize {
78  Vertex v;
79  unsigned int size;
80  inline VertexSize( const Vertex & aV, unsigned int aSize )
81  : v( aV ), size( aSize )
82  {}
83 };
84 
85 inline
86 bool operator<( const VertexSize & vs1, const VertexSize & vs2 )
87 {
88  return vs1.size < vs2.size;
89 }
91 
93 
94 int main( int argc, char** argv )
95 {
97  trace.info() << "Segments the surface at given threshold within given volume into digital planes of rational width num/den." << std::endl;
98  // Setting default options: ----------------------------------------------
99  // input file used:
100  string inputFilename = examplesPath + "samples/Al.100.vol" ;
101  trace.info() << "input file used " << inputFilename << std::endl;
102  // parameter threshold
103  unsigned int threshold = 0;
104  trace.info() << "the value that defines the isosurface in the image (an integer between 0 and 255)= " << threshold<< std::endl;
105  // parameter widthNum
106  unsigned int widthNum = 1;
107  trace.info() << "the numerator of the rational width (a non-null integer) =" << widthNum<< std::endl;
108  // parameter widthDen
109  unsigned int widthDen = 1;
110  trace.info() << "the denominator of the rational width (a non-null integer)= " << widthDen<< std::endl;
112 
114  QApplication application(argc,argv);
116  Image image = VolReader<Image>::importVol(inputFilename);
117  DigitalSet set3d (image.domain());
118  SetFromImage<DigitalSet>::append<Image>(set3d, image, threshold,255);
120 
122  trace.beginBlock( "Set up digital surface." );
123  // We initializes the cellular grid space used for defining the
124  // digital surface.
125  KSpace ks;
126  bool ok = ks.init( set3d.domain().lowerBound(),
127  set3d.domain().upperBound(), true );
128  if ( ! ok ) std::cerr << "[KSpace.init] Failed." << std::endl;
129  SurfelAdjacency<KSpace::dimension> surfAdj( true ); // interior in all directions.
130  MyDigitalSurfaceContainer* ptrSurfContainer =
131  new MyDigitalSurfaceContainer( ks, set3d, surfAdj );
132  MyDigitalSurface digSurf( ptrSurfContainer ); // acquired
133  trace.endBlock();
135 
137  Point p;
138  Dimension axis;
139  unsigned int j = 0;
140  unsigned int nb = digSurf.size();
141 
142  // First pass to find biggest planes.
143  trace.beginBlock( "1) Segmentation first pass. Computes all planes so as to sort vertices by the plane size." );
144  std::map<Vertex,unsigned int> v2size;
145  NaivePlaneComputer planeComputer;
146  std::priority_queue<VertexSize> Q;
147  for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
148  {
149  if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
150  Vertex v = *it;
151  axis = ks.sOrthDir( v );
152  planeComputer.init( axis, 500, widthNum, widthDen );
153  // The visitor takes care of all the breadth-first traversal.
154  Visitor visitor( digSurf, v );
155  while ( ! visitor.finished() )
156  {
157  Visitor::Node node = visitor.current();
158  v = node.first;
159  axis = ks.sOrthDir( v );
160  p = ks.sCoords( ks.sDirectIncident( v, axis ) );
161  bool isExtended = planeComputer.extend( p );
162  if ( isExtended )
163  // surfel is in plane.
164  visitor.expand();
165  else // surfel is not in plane and should not be used in the visit.
166  visitor.ignore();
167  }
168  Q.push( VertexSize( v, planeComputer.size() ) );
169  }
170  trace.endBlock();
171 
172  trace.beginBlock( "2) Segmentation second pass. Visits vertices from the one with biggest plane to the one with smallest plane." );
173  std::set<Vertex> processedVertices;
174  std::vector<SegmentedPlane*> segmentedPlanes;
175  std::map<Vertex,SegmentedPlane*> v2plane;
176  j = 0;
177  while ( ! Q.empty() )
178  {
179  if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
180  Vertex v = Q.top().v;
181  Q.pop();
182  if ( processedVertices.find( v ) != processedVertices.end() ) // already in set
183  continue; // process to next vertex
184 
185  SegmentedPlane* ptrSegment = new SegmentedPlane;
186  segmentedPlanes.push_back( ptrSegment ); // to delete them afterwards.
187  axis = ks.sOrthDir( v );
188  ptrSegment->plane.init( axis, 500, widthNum, widthDen );
189  // The visitor takes care of all the breadth-first traversal.
190  Visitor visitor( digSurf, v );
191  while ( ! visitor.finished() )
192  {
193  Visitor::Node node = visitor.current();
194  v = node.first;
195  if ( processedVertices.find( v ) == processedVertices.end() )
196  { // Vertex is not in processedVertices
197  axis = ks.sOrthDir( v );
198  p = ks.sCoords( ks.sDirectIncident( v, axis ) );
199  bool isExtended = ptrSegment->plane.extend( p );
200  if ( isExtended )
201  { // surfel is in plane.
202  processedVertices.insert( v );
203  v2plane[ v ] = ptrSegment;
204  visitor.expand();
205  }
206  else // surfel is not in plane and should not be used in the visit.
207  visitor.ignore();
208  }
209  else // surfel is already in some plane.
210  visitor.ignore();
211  }
212  // Assign random color for each plane.
213  ptrSegment->color = Color( rand() % 256, rand() % 256, rand() % 256, 255 );
214  }
215  trace.endBlock();
217 
219  Viewer3D<> viewer( ks );
220  viewer.show();
221  for ( std::map<Vertex,SegmentedPlane*>::const_iterator
222  it = v2plane.begin(), itE = v2plane.end();
223  it != itE; ++it )
224  {
225  viewer << CustomColors3D( it->second->color, it->second->color );
226  viewer << ks.unsigns( it->first );
227  }
228  viewer << Viewer3D<>::updateDisplay;
230 
232  for ( std::vector<SegmentedPlane*>::iterator
233  it = segmentedPlanes.begin(), itE = segmentedPlanes.end();
234  it != itE; ++it )
235  delete *it;
236  segmentedPlanes.clear();
237  v2plane.clear();
239 
240  return application.exec();
241 }
void beginBlock(const std::string &keyword="")
void init(Dimension axis, InternalInteger diameter, InternalInteger widthNumerator=NumberTraits< InternalInteger >::ONE, InternalInteger widthDenominator=NumberTraits< InternalInteger >::ONE)
void progressBar(const double currentValue, const double maximalValue)
Aim: implements association bewteen points lying in a digital domain and values.
Definition: Image.h:69
Trace trace
Definition: Common.h:137
DigitalSurfaceContainer::SurfelConstIterator ConstIterator
DGtal::uint32_t Dimension
Definition: Common.h:120
KSpace::SurfelSet SurfelSet
STL namespace.
double endBlock()
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of a given...
SCell sDirectIncident(const SCell &p, Dimension k) const
Dimension sOrthDir(const SCell &s) const
bool init(const Point &lower, const Point &upper, bool isClosed)
bool extend(const Point &p)
Aim: implements methods to read a "Vol" file format.
Definition: VolReader.h:88
Surfel Vertex
Defines the type for a vertex.
Point sCoords(const SCell &c) const
DGtal is the top-level namespace which contains all DGtal functions and types.
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
Aim: Define utilities to convert a digital set into an image.
Definition: SetFromImage.h:63
static Dimension size()
std::pair< Vertex, Data > Node
FIXME.
std::ostream & info()
Aim: This class is useful to perform a breadth-first exploration of a graph given a starting point or...
Structure representing an RGB triple with alpha component.
Definition: Color.h:66
const Domain & domain() const
Definition: Image.h:192
Cell unsigns(const SCell &p) const
boost::int64_t int64_t
signed 94-bit integer.
Definition: BasicTypes.h:74
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex...