DGtal 1.4.0
Loading...
Searching...
No Matches
DigitalConvexity.ih
1/**
2 * This program is free software: you can redistribute it and/or modify
3 * it under the terms of the GNU Lesser General Public License as
4 * published by the Free Software Foundation, either version 3 of the
5 * License, or (at your option) any later version.
6 *
7 * This program is distributed in the hope that it will be useful,
8 * but WITHOUT ANY WARRANTY; without even the implied warranty of
9 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10 * GNU General Public License for more details.
11 *
12 * You should have received a copy of the GNU General Public License
13 * along with this program. If not, see <http://www.gnu.org/licenses/>.
14 *
15 **/
16
17/**
18 * @file DigitalConvexity.ih
19 * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20 * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
21 *
22 * @date 2020/01/31
23 *
24 * Implementation of inline methods defined in DigitalConvexity.h
25 *
26 * This file is part of the DGtal library.
27 */
28
29
30//////////////////////////////////////////////////////////////////////////////
31#include <cstdlib>
32#include "DGtal/geometry/volumes/ConvexityHelper.h"
33//////////////////////////////////////////////////////////////////////////////
34
35//-----------------------------------------------------------------------------
36template <typename TKSpace>
37DGtal::DigitalConvexity<TKSpace>::
38DigitalConvexity( Clone<KSpace> K, bool safe )
39 : myK( K ), mySafe( safe )
40{
41}
42//-----------------------------------------------------------------------------
43template <typename TKSpace>
44DGtal::DigitalConvexity<TKSpace>::
45DigitalConvexity( Point lo, Point hi, bool safe )
46 : mySafe( safe )
47{
48 myK.init( lo, hi, true );
49}
50
51//-----------------------------------------------------------------------------
52template <typename TKSpace>
53const TKSpace&
54DGtal::DigitalConvexity<TKSpace>::
55space() const
56{
57 return myK;
58}
59
60//-----------------------------------------------------------------------------
61template <typename TKSpace>
62template <typename PointIterator>
63typename DGtal::DigitalConvexity<TKSpace>::LatticePolytope
64DGtal::DigitalConvexity<TKSpace>::
65makeSimplex( PointIterator itB, PointIterator itE )
66{
67 return LatticePolytope( itB, itE );
68}
69
70//-----------------------------------------------------------------------------
71template <typename TKSpace>
72typename DGtal::DigitalConvexity<TKSpace>::LatticePolytope
73DGtal::DigitalConvexity<TKSpace>::
74makeSimplex( std::initializer_list<Point> l )
75{
76 return LatticePolytope( l );
77}
78
79//-----------------------------------------------------------------------------
80template <typename TKSpace>
81template <typename PointIterator>
82typename DGtal::DigitalConvexity<TKSpace>::RationalPolytope
83DGtal::DigitalConvexity<TKSpace>::
84makeRationalSimplex( Integer d, PointIterator itB, PointIterator itE )
85{
86 return RationalPolytope( d, itB, itE );
87}
88
89//-----------------------------------------------------------------------------
90template <typename TKSpace>
91typename DGtal::DigitalConvexity<TKSpace>::RationalPolytope
92DGtal::DigitalConvexity<TKSpace>::
93makeRationalSimplex( std::initializer_list<Point> l )
94{
95 return RationalPolytope( l );
96}
97
98//-----------------------------------------------------------------------------
99template <typename TKSpace>
100template <typename PointIterator>
101bool
102DGtal::DigitalConvexity<TKSpace>::
103isSimplexFullDimensional( PointIterator itB, PointIterator itE )
104{
105 typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
106 const Dimension d = KSpace::dimension;
107 std::vector<Point> pts( d+1 );
108 Dimension k = 0;
109 for ( ; itB != itE && k <= d; ++itB, ++k ) pts[ k ] = *itB;
110 // A simplex has exactly d+1 vertices.
111 if ( k != d+1 ) return false;
112 Matrix M;
113 for ( Dimension i = 0; i < d; ++i )
114 for ( Dimension j = 0; j < d; ++j )
115 M.setComponent( i, j, pts[ i ][ j ] - pts[ d ][ j ] );
116 // A simplex has its vectors linearly independent.
117 return M.determinant() != 0;
118}
119
120//-----------------------------------------------------------------------------
121template <typename TKSpace>
122bool
123DGtal::DigitalConvexity<TKSpace>::
124isSimplexFullDimensional( std::initializer_list<Point> l )
125{
126 return isSimplexFullDimensional( l.begin(), l.end() );
127}
128
129//-----------------------------------------------------------------------------
130template <typename TKSpace>
131template <typename PointIterator>
132typename DGtal::DigitalConvexity<TKSpace>::SimplexType
133DGtal::DigitalConvexity<TKSpace>::
134simplexType( PointIterator itB, PointIterator itE )
135{
136 typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
137 const Dimension d = KSpace::dimension;
138 std::vector<Point> pts( d+1 );
139 Dimension k = 0;
140 for ( ; itB != itE && k <= d; ++itB, ++k ) pts[ k ] = *itB;
141 // A simplex has exactly d+1 vertices.
142 if ( k != d+1 ) return SimplexType::INVALID;
143 Matrix M;
144 for ( Dimension i = 0; i < d; ++i )
145 for ( Dimension j = 0; j < d; ++j )
146 M.setComponent( i, j, pts[ i ][ j ] - pts[ d ][ j ] );
147 // A simplex has its vectors linearly independent.
148 auto V = M.determinant();
149 return (V == 0) ? SimplexType::DEGENERATED
150 : ( ((V == 1) || (V==-1)) ? SimplexType::UNITARY : SimplexType::COMMON );
151}
152
153//-----------------------------------------------------------------------------
154template <typename TKSpace>
155void
156DGtal::DigitalConvexity<TKSpace>::
157displaySimplex( std::ostream& out, std::initializer_list<Point> l )
158{
159 displaySimplex( out, l.begin(), l.end() );
160}
161
162//-----------------------------------------------------------------------------
163template <typename TKSpace>
164template <typename PointIterator>
165void
166DGtal::DigitalConvexity<TKSpace>::
167displaySimplex( std::ostream& out, PointIterator itB, PointIterator itE )
168{
169 typedef SimpleMatrix<Integer,dimension,dimension> Matrix;
170 const Dimension d = KSpace::dimension;
171 std::vector<Point> pts( d+1 );
172 Dimension k = 0;
173 for ( ; itB != itE && k <= d; ++itB, ++k ) pts[ k ] = *itB;
174 // A simplex has exactly d+1 vertices.
175 if ( k != d+1 ) { out << "[SPLX INVALID]"; return; }
176 Matrix M;
177 for ( Dimension i = 0; i < d; ++i )
178 for ( Dimension kk = 0; kk < d; ++kk )
179 M.setComponent( i, kk, pts[ i ][ kk ] - pts[ d ][ kk ] );
180 // A simplex has its vectors linearly independent.
181 auto V = M.determinant();
182 out << "[SPLX V=" << V;
183 for ( Dimension i = 0; i < d; ++i ) {
184 out << " (";
185 for ( Dimension j = 0; j < d; ++j )
186 out << " " << M( i, j );
187 out << " )";
188 }
189 out << " ]";
190}
191
192//-----------------------------------------------------------------------------
193template <typename TKSpace>
194typename DGtal::DigitalConvexity<TKSpace>::SimplexType
195DGtal::DigitalConvexity<TKSpace>::
196simplexType( std::initializer_list<Point> l )
197{
198 return simplexType( l.begin(), l.end() );
199}
200
201
202//-----------------------------------------------------------------------------
203template <typename TKSpace>
204typename DGtal::DigitalConvexity<TKSpace>::PointRange
205DGtal::DigitalConvexity<TKSpace>::
206insidePoints( const LatticePolytope& polytope )
207{
208 PointRange pts;
209 polytope.getPoints( pts );
210 return pts;
211}
212//-----------------------------------------------------------------------------
213template <typename TKSpace>
214typename DGtal::DigitalConvexity<TKSpace>::PointRange
215DGtal::DigitalConvexity<TKSpace>::
216interiorPoints( const LatticePolytope& polytope )
217{
218 PointRange pts;
219 polytope.getInteriorPoints( pts );
220 return pts;
221}
222
223//-----------------------------------------------------------------------------
224template <typename TKSpace>
225typename DGtal::DigitalConvexity<TKSpace>::PointRange
226DGtal::DigitalConvexity<TKSpace>::
227insidePoints( const RationalPolytope& polytope )
228{
229 PointRange pts;
230 polytope.getPoints( pts );
231 return pts;
232}
233//-----------------------------------------------------------------------------
234template <typename TKSpace>
235typename DGtal::DigitalConvexity<TKSpace>::PointRange
236DGtal::DigitalConvexity<TKSpace>::
237interiorPoints( const RationalPolytope& polytope )
238{
239 PointRange pts;
240 polytope.getInteriorPoints( pts );
241 return pts;
242}
243
244//-----------------------------------------------------------------------------
245template <typename TKSpace>
246template <typename PointIterator>
247typename DGtal::DigitalConvexity<TKSpace>::CellGeometry
248DGtal::DigitalConvexity<TKSpace>::
249makeCellCover( PointIterator itB, PointIterator itE,
250 Dimension i, Dimension k ) const
251{
252 ASSERT( i <= k );
253 ASSERT( k <= KSpace::dimension );
254 CellGeometry cgeom( myK, i, k, false );
255 cgeom.addCellsTouchingPoints( itB, itE );
256 return cgeom;
257}
258
259//-----------------------------------------------------------------------------
260template <typename TKSpace>
261typename DGtal::DigitalConvexity<TKSpace>::CellGeometry
262DGtal::DigitalConvexity<TKSpace>::
263makeCellCover( const LatticePolytope& P,
264 Dimension i, Dimension k ) const
265{
266 ASSERT( i <= k );
267 ASSERT( k <= KSpace::dimension );
268 CellGeometry cgeom( myK, i, k, false );
269 cgeom.addCellsTouchingPolytope( P );
270 return cgeom;
271}
272
273//-----------------------------------------------------------------------------
274template <typename TKSpace>
275typename DGtal::DigitalConvexity<TKSpace>::CellGeometry
276DGtal::DigitalConvexity<TKSpace>::
277makeCellCover( const RationalPolytope& P,
278 Dimension i, Dimension k ) const
279{
280 ASSERT( i <= k );
281 ASSERT( k <= KSpace::dimension );
282 CellGeometry cgeom( myK, i, k, false );
283 cgeom.addCellsTouchingPolytope( P );
284 return cgeom;
285}
286
287//-----------------------------------------------------------------------------
288template <typename TKSpace>
289typename DGtal::DigitalConvexity<TKSpace>::LatticePolytope
290DGtal::DigitalConvexity<TKSpace>::
291makePolytope( const PointRange& X, bool make_minkowski_summable ) const
292{
293 if ( mySafe )
294 {
295 typedef typename detail::ConvexityHelperInternalInteger< Integer, true >::Type
296 InternalInteger;
297 return ConvexityHelper< dimension, Integer, InternalInteger >::
298 computeLatticePolytope( X, false, make_minkowski_summable );
299 }
300 else
301 {
302 typedef typename detail::ConvexityHelperInternalInteger< Integer, false >::Type
303 InternalInteger;
304 return ConvexityHelper< dimension, Integer, InternalInteger >::
305 computeLatticePolytope( X, false, make_minkowski_summable );
306 }
307}
308
309//-----------------------------------------------------------------------------
310template <typename TKSpace>
311typename DGtal::DigitalConvexity<TKSpace>::PointRange
312DGtal::DigitalConvexity<TKSpace>::
313U( Dimension i, const PointRange& X ) const
314{
315 PointRange Y( X );
316 PointRange Z;
317 Z.reserve( X.size() );
318 auto add_one = [&i] ( Point & p ) { p[ i ] += 1; };
319 std::for_each( Y.begin(), Y.end(), add_one );
320 std::set_union( X.cbegin(), X.cend(), Y.cbegin(), Y.cend(),
321 std::back_inserter( Z ) );
322 return Z;
323}
324
325//-----------------------------------------------------------------------------
326template <typename TKSpace>
327bool
328DGtal::DigitalConvexity<TKSpace>::
329is0Convex( const PointRange& X ) const
330{
331 if ( X.empty() ) return true;
332 const auto P = makePolytope( X );
333 return P.count() == (DGtal::BoundedLatticePolytope<DGtal::SpaceND<3>>::Integer)X.size();
334}
335
336//-----------------------------------------------------------------------------
337template <typename TKSpace>
338bool
339DGtal::DigitalConvexity<TKSpace>::
340isFullyConvex( const PointRange& Z, bool convex0 ) const
341{
342 ASSERT( dimension <= 64 );
343 typedef DGtal::int64_t Direction;
344 typedef std::vector< Direction > Directions;
345 std::array< Directions, dimension > C;
346 const bool cvx0 = convex0 ? true : is0Convex( Z );
347 if ( ! cvx0 ) return false;
348 C[ 0 ].push_back( (Direction) 0 );
349 std::map< Direction, PointRange > X;
350 X[ 0 ] = Z;
351 std::sort( X[ 0 ].begin(), X[ 0 ].end() );
352 for ( Dimension k = 1; k < dimension; k++ )
353 {
354 for ( const auto beta : C[ k-1 ] )
355 {
356 for ( Dimension j = 0; j < dimension; j++ )
357 {
358 const Direction dir_j = Direction(1) << j;
359 if ( beta < dir_j )
360 {
361 const Direction alpha = beta | dir_j;
362 C[ k ].push_back( alpha );
363 X[ alpha ] = U( j, X[ beta ] );
364 if ( ! is0Convex( X[ alpha ] ) ) return false;
365 }
366 }
367 }
368 }
369 return true;
370}
371
372//-----------------------------------------------------------------------------
373template <typename TKSpace>
374bool
375DGtal::DigitalConvexity<TKSpace>::
376isFullyConvexFast( const PointRange& Z ) const
377{
378 LatticeSet C_Z( Z.cbegin(), Z.cend(), 0 );
379 const auto nb_cells = C_Z.starOfPoints().size();
380 const auto s = sizeStarCvxH( Z );
381 return s == (Integer)nb_cells;
382}
383
384//-----------------------------------------------------------------------------
385template <typename TKSpace>
386typename DGtal::DigitalConvexity<TKSpace>::PointRange
387DGtal::DigitalConvexity<TKSpace>::
388ExtrCvxH( const PointRange& X ) const
389{
390 if ( mySafe )
391 {
392 typedef typename detail::ConvexityHelperInternalInteger< Integer, true >::Type
393 InternalInteger;
394 return ConvexityHelper< dimension, Integer, InternalInteger >::
395 computeLatticePolytopeVertices( X, false, false );
396 }
397 else
398 {
399 typedef typename detail::ConvexityHelperInternalInteger< Integer, false >::Type
400 InternalInteger;
401 return ConvexityHelper< dimension, Integer, InternalInteger >::
402 computeLatticePolytopeVertices( X, false, false );
403 }
404}
405
406//-----------------------------------------------------------------------------
407template <typename TKSpace>
408typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
409DGtal::DigitalConvexity<TKSpace>::
410StarCvxH( const PointRange& X, Dimension axis ) const
411{
412 PointRange cells;
413 // Computes Minkowski sum of Z with hypercube
414 PointRange Z = U( 0, X );
415 for ( Dimension k = 1; k < dimension; k++ )
416 Z = U( k, Z );
417 // Builds polytope
418 const auto P = makePolytope( Z );
419 // Extracts lattice points within polytope
420 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
421 Counter C( P );
422 const Dimension a = axis >= dimension ? C.longestAxis() : axis;
423 auto cellP = C.getLatticeCells( a );
424 return LatticeSet( cellP, a );
425}
426
427//-----------------------------------------------------------------------------
428template <typename TKSpace>
429typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
430DGtal::DigitalConvexity<TKSpace>::
431StarCvxH( const Point& a, const Point& b, const Point& c,
432 Dimension axis ) const
433{
434 LatticeSet LS;
435 if ( mySafe )
436 {
437 using InternalInteger
438 = typename detail::ConvexityHelperInternalInteger< Integer, true >::Type;
439 using Helper = ConvexityHelper< dimension, Integer, InternalInteger >;
440 using UnitSegment = typename Helper::LatticePolytope::UnitSegment;
441 auto P = Helper::compute3DTriangle( a, b, c, true );
442 if ( ! P.isValid() ) return LS;
443 P += UnitSegment( 0 );
444 P += UnitSegment( 1 );
445 P += UnitSegment( 2 );
446 // Extracts lattice points within polytope
447 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
448 Counter C( P );
449 if ( axis >= dimension ) axis = C.longestAxis();
450 const auto cellP = C.getLatticeCells( axis );
451 return LatticeSet( cellP, axis );
452 }
453 else
454 {
455 using InternalInteger
456 = typename detail::ConvexityHelperInternalInteger< Integer, false >::Type;
457 using Helper = ConvexityHelper< dimension, Integer, InternalInteger >;
458 using UnitSegment = typename Helper::LatticePolytope::UnitSegment;
459 auto P = Helper::compute3DTriangle( a, b, c, true );
460 if ( ! P.isValid() ) return LS;
461 P += UnitSegment( 0 );
462 P += UnitSegment( 1 );
463 P += UnitSegment( 2 );
464 // Extracts lattice points within polytope
465 // they correspond 1-1 to the d-cells intersected by Cvxh( Z )
466 Counter C( P );
467 if ( axis >= dimension ) axis = C.longestAxis();
468 const auto cellP = C.getLatticeCells( axis );
469 return LatticeSet( cellP, axis );
470 }
471}
472
473//-----------------------------------------------------------------------------
474template <typename TKSpace>
475typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
476DGtal::DigitalConvexity<TKSpace>::
477Star( 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//-----------------------------------------------------------------------------
484template <typename TKSpace>
485typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
486DGtal::DigitalConvexity<TKSpace>::
487StarCells( 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//-----------------------------------------------------------------------------
495template <typename TKSpace>
496typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
497DGtal::DigitalConvexity<TKSpace>::
498toLatticeSet( 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//-----------------------------------------------------------------------------
505template <typename TKSpace>
506typename DGtal::DigitalConvexity<TKSpace>::PointRange
507DGtal::DigitalConvexity<TKSpace>::
508toPointRange( const LatticeSet& L ) const
509{
510 return L.toPointRange();
511}
512
513//-----------------------------------------------------------------------------
514template <typename TKSpace>
515typename DGtal::DigitalConvexity<TKSpace>::Integer
516DGtal::DigitalConvexity<TKSpace>::
517sizeStarCvxH( 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//-----------------------------------------------------------------------------
544template <typename TKSpace>
545typename DGtal::DigitalConvexity<TKSpace>::PointRange
546DGtal::DigitalConvexity<TKSpace>::
547Extr( 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//-----------------------------------------------------------------------------
572template <typename TKSpace>
573typename DGtal::DigitalConvexity<TKSpace>::PointRange
574DGtal::DigitalConvexity<TKSpace>::
575Extr( const LatticeSet& C ) const
576{
577 return C.extremaOfCells();
578}
579//-----------------------------------------------------------------------------
580template <typename TKSpace>
581typename DGtal::DigitalConvexity<TKSpace>::LatticeSet
582DGtal::DigitalConvexity<TKSpace>::
583Skel( const LatticeSet& C ) const
584{
585 return C.skeletonOfCells();
586}
587//-----------------------------------------------------------------------------
588template <typename TKSpace>
589typename DGtal::DigitalConvexity<TKSpace>::PointRange
590DGtal::DigitalConvexity<TKSpace>::
591ExtrSkel( const LatticeSet& C ) const
592{
593 return C.skeletonOfCells().extremaOfCells();
594}
595
596//-----------------------------------------------------------------------------
597template <typename TKSpace>
598typename DGtal::DigitalConvexity<TKSpace>::PointRange
599DGtal::DigitalConvexity<TKSpace>::
600FC_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//-----------------------------------------------------------------------------
751template <typename TKSpace>
752typename DGtal::DigitalConvexity<TKSpace>::PointRange
753DGtal::DigitalConvexity<TKSpace>::
754FC_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//-----------------------------------------------------------------------------
762template <typename TKSpace>
763typename DGtal::DigitalConvexity<TKSpace>::PointRange
764DGtal::DigitalConvexity<TKSpace>::
765FC( 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//-----------------------------------------------------------------------------
776template <typename TKSpace>
777typename DGtal::DigitalConvexity<TKSpace>::PointRange
778DGtal::DigitalConvexity<TKSpace>::
779envelope( 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//-----------------------------------------------------------------------------
795template <typename TKSpace>
796typename DGtal::DigitalConvexity<TKSpace>::PointRange
797DGtal::DigitalConvexity<TKSpace>::
798relativeEnvelope( 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//-----------------------------------------------------------------------------
816template <typename TKSpace>
817template <typename Predicate>
818typename DGtal::DigitalConvexity<TKSpace>::PointRange
819DGtal::DigitalConvexity<TKSpace>::
820relativeEnvelope( 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//-----------------------------------------------------------------------------
835template <typename TKSpace>
836typename DGtal::DigitalConvexity<TKSpace>::Size
837DGtal::DigitalConvexity<TKSpace>::
838depthLastEnvelope() const
839{
840 return myDepthLastFCE;
841}
842
843//-----------------------------------------------------------------------------
844template <typename TKSpace>
845bool
846DGtal::DigitalConvexity<TKSpace>::
847isKConvex( 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//-----------------------------------------------------------------------------
859template <typename TKSpace>
860bool
861DGtal::DigitalConvexity<TKSpace>::
862isFullyConvex( 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//-----------------------------------------------------------------------------
878template <typename TKSpace>
879bool
880DGtal::DigitalConvexity<TKSpace>::
881isKSubconvex( 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//-----------------------------------------------------------------------------
887template <typename TKSpace>
888bool
889DGtal::DigitalConvexity<TKSpace>::
890isFullySubconvex( 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//-----------------------------------------------------------------------------
897template <typename TKSpace>
898bool
899DGtal::DigitalConvexity<TKSpace>::
900isFullySubconvex( 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//-----------------------------------------------------------------------------
912template <typename TKSpace>
913bool
914DGtal::DigitalConvexity<TKSpace>::
915isFullySubconvex( const PointRange& Y, const LatticeSet& StarX ) const
916{
917 const auto SCY = StarCvxH( Y, StarX.axis() );
918 return StarX.includes( SCY );
919}
920
921//-----------------------------------------------------------------------------
922template <typename TKSpace>
923bool
924DGtal::DigitalConvexity<TKSpace>::
925isFullySubconvex( 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//-----------------------------------------------------------------------------
934template <typename TKSpace>
935bool
936DGtal::DigitalConvexity<TKSpace>::
937isKSubconvex( 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//-----------------------------------------------------------------------------
946template <typename TKSpace>
947bool
948DGtal::DigitalConvexity<TKSpace>::
949isFullySubconvex( 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//-----------------------------------------------------------------------------
958template <typename TKSpace>
959bool
960DGtal::DigitalConvexity<TKSpace>::
961isFullySubconvex( 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//-----------------------------------------------------------------------------
997template <typename TKSpace>
998bool
999DGtal::DigitalConvexity<TKSpace>::
1000isKConvex( 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//-----------------------------------------------------------------------------
1011template <typename TKSpace>
1012bool
1013DGtal::DigitalConvexity<TKSpace>::
1014isFullyConvex( 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//-----------------------------------------------------------------------------
1029template <typename TKSpace>
1030bool
1031DGtal::DigitalConvexity<TKSpace>::
1032isKSubconvex( 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//-----------------------------------------------------------------------------
1039template <typename TKSpace>
1040bool
1041DGtal::DigitalConvexity<TKSpace>::
1042isFullySubconvex( 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//-----------------------------------------------------------------------------
1049template <typename TKSpace>
1050void
1051DGtal::DigitalConvexity<TKSpace>::
1052eraseInterval( 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//-----------------------------------------------------------------------------
1092template <typename TKSpace>
1093template <typename Predicate>
1094typename DGtal::DigitalConvexity<TKSpace>::PointRange
1095DGtal::DigitalConvexity<TKSpace>::
1096filter( 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 */
1112template <typename TKSpace>
1113inline
1114void
1115DGtal::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 */
1124template <typename TKSpace>
1125inline
1126bool
1127DGtal::DigitalConvexity<TKSpace>::isValid() const
1128{
1129 return true;
1130}
1131
1132
1133///////////////////////////////////////////////////////////////////////////////
1134// Implementation of inline functions //
1135
1136//-----------------------------------------------------------------------------
1137template <typename TKSpace>
1138inline
1139std::ostream&
1140DGtal::operator<< ( std::ostream & out,
1141 const DigitalConvexity<TKSpace> & object )
1142{
1143 object.selfDisplay( out );
1144 return out;
1145}
1146
1147// //
1148///////////////////////////////////////////////////////////////////////////////