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 DigitalConvexity.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 DigitalConvexity.h
26 * This file is part of the DGtal library.
30 //////////////////////////////////////////////////////////////////////////////
32 #include "DGtal/geometry/volumes/ConvexityHelper.h"
33 //////////////////////////////////////////////////////////////////////////////
35 //-----------------------------------------------------------------------------
36 template <typename TKSpace>
37 DGtal::DigitalConvexity<TKSpace>::
38 DigitalConvexity( Clone<KSpace> K, bool safe )
39 : myK( K ), mySafe( safe )
42 //-----------------------------------------------------------------------------
43 template <typename TKSpace>
44 DGtal::DigitalConvexity<TKSpace>::
45 DigitalConvexity( Point lo, Point hi, bool safe )
48 myK.init( lo, hi, true );
51 //-----------------------------------------------------------------------------
52 template <typename TKSpace>
54 DGtal::DigitalConvexity<TKSpace>::
60 //-----------------------------------------------------------------------------
61 template <typename TKSpace>
62 template <typename PointIterator>
63 typename DGtal::DigitalConvexity<TKSpace>::LatticePolytope
64 DGtal::DigitalConvexity<TKSpace>::
65 makeSimplex( PointIterator itB, PointIterator itE )
67 return LatticePolytope( itB, itE );
70 //-----------------------------------------------------------------------------
71 template <typename TKSpace>
72 typename DGtal::DigitalConvexity<TKSpace>::LatticePolytope
73 DGtal::DigitalConvexity<TKSpace>::
74 makeSimplex( std::initializer_list<Point> l )
76 return LatticePolytope( l );
79 //-----------------------------------------------------------------------------
80 template <typename TKSpace>
81 template <typename PointIterator>
82 typename DGtal::DigitalConvexity<TKSpace>::RationalPolytope
83 DGtal::DigitalConvexity<TKSpace>::
84 makeRationalSimplex( Integer d, PointIterator itB, PointIterator itE )
86 return RationalPolytope( d, itB, itE );
89 //-----------------------------------------------------------------------------
90 template <typename TKSpace>
91 typename DGtal::DigitalConvexity<TKSpace>::RationalPolytope
92 DGtal::DigitalConvexity<TKSpace>::
93 makeRationalSimplex( std::initializer_list<Point> l )
95 return RationalPolytope( l );
98 //-----------------------------------------------------------------------------
99 template <typename TKSpace>
100 template <typename PointIterator>
102 DGtal::DigitalConvexity<TKSpace>::
103 isSimplexFullDimensional( PointIterator itB, PointIterator itE )
105 typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
106 const Dimension d = KSpace::dimension;
107 std::vector<Point> pts( d+1 );
109 for ( ; itB != itE && k <= d; ++itB, ++k ) pts[ k ] = *itB;
110 // A simplex has exactly d+1 vertices.
111 if ( k != d+1 ) return false;
113 for ( Dimension i = 0; i < d; ++i )
114 for ( Dimension j = 0; j < d; ++j )
115 M.setComponent( i, j, pts[ i ][ j ] - pts[ d ][ j ] );
116 // A simplex has its vectors linearly independent.
117 return M.determinant() != 0;
120 //-----------------------------------------------------------------------------
121 template <typename TKSpace>
123 DGtal::DigitalConvexity<TKSpace>::
124 isSimplexFullDimensional( std::initializer_list<Point> l )
126 return isSimplexFullDimensional( l.begin(), l.end() );
129 //-----------------------------------------------------------------------------
130 template <typename TKSpace>
131 template <typename PointIterator>
132 typename DGtal::DigitalConvexity<TKSpace>::SimplexType
133 DGtal::DigitalConvexity<TKSpace>::
134 simplexType( PointIterator itB, PointIterator itE )
136 typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
137 const Dimension d = KSpace::dimension;
138 std::vector<Point> pts( d+1 );
140 for ( ; itB != itE && k <= d; ++itB, ++k ) pts[ k ] = *itB;
141 // A simplex has exactly d+1 vertices.
142 if ( k != d+1 ) return SimplexType::INVALID;
144 for ( Dimension i = 0; i < d; ++i )
145 for ( Dimension j = 0; j < d; ++j )
146 M.setComponent( i, j, pts[ i ][ j ] - pts[ d ][ j ] );
147 // A simplex has its vectors linearly independent.
148 auto V = M.determinant();
149 return (V == 0) ? SimplexType::DEGENERATED
150 : ( ((V == 1) || (V==-1)) ? SimplexType::UNITARY : SimplexType::COMMON );
153 //-----------------------------------------------------------------------------
154 template <typename TKSpace>
156 DGtal::DigitalConvexity<TKSpace>::
157 displaySimplex( std::ostream& out, std::initializer_list<Point> l )
159 displaySimplex( out, l.begin(), l.end() );
162 //-----------------------------------------------------------------------------
163 template <typename TKSpace>
164 template <typename PointIterator>
166 DGtal::DigitalConvexity<TKSpace>::
167 displaySimplex( std::ostream& out, PointIterator itB, PointIterator itE )
169 typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
170 const Dimension d = KSpace::dimension;
171 std::vector<Point> pts( d+1 );
173 for ( ; itB != itE && k <= d; ++itB, ++k ) pts[ k ] = *itB;
174 // A simplex has exactly d+1 vertices.
175 if ( k != d+1 ) { out << "[SPLX INVALID]"; return; }
177 for ( Dimension i = 0; i < d; ++i )
178 for ( Dimension kk = 0; kk < d; ++kk )
179 M.setComponent( i, kk, pts[ i ][ kk ] - pts[ d ][ kk ] );
180 // A simplex has its vectors linearly independent.
181 auto V = M.determinant();
182 out << "[SPLX V=" << V;
183 for ( Dimension i = 0; i < d; ++i ) {
185 for ( Dimension j = 0; j < d; ++j )
186 out << " " << M( i, j );
192 //-----------------------------------------------------------------------------
193 template <typename TKSpace>
194 typename DGtal::DigitalConvexity<TKSpace>::SimplexType
195 DGtal::DigitalConvexity<TKSpace>::
196 simplexType( std::initializer_list<Point> l )
198 return simplexType( l.begin(), l.end() );
202 //-----------------------------------------------------------------------------
203 template <typename TKSpace>
204 typename DGtal::DigitalConvexity<TKSpace>::PointRange
205 DGtal::DigitalConvexity<TKSpace>::
206 insidePoints( const LatticePolytope& polytope )
209 polytope.getPoints( pts );
212 //-----------------------------------------------------------------------------
213 template <typename TKSpace>
214 typename DGtal::DigitalConvexity<TKSpace>::PointRange
215 DGtal::DigitalConvexity<TKSpace>::
216 interiorPoints( const LatticePolytope& polytope )
219 polytope.getInteriorPoints( pts );
223 //-----------------------------------------------------------------------------
224 template <typename TKSpace>
225 typename DGtal::DigitalConvexity<TKSpace>::PointRange
226 DGtal::DigitalConvexity<TKSpace>::
227 insidePoints( const RationalPolytope& polytope )
230 polytope.getPoints( pts );
233 //-----------------------------------------------------------------------------
234 template <typename TKSpace>
235 typename DGtal::DigitalConvexity<TKSpace>::PointRange
236 DGtal::DigitalConvexity<TKSpace>::
237 interiorPoints( const RationalPolytope& polytope )
240 polytope.getInteriorPoints( pts );
244 //-----------------------------------------------------------------------------
245 template <typename TKSpace>
246 template <typename PointIterator>
247 typename DGtal::DigitalConvexity<TKSpace>::CellGeometry
248 DGtal::DigitalConvexity<TKSpace>::
249 makeCellCover( PointIterator itB, PointIterator itE,
250 Dimension i, Dimension k ) const
253 ASSERT( k <= KSpace::dimension );
254 CellGeometry cgeom( myK, i, k, false );
255 cgeom.addCellsTouchingPoints( itB, itE );
259 //-----------------------------------------------------------------------------
260 template <typename TKSpace>
261 typename DGtal::DigitalConvexity<TKSpace>::CellGeometry
262 DGtal::DigitalConvexity<TKSpace>::
263 makeCellCover( const LatticePolytope& P,
264 Dimension i, Dimension k ) const
267 ASSERT( k <= KSpace::dimension );
268 CellGeometry cgeom( myK, i, k, false );
269 cgeom.addCellsTouchingPolytope( P );
273 //-----------------------------------------------------------------------------
274 template <typename TKSpace>
275 typename DGtal::DigitalConvexity<TKSpace>::CellGeometry
276 DGtal::DigitalConvexity<TKSpace>::
277 makeCellCover( const RationalPolytope& P,
278 Dimension i, Dimension k ) const
281 ASSERT( k <= KSpace::dimension );
282 CellGeometry cgeom( myK, i, k, false );
283 cgeom.addCellsTouchingPolytope( P );
287 //-----------------------------------------------------------------------------
288 template <typename TKSpace>
289 typename DGtal::DigitalConvexity<TKSpace>::LatticePolytope
290 DGtal::DigitalConvexity<TKSpace>::
291 makePolytope( const PointRange& X, bool make_minkowski_summable ) const
295 typedef typename detail::ConvexityHelperInternalInteger< Integer, true >::Type
297 return ConvexityHelper< dimension, Integer, InternalInteger >::
298 computeLatticePolytope( X, false, make_minkowski_summable );
302 typedef typename detail::ConvexityHelperInternalInteger< Integer, false >::Type
304 return ConvexityHelper< dimension, Integer, InternalInteger >::
305 computeLatticePolytope( X, false, make_minkowski_summable );
309 //-----------------------------------------------------------------------------
310 template <typename TKSpace>
311 typename DGtal::DigitalConvexity<TKSpace>::PointRange
312 DGtal::DigitalConvexity<TKSpace>::
313 U( Dimension i, const PointRange& X ) const
317 Z.reserve( X.size() );
318 auto add_one = [&i] ( Point & p ) { p[ i ] += 1; };
319 std::for_each( Y.begin(), Y.end(), add_one );
320 std::set_union( X.cbegin(), X.cend(), Y.cbegin(), Y.cend(),
321 std::back_inserter( Z ) );
325 //-----------------------------------------------------------------------------
326 template <typename TKSpace>
328 DGtal::DigitalConvexity<TKSpace>::
329 is0Convex( const PointRange& X ) const
331 if ( X.empty() ) return true;
332 const auto P = makePolytope( X );
333 return P.count() == (DGtal::BoundedLatticePolytope<DGtal::SpaceND<3>>::Integer)X.size();
336 //-----------------------------------------------------------------------------
337 template <typename TKSpace>
339 DGtal::DigitalConvexity<TKSpace>::
340 isFullyConvex( const PointRange& Z, bool convex0 ) const
342 ASSERT( dimension <= 64 );
343 typedef DGtal::int64_t Direction;
344 typedef std::vector< Direction > Directions;
345 std::array< Directions, dimension > C;
346 const bool cvx0 = convex0 ? true : is0Convex( Z );
347 if ( ! cvx0 ) return false;
348 C[ 0 ].push_back( (Direction) 0 );
349 std::map< Direction, PointRange > X;
351 std::sort( X[ 0 ].begin(), X[ 0 ].end() );
352 for ( Dimension k = 1; k < dimension; k++ )
354 for ( const auto beta : C[ k-1 ] )
356 for ( Dimension j = 0; j < dimension; j++ )
358 const Direction dir_j = Direction(1) << j;
361 const Direction alpha = beta | dir_j;
362 C[ k ].push_back( alpha );
363 X[ alpha ] = U( j, X[ beta ] );
364 if ( ! is0Convex( X[ alpha ] ) ) return false;
372 //-----------------------------------------------------------------------------
373 template <typename TKSpace>
375 DGtal::DigitalConvexity<TKSpace>::
376 isFullyConvexFast( const PointRange& Z ) const
378 LatticeSet C_Z( Z.cbegin(), Z.cend(), 0 );
379 const auto nb_cells = C_Z.starOfPoints().size();
380 const auto s = sizeStarCvxH( Z );
381 return s == (Integer)nb_cells;
384 //-----------------------------------------------------------------------------
385 template <typename TKSpace>
386 typename DGtal::DigitalConvexity<TKSpace>::PointRange
387 DGtal::DigitalConvexity<TKSpace>::
388 ExtrCvxH( const PointRange& X ) const
392 typedef typename detail::ConvexityHelperInternalInteger< Integer, true >::Type
394 return ConvexityHelper< dimension, Integer, InternalInteger >::
395 computeLatticePolytopeVertices( X, false, false );
399 typedef typename detail::ConvexityHelperInternalInteger< Integer, false >::Type
401 return ConvexityHelper< dimension, Integer, InternalInteger >::
402 computeLatticePolytopeVertices( X, false, false );
406 //-----------------------------------------------------------------------------
407 template <typename TKSpace>
408 typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
409 DGtal::DigitalConvexity<TKSpace>::
410 StarCvxH( const PointRange& X, Dimension axis ) const
413 // Computes Minkowski sum of Z with hypercube
414 PointRange Z = U( 0, X );
415 for ( Dimension k = 1; k < dimension; k++ )
418 const auto P = makePolytope( Z );
419 // Extracts lattice points within polytope
420 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
422 const Dimension a = axis >= dimension ? C.longestAxis() : axis;
423 auto cellP = C.getLatticeCells( a );
424 return LatticeSet( cellP, a );
427 //-----------------------------------------------------------------------------
428 template <typename TKSpace>
429 typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
430 DGtal::DigitalConvexity<TKSpace>::
431 Star( const PointRange& X, const Dimension axis ) const
433 const Dimension a = axis >= dimension ? 0 : axis;
434 LatticeSet L( X.cbegin(), X.cend(), a );
435 return L.starOfPoints();
437 //-----------------------------------------------------------------------------
438 template <typename TKSpace>
439 typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
440 DGtal::DigitalConvexity<TKSpace>::
441 StarCells( const PointRange& C, const Dimension axis ) const
443 const Dimension a = axis >= dimension ? 0 : axis;
444 LatticeSet L( C.cbegin(), C.cend(), a );
445 return L.starOfCells();
448 //-----------------------------------------------------------------------------
449 template <typename TKSpace>
450 typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
451 DGtal::DigitalConvexity<TKSpace>::
452 toLatticeSet( const PointRange& X, const Dimension axis ) const
454 const Dimension a = axis >= dimension ? 0 : axis;
455 return LatticeSet( X.cbegin(), X.cend(), a );
458 //-----------------------------------------------------------------------------
459 template <typename TKSpace>
460 typename DGtal::DigitalConvexity<TKSpace>::PointRange
461 DGtal::DigitalConvexity<TKSpace>::
462 toPointRange( const LatticeSet& L ) const
464 return L.toPointRange();
467 //-----------------------------------------------------------------------------
468 template <typename TKSpace>
469 typename DGtal::DigitalConvexity<TKSpace>::Integer
470 DGtal::DigitalConvexity<TKSpace>::
471 sizeStarCvxH( const PointRange& X ) const
474 // Computes Minkowski sum of Z with hypercube
475 PointRange Z = U( 0, X );
476 for ( Dimension k = 1; k < dimension; k++ )
479 const auto P = makePolytope( Z );
480 // Extracts lattice points within polytope
481 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
483 const Dimension a = C.longestAxis();
484 auto cellP = C.getLatticeCells( a );
486 // Counts the number of cells
488 for ( const auto& value : cellP )
490 Point p = value.first;
491 Interval I = value.second;
492 nb += I.second - I.first + 1;
497 //-----------------------------------------------------------------------------
498 template <typename TKSpace>
499 typename DGtal::DigitalConvexity<TKSpace>::PointRange
500 DGtal::DigitalConvexity<TKSpace>::
501 Extr( const PointRange& C ) const
503 // JOL: using std::set< Point > or std::unordered_set< Point > is slightly slower.
504 // We prefer to use vector for easier vectorization.
505 std::vector<Point> E;
506 E.reserve( 2*C.size() );
507 for ( auto&& kp : C )
509 auto c = myK.uCell( kp );
510 if ( myK.uDim( c ) == 0 )
511 E.push_back( myK.uCoords( c ) );
514 auto faces = myK.uFaces( c );
515 for ( auto&& f : faces )
516 if ( myK.uDim( f ) == 0 )
517 E.push_back( myK.uCoords( f ) );
520 std::sort( E.begin(), E.end() );
521 auto last = std::unique( E.begin(), E.end() );
522 E.erase( last, E.end() );
525 //-----------------------------------------------------------------------------
526 template <typename TKSpace>
527 typename DGtal::DigitalConvexity<TKSpace>::PointRange
528 DGtal::DigitalConvexity<TKSpace>::
529 Extr( const LatticeSet& C ) const
531 return C.extremaOfCells();
533 //-----------------------------------------------------------------------------
534 template <typename TKSpace>
535 typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
536 DGtal::DigitalConvexity<TKSpace>::
537 Skel( const LatticeSet& C ) const
539 return C.skeletonOfCells();
541 //-----------------------------------------------------------------------------
542 template <typename TKSpace>
543 typename DGtal::DigitalConvexity<TKSpace>::PointRange
544 DGtal::DigitalConvexity<TKSpace>::
545 ExtrSkel( const LatticeSet& C ) const
547 return C.skeletonOfCells().extremaOfCells();
550 //-----------------------------------------------------------------------------
551 template <typename TKSpace>
552 typename DGtal::DigitalConvexity<TKSpace>::PointRange
553 DGtal::DigitalConvexity<TKSpace>::
554 FC_direct( const PointRange& Z ) const
556 typedef typename LatticePolytope::Domain Domain;
558 // Computes Minkowski sum of Z with hypercube
560 for ( Dimension k = 0; k < dimension; k++ )
563 const auto P = makePolytope( X );
564 // Extracts lattice points within polytope
565 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
567 const Dimension a = C.longestAxis();
568 Point lo = C.lowerBound();
569 Point hi = C.upperBound();
571 const Domain projD( lo, hi ); //< the projected domain of the polytope.
572 const Point One = Point::diagonal( 1 );
573 //const Size size = projD.size(); //not USED
574 std::unordered_map< Point, Interval > cellP;
576 for ( auto&& p : projD )
579 const auto I = C.intersectionIntervalAlongAxis( p, a );
580 const auto n = I.second - I.first;
583 // Now the second bound is included
584 cellP[ q ] = Interval( 2 * I.first - 1, 2 * I.second - 3 );
587 // It remains to compute all the k-cells, 0 <= k < d, intersected by Cvxh( Z )
588 for ( Dimension k = 0; k < dimension; k++ )
590 if ( k == a ) continue;
591 std::vector< Point > q_computed;
592 std::vector< Interval > I_computed;
593 for ( const auto& value : cellP )
595 Point p = value.first;
596 Interval I = value.second;
597 Point r = p; r[ k ] += 2;
598 const auto it = cellP.find( r );
599 if ( it == cellP.end() ) continue; // neighbor is empty
600 // Otherwise compute common part.
601 Interval J = it->second;
602 auto f = std::max( I.first, J.first );
603 auto s = std::min( I.second, J.second );
606 Point qq = p; qq[ k ] += 1;
607 q_computed.push_back( qq );
608 I_computed.push_back( Interval( f, s ) );
611 // Add new columns to map Point -> column
612 for ( size_t i = 0; i < q_computed.size(); ++i )
614 cellP[ q_computed[ i ] ] = I_computed[ i ];
617 // The built complex is open.
618 // Check it and fill skelP
619 std::unordered_map< Point, std::vector< Interval > > skelP;
620 for ( const auto& value : cellP )
622 Point p = value.first;
623 Interval I = value.second;
624 auto n = I.second - I.first + 1;
626 trace.error() << "Weird thickness step 1="
628 std::vector< Interval > V( 1, I );
629 auto result = skelP.insert( std::make_pair( p, V ) );
630 (void)result;//unused var;
633 // std::cout << "Extract skel" << std::endl;
634 // Now extracting implicitly its Skel
635 for ( const auto& value : cellP )
637 const Point& p = value.first;
638 const auto& I = value.second;
639 for ( Dimension k = 0; k < dimension; k++ )
641 if ( k == a ) continue;
642 if ( ( p[ k ] & 0x1 ) != 0 ) continue; // if open along axis continue
643 // if closed, check upper incident cells along direction k
644 Point qq = p; qq[ k ] -= 1;
645 Point r = p; r[ k ] += 1;
646 auto itq = skelP.find( qq );
647 if ( itq != skelP.end() )
649 auto& W = itq->second;
650 eraseInterval( I, W );
651 // if ( W.empty() ) skelP.erase( itq );
653 auto itr = skelP.find( r );
654 if ( itr != skelP.end() )
656 auto& W = itr->second;
657 eraseInterval( I, W );
658 // if ( W.empty() ) skelP.erase( itr );
662 // Extract skel along main axis
663 for ( const auto& value : cellP )
665 const Point& p = value.first;
666 auto I = value.second;
667 const auto itp = skelP.find( p );
668 if ( itp == skelP.end() ) continue;
669 if ( ( I.first & 0x1 ) != 0 ) I.first += 1;
670 if ( ( I.second & 0x1 ) != 0 ) I.second -= 1;
671 auto& W = itp->second;
672 for ( auto x = I.first; x <= I.second; x += 2 )
674 eraseInterval( Interval{ x-1,x-1} , W );
675 eraseInterval( Interval{ x+1,x+1} , W );
678 // Erase empty stacks
679 for ( auto it = skelP.begin(), itE = skelP.end(); it != itE; )
681 const auto& V = it->second;
686 skelP.erase( it_erase );
690 // Skel is constructed, now compute its Extr.
692 for ( const auto& value : skelP )
694 Point p = value.first;
695 auto W = value.second;
704 //-----------------------------------------------------------------------------
705 template <typename TKSpace>
706 typename DGtal::DigitalConvexity<TKSpace>::PointRange
707 DGtal::DigitalConvexity<TKSpace>::
708 FC_LatticeSet( const PointRange& Z ) const
710 const auto StarCvxZ = StarCvxH( Z );
711 const auto SkelStarCvxZ = StarCvxZ.skeletonOfCells().toPointRange();
712 return Extr( SkelStarCvxZ );
715 //-----------------------------------------------------------------------------
716 template <typename TKSpace>
717 typename DGtal::DigitalConvexity<TKSpace>::PointRange
718 DGtal::DigitalConvexity<TKSpace>::
719 FC( const PointRange& Z, EnvelopeAlgorithm algo ) const
721 if ( algo == EnvelopeAlgorithm::DIRECT )
722 return FC_direct( Z );
723 else if ( algo == EnvelopeAlgorithm::LATTICE_SET )
724 return FC_LatticeSet( Z );
729 //-----------------------------------------------------------------------------
730 template <typename TKSpace>
731 typename DGtal::DigitalConvexity<TKSpace>::PointRange
732 DGtal::DigitalConvexity<TKSpace>::
733 envelope( const PointRange& Z, EnvelopeAlgorithm algo ) const
738 auto card_In = In.size();
740 if ( In.size() == card_In ) return In;
743 trace.error() << "[DigitalConvexity::envelope] Should never pass here."
745 return Z; // to avoid warnings.
748 //-----------------------------------------------------------------------------
749 template <typename TKSpace>
750 typename DGtal::DigitalConvexity<TKSpace>::PointRange
751 DGtal::DigitalConvexity<TKSpace>::
752 relativeEnvelope( const PointRange& Z, const PointRange& Y,
753 EnvelopeAlgorithm algo ) const
759 std::set_intersection( In.cbegin(), In.cend(),
760 Y.cbegin(), Y.cend(),
761 std::back_inserter( Out ) );
762 In = FC( Out, algo );
763 if ( In.size() == Out.size() ) return Out;
769 //-----------------------------------------------------------------------------
770 template <typename TKSpace>
771 template <typename Predicate>
772 typename DGtal::DigitalConvexity<TKSpace>::PointRange
773 DGtal::DigitalConvexity<TKSpace>::
774 relativeEnvelope( const PointRange& Z, const Predicate& PredY,
775 EnvelopeAlgorithm algo ) const
780 auto Out = filter( In, PredY );
781 In = FC( Out, algo );
782 if ( In.size() == Out.size() ) return In;
788 //-----------------------------------------------------------------------------
789 template <typename TKSpace>
790 typename DGtal::DigitalConvexity<TKSpace>::Size
791 DGtal::DigitalConvexity<TKSpace>::
792 depthLastEnvelope() const
794 return myDepthLastFCE;
797 //-----------------------------------------------------------------------------
798 template <typename TKSpace>
800 DGtal::DigitalConvexity<TKSpace>::
801 isKConvex( const LatticePolytope& P, const Dimension k ) const
803 if ( k == 0 ) return true;
804 auto S = insidePoints( P );
805 auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
806 auto intersected_cells = makeCellCover( P, k, k );
807 return intersected_cells.nbCells() == touched_cells.nbCells();
808 // JOL: number should be enough
809 // && intersected_cells.subset( touched_cells );
812 //-----------------------------------------------------------------------------
813 template <typename TKSpace>
815 DGtal::DigitalConvexity<TKSpace>::
816 isFullyConvex( const LatticePolytope& P ) const
818 auto S = insidePoints( P );
819 for ( Dimension k = 1; k < KSpace::dimension; ++ k )
821 auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
822 auto intersected_cells = makeCellCover( P, k, k );
823 if ( ( intersected_cells.nbCells() != touched_cells.nbCells() ) )
824 // JOL: number should be enough
825 // || ( ! intersected_cells.subset( touched_cells ) ) )
831 //-----------------------------------------------------------------------------
832 template <typename TKSpace>
834 DGtal::DigitalConvexity<TKSpace>::
835 isKSubconvex( const LatticePolytope& P, const CellGeometry& C, const Dimension k ) const
837 auto intersected_cells = makeCellCover( P, k, k );
838 return intersected_cells.subset( C );
840 //-----------------------------------------------------------------------------
841 template <typename TKSpace>
843 DGtal::DigitalConvexity<TKSpace>::
844 isFullySubconvex( const LatticePolytope& P, const CellGeometry& C ) const
846 auto intersected_cells = makeCellCover( P, C.minCellDim(), C.maxCellDim() );
847 return intersected_cells.subset( C );
850 //-----------------------------------------------------------------------------
851 template <typename TKSpace>
853 DGtal::DigitalConvexity<TKSpace>::
854 isFullySubconvex( const LatticePolytope& P, const LatticeSet& StarX ) const
856 LatticePolytope Q = P + typename LatticePolytope::UnitSegment( 0 );
857 for ( Dimension k = 1; k < dimension; k++ )
858 Q = Q + typename LatticePolytope::UnitSegment( k );
860 const auto cells = C.getLatticeCells( StarX.axis() );
861 LatticeSet StarP( cells, StarX.axis() );
862 return StarX.includes( StarP );
865 //-----------------------------------------------------------------------------
866 template <typename TKSpace>
868 DGtal::DigitalConvexity<TKSpace>::
869 isFullySubconvex( const PointRange& Y, const LatticeSet& StarX ) const
871 const auto SCY = StarCvxH( Y, StarX.axis() );
872 return StarX.includes( SCY );
875 //-----------------------------------------------------------------------------
876 template <typename TKSpace>
878 DGtal::DigitalConvexity<TKSpace>::
879 isKSubconvex( const Point& a, const Point& b,
880 const CellGeometry& C, const Dimension k ) const
882 CellGeometry cgeom( myK, k, k, false );
883 cgeom.addCellsTouchingSegment( a, b );
884 return cgeom.subset( C );
887 //-----------------------------------------------------------------------------
888 template <typename TKSpace>
890 DGtal::DigitalConvexity<TKSpace>::
891 isFullySubconvex( const Point& a, const Point& b,
892 const CellGeometry& C ) const
894 CellGeometry cgeom( myK, C.minCellDim(), C.maxCellDim(), false );
895 cgeom.addCellsTouchingSegment( a, b );
896 return cgeom.subset( C );
899 //-----------------------------------------------------------------------------
900 template <typename TKSpace>
902 DGtal::DigitalConvexity<TKSpace>::
903 isFullySubconvex( const Point& a, const Point& b,
904 const LatticeSet& StarX ) const
907 LatticeSet L_ab( StarX.axis() );
908 const auto V = b - a;
910 // addCellsTouchingPoint( a );
911 for ( Dimension k = 0; k < dimension; k++ )
913 const Integer n = ( V[ k ] >= 0 ) ? V[ k ] : -V[ k ];
914 const Integer d = ( V[ k ] >= 0 ) ? 1 : -1;
915 if ( n == 0 ) continue;
917 for ( Integer i = 1; i < n; i++ )
919 for ( Dimension j = 0; j < dimension; j++ )
921 if ( j == k ) kc[ k ] = 2 * ( a[ k ] + d * i );
924 const auto v = V[ j ];
925 const auto q = ( v * i ) / n;
926 const auto r = ( v * i ) % n; // might be negative
927 kc[ j ] = 2 * ( a[ j ] + q );
928 if ( r < 0 ) kc[ j ] -= 1;
929 else if ( r > 0 ) kc[ j ] += 1;
935 if ( a != b ) L_ab.insert( 2*b );
936 LatticeSet Star_ab = L_ab.starOfCells();
937 return StarX.includes( Star_ab );
940 //-----------------------------------------------------------------------------
941 template <typename TKSpace>
943 DGtal::DigitalConvexity<TKSpace>::
944 isKConvex( const RationalPolytope& P, const Dimension k ) const
946 if ( k == 0 ) return true;
947 auto S = insidePoints( P );
948 auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
949 auto intersected_cells = makeCellCover( P, k, k );
950 return intersected_cells.nbCells() == touched_cells.nbCells()
951 && intersected_cells.subset( touched_cells );
954 //-----------------------------------------------------------------------------
955 template <typename TKSpace>
957 DGtal::DigitalConvexity<TKSpace>::
958 isFullyConvex( const RationalPolytope& P ) const
960 auto S = insidePoints( P );
961 for ( Dimension k = 1; k < KSpace::dimension; ++ k )
963 auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
964 auto intersected_cells = makeCellCover( P, k, k );
965 if ( ( intersected_cells.nbCells() != touched_cells.nbCells() )
966 || ( ! intersected_cells.subset( touched_cells ) ) )
972 //-----------------------------------------------------------------------------
973 template <typename TKSpace>
975 DGtal::DigitalConvexity<TKSpace>::
976 isKSubconvex( const RationalPolytope& P, const CellGeometry& C,
977 const Dimension k ) const
979 auto intersected_cells = makeCellCover( P, k, k );
980 return intersected_cells.subset( C );
982 //-----------------------------------------------------------------------------
983 template <typename TKSpace>
985 DGtal::DigitalConvexity<TKSpace>::
986 isFullySubconvex( const RationalPolytope& P, const CellGeometry& C ) const
988 auto intersected_cells = makeCellCover( P, C.minCellDim(), C.maxCellDim() );
989 return intersected_cells.subset( C );
992 //-----------------------------------------------------------------------------
993 template <typename TKSpace>
995 DGtal::DigitalConvexity<TKSpace>::
996 eraseInterval( Interval I, std::vector< Interval > & V )
998 for ( std::size_t i = 0; i < V.size(); )
1000 Interval& J = V[ i ];
1001 // I=[a,b], J=[a',b'], a <= b, a' <= b'
1002 if ( I.second < J.first ) { break; } // b < a' : no further intersection
1003 if ( J.second < I.first ) { ++i; continue; } // b' < a : no further intersection
1004 // a' <= b and a <= b'
1006 // a' ............... a'
1007 // b' ................. b'
1009 // a' ..................... b' => a'..a-1 b+1 b'
1010 Interval K1( J.first, I.first - 1 );
1011 Interval K2( I.second + 1, J.second );
1012 bool K1_exist = K1.second >= K1.first;
1013 bool K2_exist = K2.second >= K2.first;
1014 if ( K1_exist && K2_exist )
1017 V.insert( V.begin() + i, K1 );
1018 break; // no further intersection possible
1020 else if ( K1_exist )
1024 else if ( K2_exist )
1030 V.erase( V.begin() + i );
1035 //-----------------------------------------------------------------------------
1036 template <typename TKSpace>
1037 template <typename Predicate>
1038 typename DGtal::DigitalConvexity<TKSpace>::PointRange
1039 DGtal::DigitalConvexity<TKSpace>::
1040 filter( const PointRange& E, const Predicate& Pred )
1043 Out.reserve( E.size() );
1044 for ( auto&& p : E )
1045 if ( Pred( p ) ) Out.push_back( p );
1049 ///////////////////////////////////////////////////////////////////////////////
1050 // Interface - public :
1053 * Writes/Displays the object on an output stream.
1054 * @param out the output stream where the object is written.
1056 template <typename TKSpace>
1059 DGtal::DigitalConvexity<TKSpace>::selfDisplay ( std::ostream & out ) const
1061 out << "[DigitalConvexity]";
1065 * Checks the validity/consistency of the object.
1066 * @return 'true' if the object is valid, 'false' otherwise.
1068 template <typename TKSpace>
1071 DGtal::DigitalConvexity<TKSpace>::isValid() const
1077 ///////////////////////////////////////////////////////////////////////////////
1078 // Implementation of inline functions //
1080 //-----------------------------------------------------------------------------
1081 template <typename TKSpace>
1084 DGtal::operator<< ( std::ostream & out,
1085 const DigitalConvexity<TKSpace> & object )
1087 object.selfDisplay( out );
1092 ///////////////////////////////////////////////////////////////////////////////