DGtal  1.3.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 [74] 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;
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 );
// 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 )
: 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; // 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;
}
// //
DGtal::KhalimskySpaceND::upperBound
const Point & upperBound() const
Return the upper bound for digital points in this space.
DGtal::Trace::endBlock
double endBlock()
SH3
Shortcuts< KSpace > SH3
Definition: testArithmeticalDSSComputerOnSurfels.cpp:49
max
int max(int a, int b)
Definition: testArithmeticalDSS.cpp:1108
DGtal::NeighborhoodConvexityAnalyzer::isComplementaryFullyConvex
bool isComplementaryFullyConvex(bool with_center)
Definition: NeighborhoodConvexityAnalyzer.h:355
DGtal::NeighborhoodConvexityAnalyzer::configuration
Configuration configuration() const
Definition: NeighborhoodConvexityAnalyzer.h:231
DGtal::Color
Structure representing an RGB triple with alpha component.
Definition: Color.h:67
DGtal::trace
Trace trace
Definition: Common.h:154
DGtal::Shortcuts::Colors
std::vector< Color > Colors
Definition: Shortcuts.h:192
K
KSpace K
Definition: testCubicalComplex.cpp:62
DGtal::Dimension
DGtal::uint32_t Dimension
Definition: Common.h:137
DGtal::Shortcuts::makeDigitalSurface
static CountedPtr< DigitalSurface > makeDigitalSurface(CountedPtr< TPointPredicate > bimage, const KSpace &K, const Parameters &params=parametersDigitalSurface())
Definition: Shortcuts.h:1209
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::KhalimskySpaceND::lowerBound
const Point & lowerBound() const
Return the lower bound for digital points in this space.
DGtal::Shortcuts
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:104
DGtal::KhalimskySpaceND::sCoords
Point sCoords(const SCell &c) const
Return its digital coordinates.
DGtal::SignedKhalimskyCell
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
Definition: KhalimskySpaceND.h:208
DGtal::Shortcuts::RealVectors
std::vector< RealVector > RealVectors
Definition: Shortcuts.h:179
DGtal::KhalimskySpaceND::sDirect
bool sDirect(const SCell &p, Dimension k) const
Return 'true' if the direct orientation of [p] along [k] is in the positive coordinate direction.
DGtal::Trace::info
std::ostream & info()
DGtal::KhalimskySpaceND::dimension
static const constexpr Dimension dimension
Definition: KhalimskySpaceND.h:430
DGtal::Trace::progressBar
void progressBar(const double currentValue, const double maximalValue)
DGtal::SignedKhalimskyCell::className
std::string className() const
Return the style name used for drawing this object.
DGtal::Viewer3D
Definition: Viewer3D.h:135
DGtal::Shortcuts::defaultParameters
static Parameters defaultParameters()
Definition: Shortcuts.h:203
DGtal::NeighborhoodConvexityAnalyzer::getLocalCompX
void getLocalCompX(std::vector< Point > &localCompX, bool with_center) const
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::NeighborhoodConvexityAnalyzer::setCenter
void setCenter(Point c, const PointPredicate &X)
DGtal::CustomColors3D
Definition: DrawWithDisplay3DModifier.h:129
DGtal::NeighborhoodConvexityAnalyzer::isFullyConvex
bool isFullyConvex(bool with_center)
Definition: NeighborhoodConvexityAnalyzer.h:288
DGtal::SetMode3D
Modifier class in a Display3D stream. Useful to choose your own mode for a given class....
Definition: DrawWithDisplay3DModifier.h:73
DGtal::NeighborhoodConvexityAnalyzer
Aim: A class that models a neighborhood and that provides services to analyse the convexity properti...
Definition: NeighborhoodConvexityAnalyzer.h:94
main
int main(int argc, char **argv)
Definition: testArithmeticDSS-benchmark.cpp:147
DGtal::functions::abs
T abs(const T &a)
Definition: BasicMathFunctions.h:116
DGtal::KhalimskySpaceND::sOrthDir
Dimension sOrthDir(const SCell &s) const
Given a signed surfel [s], returns its orthogonal direction (ie, the coordinate where the surfel is c...
DGtal::Shortcuts::getKSpace
static KSpace getKSpace(const Point &low, const Point &up, Parameters params=parametersKSpace())
Definition: Shortcuts.h:332
DGtal::NeighborhoodConvexityAnalyzer::makeConfiguration
static Configuration makeConfiguration(Configuration current, bool complement, bool with_center)
Definition: NeighborhoodConvexityAnalyzer.h:473
DGtal::PointVector< dim, Integer >
DGtal::Viewer3D::show
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
DGtal::KhalimskySpaceND::sIncident
SCell sIncident(const SCell &c, Dimension k, bool up) const
Return the forward or backward signed cell incident to [c] along axis [k], depending on [up].
DGtal::Shortcuts::saveOBJ
static bool saveOBJ(CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer > > digsurf, const TCellEmbedder &embedder, const RealVectors &normals, const Colors &diffuse_colors, std::string objfile, const Color &ambient_color=Color(32, 32, 32), const Color &diffuse_color=Color(200, 200, 255), const Color &specular_color=Color::White)
Definition: Shortcuts.h:1739
DGtal::Shortcuts::getSurfelRange
static SurfelRange getSurfelRange(CountedPtr< ::DGtal::DigitalSurface< TDigitalSurfaceContainer > > surface, const Parameters &params=parametersDigitalSurface())
Definition: Shortcuts.h:1547
DGtal::Shortcuts::makeBinaryImage
static CountedPtr< BinaryImage > makeBinaryImage(Domain shapeDomain)
Definition: Shortcuts.h:561
Point
MyPointD Point
Definition: testClone2.cpp:383
NCA
NeighborhoodConvexityAnalyzer< KSpace, 1 > NCA
Definition: fullConvexityThinning3D.cpp:60
DGtal::KhalimskySpaceND
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Definition: KhalimskySpaceND.h:64