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" 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" 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" 62 #include <boost/program_options/options_description.hpp> 63 #include <boost/program_options/parsers.hpp> 64 #include <boost/program_options/variables_map.hpp> 69 using namespace DGtal;
72 namespace po = boost::program_options;
134 void missingParam ( std::string param )
136 trace.error() <<
" Parameter: "<<param<<
" is required..";
137 trace.info() <<std::endl;
142 int main (
int argc,
char**argv )
146 po::options_description general_opt (
"Allowed options are: " );
147 general_opt.add_options()
148 (
"help,h",
"display this message." )
149 (
"input,i", po::value<std::string>(),
"Input vol file." )
150 (
"output,o", po::value<string>(),
"Output filename." )
151 (
"level,l", po::value<unsigned int>()->default_value ( 0 ),
"Iso-level for the surface construction." )
152 (
"sigma,s", po::value<double>()->default_value ( 5.0 ),
"Sigma parameter of the Gaussian kernel." )
153 (
"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).")
154 (
"vectorsNorm,N", po::value<double>()->default_value(1.0),
"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.")
155 (
"neighborhood,n", po::value<unsigned int>()->default_value ( 10 ),
"Size of the neighborhood for the convolution (distance on surfel graph)." );
158 po::variables_map vm;
160 po::store(po::parse_command_line(argc, argv, general_opt), vm);
161 }
catch(
const std::exception& ex){
163 trace.info()<<
"Error checking program options: "<< ex.what()<< endl;
167 if (!parseOK || vm.count (
"help" ) ||argc<=1 )
169 trace.info() <<
"Generate normal vector field from a vol file using DGtal library."<<std::endl
170 <<
"It will output the embedded vector field (Gaussian convolution on elementary normal vectors)"<<std::endl
171 <<
"an OFF file, and a TXT normal vector file (theta, phi in degree)." 172 << std::endl <<
"Basic usage: "<<std::endl
173 <<
"\tvol2normalField[options] --input <volFileName> --o <outputFileName> "<<std::endl
174 << general_opt <<
"\n";
179 if ( ! ( vm.count (
"input" ) ) ) missingParam (
"--input" );
180 std::string filename = vm[
"input"].as<std::string>();
181 if ( ! ( vm.count (
"output" ) ) ) missingParam (
"--output" );
182 std::string outputFileName = vm[
"output"].as<std::string>();
184 unsigned int level = vm[
"level"].as<
unsigned int>();
185 double sigma = vm[
"sigma"].as<
double>();
186 unsigned int neighborhood = vm[
"neighborhood"].as<
unsigned int>();
187 double normExport = vm[
"vectorsNorm"].as<
double>();
188 typedef ImageSelector < Z3i::Domain, unsigned char>::Type Image;
189 Image image = VolReader<Image>::importVol ( filename );
191 trace.info() <<image<<std::endl;
193 functors::SimpleThresholdForegroundPredicate<Image> simplePredicate ( image, level );
196 bool space_ok = ks.init ( image.domain().lowerBound(),
197 image.domain().upperBound(), true );
200 trace.error() <<
"Error in the Khamisky space construction."<<std::endl;
204 typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
205 MySurfelAdjacency surfAdj (
true );
208 typedef LightImplicitDigitalSurface<KSpace, functors::SimpleThresholdForegroundPredicate<Image> > MyDigitalSurfaceContainer;
209 typedef DigitalSurface<MyDigitalSurfaceContainer> MyDigitalSurface;
210 SCell bel = Surfaces<KSpace>::findABel ( ks, simplePredicate );
212 MyDigitalSurfaceContainer* ptrSurfContainer =
213 new MyDigitalSurfaceContainer ( ks, simplePredicate, surfAdj, bel );
214 MyDigitalSurface digSurf ( ptrSurfContainer );
215 MyDigitalSurface::ConstIterator it = digSurf.begin();
218 typedef CanonicDigitalSurfaceEmbedder<MyDigitalSurface> SurfaceEmbedder;
219 SurfaceEmbedder surfaceEmbedder ( digSurf );
222 deprecated::GaussianConvolutionWeights < MyDigitalSurface::Size > Gkernel ( sigma );
225 typedef deprecated::LocalConvolutionNormalVectorEstimator < MyDigitalSurface,
226 deprecated::GaussianConvolutionWeights< MyDigitalSurface::Size> > MyGaussianEstimator;
227 BOOST_CONCEPT_ASSERT ( ( concepts::CNormalVectorEstimator< MyGaussianEstimator > ) );
228 MyGaussianEstimator myNormalEstimatorG ( digSurf, Gkernel );
231 typedef DigitalSurfaceEmbedderWithNormalVectorEstimator<SurfaceEmbedder,MyGaussianEstimator> SurfaceEmbedderWithGaussianNormal;
232 SurfaceEmbedderWithGaussianNormal mySurfelEmbedderG ( surfaceEmbedder, myNormalEstimatorG );
235 myNormalEstimatorG.init ( 1.0, neighborhood );
237 trace.info() <<
"Generating the NOFF surface "<< std::endl;
238 ofstream out2 ( ( outputFileName +
".off" ).c_str() );
240 digSurf.exportAs3DNOFF ( out2 ,mySurfelEmbedderG );
243 trace.info() <<
"Generating the polar coordinates file"<< std::endl;
244 ofstream out3 ( ( outputFileName +
".txt" ).c_str() );
247 MyGaussianEstimator::Quantity res;
248 for ( MyDigitalSurface::ConstIterator it =digSurf.begin(),
249 itend = digSurf.end(); it != itend; ++it )
251 res = myNormalEstimatorG.eval ( it );
253 out3<< acos ( res [2] ) *180.0/M_PI <<
" " << ( atan2 ( res [1], res [0] ) + M_PI ) *180.0/M_PI;
255 if (vm.count(
"exportOriginAndExtremity"))
258 out3 <<
" " << mySurfelEmbedderG(*it)[0]
259 <<
" " << mySurfelEmbedderG(*it)[1]
260 <<
" " << mySurfelEmbedderG(*it)[2] <<
" " 261 << mySurfelEmbedderG(*it)[0]+res[0] <<
" " 262 << mySurfelEmbedderG(*it)[1]+res[1] <<
" " 263 << mySurfelEmbedderG(*it)[2]+res[2];