Computing the area of a sphere with a 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::DigitalSurface and DGtal::BreadthFirstVisitor onto it. This method is a bit slower than using an DGtal::IndexedDigitalSurface.
#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"
int main(
int argc,
char** argv )
{
const double R = argc >= 2 ? atof( argv[ 1 ] ) : 13.0;
const unsigned int KN = argc >= 3 ? atoi( argv[ 2 ] ) : 6;
const int M = (int) ceil( R + 2.0 );
trace.
info() <<
"Sphere of radius " << R
<< ", " << KN << "-ring neighborhood" << std::endl;
typedef DigitalSetBoundary < KSpace, DigitalSet > DigitalSurfaceContainer;
typedef DigitalSurface < DigitalSurfaceContainer > DigSurface;
typedef BreadthFirstVisitor < DigSurface > BFSVisitor;
DigSurface dsurf(
new DigitalSurfaceContainer(
K, aSet ) );
trace.
info() <<
"[OK]" <<
" Sphere has " << dsurf.size() <<
" vertices/surfels"
<< std::endl;
std::map< SCell, RealPoint > v2n;
for ( auto v : dsurf )
{
int nbv = 0;
BFSVisitor bfv( dsurf, v );
while( ! bfv.finished() )
{
auto node = bfv.current();
if ( KN < node.second ) break;
auto surfel = node.first;
nv[ k ] =
K.sDirect( surfel, k ) ? 1.0 : -1.0;
normal += nv;
nbv += 1;
bfv.expand();
}
normal /= nbv;
}
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 = v;
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;
}
PointVector< dim, double, std::array< double, dim > > getNormalized() const
Aim: A utility class for constructing different shapes (balls, diamonds, and others).
void beginBlock(const std::string &keyword="")
Z3i this namespace gathers the standard of types for 3D imagery.
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
HyperRectDomain< Space > Domain
Z2i::DigitalSet DigitalSet