DGtal  1.3.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;
// 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
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;
}
// //
DGtal::TangencyComputer::Index
std::size_t Index
Definition: TangencyComputer.h:83
DGtal::HyperRectDomain< Space >
SH3
Shortcuts< KSpace > SH3
Definition: testArithmeticalDSSComputerOnSurfels.cpp:49
DGtal::Color
Structure representing an RGB triple with alpha component.
Definition: Color.h:67
DGtal::int32_t
boost::int32_t int32_t
signed 32-bit integer.
Definition: BasicTypes.h:72
Index
SMesh::Index Index
Definition: fullConvexitySphereGeodesics.cpp:117
DGtal::trace
Trace trace
Definition: Common.h:154
DGtal::TangencyComputer
Aim: A class that computes tangency to a given digital set. It provides services to compute all the c...
Definition: TangencyComputer.h:72
K
KSpace K
Definition: testCubicalComplex.cpp:62
DGtal::Dimension
DGtal::uint32_t Dimension
Definition: Common.h:137
DGtal::Shortcuts::makeDigitalSurface
static CountedPtr< DigitalSurface > makeDigitalSurface(CountedPtr< TPointPredicate > bimage, const KSpace &K, const Parameters &params=parametersDigitalSurface())
Definition: Shortcuts.h:1209
reaction
int reaction(void *viewer, DGtal::int32_t name, void *data)
Definition: fullConvexityShortestPaths3D.cpp:90
SCell
Z3i::SCell SCell
Definition: fullConvexityShortestPaths3D.cpp:83
DGtal::Shortcuts
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:104
DGtal::KhalimskySpaceND::sCoords
Point sCoords(const SCell &c) const
Return its digital coordinates.
DGtal::SignedKhalimskyCell
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
Definition: KhalimskySpaceND.h:208
DGtal::KhalimskySpaceND::sDirect
bool sDirect(const SCell &p, Dimension k) const
Return 'true' if the direct orientation of [p] along [k] is in the positive coordinate direction.
DGtal::SetName3D
Definition: DrawWithDisplay3DModifier.h:242
DGtal::SpaceND
Definition: SpaceND.h:95
KSpace
Z3i::KSpace KSpace
Definition: testArithmeticalDSSComputerOnSurfels.cpp:48
DGtal::SimpleDistanceColorMap
Aim: simple blue to red colormap for distance information for instance.
Definition: SimpleDistanceColorMap.h:64
DGtal::Trace::info
std::ostream & info()
DGtal::Viewer3D
Definition: Viewer3D.h:135
DGtal::Shortcuts::defaultParameters
static Parameters defaultParameters()
Definition: Shortcuts.h:203
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
Domain
HyperRectDomain< Space > Domain
Definition: testSimpleRandomAccessRangeFromPoint.cpp:44
Vector
FreemanChain< int >::Vector Vector
Definition: testCombinDSS.cpp:60
DGtal::SetMode3D
Modifier class in a Display3D stream. Useful to choose your own mode for a given class....
Definition: DrawWithDisplay3DModifier.h:73
main
int main(int argc, char **argv)
Definition: testArithmeticDSS-benchmark.cpp:147
DGtal::KhalimskySpaceND::sOrthDir
Dimension sOrthDir(const SCell &s) const
Given a signed surfel [s], returns its orthogonal direction (ie, the coordinate where the surfel is c...
DGtal::Shortcuts::getKSpace
static KSpace getKSpace(const Point &low, const Point &up, Parameters params=parametersKSpace())
Definition: Shortcuts.h:332
DGtal::PointVector< dim, Integer >
DGtal::Viewer3D::show
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
Space
SpaceND< 2 > Space
Definition: testSimpleRandomAccessRangeFromPoint.cpp:42
DGtal::KhalimskySpaceND::sIncident
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].
DGtal::Shortcuts::getSurfelRange
static SurfelRange getSurfelRange(CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer > > surface, const Parameters &params=parametersDigitalSurface())
Definition: Shortcuts.h:1547
DGtal::Shortcuts::makeBinaryImage
static CountedPtr< BinaryImage > makeBinaryImage(Domain shapeDomain)
Definition: Shortcuts.h:561
DGtal::SetSelectCallback3D
Definition: DrawWithDisplay3DModifier.h:257
Point
MyPointD Point
Definition: testClone2.cpp:383
RealPoint
Z2i::RealPoint RealPoint
Definition: testAstroid2D.cpp:46
DGtal::KhalimskySpaceND
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Definition: KhalimskySpaceND.h:64