testShortestPaths.cpp

Go to the documentation of this file.

#include <iostream>
#include <vector>
#include <algorithm>
#include "DGtal/base/Common.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/helpers/Shortcuts.h"
#include "DGtal/geometry/volumes/TangencyComputer.h"
#include "DGtalCatch.h"
 
using namespace std;
using namespace DGtal;
 

// Functions for testing class TangencyComputer::ShortestPaths
 
SCENARIO( "TangencyComputer::ShortestPaths 3D tests", "[shortest_paths][3d][tangency]" )
{
  typedef Z3i::Space          Space;
  typedef Z3i::KSpace         KSpace;
  typedef Shortcuts< KSpace > SH3;
  typedef Space::Point        Point;
  typedef std::size_t         Index;
 
  SECTION( "Computing shortest paths on a 3D unit sphere digitized at gridstep 0.25" )
    {
      // Make digital sphere
      const double h = 0.25;
      auto   params  = SH3::defaultParameters();
      params( "polynomial", "sphere1" )( "gridstep",  h );
      params( "minAABB", -2)( "maxAABB", 2)( "offset", 1.0 )( "closed", 1 );
      auto implicit_shape  = SH3::makeImplicitShape3D  ( params );
      auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
      auto K            = SH3::getKSpace( params );
      auto binary_image = SH3::makeBinaryImage(digitized_shape,
                                               SH3::Domain(K.lowerBound(),K.upperBound()),
                                               params );
      auto surface = SH3::makeDigitalSurface( binary_image, K, params );
      std::vector< Point >    lattice_points;
      auto pointels = SH3::getPointelRange( surface );
      for ( auto p : pointels ) lattice_points.push_back( K.uCoords( p ) );
      REQUIRE( pointels.size() == 296 );
      // Find lowest and uppest point.
      const Index nb = lattice_points.size();
      Index   lowest = 0;
      Index   uppest = 0;
      for ( Index i = 1; i < nb; i++ )
        {
          if ( lattice_points[ i ] < lattice_points[ lowest ] ) lowest = i;
          if ( lattice_points[ uppest ] < lattice_points[ i ] ) uppest = i;
        }
      // Compute shortest paths
      typedef TangencyComputer< KSpace >::Index _Index;
      TangencyComputer< KSpace > TC( K );
      TC.init( lattice_points.cbegin(), lattice_points.cend() );
      auto SP = TC.makeShortestPaths( sqrt(3.0) );
      SP.init( lowest ); //< set source
      double last_distance = 0.0;
      _Index  last = 0;
      double prev_distance = 0.0;
      std::set< _Index > V;
      unsigned int nb_multiple_pops = 0;
      unsigned int nb_decreasing_distance = 0;
      while ( ! SP.finished() )
        {
          last = std::get<0>( SP.current() );
          last_distance = std::get<2>( SP.current() );
          if ( V.count( last ) ) nb_multiple_pops += 1;
          V.insert( last );
          SP.expand();
          if ( last_distance < prev_distance ) nb_decreasing_distance += 1;
          prev_distance = last_distance;
        }
      // THEN( "No point is popped several times" )
      REQUIRE( nb_multiple_pops == 0 );
      // AND_THEN( "The sequence of popped points has non decreasing distances" )
      REQUIRE( nb_decreasing_distance == 0 );
      // AND_THEN( "The furthest point is also the antipodal point" )
      REQUIRE( last == uppest );
      // AND_THEN( "The furthest point is at distance close but lower than pi" )
      REQUIRE( last_distance*h >= 2.8 );
      REQUIRE( last_distance*h <= 3.14159265358979323844 );
      // Compute approximate shortest paths
      SP = TC.makeShortestPaths( 0 );
      SP.init( uppest ); //< set source
      double last_distance_opt = 0.0;
      while ( ! SP.finished() )
        {
          last = std::get<0>( SP.current() );
          last_distance_opt = std::get<2>( SP.current() );
          SP.expand();
        }
      // THEN( "The furthest point is also the antipodal point" )
      REQUIRE( last == lowest );
      // AND_THEN( "The furthest point is at distance close but lower than pi" )
      REQUIRE( last_distance_opt*h >= 2.8 );
      REQUIRE( last_distance_opt*h <= 3.14159265358979323844 );
      // AND_THEN( "This distance is greater or equal to the exacts shortest path" )
      REQUIRE( last_distance_opt*h >= last_distance*h );
    }
 
  SECTION( "Computing different shortest paths on a 3D unit sphere digitized at gridstep 0.125" )
    {
      // Make digital sphere
      const double h = 0.125;
      auto   params  = SH3::defaultParameters();
      params( "polynomial", "sphere1" )( "gridstep",  h );
      params( "minAABB", -2)( "maxAABB", 2)( "offset", 1.0 )( "closed", 1 );
      auto implicit_shape  = SH3::makeImplicitShape3D  ( params );
      auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
      auto K            = SH3::getKSpace( params );
      auto binary_image = SH3::makeBinaryImage(digitized_shape,
                                               SH3::Domain(K.lowerBound(),K.upperBound()),
                                               params );
      auto surface = SH3::makeDigitalSurface( binary_image, K, params );
      std::vector< Point >    lattice_points;
      auto pointels = SH3::getPointelRange( surface );
      for ( auto p : pointels ) lattice_points.push_back( K.uCoords( p ) );
      // Find lowest and uppest point.
      const Index nb = lattice_points.size();
      Index   lowest = 0;
      Index   uppest = 0;
      for ( Index i = 1; i < nb; i++ )
        {
          if ( lattice_points[ i ] < lattice_points[ lowest ] ) lowest = i;
          if ( lattice_points[ uppest ] < lattice_points[ i ] ) uppest = i;
        }
      // Compute shortest paths
      TangencyComputer< KSpace > TC( K );
      TC.init( lattice_points.cbegin(), lattice_points.cend() );
      auto SP = TC.makeShortestPaths( sqrt(3.0) );
 
      const double max_discrete_distance = 5.0;
      const int step = lattice_points.size() / 8;
      std::size_t sum_nb_visited = 0;
      std::size_t min_nb_visited = 100000;
      std::size_t max_nb_visited = 0;
      int n              = 0;
      for ( auto i = 0; i < lattice_points.size(); i += step, n += 1 )
        {
          SP.init( i );
          std::vector< std::size_t > visited;
          while ( ! SP.finished() )
            {
              auto last          = std::get<0>( SP.current() );
              auto last_distance = std::get<2>( SP.current() );
              visited.push_back( last );
              if ( last_distance >= max_discrete_distance ) break;
              SP.expand();
            }
          SP.clearVisited( visited );
          sum_nb_visited += visited.size();
          min_nb_visited  = std::min( min_nb_visited, visited.size() );
          max_nb_visited  = std::max( max_nb_visited, visited.size() );
        }
      double average = sum_nb_visited / (double) n;
      // Shortest paths from sources on a digital sphere should be of
      // approximately the same cardinal.
      REQUIRE( ( average - min_nb_visited ) / average < 0.2 );
      REQUIRE( ( max_nb_visited - average ) / average < 0.2 );
    }
}
 
SCENARIO( "TangencyComputer 3D tests", "[3d][tangency]" )
{
  typedef Z3i::Space          Space;
  typedef Z3i::KSpace         KSpace;
  typedef Shortcuts< KSpace > SH3;
  typedef Space::Point        Point;
  typedef std::size_t         Index;
 
  SECTION( "Computing shortest paths on a 3D unit sphere digitized at gridstep 0.125" )
    {
      // Make digital sphere
      const double h = 0.125;
      auto   params  = SH3::defaultParameters();
      params( "polynomial", "sphere1" )( "gridstep",  h );
      params( "minAABB", -2)( "maxAABB", 2)( "offset", 1.0 )( "closed", 1 );
      auto implicit_shape  = SH3::makeImplicitShape3D  ( params );
      auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
      auto K            = SH3::getKSpace( params );
      auto binary_image = SH3::makeBinaryImage(digitized_shape,
                                               SH3::Domain(K.lowerBound(),K.upperBound()),
                                               params );
      auto surface = SH3::makeDigitalSurface( binary_image, K, params );
      std::vector< Point >    lattice_points;
      auto pointels = SH3::getPointelRange( surface );
      for ( auto p : pointels ) lattice_points.push_back( K.uCoords( p ) );
      // Find lowest and uppest point.
      const Index nb = lattice_points.size();
      Index   lowest = 0;
      Index   uppest = 0;
      for ( Index i = 1; i < nb; i++ )
        {
          if ( lattice_points[ i ] < lattice_points[ lowest ] ) lowest = i;
          if ( lattice_points[ uppest ] < lattice_points[ i ] ) uppest = i;
        }
      // Compute shortest paths
      TangencyComputer< KSpace > TC( K );
      TC.init( lattice_points.cbegin(), lattice_points.cend() );
      const Point a = TC.point( lowest );
      std::vector< Index > V1 = TC.getCotangentPoints( a, 5.0 );
      std::vector< Index > V2 = TC.getCotangentPoints( a, 10.0 );
      std::vector< Index > V3 = TC.getCotangentPoints( a, 15.0 );
      std::sort( V1.begin(), V1.end() );
      std::sort( V2.begin(), V2.end() );
      std::sort( V3.begin(), V3.end() );
      REQUIRE( 0 < V1.size() );
      REQUIRE( V1.size() < V2.size() );
      REQUIRE( V2.size() <= V3.size() );
      REQUIRE( V3.size() < pointels.size()/2 );
      REQUIRE( std::includes( V2.begin(), V2.end(), V1.begin(), V1.end() ) );
      REQUIRE( std::includes( V3.begin(), V3.end(), V2.begin(), V2.end() ) );
      double    md1 = 0.0;
      double    md2 = 0.0;
      double    md3 = 0.0;
      for ( auto i : V1 ) md1 = std::max( md1, (TC.point( i ) - a).norm() );
      for ( auto i : V2 ) md2 = std::max( md2, (TC.point( i ) - a).norm() );
      for ( auto i : V3 ) md3 = std::max( md3, (TC.point( i ) - a).norm() );
      REQUIRE( md1 <= 5.0 );
      REQUIRE( md2 <= 10.0 );
      REQUIRE( md3 <= 15.0 );
    }
}