83int main(
int argc,
char** argv )
85 const double R = argc >= 2 ? atof( argv[ 1 ] ) : 13.0;
86 const unsigned int KN = argc >= 3 ? atoi( argv[ 2 ] ) : 6;
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 );
97 K.init( p1, p2,
true );
100 DigSurface dsurf(
new DigitalSurfaceContainer(
K, aSet ) );
101 trace.info() <<
"[OK]" <<
" Sphere has " << dsurf.size() <<
" vertices/surfels"
105 trace.beginBlock(
"Estimating normals" );
106 std::map< SCell, RealPoint > v2n;
107 for (
auto v : dsurf )
111 BFSVisitor bfv( dsurf, v );
112 while( ! bfv.finished() )
114 auto node = bfv.current();
115 if ( KN < node.second )
break;
116 auto surfel = node.first;
119 nv[ k ] =
K.sDirect( surfel, k ) ? 1.0 : -1.0;
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 )
137 area_corrected += fabs( v2n[ v ][ k ] );
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;