DGtal  0.9.2
AreaSurfaceEstimation-final.cpp
1 #include "DGtal/shapes/parametric/Ball3D.h"
4 
6 #include "DGtal/shapes/GaussDigitizer.h"
7 #include "DGtal/topology/LightImplicitDigitalSurface.h"
8 #include "DGtal/topology/DigitalSurface.h"
9 #include "DGtal/graph/DepthFirstVisitor.h"
10 #include "DGtal/graph/GraphVisitorRange.h"
11 
13 #include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h"
14 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h"
15 
16 using namespace DGtal;
17 
18 int main( int argc, char** argv )
19 {
20  const double h = 1;
21  const double radiusBall = 12.0;
22  const double radiusII = 6;
23  const double trueAreaSurface = 4.0*M_PI*radiusBall*radiusBall;
24 
25  trace.beginBlock( "Make parametric shape..." );
26 
27  typedef Ball3D< Z3i::Space > Shape;
28 
29  Z3i::RealPoint center( 0.0, 0.0, 0.0 );
30  Shape ball( center, radiusBall );
31 
32  trace.endBlock();
33 
34  trace.beginBlock( "Make digital shape..." );
35 
38 
39  DigitalShape digitalBall;
40  digitalBall.attach( ball );
41  digitalBall.init( ball.getLowerBound() - Z3i::RealPoint( 1.0, 1.0, 1.0 ),
42  ball.getUpperBound() + Z3i::RealPoint( 1.0, 1.0, 1.0 ),
43  h );
44  Domain domain = digitalBall.getDomain();
45  Z3i::KSpace kspace;
46  kspace.init( domain.lowerBound(), domain.upperBound(), true );
47 
48  trace.endBlock();
49 
50  trace.beginBlock( "Make digital surface..." );
51 
56  typedef GraphVisitorRange::ConstIterator SurfelConstIterator;
57  typedef Z3i::KSpace::Surfel Surfel;
58 
59  Surfel bel = Surfaces< Z3i::KSpace >::findABel( kspace, digitalBall, 500 );
60  SurfelAdjacency< Z3i::KSpace::dimension > surfelAdjacency( true );
61  LightDigitalSurface lightDigitalSurface( kspace, digitalBall, surfelAdjacency, bel );
62  DigitalSurface digitalSurface( lightDigitalSurface );
63 
64  GraphVisitorRange graphVisitorRange( new DepthFirstVisitor( digitalSurface, *digitalSurface.begin() ) );
65  SurfelConstIterator sbegin = graphVisitorRange.begin();
66  SurfelConstIterator send = graphVisitorRange.end();
67  std::vector< Surfel > v_border;
68  while( sbegin != send )
69  {
70  v_border.push_back( *sbegin );
71  ++sbegin;
72  }
73 
74  trace.endBlock();
75 
76  trace.beginBlock( "Computation with normal estimation ..." );
77 
78  typedef IIGeometricFunctors::IINormalDirectionFunctor< Z3i::Space > NormalFunctor;
80 
81  NormalFunctor normalFunctor;
82  IINormalEstimator normalEstimator( normalFunctor );
83  normalEstimator.attach( kspace, digitalBall );
84  normalEstimator.setParams( radiusII / h );
85  normalEstimator.init( h, v_border.begin(), v_border.end() );
86 
87  double areaSurfaceEstimated = 0.0;
88 
89  for( unsigned int i_position = 0; i_position < v_border.size(); ++i_position )
90  {
91  Z3i::RealPoint normalEstimated = normalEstimator.eval( &(v_border[i_position]) );
92  Z3i::RealPoint normalSurfel = kspace.sKCoords( kspace.sDirectIncident( v_border[i_position], kspace.sOrthDir( v_border[i_position] ))) - kspace.sKCoords( v_border[i_position] );
93  normalEstimated = normalEstimated.getNormalized();
94  areaSurfaceEstimated += std::abs( normalEstimated.dot( normalSurfel )) * h * h;
95  }
96 
97  trace.endBlock();
98 
99  trace.info() << "Area Surface estimated : " << areaSurfaceEstimated << std::endl;
100  trace.info() << "True areaSurface : " << trueAreaSurface << std::endl;
101  trace.info() << "Ratio : " << areaSurfaceEstimated / trueAreaSurface << std::endl;
102 
103  return 0;
104 }
void beginBlock(const std::string &keyword="")
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of an impl...
Trace trace
Definition: Common.h:130
const Point & sKCoords(const SCell &c) const
Aim: This class is useful to perform a depth-first exploration of a graph given a starting point or s...
Aim: Represent adjacencies between surfel elements, telling if it follows an interior to exterior ord...
double endBlock()
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
Aim: This class implement an Integral Invariant estimator which computes for each surfel the covarian...
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:141
SCell sDirectIncident(const SCell &p, Dimension k) const
static SCell findABel(const KSpace &K, const PointPredicate &pp, unsigned int nbtries=1000)
void attach(const EuclideanShape &shape)
void init(const RealPoint &xLow, const RealPoint &xUp, typename RealVector::Component gridStep)
Dimension sOrthDir(const SCell &s) const
bool init(const Point &lower, const Point &upper, bool isClosed)
Aim: A class for computing the Gauss digitization of some Euclidean shape, i.e. its intersection with...
Domain getDomain() const
const Point & upperBound() const
DGtal is the top-level namespace which contains all DGtal functions and types.
Aim: Transforms a graph visitor into a single pass input range.
std::ostream & info()
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value...
const Point & lowerBound() const
Space::RealPoint RealPoint
Definition: StdDefs.h:170
Component dot(const Self &v) const
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex...
PointVector< dim, double, std::array< double, dim > > getNormalized() const
Aim: Model of the concept StarShaped3D represents any Sphere in the space.
Definition: Ball3D.h:60