DGtal 1.3.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( 0 <= myMinCellDim );
64 ASSERT( myMinCellDim <= myMaxCellDim );
65 ASSERT( myMaxCellDim <= myK.dimension );
66}
67//-----------------------------------------------------------------------------
68template <typename TKSpace>
69void
70DGtal::CellGeometry<TKSpace>::
71init( 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//-----------------------------------------------------------------------------
85template <typename TKSpace>
86void
87DGtal::CellGeometry<TKSpace>::
88addCellsTouchingPoint( const Point& p )
89{
90 addCellsTouchingPointel( myK.uPointel( p ) );
91}
92
93//-----------------------------------------------------------------------------
94template <typename TKSpace>
95void
96DGtal::CellGeometry<TKSpace>::
97addCellsTouchingPointel( 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//-----------------------------------------------------------------------------
120template <typename TKSpace>
121void
122DGtal::CellGeometry<TKSpace>::
123addCellsTouchingCell( 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//-----------------------------------------------------------------------------
139template <typename TKSpace>
140void
141DGtal::CellGeometry<TKSpace>::
142addCellsTouchingSegment( 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//-----------------------------------------------------------------------------
174template <typename TKSpace>
175template <typename PointIterator>
176void
177DGtal::CellGeometry<TKSpace>::
178addCellsTouchingPoints( 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//-----------------------------------------------------------------------------
210template <typename TKSpace>
211template <typename PointelIterator>
212void
213DGtal::CellGeometry<TKSpace>::
214addCellsTouchingPointels( 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//-----------------------------------------------------------------------------
245template <typename TKSpace>
246void
247DGtal::CellGeometry<TKSpace>::
248addCellsTouchingPolytopePoints( const LatticePolytope& polytope )
249{
250 std::vector< Point > points;
251 polytope.getPoints( points );
252 addCellsTouchingPoints( points.begin(), points.end() );
253}
254//-----------------------------------------------------------------------------
255template <typename TKSpace>
256void
257DGtal::CellGeometry<TKSpace>::
258addCellsTouchingPolytopePoints( const RationalPolytope& polytope )
259{
260 std::vector< Point > points;
261 polytope.getPoints( points );
262 addCellsTouchingPoints( points.cbegin(), points.cend() );
263}
264//-----------------------------------------------------------------------------
265template <typename TKSpace>
266void
267DGtal::CellGeometry<TKSpace>::
268addCellsTouchingPolytope( 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//-----------------------------------------------------------------------------
277template <typename TKSpace>
278void
279DGtal::CellGeometry<TKSpace>::
280addCellsTouchingPolytope( 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//-----------------------------------------------------------------------------
289template <typename TKSpace>
290typename DGtal::CellGeometry<TKSpace>&
291DGtal::CellGeometry<TKSpace>::
292operator+=( 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//-----------------------------------------------------------------------------
304template <typename TKSpace>
305typename DGtal::CellGeometry<TKSpace>::Size
306DGtal::CellGeometry<TKSpace>::
307nbCells() const
308{
309 return myKPoints.size();
310}
311//-----------------------------------------------------------------------------
312template <typename TKSpace>
313typename DGtal::CellGeometry<TKSpace>::Size
314DGtal::CellGeometry<TKSpace>::
315computeNbCells( 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//-----------------------------------------------------------------------------
325template <typename TKSpace>
326std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
327DGtal::CellGeometry<TKSpace>::
328getKPoints( 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//-----------------------------------------------------------------------------
338template <typename TKSpace>
339typename DGtal::CellGeometry<TKSpace>::Integer
340DGtal::CellGeometry<TKSpace>::
341computeEuler() 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//-----------------------------------------------------------------------------
354template <typename TKSpace>
355DGtal::Dimension
356DGtal::CellGeometry<TKSpace>::
357minCellDim() const
358{
359 return myMinCellDim;
360}
361//-----------------------------------------------------------------------------
362template <typename TKSpace>
363DGtal::Dimension
364DGtal::CellGeometry<TKSpace>::
365maxCellDim() const
366{
367 return myMaxCellDim;
368}
369
370//-----------------------------------------------------------------------------
371template <typename TKSpace>
372bool
373DGtal::CellGeometry<TKSpace>::
374subset( const CellGeometry& other ) const
375{
376 return other.myKPoints.includes( myKPoints );
377}
378//-----------------------------------------------------------------------------
379template <typename TKSpace>
380bool
381DGtal::CellGeometry<TKSpace>::
382subset( 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//-----------------------------------------------------------------------------
392template <typename TKSpace>
393template <typename RandomIterator>
394bool
395DGtal::CellGeometry<TKSpace>::
396includes( 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//-----------------------------------------------------------------------------
426template <typename TKSpace>
427std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
428DGtal::CellGeometry<TKSpace>::
429getIntersectedKPoints( 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//-----------------------------------------------------------------------------
490template <typename TKSpace>
491std::vector< typename DGtal::CellGeometry<TKSpace>::Cell >
492DGtal::CellGeometry<TKSpace>::
493getIntersectedCells( 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//-----------------------------------------------------------------------------
553template <typename TKSpace>
554std::vector< typename DGtal::CellGeometry<TKSpace>::Cell >
555DGtal::CellGeometry<TKSpace>::
556getIntersectedCells( 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//-----------------------------------------------------------------------------
616template <typename TKSpace>
617std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
618DGtal::CellGeometry<TKSpace>::
619getIntersectedKPoints( 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//-----------------------------------------------------------------------------
680template <typename TKSpace>
681std::vector< typename DGtal::CellGeometry<TKSpace>::Cell >
682DGtal::CellGeometry<TKSpace>::
683getTouchedCells( 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//-----------------------------------------------------------------------------
710template <typename TKSpace>
711std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
712DGtal::CellGeometry<TKSpace>::
713getTouchedKPoints( 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//-----------------------------------------------------------------------------
741template <typename TKSpace>
742DGtal::Dimension
743DGtal::CellGeometry<TKSpace>::
744dim( 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 */
759template <typename TKSpace>
760inline
761void
762DGtal::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 */
781template <typename TKSpace>
782inline
783bool
784DGtal::CellGeometry<TKSpace>::isValid() const
785{
786 return true;
787}
788//-----------------------------------------------------------------------------
789template <typename TKSpace>
790inline
791std::string
792DGtal::CellGeometry<TKSpace>::className
793() const
794{
795 return "CellGeometry";
796}
797
798
799
800///////////////////////////////////////////////////////////////////////////////
801// Implementation of inline functions //
802
803//-----------------------------------------------------------------------------
804template <typename TKSpace>
805inline
806std::ostream&
807DGtal::operator<< ( std::ostream & out,
808 const CellGeometry<TKSpace> & object )
809{
810 object.selfDisplay( out );
811 return out;
812}
813
814// //
815///////////////////////////////////////////////////////////////////////////////
Generated on Sun Nov 27 2022 15:15:41 for DGtal by doxygen 1.9.5