DGtal 1.4.0
Loading...
Searching...
No Matches
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//-----------------------------------------------------------------------------
45template <typename TKSpace>
46DGtal::CellGeometry<TKSpace>::
47CellGeometry()
48 : myK(), myKPoints(),
49 myMinCellDim( 0 ), myMaxCellDim( KSpace::dimension ),
50 myVerbose( false )
51
52{
53}
54//-----------------------------------------------------------------------------
55template <typename TKSpace>
56DGtal::CellGeometry<TKSpace>::
57CellGeometry( 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//-----------------------------------------------------------------------------
67template <typename TKSpace>
68void
69DGtal::CellGeometry<TKSpace>::
70init( 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//-----------------------------------------------------------------------------
83template <typename TKSpace>
84void
85DGtal::CellGeometry<TKSpace>::
86addCellsTouchingPoint( const Point& p )
87{
88 addCellsTouchingPointel( myK.uPointel( p ) );
89}
90
91//-----------------------------------------------------------------------------
92template <typename TKSpace>
93void
94DGtal::CellGeometry<TKSpace>::
95addCellsTouchingPointel( 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//-----------------------------------------------------------------------------
118template <typename TKSpace>
119void
120DGtal::CellGeometry<TKSpace>::
121addCellsTouchingCell( 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//-----------------------------------------------------------------------------
137template <typename TKSpace>
138void
139DGtal::CellGeometry<TKSpace>::
140addCellsTouchingSegment( 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//-----------------------------------------------------------------------------
172template <typename TKSpace>
173template <typename PointIterator>
174void
175DGtal::CellGeometry<TKSpace>::
176addCellsTouchingPoints( 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//-----------------------------------------------------------------------------
208template <typename TKSpace>
209template <typename PointelIterator>
210void
211DGtal::CellGeometry<TKSpace>::
212addCellsTouchingPointels( 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//-----------------------------------------------------------------------------
243template <typename TKSpace>
244void
245DGtal::CellGeometry<TKSpace>::
246addCellsTouchingPolytopePoints( const LatticePolytope& polytope )
247{
248 std::vector< Point > points;
249 polytope.getPoints( points );
250 addCellsTouchingPoints( points.begin(), points.end() );
251}
252//-----------------------------------------------------------------------------
253template <typename TKSpace>
254void
255DGtal::CellGeometry<TKSpace>::
256addCellsTouchingPolytopePoints( const RationalPolytope& polytope )
257{
258 std::vector< Point > points;
259 polytope.getPoints( points );
260 addCellsTouchingPoints( points.cbegin(), points.cend() );
261}
262//-----------------------------------------------------------------------------
263template <typename TKSpace>
264void
265DGtal::CellGeometry<TKSpace>::
266addCellsTouchingPolytope( 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//-----------------------------------------------------------------------------
275template <typename TKSpace>
276void
277DGtal::CellGeometry<TKSpace>::
278addCellsTouchingPolytope( 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//-----------------------------------------------------------------------------
287template <typename TKSpace>
288typename DGtal::CellGeometry<TKSpace>&
289DGtal::CellGeometry<TKSpace>::
290operator+=( 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//-----------------------------------------------------------------------------
302template <typename TKSpace>
303typename DGtal::CellGeometry<TKSpace>::Size
304DGtal::CellGeometry<TKSpace>::
305nbCells() const
306{
307 return (DGtal::CellGeometry<TKSpace>::Size)myKPoints.size();
308}
309//-----------------------------------------------------------------------------
310template <typename TKSpace>
311typename DGtal::CellGeometry<TKSpace>::Size
312DGtal::CellGeometry<TKSpace>::
313computeNbCells( 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//-----------------------------------------------------------------------------
323template <typename TKSpace>
324std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
325DGtal::CellGeometry<TKSpace>::
326getKPoints( 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//-----------------------------------------------------------------------------
336template <typename TKSpace>
337typename DGtal::CellGeometry<TKSpace>::Integer
338DGtal::CellGeometry<TKSpace>::
339computeEuler() 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//-----------------------------------------------------------------------------
352template <typename TKSpace>
353DGtal::Dimension
354DGtal::CellGeometry<TKSpace>::
355minCellDim() const
356{
357 return myMinCellDim;
358}
359//-----------------------------------------------------------------------------
360template <typename TKSpace>
361DGtal::Dimension
362DGtal::CellGeometry<TKSpace>::
363maxCellDim() const
364{
365 return myMaxCellDim;
366}
367
368//-----------------------------------------------------------------------------
369template <typename TKSpace>
370bool
371DGtal::CellGeometry<TKSpace>::
372subset( const CellGeometry& other ) const
373{
374 return other.myKPoints.includes( myKPoints );
375}
376//-----------------------------------------------------------------------------
377template <typename TKSpace>
378bool
379DGtal::CellGeometry<TKSpace>::
380subset( 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//-----------------------------------------------------------------------------
390template <typename TKSpace>
391template <typename RandomIterator>
392bool
393DGtal::CellGeometry<TKSpace>::
394includes( 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//-----------------------------------------------------------------------------
424template <typename TKSpace>
425std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
426DGtal::CellGeometry<TKSpace>::
427getIntersectedKPoints( 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//-----------------------------------------------------------------------------
488template <typename TKSpace>
489std::vector< typename DGtal::CellGeometry<TKSpace>::Cell >
490DGtal::CellGeometry<TKSpace>::
491getIntersectedCells( 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//-----------------------------------------------------------------------------
551template <typename TKSpace>
552std::vector< typename DGtal::CellGeometry<TKSpace>::Cell >
553DGtal::CellGeometry<TKSpace>::
554getIntersectedCells( 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//-----------------------------------------------------------------------------
614template <typename TKSpace>
615std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
616DGtal::CellGeometry<TKSpace>::
617getIntersectedKPoints( 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//-----------------------------------------------------------------------------
678template <typename TKSpace>
679std::vector< typename DGtal::CellGeometry<TKSpace>::Cell >
680DGtal::CellGeometry<TKSpace>::
681getTouchedCells( 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//-----------------------------------------------------------------------------
708template <typename TKSpace>
709std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
710DGtal::CellGeometry<TKSpace>::
711getTouchedKPoints( 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//-----------------------------------------------------------------------------
739template <typename TKSpace>
740DGtal::Dimension
741DGtal::CellGeometry<TKSpace>::
742dim( 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 */
757template <typename TKSpace>
758inline
759void
760DGtal::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 */
779template <typename TKSpace>
780inline
781bool
782DGtal::CellGeometry<TKSpace>::isValid() const
783{
784 return true;
785}
786//-----------------------------------------------------------------------------
787template <typename TKSpace>
788inline
789std::string
790DGtal::CellGeometry<TKSpace>::className
791() const
792{
793 return "CellGeometry";
794}
795
796
797
798///////////////////////////////////////////////////////////////////////////////
799// Implementation of inline functions //
800
801//-----------------------------------------------------------------------------
802template <typename TKSpace>
803inline
804std::ostream&
805DGtal::operator<< ( std::ostream & out,
806 const CellGeometry<TKSpace> & object )
807{
808 object.selfDisplay( out );
809 return out;
810}
811
812// //
813///////////////////////////////////////////////////////////////////////////////