DGtal  1.5.beta
ConvexCellComplex.h
1 
17 #pragma once
18 
31 #if defined(ConvexCellComplex_RECURSES)
32 #error Recursive header files inclusion detected in ConvexCellComplex.h
33 #else // defined(ConvexCellComplex_RECURSES)
35 #define ConvexCellComplex_RECURSES
36 
37 #if !defined ConvexCellComplex_h
39 #define ConvexCellComplex_h
40 
42 // Inclusions
43 #include <iostream>
44 #include <string>
45 #include <vector>
46 #include "DGtal/base/Common.h"
47 #include "DGtal/kernel/PointVector.h"
48 #include "DGtal/math/linalg/SimpleMatrix.h"
49 #include "DGtal/geometry/volumes/BoundedLatticePolytope.h"
50 
51 namespace DGtal
52 {
54  // template class ConvexCellComplex
84  template < typename TPoint >
86  static const Dimension dimension = TPoint::dimension;
87 
88  typedef TPoint Point;
89  typedef std::size_t Index;
90  typedef std::size_t Size;
91 
92  static const Index INFINITE_CELL = (Index) -1;
93  typedef Index Cell;
94  typedef std::pair< Index, bool > Face;
95  typedef Index Vertex;
96  typedef std::vector< Index > IndexRange;
97  typedef std::vector< Vertex > VertexRange;
98  typedef std::vector< Face > FaceRange;
99 
100  typedef typename Point::Coordinate Scalar;
104 
105  // ----------------------- standard services ---------------------------
106  public:
109 
112  { clear(); }
113 
115  void clear()
116  {
117  cell_faces.clear();
118  cell_vertices.clear();
119  true_face_cell.clear();
120  false_face_cell.clear();
121  true_face_vertices.clear();
122  vertex_position.clear();
123  has_face_geometry = false;
124  true_face_normal.clear();
125  true_face_intercept.clear();
126  }
127 
129  Size nbCells() const
130  { return cell_faces.size(); }
132  Size nbFaces() const
133  { return true_face_cell.size(); }
135  Size nbVertices() const
136  { return vertex_position.size(); }
137 
139  // -------------------- primal structure topology services ---------------------
140  public:
143 
146  { return INFINITE_CELL; }
147 
150  bool isInfinite( const Cell c ) const
151  { return c == INFINITE_CELL; }
152 
153 
156  Face opposite( const Face f ) const
157  { return std::make_pair( f.first, ! f.second ); }
158 
161  const FaceRange& cellFaces( const Cell c ) const
162  { return cell_faces[ c ]; }
163 
169  const VertexRange& cellVertices( const Cell c ) const
170  {
171  ASSERT( c != INFINITE_CELL && c < nbCells() );
172  ASSERT( ! cell_vertices.empty() );
173  if ( cell_vertices.empty() )
174  cell_vertices.resize( nbCells() );
175  if ( cell_vertices[ c ].empty() )
176  computeCellVertices( c );
177  return cell_vertices[ c ];
178  }
179 
181  const std::vector< VertexRange > & allCellVertices() const
182  {
183  if ( cell_vertices.empty() )
184  cell_vertices.resize( nbCells() );
185  for ( Index c = 0; c < nbCells(); ++c )
186  if ( cell_vertices[ c ].empty() )
187  computeCellVertices( c );
188  return cell_vertices;
189  }
190 
193  VertexRange faceVertices( const Face f ) const
194  {
195  if ( f.second ) return true_face_vertices[ f.first ];
196  return VertexRange( true_face_vertices[ f.first ].crbegin(),
197  true_face_vertices[ f.first ].crend() );
198  }
199 
202  Cell faceCell( const Face f ) const
203  {
204  return f.second ? true_face_cell[ f.first ] : false_face_cell[ f.first ];
205  }
206 
212  {
213  VertexRange result;
214  const Cell c = faceCell( f );
215  if ( isInfinite( c ) ) return result;
216  const auto& c_vtcs = cellVertices( c );
217  auto f_vtcs = true_face_vertices[ f.first ];
218  std::sort( f_vtcs.begin(), f_vtcs.end() );
219  for ( auto v : c_vtcs )
220  if ( ! std::binary_search( f_vtcs.cbegin(), f_vtcs.cend(), v ) )
221  result.push_back( v );
222  return result;
223  }
224 
226  // -------------------- primal structure geometry services ---------------------
227  public:
230 
233  std::vector< Point > cellVertexPositions( const Cell c ) const
234  {
235  ASSERT( ! isInfinite( c ) );
236  const auto vtcs = cellVertices( c );
237  std::vector< Point > pts( vtcs.size() );
238  std::transform( vtcs.cbegin(), vtcs.cend(), pts.begin(),
239  [&] ( Vertex v ) { return this->position( v ); } );
240  return pts;
241  }
242 
245  std::vector< Point > faceVertexPositions( const Face f ) const
246  {
247  const auto vtcs = faceVertices( f );
248  std::vector< Point > pts( vtcs.size() );
249  std::transform( vtcs.cbegin(), vtcs.cend(), pts.begin(),
250  [&] ( Vertex v ) { return this->position( v ); } );
251  return pts;
252  }
253 
256  Point position( const Vertex v ) const
257  {
258  return vertex_position[ v ];
259  }
260 
263  RealPoint toReal( const Point p ) const
264  {
265  RealPoint x;
266  for ( Dimension k = 0; k < dimension; ++k )
267  x[ k ] = (double) p[ k ];
268  return x;
269  }
270 
274  template < typename LatticePoint >
275  LatticePoint toLattice( const RealPoint p, double factor = 1.0 ) const
276  {
277  LatticePoint x;
278  for ( Dimension k = 0; k < dimension; ++k )
279  x[ k ] = (typename LatticePoint::Coordinate) round( p[ k ] * factor );
280  return x;
281  }
282 
283 
286  RealPoint cellBarycenter( const Cell c ) const
287  {
288  ASSERT( ! isInfinite( c ) );
289  RealPoint b;
290  const auto& vtcs = cellVertices( c );
291  for ( auto v : vtcs )
292  b += toReal( position( v ) );
293  return b / vtcs.size();
294  }
295 
300  template < typename LatticePolytope >
301  LatticePolytope
302  cellLatticePolytope( const Cell c, const double factor = 1.0 ) const
303  {
304  ASSERT( hasFaceGeometry() );
305  ASSERT( ! isInfinite( c ) );
306  typedef typename LatticePolytope::Point LatticePoint;
307  typedef typename LatticePolytope::HalfSpace HalfSpace;
308  typedef typename LatticePolytope::Domain Domain;
309  const auto vtcs = cellVertices( c );
310  std::vector< LatticePoint > pts;
311  for ( auto v : vtcs )
312  pts.push_back( toLattice< LatticePoint >( position ( v ), factor ) );
313  LatticePoint l = pts[ 0 ];
314  LatticePoint u = pts[ 0 ];
315  for ( const auto& p : pts ) {
316  l = l.inf( p );
317  u = u.sup( p );
318  }
319  Domain domain( l, u );
320  std::vector< HalfSpace > HS;
321  for ( auto f : cellFaces( c ) ) {
322  HS.emplace_back( faceNormal( f ), faceIntercept( f ) );
323  }
324  return LatticePolytope( domain, HS.cbegin(), HS.cend(), false );
325  }
326 
330  const Vector& faceNormal( const Face f ) const
331  {
332  ASSERT( hasFaceGeometry() );
333  ASSERT( f.first < nbFaces() );
334  if ( f.second ) return true_face_normal[ f.first ];
335  else return -true_face_normal[ f.first ];
336  }
337 
341  Scalar faceIntercept( const Face f ) const
342  {
343  ASSERT( hasFaceGeometry() );
344  ASSERT( f.first < nbFaces() );
345  if ( f.second ) return true_face_intercept[ f.first ];
346  else return -true_face_intercept[ f.first ];
347  }
348 
350  bool hasFaceGeometry() const
351  { return has_face_geometry; }
352 
355  {
356  if ( ! hasFaceGeometry() )
358  }
359 
362  {
363  has_face_geometry = false;
364  true_face_normal.clear();
365  true_face_intercept.clear();
366  }
367 
368 
370  // ----------------------- datas for primal structure ----------------------
371  public:
374 
375  // for each cell, gives its faces
376  std::vector< FaceRange > cell_faces;
377  // for each cell, gives its vertices
378  mutable std::vector< VertexRange > cell_vertices;
379  // for each 'true' face, gives its cell (or INFINITE_CELL)
380  std::vector< Cell > true_face_cell;
381  // for each 'false' face, gives its cell (or INFINITE_CELL)
382  std::vector< Cell > false_face_cell;
383  // for each true face, gives its vertices in order.
384  std::vector< VertexRange > true_face_vertices;
385  // for each vertex, gives its position
386  std::vector< Point > vertex_position;
387 
391  std::vector< Vector > true_face_normal;
393  std::vector< Scalar > true_face_intercept;
394 
396  // ----------------------- Interface ----------------------
397  public:
400 
405  void selfDisplay ( std::ostream & out ) const
406  {
407  out << "[ConvexCellComplex<" << dimension << ">"
408  << " #C=" << nbCells()
409  << " #F=" << nbFaces()
410  << " #V=" << nbVertices();
411  if ( hasFaceGeometry() ) out << " hasFaceGeometry";
412  out << " ]";
413  }
414 
419  bool isValid() const
420  {
421  return nbCells() != 0;
422  }
423 
425  // ----------------------- protected services ----------------------
426  protected:
429 
432  void computeCellVertices( const Cell c ) const
433  {
434  std::set< Vertex > vertices;
435  for ( auto f : cell_faces[ c ] )
436  vertices.insert( true_face_vertices[ f.first ].cbegin(),
437  true_face_vertices[ f.first ].cend() );
438  cell_vertices[ c ] = VertexRange( vertices.cbegin(), vertices.cend() );
439  }
440 
443  {
444  true_face_normal .resize( nbFaces() );
445  true_face_intercept.resize( nbFaces() );
446  // Compute normals and intercepts
447  for ( Index f = 0; f < nbFaces(); ++f ) {
448  std::tie( true_face_normal[ f ], true_face_intercept[ f ] )
450  if ( true_face_normal[ f ] == Vector::zero )
451  trace.error() << "[ConvexCellComplex::computeFaceGeometry]"
452  << " null normal vector at face " << f << std::endl;
453  }
454  // Orient them consistently
455  for ( Index c = 0; c < nbCells(); c++ ) {
456  //const auto& c_vtcs = cellVertices( c ); //not used
457  for ( auto f : cellFaces( c ) ) {
458  if ( ! f.second ) continue; // process only 'true' faces
459  const auto ov = faceComplementVertices( f );
460  ASSERT( ! ov.empty() );
461  const Vector N = true_face_normal [ f.first ];
462  const Scalar nu = true_face_intercept[ f.first ];
463  const Scalar iv = N.dot( position( ov.front() ) );
464  if ( ( iv - nu ) > 0 ) {
465  // Should be inside
466  true_face_normal [ f.first ] = -N;
467  true_face_intercept[ f.first ] = -nu;
468  std::reverse( true_face_vertices[ f.first ].begin(),
469  true_face_vertices[ f.first ].end() );
470  }
471  reorderFaceVertices( f.first );
472  }
473  }
474  has_face_geometry = true;
475  }
476 
479  std::pair< Vector, Scalar >
480  computeHalfSpace( const VertexRange& v ) const
481  {
483  Matrix A;
484  Vector N;
485  Scalar c;
486  for ( int i = 1; i < dimension; i++ )
487  for ( int j = 0; j < dimension; j++ )
488  A.setComponent( i-1, j, position( v[ i ] )[ j ] - position( v[ 0 ] )[ j ] );
489  for ( int j = 0; j < dimension; j++ )
490  N[ j ] = A.cofactor( dimension-1, j );
491  c = N.dot( position( v[ 0 ] ) );
492  return std::make_pair( N, c );
493  }
494 
504  {
505  VertexRange& result = true_face_vertices[ f ];
506  Vector N = true_face_normal [ f ];
507  // Sort a facet such that its points, taken in order, have
508  // always the same orientation of the facet. More precisely,
509  // facets span a `dimension-1` vector space. There are thus
510  // dimension-2 fixed points, and the last ones (at least two)
511  // may be reordered.
512  VertexRange splx( dimension );
513  for ( Dimension k = 0; k < dimension-2; k++ )
514  splx[ k ] = result[ k ];
515  std::sort( result.begin()+dimension-2, result.end(),
516  [&] ( Index i, Index j )
517  {
518  splx[ dimension-2 ] = i;
519  splx[ dimension-1 ] = j;
520  const auto H = computeHalfSpace( splx );
521  const auto orient = N.dot( H.first );
522  return orient > 0;
523  } );
524  }
525 
527 
528  }; // class ConvexCellComplex
529 
536  template <typename TPoint>
537  std::ostream&
538  operator<< ( std::ostream & out, const ConvexCellComplex<TPoint> & object )
539  {
540  object.selfDisplay( out );
541  return out;
542  }
543 
544 } // namespace DGtal
545 
546 #endif // !defined ConvexCellComplex_h
547 
548 #undef ConvexCellComplex_RECURSES
549 #endif // else defined(ConvexCellComplex_RECURSES)
550 
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:593
auto dot(const PointVector< dim, OtherComponent, OtherStorage > &v) const -> decltype(DGtal::dotProduct(*this, v))
Dot product with a PointVector.
static Self zero
Static const for zero PointVector.
Definition: PointVector.h:1595
Aim: implements basic MxN Matrix services (M,N>=1).
Definition: SimpleMatrix.h:76
std::ostream & error()
DGtal is the top-level namespace which contains all DGtal functions and types.
std::ostream & operator<<(std::ostream &out, const ATu0v1< TKSpace, TLinearAlgebra > &object)
DGtal::uint32_t Dimension
Definition: Common.h:136
Trace trace
Definition: Common.h:153
Aim: represents a d-dimensional complex in a d-dimensional space with the following properties and re...
const Vector & faceNormal(const Face f) const
std::vector< Point > cellVertexPositions(const Cell c) const
std::vector< VertexRange > true_face_vertices
Tells if the face geometry has been computed.
Cell faceCell(const Face f) const
std::vector< Scalar > true_face_intercept
Contains the intercept of each 'true' face.
std::vector< Cell > false_face_cell
Tells if the face geometry has been computed.
LatticePolytope cellLatticePolytope(const Cell c, const double factor=1.0) const
std::vector< Cell > true_face_cell
Tells if the face geometry has been computed.
bool has_face_geometry
Tells if the face geometry has been computed.
void selfDisplay(std::ostream &out) const
const VertexRange & cellVertices(const Cell c) const
std::vector< Index > IndexRange
std::pair< Index, bool > Face
static const Index INFINITE_CELL
VertexRange faceComplementVertices(const Face f) const
Face opposite(const Face f) const
Point position(const Vertex v) const
void clear()
Clears the complex (as if it was just default constructed).
Scalar faceIntercept(const Face f) const
std::vector< FaceRange > cell_faces
Tells if the face geometry has been computed.
PointVector< dimension, double > RealPoint
std::vector< Vertex > VertexRange
static const Dimension dimension
std::vector< Face > FaceRange
PointVector< dimension, Scalar > Vector
VertexRange faceVertices(const Face f) const
RealPoint toReal(const Point p) const
std::vector< Point > faceVertexPositions(const Face f) const
RealPoint cellBarycenter(const Cell c) const
std::pair< Vector, Scalar > computeHalfSpace(const VertexRange &v) const
void requireFaceGeometry()
Forces the computation of face geometry.
void computeCellVertices(const Cell c) const
const FaceRange & cellFaces(const Cell c) const
bool isInfinite(const Cell c) const
ConvexCellComplex()
Defaut constructor.
LatticePoint toLattice(const RealPoint p, double factor=1.0) const
std::vector< VertexRange > cell_vertices
Tells if the face geometry has been computed.
std::vector< Vector > true_face_normal
Contains the outward oriented normal of each 'true' face.
std::vector< Point > vertex_position
Tells if the face geometry has been computed.
const std::vector< VertexRange > & allCellVertices() const
void computeFaceGeometry()
Computes for each face its outward oriented normal vector.
void unrequireFaceGeometry()
Release ressources allocated to face geometry.
PointVector< dimension, double > RealVector
MyPointD Point
Definition: testClone2.cpp:383
ch reverse()
Domain domain
HyperRectDomain< Space > Domain