Merge pull request #1717 from JacquesOlivierLachaud/AddTriangleToConv…

…exityHelper

Add triangle to convexity helper

ChangeLog.md

@@ -43,6 +43,11 @@
   - The WindingNumberShape class can output the raw winding number values
     (David Coeurjolly,[#1719](https://github.com/DGtal-team/DGtal/issues/1719))
 
+- *Geometry package*
+  - Add creation of polytopes from segments and triangles in
+    ConvexityHelper and 3-5xfaster full subconvexity tests for triangles
+    in DigitalConvexity (Jacques-Olivier Lachaud,
+    [#1717](https://github.com/DGtal-team/DGtal/pull/1717))
 
 - *Project*
   - Add CMake option DGTAL_WRAP_PYTHON (Pablo Hernandez-Cerdan,
@@ -96,6 +101,8 @@
 - *Geometry package*
   - Fix Issue #1676 in testStabbingCircleComputer (Tristan Roussillon,
     [#1688](https://github.com/DGtal-team/DGtal/pull/1688)
+  - Fix BoundedLatticePolytopeCounter::countInterior method (Jacques-Olivier Lachaud,
+    [#1717](https://github.com/DGtal-team/DGtal/pull/1717))
   - Fix const attribute that shouldn't be in FreemanChain (Colin Weill--Duflos,
     [#1723](https://github.com/DGtal-team/DGtal/pull/1723))
 

src/DGtal/geometry/doc/moduleDigitalConvexity.dox

@@ -199,6 +199,12 @@ Lattice point retrieval services:
 - BoundedLatticePolytope::getBoundaryPoints outputs the lattice points on the boundary of the polytope
 - BoundedLatticePolytope::insertPoints inserts the lattice points in the polytope into some point set
 
+The class ConvexityHelper also provides several static methods related to BoundedLatticePolytope:
+- ConvexityHelper::computeLatticePolytope computes a BoundedLatticePolytope from an arbitrary range of points (use QuickHull algorithm, \see moduleQuickHull).
+- ConvexityHelper::compute3DTriangle computes the BoundedLatticePolytope enclosing a 3D triangle (faster way than above, but limited to 3D triangles).
+- ConvexityHelper::computeDegeneratedTriangle computes the BoundedLatticePolytope enclosing a degenerated triangle (the three points are aligned or non distinct, but may lie in arbitrary dimension).
+- ConvexityHelper::computeSegment computes the BoundedLatticePolytope enclosing a segment in arbitrary dimension.
+
 
 @subsection dgtal_dconvexity_sec22 Building a set of lattice cells from digital points
 
@@ -304,8 +310,8 @@ bool ok  = simplex_cover.subset( point_cover ); // true
 
 Morphological services (since 1.3):
 - DigitalConvexity::makePolytope( const PointRange& X, bool safe ) const builds the tightest polytope enclosing the digital set \a X
-- DigitalConvexity::isFullyConvex( const PointRange& X, bool convex0, bool safe ) const checks the full convexity of the digital set \a X
-- DigitalConvexity::is0Convex( const PointRange& X, bool safe ) const checks the usual digital convexity (H-convexity or 0-convexity) of the digital set \a X
+- DigitalConvexity::isFullyConvex( const PointRange& X, bool convex0 ) const checks the full convexity of the digital set \a X
+- DigitalConvexity::is0Convex( const PointRange& X ) const checks the usual digital convexity (H-convexity or 0-convexity) of the digital set \a X
 - DigitalConvexity::U( Dimension i, const PointRange& X ) const performs the digital Minkowski sum of \a X along direction \a i
 
 The following snippet shows that a 4D ball is indeed fully convex.
@@ -334,9 +340,13 @@ Convexity services:
 - DigitalConvexity::isKConvex tells if a given polytope is k-convex
 - DigitalConvexity::isFullyConvex tells if a given polytope is fully convex
 - DigitalConvexity::isKSubconvex tells if a given polytope is k-subconvex to some cell cover (i.e. k-tangent)
-- DigitalConvexity::isKSubconvex( const Point&, const Point&, const CellGeometry&, const Dimension ) const tells if the Euclidean straight segment \f$ \lbrack a,b \rbrack \f$ is k-subconvex (i.e. k-tangent) to the given cell cover.
+- DigitalConvexity::isKSubconvex( const Point& a, const Point& b, const CellGeometry& C, const Dimension k ) const tells if the Euclidean straight segment \f$ \lbrack a,b \rbrack \f$ is k-subconvex (i.e. k-tangent) to the given cell cover.
 - DigitalConvexity::isFullySubconvex tells if a given polytope is fully subconvex to some cell cover
-- DigitalConvexity::isFullySubconvex( const Point&, const Point&, const CellGeometry& ) const tells if the Euclidean straight segment \f$ \lbrack a,b \rbrack \f$ is fully k-subconvex (i.e. tangent) to the given cell cover.
+- DigitalConvexity::isFullySubconvex( const Point& a, const Point& b, const CellGeometry& C ) const tells if the Euclidean straight segment \f$ \lbrack a,b \rbrack \f$ is fully k-subconvex (i.e. tangent) to the given cell cover.
+- 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).
 
 @subsection dgtal_dconvexity_sec25 Ehrhart polynomial of a lattice polytope
 

src/DGtal/geometry/volumes/BoundedLatticePolytopeCounter.ih

@@ -128,7 +128,8 @@ interiorIntersectionIntervalAlongAxis( 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++ )
+  // We must take into account also bounding box constraints for interior points.
+  for ( Dimension k = 0; k < A.size(); k++ )
     {
       const Integer c = A[ k ].dot( p );
       const Integer n = A[ k ][ a ];

src/DGtal/geometry/volumes/ConvexityHelper.h

@@ -379,7 +379,51 @@ namespace DGtal
     static
     PointRange
     computeDegeneratedConvexHullVertices( PointRange& input_points );
-
+
+    /// Computes the lattice polytope enclosing a 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 range of points.
+    static
+    LatticePolytope
+    compute3DTriangle( 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).
+    ///
+    /// @param a any point
+    /// @param b any point
+    /// @param c any point
+    ///
+    /// @return the tightiest bounded lattice polytope
+    /// (i.e. H-representation) including the given range of points.
+    static
+    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.
+    static
+    LatticePolytope
+    computeSegment( const Point& a, const Point& b );
 
     /// @}
 

src/DGtal/geometry/volumes/ConvexityHelper.ih

@@ -183,6 +183,129 @@ computeDegeneratedLatticePolytope
                           false, false );
 }
 
+//-----------------------------------------------------------------------------
+template < int dim, typename TInteger, typename TInternalInteger >
+typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
+DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::compute3DTriangle
+( const Point& a, const Point& b, const Point& c,
+  bool make_minkowski_summable )
+{
+  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 Vector ab   = b - a;
+  const Vector bc   = c - b;
+  const Vector ca   = a - c;
+  const Vector n    = Op::crossProduct( ab, bc );
+  if ( n == Vector::zero )
+    return computeDegeneratedTriangle( a, b, c );
+  const Point  low  = a.inf( b ).inf( c );
+  const Point  high = a.sup( b ).sup( c );
+  // Initialize polytope
+  std::vector< PolytopeHalfSpace >   PHS;
+  PHS.reserve( make_minkowski_summable ? 11 : 5 );
+  const Integer n_a = n.dot( a );
+  const Vector  u   = Op::crossProduct( ab, n );
+  const Vector  v   = Op::crossProduct( bc, n );
+  const Vector  w   = Op::crossProduct( ca, n );
+  PHS.emplace_back(  n,  n_a  );
+  PHS.emplace_back( -n, -n_a );
+  if ( ! make_minkowski_summable )
+    { // It is enough to have one constraint per edge.
+      PHS.emplace_back(  u, u.dot( a ) );
+      PHS.emplace_back(  v, v.dot( b ) );
+      PHS.emplace_back(  w, w.dot( c ) );
+    }
+  else
+    { // 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  eab   = Op::crossProduct( ab, i );
+	    const Integer eab_a = eab.dot( a );
+	    if ( eab.dot( c ) < eab_a ) // c must be below plane (a,eab)
+	      PHS.emplace_back( eab, eab_a );
+	    const Vector  ebc   = Op::crossProduct( bc, i );
+	    const Integer ebc_b = ebc.dot( b );
+	    if ( ebc.dot( a ) < ebc_b ) // a must be below plane (b,ebc)
+	      PHS.emplace_back( ebc, ebc_b );
+	    const Vector  eca   = Op::crossProduct( ca, i );
+	    const Integer eca_c = eca.dot( c );
+	    if ( eca.dot( b ) < eca_c ) // b must be below plane (c,eca)
+	      PHS.emplace_back( eca, eca_c );
+	  }
+    }
+  return LatticePolytope( Domain( low, high ),
+			  PHS.cbegin(), PHS.cend(),
+                          make_minkowski_summable, false );
+}
+
+//-----------------------------------------------------------------------------
+template < int dim, typename TInteger, typename TInternalInteger >
+typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
+DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
+computeDegeneratedTriangle
+( const Point& a, const Point& b, const Point& c )
+{
+  if ( a == b ) return computeSegment( a, c );
+  if ( ( a == c ) || ( b == c ) ) return computeSegment( a, b );
+  // The three points are distinct, hence aligned. One is in-between the two others.
+  const Point  low  = a.inf( b ).inf( c );
+  const Point  high = a.sup( b ).sup( c );
+  for ( Dimension k = 0; k < dim; k++ )
+    {
+      const auto lk = low [ k ];
+      const auto hk = high[ k ];      
+      if ( ( a[ k ] != lk ) && ( a[ k ] != hk ) ) return computeSegment( b, c );
+      if ( ( b[ k ] != lk ) && ( b[ k ] != hk ) ) return computeSegment( a, c );
+      if ( ( c[ k ] != lk ) && ( c[ k ] != hk ) ) return computeSegment( a, b );
+    }
+  trace.error() << "[ConvexityHelper::computeSegmentFromDegeneratedTriangle] "
+		<< "Should never arrive here." << std::endl;
+  return computeSegment( a, a );
+}
+
+//-----------------------------------------------------------------------------
+template < int dim, typename TInteger, typename TInternalInteger >
+typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
+DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeSegment
+( 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 );
+  // Initialize polytope
+  std::vector< PolytopeHalfSpace >   PHS;
+  if ( degenerate ) 
+    return LatticePolytope( Domain( low, high ), PHS.cbegin(), PHS.cend(), true, false );
+  PHS.reserve( 2*dim );
+  // 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 );
+	  }
+      }
+  return LatticePolytope( Domain( low, high ), PHS.cbegin(), PHS.cend(), true, false );
+}
+
+
 //-----------------------------------------------------------------------------
 template < int dim, typename TInteger, typename TInternalInteger >
 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::PointRange
@@ -306,10 +429,12 @@ computeLatticePolytope
   // Compute domain
   Point l = input_points[ 0 ];
   Point u = input_points[ 0 ];
-  for ( const auto& p : input_points ) {
-    l = l.inf( p );
-    u = u.sup( p );
-  }
+  for ( std::size_t i = 1; i < input_points.size(); i++ )
+    {
+      const auto& p = input_points[ i ];
+      l = l.inf( p );
+      u = u.sup( p );
+    }
   Domain domain( l, u );
   // Compute convex hull
   ConvexHull hull;

src/DGtal/geometry/volumes/DigitalConvexity.h

@@ -437,6 +437,24 @@ namespace DGtal
     /// generic since the two other methods require a Minkowski
     /// summable polytope, i.e. `P.canBeSummed() == true`.
     bool isFullySubconvex( const PointRange& Y, const LatticeSet& StarX ) const;
+
+
+    /// Tells if the non-degenerated 3D triangle a,b,c 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 triangles.
+    bool isFullySubconvex( 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
@@ -491,7 +509,7 @@ namespace DGtal
     /// polytope and then checking if it subconvex.
     bool isFullySubconvex( const Point& a, const Point& b,
                            const LatticeSet& StarX ) const;
-
+    
     /// Given a range of distinct points \a X, computes the tightiest
     /// polytope that enclosed it. Note that this polytope may contain
     /// more lattice points than the given input points.
@@ -537,6 +555,30 @@ namespace DGtal
     LatticeSet StarCvxH( const PointRange& X,
                          Dimension axis = dimension ) const;
 
+    /// Builds the cell complex `Star(CvxH({a,b,c}))` for `a,b,c` a
+    /// non-degenerate 3D triangle, 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 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 StarCvxH( 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.
     ///
     /// @param X any range of lattice points