DGtal  1.5.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 PointRange& input_points,
49  bool remove_duplicates )
50 {
51  PointRange X;
52  if ( remove_duplicates )
53  {
54  std::set<Point> S;
55  for ( auto&& p : input_points ) S.insert( p );
56  X = PointRange( 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 ( PointRange& input_points )
72 {
73  typedef typename LatticePolytope::Domain Domain;
74  typedef typename LatticePolytope::HalfSpace PolytopeHalfSpace;
75  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
76  typedef typename ConvexHull::HalfSpace ConvexHullHalfSpace;
77  // input_points is a range of points with no duplicates, but which
78  // seems to be not full dimensional.
79  if ( input_points.size() <= 1 )
80  return LatticePolytope( input_points.cbegin(), input_points.cend() );
81  // At least 1-dimensional
82  std::vector< Vector > basis;
83  std::vector< Integer > alpha;
84  basis.push_back( input_points[ 1 ] - input_points[ 0 ] );
85  const auto n0 = basis[ 0 ].norm();
86  alpha.push_back( Integer( 0 ) );
87  alpha.push_back( basis[ 0 ].dot( basis[ 0 ] ) );
88  Index i = 2;
89  while ( i < input_points.size() ) {
90  Vector v = input_points[ i ] - input_points[ 0 ];
91  alpha.push_back( basis[ 0 ].dot( v ) );
92  const auto ni = v.norm();
93  const double alignment =
94  fabs( fabs( NumberTraits< Integer >::castToDouble( alpha.back() ) )
95  - ( n0 * ni ) );
96  if ( alignment > 1e-8 ) break;
97  i++;
98  }
99  if ( i == input_points.size() )
100  { // 1-dimensional
101  Index a = 0, b = 0;
102  for ( i = 1; i < input_points.size(); i++ )
103  {
104  if ( alpha[ i ] < alpha[ a ] ) a = i;
105  if ( alpha[ i ] > alpha[ b ] ) b = i;
106  }
107  PointRange X( 2 );
108  X[ 0 ] = input_points[ a ];
109  X[ 1 ] = input_points[ b ];
110  return LatticePolytope( X.cbegin(), X.cend() );
111  }
112  // at least 2-dimensional
113  ASSERT( dimension > 1 );
114  if ( dimension == 2 )
115  {
116  std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
117  << " Weird error: found initial full dimensional simplex" << std::endl;
118  return LatticePolytope();
119  }
120  if ( dimension >= 4 )
121  {
122  std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
123  << "Degenerated lattice polytope in nD, n >= 4 is not implemented"
124  << std::endl;
125  return LatticePolytope();
126  }
127  basis.push_back( input_points[ i ] - input_points[ 0 ] );
128  Vector n = detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
129  ::crossProduct( basis[ 0 ], basis[ 1 ] );
130  if ( n == Vector::zero )
131  {
132  std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
133  << "Weird numerical error, u . v != |u| |v| but u x v != 0"
134  << std::endl;
135  return LatticePolytope();
136  }
137  // Now the set of input points should be full dimensional.
138  input_points.push_back( input_points[ 0 ] + n );
139  // Compute domain
140  Point l = input_points[ 0 ];
141  Point u = input_points[ 0 ];
142  for ( const auto& p : input_points ) {
143  l = l.inf( p );
144  u = u.sup( p );
145  }
146  Domain domain( l, u );
147  // Compute convex hull
148  ConvexHull hull;
149  hull.setInput( input_points, false );
150  const auto target = ConvexHull::Status::FacetsCompleted;
151  IndexRange full_splx = { 0, 1, i, input_points.size() - 1 };
152  bool ok_init = hull.setInitialSimplex( full_splx );
153  if ( ! ok_init )
154  {
155  std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
156  << "Weird error in hull.setInitialSimplex" << std::endl;
157  return LatticePolytope();
158  }
159  bool ok_hull = hull.computeConvexHull( target );
160  if ( ! ok_hull )
161  {
162  std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
163  << "Weird error in hull.computeConvexHull" << std::endl;
164  return LatticePolytope();
165  }
166  // Initialize polytope
167  std::vector< ConvexHullHalfSpace > HS;
168  std::vector< PolytopeHalfSpace > PHS;
169  hull.getFacetHalfSpaces( HS );
170  PHS.reserve( HS.size()+2 );
171  for ( auto& H : HS ) {
172  Vector N;
173  Integer nu;
174  for ( Dimension ii = 0; ii < dim; ++ii )
175  N[ ii ] = IntegerConverter< dimension, Integer >::cast( H.internalNormal()[ ii ] );
176  nu = IntegerConverter< dimension, Integer >::cast( H.internalIntercept() );
177  PHS.emplace_back( N, nu );
178  }
179  // Add top constraint.
180  Integer nu0 = input_points[ 0 ].dot( n );
181  PHS.emplace_back( n, nu0 );
182  return LatticePolytope( domain, PHS.cbegin(), PHS.cend(),
183  false, false );
184 }
185 
186 //-----------------------------------------------------------------------------
187 template < int dim, typename TInteger, typename TInternalInteger >
188 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
189 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::compute3DTriangle
190 ( const Point& a, const Point& b, const Point& c,
191  bool make_minkowski_summable )
192 {
193  if constexpr( dim != 3 ) return LatticePolytope();
194  using Op = detail::BoundedRationalPolytopeSpecializer< dimension, Integer>;
195  typedef typename LatticePolytope::Domain Domain;
196  typedef typename LatticePolytope::HalfSpace PolytopeHalfSpace;
197  // Compute domain
198  const Vector ab = b - a;
199  const Vector bc = c - b;
200  const Vector ca = a - c;
201  const Vector n = Op::crossProduct( ab, bc );
202  if ( n == Vector::zero )
203  return computeDegeneratedTriangle( a, b, c );
204  const Point low = a.inf( b ).inf( c );
205  const Point high = a.sup( b ).sup( c );
206  // Initialize polytope
207  std::vector< PolytopeHalfSpace > PHS;
208  PHS.reserve( make_minkowski_summable ? 11 : 5 );
209  const Integer n_a = n.dot( a );
210  const Vector u = Op::crossProduct( ab, n );
211  const Vector v = Op::crossProduct( bc, n );
212  const Vector w = Op::crossProduct( ca, n );
213  PHS.emplace_back( n, n_a );
214  PHS.emplace_back( -n, -n_a );
215  if ( ! make_minkowski_summable )
216  { // It is enough to have one constraint per edge.
217  PHS.emplace_back( u, u.dot( a ) );
218  PHS.emplace_back( v, v.dot( b ) );
219  PHS.emplace_back( w, w.dot( c ) );
220  }
221  else
222  { // Compute additionnal constraints on edges so that the
223  // Minkowski sum with axis-aligned edges is valid.
224  for ( Integer d = -1; d <= 1; d += 2 )
225  for ( Dimension k = 0; k < dim; k++ )
226  {
227  const Vector i = Vector::base( k, d );
228  const Vector eab = Op::crossProduct( ab, i );
229  const Integer eab_a = eab.dot( a );
230  if ( eab.dot( c ) < eab_a ) // c must be below plane (a,eab)
231  PHS.emplace_back( eab, eab_a );
232  const Vector ebc = Op::crossProduct( bc, i );
233  const Integer ebc_b = ebc.dot( b );
234  if ( ebc.dot( a ) < ebc_b ) // a must be below plane (b,ebc)
235  PHS.emplace_back( ebc, ebc_b );
236  const Vector eca = Op::crossProduct( ca, i );
237  const Integer eca_c = eca.dot( c );
238  if ( eca.dot( b ) < eca_c ) // b must be below plane (c,eca)
239  PHS.emplace_back( eca, eca_c );
240  }
241  }
242  return LatticePolytope( Domain( low, high ),
243  PHS.cbegin(), PHS.cend(),
244  make_minkowski_summable, false );
245 }
246 
247 //-----------------------------------------------------------------------------
248 template < int dim, typename TInteger, typename TInternalInteger >
249 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
250 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
251 computeDegeneratedTriangle
252 ( const Point& a, const Point& b, const Point& c )
253 {
254  if ( a == b ) return computeSegment( a, c );
255  if ( ( a == c ) || ( b == c ) ) return computeSegment( a, b );
256  // The three points are distinct, hence aligned. One is in-between the two others.
257  const Point low = a.inf( b ).inf( c );
258  const Point high = a.sup( b ).sup( c );
259  for ( Dimension k = 0; k < dim; k++ )
260  {
261  const auto lk = low [ k ];
262  const auto hk = high[ k ];
263  if ( ( a[ k ] != lk ) && ( a[ k ] != hk ) ) return computeSegment( b, c );
264  if ( ( b[ k ] != lk ) && ( b[ k ] != hk ) ) return computeSegment( a, c );
265  if ( ( c[ k ] != lk ) && ( c[ k ] != hk ) ) return computeSegment( a, b );
266  }
267  trace.error() << "[ConvexityHelper::computeSegmentFromDegeneratedTriangle] "
268  << "Should never arrive here." << std::endl;
269  return computeSegment( a, a );
270 }
271 
272 //-----------------------------------------------------------------------------
273 template < int dim, typename TInteger, typename TInternalInteger >
274 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
275 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeSegment
276 ( const Point& a, const Point& b )
277 {
278  if constexpr( dim != 3 ) return LatticePolytope();
279  using Op = detail::BoundedRationalPolytopeSpecializer< dimension, Integer>;
280  typedef typename LatticePolytope::Domain Domain;
281  typedef typename LatticePolytope::HalfSpace PolytopeHalfSpace;
282  // Compute domain
283  const Point low = a.inf( b );
284  const Point high = a.sup( b );
285  const Vector ab = b - a;
286  bool degenerate = ( ab == Vector::zero );
287  // Initialize polytope
288  std::vector< PolytopeHalfSpace > PHS;
289  if ( degenerate )
290  return LatticePolytope( Domain( low, high ), PHS.cbegin(), PHS.cend(), true, false );
291  PHS.reserve( 2*dim );
292  // Compute additionnal constraints on edges so that the
293  // Minkowski sum with axis-aligned edges is valid.
294  for ( Integer d = -1; d <= 1; d += 2 )
295  for ( Dimension k = 0; k < dim; k++ )
296  {
297  const Vector i = Vector::base( k, d );
298  const Vector e = Op::crossProduct( ab, i );
299  if ( e != Vector::zero )
300  {
301  const Integer e_a = e.dot( a );
302  PHS.emplace_back( e, e_a );
303  }
304  }
305  return LatticePolytope( Domain( low, high ), PHS.cbegin(), PHS.cend(), true, false );
306 }
307 
308 
309 //-----------------------------------------------------------------------------
310 template < int dim, typename TInteger, typename TInternalInteger >
311 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::PointRange
312 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
313 computeDegeneratedConvexHullVertices
314 ( PointRange& input_points )
315 {
316  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
317  // input_points is a range of points with no duplicates, but which
318  // seems to be not full dimensional.
319  if ( input_points.size() <= 1 )
320  return input_points;
321  // At least 1-dimensional
322  std::vector< Vector > basis;
323  std::vector< Integer > alpha;
324  basis.push_back( input_points[ 1 ] - input_points[ 0 ] );
325  const auto n0 = basis[ 0 ].norm();
326  alpha.push_back( Integer( 0 ) );
327  alpha.push_back( basis[ 0 ].dot( basis[ 0 ] ) );
328  Index i = 2;
329  while ( i < input_points.size() ) {
330  Vector v = input_points[ i ] - input_points[ 0 ];
331  alpha.push_back( basis[ 0 ].dot( v ) );
332  const auto ni = v.norm();
333  const double alignment =
334  fabs( fabs( NumberTraits< Integer >::castToDouble( alpha.back() ) )
335  - ( n0 * ni ) );
336  if ( alignment > 1e-8 ) break;
337  i++;
338  }
339  if ( i == input_points.size() )
340  { // 1-dimensional
341  Index a = 0, b = 0;
342  for ( i = 1; i < input_points.size(); i++ )
343  {
344  if ( alpha[ i ] < alpha[ a ] ) a = i;
345  if ( alpha[ i ] > alpha[ b ] ) b = i;
346  }
347  PointRange X( 2 );
348  X[ 0 ] = input_points[ a ];
349  X[ 1 ] = input_points[ b ];
350  return X;
351  }
352  // at least 2-dimensional
353  ASSERT( dimension > 1 );
354  if ( dimension == 2 )
355  {
356  std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
357  << " Weird error: found initial full dimensional simplex" << std::endl;
358  return PointRange();
359  }
360  if ( dimension >= 4 )
361  {
362  std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
363  << "Degenerated lattice polytope in nD, n >= 4 is not implemented"
364  << std::endl;
365  return PointRange();
366  }
367  basis.push_back( input_points[ i ] - input_points[ 0 ] );
368  Vector n = detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
369  ::crossProduct( basis[ 0 ], basis[ 1 ] );
370  if ( n == Vector::zero )
371  {
372  std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
373  << "Weird numerical error, u . v != |u| |v| but u x v != 0"
374  << std::endl;
375  return PointRange();
376  }
377  // Now the set of input points should be full dimensional.
378  input_points.push_back( input_points[ 0 ] + n );
379  // Compute convex hull
380  ConvexHull hull;
381  hull.setInput( input_points, false );
382  const auto target = ConvexHull::Status::VerticesCompleted;
383  IndexRange full_splx = { 0, 1, i, input_points.size() - 1 };
384  bool ok_init = hull.setInitialSimplex( full_splx );
385  if ( ! ok_init )
386  {
387  std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
388  << "Weird error in hull.setInitialSimplex" << std::endl;
389  return PointRange();
390  }
391  bool ok_hull = hull.computeConvexHull( target );
392  if ( ! ok_hull )
393  {
394  std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
395  << "Weird error in hull.computeConvexHull" << std::endl;
396  return PointRange();
397  }
398  // Get convex hull vertices and remove top point
399  PointRange X;
400  hull.getVertexPositions( X );
401  const std::size_t nX = X.size();
402  for ( std::size_t j = 0; j < nX; j++ )
403  if ( X[ j ] == input_points.back() )
404  {
405  X[ j ] = X.back();
406  break;
407  }
408  X.resize( nX - 1 );
409  return X;
410 }
411 
412 //-----------------------------------------------------------------------------
413 template < int dim, typename TInteger, typename TInternalInteger >
414 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
415 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
416 computeLatticePolytope
417 ( const PointRange& input_points,
418  bool remove_duplicates,
419  bool make_minkowski_summable )
420 {
421  typedef typename LatticePolytope::Domain Domain;
422  typedef typename LatticePolytope::HalfSpace PolytopeHalfSpace;
423  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
424  typedef typename ConvexHull::HalfSpace ConvexHullHalfSpace;
425  typedef typename ConvexHull::Ridge Ridge;
426  if ( input_points.empty() ) return LatticePolytope();
427  if ( input_points.size() <= ( dimension + 1) )
428  return computeSimplex( input_points, remove_duplicates );
429  // Compute domain
430  Point l = input_points[ 0 ];
431  Point u = input_points[ 0 ];
432  for ( std::size_t i = 1; i < input_points.size(); i++ )
433  {
434  const auto& p = input_points[ i ];
435  l = l.inf( p );
436  u = u.sup( p );
437  }
438  Domain domain( l, u );
439  // Compute convex hull
440  ConvexHull hull;
441  hull.setInput( input_points, remove_duplicates );
442  const auto target = ( make_minkowski_summable && dimension == 3 )
443  ? ConvexHull::Status::VerticesCompleted
444  : ConvexHull::Status::FacetsCompleted;
445  bool ok = hull.computeConvexHull( target );
446  if ( ! ok ) // set of points is not full dimensional
447  return computeDegeneratedLatticePolytope( hull.points );
448  // Initialize polytope
449  std::vector< ConvexHullHalfSpace > HS;
450  std::vector< PolytopeHalfSpace > PHS;
451  hull.getFacetHalfSpaces( HS );
452  PHS.reserve( HS.size() );
453  for ( auto& H : HS ) {
454  Vector N;
455  Integer nu;
456  for ( Dimension i = 0; i < dim; ++i )
457  N[ i ] = IntegerConverter< dimension, Integer >::cast( H.internalNormal()[ i ] );
458  nu = IntegerConverter< dimension, Integer >::cast( H.internalIntercept() );
459  PHS.emplace_back( N, nu );
460  }
461  if ( make_minkowski_summable && dimension >= 4 )
462  trace.warning() << "[ConvexityHelper::computeLatticePolytope]"
463  << " Not implemented starting from dimension 4."
464  << std::endl;
465  if ( make_minkowski_summable && dimension == 3 )
466  {
467  // Compute ridge vertices to add edge constraints.
468  PointRange positions;
469  std::vector< IndexRange > facet_vertices;
470  std::vector< IndexRange > ridge_vertices;
471  std::map< Ridge, Index > ridge2index;
472  hull.getVertexPositions( positions );
473  computeFacetAndRidgeVertices( hull, facet_vertices,
474  ridge2index, ridge_vertices );
475  for ( auto p : ridge2index ) {
476  const auto r = p.first;
477  // Copy by value since PHS may be reallocated during the iteration.
478  const auto U = PHS[ r.first ].N; // normal of facet 1
479  const auto V = PHS[ r.second ].N; // normal of facet 2
480  const auto& S = ridge_vertices[ p.second ]; // vertices along facets 1, 2
481  ASSERT( S.size() == 2 && "Invalid ridge" );
482  const auto& P0 = positions[ S[ 0 ] ];
483  const auto& P1 = positions[ S[ 1 ] ];
484  auto E = P1 - P0; // edge 1, 2
485  const auto UxV =
486  detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
487  ::crossProduct( U, V ); // parallel to E
488  ASSERT( E.dot( UxV ) != 0 && "Invalid E / UxV" );
489  if ( E.dot( UxV ) <= 0 ) E = -E; // force correct orientation
490  const auto E1 =
491  detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
492  ::crossProduct( U, E ); // edge on facet 1
493  const auto E2 =
494  detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
495  ::crossProduct( E, V ); // edge on facet 2
496  ASSERT( E1.dot( U ) == 0 && "Invalid E1 / U" );
497  ASSERT( E1.dot( V ) < 0 && "Invalid E1 / V" );
498  ASSERT( E2.dot( V ) == 0 && "Invalid E1 / V" );
499  ASSERT( E2.dot( U ) < 0 && "Invalid E1 / U" );
500  for ( Dimension k = 0; k < dimension; ++k ) {
501  const auto W = U[ k ] * V - V[ k ] * U;
502  const auto nn1 = W.dot( E1 );
503  const auto nn2 = W.dot( E2 );
504  if ( nn1 > 0 && nn2 > 0 ) {
505  PHS.emplace_back( -W, -W.dot( P0 ) );
506  ASSERT( E1.dot(-W ) < 0 && "Invalid E1 /-W" );
507  ASSERT( E2.dot(-W ) < 0 && "Invalid E2 /-W" );
508  }
509  else if ( nn1 < 0 && nn2 < 0 ) {
510  PHS.emplace_back( W, W.dot( P0 ) );
511  ASSERT( E1.dot( W ) < 0 && "Invalid E1 / W" );
512  ASSERT( E2.dot( W ) < 0 && "Invalid E2 / W" );
513  }
514  }
515  }
516  }
517  return LatticePolytope( domain, PHS.cbegin(), PHS.cend(),
518  make_minkowski_summable && ( dimension <= 3 ), true );
519 }
520 
521 //-----------------------------------------------------------------------------
522 template < int dim, typename TInteger, typename TInternalInteger >
523 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::PointRange
524 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
525 computeConvexHullVertices
526 ( const PointRange& input_points,
527  bool remove_duplicates )
528 {
529  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
530  PointRange positions;
531  ConvexHull hull;
532  hull.setInput( input_points, remove_duplicates );
533  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
534  if ( !ok )
535  {
536  PointRange Z( input_points );
537  return computeDegeneratedConvexHullVertices( Z );
538  }
539  hull.getVertexPositions( positions );
540  return positions;
541 }
542 
543 //-----------------------------------------------------------------------------
544 template < int dim, typename TInteger, typename TInternalInteger >
545 template < typename TSurfaceMesh >
546 bool
547 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
548 ( TSurfaceMesh& mesh,
549  const PointRange& input_points,
550  bool remove_duplicates )
551 {
552  typedef TSurfaceMesh SurfaceMesh;
553  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
554  ConvexHull hull;
555  hull.setInput( input_points, remove_duplicates );
556  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
557  if ( !ok ) return false;
558  std::vector< RealPoint > positions;
559  hull.getVertexPositions( positions );
560  std::vector< IndexRange > faces;
561  hull.getFacetVertices( faces );
562  mesh = SurfaceMesh( positions.cbegin(), positions.cend(),
563  faces.cbegin(), faces.cend() );
564  return true;
565 }
566 
567 //-----------------------------------------------------------------------------
568 template < int dim, typename TInteger, typename TInternalInteger >
569 bool
570 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
571 ( PolygonalSurface< Point >& polysurf,
572  const PointRange& input_points,
573  bool remove_duplicates )
574 {
575  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
576  ConvexHull hull;
577  hull.setInput( input_points, remove_duplicates );
578  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
579  if ( !ok ) return false;
580  PointRange positions;
581  hull.getVertexPositions( positions );
582  std::vector< IndexRange > faces;
583  hull.getFacetVertices( faces );
584  // build polygonal surface
585  polysurf.clear();
586  for ( auto p : positions ) polysurf.addVertex( p );
587  for ( auto f : faces ) polysurf.addPolygonalFace( f );
588  return polysurf.build();
589 }
590 
591 //-----------------------------------------------------------------------------
592 template < int dim, typename TInteger, typename TInternalInteger >
593 bool
594 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullCellComplex
595 ( ConvexCellComplex< Point >& cell_complex,
596  const PointRange& input_points,
597  bool remove_duplicates )
598 {
599  typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
600  typedef typename ConvexCellComplex< Point >::FaceRange FaceRange;
601  ConvexHull hull;
602  hull.setInput( input_points, remove_duplicates );
603  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
604  cell_complex.clear();
605  if ( ! ok ) return false;
606  // Build complex, only 1 finite cell and as many faces as convex hull facets.
607  // Taking care of faces for each cell (here one cell borders all faces).
608  std::vector< IndexRange > faces;
609  hull.getFacetVertices( faces );
610  FaceRange all_faces;
611  for ( Index i = 0; i < faces.size(); i++ )
612  all_faces.push_back( { i, true } );
613  cell_complex.cell_faces.push_back( all_faces );
614  // Vertices of this unique cell will be computed lazily on request.
615  // Taking care of each face.
616  for ( Index i = 0; i < faces.size(); i++ )
617  {
618  // every inner face borders cell 0
619  cell_complex.true_face_cell.push_back( 0 );
620  // every outer face borders the infinite cell
621  cell_complex.false_face_cell.push_back( cell_complex.infiniteCell() );
622  }
623  // Taking care of vertices (in consistent order) of every face
624  cell_complex.true_face_vertices.swap( faces );
625  // Taking care of vertex positions
626  hull.getVertexPositions( cell_complex.vertex_position );
627  return true;
628 }
629 
630 //-----------------------------------------------------------------------------
631 template < int dim, typename TInteger, typename TInternalInteger >
632 bool
633 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeDelaunayCellComplex
634 ( ConvexCellComplex< Point >& cell_complex,
635  const PointRange& input_points,
636  bool remove_duplicates )
637 {
638  typedef QuickHull< LatticeDelaunayKernel > Delaunay;
639  typedef typename Delaunay::Ridge Ridge;
640  typedef typename ConvexCellComplex< Point >::FaceRange FaceRange;
641 
642  Delaunay del;
643  del.setInput( input_points, remove_duplicates );
644  bool ok = del.computeConvexHull( Delaunay::Status::VerticesCompleted );
645  cell_complex.clear();
646  if ( ! ok ) return false;
647 
648  // Build complex, as many maximal cells as convex hull facets.
649  // convex hull facet -> cell of complex
650  // convex hull ridge -> face of complex
651  // (1) Get cell vertices, count ridges/faces and compute their vertices
652  std::map< Ridge, Index > r2f;
653  computeFacetAndRidgeVertices( del,
654  cell_complex.cell_vertices,
655  r2f,
656  cell_complex.true_face_vertices );
657  // (2) assign ridges/faces to cell and conversely
658  const Index nb_r = r2f.size();
659  cell_complex.true_face_cell .resize( nb_r, cell_complex.infiniteCell() );
660  cell_complex.false_face_cell.resize( nb_r, cell_complex.infiniteCell() );
661  cell_complex.true_face_vertices.resize( nb_r );
662  for ( Index cur_f = 0; cur_f < del.nbFiniteFacets(); ++cur_f ) {
663  const auto& facet = del.facets[ cur_f ];
664  FaceRange current_faces;
665  for ( auto neigh_f : facet.neighbors ) {
666  const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
667  const bool pos = cur_f < neigh_f;
668  const Index cur_r = r2f[ r ];
669  cell_complex.true_face_cell [ cur_r ] = r.first;
670  if ( r.second >= del.nbFiniteFacets() )
671  cell_complex.false_face_cell[ cur_r ] = cell_complex.infiniteCell();
672  else
673  cell_complex.false_face_cell[ cur_r ] = r.second;
674  current_faces.emplace_back( cur_r, pos );
675  }
676  cell_complex.cell_faces.push_back( current_faces );
677  }
678  // (3) Takes care of vertex positions
679  del.getVertexPositions( cell_complex.vertex_position );
680  return true;
681 }
682 
683 
684 //-----------------------------------------------------------------------------
685 template < int dim, typename TInteger, typename TInternalInteger >
686 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::RationalPolytope
687 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeRationalPolytope
688 ( const std::vector< RealPoint >& input_points,
689  Integer denominator,
690  bool remove_duplicates,
691  bool make_minkowski_summable )
692 {
693  if ( denominator < 1 )
694  trace.error() << "Invalid denominator " << denominator
695  << ". Should be greater or equal to 1." << std::endl;
696  typedef typename RationalPolytope::Domain Domain;
697  typedef typename RationalPolytope::HalfSpace PolytopeHalfSpace;
698  typedef QuickHull< RealConvexHullKernel > ConvexHull;
699  typedef typename ConvexHull::HalfSpace ConvexHullHalfSpace;
700  typedef typename ConvexHull::Ridge Ridge;
701  if ( input_points.empty() ) return RationalPolytope();
702  // Compute convex hull
703  ConvexHull hull( denominator );
704  hull.setInput( input_points, remove_duplicates );
705  const auto target = ( make_minkowski_summable && dimension == 3 )
706  ? ConvexHull::Status::VerticesCompleted
707  : ConvexHull::Status::FacetsCompleted;
708  bool ok = hull.computeConvexHull( target );
709  if ( ! ok ) return RationalPolytope();
710  // Compute domain (as a lattice domain)
711  auto l = hull.points[ 0 ];
712  auto u = hull.points[ 0 ];
713  for ( const auto& p : hull.points ) {
714  l = l.inf( p );
715  u = u.sup( p );
716  }
717  Domain domain( l, u );
718  trace.info() << "Domain l=" << l << " u=" << u << std::endl;
719  // Initialize polytope
720  std::vector< ConvexHullHalfSpace > HS;
721  std::vector< PolytopeHalfSpace > PHS;
722  hull.getFacetHalfSpaces( HS );
723  PHS.reserve( HS.size() );
724  for ( auto& H : HS ) {
725  Vector N;
726  Integer nu;
727  for ( Dimension i = 0; i < dim; ++i )
728  N[ i ] = (Integer) H.internalNormal()[ i ];
729  nu = (Integer) H.internalIntercept();
730  PHS.emplace_back( N, nu );
731  }
732  if ( make_minkowski_summable && dimension >= 4 )
733  trace.warning() << "[ConvexityHelper::computeRationalPolytope]"
734  << " Not implemented starting from dimension 4."
735  << std::endl;
736  if ( make_minkowski_summable && dimension == 3 )
737  {
738  // Compute ridge vertices to add edge constraints.
739  PointRange positions;
740  std::vector< IndexRange > facet_vertices;
741  std::vector< IndexRange > ridge_vertices;
742  std::map< Ridge, Index > ridge2index;
743  hull.getVertexPositions( positions );
744  computeFacetAndRidgeVertices( hull, facet_vertices,
745  ridge2index, ridge_vertices );
746  for ( auto p : ridge2index ) {
747  const auto r = p.first;
748  // Copy by value since PHS may be reallocated during the iteration.
749  const auto U = PHS[ r.first ].N; // normal of facet 1
750  const auto V = PHS[ r.second ].N; // normal of facet 2
751  const auto& S = ridge_vertices[ p.second ]; // vertices along facets 1, 2
752  ASSERT( S.size() == 2 && "Invalid ridge" );
753  const auto& P0 = positions[ S[ 0 ] ];
754  const auto& P1 = positions[ S[ 1 ] ];
755  auto E = P1 - P0; // edge 1, 2
756  const auto UxV =
757  detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
758  ::crossProduct( U, V ); // parallel to E
759  ASSERT( E.dot( UxV ) != 0 && "Invalid E / UxV" );
760  if ( E.dot( UxV ) <= 0 ) E = -E; // force correct orientation
761  const auto E1 =
762  detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
763  ::crossProduct( U, E ); // edge on facet 1
764  const auto E2 =
765  detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
766  ::crossProduct( E, V ); // edge on facet 2
767  ASSERT( E1.dot( U ) == 0 && "Invalid E1 / U" );
768  ASSERT( E1.dot( V ) < 0 && "Invalid E1 / V" );
769  ASSERT( E2.dot( V ) == 0 && "Invalid E1 / V" );
770  ASSERT( E2.dot( U ) < 0 && "Invalid E1 / U" );
771  for ( Dimension k = 0; k < dimension; ++k ) {
772  const auto W = U[ k ] * V - V[ k ] * U;
773  const auto nn1 = W.dot( E1 );
774  const auto nn2 = W.dot( E2 );
775  if ( nn1 > 0 && nn2 > 0 ) {
776  PHS.emplace_back( -W, -W.dot( P0 ) );
777  ASSERT( E1.dot(-W ) < 0 && "Invalid E1 /-W" );
778  ASSERT( E2.dot(-W ) < 0 && "Invalid E2 /-W" );
779  }
780  else if ( nn1 < 0 && nn2 < 0 ) {
781  PHS.emplace_back( W, W.dot( P0 ) );
782  ASSERT( E1.dot( W ) < 0 && "Invalid E1 / W" );
783  ASSERT( E2.dot( W ) < 0 && "Invalid E2 / W" );
784  }
785  }
786  }
787  }
788  return RationalPolytope( denominator, domain, PHS.cbegin(), PHS.cend(),
789  make_minkowski_summable && ( dimension <= 3 ), true );
790 }
791 
792 
793 //-----------------------------------------------------------------------------
794 template < int dim, typename TInteger, typename TInternalInteger >
795 template < typename TSurfaceMesh >
796 bool
797 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
798 ( TSurfaceMesh& mesh,
799  const std::vector< RealPoint >& input_points,
800  double precision,
801  bool remove_duplicates )
802 {
803  typedef TSurfaceMesh SurfaceMesh;
804  typedef QuickHull< RealConvexHullKernel > ConvexHull;
805  ConvexHull hull( precision );
806  hull.setInput( input_points, remove_duplicates );
807  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
808  if ( !ok ) return false;
809  std::vector< RealPoint > positions;
810  hull.getVertexPositions( positions );
811  std::vector< IndexRange > faces;
812  hull.getFacetVertices( faces );
813  mesh = SurfaceMesh( positions.cbegin(), positions.cend(),
814  faces.cbegin(), faces.cend() );
815  return true;
816 }
817 
818 
819 //-----------------------------------------------------------------------------
820 template < int dim, typename TInteger, typename TInternalInteger >
821 bool
822 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
823 ( PolygonalSurface< RealPoint >& polysurf,
824  const std::vector< RealPoint >& input_points,
825  double precision,
826  bool remove_duplicates )
827 {
828  typedef QuickHull< RealConvexHullKernel > ConvexHull;
829  ConvexHull hull( precision );
830  hull.setInput( input_points, remove_duplicates );
831  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
832  if ( !ok ) return false;
833  std::vector< RealPoint > positions;
834  hull.getVertexPositions( positions );
835  std::vector< IndexRange > faces;
836  hull.getFacetVertices( faces );
837  // build polygonal surface
838  polysurf.clear();
839  for ( auto p : positions ) polysurf.addVertex( p );
840  for ( auto f : faces ) polysurf.addPolygonalFace( f );
841  return polysurf.build();
842 }
843 
844 //-----------------------------------------------------------------------------
845 template < int dim, typename TInteger, typename TInternalInteger >
846 bool
847 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullCellComplex
848 ( ConvexCellComplex< RealPoint >& cell_complex,
849  const std::vector< RealPoint >& input_points,
850  double precision,
851  bool remove_duplicates )
852 {
853  typedef QuickHull< RealConvexHullKernel > ConvexHull;
854  typedef typename ConvexCellComplex< RealPoint >::FaceRange FaceRange;
855  ConvexHull hull( precision );
856  hull.setInput( input_points, remove_duplicates );
857  bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
858  cell_complex.clear();
859  if ( ! ok ) return false;
860  // Build complex, only 1 finite cell and as many faces as convex hull facets.
861  // Taking care of faces for each cell (here one cell borders all faces).
862  std::vector< IndexRange > faces;
863  hull.getFacetVertices( faces );
864  FaceRange all_faces;
865  for ( Index i = 0; i < faces.size(); i++ )
866  all_faces.push_back( { i, true } );
867  cell_complex.cell_faces.push_back( all_faces );
868  // Vertices of this unique cell will be computed lazily on request.
869  // Taking care of each face.
870  for ( Index i = 0; i < faces.size(); i++ )
871  {
872  // every inner face borders cell 0
873  cell_complex.true_face_cell.push_back( 0 );
874  // every outer face borders the infinite cell
875  cell_complex.false_face_cell.push_back( cell_complex.infiniteCell() );
876  }
877  // Taking care of vertices (in consistent order) of every face
878  cell_complex.true_face_vertices.swap( faces );
879  // Taking care of vertex positions
880  hull.getVertexPositions( cell_complex.vertex_position );
881  return true;
882 }
883 
884 
885 //-----------------------------------------------------------------------------
886 template < int dim, typename TInteger, typename TInternalInteger >
887 bool
888 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeDelaunayCellComplex
889 ( ConvexCellComplex< RealPoint >& cell_complex,
890  const std::vector< RealPoint >& input_points,
891  double precision,
892  bool remove_duplicates )
893 {
894  typedef QuickHull< RealDelaunayKernel > Delaunay;
895  typedef typename Delaunay::Ridge Ridge;
896  typedef typename ConvexCellComplex< RealPoint >::FaceRange FaceRange;
897 
898  Delaunay del( precision );
899  del.setInput( input_points, remove_duplicates );
900  bool ok = del.computeConvexHull( Delaunay::Status::VerticesCompleted );
901  cell_complex.clear();
902  if ( ! ok ) return false;
903 
904  // Build complex, as many maximal cells as convex hull facets.
905  // convex hull facet -> cell of complex
906  // convex hull ridge -> face of complex
907  // (1) Get cell vertices, count ridges/faces and compute their vertices
908  std::map< Ridge, Index > r2f;
909  computeFacetAndRidgeVertices( del,
910  cell_complex.cell_vertices,
911  r2f,
912  cell_complex.true_face_vertices );
913  // (2) assign ridges/faces to cell and conversely
914  const Index nb_r = r2f.size();
915  cell_complex.true_face_cell .resize( nb_r, cell_complex.infiniteCell() );
916  cell_complex.false_face_cell.resize( nb_r, cell_complex.infiniteCell() );
917  cell_complex.true_face_vertices.resize( nb_r );
918  for ( Index cur_f = 0; cur_f < del.nbFiniteFacets(); ++cur_f ) {
919  const auto& facet = del.facets[ cur_f ];
920  FaceRange current_faces;
921  for ( auto neigh_f : facet.neighbors ) {
922  const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
923  const bool pos = cur_f < neigh_f;
924  const Index cur_r = r2f[ r ];
925  cell_complex.true_face_cell [ cur_r ] = r.first;
926  if ( r.second >= del.nbFiniteFacets() )
927  cell_complex.false_face_cell[ cur_r ] = cell_complex.infiniteCell();
928  else
929  cell_complex.false_face_cell[ cur_r ] = r.second;
930  current_faces.emplace_back( cur_r, pos );
931  }
932  cell_complex.cell_faces.push_back( current_faces );
933  }
934  // (3) Takes care of vertex positions
935  del.getVertexPositions( cell_complex.vertex_position );
936  return true;
937 }
938 
939 
940 //-----------------------------------------------------------------------------
941 template < int dim, typename TInteger, typename TInternalInteger >
942 template < typename QHull >
943 void
944 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeFacetAndRidgeVertices
945 ( const QHull& hull,
946  std::vector< IndexRange >& cell_vertices,
947  std::map< typename QHull::Ridge, Index >& r2f,
948  std::vector< IndexRange >& face_vertices )
949 {
950  typedef typename QHull::Ridge Ridge;
951 
952  ASSERT( hull.status() >= QHull::Status::VerticesCompleted
953  && hull.status() <= QHull::Status::AllCompleted );
954 
955  // Get cell vertices and sort them
956  bool ok_fv = hull.getFacetVertices( cell_vertices );
957  if ( ! ok_fv )
958  trace.error() << "[ConvexityHelper::computeFacetAndRidgeVertices]"
959  << " method hull.getFacetVertices failed."
960  << " Maybe QuickHull was not computed till VerticesCompleted."
961  << std::endl;
962  std::vector< IndexRange > sorted_cell_vertices = cell_vertices;
963  for ( auto& vtcs : sorted_cell_vertices )
964  std::sort( vtcs.begin(), vtcs.end() );
965  cell_vertices.resize( hull.nbFiniteFacets() );
966 
967  // Count ridges/faces and compute their vertices
968  Index nb_r = 0;
969  face_vertices.clear();
970  for ( Index cur_f = 0; cur_f < hull.nbFiniteFacets(); ++cur_f ) {
971  const auto& facet = hull.facets[ cur_f ];
972  for ( auto neigh_f : facet.neighbors ) {
973  const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
974  auto itr = r2f.find( r );
975  if ( itr == r2f.end() ) {
976  IndexRange result;
977  std::set_intersection( sorted_cell_vertices[ cur_f ].cbegin(),
978  sorted_cell_vertices[ cur_f ].cend(),
979  sorted_cell_vertices[ neigh_f ].cbegin(),
980  sorted_cell_vertices[ neigh_f ].cend(),
981  std::back_inserter( result ) );
982  face_vertices.push_back( result );
983  r2f[ r ] = nb_r++;
984  }
985  }
986  }
987 }
988 
989 // //
990 ///////////////////////////////////////////////////////////////////////////////