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.
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.
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/>.
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
24 * Implementation of inline methods defined in CellGeometry.h
26 * This file is part of the DGtal library.
30 //////////////////////////////////////////////////////////////////////////////
34 //////////////////////////////////////////////////////////////////////////////
36 ///////////////////////////////////////////////////////////////////////////////
37 // IMPLEMENTATION of inline methods.
38 ///////////////////////////////////////////////////////////////////////////////
40 ///////////////////////////////////////////////////////////////////////////////
41 // ----------------------- Standard services ------------------------------
43 //-----------------------------------------------------------------------------
44 template <typename TKSpace>
45 DGtal::CellGeometry<TKSpace>::
48 myMinCellDim( 0 ), myMaxCellDim( KSpace::dimension ),
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 ),
62 ASSERT( 0 <= myMinCellDim );
63 ASSERT( myMinCellDim <= myMaxCellDim );
64 ASSERT( myMaxCellDim <= myK.dimension );
66 //-----------------------------------------------------------------------------
67 template <typename TKSpace>
69 DGtal::CellGeometry<TKSpace>::
70 init( const KSpace & K,
71 Dimension min_cell_dim, Dimension max_cell_dim, bool verbose )
73 ASSERT( 0 <= myMinCellDim );
74 ASSERT( myMinCellDim <= myMaxCellDim );
75 ASSERT( myMaxCellDim <= myK.dimension );
78 myMinCellDim = min_cell_dim;
79 myMaxCellDim = max_cell_dim;
83 //-----------------------------------------------------------------------------
84 template <typename TKSpace>
85 template <typename PointIterator>
87 DGtal::CellGeometry<TKSpace>::
88 addCellsTouchingPoints( PointIterator itB, PointIterator itE )
90 if ( ( myMinCellDim == 0 ) && ( myMaxCellDim == KSpace::dimension ) ) {
91 for ( auto it = itB; it != itE; ++it )
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 ) );
100 for ( auto it = itB; it != itE; ++it )
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 ) );
114 //-----------------------------------------------------------------------------
115 template <typename TKSpace>
116 template <typename PointelIterator>
118 DGtal::CellGeometry<TKSpace>::
119 addCellsTouchingPointels( PointelIterator itB, PointelIterator itE )
121 if ( ( myMinCellDim == 0 ) && ( myMaxCellDim == KSpace::dimension ) ) {
122 for ( auto it = itB; it != itE; ++it )
125 auto cofaces = myK.uCoFaces( pointel );
126 myKPoints.emplace( myK.uKCoords( pointel ) );
127 for ( auto && c : cofaces )
128 myKPoints.emplace( myK.uKCoords( c ) );
131 for ( auto it = itB; it != itE; ++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 ) );
144 //-----------------------------------------------------------------------------
145 template <typename TKSpace>
147 DGtal::CellGeometry<TKSpace>::
148 addCellsTouchingPolytopePoints( const LatticePolytope& polytope )
150 std::vector< Point > points;
151 polytope.getPoints( points );
152 addCellsTouchingPoints( points.begin(), points.end() );
154 //-----------------------------------------------------------------------------
155 template <typename TKSpace>
157 DGtal::CellGeometry<TKSpace>::
158 addCellsTouchingPolytopePoints( const RationalPolytope& polytope )
160 std::vector< Point > points;
161 polytope.getPoints( points );
162 addCellsTouchingPoints( points.begin(), points.end() );
164 //-----------------------------------------------------------------------------
165 template <typename TKSpace>
167 DGtal::CellGeometry<TKSpace>::
168 addCellsTouchingPolytope( const LatticePolytope& polytope )
170 for ( Dimension i = myMinCellDim; i <= myMaxCellDim; ++i )
172 auto kpoints = getIntersectedKPoints( polytope, i );
173 for ( auto&& c : kpoints )
174 myKPoints.insert( c );
177 //-----------------------------------------------------------------------------
178 template <typename TKSpace>
180 DGtal::CellGeometry<TKSpace>::
181 addCellsTouchingPolytope( const RationalPolytope& polytope )
183 for ( Dimension i = myMinCellDim; i <= myMaxCellDim; ++i )
185 auto kpoints = getIntersectedKPoints( polytope, i );
186 for ( auto&& c : kpoints )
187 myKPoints.insert( c );
190 //-----------------------------------------------------------------------------
191 template <typename TKSpace>
192 typename DGtal::CellGeometry<TKSpace>&
193 DGtal::CellGeometry<TKSpace>::
194 operator+=( const CellGeometry& other )
196 if ( this != &other )
198 for ( auto&& c : other.myKPoints )
199 myKPoints.insert( c );
200 myMinCellDim = std::min( myMinCellDim, other.myMinCellDim );
201 myMaxCellDim = std::max( myMaxCellDim, other.myMaxCellDim );
206 //-----------------------------------------------------------------------------
207 template <typename TKSpace>
208 typename DGtal::CellGeometry<TKSpace>::Size
209 DGtal::CellGeometry<TKSpace>::
212 return myKPoints.size();
214 //-----------------------------------------------------------------------------
215 template <typename TKSpace>
216 typename DGtal::CellGeometry<TKSpace>::Size
217 DGtal::CellGeometry<TKSpace>::
218 computeNbCells( const Dimension k ) const
220 if ( k < minCellDim() || k > maxCellDim() ) return 0;
222 for ( auto&& c : myKPoints )
223 nb += ( dim( c ) == k ) ? 1 : 0;
227 //-----------------------------------------------------------------------------
228 template <typename TKSpace>
229 typename DGtal::CellGeometry<TKSpace>::Integer
230 DGtal::CellGeometry<TKSpace>::
235 for ( Dimension k = 0; k <= maxCellDim(); ++k )
237 if ( pos ) euler += computeNbCells( k );
238 else euler -= computeNbCells( k );
243 //-----------------------------------------------------------------------------
244 template <typename TKSpace>
246 DGtal::CellGeometry<TKSpace>::
251 //-----------------------------------------------------------------------------
252 template <typename TKSpace>
254 DGtal::CellGeometry<TKSpace>::
260 //-----------------------------------------------------------------------------
261 template <typename TKSpace>
263 DGtal::CellGeometry<TKSpace>::
264 subset( const CellGeometry& other ) const
266 return other.myKPoints.includes( myKPoints );
267 // for ( auto&& c : myKPoints )
268 // if ( other.myKPoints.count( c ) == 0 )
272 //-----------------------------------------------------------------------------
273 template <typename TKSpace>
275 DGtal::CellGeometry<TKSpace>::
276 subset( const CellGeometry& other, const Dimension k ) const
278 UnorderedSetByBlock< Point, Splitter< Point, uint64_t > > k_dim_points;
279 for ( auto&& c : myKPoints )
281 k_dim_points.insert( c );
282 return other.myKPoints.includes( k_dim_points );
285 //-----------------------------------------------------------------------------
286 template <typename TKSpace>
287 template <typename RandomIterator>
289 DGtal::CellGeometry<TKSpace>::
290 includes( RandomIterator it1, RandomIterator itE1,
291 RandomIterator it2, RandomIterator itE2 )
294 for ( ; it2 != itE2; ++it1)
296 if (it1 == itE1 || *it2 < *it1) return false;
298 for ( k = 1; ( it1 < itE1 ) && ( *it1 < *it2 ); k *= 2 ) it1 += k;
300 if ( *it2 == *it1 ) ++it2; //equality
302 it1 = lower_bound( it1 - k/2, it1, *it2 );
303 if ( *it2 != *it1 ) return false;
307 it1 = lower_bound( it1 - k/2, itE1, *it2 );
308 if ( it1 == itE1 || *it2 != *it1 ) return false;
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
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 )
333 extensions[ k ] = k - 1;
334 polytopes [ k ] = polytopes[ k - 1 ] + typename LatticePolytope::UnitSegment( k - 1 );
336 // We have to build several dilated polytopes which corresponds
337 // to the binom(d,i) possible cell types.
343 std::ostringstream ostr( str );
344 ostr << "Dilating Polytope along directions {";
345 for ( Dimension k = 1; k < extensions.size(); ++k )
346 ostr << " + " << extensions[ k ];
348 trace.info() << ostr.str() << std::endl;
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 );
360 // Go to next type of cell
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 ] );
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
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 )
394 extensions[ k ] = k - 1;
395 polytopes [ k ] = polytopes[ k - 1 ] + typename LatticePolytope::UnitSegment( k - 1 );
397 // We have to build several dilated polytopes which corresponds
398 // to the binom(d,i) possible cell types.
404 std::ostringstream ostr( str );
405 ostr << "Dilating Polytope along directions {";
406 for ( Dimension k = 1; k < extensions.size(); ++k )
407 ostr << " + " << extensions[ k ];
409 trace.info() << ostr.str() << std::endl;
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 );
420 // Go to next type of cell
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 ] );
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
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 )
454 extensions[ k ] = k - 1;
455 polytopes [ k ] = polytopes[ k - 1 ]
456 + typename RationalPolytope::UnitSegment( k - 1 );
458 // We have to build several dilated polytopes which corresponds
459 // to the binom(d,i) possible cell types.
465 std::ostringstream ostr( str );
466 ostr << "Dilating Polytope along directions {";
467 for ( Dimension k = 1; k < extensions.size(); ++k )
468 ostr << " + " << extensions[ k ];
470 trace.info() << ostr.str() << std::endl;
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 );
481 // Go to next type of cell
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 ] );
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
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 )
515 extensions[ k ] = k - 1;
516 polytopes [ k ] = polytopes[ k - 1 ]
517 + typename RationalPolytope::UnitSegment( k - 1 );
519 // We have to build several dilated polytopes which corresponds
520 // to the binom(d,i) possible cell types.
526 std::ostringstream ostr( str );
527 ostr << "Dilating Polytope along directions {";
528 for ( Dimension k = 1; k < extensions.size(); ++k )
529 ostr << " + " << extensions[ k ];
531 trace.info() << ostr.str() << std::endl;
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 );
543 // Go to next type of cell
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 ] );
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
565 std::unordered_set< Cell > cells;
567 cells = CellGeometryFunctions< KSpace, 0, KSpace::dimension>
568 ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
570 cells = CellGeometryFunctions< KSpace, 1, KSpace::dimension>
571 ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
573 cells = CellGeometryFunctions< KSpace, 2, KSpace::dimension>
574 ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
576 cells = CellGeometryFunctions< KSpace, 3, KSpace::dimension>
577 ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
579 cells = CellGeometryFunctions< KSpace, 4, KSpace::dimension>
580 ::getIncidentCellsToPoints( myK, points.begin(), points.end() );
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() );
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
595 UnorderedSetByBlock< Point, Splitter< Point, uint64_t > > kpoints;
597 kpoints = CellGeometryFunctions< KSpace, 0, KSpace::dimension>
598 ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
600 kpoints = CellGeometryFunctions< KSpace, 1, KSpace::dimension>
601 ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
603 kpoints = CellGeometryFunctions< KSpace, 2, KSpace::dimension>
604 ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
606 kpoints = CellGeometryFunctions< KSpace, 3, KSpace::dimension>
607 ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
609 kpoints = CellGeometryFunctions< KSpace, 4, KSpace::dimension>
610 ::getIncidentKPointsToPoints( myK, points.begin(), points.end() );
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() );
620 //-----------------------------------------------------------------------------
621 template <typename TKSpace>
623 DGtal::CellGeometry<TKSpace>::
624 dim( const Point & kp )
627 for ( Dimension i = 0; i < KSpace::dimension; ++i )
628 d += kp[ i ] & 1 ? 1 : 0;
632 ///////////////////////////////////////////////////////////////////////////////
633 // Interface - public :
636 * Writes/Displays the object on an output stream.
637 * @param out the output stream where the object is written.
639 template <typename TKSpace>
642 DGtal::CellGeometry<TKSpace>::selfDisplay ( std::ostream & out ) const
644 out << "[CellGeometry]";
648 * Checks the validity/consistency of the object.
649 * @return 'true' if the object is valid, 'false' otherwise.
651 template <typename TKSpace>
654 DGtal::CellGeometry<TKSpace>::isValid() const
658 //-----------------------------------------------------------------------------
659 template <typename TKSpace>
662 DGtal::CellGeometry<TKSpace>::className
665 return "CellGeometry";
670 ///////////////////////////////////////////////////////////////////////////////
671 // Implementation of inline functions //
673 //-----------------------------------------------------------------------------
674 template <typename TKSpace>
677 DGtal::operator<< ( std::ostream & out,
678 const CellGeometry<TKSpace> & object )
680 object.selfDisplay( out );
685 ///////////////////////////////////////////////////////////////////////////////