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