DGtal  0.9.2
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/topology/helpers/NeighborhoodConfigurationsHelper.h>
33 //////////////////////////////////////////////////////////////////////////////
34 
35 ///////////////////////////////////////////////////////////////////////////////
36 // IMPLEMENTATION of inline methods.
37 ///////////////////////////////////////////////////////////////////////////////
38 
39 //-----------------------------------------------------------------------------
40 template <typename TKSpace, typename TCellContainer,
41  typename CellConstIterator,
42  typename CellMapIteratorPriority >
43 DGtal::uint64_t
44 DGtal::functions::
45 collapse( CubicalComplex< TKSpace, TCellContainer > & K,
46  CellConstIterator S_itB, CellConstIterator S_itE,
47  const CellMapIteratorPriority& priority,
48  bool hintIsSClosed, bool hintIsKClosed,
49  bool verbose )
50 {
51  using namespace std;
52  typedef CubicalComplex< TKSpace, TCellContainer > CC;
53  typedef typename CC::Cell Cell;
54  typedef typename CC::CellType CellType;
55  typedef typename CC::CellMapIterator CellMapIterator;
56  typedef vector< CellMapIterator > CMIVector;
57  typedef typename CMIVector::const_iterator CMIVectorConstIterator;
58  // NB : a maximal k-cell is collapsible if it has a free incident k-1-cell.
59  Dimension n = K.dim();
60  CMIVector S; // stores the cells to process
61  CMIVector Q_low; // stores the iterators on direct faces of the maximal cell.
62  CMIVector Q_collapsible;// stores collapsible cells in order to clean them at the end.
63  CellMapIterator it_cell; // generic iterator on a cell.
64  CellMapIterator it_cell_c; // points on c in the free pair (c,d)
65  CellMapIterator it_cell_d; // points on d in the free pair (c,d)
66  CMIVectorConstIterator itlow;
67  CMIVectorConstIterator itqup;
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  CellMapIterator best_free_face_it = K.end( 0 );
167  for ( CMIVectorConstIterator it = Q_low.begin(), itE = Q_low.end();
168  it != itE; ++it )
169  {
170  CellMapIterator low_ic = *it;
171  uint32_t& data = low_ic->second.data;
172  // trace.info() << " + Cell " << low_ic->first << " data=" << data << std::endl;
173  if ( ( data & CC::REMOVED ) || ! ( data & CC::COLLAPSIBLE ) ) continue;
174  const Cell& cur_d = low_ic->first;
175  CellType cur_d_type = K.computeCellType( cur_d, itup, n );
176  // trace.info() << " + Type=" << cur_d_type << std::endl;
177  if ( cur_d_type == CC::Free )
178  { // found a free n-1-face ic
179  if ( ( best_free_face_it == K.end( 0 ) )
180  || ( ! priority( low_ic, best_free_face_it ) ) )
181  best_free_face_it = low_ic;
182  }
183  }
184  if ( best_free_face_it != K.end( 0 ) )
185  {
186  // delete c and ic.
187  found_pair = true;
188  it_cell_c = itcur;
189  it_cell_d = best_free_face_it;
190  // Q_low already contains cells that should be
191  // checked again
192  }
193  }
194  else if ( cur_c_type == CC::Free )
195  { // free face... check that its 1-up-incident face is maximal.
196  CellMapIterator it_up_up;
197  const Cell& cur_d = itup->first;
198  CellType cur_d_type = K.computeCellType( cur_d, it_up_up, n );
199  if ( cur_d_type == CC::Maximal )
200  { // found a maximal face.
201  found_pair = true;
202  it_cell_c = itup;
203  it_cell_d = itcur;
204  // Q_low will contain cells that should be checked
205  // again
206  back_insert_iterator< CMIVector > back_it( Q_low );
207  K.directFacesIterators( back_it, it_cell_c->first );
208  }
209  }
210  if ( found_pair )
211  { // If found, remove pair from complex (logical removal).
212  it_cell_c->second.data |= CC::REMOVED;
213  it_cell_d->second.data |= CC::REMOVED;
214  nb_removed += 2;
215  // Incident cells have to be checked again.
216  for ( CMIVectorConstIterator it = Q_low.begin(), itE = Q_low.end();
217  it != itE; ++it )
218  {
219  it_cell = *it;
220  uint32_t& data_qlow = it_cell->second.data;
221  if ( ( ! ( data_qlow & CC::REMOVED ) )
222  && ( data_qlow & CC::COLLAPSIBLE )
223  && ( ! ( data_qlow & CC::USER1 ) ) )
224  {
225  S.push_back( it_cell );
226  }
227  }
228  }
229  Q_low.clear();
230  } // while ( ! PQ.empty() )
231  } // while ( ! S.empty() )
232 
233  if ( verbose ) trace.info() << "[CC::collapse]-+ cleaning complex." << std::endl;
234 
235  // Now clean the complex so that removed cells are effectively
236  // removed and no more cell is tagged as collapsible.
237  for ( CMIVectorConstIterator it = Q_collapsible.begin(), itE = Q_collapsible.end();
238  it != itE; ++it )
239  {
240  CellMapIterator cmIt = *it;
241  uint32_t& cur_data = cmIt->second.data;
242  if ( cur_data & CC::REMOVED ) K.eraseCell( cmIt );
243  else cur_data &= ~CC::COLLAPSIBLE;
244  // if ( (*it)->second.data & CC::REMOVED )
245  // K.eraseCell( *it );
246  // else
247  // (*it)->second.data &= ~CC::COLLAPSIBLE;
248  }
249  return nb_removed;
250 }
251 
252 
253 //-----------------------------------------------------------------------------
254 template <typename TKSpace, typename TCellContainer,
255  typename BdryCellOutputIterator,
256  typename InnerCellOutputIterator>
257 void
258 DGtal::functions::
259 filterCellsWithinBounds( const CubicalComplex< TKSpace, TCellContainer > & K,
260  const typename TKSpace::Point& kLow,
261  const typename TKSpace::Point& kUp,
262  BdryCellOutputIterator itBdry,
263  InnerCellOutputIterator itInner )
264 {
265  typedef CubicalComplex< TKSpace, TCellContainer > CC;
266  typedef typename CC::Cell Cell;
267  typedef typename CC::Point Point;
268  typedef typename CC::CellMapConstIterator CellMapConstIterator;
269  Dimension d = K.dim();
270  for ( Dimension i = 0; i <= d; ++i )
271  {
272  for ( CellMapConstIterator it = K.begin( i ), itE = K.end( i ); it != itE; ++it )
273  {
274  Cell cell = it->first;
275  Point kCell = K.space().uKCoords( cell );
276  if ( ( kLow.inf( kCell ) == kLow ) && ( kUp.sup( kCell ) == kUp ) )
277  { // Inside or on boundary.
278  bool bdry = false;
279  for ( Dimension j = 0; j < Point::dimension; ++j )
280  {
281  if ( ( kCell[ j ] == kLow[ j ] ) || ( kCell[ j ] == kUp[ j ] ) )
282  {
283  bdry = true;
284  break;
285  }
286  }
287  if ( bdry ) *itBdry++ = cell;
288  else *itInner++ = cell;
289  }
290  }
291  }
292 }
293 
294  template<typename TComplex>
295  DGtal::NeighborhoodConfiguration
296  DGtal::functions::
297  getSpelNeighborhoodConfigurationOccupancy(
298  const TComplex & input_complex,
299  const typename TComplex::Point & center,
300  const std::unordered_map<
301  typename TComplex::Point, NeighborhoodConfiguration> & mapPointToMask)
302  {
303  using Point = typename TComplex::Point;
304  using Space = SpaceND< Point::dimension > ;
305  using Domain = HyperRectDomain< Space >;
306  using DomainConstIterator = typename Domain::ConstIterator;
307 
308  Point p1 = Point::diagonal( -1 );
309  Point p2 = Point::diagonal( 1 );
310  Point c = Point::diagonal( 0 );
311  Domain domain( p1, p2 );
312  const auto & not_found( input_complex.end() );
313  // const auto & not_found( input_complex.end(3) );
314  NeighborhoodConfiguration cfg{0};
315  for ( DomainConstIterator it = domain.begin(); it != domain.end(); ++it ) {
316  if( *it != c &&
317  input_complex.find(
318  input_complex.space().uSpel(center + *it) ) != not_found )
319  // input_complex.findCell( 3,
320  // input_complex.space().uSpel(center + *it) ) != not_found )
321  cfg |= mapPointToMask.at(*it) ;
322  }
323  return cfg;
324  }
325 // //
326 ///////////////////////////////////////////////////////////////////////////////
327 
328