DGtal  1.5.beta
CellGeometry.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 CellGeometry.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 2020/01/02
23  *
24  * Implementation of inline methods defined in CellGeometry.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 #include <string>
33 #include <sstream>
34 #include <algorithm>
35 //////////////////////////////////////////////////////////////////////////////
36 
37 ///////////////////////////////////////////////////////////////////////////////
38 // IMPLEMENTATION of inline methods.
39 ///////////////////////////////////////////////////////////////////////////////
40 
41 ///////////////////////////////////////////////////////////////////////////////
42 // ----------------------- Standard services ------------------------------
43 
44 //-----------------------------------------------------------------------------
45 template <typename TKSpace>
46 DGtal::CellGeometry<TKSpace>::
47 CellGeometry()
48  : myK(), myKPoints(),
49  myMinCellDim( 0 ), myMaxCellDim( KSpace::dimension ),
50  myVerbose( false )
51 
52 {
53 }
54 //-----------------------------------------------------------------------------
55 template <typename TKSpace>
56 DGtal::CellGeometry<TKSpace>::
57 CellGeometry( const KSpace & K,
58  Dimension min_cell_dim, Dimension max_cell_dim, bool verbose )
59  : myK( K ), myKPoints(),
60  myMinCellDim( min_cell_dim ), myMaxCellDim( max_cell_dim ),
61  myVerbose( verbose )
62 {
63  ASSERT( myMinCellDim <= myMaxCellDim );
64  ASSERT( myMaxCellDim <= myK.dimension );
65 }
66 //-----------------------------------------------------------------------------
67 template <typename TKSpace>
68 void
69 DGtal::CellGeometry<TKSpace>::
70 init( const KSpace & K,
71  Dimension min_cell_dim, Dimension max_cell_dim, bool verbose )
72 {
73  ASSERT( myMinCellDim <= myMaxCellDim );
74  ASSERT( myMaxCellDim <= myK.dimension );
75  myK = K;
76  myKPoints.clear();
77  myMinCellDim = min_cell_dim;
78  myMaxCellDim = max_cell_dim;
79  myVerbose = verbose;
80 }
81 
82 //-----------------------------------------------------------------------------
83 template <typename TKSpace>
84 void
85 DGtal::CellGeometry<TKSpace>::
86 addCellsTouchingPoint( const Point& p )
87 {
88  addCellsTouchingPointel( myK.uPointel( p ) );
89 }
90 
91 //-----------------------------------------------------------------------------
92 template <typename TKSpace>
93 void
94 DGtal::CellGeometry<TKSpace>::
95 addCellsTouchingPointel( const Cell& pointel )
96 {
97  auto cofaces = myK.uCoFaces( pointel );
98  if ( ( myMinCellDim == 0 ) && ( myMaxCellDim == KSpace::dimension ) )
99  {
100  myKPoints.emplace( myK.uKCoords( pointel ) );
101  for ( auto && c : cofaces )
102  myKPoints.emplace( myK.uKCoords( c ) );
103  }
104  else
105  {
106  if ( myMinCellDim <= 0 )
107  myKPoints.emplace( myK.uKCoords( pointel ) );
108  for ( auto&& f : cofaces )
109  {
110  Dimension d = myK.uDim( f );
111  if ( ( myMinCellDim <= d ) && ( d <= myMaxCellDim ) )
112  myKPoints.emplace( myK.uKCoords( f ) );
113  }
114  }
115 }
116 
117 //-----------------------------------------------------------------------------
118 template <typename TKSpace>
119 void
120 DGtal::CellGeometry<TKSpace>::
121 addCellsTouchingCell( const Cell& c )
122 {
123  const Dimension dc = myK.uDim( c );
124  if ( myMinCellDim <= dc && dc <= myMaxCellDim )
125  myKPoints.emplace( myK.uKCoords( c ) );
126  if ( myMaxCellDim <= dc ) return;
127  const auto cofaces = myK.uCoFaces( c );
128  for ( auto&& f : cofaces )
129  {
130  Dimension d = myK.uDim( f );
131  if ( ( myMinCellDim <= d ) && ( d <= myMaxCellDim ) )
132  myKPoints.emplace( myK.uKCoords( f ) );
133  }
134 }
135 
136 //-----------------------------------------------------------------------------
137 template <typename TKSpace>
138 void
139 DGtal::CellGeometry<TKSpace>::
140 addCellsTouchingSegment( const Point& a, const Point& b )
141 {
142  const auto V = b - a;
143  addCellsTouchingPoint( a );
144  for ( Dimension k = 0; k < dimension; k++ )
145  {
146  const Integer n = ( V[ k ] >= 0 ) ? V[ k ] : -V[ k ];
147  const Integer d = ( V[ k ] >= 0 ) ? 1 : -1;
148  if ( n == 0 ) continue;
149  Point kc;
150  for ( Integer i = 1; i < n; i++ )
151  {
152  for ( Dimension j = 0; j < dimension; j++ )
153  {
154  if ( j == k ) kc[ k ] = 2 * ( a[ k ] + d * i );
155  else
156  {
157  const auto v = V[ j ];
158  const auto q = ( v * i ) / n;
159  const auto r = ( v * i ) % n; // might be negative
160  kc[ j ] = 2 * ( a[ j ] + q );
161  if ( r < 0 ) kc[ j ] -= 1;
162  else if ( r > 0 ) kc[ j ] += 1;
163  }
164  }
165  addCellsTouchingCell( myK.uCell( kc ) );
166  }
167  }
168  if ( a != b ) addCellsTouchingPoint( b );
169 }
170 
171 //-----------------------------------------------------------------------------
172 template <typename TKSpace>
173 template <typename PointIterator>
174 void
175 DGtal::CellGeometry<TKSpace>::
176 addCellsTouchingPoints( PointIterator itB, PointIterator itE )
177 {
178  if ( ( myMinCellDim == 0 ) && ( myMaxCellDim == KSpace::dimension ) )
179  {
180  for ( auto it = itB; it != itE; ++it )
181  {
182  auto pointel = myK.uPointel( *it );
183  auto cofaces = myK.uCoFaces( pointel );
184  myKPoints.emplace( myK.uKCoords( pointel ) );
185  for ( auto && c : cofaces )
186  myKPoints.emplace( myK.uKCoords( c ) );
187  }
188  }
189  else
190  {
191  for ( auto it = itB; it != itE; ++it )
192  {
193  auto pointel = myK.uPointel( *it );
194  auto cofaces = myK.uCoFaces( pointel );
195  if ( myMinCellDim <= 0 )
196  myKPoints.emplace( myK.uKCoords( pointel ) );
197  for ( auto&& f : cofaces )
198  {
199  Dimension d = myK.uDim( f );
200  if ( ( myMinCellDim <= d ) && ( d <= myMaxCellDim ) )
201  myKPoints.emplace( myK.uKCoords( f ) );
202  }
203  }
204  }
205 }
206 
207 //-----------------------------------------------------------------------------
208 template <typename TKSpace>
209 template <typename PointelIterator>
210 void
211 DGtal::CellGeometry<TKSpace>::
212 addCellsTouchingPointels( PointelIterator itB, PointelIterator itE )
213 {
214  if ( ( myMinCellDim == 0 ) && ( myMaxCellDim == KSpace::dimension ) )
215  {
216  for ( auto it = itB; it != itE; ++it )
217  {
218  auto pointel = *it;
219  auto cofaces = myK.uCoFaces( pointel );
220  myKPoints.emplace( myK.uKCoords( pointel ) );
221  for ( auto && c : cofaces )
222  myKPoints.emplace( myK.uKCoords( c ) );
223  }
224  }
225  else
226  {
227  for ( auto it = itB; it != itE; ++it )
228  {
229  auto pointel = *it;
230  auto cofaces = myK.uCoFaces( pointel );
231  if ( myMinCellDim <= 0 )
232  myKPoints.emplace( myK.uKCoords( pointel ) );
233  for ( auto&& f : cofaces )
234  {
235  Dimension d = myK.uDim( f );
236  if ( ( myMinCellDim <= d ) && ( d <= myMaxCellDim ) )
237  myKPoints.emplace( myK.uKCoords( f ) );
238  }
239  }
240  }
241 }
242 //-----------------------------------------------------------------------------
243 template <typename TKSpace>
244 void
245 DGtal::CellGeometry<TKSpace>::
246 addCellsTouchingPolytopePoints( const LatticePolytope& polytope )
247 {
248  std::vector< Point > points;
249  polytope.getPoints( points );
250  addCellsTouchingPoints( points.begin(), points.end() );
251 }
252 //-----------------------------------------------------------------------------
253 template <typename TKSpace>
254 void
255 DGtal::CellGeometry<TKSpace>::
256 addCellsTouchingPolytopePoints( const RationalPolytope& polytope )
257 {
258  std::vector< Point > points;
259  polytope.getPoints( points );
260  addCellsTouchingPoints( points.cbegin(), points.cend() );
261 }
262 //-----------------------------------------------------------------------------
263 template <typename TKSpace>
264 void
265 DGtal::CellGeometry<TKSpace>::
266 addCellsTouchingPolytope( const LatticePolytope& polytope )
267 {
268  for ( Dimension i = myMinCellDim; i <= myMaxCellDim; ++i )
269  {
270  auto kpoints = getIntersectedKPoints( polytope, i );
271  myKPoints.insert( kpoints.cbegin(), kpoints.cend() );
272  }
273 }
274 //-----------------------------------------------------------------------------
275 template <typename TKSpace>
276 void
277 DGtal::CellGeometry<TKSpace>::
278 addCellsTouchingPolytope( const RationalPolytope& polytope )
279 {
280  for ( Dimension i = myMinCellDim; i <= myMaxCellDim; ++i )
281  {
282  auto kpoints = getIntersectedKPoints( polytope, i );
283  myKPoints.insert( kpoints.cbegin(), kpoints.cend() );
284  }
285 }
286 //-----------------------------------------------------------------------------
287 template <typename TKSpace>
288 typename DGtal::CellGeometry<TKSpace>&
289 DGtal::CellGeometry<TKSpace>::
290 operator+=( const CellGeometry& other )
291 {
292  if ( this != &other )
293  {
294  myKPoints.insert( other.myKPoints.cbegin(), other.myKPoints.cend() );
295  myMinCellDim = std::min( myMinCellDim, other.myMinCellDim );
296  myMaxCellDim = std::max( myMaxCellDim, other.myMaxCellDim );
297  }
298  return *this;
299 }
300 
301 //-----------------------------------------------------------------------------
302 template <typename TKSpace>
303 typename DGtal::CellGeometry<TKSpace>::Size
304 DGtal::CellGeometry<TKSpace>::
305 nbCells() const
306 {
307  return (DGtal::CellGeometry<TKSpace>::Size)myKPoints.size();
308 }
309 //-----------------------------------------------------------------------------
310 template <typename TKSpace>
311 typename DGtal::CellGeometry<TKSpace>::Size
312 DGtal::CellGeometry<TKSpace>::
313 computeNbCells( const Dimension k ) const
314 {
315  if ( k < minCellDim() || k > maxCellDim() ) return 0;
316  Size nb = 0;
317  for ( auto&& c : myKPoints )
318  nb += ( dim( c ) == k ) ? 1 : 0;
319  return nb;
320 }
321 
322 //-----------------------------------------------------------------------------
323 template <typename TKSpace>
324 std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
325 DGtal::CellGeometry<TKSpace>::
326 getKPoints( const Dimension k ) const
327 {
328  std::vector< Point > R;
329  if ( k < minCellDim() || k > maxCellDim() ) return R;
330  for ( auto&& c : myKPoints )
331  if ( dim( c ) == k ) R.push_back( c );
332  return R;
333 }
334 
335 //-----------------------------------------------------------------------------
336 template <typename TKSpace>
337 typename DGtal::CellGeometry<TKSpace>::Integer
338 DGtal::CellGeometry<TKSpace>::
339 computeEuler() const
340 {
341  Integer euler = 0;
342  bool pos = true;
343  for ( Dimension k = 0; k <= maxCellDim(); ++k )
344  {
345  if ( pos ) euler += computeNbCells( k );
346  else euler -= computeNbCells( k );
347  pos = ! pos;
348  }
349  return euler;
350 }
351 //-----------------------------------------------------------------------------
352 template <typename TKSpace>
353 DGtal::Dimension
354 DGtal::CellGeometry<TKSpace>::
355 minCellDim() const
356 {
357  return myMinCellDim;
358 }
359 //-----------------------------------------------------------------------------
360 template <typename TKSpace>
361 DGtal::Dimension
362 DGtal::CellGeometry<TKSpace>::
363 maxCellDim() const
364 {
365  return myMaxCellDim;
366 }
367 
368 //-----------------------------------------------------------------------------
369 template <typename TKSpace>
370 bool
371 DGtal::CellGeometry<TKSpace>::
372 subset( const CellGeometry& other ) const
373 {
374  return other.myKPoints.includes( myKPoints );
375 }
376 //-----------------------------------------------------------------------------
377 template <typename TKSpace>
378 bool
379 DGtal::CellGeometry<TKSpace>::
380 subset( const CellGeometry& other, const Dimension k ) const
381 {
382  UnorderedSetByBlock< Point, Splitter< Point, uint64_t > > k_dim_points;
383  for ( auto&& c : myKPoints )
384  if ( dim( c ) == k )
385  k_dim_points.insert( c );
386  return other.myKPoints.includes( k_dim_points );
387 }
388 
389 //-----------------------------------------------------------------------------
390 template <typename TKSpace>
391 template <typename RandomIterator>
392 bool
393 DGtal::CellGeometry<TKSpace>::
394 includes( RandomIterator it1, RandomIterator itE1,
395  RandomIterator it2, RandomIterator itE2 )
396 {
397  std::size_t k;
398  for ( ; it2 != itE2; ++it1)
399  {
400  if (it1 == itE1 || *it2 < *it1) return false;
401  // exponential march
402  for ( k = 1; ( it1 < itE1 ) && ( *it1 < *it2 ); k *= 2 ) it1 += k;
403  if ( it1 < itE1 )
404  {
405  if ( *it2 == *it1 ) ++it2; //equality
406  else
407  {
408  it1 = lower_bound( it1 - k/2, it1, *it2 );
409  if ( *it2 != *it1 ) return false;
410  ++it2;
411  }
412  }
413  else
414  {
415  it1 = lower_bound( it1 - k/2, itE1, *it2 );
416  if ( it1 == itE1 || *it2 != *it1 ) return false;
417  ++it2;
418  }
419  }
420  return true;
421 }
422 
423 //-----------------------------------------------------------------------------
424 template <typename TKSpace>
425 std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
426 DGtal::CellGeometry<TKSpace>::
427 getIntersectedKPoints( const LatticePolytope& polytope,
428  const Dimension i ) const
429 {
430  ASSERT( polytope.canBeSummed() );
431  if ( ! polytope.canBeSummed() )
432  trace.warning() << "[CellGeometryFunctions::getIntersectedKPoints]"
433  << " LatticePolytope is not valid for Minkowski sums. "
434  << std::endl;
435  static const Dimension d = KSpace::dimension;
436  std::vector< Point > result;
437  std::vector< Point > points;
438  std::vector< LatticePolytope > polytopes( i+1, polytope );
439  std::vector< Dimension > extensions( i+1, 0 );
440  for ( Dimension k = 1; k < extensions.size(); ++k )
441  {
442  extensions[ k ] = k - 1;
443  polytopes [ k ] = polytopes[ k - 1 ]
444  + typename LatticePolytope::UnitSegment( k - 1 );
445  }
446  // We have to build several dilated polytopes which corresponds
447  // to the binom(d,i) possible cell types.
448  while ( true )
449  {
450  if ( myVerbose )
451  {
452  std::string str;
453  std::ostringstream ostr( str );
454  ostr << "Dilating Polytope along directions {";
455  for ( Dimension k = 1; k < extensions.size(); ++k )
456  ostr << " + " << extensions[ k ];
457  ostr << "}" ;
458  trace.info() << ostr.str() << std::endl;
459  }
460  // Intersected cells are bijective to points in a dilated polytope.
461  polytopes.back().getPoints( points );
462  // For each point, build its Khalimsky points and push it into result.
463  for ( auto p : points )
464  {
465  auto kp = myK.uKCoords( myK.uPointel( p ) );
466  for ( Dimension k = 1; k < extensions.size(); ++k )
467  // decrease Khalimsky coordinate to get incident cell
468  kp[ extensions[ k ] ] -= 1;
469  result.push_back( kp );
470  }
471  // Go to next type of cell
472  Dimension k = i;
473  extensions[ k ] += 1;
474  // will quit when k == 0
475  while ( k > 0 && extensions[ k ] >= d+k-i ) extensions[ --k ] += 1;
476  if ( k == 0 ) break; // finished
477  for ( Dimension l = k + 1; l < extensions.size(); ++l )
478  extensions[ l ] = extensions[ l - 1 ] + 1;
479  // Recomputes polytopes
480  for ( ; k < extensions.size(); ++k )
481  polytopes [ k ] = polytopes[ k - 1 ]
482  + typename LatticePolytope::UnitSegment( extensions[ k ] );
483  } // while ( true )
484  return result;
485 }
486 
487 //-----------------------------------------------------------------------------
488 template <typename TKSpace>
489 std::vector< typename DGtal::CellGeometry<TKSpace>::Cell >
490 DGtal::CellGeometry<TKSpace>::
491 getIntersectedCells( const LatticePolytope& polytope,
492  const Dimension i ) const
493 {
494  ASSERT( polytope.canBeSummed() );
495  if ( ! polytope.canBeSummed() )
496  trace.warning() << "[CellGeometryFunctions::getIntersectedCells]"
497  << " LatticePolytope is not valid for Minkowski sums. "
498  << std::endl;
499  static const Dimension d = KSpace::dimension;
500  std::vector< Cell > result;
501  std::vector< Point > points;
502  std::vector< LatticePolytope > polytopes( i+1, polytope );
503  std::vector< Dimension > extensions( i+1, 0 );
504  for ( Dimension k = 1; k < extensions.size(); ++k )
505  {
506  extensions[ k ] = k - 1;
507  polytopes [ k ] = polytopes[ k - 1 ]
508  + typename LatticePolytope::UnitSegment( k - 1 );
509  }
510  // We have to build several dilated polytopes which corresponds
511  // to the binom(d,i) possible cell types.
512  while ( true )
513  {
514  if ( myVerbose )
515  {
516  std::string str;
517  std::ostringstream ostr( str );
518  ostr << "Dilating Polytope along directions {";
519  for ( Dimension k = 1; k < extensions.size(); ++k )
520  ostr << " + " << extensions[ k ];
521  ostr << "}" ;
522  trace.info() << ostr.str() << std::endl;
523  }
524  // Intersected cells are bijective to points in a dilated polytope.
525  polytopes.back().getPoints( points );
526  // For each point, build its cell and push it into result.
527  for ( auto p : points )
528  {
529  auto cell = myK.uPointel( p );
530  for ( Dimension k = 1; k < extensions.size(); ++k )
531  cell = myK.uIncident( cell, extensions[ k ], false );
532  result.push_back( cell );
533  }
534  // Go to next type of cell
535  Dimension k = i;
536  extensions[ k ] += 1;
537  // will quit when k == 0
538  while ( k > 0 && extensions[ k ] >= d+k-i ) extensions[ --k ] += 1;
539  if ( k == 0 ) break; // finished
540  for ( Dimension l = k + 1; l < extensions.size(); ++l )
541  extensions[ l ] = extensions[ l - 1 ] + 1;
542  // Recomputes polytopes
543  for ( ; k < extensions.size(); ++k )
544  polytopes [ k ] = polytopes[ k - 1 ]
545  + typename LatticePolytope::UnitSegment( extensions[ k ] );
546  } // while ( true )
547  return result;
548 }
549 
550 //-----------------------------------------------------------------------------
551 template <typename TKSpace>
552 std::vector< typename DGtal::CellGeometry<TKSpace>::Cell >
553 DGtal::CellGeometry<TKSpace>::
554 getIntersectedCells( const RationalPolytope& polytope,
555  const Dimension i ) const
556 {
557  ASSERT( polytope.canBeSummed() );
558  if ( ! polytope.canBeSummed() )
559  trace.warning() << "[CellGeometryFunctions::getIntersectedCells]"
560  << " RationalPolytope is not valid for Minkowski sums. "
561  << std::endl;
562  static const Dimension d = KSpace::dimension;
563  std::vector< Cell > result;
564  std::vector< Point > points;
565  std::vector< RationalPolytope > polytopes( i+1, polytope );
566  std::vector< Dimension > extensions( i+1, 0 );
567  for ( Dimension k = 1; k < extensions.size(); ++k )
568  {
569  extensions[ k ] = k - 1;
570  polytopes [ k ] = polytopes[ k - 1 ]
571  + typename RationalPolytope::UnitSegment( k - 1 );
572  }
573  // We have to build several dilated polytopes which corresponds
574  // to the binom(d,i) possible cell types.
575  while ( true )
576  {
577  if ( myVerbose )
578  {
579  std::string str;
580  std::ostringstream ostr( str );
581  ostr << "Dilating Polytope along directions {";
582  for ( Dimension k = 1; k < extensions.size(); ++k )
583  ostr << " + " << extensions[ k ];
584  ostr << "}" ;
585  trace.info() << ostr.str() << std::endl;
586  }
587  // Intersected cells are bijective to points in a dilated polytope.
588  polytopes.back().getPoints( points );
589  // For each point, build its cell and push it into result.
590  for ( auto p : points )
591  {
592  auto cell = myK.uPointel( p );
593  for ( Dimension k = 1; k < extensions.size(); ++k )
594  cell = myK.uIncident( cell, extensions[ k ], false );
595  result.push_back( cell );
596  }
597  // Go to next type of cell
598  Dimension k = i;
599  extensions[ k ] += 1;
600  // will quit when k == 0
601  while ( k > 0 && extensions[ k ] >= d+k-i ) extensions[ --k ] += 1;
602  if ( k == 0 ) break; // finished
603  for ( Dimension l = k + 1; l < extensions.size(); ++l )
604  extensions[ l ] = extensions[ l - 1 ] + 1;
605  // Recomputes polytopes
606  for ( ; k < extensions.size(); ++k )
607  polytopes [ k ] = polytopes[ k - 1 ]
608  + typename RationalPolytope::UnitSegment( extensions[ k ] );
609  } // while ( true )
610  return result;
611 }
612 
613 //-----------------------------------------------------------------------------
614 template <typename TKSpace>
615 std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
616 DGtal::CellGeometry<TKSpace>::
617 getIntersectedKPoints( const RationalPolytope& polytope,
618  const Dimension i ) const
619 {
620  ASSERT( polytope.canBeSummed() );
621  if ( ! polytope.canBeSummed() )
622  trace.warning() << "[CellGeometryFunctions::getIntersectedKPoints]"
623  << " RationalPolytope is not valid for Minkowski sums. "
624  << std::endl;
625  static const Dimension d = KSpace::dimension;
626  std::vector< Point > result;
627  std::vector< Point > points;
628  std::vector< RationalPolytope > polytopes( i+1, polytope );
629  std::vector< Dimension > extensions( i+1, 0 );
630  for ( Dimension k = 1; k < extensions.size(); ++k )
631  {
632  extensions[ k ] = k - 1;
633  polytopes [ k ] = polytopes[ k - 1 ]
634  + typename RationalPolytope::UnitSegment( k - 1 );
635  }
636  // We have to build several dilated polytopes which corresponds
637  // to the binom(d,i) possible cell types.
638  while ( true )
639  {
640  if ( myVerbose )
641  {
642  std::string str;
643  std::ostringstream ostr( str );
644  ostr << "Dilating Polytope along directions {";
645  for ( Dimension k = 1; k < extensions.size(); ++k )
646  ostr << " + " << extensions[ k ];
647  ostr << "}" ;
648  trace.info() << ostr.str() << std::endl;
649  }
650  // Intersected cells are bijective to points in a dilated polytope.
651  polytopes.back().getPoints( points );
652  // For each point, build its Khalimsky points and push it into result.
653  for ( auto p : points )
654  {
655  auto kp = myK.uKCoords( myK.uPointel( p ) );
656  for ( Dimension k = 1; k < extensions.size(); ++k )
657  // decrease Khalimsky coordinate to get incident cell
658  kp[ extensions[ k ] ] -= 1;
659  result.push_back( kp );
660  }
661  // Go to next type of cell
662  Dimension k = i;
663  extensions[ k ] += 1;
664  // will quit when k == 0
665  while ( k > 0 && extensions[ k ] >= d+k-i ) extensions[ --k ] += 1;
666  if ( k == 0 ) break; // finished
667  for ( Dimension l = k + 1; l < extensions.size(); ++l )
668  extensions[ l ] = extensions[ l - 1 ] + 1;
669  // Recomputes polytopes
670  for ( ; k < extensions.size(); ++k )
671  polytopes [ k ] = polytopes[ k - 1 ]
672  + typename RationalPolytope::UnitSegment( extensions[ k ] );
673  } // while ( true )
674  return result;
675 }
676 
677 //-----------------------------------------------------------------------------
678 template <typename TKSpace>
679 std::vector< typename DGtal::CellGeometry<TKSpace>::Cell >
680 DGtal::CellGeometry<TKSpace>::
681 getTouchedCells( const std::vector< Point >& points, const Dimension i ) const
682 {
683  std::unordered_set< Cell > cells;
684  if ( i == 0 )
685  cells = CellGeometryFunctions< KSpace, 0, KSpace::dimension>
686  ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
687  else if ( i == 1 )
688  cells = CellGeometryFunctions< KSpace, 1, KSpace::dimension>
689  ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
690  else if ( i == 2 )
691  cells = CellGeometryFunctions< KSpace, 2, KSpace::dimension>
692  ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
693  else if ( i == 3 )
694  cells = CellGeometryFunctions< KSpace, 3, KSpace::dimension>
695  ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
696  else if ( i == 4 )
697  cells = CellGeometryFunctions< KSpace, 4, KSpace::dimension>
698  ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
699  else if ( i == 5 )
700  cells = CellGeometryFunctions< KSpace, 5, KSpace::dimension>
701  ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
702  else trace.error() << "[DGtal::CellGeometry<TKSpace>::getTouchedCells]"
703  << " Computation are limited to n-D, n <= 5" << std::endl;
704  return std::vector< Cell >( cells.begin(), cells.end() );
705 }
706 
707 //-----------------------------------------------------------------------------
708 template <typename TKSpace>
709 std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
710 DGtal::CellGeometry<TKSpace>::
711 getTouchedKPoints( const std::vector< Point >& points, const Dimension i ) const
712 {
713  UnorderedSetByBlock< Point, Splitter< Point, uint64_t > > kpoints;
714  if ( i == 0 )
715  kpoints = CellGeometryFunctions< KSpace, 0, KSpace::dimension>
716  ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
717  else if ( i == 1 )
718  kpoints = CellGeometryFunctions< KSpace, 1, KSpace::dimension>
719  ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
720  else if ( i == 2 )
721  kpoints = CellGeometryFunctions< KSpace, 2, KSpace::dimension>
722  ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
723  else if ( i == 3 )
724  kpoints = CellGeometryFunctions< KSpace, 3, KSpace::dimension>
725  ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
726  else if ( i == 4 )
727  kpoints = CellGeometryFunctions< KSpace, 4, KSpace::dimension>
728  ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
729  else if ( i == 5 )
730  kpoints = CellGeometryFunctions< KSpace, 5, KSpace::dimension>
731  ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
732  else trace.error() << "[DGtal::CellGeometry<TKSpace>::getTouchedKPoints]"
733  << " Computation are limited to n-D, n <= 5" << std::endl;
734  return std::vector< Cell >( kpoints.begin(), kpoints.end() );
735 }
736 
737 
738 //-----------------------------------------------------------------------------
739 template <typename TKSpace>
740 DGtal::Dimension
741 DGtal::CellGeometry<TKSpace>::
742 dim( const Point & kp )
743 {
744  Dimension d = 0;
745  for ( Dimension i = 0; i < KSpace::dimension; ++i )
746  d += kp[ i ] & 1 ? 1 : 0;
747  return d;
748 }
749 
750 ///////////////////////////////////////////////////////////////////////////////
751 // Interface - public :
752 
753 /**
754  * Writes/Displays the object on an output stream.
755  * @param out the output stream where the object is written.
756  */
757 template <typename TKSpace>
758 inline
759 void
760 DGtal::CellGeometry<TKSpace>::selfDisplay ( std::ostream & out ) const
761 {
762  out << "[CellGeometry] ";
763  std::vector< Point > X;
764  for ( auto && c : myKPoints ) X.push_back( c );
765  std::sort( X.begin(), X.end() );
766  for ( auto p : X )
767  {
768  out << "(" << p[ 0 ];
769  for ( Dimension k = 1; k < TKSpace::dimension; ++k )
770  out << "," << p[ k ];
771  out << ") ";
772  }
773 }
774 
775 /**
776  * Checks the validity/consistency of the object.
777  * @return 'true' if the object is valid, 'false' otherwise.
778  */
779 template <typename TKSpace>
780 inline
781 bool
782 DGtal::CellGeometry<TKSpace>::isValid() const
783 {
784  return true;
785 }
786 //-----------------------------------------------------------------------------
787 template <typename TKSpace>
788 inline
789 std::string
790 DGtal::CellGeometry<TKSpace>::className
791 () const
792 {
793  return "CellGeometry";
794 }
795 
796 
797 
798 ///////////////////////////////////////////////////////////////////////////////
799 // Implementation of inline functions //
800 
801 //-----------------------------------------------------------------------------
802 template <typename TKSpace>
803 inline
804 std::ostream&
805 DGtal::operator<< ( std::ostream & out,
806  const CellGeometry<TKSpace> & object )
807 {
808  object.selfDisplay( out );
809  return out;
810 }
811 
812 // //
813 ///////////////////////////////////////////////////////////////////////////////