DGtal  0.9.2
testLinearizer.cpp
1 
33 #include <cstddef>
34 #include <cmath>
35 
36 #include "DGtalCatch.h"
37 
38 #include <DGtal/kernel/SpaceND.h>
39 #include <DGtal/kernel/domains/HyperRectDomain.h>
40 #include <DGtal/kernel/domains/Linearizer.h>
41 
42 using namespace DGtal;
43 
44 // Previous version of the linearizer, used here for reference results.
45 namespace
46 {
47 
56  template < typename Domain, int dimension>
57  struct linearizer
58  {
59 
60  typedef typename Domain::Point Point;
61  typedef typename Domain::Size Size;
62 
72  static Size apply( const Point & aPoint, const Point & lowerBound,
73  const Point & extent )
74  {
75  Size pos = aPoint[ 0 ] - lowerBound[ 0 ] ;
76  Size multiplier = 1;
77  for (typename Domain::Dimension k = 1 ; k < dimension ; ++k)
78  {
79  multiplier *= extent[ k-1 ];
80  pos += multiplier * ( aPoint[ k ] - lowerBound[ k ] );
81  }
82  return pos;
83  }
84  };
85 
90  template < typename Domain >
91  struct linearizer< Domain, 1 >
92  {
93  typedef typename Domain::Point Point;
94  typedef typename Domain::Size Size;
95 
96  static Size apply( const Point & aPoint,
97  const Point & lowerBound,
98  const Point & /*extent*/ )
99  {
100  return aPoint[ 0 ] - lowerBound[ 0 ];
101  }
102  };
103 
108  template < typename Domain >
109  struct linearizer< Domain, 2 >
110  {
111  typedef typename Domain::Point Point;
112  typedef typename Domain::Size Size;
113 
114  static Size apply( const Point & aPoint,
115  const Point & lowerBound,
116  const Point & extent )
117  {
118  return ( aPoint[ 0 ] - lowerBound[ 0 ] ) + extent[ 0 ] *
119  (aPoint[ 1 ] - lowerBound[ 1 ] );
120  }
121  };
122 
127  template < typename Domain >
128  struct linearizer< Domain, 3 >
129  {
130  typedef typename Domain::Point Point;
131  typedef typename Domain::Size Size;
132 
133  static Size apply( const Point & aPoint,
134  const Point & lowerBound,
135  const Point & extent )
136  {
137  Size res = aPoint[ 0 ] - lowerBound[ 0 ];
138  Size multiplier = extent[ 0 ];
139  res += multiplier * ( aPoint[ 1 ] - lowerBound[ 1 ] );
140  multiplier *= extent[ 1 ];
141  res += multiplier * ( aPoint[ 2 ] - lowerBound[ 2 ] );
142  return res;
143  }
144  };
145 }
146 
147 
149 template < typename StorageOrder >
150 struct PointConverter;
151 
152 template <>
153 struct PointConverter<ColMajorStorage>
154 {
155  template < typename TPoint >
156  static inline
157  TPoint apply( TPoint const& aPoint )
158  {
159  return aPoint;
160  }
161 };
162 
163 template <>
164 struct PointConverter<RowMajorStorage>
165 {
166  template < typename TPoint >
167  static inline
168  TPoint apply( TPoint const& aPoint )
169  {
170  TPoint result;
171  for ( std::size_t i = 0 ; i < TPoint::dimension ; ++i )
172  result[i] = aPoint[ TPoint::dimension - i - 1 ];
173 
174  return result;
175  }
176 };
177 
178 #define TEST_LINEARIZER( N , ORDER ) \
179 TEST_CASE( "Testing Linearizer in dimension " #N " with " #ORDER, "[test][dim" #N "][" #ORDER "]" )\
180 {\
181 \
182  typedef SpaceND<N> Space;\
183  typedef HyperRectDomain<Space> Domain;\
184  typedef Space::Point Point;\
185 \
186  typedef linearizer<Domain, N> RefLinearizer;\
187  typedef Linearizer<Domain, ORDER> NewLinearizer;\
188 \
189  typedef PointConverter<ORDER> RefConverter;\
190 \
191  std::size_t size = 1e3;\
192 \
193  Point lowerBound;\
194  for ( std::size_t i = 0 ; i < N ; ++i )\
195  lowerBound[i] = 1 + 7*i;\
196 \
197  std::size_t dim_size = std::size_t( std::pow( double(size), 1./N ) + 0.5 );\
198  Point upperBound;\
199  for ( std::size_t i = 0; i < N ; ++i )\
200  upperBound[i] = lowerBound[i] + dim_size + i;\
201 \
202  Domain domain( lowerBound, upperBound );\
203  Point extent = upperBound - lowerBound + Point::diagonal(1);\
204 \
205  Point refLowerBound = RefConverter::apply(lowerBound);\
206  Point refExtent = RefConverter::apply(extent);\
207 \
208  SECTION( "Testing getIndex(Point, Point, Extent) syntax" )\
209  {\
210  for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
211  REQUIRE( RefLinearizer::apply( RefConverter::apply(*it), refLowerBound, refExtent ) == NewLinearizer::getIndex( *it, lowerBound, extent ) );\
212  }\
213 \
214  SECTION( "Testing getIndex(Point, Extent) syntax" )\
215  {\
216  for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
217  REQUIRE( RefLinearizer::apply( RefConverter::apply(*it), refLowerBound, refExtent ) == NewLinearizer::getIndex( *it - lowerBound, extent ) );\
218  }\
219 \
220  SECTION( "Testing getIndex(Point, Domain) syntax" )\
221  {\
222  for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
223  REQUIRE( RefLinearizer::apply( RefConverter::apply(*it), refLowerBound, refExtent ) == NewLinearizer::getIndex( *it, domain ) );\
224  }\
225 \
226  SECTION( "Testing getPoint(Index, Point, Extent) syntax" )\
227  {\
228  for ( std::size_t i = 0; i < domain.size(); ++i )\
229  REQUIRE( RefLinearizer::apply( RefConverter::apply( NewLinearizer::getPoint( i, lowerBound, extent ) ), refLowerBound, refExtent ) == i );\
230  }\
231 \
232  SECTION( "Testing getPoint(Index, Extent) syntax" )\
233  {\
234  for ( std::size_t i = 0; i < domain.size(); ++i )\
235  REQUIRE( RefLinearizer::apply( RefConverter::apply( NewLinearizer::getPoint( i, extent ) + lowerBound ), refLowerBound, refExtent ) == i );\
236  }\
237 \
238  SECTION( "Testing getPoint(Index, Domain) syntax" )\
239  {\
240  for ( std::size_t i = 0; i < domain.size(); ++i )\
241  REQUIRE( RefLinearizer::apply( RefConverter::apply( NewLinearizer::getPoint( i, domain ) ), refLowerBound, refExtent ) == i );\
242  }\
243 }
244 
245 #define BENCH_LINEARIZER( N , ORDER ) \
246 TEST_CASE( "Benchmarking Linearizer in dimension " #N " with " #ORDER, "[.bench][dim" #N "][" #ORDER "]" )\
247 {\
248 \
249  typedef SpaceND<N> Space;\
250  typedef HyperRectDomain<Space> Domain;\
251  typedef Space::Point Point;\
252 \
253  typedef linearizer<Domain, N> RefLinearizer;\
254  typedef Linearizer<Domain, ORDER> NewLinearizer;\
255 \
256  typedef PointConverter<ORDER> RefConverter;\
257 \
258  std::size_t size = 1e8;\
259 \
260  Point lowerBound;\
261  for ( std::size_t i = 0 ; i < N ; ++i )\
262  lowerBound[i] = 1 + 7*i;\
263 \
264  std::size_t dim_size = std::size_t( std::pow( double(size), 1./N ) + 0.5 );\
265  Point upperBound;\
266  for ( std::size_t i = 0; i < N ; ++i )\
267  upperBound[i] = lowerBound[i] + dim_size + i;\
268 \
269  Domain domain( lowerBound, upperBound );\
270  Point extent = upperBound - lowerBound + Point::diagonal(1);\
271 \
272  Point refLowerBound = RefConverter::apply(lowerBound);\
273  Point refExtent = RefConverter::apply(extent);\
274 \
275  std::size_t sum = 0;\
276 \
277  for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
278  sum += RefLinearizer::apply( RefConverter::apply(*it), refLowerBound, refExtent );\
279  REQUIRE( sum > 0 );\
280  sum = 0;\
281 \
282  SECTION( "Benchmarking reference linearizer" )\
283  {\
284  for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
285  sum += RefLinearizer::apply( RefConverter::apply(*it), refLowerBound, refExtent );\
286  }\
287 \
288  SECTION( "Benchmarking getIndex(Point, Point, Extent) syntax" )\
289  {\
290  for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
291  sum += NewLinearizer::getIndex( *it, lowerBound, extent );\
292  }\
293 \
294  SECTION( "Benchmarking getIndex(Point, Extent) syntax" )\
295  {\
296  for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
297  sum += NewLinearizer::getIndex( *it, extent );\
298  }\
299 \
300  SECTION( "Benchmarking getIndex(Point, Domain) syntax" )\
301  {\
302  for ( Domain::ConstIterator it = domain.begin(), it_end = domain.end(); it != it_end ; ++it )\
303  sum += NewLinearizer::getIndex( *it, domain );\
304  }\
305 \
306  SECTION( "Benchmarking getPoint(Index, Point, Extent) syntax" )\
307  {\
308  for ( std::size_t i = 0; i < domain.size(); ++i )\
309  sum += NewLinearizer::getPoint( i, lowerBound, extent )[N-1];\
310  }\
311 \
312  SECTION( "Benchmarking getPoint(Index, Extent) syntax" )\
313  {\
314  for ( std::size_t i = 0; i < domain.size(); ++i )\
315  sum += NewLinearizer::getPoint( i, extent )[N-1];\
316  }\
317 \
318  SECTION( "Benchmarking getPoint(Index, Domain) syntax" )\
319  {\
320  for ( std::size_t i = 0; i < domain.size(); ++i )\
321  sum += NewLinearizer::getPoint( i, domain )[N-1];\
322  }\
323 \
324  REQUIRE( sum > 0 );\
325 }
326 
327 TEST_LINEARIZER( 1, ColMajorStorage )
328 TEST_LINEARIZER( 2, ColMajorStorage )
329 TEST_LINEARIZER( 3, ColMajorStorage )
330 TEST_LINEARIZER( 4, ColMajorStorage )
331 TEST_LINEARIZER( 5, ColMajorStorage )
332 
333 TEST_LINEARIZER( 1, RowMajorStorage )
334 TEST_LINEARIZER( 2, RowMajorStorage )
335 TEST_LINEARIZER( 3, RowMajorStorage )
336 TEST_LINEARIZER( 4, RowMajorStorage )
337 TEST_LINEARIZER( 5, RowMajorStorage )
338 
339 BENCH_LINEARIZER( 1, ColMajorStorage )
340 BENCH_LINEARIZER( 2, ColMajorStorage )
341 BENCH_LINEARIZER( 3, ColMajorStorage )
342 BENCH_LINEARIZER( 4, ColMajorStorage )
343 BENCH_LINEARIZER( 5, ColMajorStorage )
344 
345 BENCH_LINEARIZER( 1, RowMajorStorage )
346 BENCH_LINEARIZER( 2, RowMajorStorage )
347 BENCH_LINEARIZER( 3, RowMajorStorage )
348 BENCH_LINEARIZER( 4, RowMajorStorage )
349 BENCH_LINEARIZER( 5, RowMajorStorage )
Tag (empty structure) specifying a row-major storage order.
Definition: Linearizer.h:54
Tag (empty structure) specifying a col-major storage order.
Definition: Linearizer.h:61
DGtal is the top-level namespace which contains all DGtal functions and types.