DGtal  1.5.beta
CubicalComplexFunctions.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 CubicalComplexFunctions.ih
19  * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20  * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
21  *
22  * @date 2015/11/24
23  *
24  * Implementation of inline methods defined in CubicalComplexFunctions.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 #include "DGtal/kernel/domains/HyperRectDomain.h"
33 #include "DGtal/topology/DigitalTopology.h"
34 #include "DGtal/topology/helpers/NeighborhoodConfigurationsHelper.h"
35 //////////////////////////////////////////////////////////////////////////////
36 
37 ///////////////////////////////////////////////////////////////////////////////
38 // IMPLEMENTATION of inline methods.
39 ///////////////////////////////////////////////////////////////////////////////
40 
41 //-----------------------------------------------------------------------------
42 template <typename TKSpace, typename TCellContainer,
43  typename CellConstIterator,
44  typename CellMapIteratorPriority >
45 DGtal::uint64_t
46 DGtal::functions::
47 collapse( CubicalComplex< TKSpace, TCellContainer > & K,
48  CellConstIterator S_itB, CellConstIterator S_itE,
49  const CellMapIteratorPriority& priority,
50  bool hintIsSClosed, bool hintIsKClosed,
51  bool verbose )
52 {
53  using namespace std;
54  typedef CubicalComplex< TKSpace, TCellContainer > CC;
55  typedef typename CC::Cell Cell;
56  typedef typename CC::CellType CellType;
57  typedef typename CC::CellMapIterator CellMapIterator;
58  typedef vector< CellMapIterator > CMIVector;
59  typedef typename CMIVector::const_iterator CMIVectorConstIterator;
60  // NB : a maximal k-cell is collapsible if it has a free incident k-1-cell.
61  Dimension n = K.dim();
62  CMIVector S; // stores the cells to process
63  CMIVector Q_low; // stores the iterators on direct faces of the maximal cell.
64  CMIVector Q_collapsible;// stores collapsible cells in order to clean them at the end.
65  CellMapIterator it_cell; // generic iterator on a cell.
66  CellMapIterator it_cell_c; // points on c in the free pair (c,d)
67  CellMapIterator it_cell_d; // points on d in the free pair (c,d)
68 
69  if ( verbose ) trace.info() << "[CC::collapse]-+ tag collapsible elements... " << flush;
70  // Restricts the set of elements that are collapsible.
71  if ( hintIsSClosed )
72  for ( CellConstIterator S_it = S_itB; S_it != S_itE; ++S_it )
73  {
74  Cell c = *S_it;
75  Dimension k = K.dim( c );
76  it_cell = K.findCell( k, c );
77  uint32_t& ccdata = it_cell->second.data;
78  ASSERT( it_cell != K.end( k ) );
79  S.push_back( it_cell );
80  if ( ! ( ccdata & (CC::FIXED | CC::COLLAPSIBLE ) ) )
81  {
82  ccdata |= CC::COLLAPSIBLE;
83  Q_collapsible.push_back( it_cell );
84  }
85  }
86  else // not ( hintIsSClosed )
87  for ( CellConstIterator S_it = S_itB; S_it != S_itE; ++S_it )
88  {
89  Cell c = *S_it;
90  Dimension k = K.dim( c );
91  it_cell = K.findCell( k, c );
92  uint32_t& ccdata = it_cell->second.data;
93  ASSERT( it_cell != K.end( k ) );
94  S.push_back( it_cell );
95  if ( ! ( ccdata & (CC::FIXED | CC::COLLAPSIBLE ) ) )
96  {
97  ccdata |= CC::COLLAPSIBLE;
98  Q_collapsible.push_back( it_cell );
99  }
100  vector<Cell> cells;
101  back_insert_iterator< vector<Cell> > back_it( cells );
102  K.faces( back_it, c, hintIsKClosed );
103  for ( typename vector<Cell>::const_iterator
104  it = cells.begin(), itE = cells.end(); it != itE; ++it )
105  {
106  it_cell = K.findCell( *it );
107  uint32_t& ccdata2 = it_cell->second.data;
108  if ( ! ( ccdata2 & (CC::FIXED | CC::COLLAPSIBLE ) ) )
109  {
110  ccdata2 |= CC::COLLAPSIBLE;
111  Q_collapsible.push_back( it_cell );
112  }
113  }
114  }
115  if ( verbose ) trace.info() << " " << Q_collapsible.size() << " found." << endl;
116 
117  // Fill queue
118  priority_queue<CellMapIterator, CMIVector, CellMapIteratorPriority> PQ( priority );
119 
120  if ( verbose ) trace.info() << "[CC::collapse]-+ entering collapsing loop. " << endl;
121  uint64_t nb_pass = 0;
122  uint64_t nb_examined = 0;
123  uint64_t nb_removed = 0;
124  while ( ! S.empty() )
125  {
126  for ( CMIVectorConstIterator it = S.begin(), itE = S.end();
127  it != itE; ++it )
128  {
129  PQ.push( *it );
130  (*it)->second.data |= CC::USER1;
131  }
132  S.clear();
133  if ( verbose ) trace.info() << "[CC::collapse]---+ Pass " << ++nb_pass
134  << ", Card(PQ)=" << PQ.size() << " elements, "
135  << "nb_exam=" << nb_examined << endl;
136 
137  // Try to collapse elements according to priority queue.
138  while ( ! PQ.empty() )
139  {
140  // Get top element.
141  CellMapIterator itcur = PQ.top();
142  uint32_t& cur_data = itcur->second.data;
143  PQ.pop();
144  ++nb_examined;
145  // Check if the cell is removable
146  if ( ( cur_data & CC::REMOVED ) || ( ! ( cur_data & CC::COLLAPSIBLE ) ) )
147  continue;
148  // Check if the cell was several time in the queue and is already processed.
149  if ( ! ( cur_data & CC::USER1 ) )
150  continue;
151  ASSERT( cur_data & CC::USER1 );
152  cur_data &= ~CC::USER1;
153 
154  // Cell may be removable.
155  // Check if it is a maximal cell
156  CellMapIterator itup;
157  const Cell & cur_c = itcur->first;
158  CellType cur_c_type = K.computeCellType( cur_c, itup, n );
159  bool found_pair = false;
160  // trace.info() << " - Cell " << cur_c << " Dim=" << dim( cur_c ) << " Type=" << cur_c_type << std::endl;
161  if ( cur_c_type == CC::Maximal )
162  { // maximal cell... must find a free face
163  // check faces to find a free face.
164  back_insert_iterator< CMIVector > back_it( Q_low );
165  K.directFacesIterators( back_it, cur_c );
166  bool best_free_face_found = false;
167  CellMapIterator best_free_face_it;
168  for ( CMIVectorConstIterator it = Q_low.begin(), itE = Q_low.end();
169  it != itE; ++it )
170  {
171  CellMapIterator low_ic = *it;
172  uint32_t& data = low_ic->second.data;
173  // trace.info() << " + Cell " << low_ic->first << " data=" << data << std::endl;
174  if ( ( data & CC::REMOVED ) || ! ( data & CC::COLLAPSIBLE ) ) continue;
175  const Cell& cur_d = low_ic->first;
176  CellType cur_d_type = K.computeCellType( cur_d, itup, n );
177  // trace.info() << " + Type=" << cur_d_type << std::endl;
178  if ( cur_d_type == CC::Free )
179  { // found a free n-1-face ic
180  if ( ( ! best_free_face_found )
181  || ( ! priority( low_ic, best_free_face_it ) ) )
182  {
183  best_free_face_it = low_ic;
184  best_free_face_found = true;
185  }
186  }
187  }
188  if ( best_free_face_found )
189  {
190  // delete c and ic.
191  found_pair = true;
192  it_cell_c = itcur;
193  it_cell_d = best_free_face_it;
194  // Q_low already contains cells that should be
195  // checked again
196  }
197  }
198  else if ( cur_c_type == CC::Free )
199  { // free face... check that its 1-up-incident face is maximal.
200  CellMapIterator it_up_up;
201  const Cell& cur_d = itup->first;
202  CellType cur_d_type = K.computeCellType( cur_d, it_up_up, n );
203  if ( cur_d_type == CC::Maximal )
204  { // found a maximal face.
205  found_pair = true;
206  it_cell_c = itup;
207  it_cell_d = itcur;
208  // Q_low will contain cells that should be checked
209  // again
210  back_insert_iterator< CMIVector > back_it( Q_low );
211  K.directFacesIterators( back_it, it_cell_c->first );
212  }
213  }
214  if ( found_pair )
215  { // If found, remove pair from complex (logical removal).
216  it_cell_c->second.data |= CC::REMOVED;
217  it_cell_d->second.data |= CC::REMOVED;
218  nb_removed += 2;
219  // Incident cells have to be checked again.
220  for ( CMIVectorConstIterator it = Q_low.begin(), itE = Q_low.end();
221  it != itE; ++it )
222  {
223  it_cell = *it;
224  uint32_t& data_qlow = it_cell->second.data;
225  if ( ( ! ( data_qlow & CC::REMOVED ) )
226  && ( data_qlow & CC::COLLAPSIBLE )
227  && ( ! ( data_qlow & CC::USER1 ) ) )
228  {
229  S.push_back( it_cell );
230  }
231  }
232  }
233  Q_low.clear();
234  } // while ( ! PQ.empty() )
235  } // while ( ! S.empty() )
236 
237  if ( verbose ) trace.info() << "[CC::collapse]-+ cleaning complex." << std::endl;
238 
239  // Now clean the complex so that removed cells are effectively
240  // removed and no more cell is tagged as collapsible.
241  for ( CMIVectorConstIterator it = Q_collapsible.begin(), itE = Q_collapsible.end();
242  it != itE; ++it )
243  {
244  CellMapIterator cmIt = *it;
245  uint32_t& cur_data = cmIt->second.data;
246  if ( cur_data & CC::REMOVED ) K.eraseCell( cmIt );
247  else cur_data &= ~CC::COLLAPSIBLE;
248  // if ( (*it)->second.data & CC::REMOVED )
249  // K.eraseCell( *it );
250  // else
251  // (*it)->second.data &= ~CC::COLLAPSIBLE;
252  }
253  return nb_removed;
254 }
255 
256 
257 //-----------------------------------------------------------------------------
258 template <typename TKSpace, typename TCellContainer,
259  typename BdryCellOutputIterator,
260  typename InnerCellOutputIterator>
261 void
262 DGtal::functions::
263 filterCellsWithinBounds( const CubicalComplex< TKSpace, TCellContainer > & K,
264  const typename TKSpace::Point& kLow,
265  const typename TKSpace::Point& kUp,
266  BdryCellOutputIterator itBdry,
267  InnerCellOutputIterator itInner )
268 {
269  typedef CubicalComplex< TKSpace, TCellContainer > CC;
270  typedef typename CC::Cell Cell;
271  typedef typename CC::Point Point;
272  typedef typename CC::CellMapConstIterator CellMapConstIterator;
273  Dimension d = K.dim();
274  for ( Dimension i = 0; i <= d; ++i )
275  {
276  for ( CellMapConstIterator it = K.begin( i ), itE = K.end( i ); it != itE; ++it )
277  {
278  Cell cell = it->first;
279  Point kCell = K.space().uKCoords( cell );
280  if ( ( kLow.inf( kCell ) == kLow ) && ( kUp.sup( kCell ) == kUp ) )
281  { // Inside or on boundary.
282  bool bdry = false;
283  for ( Dimension j = 0; j < Point::dimension; ++j )
284  {
285  if ( ( kCell[ j ] == kLow[ j ] ) || ( kCell[ j ] == kUp[ j ] ) )
286  {
287  bdry = true;
288  break;
289  }
290  }
291  if ( bdry ) *itBdry++ = cell;
292  else *itInner++ = cell;
293  }
294  }
295  }
296 }
297 
298 //-----------------------------------------------------------------------------
299 template <typename TObject, typename TKSpace, typename TCellContainer>
300 std::unique_ptr<TObject>
301 DGtal::functions::
302 objectFromSpels(
303  const CubicalComplex< TKSpace, TCellContainer > & C)
304 {
305 
306  const auto & cs = C.space();
307  // Get domain of C KSpace.
308  typename TObject::Domain domain(cs.lowerBound(), cs.upperBound());
309  typename TObject::DigitalSet dset(domain);
310  for(auto it = C.begin(TKSpace::DIM);
311  it!= C.end(TKSpace::DIM) ; ++it)
312  dset.insertNew(cs.uCoords(it->first));
313 
314  typename TObject::ForegroundAdjacency fAdj;
315  typename TObject::BackgroundAdjacency bAdj;
316  typename TObject::DigitalTopology topo(fAdj,bAdj, JORDAN_DT);
317  std::unique_ptr<TObject> obj( new TObject(topo, dset) );
318  return obj;
319 }
320 
321 template<typename TComplex>
322 DGtal::NeighborhoodConfiguration
323 DGtal::functions::
324 getSpelNeighborhoodConfigurationOccupancy
325 ( const TComplex & input_complex,
326  const typename TComplex::Point & center,
327  const std::unordered_map<
328  typename TComplex::Point, NeighborhoodConfiguration> & mapPointToMask )
329 {
330  using Point = typename TComplex::Point;
331  using Space = typename TComplex::Space;
332  using Domain = HyperRectDomain< Space >;
333 
334  Point p1 = Point::diagonal( -1 );
335  Point p2 = Point::diagonal( 1 );
336  Point c = Point::diagonal( 0 );
337  Domain domain( p1, p2 );
338 
339  using KPreSpace = typename TComplex::KSpace::PreCellularGridSpace;
340  const auto & kspace = input_complex.space();
341  const auto & not_found( input_complex.end(3) );
342  NeighborhoodConfiguration cfg{0};
343  for ( auto it = domain.begin(); it != domain.end(); ++it ) {
344  if( *it == c) continue;
345  const auto pre_cell = KPreSpace::uSpel(center + *it);
346  if (input_complex.belongs(pre_cell) &&
347  input_complex.findCell(3, kspace.uCell(pre_cell)) != not_found )
348  cfg |= mapPointToMask.at(*it) ;
349  }
350  return cfg;
351 }
352 // //
353 ///////////////////////////////////////////////////////////////////////////////
354 
355