DGtal 1.3.0
Searching...
No Matches
geometry/volumes/fullConvexityShortestPaths3D.cpp

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

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 Space::Point Point;
typedef Space::RealPoint RealPoint;
typedef Space::Vector Vector;
// 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 << 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 ] );
}
trace.info() << "length=" << TC.length( path ) << std::endl;
}
viewerCore << MViewer3D::updateDisplay;
application.exec();
}
return 0;
}
// //
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,...
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...
Z3i::SCell SCell
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:137
Trace trace
Definition: Common.h:154
boost::int32_t int32_t
signed 32-bit integer.
Definition: BasicTypes.h:72
STL namespace.
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.
Shortcuts< KSpace > SH3
Z2i::RealPoint RealPoint
int main()
Definition: testBits.cpp:56
MyPointD Point
Definition: testClone2.cpp:383
FreemanChain< int >::Vector Vector
KSpace K
HyperRectDomain< Space > Domain