DGtal  1.0.0
VoxelComplexFunctions.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file VoxelComplexFunctions.ih
19  * @author Pablo Hernandez-Cerdan (\c pablo.hernandez.cerdan@outlook.com)
20  * Institute of Fundamental Sciences. Massey University.
21  * Palmerston North, New Zealand
22  *
23  * @date 2018/01/01
24  *
25  * Implementation of inline methods defined in VoxelComplexFunctions.h
26  *
27  * This file is part of the DGtal library.
28  */
29 
30 
31 //////////////////////////////////////////////////////////////////////////////
32 #include <cstdlib>
33 #include <DGtal/topology/DigitalTopology.h>
34 #include <random>
35 //////////////////////////////////////////////////////////////////////////////
36 
37 ///////////////////////////////////////////////////////////////////////////////
38 // IMPLEMENTATION of inline methods.
39 ///////////////////////////////////////////////////////////////////////////////
40 
41 //-----------------------------------------------------------------------------
42 template < typename TComplex >
43 TComplex
44 DGtal::functions::
45 asymetricThinningScheme(
46  TComplex & vc ,
47  std::function<
48  std::pair<typename TComplex::Cell, typename TComplex::Data>(
49  const typename TComplex::Clique &)
50  > Select ,
51  std::function<
52  bool(
53  const TComplex & ,
54  const typename TComplex::Cell & )
55  > Skel,
56  bool verbose )
57 {
58  if(verbose) trace.beginBlock("Asymetric Thinning Scheme");
59 
60  using Cell = typename TComplex::Cell;
61  using Data = typename TComplex::Data;
62 
63  // Initialize K set, and VoxelComplex containers;
64  struct FirstComparator
65  {
66  bool operator() (const std::pair<Cell, Data>& lhs, const std::pair<Cell, Data>& rhs) const
67  {
68  return lhs.first < rhs.first;
69  }
70  };
71  std::set<std::pair<Cell, Data>, FirstComparator> selected_voxels_per_dim;
72  TComplex already_selected_voxels(vc.space());
73  TComplex constraint_set(vc.space());
74 
75  auto & Z = selected_voxels_per_dim;
76  auto & Y = already_selected_voxels;
77  auto & K = constraint_set;
78  auto X = vc;
79  Y.copySimplicityTable(vc);
80  K.copySimplicityTable(vc);
81  typename TComplex::Parent x_y(vc.space());
82  typename TComplex::Parent x_k(vc.space());
83 
84  bool stability{false};
85  uint64_t generation{0};
86  auto xsize = X.nbCells(3);
87  auto xsize_old = xsize;
88  typename TComplex::CliqueContainer critical_cliques;
89 
90  if(verbose){
91  trace.info() << "generation: " << generation <<
92  " ; X.nbCells(3): " << xsize << std::endl;
93  }
94  do {
95  ++generation;
96 
97  Y = K ;
98  /* Select voxels from critical d-cliques,
99  * having priority d=3 down to 0 using Select(clique)
100  * which returns the selected voxel. */
101  // X.close(); // Close K-Space from current voxels.
102  for (int d = 3 ; d >= 0 ; --d) {
103  Z.clear();
104  /* Y is the already_selected voxels (in this generation)
105  * Search critical cliques only in the cells of X-Y,
106  * but which are critical for X. */
107  x_y = X - Y;
108  // x_y.close();
109  critical_cliques = X.criticalCliquesForD(d, x_y);
110  for(auto & clique : critical_cliques){
111  // if (d!=3)
112  // if (! (clique <= (x_y)) ) continue ;
113  auto yinclude = std::find_if(
114  std::begin(clique),
115  std::end(clique),
116  [&Y]( const Cell& c) {
117  return Y.belongs(c);
118  });
119  if (yinclude != std::end(clique)) continue ;
120  Z.insert(Select(clique));
121  } // critical_cliques
122  // Y = Y union Z. From set to VoxelComplex
123  for( const auto & selected_celldata_pair : Z)
124  Y.insertVoxelCell(selected_celldata_pair);
125  } // d-loop
126  X = Y;
127  // X - K is equal to X-Y, which is equal to a Ynew - Yold
128  x_k = X - K;
129  for (auto it = x_k.begin(3), itE = x_k.end(3) ; it != itE ; ++it ){
130  auto new_voxel = it->first ;
131  if( Skel(X, new_voxel) == true){
132  K.insertVoxelCell(new_voxel);
133  // K.insertCell(3, new_voxel);
134  // K.objectSet().insert(K.space().uCoords(new_voxel));
135  }
136  }
137 
138  // Stability Update:
139  xsize = X.nbCells(3);
140  if(xsize == xsize_old) stability = true;
141  xsize_old = xsize;
142 
143  if(verbose){
144  trace.info() << "generation: " << generation <<
145  " ; X.nbCells(3): " << xsize <<
146  " ; K (constrain set): " << K.nbCells(3) <<
147  " ; Y (selected in this generation): " << Y.nbCells(3) <<
148  " ; X - K ie. Ynew - Yold : " << x_k.nbCells(3) <<
149  " ; X - Y : " << x_y.nbCells(3) << std::endl;
150  }
151  } while( !stability );
152 
153  if(verbose)
154  trace.endBlock();
155 
156  return X;
157 }
158 
159 
160 template < typename TComplex >
161 TComplex
162 DGtal::functions::
163 persistenceAsymetricThinningScheme(
164  TComplex & vc ,
165  std::function<
166  std::pair<typename TComplex::Cell, typename TComplex::Data> (
167  const typename TComplex::Clique &)
168  > Select ,
169  std::function<
170  bool(
171  const TComplex & ,
172  const typename TComplex::Cell & )
173  > Skel,
174  uint32_t persistence,
175  bool verbose )
176 {
177  if(verbose) trace.beginBlock("Persistence asymetricThinningScheme");
178 
179  using Cell = typename TComplex::Cell;
180  using Data = typename TComplex::Data;
181 
182  // Initialize Z set, and VoxelComplex containers;
183  struct FirstComparator
184  {
185  bool operator() (const std::pair<Cell, Data>& lhs, const std::pair<Cell, Data>& rhs) const
186  {
187  return lhs.first < rhs.first;
188  }
189  };
190  TComplex already_selected_voxels(vc.space());
191  TComplex constraint_set(vc.space());
192 
193  auto & Y = already_selected_voxels;
194  auto & K = constraint_set;
195  TComplex X = vc;
196  TComplex x_y(vc.space());
197  Y.copySimplicityTable(vc);
198  K.copySimplicityTable(vc);
199  x_y.copySimplicityTable(vc);
200 
201  bool stability{false};
202  uint64_t generation{0};
203  auto xsize = X.nbCells(3);
204  auto xsize_old = xsize;
205  const bool close_it = true;
206  typename TComplex::CliqueContainer critical_cliques;
207 
208  if(verbose){
209  trace.info() << "Initial Voxels at generation: " << generation <<
210  " ; X.nbCells(3): " << xsize << std::endl;
211  }
212  // cubical cell data (aka, the birth_date) is initialized to zero already.
213  do {
214  ++generation;
215  // Update birth_date for our Skel function. (isIsthmus for example)
216  for (auto it = X.begin(3), itE = X.end(3) ; it != itE ; ++it ){
217  // Ignore voxels existing in K set.(ie: X-K)
218  if (K.findCell(3, it->first) != K.end(3))
219  continue;
220  auto & voxel = it->first ;
221  auto & ccdata = it->second.data ;
222  if( Skel(X, voxel) == true && ccdata == 0)
223  ccdata = generation;
224  }
225  Y = K ;
226  x_y = X; //optimization instead of x_y = X-Y, use x_y -= Y;
227  // d-cliques: From voxels (d=3) to pointels (d=0)
228  for (int d = 3 ; d >= 0 ; --d) {
229  // Search only critical cliques of X in the set X \ Y
230  x_y -= Y; // Y is growing in this d-loop.
231  critical_cliques = X.criticalCliquesForD(d, x_y);
232 
233  for(auto & clique : critical_cliques)
234  {
235  // Find the selected voxel in x to get associated Data.
236  auto selected_voxelpair = Select(clique);
237  selected_voxelpair.second = X.findCell(3,selected_voxelpair.first)->second;
238  Y.insertVoxelCell(selected_voxelpair, close_it);
239  }
240  } // d-loop
241  // Y is now populated with some selected critical cliques for this generation,
242  // plus with the fixed points in K.
243  // Store result (final or for the next generation)
244  X = Y;
245 
246  /*** Modify K (the set of voxels to conserve) ***/
247  // Insert in K (inmutable set) only if the NEW critical cells of this generation (Y-K) are:
248  // a) Part of the Skeleton (Skel) user wants to preserve.
249  // b) persistent enough.
250 
251  // Update K
252  Y -= K;
253  for (auto it = Y.begin(3), itE = Y.end(3) ; it != itE ; ++it ){
254  auto & voxel = it->first;
255  auto & ccdata = it->second.data;
256  bool is_skel = Skel(X, voxel);
257  bool is_persistent_enough = (generation + 1 - ccdata) >= persistence;
258  if (is_skel && is_persistent_enough)
259  K.insertVoxelCell(voxel, close_it, ccdata);
260  }
261 
262  if(verbose){
263  trace.info() << "generation: " << generation <<
264  " ; X at start generation: " << xsize_old <<
265  " ; K (constraint set): " << K.nbCells(3) <<
266  " ; Y - K: " << Y.nbCells(3) <<
267  " ; Y (X at the end): " << X.nbCells(3) <<
268  std::endl;
269  }
270 
271  /***** Stability Update: ****/
272  xsize = X.nbCells(3);
273  if(xsize == xsize_old)
274  stability = true;
275  xsize_old = xsize;
276 
277  } while( !stability );
278 
279  if(verbose) trace.endBlock();
280 
281  return X;
282 }
283 
284 //////////////////////////////////////////////////////////////////////////////
285 // Select Functions
286 //////////////////////////////////////////////////////////////////////////////
287 template < typename TComplex >
288 std::pair<typename TComplex::Cell, typename TComplex::Data>
289 DGtal::functions::selectFirst(
290  const typename TComplex::Clique & clique)
291 {
292  return *clique.begin(3);
293 }
294 
295 template < typename TComplex >
296 std::pair<typename TComplex::Cell, typename TComplex::Data>
297 DGtal::functions::selectRandom(
298  const typename TComplex::Clique & clique)
299 {
300  static std::mt19937 gen( std::random_device{}() );
301  return selectRandom<TComplex>(clique, gen);
302 }
303 
304 template < typename TComplex, typename TRandomGenerator >
305 std::pair<typename TComplex::Cell, typename TComplex::Data>
306 DGtal::functions::selectRandom(
307  const typename TComplex::Clique & clique,
308  TRandomGenerator & gen
309 )
310 {
311  auto size = clique.nbCells(3);
312  std::uniform_int_distribution<> dis(0, size - 1);
313  auto it = clique.begin(3);
314  std::advance(it, dis(gen));
315  return *it;
316 }
317 
318 template < typename TDistanceTransform, typename TComplex >
319 std::pair<typename TComplex::Cell, typename TComplex::Data>
320 DGtal::functions::selectMaxValue(
321  const TDistanceTransform & dist_map,
322  const typename TComplex::Clique & clique)
323 {
324  typename TDistanceTransform::Value max_value{0};
325  typename TDistanceTransform::Value value{0};
326  using ComplexConstIterator = typename TComplex::CellMapConstIterator;
327  ComplexConstIterator selected_pair;
328 
329  for(ComplexConstIterator it = clique.begin(3), itE = clique.end(3);
330  it != itE ; ++it){
331  value = dist_map(clique.space().uCoords(it->first));
332  // trace.info() << "P: " << it->first << " V: " << value << std::endl;
333  if(value > max_value){
334  max_value = value;
335  selected_pair = it;
336  continue;
337  }
338  if(value == max_value){
339  // max_value = value;
340  // selected_pair = it;
341  continue; // TODO choose one wisely when they have same DM value.
342  }
343  }
344 
345  return *selected_pair;
346 }
347 //////////////////////////////////////////////////////////////////////////////
348 // Skeleton Functions
349 //////////////////////////////////////////////////////////////////////////////
350 
351 template < typename TComplex >
352 bool
353 DGtal::functions::skelUltimate( const TComplex &,
354  const typename TComplex::Cell &)
355 {
356  return false;
357 }
358 /*--------------------------------------------------------------------------*/
359 
360 template < typename TComplex >
361 bool
362 DGtal::functions::skelEnd( const TComplex & vc,
363  const typename TComplex::Cell & cell)
364 {
365  const auto pnsize = vc.properNeighborhoodVoxels(cell).size();
366  return (pnsize == 1);
367 }
368 
369 template < typename TComplex >
370 bool
371 DGtal::functions::skelSimple( const TComplex & vc,
372  const typename TComplex::Cell & cell)
373 {
374  return vc.isSimple(cell);
375 }
376 
377 template < typename TComplex >
378 bool
379 DGtal::functions::skelIsthmus(const TComplex & vc,
380  const typename TComplex::Cell & cell)
381 {
382  if(twoIsthmus<TComplex>(vc,cell))
383  return true;
384  else if(oneIsthmus<TComplex>(vc,cell))
385  return true;
386  else
387  return false;
388 }
389 
390 template < typename TComplex >
391 bool
392 DGtal::functions::oneIsthmus(const TComplex & vc,
393  const typename TComplex::Cell & cell)
394 {
395 
396  auto &ks = vc.space();
397  auto point_cell = ks.uCoords(cell);
398  // Create an object from the set.
399  using DigitalTopology = DGtal::Z3i::DT26_6;
400  using DigitalSet = DGtal::DigitalSetByAssociativeContainer<
401  DGtal::Z3i::Domain,
402  std::unordered_set<typename DGtal::Z3i::Domain::Point>>;
403  using Object = DGtal::Object<DigitalTopology, DigitalSet>;
404  DigitalTopology::ForegroundAdjacency adjF;
405  DigitalTopology::BackgroundAdjacency adjB;
406  auto domain = DGtal::Z3i::Domain(ks.lowerBound(), ks.upperBound());
407  DigitalTopology topo(
408  adjF, adjB, DGtal::DigitalTopologyProperties::JORDAN_DT);
409  Object object(topo, domain);
410  for(auto & properNeighborhoodVoxel : vc.properNeighborhoodVoxels(cell))
411  {
412  object.pointSet().insertNew(ks.uCoords(properNeighborhoodVoxel));
413  }
414 
415  // auto spN = vc.object().properNeighborhood(point_cell);
416  if (isZeroSurface(object)) return true;
417 
418  // else thinning
419  TComplex pre_thin(ks);
420  pre_thin.construct(object.pointSet());
421 
422  auto after_thin = asymetricThinningScheme<TComplex>( pre_thin,
423  selectFirst<TComplex>,
424  skelUltimate<TComplex>);
425 
426  object.pointSet().clear();
427  for (auto it = after_thin.begin(3), itE = after_thin.end(3) ; it != itE ; ++it ){
428  object.pointSet().insertNew(ks.uCoords(it->first));
429  }
430 
431  return isZeroSurface(object);
432 }
433 
434 template < typename TComplex >
435 bool
436 DGtal::functions::twoIsthmus( const TComplex & vc,
437  const typename TComplex::Cell & cell)
438 {
439  auto &ks = vc.space();
440  auto point_cell = ks.uCoords(cell);
441  // Create an object from the set.
442  using DigitalTopology = DGtal::Z3i::DT26_6;
443  using DigitalSet = DGtal::DigitalSetByAssociativeContainer<
444  DGtal::Z3i::Domain,
445  std::unordered_set<typename DGtal::Z3i::Domain::Point>>;
446  using Object = DGtal::Object<DigitalTopology, DigitalSet>;
447  DigitalTopology::ForegroundAdjacency adjF;
448  DigitalTopology::BackgroundAdjacency adjB;
449  auto domain = DGtal::Z3i::Domain(ks.lowerBound(), ks.upperBound());
450  DigitalTopology topo(
451  adjF, adjB, DGtal::DigitalTopologyProperties::JORDAN_DT);
452  Object object(topo, domain);
453  for(auto & properNeighborhoodVoxel : vc.properNeighborhoodVoxels(cell))
454  {
455  object.pointSet().insertNew(ks.uCoords(properNeighborhoodVoxel));
456  }
457  if (isOneSurface(object)) return true;
458 
459  // else thinning
460  TComplex pre_thin(ks);
461  pre_thin.construct(object.pointSet());
462 
463  auto after_thin = asymetricThinningScheme<TComplex>( pre_thin,
464  selectFirst<TComplex>,
465  skelUltimate<TComplex>);
466 
467  object.pointSet().clear();
468  for (auto it = after_thin.begin(3), itE = after_thin.end(3) ; it != itE ; ++it ){
469  object.pointSet().insertNew(ks.uCoords(it->first));
470  }
471 
472  return isOneSurface(object);
473 }
474 
475 template < typename TComplex >
476 bool
477 DGtal::functions::skelWithTable(
478  const boost::dynamic_bitset<> & table,
479  const std::unordered_map<typename TComplex::Point, unsigned int> & pointToMaskMap,
480  const TComplex & vc,
481  const typename TComplex::Cell & cell)
482 {
483  auto conf = functions::getSpelNeighborhoodConfigurationOccupancy(
484  vc,
485  vc.space().uCoords(cell),
486  pointToMaskMap);
487  return table[conf];
488 }
489 ///////////////////////////////////////////////////////////////////////////////
490 // Object Helpers
491 template < typename TObject >
492 bool
493 DGtal::functions::isZeroSurface(const TObject & small_obj)
494 {
495  if (small_obj.size() != 2) return false;
496  auto connectedness = small_obj.computeConnectedness();
497  if( connectedness == DGtal::Connectedness::CONNECTED ) return false;
498 
499  return true;
500 }
501 
502 template < typename TObject >
503 bool
504 DGtal::functions::isOneSurface(const TObject & small_obj)
505 {
506  auto connectedness = small_obj.computeConnectedness();
507  if( connectedness == DGtal::Connectedness::DISCONNECTED )
508  return false;
509 
510  for (const auto & p : small_obj.pointSet()){
511  auto pN = small_obj.properNeighborhood(p);
512  if(!isZeroSurface(pN) ) return false;
513  }
514 
515  return true;
516 }
517 
518 template <typename TObject >
519 std::vector< TObject >
520 DGtal::functions::
521 connectedComponents(const TObject & input_obj, bool verbose)
522 {
523 
524  // Graph related alias
525  using Graph = TObject ;
526  using vertex_descriptor =
527  typename boost::graph_traits<Graph>::vertex_descriptor ; // Object::Vertex
528  using edge_descriptor =
529  typename boost::graph_traits<Graph>::edge_descriptor ; // Object::Edge
530  using vertices_size_type =
531  typename boost::graph_traits<Graph>::vertices_size_type ; // Object::Size
532 
533  using StdColorMap = std::map< vertex_descriptor, boost::default_color_type > ;
534  StdColorMap colorMap;
535  boost::associative_property_map< StdColorMap > propColorMap( colorMap );
536 
537  using StdComponentMap = std::map< vertex_descriptor, vertices_size_type > ;
538  StdComponentMap componentMap;
539  boost::associative_property_map< StdComponentMap > propComponentMap( componentMap );
540 
541  if(verbose) trace.beginBlock( "Boost connected_components");
542  vertices_size_type nbComp =
543  boost::connected_components // boost graph connected components algorithm.
544  ( input_obj, // the graph
545  propComponentMap, // the mapping vertex -> label
546  boost::color_map( propColorMap ) // this map is used internally when computing connected components.
547  );
548  if(verbose) trace.info() << "num_components = " << nbComp << std::endl;
549 
550  if(verbose) trace.beginBlock( "Filter graph and get components");
551 
552  using ComponentGraph =
553  boost::filtered_graph<
554  Graph,
555  std::function<bool(edge_descriptor)>,
556  std::function<bool(vertex_descriptor)>
557  >;
558  auto &g = input_obj;
559 
560  std::vector<ComponentGraph> component_graphs;
561 
562  for (size_t i = 0; i < nbComp; i++)
563  component_graphs.emplace_back(g,
564  [componentMap,i,&g](edge_descriptor e) {
565  return componentMap.at(boost::source(e,g) )==i
566  || componentMap.at(boost::target(e,g))==i;
567  },
568  [componentMap,i](vertex_descriptor v) {
569  return componentMap.at(v)==i;
570  });
571  if(verbose) trace.endBlock();
572 
573  std::vector<TObject> obj_components;
574  for (auto && c : component_graphs){
575  // Create empty but valid object.
576  obj_components.emplace_back(
577  TObject( input_obj.topology(),input_obj.domain() ));
578  for (auto && vp = boost::vertices(c); vp.first != vp.second ; ++vp.first){
579  obj_components.back().pointSet().insertNew(*vp.first);
580  }
581  }
582  if(verbose) trace.endBlock();
583 
584  return obj_components;
585 }
586 ///////////////////////////////////////////////////////////////////////////////
587 
588 ///////////////////////////////////////////////////////////////////////////////