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