DGtal  1.2.0
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 //////////////////////////////////////////////////////////////////////////////
33 
34 //-----------------------------------------------------------------------------
35 template <typename TKSpace>
36 DGtal::DigitalConvexity<TKSpace>::
37 DigitalConvexity( Clone<KSpace> K )
38  : myK( K )
39 {
40 }
41 //-----------------------------------------------------------------------------
42 template <typename TKSpace>
43 DGtal::DigitalConvexity<TKSpace>::
44 DigitalConvexity( Point lo, Point hi )
45 {
46  myK.init( lo, hi, true );
47 }
48 
49 //-----------------------------------------------------------------------------
50 template <typename TKSpace>
51 const TKSpace&
52 DGtal::DigitalConvexity<TKSpace>::
53 space() const
54 {
55  return myK;
56 }
57 
58 //-----------------------------------------------------------------------------
59 template <typename TKSpace>
60 template <typename PointIterator>
61 typename DGtal::DigitalConvexity<TKSpace>::LatticePolytope
62 DGtal::DigitalConvexity<TKSpace>::
63 makeSimplex( PointIterator itB, PointIterator itE )
64 {
65  return LatticePolytope( itB, itE );
66 }
67 
68 //-----------------------------------------------------------------------------
69 template <typename TKSpace>
70 typename DGtal::DigitalConvexity<TKSpace>::LatticePolytope
71 DGtal::DigitalConvexity<TKSpace>::
72 makeSimplex( std::initializer_list<Point> l )
73 {
74  return LatticePolytope( l );
75 }
76 
77 //-----------------------------------------------------------------------------
78 template <typename TKSpace>
79 template <typename PointIterator>
80 typename DGtal::DigitalConvexity<TKSpace>::RationalPolytope
81 DGtal::DigitalConvexity<TKSpace>::
82 makeRationalSimplex( Integer d, PointIterator itB, PointIterator itE )
83 {
84  return RationalPolytope( d, itB, itE );
85 }
86 
87 //-----------------------------------------------------------------------------
88 template <typename TKSpace>
89 typename DGtal::DigitalConvexity<TKSpace>::RationalPolytope
90 DGtal::DigitalConvexity<TKSpace>::
91 makeRationalSimplex( std::initializer_list<Point> l )
92 {
93  return RationalPolytope( l );
94 }
95 
96 //-----------------------------------------------------------------------------
97 template <typename TKSpace>
98 template <typename PointIterator>
99 bool
100 DGtal::DigitalConvexity<TKSpace>::
101 isSimplexFullDimensional( PointIterator itB, PointIterator itE )
102 {
103  typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
104  const Dimension d = KSpace::dimension;
105  std::vector<Point> pts( d+1 );
106  Dimension k = 0;
107  for ( ; itB != itE && k <= d; ++itB, ++k ) pts[ k ] = *itB;
108  // A simplex has exactly d+1 vertices.
109  if ( k != d+1 ) return false;
110  Matrix M;
111  for ( Dimension i = 0; i < d; ++i )
112  for ( Dimension j = 0; j < d; ++j )
113  M.setComponent( i, j, pts[ i ][ j ] - pts[ d ][ j ] );
114  // A simplex has its vectors linearly independent.
115  return M.determinant() != 0;
116 }
117 
118 //-----------------------------------------------------------------------------
119 template <typename TKSpace>
120 bool
121 DGtal::DigitalConvexity<TKSpace>::
122 isSimplexFullDimensional( std::initializer_list<Point> l )
123 {
124  return isSimplexFullDimensional( l.begin(), l.end() );
125 }
126 
127 //-----------------------------------------------------------------------------
128 template <typename TKSpace>
129 template <typename PointIterator>
130 typename DGtal::DigitalConvexity<TKSpace>::SimplexType
131 DGtal::DigitalConvexity<TKSpace>::
132 simplexType( PointIterator itB, PointIterator itE )
133 {
134  typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
135  const Dimension d = KSpace::dimension;
136  std::vector<Point> pts( d+1 );
137  Dimension k = 0;
138  for ( ; itB != itE && k <= d; ++itB, ++k ) pts[ k ] = *itB;
139  // A simplex has exactly d+1 vertices.
140  if ( k != d+1 ) return SimplexType::INVALID;
141  Matrix M;
142  for ( Dimension i = 0; i < d; ++i )
143  for ( Dimension j = 0; j < d; ++j )
144  M.setComponent( i, j, pts[ i ][ j ] - pts[ d ][ j ] );
145  // A simplex has its vectors linearly independent.
146  auto V = M.determinant();
147  return (V == 0) ? SimplexType::DEGENERATED
148  : ( ((V == 1) || (V==-1)) ? SimplexType::UNITARY : SimplexType::COMMON );
149 }
150 
151 //-----------------------------------------------------------------------------
152 template <typename TKSpace>
153 void
154 DGtal::DigitalConvexity<TKSpace>::
155 displaySimplex( std::ostream& out, std::initializer_list<Point> l )
156 {
157  displaySimplex( out, l.begin(), l.end() );
158 }
159 
160 //-----------------------------------------------------------------------------
161 template <typename TKSpace>
162 template <typename PointIterator>
163 void
164 DGtal::DigitalConvexity<TKSpace>::
165 displaySimplex( std::ostream& out, PointIterator itB, PointIterator itE )
166 {
167  typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
168  const Dimension d = KSpace::dimension;
169  std::vector<Point> pts( d+1 );
170  Dimension k = 0;
171  for ( ; itB != itE && k <= d; ++itB, ++k ) pts[ k ] = *itB;
172  // A simplex has exactly d+1 vertices.
173  if ( k != d+1 ) { out << "[SPLX INVALID]"; return; }
174  Matrix M;
175  for ( Dimension i = 0; i < d; ++i )
176  for ( Dimension k = 0; k < d; ++k )
177  M.setComponent( i, k, pts[ i ][ k ] - pts[ d ][ k ] );
178  // A simplex has its vectors linearly independent.
179  auto V = M.determinant();
180  out << "[SPLX V=" << V;
181  for ( Dimension i = 0; i < d; ++i ) {
182  out << " (";
183  for ( Dimension j = 0; j < d; ++j )
184  out << " " << M( i, j );
185  out << " )";
186  }
187  out << " ]";
188 }
189 
190 //-----------------------------------------------------------------------------
191 template <typename TKSpace>
192 typename DGtal::DigitalConvexity<TKSpace>::SimplexType
193 DGtal::DigitalConvexity<TKSpace>::
194 simplexType( std::initializer_list<Point> l )
195 {
196  return simplexType( l.begin(), l.end() );
197 }
198 
199 
200 //-----------------------------------------------------------------------------
201 template <typename TKSpace>
202 typename DGtal::DigitalConvexity<TKSpace>::PointRange
203 DGtal::DigitalConvexity<TKSpace>::
204 insidePoints( const LatticePolytope& polytope )
205 {
206  PointRange pts;
207  polytope.getPoints( pts );
208  return pts;
209 }
210 //-----------------------------------------------------------------------------
211 template <typename TKSpace>
212 typename DGtal::DigitalConvexity<TKSpace>::PointRange
213 DGtal::DigitalConvexity<TKSpace>::
214 interiorPoints( const LatticePolytope& polytope )
215 {
216  PointRange pts;
217  polytope.getInteriorPoints( pts );
218  return pts;
219 }
220 
221 //-----------------------------------------------------------------------------
222 template <typename TKSpace>
223 typename DGtal::DigitalConvexity<TKSpace>::PointRange
224 DGtal::DigitalConvexity<TKSpace>::
225 insidePoints( const RationalPolytope& polytope )
226 {
227  PointRange pts;
228  polytope.getPoints( pts );
229  return pts;
230 }
231 //-----------------------------------------------------------------------------
232 template <typename TKSpace>
233 typename DGtal::DigitalConvexity<TKSpace>::PointRange
234 DGtal::DigitalConvexity<TKSpace>::
235 interiorPoints( const RationalPolytope& polytope )
236 {
237  PointRange pts;
238  polytope.getInteriorPoints( pts );
239  return pts;
240 }
241 
242 //-----------------------------------------------------------------------------
243 template <typename TKSpace>
244 template <typename PointIterator>
245 typename DGtal::DigitalConvexity<TKSpace>::CellGeometry
246 DGtal::DigitalConvexity<TKSpace>::
247 makeCellCover( PointIterator itB, PointIterator itE,
248  Dimension i, Dimension k ) const
249 {
250  ASSERT( 0 <= i );
251  ASSERT( i <= k );
252  ASSERT( k <= KSpace::dimension );
253  CellGeometry cgeom( myK, i, k, false );
254  cgeom.addCellsTouchingPoints( itB, itE );
255  return cgeom;
256 }
257 
258 //-----------------------------------------------------------------------------
259 template <typename TKSpace>
260 typename DGtal::DigitalConvexity<TKSpace>::CellGeometry
261 DGtal::DigitalConvexity<TKSpace>::
262 makeCellCover( const LatticePolytope& P,
263  Dimension i, Dimension k ) const
264 {
265  ASSERT( 0 <= i );
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( 0 <= i );
281  ASSERT( i <= k );
282  ASSERT( k <= KSpace::dimension );
283  CellGeometry cgeom( myK, i, k, false );
284  cgeom.addCellsTouchingPolytope( P );
285  return cgeom;
286 }
287 
288 //-----------------------------------------------------------------------------
289 template <typename TKSpace>
290 bool
291 DGtal::DigitalConvexity<TKSpace>::
292 isKConvex( const LatticePolytope& P, const Dimension k ) const
293 {
294  if ( k == 0 ) return true;
295  auto S = insidePoints( P );
296  auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
297  auto intersected_cells = makeCellCover( P, k, k );
298  return intersected_cells.nbCells() == touched_cells.nbCells()
299  && intersected_cells.subset( touched_cells );
300 }
301 
302 //-----------------------------------------------------------------------------
303 template <typename TKSpace>
304 bool
305 DGtal::DigitalConvexity<TKSpace>::
306 isFullyConvex( const LatticePolytope& P ) const
307 {
308  auto S = insidePoints( P );
309  for ( Dimension k = 1; k < KSpace::dimension; ++ k )
310  {
311  auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
312  auto intersected_cells = makeCellCover( P, k, k );
313  if ( ( intersected_cells.nbCells() != touched_cells.nbCells() )
314  || ( ! intersected_cells.subset( touched_cells ) ) )
315  return false;
316  }
317  return true;
318 }
319 
320 //-----------------------------------------------------------------------------
321 template <typename TKSpace>
322 bool
323 DGtal::DigitalConvexity<TKSpace>::
324 isKSubconvex( const LatticePolytope& P, const CellGeometry& C, const Dimension k ) const
325 {
326  auto intersected_cells = makeCellCover( P, k, k );
327  return intersected_cells.subset( C );
328 }
329 //-----------------------------------------------------------------------------
330 template <typename TKSpace>
331 bool
332 DGtal::DigitalConvexity<TKSpace>::
333 isFullySubconvex( const LatticePolytope& P, const CellGeometry& C ) const
334 {
335  auto intersected_cells = makeCellCover( P, C.minCellDim(), C.maxCellDim() );
336  return intersected_cells.subset( C );
337 }
338 
339 //-----------------------------------------------------------------------------
340 template <typename TKSpace>
341 bool
342 DGtal::DigitalConvexity<TKSpace>::
343 isKConvex( const RationalPolytope& P, const Dimension k ) const
344 {
345  if ( k == 0 ) return true;
346  auto S = insidePoints( P );
347  auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
348  auto intersected_cells = makeCellCover( P, k, k );
349  return intersected_cells.nbCells() == touched_cells.nbCells()
350  && intersected_cells.subset( touched_cells );
351 }
352 
353 //-----------------------------------------------------------------------------
354 template <typename TKSpace>
355 bool
356 DGtal::DigitalConvexity<TKSpace>::
357 isFullyConvex( const RationalPolytope& P ) const
358 {
359  auto S = insidePoints( P );
360  for ( Dimension k = 1; k < KSpace::dimension; ++ k )
361  {
362  auto touched_cells = makeCellCover( S.begin(), S.end(), k, k );
363  auto intersected_cells = makeCellCover( P, k, k );
364  if ( ( intersected_cells.nbCells() != touched_cells.nbCells() )
365  || ( ! intersected_cells.subset( touched_cells ) ) )
366  return false;
367  }
368  return true;
369 }
370 
371 //-----------------------------------------------------------------------------
372 template <typename TKSpace>
373 bool
374 DGtal::DigitalConvexity<TKSpace>::
375 isKSubconvex( const RationalPolytope& P, const CellGeometry& C,
376  const Dimension k ) const
377 {
378  auto intersected_cells = makeCellCover( P, k, k );
379  return intersected_cells.subset( C );
380 }
381 //-----------------------------------------------------------------------------
382 template <typename TKSpace>
383 bool
384 DGtal::DigitalConvexity<TKSpace>::
385 isFullySubconvex( const RationalPolytope& P, const CellGeometry& C ) const
386 {
387  auto intersected_cells = makeCellCover( P, C.minCellDim(), C.maxCellDim() );
388  return intersected_cells.subset( C );
389 }
390 
391 ///////////////////////////////////////////////////////////////////////////////
392 // Interface - public :
393 
394 /**
395  * Writes/Displays the object on an output stream.
396  * @param out the output stream where the object is written.
397  */
398 template <typename TKSpace>
399 inline
400 void
401 DGtal::DigitalConvexity<TKSpace>::selfDisplay ( std::ostream & out ) const
402 {
403  out << "[DigitalConvexity]";
404 }
405 
406 /**
407  * Checks the validity/consistency of the object.
408  * @return 'true' if the object is valid, 'false' otherwise.
409  */
410 template <typename TKSpace>
411 inline
412 bool
413 DGtal::DigitalConvexity<TKSpace>::isValid() const
414 {
415  return true;
416 }
417 
418 
419 ///////////////////////////////////////////////////////////////////////////////
420 // Implementation of inline functions //
421 
422 //-----------------------------------------------------------------------------
423 template <typename TKSpace>
424 inline
425 std::ostream&
426 DGtal::operator<< ( std::ostream & out,
427  const DigitalConvexity<TKSpace> & object )
428 {
429  object.selfDisplay( out );
430  return out;
431 }
432 
433 // //
434 ///////////////////////////////////////////////////////////////////////////////