Computation of mean and Gaussian curvatures on a mesh defined by a an implicit shape discretized as a digital surface, using constant or interpolated corrected curvature measures (based on the theory of corrected normal currents). It uses a digital normal vector estimator to improve curvature estimations. Errors with respect to true expected curvatures are also computed.
This first example uses constant per face corrected normal vector field to compute curvatures.
This second example uses vertex-interpolated corrected normal vector field to compute curvatures.
#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"
{
return gradcmap;
}
void usage(
int argc,
char* argv[] )
{
std::cout << "Usage: " << std::endl
<< "\t" << argv[ 0 ] << " <P> <B> <h> <R> <mode>" << std::endl
<< std::endl
<< "Computation of mean and Gaussian curvatures on an " << std::endl
<< "digitized implicit shape using constant or " << std::endl
<< "interpolated corrected curvature measures (based " << std::endl
<< "on the theory of corrected normal currents)." << 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
<< "- <R> is the radius of the measuring balls" << 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 produces several OBJ files to display mean and" << std::endl
<< "Gaussian curvature estimation results: `example-cnc-H.obj`" << std::endl
<< "and `example-cnc-G.obj` as well as the associated MTL file." << std::endl;
std::cout << "You may either write your own polynomial as 3*x^2*y-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 )
{
return 0;
}
std::string poly = argv[ 1 ];
const double B = argc > 2 ? atof( argv[ 2 ] ) : 1.0;
const double h = argc > 3 ? atof( argv[ 3 ] ) : 1.0;
const double R = argc > 4 ? atof( argv[ 4 ] ) : 2.0;
std::string mode = argc > 5 ? argv[ 5 ] : "Const";
bool interpolated = mode == "Interp";
if ( interpolated )
std::cout << "Using vertex-*Interpolated* Corrected Normal Current" << std::endl;
else
std::cout << "Using face-*Constant* Corrected Normal Current" << std::endl;
auto params = SH::defaultParameters() | SHG::defaultParameters();
params( "t-ring", 3 )( "surfaceTraversal", "Default" );
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 )
{
<< 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;
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 : *surface )
{
const auto primal_surfel_vtcs = SH::getPointelRange(
K, surfel );
for ( auto&& primal_vtx : primal_surfel_vtcs )
face.push_back( c2i[ primal_vtx ] );
faces.push_back( face );
}
faces.cbegin(), faces.cend() );
auto exp_H = SHG::getMeanCurvatures( shape,
K, surfels, params );
auto exp_G = SHG::getGaussianCurvatures( shape,
K, surfels, params );
CNC cnc( smesh );
auto face_normals = SHG::getCTrivialNormalVectors( surface, surfels, params );
smesh.setFaceNormals( face_normals.cbegin(), face_normals.cend() );
if ( interpolated ) smesh.computeVertexNormalsFromFaceNormals();
auto mu0 = cnc.computeMu0();
auto mu1 = cnc.computeMu1();
auto mu2 = cnc.computeMu2();
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, R, f );
H[ f ] = cnc.meanCurvature ( area, mu1.measure( b, R, f ) );
G[ f ] = cnc.GaussianCurvature( area, mu2.measure( b, R, f ) );
}
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() );
std::cout << "Expected mean curvatures:"
<< " min=" << *exp_H_min_max.first << " max=" << *exp_H_min_max.second
<< std::endl;
std::cout << "Computed mean curvatures:"
<< " min=" << *H_min_max.first << " max=" << *H_min_max.second
<< std::endl;
std::cout << "Expected Gaussian curvatures:"
<< " min=" << *exp_G_min_max.first << " max=" << *exp_G_min_max.second
<< std::endl;
std::cout << "Computed Gaussian curvatures:"
<< " min=" << *G_min_max.first << " max=" << *G_min_max.second
<< 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_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;
smesh.vertexNormals() = SH::RealVectors();
smesh.faceNormals() = SH::RealVectors();
const double Hmax =
std::max( fabs( *exp_H_min_max.first ),
fabs( *exp_H_min_max.second ) );
const double Gmax =
std::max( fabs( *exp_G_min_max.first ),
fabs( *exp_G_min_max.second ) );
auto colorsH = SMW::Colors( smesh.nbFaces() );
auto colorsG = SMW::Colors( smesh.nbFaces() );
for ( auto i = 0; i < smesh.nbFaces(); i++ )
{
colorsH[ i ] = colormapH( H[ i ] );
colorsG[ i ] = colormapG( G[ i ] );
}
SMW::writeOBJ( "example-cnc-H", smesh, colorsH );
SMW::writeOBJ( "example-cnc-G", smesh, colorsG );
return 0;
}