Add full covering of triangle and full subconvexity of open triangles in 3D

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
 

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
 

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<Dimension> dims;
+      StrictUnitCell( std::initializer_list<Dimension> 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,26 +552,42 @@ 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'.
      */
     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'.
      */
     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<TSpace> & P,
               typename BoundedLatticePolytope<TSpace>::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 <typename TSpace>
+  BoundedLatticePolytope<TSpace>
+  operator+ ( const BoundedLatticePolytope<TSpace> & P,
+              typename BoundedLatticePolytope<TSpace>::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 <typename TSpace>
+  BoundedLatticePolytope<TSpace>
+  operator+ ( const BoundedLatticePolytope<TSpace> & P,
+              typename BoundedLatticePolytope<TSpace>::StrictUnitCell c );
+
   /**
    * Minkowski sum of polytope \a P with strict unit segment \a s aligned with some axis.
    *

src/DGtal/geometry/volumes/BoundedLatticePolytope.ih

@@ -503,6 +503,28 @@ operator+=( UnitSegment s )
   return *this;
 }
 
+//-----------------------------------------------------------------------------
+template <typename TSpace>
+typename DGtal::BoundedLatticePolytope<TSpace>::Self&
+DGtal::BoundedLatticePolytope<TSpace>::
+operator+=( StrictUnitSegment s )
+{
+  for ( Dimension i = 0; i < A.size(); ++i )
+    {
+      if ( A[ i ][ s.k ] > NumberTraits<Integer>::ZERO )
+	{
+	  B[ i ] += A[ i ][ s.k ];
+	  I[ i ]  = false;
+	}
+      else if ( A[ i ][ s.k ] < NumberTraits<Integer>::ZERO )
+        I[ i ] = false;
+    }
+  Vector z = Vector::zero;
+  z[ s.k ] = NumberTraits<Integer>::ONE;
+  D = Domain( D.lowerBound(), D.upperBound() + z );
+  return *this;
+}
+
 //-----------------------------------------------------------------------------
 template <typename TSpace>
 typename DGtal::BoundedLatticePolytope<TSpace>::Self&
@@ -553,6 +575,17 @@ operator+=( UnitCell c )
   return *this;
 }
 
+//-----------------------------------------------------------------------------
+template <typename TSpace>
+typename DGtal::BoundedLatticePolytope<TSpace>::Self&
+DGtal::BoundedLatticePolytope<TSpace>::
+operator+=( StrictUnitCell c )
+{
+  for ( Dimension i = 0; i < c.dims.size(); ++i )
+    *this += StrictUnitSegment( c.dims[ i ] );
+  return *this;
+}
+
 //-----------------------------------------------------------------------------
 template <typename TSpace>
 typename DGtal::BoundedLatticePolytope<TSpace>::Self&
@@ -927,6 +960,25 @@ DGtal::BoundedLatticePolytope<TSpace>::isLarge( unsigned int i ) const
   return I[ i ];
 }
 
+//-----------------------------------------------------------------------------
+template <typename TSpace>
+void
+DGtal::BoundedLatticePolytope<TSpace>::setLarge( unsigned int i )
+{
+  ASSERT( i < nbHalfSpaces() );
+  I[ i ] = true;
+}
+
+//-----------------------------------------------------------------------------
+template <typename TSpace>
+void
+DGtal::BoundedLatticePolytope<TSpace>::setStrict( unsigned int i )
+{
+  ASSERT( i < nbHalfSpaces() );
+  I[ i ] = false;
+}
+
+
 //-----------------------------------------------------------------------------
 template <typename TSpace>
 const typename DGtal::BoundedLatticePolytope<TSpace>::InequalityMatrix&
@@ -979,7 +1031,7 @@ DGtal::BoundedLatticePolytope<TSpace>::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<TSpace> & P,
 //-----------------------------------------------------------------------------
 template <typename TSpace>
 DGtal::BoundedLatticePolytope<TSpace>
+DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
+                   typename BoundedLatticePolytope<TSpace>::StrictUnitSegment s )
+{
+  BoundedLatticePolytope<TSpace> Q = P;
+  Q += s;
+  return Q;
+}
+//-----------------------------------------------------------------------------
+template <typename TSpace>
+DGtal::BoundedLatticePolytope<TSpace>
+DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
+                   typename BoundedLatticePolytope<TSpace>::StrictUnitCell c )
+{
+  BoundedLatticePolytope<TSpace> Q = P;
+  Q += c;
+  return Q;
+}
+//-----------------------------------------------------------------------------
+template <typename TSpace>
+DGtal::BoundedLatticePolytope<TSpace>
 DGtal::operator+ ( const BoundedLatticePolytope<TSpace> & P,
                    typename BoundedLatticePolytope<TSpace>::RightStrictUnitSegment s )
 {

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;
 }
 
 //-----------------------------------------------------------------------------

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 );
 
     /// @}