fullConvexityShortestPaths3D.cpp

Go to the documentation of this file.

 
#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::KSpace         KSpace;
typedef Z3i::Domain         Domain;
typedef Z3i::SCell          SCell;
typedef Shortcuts< KSpace > SH3;
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
 
  PolyscopeViewer<> viewer;
 
  // 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 ... ";
  KSpace K;
  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.
  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
 
  typedef TangencyComputer< KSpace >::Index Index;
  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;
    Color colSurfel( 200, 200, 255, 128 );
    Color colStart( 255, 0, 0, 128 );
 
    for ( auto 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.
    
    auto surfels = SH3::getSurfelRange ( surface );
    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 ( 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.drawColor( c_s );
        viewerCore.drawBall( RealPoint( points[ i ][ 0 ],
                                       points[ i ][ 1 ],
                                       points[ i ][ 2 ] ) );
      }
 
    viewerCore.drawColor( Color::Green );
    for ( auto 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 ( auto 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 ( auto 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;
}
//                                                                           //