diff --git a/src/DGtal/geometry/doc/moduleDigitalConvexity.dox b/src/DGtal/geometry/doc/moduleDigitalConvexity.dox index 20567d66b4..de8107110f 100644 --- a/src/DGtal/geometry/doc/moduleDigitalConvexity.dox +++ b/src/DGtal/geometry/doc/moduleDigitalConvexity.dox @@ -361,7 +361,8 @@ Convexity services: - DigitalConvexity::isFullySubconvex( const LatticePolytope& P, const LatticeSet& StarX ) const tells if a given polytope is fully subconvex to a set of cells, represented as a LatticeSetByIntervals. - DigitalConvexity::isFullySubconvex( const PointRange& Y, const LatticeSet& StarX ) const tells is the given range of points is fully subconvex to a set of cells, represented as a LatticeSetByIntervals. - DigitalConvexity::isFullySubconvex( const Point& a, const Point& b, const LatticeSet& StarX ) const tells if the given pair of points is fully subconvex to a set of cells, represented as a LatticeSetByIntervals (three times faster than previous generic method). -- DigitalConvexity::isFullySubconvex( const Point& a, const Point& b, const Point& c, const LatticeSet& StarX ) const tells if the given triplets of points is fully subconvex to a set of cells, represented as a LatticeSetByIntervals (five times faster than previous generic method). +- DigitalConvexity::isFullySubconvex( const Point& a, const Point& b, const Point& c, const LatticeSet& cells ) const tells if the given triplets of points is fully subconvex to a set of cells, represented as a LatticeSetByIntervals (five times faster than previous generic method). +- DigitalConvexity::isFullyCovered( const Point& a, const Point& b, const Point& c, const LatticeSet& cells ) const tells if the given triplets of points is fully covered by the given set of cells, represented as a LatticeSetByIntervals (limited to 3D at the moment). @subsection dgtal_dconvexity_sec25 Ehrhart polynomial of a lattice polytope diff --git a/src/DGtal/geometry/doc/moduleEnvelope.dox b/src/DGtal/geometry/doc/moduleEnvelope.dox index f9c95ebc37..2f1d758d30 100644 --- a/src/DGtal/geometry/doc/moduleEnvelope.dox +++ b/src/DGtal/geometry/doc/moduleEnvelope.dox @@ -57,9 +57,10 @@ operator defined as: where: - \f$ \mathrm{CvxH} \f$ denotes the standard convex hull in Euclidean space, -- \f$ \mathrm{Star} \f$ denotes the set of cells of the cubical grid complex induced by lattice points, whose closure intersects the given set, +- \f$ \mathrm{Star} \f$ denotes the set of cells of the cubical grid complex whose closure intersects the given Euclidean set, - \f$ \mathrm{Skel} \f$ denotes the skeleton of a cell complex, i.e. the smallest cell complex whose star is equal to the input cell complex, -- \f$ \mathrm{Extr} \f$ returns the set of lattice points that are vertices of cells of the given cell complex. +- \f$ \mathrm{Extr} \f$ returns the set of lattice points that are vertices of cells of the given cell complex, +- \f$ \mathrm{Cover} \f$ denotes the set of cells of the cubical grid complex which have a non empty intersection with the given Euclidean set. The envelope operator is the limit of the iterated repetition of operator \f$ \mathrm{FC} \f$. It is proven that this sequence achieves @@ -105,6 +106,7 @@ intermediate steps, which can be useful in some contexts. - DigitalConvexity::Extr computes the extremal vertices of a range of cells. - DigitalConvexity::ExtrCvxH computes the extremal vertices of the convex hull of a range of lattice points. - DigitalConvexity::ExtrSkel computes the extremal vertices of the skeleton of a range of cells. +- DigitalConvexity::CoverCvxH computes the set of lattice cells \f$ \mathrm{Cover}(\mathrm{CvxH}(a,b,c)) \f$, i.e. the cells that have a non-empty intersection with the convex hull of the given triangle (limited to 3D at the moment). \b Representations and \b conversions diff --git a/src/DGtal/geometry/volumes/BoundedLatticePolytope.h b/src/DGtal/geometry/volumes/BoundedLatticePolytope.h index 902c1ef39d..ccecba6f3d 100644 --- a/src/DGtal/geometry/volumes/BoundedLatticePolytope.h +++ b/src/DGtal/geometry/volumes/BoundedLatticePolytope.h @@ -100,6 +100,15 @@ namespace DGtal UnitSegment( Dimension d ) : k( d ) {} }; + /** + * Represents the unit segment from (0,...,0) (excluded) to + * (0,...,1,...,0) (excluded) with the 1 at position \a k. + */ + struct StrictUnitSegment { + Dimension k; + StrictUnitSegment( Dimension d ) : k( d ) {} + }; + /** * Represents the unit segment from (0,...,0) (included) to * (0,...,1,...,0) (excluded) with the 1 at position \a k. @@ -199,6 +208,33 @@ namespace DGtal } }; + /** + * Represents the unit cell obtained by successive Minkowski sum + * of StrictUnitSegment whose dimensions are stored in \a dims. When \a + * dims is empty, it is only the point (0,...,0). + */ + struct StrictUnitCell { + std::vector dims; + StrictUnitCell( std::initializer_list l ) + : dims( l.begin(), l.end() ) {} + + /** + * Overloads 'operator<<' for displaying objects of class 'BoundedLatticePolytope::UnitCell'. + * @param out the output stream where the object is written. + * @param object the object of class 'BoundedLatticePolytope::UnitCell' to write. + * @return the output stream after the writing. + */ + friend std::ostream& + operator<< ( std::ostream & out, + const StrictUnitCell & object ) + { + out << "{"; + for ( Dimension i = 0; i < object.dims.size(); ++i ) out << object.dims[ i ]; + out << "}"; + return out; + } + }; + /// @name Standard services (construction, initialization, assignment, interior, closure) /// @{ @@ -356,6 +392,18 @@ namespace DGtal /// form `Ax <= b`, 'false' if it is of the form `Ax < b`. bool isLarge( unsigned int i ) const; + /// Sets the \a i-th half space constraint to the form `Ax <= b`. + /// + /// @param i the index of the half-space constraint between 0 and + /// `nbHalfSpaces()` (excluded). + void setLarge( unsigned int i ); + + /// Sets the \a i-th half space constraint to the form `Ax < b`. + /// + /// @param i the index of the half-space constraint between 0 and + /// `nbHalfSpaces()` (excluded). + void setStrict( unsigned int i ); + /// @return the matrix A in the polytope representation \f$ Ax \le B \f$. const InequalityMatrix& getA() const; @@ -504,10 +552,26 @@ namespace DGtal * @param s any strict unit segment. * @return a reference to 'this'. */ + Self& operator+=( StrictUnitSegment s ); + + /** + * Minkowski sum of this polytope with an axis-aligned strict unit cell. + * + * @param c any unit cell. + * @return a reference to 'this'. + */ + Self& operator+=( StrictUnitCell c ); + + /** + * Minkowski sum of this polytope with a right strict unit segment aligned with some axis. + * + * @param s any strict unit segment. + * @return a reference to 'this'. + */ Self& operator+=( RightStrictUnitSegment s ); /** - * Minkowski sum of this polytope with an axis-aligned strict unit cell. + * Minkowski sum of this polytope with an axis-aligned right strict unit cell. * * @param c any strict unit cell. * @return a reference to 'this'. @@ -515,7 +579,7 @@ namespace DGtal Self& operator+=( RightStrictUnitCell c ); /** - * Minkowski sum of this polytope with a strict unit segment aligned with some axis. + * Minkowski sum of this polytope with a left strict unit segment aligned with some axis. * * @param s any strict unit segment. * @return a reference to 'this'. @@ -523,7 +587,7 @@ namespace DGtal Self& operator+=( LeftStrictUnitSegment s ); /** - * Minkowski sum of this polytope with an axis-aligned strict unit cell. + * Minkowski sum of this polytope with an axis-aligned left strict unit cell. * * @param c any strict unit cell. * @return a reference to 'this'. @@ -1029,6 +1093,30 @@ namespace DGtal operator+ ( const BoundedLatticePolytope & P, typename BoundedLatticePolytope::UnitCell c ); + /** + * Minkowski sum of polytope \a P with strict unit segment \a s aligned with some axis. + * + * @param P any polytope. + * @param s any strict unit segment. + * @return the Polytope P + s. + */ + template + BoundedLatticePolytope + operator+ ( const BoundedLatticePolytope & P, + typename BoundedLatticePolytope::StrictUnitSegment s ); + + /** + * Minkowski sum of polytope \a P with an axis-aligned strict unit cell \a c. + * + * @param P any polytope. + * @param c any strict unit cell. + * @return the Polytope P + c. + */ + template + BoundedLatticePolytope + operator+ ( const BoundedLatticePolytope & P, + typename BoundedLatticePolytope::StrictUnitCell c ); + /** * Minkowski sum of polytope \a P with strict unit segment \a s aligned with some axis. * diff --git a/src/DGtal/geometry/volumes/BoundedLatticePolytope.ih b/src/DGtal/geometry/volumes/BoundedLatticePolytope.ih index 1787964fc3..3df4348ceb 100644 --- a/src/DGtal/geometry/volumes/BoundedLatticePolytope.ih +++ b/src/DGtal/geometry/volumes/BoundedLatticePolytope.ih @@ -503,6 +503,28 @@ operator+=( UnitSegment s ) return *this; } +//----------------------------------------------------------------------------- +template +typename DGtal::BoundedLatticePolytope::Self& +DGtal::BoundedLatticePolytope:: +operator+=( StrictUnitSegment s ) +{ + for ( Dimension i = 0; i < A.size(); ++i ) + { + if ( A[ i ][ s.k ] > NumberTraits::ZERO ) + { + B[ i ] += A[ i ][ s.k ]; + I[ i ] = false; + } + else if ( A[ i ][ s.k ] < NumberTraits::ZERO ) + I[ i ] = false; + } + Vector z = Vector::zero; + z[ s.k ] = NumberTraits::ONE; + D = Domain( D.lowerBound(), D.upperBound() + z ); + return *this; +} + //----------------------------------------------------------------------------- template typename DGtal::BoundedLatticePolytope::Self& @@ -553,6 +575,17 @@ operator+=( UnitCell c ) return *this; } +//----------------------------------------------------------------------------- +template +typename DGtal::BoundedLatticePolytope::Self& +DGtal::BoundedLatticePolytope:: +operator+=( StrictUnitCell c ) +{ + for ( Dimension i = 0; i < c.dims.size(); ++i ) + *this += StrictUnitSegment( c.dims[ i ] ); + return *this; +} + //----------------------------------------------------------------------------- template typename DGtal::BoundedLatticePolytope::Self& @@ -927,6 +960,25 @@ DGtal::BoundedLatticePolytope::isLarge( unsigned int i ) const return I[ i ]; } +//----------------------------------------------------------------------------- +template +void +DGtal::BoundedLatticePolytope::setLarge( unsigned int i ) +{ + ASSERT( i < nbHalfSpaces() ); + I[ i ] = true; +} + +//----------------------------------------------------------------------------- +template +void +DGtal::BoundedLatticePolytope::setStrict( unsigned int i ) +{ + ASSERT( i < nbHalfSpaces() ); + I[ i ] = false; +} + + //----------------------------------------------------------------------------- template const typename DGtal::BoundedLatticePolytope::InequalityMatrix& @@ -979,7 +1031,7 @@ DGtal::BoundedLatticePolytope::selfDisplay ( std::ostream & out ) const out << " ["; for ( Dimension j = 0; j < dimension; ++j ) out << " " << A[ i ][ j ]; - out << " ] . x <= " << B[ i ] << std::endl; + out << " ] . x " << ( isLarge( i ) ? "<=" : "<" ) << " " << B[ i ] << std::endl; } } @@ -1052,6 +1104,26 @@ DGtal::operator+ ( const BoundedLatticePolytope & P, //----------------------------------------------------------------------------- template DGtal::BoundedLatticePolytope +DGtal::operator+ ( const BoundedLatticePolytope & P, + typename BoundedLatticePolytope::StrictUnitSegment s ) +{ + BoundedLatticePolytope Q = P; + Q += s; + return Q; +} +//----------------------------------------------------------------------------- +template +DGtal::BoundedLatticePolytope +DGtal::operator+ ( const BoundedLatticePolytope & P, + typename BoundedLatticePolytope::StrictUnitCell c ) +{ + BoundedLatticePolytope Q = P; + Q += c; + return Q; +} +//----------------------------------------------------------------------------- +template +DGtal::BoundedLatticePolytope DGtal::operator+ ( const BoundedLatticePolytope & P, typename BoundedLatticePolytope::RightStrictUnitSegment s ) { diff --git a/src/DGtal/geometry/volumes/BoundedLatticePolytopeCounter.ih b/src/DGtal/geometry/volumes/BoundedLatticePolytopeCounter.ih index f3ab0e1a4f..4180b94b5e 100644 --- a/src/DGtal/geometry/volumes/BoundedLatticePolytopeCounter.ih +++ b/src/DGtal/geometry/volumes/BoundedLatticePolytopeCounter.ih @@ -78,7 +78,7 @@ intersectionIntervalAlongAxis( Point p, Dimension a ) const const Integer x_a = x_min; p[ a ] = x_a; bool empty = false; - for ( Dimension k = 2*dimension; k < A.size(); k++ ) + for ( Dimension k = 0 /* 2*dimension */; k < A.size(); k++ ) { const Integer c = A[ k ].dot( p ); const Integer n = A[ k ][ a ]; @@ -103,6 +103,7 @@ intersectionIntervalAlongAxis( Point p, Dimension a ) const if ( d >= 0 ) { x = I[ k ] ? ( (d-n-1) / -n ) : ( d / -n + 1 ); + // x = I[ k ] ? ( (d-n-1) / -n ) : ( (d-n-1) / -n + 1 ); ? x_min = std::max( x_min, x_a + x ); } // otherwise the constraint is true @@ -267,9 +268,11 @@ getLatticeSet( Dimension a ) const LatticeSetByInterval L; for ( auto&& p : D ) { - auto I = intersectionIntervalAlongAxis( p, a ); - L[ p ] = I; + auto I = intersectionIntervalAlongAxis( p, a ); + if ( I.second > I.first ) + L[ p ] = Interval( I.first, I.second - 1 ); } + return L; } //----------------------------------------------------------------------------- diff --git a/src/DGtal/geometry/volumes/ConvexityHelper.h b/src/DGtal/geometry/volumes/ConvexityHelper.h index c54a35e9b4..bc482508fc 100644 --- a/src/DGtal/geometry/volumes/ConvexityHelper.h +++ b/src/DGtal/geometry/volumes/ConvexityHelper.h @@ -394,12 +394,37 @@ namespace DGtal /// moduleDigitalConvexity). /// /// @return the tightiest bounded lattice polytope - /// (i.e. H-representation) including the given range of points. + /// (i.e. H-representation) including the given triangle. + /// + /// @warning Implemented only in 3D. static LatticePolytope compute3DTriangle( const Point& a, const Point& b, const Point& c, bool make_minkowski_summable = false ); + /// Computes the lattice polytope enclosing an open triangle in + /// dimension 3. Takes care of degeneracies (non distinct points + /// or alignment). + /// + /// @param a any point + /// @param b any point + /// @param c any point + /// + /// @param[in] make_minkowski_summable Other constraints are added + /// so that we can perform axis aligned Minkowski sums on this + /// polytope. Useful for checking full convexity (see + /// moduleDigitalConvexity). + /// + /// @return the tightiest bounded lattice polytope + /// (i.e. H-representation) including the given open triangle + /// (i.e. without edges and vertices). + /// + /// @warning Implemented only in 3D. + static + LatticePolytope + compute3DOpenTriangle( const Point& a, const Point& b, const Point& c, + bool make_minkowski_summable = false ); + /// Computes the lattice polytope enclosing a degenerated /// triangle. The points must be aligned (or non distinct). /// @@ -413,17 +438,34 @@ namespace DGtal LatticePolytope computeDegeneratedTriangle( const Point& a, const Point& b, const Point& c ); + /// Computes the lattice polytope enclosing a segment. /// /// @param a any point /// @param b any point /// /// @return the tightiest bounded lattice polytope - /// (i.e. H-representation) including the given range of - /// points. It is always Minkowski summable. + /// (i.e. H-representation) including the closed segment + /// `[a,b]`. It is always Minkowski summable. + /// + /// @warning Implemented only in 3D. static LatticePolytope computeSegment( const Point& a, const Point& b ); + + /// Computes the lattice polytope enclosing an open segment. + /// + /// @param a any point + /// @param b any point + /// + /// @return the tightiest bounded lattice polytope + /// (i.e. H-representation) including the open segment + /// `]a,b[`. It is always Minkowski summable. + /// + /// @warning Implemented only in 3D. + static + LatticePolytope + computeOpenSegment( const Point& a, const Point& b ); /// @} diff --git a/src/DGtal/geometry/volumes/ConvexityHelper.ih b/src/DGtal/geometry/volumes/ConvexityHelper.ih index 656a3a7b9e..cd06510c8e 100644 --- a/src/DGtal/geometry/volumes/ConvexityHelper.ih +++ b/src/DGtal/geometry/volumes/ConvexityHelper.ih @@ -244,6 +244,30 @@ DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::compute3DTriangle make_minkowski_summable, false ); } +//----------------------------------------------------------------------------- +template < int dim, typename TInteger, typename TInternalInteger > +typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope +DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::compute3DOpenTriangle +( const Point& a, const Point& b, const Point& c, + bool make_minkowski_summable ) +{ + using Op = detail::BoundedRationalPolytopeSpecializer< dimension, Integer>; + const Vector ab = b - a; + const Vector bc = c - b; + const Vector n = Op::crossProduct( ab, bc ); + if ( n == Vector::zero ) + return LatticePolytope(); // empty polytope + auto P = compute3DTriangle( a, b, c, make_minkowski_summable ); + // Compute domain + for ( auto k = 0; k < P.nbHalfSpaces(); k++ ) + { + if ( Op::crossProduct( P.getA( k ), n ) != Vector::zero ) + P.setStrict( k ); + } + return P; +} + + //----------------------------------------------------------------------------- template < int dim, typename TInteger, typename TInternalInteger > typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope @@ -269,6 +293,59 @@ computeDegeneratedTriangle return computeSegment( a, a ); } + +//----------------------------------------------------------------------------- +template < int dim, typename TInteger, typename TInternalInteger > +typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope +DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeOpenSegment +( const Point& a, const Point& b ) +{ + + if constexpr( dim != 3 ) return LatticePolytope(); + using Op = detail::BoundedRationalPolytopeSpecializer< dimension, Integer>; + typedef typename LatticePolytope::Domain Domain; + typedef typename LatticePolytope::HalfSpace PolytopeHalfSpace; + // Compute domain + const Point low = a.inf( b ); + const Point high = a.sup( b ); + const Vector ab = b - a; + bool degenerate = ( ab == Vector::zero ); + if ( degenerate ) return LatticePolytope(); + // Initialize polytope + std::vector< PolytopeHalfSpace > PHS; + PHS.reserve( 4*dim ); + // Compute additionnal constraints on domain boundary to make it open. + for ( Dimension k = 0; k < dim; k++ ) + { + const Vector bp = Vector::base( k, 1 ); + PHS.emplace_back( bp, high[ k ] ); + const Vector bn = Vector::base( k, -1 ); + PHS.emplace_back( bn, -low[ k ] ); + } + // Compute additionnal constraints on edges so that the + // Minkowski sum with axis-aligned edges is valid. + for ( Integer d = -1; d <= 1; d += 2 ) + for ( Dimension k = 0; k < dim; k++ ) + { + const Vector i = Vector::base( k, d ); + const Vector e = Op::crossProduct( ab, i ); + if ( e != Vector::zero ) + { + const Integer e_a = e.dot( a ); + PHS.emplace_back( e, e_a ); + } + } + auto P = LatticePolytope( Domain( low, high ), PHS.cbegin(), PHS.cend(), true, false ); + // Fix < inequalities + for ( auto k = 2 * dim; k < 4 * dim; k++ ) + { + auto V = P.getA( k ); + if ( V.dot( a ) != V.dot( b ) ) P.setStrict( k ); + } + return P; + +} + //----------------------------------------------------------------------------- template < int dim, typename TInteger, typename TInternalInteger > typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope diff --git a/src/DGtal/geometry/volumes/DigitalConvexity.h b/src/DGtal/geometry/volumes/DigitalConvexity.h index 5338fd1ac3..5a57b9517a 100644 --- a/src/DGtal/geometry/volumes/DigitalConvexity.h +++ b/src/DGtal/geometry/volumes/DigitalConvexity.h @@ -455,6 +455,41 @@ namespace DGtal bool isFullySubconvex( const Point& a, const Point& b, const Point& c, const LatticeSet& StarX ) const; + /// Tells if the non-degenerated 3D triangle a,b,c is + /// fully covered by some lattice set of cells \a cells. + /// + /// @param a any 3D point (distinct from the two others) + /// @param b any 3D point (distinct from the two others) + /// @param c any 3D point (distinct from the two others) + /// + /// @param cells any lattice set representing a set of cells. + /// + /// @return 'true' iff `Cvxh({a,b,c})` is fully covered by \a cells, + /// i.e. all the cells intersected by the triangle belong to \a + /// cells. + /// + /// @note This method is limited to 3D triangles. + bool isFullyCovered( const Point& a, const Point& b, const Point& c, + const LatticeSet& cells ) const; + + + /// Tells if the non-degenerated 3D open triangle a,b,c (ie + /// without its edges and vertices) is digitally fully subconvex + /// to some lattice set \a Star_X, i.e. the cell cover of some set + /// X represented by lattice points. + /// + /// @param a any 3D point (distinct from the two others) + /// @param b any 3D point (distinct from the two others) + /// @param c any 3D point (distinct from the two others) + /// + /// @param StarX any lattice set representing an open cubical complex. + /// @return 'true' iff Y is digitally fully subconvex to X. + /// + /// @note This method is supposed to be faster than the others, + /// but is limited to 3D open triangles. + bool isOpenTriangleFullySubconvex( const Point& a, const Point& b, const Point& c, + const LatticeSet& StarX ) const; + /// Tells if a given segment from \a a to \a b is digitally /// k-subconvex (i.e. k-tangent) to some cell cover \a C. The @@ -578,6 +613,32 @@ namespace DGtal /// must indeed have the same axis. LatticeSet StarCvxH( const Point& a, const Point& b, const Point& c, Dimension axis = dimension ) const; + + /// Builds the cell complex `Star(OpenTriangle({a,b,c}))` for + /// `a,b,c` the vertices of a non-degenerate open 3D triangle (ie + /// without edges or vertices), represented as a lattice set + /// (stacked row representation). + /// + /// @param a any 3D point (distinct from the two others) + /// @param b any 3D point (distinct from the two others) + /// @param c any 3D point (distinct from the two others) + /// + /// @param axis specifies the projection axis for the row + /// representation if below space dimension, otherwise chooses the + /// axis that minimizes memory/computations. + /// + /// @return the range of cells touching the open triangle `abc`, + /// represented as a lattice set (cells are represented with + /// Khalimsky coordinates). If the triangle is degenerate (a,b,c + /// not distinct or aligned), then it returns an empty range of + /// cells. + /// + /// @note It is useful to specify an axis if you wish later to + /// compare or make operations with several lattice sets. They + /// must indeed have the same axis. + LatticeSet StarOpenTriangle( const Point& a, const Point& b, const Point& c, + Dimension axis = dimension ) const; + /// Computes the number of cells in Star(CvxH(X)) for X a digital set. /// @@ -646,6 +707,35 @@ namespace DGtal /// @return the range of digital points that are the extremal /// vertices to the skeleton of the cells in \a C. PointRange ExtrSkel( const LatticeSet& C ) const; + + /// Builds the cell complex `Cover(CvxH({a,b,c}))` for `a,b,c` a + /// non-degenerate 3D triangle, represented as a lattice set + /// (stacked row representation). The cover of a Euclidean set X + /// is the set of grid cells that have a non-empty intersection + /// with X. It is always included in the star of X, which is the + /// of cells whose closure has a non-empty intersection with + /// X. However the cover of X is not necessarily closed or open. + /// + /// @param a any 3D point (distinct from the two others) + /// @param b any 3D point (distinct from the two others) + /// @param c any 3D point (distinct from the two others) + /// + /// @param axis specifies the projection axis for the row + /// representation if below space dimension, otherwise chooses the + /// axis that minimizes memory/computations. + /// + /// @return the range of cells intersecting the triangle `abc`, + /// represented as a lattice set (cells are represented with + /// Khalimsky coordinates). If the triangle is degenerate (a,b,c + /// not distinct or aligned), then it returns an empty range of + /// cells. + /// + /// @note It is useful to specify an axis if you wish later to + /// compare or make operations with several lattice sets. They + /// must indeed have the same axis. + LatticeSet CoverCvxH( const Point& a, const Point& b, const Point& c, + Dimension axis = dimension ) const; + /// Builds the lattice set (stacked row representation) associated /// to the given range of points. diff --git a/src/DGtal/geometry/volumes/DigitalConvexity.ih b/src/DGtal/geometry/volumes/DigitalConvexity.ih index ce71e3a5d9..0de363aaac 100644 --- a/src/DGtal/geometry/volumes/DigitalConvexity.ih +++ b/src/DGtal/geometry/volumes/DigitalConvexity.ih @@ -470,6 +470,52 @@ StarCvxH( const Point& a, const Point& b, const Point& c, } } +//----------------------------------------------------------------------------- +template +typename DGtal::DigitalConvexity::LatticeSet +DGtal::DigitalConvexity:: +StarOpenTriangle( const Point& a, const Point& b, const Point& c, + Dimension axis ) const +{ + LatticeSet LS; + if ( mySafe ) + { + using InternalInteger + = typename detail::ConvexityHelperInternalInteger< Integer, true >::Type; + using Helper = ConvexityHelper< dimension, Integer, InternalInteger >; + using UnitSegment = typename Helper::LatticePolytope::UnitSegment; + auto P = Helper::compute3DOpenTriangle( a, b, c, true ); + if ( ! P.isValid() ) return LS; + P += UnitSegment( 0 ); + P += UnitSegment( 1 ); + P += UnitSegment( 2 ); + // Extracts lattice points within polytope + // they correspond 1-1 to the d-cells intersected by Cvxh( Z ) + Counter C( P ); + if ( axis >= dimension ) axis = C.longestAxis(); + const auto cellP = C.getLatticeCells( axis ); + return LatticeSet( cellP, axis ); + } + else + { + using InternalInteger + = typename detail::ConvexityHelperInternalInteger< Integer, false >::Type; + using Helper = ConvexityHelper< dimension, Integer, InternalInteger >; + using UnitSegment = typename Helper::LatticePolytope::UnitSegment; + auto P = Helper::compute3DOpenTriangle( a, b, c, true ); + if ( ! P.isValid() ) return LS; + P += UnitSegment( 0 ); + P += UnitSegment( 1 ); + P += UnitSegment( 2 ); + // Extracts lattice points within polytope + // they correspond 1-1 to the d-cells intersected by Cvxh( Z ) + Counter C( P ); + if ( axis >= dimension ) axis = C.longestAxis(); + const auto cellP = C.getLatticeCells( axis ); + return LatticeSet( cellP, axis ); + } +} + //----------------------------------------------------------------------------- template typename DGtal::DigitalConvexity::LatticeSet @@ -491,6 +537,105 @@ StarCells( const PointRange& C, const Dimension axis ) const return L.starOfCells(); } +//----------------------------------------------------------------------------- +template +typename DGtal::DigitalConvexity::LatticeSet +DGtal::DigitalConvexity:: +CoverCvxH( const Point& a, const Point& b, const Point& c, + Dimension axis ) const +{ + LatticeSet LS; + if ( mySafe ) + { + using InternalInteger + = typename detail::ConvexityHelperInternalInteger< Integer, true >::Type; + using Helper = ConvexityHelper< dimension, Integer, InternalInteger >; + using StrictUnitSegment = typename Helper::LatticePolytope::StrictUnitSegment; + auto P = Helper::compute3DTriangle( a, b, c, true ); + if ( ! P.isValid() ) return LS; + Counter C( P ); + if ( axis >= dimension ) axis = C.longestAxis(); + auto V = LatticeSet( C.getLatticeSet( axis ), axis ).toPointRange(); // << 0-cells + auto P0 = P + StrictUnitSegment( 0 ); + auto P1 = P + StrictUnitSegment( 1 ); + auto P2 = P + StrictUnitSegment( 2 ); + Counter C0( P0 ); + Counter C1( P1 ); + Counter C2( P2 ); + auto V0 = LatticeSet( C0.getLatticeSet( axis ), axis ).toPointRange(); // << 1-cells + auto V1 = LatticeSet( C1.getLatticeSet( axis ), axis ).toPointRange(); // << 1-cells + auto V2 = LatticeSet( C2.getLatticeSet( axis ), axis ).toPointRange(); // << 1-cells + auto P01 = P0 + StrictUnitSegment( 1 ); + auto P02 = P0 + StrictUnitSegment( 2 ); + auto P12 = P1 + StrictUnitSegment( 2 ); + Counter C01( P01 ); + Counter C02( P02 ); + Counter C12( P12 ); + auto V01 = LatticeSet( C01.getLatticeSet( axis ), axis ).toPointRange(); // << 2-cells + auto V02 = LatticeSet( C02.getLatticeSet( axis ), axis ).toPointRange(); // << 2-cells + auto V12 = LatticeSet( C12.getLatticeSet( axis ), axis ).toPointRange(); // << 2-cells + auto P012 = P01 + StrictUnitSegment( 2 ); + Counter C012( P012 ); + auto V012 = LatticeSet( C012.getLatticeSet( axis ), axis ).toPointRange(); // << 3-cells + // Convert points to Khalimsky coordinates and insert them in LatticeSet. + LS = LatticeSet( axis ); + for ( const auto& p : V ) LS.insert( Point( 2*p[0] , 2*p[1] , 2*p[2] ) ); + for ( const auto& p : V0 ) LS.insert( Point( 2*p[0]-1, 2*p[1] , 2*p[2] ) ); + for ( const auto& p : V1 ) LS.insert( Point( 2*p[0] , 2*p[1]-1, 2*p[2] ) ); + for ( const auto& p : V2 ) LS.insert( Point( 2*p[0] , 2*p[1] , 2*p[2]-1 ) ); + for ( const auto& p : V01 ) LS.insert( Point( 2*p[0]-1, 2*p[1]-1, 2*p[2] ) ); + for ( const auto& p : V02 ) LS.insert( Point( 2*p[0]-1, 2*p[1] , 2*p[2]-1 ) ); + for ( const auto& p : V12 ) LS.insert( Point( 2*p[0] , 2*p[1]-1, 2*p[2]-1 ) ); + for ( const auto& p : V012 ) LS.insert( Point( 2*p[0]-1, 2*p[1]-1, 2*p[2]-1 ) ); + return LS; + } + else + { + using InternalInteger + = typename detail::ConvexityHelperInternalInteger< Integer, false >::Type; + using Helper = ConvexityHelper< dimension, Integer, InternalInteger >; + using StrictUnitSegment = typename Helper::LatticePolytope::StrictUnitSegment; + auto P = Helper::compute3DTriangle( a, b, c, true ); + if ( ! P.isValid() ) return LS; + Counter C( P ); + if ( axis >= dimension ) axis = C.longestAxis(); + auto V = LatticeSet( C.getLatticeSet( axis ), axis ).toPointRange(); // << 0-cells + auto P0 = P + StrictUnitSegment( 0 ); + auto P1 = P + StrictUnitSegment( 1 ); + auto P2 = P + StrictUnitSegment( 2 ); + Counter C0( P0 ); + Counter C1( P1 ); + Counter C2( P2 ); + auto V0 = LatticeSet( C0.getLatticeSet( axis ), axis ).toPointRange(); // << 1-cells + auto V1 = LatticeSet( C1.getLatticeSet( axis ), axis ).toPointRange(); // << 1-cells + auto V2 = LatticeSet( C2.getLatticeSet( axis ), axis ).toPointRange(); // << 1-cells + auto P01 = P0 + StrictUnitSegment( 1 ); + auto P02 = P0 + StrictUnitSegment( 2 ); + auto P12 = P1 + StrictUnitSegment( 2 ); + Counter C01( P01 ); + Counter C02( P02 ); + Counter C12( P12 ); + auto V01 = LatticeSet( C01.getLatticeSet( axis ), axis ).toPointRange(); // << 2-cells + auto V02 = LatticeSet( C02.getLatticeSet( axis ), axis ).toPointRange(); // << 2-cells + auto V12 = LatticeSet( C12.getLatticeSet( axis ), axis ).toPointRange(); // << 2-cells + auto P012 = P01 + StrictUnitSegment( 2 ); + Counter C012( P012 ); + auto V012 = LatticeSet( C012.getLatticeSet( axis ), axis ).toPointRange(); // << 3-cells + // Convert points to Khalimsky coordinates and insert them in LatticeSet. + LS = LatticeSet( axis ); + for ( const auto& p : V ) LS.insert( Point( 2*p[0] , 2*p[1] , 2*p[2] ) ); + for ( const auto& p : V0 ) LS.insert( Point( 2*p[0]-1, 2*p[1] , 2*p[2] ) ); + for ( const auto& p : V1 ) LS.insert( Point( 2*p[0] , 2*p[1]-1, 2*p[2] ) ); + for ( const auto& p : V2 ) LS.insert( Point( 2*p[0] , 2*p[1] , 2*p[2]-1 ) ); + for ( const auto& p : V01 ) LS.insert( Point( 2*p[0]-1, 2*p[1]-1, 2*p[2] ) ); + for ( const auto& p : V02 ) LS.insert( Point( 2*p[0]-1, 2*p[1] , 2*p[2]-1 ) ); + for ( const auto& p : V12 ) LS.insert( Point( 2*p[0] , 2*p[1]-1, 2*p[2]-1 ) ); + for ( const auto& p : V012 ) LS.insert( Point( 2*p[0]-1, 2*p[1]-1, 2*p[2]-1 ) ); + return LS; + } +} + + //----------------------------------------------------------------------------- template typename DGtal::DigitalConvexity::LatticeSet @@ -929,6 +1074,29 @@ isFullySubconvex( const Point& a, const Point& b, const Point& c, const auto SCabc = StarCvxH( a, b, c, StarX.axis() ); return StarX.includes( SCabc ); } +//----------------------------------------------------------------------------- +template +bool +DGtal::DigitalConvexity:: +isFullyCovered( const Point& a, const Point& b, const Point& c, + const LatticeSet& cells ) const +{ + ASSERT( dimension == 3 ); + const auto SCabc = CoverCvxH( a, b, c, cells.axis() ); + return cells.includes( SCabc ); +} + +//----------------------------------------------------------------------------- +template +bool +DGtal::DigitalConvexity:: +isOpenTriangleFullySubconvex( const Point& a, const Point& b, const Point& c, + const LatticeSet& StarX ) const +{ + ASSERT( dimension == 3 ); + const auto SCabc = StarOpenTriangle( a, b, c, StarX.axis() ); + return StarX.includes( SCabc ); +} //----------------------------------------------------------------------------- template diff --git a/tests/geometry/volumes/testConvexityHelper.cpp b/tests/geometry/volumes/testConvexityHelper.cpp index 94fe648a38..2b9af8c5c8 100644 --- a/tests/geometry/volumes/testConvexityHelper.cpp +++ b/tests/geometry/volumes/testConvexityHelper.cpp @@ -629,3 +629,158 @@ SCENARIO( "ConvexityHelper< 3 > degenerated triangle tests", } } } + + +SCENARIO( "ConvexityHelper< 3 > open segment tests", + "[convexity_helper][open segment][3d]" ) +{ + typedef ConvexityHelper< 3 > Helper; + typedef Helper::Point Point; + typedef Helper::Vector Vector; + typedef Helper::LatticePolytope::UnitSegment UnitSegment; + + auto nb_open_segment_smaller_than_segment = 0; + auto nb_dilated_open_segment_smaller_than_dilated_segment = 0; + auto nb_dilated_open_segment_greater_than_interior_dilated_segment = 0; + const int N = 100; + const int K = 10; + for ( auto n = 0; n < N; n++ ) + { + Point a( rand() % K, rand() % K, rand() % K ); + Point b( rand() % K, rand() % K, rand() % K ); + if ( a == b ) b[ rand() % 3 ] += 1; + Helper::LatticePolytope CS = Helper::computeSegment( a, b ); + Helper::LatticePolytope OS = Helper::computeOpenSegment( a, b ); + { + CAPTURE( CS ); + CAPTURE( OS ); + auto cs_in = CS.count(); + auto os_in = OS.count(); + REQUIRE( os_in < cs_in ); + nb_open_segment_smaller_than_segment += ( os_in < cs_in ) ? 1 : 0; + } + for ( Dimension k = 0; k < 3; k++ ) + { + CS += UnitSegment( k ); + OS += UnitSegment( k ); + } + { + CAPTURE( CS ); + CAPTURE( OS ); + auto cs_in = CS.count(); + auto cs_int = CS.countInterior(); + auto os_in = OS.count(); + REQUIRE( os_in < cs_in ); + REQUIRE( cs_int <= os_in ); + nb_dilated_open_segment_smaller_than_dilated_segment + += ( os_in < cs_in ) ? 1 : 0; + nb_dilated_open_segment_greater_than_interior_dilated_segment + += ( cs_int <= os_in ) ? 1 : 0; + } + } + WHEN( "Computing open segments" ) { + THEN( "They contain less lattice points than closed segments" ) { + REQUIRE( nb_open_segment_smaller_than_segment == N ); + } + THEN( "When dilated, they contain less lattice points than dilated closed segments" ) { + REQUIRE( nb_dilated_open_segment_smaller_than_dilated_segment == N ); + } + THEN( "When dilated, they contain more lattice points than the interior of dilated closed segments" ) { + REQUIRE( nb_dilated_open_segment_greater_than_interior_dilated_segment == N ); + } + } +} + + +SCENARIO( "ConvexityHelper< 3 > open triangle tests", + "[convexity_helper][open triangle][3d]" ) +{ + typedef ConvexityHelper< 3 > Helper; + typedef Helper::Point Point; + typedef Helper::Vector Vector; + typedef Helper::LatticePolytope::UnitSegment UnitSegment; + + auto nb_open_triangle_smaller_than_triangle = 0; + auto nb_dilated_open_triangle_smaller_than_dilated_triangle = 0; + auto nb_dilated_open_triangle_greater_than_interior_dilated_triangle = 0; + const int N = 100; + const int K = 10; + for ( auto n = 0; n < N; n++ ) + { + Point a( rand() % K, rand() % K, rand() % K ); + Point b( rand() % K, rand() % K, rand() % K ); + Point c( rand() % K, rand() % K, rand() % K ); + if ( a == b ) b[ rand() % 3 ] += 1; + if ( b == c ) c[ rand() % 3 ] += 1; + if ( c == a ) c[ rand() % 3 ] += 1; + if ( b == c ) c[ rand() % 3 ] += 1; + CAPTURE( a, b, c ); + Helper::LatticePolytope CS = Helper::compute3DTriangle( a, b, c, true ); + Helper::LatticePolytope OS = Helper::compute3DOpenTriangle( a, b, c, true ); + { + CAPTURE( CS ); + CAPTURE( OS ); + auto cs_in = CS.count(); + auto os_in = OS.count(); + REQUIRE( os_in < cs_in ); + nb_open_triangle_smaller_than_triangle += ( os_in < cs_in ) ? 1 : 0; + } + for ( Dimension k = 0; k < 3; k++ ) + { + CS += UnitSegment( k ); + OS += UnitSegment( k ); + } + { + CAPTURE( CS ); + CAPTURE( OS ); + auto cs_in = CS.count(); + auto cs_int = CS.countInterior(); + auto os_in = OS.count(); + REQUIRE( os_in < cs_in ); + REQUIRE( cs_int <= os_in ); + nb_dilated_open_triangle_smaller_than_dilated_triangle + += ( os_in < cs_in ) ? 1 : 0; + nb_dilated_open_triangle_greater_than_interior_dilated_triangle + += ( cs_int <= os_in ) ? 1 : 0; + } + } + WHEN( "Computing open triangles" ) { + THEN( "They contain less lattice points than closed triangles" ) { + REQUIRE( nb_open_triangle_smaller_than_triangle == N ); + } + THEN( "When dilated, they contain less lattice points than dilated closed triangles" ) { + REQUIRE( nb_dilated_open_triangle_smaller_than_dilated_triangle == N ); + } + THEN( "When dilated, they contain more lattice points than the interior of dilated closed triangles" ) { + REQUIRE( nb_dilated_open_triangle_greater_than_interior_dilated_triangle == N ); + } + } +} + +SCENARIO( "ConvexityHelper< 3 > open triangle unit tests", + "[convexity_helper][open triangle][3d]" ) +{ + typedef ConvexityHelper< 3 > Helper; + typedef Helper::Point Point; + typedef Helper::Vector Vector; + typedef Helper::LatticePolytope::UnitSegment UnitSegment; + + Point a( 0, 0, 0 ); + Point b( 2, 1, 1 ); + Point c( 2, 2, 1 ); + CAPTURE( a, b, c ); + Helper::LatticePolytope CS = Helper::compute3DTriangle( a, b, c, true ); + Helper::LatticePolytope OS = Helper::compute3DOpenTriangle( a, b, c, true ); + for ( Dimension k = 0; k < 3; k++ ) + { + CS += UnitSegment( k ); + OS += UnitSegment( k ); + } + std::vector< Point > PCS, POS; + CS.getPoints( PCS ); + OS.getPoints( POS ); + CAPTURE( PCS, POS ); + REQUIRE( PCS.size() == 21 ); + REQUIRE( POS.size() == 3 ); +} + diff --git a/tests/geometry/volumes/testDigitalConvexity.cpp b/tests/geometry/volumes/testDigitalConvexity.cpp index eb15ce1385..43297325fd 100644 --- a/tests/geometry/volumes/testDigitalConvexity.cpp +++ b/tests/geometry/volumes/testDigitalConvexity.cpp @@ -1024,3 +1024,113 @@ SCENARIO( "DigitalConvexity< Z3 > full subconvexity of points and triangles", "[ } } + +SCENARIO( "DigitalConvexity< Z3 > full covering of triangles", "[full_cover][3d]" ) +{ + typedef KhalimskySpaceND<3,int> KSpace; + typedef KSpace::Point Point; + typedef KSpace::Vector Vector; + typedef KSpace::Space Space; + typedef HyperRectDomain< Space > Domain; + typedef DigitalConvexity< KSpace > DConvexity; + + Domain domain( Point( -20, -20, -20 ), Point( 20, 20, 20 ) ); + DConvexity dconv ( Point( -21, -21, -21 ), Point( 21, 21, 21 ) ); + { + Point a( 0, 0, 0 ); + Point b( 2, 1, 0 ); + Point c( 2, 2, 0 ); + auto LS = dconv.CoverCvxH( a, b, c ); + auto P = LS.toPointRange(); + CAPTURE( P ); + REQUIRE( P.size() == 10 ); + } + { + Point a( 0, 0, 1 ); + Point b( 2, 1, 0 ); + Point c( 2, 2, 0 ); + auto LS = dconv.CoverCvxH( a, b, c ); + auto P = LS.toPointRange(); + CAPTURE( P ); + REQUIRE( P.size() == 10 ); + } + { + Point a( 1, 0, 0 ); + Point b( 0, 1, 0 ); + Point c( 0, 0, 1 ); + auto LS = dconv.CoverCvxH( a, b, c ); + auto P = LS.toPointRange(); + CAPTURE( P ); + REQUIRE( P.size() == 7 ); + } +} + + +SCENARIO( "DigitalConvexity< Z3 > full subconvexity and full covering of triangles", "[subconvexity][3d][full_cover]" ) +{ + typedef KhalimskySpaceND<3,int> KSpace; + typedef KSpace::Point Point; + typedef KSpace::Vector Vector; + typedef KSpace::Space Space; + typedef HyperRectDomain< Space > Domain; + typedef DigitalConvexity< KSpace > DConvexity; + + Domain domain( Point( -20, -20, -20 ), Point( 20, 20, 20 ) ); + DConvexity dconv ( Point( -21, -21, -21 ), Point( 21, 21, 21 ) ); + + WHEN( "Computing many tetrahedra" ) { + const unsigned int nb = 50; + unsigned int nb_total = 0; + unsigned int nb_ok_tri = 0; + unsigned int nb_subconvex = 0; + unsigned int nb_covered = 0; + for ( unsigned int l = 0; l < nb; ++l ) + { + const Point a { (rand() % 10 - 10), (rand() % 10 - 10), (rand() % 10 - 10) }; + const Point b { (rand() % 20 ), (rand() % 20 - 10), (rand() % 20 - 10) }; + const Point c { (rand() % 20 - 10), (rand() % 20 ), (rand() % 20 - 10) }; + const Point d { (rand() % 20 - 10), (rand() % 20 - 10), (rand() % 20 ) }; + const std::vector pts = { a, b, c, d }; + const bool fulldim = dconv.isSimplexFullDimensional(pts.cbegin(), pts.cend()); + if ( ! fulldim ) continue; + auto simplex = dconv.makeSimplex( pts.cbegin(), pts.cend() ); + auto cover = dconv.makeCellCover( simplex, 0, 3 ); + auto ls = dconv.StarCvxH( pts ); + { + for ( unsigned int i = 0; i < 100; i++ ) + { + const Point p { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) }; + const Point q { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) }; + const Point r { (rand() % 10 - 5), (rand() % 10 - 5), (rand() % 10 - 5) }; + const Vector n = ( q - p ).crossProduct( r - p ); + if ( n == Vector::zero ) continue; + auto tri1 = dconv.makeSimplex( { p, q, r } ); + bool ok1 = dconv.isFullySubconvex( p, q, r, ls ); + bool ok2 = dconv.isFullyCovered( p, q, r, ls ); + nb_subconvex += ok1 ? 1 : 0; + nb_covered += ok2 ? 1 : 0; + bool ok = ok1 == ok2; + if ( ! ok ) + { + std::cout << "***** FULL SUBCONVEXITY ERROR ON TRIANGLE ****" << std::endl; + std::cout << "splx v =" << a << b << c << d << std::endl; + std::cout << "simplex=" << simplex << std::endl; + std::cout << "tri v =" << p << q << r << std::endl; + std::cout << "tri1=" << tri1 << std::endl; + std::cout << "3 points are fully subconvex: " << ( ok1 ? "YES" : "NO" ) << std::endl; + std::cout << "3 points are fully covered by star: " << ( ok2 ? "YES" : "NO" ) << std::endl; + } + nb_ok_tri += ( ok ) ? 1 : 0; + nb_total += 1; + } + } + } + THEN( "The number of subconvex and covered triangles should be equal." ) { + REQUIRE( nb_subconvex == nb_covered ); + } + THEN( "Full subconvexity and covering should agree on every subset." ) { + REQUIRE( nb_ok_tri == nb_total ); + } + } + +}