DGtal  1.4.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)
34 
35 #define ConvexCellComplex_RECURSES
36 
37 #if !defined ConvexCellComplex_h
38 
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 );
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 
DGtal::ConvexCellComplex::cellLatticePolytope
LatticePolytope cellLatticePolytope(const Cell c, const double factor=1.0) const
Definition: ConvexCellComplex.h:302
DGtal::ConvexCellComplex::Index
std::size_t Index
Definition: ConvexCellComplex.h:89
DGtal::ConvexCellComplex::cellVertexPositions
std::vector< Point > cellVertexPositions(const Cell c) const
Definition: ConvexCellComplex.h:233
DGtal::ConvexCellComplex::position
Point position(const Vertex v) const
Definition: ConvexCellComplex.h:256
DGtal::ConvexCellComplex::allCellVertices
const std::vector< VertexRange > & allCellVertices() const
Definition: ConvexCellComplex.h:181
DGtal::ConvexCellComplex::cellVertices
const VertexRange & cellVertices(const Cell c) const
Definition: ConvexCellComplex.h:169
DGtal::PointVector::zero
static Self zero
Static const for zero PointVector.
Definition: PointVector.h:1595
DGtal::ConvexCellComplex::faceCell
Cell faceCell(const Face f) const
Definition: ConvexCellComplex.h:202
DGtal::ConvexCellComplex::Scalar
Point::Coordinate Scalar
Definition: ConvexCellComplex.h:100
DGtal::ConvexCellComplex::VertexRange
std::vector< Vertex > VertexRange
Definition: ConvexCellComplex.h:97
DGtal::ConvexCellComplex::hasFaceGeometry
bool hasFaceGeometry() const
Definition: ConvexCellComplex.h:350
DGtal::ConvexCellComplex::faceIntercept
Scalar faceIntercept(const Face f) const
Definition: ConvexCellComplex.h:341
DGtal::ConvexCellComplex::isValid
bool isValid() const
Definition: ConvexCellComplex.h:419
DGtal::HyperRectDomain< Space >
DGtal::ConvexCellComplex::unrequireFaceGeometry
void unrequireFaceGeometry()
Release ressources allocated to face geometry.
Definition: ConvexCellComplex.h:361
DGtal::Trace::error
std::ostream & error()
DGtal::ConvexCellComplex::Vector
PointVector< dimension, Scalar > Vector
Definition: ConvexCellComplex.h:101
DGtal::ConvexCellComplex::true_face_normal
std::vector< Vector > true_face_normal
Contains the outward oriented normal of each 'true' face.
Definition: ConvexCellComplex.h:391
DGtal::ConvexCellComplex::toReal
RealPoint toReal(const Point p) const
Definition: ConvexCellComplex.h:263
DGtal::trace
Trace trace
Definition: Common.h:154
DGtal::ConvexCellComplex::Point
TPoint Point
Definition: ConvexCellComplex.h:88
DGtal::Dimension
DGtal::uint32_t Dimension
Definition: Common.h:137
DGtal::ConvexCellComplex::requireFaceGeometry
void requireFaceGeometry()
Forces the computation of face geometry.
Definition: ConvexCellComplex.h:354
DGtal::ConvexCellComplex::toLattice
LatticePoint toLattice(const RealPoint p, double factor=1.0) const
Definition: ConvexCellComplex.h:275
DGtal::ConvexCellComplex::Vertex
Index Vertex
Definition: ConvexCellComplex.h:95
DGtal::ConvexCellComplex::false_face_cell
std::vector< Cell > false_face_cell
Tells if the face geometry has been computed.
Definition: ConvexCellComplex.h:382
DGtal::ConvexCellComplex
Aim: represents a d-dimensional complex in a d-dimensional space with the following properties and re...
Definition: ConvexCellComplex.h:85
DGtal::ConvexCellComplex::ConvexCellComplex
ConvexCellComplex()
Defaut constructor.
Definition: ConvexCellComplex.h:111
DGtal::ConvexCellComplex::reorderFaceVertices
void reorderFaceVertices(Index f)
Definition: ConvexCellComplex.h:503
DGtal::operator<<
std::ostream & operator<<(std::ostream &out, const ATu0v1< TKSpace, TLinearAlgebra > &object)
DGtal::ConvexCellComplex::faceVertexPositions
std::vector< Point > faceVertexPositions(const Face f) const
Definition: ConvexCellComplex.h:245
DGtal::ConvexCellComplex::faceComplementVertices
VertexRange faceComplementVertices(const Face f) const
Definition: ConvexCellComplex.h:211
DGtal::ConvexCellComplex::nbFaces
Size nbFaces() const
Definition: ConvexCellComplex.h:132
DGtal::SimpleMatrix
Aim: implements basic MxN Matrix services (M,N>=1).
Definition: SimpleMatrix.h:75
DGtal::ConvexCellComplex::Cell
Index Cell
Definition: ConvexCellComplex.h:93
DGtal::ConvexCellComplex::RealVector
PointVector< dimension, double > RealVector
Definition: ConvexCellComplex.h:103
DGtal::ConvexCellComplex::cellFaces
const FaceRange & cellFaces(const Cell c) const
Definition: ConvexCellComplex.h:161
reverse
ch reverse()
DGtal::ConvexCellComplex::vertex_position
std::vector< Point > vertex_position
Tells if the face geometry has been computed.
Definition: ConvexCellComplex.h:386
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::ConvexCellComplex::computeHalfSpace
std::pair< Vector, Scalar > computeHalfSpace(const VertexRange &v) const
Definition: ConvexCellComplex.h:480
DGtal::ConvexCellComplex::cell_vertices
std::vector< VertexRange > cell_vertices
Tells if the face geometry has been computed.
Definition: ConvexCellComplex.h:378
Domain
HyperRectDomain< Space > Domain
Definition: testSimpleRandomAccessRangeFromPoint.cpp:44
DGtal::ConvexCellComplex::true_face_intercept
std::vector< Scalar > true_face_intercept
Contains the intercept of each 'true' face.
Definition: ConvexCellComplex.h:393
DGtal::ConvexCellComplex::has_face_geometry
bool has_face_geometry
Tells if the face geometry has been computed.
Definition: ConvexCellComplex.h:389
DGtal::ConvexCellComplex::dimension
static const Dimension dimension
Definition: ConvexCellComplex.h:86
DGtal::ConvexCellComplex::infiniteCell
Cell infiniteCell() const
Definition: ConvexCellComplex.h:145
DGtal::ConvexCellComplex::true_face_cell
std::vector< Cell > true_face_cell
Tells if the face geometry has been computed.
Definition: ConvexCellComplex.h:380
DGtal::ConvexCellComplex::true_face_vertices
std::vector< VertexRange > true_face_vertices
Tells if the face geometry has been computed.
Definition: ConvexCellComplex.h:384
DGtal::ConvexCellComplex::Size
std::size_t Size
Definition: ConvexCellComplex.h:90
DGtal::ConvexCellComplex::INFINITE_CELL
static const Index INFINITE_CELL
Definition: ConvexCellComplex.h:92
DGtal::ConvexCellComplex::FaceRange
std::vector< Face > FaceRange
Definition: ConvexCellComplex.h:98
DGtal::ConvexCellComplex::computeCellVertices
void computeCellVertices(const Cell c) const
Definition: ConvexCellComplex.h:432
DGtal::ConvexCellComplex::isInfinite
bool isInfinite(const Cell c) const
Definition: ConvexCellComplex.h:150
DGtal::ConvexCellComplex::Face
std::pair< Index, bool > Face
Definition: ConvexCellComplex.h:94
DGtal::ConvexCellComplex::faceVertices
VertexRange faceVertices(const Face f) const
Definition: ConvexCellComplex.h:193
DGtal::PointVector::dot
auto dot(const PointVector< dim, OtherComponent, OtherStorage > &v) const -> decltype(DGtal::dotProduct(*this, v))
Dot product with a PointVector.
domain
Domain domain
Definition: testProjection.cpp:88
DGtal::PointVector
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:165
DGtal::ConvexCellComplex::faceNormal
const Vector & faceNormal(const Face f) const
Definition: ConvexCellComplex.h:330
DGtal::ConvexCellComplex::clear
void clear()
Clears the complex (as if it was just default constructed).
Definition: ConvexCellComplex.h:115
DGtal::ConvexCellComplex::cellBarycenter
RealPoint cellBarycenter(const Cell c) const
Definition: ConvexCellComplex.h:286
DGtal::ConvexCellComplex::nbCells
Size nbCells() const
Definition: ConvexCellComplex.h:129
DGtal::ConvexCellComplex::IndexRange
std::vector< Index > IndexRange
Definition: ConvexCellComplex.h:96
DGtal::ConvexCellComplex::nbVertices
Size nbVertices() const
Definition: ConvexCellComplex.h:135
Point
MyPointD Point
Definition: testClone2.cpp:383
DGtal::ConvexCellComplex::selfDisplay
void selfDisplay(std::ostream &out) const
Definition: ConvexCellComplex.h:405
DGtal::ConvexCellComplex::cell_faces
std::vector< FaceRange > cell_faces
Tells if the face geometry has been computed.
Definition: ConvexCellComplex.h:376
DGtal::ConvexCellComplex::RealPoint
PointVector< dimension, double > RealPoint
Definition: ConvexCellComplex.h:102
DGtal::ConvexCellComplex::opposite
Face opposite(const Face f) const
Definition: ConvexCellComplex.h:156
DGtal::ConvexCellComplex::computeFaceGeometry
void computeFaceGeometry()
Computes for each face its outward oriented normal vector.
Definition: ConvexCellComplex.h:442