DGtal  1.5.beta
shortcuts-geometry.cpp
Go to the documentation of this file.
1 
31 #include <iostream>
32 #include "ConfigExamples.h"
33 #include "DGtal/helpers/StdDefs.h"
34 #include "DGtal/helpers/Shortcuts.h"
35 #include "DGtal/helpers/ShortcutsGeometry.h"
36 #include "DGtal/base/Common.h"
38 
39 using namespace std;
40 using namespace DGtal;
41 
43 
44 int main( int /* argc */, char** /* argv */ )
45 {
46  unsigned int nb = 0, nbok = 0;
47  // 3d tests
50 
51  trace.beginBlock ( "Load vol file -> build digital surface -> estimate mean curvature -> save OBJ." );
52  {
54  auto params = SH3::defaultParameters() | SHG3::defaultParameters();
55  params( "colormap", "Tics" );
56  auto bimage = SH3::makeBinaryImage( examplesPath + "samples/Al.100.vol", params );
57  auto K = SH3::getKSpace( bimage, params );
58  auto surface = SH3::makeDigitalSurface( bimage, K, params );
59  auto surfels = SH3::getSurfelRange( surface, params );
60  auto curv = SHG3::getIIMeanCurvatures( bimage, surfels, params );
61  // To get Gauss curvatures, write instead:
62  // auto curv = SHG3::getIIGaussianCurvatures( bimage, surfels, params );
63  auto cmap = SH3::getColorMap( -0.5, 0.5, params );
64  auto colors = SH3::Colors( surfels.size() );
65  std::transform( curv.cbegin(), curv.cend(), colors.begin(), cmap );
66  bool ok = SH3::saveOBJ( surface, SH3::RealVectors(), colors,
67  "al-H-II.obj" );
69  ++nb; nbok += ok ? 1 : 0;
70  }
71  trace.endBlock();
72 
73  trace.beginBlock ( "Load vol file -> build digital surface -> estimate Gauss curvature -> save OBJ." );
74  {
75  auto params = SH3::defaultParameters() | SHG3::defaultParameters();
76  params( "colormap", "Tics" );
77  auto bimage = SH3::makeBinaryImage( examplesPath + "samples/Al.100.vol", params );
78  auto K = SH3::getKSpace( bimage, params );
79  auto surface = SH3::makeDigitalSurface( bimage, K, params );
80  auto surfels = SH3::getSurfelRange( surface, params );
81  auto curv = SHG3::getIIGaussianCurvatures( bimage, surfels, params );
82  auto cmap = SH3::getColorMap( -0.25, 0.25, params );
83  auto colors = SH3::Colors( surfels.size() );
84  std::transform( curv.cbegin(), curv.cend(), colors.begin(), cmap );
85  bool ok = SH3::saveOBJ( surface, SH3::RealVectors(), colors,
86  "al-G-II.obj" );
87  ++nb; nbok += ok ? 1 : 0;
88  }
89  trace.endBlock();
90 
91  trace.beginBlock ( "Build polynomial shape -> digitize -> extract ground-truth geometry." );
92  {
93  auto params = SH3::defaultParameters() | SHG3::defaultParameters();
95  params( "polynomial", "3*x^2+2*y^2+z^2-90" )( "gridstep", 0.25 );
96  auto implicit_shape = SH3::makeImplicitShape3D ( params );
97  auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
98  auto binary_image = SH3::makeBinaryImage ( digitized_shape, params );
99  auto K = SH3::getKSpace( params );
100  auto surface = SH3::makeLightDigitalSurface( binary_image, K, params );
101  auto surfels = SH3::getSurfelRange( surface, params );
102  auto positions = SHG3::getPositions( implicit_shape, K, surfels, params );
103  auto normals = SHG3::getNormalVectors( implicit_shape, K, surfels, params );
104  auto mean_curvs = SHG3::getMeanCurvatures( implicit_shape, K, surfels, params );
105  auto gauss_curvs = SHG3::getGaussianCurvatures( implicit_shape, K, surfels, params );
107  auto stat_mean = SHG3::getStatistic( mean_curvs );
108  auto stat_gauss = SHG3::getStatistic( gauss_curvs );
109  trace.info() << " min(H)=" << stat_mean.min()
110  << " avg(H)=" << stat_mean.mean()
111  << " max(H)=" << stat_mean.max() << std::endl;
112  trace.info() << " min(G)=" << stat_gauss.min()
113  << " avg(G)=" << stat_gauss.mean()
114  << " max(G)=" << stat_gauss.max() << std::endl;
115  ++nb; nbok += positions.size() == surfels.size() ? 1 : 0;
116  ++nb; nbok += normals.size() == surfels.size() ? 1 : 0;
117  ++nb; nbok += mean_curvs.size() == surfels.size() ? 1 : 0;
118  ++nb; nbok += gauss_curvs.size() == surfels.size() ? 1 : 0;
119  ++nb; nbok += stat_mean.min() > 0.08 ? 1 : 0;
120  ++nb; nbok += stat_gauss.min() > 0.0064 ? 1 : 0;
121  }
122  trace.endBlock();
123 
124  trace.beginBlock ( "Build polynomial shape -> digitize -> get pointels -> save projected quadrangulated surface." );
125  {
126  auto params = SH3::defaultParameters() | SHG3::defaultParameters();
128  const double h = 0.25;
129  params( "polynomial", "goursat" )( "gridstep", h );
130  auto implicit_shape = SH3::makeImplicitShape3D ( params );
131  auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
132  auto binary_image = SH3::makeBinaryImage ( digitized_shape, params );
133  auto K = SH3::getKSpace( params );
134  auto embedder = SH3::getCellEmbedder( K );
135  auto surface = SH3::makeLightDigitalSurface( binary_image, K, params );
136  SH3::Cell2Index c2i;
137  auto pointels = SH3::getPointelRange( c2i, surface );
138  SH3::RealPoints pos( pointels.size() );
139  std::transform( pointels.cbegin(), pointels.cend(), pos.begin(),
140  [&] (const SH3::Cell& c) { return h * embedder( c ); } );
141  auto ppos = SHG3::getPositions( implicit_shape, pos, params );
142  bool ok = SH3::saveOBJ( surface,
143  [&] (const SH3::Cell& c){ return ppos[ c2i[ c ] ];},
145  "goursat-quad-proj.obj" );
147  ++nb; nbok += ok ? 1 : 0;
148  }
149  trace.endBlock();
150 
151  trace.beginBlock ( "Build polynomial shape -> digitize -> extract mean curvature -> save as OBJ with colors." );
152  {
153  auto params = SH3::defaultParameters() | SHG3::defaultParameters();
155  params( "polynomial", "goursat" )( "gridstep", 0.25 )( "colormap", "Tics" );
156  auto implicit_shape = SH3::makeImplicitShape3D ( params );
157  auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
158  auto binary_image = SH3::makeBinaryImage ( digitized_shape, params );
159  auto K = SH3::getKSpace( params );
160  auto surface = SH3::makeLightDigitalSurface( binary_image, K, params );
161  auto surfels = SH3::getSurfelRange( surface, params );
162  auto mean_curv = SHG3::getMeanCurvatures( implicit_shape, K, surfels, params );
163  auto cmap = SH3::getColorMap( -0.3, 0.3, params );
164  auto colors = SH3::Colors( surfels.size() );
165  std::transform( mean_curv.cbegin(), mean_curv.cend(), colors.begin(), cmap );
166  bool ok = SH3::saveOBJ( surface, SH3::RealVectors(), colors,
167  "goursat-H.obj" );
169  ++nb; nbok += ok ? 1 : 0;
170  }
171  trace.endBlock();
172 
173  trace.beginBlock ( "Build polynomial shape -> digitize -> extract ground-truth and estimated mean curvature -> display errors in OBJ with colors." );
174  {
175  auto params = SH3::defaultParameters() | SHG3::defaultParameters();
177  params( "polynomial", "goursat" )( "gridstep", 0.25 )( "colormap", "Tics" )
178  ( "R-radius", 5.0 );
179  auto implicit_shape = SH3::makeImplicitShape3D ( params );
180  auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
181  auto bimage = SH3::makeBinaryImage ( digitized_shape, params );
182  auto K = SH3::getKSpace( params );
183  auto surface = SH3::makeLightDigitalSurface( bimage, K, params );
184  auto surfels = SH3::getSurfelRange( surface, params );
185  auto t_curv = SHG3::getMeanCurvatures( implicit_shape, K, surfels, params );
186  auto ii_curv = SHG3::getIIMeanCurvatures( bimage, surfels, params );
187  auto cmap = SH3::getColorMap( -0.5, 0.5, params );
188  auto colors = SH3::Colors( surfels.size() );
189  std::transform( t_curv.cbegin(), t_curv.cend(), colors.begin(), cmap );
190  bool ok_t = SH3::saveOBJ( surface, SH3::RealVectors(), colors, "goursat-H.obj" );
191  std::transform( ii_curv.cbegin(), ii_curv.cend(), colors.begin(), cmap );
192  bool ok_ii = SH3::saveOBJ( surface, SH3::RealVectors(), colors, "goursat-H-ii.obj" );
193  auto errors = SHG3::getScalarsAbsoluteDifference( t_curv, ii_curv );
194  auto stat_errors = SHG3::getStatistic( errors );
195  auto cmap_errors = SH3::getColorMap( 0.0, stat_errors.max(), params );
196  std::transform( errors.cbegin(), errors.cend(), colors.begin(), cmap_errors );
197  bool ok_err = SH3::saveOBJ( surface, SH3::RealVectors(), colors, "goursat-H-ii-err.obj" );
198  trace.info() << "Error Loo=" << SHG3::getScalarsNormLoo( t_curv, ii_curv )
199  << " L1=" << SHG3::getScalarsNormL1 ( t_curv, ii_curv )
200  << " L2=" << SHG3::getScalarsNormL2 ( t_curv, ii_curv )
201  << std::endl;
203  ++nb; nbok += ( ok_t && ok_ii && ok_err ) ? 1 : 0;
204  }
205  trace.endBlock();
206 
207  trace.beginBlock ( "Build polynomial shape -> digitize -> build digital surface -> save primal surface with VCM normals as obj." );
208  {
209  auto params = SH3::defaultParameters() | SHG3::defaultParameters();
211  params( "polynomial", "goursat" )( "gridstep", 0.25 )
212  ( "surfaceTraversal", "Default" );
213  auto implicit_shape = SH3::makeImplicitShape3D ( params );
214  auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
215  auto K = SH3::getKSpace( params );
216  auto binary_image = SH3::makeBinaryImage( digitized_shape, params );
217  auto surface = SH3::makeDigitalSurface( binary_image, K, params );
218  auto surfels = SH3::getSurfelRange( surface, params );
219  auto vcm_normals = SHG3::getVCMNormalVectors( surface, surfels, params );
220  bool ok = SH3::saveOBJ( surface, vcm_normals, SH3::Colors(),
221  "goursat-primal-vcm.obj" );
223  ++nb; nbok += ok ? 1 : 0;
224  }
225  trace.endBlock();
226 
227  trace.beginBlock ( "Build polynomial shape -> digitize implicitly -> estimate II normals and curvature." );
228  {
229  auto params = SH3::defaultParameters() | SHG3::defaultParameters();
231  params( "polynomial", "goursat" )( "gridstep", .25 );
232  auto implicit_shape = SH3::makeImplicitShape3D ( params );
233  auto dig_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
234  auto K = SH3::getKSpace ( params );
235  auto surface = SH3::makeDigitalSurface ( dig_shape, K, params );
236  auto surfels = SH3::getSurfelRange ( surface, params( "surfaceTraversal", "DepthFirst" ) );
237  auto def_surfels = SH3::getSurfelRange ( surface, params( "surfaceTraversal", "Default" ) );
238  auto ii_normals = SHG3::getIINormalVectors ( dig_shape, surfels, params );
239  trace.beginBlock( "II with default traversal (slower)" );
240  auto ii_mean_curv = SHG3::getIIMeanCurvatures ( dig_shape, def_surfels, params );
241  trace.endBlock();
242  trace.beginBlock( "II with depth-first traversal (faster)" );
243  auto ii_mean_curv2 = SHG3::getIIMeanCurvatures ( dig_shape, surfels, params );
244  trace.endBlock();
245  auto cmap = SH3::getColorMap ( -0.5, 0.5, params );
246  auto colors = SH3::Colors ( def_surfels.size() );
247  auto match = SH3::getRangeMatch ( def_surfels, surfels );
248  auto normals = SH3::getMatchedRange ( ii_normals, match );
249  for ( SH3::Idx i = 0; i < colors.size(); i++ )
250  colors[ i ] = cmap( ii_mean_curv[ match[ i ] ] );
251  bool ok_H = SH3::saveOBJ( surface, SH3::RealVectors(), colors, "goursat-imp-H-ii.obj" );
253  ++nb; nbok += ( ok_H && ii_mean_curv.size() == ii_mean_curv2.size() ) ? 1 : 0;
254  }
255  trace.endBlock();
256 
257  trace.beginBlock ( "Build polynomial shape -> save several projected quadrangulated surface and digitized boundaries." );
258  {
259  auto params = SH3::defaultParameters() | SHG3::defaultParameters();
260  std::vector<double> gridsteps {0.5, 0.25, 0.125};
261  for ( auto h : gridsteps ) {
262  params( "polynomial", "goursat" )( "gridstep", h );
263  auto implicit_shape = SH3::makeImplicitShape3D ( params );
264  auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
265  auto binary_image = SH3::makeBinaryImage ( digitized_shape, params );
266  auto K = SH3::getKSpace( params );
267  auto embedder = SH3::getCellEmbedder( K );
268  auto surface = SH3::makeLightDigitalSurface( binary_image, K, params );
269  SH3::Cell2Index c2i;
270  auto pointels = SH3::getPointelRange( c2i, surface );
271  SH3::RealPoints pos( pointels.size() );
272  std::transform( pointels.cbegin(), pointels.cend(), pos.begin(),
273  [&] (const SH3::Cell& c) { return h * embedder( c ); } );
274  auto ppos = SHG3::getPositions( implicit_shape, pos, params );
275  auto fname = std::string( "goursat-quad-" ) + std::to_string( h ) + std::string( ".obj" );
276  bool ok = SH3::saveOBJ( surface,
277  [&] (const SH3::Cell& c){ return pos[ c2i[ c ] ];},
279  fname );
280  auto proj_fname = std::string( "goursat-quad-proj-" ) + std::to_string( h ) + std::string( ".obj" );
281  bool proj_ok = SH3::saveOBJ( surface,
282  [&] (const SH3::Cell& c){ return ppos[ c2i[ c ] ];},
284  proj_fname );
285  ++nb; nbok += ok ? 1 : 0;
286  ++nb; nbok += proj_ok ? 1 : 0;
287  }
288  }
289  trace.endBlock();
290 
291  trace.beginBlock ( "Build polynomial shape -> digitize -> digital surface -> save primal surface and VCM normal field as obj." );
292  {
293  auto params = SH3::defaultParameters() | SHG3::defaultParameters();
295  params( "polynomial", "goursat" )( "gridstep", 0.5 )
296  ( "surfaceTraversal", "Default" );
297  auto implicit_shape = SH3::makeImplicitShape3D ( params );
298  auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
299  auto K = SH3::getKSpace( params );
300  auto binary_image = SH3::makeBinaryImage( digitized_shape, params );
301  auto surface = SH3::makeDigitalSurface( binary_image, K, params );
302  auto surfels = SH3::getSurfelRange( surface, params );
303  auto vcm_normals = SHG3::getVCMNormalVectors( surface, surfels, params );
304  auto embedder = SH3::getSCellEmbedder( K );
305  SH3::RealPoints positions( surfels.size() );
306  std::transform( surfels.cbegin(), surfels.cend(), positions.begin(),
307  [&] (const SH3::SCell& c) { return embedder( c ); } );
308  bool ok = SH3::saveOBJ( surface, vcm_normals, SH3::Colors(),
309  "goursat-primal-vcm.obj" );
310  bool ok2 = SH3::saveVectorFieldOBJ( positions, vcm_normals, 0.05, SH3::Colors(),
311  "goursat-primal-vcm-normals.obj",
312  SH3::Color( 0, 0, 0 ), SH3::Color::Red );
314  ++nb, nbok += ok ? 1 : 0;
315  ++nb, nbok += ok2 ? 1 : 0;
316  }
317  trace.endBlock();
318 
319  trace.beginBlock ( "Build polynomial shape -> digitize -> extract ground-truth curvatures -> display in OBJ." );
320  {
321  auto params = SH3::defaultParameters() | SHG3::defaultParameters();
323  params( "polynomial", "goursat" )( "gridstep", 0.25 )( "colormap", "Tics" );
324  auto implicit_shape = SH3::makeImplicitShape3D ( params );
325  auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
326  auto bimage = SH3::makeBinaryImage ( digitized_shape, params );
327  auto K = SH3::getKSpace( params );
328  auto surface = SH3::makeLightDigitalSurface( bimage, K, params );
329  auto surfels = SH3::getSurfelRange( surface, params );
330  auto k1 = SHG3::getFirstPrincipalCurvatures( implicit_shape, K, surfels, params );
331  auto k2 = SHG3::getSecondPrincipalCurvatures( implicit_shape, K, surfels, params );
332  auto d1 = SHG3::getFirstPrincipalDirections( implicit_shape, K, surfels, params );
333  auto d2 = SHG3::getSecondPrincipalDirections( implicit_shape, K, surfels, params );
334  auto embedder = SH3::getSCellEmbedder( K );
335  SH3::RealPoints positions( surfels.size() );
336  std::transform( surfels.cbegin(), surfels.cend(), positions.begin(),
337  [&] (const SH3::SCell& c) { return embedder( c ); } );
338  SH3::saveOBJ( surface, SH3::RealVectors(), SH3::Colors(),
339  "goursat-primal.obj" );
340  // output principal curvatures and directions
341  auto cmap = SH3::getColorMap( -0.5, 0.5, params );
342  auto colors= SH3::Colors( surfels.size() );
343  std::transform( k1.cbegin(), k1.cend(), colors.begin(), cmap );
344  bool ok_k1 = SH3::saveOBJ( surface, SH3::RealVectors(), colors, "goursat-primal-k1.obj" );
345  bool ok_d1 = SH3::saveVectorFieldOBJ( positions, d1, 0.05, colors,
346  "goursat-primal-d1.obj", SH3::Color::Black );
347  std::transform( k2.cbegin(), k2.cend(), colors.begin(), cmap );
348  bool ok_k2 = SH3::saveOBJ( surface, SH3::RealVectors(), colors, "goursat-primal-k2.obj" );
349  bool ok_d2 = SH3::saveVectorFieldOBJ( positions, d2, 0.05, colors,
350  "goursat-primal-d2.obj", SH3::Color::Black );
351  ASSERT(ok_k1 && ok_d1 && ok_k2 && ok_d2);
353  }
354  trace.endBlock();
355 
356 
357 #if defined(WITH_EIGEN)
358 
359  trace.beginBlock ( "Load vol file -> build main digital surface -> II normals -> AT regularization -> save OBJ with colored normals." );
360  {
361  auto params = SH3::defaultParameters() | SHG3::defaultParameters();
363  auto al_capone = SH3::makeBinaryImage( examplesPath + "samples/Al.100.vol", params );
364  auto K = SH3::getKSpace( al_capone );
365  auto surface = SH3::makeLightDigitalSurface( al_capone, K, params );
366  auto surfels = SH3::getSurfelRange( surface, params );
367  auto ii_normals = SHG3::getIINormalVectors( al_capone, surfels, params );
368  auto linels = SH3::getCellRange( surface, 1 );
369  auto uembedder = SH3::getCellEmbedder( K );
370  SH3::Scalars features( linels.size() );
371  auto at_normals = SHG3::getATVectorFieldApproximation( features, linels.cbegin(), linels.cend(),
372  surface, surfels,
373  ii_normals, params );
374  // Output normals as colors depending on directions
375  SH3::Colors colors( surfels.size() );
376  for ( size_t i = 0; i < surfels.size(); i++ )
377  colors[ i ] = SH3::Color( (unsigned char) 255.0*fabs( at_normals[ i ][ 0 ] ),
378  (unsigned char) 255.0*fabs( at_normals[ i ][ 1 ] ),
379  (unsigned char) 255.0*fabs( at_normals[ i ][ 2 ] ) );
380  bool ok1 = SH3::saveOBJ( surface, SH3::RealVectors(), SH3::Colors(), "al-surface.obj" );
381  bool ok2 = SH3::saveOBJ( surface, at_normals, colors, "al-colored-at-normals.obj" );
382  // Output discontinuities as sticks on linels.
383  SH3::RealPoints f0;
384  SH3::RealVectors f1;
385  for ( size_t i = 0; i < linels.size(); i++ )
386  {
387  if ( features[ i ] < 0.5 )
388  {
389  const SH3::Cell linel = linels[ i ];
390  const Dimension d = * K.uDirs( linel );
391  const SH3::Cell p0 = K.uIncident( linel, d, false );
392  const SH3::Cell p1 = K.uIncident( linel, d, true );
393  f0.push_back( uembedder( p0 ) );
394  f1.push_back( uembedder( p1 ) - uembedder( p0 ) );
395  }
396  }
397  bool ok3 = SH3::saveVectorFieldOBJ( f0, f1, 0.1, SH3::Colors(),
398  "al-features.obj",
399  SH3::Color( 0, 0, 0 ), SH3::Color::Red );
401  ++nb; nbok += ok1 ? 1 : 0;
402  ++nb; nbok += ok2 ? 1 : 0;
403  ++nb; nbok += ok3 ? 1 : 0;
404  }
405  trace.endBlock();
406 
407 #endif // defined(WITH_EIGEN)
408 
409  trace.info() << nbok << "/" << nb << " passed tests." << std::endl;
410  return 0;
411 }
412 // //
Structure representing an RGB triple with alpha component.
Definition: Color.h:68
static const Color Red
Definition: Color.h:416
static const Color Black
Definition: Color.h:413
Cell uIncident(const Cell &c, Dimension k, bool up) const
Return the forward or backward unsigned cell incident to [c] along axis [k], depending on [up].
DirIterator uDirs(const Cell &p) const
Given an unsigned cell [p], returns an iterator to iterate over each coordinate the cell spans.
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...
Definition: Shortcuts.h:105
std::map< Cell, IdxVertex > Cell2Index
Definition: Shortcuts.h:189
std::vector< Color > Colors
Definition: Shortcuts.h:192
std::vector< RealPoint > RealPoints
Definition: Shortcuts.h:180
std::vector< RealVector > RealVectors
Definition: Shortcuts.h:179
std::vector< Scalar > Scalars
Definition: Shortcuts.h:178
IdxVertex Idx
Definition: Shortcuts.h:181
LightDigitalSurface::SCell SCell
Definition: Shortcuts.h:163
LightDigitalSurface::Cell Cell
Definition: Shortcuts.h:162
void beginBlock(const std::string &keyword="")
std::ostream & info()
double endBlock()
CountedPtr< SH3::DigitalSurface > surface
CountedPtr< SH3::BinaryImage > binary_image
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition: Common.h:136
Trace trace
Definition: Common.h:153
int main(int, char **)
Shortcuts< KSpace > SH3
KSpace K
ShortcutsGeometry< Z3i::KSpace > SHG3