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