DGtal  1.5.beta
geometry/volumes/fullConvexityShortestPaths3D.cpp

This example shows how to use tangency to compute shortest paths on 3D digital objects

See also
Tangency and shortest paths

For instance, you may call it on object "cube+sphere" as

fullConvexityShortestPaths3D cps.vol 0 255 0.0

The user selects two surfels (with shift + left click), and then shortest paths are computed and displayed.

Geodesic distances and geodesics on cube+sphere shape
Geodesic distances on cube+sphere shape
Shortest path between two points
#include <iostream>
#include <queue>
#include "DGtal/base/Common.h"
#include "DGtal/io/viewers/Viewer3D.h"
#include "DGtal/io/DrawWithDisplay3DModifier.h"
#include "DGtal/io/Color.h"
#include "DGtal/io/colormaps/SimpleDistanceColorMap.h"
#include "DGtal/shapes/Shapes.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/helpers/Shortcuts.h"
#include "DGtal/images/ImageContainerBySTLVector.h"
#include "DGtal/geometry/volumes/TangencyComputer.h"
#include "ConfigExamples.h"
using namespace std;
using namespace DGtal;
typedef Z3i::Space Space;
typedef Z3i::SCell SCell;
typedef Shortcuts< KSpace > SH3;
// Called when an user clicks on a surfel.
int reaction( void* viewer, DGtal::int32_t name, void* data )
{
DGtal::int32_t* selected = (DGtal::int32_t*) data;
*selected = name;
std::cout << "Selected surfel=" << *selected << std::endl;
return 0;
}
int main( int argc, char** argv )
{
trace.info() << "Usage: " << argv[ 0 ] << " <input.vol> <m> <M> <opt>" << std::endl;
trace.info() << "\tComputes shortest paths to a source point" << std::endl;
trace.info() << "\t- input.vol: choose your favorite shape" << std::endl;
trace.info() << "\t- m [==0], M [==255]: used to threshold input vol image" << std::endl;
trace.info() << "\t- opt >= sqrt(3): secure shortest paths, 0: fast" << std::endl;
string inputFilename = examplesPath + "samples/Al.100.vol";
std::string fn= argc > 1 ? argv[ 1 ] : inputFilename; //< vol filename
int m = argc > 2 ? atoi( argv[ 2 ] ) : 0; //< low for thresholding
int M = argc > 3 ? atoi( argv[ 3 ] ) : 255; //< up for thresholding
double opt = argc > 4 ? atof( argv[ 4 ] ) : sqrt(3.0); //< exact (sqrt(3)) or inexact (0) computations
QApplication application(argc,argv);
Viewer3D<> viewer;
viewer.setWindowTitle("fullConvexityShortestPaths3D");
viewer.show();
// Set up shortcuts parameters.
auto params = SH3::defaultParameters();
params( "thresholdMin", m )( "thresholdMax", M );
params( "surfaceComponents" , "All" );
// Domain creation from two bounding points.
trace.info() << "Building set or importing vol ... ";
auto bimage = SH3::makeBinaryImage( fn, params );
K = SH3::getKSpace( bimage );
trace.info() << " [Done]" << std::endl;
// Compute surface
auto surface = SH3::makeDigitalSurface( bimage, K, params );
// Compute interior boundary points
// They are less immediate interior points than surfels.
std::vector< Point > points;
std::map< SCell, int > surfel2idx;
std::map< Point, int > point2idx;
int idx = 0;
for ( auto s : (*surface) )
{
// get inside point on the border of the shape.
Dimension k = K.sOrthDir( s );
auto voxel = K.sIncident( s, k, K.sDirect( s, k ) );
Point p = K.sCoords( voxel );
auto it = point2idx.find( p );
if ( it == point2idx.end() )
{
points.push_back( p );
surfel2idx[ s ] = idx;
point2idx [ p ] = idx++;
}
else
surfel2idx[ s ] = it->second;
}
trace.info() << "Shape has " << points.size() << " interior boundary points"
<< std::endl;
// Select a starting point.
typedef Viewer3D<> MViewer3D;
DGtal::int32_t selected_surfels[ 2 ] = { 0, 0 };
auto surfels = SH3::getSurfelRange ( surface );
for ( int i = 0; i < 2; i++ )
{
MViewer3D viewerCore( K );
viewerCore.show();
Color colSurfel( 200, 200, 255, 255 );
Color colStart( 255, 0, 0, 255 );
DGtal::int32_t name = 0;
viewerCore << SetMode3D( surfels[ 0 ].className(), "Basic");
viewerCore.setFillColor( colSurfel );
for ( auto && s : surfels ) viewerCore << SetName3D( name++ ) << s;
viewerCore << SetSelectCallback3D( reaction, &selected_surfels[ i ],
0, surfels.size() - 1 );
viewerCore << MViewer3D::updateDisplay;
application.exec();
}
// Get selected surfel/point
const auto s0 = surfels[ selected_surfels[ 0 ] ];
Dimension k0 = K.sOrthDir( s0 );
auto voxel0 = K.sIncident( s0, k0, K.sDirect( s0, k0 ) );
Point p0 = K.sCoords( voxel0 );
auto start0 = point2idx[ p0 ];
std::cout << "Start0 index is " << start0 << std::endl;
const auto s1 = surfels[ selected_surfels[ 1 ] ];
Dimension k1 = K.sOrthDir( s1 );
auto voxel1 = K.sIncident( s1, k1, K.sDirect( s1, k1 ) );
Point p1 = K.sCoords( voxel1 );
auto start1 = point2idx[ p1 ];
std::cout << "Start1 index is " << start1 << std::endl;
// (I) Extracts shortest paths to a target
TangencyComputer< KSpace > TC( K );
TC.init( points.cbegin(), points.cend() );
auto SP = TC.makeShortestPaths( opt );
SP.init( start0 ); //< set source
double last_distance = 0.0;
while ( ! SP.finished() )
{
last_distance = std::get<2>( SP.current() );
SP.expand();
}
std::cout << "Max distance is " << last_distance << std::endl;
{
const int nb_repetitions = 10;
const double period = last_distance / nb_repetitions;
SimpleDistanceColorMap< double > cmap( 0.0, period, false );
MViewer3D viewerCore;
viewerCore.show();
Color colSurfel( 200, 200, 255, 128 );
Color colStart( 255, 0, 0, 128 );
viewerCore.setUseGLPointForBalls(true);
for ( auto i = 0; i < points.size(); ++i )
{
const double d_s = SP.distance( i );
Color c_s = cmap( fmod( d_s, period ) );
viewerCore.setFillColor( c_s );
viewerCore.addBall( RealPoint( points[ i ][ 0 ],
points[ i ][ 1 ],
points[ i ][ 2 ] ), 12.0 );
}
// JOL: Left if you wish to display it as a surface, and to display paths.
// auto surfels = SH3::getSurfelRange ( surface );
// viewerCore << SetMode3D( surfels[ 0 ].className(), "Basic");
// for ( auto && s : surfels )
// {
// const double d_s = SP.distance( surfel2idx[ s ] );
// Color c_s = cmap( fmod( d_s, period ) );
// viewerCore.setFillColor( c_s );
// viewerCore << s;
// }
// viewerCore.setFillColor( colStart );
// viewerCore.setLineColor( Color::Green );
// viewerCore << s;
// for ( Index i = 0; i < SP.size(); i++ ) {
// Point p1 = SP.point( i );
// Point p2 = SP.point( SP.ancestor( i ) );
// viewerCore.addLine( p1, p2, 1.0 );
// }
viewerCore << MViewer3D::updateDisplay;
application.exec();
}
// (II) Extracts a shortest path between two points.
auto SP0 = TC.makeShortestPaths( opt );
auto SP1 = TC.makeShortestPaths( opt );
SP0.init( start0 ); //< source point
SP1.init( start1 ); //< target point
std::vector< Index > Q; //< the output shortest path
while ( ! SP0.finished() && ! SP1.finished() )
{
auto n0 = SP0.current();
auto n1 = SP1.current();
auto p0 = std::get<0>( n0 );
auto p1 = std::get<0>( n1 );
SP0.expand();
SP1.expand();
if ( SP0.isVisited( p1 ) )
{
auto c0 = SP0.pathToSource( p1 );
auto c1 = SP1.pathToSource( p1 );
std::copy(c0.rbegin(), c0.rend(), std::back_inserter(Q));
Q.pop_back();
std::copy(c1.begin(), c1.end(), std::back_inserter(Q));
break;
}
}
// Q is empty if there is no path.
// Display computed distances and shortest path
{
const int nb_repetitions = 10;
const double period = last_distance / nb_repetitions;
SimpleDistanceColorMap< double > cmap( 0.0, period, false );
MViewer3D viewerCore;
viewerCore.show();
Color colSurfel( 200, 200, 255, 128 );
Color colStart( 255, 0, 0, 128 );
viewerCore.setUseGLPointForBalls(true);
for ( auto i = 0; i < points.size(); ++i )
{
const double d_s0 = SP0.isVisited( i ) ? SP0.distance( i ) : SP0.infinity();
const double d_s1 = SP1.isVisited( i ) ? SP1.distance( i ) : SP1.infinity();
const double d_s = std::min( d_s0, d_s1 );
Color c_s = ( d_s != SP0.infinity() )
? cmap( fmod( d_s, period ) )
: Color::Black;
viewerCore.setFillColor( c_s );
viewerCore.addBall( RealPoint( points[ i ][ 0 ],
points[ i ][ 1 ],
points[ i ][ 2 ] ), 12 );
}
viewerCore.setLineColor( Color::Green );
for ( auto i = 1; i < Q.size(); i++ )
{
Point p1 = TC.point( Q[ i-1 ] );
Point p2 = TC.point( Q[ i ] );
viewerCore.addLine( p1, p2, 18.0 );
}
viewerCore << MViewer3D::updateDisplay;
application.exec();
}
// (III) Extracts multiple shortest paths between many sources and two targets.
std::vector< Index > sources;
std::vector< Index > dests;
for ( int i = 0; i < 20; i++ )
sources.push_back( rand() % TC.size() );
dests.push_back( start0 );
dests.push_back( start1 );
auto paths = TC.shortestPaths( sources, dests, opt );
// Display them.
{
MViewer3D viewerCore;
viewerCore.show();
Color colSurfel( 200, 200, 255, 128 );
Color colStart( 255, 0, 0, 128 );
viewerCore.setUseGLPointForBalls(true);
for ( auto i = 0; i < points.size(); ++i )
{
viewerCore.setFillColor( Color( 150, 150, 150, 255 ) );
viewerCore.addBall( RealPoint( points[ i ][ 0 ],
points[ i ][ 1 ],
points[ i ][ 2 ] ), 12 );
}
viewerCore.setLineColor( Color::Green );
for ( auto path : paths )
{
for ( auto i = 1; i < path.size(); i++ )
{
Point p1 = TC.point( path[ i-1 ] );
Point p2 = TC.point( path[ i ] );
viewerCore.addLine( p1, p2, 18.0 );
}
trace.info() << "length=" << TC.length( path ) << std::endl;
}
viewerCore << MViewer3D::updateDisplay;
application.exec();
}
return 0;
}
// //
std::ostream & info()
CountedPtr< SH3::DigitalSurface > surface
Z3i::SCell SCell
DigitalPlane::Point Vector
int reaction(void *viewer, DGtal::int32_t name, void *data)
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
int main(int argc, char **argv)
Shortcuts< KSpace > SH3
MyPointD Point
Definition: testClone2.cpp:383
KSpace K
HyperRectDomain< Space > Domain
PointVector< 3, double > RealPoint