DGtal 1.3.0
Loading...
Searching...
No Matches
fullConvexityAnalysis3D.cpp
Go to the documentation of this file.
1
71#include <iostream>
72#include <queue>
73#include "DGtal/base/Common.h"
74#include "DGtal/io/viewers/Viewer3D.h"
75#include "DGtal/io/DrawWithDisplay3DModifier.h"
76#include "DGtal/io/Color.h"
77#include "DGtal/shapes/Shapes.h"
78#include "DGtal/helpers/StdDefs.h"
79#include "DGtal/helpers/Shortcuts.h"
80#include "DGtal/images/ImageContainerBySTLVector.h"
81#include "DGtal/geometry/volumes/NeighborhoodConvexityAnalyzer.h"
82
84
85using namespace std;
86using namespace DGtal;
87using namespace Z3i;
89
91
92
93template < typename KSpace, int N >
94struct Analyzer {
95 typedef typename KSpace::Point Point;
97
98 template < typename ImagePtr >
99 static
100 std::vector<Point>
101 debug_one( const KSpace& aK, Point p, ImagePtr bimage )
102 {
103 NCA nca( aK.lowerBound(), aK.upperBound(),
104 KSpace::dimension <= 2 ? 0 : 10000*KSpace::dimension*N );
105 auto& image = *bimage;
106 int geom = 0;
107 nca.setCenter( p, image );
108 bool cvx = nca.isFullyConvex( true );
109 bool ccvx = nca.isComplementaryFullyConvex( false );
110 auto cfg = nca.makeConfiguration( nca.configuration(), true, false );
111 std::vector< Point > localCompX;
112 nca.getLocalCompX( localCompX, false );
113 std::cout << "InC=" << nca.configuration() << std::endl;
114 std::cout << "Cfg=" << cfg << std::endl;
115 for ( auto q : localCompX ) std::cout << q;
116 std::cout << std::endl;
117 geom = ( cvx ? 0x1 : 0x0 ) | ( ccvx ? 0x2 : 0x0 );
118 std::cout << "cvx=" << cvx << " ccvx=" << ccvx << std::endl;
119 std::cout << "geom=" << geom << std::endl;
120 return localCompX;
121 }
122
123 template < typename ImagePtr >
124 static
125 std::vector<int>
126 run( const KSpace& aK, std::vector<Point> pts, ImagePtr bimage )
127 {
128 NCA nca( aK.lowerBound(), aK.upperBound(), 0 );
129 // KSpace::dimension <= 2 ? 0 : 10000*KSpace::dimension*N );
130 auto& image = *bimage;
131 std::vector<int> result;
132 std::map< Point, int > computed;
133 int geom;
134 int i = 0;
135 int nb = pts.size();
136 int nb_cvx = 0;
137 int nb_ccvx = 0;
138 for ( auto p : pts )
139 {
140 if ( i % 100 == 0 ) trace.progressBar( i, nb );
141 auto it = computed.find( p );
142 if ( it == computed.end() )
143 {
144 nca.setCenter( p, image );
145 bool cvx = nca.isFullyConvex( true );
146 bool ccvx = nca.isComplementaryFullyConvex( false );
147 if ( cvx ) nb_cvx += 1;
148 if ( ccvx ) nb_ccvx += 1;
149 geom = ( cvx ? 0x1 : 0x0 ) | ( ccvx ? 0x2 : 0x0 );
150 computed[ p ] = geom;
151 }
152 else geom = it->second;
153 result.push_back( geom );
154 i++;
155 }
156 trace.info() << "nb_cvx=" << nb_cvx << " nb_ccvx=" << nb_ccvx << std::endl;
157 return result;
158 }
159
160 template < typename ImagePtr >
161 static
162 void
163 run( std::vector<int> & to_update,
164 const KSpace& aK, std::vector<Point> pts, ImagePtr bimage )
165 {
166 NCA nca( aK.lowerBound(), aK.upperBound() );
167 // KSpace::dimension <= 2 ? 0 : 10000*KSpace::dimension*N );
168 auto& image = *bimage;
169 std::map< Point, int > computed;
170 int geom;
171 int i = 0;
172 int nb = pts.size();
173 for ( auto p : pts )
174 {
175 if ( i % 100 == 0 ) trace.progressBar( i, nb );
176 auto it = computed.find( p );
177 if ( it == computed.end() )
178 {
179 nca.setCenter( p, image );
180 bool cvx = ( to_update[ i ] & 0x1 )
181 ? nca.isFullyConvex( true )
182 : false;
183 bool ccvx = ( to_update[ i ] & 0x2 )
184 ? nca.isComplementaryFullyConvex( false )
185 : false;
186 geom = ( cvx ? 0x1 : 0x0 ) | ( ccvx ? 0x2 : 0x0 );
187 computed[ p ] = geom;
188 }
189 else geom = it->second;
190 to_update[ i++ ] = geom;
191 }
192 }
193};
194
195template < typename KSpace, int N >
196struct MultiScaleAnalyzer {
197 typedef typename KSpace::Point Point;
199 typedef std::pair< int, int > Geometry; // convex, concave
200
201 template < typename ImagePtr >
202 static
203 std::vector< Geometry >
204 multiscale_run( const KSpace& aK,
205 std::vector<Point> pts,
206 ImagePtr bimage )
207 {
208 auto prev_geometry
209 = MultiScaleAnalyzer< KSpace, N-1>::multiscale_run( aK, pts, bimage );
210 trace.info() << "------- Analyzing scale " << N << " --------" << std::endl;
211 std::vector< int > geom( prev_geometry.size() );
212 for ( int i = 0; i < geom.size(); i++ )
213 geom[ i ] = ( prev_geometry[ i ].first == N-1 ? 0x1 : 0x0 )
214 | ( prev_geometry[ i ].second == N-1 ? 0x2 : 0x0 );
215 Analyzer< KSpace, N>::run( geom, aK, pts, bimage );
216 for ( int i = 0; i < geom.size(); i++ ) {
217 prev_geometry[ i ].first += ( geom[ i ] & 0x1 ) ? 1 : 0;
218 prev_geometry[ i ].second += ( geom[ i ] & 0x2 ) ? 1 : 0;
219 }
220 return prev_geometry;
221 }
222};
223
225template < typename KSpace>
226struct MultiScaleAnalyzer< KSpace, 0 > {
227 typedef typename KSpace::Point Point;
228 typedef std::pair< int, int > Geometry; // convex, concave
229
230 template < typename ImagePtr >
231 static
232 std::vector< Geometry >
233 multiscale_run( const KSpace& aK,
234 std::vector<Point> pts,
235 ImagePtr bimage )
236 {
237 return std::vector< Geometry >( pts.size(), std::make_pair( 0, 0 ) );
238 }
239};
240
241int main( int argc, char** argv )
242{
243 if ( argc <= 2 )
244 {
245 trace.info() << "Usage: " << argv[ 0 ] << " <K> <input.vol> <m> <M>" << std::endl;
246 trace.info() << "\tAnalyze the shape with local full convexity" << std::endl;
247 trace.info() << "\t- 1 <= K <= 5: analysis at scale K" << std::endl;
248 trace.info() << "\t- K == 0: multiscale analysis (using scales 1-5)" << std::endl;
249 trace.info() << "\t- input.vol: choose your favorite shape" << std::endl;
250 trace.info() << "\t- m [==0], M [==255]: used to threshold input vol image" << std::endl;
251 return 1;
252 }
253 int N = atoi( argv[ 1 ] );
254 std::string fn= argv[ 2 ];
255 int m = argc > 3 ? atoi( argv[ 3 ] ) : 0;
256 int M = argc > 4 ? atoi( argv[ 4 ] ) : 255;
257
258 QApplication application(argc,argv);
259
260 auto params = SH3::defaultParameters();
261
262 // Domain creation from two bounding points.
263 trace.info() << "Building set or importing vol ... ";
264 KSpace K;
265 params( "thresholdMin", m );
266 params( "thresholdMax", M );
267 auto bimage = SH3::makeBinaryImage( fn, params );
268 K = SH3::getKSpace( bimage );
269 trace.info() << " [Done]" << std::endl;
270 // Compute surface
271 params( "surfaceComponents" , "All" );
272 auto surface = SH3::makeDigitalSurface( bimage, K, params );
273
274 // Compute interior boundary points
275 // They are less immediate interior points than surfels.
276 std::vector< Point > points;
277 std::map< SCell, int > surfel2idx;
278 std::map< Point, int > point2idx;
279 int idx = 0;
280 for ( auto s : (*surface) )
281 {
282 // get inside point on the border of the shape.
283 Dimension k = K.sOrthDir( s );
284 auto voxel = K.sIncident( s, k, K.sDirect( s, k ) );
285 Point p = K.sCoords( voxel );
286 auto it = point2idx.find( p );
287 if ( it == point2idx.end() )
288 {
289 points.push_back( p );
290 surfel2idx[ s ] = idx;
291 point2idx [ p ] = idx++;
292 }
293 else
294 surfel2idx[ s ] = it->second;
295 }
296 trace.info() << "Shape has " << points.size() << " interior boundary points"
297 << std::endl;
298 if ( N != 0 )
299 {
300 std::vector< int > result;
301 trace.beginBlock ( "Single scale analysis" );
302 if ( N == 1 ) result = Analyzer< KSpace, 1 >::run( K, points, bimage );
303 if ( N == 2 ) result = Analyzer< KSpace, 2 >::run( K, points, bimage );
304 if ( N == 3 ) result = Analyzer< KSpace, 3 >::run( K, points, bimage );
305 if ( N == 4 ) result = Analyzer< KSpace, 4 >::run( K, points, bimage );
306 if ( N == 5 ) result = Analyzer< KSpace, 5 >::run( K, points, bimage );
307 trace.endBlock();
308 SCell dummy;
309 Color colors[ 4 ] =
310 { Color( 255, 0, 0, 255 ), Color( 0, 255, 0, 255 ),
311 Color( 0, 0, 255, 255 ), Color( 255, 255, 255, 255 ) };
312 auto surfels = SH3::getSurfelRange( surface, params );
313 SH3::Colors all_colors( surfels.size() );
314 for ( int i = 0; i < surfels.size(); i++ )
315 {
316 const auto j = surfel2idx[ surfels[ i ] ];
317 all_colors[ i ] = colors[ result[ j ] ];
318 }
319 SH3::saveOBJ( surface, SH3::RealVectors(), all_colors, "geom-cvx.obj" );
320 Viewer3D<> viewer;
321 viewer.setWindowTitle("fullConvexityAnalysis3D");
322 viewer.show();
323 int i = 0;
324 viewer << SetMode3D( dummy.className(), "Basic" );
325 for ( auto s : (*surface) )
326 {
327 viewer << CustomColors3D( all_colors[ i ], all_colors[ i ] )
328 << s;
329 i++;
330 }
331 viewer<< Viewer3D<>::updateDisplay;
332 application.exec();
333 }
334 else
335 {
336 trace.beginBlock ( "Multiscale analysis" );
337 auto geometry =
338 MultiScaleAnalyzer< KSpace, 5 >::multiscale_run( K, points, bimage );
339 trace.endBlock();
340 Color colors_planar[ 6 ] =
341 { Color( 0, 255, 255, 255),
342 Color( 50, 255, 255, 255), Color( 100, 255, 255, 255),
343 Color( 150, 255, 255, 255), Color( 200, 255, 255, 255 ),
344 Color( 255, 255, 255, 255 ) };
345 Color color_atypical( 255, 0, 0, 255 );
346 Color colors_cvx[ 5 ] =
347 { Color( 0, 255, 0, 255 ), Color( 50, 255, 50, 255 ),
348 Color( 100, 255, 100, 255 ), Color( 150, 255, 150, 255 ),
349 Color( 200, 255, 200, 255 ) };
350 Color colors_ccv[ 5 ] =
351 { Color( 0, 0, 255, 255 ), Color( 50, 50, 255, 255 ),
352 Color( 100, 100, 255, 255 ), Color( 150, 150, 255, 255 ),
353 Color( 200, 200, 255, 255 ) };
354 auto surfels = SH3::getSurfelRange( surface, params );
355 SH3::Colors all_colors( surfels.size() );
356 for ( int i = 0; i < surfels.size(); i++ ) {
357 const auto j = surfel2idx[ surfels[ i ] ];
358 int m0 = std::min( geometry[ j ].first, geometry[ j ].second );
359 int m1 = std::max( geometry[ j ].first, geometry[ j ].second );
360 if ( m1 == 0 ) all_colors[ i ] = color_atypical;
361 else if ( m0 == m1 ) all_colors[ i ] = colors_planar[ 5 ];
362 else if ( geometry[ j ].first > geometry[ j ].second )
363 all_colors[ i ] = colors_cvx[ 5 - abs( m0 - m1 ) ];
364 else
365 all_colors[ i ] = colors_ccv[ 5 - abs( m0 - m1 ) ];
366 }
367 SH3::saveOBJ( surface, SH3::RealVectors(), all_colors, "geom-scale-cvx.obj" );
368 SCell dummy;
369 int i = 0;
370 Viewer3D<> viewer;
371 viewer.setWindowTitle("fullConvexityAnalysis3D");
372 viewer.show();
373 viewer << SetMode3D( dummy.className(), "Basic" );
374 for ( auto s : (*surface) )
375 {
376 viewer << CustomColors3D( all_colors[ i ], all_colors[ i ] )
377 << s;
378 i++;
379 }
380 viewer<< Viewer3D<>::updateDisplay;
381 application.exec();
382 }
383 return 0;
384}
385// //
387
Structure representing an RGB triple with alpha component.
Definition: Color.h:68
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
const Point & lowerBound() const
Return the lower bound for digital points 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...
Point sCoords(const SCell &c) const
Return its digital coordinates.
bool sDirect(const SCell &p, Dimension k) const
Return 'true' if the direct orientation of [p] along [k] is in the positive coordinate direction.
const Point & upperBound() const
Return the upper bound for digital points in this space.
SCell sIncident(const SCell &c, Dimension k, bool up) const
Return the forward or backward signed cell incident to [c] along axis [k], depending on [up].
static const constexpr Dimension dimension
Aim: A class that models a neighborhood and that provides services to analyse the convexity properti...
static Configuration makeConfiguration(Configuration current, bool complement, bool with_center)
void setCenter(Point c, const PointPredicate &X)
void getLocalCompX(std::vector< Point > &localCompX, bool with_center) const
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:105
static KSpace getKSpace(const Point &low, const Point &up, Parameters params=parametersKSpace())
Definition: Shortcuts.h:332
std::vector< Color > Colors
Definition: Shortcuts.h:192
static SurfelRange getSurfelRange(CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer > > surface, const Parameters &params=parametersDigitalSurface())
Definition: Shortcuts.h:1547
static CountedPtr< DigitalSurface > makeDigitalSurface(CountedPtr< TPointPredicate > bimage, const KSpace &K, const Parameters &params=parametersDigitalSurface())
Definition: Shortcuts.h:1209
std::vector< RealVector > RealVectors
Definition: Shortcuts.h:179
static Parameters defaultParameters()
Definition: Shortcuts.h:203
static bool saveOBJ(CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer > > digsurf, const TCellEmbedder &embedder, const RealVectors &normals, const Colors &diffuse_colors, std::string objfile, const Color &ambient_color=Color(32, 32, 32), const Color &diffuse_color=Color(200, 200, 255), const Color &specular_color=Color::White)
Definition: Shortcuts.h:1739
static CountedPtr< BinaryImage > makeBinaryImage(Domain shapeDomain)
Definition: Shortcuts.h:561
void beginBlock(const std::string &keyword="")
std::ostream & info()
void progressBar(const double currentValue, const double maximalValue)
double endBlock()
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
Shortcuts< KSpace > SH3
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition: Common.h:137
Trace trace
Definition: Common.h:154
STL namespace.
Modifier class in a Display3D stream. Useful to choose your own mode for a given class....
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
std::string className() const
Return the style name used for drawing this object.
int main()
Definition: testBits.cpp:56
KSpace K