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