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