99int main(
int argc,
char** argv )
101 trace.
info() <<
"Usage: " << argv[ 0 ] <<
" <input.vol> <m> <M> <opt>" << std::endl;
102 trace.
info() <<
"\tComputes shortest paths to a source point" << std::endl;
103 trace.
info() <<
"\t- input.vol: choose your favorite shape" << std::endl;
104 trace.
info() <<
"\t- m [==0], M [==255]: used to threshold input vol image" << std::endl;
105 trace.
info() <<
"\t- opt >= sqrt(3): secure shortest paths, 0: fast" << std::endl;
106 string inputFilename = examplesPath +
"samples/Al.100.vol";
107 std::string fn= argc > 1 ? argv[ 1 ] : inputFilename;
108 int m = argc > 2 ? atoi( argv[ 2 ] ) : 0;
109 int M = argc > 3 ? atoi( argv[ 3 ] ) : 255;
110 double opt = argc > 4 ? atof( argv[ 4 ] ) : sqrt(3.0);
115 auto params = SH3::defaultParameters();
116 params(
"thresholdMin", m )(
"thresholdMax", M );
117 params(
"surfaceComponents" ,
"All" );
120 trace.
info() <<
"Building set or importing vol ... ";
122 auto bimage = SH3::makeBinaryImage( fn, params );
123 K = SH3::getKSpace( bimage );
127 auto surface = SH3::makeDigitalSurface( bimage,
K, params );
131 std::vector< Point > points;
132 std::map< SCell, int > surfel2idx;
133 std::map< Point, int > point2idx;
135 for (
auto s : (*surface) )
141 auto it = point2idx.find( p );
142 if ( it == point2idx.end() )
144 points.push_back( p );
145 surfel2idx[ s ] = idx;
146 point2idx [ p ] = idx++;
149 surfel2idx[ s ] = it->second;
151 trace.
info() <<
"Shape has " << points.size() <<
" interior boundary points"
157 auto surfels = SH3::getSurfelRange (
surface );
158 for (
int i = 0; i < 2; i++ )
160 MViewer3D viewerCore(
K );
161 Color colSurfel( 200, 200, 255, 255 );
162 Color colStart( 255, 0, 0, 255 );
163 viewerCore.drawColor( colSurfel );
164 for (
auto && s : surfels ) viewerCore << s;
170 const auto s0 = surfels[ selected_surfels[ 0 ] ];
174 auto start0 = point2idx[ p0 ];
175 std::cout <<
"Start0 index is " << start0 << std::endl;
176 const auto s1 = surfels[ selected_surfels[ 1 ] ];
180 auto start1 = point2idx[ p1 ];
181 std::cout <<
"Start1 index is " << start1 << std::endl;
188 TC.init( points.cbegin(), points.cend() );
189 auto SP = TC.makeShortestPaths( opt );
191 double last_distance = 0.0;
192 while ( ! SP.finished() )
194 last_distance = std::get<2>( SP.current() );
197 std::cout <<
"Max distance is " << last_distance << std::endl;
201 const int nb_repetitions = 10;
202 const double period = last_distance / nb_repetitions;
204 MViewer3D viewerCore;
205 Color colSurfel( 200, 200, 255, 128 );
206 Color colStart( 255, 0, 0, 128 );
208 for (
auto i = 0; i < points.size(); ++i )
210 const double d_s = SP.distance( i );
211 Color c_s = cmap( fmod( d_s, period ) );
212 viewerCore.drawColor( c_s );
213 viewerCore.drawBall(
RealPoint( points[ i ][ 0 ],
215 points[ i ][ 2 ] ) );
220 auto surfels = SH3::getSurfelRange (
surface );
221 for (
auto && s : surfels )
223 const double d_s = SP.distance( surfel2idx[ s ] );
224 Color c_s = cmap( fmod( d_s, period ) );
225 viewerCore.drawColor( c_s );
228 viewerCore.drawColor( colStart );
230 for (
Index i = 0; i < SP.size(); i++ ) {
231 Point p1 = SP.point( i );
232 Point p2 = SP.point( SP.ancestor( i ) );
233 viewerCore.drawLine( p1, p2 );
241 auto SP0 = TC.makeShortestPaths( opt );
242 auto SP1 = TC.makeShortestPaths( opt );
245 std::vector< Index > Q;
246 while ( ! SP0.finished() && ! SP1.finished() )
248 auto n0_ = SP0.current(); ((void) n0_);
249 auto n1_ = SP1.current();
251 auto p1_ = std::get<0>( n1_ );
254 if ( SP0.isVisited( p1_ ) )
256 auto c0 = SP0.pathToSource( p1_ );
257 auto c1 = SP1.pathToSource( p1_ );
258 std::copy(c0.rbegin(), c0.rend(), std::back_inserter(Q));
260 std::copy(c1.begin(), c1.end(), std::back_inserter(Q));
269 const int nb_repetitions = 10;
270 const double period = last_distance / nb_repetitions;
272 MViewer3D viewerCore;
273 Color colSurfel( 200, 200, 255, 128 );
274 Color colStart( 255, 0, 0, 128 );
275 for (
auto i = 0; i < points.size(); ++i )
277 const double d_s0 = SP0.isVisited( i ) ? SP0.distance( i ) : SP0.infinity();
278 const double d_s1 = SP1.isVisited( i ) ? SP1.distance( i ) : SP1.infinity();
279 const double d_s = std::min( d_s0, d_s1 );
280 Color c_s = ( d_s != SP0.infinity() )
281 ? cmap( fmod( d_s, period ) )
283 viewerCore.drawColor( c_s );
284 viewerCore.drawBall(
RealPoint( points[ i ][ 0 ],
286 points[ i ][ 2 ] ) );
290 for (
auto i = 1; i < Q.size(); i++ )
292 Point p1_ = TC.point( Q[ i-1 ] );
293 Point p2_ = TC.point( Q[ i ] );
294 viewerCore.drawLine( p1_, p2_ );
301 std::vector< Index > sources;
302 std::vector< Index > dests;
303 for (
int i = 0; i < 20; i++ )
304 sources.push_back( rand() % TC.size() );
305 dests.push_back( start0 );
306 dests.push_back( start1 );
307 auto paths = TC.shortestPaths( sources, dests, opt );
311 MViewer3D viewerCore;
312 Color colSurfel( 200, 200, 255, 128 );
313 Color colStart( 255, 0, 0, 128 );
314 for (
auto i = 0; i < points.size(); ++i )
316 viewerCore.drawColor(
Color( 150, 150, 150, 255 ) );
317 viewerCore.drawBall(
RealPoint( points[ i ][ 0 ],
319 points[ i ][ 2 ] ) );
322 for (
auto path : paths )
324 for (
auto i = 1; i < path.size(); i++ )
326 Point p1_ = TC.point( path[ i-1 ] );
327 Point p2_ = TC.point( path[ i ] );
328 viewerCore.drawLine( p1_, p2_ );
330 trace.
info() <<
"length=" << TC.length( path ) << std::endl;