DGtal 1.4.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::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//-----------------------------------------------------------------------------
187template < int dim, typename TInteger, typename TInternalInteger >
188typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
189DGtal::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//-----------------------------------------------------------------------------
248template < int dim, typename TInteger, typename TInternalInteger >
249typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
250DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
251computeDegeneratedTriangle
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//-----------------------------------------------------------------------------
273template < int dim, typename TInteger, typename TInternalInteger >
274typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
275DGtal::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//-----------------------------------------------------------------------------
310template < int dim, typename TInteger, typename TInternalInteger >
311typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::PointRange
312DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
313computeDegeneratedConvexHullVertices
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//-----------------------------------------------------------------------------
413template < int dim, typename TInteger, typename TInternalInteger >
414typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
415DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
416computeLatticePolytope
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//-----------------------------------------------------------------------------
522template < int dim, typename TInteger, typename TInternalInteger >
523typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::PointRange
524DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
525computeConvexHullVertices
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//-----------------------------------------------------------------------------
544template < int dim, typename TInteger, typename TInternalInteger >
545template < typename TSurfaceMesh >
546bool
547DGtal::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//-----------------------------------------------------------------------------
568template < int dim, typename TInteger, typename TInternalInteger >
569bool
570DGtal::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//-----------------------------------------------------------------------------
592template < int dim, typename TInteger, typename TInternalInteger >
593bool
594DGtal::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//-----------------------------------------------------------------------------
631template < int dim, typename TInteger, typename TInternalInteger >
632bool
633DGtal::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//-----------------------------------------------------------------------------
685template < int dim, typename TInteger, typename TInternalInteger >
686typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::RationalPolytope
687DGtal::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//-----------------------------------------------------------------------------
794template < int dim, typename TInteger, typename TInternalInteger >
795template < typename TSurfaceMesh >
796bool
797DGtal::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//-----------------------------------------------------------------------------
820template < int dim, typename TInteger, typename TInternalInteger >
821bool
822DGtal::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//-----------------------------------------------------------------------------
845template < int dim, typename TInteger, typename TInternalInteger >
846bool
847DGtal::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//-----------------------------------------------------------------------------
886template < int dim, typename TInteger, typename TInternalInteger >
887bool
888DGtal::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//-----------------------------------------------------------------------------
941template < int dim, typename TInteger, typename TInternalInteger >
942template < typename QHull >
943void
944DGtal::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///////////////////////////////////////////////////////////////////////////////