DGtal  1.3.beta
geometry/meshes/vol-curvature-measures-icnc-XY-3d.cpp

Computation of principal curvatures and directions on a mesh defined by a VOL digital file, using interpolated corrected curvature measures (based on the theory of corrected normal currents). It uses a digital normal vector estimator to improve curvature estimations.

# "Al" vol file
./examples/geometry/meshes/vol-curvature-measures-icnc-XY-3d ../examples/samples/Al.100.vol 2.0 0 1 0.33

outputs

- Domain size is 100 x 100 x 100
- digital shape has 70413 voxels.
- surface has 21239 surfels.
[SurfaceMesh (OK) #V=21278 #VN=0 #E=42522 #F=21239 #FN=0 E[IF]=4.00151 E[IV]=3.9968 E[IFE]=1.99793]
- CTrivial normal t-ring=3 (discrete)
Computed k1 curvatures: min=-0.591355 max=0.332964
Computed k2 curvatures: min=-0.0853815 max=0.559684

It also produces several OBJ files to display curvature estimation results, example-cnc-K1.obj, example-cnc-D1.obj, example-cnc-K2.obj, and example-cnc-D2.obj as well as the associated MTL file.

Interpolated corrected smallest principal curvature and direction, r=2
Interpolated corrected greatest principal curvature and direction, r=2
See also
Curvature measures on meshes and digital surfaces
Note
In opposition with Normal Cycle curvature measures, (interpolated) corrected curvature measures can take into account an external normal vector field to estimate curvatures with better accuracy.
#include <iostream>
#include <algorithm>
#include "DGtal/base/Common.h"
#include "DGtal/math/linalg/EigenDecomposition.h"
#include "DGtal/shapes/SurfaceMesh.h"
#include "DGtal/geometry/meshes/CorrectedNormalCurrentComputer.h"
#include "DGtal/io/writers/SurfaceMeshWriter.h"
#include "DGtal/io/colormaps/GradientColorMap.h"
#include "DGtal/helpers/Shortcuts.h"
#include "DGtal/helpers/ShortcutsGeometry.h"
#include "DGtal/io/readers/SurfaceMeshReader.h"
#include "DGtal/io/colormaps/GradientColorMap.h"
#include "DGtal/io/colormaps/QuantifiedColorMap.h"
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[] )
{
std::cout << "Usage: " << std::endl
<< "\t" << argv[ 0 ] << " <filename.vol> <R> <m> <M> <Kmax>" << std::endl
<< std::endl
<< "Computation of principal curvatures and directions on a vol file, " << std::endl
<< "using interpolated corrected curvature measures (based " << std::endl
<< "on the theory of corrected normal currents)." << std::endl
<< "- builds the surface mesh from file <filename.obj>" << std::endl
<< "- <R> is the radius of the measuring balls." << std::endl
<< "- <m> is the min threshold value for the vol file" << std::endl
<< "- <M> is the max threshold value for the vol file" << std::endl
<< "- <Kmax> gives the colormap range [-Kmax,Kmax] for" << std::endl
<< " the output of principal curvatures estimates" << std::endl
<< "It produces several OBJ files to display principal " << std::endl
<< "curvatures and directions estimations: `example-cnc-K1.obj`" << std::endl
<< "`example-cnc-K2.obj`, `example-cnc-D1.obj`, and" << std::endl
<< "`example-cnc-D2.obj` as well as associated MTL files." << std::endl;
}
int main( int argc, char* argv[] )
{
if ( argc <= 1 )
{
usage( argc, argv );
return 0;
}
using namespace DGtal;
using namespace DGtal::Z3i;
typedef Shortcuts< KSpace > SH;
// VOL file
std::string input = argv[ 1 ];
const double R = argc > 2 ? atof( argv[ 2 ] ) : 2.0; // radius of measuring ball
const int m = argc > 3 ? atoi( argv[ 3 ] ) : 0; // min threshold
const int M = argc > 4 ? atoi( argv[ 4 ] ) : 1; // max threshold
const double Kmax = argc > 5 ? atof( argv[ 5 ] ) : 0.33; // range mean curvature colormap
// Read VOL file and build digital surface
auto params = SH::defaultParameters() | SHG::defaultParameters();
params( "thresholdMin", m )( "thresholdMax", M )( "closed", 1);
params( "t-ring", 3 )( "surfaceTraversal", "Default" );
auto bimage = SH::makeBinaryImage( input.c_str(), params );
if ( bimage == nullptr )
{
trace.error() << "Unable to read file <" << input.c_str() << ">" << std::endl;
return 1;
}
auto K = SH::getKSpace( bimage, params );
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 embedder( c ); } );
for ( auto&& surfel : *surface )
{
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;
// Builds a CorrectedNormalCurrentComputer object onto the SurfaceMesh object
CNC cnc( smesh );
// Estimates normal vectors using Convolved Trivial Normal estimator
auto face_normals = SHG::getCTrivialNormalVectors( surface, surfels, params );
smesh.setFaceNormals( face_normals.cbegin(), face_normals.cend() );
if ( smesh.vertexNormals().empty() )
smesh.computeVertexNormalsFromFaceNormals();
// computes area, anisotropic XY curvature measures
auto mu0 = cnc.computeMu0();
auto muXY = cnc.computeMuXY();
// estimates principal curvatures (K1,K2) and directions (D1,D2) by
// measure normalization.
std::vector< double > K1( smesh.nbFaces() );
std::vector< double > K2( smesh.nbFaces() );
std::vector< RealVector > D1( smesh.nbFaces() );
std::vector< RealVector > D2( smesh.nbFaces() );
for ( auto f = 0; f < smesh.nbFaces(); ++f )
{
const auto b = smesh.faceCentroid( f );
const auto N = smesh.faceNormals()[ f ];
const auto area = mu0 .measure( b, R, f );
const auto M = muXY.measure( b, R, f );
std::tie( K1[ f ], K2[ f ], D1[ f ], D2[ f ] )
= cnc.principalCurvatures( area, M, N );
}
auto K1_min_max = std::minmax_element( K1.cbegin(), K1.cend() );
auto K2_min_max = std::minmax_element( K2.cbegin(), K2.cend() );
std::cout << "Computed k1 curvatures:"
<< " min=" << *K1_min_max.first << " max=" << *K1_min_max.second
<< std::endl;
std::cout << "Computed k2 curvatures:"
<< " min=" << *K2_min_max.first << " max=" << *K2_min_max.second
<< std::endl;
// Remove normals for better blocky display.
smesh.vertexNormals() = SH::RealVectors();
smesh.faceNormals() = SH::RealVectors();
const auto colormapK1 = makeQuantifiedColorMap( makeColorMap( -Kmax, Kmax ) );
const auto colormapK2 = makeQuantifiedColorMap( makeColorMap( -Kmax, Kmax ) );
auto colorsK1 = SMW::Colors( smesh.nbFaces() );
auto colorsK2 = SMW::Colors( smesh.nbFaces() );
for ( auto i = 0; i < smesh.nbFaces(); i++ )
{
colorsK1[ i ] = colormapK1( K1[ i ] );
colorsK2[ i ] = colormapK2( K2[ i ] );
}
SMW::writeOBJ( "example-cnc-K1", smesh, colorsK1 );
SMW::writeOBJ( "example-cnc-K2", smesh, colorsK2 );
const auto avg_e = smesh.averageEdgeLength();
SH::RealPoints positions( smesh.nbFaces() );
for ( auto f = 0; f < positions.size(); ++f )
{
D1[ f ] *= smesh.localWindow( f );
positions[ f ] = smesh.faceCentroid( f ) - 0.5 * D1[ f ];
}
SH::saveVectorFieldOBJ( positions, D1, 0.05 * avg_e, SH::Colors(),
"example-cnc-D1",
SH::Color::Black, SH::Color( 0, 128, 0 ) );
for ( auto f = 0; f < positions.size(); ++f )
{
D2[ f ] *= smesh.localWindow( f );
positions[ f ] = smesh.faceCentroid( f ) - 0.5 * D2[ f ];
}
SH::saveVectorFieldOBJ( positions, D2, 0.05 * avg_e, SH::Colors(),
"example-cnc-D2",
SH::Color::Black, SH::Color(128, 0,128 ) );
return 0;
}
DGtal::Trace::error
std::ostream & error()
DGtal::Color
Structure representing an RGB triple with alpha component.
Definition: Color.h:67
DGtal::trace
Trace trace
Definition: Common.h:154
K
KSpace K
Definition: testCubicalComplex.cpp:62
DGtal::SurfaceMesh
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:91
DGtal::SurfaceMeshWriter
Aim: An helper class for writing mesh file formats (Waverfront OBJ at this point) and creating a Surf...
Definition: SurfaceMeshWriter.h:64
makeColorMap
DGtal::GradientColorMap< double > makeColorMap(double min_value, double max_value)
[curvature-comparator-Includes]
Definition: curvature-comparator-ii-cnc-3d.cpp:89
DGtal::Shortcuts
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:104
DGtal::ProbingMode::R
@ R
DGtal::CorrectedNormalCurrentComputer
Aim: Utility class to compute curvature measures induced by (1) a corrected normal current defined by...
Definition: CorrectedNormalCurrentComputer.h:69
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)
DGtal::Trace::info
std::ostream & info()
Vertices
SMesh::Vertices Vertices
Definition: fullConvexitySphereGeodesics.cpp:118
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::Z3i
Z3i this namespace gathers the standard of types for 3D imagery.
DGtal::ShortcutsGeometry
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: ShortcutsGeometry.h:76
DGtal::SurfaceMeshReader
Aim: An helper class for reading mesh files (Wavefront OBJ at this point) and creating a SurfaceMesh.
Definition: SurfaceMeshReader.h:63
main
int main(int argc, char **argv)
Definition: testArithmeticDSS-benchmark.cpp:147
DGtal::GradientColorMap::addColor
void addColor(const Color &color)
usage
void usage(int, char **argv)
Definition: approximation.cpp:64
Cell
KSpace::Cell Cell
Definition: testCubicalComplex.cpp:56
DGtal::Z2i::K2
KSpace K2
Definition: StdDefs.h:78
DGtal::makeQuantifiedColorMap
QuantifiedColorMap< TColorMap > makeQuantifiedColorMap(TColorMap colormap, int nb=50)
Definition: QuantifiedColorMap.h:113
DGtal::GradientColorMap
Aim: This class template may be used to (linearly) convert scalar values in a given range into a colo...
Definition: GradientColorMap.h:119