DGtal 1.3.0
Searching...
No Matches
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 [75] 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 )
? 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; // 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;
}
// //
Structure representing an RGB triple with alpha component.
Definition: Color.h:68
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...
Definition: Shortcuts.h:105
void beginBlock(const std::string &keyword="")
std::ostream & info()
void progressBar(const double currentValue, const double maximalValue)
double endBlock()
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
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:137
Trace trace
Definition: Common.h:154
STL namespace.
Modifier class in a Display3D stream. Useful to choose your own mode for a given class....
Shortcuts< KSpace > SH3
int main()
Definition: testBits.cpp:56
MyPointD Point
Definition: testClone2.cpp:383
KSpace K