DGtal  1.4.beta
fullConvexityShortestPaths3D.cpp
Go to the documentation of this file.
1 
60 #include <iostream>
61 #include <queue>
62 #include "DGtal/base/Common.h"
63 #include "DGtal/io/viewers/Viewer3D.h"
64 #include "DGtal/io/DrawWithDisplay3DModifier.h"
65 #include "DGtal/io/Color.h"
66 #include "DGtal/io/colormaps/SimpleDistanceColorMap.h"
67 #include "DGtal/shapes/Shapes.h"
68 #include "DGtal/helpers/StdDefs.h"
69 #include "DGtal/helpers/Shortcuts.h"
70 #include "DGtal/images/ImageContainerBySTLVector.h"
72 #include "DGtal/geometry/volumes/TangencyComputer.h"
74 #include "ConfigExamples.h"
75 
77 
78 using namespace std;
79 using namespace DGtal;
80 typedef Z3i::Space Space;
83 typedef Z3i::SCell SCell;
88 
89 // Called when an user clicks on a surfel.
90 int reaction( void* viewer, DGtal::int32_t name, void* data )
91 {
92  DGtal::int32_t* selected = (DGtal::int32_t*) data;
93  *selected = name;
94  std::cout << "Selected surfel=" << *selected << std::endl;
95  return 0;
96 }
97 
98 int main( int argc, char** argv )
99 {
100  trace.info() << "Usage: " << argv[ 0 ] << " <input.vol> <m> <M> <opt>" << std::endl;
101  trace.info() << "\tComputes shortest paths to a source point" << std::endl;
102  trace.info() << "\t- input.vol: choose your favorite shape" << std::endl;
103  trace.info() << "\t- m [==0], M [==255]: used to threshold input vol image" << std::endl;
104  trace.info() << "\t- opt >= sqrt(3): secure shortest paths, 0: fast" << std::endl;
105  string inputFilename = examplesPath + "samples/Al.100.vol";
106  std::string fn= argc > 1 ? argv[ 1 ] : inputFilename; //< vol filename
107  int m = argc > 2 ? atoi( argv[ 2 ] ) : 0; //< low for thresholding
108  int M = argc > 3 ? atoi( argv[ 3 ] ) : 255; //< up for thresholding
109  double opt = argc > 4 ? atof( argv[ 4 ] ) : sqrt(3.0); //< exact (sqrt(3)) or inexact (0) computations
110 
111  QApplication application(argc,argv);
112  Viewer3D<> viewer;
113  viewer.setWindowTitle("fullConvexityShortestPaths3D");
114  viewer.show();
115 
116  // Set up shortcuts parameters.
117  auto params = SH3::defaultParameters();
118  params( "thresholdMin", m )( "thresholdMax", M );
119  params( "surfaceComponents" , "All" );
120 
121  // Domain creation from two bounding points.
122  trace.info() << "Building set or importing vol ... ";
123  KSpace K;
124  auto bimage = SH3::makeBinaryImage( fn, params );
125  K = SH3::getKSpace( bimage );
126  trace.info() << " [Done]" << std::endl;
127 
128  // Compute surface
129  auto surface = SH3::makeDigitalSurface( bimage, K, params );
130 
131  // Compute interior boundary points
132  // They are less immediate interior points than surfels.
133  std::vector< Point > points;
134  std::map< SCell, int > surfel2idx;
135  std::map< Point, int > point2idx;
136  int idx = 0;
137  for ( auto s : (*surface) )
138  {
139  // get inside point on the border of the shape.
140  Dimension k = K.sOrthDir( s );
141  auto voxel = K.sIncident( s, k, K.sDirect( s, k ) );
142  Point p = K.sCoords( voxel );
143  auto it = point2idx.find( p );
144  if ( it == point2idx.end() )
145  {
146  points.push_back( p );
147  surfel2idx[ s ] = idx;
148  point2idx [ p ] = idx++;
149  }
150  else
151  surfel2idx[ s ] = it->second;
152  }
153  trace.info() << "Shape has " << points.size() << " interior boundary points"
154  << std::endl;
155 
156  // Select a starting point.
157  typedef Viewer3D<> MViewer3D;
158  DGtal::int32_t selected_surfels[ 2 ] = { 0, 0 };
159  auto surfels = SH3::getSurfelRange ( surface );
160  for ( int i = 0; i < 2; i++ )
161  {
162  MViewer3D viewerCore( K );
163  viewerCore.show();
164  Color colSurfel( 200, 200, 255, 255 );
165  Color colStart( 255, 0, 0, 255 );
166  DGtal::int32_t name = 0;
167  viewerCore << SetMode3D( surfels[ 0 ].className(), "Basic");
168  viewerCore.setFillColor( colSurfel );
169  for ( auto && s : surfels ) viewerCore << SetName3D( name++ ) << s;
170  viewerCore << SetSelectCallback3D( reaction, &selected_surfels[ i ],
171  0, surfels.size() - 1 );
172  viewerCore << MViewer3D::updateDisplay;
173  application.exec();
174  }
175 
176  // Get selected surfel/point
177  const auto s0 = surfels[ selected_surfels[ 0 ] ];
178  Dimension k0 = K.sOrthDir( s0 );
179  auto voxel0 = K.sIncident( s0, k0, K.sDirect( s0, k0 ) );
180  Point p0 = K.sCoords( voxel0 );
181  auto start0 = point2idx[ p0 ];
182  std::cout << "Start0 index is " << start0 << std::endl;
183  const auto s1 = surfels[ selected_surfels[ 1 ] ];
184  Dimension k1 = K.sOrthDir( s1 );
185  auto voxel1 = K.sIncident( s1, k1, K.sDirect( s1, k1 ) );
186  Point p1 = K.sCoords( voxel1 );
187  auto start1 = point2idx[ p1 ];
188  std::cout << "Start1 index is " << start1 << std::endl;
189 
190  // (I) Extracts shortest paths to a target
191 
195  TC.init( points.cbegin(), points.cend() );
196  auto SP = TC.makeShortestPaths( opt );
197  SP.init( start0 ); //< set source
198  double last_distance = 0.0;
199  while ( ! SP.finished() )
200  {
201  last_distance = std::get<2>( SP.current() );
202  SP.expand();
203  }
204  std::cout << "Max distance is " << last_distance << std::endl;
206 
207  {
208  const int nb_repetitions = 10;
209  const double period = last_distance / nb_repetitions;
210  SimpleDistanceColorMap< double > cmap( 0.0, period, false );
211  MViewer3D viewerCore;
212  viewerCore.show();
213  Color colSurfel( 200, 200, 255, 128 );
214  Color colStart( 255, 0, 0, 128 );
215 
216  viewerCore.setUseGLPointForBalls(true);
217  for ( auto i = 0; i < points.size(); ++i )
218  {
219  const double d_s = SP.distance( i );
220  Color c_s = cmap( fmod( d_s, period ) );
221  viewerCore.setFillColor( c_s );
222  viewerCore.addBall( RealPoint( points[ i ][ 0 ],
223  points[ i ][ 1 ],
224  points[ i ][ 2 ] ), 12.0 );
225  }
226 
227  // JOL: Left if you wish to display it as a surface, and to display paths.
228 
229  // auto surfels = SH3::getSurfelRange ( surface );
230  // viewerCore << SetMode3D( surfels[ 0 ].className(), "Basic");
231  // for ( auto && s : surfels )
232  // {
233  // const double d_s = SP.distance( surfel2idx[ s ] );
234  // Color c_s = cmap( fmod( d_s, period ) );
235  // viewerCore.setFillColor( c_s );
236  // viewerCore << s;
237  // }
238  // viewerCore.setFillColor( colStart );
239  // viewerCore.setLineColor( Color::Green );
240  // viewerCore << s;
241  // for ( Index i = 0; i < SP.size(); i++ ) {
242  // Point p1 = SP.point( i );
243  // Point p2 = SP.point( SP.ancestor( i ) );
244  // viewerCore.addLine( p1, p2, 1.0 );
245  // }
246  viewerCore << MViewer3D::updateDisplay;
247  application.exec();
248  }
249 
250  // (II) Extracts a shortest path between two points.
251 
253  auto SP0 = TC.makeShortestPaths( opt );
254  auto SP1 = TC.makeShortestPaths( opt );
255  SP0.init( start0 ); //< source point
256  SP1.init( start1 ); //< target point
257  std::vector< Index > Q; //< the output shortest path
258  while ( ! SP0.finished() && ! SP1.finished() )
259  {
260  auto n0 = SP0.current();
261  auto n1 = SP1.current();
262  auto p0 = std::get<0>( n0 );
263  auto p1 = std::get<0>( n1 );
264  SP0.expand();
265  SP1.expand();
266  if ( SP0.isVisited( p1 ) )
267  {
268  auto c0 = SP0.pathToSource( p1 );
269  auto c1 = SP1.pathToSource( p1 );
270  std::copy(c0.rbegin(), c0.rend(), std::back_inserter(Q));
271  Q.pop_back();
272  std::copy(c1.begin(), c1.end(), std::back_inserter(Q));
273  break;
274  }
275  }
276  // Q is empty if there is no path.
278 
279  // Display computed distances and shortest path
280  {
281  const int nb_repetitions = 10;
282  const double period = last_distance / nb_repetitions;
283  SimpleDistanceColorMap< double > cmap( 0.0, period, false );
284  MViewer3D viewerCore;
285  viewerCore.show();
286  Color colSurfel( 200, 200, 255, 128 );
287  Color colStart( 255, 0, 0, 128 );
288  viewerCore.setUseGLPointForBalls(true);
289  for ( auto i = 0; i < points.size(); ++i )
290  {
291  const double d_s0 = SP0.isVisited( i ) ? SP0.distance( i ) : SP0.infinity();
292  const double d_s1 = SP1.isVisited( i ) ? SP1.distance( i ) : SP1.infinity();
293  const double d_s = std::min( d_s0, d_s1 );
294  Color c_s = ( d_s != SP0.infinity() )
295  ? cmap( fmod( d_s, period ) )
296  : Color::Black;
297  viewerCore.setFillColor( c_s );
298  viewerCore.addBall( RealPoint( points[ i ][ 0 ],
299  points[ i ][ 1 ],
300  points[ i ][ 2 ] ), 12 );
301  }
302  viewerCore.setLineColor( Color::Green );
303  for ( auto i = 1; i < Q.size(); i++ )
304  {
305  Point p1 = TC.point( Q[ i-1 ] );
306  Point p2 = TC.point( Q[ i ] );
307  viewerCore.addLine( p1, p2, 18.0 );
308  }
309  viewerCore << MViewer3D::updateDisplay;
310  application.exec();
311  }
312 
313  // (III) Extracts multiple shortest paths between many sources and two targets.
314 
315  std::vector< Index > sources;
316  std::vector< Index > dests;
317  for ( int i = 0; i < 20; i++ )
318  sources.push_back( rand() % TC.size() );
319  dests.push_back( start0 );
320  dests.push_back( start1 );
321  auto paths = TC.shortestPaths( sources, dests, opt );
322 
323  // Display them.
324  {
325  MViewer3D viewerCore;
326  viewerCore.show();
327  Color colSurfel( 200, 200, 255, 128 );
328  Color colStart( 255, 0, 0, 128 );
329  viewerCore.setUseGLPointForBalls(true);
330  for ( auto i = 0; i < points.size(); ++i )
331  {
332  viewerCore.setFillColor( Color( 150, 150, 150, 255 ) );
333  viewerCore.addBall( RealPoint( points[ i ][ 0 ],
334  points[ i ][ 1 ],
335  points[ i ][ 2 ] ), 12 );
336  }
337  viewerCore.setLineColor( Color::Green );
338  for ( auto path : paths )
339  {
340  for ( auto i = 1; i < path.size(); i++ )
341  {
342  Point p1 = TC.point( path[ i-1 ] );
343  Point p2 = TC.point( path[ i ] );
344  viewerCore.addLine( p1, p2, 18.0 );
345  }
346  trace.info() << "length=" << TC.length( path ) << std::endl;
347  }
348  viewerCore << MViewer3D::updateDisplay;
349  application.exec();
350  }
351 
352  return 0;
353 }
354 // //
356 
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.
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].
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:105
Aim: simple blue to red colormap for distance information for instance.
Aim: A class that computes tangency to a given digital set. It provides services to compute all the c...
std::ostream & info()
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
Space::Point Point
int main(int argc, char **argv)
Z3i::KSpace KSpace
int reaction(void *viewer, DGtal::int32_t name, void *data)
Z3i::SCell SCell
Space::Vector Vector
Shortcuts< KSpace > SH3
Z3i::Space Space
Space::RealPoint RealPoint
Z3i::Domain Domain
SMesh::Index Index
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
boost::int32_t int32_t
signed 32-bit integer.
Definition: BasicTypes.h:72
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.
MyPointD Point
Definition: testClone2.cpp:383
KSpace K