DGtalTools  1.2.0
criticalKernelsThinning3D.cpp
1 
75 #include <iostream>
76 #include <chrono>
77 #include <unordered_map>
78 
79 #include <DGtal/io/viewers/Viewer3D.h>
80 #include <DGtal/base/Common.h>
81 #include <DGtal/helpers/StdDefs.h>
82 #include <DGtal/io/readers/GenericReader.h>
83 #include <DGtal/io/writers/GenericWriter.h>
84 #include "DGtal/images/imagesSetsUtils/SetFromImage.h"
85 #include "DGtal/images/SimpleThresholdForegroundPredicate.h"
86 #include "DGtal/images/ImageSelector.h"
87 #include "DGtal/images/imagesSetsUtils/ImageFromSet.h"
88 
89 #include <DGtal/topology/SurfelAdjacency.h>
90 #include <DGtal/io/boards/Board2D.h>
91 #include <DGtal/topology/CubicalComplex.h>
92 #include <DGtal/topology/CubicalComplexFunctions.h>
93 
94 #include <DGtal/topology/VoxelComplex.h>
95 #include <DGtal/topology/VoxelComplexFunctions.h>
96 #include "DGtal/topology/NeighborhoodConfigurations.h"
97 #include "DGtal/topology/tables/NeighborhoodTables.h"
98 
99 // Distance Map
100 #include "DGtal/geometry/volumes/distance/ExactPredicateLpSeparableMetric.h"
101 #include "DGtal/geometry/volumes/distance/VoronoiMap.h"
102 #include "DGtal/geometry/volumes/distance/DistanceTransformation.h"
103 
104 #include "CLI11.hpp"
105 
106 
107 using namespace DGtal;
108 using namespace std;
109 using namespace DGtal::Z3i;
110 
111 int main(int argc, char* const argv[]){
112 
113  // parse command line using CLI ----------------------------------------------
114  CLI::App app;
115  std::string inputFileName;
116  std::string sk_string {"1isthmus"};
117  string select_string {"dmax"};
118  string foreground {"black"};
119  string outputFilenameImg;
120  string outputFilenameSDP;
121 
122  int thresholdMin {0};
123  int thresholdMax {255};
124  int persistence {0};
125  bool profile {false};
126  bool verbose {false};
127  bool visualize {false};
128 
129  app.description("Compute the thinning of a volume using the CriticalKernels framework\nBasic usage: criticalKernelsThinning3D --input <volFileName> --skel <ulti,end, 1isthmus, isthmus> --select [ -f <white,black> -m <minlevel> -M <maxlevel> -v ] [--persistence <value> ] --persistence <value> ] \n options for --skel {ulti end 1isthmus isthmus} \n options for --select = {dmax random first} \n Example: \n criticalKernelsThinning3D --input ${DGtal}/examples/samples/Al.100.vol --select dmax --skel 1isthmus --persistence 1 --visualize --verbose --exportImage ./Al100_dmax_1isthmus_p1.vol \n");
130  app.add_option("-i,--input,1", inputFileName, "Input vol file." )
131  ->required()
132  ->check(CLI::ExistingFile);
133 
134  app.add_option("--skel,-s", sk_string,"Type of skeletonization. Options: 1isthmus, isthmus, end, ulti.", true )
135  -> check(CLI::IsMember({"ulti", "end","isthmus", "1isthmus"}));
136  app.add_option("--select,-c", select_string, "Select the ordering for skeletonization. Options: dmax, random, first", true)
137  -> check(CLI::IsMember({"random", "dmax", "first"}));
138  app.add_option("--foreground,-f",foreground, "Foreground color in binary image", true )
139  -> check(CLI::IsMember({"white", "black"}));
140  app.add_option("--thresholdMin,-m", thresholdMin, "Threshold min (excluded) to define binary shape", true );
141  app.add_option("--thresholdMax,-M", thresholdMax, "Threshold max (included) to define binary shape", true );
142  app.add_option("--persistence,-p",persistence,"Persistence value, implies use of persistence algorithm if p>=1", true )
143  ->check(CLI::PositiveNumber);
144  app.add_flag("--profile", profile, "Profile algorithm");
145  app.add_flag("--verbose,-v",verbose, "Verbose output");
146  app.add_option("--exportImage,-o",outputFilenameImg, "Export the resulting set of points to a image compatible with GenericWriter.");
147  app.add_option("--exportSDP,-e",outputFilenameSDP, "Export the resulting set of points in a simple (sequence of discrete point (sdp))." );
148  app.add_flag("--visualize,-t", visualize, "Visualize result in viewer");
149 
150  app.get_formatter()->column_width(40);
151  CLI11_PARSE(app, argc, argv);
152  // END parse command line using CLI ----------------------------------------------
153 
154 
155 
156  if(verbose){
157  std::cout << "Skel: " << sk_string << std::endl;
158  std::cout << "Select: " << select_string << std::endl;
159  std::cout << "Persistence: " << persistence << std::endl;
160  std::cout << "Input: " << inputFileName << std::endl;
161  }
162  trace.beginBlock("Reading input");
163  using Domain = Z3i::Domain ;
164  using Image = ImageSelector < Z3i::Domain, unsigned char>::Type ;
165  Image image = GenericReader<Image>::import(inputFileName);
166  trace.endBlock();
167 
168  DigitalSet image_set (image.domain());
169  SetFromImage<Z3i::DigitalSet>::append<Image>(
170  image_set, image,
171  thresholdMin, thresholdMax);
172 
173 
174  // Create a VoxelComplex from the set
175 
176  using DigitalTopology = DT26_6;
177  using DigitalSet =
178  DGtal::DigitalSetByAssociativeContainer<Domain ,
179  std::unordered_set< typename Domain::Point> >;
180  using Complex = DGtal::VoxelComplex<KSpace>;
181 
182  auto & sk = sk_string;
183  KSpace ks;
184  KSpace::Point d1( KSpace::Point::diagonal( 1 ) );
185  ks.init(image.domain().lowerBound() - d1 ,
186  image.domain().upperBound() + d1 , true);
187 
188  trace.beginBlock("construct with table");
189  Complex vc(ks);
190  vc.construct(image_set, functions::loadTable(simplicity::tableSimple26_6 ));
191  trace.endBlock();
192  trace.beginBlock("load isthmus table");
193  boost::dynamic_bitset<> isthmus_table;
194  if (sk == "isthmus")
195  isthmus_table = *functions::loadTable(isthmusicity::tableIsthmus);
196  else if (sk == "1isthmus")
197  isthmus_table = *functions::loadTable(isthmusicity::tableOneIsthmus);
198  auto pointMap = *functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
199 
200  trace.endBlock();
201  using namespace DGtal::functions ;
202  // SKEL FUNCTION:
203  std::function< bool(const Complex&, const Cell&) > Skel ;
204  if (sk == "ulti") Skel = skelUltimate<Complex>;
205  else if (sk == "end") Skel = skelEnd<Complex>;
206  else if (sk == "isthmus" || sk == "1isthmus")
207  Skel = [&isthmus_table, &pointMap](const Complex & fc,
208  const Complex::Cell & c){
209  return skelWithTable(isthmus_table, pointMap, fc, c);
210  };
211  else throw std::runtime_error("Invalid skel string");
212  auto start = std::chrono::system_clock::now();
213 
214  // SELECT FUNCTION
215  /*
216  * Calculate distance map even if not requested:
217  */
218  trace.beginBlock("Create Distance Map");
219  using L3Metric = ExactPredicateLpSeparableMetric<Z3i::Space, 3>;
220  using DT = DistanceTransformation<Z3i::Space, DigitalSet, L3Metric>;
221  L3Metric l3;
222  DT dt(image.domain(), image_set, l3);
223  trace.endBlock();
224 
225  std::function< std::pair<typename Complex::Cell, typename Complex::Data>(const Complex::Clique&) > Select ;
226  auto & sel = select_string;
227  if (sel == "random") Select = selectRandom<Complex>;
228  else if (sel == "first") Select = selectFirst<Complex>;
229  else if (sel == "dmax"){
230  Select =
231  [&dt](const Complex::Clique & clique){
232  return selectMaxValue<DT, Complex>(dt,clique);
233  };
234  } else throw std::runtime_error("Invalid skel string");
235 
236  trace.beginBlock("Thinning");
237  Complex vc_new(ks);
238  if (persistence == 0)
239  vc_new = asymetricThinningScheme< Complex >(
240  vc, Select, Skel, verbose);
241  else
242  vc_new = persistenceAsymetricThinningScheme< Complex >(
243  vc, Select, Skel, persistence, verbose);
244  trace.endBlock();
245 
246  auto end = std::chrono::system_clock::now();
247  auto elapsed = std::chrono::duration_cast<std::chrono::seconds> (end - start) ;
248  if (profile) std::cout <<"Time elapsed: " << elapsed.count() << std::endl;
249 
250 
251  DigitalSet thin_set(image.domain());
252  vc_new.dumpVoxels(thin_set);
253  const auto & all_set = image_set;
254 
255  if (outputFilenameSDP != "")
256  {
257  std::ofstream out;
258  out.open(outputFilenameSDP.c_str());
259  for (auto &p : thin_set)
260  {
261  out << p[0] << " " << p[1] << " " << p[2] << std::endl;
262  }
263  }
264 
265  if (outputFilenameImg != "")
266  {
267  if(verbose)
268  std::cout << "outputFilename" << outputFilenameImg << std::endl;
269 
270  unsigned int foreground_value = 255;
271  auto thin_image = ImageFromSet<Image>::create(thin_set, foreground_value);
272  thin_image >> outputFilenameImg;
273  }
274 
275  if(visualize)
276  {
277  int argc(1);
278  char** argv(nullptr);
279  QApplication app(argc, argv);
280  Viewer3D<> viewer;
281  viewer.setWindowTitle("criticalKernelsThinning3D");
282  viewer.show();
283 
284  viewer.setFillColor(Color(255, 255, 255, 255));
285  viewer << thin_set;
286 
287  // All kspace voxels
288  viewer.setFillColor(Color(40, 200, 55, 10));
289  viewer << all_set;
290 
291  viewer << Viewer3D<>::updateDisplay;
292 
293  app.exec();
294  }
295 }
Definition: ATu0v1.h:57