DGtal  1.0.0
Functions
area-estimation-with-indexed-digital-surface.cpp File Reference
#include <cstdlib>
#include <iostream>
#include <map>
#include "DGtal/base/Common.h"
#include "DGtal/base/CountedPtr.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/kernel/PointVector.h"
#include "DGtal/graph/CUndirectedSimpleGraph.h"
#include "DGtal/graph/BreadthFirstVisitor.h"
#include "DGtal/topology/DigitalSetBoundary.h"
#include "DGtal/topology/IndexedDigitalSurface.h"
#include "DGtal/shapes/Shapes.h"
Include dependency graph for area-estimation-with-indexed-digital-surface.cpp:

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

This program is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.

Author
Jacques-Olivier Lachaud (jacqu.nosp@m.es-o.nosp@m.livie.nosp@m.r.la.nosp@m.chaud.nosp@m.@uni.nosp@m.v-sav.nosp@m.oie..nosp@m.fr ) Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
Date
2015/08/28

An example file named area-estimation-with-indexed-digital-surface.cpp.

This file is part of the DGtal library.

Definition in file area-estimation-with-indexed-digital-surface.cpp.

Function Documentation

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 83 of file area-estimation-with-indexed-digital-surface.cpp.

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;
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  auto dsc = CountedPtr< DigitalSurfaceContainer >( new DigitalSurfaceContainer( K, aSet ) );
101  DigSurface dsurf( dsc );
102  trace.info() << "[OK]"
103  << " Sphere has " << dsurf.size() << " vertices/surfels"
104  << " with euler " << dsurf.Euler() << std::endl;
105  trace.endBlock();
106 
107  trace.beginBlock( "Estimating normals" );
108  auto v2n = dsurf.makeVertexMap( RealPoint() );
109  for ( auto v : dsurf )
110  {
111  int nbv = 0;
112  RealPoint normal = RealPoint::zero;
113  BFSVisitor bfv( dsurf, v );
114  while( ! bfv.finished() )
115  { // Vertex are colored according to the distance to initial seed.
116  auto node = bfv.current();
117  if ( KN < node.second ) break;
118  auto surfel = dsurf.surfel( node.first );
120  Dimension k = K.sOrthDir( surfel );
121  nv[ k ] = K.sDirect( surfel, k ) ? 1.0 : -1.0;
122  normal += nv;
123  nbv += 1;
124  bfv.expand();
125  }
126  normal /= nbv;
127  v2n[ v ] = normal.getNormalized();
128  }
129  trace.endBlock();
130 
131  trace.beginBlock( "Estimating area" );
132  const double area_true = 4.0 * M_PI * R * R;
133  double area_averaged = 0.0;
134  double area_corrected = 0.0;
135  for ( auto v : dsurf )
136  {
137  auto surfel = dsurf.surfel( v );
138  Dimension k = K.sOrthDir( surfel );
139  area_corrected += fabs( v2n[ v ][ k ] );
140  area_averaged += 1.0 / v2n[ v ].norm( RealPoint::L_1 );
141  }
142  trace.info() << "- true area = " << area_true << std::endl;
143  trace.info() << "- corrected area = " << area_corrected << std::endl;
144  trace.info() << "- averaged area = " << area_averaged << std::endl;
145  trace.endBlock();
146 
147  return 0;
148 }
void beginBlock(const std::string &keyword="")
HyperRectDomain< Space > Domain
Trace trace
Definition: Common.h:144
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::uint32_t Dimension
Definition: Common.h:127
PointVector< dim, double, std::array< double, dim > > getNormalized() const
double endBlock()
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of a given...
PointVector< 3, double > RealPoint
Dimension sOrthDir(const SCell &s) const
Given a signed surfel [s], returns its orthogonal direction (ie, the coordinate where the surfel is c...
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
MyPointD Point
Definition: testClone2.cpp:383
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
std::ostream & info()
Aim: Represents a digital surface with the topology of its dual surface. Its aim is to mimick the sta...
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:1582
KSpace K
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,...

References DGtal::Trace::beginBlock(), DGtal::Trace::endBlock(), DGtal::PointVector< dim, TEuclideanRing, TContainer >::getNormalized(), DGtal::Trace::info(), DGtal::KhalimskySpaceND< dim, TInteger >::init(), K, DGtal::PointVector< 3, double >::L_1, DGtal::KhalimskySpaceND< dim, TInteger >::sDirect(), DGtal::KhalimskySpaceND< dim, TInteger >::sOrthDir(), DGtal::trace, and DGtal::PointVector< 3, double >::zero.