DGtal  1.5.beta
geometry/meshes/curvature-comparator-ii-cnc-3d.cpp

Computation of mean and Gaussian curvatures on a mesh defined by a an implicit shape discretized as a digital surface, using Integral Invariant (II) method or constant (CCNC) or interpolated corrected curvature measures (ICNC) (based on the theory of corrected normal currents). For CCNC and ICNC, we use also II normal vector estimations. Timings and errors with respect to true expected curvatures are computed.

This first example compares II with CCNC for estimating curvatures on "goursat" polynomial, for gridstep 0.25

./examples/geometry/meshes/curvature-comparator-ii-cnc-3d goursat 10 0.25 Const

outputs

This second example compares II with ICNC for estimating curvatures on "goursat" polynomial, for gridstep 0.25


./examples/geometry/meshes/curvature-comparator-ii-cnc-3d goursat 10 0.25 Interp

outputs


See also
Curvature measures on meshes and digital surfaces


Note
In opposition with Normal Cycle curvature measures, constant or interpolated corrected curvature measures can take into account an external normal vector field to estimate curvatures with better accuracy.


 


#include <iostream>


#include <fstream>


#include <algorithm>


#include "DGtal/base/Common.h"


#include "DGtal/shapes/SurfaceMesh.h"


#include "DGtal/geometry/meshes/CorrectedNormalCurrentComputer.h"


 


#include "DGtal/helpers/Shortcuts.h"


#include "DGtal/helpers/ShortcutsGeometry.h"


#include "DGtal/io/writers/SurfaceMeshWriter.h"


#include "DGtal/io/colormaps/GradientColorMap.h"


#include "DGtal/io/colormaps/QuantifiedColorMap.h"


 


DGtal::GradientColorMap< double >


makeColorMap( double min_value, double max_value )


{


  DGtal::GradientColorMap< double > gradcmap( min_value, max_value );


  gradcmap.addColor( DGtal::Color( 0, 0, 255 ) );


  gradcmap.addColor( DGtal::Color( 0, 255, 255 ) );


  gradcmap.addColor( DGtal::Color( 255, 255, 255 ) );


  gradcmap.addColor( DGtal::Color( 255, 255, 0 ) );


  gradcmap.addColor( DGtal::Color( 255, 0, 0 ) );


  return gradcmap;


}


 


void usage( int argc, char* argv[] )


{


  using namespace DGtal;


  using namespace DGtal::Z3i;


  typedef Shortcuts< KSpace >          SH;


  std::cout << "Usage: " << std::endl


            << "\t" << argv[ 0 ] << " <P> <B> <h> <mode>" << std::endl


            << std::endl


            << "Compare Integral Invariant (II) curvature estimations "  << std::endl


            << "with Corrected Normal Current (CNC) estimations, either "<< std::endl


            << "interpolated (Interp) or constant per face (Const)."     << std::endl


            << "- builds the surface mesh from polynomial <P>"           << std::endl


            << "- <B> defines the digitization space size [-B,B]^3"      << std::endl


            << "- <h> is the gridstep digitization"                      << std::endl


            << "- <mode> is either Const for constant corrected normal"  << std::endl


            << "  vector field or Interp for interpolated corrected"     << std::endl


            << "  normal vector field."                                  << std::endl


            << "It outputs timings and accuracy statistics for both "    << std::endl


            << "methods."                                                << std::endl;


  std::cout << "You may either write your own polynomial as 3*x^2-z^2*x*y+1" << std::endl


            <<"or use a predefined polynomial in the following list:" << std::endl;


  auto L = SH::getPolynomialList();


  for ( const auto& p : L )


    std::cout << p.first << " : " << p.second << std::endl;


}


 


int main( int argc, char* argv[] )


{


  if ( argc <= 1 )


    {


      usage( argc, argv );


      return 0;


    }


  using namespace DGtal;


  using namespace DGtal::Z3i;


  typedef SurfaceMesh< RealPoint, RealVector >                    SM;


  typedef CorrectedNormalCurrentComputer< RealPoint, RealVector > CNC;


  typedef Shortcuts< KSpace >          SH;


  typedef ShortcutsGeometry< KSpace >  SHG;


  std::string  poly = argv[ 1 ]; // polynomial


  const double    B = argc > 2 ? atof( argv[ 2 ] ) : 1.0; // max ||_oo bbox


  const double    h = argc > 3 ? atof( argv[ 3 ] ) : 1.0; // gridstep  


  std::string  mode = argc > 4 ? argv[ 4 ] : "Const"; // either Const or Interp


  bool interpolated = mode == "Interp";


  if ( interpolated )


    trace.info() << "Using vertex-*Interpolated* Corrected Normal Current" << std::endl;


  else


    trace.info() << "Using face-*Constant* Corrected Normal Current" << std::endl;


  // Read polynomial and build digital surface


  auto params = SH::defaultParameters() | SHG::defaultParameters();


  // Choose depth-first traversal to speed-up II computations.


  params( "surfaceTraversal", "DepthFirst" );


  params( "t-ring", 3 );


  params( "polynomial", poly )( "gridstep", h ); 


  params( "minAABB", -B )( "maxAABB", B );


  params( "offset", 3.0 );


  auto shape       = SH::makeImplicitShape3D( params );


  auto K           = SH::getKSpace( params );


  auto dshape      = SH::makeDigitizedImplicitShape3D( shape, params );


  auto bimage      = SH::makeBinaryImage( dshape, params );


  if ( bimage == nullptr ) 


    {


      trace.error() <<  "Unable to read polynomial <"


                    << poly.c_str() << ">" << std::endl;


      return 1;


    }


  auto sembedder   = SH::getSCellEmbedder( K );


  auto embedder    = SH::getCellEmbedder( K );


  auto surface     = SH::makeDigitalSurface( bimage, K, params );


  auto surfels     = SH::getSurfelRange( surface, params );


  trace.info() << "- surface has " << surfels.size()<< " surfels." << std::endl;


 


  params( "r-radius", 3.0 );


  params( "alpha",    0.33 );


  double ii_r = 3.0 * pow( h, 0.33 );


  Clock c;


  trace.beginBlock( "Computing II curvatures" );


  std::vector< double > HII  = SHG::getIIMeanCurvatures    ( bimage, surfels, params );


  std::vector< double > GII  = SHG::getIIGaussianCurvatures( bimage, surfels, params );


  auto ii_t = trace.endBlock();


  


  SM smesh;


  std::vector< SM::Vertices > faces;


  SH::Cell2Index c2i;


  auto pointels = SH::getPointelRange( c2i, surface );


  auto vertices = SH::RealPoints( pointels.size() );


  std::transform( pointels.cbegin(), pointels.cend(), vertices.begin(),


                  [&] (const SH::Cell& c) { return h * embedder( c ); } ); 


  for ( auto&& surfel : surfels )


    {


      const auto primal_surfel_vtcs = SH::getPointelRange( K, surfel );


      SM::Vertices face;              


      for ( auto&& primal_vtx : primal_surfel_vtcs )


        face.push_back( c2i[ primal_vtx ] );


      faces.push_back( face );


    }


  smesh.init( vertices.cbegin(), vertices.cend(),


              faces.cbegin(),    faces.cend() );


  trace.info() << smesh << std::endl;


 


  trace.beginBlock( "Computing True curvatures" );


  auto exp_H = SHG::getMeanCurvatures( shape, K, surfels, params );


  auto exp_G = SHG::getGaussianCurvatures( shape, K, surfels, params );


  trace.endBlock();


 


  trace.beginBlock( "Computing CNC curvatures" );


  // Builds a CorrectedNormalCurrentComputer object onto the SurfaceMesh object


  CNC cnc( smesh );


  // Estimates normal vectors using Integral Invariant Normal estimator 


  trace.beginBlock( "Computing II normal vectors" );


  auto face_normals = SHG::getIINormalVectors( bimage, surfels, params );


  double cnc_tn = trace.endBlock();


  // Set corrected face normals => Corrected Normal Current with


  // constant per face corrected vector field.


  smesh.setFaceNormals( face_normals.cbegin(), face_normals.cend() ); // CCNC


  // Set corrected vertex normals => Corrected Normal Current with


  // smooth linearly interpolated per face corrected vector field.


  if ( interpolated ) smesh.computeVertexNormalsFromFaceNormals();    // ICNC


  // computes area, mean and Gaussian curvature measures


  auto mu0 = cnc.computeMu0();


  auto mu1 = cnc.computeMu1();


  auto mu2 = cnc.computeMu2();


 


  // estimates mean (H) and Gaussian (G) curvatures by measure normalization.


  double cnc_mr = 1.0 * sqrt( h );


  trace.info() << "CNC measuring radius = " << cnc_mr << std::endl;


  std::vector< double > H( smesh.nbFaces() );


  std::vector< double > G( smesh.nbFaces() );


  for ( auto f = 0; f < smesh.nbFaces(); ++f )


    {


      const auto b    = smesh.faceCentroid( f );


      const auto area = mu0.measure( b, cnc_mr, f );


      H[ f ] = cnc.meanCurvature    ( area, mu1.measure( b, cnc_mr, f ) );


      G[ f ] = cnc.GaussianCurvature( area, mu2.measure( b, cnc_mr, f ) );


    }


  double cnc_t = trace.endBlock();


  


  auto HII_min_max = std::minmax_element( HII.cbegin(), HII.cend() );


  auto GII_min_max = std::minmax_element( GII.cbegin(), GII.cend() );


  auto H_min_max = std::minmax_element( H.cbegin(), H.cend() );


  auto G_min_max = std::minmax_element( G.cbegin(), G.cend() );


  auto exp_H_min_max = std::minmax_element( exp_H.cbegin(), exp_H.cend() );


  auto exp_G_min_max = std::minmax_element( exp_G.cbegin(), exp_G.cend() );


  trace.info() << "Expected mean curvatures:"


            << " min=" << *exp_H_min_max.first << " max=" << *exp_H_min_max.second


            << std::endl;


  trace.info() << "Computed II mean curvatures:"


            << " min=" << *HII_min_max.first << " max=" << *HII_min_max.second


            << std::endl;


  trace.info() << "Computed CNC mean curvatures:"


            << " min=" << *H_min_max.first << " max=" << *H_min_max.second


            << std::endl;


  trace.info() << "Expected Gaussian curvatures:"


            << " min=" << *exp_G_min_max.first << " max=" << *exp_G_min_max.second


            << std::endl;


  trace.info() << "Computed II Gaussian curvatures:"


            << " min=" << *GII_min_max.first << " max=" << *GII_min_max.second


            << std::endl;


  trace.info() << "Computed CNC Gaussian curvatures:"


            << " min=" << *G_min_max.first << " max=" << *G_min_max.second


            << std::endl;


  const auto      error_HII = SHG::getScalarsAbsoluteDifference( HII, exp_H );


  const auto stat_error_HII = SHG::getStatistic( error_HII );


  const auto   error_HII_l2 = SHG::getScalarsNormL2( HII, exp_H );


  trace.info() << "|H-H_II|_oo = " << stat_error_HII.max() << std::endl;


  trace.info() << "|H-H_II|_2  = " << error_HII_l2 << std::endl;


  const auto      error_H = SHG::getScalarsAbsoluteDifference( H, exp_H );


  const auto stat_error_H = SHG::getStatistic( error_H );


  const auto   error_H_l2 = SHG::getScalarsNormL2( H, exp_H );


  trace.info() << "|H-H_CNC|_oo = " << stat_error_H.max() << std::endl;


  trace.info() << "|H-H_CNC|_2  = " << error_H_l2 << std::endl;


  const auto      error_GII = SHG::getScalarsAbsoluteDifference( GII, exp_G );


  const auto stat_error_GII = SHG::getStatistic( error_GII );


  const auto   error_GII_l2 = SHG::getScalarsNormL2( GII, exp_G );


  trace.info() << "|G-G_II|_oo = " << stat_error_GII.max() << std::endl;


  trace.info() << "|G-G_II|_2  = " << error_GII_l2 << std::endl;


  const auto      error_G = SHG::getScalarsAbsoluteDifference( G, exp_G );


  const auto stat_error_G = SHG::getStatistic( error_G );


  const auto   error_G_l2 = SHG::getScalarsNormL2( G, exp_G );


  trace.info() << "|G-G_CNC|_oo = " << stat_error_G.max() << std::endl;


  trace.info() << "|G-G_CNC|_2  = " << error_G_l2 << std::endl;


 


  std::cout << "# " << argv[ 0 ] << std::endl


            << "# polynomial: " << poly << std::endl


            << "# CNC mode:   " << mode << std::endl;


  //              1 2          3    4    5      6     7      8


  std::cout << "# h nb_surfels ii_t ii_r ii_Hoo ii_H2 ii_Goo ii_G2 ";


  //            9      10    11     12      13     14      15


  std::cout << "cnc_tn cnc_t cnc_mr cnc_Hoo cnc_H2 cnc_Goo cnc_G2" << std::endl;


  std::cout << h << " " << surfels.size() << " " << ii_t << " " << ii_r


            << " " << stat_error_HII.max() << " " << error_HII_l2


            << " " << stat_error_GII.max() << " " << error_GII_l2


            << " " << cnc_tn << " " << cnc_t << " " << cnc_mr


            << " " << stat_error_H.max() << " " << error_H_l2


            << " " << stat_error_G.max() << " " << error_G_l2 << std::endl;


  return 0;


}


usage
void usage(int, char **argv)
Definition: approximation.cpp:64


DGtal::Color
Structure representing an RGB triple with alpha component.
Definition: Color.h:68


DGtal::GradientColorMap
Aim: This class template may be used to (linearly) convert scalar values in a given range into a colo...
Definition: GradientColorMap.h:120


DGtal::GradientColorMap::addColor
void addColor(const Color &color)


DGtal::Trace::error
std::ostream & error()


DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")


DGtal::Trace::info
std::ostream & info()


DGtal::Trace::endBlock
double endBlock()


makeColorMap
DGtal::GradientColorMap< double > makeColorMap(double min_value, double max_value)
[curvature-comparator-Includes]
Definition: curvature-comparator-ii-cnc-3d.cpp:89


surface
CountedPtr< SH3::DigitalSurface > surface
Definition: dgtalCalculus-geodesic.cpp:66


Vertices
SMesh::Vertices Vertices
Definition: fullConvexitySphereGeodesics.cpp:118


DGtal::Z3i
Z3i this namespace gathers the standard of types for 3D imagery.


DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.


DGtal::trace
Trace trace
Definition: Common.h:153




boost::vertices
std::pair< typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator, typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator > vertices(const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)


main
int main(int argc, char **argv)
Definition: testArithmeticDSS-benchmark.cpp:147


K
KSpace K
Definition: testCubicalComplex.cpp:62


Cell
KSpace::Cell Cell
Definition: testCubicalComplex.cpp:56