DGtal  1.5.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 StarCvxH( 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 //-----------------------------------------------------------------------------
474 template <typename TKSpace>
475 typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
476 DGtal::DigitalConvexity<TKSpace>::
477 Star( const PointRange& X, const Dimension axis ) const
478 {
479  const Dimension a = axis >= dimension ? 0 : axis;
480  LatticeSet L( X.cbegin(), X.cend(), a );
481  return L.starOfPoints();
482 }
483 //-----------------------------------------------------------------------------
484 template <typename TKSpace>
485 typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
486 DGtal::DigitalConvexity<TKSpace>::
487 StarCells( const PointRange& C, const Dimension axis ) const
488 {
489  const Dimension a = axis >= dimension ? 0 : axis;
490  LatticeSet L( C.cbegin(), C.cend(), a );
491  return L.starOfCells();
492 }
493 
494 //-----------------------------------------------------------------------------
495 template <typename TKSpace>
496 typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
497 DGtal::DigitalConvexity<TKSpace>::
498 toLatticeSet( const PointRange& X, const Dimension axis ) const
499 {
500  const Dimension a = axis >= dimension ? 0 : axis;
501  return LatticeSet( X.cbegin(), X.cend(), a );
502 }
503 
504 //-----------------------------------------------------------------------------
505 template <typename TKSpace>
506 typename DGtal::DigitalConvexity<TKSpace>::PointRange
507 DGtal::DigitalConvexity<TKSpace>::
508 toPointRange( const LatticeSet& L ) const
509 {
510  return L.toPointRange();
511 }
512 
513 //-----------------------------------------------------------------------------
514 template <typename TKSpace>
515 typename DGtal::DigitalConvexity<TKSpace>::Integer
516 DGtal::DigitalConvexity<TKSpace>::
517 sizeStarCvxH( const PointRange& X ) const
518 {
519  PointRange cells;
520  // Computes Minkowski sum of Z with hypercube
521  PointRange Z = U( 0, X );
522  for ( Dimension k = 1; k < dimension; k++ )
523  Z = U( k, Z );
524  // Builds polytope
525  const auto P = makePolytope( Z );
526  // Extracts lattice points within polytope
527  // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
528  Counter C( P );
529  const Dimension a = C.longestAxis();
530  auto cellP = C.getLatticeCells( a );
531 
532  // Counts the number of cells
533  Integer nb = 0;
534  for ( const auto& value : cellP )
535  {
536  Point p = value.first;
537  Interval I = value.second;
538  nb += I.second - I.first + 1;
539  }
540  return nb;
541 }
542 
543 //-----------------------------------------------------------------------------
544 template <typename TKSpace>
545 typename DGtal::DigitalConvexity<TKSpace>::PointRange
546 DGtal::DigitalConvexity<TKSpace>::
547 Extr( const PointRange& C ) const
548 {
549  // JOL: using std::set< Point > or std::unordered_set< Point > is slightly slower.
550  // We prefer to use vector for easier vectorization.
551  std::vector<Point> E;
552  E.reserve( 2*C.size() );
553  for ( auto&& kp : C )
554  {
555  auto c = myK.uCell( kp );
556  if ( myK.uDim( c ) == 0 )
557  E.push_back( myK.uCoords( c ) );
558  else
559  {
560  auto faces = myK.uFaces( c );
561  for ( auto&& f : faces )
562  if ( myK.uDim( f ) == 0 )
563  E.push_back( myK.uCoords( f ) );
564  }
565  }
566  std::sort( E.begin(), E.end() );
567  auto last = std::unique( E.begin(), E.end() );
568  E.erase( last, E.end() );
569  return E;
570 }
571 //-----------------------------------------------------------------------------
572 template <typename TKSpace>
573 typename DGtal::DigitalConvexity<TKSpace>::PointRange
574 DGtal::DigitalConvexity<TKSpace>::
575 Extr( const LatticeSet& C ) const
576 {
577  return C.extremaOfCells();
578 }
579 //-----------------------------------------------------------------------------
580 template <typename TKSpace>
581 typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
582 DGtal::DigitalConvexity<TKSpace>::
583 Skel( const LatticeSet& C ) const
584 {
585  return C.skeletonOfCells();
586 }
587 //-----------------------------------------------------------------------------
588 template <typename TKSpace>
589 typename DGtal::DigitalConvexity<TKSpace>::PointRange
590 DGtal::DigitalConvexity<TKSpace>::
591 ExtrSkel( const LatticeSet& C ) const
592 {
593  return C.skeletonOfCells().extremaOfCells();
594 }
595 
596 //-----------------------------------------------------------------------------
597 template <typename TKSpace>
598 typename DGtal::DigitalConvexity<TKSpace>::PointRange
599 DGtal::DigitalConvexity<TKSpace>::
600 FC_direct( const PointRange& Z ) const
601 {
602  typedef typename LatticePolytope::Domain Domain;
603  PointRange cells;
604  // Computes Minkowski sum of Z with hypercube
605  PointRange X( Z );
606  for ( Dimension k = 0; k < dimension; k++ )
607  X = U( k, X );
608  // Builds polytope
609  const auto P = makePolytope( X );
610  // Extracts lattice points within polytope
611  // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
612  Counter C( P );
613  const Dimension a = C.longestAxis();
614  Point lo = C.lowerBound();
615  Point hi = C.upperBound();
616  hi[ a ] = lo[ a ];
617  const Domain projD( lo, hi ); //< the projected domain of the polytope.
618  const Point One = Point::diagonal( 1 );
619  //const Size size = projD.size(); //not USED
620  std::unordered_map< Point, Interval > cellP;
621  Point q;
622  for ( auto&& p : projD )
623  {
624  q = 2*p - One;
625  const auto I = C.intersectionIntervalAlongAxis( p, a );
626  const auto n = I.second - I.first;
627  if ( n != 0 )
628  {
629  // Now the second bound is included
630  cellP[ q ] = Interval( 2 * I.first - 1, 2 * I.second - 3 );
631  }
632  }
633  // It remains to compute all the k-cells, 0 <= k < d, intersected by Cvxh( Z )
634  for ( Dimension k = 0; k < dimension; k++ )
635  {
636  if ( k == a ) continue;
637  std::vector< Point > q_computed;
638  std::vector< Interval > I_computed;
639  for ( const auto& value : cellP )
640  {
641  Point p = value.first;
642  Interval I = value.second;
643  Point r = p; r[ k ] += 2;
644  const auto it = cellP.find( r );
645  if ( it == cellP.end() ) continue; // neighbor is empty
646  // Otherwise compute common part.
647  Interval J = it->second;
648  auto f = std::max( I.first, J.first );
649  auto s = std::min( I.second, J.second );
650  if ( f <= s )
651  {
652  Point qq = p; qq[ k ] += 1;
653  q_computed.push_back( qq );
654  I_computed.push_back( Interval( f, s ) );
655  }
656  }
657  // Add new columns to map Point -> column
658  for ( size_t i = 0; i < q_computed.size(); ++i )
659  {
660  cellP[ q_computed[ i ] ] = I_computed[ i ];
661  }
662  }
663  // The built complex is open.
664  // Check it and fill skelP
665  std::unordered_map< Point, std::vector< Interval > > skelP;
666  for ( const auto& value : cellP )
667  {
668  Point p = value.first;
669  Interval I = value.second;
670  auto n = I.second - I.first + 1;
671  if ( n % 2 == 0 )
672  trace.error() << "Weird thickness step 1="
673  << n << std::endl;
674  std::vector< Interval > V( 1, I );
675  auto result = skelP.insert( std::make_pair( p, V ) );
676  (void)result;//unused var;
677  }
678 
679  // std::cout << "Extract skel" << std::endl;
680  // Now extracting implicitly its Skel
681  for ( const auto& value : cellP )
682  {
683  const Point& p = value.first;
684  const auto& I = value.second;
685  for ( Dimension k = 0; k < dimension; k++ )
686  {
687  if ( k == a ) continue;
688  if ( ( p[ k ] & 0x1 ) != 0 ) continue; // if open along axis continue
689  // if closed, check upper incident cells along direction k
690  Point qq = p; qq[ k ] -= 1;
691  Point r = p; r[ k ] += 1;
692  auto itq = skelP.find( qq );
693  if ( itq != skelP.end() )
694  {
695  auto& W = itq->second;
696  eraseInterval( I, W );
697  // if ( W.empty() ) skelP.erase( itq );
698  }
699  auto itr = skelP.find( r );
700  if ( itr != skelP.end() )
701  {
702  auto& W = itr->second;
703  eraseInterval( I, W );
704  // if ( W.empty() ) skelP.erase( itr );
705  }
706  }
707  }
708  // Extract skel along main axis
709  for ( const auto& value : cellP )
710  {
711  const Point& p = value.first;
712  auto I = value.second;
713  const auto itp = skelP.find( p );
714  if ( itp == skelP.end() ) continue;
715  if ( ( I.first & 0x1 ) != 0 ) I.first += 1;
716  if ( ( I.second & 0x1 ) != 0 ) I.second -= 1;
717  auto& W = itp->second;
718  for ( auto x = I.first; x <= I.second; x += 2 )
719  {
720  eraseInterval( Interval{ x-1,x-1} , W );
721  eraseInterval( Interval{ x+1,x+1} , W );
722  }
723  }
724  // Erase empty stacks
725  for ( auto it = skelP.begin(), itE = skelP.end(); it != itE; )
726  {
727  const auto& V = it->second;
728  if ( V.empty() )
729  {
730  auto it_erase = it;
731  ++it;
732  skelP.erase( it_erase );
733  }
734  else ++it;
735  }
736  // Skel is constructed, now compute its Extr.
737  PointRange O;
738  for ( const auto& value : skelP )
739  {
740  Point p = value.first;
741  auto W = value.second;
742  for ( auto&& I : W )
743  {
744  p[ a ] = I.first;
745  O.push_back( p );
746  }
747  }
748  return Extr( O );
749 }
750 //-----------------------------------------------------------------------------
751 template <typename TKSpace>
752 typename DGtal::DigitalConvexity<TKSpace>::PointRange
753 DGtal::DigitalConvexity<TKSpace>::
754 FC_LatticeSet( const PointRange& Z ) const
755 {
756  const auto StarCvxZ = StarCvxH( Z );
757  const auto SkelStarCvxZ = StarCvxZ.skeletonOfCells().toPointRange();
758  return Extr( SkelStarCvxZ );
759 }
760 
761 //-----------------------------------------------------------------------------
762 template <typename TKSpace>
763 typename DGtal::DigitalConvexity<TKSpace>::PointRange
764 DGtal::DigitalConvexity<TKSpace>::
765 FC( const PointRange& Z, EnvelopeAlgorithm algo ) const
766 {
767  if ( algo == EnvelopeAlgorithm::DIRECT )
768  return FC_direct( Z );
769  else if ( algo == EnvelopeAlgorithm::LATTICE_SET )
770  return FC_LatticeSet( Z );
771  else
772  return Z;
773 }
774 
775 //-----------------------------------------------------------------------------
776 template <typename TKSpace>
777 typename DGtal::DigitalConvexity<TKSpace>::PointRange
778 DGtal::DigitalConvexity<TKSpace>::
779 envelope( const PointRange& Z, EnvelopeAlgorithm algo ) const
780 {
781  myDepthLastFCE = 0;
782  auto In = Z;
783  while (true) {
784  auto card_In = In.size();
785  In = FC( In, algo );
786  if ( In.size() == card_In ) return In;
787  myDepthLastFCE++;
788  }
789  trace.error() << "[DigitalConvexity::envelope] Should never pass here."
790  << std::endl;
791  return Z; // to avoid warnings.
792 }
793 
794 //-----------------------------------------------------------------------------
795 template <typename TKSpace>
796 typename DGtal::DigitalConvexity<TKSpace>::PointRange
797 DGtal::DigitalConvexity<TKSpace>::
798 relativeEnvelope( const PointRange& Z, const PointRange& Y,
799  EnvelopeAlgorithm algo ) const
800 {
801  myDepthLastFCE = 0;
802  auto In = Z;
803  while (true) {
804  PointRange Out;
805  std::set_intersection( In.cbegin(), In.cend(),
806  Y.cbegin(), Y.cend(),
807  std::back_inserter( Out ) );
808  In = FC( Out, algo );
809  if ( In.size() == Out.size() ) return Out;
810  myDepthLastFCE++;
811  }
812  return Z;
813 }
814 
815 //-----------------------------------------------------------------------------
816 template <typename TKSpace>
817 template <typename Predicate>
818 typename DGtal::DigitalConvexity<TKSpace>::PointRange
819 DGtal::DigitalConvexity<TKSpace>::
820 relativeEnvelope( const PointRange& Z, const Predicate& PredY,
821  EnvelopeAlgorithm algo ) const
822 {
823  myDepthLastFCE = 0;
824  auto In = Z;
825  while (true) {
826  auto Out = filter( In, PredY );
827  In = FC( Out, algo );
828  if ( In.size() == Out.size() ) return In;
829  myDepthLastFCE++;
830  }
831  return Z;
832 }
833 
834 //-----------------------------------------------------------------------------
835 template <typename TKSpace>
836 typename DGtal::DigitalConvexity<TKSpace>::Size
837 DGtal::DigitalConvexity<TKSpace>::
838 depthLastEnvelope() const
839 {
840  return myDepthLastFCE;
841 }
842 
843 //-----------------------------------------------------------------------------
844 template <typename TKSpace>
845 bool
846 DGtal::DigitalConvexity<TKSpace>::
847 isKConvex( const LatticePolytope& P, const Dimension k ) const
848 {
849  if ( k == 0 ) return true;
850  auto S = insidePoints( P );
851  auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
852  auto intersected_cells = makeCellCover( P, k, k );
853  return intersected_cells.nbCells() == touched_cells.nbCells();
854  // JOL: number should be enough
855  // && intersected_cells.subset( touched_cells );
856 }
857 
858 //-----------------------------------------------------------------------------
859 template <typename TKSpace>
860 bool
861 DGtal::DigitalConvexity<TKSpace>::
862 isFullyConvex( const LatticePolytope& P ) const
863 {
864  auto S = insidePoints( P );
865  for ( Dimension k = 1; k < KSpace::dimension; ++ k )
866  {
867  auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
868  auto intersected_cells = makeCellCover( P, k, k );
869  if ( ( intersected_cells.nbCells() != touched_cells.nbCells() ) )
870  // JOL: number should be enough
871  // || ( ! intersected_cells.subset( touched_cells ) ) )
872  return false;
873  }
874  return true;
875 }
876 
877 //-----------------------------------------------------------------------------
878 template <typename TKSpace>
879 bool
880 DGtal::DigitalConvexity<TKSpace>::
881 isKSubconvex( const LatticePolytope& P, const CellGeometry& C, const Dimension k ) const
882 {
883  auto intersected_cells = makeCellCover( P, k, k );
884  return intersected_cells.subset( C );
885 }
886 //-----------------------------------------------------------------------------
887 template <typename TKSpace>
888 bool
889 DGtal::DigitalConvexity<TKSpace>::
890 isFullySubconvex( const LatticePolytope& P, const CellGeometry& C ) const
891 {
892  auto intersected_cells = makeCellCover( P, C.minCellDim(), C.maxCellDim() );
893  return intersected_cells.subset( C );
894 }
895 
896 //-----------------------------------------------------------------------------
897 template <typename TKSpace>
898 bool
899 DGtal::DigitalConvexity<TKSpace>::
900 isFullySubconvex( const LatticePolytope& P, const LatticeSet& StarX ) const
901 {
902  LatticePolytope Q = P + typename LatticePolytope::UnitSegment( 0 );
903  for ( Dimension k = 1; k < dimension; k++ )
904  Q = Q + typename LatticePolytope::UnitSegment( k );
905  Counter C( Q );
906  const auto cells = C.getLatticeCells( StarX.axis() );
907  LatticeSet StarP( cells, StarX.axis() );
908  return StarX.includes( StarP );
909 }
910 
911 //-----------------------------------------------------------------------------
912 template <typename TKSpace>
913 bool
914 DGtal::DigitalConvexity<TKSpace>::
915 isFullySubconvex( const PointRange& Y, const LatticeSet& StarX ) const
916 {
917  const auto SCY = StarCvxH( Y, StarX.axis() );
918  return StarX.includes( SCY );
919 }
920 
921 //-----------------------------------------------------------------------------
922 template <typename TKSpace>
923 bool
924 DGtal::DigitalConvexity<TKSpace>::
925 isFullySubconvex( const Point& a, const Point& b, const Point& c,
926  const LatticeSet& StarX ) const
927 {
928  ASSERT( dimension == 3 );
929  const auto SCabc = StarCvxH( a, b, c, StarX.axis() );
930  return StarX.includes( SCabc );
931 }
932 
933 //-----------------------------------------------------------------------------
934 template <typename TKSpace>
935 bool
936 DGtal::DigitalConvexity<TKSpace>::
937 isKSubconvex( const Point& a, const Point& b,
938  const CellGeometry& C, const Dimension k ) const
939 {
940  CellGeometry cgeom( myK, k, k, false );
941  cgeom.addCellsTouchingSegment( a, b );
942  return cgeom.subset( C );
943 }
944 
945 //-----------------------------------------------------------------------------
946 template <typename TKSpace>
947 bool
948 DGtal::DigitalConvexity<TKSpace>::
949 isFullySubconvex( const Point& a, const Point& b,
950  const CellGeometry& C ) const
951 {
952  CellGeometry cgeom( myK, C.minCellDim(), C.maxCellDim(), false );
953  cgeom.addCellsTouchingSegment( a, b );
954  return cgeom.subset( C );
955 }
956 
957 //-----------------------------------------------------------------------------
958 template <typename TKSpace>
959 bool
960 DGtal::DigitalConvexity<TKSpace>::
961 isFullySubconvex( const Point& a, const Point& b,
962  const LatticeSet& StarX ) const
963 {
964  LatticeSet L_ab( StarX.axis() );
965  const auto V = b - a;
966  L_ab.insert( 2*a );
967  for ( Dimension k = 0; k < dimension; k++ )
968  {
969  const Integer n = ( V[ k ] >= 0 ) ? V[ k ] : -V[ k ];
970  const Integer d = ( V[ k ] >= 0 ) ? 1 : -1;
971  if ( n == 0 ) continue;
972  Point kc;
973  for ( Integer i = 1; i < n; i++ )
974  {
975  for ( Dimension j = 0; j < dimension; j++ )
976  {
977  if ( j == k ) kc[ k ] = 2 * ( a[ k ] + d * i );
978  else
979  {
980  const auto v = V[ j ];
981  const auto q = ( v * i ) / n;
982  const auto r = ( v * i ) % n; // might be negative
983  kc[ j ] = 2 * ( a[ j ] + q );
984  if ( r < 0 ) kc[ j ] -= 1;
985  else if ( r > 0 ) kc[ j ] += 1;
986  }
987  }
988  L_ab.insert( kc );
989  }
990  }
991  if ( a != b ) L_ab.insert( 2*b );
992  LatticeSet Star_ab = L_ab.starOfCells();
993  return StarX.includes( Star_ab );
994 }
995 
996 //-----------------------------------------------------------------------------
997 template <typename TKSpace>
998 bool
999 DGtal::DigitalConvexity<TKSpace>::
1000 isKConvex( const RationalPolytope& P, const Dimension k ) const
1001 {
1002  if ( k == 0 ) return true;
1003  auto S = insidePoints( P );
1004  auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
1005  auto intersected_cells = makeCellCover( P, k, k );
1006  return intersected_cells.nbCells() == touched_cells.nbCells()
1007  && intersected_cells.subset( touched_cells );
1008 }
1009 
1010 //-----------------------------------------------------------------------------
1011 template <typename TKSpace>
1012 bool
1013 DGtal::DigitalConvexity<TKSpace>::
1014 isFullyConvex( const RationalPolytope& 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  || ( ! intersected_cells.subset( touched_cells ) ) )
1023  return false;
1024  }
1025  return true;
1026 }
1027 
1028 //-----------------------------------------------------------------------------
1029 template <typename TKSpace>
1030 bool
1031 DGtal::DigitalConvexity<TKSpace>::
1032 isKSubconvex( const RationalPolytope& P, const CellGeometry& C,
1033  const Dimension k ) const
1034 {
1035  auto intersected_cells = makeCellCover( P, k, k );
1036  return intersected_cells.subset( C );
1037 }
1038 //-----------------------------------------------------------------------------
1039 template <typename TKSpace>
1040 bool
1041 DGtal::DigitalConvexity<TKSpace>::
1042 isFullySubconvex( const RationalPolytope& 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 //-----------------------------------------------------------------------------
1049 template <typename TKSpace>
1050 void
1051 DGtal::DigitalConvexity<TKSpace>::
1052 eraseInterval( Interval I, std::vector< Interval > & V )
1053 {
1054  for ( std::size_t i = 0; i < V.size(); )
1055  {
1056  Interval& J = V[ i ];
1057  // I=[a,b], J=[a',b'], a <= b, a' <= b'
1058  if ( I.second < J.first ) { break; } // b < a' : no further intersection
1059  if ( J.second < I.first ) { ++i; continue; } // b' < a : no further intersection
1060  // a' <= b and a <= b'
1061  // a ---------- b
1062  // a' ............... a'
1063  // b' ................. b'
1064 
1065  // a' ..................... b' => a'..a-1 b+1 b'
1066  Interval K1( J.first, I.first - 1 );
1067  Interval K2( I.second + 1, J.second );
1068  bool K1_exist = K1.second >= K1.first;
1069  bool K2_exist = K2.second >= K2.first;
1070  if ( K1_exist && K2_exist )
1071  {
1072  V[ i ] = K2;
1073  V.insert( V.begin() + i, K1 );
1074  break; // no further intersection possible
1075  }
1076  else if ( K1_exist )
1077  {
1078  V[ i ] = K1; i++;
1079  }
1080  else if ( K2_exist )
1081  {
1082  V[ i ] = K2; break;
1083  }
1084  else
1085  {
1086  V.erase( V.begin() + i );
1087  }
1088  }
1089 }
1090 
1091 //-----------------------------------------------------------------------------
1092 template <typename TKSpace>
1093 template <typename Predicate>
1094 typename DGtal::DigitalConvexity<TKSpace>::PointRange
1095 DGtal::DigitalConvexity<TKSpace>::
1096 filter( const PointRange& E, const Predicate& Pred )
1097 {
1098  PointRange Out;
1099  Out.reserve( E.size() );
1100  for ( auto&& p : E )
1101  if ( Pred( p ) ) Out.push_back( p );
1102  return Out;
1103 }
1104 
1105 ///////////////////////////////////////////////////////////////////////////////
1106 // Interface - public :
1107 
1108 /**
1109  * Writes/Displays the object on an output stream.
1110  * @param out the output stream where the object is written.
1111  */
1112 template <typename TKSpace>
1113 inline
1114 void
1115 DGtal::DigitalConvexity<TKSpace>::selfDisplay ( std::ostream & out ) const
1116 {
1117  out << "[DigitalConvexity]";
1118 }
1119 
1120 /**
1121  * Checks the validity/consistency of the object.
1122  * @return 'true' if the object is valid, 'false' otherwise.
1123  */
1124 template <typename TKSpace>
1125 inline
1126 bool
1127 DGtal::DigitalConvexity<TKSpace>::isValid() const
1128 {
1129  return true;
1130 }
1131 
1132 
1133 ///////////////////////////////////////////////////////////////////////////////
1134 // Implementation of inline functions //
1135 
1136 //-----------------------------------------------------------------------------
1137 template <typename TKSpace>
1138 inline
1139 std::ostream&
1140 DGtal::operator<< ( std::ostream & out,
1141  const DigitalConvexity<TKSpace> & object )
1142 {
1143  object.selfDisplay( out );
1144  return out;
1145 }
1146 
1147 // //
1148 ///////////////////////////////////////////////////////////////////////////////