DGtal 2.1.0
Loading...
Searching...
No Matches
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/PolyscopeViewer.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 )
{
((void) viewer);
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
// 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
const 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.
DGtal::int32_t selected_surfels[ 2 ] = { 0, 0 };
typedef PolyscopeViewer<> MViewer3D;
auto surfels = SH3::getSurfelRange ( surface );
for ( int i = 0; i < 2; i++ )
{
MViewer3D viewerCore( K );
Color colSurfel( 200, 200, 255, 255 );
Color colStart( 255, 0, 0, 255 );
viewerCore.drawColor( colSurfel );
for ( auto && s : surfels ) viewerCore << s;
viewerCore.show();
}
// 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;
Color colSurfel( 200, 200, 255, 128 );
Color colStart( 255, 0, 0, 128 );
for ( size_t i = 0; i < points.size(); ++i )
{
const double d_s = SP.distance( i );
Color c_s = cmap( fmod( d_s, period ) );
viewerCore.drawColor( c_s );
viewerCore.drawBall( RealPoint( points[ i ][ 0 ],
points[ i ][ 1 ],
points[ i ][ 2 ] ) );
}
// JOL: Left if you wish to display it as a surface, and to display paths.
for ( auto && s : surfels )
{
const double d_s = SP.distance( surfel2idx[ s ] );
Color c_s = cmap( fmod( d_s, period ) );
viewerCore.drawColor( c_s );
viewerCore << s;
}
viewerCore.drawColor( colStart );
viewerCore.drawColor( Color::Green );
for ( Index i = 0; i < SP.size(); i++ ) {
Point p1_ = SP.point( i );
Point p2_ = SP.point( SP.ancestor( i ) );
viewerCore.drawLine( p1_, p2_ );
}
viewerCore.show();
}
// (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(); ((void) n0_);
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;
Color colSurfel( 200, 200, 255, 128 );
Color colStart( 255, 0, 0, 128 );
for ( size_t 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.drawColor( c_s );
viewerCore.drawBall( RealPoint( points[ i ][ 0 ],
points[ i ][ 1 ],
points[ i ][ 2 ] ) );
}
viewerCore.drawColor( Color::Green );
for ( size_t i = 1; i < Q.size(); i++ )
{
Point p1_ = TC.point( Q[ i-1 ] );
Point p2_ = TC.point( Q[ i ] );
viewerCore.drawLine( p1_, p2_ );
}
viewerCore.show();
}
// (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;
Color colSurfel( 200, 200, 255, 128 );
Color colStart( 255, 0, 0, 128 );
for ( size_t i = 0; i < points.size(); ++i )
{
viewerCore.drawColor( Color( 150, 150, 150, 255 ) );
viewerCore.drawBall( RealPoint( points[ i ][ 0 ],
points[ i ][ 1 ],
points[ i ][ 2 ] ) );
}
viewerCore.drawColor( Color::Green );
for ( auto path : paths )
{
for ( size_t i = 1; i < path.size(); i++ )
{
Point p1_ = TC.point( path[ i-1 ] );
Point p2_ = TC.point( path[ i ] );
viewerCore.drawLine( p1_, p2_ );
}
trace.info() << "length=" << TC.length( path ) << std::endl;
}
viewerCore.show();
}
return 0;
}
// //
Structure representing an RGB triple with alpha component.
Definition Color.h:77
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:102
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()
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.
std::int32_t int32_t
signed 32-bit integer.
Definition BasicTypes.h:71
DGtal::uint32_t Dimension
Definition Common.h:119
Trace trace
STL namespace.
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
Shortcuts< KSpace > SH3
int main()
Definition testBits.cpp:56
MyPointD Point
KSpace K
HyperRectDomain< Space > Domain
PointVector< 3, double > RealPoint