DGtal 1.3.0
Loading...
Searching...
No Matches
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//-----------------------------------------------------------------------------
45template < int dim, typename TInteger, typename TInternalInteger >
46typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
47DGtal::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//-----------------------------------------------------------------------------
67template < int dim, typename TInteger, typename TInternalInteger >
68typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
69DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
70computeDegeneratedLatticePolytope
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::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 PointRange 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//-----------------------------------------------------------------------------
188template < int dim, typename TInteger, typename TInternalInteger >
189typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::PointRange
190DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
191computeDegeneratedConvexHullVertices
192( PointRange& input_points )
193{
194 typedef typename LatticePolytope::Domain Domain;
195 typedef typename LatticePolytope::HalfSpace PolytopeHalfSpace;
196 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
197 typedef typename ConvexHull::IndexRange IndexRange;
198 typedef typename ConvexHull::HalfSpace ConvexHullHalfSpace;
199 // input_points is a range of points with no duplicates, but which
200 // seems to be not full dimensional.
201 if ( input_points.size() <= 1 )
202 return input_points;
203 // At least 1-dimensional
204 std::vector< Vector > basis;
205 std::vector< Integer > alpha;
206 basis.push_back( input_points[ 1 ] - input_points[ 0 ] );
207 const auto n0 = basis[ 0 ].norm();
208 alpha.push_back( Integer( 0 ) );
209 alpha.push_back( basis[ 0 ].dot( basis[ 0 ] ) );
210 Index i = 2;
211 while ( i < input_points.size() ) {
212 Vector v = input_points[ i ] - input_points[ 0 ];
213 alpha.push_back( basis[ 0 ].dot( v ) );
214 const auto ni = v.norm();
215 const double alignment =
216 fabs( fabs( NumberTraits< Integer >::castToDouble( alpha.back() ) )
217 - ( n0 * ni ) );
218 if ( alignment > 1e-8 ) break;
219 i++;
220 }
221 if ( i == input_points.size() )
222 { // 1-dimensional
223 Index a = 0, b = 0;
224 for ( i = 1; i < input_points.size(); i++ )
225 {
226 if ( alpha[ i ] < alpha[ a ] ) a = i;
227 if ( alpha[ i ] > alpha[ b ] ) b = i;
228 }
229 PointRange X( 2 );
230 X[ 0 ] = input_points[ a ];
231 X[ 1 ] = input_points[ b ];
232 return X;
233 }
234 // at least 2-dimensional
235 ASSERT( dimension > 1 );
236 if ( dimension == 2 )
237 {
238 std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
239 << " Weird error: found initial full dimensional simplex" << std::endl;
240 return PointRange();
241 }
242 if ( dimension >= 4 )
243 {
244 std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
245 << "Degenerated lattice polytope in nD, n >= 4 is not implemented"
246 << std::endl;
247 return PointRange();
248 }
249 basis.push_back( input_points[ i ] - input_points[ 0 ] );
250 Vector n = detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
251 ::crossProduct( basis[ 0 ], basis[ 1 ] );
252 if ( n == Vector::zero )
253 {
254 std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
255 << "Weird numerical error, u . v != |u| |v| but u x v != 0"
256 << std::endl;
257 return PointRange();
258 }
259 // Now the set of input points should be full dimensional.
260 input_points.push_back( input_points[ 0 ] + n );
261 // Compute convex hull
262 ConvexHull hull;
263 hull.setInput( input_points, false );
264 const auto target = ConvexHull::Status::VerticesCompleted;
265 IndexRange full_splx = { 0, 1, i, input_points.size() - 1 };
266 bool ok_init = hull.setInitialSimplex( full_splx );
267 if ( ! ok_init )
268 {
269 std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
270 << "Weird error in hull.setInitialSimplex" << std::endl;
271 return PointRange();
272 }
273 bool ok_hull = hull.computeConvexHull( target );
274 if ( ! ok_hull )
275 {
276 std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
277 << "Weird error in hull.computeConvexHull" << std::endl;
278 return PointRange();
279 }
280 // Get convex hull vertices and remove top point
281 PointRange X;
282 hull.getVertexPositions( X );
283 const std::size_t nX = X.size();
284 for ( std::size_t j = 0; j < nX; j++ )
285 if ( X[ j ] == input_points.back() )
286 {
287 X[ j ] = X.back();
288 break;
289 }
290 X.resize( nX - 1 );
291 return X;
292}
293
294//-----------------------------------------------------------------------------
295template < int dim, typename TInteger, typename TInternalInteger >
296typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
297DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
298computeLatticePolytope
299( const PointRange& input_points,
300 bool remove_duplicates,
301 bool make_minkowski_summable )
302{
303 typedef typename LatticePolytope::Domain Domain;
304 typedef typename LatticePolytope::HalfSpace PolytopeHalfSpace;
305 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
306 typedef typename ConvexHull::HalfSpace ConvexHullHalfSpace;
307 typedef typename ConvexHull::Ridge Ridge;
308 if ( input_points.empty() ) return LatticePolytope();
309 if ( input_points.size() <= ( dimension + 1) )
310 return computeSimplex( input_points, remove_duplicates );
311 // Compute domain
312 Point l = input_points[ 0 ];
313 Point u = input_points[ 0 ];
314 for ( const auto& p : input_points ) {
315 l = l.inf( p );
316 u = u.sup( p );
317 }
318 Domain domain( l, u );
319 // Compute convex hull
320 ConvexHull hull;
321 hull.setInput( input_points, remove_duplicates );
322 const auto target = ( make_minkowski_summable && dimension == 3 )
323 ? ConvexHull::Status::VerticesCompleted
324 : ConvexHull::Status::FacetsCompleted;
325 bool ok = hull.computeConvexHull( target );
326 if ( ! ok ) // set of points is not full dimensional
327 return computeDegeneratedLatticePolytope( hull.points );
328 // Initialize polytope
329 std::vector< ConvexHullHalfSpace > HS;
330 std::vector< PolytopeHalfSpace > PHS;
331 hull.getFacetHalfSpaces( HS );
332 PHS.reserve( HS.size() );
333 for ( auto& H : HS ) {
334 Vector N;
335 Integer nu;
336 for ( Dimension i = 0; i < dim; ++i )
337 N[ i ] = IntegerConverter< dimension, Integer >::cast( H.internalNormal()[ i ] );
338 nu = IntegerConverter< dimension, Integer >::cast( H.internalIntercept() );
339 PHS.emplace_back( N, nu );
340 }
341 if ( make_minkowski_summable && dimension >= 4 )
342 trace.warning() << "[ConvexityHelper::computeLatticePolytope]"
343 << " Not implemented starting from dimension 4."
344 << std::endl;
345 if ( make_minkowski_summable && dimension == 3 )
346 {
347 // Compute ridge vertices to add edge constraints.
348 PointRange positions;
349 std::vector< IndexRange > facet_vertices;
350 std::vector< IndexRange > ridge_vertices;
351 std::map< Ridge, Index > ridge2index;
352 hull.getVertexPositions( positions );
353 computeFacetAndRidgeVertices( hull, facet_vertices,
354 ridge2index, ridge_vertices );
355 for ( auto p : ridge2index ) {
356 const auto r = p.first;
357 // Copy by value since PHS may be reallocated during the iteration.
358 const auto U = PHS[ r.first ].N; // normal of facet 1
359 const auto V = PHS[ r.second ].N; // normal of facet 2
360 const auto& S = ridge_vertices[ p.second ]; // vertices along facets 1, 2
361 ASSERT( S.size() == 2 && "Invalid ridge" );
362 const auto& P0 = positions[ S[ 0 ] ];
363 const auto& P1 = positions[ S[ 1 ] ];
364 auto E = P1 - P0; // edge 1, 2
365 const auto UxV =
366 detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
367 ::crossProduct( U, V ); // parallel to E
368 ASSERT( E.dot( UxV ) != 0 && "Invalid E / UxV" );
369 if ( E.dot( UxV ) <= 0 ) E = -E; // force correct orientation
370 const auto E1 =
371 detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
372 ::crossProduct( U, E ); // edge on facet 1
373 const auto E2 =
374 detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
375 ::crossProduct( E, V ); // edge on facet 2
376 ASSERT( E1.dot( U ) == 0 && "Invalid E1 / U" );
377 ASSERT( E1.dot( V ) < 0 && "Invalid E1 / V" );
378 ASSERT( E2.dot( V ) == 0 && "Invalid E1 / V" );
379 ASSERT( E2.dot( U ) < 0 && "Invalid E1 / U" );
380 for ( Dimension k = 0; k < dimension; ++k ) {
381 const auto W = U[ k ] * V - V[ k ] * U;
382 const auto nn1 = W.dot( E1 );
383 const auto nn2 = W.dot( E2 );
384 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" );
388 }
389 else if ( nn1 < 0 && nn2 < 0 ) {
390 PHS.emplace_back( W, W.dot( P0 ) );
391 ASSERT( E1.dot( W ) < 0 && "Invalid E1 / W" );
392 ASSERT( E2.dot( W ) < 0 && "Invalid E2 / W" );
393 }
394 }
395 }
396 }
397 return LatticePolytope( domain, PHS.cbegin(), PHS.cend(),
398 make_minkowski_summable && ( dimension <= 3 ), true );
399}
400
401//-----------------------------------------------------------------------------
402template < int dim, typename TInteger, typename TInternalInteger >
403typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::PointRange
404DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
405computeConvexHullVertices
406( const PointRange& input_points,
407 bool remove_duplicates )
408{
409 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
410 typedef typename ConvexHull::IndexRange IndexRange;
411 PointRange positions;
412 ConvexHull hull;
413 hull.setInput( input_points, remove_duplicates );
414 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
415 if ( !ok )
416 {
417 PointRange Z( input_points );
418 return computeDegeneratedConvexHullVertices( Z );
419 }
420 hull.getVertexPositions( positions );
421 return positions;
422}
423
424//-----------------------------------------------------------------------------
425template < int dim, typename TInteger, typename TInternalInteger >
426template < typename TSurfaceMesh >
427bool
428DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
429( TSurfaceMesh& mesh,
430 const PointRange& input_points,
431 bool remove_duplicates )
432{
433 typedef TSurfaceMesh SurfaceMesh;
434 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
435 typedef typename ConvexHull::IndexRange IndexRange;
436 ConvexHull hull;
437 hull.setInput( input_points, remove_duplicates );
438 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
439 if ( !ok ) return false;
440 std::vector< RealPoint > positions;
441 hull.getVertexPositions( positions );
442 std::vector< IndexRange > faces;
443 hull.getFacetVertices( faces );
444 mesh = SurfaceMesh( positions.cbegin(), positions.cend(),
445 faces.cbegin(), faces.cend() );
446 return true;
447}
448
449//-----------------------------------------------------------------------------
450template < int dim, typename TInteger, typename TInternalInteger >
451bool
452DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
453( PolygonalSurface< Point >& polysurf,
454 const PointRange& input_points,
455 bool remove_duplicates )
456{
457 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
458 typedef typename ConvexHull::IndexRange IndexRange;
459 ConvexHull hull;
460 hull.setInput( input_points, remove_duplicates );
461 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
462 if ( !ok ) return false;
463 PointRange positions;
464 hull.getVertexPositions( positions );
465 std::vector< IndexRange > faces;
466 hull.getFacetVertices( faces );
467 // build polygonal surface
468 polysurf.clear();
469 for ( auto p : positions ) polysurf.addVertex( p );
470 for ( auto f : faces ) polysurf.addPolygonalFace( f );
471 return polysurf.build();
472}
473
474//-----------------------------------------------------------------------------
475template < int dim, typename TInteger, typename TInternalInteger >
476bool
477DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullCellComplex
478( ConvexCellComplex< Point >& cell_complex,
479 const PointRange& input_points,
480 bool remove_duplicates )
481{
482 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
483 typedef typename ConvexHull::IndexRange IndexRange;
484 typedef typename ConvexCellComplex< Point >::FaceRange FaceRange;
485 ConvexHull hull;
486 hull.setInput( input_points, remove_duplicates );
487 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
488 cell_complex.clear();
489 if ( ! ok ) return false;
490 // Build complex, only 1 finite cell and as many faces as convex hull facets.
491 // Taking care of faces for each cell (here one cell borders all faces).
492 std::vector< IndexRange > faces;
493 hull.getFacetVertices( faces );
494 FaceRange all_faces;
495 for ( Index i = 0; i < faces.size(); i++ )
496 all_faces.push_back( { i, true } );
497 cell_complex.cell_faces.push_back( all_faces );
498 // Vertices of this unique cell will be computed lazily on request.
499 // Taking care of each face.
500 for ( const auto& vtcs : faces )
501 {
502 // every inner face borders cell 0
503 cell_complex.true_face_cell.push_back( 0 );
504 // every outer face borders the infinite cell
505 cell_complex.false_face_cell.push_back( cell_complex.infiniteCell() );
506 }
507 // Taking care of vertices (in consistent order) of every face
508 cell_complex.true_face_vertices.swap( faces );
509 // Taking care of vertex positions
510 hull.getVertexPositions( cell_complex.vertex_position );
511 return true;
512}
513
514//-----------------------------------------------------------------------------
515template < int dim, typename TInteger, typename TInternalInteger >
516bool
517DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeDelaunayCellComplex
518( ConvexCellComplex< Point >& cell_complex,
519 const PointRange& input_points,
520 bool remove_duplicates )
521{
522 typedef QuickHull< LatticeDelaunayKernel > Delaunay;
523 typedef typename Delaunay::Ridge Ridge;
524 typedef typename ConvexCellComplex< Point >::FaceRange FaceRange;
525
526 Delaunay del;
527 del.setInput( input_points, remove_duplicates );
528 bool ok = del.computeConvexHull( Delaunay::Status::VerticesCompleted );
529 cell_complex.clear();
530 if ( ! ok ) return false;
531
532 // Build complex, as many maximal cells as convex hull facets.
533 // convex hull facet -> cell of complex
534 // convex hull ridge -> face of complex
535 // (1) Get cell vertices, count ridges/faces and compute their vertices
536 std::map< Ridge, Index > r2f;
537 computeFacetAndRidgeVertices( del,
538 cell_complex.cell_vertices,
539 r2f,
540 cell_complex.true_face_vertices );
541 // (2) assign ridges/faces to cell and conversely
542 const Index nb_r = r2f.size();
543 cell_complex.true_face_cell .resize( nb_r, cell_complex.infiniteCell() );
544 cell_complex.false_face_cell.resize( nb_r, cell_complex.infiniteCell() );
545 cell_complex.true_face_vertices.resize( nb_r );
546 for ( Index cur_f = 0; cur_f < del.nbFiniteFacets(); ++cur_f ) {
547 const auto& facet = del.facets[ cur_f ];
548 FaceRange current_faces;
549 for ( auto neigh_f : facet.neighbors ) {
550 const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
551 const bool pos = cur_f < neigh_f;
552 const Index cur_r = r2f[ r ];
553 cell_complex.true_face_cell [ cur_r ] = r.first;
554 if ( r.second >= del.nbFiniteFacets() )
555 cell_complex.false_face_cell[ cur_r ] = cell_complex.infiniteCell();
556 else
557 cell_complex.false_face_cell[ cur_r ] = r.second;
558 current_faces.emplace_back( cur_r, pos );
559 }
560 cell_complex.cell_faces.push_back( current_faces );
561 }
562 // (3) Takes care of vertex positions
563 del.getVertexPositions( cell_complex.vertex_position );
564 return true;
565}
566
567
568//-----------------------------------------------------------------------------
569template < int dim, typename TInteger, typename TInternalInteger >
570typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::RationalPolytope
571DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeRationalPolytope
572( const std::vector< RealPoint >& input_points,
573 Integer denominator,
574 bool remove_duplicates,
575 bool make_minkowski_summable )
576{
577 if ( denominator < 1 )
578 trace.error() << "Invalid denominator " << denominator
579 << ". Should be greater or equal to 1." << std::endl;
580 typedef typename RationalPolytope::Domain Domain;
581 typedef typename RationalPolytope::HalfSpace PolytopeHalfSpace;
582 typedef QuickHull< RealConvexHullKernel > ConvexHull;
583 typedef typename ConvexHull::HalfSpace ConvexHullHalfSpace;
584 typedef typename ConvexHull::Ridge Ridge;
585 if ( input_points.empty() ) return RationalPolytope();
586 // Compute convex hull
587 ConvexHull hull( denominator );
588 hull.setInput( input_points, remove_duplicates );
589 const auto target = ( make_minkowski_summable && dimension == 3 )
590 ? ConvexHull::Status::VerticesCompleted
591 : ConvexHull::Status::FacetsCompleted;
592 bool ok = hull.computeConvexHull( target );
593 if ( ! ok ) return RationalPolytope();
594 // Compute domain (as a lattice domain)
595 auto l = hull.points[ 0 ];
596 auto u = hull.points[ 0 ];
597 for ( const auto& p : hull.points ) {
598 l = l.inf( p );
599 u = u.sup( p );
600 }
601 Domain domain( l, u );
602 trace.info() << "Domain l=" << l << " u=" << u << std::endl;
603 // Initialize polytope
604 std::vector< ConvexHullHalfSpace > HS;
605 std::vector< PolytopeHalfSpace > PHS;
606 hull.getFacetHalfSpaces( HS );
607 PHS.reserve( HS.size() );
608 for ( auto& H : HS ) {
609 Vector N;
610 Integer nu;
611 for ( Dimension i = 0; i < dim; ++i )
612 N[ i ] = (Integer) H.internalNormal()[ i ];
613 nu = (Integer) H.internalIntercept();
614 PHS.emplace_back( N, nu );
615 }
616 if ( make_minkowski_summable && dimension >= 4 )
617 trace.warning() << "[ConvexityHelper::computeRationalPolytope]"
618 << " Not implemented starting from dimension 4."
619 << std::endl;
620 if ( make_minkowski_summable && dimension == 3 )
621 {
622 // Compute ridge vertices to add edge constraints.
623 PointRange positions;
624 std::vector< IndexRange > facet_vertices;
625 std::vector< IndexRange > ridge_vertices;
626 std::map< Ridge, Index > ridge2index;
627 hull.getVertexPositions( positions );
628 computeFacetAndRidgeVertices( hull, facet_vertices,
629 ridge2index, ridge_vertices );
630 for ( auto p : ridge2index ) {
631 const auto r = p.first;
632 // Copy by value since PHS may be reallocated during the iteration.
633 const auto U = PHS[ r.first ].N; // normal of facet 1
634 const auto V = PHS[ r.second ].N; // normal of facet 2
635 const auto& S = ridge_vertices[ p.second ]; // vertices along facets 1, 2
636 ASSERT( S.size() == 2 && "Invalid ridge" );
637 const auto& P0 = positions[ S[ 0 ] ];
638 const auto& P1 = positions[ S[ 1 ] ];
639 auto E = P1 - P0; // edge 1, 2
640 const auto UxV =
641 detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
642 ::crossProduct( U, V ); // parallel to E
643 ASSERT( E.dot( UxV ) != 0 && "Invalid E / UxV" );
644 if ( E.dot( UxV ) <= 0 ) E = -E; // force correct orientation
645 const auto E1 =
646 detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
647 ::crossProduct( U, E ); // edge on facet 1
648 const auto E2 =
649 detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
650 ::crossProduct( E, V ); // edge on facet 2
651 ASSERT( E1.dot( U ) == 0 && "Invalid E1 / U" );
652 ASSERT( E1.dot( V ) < 0 && "Invalid E1 / V" );
653 ASSERT( E2.dot( V ) == 0 && "Invalid E1 / V" );
654 ASSERT( E2.dot( U ) < 0 && "Invalid E1 / U" );
655 for ( Dimension k = 0; k < dimension; ++k ) {
656 const auto W = U[ k ] * V - V[ k ] * U;
657 const auto nn1 = W.dot( E1 );
658 const auto nn2 = W.dot( E2 );
659 if ( nn1 > 0 && nn2 > 0 ) {
660 PHS.emplace_back( -W, -W.dot( P0 ) );
661 ASSERT( E1.dot(-W ) < 0 && "Invalid E1 /-W" );
662 ASSERT( E2.dot(-W ) < 0 && "Invalid E2 /-W" );
663 }
664 else if ( nn1 < 0 && nn2 < 0 ) {
665 PHS.emplace_back( W, W.dot( P0 ) );
666 ASSERT( E1.dot( W ) < 0 && "Invalid E1 / W" );
667 ASSERT( E2.dot( W ) < 0 && "Invalid E2 / W" );
668 }
669 }
670 }
671 }
672 return RationalPolytope( denominator, domain, PHS.cbegin(), PHS.cend(),
673 make_minkowski_summable && ( dimension <= 3 ), true );
674}
675
676
677//-----------------------------------------------------------------------------
678template < int dim, typename TInteger, typename TInternalInteger >
679template < typename TSurfaceMesh >
680bool
681DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
682( TSurfaceMesh& mesh,
683 const std::vector< RealPoint >& input_points,
684 double precision,
685 bool remove_duplicates )
686{
687 typedef TSurfaceMesh SurfaceMesh;
688 typedef QuickHull< RealConvexHullKernel > ConvexHull;
689 typedef typename ConvexHull::IndexRange IndexRange;
690 ConvexHull hull( precision );
691 hull.setInput( input_points, remove_duplicates );
692 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
693 if ( !ok ) return false;
694 std::vector< RealPoint > positions;
695 hull.getVertexPositions( positions );
696 std::vector< IndexRange > faces;
697 hull.getFacetVertices( faces );
698 mesh = SurfaceMesh( positions.cbegin(), positions.cend(),
699 faces.cbegin(), faces.cend() );
700 return true;
701}
702
703
704//-----------------------------------------------------------------------------
705template < int dim, typename TInteger, typename TInternalInteger >
706bool
707DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
708( PolygonalSurface< RealPoint >& polysurf,
709 const std::vector< RealPoint >& input_points,
710 double precision,
711 bool remove_duplicates )
712{
713 typedef QuickHull< RealConvexHullKernel > ConvexHull;
714 typedef typename ConvexHull::IndexRange IndexRange;
715 ConvexHull hull( precision );
716 hull.setInput( input_points, remove_duplicates );
717 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
718 if ( !ok ) return false;
719 std::vector< RealPoint > positions;
720 hull.getVertexPositions( positions );
721 std::vector< IndexRange > faces;
722 hull.getFacetVertices( faces );
723 // build polygonal surface
724 polysurf.clear();
725 for ( auto p : positions ) polysurf.addVertex( p );
726 for ( auto f : faces ) polysurf.addPolygonalFace( f );
727 return polysurf.build();
728}
729
730//-----------------------------------------------------------------------------
731template < int dim, typename TInteger, typename TInternalInteger >
732bool
733DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullCellComplex
734( ConvexCellComplex< RealPoint >& cell_complex,
735 const std::vector< RealPoint >& input_points,
736 double precision,
737 bool remove_duplicates )
738{
739 typedef QuickHull< RealConvexHullKernel > ConvexHull;
740 typedef typename ConvexHull::IndexRange IndexRange;
741 typedef typename ConvexCellComplex< RealPoint >::FaceRange FaceRange;
742 ConvexHull hull( precision );
743 hull.setInput( input_points, remove_duplicates );
744 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
745 cell_complex.clear();
746 if ( ! ok ) return false;
747 // Build complex, only 1 finite cell and as many faces as convex hull facets.
748 // Taking care of faces for each cell (here one cell borders all faces).
749 std::vector< IndexRange > faces;
750 hull.getFacetVertices( faces );
751 FaceRange all_faces;
752 for ( Index i = 0; i < faces.size(); i++ )
753 all_faces.push_back( { i, true } );
754 cell_complex.cell_faces.push_back( all_faces );
755 // Vertices of this unique cell will be computed lazily on request.
756 // Taking care of each face.
757 for ( const auto& vtcs : faces )
758 {
759 // every inner face borders cell 0
760 cell_complex.true_face_cell.push_back( 0 );
761 // every outer face borders the infinite cell
762 cell_complex.false_face_cell.push_back( cell_complex.infiniteCell() );
763 }
764 // Taking care of vertices (in consistent order) of every face
765 cell_complex.true_face_vertices.swap( faces );
766 // Taking care of vertex positions
767 hull.getVertexPositions( cell_complex.vertex_position );
768 return true;
769}
770
771
772//-----------------------------------------------------------------------------
773template < int dim, typename TInteger, typename TInternalInteger >
774bool
775DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeDelaunayCellComplex
776( ConvexCellComplex< RealPoint >& cell_complex,
777 const std::vector< RealPoint >& input_points,
778 double precision,
779 bool remove_duplicates )
780{
781 typedef QuickHull< RealDelaunayKernel > Delaunay;
782 typedef typename Delaunay::Ridge Ridge;
783 typedef typename ConvexCellComplex< RealPoint >::FaceRange FaceRange;
784
785 Delaunay del( precision );
786 del.setInput( input_points, remove_duplicates );
787 bool ok = del.computeConvexHull( Delaunay::Status::VerticesCompleted );
788 cell_complex.clear();
789 if ( ! ok ) return false;
790
791 // Build complex, as many maximal cells as convex hull facets.
792 // convex hull facet -> cell of complex
793 // convex hull ridge -> face of complex
794 // (1) Get cell vertices, count ridges/faces and compute their vertices
795 std::map< Ridge, Index > r2f;
796 computeFacetAndRidgeVertices( del,
797 cell_complex.cell_vertices,
798 r2f,
799 cell_complex.true_face_vertices );
800 // (2) assign ridges/faces to cell and conversely
801 const Index nb_r = r2f.size();
802 cell_complex.true_face_cell .resize( nb_r, cell_complex.infiniteCell() );
803 cell_complex.false_face_cell.resize( nb_r, cell_complex.infiniteCell() );
804 cell_complex.true_face_vertices.resize( nb_r );
805 for ( Index cur_f = 0; cur_f < del.nbFiniteFacets(); ++cur_f ) {
806 const auto& facet = del.facets[ cur_f ];
807 FaceRange current_faces;
808 for ( auto neigh_f : facet.neighbors ) {
809 const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
810 const bool pos = cur_f < neigh_f;
811 const Index cur_r = r2f[ r ];
812 cell_complex.true_face_cell [ cur_r ] = r.first;
813 if ( r.second >= del.nbFiniteFacets() )
814 cell_complex.false_face_cell[ cur_r ] = cell_complex.infiniteCell();
815 else
816 cell_complex.false_face_cell[ cur_r ] = r.second;
817 current_faces.emplace_back( cur_r, pos );
818 }
819 cell_complex.cell_faces.push_back( current_faces );
820 }
821 // (3) Takes care of vertex positions
822 del.getVertexPositions( cell_complex.vertex_position );
823 return true;
824}
825
826
827//-----------------------------------------------------------------------------
828template < int dim, typename TInteger, typename TInternalInteger >
829template < typename QHull >
830void
831DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeFacetAndRidgeVertices
832( const QHull& hull,
833 std::vector< IndexRange >& cell_vertices,
834 std::map< typename QHull::Ridge, Index >& r2f,
835 std::vector< IndexRange >& face_vertices )
836{
837 typedef typename QHull::Ridge Ridge;
838
839 ASSERT( hull.status() >= QHull::Status::VerticesCompleted
840 && hull.status() <= QHull::Status::AllCompleted );
841
842 // Get cell vertices and sort them
843 bool ok_fv = hull.getFacetVertices( cell_vertices );
844 if ( ! ok_fv )
845 trace.error() << "[ConvexityHelper::computeFacetAndRidgeVertices]"
846 << " method hull.getFacetVertices failed."
847 << " Maybe QuickHull was not computed till VerticesCompleted."
848 << std::endl;
849 std::vector< IndexRange > sorted_cell_vertices = cell_vertices;
850 for ( auto& vtcs : sorted_cell_vertices )
851 std::sort( vtcs.begin(), vtcs.end() );
852 cell_vertices.resize( hull.nbFiniteFacets() );
853
854 // Count ridges/faces and compute their vertices
855 Index nb_r = 0;
856 face_vertices.clear();
857 for ( Index cur_f = 0; cur_f < hull.nbFiniteFacets(); ++cur_f ) {
858 const auto& facet = hull.facets[ cur_f ];
859 for ( auto neigh_f : facet.neighbors ) {
860 const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
861 auto itr = r2f.find( r );
862 if ( itr == r2f.end() ) {
863 IndexRange result;
864 std::set_intersection( sorted_cell_vertices[ cur_f ].cbegin(),
865 sorted_cell_vertices[ cur_f ].cend(),
866 sorted_cell_vertices[ neigh_f ].cbegin(),
867 sorted_cell_vertices[ neigh_f ].cend(),
868 std::back_inserter( result ) );
869 face_vertices.push_back( result );
870 r2f[ r ] = nb_r++;
871 }
872 }
873 }
874}
875
876// //
877///////////////////////////////////////////////////////////////////////////////
Generated on Sun Nov 27 2022 15:15:41 for DGtal by doxygen 1.9.5