DGtal  1.5.beta
geometry/volumes/fullConvexityAnalysis3D.cpp

This example shows how to analyze the local geometry of 3D digital sets with full convexity over cubical neighborhoods.

See also
Local convexity analysis

For instance, you may call it to analyse image Al.100.vol at scale 2 as

fullConvexityAnalysis3D 2 ${DGTAL}/examples/samples/Al.100.vol

Results are displayed, then saved in 'geom-cvx.obj'. You may also analyse the same shape in multiscale fashion with

fullConvexityAnalysis3D 2 ${DGTAL}/examples/samples/Al.100.vol

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 [79] for details).

Full convexity analysis at scale 1
Full convexity analysis at scale 2
Full convexity analysis at scale 3
Full convexity smooth multiscale analysis (scales 1-5)
#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 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, 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, 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,
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,
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();
// 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" );
Viewer3D<> viewer;
viewer.setWindowTitle("fullConvexityAnalysis3D");
viewer.show();
int i = 0;
viewer << SetMode3D( dummy.className(), "Basic" );
for ( auto s : (*surface) )
{
viewer << CustomColors3D( all_colors[ i ], all_colors[ i ] )
<< s;
i++;
}
viewer<< Viewer3D<>::updateDisplay;
application.exec();
}
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;
Viewer3D<> viewer;
viewer.setWindowTitle("fullConvexityAnalysis3D");
viewer.show();
viewer << SetMode3D( dummy.className(), "Basic" );
for ( auto s : (*surface) )
{
viewer << CustomColors3D( all_colors[ i ], all_colors[ i ] )
<< s;
i++;
}
viewer<< Viewer3D<>::updateDisplay;
application.exec();
}
return 0;
}
// //
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::Trace::info
std::ostream & info()
DGtal::Trace::progressBar
void progressBar(const double currentValue, const double maximalValue)
DGtal::Trace::endBlock
double endBlock()
surface
CountedPtr< SH3::DigitalSurface > surface
Definition: dgtalCalculus-geodesic.cpp:66
SCell
Z3i::SCell SCell
Definition: digitalPolyhedronBuilder3D.cpp:79
NCA
NeighborhoodConvexityAnalyzer< KSpace, 1 > NCA
Definition: fullConvexityThinning3D.cpp:60
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::Dimension
DGtal::uint32_t Dimension
Definition: Common.h:136
DGtal::trace
Trace trace
Definition: Common.h:153
main
int main(int argc, char **argv)
Definition: testArithmeticDSS-benchmark.cpp:147
SH3
Shortcuts< KSpace > SH3
Definition: testArithmeticalDSSComputerOnSurfels.cpp:49
max
int max(int a, int b)
Definition: testArithmeticalDSS.cpp:1108
Point
MyPointD Point
Definition: testClone2.cpp:383
K
KSpace K
Definition: testCubicalComplex.cpp:62