DGtal  1.3.beta
topology/area-estimation-with-indexed-digital-surface.cpp

Computing the area of a sphere with an indexed digital surface. First normals are estimated by averaging the trivial normals in a k-neighborhood of each surfel. Then both the corrected area (scalar product of estimated normal with trivial normal) and the averaged area (inverse of L1-norm of estimated normal) is computed per surfel, and summed in order to get area estimation. This example illustrates the use of a DGtal::IndexedDigitalSurface and DGtal::BreadthFirstVisitor onto it. This method is slightly faster than using a DGtal::DigitalSurface.

Precomputed 3D digital surface with IndexedDigitalSurface
\$ ./examples/topology/area-estimation-with-indexed-digital-surface 13 6
New Block [Creating surface]
Sphere of radius 13, 6-ring neighborhood
[OK] Sphere has 3174 vertices/surfels with euler 2
EndBlock [Creating surface] (55.999 ms)
New Block [Estimating normals]
EndBlock [Estimating normals] (274.712 ms)
New Block [Estimating area]
- true area      = 2123.72
- corrected area = 2109.58
- averaged area  = 2150.34
EndBlock [Estimating area] (0.023 ms)

#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/topology/DigitalSetBoundary.h"
#include "DGtal/topology/IndexedDigitalSurface.h"
#include "DGtal/shapes/Shapes.h"
using namespace std;
using namespace DGtal;
using namespace DGtal::Z3i;
int main( int argc, char** argv )
{
const double R = argc >= 2 ? atof( argv[ 1 ] ) : 13.0; // radius of ball
const unsigned int KN = argc >= 3 ? atoi( argv[ 2 ] ) : 6; // size of neighborhood
const int M = (int) ceil( R + 2.0 );
trace.beginBlock( "Creating surface" );
trace.info() << "Sphere of radius " << R
<< ", " << KN << "-ring neighborhood" << std::endl;
typedef DigitalSetBoundary < KSpace, DigitalSet > DigitalSurfaceContainer;
Point p1( -M, -M, -M );
Point p2( M, M, M );
K.init( p1, p2, true );
DigitalSet aSet( Domain( p1, p2 ) );
Shapes<Domain>::addNorm2Ball( aSet, Point( 0, 0, 0 ), R );
auto dsc = CountedPtr< DigitalSurfaceContainer >( new DigitalSurfaceContainer( K, aSet ) );
DigSurface dsurf( dsc );
trace.info() << "[OK]"
<< " Sphere has " << dsurf.size() << " vertices/surfels"
<< " with euler " << dsurf.Euler() << std::endl;
trace.beginBlock( "Estimating normals" );
auto v2n = dsurf.makeVertexMap( RealPoint() );
for ( auto v : dsurf )
{
int nbv = 0;
BFSVisitor bfv( dsurf, v );
while( ! bfv.finished() )
{ // Vertex are colored according to the distance to initial seed.
auto node = bfv.current();
if ( KN < node.second ) break;
auto surfel = dsurf.surfel( node.first );
Dimension k = K.sOrthDir( surfel );
nv[ k ] = K.sDirect( surfel, k ) ? 1.0 : -1.0;
normal += nv;
nbv += 1;
bfv.expand();
}
normal /= nbv;
v2n[ v ] = normal.getNormalized();
}
trace.beginBlock( "Estimating area" );
const double area_true = 4.0 * M_PI * R * R;
double area_averaged = 0.0;
double area_corrected = 0.0;
for ( auto v : dsurf )
{
auto surfel = dsurf.surfel( v );
Dimension k = K.sOrthDir( surfel );
area_corrected += fabs( v2n[ v ][ k ] );
area_averaged += 1.0 / v2n[ v ].norm( RealPoint::L_1 );
}
trace.info() << "- true area = " << area_true << std::endl;
trace.info() << "- corrected area = " << area_corrected << std::endl;
trace.info() << "- averaged area = " << area_averaged << std::endl;
return 0;
}
DGtal::PointVector::zero
static Self zero
Static const for zero PointVector.
Definition: PointVector.h:1595
DGtal::Trace::endBlock
double endBlock()
DGtal::IndexedDigitalSurface
Aim: Represents a digital surface with the topology of its dual surface. Its aim is to mimick the sta...
Definition: IndexedDigitalSurface.h:96
DGtal::KhalimskySpaceND::init
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
DGtal::trace
Trace trace
Definition: Common.h:154
DGtal::DigitalSetBoundary
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as the boundary of a given...
Definition: DigitalSetBoundary.h:69
Aim: This class is useful to perform a breadth-first exploration of a graph given a starting point or...
K
KSpace K
Definition: testCubicalComplex.cpp:62
DGtal::Dimension
DGtal::uint32_t Dimension
Definition: Common.h:137
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::PointVector::L_1
@ L_1
Definition: PointVector.h:1493
DGtal::KhalimskySpaceND::sDirect
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::Trace::info
std::ostream & info()
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
Domain
HyperRectDomain< Space > Domain
Definition: testSimpleRandomAccessRangeFromPoint.cpp:44
DGtal::PointVector::getNormalized
PointVector< dim, double, std::array< double, dim > > getNormalized() const
DGtal::Z3i
Z3i this namespace gathers the standard of types for 3D imagery.
DGtal::CountedPtr< DigitalSurfaceContainer >
main
int main(int argc, char **argv)
Definition: testArithmeticDSS-benchmark.cpp:147
DGtal::KhalimskySpaceND::sOrthDir
Dimension sOrthDir(const SCell &s) const
Given a signed surfel [s], returns its orthogonal direction (ie, the coordinate where the surfel is c...
DGtal::Shapes
Aim: A utility class for constructing different shapes (balls, diamonds, and others).
Definition: DGtal/shapes/Shapes.h:71
DGtal::PointVector< dim, Integer >
Point
MyPointD Point
Definition: testClone2.cpp:383
RealPoint
Z2i::RealPoint RealPoint
Definition: testAstroid2D.cpp:46
DGtal::DigitalSetByAssociativeContainer
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
Definition: DigitalSetByAssociativeContainer.h:89
DGtal::KhalimskySpaceND
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Definition: KhalimskySpaceND.h:64