File failed to load: https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.3/config/TeX-MML-AM_CHTML/MathJax.js
DGtal 2.0.0
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/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;
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, 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 >
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 )
{
auto prev_geometry
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 ... ";
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 =
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;
}
// //
Structure representing an RGB triple with alpha component.
Definition Color.h:77
PointVector< dim, Integer > Point
static const constexpr Dimension dimension
void show() override
Starts the event loop and display of elements.
CountedPtr< SH3::DigitalSurface > surface
Z3i::SCell SCell
NeighborhoodConvexityAnalyzer< KSpace, 1 > NCA
Z3i this namespace gathers the standard of types for 3D imagery.
KhalimskySpaceND< 3, Integer > KSpace
Definition StdDefs.h:146
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition Common.h:119
Trace trace
STL namespace.
static std::vector< Point > debug_one(const KSpace &aK, Point p, ImagePtr bimage)
static std::vector< int > run(const KSpace &aK, const std::vector< Point > &pts, ImagePtr bimage)
static std::vector< Geometry > multiscale_run(const KSpace &aK, const std::vector< Point > &pts, ImagePtr bimage)
std::pair< int, int > Geometry
Shortcuts< KSpace > SH3
int main()
Definition testBits.cpp:56
MyPointD Point
KSpace K
Image image(domain)