Specializes DigitalSurface::facesAroundVertex in the 3D case

ChangeLog.md

@@ -4,7 +4,11 @@
 ## New Features / Critical Changes
 
 - *Helpers*
-  - Classes Shortcuts and ShortcutsGeometry to simplify coding with DGtal. Integrate a lot of volume, digital surfaces, mesh, surface, geometry, estimators functions, with many conversion and input/output tools. (Jacques-Olivier Lachaud, [#1357](https://github.com/DGtal-team/DGtal/pull/1357))
+  - Classes Shortcuts and ShortcutsGeometry to simplify coding with
+    DGtal. Integrate a lot of volume, digital surfaces, mesh,
+    surface, geometry, estimators functions, with many conversion
+    and input/output tools. (Jacques-Olivier Lachaud,
+    [#1357](https://github.com/DGtal-team/DGtal/pull/1357))
 
 ## Changes
 
@@ -90,6 +94,13 @@
   - Small fixes in Shortcuts and ShortcutsGeometry, doc, and colormaps.
     (Jacques-Olivier Lachaud, [#1364](https://github.com/DGtal-team/DGtal/pull/1364))
 
+- *Topology*
+  - Specializes the method DigitalSurface::facesAroundVertex in the
+    3D case, such that faces (ie pointels) are ordered
+    counterclockwise with respect of the vertex (ie surfel) seen from
+    the exterior. (Jacques-Olivier Lachaud,
+    [#1377](https://github.com/DGtal-team/DGtal/pull/1377))
+
 - *Miscellaneous*
   - Fix Small bug in Integral Invariant Volume Estimator in 2D
     (Thomas Caissard, [#1316](https://github.com/DGtal-team/DGtal/pull/1316))

src/DGtal/topology/DigitalSurface.h

@@ -82,6 +82,24 @@ namespace DGtal
   manifold (i.e. a kind of surface), except when it has the
   property of being well-composed.
 
+  @note On \b interior / \b exterior side of a surface. A digital
+  surface is composed of oriented n-1-cells, and it is thus
+  oriented. By convention, the interior side of a digital surface
+  (globally) or of a surfel (locally) is in its direct incidence,
+  while the exterior side is in its indirect incidence. You can obtain
+  the normal vector to a surfel pointing outside as follows:
+
+  \code
+  // S is the current digital surface
+  // v is a vertex / surfel
+  auto& K = S.container().space(); 
+  auto  j = K.sOrthDir( v );
+  bool  d = K.sDirect( v, j );
+  // N is the normal vector pointing outside.
+  auto  N = K.sCoords( K.sIncident( v, j, ! d ) )   // outside voxel
+            - K.sCoords( K.sIncident( v, j, d ) );  // inside voxel
+  \endcode
+
   For geometric analysis or visualization, it is often interesting
   to look at the "dual" of the digital surface. n-1-cells form now
   vertices, n-2-cells are edges, n-3-cells are faces, and so on.  A
@@ -115,7 +133,7 @@ namespace DGtal
   CDigitalSurfaceContainer: the concrete representation chosen for
   the digital surface.
 
-  ee \ref moduleDigitalSurfaces
+  \ref moduleDigitalSurfaces
    */
   template <typename TDigitalSurfaceContainer>
   class DigitalSurface
@@ -144,7 +162,8 @@ namespace DGtal
     typedef typename DigitalSurfaceContainer::Surfel Surfel;
     typedef typename DigitalSurfaceContainer::SurfelConstIterator ConstIterator;
     typedef typename DigitalSurfaceContainer::DigitalSurfaceTracker DigitalSurfaceTracker; 
-    typedef typename KSpace::Point Point;
+    typedef typename KSpace::Point     Point;
+    typedef typename KSpace::Vector    Vector;
     typedef typename KSpace::SurfelSet SurfelSet;
     /// Template rebinding for defining the type that is a mapping
     /// SCell -> Value.
@@ -467,11 +486,23 @@ namespace DGtal
 
     /**
        @param v any vertex (surfel) of the surface.
+
+       @param order_ccw_in_3d when 'true', orders faces
+       counterclockwise around vertex (solely in 3d). It corresponds
+       to ordering the pointels counterclockwise around the given
+       surfel \a v, seen from the exterior side of the surfel.
+
        @return the faces containing this vertex [v]: 0 in 2D, 4 in 3D,
-       12 in 4D, 2(n-1)(n-2) in nD.
-    */
-    FaceRange facesAroundVertex( const Vertex & v ) const;
+       12 in 4D, 2(n-1)(n-2) in nD, in no specific order.
 
+       @note By convention, the interior side of a digital surface
+       (globally) or of a surfel (locally) is in its direct incidence,
+       while the exterior side is in its indirect incidence (see text
+       in \ref DigitalSurface description).
+    */
+    FaceRange facesAroundVertex( const Vertex & v,
+				 bool order_ccw_in_3d = false ) const;
+
     /**
       @param a any arc (s,t)
       @return the vertex t

src/DGtal/topology/DigitalSurface.ih

@@ -268,23 +268,45 @@ template <typename TDigitalSurfaceContainer>
 inline
 typename DGtal::DigitalSurface<TDigitalSurfaceContainer>::FaceRange
 DGtal::DigitalSurface<TDigitalSurfaceContainer>::
-facesAroundVertex( const Vertex & v ) const
+facesAroundVertex( const Vertex & v, bool order_ccw_in_3d ) const
 {
   typedef typename ArcRange::const_iterator ArcRangeConstIterator;
   // std::cerr << "  - facesAroundVertex(" << v << ")" << std::endl;
   ArcRange arcs = outArcs( v );
-  FaceRange faces;
-  std::back_insert_iterator<FaceRange> output_it = std::back_inserter( faces );
-  for ( ArcRangeConstIterator it = arcs.begin(), it_end = arcs.end();
-	it != it_end; ++it )
-    {
-      // std::cerr << "     + arc " << tail( *it )
-      //           << " -> " << head( *it ) << std::endl;
-      FaceRange faces_of_arc = facesAroundArc( *it );
-      output_it = 
-	std::copy( faces_of_arc.begin(), faces_of_arc.end(), output_it );
+  if ( order_ccw_in_3d && ( arcs.size() == 4 ) )
+    { // 3D method to order faces/pointels CCW around vertices/surfels
+      FaceRange faces;
+      faces.push_back( facesAroundArc( arcs[ 0 ] )[ 0 ] );
+      faces.push_back( facesAroundArc( arcs[ 2 ] )[ 0 ] );
+      faces.push_back( facesAroundArc( arcs[ 1 ] )[ 0 ] ); //< to change order.
+      faces.push_back( facesAroundArc( arcs[ 3 ] )[ 0 ] ); //< to change order.
+      const KSpace& K = container().space();
+      auto  orth_dir = K.sOrthDir( v );
+      auto    direct = K.sDirect( v, orth_dir ); // true: ccw, false: ! ccw
+      Point   vtcs[] = { K.sCoords( pivot( faces[ 0 ] ) ),
+			 K.sCoords( pivot( faces[ 1 ] ) ),
+			 K.sCoords( pivot( faces[ 2 ] ) ) };
+      Vector    s0s1 = vtcs[ 1 ] - vtcs[ 0 ];
+      Vector    s0s2 = vtcs[ 2 ] - vtcs[ 0 ];
+      Vector       t = s0s1.crossProduct( s0s2 );
+      if ( ( ( t[ orth_dir ] > 0 ) && direct )
+	   || ( ( t[ orth_dir ] < 0 ) && ! direct ) )
+	std::reverse( faces.begin(), faces.end() );
+      return faces;
+    }
+  else
+    { // nD method
+      FaceRange faces;
+      std::back_insert_iterator<FaceRange> output_it = std::back_inserter( faces );
+      for ( ArcRangeConstIterator it = arcs.begin(), it_end = arcs.end();
+	    it != it_end; ++it )
+	{
+	  FaceRange faces_of_arc = facesAroundArc( *it );
+	  output_it = 
+	    std::copy( faces_of_arc.begin(), faces_of_arc.end(), output_it );
+	}
+      return faces;
     }
-  return faces;
 }
 //-----------------------------------------------------------------------------
 template <typename TDigitalSurfaceContainer>

src/DGtal/topology/doc/moduleDigitalSurfaces.dox

@@ -563,8 +563,10 @@ edges or vertices and reciprocally.
 - DigitalSurface::inArcs (Vertex): the range of ingoing arcs from
   the given vertex.
 - DigitalSurface::facesAroundVertex (Vertex): the range of faces
-  containing this vertex [v]: 0 in 2D, 4 in 3D, 12 in 4D,
-  2(n-1)(n-2) in nD.
+  containing this vertex [v]: 0 in 2D, 4 in 3D, 12 in 4D, 2(n-1)(n-2)
+  in nD. Since DGtal 1.0, in 3D, the 4 faces may be ordered
+  counterclockwise around the vertex (seen from the exterior), if the
+  user specifies it. 
 - DigitalSurface::head (Arc): the head vertex of the given arc
 - DigitalSurface::tail (Arc): the tail vertex of the given arc
 - DigitalSurface::opposite (Arc): the opposite arc of the given arc

tests/topology/testDigitalSurface.cpp

@@ -585,6 +585,83 @@ bool testDigitalSurface()
   return nbok == nb;
 }
 
+bool testOrderingDigitalSurfaceFacesAroundVertex()
+{
+  typedef KhalimskySpaceND<3>     KSpace;
+  typedef typename KSpace::Space  Space;
+  typedef typename KSpace::Size   Size;
+  typedef typename Space::Point   Point;
+  typedef HyperRectDomain<Space>  Domain;
+  typedef typename DigitalSetSelector < Domain, BIG_DS + HIGH_ITER_DS + HIGH_BEL_DS >::Type
+                                  DigitalSet;
+  typedef DigitalSetBoundary<KSpace,DigitalSet> DSContainer;
+  typedef DigitalSurface<DSContainer>           MyDS;
+  typedef typename MyDS::Vertex                 Vertex;
+  typedef typename MyDS::SCell                  SCell;
+
+  unsigned int nbok = 0;
+  unsigned int nb = 0;
+  trace.beginBlock ( "Creating surface around one voxel" );
+  Point p0;
+  Point p1 = Point::diagonal( -1 );
+  Point p2 = Point::diagonal( 1 );
+  Domain domain( p1, p2 );
+  DigitalSet dig_set( domain );
+  dig_set.insert( p0 );
+  KSpace K;
+  nbok += K.init( domain.lowerBound(), domain.upperBound(), true ) ? 1 : 0; 
+  nb++;
+  trace.info() << "(" << nbok << "/" << nb << ") "
+	       << "K.init() is ok" << std::endl;
+  DSContainer* ptrBdry = new DSContainer( K, dig_set );
+  MyDS digsurf( ptrBdry ); // acquired
+  ++nb; nbok += digsurf.size() == 6 ? 1 : 0;
+  trace.info() << "(" << nbok << "/" << nb << ") "
+	       << "digsurf.size() = " << digsurf.size() 
+               << " == " << 6 << std::endl;
+  trace.endBlock();
+  trace.beginBlock ( "Check ordering" );
+  std::map< std::pair< SCell, SCell >, unsigned int > nb_arcs;
+  for ( Vertex v : digsurf )
+    {
+      auto faces = digsurf.facesAroundVertex( v, true );
+      for ( unsigned int i = 0; i < faces.size(); ++i )
+	{
+	  auto p = std::make_pair( digsurf.pivot( faces[ i ] ),
+				   digsurf.pivot( faces[ (i+1)%4 ] ) );
+	  trace.info() << "Arc " << p.first << " --- " << p.second << std::endl;
+	  if ( nb_arcs.find( p ) == nb_arcs.end() )
+	    nb_arcs[ p ] = 1;
+	  else
+	    nb_arcs[ p ] += 1;
+	}
+    }
+  ++nb; nbok += nb_arcs.size() == 24 ? 1 : 0;
+  trace.info() << "(" << nbok << "/" << nb << ") "
+	       << "nb_arcs.size() = " << nb_arcs.size() 
+               << " == " << 24 << " (cube has 24 arcs)" << std::endl;
+  for ( auto arc : nb_arcs )
+    {
+      SCell p1 = arc.first.first;
+      SCell p2 = arc.first.second;
+      ++nb; nbok += (K.sCoords( p1 ) - K.sCoords( p2 )).norm1() == 1 ? 1 : 0;
+    }
+  trace.info() << "(" << nbok << "/" << nb << ") "
+	       << "arcs have all norm 1." << std::endl;
+  for ( auto arc : nb_arcs )
+    {
+      SCell p1 = arc.first.first;
+      SCell p2 = arc.first.second;
+      auto  it = nb_arcs.find( std::make_pair( p2, p1 ) );
+      ++nb; nbok += arc.second == 1 ? 1 : 0;
+      ++nb; nbok += ( it != nb_arcs.end() ) && it->second == 1 ? 1 : 0;
+    }
+  trace.info() << "(" << nbok << "/" << nb << ") "
+	       << "arcs are unique and their inverse are unique." << std::endl;
+  trace.endBlock();
+  return nb == nbok;
+}
+
 ///////////////////////////////////////////////////////////////////////////////
 // Standard services - public :
 
@@ -603,7 +680,8 @@ int main( int argc, char** argv )
     && testLightExplicitDigitalSurface()
     && testDigitalSurface<KhalimskySpaceND<2> >()
     && testDigitalSurface<KhalimskySpaceND<3> >()
-    && testDigitalSurface<KhalimskySpaceND<4> >();
+    && testDigitalSurface<KhalimskySpaceND<4> >()
+    && testOrderingDigitalSurfaceFacesAroundVertex();
   trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
   trace.endBlock();
   return res ? 0 : 1;