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