DGtal 2.1.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 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 ) {
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( i <= k );
253 ASSERT( k <= KSpace::dimension );
254 CellGeometry cgeom( myK, i, k, false );
255 cgeom.addCellsTouchingPoints( itB, itE );
256 return cgeom;
257}
258
259//-----------------------------------------------------------------------------
260template <typename TKSpace>
261typename DGtal::DigitalConvexity<TKSpace>::CellGeometry
262DGtal::DigitalConvexity<TKSpace>::
263makeCellCover( const LatticePolytope& P,
264 Dimension i, Dimension k ) const
265{
266 ASSERT( i <= k );
267 ASSERT( k <= KSpace::dimension );
268 CellGeometry cgeom( myK, i, k, false );
269 cgeom.addCellsTouchingPolytope( P );
270 return cgeom;
271}
272
273//-----------------------------------------------------------------------------
274template <typename TKSpace>
275typename DGtal::DigitalConvexity<TKSpace>::CellGeometry
276DGtal::DigitalConvexity<TKSpace>::
277makeCellCover( const RationalPolytope& P,
278 Dimension i, Dimension k ) const
279{
280 ASSERT( i <= k );
281 ASSERT( k <= KSpace::dimension );
282 CellGeometry cgeom( myK, i, k, false );
283 cgeom.addCellsTouchingPolytope( P );
284 return cgeom;
285}
286
287//-----------------------------------------------------------------------------
288template <typename TKSpace>
289typename DGtal::DigitalConvexity<TKSpace>::LatticePolytope
290DGtal::DigitalConvexity<TKSpace>::
291makePolytope( const PointRange& X, bool make_minkowski_summable ) const
292{
293 if ( mySafe )
294 {
295 typedef typename detail::ConvexityHelperInternalInteger< Integer, true >::Type
296 InternalInteger;
297 return ConvexityHelper< dimension, Integer, InternalInteger >::
298 computeLatticePolytope( X, false, make_minkowski_summable );
299 }
300 else
301 {
302 typedef typename detail::ConvexityHelperInternalInteger< Integer, false >::Type
303 InternalInteger;
304 return ConvexityHelper< dimension, Integer, InternalInteger >::
305 computeLatticePolytope( X, false, make_minkowski_summable );
306 }
307}
308
309//-----------------------------------------------------------------------------
310template <typename TKSpace>
311typename DGtal::DigitalConvexity<TKSpace>::PointRange
312DGtal::DigitalConvexity<TKSpace>::
313U( Dimension i, const PointRange& X ) const
314{
315 PointRange Y( X );
316 PointRange Z;
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 ) );
322 return Z;
323}
324
325//-----------------------------------------------------------------------------
326template <typename TKSpace>
327bool
328DGtal::DigitalConvexity<TKSpace>::
329is0Convex( const PointRange& X ) const
330{
331 if ( X.empty() ) return true;
332 const auto P = makePolytope( X );
333 return P.count() == (DGtal::BoundedLatticePolytope<DGtal::SpaceND<3>>::Integer)X.size();
334}
335
336//-----------------------------------------------------------------------------
337template <typename TKSpace>
338bool
339DGtal::DigitalConvexity<TKSpace>::
340isFullyConvex( const PointRange& Z, bool convex0 ) const
341{
342 if ( Z.size() <= 1 ) return true;
343 ASSERT( dimension <= 64 );
344 typedef DGtal::int64_t Direction;
345 typedef std::vector< Direction > Directions;
346 std::array< Directions, dimension > C;
347 const bool cvx0 = convex0 ? true : is0Convex( Z );
348 if ( ! cvx0 ) return false;
349 C[ 0 ].push_back( (Direction) 0 );
350 std::map< Direction, PointRange > X;
351 X[ 0 ] = Z;
352 std::sort( X[ 0 ].begin(), X[ 0 ].end() );
353 for ( Dimension k = 1; k < dimension; k++ )
354 {
355 for ( const auto beta : C[ k-1 ] )
356 {
357 for ( Dimension j = 0; j < dimension; j++ )
358 {
359 const Direction dir_j = Direction(1) << j;
360 if ( beta < dir_j )
361 {
362 const Direction alpha = beta | dir_j;
363 C[ k ].push_back( alpha );
364 X[ alpha ] = U( j, X[ beta ] );
365 if ( ! is0Convex( X[ alpha ] ) ) return false;
366 }
367 }
368 }
369 }
370 return true;
371}
372
373//-----------------------------------------------------------------------------
374template <typename TKSpace>
375bool
376DGtal::DigitalConvexity<TKSpace>::
377isFullyConvexFast( const PointRange& Z ) const
378{
379 if ( Z.size() <= 1 ) return true;
380 LatticeSet C_Z( Z.cbegin(), Z.cend(), 0 );
381 const auto nb_cells = C_Z.starOfPoints().size();
382 const auto s = sizeStarCvxH( Z );
383 return s == (Integer)nb_cells;
384}
385
386//-----------------------------------------------------------------------------
387template <typename TKSpace>
388typename DGtal::DigitalConvexity<TKSpace>::PointRange
389DGtal::DigitalConvexity<TKSpace>::
390ExtrCvxH( const PointRange& X ) const
391{
392 if ( mySafe )
393 {
394 typedef typename detail::ConvexityHelperInternalInteger< Integer, true >::Type
395 InternalInteger;
396 return ConvexityHelper< dimension, Integer, InternalInteger >::
397 computeLatticePolytopeVertices( X, false, false );
398 }
399 else
400 {
401 typedef typename detail::ConvexityHelperInternalInteger< Integer, false >::Type
402 InternalInteger;
403 return ConvexityHelper< dimension, Integer, InternalInteger >::
404 computeLatticePolytopeVertices( X, false, false );
405 }
406}
407
408//-----------------------------------------------------------------------------
409template <typename TKSpace>
410typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
411DGtal::DigitalConvexity<TKSpace>::
412StarCvxH( const PointRange& X, Dimension axis ) const
413{
414 PointRange cells;
415 // Computes Minkowski sum of Z with hypercube
416 PointRange Z = U( 0, X );
417 for ( Dimension k = 1; k < dimension; k++ )
418 Z = U( k, Z );
419 // Builds polytope
420 const auto P = makePolytope( Z );
421 // Extracts lattice points within polytope
422 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
423 Counter C( P );
424 const Dimension a = axis >= dimension ? C.longestAxis() : axis;
425 auto cellP = C.getLatticeCells( a );
426 return LatticeSet( cellP, a );
427}
428
429//-----------------------------------------------------------------------------
430template <typename TKSpace>
431typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
432DGtal::DigitalConvexity<TKSpace>::
433StarCvxH( const Point& a, const Point& b, const Point& c,
434 Dimension axis ) const
435{
436 LatticeSet LS;
437 if ( mySafe )
438 {
439 using InternalInteger
440 = typename detail::ConvexityHelperInternalInteger< Integer, true >::Type;
441 using Helper = ConvexityHelper< dimension, Integer, InternalInteger >;
442 using UnitSegment = typename Helper::LatticePolytope::UnitSegment;
443 auto P = Helper::compute3DTriangle( a, b, c, true );
444 if ( ! P.isValid() ) return LS;
445 P += UnitSegment( 0 );
446 P += UnitSegment( 1 );
447 P += UnitSegment( 2 );
448 // Extracts lattice points within polytope
449 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
450 Counter C( P );
451 if ( axis >= dimension ) axis = C.longestAxis();
452 const auto cellP = C.getLatticeCells( axis );
453 return LatticeSet( cellP, axis );
454 }
455 else
456 {
457 using InternalInteger
458 = typename detail::ConvexityHelperInternalInteger< Integer, false >::Type;
459 using Helper = ConvexityHelper< dimension, Integer, InternalInteger >;
460 using UnitSegment = typename Helper::LatticePolytope::UnitSegment;
461 auto P = Helper::compute3DTriangle( a, b, c, true );
462 if ( ! P.isValid() ) return LS;
463 P += UnitSegment( 0 );
464 P += UnitSegment( 1 );
465 P += UnitSegment( 2 );
466 // Extracts lattice points within polytope
467 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
468 Counter C( P );
469 if ( axis >= dimension ) axis = C.longestAxis();
470 const auto cellP = C.getLatticeCells( axis );
471 return LatticeSet( cellP, axis );
472 }
473}
474
475//-----------------------------------------------------------------------------
476template <typename TKSpace>
477typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
478DGtal::DigitalConvexity<TKSpace>::
479StarOpenTriangle( const Point& a, const Point& b, const Point& c,
480 Dimension axis ) const
481{
482 LatticeSet LS;
483 if ( mySafe )
484 {
485 using InternalInteger
486 = typename detail::ConvexityHelperInternalInteger< Integer, true >::Type;
487 using Helper = ConvexityHelper< dimension, Integer, InternalInteger >;
488 using UnitSegment = typename Helper::LatticePolytope::UnitSegment;
489 auto P = Helper::compute3DOpenTriangle( a, b, c, true );
490 if ( ! P.isValid() ) return LS;
491 P += UnitSegment( 0 );
492 P += UnitSegment( 1 );
493 P += UnitSegment( 2 );
494 // Extracts lattice points within polytope
495 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
496 Counter C( P );
497 if ( axis >= dimension ) axis = C.longestAxis();
498 const auto cellP = C.getLatticeCells( axis );
499 return LatticeSet( cellP, axis );
500 }
501 else
502 {
503 using InternalInteger
504 = typename detail::ConvexityHelperInternalInteger< Integer, false >::Type;
505 using Helper = ConvexityHelper< dimension, Integer, InternalInteger >;
506 using UnitSegment = typename Helper::LatticePolytope::UnitSegment;
507 auto P = Helper::compute3DOpenTriangle( a, b, c, true );
508 if ( ! P.isValid() ) return LS;
509 P += UnitSegment( 0 );
510 P += UnitSegment( 1 );
511 P += UnitSegment( 2 );
512 // Extracts lattice points within polytope
513 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
514 Counter C( P );
515 if ( axis >= dimension ) axis = C.longestAxis();
516 const auto cellP = C.getLatticeCells( axis );
517 return LatticeSet( cellP, axis );
518 }
519}
520
521//-----------------------------------------------------------------------------
522template <typename TKSpace>
523typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
524DGtal::DigitalConvexity<TKSpace>::
525Star( const PointRange& X, const Dimension axis ) const
526{
527 const Dimension a = axis >= dimension ? 0 : axis;
528 LatticeSet L( X.cbegin(), X.cend(), a );
529 return L.starOfPoints();
530}
531//-----------------------------------------------------------------------------
532template <typename TKSpace>
533typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
534DGtal::DigitalConvexity<TKSpace>::
535StarCells( const PointRange& C, const Dimension axis ) const
536{
537 const Dimension a = axis >= dimension ? 0 : axis;
538 LatticeSet L( C.cbegin(), C.cend(), a );
539 return L.starOfCells();
540}
541
542//-----------------------------------------------------------------------------
543template <typename TKSpace>
544template <typename TPolytope>
545typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
546DGtal::DigitalConvexity<TKSpace>::
547CoverPolytope( const TPolytope& P, Dimension axis ) const
548{
549 using StrictUnitSegment = typename Polytope::StrictUnitSegment;
550 LatticeSet LS;
551 if ( ! P.isValid() ) return LS;
552 if ( ! P.canBeSummed() ) {
553 trace.error() << "[DigitalConvexity::CoverPolytope] Polytope should be Minkowski summabel." << std::endl;
554 return LS;
555 }
556 Counter C( P );
557 if ( axis >= dimension ) axis = C.longestAxis();
558 auto V = LatticeSet( C.getLatticeSet( axis ), axis ).toPointRange(); // << 0-cells
559 auto P0 = P + StrictUnitSegment( 0 );
560 auto P1 = P + StrictUnitSegment( 1 );
561 auto P2 = P + StrictUnitSegment( 2 );
562 Counter C0( P0 );
563 Counter C1( P1 );
564 Counter C2( P2 );
565 auto V0 = LatticeSet( C0.getLatticeSet( axis ), axis ).toPointRange(); // << 1-cells
566 auto V1 = LatticeSet( C1.getLatticeSet( axis ), axis ).toPointRange(); // << 1-cells
567 auto V2 = LatticeSet( C2.getLatticeSet( axis ), axis ).toPointRange(); // << 1-cells
568 auto P01 = P0 + StrictUnitSegment( 1 );
569 auto P02 = P0 + StrictUnitSegment( 2 );
570 auto P12 = P1 + StrictUnitSegment( 2 );
571 Counter C01( P01 );
572 Counter C02( P02 );
573 Counter C12( P12 );
574 auto V01 = LatticeSet( C01.getLatticeSet( axis ), axis ).toPointRange(); // << 2-cells
575 auto V02 = LatticeSet( C02.getLatticeSet( axis ), axis ).toPointRange(); // << 2-cells
576 auto V12 = LatticeSet( C12.getLatticeSet( axis ), axis ).toPointRange(); // << 2-cells
577 auto P012 = P01 + StrictUnitSegment( 2 );
578 Counter C012( P012 );
579 auto V012 = LatticeSet( C012.getLatticeSet( axis ), axis ).toPointRange(); // << 3-cells
580 // Convert points to Khalimsky coordinates and insert them in LatticeSet.
581 LS = LatticeSet( axis );
582 for ( const auto& p : V ) LS.insert( Point( 2*p[0] , 2*p[1] , 2*p[2] ) );
583 for ( const auto& p : V0 ) LS.insert( Point( 2*p[0]-1, 2*p[1] , 2*p[2] ) );
584 for ( const auto& p : V1 ) LS.insert( Point( 2*p[0] , 2*p[1]-1, 2*p[2] ) );
585 for ( const auto& p : V2 ) LS.insert( Point( 2*p[0] , 2*p[1] , 2*p[2]-1 ) );
586 for ( const auto& p : V01 ) LS.insert( Point( 2*p[0]-1, 2*p[1]-1, 2*p[2] ) );
587 for ( const auto& p : V02 ) LS.insert( Point( 2*p[0]-1, 2*p[1] , 2*p[2]-1 ) );
588 for ( const auto& p : V12 ) LS.insert( Point( 2*p[0] , 2*p[1]-1, 2*p[2]-1 ) );
589 for ( const auto& p : V012 ) LS.insert( Point( 2*p[0]-1, 2*p[1]-1, 2*p[2]-1 ) );
590 return LS;
591}
592
593//-----------------------------------------------------------------------------
594template <typename TKSpace>
595typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
596DGtal::DigitalConvexity<TKSpace>::
597CoverCvxH( const Point& a, const Point& b, const Point& c,
598 Dimension axis ) const
599{
600 LatticeSet LS;
601 if ( mySafe )
602 {
603 using InternalInteger
604 = typename detail::ConvexityHelperInternalInteger< Integer, true >::Type;
605 using Helper = ConvexityHelper< dimension, Integer, InternalInteger >;
606 // using StrictUnitSegment = typename Helper::LatticePolytope::StrictUnitSegment;
607 auto P = Helper::compute3DTriangle( a, b, c, true );
608 return CoverPolytope( P, axis );
609 }
610 else
611 {
612 using InternalInteger
613 = typename detail::ConvexityHelperInternalInteger< Integer, false >::Type;
614 using Helper = ConvexityHelper< dimension, Integer, InternalInteger >;
615 // using StrictUnitSegment = typename Helper::LatticePolytope::StrictUnitSegment;
616 auto P = Helper::compute3DTriangle( a, b, c, true );
617 return CoverPolytope( P, axis );
618 }
619}
620
621//-----------------------------------------------------------------------------
622template <typename TKSpace>
623typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
624DGtal::DigitalConvexity<TKSpace>::
625CoverCvxH( const Point& a, const Point& b,
626 Dimension axis ) const
627{
628 LatticeSet LS;
629 if ( mySafe )
630 {
631 using InternalInteger
632 = typename detail::ConvexityHelperInternalInteger< Integer, true >::Type;
633 using Helper = ConvexityHelper< dimension, Integer, InternalInteger >;
634 auto P = Helper::computeSegment( a, b );
635 return CoverPolytope( P, axis );
636 }
637 else
638 {
639 using InternalInteger
640 = typename detail::ConvexityHelperInternalInteger< Integer, false >::Type;
641 using Helper = ConvexityHelper< dimension, Integer, InternalInteger >;
642 auto P = Helper::computeSegment( a, b );
643 return CoverPolytope( P, axis );
644 }
645}
646
647
648//-----------------------------------------------------------------------------
649template <typename TKSpace>
650typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
651DGtal::DigitalConvexity<TKSpace>::
652toLatticeSet( const PointRange& X, const Dimension axis ) const
653{
654 const Dimension a = axis >= dimension ? 0 : axis;
655 return LatticeSet( X.cbegin(), X.cend(), a );
656}
657
658//-----------------------------------------------------------------------------
659template <typename TKSpace>
660typename DGtal::DigitalConvexity<TKSpace>::PointRange
661DGtal::DigitalConvexity<TKSpace>::
662toPointRange( const LatticeSet& L ) const
663{
664 return L.toPointRange();
665}
666
667//-----------------------------------------------------------------------------
668template <typename TKSpace>
669typename DGtal::DigitalConvexity<TKSpace>::Integer
670DGtal::DigitalConvexity<TKSpace>::
671sizeStarCvxH( const PointRange& X ) const
672{
673 PointRange cells;
674 // Computes Minkowski sum of Z with hypercube
675 PointRange Z = U( 0, X );
676 for ( Dimension k = 1; k < dimension; k++ )
677 Z = U( k, Z );
678 // Builds polytope
679 const auto P = makePolytope( Z );
680 // Extracts lattice points within polytope
681 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
682 Counter C( P );
683 const Dimension a = C.longestAxis();
684 auto cellP = C.getLatticeCells( a );
685
686 // Counts the number of cells
687 Integer nb = 0;
688 for ( const auto& value : cellP )
689 {
690 Point p = value.first;
691 Interval I = value.second;
692 nb += I.second - I.first + 1;
693 }
694 return nb;
695}
696
697//-----------------------------------------------------------------------------
698template <typename TKSpace>
699typename DGtal::DigitalConvexity<TKSpace>::PointRange
700DGtal::DigitalConvexity<TKSpace>::
701Extr( const PointRange& C ) const
702{
703 // JOL: using std::set< Point > or std::unordered_set< Point > is slightly slower.
704 // We prefer to use vector for easier vectorization.
705 std::vector<Point> E;
706 E.reserve( 2*C.size() );
707 for ( auto&& kp : C )
708 {
709 auto c = myK.uCell( kp );
710 if ( myK.uDim( c ) == 0 )
711 E.push_back( myK.uCoords( c ) );
712 else
713 {
714 auto faces = myK.uFaces( c );
715 for ( auto&& f : faces )
716 if ( myK.uDim( f ) == 0 )
717 E.push_back( myK.uCoords( f ) );
718 }
719 }
720 std::sort( E.begin(), E.end() );
721 auto last = std::unique( E.begin(), E.end() );
722 E.erase( last, E.end() );
723 return E;
724}
725//-----------------------------------------------------------------------------
726template <typename TKSpace>
727typename DGtal::DigitalConvexity<TKSpace>::PointRange
728DGtal::DigitalConvexity<TKSpace>::
729Extr( const LatticeSet& C ) const
730{
731 return C.extremaOfCells();
732}
733//-----------------------------------------------------------------------------
734template <typename TKSpace>
735typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
736DGtal::DigitalConvexity<TKSpace>::
737Skel( const LatticeSet& C ) const
738{
739 return C.skeletonOfCells();
740}
741//-----------------------------------------------------------------------------
742template <typename TKSpace>
743typename DGtal::DigitalConvexity<TKSpace>::PointRange
744DGtal::DigitalConvexity<TKSpace>::
745ExtrSkel( const LatticeSet& C ) const
746{
747 return C.skeletonOfCells().extremaOfCells();
748}
749
750//-----------------------------------------------------------------------------
751template <typename TKSpace>
752typename DGtal::DigitalConvexity<TKSpace>::PointRange
753DGtal::DigitalConvexity<TKSpace>::
754FC_direct( const PointRange& Z ) const
755{
756 typedef typename LatticePolytope::Domain Domain;
757 PointRange cells;
758 // Computes Minkowski sum of Z with hypercube
759 PointRange X( Z );
760 for ( Dimension k = 0; k < dimension; k++ )
761 X = U( k, X );
762 // Builds polytope
763 const auto P = makePolytope( X );
764 // Extracts lattice points within polytope
765 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
766 Counter C( P );
767 const Dimension a = C.longestAxis();
768 Point lo = C.lowerBound();
769 Point hi = C.upperBound();
770 hi[ a ] = lo[ a ];
771 const Domain projD( lo, hi ); //< the projected domain of the polytope.
772 const Point One = Point::diagonal( 1 );
773 //const Size size = projD.size(); //not USED
774 std::unordered_map< Point, Interval > cellP;
775 Point q;
776 for ( auto&& p : projD )
777 {
778 q = 2*p - One;
779 const auto I = C.intersectionIntervalAlongAxis( p, a );
780 const auto n = I.second - I.first;
781 if ( n != 0 )
782 {
783 // Now the second bound is included
784 cellP[ q ] = Interval( 2 * I.first - 1, 2 * I.second - 3 );
785 }
786 }
787 // It remains to compute all the k-cells, 0 <= k < d, intersected by Cvxh( Z )
788 for ( Dimension k = 0; k < dimension; k++ )
789 {
790 if ( k == a ) continue;
791 std::vector< Point > q_computed;
792 std::vector< Interval > I_computed;
793 for ( const auto& value : cellP )
794 {
795 Point p = value.first;
796 Interval I = value.second;
797 Point r = p; r[ k ] += 2;
798 const auto it = cellP.find( r );
799 if ( it == cellP.end() ) continue; // neighbor is empty
800 // Otherwise compute common part.
801 Interval J = it->second;
802 auto f = std::max( I.first, J.first );
803 auto s = std::min( I.second, J.second );
804 if ( f <= s )
805 {
806 Point qq = p; qq[ k ] += 1;
807 q_computed.push_back( qq );
808 I_computed.push_back( Interval( f, s ) );
809 }
810 }
811 // Add new columns to map Point -> column
812 for ( size_t i = 0; i < q_computed.size(); ++i )
813 {
814 cellP[ q_computed[ i ] ] = I_computed[ i ];
815 }
816 }
817 // The built complex is open.
818 // Check it and fill skelP
819 std::unordered_map< Point, std::vector< Interval > > skelP;
820 for ( const auto& value : cellP )
821 {
822 Point p = value.first;
823 Interval I = value.second;
824 auto n = I.second - I.first + 1;
825 if ( n % 2 == 0 )
826 trace.error() << "Weird thickness step 1="
827 << n << std::endl;
828 std::vector< Interval > V( 1, I );
829 auto result = skelP.insert( std::make_pair( p, V ) );
830 (void)result;//unused var;
831 }
832
833 // std::cout << "Extract skel" << std::endl;
834 // Now extracting implicitly its Skel
835 for ( const auto& value : cellP )
836 {
837 const Point& p = value.first;
838 const auto& I = value.second;
839 for ( Dimension k = 0; k < dimension; k++ )
840 {
841 if ( k == a ) continue;
842 if ( ( p[ k ] & 0x1 ) != 0 ) continue; // if open along axis continue
843 // if closed, check upper incident cells along direction k
844 Point qq = p; qq[ k ] -= 1;
845 Point r = p; r[ k ] += 1;
846 auto itq = skelP.find( qq );
847 if ( itq != skelP.end() )
848 {
849 auto& W = itq->second;
850 eraseInterval( I, W );
851 // if ( W.empty() ) skelP.erase( itq );
852 }
853 auto itr = skelP.find( r );
854 if ( itr != skelP.end() )
855 {
856 auto& W = itr->second;
857 eraseInterval( I, W );
858 // if ( W.empty() ) skelP.erase( itr );
859 }
860 }
861 }
862 // Extract skel along main axis
863 for ( const auto& value : cellP )
864 {
865 const Point& p = value.first;
866 auto I = value.second;
867 const auto itp = skelP.find( p );
868 if ( itp == skelP.end() ) continue;
869 if ( ( I.first & 0x1 ) != 0 ) I.first += 1;
870 if ( ( I.second & 0x1 ) != 0 ) I.second -= 1;
871 auto& W = itp->second;
872 for ( auto x = I.first; x <= I.second; x += 2 )
873 {
874 eraseInterval( Interval{ x-1,x-1} , W );
875 eraseInterval( Interval{ x+1,x+1} , W );
876 }
877 }
878 // Erase empty stacks
879 for ( auto it = skelP.begin(), itE = skelP.end(); it != itE; )
880 {
881 const auto& V = it->second;
882 if ( V.empty() )
883 {
884 auto it_erase = it;
885 ++it;
886 skelP.erase( it_erase );
887 }
888 else ++it;
889 }
890 // Skel is constructed, now compute its Extr.
891 PointRange O;
892 for ( const auto& value : skelP )
893 {
894 Point p = value.first;
895 auto W = value.second;
896 for ( auto&& I : W )
897 {
898 p[ a ] = I.first;
899 O.push_back( p );
900 }
901 }
902 return Extr( O );
903}
904//-----------------------------------------------------------------------------
905template <typename TKSpace>
906typename DGtal::DigitalConvexity<TKSpace>::PointRange
907DGtal::DigitalConvexity<TKSpace>::
908FC_LatticeSet( const PointRange& Z ) const
909{
910 const auto StarCvxZ = StarCvxH( Z );
911 const auto SkelStarCvxZ = StarCvxZ.skeletonOfCells().toPointRange();
912 return Extr( SkelStarCvxZ );
913}
914
915//-----------------------------------------------------------------------------
916template <typename TKSpace>
917typename DGtal::DigitalConvexity<TKSpace>::PointRange
918DGtal::DigitalConvexity<TKSpace>::
919FC( const PointRange& Z, EnvelopeAlgorithm algo ) const
920{
921 if ( algo == EnvelopeAlgorithm::DIRECT )
922 return FC_direct( Z );
923 else if ( algo == EnvelopeAlgorithm::LATTICE_SET )
924 return FC_LatticeSet( Z );
925 else
926 return Z;
927}
928
929//-----------------------------------------------------------------------------
930template <typename TKSpace>
931typename DGtal::DigitalConvexity<TKSpace>::PointRange
932DGtal::DigitalConvexity<TKSpace>::
933envelope( const PointRange& Z, EnvelopeAlgorithm algo ) const
934{
935 myDepthLastFCE = 0;
936 auto In = Z;
937 while (true) {
938 auto card_In = In.size();
939 In = FC( In, algo );
940 if ( In.size() == card_In ) return In;
941 myDepthLastFCE++;
942 }
943 trace.error() << "[DigitalConvexity::envelope] Should never pass here."
944 << std::endl;
945 return Z; // to avoid warnings.
946}
947
948//-----------------------------------------------------------------------------
949template <typename TKSpace>
950typename DGtal::DigitalConvexity<TKSpace>::PointRange
951DGtal::DigitalConvexity<TKSpace>::
952relativeEnvelope( const PointRange& Z, const PointRange& Y,
953 EnvelopeAlgorithm algo ) const
954{
955 myDepthLastFCE = 0;
956 auto In = Z;
957 while (true) {
958 PointRange Out;
959 std::set_intersection( In.cbegin(), In.cend(),
960 Y.cbegin(), Y.cend(),
961 std::back_inserter( Out ) );
962 In = FC( Out, algo );
963 if ( In.size() == Out.size() ) return Out;
964 myDepthLastFCE++;
965 }
966 return Z;
967}
968
969//-----------------------------------------------------------------------------
970template <typename TKSpace>
971template <typename Predicate>
972typename DGtal::DigitalConvexity<TKSpace>::PointRange
973DGtal::DigitalConvexity<TKSpace>::
974relativeEnvelope( const PointRange& Z, const Predicate& PredY,
975 EnvelopeAlgorithm algo ) const
976{
977 myDepthLastFCE = 0;
978 auto In = Z;
979 while (true) {
980 auto Out = filter( In, PredY );
981 In = FC( Out, algo );
982 if ( In.size() == Out.size() ) return In;
983 myDepthLastFCE++;
984 }
985 return Z;
986}
987
988//-----------------------------------------------------------------------------
989template <typename TKSpace>
990typename DGtal::DigitalConvexity<TKSpace>::Size
991DGtal::DigitalConvexity<TKSpace>::
992depthLastEnvelope() const
993{
994 return myDepthLastFCE;
995}
996
997//-----------------------------------------------------------------------------
998template <typename TKSpace>
999bool
1000DGtal::DigitalConvexity<TKSpace>::
1001isKConvex( const LatticePolytope& P, const Dimension k ) const
1002{
1003 if ( k == 0 ) return true;
1004 auto S = insidePoints( P );
1005 auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
1006 auto intersected_cells = makeCellCover( P, k, k );
1007 return intersected_cells.nbCells() == touched_cells.nbCells();
1008 // JOL: number should be enough
1009 // && intersected_cells.subset( touched_cells );
1010}
1011
1012//-----------------------------------------------------------------------------
1013template <typename TKSpace>
1014bool
1015DGtal::DigitalConvexity<TKSpace>::
1016isFullyConvex( const LatticePolytope& P ) const
1017{
1018 auto S = insidePoints( P );
1019 if ( S.size() <= 1 ) return true;
1020 for ( Dimension k = 1; k < KSpace::dimension; ++ k )
1021 {
1022 auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
1023 auto intersected_cells = makeCellCover( P, k, k );
1024 if ( ( intersected_cells.nbCells() != touched_cells.nbCells() ) )
1025 // JOL: number should be enough
1026 // || ( ! intersected_cells.subset( touched_cells ) ) )
1027 return false;
1028 }
1029 return true;
1030}
1031
1032//-----------------------------------------------------------------------------
1033template <typename TKSpace>
1034bool
1035DGtal::DigitalConvexity<TKSpace>::
1036isKSubconvex( const LatticePolytope& P, const CellGeometry& C, const Dimension k ) const
1037{
1038 auto intersected_cells = makeCellCover( P, k, k );
1039 return intersected_cells.subset( C );
1040}
1041//-----------------------------------------------------------------------------
1042template <typename TKSpace>
1043bool
1044DGtal::DigitalConvexity<TKSpace>::
1045isFullySubconvex( const LatticePolytope& P, const CellGeometry& C ) const
1046{
1047 auto intersected_cells = makeCellCover( P, C.minCellDim(), C.maxCellDim() );
1048 return intersected_cells.subset( C );
1049}
1050
1051//-----------------------------------------------------------------------------
1052template <typename TKSpace>
1053bool
1054DGtal::DigitalConvexity<TKSpace>::
1055isFullySubconvex( const LatticePolytope& P, const LatticeSet& StarX ) const
1056{
1057 LatticePolytope Q = P + typename LatticePolytope::UnitSegment( 0 );
1058 for ( Dimension k = 1; k < dimension; k++ )
1059 Q = Q + typename LatticePolytope::UnitSegment( k );
1060 Counter C( Q );
1061 const auto cells = C.getLatticeCells( StarX.axis() );
1062 LatticeSet StarP( cells, StarX.axis() );
1063 return StarX.includes( StarP );
1064}
1065
1066//-----------------------------------------------------------------------------
1067template <typename TKSpace>
1068bool
1069DGtal::DigitalConvexity<TKSpace>::
1070isFullySubconvex( const PointRange& Y, const LatticeSet& StarX ) const
1071{
1072 const auto SCY = StarCvxH( Y, StarX.axis() );
1073 return StarX.includes( SCY );
1074}
1075
1076//-----------------------------------------------------------------------------
1077template <typename TKSpace>
1078bool
1079DGtal::DigitalConvexity<TKSpace>::
1080isFullySubconvex( const Point& a, const Point& b, const Point& c,
1081 const LatticeSet& StarX ) const
1082{
1083 ASSERT( dimension == 3 );
1084 const auto SCabc = StarCvxH( a, b, c, StarX.axis() );
1085 return StarX.includes( SCabc );
1086}
1087//-----------------------------------------------------------------------------
1088template <typename TKSpace>
1089bool
1090DGtal::DigitalConvexity<TKSpace>::
1091isFullyCovered( const Point& a, const Point& b, const Point& c,
1092 const LatticeSet& cells ) const
1093{
1094 ASSERT( dimension == 3 );
1095 const auto SCabc = CoverCvxH( a, b, c, cells.axis() );
1096 return cells.includes( SCabc );
1097}
1098//-----------------------------------------------------------------------------
1099template <typename TKSpace>
1100bool
1101DGtal::DigitalConvexity<TKSpace>::
1102isFullyCovered( const Point& a, const Point& b,
1103 const LatticeSet& cells ) const
1104{
1105 ASSERT( dimension == 3 );
1106 LatticeSet L_ab( cells.axis() );
1107 const auto V = b - a;
1108 L_ab.insert( 2*a );
1109 for ( Dimension k = 0; k < dimension; k++ )
1110 {
1111 const Integer n = ( V[ k ] >= 0 ) ? V[ k ] : -V[ k ];
1112 const Integer d = ( V[ k ] >= 0 ) ? 1 : -1;
1113 if ( n == 0 ) continue;
1114 Point kc;
1115 for ( Integer i = 1; i < n; i++ )
1116 {
1117 for ( Dimension j = 0; j < dimension; j++ )
1118 {
1119 if ( j == k ) kc[ k ] = 2 * ( a[ k ] + d * i );
1120 else
1121 {
1122 const auto v = V[ j ];
1123 const auto q = ( v * i ) / n;
1124 const auto r = ( v * i ) % n; // might be negative
1125 kc[ j ] = 2 * ( a[ j ] + q );
1126 if ( r < 0 ) kc[ j ] -= 1;
1127 else if ( r > 0 ) kc[ j ] += 1;
1128 }
1129 }
1130 L_ab.insert( kc );
1131 }
1132 }
1133 if ( a != b ) L_ab.insert( 2*b );
1134 // LatticeSet Star_ab = L_ab.starOfCells();
1135 return cells.includes( L_ab );
1136
1137 // ASSERT( dimension == 3 );
1138 // const auto SCab = CoverCvxH( a, b, cells.axis() );
1139 // return cells.includes( SCab );
1140}
1141
1142//-----------------------------------------------------------------------------
1143template <typename TKSpace>
1144bool
1145DGtal::DigitalConvexity<TKSpace>::
1146isOpenTriangleFullySubconvex( const Point& a, const Point& b, const Point& c,
1147 const LatticeSet& StarX ) const
1148{
1149 ASSERT( dimension == 3 );
1150 const auto SCabc = StarOpenTriangle( a, b, c, StarX.axis() );
1151 return StarX.includes( SCabc );
1152}
1153
1154//-----------------------------------------------------------------------------
1155template <typename TKSpace>
1156bool
1157DGtal::DigitalConvexity<TKSpace>::
1158isKSubconvex( const Point& a, const Point& b,
1159 const CellGeometry& C, const Dimension k ) const
1160{
1161 CellGeometry cgeom( myK, k, k, false );
1162 cgeom.addCellsTouchingSegment( a, b );
1163 return cgeom.subset( C );
1164}
1165
1166//-----------------------------------------------------------------------------
1167template <typename TKSpace>
1168bool
1169DGtal::DigitalConvexity<TKSpace>::
1170isFullySubconvex( const Point& a, const Point& b,
1171 const CellGeometry& C ) const
1172{
1173 CellGeometry cgeom( myK, C.minCellDim(), C.maxCellDim(), false );
1174 cgeom.addCellsTouchingSegment( a, b );
1175 return cgeom.subset( C );
1176}
1177
1178//-----------------------------------------------------------------------------
1179template <typename TKSpace>
1180bool
1181DGtal::DigitalConvexity<TKSpace>::
1182isFullySubconvex( const Point& a, const Point& b,
1183 const LatticeSet& StarX ) const
1184{
1185 LatticeSet L_ab( StarX.axis() );
1186 const auto V = b - a;
1187 L_ab.insert( 2*a );
1188 for ( Dimension k = 0; k < dimension; k++ )
1189 {
1190 const Integer n = ( V[ k ] >= 0 ) ? V[ k ] : -V[ k ];
1191 const Integer d = ( V[ k ] >= 0 ) ? 1 : -1;
1192 if ( n == 0 ) continue;
1193 Point kc;
1194 for ( Integer i = 1; i < n; i++ )
1195 {
1196 for ( Dimension j = 0; j < dimension; j++ )
1197 {
1198 if ( j == k ) kc[ k ] = 2 * ( a[ k ] + d * i );
1199 else
1200 {
1201 const auto v = V[ j ];
1202 const auto q = ( v * i ) / n;
1203 const auto r = ( v * i ) % n; // might be negative
1204 kc[ j ] = 2 * ( a[ j ] + q );
1205 if ( r < 0 ) kc[ j ] -= 1;
1206 else if ( r > 0 ) kc[ j ] += 1;
1207 }
1208 }
1209 L_ab.insert( kc );
1210 }
1211 }
1212 if ( a != b ) L_ab.insert( 2*b );
1213 LatticeSet Star_ab = L_ab.starOfCells();
1214 return StarX.includes( Star_ab );
1215}
1216
1217//-----------------------------------------------------------------------------
1218template <typename TKSpace>
1219bool
1220DGtal::DigitalConvexity<TKSpace>::
1221isKConvex( const RationalPolytope& P, const Dimension k ) const
1222{
1223 if ( k == 0 ) return true;
1224 auto Q = P.denominator() * P;
1225 auto S = insidePoints( Q );
1226 auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
1227 auto intersected_cells = makeCellCover( Q, k, k );
1228 return intersected_cells.nbCells() == touched_cells.nbCells()
1229 && intersected_cells.subset( touched_cells );
1230}
1231
1232//-----------------------------------------------------------------------------
1233template <typename TKSpace>
1234bool
1235DGtal::DigitalConvexity<TKSpace>::
1236isFullyConvex( const RationalPolytope& P ) const
1237{
1238 auto Q = P.denominator() * P;
1239 auto S = insidePoints( Q );
1240 if ( S.size() <= 1 ) return true;
1241 for ( Dimension k = 1; k < KSpace::dimension; ++ k )
1242 {
1243 auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
1244 auto intersected_cells = makeCellCover( Q, k, k );
1245 if ( ( intersected_cells.nbCells() != touched_cells.nbCells() )
1246 || ( ! intersected_cells.subset( touched_cells ) ) )
1247 return false;
1248 }
1249 return true;
1250}
1251
1252//-----------------------------------------------------------------------------
1253template <typename TKSpace>
1254bool
1255DGtal::DigitalConvexity<TKSpace>::
1256isKSubconvex( const RationalPolytope& P, const CellGeometry& C,
1257 const Dimension k ) const
1258{
1259 auto intersected_cells = makeCellCover( P, k, k );
1260 return intersected_cells.subset( C );
1261}
1262//-----------------------------------------------------------------------------
1263template <typename TKSpace>
1264bool
1265DGtal::DigitalConvexity<TKSpace>::
1266isFullySubconvex( const RationalPolytope& P, const CellGeometry& C ) const
1267{
1268 auto intersected_cells = makeCellCover( P, C.minCellDim(), C.maxCellDim() );
1269 return intersected_cells.subset( C );
1270}
1271
1272//-----------------------------------------------------------------------------
1273template <typename TKSpace>
1274void
1275DGtal::DigitalConvexity<TKSpace>::
1276eraseInterval( Interval I, std::vector< Interval > & V )
1277{
1278 for ( std::size_t i = 0; i < V.size(); )
1279 {
1280 Interval& J = V[ i ];
1281 // I=[a,b], J=[a',b'], a <= b, a' <= b'
1282 if ( I.second < J.first ) { break; } // b < a' : no further intersection
1283 if ( J.second < I.first ) { ++i; continue; } // b' < a : no further intersection
1284 // a' <= b and a <= b'
1285 // a ---------- b
1286 // a' ............... a'
1287 // b' ................. b'
1288
1289 // a' ..................... b' => a'..a-1 b+1 b'
1290 Interval K1( J.first, I.first - 1 );
1291 Interval K2( I.second + 1, J.second );
1292 bool K1_exist = K1.second >= K1.first;
1293 bool K2_exist = K2.second >= K2.first;
1294 if ( K1_exist && K2_exist )
1295 {
1296 V[ i ] = K2;
1297 V.insert( V.begin() + i, K1 );
1298 break; // no further intersection possible
1299 }
1300 else if ( K1_exist )
1301 {
1302 V[ i ] = K1; i++;
1303 }
1304 else if ( K2_exist )
1305 {
1306 V[ i ] = K2; break;
1307 }
1308 else
1309 {
1310 V.erase( V.begin() + i );
1311 }
1312 }
1313}
1314
1315//-----------------------------------------------------------------------------
1316template <typename TKSpace>
1317template <typename Predicate>
1318typename DGtal::DigitalConvexity<TKSpace>::PointRange
1319DGtal::DigitalConvexity<TKSpace>::
1320filter( const PointRange& E, const Predicate& Pred )
1321{
1322 PointRange Out;
1323 Out.reserve( E.size() );
1324 for ( auto&& p : E )
1325 if ( Pred( p ) ) Out.push_back( p );
1326 return Out;
1327}
1328
1329///////////////////////////////////////////////////////////////////////////////
1330// Interface - public :
1331
1332/**
1333 * Writes/Displays the object on an output stream.
1334 * @param out the output stream where the object is written.
1335 */
1336template <typename TKSpace>
1337inline
1338void
1339DGtal::DigitalConvexity<TKSpace>::selfDisplay ( std::ostream & out ) const
1340{
1341 out << "[DigitalConvexity]";
1342}
1343
1344/**
1345 * Checks the validity/consistency of the object.
1346 * @return 'true' if the object is valid, 'false' otherwise.
1347 */
1348template <typename TKSpace>
1349inline
1350bool
1351DGtal::DigitalConvexity<TKSpace>::isValid() const
1352{
1353 return true;
1354}
1355
1356
1357///////////////////////////////////////////////////////////////////////////////
1358// Implementation of inline functions //
1359
1360//-----------------------------------------------------------------------------
1361template <typename TKSpace>
1362inline
1363std::ostream&
1364DGtal::operator<< ( std::ostream & out,
1365 const DigitalConvexity<TKSpace> & object )
1366{
1367 object.selfDisplay( out );
1368 return out;
1369}
1370
1371// //
1372///////////////////////////////////////////////////////////////////////////////