DGtal 1.3.0
Loading...
Searching...
No Matches
QuickHullKernels.h
1
17#pragma once
18
31#if defined(QuickHullKernels_RECURSES)
32#error Recursive header files inclusion detected in QuickHullKernels.h
33#else // defined(QuickHullKernels_RECURSES)
35#define QuickHullKernels_RECURSES
36
37#if !defined QuickHullKernels_h
39#define QuickHullKernels_h
40
42// Inclusions
43#include <iostream>
44#include <string>
45#include <vector>
46#include <array>
47#include "DGtal/base/Common.h"
48#include "DGtal/kernel/CInteger.h"
49#include "DGtal/kernel/NumberTraits.h"
50#include "DGtal/kernel/PointVector.h"
51#include "DGtal/kernel/IntegerConverter.h"
52#include "DGtal/math/linalg/SimpleMatrix.h"
53
54namespace DGtal
55{
56 namespace detail {
57
58 // ------------------------ POINT RELATED SERVICES -----------------------
59
65 template < typename Point >
66 Point center( const std::vector< Point >& points )
67 {
68 if ( points.empty() ) return Point::zero;
69 Point l = points[ 0 ];
70 Point u = l;
71 for ( const auto& p : points ) {
72 l = l.inf( p );
73 u = u.sup( p );
74 }
75 return Point( ( l + u ) / 2 );
76 }
77
106 template < typename OutputValue,
107 typename ForwardIterator,
108 typename ConversionFct >
109 void transform( std::vector< OutputValue >& output_values,
110 std::vector< std::size_t >& input2output,
111 std::vector< std::size_t >& output2input,
112 ForwardIterator itb, ForwardIterator ite,
113 const ConversionFct& F,
114 bool remove_duplicates )
115 {
116 typedef std::size_t Size;
117 std::vector< OutputValue > input;
118 while ( itb != ite ) {
119 const auto ip = *itb++;
120 input.push_back( F( ip ) );
121 }
122 if ( ! remove_duplicates ) {
123 output_values.swap( input );
124 input2output.resize( input.size() );
125 output2input.resize( input.size() );
126 for ( Size i = 0; i < input.size(); ++i )
127 input2output[ i ] = output2input[ i ] = i;
128 }
129 else {
130 output_values.clear();
131 std::vector< std::size_t > i2c_sort( input.size() );
132 input2output.resize( input.size() );
133 for ( Size i = 0; i < input.size(); i++ ) i2c_sort[ i ] = i;
134 // indirect sort
135 std::sort( i2c_sort.begin(), i2c_sort.end(),
136 [&input] ( Size i, Size j ) { return input[ i ] < input[ j ]; } );
137 output_values.resize( input.size() );
138 output_values[ 0 ] = input[ i2c_sort[ 0 ] ];
139 input2output[ i2c_sort[ 0 ] ] = 0;
140 Size j = 0;
141 for ( Size i = 1; i < input.size(); i++ ) {
142 if ( input[ i2c_sort[ i-1 ] ] != input[ i2c_sort[ i ] ] )
143 output_values[ ++j ] = input[ i2c_sort[ i ] ];
144 input2output[ i2c_sort[ i ] ] = j;
145 }
146 output_values.resize( j+1 );
147 output2input.resize( output_values.size() );
148 for ( Size i = 0; i < input2output.size(); i++ )
149 output2input[ input2output[ i ] ] = i;
150 }
151 }
152
153 } // namespace detail {
154
155
157 // template class ConvexHullCommonKernel
175 template < Dimension dim,
176 typename TCoordinateInteger = DGtal::int64_t,
177 typename TInternalInteger = DGtal::int64_t >
181 typedef TCoordinateInteger CoordinateInteger;
182 typedef TInternalInteger InternalInteger;
183 //typedef CoordinateInteger Scalar;
186 //typedef DGtal::PointVector< dim, CoordinateInteger > Point;
187 //typedef DGtal::PointVector< dim, CoordinateInteger > Vector;
192 typedef std::size_t Size;
193 typedef Size Index;
194 typedef std::vector< Index > IndexRange;
195 typedef std::array< Index, dim > CombinatorialPlaneSimplex;
196 static const Dimension dimension = dim;
197
202
203 class HalfSpace {
208 : N( aN ), c( aC ) {}
209 public:
210 HalfSpace() = default;
211 const InternalVector& internalNormal() const { return N; }
213 };
214
217
230 compute( const std::vector< CoordinatePoint >& vpoints,
231 const CombinatorialPlaneSimplex& simplex,
232 Index idx_below )
233 {
234 HalfSpace hs = compute( vpoints, simplex );
235 if ( hs.N != InternalVector::zero )
236 {
237 const InternalPoint ip = Inner::cast( vpoints[ idx_below ] );
238 const InternalScalar nu = hs.N.dot( ip );
239 //const Scalar nu = hs.N.dot( vpoints[ idx_below ] );
240 if ( nu > hs.c ) { hs.N = -hs.N; hs.c = -hs.c; }
241 }
242 return hs;
243 }
244
256 HalfSpace
257 compute( const std::vector< CoordinatePoint >& vpoints,
258 const CombinatorialPlaneSimplex& simplex )
259 {
261 Matrix A;
264 for ( Dimension i = 1; i < dimension; i++ )
265 for ( Dimension j = 0; j < dimension; j++ )
266 A.setComponent( i-1, j,
267 Inner::cast( vpoints[ simplex[ i ] ][ j ]
268 - vpoints[ simplex[ 0 ] ][ j ] ) );
269 for ( Dimension j = 0; j < dimension; j++ )
270 N[ j ] = A.cofactor( dimension-1, j );
271 const InternalPoint ip = Inner::cast( vpoints[ simplex[ 0 ] ] );
272 // c = N.dot( vpoints[ simplex[ 0 ] ] );
273 return HalfSpace { N, N.dot( ip ) };
274 }
275
279 {
280 return Outer::cast( H.N );
281 }
282
286 {
287 return Outer::cast( H.c );
288 }
289
298 InternalScalar dot( const HalfSpace& H1, const HalfSpace& H2 ) const
299 {
300 return H1.N.dot( H2.N );
301 }
302
311 bool equal( const HalfSpace& H1, const HalfSpace& H2 ) const
312 {
313 return H1.c == H2.c && H1.N == H2.N;
314 }
315
320 { return H.N.dot( Inner::cast( p ) ) - H.c; }
321
326 {
327 InternalScalar v = height( H, p );
328 return v < InternalScalar( 0 ) ? -v : v;
329 }
330
334 bool above( const HalfSpace& H, const CoordinatePoint& p ) const
335 { return height( H, p ) > 0; }
336
340 bool aboveOrOn( const HalfSpace& H, const CoordinatePoint& p ) const
341 { return height( H, p ) >= 0; }
342
346 bool on( const HalfSpace& H, const CoordinatePoint& p ) const
347 { return height( H, p ) == 0; }
348
349
350 }; // template < Dimension dim > struct ConvexHullIntegralKernel {
351
352
353
355 // template class ConvexHullIntegralKernel
373 template < Dimension dim,
374 typename TCoordinateInteger = DGtal::int64_t,
375 typename TInternalInteger = DGtal::int64_t >
377 : public ConvexHullCommonKernel< dim, TCoordinateInteger, TInternalInteger >
378 {
380 // inheriting types
381 // using typename Base::Point;
382 // using typename Base::Vector;
383 // using typename Base::Scalar;
384 using typename Base::CoordinatePoint;
385 using typename Base::CoordinateVector;
386 using typename Base::CoordinateScalar;
387 using typename Base::InternalPoint;
388 using typename Base::InternalVector;
389 using typename Base::InternalScalar;
390 using typename Base::Size;
391 using typename Base::Index;
392 using typename Base::IndexRange;
393 using typename Base::CombinatorialPlaneSimplex;
394 using typename Base::HalfSpace;
395 // inheriting constants
396 using Base::dimension;
397 // inheriting methods
398 using Base::compute;
399 using Base::normal;
400 using Base::intercept;
401 using Base::dot;
402 using Base::equal;
403 using Base::height;
404 using Base::volume;
405 using Base::above;
406 using Base::aboveOrOn;
407 using Base::on;
408
411
415 bool hasInfiniteFacets() const
416 { return false; }
417
422 bool isHalfSpaceFacetInfinite( const HalfSpace& hs ) const
423 {
424 (void) hs; // unused parameter
425 return false;
426 }
427
453 template < typename InputPoint>
454 void makeInput( std::vector< CoordinatePoint >& processed_points,
455 IndexRange& input2comp, IndexRange& comp2input,
456 const std::vector< InputPoint >& input_points,
457 bool remove_duplicates )
458 {
459 const auto F = [&] ( InputPoint input ) -> CoordinatePoint
460 {
462 for ( Dimension i = 0; i < dimension; i++ )
463 p[ i ] = CoordinateScalar( input[ i ] );
464 return p;
465 };
466 DGtal::detail::transform( processed_points, input2comp, comp2input,
467 input_points.cbegin(), input_points.cend(),
468 F, remove_duplicates );
469 }
470
473 template < typename OutputPoint>
474 void convertPointTo( const CoordinatePoint& p, OutputPoint& out_p ) const
475 {
476 for ( Dimension k = 0; k < dimension; k++ )
477 out_p[ k ] = p[ k ];
478 }
479
480 }; // template < Dimension dim > struct ConvexHullIntegralKernel {
481
482
484 // template class DelaunayIntegralKernel
505 template < Dimension dim,
506 typename TCoordinateInteger = DGtal::int64_t,
507 typename TInternalInteger = DGtal::int64_t >
509 : public ConvexHullCommonKernel< dim+1, TCoordinateInteger, TInternalInteger >
510 {
512 // inheriting types
513 // using typename Base::Point;
514 // using typename Base::Vector;
515 // using typename Base::Scalar;
516 using typename Base::CoordinatePoint;
517 using typename Base::CoordinateVector;
518 using typename Base::CoordinateScalar;
519 using typename Base::InternalPoint;
520 using typename Base::InternalVector;
521 using typename Base::InternalScalar;
522 using typename Base::Size;
523 using typename Base::Index;
524 using typename Base::IndexRange;
525 using typename Base::CombinatorialPlaneSimplex;
526 using typename Base::HalfSpace;
527 // inheriting constants
528 using Base::dimension;
529 // inheriting methods
530 using Base::compute;
531 using Base::normal;
532 using Base::intercept;
533 using Base::dot;
534 using Base::equal;
535 using Base::height;
536 using Base::volume;
537 using Base::above;
538 using Base::aboveOrOn;
539 using Base::on;
540
543
547 bool hasInfiniteFacets() const
548 { return true; }
549
554 bool isHalfSpaceFacetInfinite( const HalfSpace& hs ) const
555 {
556 return hs.internalNormal()[ dimension - 1 ] >= InternalScalar( 0 );
557 }
558
588 template < typename InputPoint>
589 void makeInput( std::vector< CoordinatePoint >& processed_points,
590 IndexRange& input2comp, IndexRange& comp2input,
591 const std::vector< InputPoint >& input_points,
592 bool remove_duplicates )
593 {
594 const auto F = [&] ( InputPoint input ) -> CoordinatePoint
595 {
597 CoordinateScalar z = 0;
598 for ( Dimension i = 0; i < dimension-1; i++ ) {
599 const CoordinateScalar x = CoordinateScalar( input[ i ] );
600 p[ i ] = x;
601 z += x*x;
602 }
603 p[ dimension-1 ] = z;
604 return p;
605 };
606 DGtal::detail::transform( processed_points, input2comp, comp2input,
607 input_points.cbegin(), input_points.cend(),
608 F, remove_duplicates );
609 }
610
613 template < typename OutputPoint>
614 void convertPointTo( const CoordinatePoint& p, OutputPoint& out_p ) const
615 {
616 for ( Dimension k = 0; k < dimension-1; k++ )
617 out_p[ k ] = p[ k ];
618 }
619
620 }; // template < Dimension dim > struct DelaunayIntegralKernel {
621
622
624 // template class ConvexHullRationalKernel
655 template < Dimension dim,
656 typename TCoordinateInteger = DGtal::int64_t,
657 typename TInternalInteger = DGtal::int64_t >
659 : public ConvexHullCommonKernel< dim, TCoordinateInteger, TInternalInteger >
660 {
662 // inheriting types
663 // using typename Base::Point;
664 // using typename Base::Vector;
665 // using typename Base::Scalar;
666 using typename Base::CoordinatePoint;
667 using typename Base::CoordinateVector;
668 using typename Base::CoordinateScalar;
669 using typename Base::InternalPoint;
670 using typename Base::InternalVector;
671 using typename Base::InternalScalar;
672 using typename Base::Size;
673 using typename Base::Index;
674 using typename Base::IndexRange;
675 using typename Base::CombinatorialPlaneSimplex;
676 using typename Base::HalfSpace;
677 // inheriting constants
678 using Base::dimension;
679 // inheriting methods
680 using Base::compute;
681 using Base::normal;
682 using Base::intercept;
683 using Base::dot;
684 using Base::equal;
685 using Base::height;
686 using Base::volume;
687 using Base::above;
688 using Base::aboveOrOn;
689 using Base::on;
690
692 double precision;
693
698 ConvexHullRationalKernel( double aPrecision = 1024. )
699 : precision( aPrecision ) {}
700
704 bool hasInfiniteFacets() const
705 { return false; }
706
711 bool isHalfSpaceFacetInfinite( const HalfSpace& hs ) const
712 {
713 (void) hs; // unused parameter
714 return false;
715 }
716
747 template < typename InputPoint>
748 void makeInput( std::vector< CoordinatePoint >& processed_points,
749 IndexRange& input2comp, IndexRange& comp2input,
750 const std::vector< InputPoint >& input_points,
751 bool remove_duplicates )
752 {
753 const auto F = [&] ( InputPoint input ) -> CoordinatePoint
754 {
756 for ( Dimension i = 0; i < dimension; i++ )
757 p[ i ] = CoordinateScalar( round( input[ i ] * precision ) );
758 return p;
759 };
760 DGtal::detail::transform( processed_points, input2comp, comp2input,
761 input_points.cbegin(), input_points.cend(),
762 F, remove_duplicates );
763 }
764
781 template < typename OutputPoint>
782 void convertPointTo( const CoordinatePoint& p, OutputPoint& out_p ) const
783 {
784 for ( Dimension k = 0; k < dimension; k++ )
785 out_p[ k ] = ( (double) p[ k ] ) / precision;
786 }
787
788
789 }; // template < Dimension dim > struct ConvexHullRationalKernel {
790
791
793 // template class DelaunayRationalKernel
827 template < Dimension dim,
828 typename TCoordinateInteger = DGtal::int64_t,
829 typename TInternalInteger = DGtal::int64_t >
831 : public ConvexHullCommonKernel< dim+1, TCoordinateInteger, TInternalInteger >
832 {
834 // inheriting types
835 // using typename Base::Point;
836 // using typename Base::Vector;
837 // using typename Base::Scalar;
838 using typename Base::CoordinatePoint;
839 using typename Base::CoordinateVector;
840 using typename Base::CoordinateScalar;
841 using typename Base::InternalPoint;
842 using typename Base::InternalVector;
843 using typename Base::InternalScalar;
844 using typename Base::Size;
845 using typename Base::Index;
846 using typename Base::IndexRange;
847 using typename Base::CombinatorialPlaneSimplex;
848 using typename Base::HalfSpace;
849 // inheriting constants
850 using Base::dimension;
851 // inheriting methods
852 using Base::compute;
853 using Base::normal;
854 using Base::intercept;
855 using Base::dot;
856 using Base::equal;
857 using Base::height;
858 using Base::volume;
859 using Base::above;
860 using Base::aboveOrOn;
861 using Base::on;
862
864 double precision;
865
870 DelaunayRationalKernel( double aPrecision = 1024. )
871 : precision( aPrecision ) {}
872
876 bool hasInfiniteFacets() const
877 { return true; }
878
883 bool isHalfSpaceFacetInfinite( const HalfSpace& hs ) const
884 {
885 return hs.internalNormal()[ dimension - 1 ] >= InternalScalar( 0 );
886 }
887
923 template < typename InputPoint>
924 void makeInput( std::vector< CoordinatePoint >& processed_points,
925 IndexRange& input2comp, IndexRange& comp2input,
926 const std::vector< InputPoint >& input_points,
927 bool remove_duplicates )
928 {
929 const auto F = [&] ( InputPoint input ) -> CoordinatePoint
930 {
932 CoordinateScalar z = 0;
933 for ( Dimension i = 0; i < dimension - 1; i++ ) {
934 const CoordinateScalar x
935 = CoordinateScalar( round( input[ i ] * precision ) );
936 p[ i ] = x;
937 z += x*x;
938 }
939 p[ dimension-1 ] = z;
940 return p;
941 };
942 DGtal::detail::transform( processed_points, input2comp, comp2input,
943 input_points.cbegin(), input_points.cend(),
944 F, remove_duplicates );
945 }
946
966 template < typename OutputPoint>
967 void convertPointTo( const CoordinatePoint& p, OutputPoint& out_p ) const
968 {
969 for ( Dimension k = 0; k < dimension - 1; k++ )
970 out_p[ k ] = ( (double) p[ k ] ) / precision;
971 }
972
973 }; // template < Dimension dim > struct DelaunayRationalKernel {
974
975
976
977} // namespace DGtal {
978
979#endif // !defined QuickHullKernels_h
980
981#undef QuickHullKernels_RECURSES
982#endif // else defined(QuickHullKernels_RECURSES)
HalfSpace(const InternalVector &aN, const InternalScalar aC)
InternalVector N
the normal vector
const InternalVector & internalNormal() const
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:593
auto dot(const PointVector< dim, OtherComponent, OtherStorage > &v) const -> decltype(DGtal::dotProduct(*this, v))
Dot product with a PointVector.
static Self zero
Static const for zero PointVector.
Definition: PointVector.h:1595
Aim: implements basic MxN Matrix services (M,N>=1).
Definition: SimpleMatrix.h:76
void transform(std::vector< OutputValue > &output_values, std::vector< std::size_t > &input2output, std::vector< std::size_t > &output2input, ForwardIterator itb, ForwardIterator ite, const ConversionFct &F, bool remove_duplicates)
Point center(const std::vector< Point > &points)
DGtal is the top-level namespace which contains all DGtal functions and types.
boost::int64_t int64_t
signed 94-bit integer.
Definition: BasicTypes.h:74
DGtal::uint32_t Dimension
Definition: Common.h:137
Aim: the common part of all geometric kernels for computing the convex hull or Delaunay triangulation...
DGtal::PointVector< dim, InternalInteger > InternalVector
BOOST_CONCEPT_ASSERT((concepts::CInteger< TCoordinateInteger >))
std::array< Index, dim > CombinatorialPlaneSimplex
bool above(const HalfSpace &H, const CoordinatePoint &p) const
InternalScalar height(const HalfSpace &H, const CoordinatePoint &p) const
CoordinateInteger CoordinateScalar
HalfSpace compute(const std::vector< CoordinatePoint > &vpoints, const CombinatorialPlaneSimplex &simplex)
TCoordinateInteger CoordinateInteger
CoordinateVector normal(const HalfSpace &H) const
InternalScalar volume(const HalfSpace &H, const CoordinatePoint &p) const
InternalScalar dot(const HalfSpace &H1, const HalfSpace &H2) const
HalfSpace compute(const std::vector< CoordinatePoint > &vpoints, const CombinatorialPlaneSimplex &simplex, Index idx_below)
DGtal::PointVector< dim, CoordinateInteger > CoordinatePoint
std::vector< Index > IndexRange
bool equal(const HalfSpace &H1, const HalfSpace &H2) const
BOOST_CONCEPT_ASSERT((concepts::CInteger< TInternalInteger >))
DGtal::PointVector< dim, InternalInteger > InternalPoint
IntegerConverter< dim, InternalInteger > Inner
Converter to inner internal integers or lattice points / vector.
bool on(const HalfSpace &H, const CoordinatePoint &p) const
static const Dimension dimension
bool aboveOrOn(const HalfSpace &H, const CoordinatePoint &p) const
IntegerConverter< dim, CoordinateInteger > Outer
Converter to outer coordinate integers or lattice points / vector.
ConvexHullCommonKernel()=default
Default constructor.
CoordinateScalar intercept(const HalfSpace &H) const
DGtal::PointVector< dim, CoordinateInteger > CoordinateVector
Aim: a geometric kernel to compute the convex hull of digital points with integer-only arithmetic.
ConvexHullCommonKernel< dim, TCoordinateInteger, TInternalInteger > Base
void convertPointTo(const CoordinatePoint &p, OutputPoint &out_p) const
void makeInput(std::vector< CoordinatePoint > &processed_points, IndexRange &input2comp, IndexRange &comp2input, const std::vector< InputPoint > &input_points, bool remove_duplicates)
bool isHalfSpaceFacetInfinite(const HalfSpace &hs) const
static const Dimension dimension
ConvexHullIntegralKernel()=default
Default constructor.
Aim: a geometric kernel to compute the convex hull of floating points with integer-only arithmetic....
void convertPointTo(const CoordinatePoint &p, OutputPoint &out_p) const
void makeInput(std::vector< CoordinatePoint > &processed_points, IndexRange &input2comp, IndexRange &comp2input, const std::vector< InputPoint > &input_points, bool remove_duplicates)
double precision
The precision as the common denominator for all rational points.
bool isHalfSpaceFacetInfinite(const HalfSpace &hs) const
ConvexHullRationalKernel(double aPrecision=1024.)
static const Dimension dimension
ConvexHullCommonKernel< dim, TCoordinateInteger, TInternalInteger > Base
Aim: a geometric kernel to compute the Delaunay triangulation of digital points with integer-only ari...
CoordinateInteger CoordinateScalar
std::vector< Index > IndexRange
ConvexHullCommonKernel< dim+1, TCoordinateInteger, TInternalInteger > Base
void convertPointTo(const CoordinatePoint &p, OutputPoint &out_p) const
DelaunayIntegralKernel()=default
Default constructor.
void makeInput(std::vector< CoordinatePoint > &processed_points, IndexRange &input2comp, IndexRange &comp2input, const std::vector< InputPoint > &input_points, bool remove_duplicates)
static const Dimension dimension
bool isHalfSpaceFacetInfinite(const HalfSpace &hs) const
Aim: a geometric kernel to compute the Delaunay triangulation of a range of floating points with inte...
CoordinateInteger CoordinateScalar
std::vector< Index > IndexRange
double precision
The precision as the common denominator for all rational points.
ConvexHullCommonKernel< dim+1, TCoordinateInteger, TInternalInteger > Base
DelaunayRationalKernel(double aPrecision=1024.)
static const Dimension dimension
void makeInput(std::vector< CoordinatePoint > &processed_points, IndexRange &input2comp, IndexRange &comp2input, const std::vector< InputPoint > &input_points, bool remove_duplicates)
void convertPointTo(const CoordinatePoint &p, OutputPoint &out_p) const
bool isHalfSpaceFacetInfinite(const HalfSpace &hs) const
----------— INTEGER/POINT CONVERSION SERVICES -----------------—
static Integer cast(Integer i)
Aim: Concept checking for Integer Numbers. More precisely, this concept is a refinement of both CEucl...
Definition: CInteger.h:88
MyPointD Point
Definition: testClone2.cpp:383
HalfEdgeDataStructure::Size Size