DGtal  1.3.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 )
39  : myK( K )
40 {
41 }
42 //-----------------------------------------------------------------------------
43 template <typename TKSpace>
44 DGtal::DigitalConvexity<TKSpace>::
45 DigitalConvexity( Point lo, Point hi )
46 {
47  myK.init( lo, hi, true );
48 }
49 
50 //-----------------------------------------------------------------------------
51 template <typename TKSpace>
52 const TKSpace&
53 DGtal::DigitalConvexity<TKSpace>::
54 space() const
55 {
56  return myK;
57 }
58 
59 //-----------------------------------------------------------------------------
60 template <typename TKSpace>
61 template <typename PointIterator>
62 typename DGtal::DigitalConvexity<TKSpace>::LatticePolytope
63 DGtal::DigitalConvexity<TKSpace>::
64 makeSimplex( PointIterator itB, PointIterator itE )
65 {
66  return LatticePolytope( itB, itE );
67 }
68 
69 //-----------------------------------------------------------------------------
70 template <typename TKSpace>
71 typename DGtal::DigitalConvexity<TKSpace>::LatticePolytope
72 DGtal::DigitalConvexity<TKSpace>::
73 makeSimplex( std::initializer_list<Point> l )
74 {
75  return LatticePolytope( l );
76 }
77 
78 //-----------------------------------------------------------------------------
79 template <typename TKSpace>
80 template <typename PointIterator>
81 typename DGtal::DigitalConvexity<TKSpace>::RationalPolytope
82 DGtal::DigitalConvexity<TKSpace>::
83 makeRationalSimplex( Integer d, PointIterator itB, PointIterator itE )
84 {
85  return RationalPolytope( d, itB, itE );
86 }
87 
88 //-----------------------------------------------------------------------------
89 template <typename TKSpace>
90 typename DGtal::DigitalConvexity<TKSpace>::RationalPolytope
91 DGtal::DigitalConvexity<TKSpace>::
92 makeRationalSimplex( std::initializer_list<Point> l )
93 {
94  return RationalPolytope( l );
95 }
96 
97 //-----------------------------------------------------------------------------
98 template <typename TKSpace>
99 template <typename PointIterator>
100 bool
101 DGtal::DigitalConvexity<TKSpace>::
102 isSimplexFullDimensional( PointIterator itB, PointIterator itE )
103 {
104  typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
105  const Dimension d = KSpace::dimension;
106  std::vector<Point> pts( d+1 );
107  Dimension k = 0;
108  for ( ; itB != itE && k <= d; ++itB, ++k ) pts[ k ] = *itB;
109  // A simplex has exactly d+1 vertices.
110  if ( k != d+1 ) return false;
111  Matrix M;
112  for ( Dimension i = 0; i < d; ++i )
113  for ( Dimension j = 0; j < d; ++j )
114  M.setComponent( i, j, pts[ i ][ j ] - pts[ d ][ j ] );
115  // A simplex has its vectors linearly independent.
116  return M.determinant() != 0;
117 }
118 
119 //-----------------------------------------------------------------------------
120 template <typename TKSpace>
121 bool
122 DGtal::DigitalConvexity<TKSpace>::
123 isSimplexFullDimensional( std::initializer_list<Point> l )
124 {
125  return isSimplexFullDimensional( l.begin(), l.end() );
126 }
127 
128 //-----------------------------------------------------------------------------
129 template <typename TKSpace>
130 template <typename PointIterator>
131 typename DGtal::DigitalConvexity<TKSpace>::SimplexType
132 DGtal::DigitalConvexity<TKSpace>::
133 simplexType( PointIterator itB, PointIterator itE )
134 {
135  typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
136  const Dimension d = KSpace::dimension;
137  std::vector<Point> pts( d+1 );
138  Dimension k = 0;
139  for ( ; itB != itE && k <= d; ++itB, ++k ) pts[ k ] = *itB;
140  // A simplex has exactly d+1 vertices.
141  if ( k != d+1 ) return SimplexType::INVALID;
142  Matrix M;
143  for ( Dimension i = 0; i < d; ++i )
144  for ( Dimension j = 0; j < d; ++j )
145  M.setComponent( i, j, pts[ i ][ j ] - pts[ d ][ j ] );
146  // A simplex has its vectors linearly independent.
147  auto V = M.determinant();
148  return (V == 0) ? SimplexType::DEGENERATED
149  : ( ((V == 1) || (V==-1)) ? SimplexType::UNITARY : SimplexType::COMMON );
150 }
151 
152 //-----------------------------------------------------------------------------
153 template <typename TKSpace>
154 void
155 DGtal::DigitalConvexity<TKSpace>::
156 displaySimplex( std::ostream& out, std::initializer_list<Point> l )
157 {
158  displaySimplex( out, l.begin(), l.end() );
159 }
160 
161 //-----------------------------------------------------------------------------
162 template <typename TKSpace>
163 template <typename PointIterator>
164 void
165 DGtal::DigitalConvexity<TKSpace>::
166 displaySimplex( std::ostream& out, PointIterator itB, PointIterator itE )
167 {
168  typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
169  const Dimension d = KSpace::dimension;
170  std::vector<Point> pts( d+1 );
171  Dimension k = 0;
172  for ( ; itB != itE && k <= d; ++itB, ++k ) pts[ k ] = *itB;
173  // A simplex has exactly d+1 vertices.
174  if ( k != d+1 ) { out << "[SPLX INVALID]"; return; }
175  Matrix M;
176  for ( Dimension i = 0; i < d; ++i )
177  for ( Dimension k = 0; k < d; ++k )
178  M.setComponent( i, k, pts[ i ][ k ] - pts[ d ][ k ] );
179  // A simplex has its vectors linearly independent.
180  auto V = M.determinant();
181  out << "[SPLX V=" << V;
182  for ( Dimension i = 0; i < d; ++i ) {
183  out << " (";
184  for ( Dimension j = 0; j < d; ++j )
185  out << " " << M( i, j );
186  out << " )";
187  }
188  out << " ]";
189 }
190 
191 //-----------------------------------------------------------------------------
192 template <typename TKSpace>
193 typename DGtal::DigitalConvexity<TKSpace>::SimplexType
194 DGtal::DigitalConvexity<TKSpace>::
195 simplexType( std::initializer_list<Point> l )
196 {
197  return simplexType( l.begin(), l.end() );
198 }
199 
200 
201 //-----------------------------------------------------------------------------
202 template <typename TKSpace>
203 typename DGtal::DigitalConvexity<TKSpace>::PointRange
204 DGtal::DigitalConvexity<TKSpace>::
205 insidePoints( const LatticePolytope& polytope )
206 {
207  PointRange pts;
208  polytope.getPoints( pts );
209  return pts;
210 }
211 //-----------------------------------------------------------------------------
212 template <typename TKSpace>
213 typename DGtal::DigitalConvexity<TKSpace>::PointRange
214 DGtal::DigitalConvexity<TKSpace>::
215 interiorPoints( const LatticePolytope& polytope )
216 {
217  PointRange pts;
218  polytope.getInteriorPoints( pts );
219  return pts;
220 }
221 
222 //-----------------------------------------------------------------------------
223 template <typename TKSpace>
224 typename DGtal::DigitalConvexity<TKSpace>::PointRange
225 DGtal::DigitalConvexity<TKSpace>::
226 insidePoints( const RationalPolytope& polytope )
227 {
228  PointRange pts;
229  polytope.getPoints( pts );
230  return pts;
231 }
232 //-----------------------------------------------------------------------------
233 template <typename TKSpace>
234 typename DGtal::DigitalConvexity<TKSpace>::PointRange
235 DGtal::DigitalConvexity<TKSpace>::
236 interiorPoints( const RationalPolytope& polytope )
237 {
238  PointRange pts;
239  polytope.getInteriorPoints( pts );
240  return pts;
241 }
242 
243 //-----------------------------------------------------------------------------
244 template <typename TKSpace>
245 template <typename PointIterator>
246 typename DGtal::DigitalConvexity<TKSpace>::CellGeometry
247 DGtal::DigitalConvexity<TKSpace>::
248 makeCellCover( PointIterator itB, PointIterator itE,
249  Dimension i, Dimension k ) const
250 {
251  ASSERT( 0 <= i );
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( 0 <= i );
267  ASSERT( i <= k );
268  ASSERT( k <= KSpace::dimension );
269  CellGeometry cgeom( myK, i, k, false );
270  cgeom.addCellsTouchingPolytope( P );
271  return cgeom;
272 }
273 
274 //-----------------------------------------------------------------------------
275 template <typename TKSpace>
276 typename DGtal::DigitalConvexity<TKSpace>::CellGeometry
277 DGtal::DigitalConvexity<TKSpace>::
278 makeCellCover( const RationalPolytope& P,
279  Dimension i, Dimension k ) const
280 {
281  ASSERT( 0 <= i );
282  ASSERT( i <= k );
283  ASSERT( k <= KSpace::dimension );
284  CellGeometry cgeom( myK, i, k, false );
285  cgeom.addCellsTouchingPolytope( P );
286  return cgeom;
287 }
288 
289 //-----------------------------------------------------------------------------
290 template <typename TKSpace>
291 typename DGtal::DigitalConvexity<TKSpace>::LatticePolytope
292 DGtal::DigitalConvexity<TKSpace>::
293 makePolytope( const PointRange& X, bool safe ) const
294 {
295  if ( safe )
296  {
297  typedef typename detail::ConvexityHelperInternalInteger< Integer, true >::Type
298  InternalInteger;
299  return ConvexityHelper< dimension, Integer, InternalInteger >::
300  computeLatticePolytope( X, false, false );
301  }
302  else
303  {
304  typedef typename detail::ConvexityHelperInternalInteger< Integer, false >::Type
305  InternalInteger;
306  return ConvexityHelper< dimension, Integer, InternalInteger >::
307  computeLatticePolytope( X, false, false );
308  }
309 }
310 
311 //-----------------------------------------------------------------------------
312 template <typename TKSpace>
313 typename DGtal::DigitalConvexity<TKSpace>::PointRange
314 DGtal::DigitalConvexity<TKSpace>::
315 U( Dimension i, const PointRange& X ) const
316 {
317  PointRange Y( X );
318  PointRange Z;
319  Z.reserve( X.size() );
320  auto add_one = [&i] ( Point & p ) { p[ i ] += 1; };
321  std::for_each( Y.begin(), Y.end(), add_one );
322  std::set_union( X.cbegin(), X.cend(), Y.cbegin(), Y.cend(),
323  std::back_inserter( Z ) );
324  return Z;
325 }
326 
327 //-----------------------------------------------------------------------------
328 template <typename TKSpace>
329 bool
330 DGtal::DigitalConvexity<TKSpace>::
331 is0Convex( const PointRange& X, bool safe ) const
332 {
333  if ( X.empty() ) return true;
334  const auto P = makePolytope( X, safe );
335  return P.count() == X.size();
336 }
337 
338 //-----------------------------------------------------------------------------
339 template <typename TKSpace>
340 bool
341 DGtal::DigitalConvexity<TKSpace>::
342 isFullyConvex( const PointRange& Z, bool convex0, bool safe ) const
343 {
344  ASSERT( dimension <= 64 );
345  typedef DGtal::int64_t Direction;
346  typedef std::vector< Direction > Directions;
347  std::array< Directions, dimension > C;
348  const bool cvx0 = convex0 ? true : is0Convex( Z, safe );
349  if ( ! cvx0 ) return false;
350  C[ 0 ].push_back( (Direction) 0 );
351  std::map< Direction, PointRange > X;
352  X[ 0 ] = Z;
353  std::sort( X[ 0 ].begin(), X[ 0 ].end() );
354  for ( Dimension k = 1; k < dimension; k++ )
355  {
356  for ( const auto beta : C[ k-1 ] )
357  {
358  for ( Dimension j = 0; j < dimension; j++ )
359  {
360  const Direction dir_j = Direction(1) << j;
361  if ( beta < dir_j )
362  {
363  const Direction alpha = beta | dir_j;
364  C[ k ].push_back( alpha );
365  X[ alpha ] = U( j, X[ beta ] );
366  if ( ! is0Convex( X[ alpha ] ) ) return false;
367  }
368  }
369  }
370  }
371  return true;
372 }
373 
374 //-----------------------------------------------------------------------------
375 template <typename TKSpace>
376 bool
377 DGtal::DigitalConvexity<TKSpace>::
378 isKConvex( const LatticePolytope& P, const Dimension k ) const
379 {
380  if ( k == 0 ) return true;
381  auto S = insidePoints( P );
382  auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
383  auto intersected_cells = makeCellCover( P, k, k );
384  return intersected_cells.nbCells() == touched_cells.nbCells()
385  && intersected_cells.subset( touched_cells );
386 }
387 
388 //-----------------------------------------------------------------------------
389 template <typename TKSpace>
390 bool
391 DGtal::DigitalConvexity<TKSpace>::
392 isFullyConvex( const LatticePolytope& P ) const
393 {
394  auto S = insidePoints( P );
395  for ( Dimension k = 1; k < KSpace::dimension; ++ k )
396  {
397  auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
398  auto intersected_cells = makeCellCover( P, k, k );
399  if ( ( intersected_cells.nbCells() != touched_cells.nbCells() )
400  || ( ! intersected_cells.subset( touched_cells ) ) )
401  return false;
402  }
403  return true;
404 }
405 
406 //-----------------------------------------------------------------------------
407 template <typename TKSpace>
408 bool
409 DGtal::DigitalConvexity<TKSpace>::
410 isKSubconvex( const LatticePolytope& P, const CellGeometry& C, const Dimension k ) const
411 {
412  auto intersected_cells = makeCellCover( P, k, k );
413  return intersected_cells.subset( C );
414 }
415 //-----------------------------------------------------------------------------
416 template <typename TKSpace>
417 bool
418 DGtal::DigitalConvexity<TKSpace>::
419 isFullySubconvex( const LatticePolytope& P, const CellGeometry& C ) const
420 {
421  auto intersected_cells = makeCellCover( P, C.minCellDim(), C.maxCellDim() );
422  return intersected_cells.subset( C );
423 }
424 
425 //-----------------------------------------------------------------------------
426 template <typename TKSpace>
427 bool
428 DGtal::DigitalConvexity<TKSpace>::
429 isKSubconvex( const Point& a, const Point& b,
430  const CellGeometry& C, const Dimension k ) const
431 {
432  CellGeometry cgeom( myK, k, k, false );
433  cgeom.addCellsTouchingSegment( a, b );
434  return cgeom.subset( C );
435 }
436 
437 //-----------------------------------------------------------------------------
438 template <typename TKSpace>
439 bool
440 DGtal::DigitalConvexity<TKSpace>::
441 isFullySubconvex( const Point& a, const Point& b,
442  const CellGeometry& C ) const
443 {
444  CellGeometry cgeom( myK, C.minCellDim(), C.maxCellDim(), false );
445  cgeom.addCellsTouchingSegment( a, b );
446  return cgeom.subset( C );
447 }
448 
449 //-----------------------------------------------------------------------------
450 template <typename TKSpace>
451 bool
452 DGtal::DigitalConvexity<TKSpace>::
453 isKConvex( const RationalPolytope& P, const Dimension k ) const
454 {
455  if ( k == 0 ) return true;
456  auto S = insidePoints( P );
457  auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
458  auto intersected_cells = makeCellCover( P, k, k );
459  return intersected_cells.nbCells() == touched_cells.nbCells()
460  && intersected_cells.subset( touched_cells );
461 }
462 
463 //-----------------------------------------------------------------------------
464 template <typename TKSpace>
465 bool
466 DGtal::DigitalConvexity<TKSpace>::
467 isFullyConvex( const RationalPolytope& P ) const
468 {
469  auto S = insidePoints( P );
470  for ( Dimension k = 1; k < KSpace::dimension; ++ k )
471  {
472  auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
473  auto intersected_cells = makeCellCover( P, k, k );
474  if ( ( intersected_cells.nbCells() != touched_cells.nbCells() )
475  || ( ! intersected_cells.subset( touched_cells ) ) )
476  return false;
477  }
478  return true;
479 }
480 
481 //-----------------------------------------------------------------------------
482 template <typename TKSpace>
483 bool
484 DGtal::DigitalConvexity<TKSpace>::
485 isKSubconvex( const RationalPolytope& P, const CellGeometry& C,
486  const Dimension k ) const
487 {
488  auto intersected_cells = makeCellCover( P, k, k );
489  return intersected_cells.subset( C );
490 }
491 //-----------------------------------------------------------------------------
492 template <typename TKSpace>
493 bool
494 DGtal::DigitalConvexity<TKSpace>::
495 isFullySubconvex( const RationalPolytope& P, const CellGeometry& C ) const
496 {
497  auto intersected_cells = makeCellCover( P, C.minCellDim(), C.maxCellDim() );
498  return intersected_cells.subset( C );
499 }
500 
501 ///////////////////////////////////////////////////////////////////////////////
502 // Interface - public :
503 
504 /**
505  * Writes/Displays the object on an output stream.
506  * @param out the output stream where the object is written.
507  */
508 template <typename TKSpace>
509 inline
510 void
511 DGtal::DigitalConvexity<TKSpace>::selfDisplay ( std::ostream & out ) const
512 {
513  out << "[DigitalConvexity]";
514 }
515 
516 /**
517  * Checks the validity/consistency of the object.
518  * @return 'true' if the object is valid, 'false' otherwise.
519  */
520 template <typename TKSpace>
521 inline
522 bool
523 DGtal::DigitalConvexity<TKSpace>::isValid() const
524 {
525  return true;
526 }
527 
528 
529 ///////////////////////////////////////////////////////////////////////////////
530 // Implementation of inline functions //
531 
532 //-----------------------------------------------------------------------------
533 template <typename TKSpace>
534 inline
535 std::ostream&
536 DGtal::operator<< ( std::ostream & out,
537  const DigitalConvexity<TKSpace> & object )
538 {
539  object.selfDisplay( out );
540  return out;
541 }
542 
543 // //
544 ///////////////////////////////////////////////////////////////////////////////
Generated on Mon Jun 20 2022 18:23:28 for DGtal by   doxygen 1.8.17