DGtalTools  1.2.0
homotopicThinning3D.cpp
1 
30 #include <iostream>
31 #include <queue>
32 #include <QImageReader>
33 
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"
40 
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"
45 
46 #include "CLI11.hpp"
47 
48 #include "DGtal/geometry/volumes/distance/DistanceTransformation.h"
49 
51 
52 using namespace std;
53 using namespace DGtal;
54 using namespace Z3i;
55 
56 
58 
107 void missingParam ( std::string param )
108 {
109  trace.error() <<" Parameter: "<<param<<" is required..";
110  trace.info() <<std::endl;
111  exit ( 1 );
112 }
113 
114 
115 int main( int argc, char** argv )
116 {
117  // parse command line using CLI ----------------------------------------------
118  CLI::App app;
119  std::string inputFileName;
120  std::string exportSDPName;
121  std::string fixedPointSDPName;
122  double min {0};
123  double max {255};
124 
125  std::vector<int> vectFixedPoints;
126 
127  app.description("Applies an homotopic thinning of a 3d image file (vol,longvol,pgm3d...) with 3D viewer. \n Basic usage: \t homotopicThinning3d [options] --input <3dImageFileName> {vol,longvol,pgm3d...} \n Usage by forcing point to be left by the thinning: \n 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");
128  app.add_option("-i,--input,1", inputFileName, "Input volumetric file (.vol, .pgm3d or p3d)" )
129  ->required()
130  ->check(CLI::ExistingFile);
131  app.add_option("--min,-m", min, "Minimum (excluded) value for threshold.", true);
132  app.add_option("--max,-M", max, "Maximum (included) value for threshold.", true);
133  app.add_option("--exportSDP,-e",exportSDPName, "Export the resulting set of points in a simple (sequence of discrete point (sdp)).");
134  app.add_option("--fixedPoints",vectFixedPoints, "defines the coordinates of points which should not be removed.");
135  app.add_option("--fixedPointSDP,-s",fixedPointSDPName, "use fixed points from a file.")
136  ->check(CLI::ExistingFile);
137 
138 
139  app.get_formatter()->column_width(40);
140  CLI11_PARSE(app, argc, argv);
141  // END parse command line using CLI ----------------------------------------------
142 
143 
144  typedef ImageSelector < Z3i::Domain, unsigned char>::Type Image;
145  Image image = GenericReader<Image>::import ( inputFileName );
146 
147  trace.beginBlock("DT Computation");
148  typedef functors::IntervalForegroundPredicate<Image> Predicate;
149  Predicate aPredicate(image, min, max );
150 
151  const Z3i::L2Metric aMetric{};
152  DistanceTransformation<Z3i::Space, Predicate , Z3i::L2Metric> dt(image.domain(), aPredicate, aMetric );
153  trace.endBlock();
154  trace.info() <<image<<std::endl;
155 
156  // Domain creation from two bounding points.
157  trace.beginBlock("Constructing Set");
158  DigitalSet shape_set( image.domain() );
159  DigitalSet fixedSet( image.domain() );
160 
161  // Get the optional fixed points
162  if(vectFixedPoints.size() > 0){
163  if(vectFixedPoints.size()%3==0){
164  for( unsigned int i=0; i < vectFixedPoints.size()-2; i=i+3)
165  {
166  Z3i::Point pt(vectFixedPoints.at(i), vectFixedPoints.at(i+1), vectFixedPoints.at(i+2));
167  fixedSet.insertNew(pt);
168  }
169  }
170  else
171  {
172  trace.error()<< " The coordinates should be 3d coordinates, ignoring fixedPoints option." << std::endl;
173  }
174  }
175 
176  if(fixedPointSDPName != "")
177  {
178  std::vector<Z3i::Point> vPt = PointListReader<Z3i::Point>::getPointsFromFile(fixedPointSDPName);
179  for( auto &p: vPt)
180  {
181  fixedSet.insert(p);
182  }
183  }
184 
185  SetFromImage<DigitalSet>::append<Image>(shape_set, image, min, max );
186  trace.info() << shape_set<<std::endl;
187  trace.endBlock();
188 
189  trace.beginBlock("Computing skeleton");
190  // (6,18), (18,6), (26,6) seem ok.
191  // (6,26) gives sometimes weird results (but perhaps ok !).
192  Object26_6 shape( dt26_6, shape_set );
193  int nb_simple=0;
194  int layer = 1;
195  std::queue<DigitalSet::Iterator> Q;
196  do
197  {
198  trace.info() << "Layer: "<< layer << std::endl;
199  int nb=0;
200  DigitalSet & S = shape.pointSet();
201 
202  trace.progressBar(0, (double)S.size());
203  for ( DigitalSet::Iterator it = S.begin(); it != S.end(); ++it )
204  {
205  if ( nb % 100 == 0 ) trace.progressBar((double)nb, (double)S.size());
206  nb++;
207  if (dt( *it ) <= layer)
208  {
209  if ( shape.isSimple( *it ) )
210  Q.push( it );
211  }
212  }
213  trace.progressBar( (double)S.size(), (double)S.size() );
214  nb_simple = 0;
215  while ( ! Q.empty() )
216  {
217  DigitalSet::Iterator it = Q.front();
218  Q.pop();
219  if ( shape.isSimple( *it ) && fixedSet.find(*it) == fixedSet.end() )
220  {
221  S.erase( *it );
222  ++nb_simple;
223  }
224  }
225  trace.info() << "Nb simple points : "<<nb_simple<< " " << std::endl;
226  ++layer;
227  }
228  while ( nb_simple != 0 );
229  trace.endBlock();
230 
231  DigitalSet & S = shape.pointSet();
232 
233  trace.info() << "Skeleton--> "<<S<<std::endl;
234 
235  // Display by using two different list to manage OpenGL transparency.
236  QApplication application(argc,argv);
237  Viewer3D<> viewer;
238  viewer.setWindowTitle("homotopicThinning3D");
239  viewer.show();
240 
241  viewer << SetMode3D( shape_set.className(), "Paving" );
242  viewer << CustomColors3D(Color(25,25,255, 255), Color(25,25,255, 255));
243  viewer << S ;
244  viewer << CustomColors3D(Color(255,25,255, 255), Color(255,25,255, 255));
245  viewer << fixedSet;
246  viewer << SetMode3D( shape_set.className(), "PavingTransp" );
247  viewer << CustomColors3D(Color(250, 0,0, 25), Color(250, 0,0, 5));
248  viewer << shape_set;
249  viewer<< Viewer3D<>::updateDisplay;
250 
251  if (exportSDPName != "")
252  {
253  std::ofstream out;
254  out.open(exportSDPName.c_str());
255  for (auto &p : S)
256  {
257  out << p[0] << " " << p[1] << " " << p[2] << std::endl;
258  }
259  }
260  return application.exec();
261 
262 }
263 // //
265 
266 
Definition: ATu0v1.h:57