134 using namespace DGtal;
139 typedef SH::DigitalSurface DigSurf;
142 std::string poly = argv[ 1 ];
143 const double B = argc > 2 ? atof( argv[ 2 ] ) : 1.0;
144 const double h = argc > 3 ? atof( argv[ 3 ] ) : 1.0;
145 std::string mode = argc > 4 ? argv[ 4 ] :
"Const";
146 bool interpolated = mode ==
"Interp";
148 trace.
info() <<
"Using vertex-*Interpolated* Corrected Normal Current" << std::endl;
150 trace.
info() <<
"Using face-*Constant* Corrected Normal Current" << std::endl;
153 auto params = SH::defaultParameters() | SHG::defaultParameters();
155 params(
"surfaceTraversal",
"DepthFirst" );
156 params(
"t-ring", 3 );
157 params(
"polynomial", poly )(
"gridstep", h );
158 params(
"minAABB", -B )(
"maxAABB", B );
159 params(
"offset", 3.0 );
160 auto shape = SH::makeImplicitShape3D( params );
161 auto K = SH::getKSpace( params );
162 auto dshape = SH::makeDigitizedImplicitShape3D( shape, params );
163 auto bimage = SH::makeBinaryImage( dshape, params );
164 if ( bimage ==
nullptr )
166 trace.
error() <<
"Unable to read polynomial <"
167 << poly.c_str() <<
">" << std::endl;
170 auto sembedder = SH::getSCellEmbedder(
K );
171 auto embedder = SH::getCellEmbedder(
K );
172 auto surface = SH::makeDigitalSurface( bimage,
K, params );
173 auto surfels = SH::getSurfelRange( surface, params );
174 trace.
info() <<
"- surface has " << surfels.size()<<
" surfels." << std::endl;
178 params(
"r-radius", 3.0 );
179 params(
"alpha", 0.33 );
180 double ii_r = 3.0 * pow( h, 0.33 );
183 std::vector< double > HII = SHG::getIIMeanCurvatures ( bimage, surfels, params );
184 std::vector< double > GII = SHG::getIIGaussianCurvatures( bimage, surfels, params );
190 std::vector< SM::Vertices > faces;
192 auto pointels = SH::getPointelRange( c2i, surface );
193 auto vertices = SH::RealPoints( pointels.size() );
194 std::transform( pointels.cbegin(), pointels.cend(),
vertices.begin(),
195 [&] (
const SH::Cell& c) { return h * embedder( c ); } );
196 for (
auto&& surfel : surfels )
198 const auto primal_surfel_vtcs = SH::getPointelRange(
K, surfel );
200 for (
auto&& primal_vtx : primal_surfel_vtcs )
201 face.push_back( c2i[ primal_vtx ] );
202 faces.push_back( face );
205 faces.cbegin(), faces.cend() );
211 auto exp_H = SHG::getMeanCurvatures( shape,
K, surfels, params );
212 auto exp_G = SHG::getGaussianCurvatures( shape,
K, surfels, params );
222 auto face_normals = SHG::getIINormalVectors( bimage, surfels, params );
226 smesh.setFaceNormals( face_normals.cbegin(), face_normals.cend() );
229 if ( interpolated ) smesh.computeVertexNormalsFromFaceNormals();
231 auto mu0 = cnc.computeMu0();
232 auto mu1 = cnc.computeMu1();
233 auto mu2 = cnc.computeMu2();
238 double cnc_mr = 1.0 * sqrt( h );
239 trace.
info() <<
"CNC measuring radius = " << cnc_mr << std::endl;
240 std::vector< double >
H( smesh.nbFaces() );
241 std::vector< double > G( smesh.nbFaces() );
242 for (
auto f = 0; f < smesh.nbFaces(); ++f )
244 const auto b = smesh.faceCentroid( f );
245 const auto area = mu0.measure( b, cnc_mr, f );
246 H[ f ] = cnc.meanCurvature ( area, mu1.measure( b, cnc_mr, f ) );
247 G[ f ] = cnc.GaussianCurvature( area, mu2.measure( b, cnc_mr, f ) );
253 auto HII_min_max = std::minmax_element( HII.cbegin(), HII.cend() );
254 auto GII_min_max = std::minmax_element( GII.cbegin(), GII.cend() );
255 auto H_min_max = std::minmax_element(
H.cbegin(),
H.cend() );
256 auto G_min_max = std::minmax_element( G.cbegin(), G.cend() );
257 auto exp_H_min_max = std::minmax_element( exp_H.cbegin(), exp_H.cend() );
258 auto exp_G_min_max = std::minmax_element( exp_G.cbegin(), exp_G.cend() );
259 trace.
info() <<
"Expected mean curvatures:"
260 <<
" min=" << *exp_H_min_max.first <<
" max=" << *exp_H_min_max.second
262 trace.
info() <<
"Computed II mean curvatures:"
263 <<
" min=" << *HII_min_max.first <<
" max=" << *HII_min_max.second
265 trace.
info() <<
"Computed CNC mean curvatures:"
266 <<
" min=" << *H_min_max.first <<
" max=" << *H_min_max.second
268 trace.
info() <<
"Expected Gaussian curvatures:"
269 <<
" min=" << *exp_G_min_max.first <<
" max=" << *exp_G_min_max.second
271 trace.
info() <<
"Computed II Gaussian curvatures:"
272 <<
" min=" << *GII_min_max.first <<
" max=" << *GII_min_max.second
274 trace.
info() <<
"Computed CNC Gaussian curvatures:"
275 <<
" min=" << *G_min_max.first <<
" max=" << *G_min_max.second
277 const auto error_HII = SHG::getScalarsAbsoluteDifference( HII, exp_H );
278 const auto stat_error_HII = SHG::getStatistic( error_HII );
279 const auto error_HII_l2 = SHG::getScalarsNormL2( HII, exp_H );
280 trace.
info() <<
"|H-H_II|_oo = " << stat_error_HII.max() << std::endl;
281 trace.
info() <<
"|H-H_II|_2 = " << error_HII_l2 << std::endl;
282 const auto error_H = SHG::getScalarsAbsoluteDifference( H, exp_H );
283 const auto stat_error_H = SHG::getStatistic( error_H );
284 const auto error_H_l2 = SHG::getScalarsNormL2( H, exp_H );
285 trace.
info() <<
"|H-H_CNC|_oo = " << stat_error_H.max() << std::endl;
286 trace.
info() <<
"|H-H_CNC|_2 = " << error_H_l2 << std::endl;
287 const auto error_GII = SHG::getScalarsAbsoluteDifference( GII, exp_G );
288 const auto stat_error_GII = SHG::getStatistic( error_GII );
289 const auto error_GII_l2 = SHG::getScalarsNormL2( GII, exp_G );
290 trace.
info() <<
"|G-G_II|_oo = " << stat_error_GII.max() << std::endl;
291 trace.
info() <<
"|G-G_II|_2 = " << error_GII_l2 << std::endl;
292 const auto error_G = SHG::getScalarsAbsoluteDifference( G, exp_G );
293 const auto stat_error_G = SHG::getStatistic( error_G );
294 const auto error_G_l2 = SHG::getScalarsNormL2( G, exp_G );
295 trace.
info() <<
"|G-G_CNC|_oo = " << stat_error_G.max() << std::endl;
296 trace.
info() <<
"|G-G_CNC|_2 = " << error_G_l2 << std::endl;
300 std::cout <<
"# " << argv[ 0 ] << std::endl
301 <<
"# polynomial: " << poly << std::endl
302 <<
"# CNC mode: " << mode << std::endl;
304 std::cout <<
"# h nb_surfels ii_t ii_r ii_Hoo ii_H2 ii_Goo ii_G2 ";
306 std::cout <<
"cnc_tn cnc_t cnc_mr cnc_Hoo cnc_H2 cnc_Goo cnc_G2" << std::endl;
307 std::cout << h <<
" " << surfels.size() <<
" " << ii_t <<
" " << ii_r
308 <<
" " << stat_error_HII.max() <<
" " << error_HII_l2
309 <<
" " << stat_error_GII.max() <<
" " << error_GII_l2
310 <<
" " << cnc_tn <<
" " << cnc_t <<
" " << cnc_mr
311 <<
" " << stat_error_H.max() <<
" " << error_H_l2
312 <<
" " << stat_error_G.max() <<
" " << error_G_l2 << std::endl;