DGtal  1.0.0
shortcuts-geometry.cpp
Go to the documentation of this file.
1 
30 #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
48  typedef Shortcuts< Z3i::KSpace > SH3;
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 ] ];},
144  SH3::RealVectors(), SH3::Colors(),
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  ( "Traversal", "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 ] ];},
278  SH3::RealVectors(), SH3::Colors(),
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 ] ];},
283  SH3::RealVectors(), SH3::Colors(),
284  proj_fname );
285  ++nb; nbok += ok ? 1 : 0;
286  ++nb; nbok += proj_ok ? 1 : 0;
287  }
288  }
289  trace.endBlock();
290 
291 
292  trace.info() << nbok << "/" << nb << " passed tests." << std::endl;
293  return 0;
294 }
295 // //
void beginBlock(const std::string &keyword="")
Trace trace
Definition: Common.h:144
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:101
double endBlock()
int main(int, char **)
DGtal is the top-level namespace which contains all DGtal functions and types.
std::ostream & info()
KSpace K
KSpace::Cell Cell
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...