DGtal  1.2.0
ConvexityHelper.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file ConvexityHelper.ih
19  * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20  * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
21  *
22  * @date 2020/01/02
23  *
24  * Implementation of inline methods defined in ConvexityHelper.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 #include <string>
33 #include <sstream>
34 //////////////////////////////////////////////////////////////////////////////
35 
36 ///////////////////////////////////////////////////////////////////////////////
37 // IMPLEMENTATION of inline methods.
38 ///////////////////////////////////////////////////////////////////////////////
39 
40 ///////////////////////////////////////////////////////////////////////////////
41 // ----------------------- Standard services ------------------------------
42 
43 //-----------------------------------------------------------------------------
44 template < int dim, typename TInteger, typename TInternalInteger >
45 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
46 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeLatticePolytope
47 ( const std::vector< Point >& input_points,
48  bool remove_duplicates,
49  bool make_minkowski_summable )
50 {
51  typedef typename LatticePolytope::Domain Domain;
52  typedef typename LatticePolytope::HalfSpace PolytopeHalfSpace;
53  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
54  typedef typename ConvexHull::HalfSpace ConvexHullHalfSpace;
55  typedef typename ConvexHull::Ridge Ridge;
56  if ( input_points.empty() ) return LatticePolytope();
57  // Compute domain
58  Point l = input_points[ 0 ];
59  Point u = input_points[ 0 ];
60  for ( const auto& p : input_points ) {
61  l = l.inf( p );
62  u = u.sup( p );
63  }
64  Domain domain( l, u );
65  // Compute convex hull
66  ConvexHull hull;
67  hull.setInput( input_points, remove_duplicates );
68  const auto target = ( make_minkowski_summable && dimension == 3 )
69  ? ConvexHull::Status::VerticesCompleted
70  : ConvexHull::Status::FacetsCompleted;
71  bool ok = hull.computeConvexHull( target );
72  if ( ! ok ) return LatticePolytope();
73  // Initialize polytope
74  std::vector< ConvexHullHalfSpace > HS;
75  std::vector< PolytopeHalfSpace > PHS;
76  hull.getFacetHalfSpaces( HS );
77  PHS.reserve( HS.size() );
78  for ( auto& H : HS ) {
79  Vector N;
80  Integer nu;
81  for ( Dimension i = 0; i < dim; ++i )
82  N[ i ] = (Integer) H.internalNormal()[ i ];
83  nu = (Integer) H.internalIntercept();
84  PHS.emplace_back( N, nu );
85  }
86  if ( make_minkowski_summable && dimension >= 4 )
87  trace.warning() << "[ConvexityHelper::computeLatticePolytope]"
88  << " Not implemented starting from dimension 4."
89  << std::endl;
90  if ( make_minkowski_summable && dimension == 3 )
91  {
92  // Compute ridge vertices to add edge constraints.
93  std::vector< Point > positions;
94  std::vector< IndexRange > facet_vertices;
95  std::vector< IndexRange > ridge_vertices;
96  std::map< Ridge, Index > ridge2index;
97  hull.getVertexPositions( positions );
98  computeFacetAndRidgeVertices( hull, facet_vertices,
99  ridge2index, ridge_vertices );
100  for ( auto p : ridge2index ) {
101  const auto r = p.first;
102  // Copy by value since PHS may be reallocated during the iteration.
103  const auto U = PHS[ r.first ].N; // normal of facet 1
104  const auto V = PHS[ r.second ].N; // normal of facet 2
105  const auto& S = ridge_vertices[ p.second ]; // vertices along facets 1, 2
106  ASSERT( S.size() == 2 && "Invalid ridge" );
107  const auto& P0 = positions[ S[ 0 ] ];
108  const auto& P1 = positions[ S[ 1 ] ];
109  auto E = P1 - P0; // edge 1, 2
110  const auto UxV =
111  detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
112  ::crossProduct( U, V ); // parallel to E
113  ASSERT( E.dot( UxV ) != 0 && "Invalid E / UxV" );
114  if ( E.dot( UxV ) <= 0 ) E = -E; // force correct orientation
115  const auto E1 =
116  detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
117  ::crossProduct( U, E ); // edge on facet 1
118  const auto E2 =
119  detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
120  ::crossProduct( E, V ); // edge on facet 2
121  ASSERT( E1.dot( U ) == 0 && "Invalid E1 / U" );
122  ASSERT( E1.dot( V ) < 0 && "Invalid E1 / V" );
123  ASSERT( E2.dot( V ) == 0 && "Invalid E1 / V" );
124  ASSERT( E2.dot( U ) < 0 && "Invalid E1 / U" );
125  for ( Dimension k = 0; k < dimension; ++k ) {
126  const auto W = U[ k ] * V - V[ k ] * U;
127  const auto nn1 = W.dot( E1 );
128  const auto nn2 = W.dot( E2 );
129  if ( nn1 > 0 && nn2 > 0 ) {
130  PHS.emplace_back( -W, -W.dot( P0 ) );
131  ASSERT( E1.dot(-W ) < 0 && "Invalid E1 /-W" );
132  ASSERT( E2.dot(-W ) < 0 && "Invalid E2 /-W" );
133  }
134  else if ( nn1 < 0 && nn2 < 0 ) {
135  PHS.emplace_back( W, W.dot( P0 ) );
136  ASSERT( E1.dot( W ) < 0 && "Invalid E1 / W" );
137  ASSERT( E2.dot( W ) < 0 && "Invalid E2 / W" );
138  }
139  }
140  }
141  }
142  return LatticePolytope( domain, PHS.cbegin(), PHS.cend(),
143  make_minkowski_summable && ( dimension <= 3 ), true );
144 }
145 
146 //-----------------------------------------------------------------------------
147 template < int dim, typename TInteger, typename TInternalInteger >
148 template < typename TSurfaceMesh >
149 bool
150 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
151 ( TSurfaceMesh& mesh,
152  const std::vector< Point >& input_points,
153  bool remove_duplicates )
154 {
155  typedef TSurfaceMesh SurfaceMesh;
156  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
157  typedef typename ConvexHull::IndexRange IndexRange;
158  ConvexHull hull;
159  hull.setInput( input_points, remove_duplicates );
160  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
161  if ( !ok ) return false;
162  std::vector< RealPoint > positions;
163  hull.getVertexPositions( positions );
164  std::vector< IndexRange > faces;
165  hull.getFacetVertices( faces );
166  mesh = SurfaceMesh( positions.cbegin(), positions.cend(),
167  faces.cbegin(), faces.cend() );
168  return true;
169 }
170 
171 //-----------------------------------------------------------------------------
172 template < int dim, typename TInteger, typename TInternalInteger >
173 bool
174 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
175 ( PolygonalSurface< Point >& polysurf,
176  const std::vector< Point >& input_points,
177  bool remove_duplicates )
178 {
179  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
180  typedef typename ConvexHull::IndexRange IndexRange;
181  ConvexHull hull;
182  hull.setInput( input_points, remove_duplicates );
183  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
184  if ( !ok ) return false;
185  std::vector< Point > positions;
186  hull.getVertexPositions( positions );
187  std::vector< IndexRange > faces;
188  hull.getFacetVertices( faces );
189  // build polygonal surface
190  polysurf.clear();
191  for ( auto p : positions ) polysurf.addVertex( p );
192  for ( auto f : faces ) polysurf.addPolygonalFace( f );
193  return polysurf.build();
194 }
195 
196 //-----------------------------------------------------------------------------
197 template < int dim, typename TInteger, typename TInternalInteger >
198 bool
199 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullCellComplex
200 ( ConvexCellComplex< Point >& cell_complex,
201  const std::vector< Point >& input_points,
202  bool remove_duplicates )
203 {
204  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
205  typedef typename ConvexHull::IndexRange IndexRange;
206  typedef typename ConvexCellComplex< Point >::FaceRange FaceRange;
207  ConvexHull hull;
208  hull.setInput( input_points, remove_duplicates );
209  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
210  cell_complex.clear();
211  if ( ! ok ) return false;
212  // Build complex, only 1 finite cell and as many faces as convex hull facets.
213  // Taking care of faces for each cell (here one cell borders all faces).
214  std::vector< IndexRange > faces;
215  hull.getFacetVertices( faces );
216  FaceRange all_faces;
217  for ( Index i = 0; i < faces.size(); i++ )
218  all_faces.push_back( { i, true } );
219  cell_complex.cell_faces.push_back( all_faces );
220  // Vertices of this unique cell will be computed lazily on request.
221  // Taking care of each face.
222  for ( const auto& vtcs : faces )
223  {
224  // every inner face borders cell 0
225  cell_complex.true_face_cell.push_back( 0 );
226  // every outer face borders the infinite cell
227  cell_complex.false_face_cell.push_back( cell_complex.infiniteCell() );
228  }
229  // Taking care of vertices (in consistent order) of every face
230  cell_complex.true_face_vertices.swap( faces );
231  // Taking care of vertex positions
232  hull.getVertexPositions( cell_complex.vertex_position );
233  return true;
234 }
235 
236 //-----------------------------------------------------------------------------
237 template < int dim, typename TInteger, typename TInternalInteger >
238 bool
239 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeDelaunayCellComplex
240 ( ConvexCellComplex< Point >& cell_complex,
241  const std::vector< Point >& input_points,
242  bool remove_duplicates )
243 {
244  typedef QuickHull< LatticeDelaunayKernel > Delaunay;
245  typedef typename Delaunay::Ridge Ridge;
246  typedef typename ConvexCellComplex< Point >::FaceRange FaceRange;
247 
248  Delaunay del;
249  del.setInput( input_points, remove_duplicates );
250  bool ok = del.computeConvexHull( Delaunay::Status::VerticesCompleted );
251  cell_complex.clear();
252  if ( ! ok ) return false;
253 
254  // Build complex, as many maximal cells as convex hull facets.
255  // convex hull facet -> cell of complex
256  // convex hull ridge -> face of complex
257  // (1) Get cell vertices, count ridges/faces and compute their vertices
258  std::map< Ridge, Index > r2f;
259  computeFacetAndRidgeVertices( del,
260  cell_complex.cell_vertices,
261  r2f,
262  cell_complex.true_face_vertices );
263  // (2) assign ridges/faces to cell and conversely
264  const Index nb_r = r2f.size();
265  cell_complex.true_face_cell .resize( nb_r, cell_complex.infiniteCell() );
266  cell_complex.false_face_cell.resize( nb_r, cell_complex.infiniteCell() );
267  cell_complex.true_face_vertices.resize( nb_r );
268  for ( Index cur_f = 0; cur_f < del.nbFiniteFacets(); ++cur_f ) {
269  const auto& facet = del.facets[ cur_f ];
270  FaceRange current_faces;
271  for ( auto neigh_f : facet.neighbors ) {
272  const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
273  const bool pos = cur_f < neigh_f;
274  const Index cur_r = r2f[ r ];
275  cell_complex.true_face_cell [ cur_r ] = r.first;
276  if ( r.second >= del.nbFiniteFacets() )
277  cell_complex.false_face_cell[ cur_r ] = cell_complex.infiniteCell();
278  else
279  cell_complex.false_face_cell[ cur_r ] = r.second;
280  current_faces.emplace_back( cur_r, pos );
281  }
282  cell_complex.cell_faces.push_back( current_faces );
283  }
284  // (3) Takes care of vertex positions
285  del.getVertexPositions( cell_complex.vertex_position );
286  return true;
287 }
288 
289 
290 //-----------------------------------------------------------------------------
291 template < int dim, typename TInteger, typename TInternalInteger >
292 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::RationalPolytope
293 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeRationalPolytope
294 ( const std::vector< RealPoint >& input_points,
295  Integer denominator,
296  bool remove_duplicates,
297  bool make_minkowski_summable )
298 {
299  if ( denominator < 1 )
300  trace.error() << "Invalid denominator " << denominator
301  << ". Should be greater or equal to 1." << std::endl;
302  typedef typename RationalPolytope::Domain Domain;
303  typedef typename RationalPolytope::HalfSpace PolytopeHalfSpace;
304  typedef QuickHull< RealConvexHullKernel > ConvexHull;
305  typedef typename ConvexHull::HalfSpace ConvexHullHalfSpace;
306  typedef typename ConvexHull::Ridge Ridge;
307  if ( input_points.empty() ) return RationalPolytope();
308  // Compute convex hull
309  ConvexHull hull( denominator );
310  hull.setInput( input_points, remove_duplicates );
311  const auto target = ( make_minkowski_summable && dimension == 3 )
312  ? ConvexHull::Status::VerticesCompleted
313  : ConvexHull::Status::FacetsCompleted;
314  bool ok = hull.computeConvexHull( target );
315  if ( ! ok ) return RationalPolytope();
316  // Compute domain (as a lattice domain)
317  auto l = hull.points[ 0 ];
318  auto u = hull.points[ 0 ];
319  for ( const auto& p : hull.points ) {
320  l = l.inf( p );
321  u = u.sup( p );
322  }
323  Domain domain( l, u );
324  trace.info() << "Domain l=" << l << " u=" << u << std::endl;
325  // Initialize polytope
326  std::vector< ConvexHullHalfSpace > HS;
327  std::vector< PolytopeHalfSpace > PHS;
328  hull.getFacetHalfSpaces( HS );
329  PHS.reserve( HS.size() );
330  for ( auto& H : HS ) {
331  Vector N;
332  Integer nu;
333  for ( Dimension i = 0; i < dim; ++i )
334  N[ i ] = (Integer) H.internalNormal()[ i ];
335  nu = (Integer) H.internalIntercept();
336  PHS.emplace_back( N, nu );
337  }
338  if ( make_minkowski_summable && dimension >= 4 )
339  trace.warning() << "[ConvexityHelper::computeRationalPolytope]"
340  << " Not implemented starting from dimension 4."
341  << std::endl;
342  if ( make_minkowski_summable && dimension == 3 )
343  {
344  // Compute ridge vertices to add edge constraints.
345  std::vector< Point > positions;
346  std::vector< IndexRange > facet_vertices;
347  std::vector< IndexRange > ridge_vertices;
348  std::map< Ridge, Index > ridge2index;
349  hull.getVertexPositions( positions );
350  computeFacetAndRidgeVertices( hull, facet_vertices,
351  ridge2index, ridge_vertices );
352  for ( auto p : ridge2index ) {
353  const auto r = p.first;
354  // Copy by value since PHS may be reallocated during the iteration.
355  const auto U = PHS[ r.first ].N; // normal of facet 1
356  const auto V = PHS[ r.second ].N; // normal of facet 2
357  const auto& S = ridge_vertices[ p.second ]; // vertices along facets 1, 2
358  ASSERT( S.size() == 2 && "Invalid ridge" );
359  const auto& P0 = positions[ S[ 0 ] ];
360  const auto& P1 = positions[ S[ 1 ] ];
361  auto E = P1 - P0; // edge 1, 2
362  const auto UxV =
363  detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
364  ::crossProduct( U, V ); // parallel to E
365  ASSERT( E.dot( UxV ) != 0 && "Invalid E / UxV" );
366  if ( E.dot( UxV ) <= 0 ) E = -E; // force correct orientation
367  const auto E1 =
368  detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
369  ::crossProduct( U, E ); // edge on facet 1
370  const auto E2 =
371  detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
372  ::crossProduct( E, V ); // edge on facet 2
373  ASSERT( E1.dot( U ) == 0 && "Invalid E1 / U" );
374  ASSERT( E1.dot( V ) < 0 && "Invalid E1 / V" );
375  ASSERT( E2.dot( V ) == 0 && "Invalid E1 / V" );
376  ASSERT( E2.dot( U ) < 0 && "Invalid E1 / U" );
377  for ( Dimension k = 0; k < dimension; ++k ) {
378  const auto W = U[ k ] * V - V[ k ] * U;
379  const auto nn1 = W.dot( E1 );
380  const auto nn2 = W.dot( E2 );
381  if ( nn1 > 0 && nn2 > 0 ) {
382  PHS.emplace_back( -W, -W.dot( P0 ) );
383  ASSERT( E1.dot(-W ) < 0 && "Invalid E1 /-W" );
384  ASSERT( E2.dot(-W ) < 0 && "Invalid E2 /-W" );
385  }
386  else if ( nn1 < 0 && nn2 < 0 ) {
387  PHS.emplace_back( W, W.dot( P0 ) );
388  ASSERT( E1.dot( W ) < 0 && "Invalid E1 / W" );
389  ASSERT( E2.dot( W ) < 0 && "Invalid E2 / W" );
390  }
391  }
392  }
393  }
394  return RationalPolytope( denominator, domain, PHS.cbegin(), PHS.cend(),
395  make_minkowski_summable && ( dimension <= 3 ), true );
396 }
397 
398 
399 //-----------------------------------------------------------------------------
400 template < int dim, typename TInteger, typename TInternalInteger >
401 template < typename TSurfaceMesh >
402 bool
403 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
404 ( TSurfaceMesh& mesh,
405  const std::vector< RealPoint >& input_points,
406  double precision,
407  bool remove_duplicates )
408 {
409  typedef TSurfaceMesh SurfaceMesh;
410  typedef QuickHull< RealConvexHullKernel > ConvexHull;
411  typedef typename ConvexHull::IndexRange IndexRange;
412  ConvexHull hull( precision );
413  hull.setInput( input_points, remove_duplicates );
414  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
415  if ( !ok ) return false;
416  std::vector< RealPoint > positions;
417  hull.getVertexPositions( positions );
418  std::vector< IndexRange > faces;
419  hull.getFacetVertices( faces );
420  mesh = SurfaceMesh( positions.cbegin(), positions.cend(),
421  faces.cbegin(), faces.cend() );
422  return true;
423 }
424 
425 
426 //-----------------------------------------------------------------------------
427 template < int dim, typename TInteger, typename TInternalInteger >
428 bool
429 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
430 ( PolygonalSurface< RealPoint >& polysurf,
431  const std::vector< RealPoint >& input_points,
432  double precision,
433  bool remove_duplicates )
434 {
435  typedef QuickHull< RealConvexHullKernel > ConvexHull;
436  typedef typename ConvexHull::IndexRange IndexRange;
437  ConvexHull hull( precision );
438  hull.setInput( input_points, remove_duplicates );
439  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
440  if ( !ok ) return false;
441  std::vector< RealPoint > positions;
442  hull.getVertexPositions( positions );
443  std::vector< IndexRange > faces;
444  hull.getFacetVertices( faces );
445  // build polygonal surface
446  polysurf.clear();
447  for ( auto p : positions ) polysurf.addVertex( p );
448  for ( auto f : faces ) polysurf.addPolygonalFace( f );
449  return polysurf.build();
450 }
451 
452 //-----------------------------------------------------------------------------
453 template < int dim, typename TInteger, typename TInternalInteger >
454 bool
455 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullCellComplex
456 ( ConvexCellComplex< RealPoint >& cell_complex,
457  const std::vector< RealPoint >& input_points,
458  double precision,
459  bool remove_duplicates )
460 {
461  typedef QuickHull< RealConvexHullKernel > ConvexHull;
462  typedef typename ConvexHull::IndexRange IndexRange;
463  typedef typename ConvexCellComplex< RealPoint >::FaceRange FaceRange;
464  ConvexHull hull( precision );
465  hull.setInput( input_points, remove_duplicates );
466  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
467  cell_complex.clear();
468  if ( ! ok ) return false;
469  // Build complex, only 1 finite cell and as many faces as convex hull facets.
470  // Taking care of faces for each cell (here one cell borders all faces).
471  std::vector< IndexRange > faces;
472  hull.getFacetVertices( faces );
473  FaceRange all_faces;
474  for ( Index i = 0; i < faces.size(); i++ )
475  all_faces.push_back( { i, true } );
476  cell_complex.cell_faces.push_back( all_faces );
477  // Vertices of this unique cell will be computed lazily on request.
478  // Taking care of each face.
479  for ( const auto& vtcs : faces )
480  {
481  // every inner face borders cell 0
482  cell_complex.true_face_cell.push_back( 0 );
483  // every outer face borders the infinite cell
484  cell_complex.false_face_cell.push_back( cell_complex.infiniteCell() );
485  }
486  // Taking care of vertices (in consistent order) of every face
487  cell_complex.true_face_vertices.swap( faces );
488  // Taking care of vertex positions
489  hull.getVertexPositions( cell_complex.vertex_position );
490  return true;
491 }
492 
493 
494 //-----------------------------------------------------------------------------
495 template < int dim, typename TInteger, typename TInternalInteger >
496 bool
497 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeDelaunayCellComplex
498 ( ConvexCellComplex< RealPoint >& cell_complex,
499  const std::vector< RealPoint >& input_points,
500  double precision,
501  bool remove_duplicates )
502 {
503  typedef QuickHull< RealDelaunayKernel > Delaunay;
504  typedef typename Delaunay::Ridge Ridge;
505  typedef typename ConvexCellComplex< RealPoint >::FaceRange FaceRange;
506 
507  Delaunay del( precision );
508  del.setInput( input_points, remove_duplicates );
509  bool ok = del.computeConvexHull( Delaunay::Status::VerticesCompleted );
510  cell_complex.clear();
511  if ( ! ok ) return false;
512 
513  // Build complex, as many maximal cells as convex hull facets.
514  // convex hull facet -> cell of complex
515  // convex hull ridge -> face of complex
516  // (1) Get cell vertices, count ridges/faces and compute their vertices
517  std::map< Ridge, Index > r2f;
518  computeFacetAndRidgeVertices( del,
519  cell_complex.cell_vertices,
520  r2f,
521  cell_complex.true_face_vertices );
522  // (2) assign ridges/faces to cell and conversely
523  const Index nb_r = r2f.size();
524  cell_complex.true_face_cell .resize( nb_r, cell_complex.infiniteCell() );
525  cell_complex.false_face_cell.resize( nb_r, cell_complex.infiniteCell() );
526  cell_complex.true_face_vertices.resize( nb_r );
527  for ( Index cur_f = 0; cur_f < del.nbFiniteFacets(); ++cur_f ) {
528  const auto& facet = del.facets[ cur_f ];
529  FaceRange current_faces;
530  for ( auto neigh_f : facet.neighbors ) {
531  const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
532  const bool pos = cur_f < neigh_f;
533  const Index cur_r = r2f[ r ];
534  cell_complex.true_face_cell [ cur_r ] = r.first;
535  if ( r.second >= del.nbFiniteFacets() )
536  cell_complex.false_face_cell[ cur_r ] = cell_complex.infiniteCell();
537  else
538  cell_complex.false_face_cell[ cur_r ] = r.second;
539  current_faces.emplace_back( cur_r, pos );
540  }
541  cell_complex.cell_faces.push_back( current_faces );
542  }
543  // (3) Takes care of vertex positions
544  del.getVertexPositions( cell_complex.vertex_position );
545  return true;
546 }
547 
548 
549 //-----------------------------------------------------------------------------
550 template < int dim, typename TInteger, typename TInternalInteger >
551 template < typename QHull >
552 void
553 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeFacetAndRidgeVertices
554 ( const QHull& hull,
555  std::vector< IndexRange >& cell_vertices,
556  std::map< typename QHull::Ridge, Index >& r2f,
557  std::vector< IndexRange >& face_vertices )
558 {
559  typedef typename QHull::Ridge Ridge;
560 
561  ASSERT( hull.status() >= QHull::Status::VerticesCompleted
562  && hull.status() <= QHull::Status::AllCompleted );
563 
564  // Get cell vertices and sort them
565  bool ok_fv = hull.getFacetVertices( cell_vertices );
566  if ( ! ok_fv )
567  trace.error() << "[ConvexityHelper::computeFacetAndRidgeVertices]"
568  << " method hull.getFacetVertices failed."
569  << " Maybe QuickHull was not computed till VerticesCompleted."
570  << std::endl;
571  std::vector< IndexRange > sorted_cell_vertices = cell_vertices;
572  for ( auto& vtcs : sorted_cell_vertices )
573  std::sort( vtcs.begin(), vtcs.end() );
574  cell_vertices.resize( hull.nbFiniteFacets() );
575 
576  // Count ridges/faces and compute their vertices
577  Index nb_r = 0;
578  face_vertices.clear();
579  for ( Index cur_f = 0; cur_f < hull.nbFiniteFacets(); ++cur_f ) {
580  const auto& facet = hull.facets[ cur_f ];
581  for ( auto neigh_f : facet.neighbors ) {
582  const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
583  auto itr = r2f.find( r );
584  if ( itr == r2f.end() ) {
585  IndexRange result;
586  std::set_intersection( sorted_cell_vertices[ cur_f ].cbegin(),
587  sorted_cell_vertices[ cur_f ].cend(),
588  sorted_cell_vertices[ neigh_f ].cbegin(),
589  sorted_cell_vertices[ neigh_f ].cend(),
590  std::back_inserter( result ) );
591  face_vertices.push_back( result );
592  r2f[ r ] = nb_r++;
593  }
594  }
595  }
596 }
597 
598 // //
599 ///////////////////////////////////////////////////////////////////////////////
Generated on Tue Jun 1 2021 17:04:59 for DGtal by doxygen 1.9.1