fullConvexityAnalysis3D.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/shapes/Shapes.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/helpers/Shortcuts.h"
#include "DGtal/images/ImageContainerBySTLVector.h"
#include "DGtal/geometry/volumes/NeighborhoodConvexityAnalyzer.h"

 
using namespace std;
using namespace DGtal;
using namespace Z3i;
typedef Shortcuts< KSpace > SH3;

 
 
template < typename KSpace, int N >
struct Analyzer {
  typedef typename KSpace::Point Point;
  typedef NeighborhoodConvexityAnalyzer< KSpace, N > NCA;
 
  template < typename ImagePtr >
  static
  std::vector<Point>
  debug_one( const KSpace& aK, Point p, ImagePtr bimage )
  {
    NCA nca( aK.lowerBound(), aK.upperBound(),
             KSpace::dimension <= 2 ? 0 : 10000*KSpace::dimension*N );
    auto& image = *bimage;
    int    geom = 0;
    nca.setCenter( p, image );
    bool  cvx = nca.isFullyConvex( true );
    bool ccvx = nca.isComplementaryFullyConvex( false );
    auto cfg  = nca.makeConfiguration( nca.configuration(), true, false );
    std::vector< Point > localCompX;
    nca.getLocalCompX( localCompX, false );
    std::cout << "InC=" << nca.configuration() << std::endl;
    std::cout << "Cfg=" << cfg << std::endl;
    for ( auto q : localCompX ) std::cout << q;
    std::cout << std::endl;
    geom = ( cvx ? 0x1 : 0x0 ) | ( ccvx ? 0x2 : 0x0 );
    std::cout << "cvx=" << cvx << " ccvx=" << ccvx << std::endl;
    std::cout << "geom=" << geom << std::endl;
    return localCompX;
  }
 
  template < typename ImagePtr >
  static
  std::vector<int>
  run( const KSpace& aK, const std::vector<Point>& pts, ImagePtr bimage )
  {
    NCA nca( aK.lowerBound(), aK.upperBound(), 0 );
    // KSpace::dimension <= 2 ? 0 : 10000*KSpace::dimension*N );
    auto& image = *bimage;
    std::vector<int> result;
    std::map< Point, int > computed;
    int geom;
    int i  = 0;
    int nb = pts.size();
    int  nb_cvx = 0;
    int nb_ccvx = 0;
    for ( auto p : pts )
      {
        if ( i % 100 == 0 ) trace.progressBar( i, nb );
        auto it = computed.find( p );
        if ( it == computed.end() )
          {
            nca.setCenter( p, image );
            bool  cvx = nca.isFullyConvex( true );
            bool ccvx = nca.isComplementaryFullyConvex( false );
            if ( cvx  ) nb_cvx  += 1;
            if ( ccvx ) nb_ccvx += 1;
            geom = ( cvx ? 0x1 : 0x0 ) | ( ccvx ? 0x2 : 0x0 );
            computed[ p ] = geom;
          }
        else geom = it->second;
        result.push_back( geom );
        i++; 
      }
    trace.info() << "nb_cvx=" << nb_cvx << " nb_ccvx=" << nb_ccvx << std::endl;
    return result;    
  }
 
  template < typename ImagePtr >
  static
  void
  run( std::vector<int> & to_update,
       const KSpace& aK, const std::vector<Point>& pts, ImagePtr bimage )
  {
    NCA nca( aK.lowerBound(), aK.upperBound() );
    // KSpace::dimension <= 2 ? 0 : 10000*KSpace::dimension*N );
    auto& image = *bimage;
    std::map< Point, int > computed;
    int geom;
    int i  = 0;
    int nb = pts.size();
    for ( auto p : pts )
      {
        if ( i % 100 == 0 ) trace.progressBar( i, nb );
        auto it = computed.find( p );
        if ( it == computed.end() )
          {
            nca.setCenter( p, image );
            bool  cvx = ( to_update[ i ] & 0x1 )
              ? nca.isFullyConvex( true )
              : false;
            bool ccvx = ( to_update[ i ] & 0x2 )
              ? nca.isComplementaryFullyConvex( false )
              : false;
            geom = ( cvx ? 0x1 : 0x0 ) | ( ccvx ? 0x2 : 0x0 );
            computed[ p ] = geom;
          }
        else geom = it->second;
        to_update[ i++ ] = geom;
      }
  }
};
 
template < typename KSpace, int N >
struct MultiScaleAnalyzer {
  typedef typename KSpace::Point Point;
  typedef NeighborhoodConvexityAnalyzer< KSpace, N > NCA;
  typedef std::pair< int, int > Geometry; // convex, concave
 
  template < typename ImagePtr >
  static
  std::vector< Geometry >
  multiscale_run( const KSpace& aK,
                  const std::vector<Point>& pts,
                  ImagePtr bimage )
  {
    auto prev_geometry
      = MultiScaleAnalyzer< KSpace, N-1>::multiscale_run( aK, pts, bimage );
    trace.info() << "------- Analyzing scale " << N << " --------" << std::endl;
    std::vector< int > geom( prev_geometry.size() );
    for ( int i = 0; i < geom.size(); i++ )
      geom[ i ] = ( prev_geometry[ i ].first == N-1 ? 0x1 : 0x0 )
        |  ( prev_geometry[ i ].second == N-1 ? 0x2 : 0x0 );
    Analyzer< KSpace, N>::run( geom, aK, pts, bimage );
    for ( int i = 0; i < geom.size(); i++ ) {
      prev_geometry[ i ].first  += ( geom[ i ] & 0x1 ) ? 1 : 0;
      prev_geometry[ i ].second += ( geom[ i ] & 0x2 ) ? 1 : 0;
    }
    return prev_geometry;
  }
};

template < typename KSpace>
struct MultiScaleAnalyzer< KSpace, 0 > {
  typedef typename KSpace::Point Point;
  typedef std::pair< int, int > Geometry; // convex, concave
 
  template < typename ImagePtr >
  static
  std::vector< Geometry >
  multiscale_run( const KSpace& aK,
                  const std::vector<Point>& pts,
                  ImagePtr bimage )
  {
    ((void) aK);
    ((void) bimage);
    return std::vector< Geometry >( pts.size(), std::make_pair( 0, 0 ) );
  }
};
 
int main( int argc, char** argv )
{
  if ( argc <= 2 )
    {
      trace.info() << "Usage: " << argv[ 0 ] << " <K> <input.vol> <m> <M>" << std::endl;
      trace.info() << "\tAnalyze the  shape with local full convexity" << std::endl;
      trace.info() << "\t- 1 <= K <= 5: analysis at scale K" << std::endl;
      trace.info() << "\t- K == 0: multiscale analysis (using scales 1-5)" << 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;
      return 1;
    }
  int         N = atoi( argv[ 1 ] );
  std::string fn= argv[ 2 ];
  int         m = argc > 3 ? atoi( argv[ 3 ] ) : 0;
  int         M = argc > 4 ? atoi( argv[ 4 ] ) : 255;
 
  auto   params  = SH3::defaultParameters();
  
  // Domain creation from two bounding points.
  trace.info() << "Building set or importing vol ... ";
  KSpace K;
  params( "thresholdMin", m );
  params( "thresholdMax", M );
  auto bimage = SH3::makeBinaryImage( fn, params );
  K = SH3::getKSpace( bimage );
  trace.info() << "  [Done]" << std::endl;
  // Compute surface
  params( "surfaceComponents" , "All" );
  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;
  if ( N != 0 )
    {
      std::vector< int > result;
      trace.beginBlock ( "Single scale analysis" );
      if ( N == 1 ) result = Analyzer< KSpace, 1 >::run( K, points, bimage );
      if ( N == 2 ) result = Analyzer< KSpace, 2 >::run( K, points, bimage );
      if ( N == 3 ) result = Analyzer< KSpace, 3 >::run( K, points, bimage );
      if ( N == 4 ) result = Analyzer< KSpace, 4 >::run( K, points, bimage );
      if ( N == 5 ) result = Analyzer< KSpace, 5 >::run( K, points, bimage );
      trace.endBlock();
      SCell dummy;
      Color colors[ 4 ] =
        { Color( 255, 0, 0, 255 ), Color( 0, 255, 0, 255 ),
          Color( 0, 0, 255, 255 ), Color( 255, 255, 255, 255 ) };
      auto surfels   = SH3::getSurfelRange( surface, params );
      SH3::Colors all_colors( surfels.size() );
      for ( int i = 0; i < surfels.size(); i++ )
        {
          const auto    j = surfel2idx[ surfels[ i ] ];
          all_colors[ i ] = colors[ result[ j ] ];
        }
      SH3::saveOBJ( surface, SH3::RealVectors(), all_colors, "geom-cvx.obj" );
      PolyscopeViewer viewer;
      int i = 0;
      for ( auto s : (*surface) )
        {
          viewer << all_colors[ i ]
                 << s;
          i++;
        }
      viewer.show();  
    }
  else
    {
      trace.beginBlock ( "Multiscale analysis" );
      auto geometry =
        MultiScaleAnalyzer< KSpace, 5 >::multiscale_run( K, points, bimage );
      trace.endBlock();
      Color colors_planar[ 6 ] =
        { Color( 0, 255, 255, 255),
          Color( 50, 255, 255, 255), Color( 100, 255, 255, 255),
          Color( 150, 255, 255, 255), Color( 200, 255, 255, 255 ),
          Color( 255, 255, 255, 255 ) };
      Color color_atypical( 255, 0, 0, 255 ); 
      Color colors_cvx[ 5 ] =
        { Color( 0, 255, 0, 255 ), Color( 50, 255, 50, 255 ),
          Color( 100, 255, 100, 255 ), Color( 150, 255, 150, 255 ),
          Color( 200, 255, 200, 255 ) };
      Color colors_ccv[ 5 ] =
        { Color( 0, 0, 255, 255 ), Color( 50, 50, 255, 255 ),
          Color( 100, 100, 255, 255 ), Color( 150, 150, 255, 255 ),
          Color( 200, 200, 255, 255 ) };
      auto surfels   = SH3::getSurfelRange( surface, params );
      SH3::Colors all_colors( surfels.size() );
      for ( int i = 0; i < surfels.size(); i++ ) {
        const auto j = surfel2idx[ surfels[ i ] ];
        int m0 = std::min( geometry[ j ].first, geometry[ j ].second );
        int m1 = std::max( geometry[ j ].first, geometry[ j ].second );
        if ( m1 == 0 ) all_colors[ i ] = color_atypical;
        else if ( m0 == m1 ) all_colors[ i ] = colors_planar[ 5 ];
        else if ( geometry[ j ].first > geometry[ j ].second )
          all_colors[ i ] = colors_cvx[ 5 - abs( m0 - m1 ) ];
        else
          all_colors[ i ] = colors_ccv[ 5 - abs( m0 - m1 ) ];
      }
      SH3::saveOBJ( surface, SH3::RealVectors(), all_colors, "geom-scale-cvx.obj" );
      SCell dummy;
      int i = 0;
      PolyscopeViewer<> viewer;
      for ( auto s : (*surface) )
        {
          viewer << all_colors[ i ]
                 << s;
          i++;
        }
      viewer.show();  
    }      
  return 0;
}
//                                                                           //