68#include "DGtal/base/Common.h"
69#include "DGtal/shapes/SurfaceMesh.h"
71#include "DGtal/geometry/meshes/CorrectedNormalCurrentComputer.h"
73#include "DGtal/io/writers/SurfaceMeshWriter.h"
74#include "DGtal/io/colormaps/GradientColorMap.h"
75#include "DGtal/helpers/Shortcuts.h"
76#include "DGtal/io/readers/SurfaceMeshReader.h"
77#include "DGtal/io/colormaps/GradientColorMap.h"
78#include "DGtal/io/colormaps/QuantifiedColorMap.h"
92void usage(
int argc,
char* argv[] )
94 std::cout <<
"Usage: " << std::endl
95 <<
"\t" << argv[ 0 ] <<
" <filename.obj> <R> <Kmax>" << std::endl
97 <<
"Computation of principal curvatures and directions on a mesh, " << std::endl
98 <<
"using interpolated corrected curvature measures (based " << std::endl
99 <<
"on the theory of corrected normal currents)." << std::endl
100 <<
"- builds the surface mesh from file <filename.obj>" << std::endl
101 <<
"- <R> is the radius of the measuring balls." << std::endl
102 <<
"- <Kmax> gives the colormap range [-Kmax,Kmax] for" << std::endl
103 <<
" the output of principal curvatures estimates" << std::endl
104 <<
"It produces several OBJ files to display principal " << std::endl
105 <<
"curvatures and directions estimations: `example-cnc-K1.obj`" << std::endl
106 <<
"`example-cnc-K2.obj`, `example-cnc-D1.obj`, and" << std::endl
107 <<
"`example-cnc-D2.obj` as well as associated MTL files." << std::endl;
110int main(
int argc,
char* argv[] )
118 using namespace DGtal;
125 std::string input = argv[ 1 ];
126 const double R = argc > 2 ? atof( argv[ 2 ] ) : 0.0;
127 const double Kmax = argc > 3 ? atof( argv[ 3 ] ) : 5.0;
131 std::ifstream obj_stream( input.c_str() );
132 bool ok = SMR::readOBJ( obj_stream, smesh );
135 trace.
error() <<
"Unable to read file <" << input.c_str() <<
">" << std::endl;
140 for (
const auto& p : smesh.positions() )
141 lo = lo.
inf( p ), up = up.
sup( p );
142 const auto diameter = (up - lo).norm();
143 trace.
info() <<
"Mesh=" << smesh
144 <<
" diameter=" << diameter
145 <<
" radius=" << R << std::endl;
152 if ( smesh.vertexNormals().empty() )
154 if ( smesh.faceNormals().empty() )
155 smesh.computeFaceNormalsFromPositions();
156 smesh.computeVertexNormalsFromFaceNormals();
159 auto mu0 = cnc.computeMu0();
160 auto muXY = cnc.computeMuXY();
165 std::vector< double > K1( smesh.nbFaces() );
166 std::vector< double > K2( smesh.nbFaces() );
167 std::vector< RealVector > D1( smesh.nbFaces() );
168 std::vector< RealVector > D2( smesh.nbFaces() );
169 smesh.computeFaceNormalsFromPositions();
170 for (
auto f = 0; f < smesh.nbFaces(); ++f )
172 const auto b = smesh.faceCentroid( f );
173 const auto N = smesh.faceNormals()[ f ];
174 const auto area = mu0 .measure( b, R, f );
175 const auto M = muXY.measure( b, R, f );
176 std::tie( K1[ f ], K2[ f ], D1[ f ], D2[ f ] )
177 = cnc.principalCurvatures( area, M, N );
182 auto K1_min_max = std::minmax_element( K1.cbegin(), K1.cend() );
183 auto K2_min_max = std::minmax_element( K2.cbegin(), K2.cend() );
184 std::cout <<
"Computed k1 curvatures:"
185 <<
" min=" << *K1_min_max.first <<
" max=" << *K1_min_max.second
187 std::cout <<
"Computed k2 curvatures:"
188 <<
" min=" << *K2_min_max.first <<
" max=" << *K2_min_max.second
195 const auto colormapK1 = makeQuantifiedColorMap(
makeColorMap( -Kmax, Kmax ) );
196 const auto colormapK2 = makeQuantifiedColorMap(
makeColorMap( -Kmax, Kmax ) );
197 auto colorsK1 = SMW::Colors( smesh.nbFaces() );
198 auto colorsK2 = SMW::Colors( smesh.nbFaces() );
199 for (
auto i = 0; i < smesh.nbFaces(); i++ )
201 colorsK1[ i ] = colormapK1( K1[ i ] );
202 colorsK2[ i ] = colormapK2( K2[ i ] );
204 SMW::writeOBJ(
"example-cnc-K1", smesh, colorsK1 );
205 SMW::writeOBJ(
"example-cnc-K2", smesh, colorsK2 );
206 const auto avg_e = smesh.averageEdgeLength();
207 SH::RealPoints positions( smesh.nbFaces() );
208 for (
auto f = 0; f < positions.size(); ++f )
210 D1[ f ] *= smesh.localWindow( f );
211 positions[ f ] = smesh.faceCentroid( f ) - 0.5 * D1[ f ];
213 SH::saveVectorFieldOBJ( positions, D1, 0.05 * avg_e, SH::Colors(),
215 SH::Color::Black, SH::Color( 0, 128, 0 ) );
216 for (
auto f = 0; f < positions.size(); ++f )
218 D2[ f ] *= smesh.localWindow( f );
219 positions[ f ] = smesh.faceCentroid( f ) - 0.5 * D2[ f ];
221 SH::saveVectorFieldOBJ( positions, D2, 0.05 * avg_e, SH::Colors(),
223 SH::Color::Black, SH::Color(128, 0,128 ) );
Structure representing an RGB triple with alpha component.
Aim: This class template may be used to (linearly) convert scalar values in a given range into a colo...
void addColor(const Color &color)
auto inf(const PointVector< dim, OtherComponent, OtherStorage > &aPoint) const -> decltype(DGtal::inf(*this, aPoint))
Implements the infimum (or greatest lower bound).
auto sup(const PointVector< dim, OtherComponent, OtherStorage > &aPoint) const -> decltype(DGtal::sup(*this, aPoint))
Implements the supremum (or least upper bound).
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
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::GradientColorMap< double > makeColorMap(double min_value, double max_value)
[curvature-measures-Includes]
Aim: Utility class to compute curvature measures induced by (1) a corrected normal current defined by...
Aim: An helper class for reading mesh files (Wavefront OBJ at this point) and creating a SurfaceMesh.
Aim: An helper class for writing mesh file formats (Waverfront OBJ at this point) and creating a Surf...
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...