DGtal  1.3.beta
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 #include "DGtal/kernel/IntegerConverter.h"
35 //////////////////////////////////////////////////////////////////////////////
36 
37 ///////////////////////////////////////////////////////////////////////////////
38 // IMPLEMENTATION of inline methods.
39 ///////////////////////////////////////////////////////////////////////////////
40 
41 ///////////////////////////////////////////////////////////////////////////////
42 // ----------------------- Standard services ------------------------------
43 
44 //-----------------------------------------------------------------------------
45 template < int dim, typename TInteger, typename TInternalInteger >
46 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
47 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeSimplex
48 ( const std::vector< Point >& input_points,
49  bool remove_duplicates )
50 {
51  std::vector< Point > X;
52  if ( remove_duplicates )
53  {
54  std::set<Point> S;
55  for ( auto&& p : input_points ) S.insert( p );
56  X = std::vector< Point >( S.cbegin(), S.cend() );
57  }
58  else X = input_points;
59  LatticePolytope P( X.cbegin(), X.cend() );
60  if ( P.nbHalfSpaces() != 0 )
61  return P;
62  else
63  return computeDegeneratedLatticePolytope( X );
64 }
65 
66 //-----------------------------------------------------------------------------
67 template < int dim, typename TInteger, typename TInternalInteger >
68 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
69 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
70 computeDegeneratedLatticePolytope
71 ( std::vector< Point >& input_points )
72 {
73  typedef typename LatticePolytope::Domain Domain;
74  typedef typename LatticePolytope::HalfSpace PolytopeHalfSpace;
75  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
76  typedef typename ConvexHull::IndexRange IndexRange;
77  typedef typename ConvexHull::HalfSpace ConvexHullHalfSpace;
78  // input_points is a range of points with no duplicates, but which
79  // seems to be not full dimensional.
80  if ( input_points.size() <= 1 )
81  return LatticePolytope( input_points.cbegin(), input_points.cend() );
82  // At least 1-dimensional
83  std::vector< Vector > basis;
84  std::vector< Integer > alpha;
85  basis.push_back( input_points[ 1 ] - input_points[ 0 ] );
86  const auto n0 = basis[ 0 ].norm();
87  alpha.push_back( Integer( 0 ) );
88  alpha.push_back( basis[ 0 ].dot( basis[ 0 ] ) );
89  Index i = 2;
90  while ( i < input_points.size() ) {
91  Vector v = input_points[ i ] - input_points[ 0 ];
92  alpha.push_back( basis[ 0 ].dot( v ) );
93  const auto ni = v.norm();
94  const double alignment =
95  fabs( fabs( NumberTraits< Integer >::castToDouble( alpha.back() ) )
96  - ( n0 * ni ) );
97  if ( alignment > 1e-8 ) break;
98  i++;
99  }
100  if ( i == input_points.size() )
101  { // 1-dimensional
102  Index a = 0, b = 0;
103  for ( i = 1; i < input_points.size(); i++ )
104  {
105  if ( alpha[ i ] < alpha[ a ] ) a = i;
106  if ( alpha[ i ] > alpha[ b ] ) b = i;
107  }
108  std::vector< Point > X( 2 );
109  X[ 0 ] = input_points[ a ];
110  X[ 1 ] = input_points[ b ];
111  return LatticePolytope( X.cbegin(), X.cend() );
112  }
113  // at least 2-dimensional
114  ASSERT( dimension > 1 );
115  if ( dimension == 2 )
116  {
117  std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
118  << " Weird error: found initial full dimensional simplex" << std::endl;
119  return LatticePolytope();
120  }
121  if ( dimension >= 4 )
122  {
123  std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
124  << "Degenerated lattice polytope in nD, n >= 4 is not implemented"
125  << std::endl;
126  return LatticePolytope();
127  }
128  basis.push_back( input_points[ i ] - input_points[ 0 ] );
129  Vector n = detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
130  ::crossProduct( basis[ 0 ], basis[ 1 ] );
131  if ( n == Vector::zero )
132  {
133  std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
134  << "Weird numerical error, u . v != |u| |v| but u x v != 0"
135  << std::endl;
136  return LatticePolytope();
137  }
138  // Now the set of input points should be full dimensional.
139  input_points.push_back( input_points[ 0 ] + n );
140  // Compute domain
141  Point l = input_points[ 0 ];
142  Point u = input_points[ 0 ];
143  for ( const auto& p : input_points ) {
144  l = l.inf( p );
145  u = u.sup( p );
146  }
147  Domain domain( l, u );
148  // Compute convex hull
149  ConvexHull hull;
150  hull.setInput( input_points, false );
151  const auto target = ConvexHull::Status::FacetsCompleted;
152  IndexRange full_splx = { 0, 1, i, input_points.size() - 1 };
153  bool ok_init = hull.setInitialSimplex( full_splx );
154  if ( ! ok_init )
155  {
156  std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
157  << "Weird error in hull.setInitialSimplex" << std::endl;
158  return LatticePolytope();
159  }
160  bool ok_hull = hull.computeConvexHull( target );
161  if ( ! ok_hull )
162  {
163  std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
164  << "Weird error in hull.computeConvexHull" << std::endl;
165  return LatticePolytope();
166  }
167  // Initialize polytope
168  std::vector< ConvexHullHalfSpace > HS;
169  std::vector< PolytopeHalfSpace > PHS;
170  hull.getFacetHalfSpaces( HS );
171  PHS.reserve( HS.size()+2 );
172  for ( auto& H : HS ) {
173  Vector N;
174  Integer nu;
175  for ( Dimension i = 0; i < dim; ++i )
176  N[ i ] = IntegerConverter< dimension, Integer >::cast( H.internalNormal()[ i ] );
177  nu = IntegerConverter< dimension, Integer >::cast( H.internalIntercept() );
178  PHS.emplace_back( N, nu );
179  }
180  // Add top constraint.
181  Integer nu0 = input_points[ 0 ].dot( n );
182  PHS.emplace_back( n, nu0 );
183  return LatticePolytope( domain, PHS.cbegin(), PHS.cend(),
184  false, false );
185 }
186 
187 //-----------------------------------------------------------------------------
188 template < int dim, typename TInteger, typename TInternalInteger >
189 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
190 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
191 computeLatticePolytope
192 ( const std::vector< Point >& input_points,
193  bool remove_duplicates,
194  bool make_minkowski_summable )
195 {
196  typedef typename LatticePolytope::Domain Domain;
197  typedef typename LatticePolytope::HalfSpace PolytopeHalfSpace;
198  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
199  typedef typename ConvexHull::HalfSpace ConvexHullHalfSpace;
200  typedef typename ConvexHull::Ridge Ridge;
201  if ( input_points.empty() ) return LatticePolytope();
202  if ( input_points.size() <= ( dimension + 1) )
203  return computeSimplex( input_points, remove_duplicates );
204  // Compute domain
205  Point l = input_points[ 0 ];
206  Point u = input_points[ 0 ];
207  for ( const auto& p : input_points ) {
208  l = l.inf( p );
209  u = u.sup( p );
210  }
211  Domain domain( l, u );
212  // Compute convex hull
213  ConvexHull hull;
214  hull.setInput( input_points, remove_duplicates );
215  const auto target = ( make_minkowski_summable && dimension == 3 )
216  ? ConvexHull::Status::VerticesCompleted
217  : ConvexHull::Status::FacetsCompleted;
218  bool ok = hull.computeConvexHull( target );
219  if ( ! ok ) // set of points is not full dimensional
220  return computeDegeneratedLatticePolytope( hull.points );
221  // Initialize polytope
222  std::vector< ConvexHullHalfSpace > HS;
223  std::vector< PolytopeHalfSpace > PHS;
224  hull.getFacetHalfSpaces( HS );
225  PHS.reserve( HS.size() );
226  for ( auto& H : HS ) {
227  Vector N;
228  Integer nu;
229  for ( Dimension i = 0; i < dim; ++i )
230  N[ i ] = IntegerConverter< dimension, Integer >::cast( H.internalNormal()[ i ] );
231  nu = IntegerConverter< dimension, Integer >::cast( H.internalIntercept() );
232  PHS.emplace_back( N, nu );
233  }
234  if ( make_minkowski_summable && dimension >= 4 )
235  trace.warning() << "[ConvexityHelper::computeLatticePolytope]"
236  << " Not implemented starting from dimension 4."
237  << std::endl;
238  if ( make_minkowski_summable && dimension == 3 )
239  {
240  // Compute ridge vertices to add edge constraints.
241  std::vector< Point > positions;
242  std::vector< IndexRange > facet_vertices;
243  std::vector< IndexRange > ridge_vertices;
244  std::map< Ridge, Index > ridge2index;
245  hull.getVertexPositions( positions );
246  computeFacetAndRidgeVertices( hull, facet_vertices,
247  ridge2index, ridge_vertices );
248  for ( auto p : ridge2index ) {
249  const auto r = p.first;
250  // Copy by value since PHS may be reallocated during the iteration.
251  const auto U = PHS[ r.first ].N; // normal of facet 1
252  const auto V = PHS[ r.second ].N; // normal of facet 2
253  const auto& S = ridge_vertices[ p.second ]; // vertices along facets 1, 2
254  ASSERT( S.size() == 2 && "Invalid ridge" );
255  const auto& P0 = positions[ S[ 0 ] ];
256  const auto& P1 = positions[ S[ 1 ] ];
257  auto E = P1 - P0; // edge 1, 2
258  const auto UxV =
259  detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
260  ::crossProduct( U, V ); // parallel to E
261  ASSERT( E.dot( UxV ) != 0 && "Invalid E / UxV" );
262  if ( E.dot( UxV ) <= 0 ) E = -E; // force correct orientation
263  const auto E1 =
264  detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
265  ::crossProduct( U, E ); // edge on facet 1
266  const auto E2 =
267  detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
268  ::crossProduct( E, V ); // edge on facet 2
269  ASSERT( E1.dot( U ) == 0 && "Invalid E1 / U" );
270  ASSERT( E1.dot( V ) < 0 && "Invalid E1 / V" );
271  ASSERT( E2.dot( V ) == 0 && "Invalid E1 / V" );
272  ASSERT( E2.dot( U ) < 0 && "Invalid E1 / U" );
273  for ( Dimension k = 0; k < dimension; ++k ) {
274  const auto W = U[ k ] * V - V[ k ] * U;
275  const auto nn1 = W.dot( E1 );
276  const auto nn2 = W.dot( E2 );
277  if ( nn1 > 0 && nn2 > 0 ) {
278  PHS.emplace_back( -W, -W.dot( P0 ) );
279  ASSERT( E1.dot(-W ) < 0 && "Invalid E1 /-W" );
280  ASSERT( E2.dot(-W ) < 0 && "Invalid E2 /-W" );
281  }
282  else if ( nn1 < 0 && nn2 < 0 ) {
283  PHS.emplace_back( W, W.dot( P0 ) );
284  ASSERT( E1.dot( W ) < 0 && "Invalid E1 / W" );
285  ASSERT( E2.dot( W ) < 0 && "Invalid E2 / W" );
286  }
287  }
288  }
289  }
290  return LatticePolytope( domain, PHS.cbegin(), PHS.cend(),
291  make_minkowski_summable && ( dimension <= 3 ), true );
292 }
293 
294 //-----------------------------------------------------------------------------
295 template < int dim, typename TInteger, typename TInternalInteger >
296 template < typename TSurfaceMesh >
297 bool
298 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
299 ( TSurfaceMesh& mesh,
300  const std::vector< Point >& input_points,
301  bool remove_duplicates )
302 {
303  typedef TSurfaceMesh SurfaceMesh;
304  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
305  typedef typename ConvexHull::IndexRange IndexRange;
306  ConvexHull hull;
307  hull.setInput( input_points, remove_duplicates );
308  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
309  if ( !ok ) return false;
310  std::vector< RealPoint > positions;
311  hull.getVertexPositions( positions );
312  std::vector< IndexRange > faces;
313  hull.getFacetVertices( faces );
314  mesh = SurfaceMesh( positions.cbegin(), positions.cend(),
315  faces.cbegin(), faces.cend() );
316  return true;
317 }
318 
319 //-----------------------------------------------------------------------------
320 template < int dim, typename TInteger, typename TInternalInteger >
321 bool
322 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
323 ( PolygonalSurface< Point >& polysurf,
324  const std::vector< Point >& input_points,
325  bool remove_duplicates )
326 {
327  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
328  typedef typename ConvexHull::IndexRange IndexRange;
329  ConvexHull hull;
330  hull.setInput( input_points, remove_duplicates );
331  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
332  if ( !ok ) return false;
333  std::vector< Point > positions;
334  hull.getVertexPositions( positions );
335  std::vector< IndexRange > faces;
336  hull.getFacetVertices( faces );
337  // build polygonal surface
338  polysurf.clear();
339  for ( auto p : positions ) polysurf.addVertex( p );
340  for ( auto f : faces ) polysurf.addPolygonalFace( f );
341  return polysurf.build();
342 }
343 
344 //-----------------------------------------------------------------------------
345 template < int dim, typename TInteger, typename TInternalInteger >
346 bool
347 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullCellComplex
348 ( ConvexCellComplex< Point >& cell_complex,
349  const std::vector< Point >& input_points,
350  bool remove_duplicates )
351 {
352  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
353  typedef typename ConvexHull::IndexRange IndexRange;
354  typedef typename ConvexCellComplex< Point >::FaceRange FaceRange;
355  ConvexHull hull;
356  hull.setInput( input_points, remove_duplicates );
357  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
358  cell_complex.clear();
359  if ( ! ok ) return false;
360  // Build complex, only 1 finite cell and as many faces as convex hull facets.
361  // Taking care of faces for each cell (here one cell borders all faces).
362  std::vector< IndexRange > faces;
363  hull.getFacetVertices( faces );
364  FaceRange all_faces;
365  for ( Index i = 0; i < faces.size(); i++ )
366  all_faces.push_back( { i, true } );
367  cell_complex.cell_faces.push_back( all_faces );
368  // Vertices of this unique cell will be computed lazily on request.
369  // Taking care of each face.
370  for ( const auto& vtcs : faces )
371  {
372  // every inner face borders cell 0
373  cell_complex.true_face_cell.push_back( 0 );
374  // every outer face borders the infinite cell
375  cell_complex.false_face_cell.push_back( cell_complex.infiniteCell() );
376  }
377  // Taking care of vertices (in consistent order) of every face
378  cell_complex.true_face_vertices.swap( faces );
379  // Taking care of vertex positions
380  hull.getVertexPositions( cell_complex.vertex_position );
381  return true;
382 }
383 
384 //-----------------------------------------------------------------------------
385 template < int dim, typename TInteger, typename TInternalInteger >
386 bool
387 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeDelaunayCellComplex
388 ( ConvexCellComplex< Point >& cell_complex,
389  const std::vector< Point >& input_points,
390  bool remove_duplicates )
391 {
392  typedef QuickHull< LatticeDelaunayKernel > Delaunay;
393  typedef typename Delaunay::Ridge Ridge;
394  typedef typename ConvexCellComplex< Point >::FaceRange FaceRange;
395 
396  Delaunay del;
397  del.setInput( input_points, remove_duplicates );
398  bool ok = del.computeConvexHull( Delaunay::Status::VerticesCompleted );
399  cell_complex.clear();
400  if ( ! ok ) return false;
401 
402  // Build complex, as many maximal cells as convex hull facets.
403  // convex hull facet -> cell of complex
404  // convex hull ridge -> face of complex
405  // (1) Get cell vertices, count ridges/faces and compute their vertices
406  std::map< Ridge, Index > r2f;
407  computeFacetAndRidgeVertices( del,
408  cell_complex.cell_vertices,
409  r2f,
410  cell_complex.true_face_vertices );
411  // (2) assign ridges/faces to cell and conversely
412  const Index nb_r = r2f.size();
413  cell_complex.true_face_cell .resize( nb_r, cell_complex.infiniteCell() );
414  cell_complex.false_face_cell.resize( nb_r, cell_complex.infiniteCell() );
415  cell_complex.true_face_vertices.resize( nb_r );
416  for ( Index cur_f = 0; cur_f < del.nbFiniteFacets(); ++cur_f ) {
417  const auto& facet = del.facets[ cur_f ];
418  FaceRange current_faces;
419  for ( auto neigh_f : facet.neighbors ) {
420  const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
421  const bool pos = cur_f < neigh_f;
422  const Index cur_r = r2f[ r ];
423  cell_complex.true_face_cell [ cur_r ] = r.first;
424  if ( r.second >= del.nbFiniteFacets() )
425  cell_complex.false_face_cell[ cur_r ] = cell_complex.infiniteCell();
426  else
427  cell_complex.false_face_cell[ cur_r ] = r.second;
428  current_faces.emplace_back( cur_r, pos );
429  }
430  cell_complex.cell_faces.push_back( current_faces );
431  }
432  // (3) Takes care of vertex positions
433  del.getVertexPositions( cell_complex.vertex_position );
434  return true;
435 }
436 
437 
438 //-----------------------------------------------------------------------------
439 template < int dim, typename TInteger, typename TInternalInteger >
440 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::RationalPolytope
441 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeRationalPolytope
442 ( const std::vector< RealPoint >& input_points,
443  Integer denominator,
444  bool remove_duplicates,
445  bool make_minkowski_summable )
446 {
447  if ( denominator < 1 )
448  trace.error() << "Invalid denominator " << denominator
449  << ". Should be greater or equal to 1." << std::endl;
450  typedef typename RationalPolytope::Domain Domain;
451  typedef typename RationalPolytope::HalfSpace PolytopeHalfSpace;
452  typedef QuickHull< RealConvexHullKernel > ConvexHull;
453  typedef typename ConvexHull::HalfSpace ConvexHullHalfSpace;
454  typedef typename ConvexHull::Ridge Ridge;
455  if ( input_points.empty() ) return RationalPolytope();
456  // Compute convex hull
457  ConvexHull hull( denominator );
458  hull.setInput( input_points, remove_duplicates );
459  const auto target = ( make_minkowski_summable && dimension == 3 )
460  ? ConvexHull::Status::VerticesCompleted
461  : ConvexHull::Status::FacetsCompleted;
462  bool ok = hull.computeConvexHull( target );
463  if ( ! ok ) return RationalPolytope();
464  // Compute domain (as a lattice domain)
465  auto l = hull.points[ 0 ];
466  auto u = hull.points[ 0 ];
467  for ( const auto& p : hull.points ) {
468  l = l.inf( p );
469  u = u.sup( p );
470  }
471  Domain domain( l, u );
472  trace.info() << "Domain l=" << l << " u=" << u << std::endl;
473  // Initialize polytope
474  std::vector< ConvexHullHalfSpace > HS;
475  std::vector< PolytopeHalfSpace > PHS;
476  hull.getFacetHalfSpaces( HS );
477  PHS.reserve( HS.size() );
478  for ( auto& H : HS ) {
479  Vector N;
480  Integer nu;
481  for ( Dimension i = 0; i < dim; ++i )
482  N[ i ] = (Integer) H.internalNormal()[ i ];
483  nu = (Integer) H.internalIntercept();
484  PHS.emplace_back( N, nu );
485  }
486  if ( make_minkowski_summable && dimension >= 4 )
487  trace.warning() << "[ConvexityHelper::computeRationalPolytope]"
488  << " Not implemented starting from dimension 4."
489  << std::endl;
490  if ( make_minkowski_summable && dimension == 3 )
491  {
492  // Compute ridge vertices to add edge constraints.
493  std::vector< Point > positions;
494  std::vector< IndexRange > facet_vertices;
495  std::vector< IndexRange > ridge_vertices;
496  std::map< Ridge, Index > ridge2index;
497  hull.getVertexPositions( positions );
498  computeFacetAndRidgeVertices( hull, facet_vertices,
499  ridge2index, ridge_vertices );
500  for ( auto p : ridge2index ) {
501  const auto r = p.first;
502  // Copy by value since PHS may be reallocated during the iteration.
503  const auto U = PHS[ r.first ].N; // normal of facet 1
504  const auto V = PHS[ r.second ].N; // normal of facet 2
505  const auto& S = ridge_vertices[ p.second ]; // vertices along facets 1, 2
506  ASSERT( S.size() == 2 && "Invalid ridge" );
507  const auto& P0 = positions[ S[ 0 ] ];
508  const auto& P1 = positions[ S[ 1 ] ];
509  auto E = P1 - P0; // edge 1, 2
510  const auto UxV =
511  detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
512  ::crossProduct( U, V ); // parallel to E
513  ASSERT( E.dot( UxV ) != 0 && "Invalid E / UxV" );
514  if ( E.dot( UxV ) <= 0 ) E = -E; // force correct orientation
515  const auto E1 =
516  detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
517  ::crossProduct( U, E ); // edge on facet 1
518  const auto E2 =
519  detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
520  ::crossProduct( E, V ); // edge on facet 2
521  ASSERT( E1.dot( U ) == 0 && "Invalid E1 / U" );
522  ASSERT( E1.dot( V ) < 0 && "Invalid E1 / V" );
523  ASSERT( E2.dot( V ) == 0 && "Invalid E1 / V" );
524  ASSERT( E2.dot( U ) < 0 && "Invalid E1 / U" );
525  for ( Dimension k = 0; k < dimension; ++k ) {
526  const auto W = U[ k ] * V - V[ k ] * U;
527  const auto nn1 = W.dot( E1 );
528  const auto nn2 = W.dot( E2 );
529  if ( nn1 > 0 && nn2 > 0 ) {
530  PHS.emplace_back( -W, -W.dot( P0 ) );
531  ASSERT( E1.dot(-W ) < 0 && "Invalid E1 /-W" );
532  ASSERT( E2.dot(-W ) < 0 && "Invalid E2 /-W" );
533  }
534  else if ( nn1 < 0 && nn2 < 0 ) {
535  PHS.emplace_back( W, W.dot( P0 ) );
536  ASSERT( E1.dot( W ) < 0 && "Invalid E1 / W" );
537  ASSERT( E2.dot( W ) < 0 && "Invalid E2 / W" );
538  }
539  }
540  }
541  }
542  return RationalPolytope( denominator, domain, PHS.cbegin(), PHS.cend(),
543  make_minkowski_summable && ( dimension <= 3 ), true );
544 }
545 
546 
547 //-----------------------------------------------------------------------------
548 template < int dim, typename TInteger, typename TInternalInteger >
549 template < typename TSurfaceMesh >
550 bool
551 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
552 ( TSurfaceMesh& mesh,
553  const std::vector< RealPoint >& input_points,
554  double precision,
555  bool remove_duplicates )
556 {
557  typedef TSurfaceMesh SurfaceMesh;
558  typedef QuickHull< RealConvexHullKernel > ConvexHull;
559  typedef typename ConvexHull::IndexRange IndexRange;
560  ConvexHull hull( precision );
561  hull.setInput( input_points, remove_duplicates );
562  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
563  if ( !ok ) return false;
564  std::vector< RealPoint > positions;
565  hull.getVertexPositions( positions );
566  std::vector< IndexRange > faces;
567  hull.getFacetVertices( faces );
568  mesh = SurfaceMesh( positions.cbegin(), positions.cend(),
569  faces.cbegin(), faces.cend() );
570  return true;
571 }
572 
573 
574 //-----------------------------------------------------------------------------
575 template < int dim, typename TInteger, typename TInternalInteger >
576 bool
577 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
578 ( PolygonalSurface< RealPoint >& polysurf,
579  const std::vector< RealPoint >& input_points,
580  double precision,
581  bool remove_duplicates )
582 {
583  typedef QuickHull< RealConvexHullKernel > ConvexHull;
584  typedef typename ConvexHull::IndexRange IndexRange;
585  ConvexHull hull( precision );
586  hull.setInput( input_points, remove_duplicates );
587  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
588  if ( !ok ) return false;
589  std::vector< RealPoint > positions;
590  hull.getVertexPositions( positions );
591  std::vector< IndexRange > faces;
592  hull.getFacetVertices( faces );
593  // build polygonal surface
594  polysurf.clear();
595  for ( auto p : positions ) polysurf.addVertex( p );
596  for ( auto f : faces ) polysurf.addPolygonalFace( f );
597  return polysurf.build();
598 }
599 
600 //-----------------------------------------------------------------------------
601 template < int dim, typename TInteger, typename TInternalInteger >
602 bool
603 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullCellComplex
604 ( ConvexCellComplex< RealPoint >& cell_complex,
605  const std::vector< RealPoint >& input_points,
606  double precision,
607  bool remove_duplicates )
608 {
609  typedef QuickHull< RealConvexHullKernel > ConvexHull;
610  typedef typename ConvexHull::IndexRange IndexRange;
611  typedef typename ConvexCellComplex< RealPoint >::FaceRange FaceRange;
612  ConvexHull hull( precision );
613  hull.setInput( input_points, remove_duplicates );
614  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
615  cell_complex.clear();
616  if ( ! ok ) return false;
617  // Build complex, only 1 finite cell and as many faces as convex hull facets.
618  // Taking care of faces for each cell (here one cell borders all faces).
619  std::vector< IndexRange > faces;
620  hull.getFacetVertices( faces );
621  FaceRange all_faces;
622  for ( Index i = 0; i < faces.size(); i++ )
623  all_faces.push_back( { i, true } );
624  cell_complex.cell_faces.push_back( all_faces );
625  // Vertices of this unique cell will be computed lazily on request.
626  // Taking care of each face.
627  for ( const auto& vtcs : faces )
628  {
629  // every inner face borders cell 0
630  cell_complex.true_face_cell.push_back( 0 );
631  // every outer face borders the infinite cell
632  cell_complex.false_face_cell.push_back( cell_complex.infiniteCell() );
633  }
634  // Taking care of vertices (in consistent order) of every face
635  cell_complex.true_face_vertices.swap( faces );
636  // Taking care of vertex positions
637  hull.getVertexPositions( cell_complex.vertex_position );
638  return true;
639 }
640 
641 
642 //-----------------------------------------------------------------------------
643 template < int dim, typename TInteger, typename TInternalInteger >
644 bool
645 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeDelaunayCellComplex
646 ( ConvexCellComplex< RealPoint >& cell_complex,
647  const std::vector< RealPoint >& input_points,
648  double precision,
649  bool remove_duplicates )
650 {
651  typedef QuickHull< RealDelaunayKernel > Delaunay;
652  typedef typename Delaunay::Ridge Ridge;
653  typedef typename ConvexCellComplex< RealPoint >::FaceRange FaceRange;
654 
655  Delaunay del( precision );
656  del.setInput( input_points, remove_duplicates );
657  bool ok = del.computeConvexHull( Delaunay::Status::VerticesCompleted );
658  cell_complex.clear();
659  if ( ! ok ) return false;
660 
661  // Build complex, as many maximal cells as convex hull facets.
662  // convex hull facet -> cell of complex
663  // convex hull ridge -> face of complex
664  // (1) Get cell vertices, count ridges/faces and compute their vertices
665  std::map< Ridge, Index > r2f;
666  computeFacetAndRidgeVertices( del,
667  cell_complex.cell_vertices,
668  r2f,
669  cell_complex.true_face_vertices );
670  // (2) assign ridges/faces to cell and conversely
671  const Index nb_r = r2f.size();
672  cell_complex.true_face_cell .resize( nb_r, cell_complex.infiniteCell() );
673  cell_complex.false_face_cell.resize( nb_r, cell_complex.infiniteCell() );
674  cell_complex.true_face_vertices.resize( nb_r );
675  for ( Index cur_f = 0; cur_f < del.nbFiniteFacets(); ++cur_f ) {
676  const auto& facet = del.facets[ cur_f ];
677  FaceRange current_faces;
678  for ( auto neigh_f : facet.neighbors ) {
679  const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
680  const bool pos = cur_f < neigh_f;
681  const Index cur_r = r2f[ r ];
682  cell_complex.true_face_cell [ cur_r ] = r.first;
683  if ( r.second >= del.nbFiniteFacets() )
684  cell_complex.false_face_cell[ cur_r ] = cell_complex.infiniteCell();
685  else
686  cell_complex.false_face_cell[ cur_r ] = r.second;
687  current_faces.emplace_back( cur_r, pos );
688  }
689  cell_complex.cell_faces.push_back( current_faces );
690  }
691  // (3) Takes care of vertex positions
692  del.getVertexPositions( cell_complex.vertex_position );
693  return true;
694 }
695 
696 
697 //-----------------------------------------------------------------------------
698 template < int dim, typename TInteger, typename TInternalInteger >
699 template < typename QHull >
700 void
701 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeFacetAndRidgeVertices
702 ( const QHull& hull,
703  std::vector< IndexRange >& cell_vertices,
704  std::map< typename QHull::Ridge, Index >& r2f,
705  std::vector< IndexRange >& face_vertices )
706 {
707  typedef typename QHull::Ridge Ridge;
708 
709  ASSERT( hull.status() >= QHull::Status::VerticesCompleted
710  && hull.status() <= QHull::Status::AllCompleted );
711 
712  // Get cell vertices and sort them
713  bool ok_fv = hull.getFacetVertices( cell_vertices );
714  if ( ! ok_fv )
715  trace.error() << "[ConvexityHelper::computeFacetAndRidgeVertices]"
716  << " method hull.getFacetVertices failed."
717  << " Maybe QuickHull was not computed till VerticesCompleted."
718  << std::endl;
719  std::vector< IndexRange > sorted_cell_vertices = cell_vertices;
720  for ( auto& vtcs : sorted_cell_vertices )
721  std::sort( vtcs.begin(), vtcs.end() );
722  cell_vertices.resize( hull.nbFiniteFacets() );
723 
724  // Count ridges/faces and compute their vertices
725  Index nb_r = 0;
726  face_vertices.clear();
727  for ( Index cur_f = 0; cur_f < hull.nbFiniteFacets(); ++cur_f ) {
728  const auto& facet = hull.facets[ cur_f ];
729  for ( auto neigh_f : facet.neighbors ) {
730  const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
731  auto itr = r2f.find( r );
732  if ( itr == r2f.end() ) {
733  IndexRange result;
734  std::set_intersection( sorted_cell_vertices[ cur_f ].cbegin(),
735  sorted_cell_vertices[ cur_f ].cend(),
736  sorted_cell_vertices[ neigh_f ].cbegin(),
737  sorted_cell_vertices[ neigh_f ].cend(),
738  std::back_inserter( result ) );
739  face_vertices.push_back( result );
740  r2f[ r ] = nb_r++;
741  }
742  }
743  }
744 }
745 
746 // //
747 ///////////////////////////////////////////////////////////////////////////////