DGtalTools  1.3.beta
3dVolMarchingCubes.cpp
1 
29 #include <iostream>
32 #include <queue>
33 #include <boost/program_options/options_description.hpp>
34 #include <boost/program_options/parsers.hpp>
35 #include <boost/program_options/variables_map.hpp>
36 #include <DGtal/kernel/sets/SetPredicate.h>
37 #include <DGtal/io/readers/VolReader.h>
38 #include <DGtal/images/ImageSelector.h>
39 #include <DGtal/images/SimpleThresholdForegroundPredicate.h>
40 #include <DGtal/images/ImageLinearCellEmbedder.h>
41 #include <DGtal/shapes/Shapes.h>
42 #include <DGtal/kernel/CanonicEmbedder.h>
43 #include <DGtal/helpers/StdDefs.h>
44 #include <DGtal/topology/helpers/Surfaces.h>
45 #include <DGtal/topology/DigitalSurface.h>
46 #include <DGtal/topology/SetOfSurfels.h>
47 #include <DGtal/geometry/volumes/KanungoNoise.h>
49 
50 
51 #include "CLI11.hpp"
52 
54 
55 using namespace DGtal;
56 using namespace Z3i;
57 
59 
116 int main( int argc, char** argv )
117 {
119 
120  // parse command line using CLI ----------------------------------------------
121  CLI::App app;
122  std::string inputFileName;
123  std::string outputFileName {"marching-cubes.off"};
124  double noise;
125  double threshold {1.0};
126  unsigned int intAdjacency = 0;
127 
128  app.description("Outputs the isosurface of value <threshold> of the volume <fileName.vol> as an OFF file <output.off>. The <adjacency> (0/1) allows to choose between interior (6,18) and exterior (18,6) adjacency.");
129  app.add_option("-i,--input,1", inputFileName, "the volume file (.vol)." )
130  ->required()
131  ->check(CLI::ExistingFile);
132  app.add_option("--threshold,-t",threshold, "the value that defines the isosurface in the image (an integer between 0 and 255).", true);
133  app.add_option("--adjacency,-a",intAdjacency, "0: interior adjacency, 1: exterior adjacency", true);
134  app.add_option("-o,--output,2", outputFileName, "the output OFF file that represents the geometry of the isosurface", true );
135  auto noiseOpt = app.add_option("--noise,-n", noise, "Kanungo noise level in ]0,1[. Note that only the largest connected component is considered and that no specific embedder is used.");
136 
137  app.get_formatter()->column_width(40);
138  CLI11_PARSE(app, argc, argv);
139  // END parse command line using CLI ----------------------------------------------
140 
141 
142 
143  bool addNoise=false;
144  if (noiseOpt->count() > 0 )
145  {
146  addNoise=true;
147  }
149 
151  trace.beginBlock( "Reading vol file into an image." );
153  Image image = VolReader<Image>::importVol(inputFileName);
154 
156  ThresholdedImage thresholdedImage( image, threshold );
157  trace.endBlock();
159 
160 
162  trace.beginBlock( "Construct the Khalimsky space from the image domain." );
163  KSpace ks;
164  bool space_ok =
165  ks.init( image.domain().lowerBound(), image.domain().upperBound(), true );
166  if (!space_ok)
167  {
168  trace.error() << "Error in the Khamisky space construction."<<std::endl;
169  return 2;
170  }
171  trace.endBlock();
173 
175  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
176  MySurfelAdjacency surfAdj( intAdjacency ); // interior in all directions.
178 
180  trace.beginBlock( "Extracting boundary by scanning the space. " );
181  typedef KSpace::SurfelSet SurfelSet;
182  typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
183  typedef DigitalSurface< MySetOfSurfels > MyDigitalSurface;
184  MySetOfSurfels theSetOfSurfels( ks, surfAdj );
185 
186  if (addNoise)
187  {
188  trace.info()<<"Adding some noise."<<std::endl;
189  KanungoNoise<ThresholdedImage, Z3i::Domain> kanungo(thresholdedImage, image.domain(), noise);
190  std::vector< std::vector<SCell > > vectConnectedSCell;
191  trace.info()<<"Extracting all connected components."<<std::endl;
192  Surfaces<KSpace>::extractAllConnectedSCell( vectConnectedSCell, ks, surfAdj,
193  kanungo, false);
194  if( vectConnectedSCell.size() == 0 )
195  {
196  trace.error()<< "No surface component exists. Please check the vol file --min and --max parameters." << std::endl;
197  return 3;
198  }
199 
200  trace.info()<<vectConnectedSCell.size()<<" components."<<std::endl;
201 
202  trace.info()<<"Extracting the largest one."<<std::endl;
203  int cc_max_size_idx = -1;
204  auto it_max = std::max_element( vectConnectedSCell.begin(), vectConnectedSCell.end(),
205  [] (std::vector<SCell >& v1, std::vector<SCell >& v2)
206  { return v1.size() < v2.size(); } );
207  cc_max_size_idx = std::distance( vectConnectedSCell.begin(), it_max );
208  theSetOfSurfels.surfelSet().insert( vectConnectedSCell[ cc_max_size_idx ].begin(),
209  vectConnectedSCell[ cc_max_size_idx ].end() );
210  }
211  else
213  theSetOfSurfels.surfelSet(), ks, thresholdedImage,
214  image.domain().lowerBound(), image.domain().upperBound() );
215 
216  MyDigitalSurface digSurf( theSetOfSurfels );
217  trace.info() << "Digital surface has " << digSurf.size() << " surfels."
218  << std::endl;
219  trace.endBlock();
221 
222 
224  trace.beginBlock( "Making OFF surface. " );
225  // Describes how voxels are embedded into Euclidean space.
226  typedef CanonicEmbedder< Space > MyEmbedder;
227  // Describes how the centroid surface elements is placed in-between embedded voxels.
229  CellEmbedder cellEmbedder;
230  MyEmbedder trivialEmbedder;
231 
232  // The +0.5 is to avoid isosurface going exactly through a voxel
233  // center, especially for binary volumes.
234 
235  //Making sure that we probe inside the domain.
236  Image largerImage( Domain( image.domain().lowerBound() - KSpace::Vector::diagonal(2),
237  image.domain().upperBound() + KSpace::Vector::diagonal(2)));
238  for(auto p: image.domain())
239  largerImage.setValue( p, image(p));
240 
241  std::ofstream out( outputFileName.c_str() );
242  if (addNoise)
243  {
244  if ( out.good() )
245  digSurf.exportSurfaceAs3DOFF( out );
246  else
247  {
248  trace.error()<<"IO error while opening the output file"<<std::endl;
249  exit(2);
250  }
251  }
252  else
253  {
254  cellEmbedder.init( ks, image, trivialEmbedder,
255  ( (double) threshold ) + 0.5 );
256 
257  if ( out.good() )
258  digSurf.exportEmbeddedSurfaceAs3DOFF( out, cellEmbedder );
259  else
260  {
261  trace.error()<<"IO error while opening the output file"<<std::endl;
262  exit(2);
263  }
264  }
265 
266  out.close();
267  trace.endBlock();
269  return 0;
270 }
271 
void beginBlock(const std::string &keyword="")
static Self diagonal(Component val=1)
HyperRectDomain< Space > Domain
std::set< SCell > SurfelSet
double endBlock()
int main(int argc, char **argv)
static void sMakeBoundary(SCellSet &aBoundary, const KSpace &aKSpace, const PointPredicate &pp, const Point &aLowerBound, const Point &aUpperBound)
static ImageContainer importVol(const std::string &filename, const Functor &aFunctor=Functor())
bool init(const Point &lower, const Point &upper, bool isClosed)
static void extractAllConnectedSCell(std::vector< std::vector< SCell > > &aVectConnectedSCell, const KSpace &aKSpace, const SurfelAdjacency< KSpace::dimension > &aSurfelAdj, const PointPredicate &pp, bool forceOrientCellExterior=false)
const Point & upperBound() const
Trace trace(traceWriterTerm)
std::ostream & info()
const Domain & domain() const
const Point & lowerBound() const
std::ostream & error()