DGtal 2.1.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 if ( polytope.denominator() == 1 )
260 {
261 polytope.getPoints( points );
262 addCellsTouchingPoints( points.cbegin(), points.cend() );
263 }
264 else
265 {
266 auto Q = polytope.denominator() * polytope;
267 Q.getPoints( points );
268 addCellsTouchingPoints( points.cbegin(), points.cend() );
269 }
270}
271//-----------------------------------------------------------------------------
272template <typename TKSpace>
273void
274DGtal::CellGeometry<TKSpace>::
275addCellsTouchingPolytope( const LatticePolytope& polytope )
276{
277 for ( Dimension i = myMinCellDim; i <= myMaxCellDim; ++i )
278 {
279 auto kpoints = getIntersectedKPoints( polytope, i );
280 myKPoints.insert( kpoints.cbegin(), kpoints.cend() );
281 }
282}
283//-----------------------------------------------------------------------------
284template <typename TKSpace>
285void
286DGtal::CellGeometry<TKSpace>::
287addCellsTouchingPolytope( const RationalPolytope& polytope )
288{
289 if ( polytope.denominator() == 1 )
290 for ( Dimension i = myMinCellDim; i <= myMaxCellDim; ++i )
291 {
292 auto kpoints = getIntersectedKPoints( polytope, i );
293 myKPoints.insert( kpoints.cbegin(), kpoints.cend() );
294 }
295 else
296 {
297 auto Q = polytope.denominator() * polytope;
298 for ( Dimension i = myMinCellDim; i <= myMaxCellDim; ++i )
299 {
300 auto kpoints = getIntersectedKPoints( Q, i );
301 myKPoints.insert( kpoints.cbegin(), kpoints.cend() );
302 }
303 }
304}
305//-----------------------------------------------------------------------------
306template <typename TKSpace>
307typename DGtal::CellGeometry<TKSpace>&
308DGtal::CellGeometry<TKSpace>::
309operator+=( const CellGeometry& other )
310{
311 if ( this != &other )
312 {
313 myKPoints.insert( other.myKPoints.cbegin(), other.myKPoints.cend() );
314 myMinCellDim = std::min( myMinCellDim, other.myMinCellDim );
315 myMaxCellDim = std::max( myMaxCellDim, other.myMaxCellDim );
316 }
317 return *this;
318}
319
320//-----------------------------------------------------------------------------
321template <typename TKSpace>
322typename DGtal::CellGeometry<TKSpace>::Size
323DGtal::CellGeometry<TKSpace>::
324nbCells() const
325{
326 return (DGtal::CellGeometry<TKSpace>::Size)myKPoints.size();
327}
328//-----------------------------------------------------------------------------
329template <typename TKSpace>
330typename DGtal::CellGeometry<TKSpace>::Size
331DGtal::CellGeometry<TKSpace>::
332computeNbCells( const Dimension k ) const
333{
334 if ( k < minCellDim() || k > maxCellDim() ) return 0;
335 Size nb = 0;
336 for ( auto&& c : myKPoints )
337 nb += ( dim( c ) == k ) ? 1 : 0;
338 return nb;
339}
340
341//-----------------------------------------------------------------------------
342template <typename TKSpace>
343std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
344DGtal::CellGeometry<TKSpace>::
345getKPoints( const Dimension k ) const
346{
347 std::vector< Point > R;
348 if ( k < minCellDim() || k > maxCellDim() ) return R;
349 for ( auto&& c : myKPoints )
350 if ( dim( c ) == k ) R.push_back( c );
351 return R;
352}
353
354//-----------------------------------------------------------------------------
355template <typename TKSpace>
356typename DGtal::CellGeometry<TKSpace>::Integer
357DGtal::CellGeometry<TKSpace>::
358computeEuler() const
359{
360 Integer euler = 0;
361 bool pos = true;
362 for ( Dimension k = 0; k <= maxCellDim(); ++k )
363 {
364 if ( pos ) euler += computeNbCells( k );
365 else euler -= computeNbCells( k );
366 pos = ! pos;
367 }
368 return euler;
369}
370//-----------------------------------------------------------------------------
371template <typename TKSpace>
372DGtal::Dimension
373DGtal::CellGeometry<TKSpace>::
374minCellDim() const
375{
376 return myMinCellDim;
377}
378//-----------------------------------------------------------------------------
379template <typename TKSpace>
380DGtal::Dimension
381DGtal::CellGeometry<TKSpace>::
382maxCellDim() const
383{
384 return myMaxCellDim;
385}
386
387//-----------------------------------------------------------------------------
388template <typename TKSpace>
389bool
390DGtal::CellGeometry<TKSpace>::
391subset( const CellGeometry& other ) const
392{
393 return other.myKPoints.includes( myKPoints );
394}
395//-----------------------------------------------------------------------------
396template <typename TKSpace>
397bool
398DGtal::CellGeometry<TKSpace>::
399subset( const CellGeometry& other, const Dimension k ) const
400{
401 UnorderedSetByBlock< Point, Splitter< Point, uint64_t > > k_dim_points;
402 for ( auto&& c : myKPoints )
403 if ( dim( c ) == k )
404 k_dim_points.insert( c );
405 return other.myKPoints.includes( k_dim_points );
406}
407
408//-----------------------------------------------------------------------------
409template <typename TKSpace>
410template <typename RandomIterator>
411bool
412DGtal::CellGeometry<TKSpace>::
413includes( RandomIterator it1, RandomIterator itE1,
414 RandomIterator it2, RandomIterator itE2 )
415{
416 std::size_t k;
417 for ( ; it2 != itE2; ++it1)
418 {
419 if (it1 == itE1 || *it2 < *it1) return false;
420 // exponential march
421 for ( k = 1; ( it1 < itE1 ) && ( *it1 < *it2 ); k *= 2 ) it1 += k;
422 if ( it1 < itE1 )
423 {
424 if ( *it2 == *it1 ) ++it2; //equality
425 else
426 {
427 it1 = lower_bound( it1 - k/2, it1, *it2 );
428 if ( *it2 != *it1 ) return false;
429 ++it2;
430 }
431 }
432 else
433 {
434 it1 = lower_bound( it1 - k/2, itE1, *it2 );
435 if ( it1 == itE1 || *it2 != *it1 ) return false;
436 ++it2;
437 }
438 }
439 return true;
440}
441
442//-----------------------------------------------------------------------------
443template <typename TKSpace>
444std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
445DGtal::CellGeometry<TKSpace>::
446getIntersectedKPoints( const LatticePolytope& polytope,
447 const Dimension i ) const
448{
449 ASSERT( polytope.canBeSummed() );
450 if ( ! polytope.canBeSummed() )
451 trace.warning() << "[CellGeometryFunctions::getIntersectedKPoints]"
452 << " LatticePolytope is not valid for Minkowski sums. "
453 << std::endl;
454 static const Dimension d = KSpace::dimension;
455 std::vector< Point > result;
456 std::vector< Point > points;
457 std::vector< LatticePolytope > polytopes( i+1, polytope );
458 std::vector< Dimension > extensions( i+1, 0 );
459 for ( Dimension k = 1; k < extensions.size(); ++k )
460 {
461 extensions[ k ] = k - 1;
462 polytopes [ k ] = polytopes[ k - 1 ]
463 + typename LatticePolytope::UnitSegment( k - 1 );
464 }
465 // We have to build several dilated polytopes which corresponds
466 // to the binom(d,i) possible cell types.
467 while ( true )
468 {
469 if ( myVerbose )
470 {
471 std::string str;
472 std::ostringstream ostr( str );
473 ostr << "Dilating Polytope along directions {";
474 for ( Dimension k = 1; k < extensions.size(); ++k )
475 ostr << " + " << extensions[ k ];
476 ostr << "}" ;
477 trace.info() << ostr.str() << std::endl;
478 }
479 // Intersected cells are bijective to points in a dilated polytope.
480 polytopes.back().getPoints( points );
481 // For each point, build its Khalimsky points and push it into result.
482 for ( auto p : points )
483 {
484 auto kp = myK.uKCoords( myK.uPointel( p ) );
485 for ( Dimension k = 1; k < extensions.size(); ++k )
486 // decrease Khalimsky coordinate to get incident cell
487 kp[ extensions[ k ] ] -= 1;
488 result.push_back( kp );
489 }
490 // Go to next type of cell
491 Dimension k = i;
492 extensions[ k ] += 1;
493 // will quit when k == 0
494 while ( k > 0 && extensions[ k ] >= d+k-i ) extensions[ --k ] += 1;
495 if ( k == 0 ) break; // finished
496 for ( Dimension l = k + 1; l < extensions.size(); ++l )
497 extensions[ l ] = extensions[ l - 1 ] + 1;
498 // Recomputes polytopes
499 for ( ; k < extensions.size(); ++k )
500 polytopes [ k ] = polytopes[ k - 1 ]
501 + typename LatticePolytope::UnitSegment( extensions[ k ] );
502 } // while ( true )
503 return result;
504}
505
506//-----------------------------------------------------------------------------
507template <typename TKSpace>
508std::vector< typename DGtal::CellGeometry<TKSpace>::Cell >
509DGtal::CellGeometry<TKSpace>::
510getIntersectedCells( const LatticePolytope& polytope,
511 const Dimension i ) const
512{
513 ASSERT( polytope.canBeSummed() );
514 if ( ! polytope.canBeSummed() )
515 trace.warning() << "[CellGeometryFunctions::getIntersectedCells]"
516 << " LatticePolytope is not valid for Minkowski sums. "
517 << std::endl;
518 static const Dimension d = KSpace::dimension;
519 std::vector< Cell > result;
520 std::vector< Point > points;
521 std::vector< LatticePolytope > polytopes( i+1, polytope );
522 std::vector< Dimension > extensions( i+1, 0 );
523 for ( Dimension k = 1; k < extensions.size(); ++k )
524 {
525 extensions[ k ] = k - 1;
526 polytopes [ k ] = polytopes[ k - 1 ]
527 + typename LatticePolytope::UnitSegment( k - 1 );
528 }
529 // We have to build several dilated polytopes which corresponds
530 // to the binom(d,i) possible cell types.
531 while ( true )
532 {
533 if ( myVerbose )
534 {
535 std::string str;
536 std::ostringstream ostr( str );
537 ostr << "Dilating Polytope along directions {";
538 for ( Dimension k = 1; k < extensions.size(); ++k )
539 ostr << " + " << extensions[ k ];
540 ostr << "}" ;
541 trace.info() << ostr.str() << std::endl;
542 }
543 // Intersected cells are bijective to points in a dilated polytope.
544 polytopes.back().getPoints( points );
545 // For each point, build its cell and push it into result.
546 for ( auto p : points )
547 {
548 auto cell = myK.uPointel( p );
549 for ( Dimension k = 1; k < extensions.size(); ++k )
550 cell = myK.uIncident( cell, extensions[ k ], false );
551 result.push_back( cell );
552 }
553 // Go to next type of cell
554 Dimension k = i;
555 extensions[ k ] += 1;
556 // will quit when k == 0
557 while ( k > 0 && extensions[ k ] >= d+k-i ) extensions[ --k ] += 1;
558 if ( k == 0 ) break; // finished
559 for ( Dimension l = k + 1; l < extensions.size(); ++l )
560 extensions[ l ] = extensions[ l - 1 ] + 1;
561 // Recomputes polytopes
562 for ( ; k < extensions.size(); ++k )
563 polytopes [ k ] = polytopes[ k - 1 ]
564 + typename LatticePolytope::UnitSegment( extensions[ k ] );
565 } // while ( true )
566 return result;
567}
568
569//-----------------------------------------------------------------------------
570template <typename TKSpace>
571std::vector< typename DGtal::CellGeometry<TKSpace>::Cell >
572DGtal::CellGeometry<TKSpace>::
573getIntersectedCells( const RationalPolytope& polytope,
574 const Dimension i ) const
575{
576 ASSERT( polytope.canBeSummed() );
577 if ( ! polytope.canBeSummed() )
578 trace.warning() << "[CellGeometryFunctions::getIntersectedCells]"
579 << " RationalPolytope is not valid for Minkowski sums. "
580 << std::endl;
581 static const Dimension d = KSpace::dimension;
582 std::vector< Cell > result;
583 if ( polytope.denominator() != 1 )
584 {
585 trace.warning() << "[CellGeometryFunctions::getIntersectedCells]"
586 << "only valid for rational polytopes with denominator 1, i.e. lattice polytope."
587 << std::endl;
588 return result;
589 }
590
591 std::vector< Point > points;
592 std::vector< RationalPolytope > polytopes( i+1, polytope );
593 std::vector< Dimension > extensions( i+1, 0 );
594 for ( Dimension k = 1; k < extensions.size(); ++k )
595 {
596 extensions[ k ] = k - 1;
597 polytopes [ k ] = polytopes[ k - 1 ]
598 + typename RationalPolytope::UnitSegment( k - 1 );
599 }
600 // We have to build several dilated polytopes which corresponds
601 // to the binom(d,i) possible cell types.
602 while ( true )
603 {
604 if ( myVerbose )
605 {
606 std::string str;
607 std::ostringstream ostr( str );
608 ostr << "Dilating Polytope along directions {";
609 for ( Dimension k = 1; k < extensions.size(); ++k )
610 ostr << " + " << extensions[ k ];
611 ostr << "}" ;
612 trace.info() << ostr.str() << std::endl;
613 }
614 // Intersected cells are bijective to points in a dilated polytope.
615 polytopes.back().getPoints( points );
616 // For each point, build its cell and push it into result.
617 for ( auto p : points )
618 {
619 auto cell = myK.uPointel( p );
620 for ( Dimension k = 1; k < extensions.size(); ++k )
621 cell = myK.uIncident( cell, extensions[ k ], false );
622 result.push_back( cell );
623 }
624 // Go to next type of cell
625 Dimension k = i;
626 extensions[ k ] += 1;
627 // will quit when k == 0
628 while ( k > 0 && extensions[ k ] >= d+k-i ) extensions[ --k ] += 1;
629 if ( k == 0 ) break; // finished
630 for ( Dimension l = k + 1; l < extensions.size(); ++l )
631 extensions[ l ] = extensions[ l - 1 ] + 1;
632 // Recomputes polytopes
633 for ( ; k < extensions.size(); ++k )
634 polytopes [ k ] = polytopes[ k - 1 ]
635 + typename RationalPolytope::UnitSegment( extensions[ k ] );
636 } // while ( true )
637 return result;
638}
639
640//-----------------------------------------------------------------------------
641template <typename TKSpace>
642std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
643DGtal::CellGeometry<TKSpace>::
644getIntersectedKPoints( const RationalPolytope& polytope,
645 const Dimension i ) const
646{
647 ASSERT( polytope.canBeSummed() );
648 if ( ! polytope.canBeSummed() )
649 trace.warning() << "[CellGeometryFunctions::getIntersectedKPoints]"
650 << " RationalPolytope is not valid for Minkowski sums. "
651 << std::endl;
652 static const Dimension d = KSpace::dimension;
653 std::vector< Point > result;
654 if ( polytope.denominator() != 1 )
655 {
656 trace.warning() << "[CellGeometryFunctions::getIntersectedKPoints]"
657 << "only valid for rational polytopes with denominator 1, i.e. lattice polytope."
658 << std::endl;
659 return result;
660 }
661 std::vector< Point > points;
662 std::vector< RationalPolytope > polytopes( i+1, polytope );
663 std::vector< Dimension > extensions( i+1, 0 );
664 for ( Dimension k = 1; k < extensions.size(); ++k )
665 {
666 extensions[ k ] = k - 1;
667 polytopes [ k ] = polytopes[ k - 1 ]
668 + typename RationalPolytope::UnitSegment( k - 1 );
669 }
670 // We have to build several dilated polytopes which corresponds
671 // to the binom(d,i) possible cell types.
672 while ( true )
673 {
674 if ( myVerbose )
675 {
676 std::string str;
677 std::ostringstream ostr( str );
678 ostr << "Dilating Polytope along directions {";
679 for ( Dimension k = 1; k < extensions.size(); ++k )
680 ostr << " + " << extensions[ k ];
681 ostr << "}" ;
682 trace.info() << ostr.str() << std::endl;
683 }
684 // Intersected cells are bijective to points in a dilated polytope.
685 polytopes.back().getPoints( points );
686 // For each point, build its Khalimsky points and push it into result.
687 for ( auto p : points )
688 {
689 auto kp = myK.uKCoords( myK.uPointel( p ) );
690 for ( Dimension k = 1; k < extensions.size(); ++k )
691 // decrease Khalimsky coordinate to get incident cell
692 kp[ extensions[ k ] ] -= 1;
693 result.push_back( kp );
694 }
695 // Go to next type of cell
696 Dimension k = i;
697 extensions[ k ] += 1;
698 // will quit when k == 0
699 while ( k > 0 && extensions[ k ] >= d+k-i ) extensions[ --k ] += 1;
700 if ( k == 0 ) break; // finished
701 for ( Dimension l = k + 1; l < extensions.size(); ++l )
702 extensions[ l ] = extensions[ l - 1 ] + 1;
703 // Recomputes polytopes
704 for ( ; k < extensions.size(); ++k )
705 polytopes [ k ] = polytopes[ k - 1 ]
706 + typename RationalPolytope::UnitSegment( extensions[ k ] );
707 } // while ( true )
708 return result;
709}
710
711//-----------------------------------------------------------------------------
712template <typename TKSpace>
713std::vector< typename DGtal::CellGeometry<TKSpace>::Cell >
714DGtal::CellGeometry<TKSpace>::
715getTouchedCells( const std::vector< Point >& points, const Dimension i ) const
716{
717 std::unordered_set< Cell > cells;
718 if ( i == 0 )
719 cells = CellGeometryFunctions< KSpace, 0, KSpace::dimension>
720 ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
721 else if ( i == 1 )
722 cells = CellGeometryFunctions< KSpace, 1, KSpace::dimension>
723 ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
724 else if ( i == 2 )
725 cells = CellGeometryFunctions< KSpace, 2, KSpace::dimension>
726 ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
727 else if ( i == 3 )
728 cells = CellGeometryFunctions< KSpace, 3, KSpace::dimension>
729 ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
730 else if ( i == 4 )
731 cells = CellGeometryFunctions< KSpace, 4, KSpace::dimension>
732 ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
733 else if ( i == 5 )
734 cells = CellGeometryFunctions< KSpace, 5, KSpace::dimension>
735 ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
736 else trace.error() << "[DGtal::CellGeometry<TKSpace>::getTouchedCells]"
737 << " Computation are limited to n-D, n <= 5" << std::endl;
738 return std::vector< Cell >( cells.begin(), cells.end() );
739}
740
741//-----------------------------------------------------------------------------
742template <typename TKSpace>
743std::vector< typename DGtal::CellGeometry<TKSpace>::Point >
744DGtal::CellGeometry<TKSpace>::
745getTouchedKPoints( const std::vector< Point >& points, const Dimension i ) const
746{
747 UnorderedSetByBlock< Point, Splitter< Point, uint64_t > > kpoints;
748 if ( i == 0 )
749 kpoints = CellGeometryFunctions< KSpace, 0, KSpace::dimension>
750 ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
751 else if ( i == 1 )
752 kpoints = CellGeometryFunctions< KSpace, 1, KSpace::dimension>
753 ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
754 else if ( i == 2 )
755 kpoints = CellGeometryFunctions< KSpace, 2, KSpace::dimension>
756 ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
757 else if ( i == 3 )
758 kpoints = CellGeometryFunctions< KSpace, 3, KSpace::dimension>
759 ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
760 else if ( i == 4 )
761 kpoints = CellGeometryFunctions< KSpace, 4, KSpace::dimension>
762 ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
763 else if ( i == 5 )
764 kpoints = CellGeometryFunctions< KSpace, 5, KSpace::dimension>
765 ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
766 else trace.error() << "[DGtal::CellGeometry<TKSpace>::getTouchedKPoints]"
767 << " Computation are limited to n-D, n <= 5" << std::endl;
768 return std::vector< Cell >( kpoints.begin(), kpoints.end() );
769}
770
771
772//-----------------------------------------------------------------------------
773template <typename TKSpace>
774DGtal::Dimension
775DGtal::CellGeometry<TKSpace>::
776dim( const Point & kp )
777{
778 Dimension d = 0;
779 for ( Dimension i = 0; i < KSpace::dimension; ++i )
780 d += kp[ i ] & 1 ? 1 : 0;
781 return d;
782}
783
784///////////////////////////////////////////////////////////////////////////////
785// Interface - public :
786
787/**
788 * Writes/Displays the object on an output stream.
789 * @param out the output stream where the object is written.
790 */
791template <typename TKSpace>
792inline
793void
794DGtal::CellGeometry<TKSpace>::selfDisplay ( std::ostream & out ) const
795{
796 out << "[CellGeometry] ";
797 std::vector< Point > X;
798 for ( auto && c : myKPoints ) X.push_back( c );
799 std::sort( X.begin(), X.end() );
800 for ( auto p : X )
801 {
802 out << "(" << p[ 0 ];
803 for ( Dimension k = 1; k < TKSpace::dimension; ++k )
804 out << "," << p[ k ];
805 out << ") ";
806 }
807}
808
809/**
810 * Checks the validity/consistency of the object.
811 * @return 'true' if the object is valid, 'false' otherwise.
812 */
813template <typename TKSpace>
814inline
815bool
816DGtal::CellGeometry<TKSpace>::isValid() const
817{
818 return true;
819}
820//-----------------------------------------------------------------------------
821template <typename TKSpace>
822inline
823std::string
824DGtal::CellGeometry<TKSpace>::className
825() const
826{
827 return "CellGeometry";
828}
829
830
831
832///////////////////////////////////////////////////////////////////////////////
833// Implementation of inline functions //
834
835//-----------------------------------------------------------------------------
836template <typename TKSpace>
837inline
838std::ostream&
839DGtal::operator<< ( std::ostream & out,
840 const CellGeometry<TKSpace> & object )
841{
842 object.selfDisplay( out );
843 return out;
844}
845
846// //
847///////////////////////////////////////////////////////////////////////////////