DGtal 1.3.0
Loading...
Searching...
No Matches
DigitalConvexity.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 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
21 *
22 * @date 2020/01/31
23 *
24 * Implementation of inline methods defined in DigitalConvexity.h
25 *
26 * This file is part of the DGtal library.
27 */
28
29
30//////////////////////////////////////////////////////////////////////////////
31#include <cstdlib>
32#include "DGtal/geometry/volumes/ConvexityHelper.h"
33//////////////////////////////////////////////////////////////////////////////
34
35//-----------------------------------------------------------------------------
36template <typename TKSpace>
37DGtal::DigitalConvexity<TKSpace>::
38DigitalConvexity( Clone<KSpace> K, bool safe )
39 : myK( K ), mySafe( safe )
40{
41}
42//-----------------------------------------------------------------------------
43template <typename TKSpace>
44DGtal::DigitalConvexity<TKSpace>::
45DigitalConvexity( Point lo, Point hi, bool safe )
46 : mySafe( safe )
47{
48 myK.init( lo, hi, true );
49}
50
51//-----------------------------------------------------------------------------
52template <typename TKSpace>
53const TKSpace&
54DGtal::DigitalConvexity<TKSpace>::
55space() const
56{
57 return myK;
58}
59
60//-----------------------------------------------------------------------------
61template <typename TKSpace>
62template <typename PointIterator>
63typename DGtal::DigitalConvexity<TKSpace>::LatticePolytope
64DGtal::DigitalConvexity<TKSpace>::
65makeSimplex( PointIterator itB, PointIterator itE )
66{
67 return LatticePolytope( itB, itE );
68}
69
70//-----------------------------------------------------------------------------
71template <typename TKSpace>
72typename DGtal::DigitalConvexity<TKSpace>::LatticePolytope
73DGtal::DigitalConvexity<TKSpace>::
74makeSimplex( std::initializer_list<Point> l )
75{
76 return LatticePolytope( l );
77}
78
79//-----------------------------------------------------------------------------
80template <typename TKSpace>
81template <typename PointIterator>
82typename DGtal::DigitalConvexity<TKSpace>::RationalPolytope
83DGtal::DigitalConvexity<TKSpace>::
84makeRationalSimplex( Integer d, PointIterator itB, PointIterator itE )
85{
86 return RationalPolytope( d, itB, itE );
87}
88
89//-----------------------------------------------------------------------------
90template <typename TKSpace>
91typename DGtal::DigitalConvexity<TKSpace>::RationalPolytope
92DGtal::DigitalConvexity<TKSpace>::
93makeRationalSimplex( std::initializer_list<Point> l )
94{
95 return RationalPolytope( l );
96}
97
98//-----------------------------------------------------------------------------
99template <typename TKSpace>
100template <typename PointIterator>
101bool
102DGtal::DigitalConvexity<TKSpace>::
103isSimplexFullDimensional( PointIterator itB, PointIterator itE )
104{
105 typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
106 const Dimension d = KSpace::dimension;
107 std::vector<Point> pts( d+1 );
108 Dimension k = 0;
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;
112 Matrix M;
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;
118}
119
120//-----------------------------------------------------------------------------
121template <typename TKSpace>
122bool
123DGtal::DigitalConvexity<TKSpace>::
124isSimplexFullDimensional( std::initializer_list<Point> l )
125{
126 return isSimplexFullDimensional( l.begin(), l.end() );
127}
128
129//-----------------------------------------------------------------------------
130template <typename TKSpace>
131template <typename PointIterator>
132typename DGtal::DigitalConvexity<TKSpace>::SimplexType
133DGtal::DigitalConvexity<TKSpace>::
134simplexType( PointIterator itB, PointIterator itE )
135{
136 typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
137 const Dimension d = KSpace::dimension;
138 std::vector<Point> pts( d+1 );
139 Dimension k = 0;
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;
143 Matrix M;
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 );
151}
152
153//-----------------------------------------------------------------------------
154template <typename TKSpace>
155void
156DGtal::DigitalConvexity<TKSpace>::
157displaySimplex( std::ostream& out, std::initializer_list<Point> l )
158{
159 displaySimplex( out, l.begin(), l.end() );
160}
161
162//-----------------------------------------------------------------------------
163template <typename TKSpace>
164template <typename PointIterator>
165void
166DGtal::DigitalConvexity<TKSpace>::
167displaySimplex( std::ostream& out, PointIterator itB, PointIterator itE )
168{
169 typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
170 const Dimension d = KSpace::dimension;
171 std::vector<Point> pts( d+1 );
172 Dimension k = 0;
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; }
176 Matrix M;
177 for ( Dimension i = 0; i < d; ++i )
178 for ( Dimension k = 0; k < d; ++k )
179 M.setComponent( i, k, pts[ i ][ k ] - pts[ d ][ k ] );
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 ) {
184 out << " (";
185 for ( Dimension j = 0; j < d; ++j )
186 out << " " << M( i, j );
187 out << " )";
188 }
189 out << " ]";
190}
191
192//-----------------------------------------------------------------------------
193template <typename TKSpace>
194typename DGtal::DigitalConvexity<TKSpace>::SimplexType
195DGtal::DigitalConvexity<TKSpace>::
196simplexType( std::initializer_list<Point> l )
197{
198 return simplexType( l.begin(), l.end() );
199}
200
201
202//-----------------------------------------------------------------------------
203template <typename TKSpace>
204typename DGtal::DigitalConvexity<TKSpace>::PointRange
205DGtal::DigitalConvexity<TKSpace>::
206insidePoints( const LatticePolytope& polytope )
207{
208 PointRange pts;
209 polytope.getPoints( pts );
210 return pts;
211}
212//-----------------------------------------------------------------------------
213template <typename TKSpace>
214typename DGtal::DigitalConvexity<TKSpace>::PointRange
215DGtal::DigitalConvexity<TKSpace>::
216interiorPoints( const LatticePolytope& polytope )
217{
218 PointRange pts;
219 polytope.getInteriorPoints( pts );
220 return pts;
221}
222
223//-----------------------------------------------------------------------------
224template <typename TKSpace>
225typename DGtal::DigitalConvexity<TKSpace>::PointRange
226DGtal::DigitalConvexity<TKSpace>::
227insidePoints( const RationalPolytope& polytope )
228{
229 PointRange pts;
230 polytope.getPoints( pts );
231 return pts;
232}
233//-----------------------------------------------------------------------------
234template <typename TKSpace>
235typename DGtal::DigitalConvexity<TKSpace>::PointRange
236DGtal::DigitalConvexity<TKSpace>::
237interiorPoints( const RationalPolytope& polytope )
238{
239 PointRange pts;
240 polytope.getInteriorPoints( pts );
241 return pts;
242}
243
244//-----------------------------------------------------------------------------
245template <typename TKSpace>
246template <typename PointIterator>
247typename DGtal::DigitalConvexity<TKSpace>::CellGeometry
248DGtal::DigitalConvexity<TKSpace>::
249makeCellCover( PointIterator itB, PointIterator itE,
250 Dimension i, Dimension k ) const
251{
252 ASSERT( 0 <= i );
253 ASSERT( i <= k );
254 ASSERT( k <= KSpace::dimension );
255 CellGeometry cgeom( myK, i, k, false );
256 cgeom.addCellsTouchingPoints( itB, itE );
257 return cgeom;
258}
259
260//-----------------------------------------------------------------------------
261template <typename TKSpace>
262typename DGtal::DigitalConvexity<TKSpace>::CellGeometry
263DGtal::DigitalConvexity<TKSpace>::
264makeCellCover( const LatticePolytope& P,
265 Dimension i, Dimension k ) const
266{
267 ASSERT( 0 <= i );
268 ASSERT( i <= k );
269 ASSERT( k <= KSpace::dimension );
270 CellGeometry cgeom( myK, i, k, false );
271 cgeom.addCellsTouchingPolytope( P );
272 return cgeom;
273}
274
275//-----------------------------------------------------------------------------
276template <typename TKSpace>
277typename DGtal::DigitalConvexity<TKSpace>::CellGeometry
278DGtal::DigitalConvexity<TKSpace>::
279makeCellCover( const RationalPolytope& P,
280 Dimension i, Dimension k ) const
281{
282 ASSERT( 0 <= i );
283 ASSERT( i <= k );
284 ASSERT( k <= KSpace::dimension );
285 CellGeometry cgeom( myK, i, k, false );
286 cgeom.addCellsTouchingPolytope( P );
287 return cgeom;
288}
289
290//-----------------------------------------------------------------------------
291template <typename TKSpace>
292typename DGtal::DigitalConvexity<TKSpace>::LatticePolytope
293DGtal::DigitalConvexity<TKSpace>::
294makePolytope( const PointRange& X, bool make_minkowski_summable ) const
295{
296 if ( mySafe )
297 {
298 typedef typename detail::ConvexityHelperInternalInteger< Integer, true >::Type
299 InternalInteger;
300 return ConvexityHelper< dimension, Integer, InternalInteger >::
301 computeLatticePolytope( X, false, make_minkowski_summable );
302 }
303 else
304 {
305 typedef typename detail::ConvexityHelperInternalInteger< Integer, false >::Type
306 InternalInteger;
307 return ConvexityHelper< dimension, Integer, InternalInteger >::
308 computeLatticePolytope( X, false, make_minkowski_summable );
309 }
310}
311
312//-----------------------------------------------------------------------------
313template <typename TKSpace>
314typename DGtal::DigitalConvexity<TKSpace>::PointRange
315DGtal::DigitalConvexity<TKSpace>::
316U( Dimension i, const PointRange& X ) const
317{
318 PointRange Y( X );
319 PointRange Z;
320 Z.reserve( X.size() );
321 auto add_one = [&i] ( Point & p ) { p[ i ] += 1; };
322 std::for_each( Y.begin(), Y.end(), add_one );
323 std::set_union( X.cbegin(), X.cend(), Y.cbegin(), Y.cend(),
324 std::back_inserter( Z ) );
325 return Z;
326}
327
328//-----------------------------------------------------------------------------
329template <typename TKSpace>
330bool
331DGtal::DigitalConvexity<TKSpace>::
332is0Convex( const PointRange& X ) const
333{
334 if ( X.empty() ) return true;
335 const auto P = makePolytope( X );
336 return P.count() == X.size();
337}
338
339//-----------------------------------------------------------------------------
340template <typename TKSpace>
341bool
342DGtal::DigitalConvexity<TKSpace>::
343isFullyConvex( const PointRange& Z, bool convex0 ) const
344{
345 ASSERT( dimension <= 64 );
346 typedef DGtal::int64_t Direction;
347 typedef std::vector< Direction > Directions;
348 std::array< Directions, dimension > C;
349 const bool cvx0 = convex0 ? true : is0Convex( Z );
350 if ( ! cvx0 ) return false;
351 C[ 0 ].push_back( (Direction) 0 );
352 std::map< Direction, PointRange > X;
353 X[ 0 ] = Z;
354 std::sort( X[ 0 ].begin(), X[ 0 ].end() );
355 for ( Dimension k = 1; k < dimension; k++ )
356 {
357 for ( const auto beta : C[ k-1 ] )
358 {
359 for ( Dimension j = 0; j < dimension; j++ )
360 {
361 const Direction dir_j = Direction(1) << j;
362 if ( beta < dir_j )
363 {
364 const Direction alpha = beta | dir_j;
365 C[ k ].push_back( alpha );
366 X[ alpha ] = U( j, X[ beta ] );
367 if ( ! is0Convex( X[ alpha ] ) ) return false;
368 }
369 }
370 }
371 }
372 return true;
373}
374
375//-----------------------------------------------------------------------------
376template <typename TKSpace>
377bool
378DGtal::DigitalConvexity<TKSpace>::
379isFullyConvexFast( const PointRange& Z ) const
380{
381 LatticeSet C_Z( Z.cbegin(), Z.cend(), 0 );
382 const auto nb_cells = C_Z.starOfPoints().size();
383 const auto s = sizeStarCvxH( Z );
384 return s == nb_cells;
385}
386
387//-----------------------------------------------------------------------------
388template <typename TKSpace>
389typename DGtal::DigitalConvexity<TKSpace>::PointRange
390DGtal::DigitalConvexity<TKSpace>::
391ExtrCvxH( const PointRange& X ) const
392{
393 if ( mySafe )
394 {
395 typedef typename detail::ConvexityHelperInternalInteger< Integer, true >::Type
396 InternalInteger;
397 return ConvexityHelper< dimension, Integer, InternalInteger >::
398 computeLatticePolytopeVertices( X, false, false );
399 }
400 else
401 {
402 typedef typename detail::ConvexityHelperInternalInteger< Integer, false >::Type
403 InternalInteger;
404 return ConvexityHelper< dimension, Integer, InternalInteger >::
405 computeLatticePolytopeVertices( X, false, false );
406 }
407}
408
409//-----------------------------------------------------------------------------
410template <typename TKSpace>
411typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
412DGtal::DigitalConvexity<TKSpace>::
413StarCvxH( const PointRange& X, Dimension axis ) const
414{
415 PointRange cells;
416 typedef BoundedLatticePolytopeCounter<Space> Counter;
417 typedef typename Counter::Interval Interval;
418 // Computes Minkowski sum of Z with hypercube
419 PointRange Z = U( 0, X );
420 for ( Dimension k = 1; k < dimension; k++ )
421 Z = U( k, Z );
422 // Builds polytope
423 const auto P = makePolytope( Z );
424 // Extracts lattice points within polytope
425 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
426 Counter C( P );
427 const Dimension a = axis >= dimension ? C.longestAxis() : axis;
428 auto cellP = C.getLatticeCells( a );
429 return LatticeSet( cellP, a );
430}
431
432//-----------------------------------------------------------------------------
433template <typename TKSpace>
434typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
435DGtal::DigitalConvexity<TKSpace>::
436Star( const PointRange& X, const Dimension axis ) const
437{
438 const Dimension a = axis >= dimension ? 0 : axis;
439 LatticeSet L( X.cbegin(), X.cend(), a );
440 return L.starOfPoints();
441}
442//-----------------------------------------------------------------------------
443template <typename TKSpace>
444typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
445DGtal::DigitalConvexity<TKSpace>::
446StarCells( const PointRange& C, const Dimension axis ) const
447{
448 const Dimension a = axis >= dimension ? 0 : axis;
449 LatticeSet L( C.cbegin(), C.cend(), a );
450 return L.starOfCells();
451}
452
453//-----------------------------------------------------------------------------
454template <typename TKSpace>
455typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
456DGtal::DigitalConvexity<TKSpace>::
457toLatticeSet( const PointRange& X, const Dimension axis ) const
458{
459 const Dimension a = axis >= dimension ? 0 : axis;
460 return LatticeSet( X.cbegin(), X.cend(), a );
461}
462
463//-----------------------------------------------------------------------------
464template <typename TKSpace>
465typename DGtal::DigitalConvexity<TKSpace>::PointRange
466DGtal::DigitalConvexity<TKSpace>::
467toPointRange( const LatticeSet& L ) const
468{
469 return L.toPointRange();
470}
471
472//-----------------------------------------------------------------------------
473template <typename TKSpace>
474typename DGtal::DigitalConvexity<TKSpace>::Integer
475DGtal::DigitalConvexity<TKSpace>::
476sizeStarCvxH( const PointRange& X ) const
477{
478 PointRange cells;
479 typedef BoundedLatticePolytopeCounter<Space> Counter;
480 typedef typename Counter::Interval Interval;
481 // Computes Minkowski sum of Z with hypercube
482 PointRange Z = U( 0, X );
483 for ( Dimension k = 1; k < dimension; k++ )
484 Z = U( k, Z );
485 // Builds polytope
486 const auto P = makePolytope( Z );
487 // Extracts lattice points within polytope
488 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
489 Counter C( P );
490 const Dimension a = C.longestAxis();
491 auto cellP = C.getLatticeCells( a );
492
493 // Counts the number of cells
494 Integer nb = 0;
495 for ( const auto& value : cellP )
496 {
497 Point p = value.first;
498 Interval I = value.second;
499 nb += I.second - I.first + 1;
500 }
501 return nb;
502}
503
504//-----------------------------------------------------------------------------
505template <typename TKSpace>
506typename DGtal::DigitalConvexity<TKSpace>::PointRange
507DGtal::DigitalConvexity<TKSpace>::
508Extr( const PointRange& C ) const
509{
510 // JOL: using std::set< Point > or std::unordered_set< Point > is slightly slower.
511 // We prefer to use vector for easier vectorization.
512 std::vector<Point> E;
513 E.reserve( 2*C.size() );
514 for ( auto&& kp : C )
515 {
516 auto c = myK.uCell( kp );
517 if ( myK.uDim( c ) == 0 )
518 E.push_back( myK.uCoords( c ) );
519 else
520 {
521 auto faces = myK.uFaces( c );
522 for ( auto&& f : faces )
523 if ( myK.uDim( f ) == 0 )
524 E.push_back( myK.uCoords( f ) );
525 }
526 }
527 std::sort( E.begin(), E.end() );
528 auto last = std::unique( E.begin(), E.end() );
529 E.erase( last, E.end() );
530 return E;
531}
532//-----------------------------------------------------------------------------
533template <typename TKSpace>
534typename DGtal::DigitalConvexity<TKSpace>::PointRange
535DGtal::DigitalConvexity<TKSpace>::
536Extr( const LatticeSet& C ) const
537{
538 return C.extremaOfCells();
539}
540//-----------------------------------------------------------------------------
541template <typename TKSpace>
542typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
543DGtal::DigitalConvexity<TKSpace>::
544Skel( const LatticeSet& C ) const
545{
546 return C.skeletonOfCells();
547}
548//-----------------------------------------------------------------------------
549template <typename TKSpace>
550typename DGtal::DigitalConvexity<TKSpace>::PointRange
551DGtal::DigitalConvexity<TKSpace>::
552ExtrSkel( const LatticeSet& C ) const
553{
554 return C.skeletonOfCells().extremaOfCells();
555}
556
557//-----------------------------------------------------------------------------
558template <typename TKSpace>
559typename DGtal::DigitalConvexity<TKSpace>::PointRange
560DGtal::DigitalConvexity<TKSpace>::
561FC_direct( const PointRange& Z ) const
562{
563 typedef std::size_t Size;
564 typedef typename LatticePolytope::Domain Domain;
565 typedef typename Counter::Interval Interval;
566 PointRange cells;
567 // Computes Minkowski sum of Z with hypercube
568 PointRange X( Z );
569 for ( Dimension k = 0; k < dimension; k++ )
570 X = U( k, X );
571 // Builds polytope
572 const auto P = makePolytope( X );
573 // Extracts lattice points within polytope
574 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
575 Counter C( P );
576 const Dimension a = C.longestAxis();
577 Point lo = C.lowerBound();
578 Point hi = C.upperBound();
579 hi[ a ] = lo[ a ];
580 const Domain projD( lo, hi ); //< the projected domain of the polytope.
581 const Point One = Point::diagonal( 1 );
582 const Size size = projD.size();
583 std::unordered_map< Point, Interval > cellP;
584 Point q;
585 for ( auto&& p : projD )
586 {
587 q = 2*p - One;
588 const auto I = C.intersectionIntervalAlongAxis( p, a );
589 const auto n = I.second - I.first;
590 if ( n != 0 )
591 {
592 // Now the second bound is included
593 cellP[ q ] = Interval( 2 * I.first - 1, 2 * I.second - 3 );
594 }
595 }
596 // It remains to compute all the k-cells, 0 <= k < d, intersected by Cvxh( Z )
597 for ( Dimension k = 0; k < dimension; k++ )
598 {
599 if ( k == a ) continue;
600 std::vector< Point > q_computed;
601 std::vector< Interval > I_computed;
602 for ( const auto& value : cellP )
603 {
604 Point p = value.first;
605 Interval I = value.second;
606 Point r = p; r[ k ] += 2;
607 const auto it = cellP.find( r );
608 if ( it == cellP.end() ) continue; // neighbor is empty
609 // Otherwise compute common part.
610 Interval J = it->second;
611 auto f = std::max( I.first, J.first );
612 auto s = std::min( I.second, J.second );
613 if ( f <= s )
614 {
615 Point q = p; q[ k ] += 1;
616 q_computed.push_back( q );
617 I_computed.push_back( Interval( f, s ) );
618 }
619 }
620 // Add new columns to map Point -> column
621 for ( auto i = 0; i < q_computed.size(); ++i )
622 {
623 cellP[ q_computed[ i ] ] = I_computed[ i ];
624 }
625 }
626 // The built complex is open.
627 // Check it and fill skelP
628 std::unordered_map< Point, std::vector< Interval > > skelP;
629 for ( const auto& value : cellP )
630 {
631 Point p = value.first;
632 Interval I = value.second;
633 auto n = I.second - I.first + 1;
634 if ( n % 2 == 0 )
635 trace.error() << "Weird thickness step 1="
636 << n << std::endl;
637 std::vector< Interval > V( 1, I );
638 auto result = skelP.insert( std::make_pair( p, V ) );
639 }
640
641 // std::cout << "Extract skel" << std::endl;
642 // Now extracting implicitly its Skel
643 for ( const auto& value : cellP )
644 {
645 const Point& p = value.first;
646 const auto& I = value.second;
647 for ( Dimension k = 0; k < dimension; k++ )
648 {
649 if ( k == a ) continue;
650 if ( ( p[ k ] & 0x1 ) != 0 ) continue; // if open along axis continue
651 // if closed, check upper incident cells along direction k
652 Point q = p; q[ k ] -= 1;
653 Point r = p; r[ k ] += 1;
654 auto itq = skelP.find( q );
655 if ( itq != skelP.end() )
656 {
657 auto& W = itq->second;
658 eraseInterval( I, W );
659 // if ( W.empty() ) skelP.erase( itq );
660 }
661 auto itr = skelP.find( r );
662 if ( itr != skelP.end() )
663 {
664 auto& W = itr->second;
665 eraseInterval( I, W );
666 // if ( W.empty() ) skelP.erase( itr );
667 }
668 }
669 }
670 // Extract skel along main axis
671 for ( const auto& value : cellP )
672 {
673 const Point& p = value.first;
674 auto I = value.second;
675 const auto itp = skelP.find( p );
676 if ( itp == skelP.end() ) continue;
677 if ( ( I.first & 0x1 ) != 0 ) I.first += 1;
678 if ( ( I.second & 0x1 ) != 0 ) I.second -= 1;
679 auto& W = itp->second;
680 for ( auto x = I.first; x <= I.second; x += 2 )
681 {
682 eraseInterval( Interval{ x-1,x-1} , W );
683 eraseInterval( Interval{ x+1,x+1} , W );
684 }
685 }
686 // Erase empty stacks
687 for ( auto it = skelP.begin(), itE = skelP.end(); it != itE; )
688 {
689 const auto& V = it->second;
690 if ( V.empty() )
691 {
692 auto it_erase = it;
693 ++it;
694 skelP.erase( it_erase );
695 }
696 else ++it;
697 }
698 // Skel is constructed, now compute its Extr.
699 PointRange O;
700 for ( const auto& value : skelP )
701 {
702 Point p = value.first;
703 auto W = value.second;
704 for ( auto&& I : W )
705 {
706 p[ a ] = I.first;
707 O.push_back( p );
708 }
709 }
710 return Extr( O );
711}
712//-----------------------------------------------------------------------------
713template <typename TKSpace>
714typename DGtal::DigitalConvexity<TKSpace>::PointRange
715DGtal::DigitalConvexity<TKSpace>::
716FC_LatticeSet( const PointRange& Z ) const
717{
718 const auto StarCvxZ = StarCvxH( Z );
719 const auto SkelStarCvxZ = StarCvxZ.skeletonOfCells().toPointRange();
720 return Extr( SkelStarCvxZ );
721}
722
723//-----------------------------------------------------------------------------
724template <typename TKSpace>
725typename DGtal::DigitalConvexity<TKSpace>::PointRange
726DGtal::DigitalConvexity<TKSpace>::
727FC( const PointRange& Z, EnvelopeAlgorithm algo ) const
728{
729 if ( algo == EnvelopeAlgorithm::DIRECT )
730 return FC_direct( Z );
731 else if ( algo == EnvelopeAlgorithm::LATTICE_SET )
732 return FC_LatticeSet( Z );
733 else
734 return Z;
735}
736
737//-----------------------------------------------------------------------------
738template <typename TKSpace>
739typename DGtal::DigitalConvexity<TKSpace>::PointRange
740DGtal::DigitalConvexity<TKSpace>::
741envelope( const PointRange& Z, EnvelopeAlgorithm algo ) const
742{
743 myDepthLastFCE = 0;
744 auto In = Z;
745 while (true) {
746 auto card_In = In.size();
747 In = FC( In, algo );
748 if ( In.size() == card_In ) return In;
749 myDepthLastFCE++;
750 }
751 trace.error() << "[DigitalConvexity::envelope] Should never pass here."
752 << std::endl;
753 return Z; // to avoid warnings.
754}
755
756//-----------------------------------------------------------------------------
757template <typename TKSpace>
758typename DGtal::DigitalConvexity<TKSpace>::PointRange
759DGtal::DigitalConvexity<TKSpace>::
760relativeEnvelope( const PointRange& Z, const PointRange& Y,
761 EnvelopeAlgorithm algo ) const
762{
763 myDepthLastFCE = 0;
764 auto In = Z;
765 while (true) {
766 PointRange Out;
767 std::set_intersection( In.cbegin(), In.cend(),
768 Y.cbegin(), Y.cend(),
769 std::back_inserter( Out ) );
770 In = FC( Out, algo );
771 if ( In.size() == Out.size() ) return Out;
772 myDepthLastFCE++;
773 }
774 return Z;
775}
776
777//-----------------------------------------------------------------------------
778template <typename TKSpace>
779template <typename Predicate>
780typename DGtal::DigitalConvexity<TKSpace>::PointRange
781DGtal::DigitalConvexity<TKSpace>::
782relativeEnvelope( const PointRange& Z, const Predicate& PredY,
783 EnvelopeAlgorithm algo ) const
784{
785 myDepthLastFCE = 0;
786 auto In = Z;
787 while (true) {
788 auto Out = filter( In, PredY );
789 In = FC( Out, algo );
790 if ( In.size() == Out.size() ) return In;
791 myDepthLastFCE++;
792 }
793 return Z;
794}
795
796//-----------------------------------------------------------------------------
797template <typename TKSpace>
798typename DGtal::DigitalConvexity<TKSpace>::Size
799DGtal::DigitalConvexity<TKSpace>::
800depthLastEnvelope() const
801{
802 return myDepthLastFCE;
803}
804
805//-----------------------------------------------------------------------------
806template <typename TKSpace>
807bool
808DGtal::DigitalConvexity<TKSpace>::
809isKConvex( const LatticePolytope& P, const Dimension k ) const
810{
811 if ( k == 0 ) return true;
812 auto S = insidePoints( P );
813 auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
814 auto intersected_cells = makeCellCover( P, k, k );
815 return intersected_cells.nbCells() == touched_cells.nbCells();
816 // JOL: number should be enough
817 // && intersected_cells.subset( touched_cells );
818}
819
820//-----------------------------------------------------------------------------
821template <typename TKSpace>
822bool
823DGtal::DigitalConvexity<TKSpace>::
824isFullyConvex( const LatticePolytope& P ) const
825{
826 auto S = insidePoints( P );
827 for ( Dimension k = 1; k < KSpace::dimension; ++ k )
828 {
829 auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
830 auto intersected_cells = makeCellCover( P, k, k );
831 if ( ( intersected_cells.nbCells() != touched_cells.nbCells() ) )
832 // JOL: number should be enough
833 // || ( ! intersected_cells.subset( touched_cells ) ) )
834 return false;
835 }
836 return true;
837}
838
839//-----------------------------------------------------------------------------
840template <typename TKSpace>
841bool
842DGtal::DigitalConvexity<TKSpace>::
843isKSubconvex( const LatticePolytope& P, const CellGeometry& C, const Dimension k ) const
844{
845 auto intersected_cells = makeCellCover( P, k, k );
846 return intersected_cells.subset( C );
847}
848//-----------------------------------------------------------------------------
849template <typename TKSpace>
850bool
851DGtal::DigitalConvexity<TKSpace>::
852isFullySubconvex( const LatticePolytope& P, const CellGeometry& C ) const
853{
854 auto intersected_cells = makeCellCover( P, C.minCellDim(), C.maxCellDim() );
855 return intersected_cells.subset( C );
856}
857
858//-----------------------------------------------------------------------------
859template <typename TKSpace>
860bool
861DGtal::DigitalConvexity<TKSpace>::
862isFullySubconvex( const LatticePolytope& P, const LatticeSet& StarX ) const
863{
864 LatticePolytope Q = P + typename LatticePolytope::UnitSegment( 0 );
865 for ( Dimension k = 1; k < dimension; k++ )
866 Q = Q + typename LatticePolytope::UnitSegment( k );
867 typedef BoundedLatticePolytopeCounter<Space> Counter;
868 Counter C( Q );
869 const auto cells = C.getLatticeCells( StarX.axis() );
870 LatticeSet StarP( cells, StarX.axis() );
871 return StarX.includes( StarP );
872}
873
874//-----------------------------------------------------------------------------
875template <typename TKSpace>
876bool
877DGtal::DigitalConvexity<TKSpace>::
878isFullySubconvex( const PointRange& Y, const LatticeSet& StarX ) const
879{
880 const auto SCY = StarCvxH( Y, StarX.axis() );
881 return StarX.includes( SCY );
882}
883
884//-----------------------------------------------------------------------------
885template <typename TKSpace>
886bool
887DGtal::DigitalConvexity<TKSpace>::
888isKSubconvex( const Point& a, const Point& b,
889 const CellGeometry& C, const Dimension k ) const
890{
891 CellGeometry cgeom( myK, k, k, false );
892 cgeom.addCellsTouchingSegment( a, b );
893 return cgeom.subset( C );
894}
895
896//-----------------------------------------------------------------------------
897template <typename TKSpace>
898bool
899DGtal::DigitalConvexity<TKSpace>::
900isFullySubconvex( const Point& a, const Point& b,
901 const CellGeometry& C ) const
902{
903 CellGeometry cgeom( myK, C.minCellDim(), C.maxCellDim(), false );
904 cgeom.addCellsTouchingSegment( a, b );
905 return cgeom.subset( C );
906}
907
908//-----------------------------------------------------------------------------
909template <typename TKSpace>
910bool
911DGtal::DigitalConvexity<TKSpace>::
912isFullySubconvex( const Point& a, const Point& b,
913 const LatticeSet& StarX ) const
914{
915 // TODO
916 LatticeSet L_ab( StarX.axis() );
917 const auto V = b - a;
918 L_ab.insert( 2*a );
919 // addCellsTouchingPoint( a );
920 for ( Integer k = 0; k < dimension; k++ )
921 {
922 const Integer n = ( V[ k ] >= 0 ) ? V[ k ] : -V[ k ];
923 const Integer d = ( V[ k ] >= 0 ) ? 1 : -1;
924 if ( n == 0 ) continue;
925 Point kc;
926 for ( Integer i = 1; i < n; i++ )
927 {
928 for ( Dimension j = 0; j < dimension; j++ )
929 {
930 if ( j == k ) kc[ k ] = 2 * ( a[ k ] + d * i );
931 else
932 {
933 const auto v = V[ j ];
934 const auto q = ( v * i ) / n;
935 const auto r = ( v * i ) % n; // might be negative
936 kc[ j ] = 2 * ( a[ j ] + q );
937 if ( r < 0 ) kc[ j ] -= 1;
938 else if ( r > 0 ) kc[ j ] += 1;
939 }
940 }
941 L_ab.insert( kc );
942 }
943 }
944 if ( a != b ) L_ab.insert( 2*b );
945 LatticeSet Star_ab = L_ab.starOfCells();
946 return StarX.includes( Star_ab );
947}
948
949//-----------------------------------------------------------------------------
950template <typename TKSpace>
951bool
952DGtal::DigitalConvexity<TKSpace>::
953isKConvex( const RationalPolytope& P, const Dimension k ) const
954{
955 if ( k == 0 ) return true;
956 auto S = insidePoints( P );
957 auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
958 auto intersected_cells = makeCellCover( P, k, k );
959 return intersected_cells.nbCells() == touched_cells.nbCells()
960 && intersected_cells.subset( touched_cells );
961}
962
963//-----------------------------------------------------------------------------
964template <typename TKSpace>
965bool
966DGtal::DigitalConvexity<TKSpace>::
967isFullyConvex( const RationalPolytope& P ) const
968{
969 auto S = insidePoints( P );
970 for ( Dimension k = 1; k < KSpace::dimension; ++ k )
971 {
972 auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
973 auto intersected_cells = makeCellCover( P, k, k );
974 if ( ( intersected_cells.nbCells() != touched_cells.nbCells() )
975 || ( ! intersected_cells.subset( touched_cells ) ) )
976 return false;
977 }
978 return true;
979}
980
981//-----------------------------------------------------------------------------
982template <typename TKSpace>
983bool
984DGtal::DigitalConvexity<TKSpace>::
985isKSubconvex( const RationalPolytope& P, const CellGeometry& C,
986 const Dimension k ) const
987{
988 auto intersected_cells = makeCellCover( P, k, k );
989 return intersected_cells.subset( C );
990}
991//-----------------------------------------------------------------------------
992template <typename TKSpace>
993bool
994DGtal::DigitalConvexity<TKSpace>::
995isFullySubconvex( const RationalPolytope& P, const CellGeometry& C ) const
996{
997 auto intersected_cells = makeCellCover( P, C.minCellDim(), C.maxCellDim() );
998 return intersected_cells.subset( C );
999}
1000
1001//-----------------------------------------------------------------------------
1002template <typename TKSpace>
1003void
1004DGtal::DigitalConvexity<TKSpace>::
1005eraseInterval( Interval I, std::vector< Interval > & V )
1006{
1007 for ( std::size_t i = 0; i < V.size(); )
1008 {
1009 Interval& J = V[ i ];
1010 // I=[a,b], J=[a',b'], a <= b, a' <= b'
1011 if ( I.second < J.first ) { break; } // b < a' : no further intersection
1012 if ( J.second < I.first ) { ++i; continue; } // b' < a : no further intersection
1013 // a' <= b and a <= b'
1014 // a ---------- b
1015 // a' ............... a'
1016 // b' ................. b'
1017
1018 // a' ..................... b' => a'..a-1 b+1 b'
1019 Interval K1( J.first, I.first - 1 );
1020 Interval K2( I.second + 1, J.second );
1021 bool K1_exist = K1.second >= K1.first;
1022 bool K2_exist = K2.second >= K2.first;
1023 if ( K1_exist && K2_exist )
1024 {
1025 V[ i ] = K2;
1026 V.insert( V.begin() + i, K1 );
1027 break; // no further intersection possible
1028 }
1029 else if ( K1_exist )
1030 {
1031 V[ i ] = K1; i++;
1032 }
1033 else if ( K2_exist )
1034 {
1035 V[ i ] = K2; break;
1036 }
1037 else
1038 {
1039 V.erase( V.begin() + i );
1040 }
1041 }
1042}
1043
1044//-----------------------------------------------------------------------------
1045template <typename TKSpace>
1046template <typename Predicate>
1047typename DGtal::DigitalConvexity<TKSpace>::PointRange
1048DGtal::DigitalConvexity<TKSpace>::
1049filter( const PointRange& E, const Predicate& Pred )
1050{
1051 PointRange Out;
1052 Out.reserve( E.size() );
1053 for ( auto&& p : E )
1054 if ( Pred( p ) ) Out.push_back( p );
1055 return Out;
1056}
1057
1058///////////////////////////////////////////////////////////////////////////////
1059// Interface - public :
1060
1061/**
1062 * Writes/Displays the object on an output stream.
1063 * @param out the output stream where the object is written.
1064 */
1065template <typename TKSpace>
1066inline
1067void
1068DGtal::DigitalConvexity<TKSpace>::selfDisplay ( std::ostream & out ) const
1069{
1070 out << "[DigitalConvexity]";
1071}
1072
1073/**
1074 * Checks the validity/consistency of the object.
1075 * @return 'true' if the object is valid, 'false' otherwise.
1076 */
1077template <typename TKSpace>
1078inline
1079bool
1080DGtal::DigitalConvexity<TKSpace>::isValid() const
1081{
1082 return true;
1083}
1084
1085
1086///////////////////////////////////////////////////////////////////////////////
1087// Implementation of inline functions //
1088
1089//-----------------------------------------------------------------------------
1090template <typename TKSpace>
1091inline
1092std::ostream&
1093DGtal::operator<< ( std::ostream & out,
1094 const DigitalConvexity<TKSpace> & object )
1095{
1096 object.selfDisplay( out );
1097 return out;
1098}
1099
1100// //
1101///////////////////////////////////////////////////////////////////////////////
Generated on Sun Nov 27 2022 15:15:41 for DGtal by doxygen 1.9.5