DGtalTools  1.2.0
vol2normalField.cpp
1 
31 #include <iostream>
32 #include <iterator>
33 #include "DGtal/base/Common.h"
34 #include "DGtal/topology/CanonicDigitalSurfaceEmbedder.h"
35 #include "DGtal/topology/DigitalSurface.h"
36 #include "DGtal/topology/DigitalSetBoundary.h"
37 #include "DGtal/topology/ImplicitDigitalSurface.h"
38 #include "DGtal/topology/LightImplicitDigitalSurface.h"
39 #include "DGtal/topology/ExplicitDigitalSurface.h"
40 #include "DGtal/topology/LightExplicitDigitalSurface.h"
41 #include "DGtal/graph/BreadthFirstVisitor.h"
42 #include "DGtal/topology/helpers/FrontierPredicate.h"
43 #include "DGtal/topology/helpers/BoundaryPredicate.h"
44 #include "DGtal/graph/CUndirectedSimpleLocalGraph.h"
45 #include "DGtal/graph/CUndirectedSimpleGraph.h"
46 
47 #include "DGtal/io/readers/VolReader.h"
48 #include "DGtal/images/imagesSetsUtils/SetFromImage.h"
49 #include "DGtal/images/SimpleThresholdForegroundPredicate.h"
50 #include "DGtal/images/ImageSelector.h"
51 #include "DGtal/shapes/Shapes.h"
52 #include "DGtal/helpers/StdDefs.h"
53 #include "DGtal/kernel/CanonicEmbedder.h"
54 
55 #include "DGtal/geometry/surfaces/estimation/CNormalVectorEstimator.h"
56 #include "DGtal/geometry/surfaces/estimation/BasicConvolutionWeights.h"
57 #include "DGtal/geometry/surfaces/estimation/LocalConvolutionNormalVectorEstimator.h"
58 #include "DGtal/geometry/surfaces/estimation/DigitalSurfaceEmbedderWithNormalVectorEstimator.h"
59 
60 #include "CLI11.hpp"
61 
63 
64 using namespace std;
65 using namespace DGtal;
66 using namespace Z3i;
67 
125 void missingParam ( std::string param )
126 {
127  trace.error() <<" Parameter: "<<param<<" is required..";
128  trace.info() <<std::endl;
129  exit ( 1 );
130 }
131 
132 
133 int main ( int argc, char**argv )
134 {
135 
136  // parse command line CLI ----------------------------------------------
137  CLI::App app;
138  std::string filename;
139  std::string outputFileName;
140  unsigned int level {0};
141  double sigma {5.0};
142  unsigned int neighborhood {10};
143  double normExport {1.0};
144 
145  app.description("Generates normal vector field from a vol file using DGtal library.\n Typical use example:\n \t vol2normalField[options] --input <volFileName> --o <outputFileName>\n");
146  app.add_option("-i,--input,1",filename,"Input vol file.")->required()->check(CLI::ExistingFile);
147  app.add_option("-o,--output,2",outputFileName,"Output file.")->required();
148  app.add_option("--level,-l",level,"Iso-level for the surface construction (default 0).",true);
149  app.add_option("--sigma,-s", sigma,"Sigma parameter of the Gaussian kernel (default 5.0).",true);
150  auto expOpt = app.add_flag("--exportOriginAndExtremity", "exports the origin and extremity of the vector fields when exporting the vector field in TXT format (useful to be displayed in other viewer like meshViewer).");
151  app.add_option("--vectorsNorm,-N", normExport, "set the norm of the exported vectors in TXT format (when the extremity points are exported with --exportOriginAndExtremity). By using a negative value you will invert the direction of the vectors (default 1.0).",true);
152  app.add_option("--neighborhood,-n", neighborhood,"Size of the neighborhood for the convolution (distance on surfel graph, default 10).",true);
153 
154  app.get_formatter()->column_width(40);
155  CLI11_PARSE(app, argc, argv);
156  // END parse command line using CLI ----------------------------------------------
157 
158  typedef ImageSelector < Z3i::Domain, unsigned char>::Type Image;
159  Image image = VolReader<Image>::importVol ( filename );
160 
161  trace.info() <<image<<std::endl;
162 
163  functors::SimpleThresholdForegroundPredicate<Image> simplePredicate ( image, level );
164 
165  KSpace ks;
166  bool space_ok = ks.init ( image.domain().lowerBound(),
167  image.domain().upperBound(), true );
168  if ( !space_ok )
169  {
170  trace.error() << "Error in the Khamisky space construction."<<std::endl;
171  return 2;
172  }
173 
174  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
175  MySurfelAdjacency surfAdj ( true ); // interior in all directions.
176 
177  //Set up digital surface.
178  typedef LightImplicitDigitalSurface<KSpace, functors::SimpleThresholdForegroundPredicate<Image> > MyDigitalSurfaceContainer;
179  typedef DigitalSurface<MyDigitalSurfaceContainer> MyDigitalSurface;
180  SCell bel = Surfaces<KSpace>::findABel ( ks, simplePredicate );
181 
182  MyDigitalSurfaceContainer* ptrSurfContainer =
183  new MyDigitalSurfaceContainer ( ks, simplePredicate, surfAdj, bel );
184  MyDigitalSurface digSurf ( ptrSurfContainer ); // acquired
185  MyDigitalSurface::ConstIterator it = digSurf.begin();
186 
187  // Embedder definition
188  typedef CanonicDigitalSurfaceEmbedder<MyDigitalSurface> SurfaceEmbedder;
189  SurfaceEmbedder surfaceEmbedder ( digSurf );
190 
191  //Convolution kernel
192  deprecated::GaussianConvolutionWeights < MyDigitalSurface::Size > Gkernel ( sigma );
193 
194  //Estimator definition
195  typedef deprecated::LocalConvolutionNormalVectorEstimator < MyDigitalSurface,
196  deprecated::GaussianConvolutionWeights< MyDigitalSurface::Size> > MyGaussianEstimator;
197  BOOST_CONCEPT_ASSERT ( ( concepts::CNormalVectorEstimator< MyGaussianEstimator > ) );
198  MyGaussianEstimator myNormalEstimatorG ( digSurf, Gkernel );
199 
200  // Embedder definition
201  typedef DigitalSurfaceEmbedderWithNormalVectorEstimator<SurfaceEmbedder,MyGaussianEstimator> SurfaceEmbedderWithGaussianNormal;
202  SurfaceEmbedderWithGaussianNormal mySurfelEmbedderG ( surfaceEmbedder, myNormalEstimatorG );
203 
204  // Compute normal vector field and displays it.
205  myNormalEstimatorG.init ( 1.0, neighborhood );
206 
207  trace.info() << "Generating the NOFF surface "<< std::endl;
208  ofstream out2 ( ( outputFileName + ".off" ).c_str() );
209  if ( out2.good() )
210  digSurf.exportAs3DNOFF ( out2 ,mySurfelEmbedderG );
211  out2.close();
212 
213  trace.info() << "Generating the polar coordinates file"<< std::endl;
214  ofstream out3 ( ( outputFileName + ".txt" ).c_str() );
215  if ( out3.good() )
216  {
217  MyGaussianEstimator::Quantity res;
218  for ( MyDigitalSurface::ConstIterator it =digSurf.begin(),
219  itend = digSurf.end(); it != itend; ++it )
220  {
221  res = myNormalEstimatorG.eval ( it );
222  //We output Theta - Phi
223  out3<< acos ( res [2] ) *180.0/M_PI <<" " << ( atan2 ( res [1], res [0] ) + M_PI ) *180.0/M_PI;
224 
225  if (expOpt->count()>0)
226  {
227  res *= normExport;
228  out3 << " " << mySurfelEmbedderG(*it)[0]
229  << " " << mySurfelEmbedderG(*it)[1]
230  << " " << mySurfelEmbedderG(*it)[2] << " "
231  << mySurfelEmbedderG(*it)[0]+res[0] << " "
232  << mySurfelEmbedderG(*it)[1]+res[1] << " "
233  << mySurfelEmbedderG(*it)[2]+res[2];
234 
235  }
236 
237  out3 <<std::endl;
238  }
239  }
240  out3.close();
241 
242  return 0;
243 }
244 
245 
246 // //
248 
Definition: ATu0v1.h:57