DGtal 2.1.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>::compute3DOpenTriangle
251( const Point& a, const Point& b, const Point& c,
252 bool make_minkowski_summable )
253{
254 using Op = detail::BoundedRationalPolytopeSpecializer< dimension, Integer>;
255 const Vector ab = b - a;
256 const Vector bc = c - b;
257 const Vector n = Op::crossProduct( ab, bc );
258 if ( n == Vector::zero )
259 return LatticePolytope(); // empty polytope
260 auto P = compute3DTriangle( a, b, c, make_minkowski_summable );
261 // Compute domain
262 for ( auto k = 0; k < P.nbHalfSpaces(); k++ )
263 {
264 if ( Op::crossProduct( P.getA( k ), n ) != Vector::zero )
265 P.setStrict( k );
266 }
267 return P;
268}
269
270
271//-----------------------------------------------------------------------------
272template < int dim, typename TInteger, typename TInternalInteger >
273typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
274DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
275computeDegeneratedTriangle
276( const Point& a, const Point& b, const Point& c )
277{
278 if ( a == b ) return computeSegment( a, c );
279 if ( ( a == c ) || ( b == c ) ) return computeSegment( a, b );
280 // The three points are distinct, hence aligned. One is in-between the two others.
281 const Point low = a.inf( b ).inf( c );
282 const Point high = a.sup( b ).sup( c );
283 for ( Dimension k = 0; k < dim; k++ )
284 {
285 const auto lk = low [ k ];
286 const auto hk = high[ k ];
287 if ( ( a[ k ] != lk ) && ( a[ k ] != hk ) ) return computeSegment( b, c );
288 if ( ( b[ k ] != lk ) && ( b[ k ] != hk ) ) return computeSegment( a, c );
289 if ( ( c[ k ] != lk ) && ( c[ k ] != hk ) ) return computeSegment( a, b );
290 }
291 trace.error() << "[ConvexityHelper::computeSegmentFromDegeneratedTriangle] "
292 << "Should never arrive here." << std::endl;
293 return computeSegment( a, a );
294}
295
296
297//-----------------------------------------------------------------------------
298template < int dim, typename TInteger, typename TInternalInteger >
299typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
300DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeOpenSegment
301( const Point& a, const Point& b )
302{
303
304 if constexpr( dim != 3 ) return LatticePolytope();
305 using Op = detail::BoundedRationalPolytopeSpecializer< dimension, Integer>;
306 typedef typename LatticePolytope::Domain Domain;
307 typedef typename LatticePolytope::HalfSpace PolytopeHalfSpace;
308 // Compute domain
309 const Point low = a.inf( b );
310 const Point high = a.sup( b );
311 const Vector ab = b - a;
312 bool degenerate = ( ab == Vector::zero );
313 if ( degenerate ) return LatticePolytope();
314 // Initialize polytope
315 std::vector< PolytopeHalfSpace > PHS;
316 PHS.reserve( 4*dim );
317 // Compute additionnal constraints on domain boundary to make it open.
318 for ( Dimension k = 0; k < dim; k++ )
319 {
320 const Vector bp = Vector::base( k, 1 );
321 PHS.emplace_back( bp, high[ k ] );
322 const Vector bn = Vector::base( k, -1 );
323 PHS.emplace_back( bn, -low[ k ] );
324 }
325 // Compute additionnal constraints on edges so that the
326 // Minkowski sum with axis-aligned edges is valid.
327 for ( Integer d = -1; d <= 1; d += 2 )
328 for ( Dimension k = 0; k < dim; k++ )
329 {
330 const Vector i = Vector::base( k, d );
331 const Vector e = Op::crossProduct( ab, i );
332 if ( e != Vector::zero )
333 {
334 const Integer e_a = e.dot( a );
335 PHS.emplace_back( e, e_a );
336 }
337 }
338 auto P = LatticePolytope( Domain( low, high ), PHS.cbegin(), PHS.cend(), true, false );
339 // Fix < inequalities
340 for ( auto k = 2 * dim; k < 4 * dim; k++ )
341 {
342 auto V = P.getA( k );
343 if ( V.dot( a ) != V.dot( b ) ) P.setStrict( k );
344 }
345 return P;
346
347}
348
349//-----------------------------------------------------------------------------
350template < int dim, typename TInteger, typename TInternalInteger >
351typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
352DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeSegment
353( const Point& a, const Point& b )
354{
355 if constexpr( dim != 3 ) return LatticePolytope();
356 using Op = detail::BoundedRationalPolytopeSpecializer< dimension, Integer>;
357 typedef typename LatticePolytope::Domain Domain;
358 typedef typename LatticePolytope::HalfSpace PolytopeHalfSpace;
359 // Compute domain
360 const Point low = a.inf( b );
361 const Point high = a.sup( b );
362 const Vector ab = b - a;
363 bool degenerate = ( ab == Vector::zero );
364 // Initialize polytope
365 std::vector< PolytopeHalfSpace > PHS;
366 if ( degenerate )
367 return LatticePolytope( Domain( low, high ), PHS.cbegin(), PHS.cend(), true, false );
368 PHS.reserve( 2*dim );
369 // Compute additionnal constraints on edges so that the
370 // Minkowski sum with axis-aligned edges is valid.
371 for ( Integer d = -1; d <= 1; d += 2 )
372 for ( Dimension k = 0; k < dim; k++ )
373 {
374 const Vector i = Vector::base( k, d );
375 const Vector e = Op::crossProduct( ab, i );
376 if ( e != Vector::zero )
377 {
378 const Integer e_a = e.dot( a );
379 PHS.emplace_back( e, e_a );
380 }
381 }
382 return LatticePolytope( Domain( low, high ), PHS.cbegin(), PHS.cend(), true, false );
383}
384
385
386//-----------------------------------------------------------------------------
387template < int dim, typename TInteger, typename TInternalInteger >
388typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::PointRange
389DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
390computeDegeneratedConvexHullVertices
391( PointRange& input_points )
392{
393 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
394 // input_points is a range of points with no duplicates, but which
395 // seems to be not full dimensional.
396 if ( input_points.size() <= 1 )
397 return input_points;
398 // At least 1-dimensional
399 std::vector< Vector > basis;
400 std::vector< Integer > alpha;
401 basis.push_back( input_points[ 1 ] - input_points[ 0 ] );
402 const auto n0 = basis[ 0 ].norm();
403 alpha.push_back( Integer( 0 ) );
404 alpha.push_back( basis[ 0 ].dot( basis[ 0 ] ) );
405 Index i = 2;
406 while ( i < input_points.size() ) {
407 Vector v = input_points[ i ] - input_points[ 0 ];
408 alpha.push_back( basis[ 0 ].dot( v ) );
409 const auto ni = v.norm();
410 const double alignment =
411 fabs( fabs( NumberTraits< Integer >::castToDouble( alpha.back() ) )
412 - ( n0 * ni ) );
413 if ( alignment > 1e-8 ) break;
414 i++;
415 }
416 if ( i == input_points.size() )
417 { // 1-dimensional
418 Index a = 0, b = 0;
419 for ( i = 1; i < input_points.size(); i++ )
420 {
421 if ( alpha[ i ] < alpha[ a ] ) a = i;
422 if ( alpha[ i ] > alpha[ b ] ) b = i;
423 }
424 PointRange X( 2 );
425 X[ 0 ] = input_points[ a ];
426 X[ 1 ] = input_points[ b ];
427 return X;
428 }
429 // at least 2-dimensional
430 ASSERT( dimension > 1 );
431 if ( dimension == 2 )
432 {
433 std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
434 << " Weird error: found initial full dimensional simplex" << std::endl;
435 return PointRange();
436 }
437 if ( dimension >= 4 )
438 {
439 std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
440 << "Degenerated lattice polytope in nD, n >= 4 is not implemented"
441 << std::endl;
442 return PointRange();
443 }
444 basis.push_back( input_points[ i ] - input_points[ 0 ] );
445 Vector n = detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
446 ::crossProduct( basis[ 0 ], basis[ 1 ] );
447 if ( n == Vector::zero )
448 {
449 std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
450 << "Weird numerical error, u . v != |u| |v| but u x v != 0"
451 << std::endl;
452 return PointRange();
453 }
454 // Now the set of input points should be full dimensional.
455 input_points.push_back( input_points[ 0 ] + n );
456 // Compute convex hull
457 ConvexHull hull;
458 hull.setInput( input_points, false );
459 const auto target = ConvexHull::Status::VerticesCompleted;
460 IndexRange full_splx = { 0, 1, i, input_points.size() - 1 };
461 bool ok_init = hull.setInitialSimplex( full_splx );
462 if ( ! ok_init )
463 {
464 std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
465 << "Weird error in hull.setInitialSimplex" << std::endl;
466 return PointRange();
467 }
468 bool ok_hull = hull.computeConvexHull( target );
469 if ( ! ok_hull )
470 {
471 std::cerr << "[ConvexityHelper::computeDegeneratedConvexHullVertices]"
472 << "Weird error in hull.computeConvexHull" << std::endl;
473 return PointRange();
474 }
475 // Get convex hull vertices and remove top point
476 PointRange X;
477 hull.getVertexPositions( X );
478 const std::size_t nX = X.size();
479 for ( std::size_t j = 0; j < nX; j++ )
480 if ( X[ j ] == input_points.back() )
481 {
482 X[ j ] = X.back();
483 break;
484 }
485 X.resize( nX - 1 );
486 return X;
487}
488
489//-----------------------------------------------------------------------------
490template < int dim, typename TInteger, typename TInternalInteger >
491typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::LatticePolytope
492DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
493computeLatticePolytope
494( const PointRange& input_points,
495 bool remove_duplicates,
496 bool make_minkowski_summable )
497{
498 typedef typename LatticePolytope::Domain Domain;
499 typedef typename LatticePolytope::HalfSpace PolytopeHalfSpace;
500 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
501 typedef typename ConvexHull::HalfSpace ConvexHullHalfSpace;
502 typedef typename ConvexHull::Ridge Ridge;
503 if ( input_points.empty() ) return LatticePolytope();
504 if ( input_points.size() <= ( dimension + 1) )
505 return computeSimplex( input_points, remove_duplicates );
506 // Compute domain
507 Point l = input_points[ 0 ];
508 Point u = input_points[ 0 ];
509 for ( std::size_t i = 1; i < input_points.size(); i++ )
510 {
511 const auto& p = input_points[ i ];
512 l = l.inf( p );
513 u = u.sup( p );
514 }
515 Domain domain( l, u );
516 // Compute convex hull
517 ConvexHull hull;
518 hull.setInput( input_points, remove_duplicates );
519 const auto target = ( make_minkowski_summable && dimension == 3 )
520 ? ConvexHull::Status::VerticesCompleted
521 : ConvexHull::Status::FacetsCompleted;
522 bool ok = hull.computeConvexHull( target );
523 if ( ! ok ) // set of points is not full dimensional
524 return computeDegeneratedLatticePolytope( hull.points );
525 // Initialize polytope
526 std::vector< ConvexHullHalfSpace > HS;
527 std::vector< PolytopeHalfSpace > PHS;
528 hull.getFacetHalfSpaces( HS );
529 PHS.reserve( HS.size() );
530 for ( auto& H : HS ) {
531 Vector N;
532 Integer nu;
533 for ( Dimension i = 0; i < dim; ++i )
534 N[ i ] = IntegerConverter< dimension, Integer >::cast( H.internalNormal()[ i ] );
535 nu = IntegerConverter< dimension, Integer >::cast( H.internalIntercept() );
536 PHS.emplace_back( N, nu );
537 }
538 if ( make_minkowski_summable && dimension >= 4 )
539 trace.warning() << "[ConvexityHelper::computeLatticePolytope]"
540 << " Not implemented starting from dimension 4."
541 << std::endl;
542 if ( make_minkowski_summable && dimension == 3 )
543 {
544 // Compute ridge vertices to add edge constraints.
545 PointRange positions;
546 std::vector< IndexRange > facet_vertices;
547 std::vector< IndexRange > ridge_vertices;
548 std::map< Ridge, Index > ridge2index;
549 hull.getVertexPositions( positions );
550 computeFacetAndRidgeVertices( hull, facet_vertices,
551 ridge2index, ridge_vertices );
552 for ( auto p : ridge2index ) {
553 const auto r = p.first;
554 // Copy by value since PHS may be reallocated during the iteration.
555 const auto U = PHS[ r.first ].N; // normal of facet 1
556 const auto V = PHS[ r.second ].N; // normal of facet 2
557 const auto& S = ridge_vertices[ p.second ]; // vertices along facets 1, 2
558 ASSERT( S.size() == 2 && "Invalid ridge" );
559 const auto& P0 = positions[ S[ 0 ] ];
560 const auto& P1 = positions[ S[ 1 ] ];
561 auto E = P1 - P0; // edge 1, 2
562 const auto UxV =
563 detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
564 ::crossProduct( U, V ); // parallel to E
565 ASSERT( E.dot( UxV ) != 0 && "Invalid E / UxV" );
566 if ( E.dot( UxV ) <= 0 ) E = -E; // force correct orientation
567 const auto E1 =
568 detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
569 ::crossProduct( U, E ); // edge on facet 1
570 const auto E2 =
571 detail::BoundedLatticePolytopeSpecializer< dimension, Integer>
572 ::crossProduct( E, V ); // edge on facet 2
573 ASSERT( E1.dot( U ) == 0 && "Invalid E1 / U" );
574 ASSERT( E1.dot( V ) < 0 && "Invalid E1 / V" );
575 ASSERT( E2.dot( V ) == 0 && "Invalid E1 / V" );
576 ASSERT( E2.dot( U ) < 0 && "Invalid E1 / U" );
577 for ( Dimension k = 0; k < dimension; ++k ) {
578 const auto W = U[ k ] * V - V[ k ] * U;
579 const auto nn1 = W.dot( E1 );
580 const auto nn2 = W.dot( E2 );
581 if ( nn1 > 0 && nn2 > 0 ) {
582 PHS.emplace_back( -W, -W.dot( P0 ) );
583 ASSERT( E1.dot(-W ) < 0 && "Invalid E1 /-W" );
584 ASSERT( E2.dot(-W ) < 0 && "Invalid E2 /-W" );
585 }
586 else if ( nn1 < 0 && nn2 < 0 ) {
587 PHS.emplace_back( W, W.dot( P0 ) );
588 ASSERT( E1.dot( W ) < 0 && "Invalid E1 / W" );
589 ASSERT( E2.dot( W ) < 0 && "Invalid E2 / W" );
590 }
591 }
592 }
593 }
594 return LatticePolytope( domain, PHS.cbegin(), PHS.cend(),
595 make_minkowski_summable && ( dimension <= 3 ), true );
596}
597
598//-----------------------------------------------------------------------------
599template < int dim, typename TInteger, typename TInternalInteger >
600typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::PointRange
601DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::
602computeConvexHullVertices
603( const PointRange& input_points,
604 bool remove_duplicates )
605{
606 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
607 PointRange positions;
608 ConvexHull hull;
609 hull.setInput( input_points, remove_duplicates );
610 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
611 if ( !ok )
612 {
613 PointRange Z( input_points );
614 return computeDegeneratedConvexHullVertices( Z );
615 }
616 hull.getVertexPositions( positions );
617 return positions;
618}
619
620//-----------------------------------------------------------------------------
621template < int dim, typename TInteger, typename TInternalInteger >
622template < typename TSurfaceMesh >
623bool
624DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
625( TSurfaceMesh& mesh,
626 const PointRange& input_points,
627 bool remove_duplicates )
628{
629 typedef TSurfaceMesh SurfaceMesh;
630 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
631 ConvexHull hull;
632 hull.setInput( input_points, remove_duplicates );
633 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
634 if ( !ok ) return false;
635 std::vector< RealPoint > positions;
636 hull.getVertexPositions( positions );
637 std::vector< IndexRange > faces;
638 hull.getFacetVertices( faces );
639 mesh = SurfaceMesh( positions.cbegin(), positions.cend(),
640 faces.cbegin(), faces.cend() );
641 return true;
642}
643
644//-----------------------------------------------------------------------------
645template < int dim, typename TInteger, typename TInternalInteger >
646bool
647DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
648( PolygonalSurface< Point >& polysurf,
649 const PointRange& input_points,
650 bool remove_duplicates )
651{
652 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
653 ConvexHull hull;
654 hull.setInput( input_points, remove_duplicates );
655 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
656 if ( !ok ) return false;
657 PointRange positions;
658 hull.getVertexPositions( positions );
659 std::vector< IndexRange > faces;
660 hull.getFacetVertices( faces );
661 // build polygonal surface
662 polysurf.clear();
663 for ( auto p : positions ) polysurf.addVertex( p );
664 for ( auto f : faces ) polysurf.addPolygonalFace( f );
665 return polysurf.build();
666}
667
668//-----------------------------------------------------------------------------
669template < int dim, typename TInteger, typename TInternalInteger >
670bool
671DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullCellComplex
672( ConvexCellComplex< Point >& cell_complex,
673 const PointRange& input_points,
674 bool remove_duplicates )
675{
676 typedef QuickHull< LatticeConvexHullKernel > ConvexHull;
677 typedef typename ConvexCellComplex< Point >::FaceRange FaceRange;
678 ConvexHull hull;
679 hull.setInput( input_points, remove_duplicates );
680 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
681 cell_complex.clear();
682 if ( ! ok ) return false;
683 // Build complex, only 1 finite cell and as many faces as convex hull facets.
684 // Taking care of faces for each cell (here one cell borders all faces).
685 std::vector< IndexRange > faces;
686 hull.getFacetVertices( faces );
687 FaceRange all_faces;
688 for ( Index i = 0; i < faces.size(); i++ )
689 all_faces.push_back( { i, true } );
690 cell_complex.cell_faces.push_back( all_faces );
691 // Vertices of this unique cell will be computed lazily on request.
692 // Taking care of each face.
693 for ( Index i = 0; i < faces.size(); i++ )
694 {
695 // every inner face borders cell 0
696 cell_complex.true_face_cell.push_back( 0 );
697 // every outer face borders the infinite cell
698 cell_complex.false_face_cell.push_back( cell_complex.infiniteCell() );
699 }
700 // Taking care of vertices (in consistent order) of every face
701 cell_complex.true_face_vertices.swap( faces );
702 // Taking care of vertex positions
703 hull.getVertexPositions( cell_complex.vertex_position );
704 return true;
705}
706
707//-----------------------------------------------------------------------------
708template < int dim, typename TInteger, typename TInternalInteger >
709bool
710DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeDelaunayCellComplex
711( ConvexCellComplex< Point >& cell_complex,
712 const PointRange& input_points,
713 bool remove_duplicates )
714{
715 typedef QuickHull< LatticeDelaunayKernel > Delaunay;
716 typedef typename Delaunay::Ridge Ridge;
717 typedef typename ConvexCellComplex< Point >::FaceRange FaceRange;
718
719 Delaunay del;
720 del.setInput( input_points, remove_duplicates );
721 bool ok = del.computeConvexHull( Delaunay::Status::VerticesCompleted );
722 cell_complex.clear();
723 if ( ! ok ) return false;
724
725 // Build complex, as many maximal cells as convex hull facets.
726 // convex hull facet -> cell of complex
727 // convex hull ridge -> face of complex
728 // (1) Get cell vertices, count ridges/faces and compute their vertices
729 std::map< Ridge, Index > r2f;
730 computeFacetAndRidgeVertices( del,
731 cell_complex.cell_vertices,
732 r2f,
733 cell_complex.true_face_vertices );
734 // (2) assign ridges/faces to cell and conversely
735 const Index nb_r = r2f.size();
736 cell_complex.true_face_cell .resize( nb_r, cell_complex.infiniteCell() );
737 cell_complex.false_face_cell.resize( nb_r, cell_complex.infiniteCell() );
738 cell_complex.true_face_vertices.resize( nb_r );
739 for ( Index cur_f = 0; cur_f < del.nbFiniteFacets(); ++cur_f ) {
740 const auto& facet = del.facets[ cur_f ];
741 FaceRange current_faces;
742 for ( auto neigh_f : facet.neighbors ) {
743 const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
744 const bool pos = cur_f < neigh_f;
745 const Index cur_r = r2f[ r ];
746 cell_complex.true_face_cell [ cur_r ] = r.first;
747 if ( r.second >= del.nbFiniteFacets() )
748 cell_complex.false_face_cell[ cur_r ] = cell_complex.infiniteCell();
749 else
750 cell_complex.false_face_cell[ cur_r ] = r.second;
751 current_faces.emplace_back( cur_r, pos );
752 }
753 cell_complex.cell_faces.push_back( current_faces );
754 }
755 // (3) Takes care of vertex positions
756 del.getVertexPositions( cell_complex.vertex_position );
757 return true;
758}
759
760
761//-----------------------------------------------------------------------------
762template < int dim, typename TInteger, typename TInternalInteger >
763typename DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::RationalPolytope
764DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeRationalPolytope
765( const std::vector< RealPoint >& input_points,
766 Integer denominator,
767 bool remove_duplicates,
768 bool make_minkowski_summable )
769{
770 if ( denominator < 1 )
771 trace.error() << "Invalid denominator " << denominator
772 << ". Should be greater or equal to 1." << std::endl;
773 typedef typename RationalPolytope::Domain Domain;
774 typedef typename RationalPolytope::HalfSpace PolytopeHalfSpace;
775 typedef QuickHull< RealConvexHullKernel > ConvexHull;
776 typedef typename ConvexHull::HalfSpace ConvexHullHalfSpace;
777 typedef typename ConvexHull::Ridge Ridge;
778 if ( input_points.empty() ) return RationalPolytope();
779 // Compute convex hull
780 ConvexHull hull( denominator );
781 hull.setInput( input_points, remove_duplicates );
782 const auto target = ( make_minkowski_summable && dimension == 3 )
783 ? ConvexHull::Status::VerticesCompleted
784 : ConvexHull::Status::FacetsCompleted;
785 bool ok = hull.computeConvexHull( target );
786 if ( ! ok ) return RationalPolytope();
787 // Compute domain (as a lattice domain)
788 auto l = hull.points[ 0 ];
789 auto u = hull.points[ 0 ];
790 for ( const auto& p : hull.points ) {
791 l = l.inf( p );
792 u = u.sup( p );
793 }
794 Domain domain( l, u );
795 trace.info() << "Domain l=" << l << " u=" << u << std::endl;
796 // Initialize polytope
797 std::vector< ConvexHullHalfSpace > HS;
798 std::vector< PolytopeHalfSpace > PHS;
799 hull.getFacetHalfSpaces( HS );
800 PHS.reserve( HS.size() );
801 for ( auto& H : HS ) {
802 Vector N;
803 Integer nu;
804 for ( Dimension i = 0; i < dim; ++i )
805 N[ i ] = (Integer) H.internalNormal()[ i ];
806 nu = (Integer) H.internalIntercept();
807 PHS.emplace_back( N, nu );
808 }
809 if ( make_minkowski_summable && dimension >= 4 )
810 trace.warning() << "[ConvexityHelper::computeRationalPolytope]"
811 << " Not implemented starting from dimension 4."
812 << std::endl;
813 if ( make_minkowski_summable && dimension == 3 )
814 {
815 // Compute ridge vertices to add edge constraints.
816 PointRange positions;
817 std::vector< IndexRange > facet_vertices;
818 std::vector< IndexRange > ridge_vertices;
819 std::map< Ridge, Index > ridge2index;
820 hull.getVertexPositions( positions );
821 computeFacetAndRidgeVertices( hull, facet_vertices,
822 ridge2index, ridge_vertices );
823 for ( auto p : ridge2index ) {
824 const auto r = p.first;
825 // Copy by value since PHS may be reallocated during the iteration.
826 const auto U = PHS[ r.first ].N; // normal of facet 1
827 const auto V = PHS[ r.second ].N; // normal of facet 2
828 const auto& S = ridge_vertices[ p.second ]; // vertices along facets 1, 2
829 ASSERT( S.size() == 2 && "Invalid ridge" );
830 const auto& P0 = positions[ S[ 0 ] ];
831 const auto& P1 = positions[ S[ 1 ] ];
832 auto E = P1 - P0; // edge 1, 2
833 const auto UxV =
834 detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
835 ::crossProduct( U, V ); // parallel to E
836 ASSERT( E.dot( UxV ) != 0 && "Invalid E / UxV" );
837 if ( E.dot( UxV ) <= 0 ) E = -E; // force correct orientation
838 const auto E1 =
839 detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
840 ::crossProduct( U, E ); // edge on facet 1
841 const auto E2 =
842 detail::BoundedRationalPolytopeSpecializer< dimension, Integer>
843 ::crossProduct( E, V ); // edge on facet 2
844 ASSERT( E1.dot( U ) == 0 && "Invalid E1 / U" );
845 ASSERT( E1.dot( V ) < 0 && "Invalid E1 / V" );
846 ASSERT( E2.dot( V ) == 0 && "Invalid E1 / V" );
847 ASSERT( E2.dot( U ) < 0 && "Invalid E1 / U" );
848 for ( Dimension k = 0; k < dimension; ++k ) {
849 const auto W = U[ k ] * V - V[ k ] * U;
850 const auto nn1 = W.dot( E1 );
851 const auto nn2 = W.dot( E2 );
852 if ( nn1 > 0 && nn2 > 0 ) {
853 PHS.emplace_back( -W, -W.dot( P0 ) );
854 ASSERT( E1.dot(-W ) < 0 && "Invalid E1 /-W" );
855 ASSERT( E2.dot(-W ) < 0 && "Invalid E2 /-W" );
856 }
857 else if ( nn1 < 0 && nn2 < 0 ) {
858 PHS.emplace_back( W, W.dot( P0 ) );
859 ASSERT( E1.dot( W ) < 0 && "Invalid E1 / W" );
860 ASSERT( E2.dot( W ) < 0 && "Invalid E2 / W" );
861 }
862 }
863 }
864 }
865 return RationalPolytope( denominator, domain, PHS.cbegin(), PHS.cend(),
866 make_minkowski_summable && ( dimension <= 3 ), true );
867}
868
869
870//-----------------------------------------------------------------------------
871template < int dim, typename TInteger, typename TInternalInteger >
872template < typename TSurfaceMesh >
873bool
874DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
875( TSurfaceMesh& mesh,
876 const std::vector< RealPoint >& input_points,
877 double precision,
878 bool remove_duplicates )
879{
880 typedef TSurfaceMesh SurfaceMesh;
881 typedef QuickHull< RealConvexHullKernel > ConvexHull;
882 ConvexHull hull( precision );
883 hull.setInput( input_points, remove_duplicates );
884 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
885 if ( !ok ) return false;
886 std::vector< RealPoint > positions;
887 hull.getVertexPositions( positions );
888 std::vector< IndexRange > faces;
889 hull.getFacetVertices( faces );
890 mesh = SurfaceMesh( positions.cbegin(), positions.cend(),
891 faces.cbegin(), faces.cend() );
892 return true;
893}
894
895
896//-----------------------------------------------------------------------------
897template < int dim, typename TInteger, typename TInternalInteger >
898bool
899DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullBoundary
900( PolygonalSurface< RealPoint >& polysurf,
901 const std::vector< RealPoint >& input_points,
902 double precision,
903 bool remove_duplicates )
904{
905 typedef QuickHull< RealConvexHullKernel > ConvexHull;
906 ConvexHull hull( precision );
907 hull.setInput( input_points, remove_duplicates );
908 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
909 if ( !ok ) return false;
910 std::vector< RealPoint > positions;
911 hull.getVertexPositions( positions );
912 std::vector< IndexRange > faces;
913 hull.getFacetVertices( faces );
914 // build polygonal surface
915 polysurf.clear();
916 for ( auto p : positions ) polysurf.addVertex( p );
917 for ( auto f : faces ) polysurf.addPolygonalFace( f );
918 return polysurf.build();
919}
920
921//-----------------------------------------------------------------------------
922template < int dim, typename TInteger, typename TInternalInteger >
923bool
924DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeConvexHullCellComplex
925( ConvexCellComplex< RealPoint >& cell_complex,
926 const std::vector< RealPoint >& input_points,
927 double precision,
928 bool remove_duplicates )
929{
930 typedef QuickHull< RealConvexHullKernel > ConvexHull;
931 typedef typename ConvexCellComplex< RealPoint >::FaceRange FaceRange;
932 ConvexHull hull( precision );
933 hull.setInput( input_points, remove_duplicates );
934 bool ok = hull.computeConvexHull( ConvexHull::Status::VerticesCompleted );
935 cell_complex.clear();
936 if ( ! ok ) return false;
937 // Build complex, only 1 finite cell and as many faces as convex hull facets.
938 // Taking care of faces for each cell (here one cell borders all faces).
939 std::vector< IndexRange > faces;
940 hull.getFacetVertices( faces );
941 FaceRange all_faces;
942 for ( Index i = 0; i < faces.size(); i++ )
943 all_faces.push_back( { i, true } );
944 cell_complex.cell_faces.push_back( all_faces );
945 // Vertices of this unique cell will be computed lazily on request.
946 // Taking care of each face.
947 for ( Index i = 0; i < faces.size(); i++ )
948 {
949 // every inner face borders cell 0
950 cell_complex.true_face_cell.push_back( 0 );
951 // every outer face borders the infinite cell
952 cell_complex.false_face_cell.push_back( cell_complex.infiniteCell() );
953 }
954 // Taking care of vertices (in consistent order) of every face
955 cell_complex.true_face_vertices.swap( faces );
956 // Taking care of vertex positions
957 hull.getVertexPositions( cell_complex.vertex_position );
958 return true;
959}
960
961
962//-----------------------------------------------------------------------------
963template < int dim, typename TInteger, typename TInternalInteger >
964bool
965DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeDelaunayCellComplex
966( ConvexCellComplex< RealPoint >& cell_complex,
967 const std::vector< RealPoint >& input_points,
968 double precision,
969 bool remove_duplicates )
970{
971 typedef QuickHull< RealDelaunayKernel > Delaunay;
972 typedef typename Delaunay::Ridge Ridge;
973 typedef typename ConvexCellComplex< RealPoint >::FaceRange FaceRange;
974
975 Delaunay del( precision );
976 del.setInput( input_points, remove_duplicates );
977 bool ok = del.computeConvexHull( Delaunay::Status::VerticesCompleted );
978 cell_complex.clear();
979 if ( ! ok ) return false;
980
981 // Build complex, as many maximal cells as convex hull facets.
982 // convex hull facet -> cell of complex
983 // convex hull ridge -> face of complex
984 // (1) Get cell vertices, count ridges/faces and compute their vertices
985 std::map< Ridge, Index > r2f;
986 computeFacetAndRidgeVertices( del,
987 cell_complex.cell_vertices,
988 r2f,
989 cell_complex.true_face_vertices );
990 // (2) assign ridges/faces to cell and conversely
991 const Index nb_r = r2f.size();
992 cell_complex.true_face_cell .resize( nb_r, cell_complex.infiniteCell() );
993 cell_complex.false_face_cell.resize( nb_r, cell_complex.infiniteCell() );
994 cell_complex.true_face_vertices.resize( nb_r );
995 for ( Index cur_f = 0; cur_f < del.nbFiniteFacets(); ++cur_f ) {
996 const auto& facet = del.facets[ cur_f ];
997 FaceRange current_faces;
998 for ( auto neigh_f : facet.neighbors ) {
999 const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
1000 const bool pos = cur_f < neigh_f;
1001 const Index cur_r = r2f[ r ];
1002 cell_complex.true_face_cell [ cur_r ] = r.first;
1003 if ( r.second >= del.nbFiniteFacets() )
1004 cell_complex.false_face_cell[ cur_r ] = cell_complex.infiniteCell();
1005 else
1006 cell_complex.false_face_cell[ cur_r ] = r.second;
1007 current_faces.emplace_back( cur_r, pos );
1008 }
1009 cell_complex.cell_faces.push_back( current_faces );
1010 }
1011 // (3) Takes care of vertex positions
1012 del.getVertexPositions( cell_complex.vertex_position );
1013 return true;
1014}
1015
1016
1017//-----------------------------------------------------------------------------
1018template < int dim, typename TInteger, typename TInternalInteger >
1019template < typename QHull >
1020void
1021DGtal::ConvexityHelper< dim, TInteger, TInternalInteger>::computeFacetAndRidgeVertices
1022( const QHull& hull,
1023 std::vector< IndexRange >& cell_vertices,
1024 std::map< typename QHull::Ridge, Index >& r2f,
1025 std::vector< IndexRange >& face_vertices )
1026{
1027 typedef typename QHull::Ridge Ridge;
1028
1029 ASSERT( hull.status() >= QHull::Status::VerticesCompleted
1030 && hull.status() <= QHull::Status::AllCompleted );
1031
1032 // Get cell vertices and sort them
1033 bool ok_fv = hull.getFacetVertices( cell_vertices );
1034 if ( ! ok_fv )
1035 trace.error() << "[ConvexityHelper::computeFacetAndRidgeVertices]"
1036 << " method hull.getFacetVertices failed."
1037 << " Maybe QuickHull was not computed till VerticesCompleted."
1038 << std::endl;
1039 std::vector< IndexRange > sorted_cell_vertices = cell_vertices;
1040 for ( auto& vtcs : sorted_cell_vertices )
1041 std::sort( vtcs.begin(), vtcs.end() );
1042 cell_vertices.resize( hull.nbFiniteFacets() );
1043
1044 // Count ridges/faces and compute their vertices
1045 Index nb_r = 0;
1046 face_vertices.clear();
1047 for ( Index cur_f = 0; cur_f < hull.nbFiniteFacets(); ++cur_f ) {
1048 const auto& facet = hull.facets[ cur_f ];
1049 for ( auto neigh_f : facet.neighbors ) {
1050 const Ridge r { std::min( cur_f, neigh_f ), std::max( cur_f, neigh_f ) };
1051 auto itr = r2f.find( r );
1052 if ( itr == r2f.end() ) {
1053 IndexRange result;
1054 std::set_intersection( sorted_cell_vertices[ cur_f ].cbegin(),
1055 sorted_cell_vertices[ cur_f ].cend(),
1056 sorted_cell_vertices[ neigh_f ].cbegin(),
1057 sorted_cell_vertices[ neigh_f ].cend(),
1058 std::back_inserter( result ) );
1059 face_vertices.push_back( result );
1060 r2f[ r ] = nb_r++;
1061 }
1062 }
1063 }
1064}
1065
1066// //
1067///////////////////////////////////////////////////////////////////////////////