DGtal  1.4.beta
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
54 namespace 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 {
207  HalfSpace( const InternalVector& aN, const InternalScalar aC )
208  : N( aN ), c( aC ) {}
209  public:
210  HalfSpace() = default;
211  const InternalVector& internalNormal() const { return N; }
212  InternalScalar internalIntercept() const { return c; }
213  };
214
217
229  HalfSpace
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;
262  InternalVector N;
263  for ( Dimension i = 1; i < dimension; i++ )
264  for ( Dimension j = 0; j < dimension; j++ )
265  A.setComponent( i-1, j,
266  Inner::cast( vpoints[ simplex[ i ] ][ j ]
267  - vpoints[ simplex[ 0 ] ][ j ] ) );
268  for ( Dimension j = 0; j < dimension; j++ )
269  N[ j ] = A.cofactor( dimension-1, j );
270  const InternalPoint ip = Inner::cast( vpoints[ simplex[ 0 ] ] );
271  // c = N.dot( vpoints[ simplex[ 0 ] ] );
272  return HalfSpace { N, N.dot( ip ) };
273  }
274
278  {
279  return Outer::cast( H.N );
280  }
281
285  {
286  return Outer::cast( H.c );
287  }
288
297  InternalScalar dot( const HalfSpace& H1, const HalfSpace& H2 ) const
298  {
299  return H1.N.dot( H2.N );
300  }
301
310  bool equal( const HalfSpace& H1, const HalfSpace& H2 ) const
311  {
312  return H1.c == H2.c && H1.N == H2.N;
313  }
314
318  InternalScalar height( const HalfSpace& H, const CoordinatePoint& p ) const
319  { return H.N.dot( Inner::cast( p ) ) - H.c; }
320
324  InternalScalar volume( const HalfSpace& H, const CoordinatePoint& p ) const
325  {
326  InternalScalar v = height( H, p );
327  return v < InternalScalar( 0 ) ? -v : v;
328  }
329
333  bool above( const HalfSpace& H, const CoordinatePoint& p ) const
334  { return height( H, p ) > 0; }
335
339  bool aboveOrOn( const HalfSpace& H, const CoordinatePoint& p ) const
340  { return height( H, p ) >= 0; }
341
345  bool on( const HalfSpace& H, const CoordinatePoint& p ) const
346  { return height( H, p ) == 0; }
347
348
349  }; // template < Dimension dim > struct ConvexHullIntegralKernel {
350
351
352
354  // template class ConvexHullIntegralKernel
372  template < Dimension dim,
373  typename TCoordinateInteger = DGtal::int64_t,
374  typename TInternalInteger = DGtal::int64_t >
376  : public ConvexHullCommonKernel< dim, TCoordinateInteger, TInternalInteger >
377  {
379  // inheriting types
380  // using typename Base::Point;
381  // using typename Base::Vector;
382  // using typename Base::Scalar;
383  using typename Base::CoordinatePoint;
384  using typename Base::CoordinateVector;
385  using typename Base::CoordinateScalar;
386  using typename Base::InternalPoint;
387  using typename Base::InternalVector;
388  using typename Base::InternalScalar;
389  using typename Base::Size;
390  using typename Base::Index;
391  using typename Base::IndexRange;
392  using typename Base::CombinatorialPlaneSimplex;
393  using typename Base::HalfSpace;
394  // inheriting constants
395  using Base::dimension;
396  // inheriting methods
397  using Base::compute;
398  using Base::normal;
399  using Base::intercept;
400  using Base::dot;
401  using Base::equal;
402  using Base::height;
403  using Base::volume;
404  using Base::above;
405  using Base::aboveOrOn;
406  using Base::on;
407
410
414  bool hasInfiniteFacets() const
415  { return false; }
416
421  bool isHalfSpaceFacetInfinite( const HalfSpace& hs ) const
422  {
423  (void) hs; // unused parameter
424  return false;
425  }
426
452  template < typename InputPoint>
453  void makeInput( std::vector< CoordinatePoint >& processed_points,
454  IndexRange& input2comp, IndexRange& comp2input,
455  const std::vector< InputPoint >& input_points,
456  bool remove_duplicates )
457  {
458  const auto F = [&] ( InputPoint input ) -> CoordinatePoint
459  {
460  CoordinatePoint p;
461  for ( Dimension i = 0; i < dimension; i++ )
462  p[ i ] = CoordinateScalar( input[ i ] );
463  return p;
464  };
465  DGtal::detail::transform( processed_points, input2comp, comp2input,
466  input_points.cbegin(), input_points.cend(),
467  F, remove_duplicates );
468  }
469
472  template < typename OutputPoint>
473  void convertPointTo( const CoordinatePoint& p, OutputPoint& out_p ) const
474  {
475  for ( Dimension k = 0; k < dimension; k++ )
476  out_p[ k ] = static_cast<typename OutputPoint::Component>(p[ k ]);
477  }
478
479  }; // template < Dimension dim > struct ConvexHullIntegralKernel {
480
481
483  // template class DelaunayIntegralKernel
504  template < Dimension dim,
505  typename TCoordinateInteger = DGtal::int64_t,
506  typename TInternalInteger = DGtal::int64_t >
508  : public ConvexHullCommonKernel< dim+1, TCoordinateInteger, TInternalInteger >
509  {
511  // inheriting types
512  // using typename Base::Point;
513  // using typename Base::Vector;
514  // using typename Base::Scalar;
515  using typename Base::CoordinatePoint;
516  using typename Base::CoordinateVector;
517  using typename Base::CoordinateScalar;
518  using typename Base::InternalPoint;
519  using typename Base::InternalVector;
520  using typename Base::InternalScalar;
521  using typename Base::Size;
522  using typename Base::Index;
523  using typename Base::IndexRange;
524  using typename Base::CombinatorialPlaneSimplex;
525  using typename Base::HalfSpace;
526  // inheriting constants
527  using Base::dimension;
528  // inheriting methods
529  using Base::compute;
530  using Base::normal;
531  using Base::intercept;
532  using Base::dot;
533  using Base::equal;
534  using Base::height;
535  using Base::volume;
536  using Base::above;
537  using Base::aboveOrOn;
538  using Base::on;
539
542
546  bool hasInfiniteFacets() const
547  { return true; }
548
553  bool isHalfSpaceFacetInfinite( const HalfSpace& hs ) const
554  {
555  return hs.internalNormal()[ dimension - 1 ] >= InternalScalar( 0 );
556  }
557
587  template < typename InputPoint>
588  void makeInput( std::vector< CoordinatePoint >& processed_points,
589  IndexRange& input2comp, IndexRange& comp2input,
590  const std::vector< InputPoint >& input_points,
591  bool remove_duplicates )
592  {
593  const auto F = [&] ( InputPoint input ) -> CoordinatePoint
594  {
595  CoordinatePoint p;
596  CoordinateScalar z = 0;
597  for ( Dimension i = 0; i < dimension-1; i++ ) {
598  const CoordinateScalar x = CoordinateScalar( input[ i ] );
599  p[ i ] = x;
600  z += x*x;
601  }
602  p[ dimension-1 ] = z;
603  return p;
604  };
605  DGtal::detail::transform( processed_points, input2comp, comp2input,
606  input_points.cbegin(), input_points.cend(),
607  F, remove_duplicates );
608  }
609
612  template < typename OutputPoint>
613  void convertPointTo( const CoordinatePoint& p, OutputPoint& out_p ) const
614  {
615  for ( Dimension k = 0; k < dimension-1; k++ )
616  out_p[ k ] = p[ k ];
617  }
618
619  }; // template < Dimension dim > struct DelaunayIntegralKernel {
620
621
623  // template class ConvexHullRationalKernel
654  template < Dimension dim,
655  typename TCoordinateInteger = DGtal::int64_t,
656  typename TInternalInteger = DGtal::int64_t >
658  : public ConvexHullCommonKernel< dim, TCoordinateInteger, TInternalInteger >
659  {
661  // inheriting types
662  // using typename Base::Point;
663  // using typename Base::Vector;
664  // using typename Base::Scalar;
665  using typename Base::CoordinatePoint;
666  using typename Base::CoordinateVector;
667  using typename Base::CoordinateScalar;
668  using typename Base::InternalPoint;
669  using typename Base::InternalVector;
670  using typename Base::InternalScalar;
671  using typename Base::Size;
672  using typename Base::Index;
673  using typename Base::IndexRange;
674  using typename Base::CombinatorialPlaneSimplex;
675  using typename Base::HalfSpace;
676  // inheriting constants
677  using Base::dimension;
678  // inheriting methods
679  using Base::compute;
680  using Base::normal;
681  using Base::intercept;
682  using Base::dot;
683  using Base::equal;
684  using Base::height;
685  using Base::volume;
686  using Base::above;
687  using Base::aboveOrOn;
688  using Base::on;
689
691  double precision;
692
697  ConvexHullRationalKernel( double aPrecision = 1024. )
698  : precision( aPrecision ) {}
699
703  bool hasInfiniteFacets() const
704  { return false; }
705
710  bool isHalfSpaceFacetInfinite( const HalfSpace& hs ) const
711  {
712  (void) hs; // unused parameter
713  return false;
714  }
715
746  template < typename InputPoint>
747  void makeInput( std::vector< CoordinatePoint >& processed_points,
748  IndexRange& input2comp, IndexRange& comp2input,
749  const std::vector< InputPoint >& input_points,
750  bool remove_duplicates )
751  {
752  const auto F = [&] ( InputPoint input ) -> CoordinatePoint
753  {
754  CoordinatePoint p;
755  for ( Dimension i = 0; i < dimension; i++ )
756  p[ i ] = CoordinateScalar( round( input[ i ] * precision ) );
757  return p;
758  };
759  DGtal::detail::transform( processed_points, input2comp, comp2input,
760  input_points.cbegin(), input_points.cend(),
761  F, remove_duplicates );
762  }
763
780  template < typename OutputPoint>
781  void convertPointTo( const CoordinatePoint& p, OutputPoint& out_p ) const
782  {
783  for ( Dimension k = 0; k < dimension; k++ )
784  out_p[ k ] = ( (double) p[ k ] ) / precision;
785  }
786
787
788  }; // template < Dimension dim > struct ConvexHullRationalKernel {
789
790
792  // template class DelaunayRationalKernel
826  template < Dimension dim,
827  typename TCoordinateInteger = DGtal::int64_t,
828  typename TInternalInteger = DGtal::int64_t >
830  : public ConvexHullCommonKernel< dim+1, TCoordinateInteger, TInternalInteger >
831  {
833  // inheriting types
834  // using typename Base::Point;
835  // using typename Base::Vector;
836  // using typename Base::Scalar;
837  using typename Base::CoordinatePoint;
838  using typename Base::CoordinateVector;
839  using typename Base::CoordinateScalar;
840  using typename Base::InternalPoint;
841  using typename Base::InternalVector;
842  using typename Base::InternalScalar;
843  using typename Base::Size;
844  using typename Base::Index;
845  using typename Base::IndexRange;
846  using typename Base::CombinatorialPlaneSimplex;
847  using typename Base::HalfSpace;
848  // inheriting constants
849  using Base::dimension;
850  // inheriting methods
851  using Base::compute;
852  using Base::normal;
853  using Base::intercept;
854  using Base::dot;
855  using Base::equal;
856  using Base::height;
857  using Base::volume;
858  using Base::above;
859  using Base::aboveOrOn;
860  using Base::on;
861
863  double precision;
864
869  DelaunayRationalKernel( double aPrecision = 1024. )
870  : precision( aPrecision ) {}
871
875  bool hasInfiniteFacets() const
876  { return true; }
877
882  bool isHalfSpaceFacetInfinite( const HalfSpace& hs ) const
883  {
884  return hs.internalNormal()[ dimension - 1 ] >= InternalScalar( 0 );
885  }
886
922  template < typename InputPoint>
923  void makeInput( std::vector< CoordinatePoint >& processed_points,
924  IndexRange& input2comp, IndexRange& comp2input,
925  const std::vector< InputPoint >& input_points,
926  bool remove_duplicates )
927  {
928  const auto F = [&] ( InputPoint input ) -> CoordinatePoint
929  {
930  CoordinatePoint p;
931  CoordinateScalar z = 0;
932  for ( Dimension i = 0; i < dimension - 1; i++ ) {
933  const CoordinateScalar x
934  = CoordinateScalar( round( input[ i ] * precision ) );
935  p[ i ] = x;
936  z += x*x;
937  }
938  p[ dimension-1 ] = z;
939  return p;
940  };
941  DGtal::detail::transform( processed_points, input2comp, comp2input,
942  input_points.cbegin(), input_points.cend(),
943  F, remove_duplicates );
944  }
945
965  template < typename OutputPoint>
966  void convertPointTo( const CoordinatePoint& p, OutputPoint& out_p ) const
967  {
968  for ( Dimension k = 0; k < dimension - 1; k++ )
969  out_p[ k ] = ( (double) p[ k ] ) / precision;
970  }
971
972  }; // template < Dimension dim > struct DelaunayRationalKernel {
973
974
975
976 } // namespace DGtal {
977
978 #endif // !defined QuickHullKernels_h
979
980 #undef QuickHullKernels_RECURSES
981 #endif // else defined(QuickHullKernels_RECURSES)
const InternalVector & internalNormal() const
HalfSpace(const InternalVector &aN, const InternalScalar aC)
InternalVector N
the normal vector
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:136
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
unsigned int dim(const Vector &z)