DGtal  1.5.beta
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 
85 using namespace std;
86 using namespace DGtal;
87 using namespace Z3i;
89 
91 
92 
93 template < typename KSpace, int N >
94 struct 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 
195 template < typename KSpace, int N >
196 struct 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 
225 template < typename KSpace>
226 struct 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 
241 int 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,...
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 & lowerBound() const
Return the lower 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].
const Point & upperBound() const
Return the upper bound for digital points in this space.
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
std::vector< Color > Colors
Definition: Shortcuts.h:192
std::vector< RealVector > RealVectors
Definition: Shortcuts.h:179
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...
CountedPtr< SH3::DigitalSurface > surface
int main(int argc, char **argv)
Shortcuts< KSpace > SH3
NeighborhoodConvexityAnalyzer< KSpace, 1 > NCA
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition: Common.h:136
Trace trace
Definition: Common.h:153
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 max(int a, int b)
MyPointD Point
Definition: testClone2.cpp:383
KSpace K