DGtal  0.9.4beta
area-estimation-with-digital-surface.cpp
1 
63 #include <cstdlib>
65 #include <iostream>
66 #include <map>
67 #include "DGtal/base/Common.h"
68 #include "DGtal/base/CountedPtr.h"
69 #include "DGtal/helpers/StdDefs.h"
70 #include "DGtal/kernel/PointVector.h"
71 #include "DGtal/graph/CUndirectedSimpleGraph.h"
72 #include "DGtal/graph/BreadthFirstVisitor.h"
73 #include "DGtal/topology/DigitalSetBoundary.h"
74 #include "DGtal/topology/DigitalSurface.h"
75 #include "DGtal/shapes/Shapes.h"
76 
78 
79 using namespace std;
80 using namespace DGtal;
81 using namespace DGtal::Z3i;
82 
83 int main( int argc, char** argv )
84 {
85  const double R = argc >= 2 ? atof( argv[ 1 ] ) : 13.0; // radius of ball
86  const unsigned int KN = argc >= 3 ? atoi( argv[ 2 ] ) : 6; // size of neighborhood
87  const int M = (int) ceil( R + 2.0 );
88  trace.beginBlock( "Creating surface" );
89  trace.info() << "Sphere of radius " << R
90  << ", " << KN << "-ring neighborhood" << std::endl;
91  typedef DigitalSetBoundary < KSpace, DigitalSet > DigitalSurfaceContainer;
92  typedef DigitalSurface < DigitalSurfaceContainer > DigSurface;
93  typedef BreadthFirstVisitor < DigSurface > BFSVisitor;
94  Point p1( -M, -M, -M );
95  Point p2( M, M, M );
96  KSpace K;
97  K.init( p1, p2, true );
98  DigitalSet aSet( Domain( p1, p2 ) );
99  Shapes<Domain>::addNorm2Ball( aSet, Point( 0, 0, 0 ), R );
100  DigSurface dsurf( new DigitalSurfaceContainer( K, aSet ) );
101  trace.info() << "[OK]" << " Sphere has " << dsurf.size() << " vertices/surfels"
102  << std::endl;
103  trace.endBlock();
104 
105  trace.beginBlock( "Estimating normals" );
106  std::map< SCell, RealPoint > v2n;
107  for ( auto v : dsurf )
108  {
109  int nbv = 0;
110  RealPoint normal = RealPoint::zero;
111  BFSVisitor bfv( dsurf, v );
112  while( ! bfv.finished() )
113  { // Vertex are colored according to the distance to initial seed.
114  auto node = bfv.current();
115  if ( KN < node.second ) break;
116  auto surfel = node.first;
118  Dimension k = K.sOrthDir( surfel );
119  nv[ k ] = K.sDirect( surfel, k ) ? 1.0 : -1.0;
120  normal += nv;
121  nbv += 1;
122  bfv.expand();
123  }
124  normal /= nbv;
125  v2n[ v ] = normal.getNormalized();
126  }
127  trace.endBlock();
128 
129  trace.beginBlock( "Estimating area" );
130  const double area_true = 4.0 * M_PI * R * R;
131  double area_averaged = 0.0;
132  double area_corrected = 0.0;
133  for ( auto v : dsurf )
134  {
135  auto surfel = v;
136  Dimension k = K.sOrthDir( surfel );
137  area_corrected += fabs( v2n[ v ][ k ] );
138  area_averaged += 1.0 / v2n[ v ].norm( RealPoint::L_1 );
139  }
140  trace.info() << "- true area = " << area_true << std::endl;
141  trace.info() << "- corrected area = " << area_corrected << std::endl;
142  trace.info() << "- averaged area = " << area_averaged << std::endl;
143  trace.endBlock();
144 
145  return 0;
146 }
void beginBlock(const std::string &keyword="")
Z3i this namespace gathers the standard of types for 3D imagery.
Trace trace
Definition: Common.h:137
DGtal::uint32_t Dimension
Definition: Common.h:120
STL namespace.
double endBlock()
Dimension sOrthDir(const SCell &s) const
bool init(const Point &lower, const Point &upper, bool isClosed)
DGtal is the top-level namespace which contains all DGtal functions and types.
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
std::ostream & info()
Aim: This class is useful to perform a breadth-first exploration of a graph given a starting point or...
static Self zero
Static const for zero PointVector.
Definition: PointVector.h:798
Aim: A utility class for constructing different shapes (balls, diamonds, and others).
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex...
bool sDirect(const SCell &p, Dimension k) const
PointVector< dim, double, std::array< double, dim > > getNormalized() const