Loading [MathJax]/extensions/TeX/AMSsymbols.js
DGtal 2.0.0
fullConvexityAnalysis3D.cpp
Go to the documentation of this file.
1
16
29
30
68
69
71#include <iostream>
72#include <queue>
73#include "DGtal/base/Common.h"
74#include "DGtal/io/viewers/PolyscopeViewer.h"
75#include "DGtal/io/Color.h"
76#include "DGtal/shapes/Shapes.h"
77#include "DGtal/helpers/StdDefs.h"
78#include "DGtal/helpers/Shortcuts.h"
79#include "DGtal/images/ImageContainerBySTLVector.h"
80#include "DGtal/geometry/volumes/NeighborhoodConvexityAnalyzer.h"
81
83
84using namespace std;
85using namespace DGtal;
86using namespace Z3i;
88
90
91
92template < typename KSpace, int N >
93struct Analyzer {
94 typedef typename KSpace::Point Point;
96
97 template < typename ImagePtr >
98 static
99 std::vector<Point>
100 debug_one( const KSpace& aK, Point p, ImagePtr bimage )
101 {
102 NCA nca( aK.lowerBound(), aK.upperBound(),
103 KSpace::dimension <= 2 ? 0 : 10000*KSpace::dimension*N );
104 auto& image = *bimage;
105 int geom = 0;
106 nca.setCenter( p, image );
107 bool cvx = nca.isFullyConvex( true );
108 bool ccvx = nca.isComplementaryFullyConvex( false );
109 auto cfg = nca.makeConfiguration( nca.configuration(), true, false );
110 std::vector< Point > localCompX;
111 nca.getLocalCompX( localCompX, false );
112 std::cout << "InC=" << nca.configuration() << std::endl;
113 std::cout << "Cfg=" << cfg << std::endl;
114 for ( auto q : localCompX ) std::cout << q;
115 std::cout << std::endl;
116 geom = ( cvx ? 0x1 : 0x0 ) | ( ccvx ? 0x2 : 0x0 );
117 std::cout << "cvx=" << cvx << " ccvx=" << ccvx << std::endl;
118 std::cout << "geom=" << geom << std::endl;
119 return localCompX;
120 }
121
122 template < typename ImagePtr >
123 static
124 std::vector<int>
125 run( const KSpace& aK, const std::vector<Point>& pts, ImagePtr bimage )
126 {
127 NCA nca( aK.lowerBound(), aK.upperBound(), 0 );
128 // KSpace::dimension <= 2 ? 0 : 10000*KSpace::dimension*N );
129 auto& image = *bimage;
130 std::vector<int> result;
131 std::map< Point, int > computed;
132 int geom;
133 int i = 0;
134 int nb = pts.size();
135 int nb_cvx = 0;
136 int nb_ccvx = 0;
137 for ( auto p : pts )
138 {
139 if ( i % 100 == 0 ) trace.progressBar( i, nb );
140 auto it = computed.find( p );
141 if ( it == computed.end() )
142 {
143 nca.setCenter( p, image );
144 bool cvx = nca.isFullyConvex( true );
145 bool ccvx = nca.isComplementaryFullyConvex( false );
146 if ( cvx ) nb_cvx += 1;
147 if ( ccvx ) nb_ccvx += 1;
148 geom = ( cvx ? 0x1 : 0x0 ) | ( ccvx ? 0x2 : 0x0 );
149 computed[ p ] = geom;
150 }
151 else geom = it->second;
152 result.push_back( geom );
153 i++;
154 }
155 trace.info() << "nb_cvx=" << nb_cvx << " nb_ccvx=" << nb_ccvx << std::endl;
156 return result;
157 }
158
159 template < typename ImagePtr >
160 static
161 void
162 run( std::vector<int> & to_update,
163 const KSpace& aK, const std::vector<Point>& pts, ImagePtr bimage )
164 {
165 NCA nca( aK.lowerBound(), aK.upperBound() );
166 // KSpace::dimension <= 2 ? 0 : 10000*KSpace::dimension*N );
167 auto& image = *bimage;
168 std::map< Point, int > computed;
169 int geom;
170 int i = 0;
171 int nb = pts.size();
172 for ( auto p : pts )
173 {
174 if ( i % 100 == 0 ) trace.progressBar( i, nb );
175 auto it = computed.find( p );
176 if ( it == computed.end() )
177 {
178 nca.setCenter( p, image );
179 bool cvx = ( to_update[ i ] & 0x1 )
180 ? nca.isFullyConvex( true )
181 : false;
182 bool ccvx = ( to_update[ i ] & 0x2 )
183 ? nca.isComplementaryFullyConvex( false )
184 : false;
185 geom = ( cvx ? 0x1 : 0x0 ) | ( ccvx ? 0x2 : 0x0 );
186 computed[ p ] = geom;
187 }
188 else geom = it->second;
189 to_update[ i++ ] = geom;
190 }
191 }
192};
193
194template < typename KSpace, int N >
196 typedef typename KSpace::Point Point;
198 typedef std::pair< int, int > Geometry; // convex, concave
199
200 template < typename ImagePtr >
201 static
202 std::vector< Geometry >
204 const std::vector<Point>& pts,
205 ImagePtr bimage )
206 {
207 auto prev_geometry
209 trace.info() << "------- Analyzing scale " << N << " --------" << std::endl;
210 std::vector< int > geom( prev_geometry.size() );
211 for ( int i = 0; i < geom.size(); i++ )
212 geom[ i ] = ( prev_geometry[ i ].first == N-1 ? 0x1 : 0x0 )
213 | ( prev_geometry[ i ].second == N-1 ? 0x2 : 0x0 );
214 Analyzer< KSpace, N>::run( geom, aK, pts, bimage );
215 for ( int i = 0; i < geom.size(); i++ ) {
216 prev_geometry[ i ].first += ( geom[ i ] & 0x1 ) ? 1 : 0;
217 prev_geometry[ i ].second += ( geom[ i ] & 0x2 ) ? 1 : 0;
218 }
219 return prev_geometry;
220 }
221};
222
224template < typename KSpace>
226 typedef typename KSpace::Point Point;
227 typedef std::pair< int, int > Geometry; // convex, concave
228
229 template < typename ImagePtr >
230 static
231 std::vector< Geometry >
233 const std::vector<Point>& pts,
234 ImagePtr bimage )
235 {
236 ((void) aK);
237 ((void) bimage);
238 return std::vector< Geometry >( pts.size(), std::make_pair( 0, 0 ) );
239 }
240};
241
242int main( int argc, char** argv )
243{
244 if ( argc <= 2 )
245 {
246 trace.info() << "Usage: " << argv[ 0 ] << " <K> <input.vol> <m> <M>" << std::endl;
247 trace.info() << "\tAnalyze the shape with local full convexity" << std::endl;
248 trace.info() << "\t- 1 <= K <= 5: analysis at scale K" << std::endl;
249 trace.info() << "\t- K == 0: multiscale analysis (using scales 1-5)" << std::endl;
250 trace.info() << "\t- input.vol: choose your favorite shape" << std::endl;
251 trace.info() << "\t- m [==0], M [==255]: used to threshold input vol image" << std::endl;
252 return 1;
253 }
254 int N = atoi( argv[ 1 ] );
255 std::string fn= argv[ 2 ];
256 int m = argc > 3 ? atoi( argv[ 3 ] ) : 0;
257 int M = argc > 4 ? atoi( argv[ 4 ] ) : 255;
258
259 auto params = SH3::defaultParameters();
260
261 // Domain creation from two bounding points.
262 trace.info() << "Building set or importing vol ... ";
263 KSpace K;
264 params( "thresholdMin", m );
265 params( "thresholdMax", M );
266 auto bimage = SH3::makeBinaryImage( fn, params );
267 K = SH3::getKSpace( bimage );
268 trace.info() << " [Done]" << std::endl;
269 // Compute surface
270 params( "surfaceComponents" , "All" );
271 auto surface = SH3::makeDigitalSurface( bimage, K, params );
272
273 // Compute interior boundary points
274 // They are less immediate interior points than surfels.
275 std::vector< Point > points;
276 std::map< SCell, int > surfel2idx;
277 std::map< Point, int > point2idx;
278 int idx = 0;
279 for ( auto s : (*surface) )
280 {
281 // get inside point on the border of the shape.
282 Dimension k = K.sOrthDir( s );
283 auto voxel = K.sIncident( s, k, K.sDirect( s, k ) );
284 Point p = K.sCoords( voxel );
285 auto it = point2idx.find( p );
286 if ( it == point2idx.end() )
287 {
288 points.push_back( p );
289 surfel2idx[ s ] = idx;
290 point2idx [ p ] = idx++;
291 }
292 else
293 surfel2idx[ s ] = it->second;
294 }
295 trace.info() << "Shape has " << points.size() << " interior boundary points"
296 << std::endl;
297 if ( N != 0 )
298 {
299 std::vector< int > result;
300 trace.beginBlock ( "Single scale analysis" );
301 if ( N == 1 ) result = Analyzer< KSpace, 1 >::run( K, points, bimage );
302 if ( N == 2 ) result = Analyzer< KSpace, 2 >::run( K, points, bimage );
303 if ( N == 3 ) result = Analyzer< KSpace, 3 >::run( K, points, bimage );
304 if ( N == 4 ) result = Analyzer< KSpace, 4 >::run( K, points, bimage );
305 if ( N == 5 ) result = Analyzer< KSpace, 5 >::run( K, points, bimage );
306 trace.endBlock();
307 SCell dummy;
308 Color colors[ 4 ] =
309 { Color( 255, 0, 0, 255 ), Color( 0, 255, 0, 255 ),
310 Color( 0, 0, 255, 255 ), Color( 255, 255, 255, 255 ) };
311 auto surfels = SH3::getSurfelRange( surface, params );
312 SH3::Colors all_colors( surfels.size() );
313 for ( int i = 0; i < surfels.size(); i++ )
314 {
315 const auto j = surfel2idx[ surfels[ i ] ];
316 all_colors[ i ] = colors[ result[ j ] ];
317 }
318 SH3::saveOBJ( surface, SH3::RealVectors(), all_colors, "geom-cvx.obj" );
319 PolyscopeViewer viewer;
320 int i = 0;
321 for ( auto s : (*surface) )
322 {
323 viewer << all_colors[ i ]
324 << s;
325 i++;
326 }
327 viewer.show();
328 }
329 else
330 {
331 trace.beginBlock ( "Multiscale analysis" );
332 auto geometry =
334 trace.endBlock();
335 Color colors_planar[ 6 ] =
336 { Color( 0, 255, 255, 255),
337 Color( 50, 255, 255, 255), Color( 100, 255, 255, 255),
338 Color( 150, 255, 255, 255), Color( 200, 255, 255, 255 ),
339 Color( 255, 255, 255, 255 ) };
340 Color color_atypical( 255, 0, 0, 255 );
341 Color colors_cvx[ 5 ] =
342 { Color( 0, 255, 0, 255 ), Color( 50, 255, 50, 255 ),
343 Color( 100, 255, 100, 255 ), Color( 150, 255, 150, 255 ),
344 Color( 200, 255, 200, 255 ) };
345 Color colors_ccv[ 5 ] =
346 { Color( 0, 0, 255, 255 ), Color( 50, 50, 255, 255 ),
347 Color( 100, 100, 255, 255 ), Color( 150, 150, 255, 255 ),
348 Color( 200, 200, 255, 255 ) };
349 auto surfels = SH3::getSurfelRange( surface, params );
350 SH3::Colors all_colors( surfels.size() );
351 for ( int i = 0; i < surfels.size(); i++ ) {
352 const auto j = surfel2idx[ surfels[ i ] ];
353 int m0 = std::min( geometry[ j ].first, geometry[ j ].second );
354 int m1 = std::max( geometry[ j ].first, geometry[ j ].second );
355 if ( m1 == 0 ) all_colors[ i ] = color_atypical;
356 else if ( m0 == m1 ) all_colors[ i ] = colors_planar[ 5 ];
357 else if ( geometry[ j ].first > geometry[ j ].second )
358 all_colors[ i ] = colors_cvx[ 5 - abs( m0 - m1 ) ];
359 else
360 all_colors[ i ] = colors_ccv[ 5 - abs( m0 - m1 ) ];
361 }
362 SH3::saveOBJ( surface, SH3::RealVectors(), all_colors, "geom-scale-cvx.obj" );
363 SCell dummy;
364 int i = 0;
365 PolyscopeViewer<> viewer;
366 for ( auto s : (*surface) )
367 {
368 viewer << all_colors[ i ]
369 << s;
370 i++;
371 }
372 viewer.show();
373 }
374 return 0;
375}
376// //
378
Structure representing an RGB triple with alpha component.
Definition Color.h:77
const Point & lowerBound() const
Return the lower bound for digital points in this space.
PointVector< dim, Integer > Point
const Point & upperBound() const
Return the upper bound for digital points in this space.
static const constexpr Dimension dimension
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
void show() override
Starts the event loop and display of elements.
static KSpace getKSpace(const Point &low, const Point &up, Parameters params=parametersKSpace())
Definition Shortcuts.h:333
std::vector< Color > Colors
Definition Shortcuts.h:193
static SurfelRange getSurfelRange(CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer > > surface, const Parameters &params=parametersDigitalSurface())
Definition Shortcuts.h:1548
static CountedPtr< DigitalSurface > makeDigitalSurface(CountedPtr< TPointPredicate > bimage, const KSpace &K, const Parameters &params=parametersDigitalSurface())
Definition Shortcuts.h:1210
std::vector< RealVector > RealVectors
Definition Shortcuts.h:180
static Parameters defaultParameters()
Definition Shortcuts.h:204
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:1740
static CountedPtr< BinaryImage > makeBinaryImage(Domain shapeDomain)
Definition Shortcuts.h:562
CountedPtr< SH3::DigitalSurface > surface
Z3i this namespace gathers the standard of types for 3D imagery.
KhalimskySpaceND< 3, Integer > KSpace
Definition StdDefs.h:146
KSpace::SCell SCell
Definition StdDefs.h:149
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition Common.h:119
Trace trace
STL namespace.
static std::vector< Point > debug_one(const KSpace &aK, Point p, ImagePtr bimage)
NeighborhoodConvexityAnalyzer< KSpace, N > NCA
static void run(std::vector< int > &to_update, const KSpace &aK, const std::vector< Point > &pts, ImagePtr bimage)
static std::vector< int > run(const KSpace &aK, const std::vector< Point > &pts, ImagePtr bimage)
KSpace::Point Point
static std::vector< Geometry > multiscale_run(const KSpace &aK, const std::vector< Point > &pts, ImagePtr bimage)
NeighborhoodConvexityAnalyzer< KSpace, N > NCA
static std::vector< Geometry > multiscale_run(const KSpace &aK, const std::vector< Point > &pts, ImagePtr bimage)
std::pair< int, int > Geometry
Shortcuts< KSpace > SH3
int main()
Definition testBits.cpp:56
KSpace K
Image image(domain)