32 #include <QImageReader> 34 #include "DGtal/io/viewers/Viewer3D.h" 35 #include "DGtal/io/DrawWithDisplay3DModifier.h" 36 #include "DGtal/io/Color.h" 37 #include "DGtal/shapes/Shapes.h" 38 #include "DGtal/helpers/StdDefs.h" 39 #include "DGtal/io/readers/PointListReader.h" 41 #include "DGtal/io/readers/GenericReader.h" 42 #include "DGtal/images/imagesSetsUtils/SetFromImage.h" 43 #include "DGtal/images/SimpleThresholdForegroundPredicate.h" 44 #include "DGtal/images/ImageSelector.h" 46 #include <boost/program_options/options_description.hpp> 47 #include <boost/program_options/parsers.hpp> 48 #include <boost/program_options/variables_map.hpp> 50 #include "DGtal/geometry/volumes/distance/DistanceTransformation.h" 55 using namespace DGtal;
58 namespace po = boost::program_options;
106 void missingParam ( std::string param )
108 trace.error() <<
" Parameter: "<<param<<
" is required..";
109 trace.info() <<std::endl;
115 int main(
int argc,
char** argv )
119 po::options_description general_opt (
"Allowed options are: " );
120 general_opt.add_options()
121 (
"help,h",
"display this message." )
122 (
"input,i", po::value<std::string>(),
"Input volumetric file (.vol, .pgm3d or p3d)" )
123 (
"min,m", po::value<int>()->default_value( 0 ),
"Minimum (excluded) value for threshold.")
124 (
"max,M", po::value<int>()->default_value( 255 ),
"Maximum (included) value for threshold.")
125 (
"exportSDP,e", po::value<std::string>(),
"Export the resulting set of points in a simple (sequence of discrete point (sdp)).")
126 (
"fixedPoints", po::value<std::vector <int> >()->multitoken(),
"defines the coordinates of points which should not be removed." )
127 (
"fixedPointSDP,s", po::value<std::string>(),
"use fixed points from a file.");
131 po::variables_map vm;
133 po::store(po::parse_command_line(argc, argv, general_opt), vm);
134 }
catch(
const std::exception& ex){
136 trace.info()<<
"Error checking program options: "<< ex.what()<< endl;
139 if ( !parseOK || vm.count (
"help" ) ||argc<=1 )
141 trace.info() <<
"Illustration of homotopic thinning of a 3d image file (vol,longvol,pgm3d...) with 3D viewer."<<std::endl
142 << std::endl <<
"Basic usage: "<<std::endl
143 <<
"\thomotopicThinning3d [options] --input <3dImageFileName> {vol,longvol,pgm3d...} "<<std::endl
144 << general_opt <<
"\n" 145 <<
" Usage by forcing point to be left by the thinning: \n" 146 <<
"homotopicThinning3D --input ${DGtal}/examples/samples/Al.100.vol --fixedPoints 56 35 5 56 61 5 57 91 38 58 8 38 45 50 97 \n";
155 if ( ! ( vm.count (
"input" ) ) ) missingParam (
"--input" );
156 std::string filename = vm[
"input"].as<std::string>();
159 typedef ImageSelector < Z3i::Domain, unsigned char>::Type Image;
160 Image image = GenericReader<Image>::import ( filename );
162 trace.beginBlock(
"DT Computation");
163 typedef functors::IntervalForegroundPredicate<Image> Predicate;
164 Predicate aPredicate(image, vm[
"min" ].as<int>(), vm[
"max" ].as<int>() );
166 const Z3i::L2Metric aMetric{};
167 DistanceTransformation<Z3i::Space, Predicate , Z3i::L2Metric> dt(image.domain(), aPredicate, aMetric );
169 trace.info() <<image<<std::endl;
173 trace.beginBlock(
"Constructing Set");
174 DigitalSet shape_set( image.domain() );
175 DigitalSet fixedSet( image.domain() );
178 if( vm.count(
"fixedPoints")){
179 std::vector<int> vectC = vm[
"fixedPoints"].as<std::vector<int> >();
180 if(vectC.size()%3==0){
181 for(
unsigned int i=0; i < vectC.size()-2; i=i+3)
183 Z3i::Point pt(vectC.at(i), vectC.at(i+1), vectC.at(i+2));
184 fixedSet.insertNew(pt);
187 trace.error()<<
" The coordinates should be 3d coordinates, ignoring fixedPoints option." << std::endl;
190 if(vm.count(
"fixedPointSDP"))
192 std::vector<Z3i::Point> vPt = PointListReader<Z3i::Point>::getPointsFromFile(vm[
"fixedPointSDP"].as<std::string>());
199 SetFromImage<DigitalSet>::append<Image>(shape_set, image,
200 vm[
"min" ].as<
int>(), vm[
"max" ].as<int>() );
201 trace.info() << shape_set<<std::endl;
207 trace.beginBlock(
"Computing skeleton");
210 Object26_6 shape( dt26_6, shape_set );
213 std::queue<DigitalSet::Iterator> Q;
216 trace.info() <<
"Layer: "<< layer << std::endl;
218 DigitalSet & S = shape.pointSet();
220 trace.progressBar(0, (
double)S.size());
221 for ( DigitalSet::Iterator it = S.begin(); it != S.end(); ++it )
223 if ( nb % 100 == 0 ) trace.progressBar((
double)nb, (
double)S.size());
225 if (dt( *it ) <= layer)
227 if ( shape.isSimple( *it ) )
231 trace.progressBar( (
double)S.size(), (double)S.size() );
233 while ( ! Q.empty() )
235 DigitalSet::Iterator it = Q.front();
237 if ( shape.isSimple( *it ) && fixedSet.find(*it) == fixedSet.end() )
243 trace.info() <<
"Nb simple points : "<<nb_simple<<
" " << std::endl;
246 while ( nb_simple != 0 );
249 DigitalSet & S = shape.pointSet();
251 trace.info() <<
"Skeleton--> "<<S<<std::endl;
254 QApplication application(argc,argv);
256 viewer.setWindowTitle(
"homotopicThinning3D");
259 viewer << SetMode3D( shape_set.className(),
"Paving" );
260 viewer << CustomColors3D(Color(25,25,255, 255), Color(25,25,255, 255));
262 viewer << CustomColors3D(Color(255,25,255, 255), Color(255,25,255, 255));
264 viewer << SetMode3D( shape_set.className(),
"PavingTransp" );
265 viewer << CustomColors3D(Color(250, 0,0, 25), Color(250, 0,0, 5));
268 viewer<< Viewer3D<>::updateDisplay;
270 if (vm.count(
"exportSDP"))
273 out.open(vm[
"exportSDP"].as<std::string>().c_str());
276 out << p[0] <<
" " << p[1] <<
" " << p[2] << std::endl;
279 return application.exec();