DGtal  1.3.beta
VoxelComplexThinning.h
1 
17 #pragma once
18 
30 // Inclusions
32 #include <iostream>
33 #include <chrono>
34 #include "DGtal/base/Common.h"
35 #include "DGtal/topology/VoxelComplex.h"
36 #include "DGtal/topology/VoxelComplexFunctions.h"
37 #include "DGtal/geometry/volumes/distance/DistanceTransformation.h"
38 #include "DGtal/geometry/volumes/distance/ExactPredicateLpSeparableMetric.h"
40 namespace DGtal {
41 namespace functions {
42 /*
43 Get a skeletonized or thinned image from a binary image set.
44 
45 Parameters:
46 ----------
47 complex: TComplex
48  input complex to thin
49 
50 skel_type: str
51  Voxels to keep in the skeletonization process.
52 
53  [end, ulti, isthmus]
54  - end: keep end voxels.
55  - ulti: don't keep extra voxels, ultimate skeleton.
56  - isthmus: keep voxels that are isthmuses.
57 
58 select_type: str
59  [first, random, dmax]
60  - first: the first voxel in the set (no criteria)
61  - random: random voxel in the set.
62  - dmax: choose voxel with greatest distance map value.
63  Strategy to choose voxels in the asymmetric process.
64 
65 table_folder: str
66  Location of the DGtal look-up-tables for simplicity and isthmusicity,
67  for example simplicity_table26_6.zlib.
68  These tables are distributed with the sgext package. Use the variable
69  'sgext.tables_folder'.
70 
71 persistence: int
72  if >0, performs a persistence algorithm that prunes
73  branches that are not persistant (less important).
74 
75 distance_transform: str
76  file holding a distance map. Required for select_type dmax option.
77  This option provides a centered skeleton.
78 
79 profile: bool
80  time the algorithm
81 
82 verbose: bool
83  extra information displayed during the algorithm.
84 
85 Return
86 ------
87 A new thinned voxel complex.
88 */
89 template<typename TComplex,
90  typename TDistanceTransform =
91  DistanceTransformation<Z3i::Space, Z3i::DigitalSet, ExactPredicateLpSeparableMetric<Z3i::Space, 3>> >
93  TComplex & vc,
94  const std::string & skel_type_str,
95  const std::string & skel_select_type_str,
96  const std::string & tables_folder,
97  const int & persistence = 0,
98  const TDistanceTransform * distance_transform = nullptr,
99  const bool profile = false,
100  const bool verbose = false)
101 {
102  if(verbose) {
103  using DGtal::trace;
104  trace.beginBlock("thin_function parameters:");
105  trace.info() << "skel_type_str: " << skel_type_str << std::endl;
106  trace.info() << "skel_select_type_str: " << skel_select_type_str << std::endl;
107  if(distance_transform) {
108  trace.info() << " -- provided distance_transform." << std::endl;
109  }
110  trace.info() << "persistence: " << persistence << std::endl;
111  trace.info() << "profile: " << profile << std::endl;
112  trace.info() << "verbose: " << verbose << std::endl;
113  trace.info() << "----------" << std::endl;
114  trace.endBlock();
115  }
116 
117  // Validate input skel method and skel_select
118  const bool skel_type_str_is_valid =
119  skel_type_str == "ultimate" ||
120  skel_type_str == "end" ||
121  skel_type_str == "isthmus" ||
122  skel_type_str == "1isthmus" ||
123  skel_type_str == "isthmus1";
124  if(!skel_type_str_is_valid) {
125  throw std::runtime_error("skel_type_str is not valid: \"" + skel_type_str + "\"");
126  }
127 
128  const bool skel_select_type_str_is_valid =
129  skel_select_type_str == "first" ||
130  skel_select_type_str == "random" ||
131  skel_select_type_str == "dmax";
132  if(!skel_select_type_str_is_valid) {
133  throw std::runtime_error("skel_select_type_str is not valid: \"" + skel_select_type_str + "\"");
134  }
135  // // No filesystem to check the tables folder exist. If incorrect, it will fail loading the table.
136  // const fs::path tables_folder_path{tables_folder};
137  // if(!fs::exists(tables_folder_path)) {
138  // throw std::runtime_error("tables_folder " + tables_folder_path.string() +
139  // " doesn't exist in the filesystem.\n"
140  // "tables_folder should point to the folder "
141  // "where DGtal tables are: i.e simplicity_table26_6.zlib");
142  // }
143 
144  // Create a VoxelComplex from the set
145  using KSpace = DGtal::Z3i::KSpace;
146  using Complex = TComplex;
147  using ComplexCell = typename Complex::Cell;
148  using ComplexClique = typename Complex::Clique;
149  using Point = DGtal::Z3i::Point;
150 
151  if(verbose) { DGtal::trace.beginBlock("load isthmus table"); }
152  boost::dynamic_bitset<> isthmus_table;
153  auto &sk = skel_type_str;
154  if(sk == "isthmus") {
155  const std::string tableIsthmus = tables_folder + "/isthmusicity_table26_6.zlib";
156  isthmus_table = *DGtal::functions::loadTable(tableIsthmus);
157  } else if(sk == "isthmus1" || sk == "1ishtmus") {
158  const std::string tableOneIsthmus = tables_folder + "/isthmusicityOne_table26_6.zlib";
159  isthmus_table = *DGtal::functions::loadTable(tableOneIsthmus);
160  }
161  if(verbose) { DGtal::trace.endBlock(); }
162 
163 
164  // SKEL FUNCTION:
165  // Load a look-up-table for the neighborgood of a point
166  auto pointMap =
167  *DGtal::functions::mapZeroPointNeighborhoodToConfigurationMask<Point>();
168  std::function<bool(const Complex&, const ComplexCell&)> Skel;
169  if(sk == "ultimate") {
170  Skel = DGtal::functions::skelUltimate<Complex>;
171  } else if(sk == "end") {
172  Skel = DGtal::functions::skelEnd<Complex>;
173  // else if (sk == "1is") Skel = oneIsthmus<Complex>;
174  // else if (sk == "is") Skel = skelIsthmus<Complex>;
175  } else if(sk == "isthmus1" || sk == "1ishtmus") {
176  Skel = [&isthmus_table, &pointMap](const Complex& fc,
177  const ComplexCell& c) {
178  return DGtal::functions::skelWithTable(isthmus_table, pointMap, fc, c);
179  };
180  } else if(sk == "isthmus") {
181  Skel = [&isthmus_table, &pointMap](const Complex& fc,
182  const ComplexCell& c) {
183  return DGtal::functions::skelWithTable(isthmus_table, pointMap, fc, c);
184  };
185  } else {
186  throw std::runtime_error("Invalid skel string");
187  }
188 
189  // SELECT FUNCTION
190  std::function<std::pair<typename Complex::Cell, typename Complex::Data>(
191  const ComplexClique&)>
192  Select;
193 
194  // profile
195  auto start = std::chrono::system_clock::now();
196 
197  // If dmax is chosen, but no distance_transform is provided, create one.
198  using DT = TDistanceTransform;
199  using DTDigitalSet = typename DT::PointPredicate;
200  using DTDigitalSetDomain = typename DTDigitalSet::Domain;
201  using Metric = typename DT::SeparableMetric;
202  Metric l3;
203  auto &sel = skel_select_type_str;
204  const bool compute_distance_transform = !distance_transform && sel == "dmax";
205  DTDigitalSetDomain vc_domain = compute_distance_transform ?
206  DTDigitalSetDomain(vc.space().lowerBound(), vc.space().upperBound()) :
207  // dummy domain
208  DTDigitalSetDomain(Z3i::Point::zero, Z3i::Point::diagonal(1));
209 
210  DTDigitalSet image_set = DTDigitalSet(vc_domain);
211  if(compute_distance_transform) {
212  vc.dumpVoxels(image_set);
213  }
214  // Create the distance map here (computationally expensive if not dummy).
215  DT dist_map(vc_domain, image_set, l3);
216  if(compute_distance_transform) {
217  distance_transform = &dist_map;
218  }
219 
220  if(sel == "random") {
221  Select = DGtal::functions::selectRandom<Complex>;
222  } else if(sel == "first") {
223  Select = DGtal::functions::selectFirst<Complex>;
224  } else if(sel == "dmax") {
225  Select = [&distance_transform](const ComplexClique& clique) {
226  return selectMaxValue<TDistanceTransform, Complex>(*distance_transform, clique);
227  };
228  } else {
229  throw std::runtime_error("Invalid skel select type");
230  }
231 
232  // Perform the thin/skeletonization
233  Complex vc_new(vc.space());
234  if(persistence == 0) {
235  vc_new = DGtal::functions::asymetricThinningScheme<Complex>(vc, Select, Skel, verbose);
236  } else {
237  vc_new = DGtal::functions::persistenceAsymetricThinningScheme<Complex>(vc, Select, Skel,
238  persistence, verbose);
239  }
240 
241  // profile
242  auto end = std::chrono::system_clock::now();
243  auto elapsed = std::chrono::duration_cast<std::chrono::seconds>(end - start);
244  if(profile) {
245  std::cout << "Time elapsed: " << elapsed.count() << std::endl;
246  }
247 
248  return vc_new;
249 }
250 
251 } // namespace functions
252 } // namespace DGtal
DGtal::PointVector< dim, Integer >::zero
static Self zero
Static const for zero PointVector.
Definition: PointVector.h:1595
DGtal::Trace::endBlock
double endBlock()
DGtal::functions::thinningVoxelComplex
TComplex thinningVoxelComplex(TComplex &vc, const std::string &skel_type_str, const std::string &skel_select_type_str, const std::string &tables_folder, const int &persistence=0, const TDistanceTransform *distance_transform=nullptr, const bool profile=false, const bool verbose=false)
Definition: VoxelComplexThinning.h:92
DGtal::trace
Trace trace
Definition: Common.h:154
DGtal::Z3i::Point
Space::Point Point
Definition: StdDefs.h:168
DGtal::Z3i::KSpace
KhalimskySpaceND< 3, Integer > KSpace
Definition: StdDefs.h:146
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::PointVector< dim, Integer >::diagonal
static Self diagonal(Component val=1)
DGtal::functions::loadTable
DGtal::CountedPtr< boost::dynamic_bitset<> > loadTable(const std::string &input_filename, const unsigned int known_size, const bool compressed=true)
DGtal::Trace::info
std::ostream & info()
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
Domain
HyperRectDomain< Space > Domain
Definition: testSimpleRandomAccessRangeFromPoint.cpp:44
DGtal::functions::skelWithTable
bool skelWithTable(const boost::dynamic_bitset<> &table, const std::unordered_map< typename TComplex::Point, unsigned int > &pointToMaskMap, const TComplex &vc, const typename TComplex::Cell &cell)
Cell
KSpace::Cell Cell
Definition: testCubicalComplex.cpp:56
Point
MyPointD Point
Definition: testClone2.cpp:383
DGtal::KhalimskySpaceND
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
Definition: KhalimskySpaceND.h:64