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.

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 ... ";
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 );
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 );
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;
}
// //
void beginBlock(const std::string &keyword="")
std::ostream & info()
void progressBar(const double currentValue, const double maximalValue)
double endBlock()
CountedPtr< SH3::DigitalSurface > surface
Z3i::SCell SCell
NeighborhoodConvexityAnalyzer< KSpace, 1 > NCA
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition: Common.h:136
Trace trace
Definition: Common.h:153
int main(int argc, char **argv)
Shortcuts< KSpace > SH3
int max(int a, int b)
MyPointD Point
Definition: testClone2.cpp:383
KSpace K