QuickHullKernels.h

 
#pragma once
 
#if defined(QuickHullKernels_RECURSES)
#error Recursive header files inclusion detected in QuickHullKernels.h
#else // defined(QuickHullKernels_RECURSES)
#define QuickHullKernels_RECURSES
 
#if !defined QuickHullKernels_h
#define QuickHullKernels_h
 
// Inclusions
#include <iostream>
#include <string>
#include <vector>
#include <array>
#include "DGtal/base/Common.h"
#include "DGtal/kernel/CInteger.h"
#include "DGtal/kernel/NumberTraits.h"
#include "DGtal/kernel/PointVector.h"
#include "DGtal/kernel/IntegerConverter.h"
#include "DGtal/math/linalg/SimpleMatrix.h"
 
namespace DGtal
{
  namespace detail {
    
    // ------------------------ POINT RELATED SERVICES -----------------------
    
    template < typename Point >
    Point center( const std::vector< Point >& points )
    {
      if ( points.empty() ) return Point::zero;
      Point l = points[ 0 ];
      Point u = l;
      for ( const auto& p : points ) {
        l = l.inf( p );
        u = u.sup( p );
      }
      return Point( ( l + u ) / 2 );
    }
 
    template < typename OutputValue,
               typename ForwardIterator,
               typename ConversionFct >
    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 )
    {
      typedef std::size_t Size;
      std::vector< OutputValue > input;
      while ( itb != ite ) {
        const auto ip = *itb++;
        input.push_back( F( ip ) );
      }
      if ( ! remove_duplicates ) {
        output_values.swap( input );
        input2output.resize( input.size() );
        output2input.resize( input.size() );
        for ( Size i = 0; i < input.size(); ++i )
          input2output[ i ] = output2input[ i ] = i;
      }
      else {
        output_values.clear();
        std::vector< std::size_t > i2c_sort( input.size() );
        input2output.resize( input.size() );
        for ( Size i = 0; i < input.size(); i++ ) i2c_sort[ i ] = i;
        // indirect sort
        std::sort( i2c_sort.begin(), i2c_sort.end(),
                   [&input] ( Size i, Size j ) { return input[ i ] < input[ j ]; } );
        output_values.resize( input.size() );
        output_values[ 0 ] = input[ i2c_sort[ 0 ] ];
        input2output[ i2c_sort[ 0 ] ] = 0;
        Size j = 0;
        for ( Size i = 1; i < input.size(); i++ ) {
          if ( input[ i2c_sort[ i-1 ] ] != input[ i2c_sort[ i ] ] )
            output_values[ ++j ] = input[ i2c_sort[ i ] ];
          input2output[ i2c_sort[ i ] ] = j;
        }
        output_values.resize( j+1 );
        output2input.resize( output_values.size() );
        for ( Size i = 0; i < input2output.size(); i++ )
          output2input[ input2output[ i ] ] = i;
      }      
    }
 
  } // namespace detail {
 
 
  // template class ConvexHullCommonKernel
  template < Dimension dim,
             typename TCoordinateInteger  = DGtal::int64_t,
             typename TInternalInteger = DGtal::int64_t >
  struct ConvexHullCommonKernel {
    BOOST_CONCEPT_ASSERT(( concepts::CInteger<TCoordinateInteger> ));
    BOOST_CONCEPT_ASSERT(( concepts::CInteger<TInternalInteger> ));
    typedef TCoordinateInteger         CoordinateInteger;
    typedef TInternalInteger           InternalInteger;
    //typedef CoordinateInteger          Scalar;
    typedef CoordinateInteger          CoordinateScalar;
    typedef InternalInteger            InternalScalar;
    //typedef DGtal::PointVector< dim, CoordinateInteger > Point;
    //typedef DGtal::PointVector< dim, CoordinateInteger > Vector;
    typedef DGtal::PointVector< dim, CoordinateInteger > CoordinatePoint;
    typedef DGtal::PointVector< dim, CoordinateInteger > CoordinateVector;
    typedef DGtal::PointVector< dim, InternalInteger >   InternalPoint;
    typedef DGtal::PointVector< dim, InternalInteger >   InternalVector;
    typedef std::size_t                Size;
    typedef Size                       Index;
    typedef std::vector< Index >       IndexRange;
    typedef std::array< Index, dim >   CombinatorialPlaneSimplex;
    static const Dimension dimension = dim;
 
    typedef IntegerConverter< dim, CoordinateInteger > Outer; 
    typedef IntegerConverter< dim, InternalInteger >   Inner;
    
    class HalfSpace {
      friend class ConvexHullCommonKernel< dim, CoordinateInteger, InternalInteger >;
      InternalVector N; 
      InternalScalar c; 
      HalfSpace( const InternalVector& aN, const InternalScalar aC )
        : N( aN ), c( aC ) {}
    public:
      HalfSpace() = default;
      const InternalVector& internalNormal() const    { return N; }
      InternalScalar internalIntercept() const { return c; }
    };
    
    ConvexHullCommonKernel() = default;
 
    HalfSpace
    compute( const std::vector< CoordinatePoint >& vpoints,
             const CombinatorialPlaneSimplex& simplex,
             Index idx_below )
    {
      HalfSpace hs = compute( vpoints, simplex );
      if ( hs.N != InternalVector::zero ) 
        {
          const InternalPoint  ip = Inner::cast( vpoints[ idx_below ] );
          const InternalScalar nu = hs.N.dot( ip );
          //const Scalar nu = hs.N.dot( vpoints[ idx_below ] );
          if ( nu > hs.c ) { hs.N = -hs.N; hs.c = -hs.c; }
        }
      return hs;
    }
 
    HalfSpace
    compute( const std::vector< CoordinatePoint >& vpoints,
             const CombinatorialPlaneSimplex& simplex )
    {
      typedef DGtal::SimpleMatrix< InternalScalar, dimension, dimension > Matrix;
      Matrix A;
      InternalVector N;
      InternalScalar c;
      for ( Dimension i = 1; i < dimension; i++ )
        for ( Dimension j = 0; j < dimension; j++ )
          A.setComponent( i-1, j,
                          Inner::cast( vpoints[ simplex[ i ] ][ j ]
                                       - vpoints[ simplex[ 0 ] ][ j ] ) );
      for ( Dimension j = 0; j < dimension; j++ )
        N[ j ] = A.cofactor( dimension-1, j );
      const InternalPoint ip = Inner::cast( vpoints[ simplex[ 0 ] ] );
      // c = N.dot( vpoints[ simplex[ 0 ] ] );
      return HalfSpace { N, N.dot( ip ) };
    }
    
    CoordinateVector normal( const HalfSpace& H ) const
    {
      return Outer::cast( H.N );
    }
 
    CoordinateScalar intercept( const HalfSpace& H ) const
    {
      return Outer::cast( H.c );
    }
    
    InternalScalar dot( const HalfSpace& H1, const HalfSpace& H2 ) const
    {
      return H1.N.dot( H2.N );
    }
 
    bool equal( const HalfSpace& H1, const HalfSpace& H2 ) const
    {
      return H1.c == H2.c && H1.N == H2.N;
    }
    
    InternalScalar height( const HalfSpace& H, const CoordinatePoint& p ) const
    { return H.N.dot( Inner::cast( p ) ) - H.c; }
 
    InternalScalar volume( const HalfSpace& H, const CoordinatePoint& p ) const
    {
      InternalScalar v = height( H, p );
      return v < InternalScalar( 0 ) ? -v : v;
    }
 
    bool above( const HalfSpace& H, const CoordinatePoint& p ) const
    { return height( H, p ) > 0; }
 
    bool aboveOrOn( const HalfSpace& H, const CoordinatePoint& p ) const
    { return height( H, p ) >= 0; }
 
    bool on( const HalfSpace& H, const CoordinatePoint& p ) const
    { return height( H, p ) == 0; } 
    
    
  }; //   template < Dimension dim >  struct ConvexHullIntegralKernel {
 
 
  
  // template class ConvexHullIntegralKernel
  template < Dimension dim,
             typename TCoordinateInteger  = DGtal::int64_t,
             typename TInternalInteger = DGtal::int64_t >
  struct ConvexHullIntegralKernel
    : public ConvexHullCommonKernel< dim, TCoordinateInteger, TInternalInteger >
  {
    typedef ConvexHullCommonKernel< dim, TCoordinateInteger, TInternalInteger > Base;
    // inheriting types
    // using typename Base::Point;
    // using typename Base::Vector;
    // using typename Base::Scalar;
    using typename Base::CoordinatePoint;
    using typename Base::CoordinateVector;
    using typename Base::CoordinateScalar;
    using typename Base::InternalPoint;
    using typename Base::InternalVector;
    using typename Base::InternalScalar;
    using typename Base::Size;
    using typename Base::Index;
    using typename Base::IndexRange;
    using typename Base::CombinatorialPlaneSimplex;
    using typename Base::HalfSpace;
    // inheriting constants
    using Base::dimension;
    // inheriting methods
    using Base::compute;
    using Base::normal;
    using Base::intercept;
    using Base::dot;
    using Base::equal;
    using Base::height;
    using Base::volume;
    using Base::above;
    using Base::aboveOrOn;
    using Base::on;
    
    ConvexHullIntegralKernel() = default;
 
    bool hasInfiniteFacets() const
    { return false; }
 
    bool isHalfSpaceFacetInfinite( const HalfSpace& hs ) const
    {
      (void) hs; // unused parameter
      return false;
    }
    
    template < typename InputPoint>
    void makeInput( std::vector< CoordinatePoint >& processed_points,
                    IndexRange& input2comp, IndexRange& comp2input,
                    const std::vector< InputPoint >& input_points,
                    bool remove_duplicates )
    {
      const auto F = [&] ( InputPoint input ) -> CoordinatePoint
      {
        CoordinatePoint p;
        for ( Dimension i = 0; i < dimension; i++ )
          p[ i ] = CoordinateScalar( input[ i ] );
        return p;
      };
      DGtal::detail::transform( processed_points, input2comp, comp2input,
                                input_points.cbegin(), input_points.cend(),
                                F, remove_duplicates );
    }
 
    template < typename OutputPoint>
    void convertPointTo( const CoordinatePoint& p, OutputPoint& out_p ) const
    {
      for ( Dimension k = 0; k < dimension; k++ )
        out_p[ k ] = p[ k ];
    }
 
  }; //   template < Dimension dim >  struct ConvexHullIntegralKernel {
 
 
  // template class DelaunayIntegralKernel
  template < Dimension dim,
             typename TCoordinateInteger  = DGtal::int64_t,
             typename TInternalInteger = DGtal::int64_t >
  struct DelaunayIntegralKernel
    : public ConvexHullCommonKernel< dim+1, TCoordinateInteger, TInternalInteger >
  {
    typedef ConvexHullCommonKernel< dim+1, TCoordinateInteger, TInternalInteger > Base;
    // inheriting types
    // using typename Base::Point;
    // using typename Base::Vector;
    // using typename Base::Scalar;
    using typename Base::CoordinatePoint;
    using typename Base::CoordinateVector;
    using typename Base::CoordinateScalar;
    using typename Base::InternalPoint;
    using typename Base::InternalVector;
    using typename Base::InternalScalar;
    using typename Base::Size;
    using typename Base::Index;
    using typename Base::IndexRange;
    using typename Base::CombinatorialPlaneSimplex;
    using typename Base::HalfSpace;
    // inheriting constants
    using Base::dimension;
    // inheriting methods
    using Base::compute;
    using Base::normal;
    using Base::intercept;
    using Base::dot;
    using Base::equal;
    using Base::height;
    using Base::volume;
    using Base::above;
    using Base::aboveOrOn;
    using Base::on;
 
    DelaunayIntegralKernel() = default;
 
    bool hasInfiniteFacets() const
    { return true; }
 
    bool isHalfSpaceFacetInfinite( const HalfSpace& hs ) const
    {
      return hs.internalNormal()[ dimension - 1 ] >= InternalScalar( 0 );
    }
    
    template < typename InputPoint>
    void makeInput( std::vector< CoordinatePoint >& processed_points,
                    IndexRange& input2comp, IndexRange& comp2input,
                    const std::vector< InputPoint >& input_points,
                    bool remove_duplicates )
    {
      const auto F = [&] ( InputPoint input ) -> CoordinatePoint
        {
          CoordinatePoint p;
          CoordinateScalar z = 0;
          for ( Dimension i = 0; i < dimension-1; i++ ) {
            const CoordinateScalar x = CoordinateScalar( input[ i ] );
            p[ i ]   = x;
            z       += x*x;
          }
          p[ dimension-1 ] = z;
          return p;
        };
      DGtal::detail::transform( processed_points, input2comp, comp2input,
                                input_points.cbegin(), input_points.cend(),
                                F, remove_duplicates );
    }
 
    template < typename OutputPoint>
    void convertPointTo( const CoordinatePoint& p, OutputPoint& out_p ) const
    {
      for ( Dimension k = 0; k < dimension-1; k++ )
        out_p[ k ] = p[ k ];
    }
        
  }; //   template < Dimension dim >  struct DelaunayIntegralKernel {
 
 
  // template class ConvexHullRationalKernel
  template < Dimension dim,
             typename TCoordinateInteger  = DGtal::int64_t,
             typename TInternalInteger = DGtal::int64_t >
  struct ConvexHullRationalKernel
    : public ConvexHullCommonKernel< dim, TCoordinateInteger, TInternalInteger >
  {
    typedef ConvexHullCommonKernel< dim, TCoordinateInteger, TInternalInteger > Base;
    // inheriting types
    // using typename Base::Point;
    // using typename Base::Vector;
    // using typename Base::Scalar;
    using typename Base::CoordinatePoint;
    using typename Base::CoordinateVector;
    using typename Base::CoordinateScalar;
    using typename Base::InternalPoint;
    using typename Base::InternalVector;
    using typename Base::InternalScalar;
    using typename Base::Size;
    using typename Base::Index;
    using typename Base::IndexRange;
    using typename Base::CombinatorialPlaneSimplex;
    using typename Base::HalfSpace;
    // inheriting constants
    using Base::dimension;
    // inheriting methods
    using Base::compute;
    using Base::normal;
    using Base::intercept;
    using Base::dot;
    using Base::equal;
    using Base::height;
    using Base::volume;
    using Base::above;
    using Base::aboveOrOn;
    using Base::on;
 
    double precision;
    
    ConvexHullRationalKernel( double aPrecision = 1024. )
      : precision( aPrecision ) {}
 
    bool hasInfiniteFacets() const
    { return false; }
 
    bool isHalfSpaceFacetInfinite( const HalfSpace& hs ) const
    {
      (void) hs; // unused parameter      
      return false;
    }
    
    template < typename InputPoint>
    void makeInput( std::vector< CoordinatePoint >& processed_points,
                    IndexRange& input2comp, IndexRange& comp2input,
                    const std::vector< InputPoint >& input_points,
                    bool remove_duplicates )
    {
      const auto F = [&] ( InputPoint input ) -> CoordinatePoint
        {
          CoordinatePoint p;
          for ( Dimension i = 0; i < dimension; i++ )
            p[ i ] = CoordinateScalar( round( input[ i ] * precision ) );
          return p;
        };
      DGtal::detail::transform( processed_points, input2comp, comp2input,
                                input_points.cbegin(), input_points.cend(),
                                F, remove_duplicates );
    }
 
    template < typename OutputPoint>
    void convertPointTo( const CoordinatePoint& p, OutputPoint& out_p ) const
    {
      for ( Dimension k = 0; k < dimension; k++ )
        out_p[ k ] = ( (double) p[ k ] ) / precision;
    }
    
    
  }; //   template < Dimension dim >  struct ConvexHullRationalKernel {
 
 
  // template class DelaunayRationalKernel
  template < Dimension dim,
             typename TCoordinateInteger  = DGtal::int64_t,
             typename TInternalInteger = DGtal::int64_t >
  struct DelaunayRationalKernel
    : public ConvexHullCommonKernel< dim+1, TCoordinateInteger, TInternalInteger >
  {
    typedef ConvexHullCommonKernel< dim+1, TCoordinateInteger, TInternalInteger > Base;
    // inheriting types
    // using typename Base::Point;
    // using typename Base::Vector;
    // using typename Base::Scalar;
    using typename Base::CoordinatePoint;
    using typename Base::CoordinateVector;
    using typename Base::CoordinateScalar;
    using typename Base::InternalPoint;
    using typename Base::InternalVector;
    using typename Base::InternalScalar;
    using typename Base::Size;
    using typename Base::Index;
    using typename Base::IndexRange;
    using typename Base::CombinatorialPlaneSimplex;
    using typename Base::HalfSpace;
    // inheriting constants
    using Base::dimension;
    // inheriting methods
    using Base::compute;
    using Base::normal;
    using Base::intercept;
    using Base::dot;
    using Base::equal;
    using Base::height;
    using Base::volume;
    using Base::above;
    using Base::aboveOrOn;
    using Base::on;
 
    double precision;
    
    DelaunayRationalKernel( double aPrecision = 1024. )
      : precision( aPrecision ) {}
 
    bool hasInfiniteFacets() const
    { return true; }
 
    bool isHalfSpaceFacetInfinite( const HalfSpace& hs ) const
    {
      return hs.internalNormal()[ dimension - 1 ] >= InternalScalar( 0 );
    }
 
    template < typename InputPoint>
    void makeInput( std::vector< CoordinatePoint >& processed_points,
                    IndexRange& input2comp, IndexRange& comp2input,
                    const std::vector< InputPoint >& input_points,
                    bool remove_duplicates )
    {
      const auto F = [&] ( InputPoint input ) -> CoordinatePoint
        {
          CoordinatePoint p;
          CoordinateScalar z = 0;
          for ( Dimension i = 0; i < dimension - 1; i++ ) {
            const CoordinateScalar x
              = CoordinateScalar( round( input[ i ] * precision ) );
            p[ i ] = x;
            z     += x*x;
          }
          p[ dimension-1 ] = z;
          return p;
        };
      DGtal::detail::transform( processed_points, input2comp, comp2input,
                                input_points.cbegin(), input_points.cend(),
                                F, remove_duplicates );
    }
 
    template < typename OutputPoint>
    void convertPointTo( const CoordinatePoint& p, OutputPoint& out_p ) const
    {
      for ( Dimension k = 0; k < dimension - 1; k++ )
        out_p[ k ] = ( (double) p[ k ] ) / precision;
    }
    
  }; //   template < Dimension dim >  struct DelaunayRationalKernel {
 
  
  
} // namespace DGtal {
 
#endif // !defined QuickHullKernels_h
 
#undef QuickHullKernels_RECURSES
#endif // else defined(QuickHullKernels_RECURSES)