Loading [MathJax]/extensions/TeX/AMSsymbols.js
DGtal 2.0.0
greedy-plane-segmentation-ex3.cpp
Go to the documentation of this file.
1
16
29
31#include <iostream>
32#include <vector>
33#include <set>
34#include <map>
35#include <queue>
36
37#include "DGtal/base/Common.h"
38#include "DGtal/io/readers/VolReader.h"
39
40#include "DGtal/io/viewers/PolyscopeViewer.h"
41#include "DGtal/images/ImageSelector.h"
42#include "DGtal/images/imagesSetsUtils/SetFromImage.h"
43#include "DGtal/topology/DigitalSurface.h"
44#include "DGtal/topology/DigitalSetBoundary.h"
45#include "DGtal/graph/BreadthFirstVisitor.h"
46#include "DGtal/geometry/surfaces/COBANaivePlaneComputer.h"
47#include "DGtal/helpers/StdDefs.h"
48#include "ConfigExamples.h"
49
51
52using namespace std;
53using namespace DGtal;
54
56using namespace Z3i;
59// We choose the DigitalSetBoundary surface container in order to
60// segment connected or unconnected surfaces.
66typedef SurfelSet::iterator SurfelSetIterator;
72struct SegmentedPlane {
75};
76struct VertexSize {
77 Vertex v;
78 unsigned int size;
79 inline VertexSize( const Vertex & aV, unsigned int aSize )
80 : v( aV ), size( aSize )
81 {}
82};
83
84inline
85bool operator<( const VertexSize & vs1, const VertexSize & vs2 )
86{
87 return vs1.size < vs2.size;
88}
89
90
92
93int main( int argc, char** argv )
94{
96 trace.info() << "Segments the surface at given threshold within given volume into digital planes of rational width num/den." << std::endl;
97 // Setting default options: ----------------------------------------------
98 // input file used:
99 string inputFilename = examplesPath + "samples/Al.100.vol" ;
100 trace.info() << "input file used " << inputFilename << std::endl;
101 // parameter threshold
102 unsigned int threshold = 0;
103 trace.info() << "the value that defines the isosurface in the image (an integer between 0 and 255)= " << threshold<< std::endl;
104 // parameter widthNum
105 unsigned int widthNum = 1;
106 trace.info() << "the numerator of the rational width (a non-null integer) =" << widthNum<< std::endl;
107 // parameter widthDen
108 unsigned int widthDen = 1;
109 trace.info() << "the denominator of the rational width (a non-null integer)= " << widthDen<< std::endl;
111
114 Image image = VolReader<Image>::importVol(inputFilename);
115 DigitalSet set3d (image.domain());
116 SetFromImage<DigitalSet>::append<Image>(set3d, image, threshold,255);
118
120 trace.beginBlock( "Set up digital surface." );
121 // We initializes the cellular grid space used for defining the
122 // digital surface.
123 KSpace ks;
124 bool ok = ks.init( set3d.domain().lowerBound(),
125 set3d.domain().upperBound(), true );
126 if ( ! ok ) std::cerr << "[KSpace.init] Failed." << std::endl;
127 SurfelAdjacency<KSpace::dimension> surfAdj( true ); // interior in all directions.
128 MyDigitalSurfaceContainer* ptrSurfContainer =
129 new MyDigitalSurfaceContainer( ks, set3d, surfAdj );
130 MyDigitalSurface digSurf( ptrSurfContainer ); // acquired
131 trace.endBlock();
133
135 Point p;
136 Dimension axis;
137 unsigned int j = 0;
138 unsigned int nb = digSurf.size();
139
140 // First pass to find biggest planes.
141 trace.beginBlock( "1) Segmentation first pass. Computes all planes so as to sort vertices by the plane size." );
142 std::map<Vertex,unsigned int> v2size;
143 NaivePlaneComputer planeComputer;
144 std::priority_queue<VertexSize> Q;
145 std::vector<Point> layer;
146 for ( ConstIterator it = digSurf.begin(), itE= digSurf.end(); it != itE; ++it )
147 {
148 if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
149 Vertex v = *it;
150 axis = ks.sOrthDir( v );
151 planeComputer.init( axis, 500, widthNum, widthDen );
152 // The visitor takes care of all the breadth-first traversal.
153 Visitor visitor( digSurf, v );
154 layer.clear();
155 Visitor::Size currentSize = visitor.current().second;
156 while ( ! visitor.finished() )
157 {
158 Visitor::Node node = visitor.current();
159 v = node.first;
160 axis = ks.sOrthDir( v );
161 p = ks.sCoords( ks.sDirectIncident( v, axis ) );
162 if ( node.second != currentSize )
163 {
164 // std::cerr << "Layer " << currentSize << ", size=" << layer.size() << std::endl;
165 bool isExtended = planeComputer.extend( layer.begin(), layer.end() );
166 if ( ! isExtended )
167 break;
168 layer.clear();
169 currentSize = node.second;
170 }
171 layer.push_back( p );
172 visitor.expand();
173 }
174 // std::cerr << v << " -> " << planeComputer.size() << std::endl;
175 Q.push( VertexSize( v, planeComputer.size() ) );
176 }
177 trace.endBlock();
178
179 trace.beginBlock( "2) Segmentation second pass. Visits vertices from the one with biggest plane to the one with smallest plane." );
180 std::set<Vertex> processedVertices;
181 std::map<Vertex,SegmentedPlane*> v2plane;
182 std::vector<SegmentedPlane*> segmentedPlanes;
183 j = 0;
184 while ( ! Q.empty() )
185 {
186 if ( ( (++j) % 50 == 0 ) || ( j == nb ) ) trace.progressBar( j, nb );
187 Vertex v = Q.top().v;
188 Q.pop();
189 if ( processedVertices.find( v ) != processedVertices.end() ) // already in set
190 continue; // process to next vertex
191
192 SegmentedPlane* ptrSegment = new SegmentedPlane;
193 segmentedPlanes.push_back( ptrSegment ); // to delete them afterwards.
194 axis = ks.sOrthDir( v );
195 ptrSegment->plane.init( axis, 500, widthNum, widthDen );
196 // The visitor takes care of all the breadth-first traversal.
197 Visitor visitor( digSurf, v );
198 while ( ! visitor.finished() )
199 {
200 Visitor::Node node = visitor.current();
201 v = node.first;
202 if ( processedVertices.find( v ) == processedVertices.end() )
203 { // Vertex is not in processedVertices
204 axis = ks.sOrthDir( v );
205 p = ks.sCoords( ks.sDirectIncident( v, axis ) );
206 bool isExtended = ptrSegment->plane.extend( p );
207 if ( isExtended )
208 { // surfel is in plane.
209 processedVertices.insert( v );
210 v2plane[ v ] = ptrSegment;
211 visitor.expand();
212 }
213 else // surfel is not in plane and should not be used in the visit.
214 visitor.ignore();
215 }
216 else // surfel is already in some plane.
217 visitor.ignore();
218 }
219 // Assign random color for each plane.
220 ptrSegment->color = Color( rand() % 256, rand() % 256, rand() % 256, 255 );
221 }
222 trace.endBlock();
224
226 PolyscopeViewer<> viewer( ks );
227 Color col( 255, 255, 120 );
228 for ( std::map<Vertex,SegmentedPlane*>::const_iterator
229 it = v2plane.begin(), itE = v2plane.end();
230 it != itE; ++it )
231 {
232 viewer << it->second->color;
233 viewer << ks.unsigns( it->first );
234 }
236
238 for ( std::vector<SegmentedPlane*>::iterator
239 it = segmentedPlanes.begin(), itE = segmentedPlanes.end();
240 it != itE; ++it )
241 delete *it;
242 segmentedPlanes.clear();
243 v2plane.clear();
245
246 viewer.show();
247 return 0;
248}
249
Aim: This class is useful to perform a breadth-first exploration of a graph given a starting point or...
const Node & current() const
Aim: A class that contains the COBA algorithm (Emilie Charrier, Lilian Buzer, DGCI2008) for recognizi...
void init(Dimension axis, InternalInteger diameter, InternalInteger widthNumerator=NumberTraits< InternalInteger >::ONE, InternalInteger widthDenominator=NumberTraits< InternalInteger >::ONE)
bool extend(const Point &p)
Structure representing an RGB triple with alpha component.
Definition Color.h:77
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of a given...
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
DigitalSurfaceContainer::SurfelConstIterator ConstIterator
ConstIterator begin() const
ConstIterator end() const
Aim: implements association bewteen points lying in a digital domain and values.
Definition Image.h:70
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
Dimension sOrthDir(const SCell &s) const
Given a signed surfel [s], returns its orthogonal direction (ie, the coordinate where the surfel is c...
Cell unsigns(const SCell &p) const
Creates an unsigned cell from a signed one.
Point sCoords(const SCell &c) const
Return its digital coordinates.
SCell sDirectIncident(const SCell &p, Dimension k) const
Return the direct incident cell of [p] along [k] (the incident cell along [k])
void show() override
Starts the event loop and display of elements.
Aim: Represent adjacencies between surfel elements, telling if it follows an interior to exterior ord...
DigitalSurface< MyDigitalSurfaceContainer > MyDigitalSurface
SurfelSet::iterator SurfelSetIterator
COBANaivePlaneComputer< Z3, InternalInteger > NaivePlaneComputer
MyDigitalSurface::SurfelSet SurfelSet
BreadthFirstVisitor< MyDigitalSurface > Visitor
DigitalSetBoundary< KSpace, DigitalSet > MyDigitalSurfaceContainer
DGtal::int64_t InternalInteger
DigitalSetBoundary< KSpace, DigitalSet > MyDigitalSurfaceContainer
Z3i this namespace gathers the standard of types for 3D imagery.
KhalimskySpaceND< 3, Integer > KSpace
Definition StdDefs.h:146
DigitalSetSelector< Domain, BIG_DS+HIGH_BEL_DS >::Type DigitalSet
Definition StdDefs.h:173
DGtal is the top-level namespace which contains all DGtal functions and types.
std::int64_t int64_t
signed 94-bit integer.
Definition BasicTypes.h:73
bool operator<(PointVector< ptDim, LeftEuclideanRing, LeftContainer > const &lhs, PointVector< ptDim, RightEuclideanRing, RightContainer > const &rhs)
Comparison operator on Points/Vectors (LesserThan).
DGtal::uint32_t Dimension
Definition Common.h:119
Trace trace
STL namespace.
ImageContainerBySTLVector< Domain, Value > Type
static void append(Set &aSet, const ForegroundPredicate &isForeground, typename Image::Domain::ConstIterator itBegin, typename Image::Domain::ConstIterator itEnd)
static ImageContainer importVol(const std::string &filename, const Functor &aFunctor=Functor())
VertexSize(const Vertex &aV, unsigned int aSize)
int main()
Definition testBits.cpp:56
Image image(domain)
TriMesh::Vertex Vertex