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.
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.
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/>.
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
24 * Implementation of inline methods defined in ConvexityHelper.h
26 * This file is part of the DGtal library.
30 //////////////////////////////////////////////////////////////////////////////
34 #include "DGtal/kernel/IntegerConverter.h"
35 //////////////////////////////////////////////////////////////////////////////
37 ///////////////////////////////////////////////////////////////////////////////
38 // IMPLEMENTATION of inline methods.
39 ///////////////////////////////////////////////////////////////////////////////
41 ///////////////////////////////////////////////////////////////////////////////
42 // ----------------------- Standard services ------------------------------
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 )
52 if ( remove_duplicates )
55 for ( auto&& p : input_points ) S.insert( p );
56 X = PointRange( S.cbegin(), S.cend() );
58 else X = input_points;
59 LatticePolytope P( X.cbegin(), X.cend() );
60 if ( P.nbHalfSpaces() != 0 )
63 return computeDegeneratedLatticePolytope( X );
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 )
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 ] ) );
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() ) )
96 if ( alignment > 1e-8 ) break;
99 if ( i == input_points.size() )
102 for ( i = 1; i < input_points.size(); i++ )
104 if ( alpha[ i ] < alpha[ a ] ) a = i;
105 if ( alpha[ i ] > alpha[ b ] ) b = i;
108 X[ 0 ] = input_points[ a ];
109 X[ 1 ] = input_points[ b ];
110 return LatticePolytope( X.cbegin(), X.cend() );
112 // at least 2-dimensional
113 ASSERT( dimension > 1 );
114 if ( dimension == 2 )
116 std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
117 << " Weird error: found initial full dimensional simplex" << std::endl;
118 return LatticePolytope();
120 if ( dimension >= 4 )
122 std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
123 << "Degenerated lattice polytope in nD, n >= 4 is not implemented"
125 return LatticePolytope();
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 )
132 std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
133 << "Weird numerical error, u . v != |u| |v| but u x v != 0"
135 return LatticePolytope();
137 // Now the set of input points should be full dimensional.
138 input_points.push_back( input_points[ 0 ] + n );
140 Point l = input_points[ 0 ];
141 Point u = input_points[ 0 ];
142 for ( const auto& p : input_points ) {
146 Domain domain( l, u );
147 // Compute convex 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 );
155 std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
156 << "Weird error in hull.setInitialSimplex" << std::endl;
157 return LatticePolytope();
159 bool ok_hull = hull.computeConvexHull( target );
162 std::cerr << "[ConvexityHelper::computeDegeneratedLatticePolytope]"
163 << "Weird error in hull.computeConvexHull" << std::endl;
164 return LatticePolytope();
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 ) {
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 );
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(),
186 //-----------------------------------------------------------------------------
187 template < int dim, typename TInteger, typename TInternalInteger >
188 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::PointRange
189 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
190 computeDegeneratedConvexHullVertices
191 ( PointRange& input_points )
193 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
194 // input_points is a range of points with no duplicates, but which
195 // seems to be not full dimensional.
196 if ( input_points.size() <= 1 )
198 // At least 1-dimensional
199 std::vector< Vector > basis;
200 std::vector< Integer > alpha;
201 basis.push_back( input_points[ 1 ] - input_points[ 0 ] );
202 const auto n0 = basis[ 0 ].norm();
203 alpha.push_back( Integer( 0 ) );
204 alpha.push_back( basis[ 0 ].dot( basis[ 0 ] ) );
206 while ( i < input_points.size() ) {
207 Vector v = input_points[ i ] - input_points[ 0 ];
208 alpha.push_back( basis[ 0 ].dot( v ) );
209 const auto ni = v.norm();
210 const double alignment =
211 fabs( fabs( NumberTraits< Integer >::castToDouble( alpha.back() ) )
213 if ( alignment > 1e-8 ) break;
216 if ( i == input_points.size() )
219 for ( i = 1; i < input_points.size(); i++ )
221 if ( alpha[ i ] < alpha[ a ] ) a = i;
222 if ( alpha[ i ] > alpha[ b ] ) b = i;
225 X[ 0 ] = input_points[ a ];
226 X[ 1 ] = input_points[ b ];
229 // at least 2-dimensional
230 ASSERT( dimension > 1 );
231 if ( dimension == 2 )
233 std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
234 << " Weird error: found initial full dimensional simplex" << std::endl;
237 if ( dimension >= 4 )
239 std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
240 << "Degenerated lattice polytope in nD, n >= 4 is not implemented"
244 basis.push_back( input_points[ i ] - input_points[ 0 ] );
245 Vector n = detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
246 ::crossProduct( basis[ 0 ], basis[ 1 ] );
247 if ( n == Vector::zero )
249 std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
250 << "Weird numerical error, u . v != |u| |v| but u x v != 0"
254 // Now the set of input points should be full dimensional.
255 input_points.push_back( input_points[ 0 ] + n );
256 // Compute convex hull
258 hull.setInput( input_points, false );
259 const auto target = ConvexHull::Status::VerticesCompleted;
260 IndexRange full_splx = { 0, 1, i, input_points.size() - 1 };
261 bool ok_init = hull.setInitialSimplex( full_splx );
264 std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
265 << "Weird error in hull.setInitialSimplex" << std::endl;
268 bool ok_hull = hull.computeConvexHull( target );
271 std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
272 << "Weird error in hull.computeConvexHull" << std::endl;
275 // Get convex hull vertices and remove top point
277 hull.getVertexPositions( X );
278 const std::size_t nX = X.size();
279 for ( std::size_t j = 0; j < nX; j++ )
280 if ( X[ j ] == input_points.back() )
289 //-----------------------------------------------------------------------------
290 template < int dim, typename TInteger, typename TInternalInteger >
291 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
292 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
293 computeLatticePolytope
294 ( const PointRange& input_points,
295 bool remove_duplicates,
296 bool make_minkowski_summable )
298 typedef typename LatticePolytope::Domain Domain;
299 typedef typename LatticePolytope::HalfSpace PolytopeHalfSpace;
300 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
301 typedef typename ConvexHull::HalfSpace ConvexHullHalfSpace;
302 typedef typename ConvexHull::Ridge Ridge;
303 if ( input_points.empty() ) return LatticePolytope();
304 if ( input_points.size() <= ( dimension + 1) )
305 return computeSimplex( input_points, remove_duplicates );
307 Point l = input_points[ 0 ];
308 Point u = input_points[ 0 ];
309 for ( const auto& p : input_points ) {
313 Domain domain( l, u );
314 // Compute convex hull
316 hull.setInput( input_points, remove_duplicates );
317 const auto target = ( make_minkowski_summable && dimension == 3 )
318 ? ConvexHull::Status::VerticesCompleted
319 : ConvexHull::Status::FacetsCompleted;
320 bool ok = hull.computeConvexHull( target );
321 if ( ! ok ) // set of points is not full dimensional
322 return computeDegeneratedLatticePolytope( hull.points );
323 // Initialize polytope
324 std::vector< ConvexHullHalfSpace > HS;
325 std::vector< PolytopeHalfSpace > PHS;
326 hull.getFacetHalfSpaces( HS );
327 PHS.reserve( HS.size() );
328 for ( auto& H : HS ) {
331 for ( Dimension i = 0; i < dim; ++i )
332 N[ i ] = IntegerConverter< dimension, Integer >::cast( H.internalNormal()[ i ] );
333 nu = IntegerConverter< dimension, Integer >::cast( H.internalIntercept() );
334 PHS.emplace_back( N, nu );
336 if ( make_minkowski_summable && dimension >= 4 )
337 trace.warning() << "[ConvexityHelper::computeLatticePolytope]"
338 << " Not implemented starting from dimension 4."
340 if ( make_minkowski_summable && dimension == 3 )
342 // Compute ridge vertices to add edge constraints.
343 PointRange positions;
344 std::vector< IndexRange > facet_vertices;
345 std::vector< IndexRange > ridge_vertices;
346 std::map< Ridge, Index > ridge2index;
347 hull.getVertexPositions( positions );
348 computeFacetAndRidgeVertices( hull, facet_vertices,
349 ridge2index, ridge_vertices );
350 for ( auto p : ridge2index ) {
351 const auto r = p.first;
352 // Copy by value since PHS may be reallocated during the iteration.
353 const auto U = PHS[ r.first ].N; // normal of facet 1
354 const auto V = PHS[ r.second ].N; // normal of facet 2
355 const auto& S = ridge_vertices[ p.second ]; // vertices along facets 1, 2
356 ASSERT( S.size() == 2 && "Invalid ridge" );
357 const auto& P0 = positions[ S[ 0 ] ];
358 const auto& P1 = positions[ S[ 1 ] ];
359 auto E = P1 - P0; // edge 1, 2
361 detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
362 ::crossProduct( U, V ); // parallel to E
363 ASSERT( E.dot( UxV ) != 0 && "Invalid E / UxV" );
364 if ( E.dot( UxV ) <= 0 ) E = -E; // force correct orientation
366 detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
367 ::crossProduct( U, E ); // edge on facet 1
369 detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
370 ::crossProduct( E, V ); // edge on facet 2
371 ASSERT( E1.dot( U ) == 0 && "Invalid E1 / U" );
372 ASSERT( E1.dot( V ) < 0 && "Invalid E1 / V" );
373 ASSERT( E2.dot( V ) == 0 && "Invalid E1 / V" );
374 ASSERT( E2.dot( U ) < 0 && "Invalid E1 / U" );
375 for ( Dimension k = 0; k < dimension; ++k ) {
376 const auto W = U[ k ] * V - V[ k ] * U;
377 const auto nn1 = W.dot( E1 );
378 const auto nn2 = W.dot( E2 );
379 if ( nn1 > 0 && nn2 > 0 ) {
380 PHS.emplace_back( -W, -W.dot( P0 ) );
381 ASSERT( E1.dot(-W ) < 0 && "Invalid E1 /-W" );
382 ASSERT( E2.dot(-W ) < 0 && "Invalid E2 /-W" );
384 else if ( nn1 < 0 && nn2 < 0 ) {
385 PHS.emplace_back( W, W.dot( P0 ) );
386 ASSERT( E1.dot( W ) < 0 && "Invalid E1 / W" );
387 ASSERT( E2.dot( W ) < 0 && "Invalid E2 / W" );
392 return LatticePolytope( domain, PHS.cbegin(), PHS.cend(),
393 make_minkowski_summable && ( dimension <= 3 ), true );
396 //-----------------------------------------------------------------------------
397 template < int dim, typename TInteger, typename TInternalInteger >
398 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::PointRange
399 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
400 computeConvexHullVertices
401 ( const PointRange& input_points,
402 bool remove_duplicates )
404 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
405 PointRange positions;
407 hull.setInput( input_points, remove_duplicates );
408 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
411 PointRange Z( input_points );
412 return computeDegeneratedConvexHullVertices( Z );
414 hull.getVertexPositions( positions );
418 //-----------------------------------------------------------------------------
419 template < int dim, typename TInteger, typename TInternalInteger >
420 template < typename TSurfaceMesh >
422 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
423 ( TSurfaceMesh& mesh,
424 const PointRange& input_points,
425 bool remove_duplicates )
427 typedef TSurfaceMesh SurfaceMesh;
428 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
430 hull.setInput( input_points, remove_duplicates );
431 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
432 if ( !ok ) return false;
433 std::vector< RealPoint > positions;
434 hull.getVertexPositions( positions );
435 std::vector< IndexRange > faces;
436 hull.getFacetVertices( faces );
437 mesh = SurfaceMesh( positions.cbegin(), positions.cend(),
438 faces.cbegin(), faces.cend() );
442 //-----------------------------------------------------------------------------
443 template < int dim, typename TInteger, typename TInternalInteger >
445 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
446 ( PolygonalSurface< Point >& polysurf,
447 const PointRange& input_points,
448 bool remove_duplicates )
450 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
452 hull.setInput( input_points, remove_duplicates );
453 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
454 if ( !ok ) return false;
455 PointRange positions;
456 hull.getVertexPositions( positions );
457 std::vector< IndexRange > faces;
458 hull.getFacetVertices( faces );
459 // build polygonal surface
461 for ( auto p : positions ) polysurf.addVertex( p );
462 for ( auto f : faces ) polysurf.addPolygonalFace( f );
463 return polysurf.build();
466 //-----------------------------------------------------------------------------
467 template < int dim, typename TInteger, typename TInternalInteger >
469 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullCellComplex
470 ( ConvexCellComplex< Point >& cell_complex,
471 const PointRange& input_points,
472 bool remove_duplicates )
474 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
475 typedef typename ConvexCellComplex< Point >::FaceRange FaceRange;
477 hull.setInput( input_points, remove_duplicates );
478 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
479 cell_complex.clear();
480 if ( ! ok ) return false;
481 // Build complex, only 1 finite cell and as many faces as convex hull facets.
482 // Taking care of faces for each cell (here one cell borders all faces).
483 std::vector< IndexRange > faces;
484 hull.getFacetVertices( faces );
486 for ( Index i = 0; i < faces.size(); i++ )
487 all_faces.push_back( { i, true } );
488 cell_complex.cell_faces.push_back( all_faces );
489 // Vertices of this unique cell will be computed lazily on request.
490 // Taking care of each face.
491 for ( Index i = 0; i < faces.size(); i++ )
493 // every inner face borders cell 0
494 cell_complex.true_face_cell.push_back( 0 );
495 // every outer face borders the infinite cell
496 cell_complex.false_face_cell.push_back( cell_complex.infiniteCell() );
498 // Taking care of vertices (in consistent order) of every face
499 cell_complex.true_face_vertices.swap( faces );
500 // Taking care of vertex positions
501 hull.getVertexPositions( cell_complex.vertex_position );
505 //-----------------------------------------------------------------------------
506 template < int dim, typename TInteger, typename TInternalInteger >
508 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeDelaunayCellComplex
509 ( ConvexCellComplex< Point >& cell_complex,
510 const PointRange& input_points,
511 bool remove_duplicates )
513 typedef QuickHull< LatticeDelaunayKernel > Delaunay;
514 typedef typename Delaunay::Ridge Ridge;
515 typedef typename ConvexCellComplex< Point >::FaceRange FaceRange;
518 del.setInput( input_points, remove_duplicates );
519 bool ok = del.computeConvexHull( Delaunay::Status::VerticesCompleted );
520 cell_complex.clear();
521 if ( ! ok ) return false;
523 // Build complex, as many maximal cells as convex hull facets.
524 // convex hull facet -> cell of complex
525 // convex hull ridge -> face of complex
526 // (1) Get cell vertices, count ridges/faces and compute their vertices
527 std::map< Ridge, Index > r2f;
528 computeFacetAndRidgeVertices( del,
529 cell_complex.cell_vertices,
531 cell_complex.true_face_vertices );
532 // (2) assign ridges/faces to cell and conversely
533 const Index nb_r = r2f.size();
534 cell_complex.true_face_cell .resize( nb_r, cell_complex.infiniteCell() );
535 cell_complex.false_face_cell.resize( nb_r, cell_complex.infiniteCell() );
536 cell_complex.true_face_vertices.resize( nb_r );
537 for ( Index cur_f = 0; cur_f < del.nbFiniteFacets(); ++cur_f ) {
538 const auto& facet = del.facets[ cur_f ];
539 FaceRange current_faces;
540 for ( auto neigh_f : facet.neighbors ) {
541 const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
542 const bool pos = cur_f < neigh_f;
543 const Index cur_r = r2f[ r ];
544 cell_complex.true_face_cell [ cur_r ] = r.first;
545 if ( r.second >= del.nbFiniteFacets() )
546 cell_complex.false_face_cell[ cur_r ] = cell_complex.infiniteCell();
548 cell_complex.false_face_cell[ cur_r ] = r.second;
549 current_faces.emplace_back( cur_r, pos );
551 cell_complex.cell_faces.push_back( current_faces );
553 // (3) Takes care of vertex positions
554 del.getVertexPositions( cell_complex.vertex_position );
559 //-----------------------------------------------------------------------------
560 template < int dim, typename TInteger, typename TInternalInteger >
561 typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::RationalPolytope
562 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeRationalPolytope
563 ( const std::vector< RealPoint >& input_points,
565 bool remove_duplicates,
566 bool make_minkowski_summable )
568 if ( denominator < 1 )
569 trace.error() << "Invalid denominator " << denominator
570 << ". Should be greater or equal to 1." << std::endl;
571 typedef typename RationalPolytope::Domain Domain;
572 typedef typename RationalPolytope::HalfSpace PolytopeHalfSpace;
573 typedef QuickHull< RealConvexHullKernel > ConvexHull;
574 typedef typename ConvexHull::HalfSpace ConvexHullHalfSpace;
575 typedef typename ConvexHull::Ridge Ridge;
576 if ( input_points.empty() ) return RationalPolytope();
577 // Compute convex hull
578 ConvexHull hull( denominator );
579 hull.setInput( input_points, remove_duplicates );
580 const auto target = ( make_minkowski_summable && dimension == 3 )
581 ? ConvexHull::Status::VerticesCompleted
582 : ConvexHull::Status::FacetsCompleted;
583 bool ok = hull.computeConvexHull( target );
584 if ( ! ok ) return RationalPolytope();
585 // Compute domain (as a lattice domain)
586 auto l = hull.points[ 0 ];
587 auto u = hull.points[ 0 ];
588 for ( const auto& p : hull.points ) {
592 Domain domain( l, u );
593 trace.info() << "Domain l=" << l << " u=" << u << std::endl;
594 // Initialize polytope
595 std::vector< ConvexHullHalfSpace > HS;
596 std::vector< PolytopeHalfSpace > PHS;
597 hull.getFacetHalfSpaces( HS );
598 PHS.reserve( HS.size() );
599 for ( auto& H : HS ) {
602 for ( Dimension i = 0; i < dim; ++i )
603 N[ i ] = (Integer) H.internalNormal()[ i ];
604 nu = (Integer) H.internalIntercept();
605 PHS.emplace_back( N, nu );
607 if ( make_minkowski_summable && dimension >= 4 )
608 trace.warning() << "[ConvexityHelper::computeRationalPolytope]"
609 << " Not implemented starting from dimension 4."
611 if ( make_minkowski_summable && dimension == 3 )
613 // Compute ridge vertices to add edge constraints.
614 PointRange positions;
615 std::vector< IndexRange > facet_vertices;
616 std::vector< IndexRange > ridge_vertices;
617 std::map< Ridge, Index > ridge2index;
618 hull.getVertexPositions( positions );
619 computeFacetAndRidgeVertices( hull, facet_vertices,
620 ridge2index, ridge_vertices );
621 for ( auto p : ridge2index ) {
622 const auto r = p.first;
623 // Copy by value since PHS may be reallocated during the iteration.
624 const auto U = PHS[ r.first ].N; // normal of facet 1
625 const auto V = PHS[ r.second ].N; // normal of facet 2
626 const auto& S = ridge_vertices[ p.second ]; // vertices along facets 1, 2
627 ASSERT( S.size() == 2 && "Invalid ridge" );
628 const auto& P0 = positions[ S[ 0 ] ];
629 const auto& P1 = positions[ S[ 1 ] ];
630 auto E = P1 - P0; // edge 1, 2
632 detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
633 ::crossProduct( U, V ); // parallel to E
634 ASSERT( E.dot( UxV ) != 0 && "Invalid E / UxV" );
635 if ( E.dot( UxV ) <= 0 ) E = -E; // force correct orientation
637 detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
638 ::crossProduct( U, E ); // edge on facet 1
640 detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
641 ::crossProduct( E, V ); // edge on facet 2
642 ASSERT( E1.dot( U ) == 0 && "Invalid E1 / U" );
643 ASSERT( E1.dot( V ) < 0 && "Invalid E1 / V" );
644 ASSERT( E2.dot( V ) == 0 && "Invalid E1 / V" );
645 ASSERT( E2.dot( U ) < 0 && "Invalid E1 / U" );
646 for ( Dimension k = 0; k < dimension; ++k ) {
647 const auto W = U[ k ] * V - V[ k ] * U;
648 const auto nn1 = W.dot( E1 );
649 const auto nn2 = W.dot( E2 );
650 if ( nn1 > 0 && nn2 > 0 ) {
651 PHS.emplace_back( -W, -W.dot( P0 ) );
652 ASSERT( E1.dot(-W ) < 0 && "Invalid E1 /-W" );
653 ASSERT( E2.dot(-W ) < 0 && "Invalid E2 /-W" );
655 else if ( nn1 < 0 && nn2 < 0 ) {
656 PHS.emplace_back( W, W.dot( P0 ) );
657 ASSERT( E1.dot( W ) < 0 && "Invalid E1 / W" );
658 ASSERT( E2.dot( W ) < 0 && "Invalid E2 / W" );
663 return RationalPolytope( denominator, domain, PHS.cbegin(), PHS.cend(),
664 make_minkowski_summable && ( dimension <= 3 ), true );
668 //-----------------------------------------------------------------------------
669 template < int dim, typename TInteger, typename TInternalInteger >
670 template < typename TSurfaceMesh >
672 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
673 ( TSurfaceMesh& mesh,
674 const std::vector< RealPoint >& input_points,
676 bool remove_duplicates )
678 typedef TSurfaceMesh SurfaceMesh;
679 typedef QuickHull< RealConvexHullKernel > ConvexHull;
680 ConvexHull hull( precision );
681 hull.setInput( input_points, remove_duplicates );
682 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
683 if ( !ok ) return false;
684 std::vector< RealPoint > positions;
685 hull.getVertexPositions( positions );
686 std::vector< IndexRange > faces;
687 hull.getFacetVertices( faces );
688 mesh = SurfaceMesh( positions.cbegin(), positions.cend(),
689 faces.cbegin(), faces.cend() );
694 //-----------------------------------------------------------------------------
695 template < int dim, typename TInteger, typename TInternalInteger >
697 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
698 ( PolygonalSurface< RealPoint >& polysurf,
699 const std::vector< RealPoint >& input_points,
701 bool remove_duplicates )
703 typedef QuickHull< RealConvexHullKernel > ConvexHull;
704 ConvexHull hull( precision );
705 hull.setInput( input_points, remove_duplicates );
706 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
707 if ( !ok ) return false;
708 std::vector< RealPoint > positions;
709 hull.getVertexPositions( positions );
710 std::vector< IndexRange > faces;
711 hull.getFacetVertices( faces );
712 // build polygonal surface
714 for ( auto p : positions ) polysurf.addVertex( p );
715 for ( auto f : faces ) polysurf.addPolygonalFace( f );
716 return polysurf.build();
719 //-----------------------------------------------------------------------------
720 template < int dim, typename TInteger, typename TInternalInteger >
722 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullCellComplex
723 ( ConvexCellComplex< RealPoint >& cell_complex,
724 const std::vector< RealPoint >& input_points,
726 bool remove_duplicates )
728 typedef QuickHull< RealConvexHullKernel > ConvexHull;
729 typedef typename ConvexCellComplex< RealPoint >::FaceRange FaceRange;
730 ConvexHull hull( precision );
731 hull.setInput( input_points, remove_duplicates );
732 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
733 cell_complex.clear();
734 if ( ! ok ) return false;
735 // Build complex, only 1 finite cell and as many faces as convex hull facets.
736 // Taking care of faces for each cell (here one cell borders all faces).
737 std::vector< IndexRange > faces;
738 hull.getFacetVertices( faces );
740 for ( Index i = 0; i < faces.size(); i++ )
741 all_faces.push_back( { i, true } );
742 cell_complex.cell_faces.push_back( all_faces );
743 // Vertices of this unique cell will be computed lazily on request.
744 // Taking care of each face.
745 for ( Index i = 0; i < faces.size(); i++ )
747 // every inner face borders cell 0
748 cell_complex.true_face_cell.push_back( 0 );
749 // every outer face borders the infinite cell
750 cell_complex.false_face_cell.push_back( cell_complex.infiniteCell() );
752 // Taking care of vertices (in consistent order) of every face
753 cell_complex.true_face_vertices.swap( faces );
754 // Taking care of vertex positions
755 hull.getVertexPositions( cell_complex.vertex_position );
760 //-----------------------------------------------------------------------------
761 template < int dim, typename TInteger, typename TInternalInteger >
763 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeDelaunayCellComplex
764 ( ConvexCellComplex< RealPoint >& cell_complex,
765 const std::vector< RealPoint >& input_points,
767 bool remove_duplicates )
769 typedef QuickHull< RealDelaunayKernel > Delaunay;
770 typedef typename Delaunay::Ridge Ridge;
771 typedef typename ConvexCellComplex< RealPoint >::FaceRange FaceRange;
773 Delaunay del( precision );
774 del.setInput( input_points, remove_duplicates );
775 bool ok = del.computeConvexHull( Delaunay::Status::VerticesCompleted );
776 cell_complex.clear();
777 if ( ! ok ) return false;
779 // Build complex, as many maximal cells as convex hull facets.
780 // convex hull facet -> cell of complex
781 // convex hull ridge -> face of complex
782 // (1) Get cell vertices, count ridges/faces and compute their vertices
783 std::map< Ridge, Index > r2f;
784 computeFacetAndRidgeVertices( del,
785 cell_complex.cell_vertices,
787 cell_complex.true_face_vertices );
788 // (2) assign ridges/faces to cell and conversely
789 const Index nb_r = r2f.size();
790 cell_complex.true_face_cell .resize( nb_r, cell_complex.infiniteCell() );
791 cell_complex.false_face_cell.resize( nb_r, cell_complex.infiniteCell() );
792 cell_complex.true_face_vertices.resize( nb_r );
793 for ( Index cur_f = 0; cur_f < del.nbFiniteFacets(); ++cur_f ) {
794 const auto& facet = del.facets[ cur_f ];
795 FaceRange current_faces;
796 for ( auto neigh_f : facet.neighbors ) {
797 const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
798 const bool pos = cur_f < neigh_f;
799 const Index cur_r = r2f[ r ];
800 cell_complex.true_face_cell [ cur_r ] = r.first;
801 if ( r.second >= del.nbFiniteFacets() )
802 cell_complex.false_face_cell[ cur_r ] = cell_complex.infiniteCell();
804 cell_complex.false_face_cell[ cur_r ] = r.second;
805 current_faces.emplace_back( cur_r, pos );
807 cell_complex.cell_faces.push_back( current_faces );
809 // (3) Takes care of vertex positions
810 del.getVertexPositions( cell_complex.vertex_position );
815 //-----------------------------------------------------------------------------
816 template < int dim, typename TInteger, typename TInternalInteger >
817 template < typename QHull >
819 DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeFacetAndRidgeVertices
821 std::vector< IndexRange >& cell_vertices,
822 std::map< typename QHull::Ridge, Index >& r2f,
823 std::vector< IndexRange >& face_vertices )
825 typedef typename QHull::Ridge Ridge;
827 ASSERT( hull.status() >= QHull::Status::VerticesCompleted
828 && hull.status() <= QHull::Status::AllCompleted );
830 // Get cell vertices and sort them
831 bool ok_fv = hull.getFacetVertices( cell_vertices );
833 trace.error() << "[ConvexityHelper::computeFacetAndRidgeVertices]"
834 << " method hull.getFacetVertices failed."
835 << " Maybe QuickHull was not computed till VerticesCompleted."
837 std::vector< IndexRange > sorted_cell_vertices = cell_vertices;
838 for ( auto& vtcs : sorted_cell_vertices )
839 std::sort( vtcs.begin(), vtcs.end() );
840 cell_vertices.resize( hull.nbFiniteFacets() );
842 // Count ridges/faces and compute their vertices
844 face_vertices.clear();
845 for ( Index cur_f = 0; cur_f < hull.nbFiniteFacets(); ++cur_f ) {
846 const auto& facet = hull.facets[ cur_f ];
847 for ( auto neigh_f : facet.neighbors ) {
848 const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
849 auto itr = r2f.find( r );
850 if ( itr == r2f.end() ) {
852 std::set_intersection( sorted_cell_vertices[ cur_f ].cbegin(),
853 sorted_cell_vertices[ cur_f ].cend(),
854 sorted_cell_vertices[ neigh_f ].cbegin(),
855 sorted_cell_vertices[ neigh_f ].cend(),
856 std::back_inserter( result ) );
857 face_vertices.push_back( result );
865 ///////////////////////////////////////////////////////////////////////////////