This example shows how to analyze the local geometry of 3D digital sets with full convexity over cubical neighborhoods.
Results are displayed, then saved in 'geom-cvx.obj'. You may also analyse the same shape in multiscale fashion with
The result is saved in 'geom-scale-cvx.obj'. You will obtain images like below, where green means convex, blue means concave, white is planar and red is atypical (see [75] for details).
#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/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 Z3i;
template < typename KSpace, int N >
struct Analyzer {
typedef typename KSpace::Point
Point;
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, std::vector<Point> pts, ImagePtr bimage )
{
NCA nca( aK.lowerBound(), aK.upperBound(), 0 );
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 )
{
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, std::vector<Point> pts, ImagePtr bimage )
{
NCA nca( aK.lowerBound(), aK.upperBound() );
auto& image = *bimage;
std::map< Point, int > computed;
int geom;
int i = 0;
int nb = pts.size();
for ( auto p : pts )
{
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 std::pair< int, int > Geometry;
template < typename ImagePtr >
static
std::vector< Geometry >
multiscale_run(
const KSpace& aK,
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;
template < typename ImagePtr >
static
std::vector< Geometry >
multiscale_run(
const KSpace& aK,
std::vector<Point> pts,
ImagePtr 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;
QApplication application(argc,argv);
auto params = SH3::defaultParameters();
trace.
info() <<
"Building set or importing vol ... ";
params( "thresholdMin", m );
params( "thresholdMax", M );
auto bimage = SH3::makeBinaryImage( fn, params );
K = SH3::getKSpace( bimage );
params( "surfaceComponents" , "All" );
auto surface = SH3::makeDigitalSurface( bimage,
K, params );
std::vector< Point > points;
std::map< SCell, int > surfel2idx;
std::map< Point, int > point2idx;
int idx = 0;
for ( auto s : (*surface) )
{
auto voxel =
K.sIncident( s, k,
K.sDirect( s, k ) );
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;
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 );
{
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" );
viewer.setWindowTitle("fullConvexityAnalysis3D");
int i = 0;
viewer <<
SetMode3D( dummy.className(),
"Basic" );
for ( auto s : (*surface) )
{
<< s;
i++;
}
viewer<< Viewer3D<>::updateDisplay;
application.exec();
}
else
{
auto geometry =
MultiScaleAnalyzer< KSpace, 5 >::multiscale_run(
K, points, bimage );
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( 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( 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" );
int i = 0;
viewer.setWindowTitle("fullConvexityAnalysis3D");
viewer <<
SetMode3D( dummy.className(),
"Basic" );
for ( auto s : (*surface) )
{
<< s;
i++;
}
viewer<< Viewer3D<>::updateDisplay;
application.exec();
}
return 0;
}
Structure representing an RGB triple with alpha component.
Aim: A class that models a neighborhood and that provides services to analyse the convexity properti...
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
void beginBlock(const std::string &keyword="")
void progressBar(const double currentValue, const double maximalValue)
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
NeighborhoodConvexityAnalyzer< KSpace, 1 > NCA
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Modifier class in a Display3D stream. Useful to choose your own mode for a given class....