File failed to load: https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.3/config/TeX-MML-AM_CHTML/MathJax.js
DGtal 2.0.0
area-estimation-with-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/DigitalSurface.h"
#include "DGtal/shapes/Shapes.h"
Include dependency graph for area-estimation-with-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-digital-surface.cpp.

This file is part of the DGtal library.

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

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 83 of file area-estimation-with-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;
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}
PointVector< dim, double, std::array< double, dim > > getNormalized() const
static void addNorm2Ball(TDigitalSet &aSet, const Point &aCenter, UnsignedInteger aRadius)
HyperRectDomain< Space > Domain
Definition StdDefs.h:172
KhalimskySpaceND< 3, Integer > KSpace
Definition StdDefs.h:146
Space::RealPoint RealPoint
Definition StdDefs.h:170
Space::Point Point
Definition StdDefs.h:168
DigitalSetSelector< Domain, BIG_DS+HIGH_BEL_DS >::Type DigitalSet
Definition StdDefs.h:173
DGtal::uint32_t Dimension
Definition Common.h:119
Trace trace
KSpace K

References DGtal::Shapes< TDomain >::addNorm2Ball(), DGtal::PointVector< dim, TEuclideanRing, TContainer >::getNormalized(), K, DGtal::PointVector< 3, double >::L_1, DGtal::trace, and DGtal::PointVector< 3, double >::zero.