testLinearizer.cpp

Go to the documentation of this file.

 
#include <cstddef>
#include <cmath>
 
#include "DGtalCatch.h"
 
#include <DGtal/kernel/SpaceND.h>
#include <DGtal/kernel/domains/HyperRectDomain.h>
#include <DGtal/kernel/domains/Linearizer.h>
 
using namespace DGtal;
 
// Previous version of the linearizer, used here for reference results.
namespace
{
 
  template < typename Domain, int dimension>
  struct linearizer
  {
 
    typedef typename Domain::Point Point;
    typedef typename Domain::Size Size;
 
    static Size apply( const Point & aPoint, const Point & lowerBound,
        const Point & extent )
    {
      Size pos = aPoint[ 0 ] - lowerBound[ 0 ] ;
      Size multiplier = 1;
      for (typename Domain::Dimension k = 1 ; k < dimension ; ++k)
      {
        multiplier *= extent[ k-1  ];
        pos += multiplier * ( aPoint[ k ] - lowerBound[ k ] );
      }
      return pos;
    }
  };
 
  template < typename Domain >
  struct linearizer< Domain, 1 >
  {
    typedef typename Domain::Point Point;
    typedef typename Domain::Size Size;
 
    static Size apply( const Point & aPoint,
        const Point & lowerBound,
        const Point & /*extent*/ )
    {
      return aPoint[ 0 ] - lowerBound[ 0 ];
    }
  };
 
  template < typename Domain >
  struct linearizer< Domain, 2 >
  {
    typedef typename Domain::Point Point;
    typedef typename Domain::Size Size;
 
    static Size apply( const Point & aPoint,
        const Point & lowerBound,
        const Point & extent )
    {
      return ( aPoint[ 0 ] - lowerBound[ 0 ] ) + extent[ 0 ] *
        (aPoint[ 1 ] - lowerBound[ 1 ] );
    }
  };
 
  template < typename Domain >
  struct linearizer< Domain, 3 >
  {
    typedef typename Domain::Point Point;
    typedef typename Domain::Size Size;
 
    static Size apply( const Point & aPoint,
        const Point & lowerBound,
        const Point & extent )
    {
      Size res = aPoint[ 0 ] - lowerBound[ 0 ];
      Size multiplier = extent[ 0 ];
      res += multiplier * ( aPoint[ 1 ] - lowerBound[ 1 ] );
      multiplier *= extent[ 1 ];
      res += multiplier * ( aPoint[ 2 ] - lowerBound[ 2 ] );
      return res;
    }
  };
}
 
 
template < typename StorageOrder >
struct PointConverter;
 
template <>
struct PointConverter<ColMajorStorage>
{
  template < typename TPoint >
  static inline
  TPoint apply( TPoint const& aPoint )
    {
      return aPoint;
    }
};
 
template <>
struct PointConverter<RowMajorStorage>
{
  template < typename TPoint >
  static inline
  TPoint apply( TPoint const& aPoint )
    {
      TPoint result;
      for ( auto i = 0 ; i < (int)TPoint::dimension ; ++i )
        result[i] = aPoint[ TPoint::dimension - i - 1 ];
 
      return result;
    }
};
 
#define TEST_LINEARIZER( N , ORDER ) \
TEST_CASE( "Testing Linearizer in dimension " #N " with " #ORDER, "[test][dim" #N "][" #ORDER "]" )\
{\
\
  typedef SpaceND<N>              Space;\
  typedef HyperRectDomain<Space>  Domain;\
  typedef Space::Point   Point;\
\
  typedef linearizer<Domain, N>   RefLinearizer;\
  typedef Linearizer<Domain, ORDER>   NewLinearizer;\
\
  typedef PointConverter<ORDER>   RefConverter;\
\
  std::size_t size = 1e3;\
\
  Point lowerBound;\
  for ( std::size_t i = 0 ; i < N ; ++i )\
    lowerBound[i] = 1 + 7*i;\
\
  std::size_t dim_size = std::size_t( std::pow( double(size), 1./N ) + 0.5 );\
  Point upperBound;\
  for ( std::size_t i = 0; i < N ; ++i )\
    upperBound[i] = lowerBound[i] + dim_size + i;\
\
  Domain domain( lowerBound, upperBound );\
  Point extent = upperBound - lowerBound + Point::diagonal(1);\
\
  Point refLowerBound = RefConverter::apply(lowerBound);\
  Point refExtent     = RefConverter::apply(extent);\
\
  SECTION( "Testing getIndex(Point, Point, Extent) syntax" )\
    {\
      for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
          REQUIRE(  RefLinearizer::apply( RefConverter::apply(*it), refLowerBound, refExtent ) == NewLinearizer::getIndex( *it, lowerBound, extent )  );\
    }\
\
  SECTION( "Testing getIndex(Point, Extent) syntax" )\
    {\
      for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
          REQUIRE(  RefLinearizer::apply( RefConverter::apply(*it), refLowerBound, refExtent ) == NewLinearizer::getIndex( *it - lowerBound, extent )  );\
    }\
\
  SECTION( "Testing getIndex(Point, Domain) syntax" )\
    {\
      for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
          REQUIRE(  RefLinearizer::apply( RefConverter::apply(*it), refLowerBound, refExtent ) == NewLinearizer::getIndex( *it, domain )  );\
    }\
\
  SECTION( "Testing getPoint(Index, Point, Extent) syntax" )\
    {\
      for ( std::size_t i = 0; i < domain.size(); ++i )\
          REQUIRE(  RefLinearizer::apply( RefConverter::apply( NewLinearizer::getPoint( i, lowerBound, extent ) ), refLowerBound, refExtent ) == i  );\
    }\
\
  SECTION( "Testing getPoint(Index, Extent) syntax" )\
    {\
      for ( std::size_t i = 0; i < domain.size(); ++i )\
          REQUIRE(  RefLinearizer::apply( RefConverter::apply( NewLinearizer::getPoint( i, extent ) + lowerBound ), refLowerBound, refExtent ) == i  );\
    }\
\
  SECTION( "Testing getPoint(Index, Domain) syntax" )\
    {\
      for ( std::size_t i = 0; i < domain.size(); ++i )\
          REQUIRE(  RefLinearizer::apply( RefConverter::apply( NewLinearizer::getPoint( i, domain ) ), refLowerBound, refExtent ) == i  );\
    }\
}
 
#define BENCH_LINEARIZER( N , ORDER ) \
TEST_CASE( "Benchmarking Linearizer in dimension " #N " with " #ORDER, "[.bench][dim" #N "][" #ORDER "]" )\
{\
\
  typedef SpaceND<N>              Space;\
  typedef HyperRectDomain<Space>  Domain;\
  typedef Space::Point   Point;\
\
  typedef linearizer<Domain, N>   RefLinearizer;\
  typedef Linearizer<Domain, ORDER>   NewLinearizer;\
\
  typedef PointConverter<ORDER>   RefConverter;\
\
  std::size_t size = 1e8;\
\
  Point lowerBound;\
  for ( std::size_t i = 0 ; i < N ; ++i )\
    lowerBound[i] = 1 + 7*i;\
\
  std::size_t dim_size = std::size_t( std::pow( double(size), 1./N ) + 0.5 );\
  Point upperBound;\
  for ( std::size_t i = 0; i < N ; ++i )\
    upperBound[i] = lowerBound[i] + dim_size + i;\
\
  Domain domain( lowerBound, upperBound );\
  Point extent = upperBound - lowerBound + Point::diagonal(1);\
\
  Point refLowerBound = RefConverter::apply(lowerBound);\
  Point refExtent     = RefConverter::apply(extent);\
\
  std::size_t sum = 0;\
\
  for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
      sum += RefLinearizer::apply( RefConverter::apply(*it), refLowerBound, refExtent );\
  REQUIRE( sum > 0 );\
  sum = 0;\
\
  SECTION( "Benchmarking reference linearizer" )\
    {\
      for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
          sum += RefLinearizer::apply( RefConverter::apply(*it), refLowerBound, refExtent );\
    }\
\
  SECTION( "Benchmarking getIndex(Point, Point, Extent) syntax" )\
    {\
      for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
          sum += NewLinearizer::getIndex( *it, lowerBound, extent );\
    }\
\
  SECTION( "Benchmarking getIndex(Point, Extent) syntax" )\
    {\
      for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
          sum += NewLinearizer::getIndex( *it, extent );\
    }\
\
  SECTION( "Benchmarking getIndex(Point, Domain) syntax" )\
    {\
      for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
          sum += NewLinearizer::getIndex( *it, domain );\
    }\
\
  SECTION( "Benchmarking getPoint(Index, Point, Extent) syntax" )\
    {\
      for ( std::size_t i = 0; i < domain.size(); ++i )\
          sum += NewLinearizer::getPoint( i, lowerBound, extent )[N-1];\
    }\
\
  SECTION( "Benchmarking getPoint(Index, Extent) syntax" )\
    {\
      for ( std::size_t i = 0; i < domain.size(); ++i )\
          sum += NewLinearizer::getPoint( i, extent )[N-1];\
    }\
\
  SECTION( "Benchmarking getPoint(Index, Domain) syntax" )\
    {\
      for ( std::size_t i = 0; i < domain.size(); ++i )\
          sum += NewLinearizer::getPoint( i, domain )[N-1];\
    }\
\
  REQUIRE( sum > 0 );\
}
 
TEST_LINEARIZER( 1, ColMajorStorage )
TEST_LINEARIZER( 2, ColMajorStorage )
TEST_LINEARIZER( 3, ColMajorStorage )
TEST_LINEARIZER( 4, ColMajorStorage )
TEST_LINEARIZER( 5, ColMajorStorage )
 
TEST_LINEARIZER( 1, RowMajorStorage )
TEST_LINEARIZER( 2, RowMajorStorage )
TEST_LINEARIZER( 3, RowMajorStorage )
TEST_LINEARIZER( 4, RowMajorStorage )
TEST_LINEARIZER( 5, RowMajorStorage )
 
BENCH_LINEARIZER( 1, ColMajorStorage )
BENCH_LINEARIZER( 2, ColMajorStorage )
BENCH_LINEARIZER( 3, ColMajorStorage )
BENCH_LINEARIZER( 4, ColMajorStorage )
BENCH_LINEARIZER( 5, ColMajorStorage )
 
BENCH_LINEARIZER( 1, RowMajorStorage )
BENCH_LINEARIZER( 2, RowMajorStorage )
BENCH_LINEARIZER( 3, RowMajorStorage )
BENCH_LINEARIZER( 4, RowMajorStorage )
BENCH_LINEARIZER( 5, RowMajorStorage )