78 #include "DGtal/base/Common.h"
79 #include "DGtal/math/linalg/EigenDecomposition.h"
80 #include "DGtal/shapes/SurfaceMesh.h"
82 #include "DGtal/geometry/meshes/CorrectedNormalCurrentComputer.h"
84 #include "DGtal/io/writers/SurfaceMeshWriter.h"
85 #include "DGtal/io/colormaps/GradientColorMap.h"
86 #include "DGtal/helpers/Shortcuts.h"
87 #include "DGtal/helpers/ShortcutsGeometry.h"
88 #include "DGtal/io/readers/SurfaceMeshReader.h"
89 #include "DGtal/io/colormaps/GradientColorMap.h"
90 #include "DGtal/io/colormaps/QuantifiedColorMap.h"
104 void usage(
int argc,
char* argv[] )
106 std::cout <<
"Usage: " << std::endl
107 <<
"\t" << argv[ 0 ] <<
" <filename.vol> <R> <m> <M> <Kmax>" << std::endl
109 <<
"Computation of principal curvatures and directions on a vol file, " << std::endl
110 <<
"using interpolated corrected curvature measures (based " << std::endl
111 <<
"on the theory of corrected normal currents)." << std::endl
112 <<
"- builds the surface mesh from file <filename.obj>" << std::endl
113 <<
"- <R> is the radius of the measuring balls." << std::endl
114 <<
"- <m> is the min threshold value for the vol file" << std::endl
115 <<
"- <M> is the max threshold value for the vol file" << std::endl
116 <<
"- <Kmax> gives the colormap range [-Kmax,Kmax] for" << std::endl
117 <<
" the output of principal curvatures estimates" << std::endl
118 <<
"It produces several OBJ files to display principal " << std::endl
119 <<
"curvatures and directions estimations: `example-cnc-K1.obj`" << std::endl
120 <<
"`example-cnc-K2.obj`, `example-cnc-D1.obj`, and" << std::endl
121 <<
"`example-cnc-D2.obj` as well as associated MTL files." << std::endl;
124 int main(
int argc,
char* argv[] )
132 using namespace DGtal;
140 std::string input = argv[ 1 ];
141 const double R = argc > 2 ? atof( argv[ 2 ] ) : 2.0;
142 const int m = argc > 3 ? atoi( argv[ 3 ] ) : 0;
143 const int M = argc > 4 ? atoi( argv[ 4 ] ) : 1;
144 const double Kmax = argc > 5 ? atof( argv[ 5 ] ) : 0.33;
148 auto params = SH::defaultParameters() | SHG::defaultParameters();
149 params(
"thresholdMin", m )(
"thresholdMax", M )(
"closed", 1);
150 params(
"t-ring", 3 )(
"surfaceTraversal",
"Default" );
151 auto bimage = SH::makeBinaryImage( input.c_str(), params );
152 if ( bimage ==
nullptr )
154 trace.
error() <<
"Unable to read file <" << input.c_str() <<
">" << std::endl;
157 auto K = SH::getKSpace( bimage, params );
158 auto sembedder = SH::getSCellEmbedder(
K );
159 auto embedder = SH::getCellEmbedder(
K );
160 auto surface = SH::makeDigitalSurface( bimage,
K, params );
161 auto surfels = SH::getSurfelRange( surface, params );
162 trace.
info() <<
"- surface has " << surfels.size()<<
" surfels." << std::endl;
167 std::vector< SM::Vertices > faces;
169 auto pointels = SH::getPointelRange( c2i, surface );
170 auto vertices = SH::RealPoints( pointels.size() );
171 std::transform( pointels.cbegin(), pointels.cend(), vertices.begin(),
172 [&] (
const SH::Cell& c) { return embedder( c ); } );
173 for (
auto&& surfel : *surface )
175 const auto primal_surfel_vtcs = SH::getPointelRange(
K, surfel );
177 for (
auto&& primal_vtx : primal_surfel_vtcs )
178 face.push_back( c2i[ primal_vtx ] );
179 faces.push_back( face );
181 smesh.init( vertices.cbegin(), vertices.cend(),
182 faces.cbegin(), faces.cend() );
190 auto face_normals = SHG::getCTrivialNormalVectors( surface, surfels, params );
191 smesh.setFaceNormals( face_normals.cbegin(), face_normals.cend() );
192 if ( smesh.vertexNormals().empty() )
193 smesh.computeVertexNormalsFromFaceNormals();
195 auto mu0 = cnc.computeMu0();
196 auto muXY = cnc.computeMuXY();
202 std::vector< double > K1( smesh.nbFaces() );
203 std::vector< double > K2( smesh.nbFaces() );
204 std::vector< RealVector > D1( smesh.nbFaces() );
205 std::vector< RealVector > D2( smesh.nbFaces() );
206 for (
auto f = 0; f < smesh.nbFaces(); ++f )
208 const auto b = smesh.faceCentroid( f );
209 const auto N = smesh.faceNormals()[ f ];
210 const auto area = mu0 .measure( b, R, f );
211 const auto M = muXY.measure( b, R, f );
212 std::tie( K1[ f ], K2[ f ], D1[ f ], D2[ f ] )
213 = cnc.principalCurvatures( area, M, N );
218 auto K1_min_max = std::minmax_element( K1.cbegin(), K1.cend() );
219 auto K2_min_max = std::minmax_element( K2.cbegin(), K2.cend() );
220 std::cout <<
"Computed k1 curvatures:"
221 <<
" min=" << *K1_min_max.first <<
" max=" << *K1_min_max.second
223 std::cout <<
"Computed k2 curvatures:"
224 <<
" min=" << *K2_min_max.first <<
" max=" << *K2_min_max.second
230 smesh.vertexNormals() = SH::RealVectors();
231 smesh.faceNormals() = SH::RealVectors();
235 auto colorsK1 = SMW::Colors( smesh.nbFaces() );
236 auto colorsK2 = SMW::Colors( smesh.nbFaces() );
237 for (
auto i = 0; i < smesh.nbFaces(); i++ )
239 colorsK1[ i ] = colormapK1( K1[ i ] );
240 colorsK2[ i ] = colormapK2( K2[ i ] );
242 SMW::writeOBJ(
"example-cnc-K1", smesh, colorsK1 );
243 SMW::writeOBJ(
"example-cnc-K2", smesh, colorsK2 );
244 const auto avg_e = smesh.averageEdgeLength();
245 SH::RealPoints positions( smesh.nbFaces() );
246 for (
auto f = 0; f < positions.size(); ++f )
248 D1[ f ] *= smesh.localWindow( f );
249 positions[ f ] = smesh.faceCentroid( f ) - 0.5 * D1[ f ];
251 SH::saveVectorFieldOBJ( positions, D1, 0.05 * avg_e, SH::Colors(),
253 SH::Color::Black, SH::Color( 0, 128, 0 ) );
254 for (
auto f = 0; f < positions.size(); ++f )
256 D2[ f ] *= smesh.localWindow( f );
257 positions[ f ] = smesh.faceCentroid( f ) - 0.5 * D2[ f ];
259 SH::saveVectorFieldOBJ( positions, D2, 0.05 * avg_e, SH::Colors(),
261 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)
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
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.
QuantifiedColorMap< TColorMap > makeQuantifiedColorMap(TColorMap colormap, int nb=50)
Aim: Utility class to compute curvature measures induced by (1) a corrected normal current defined by...
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 ...
int main(int argc, char *argv[])
DGtal::GradientColorMap< double > makeColorMap(double min_value, double max_value)
[curvature-measures-Includes]
void usage(int argc, char *argv[])