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