[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

multi_array_chunked.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 2012-2014 by Ullrich Koethe and Thorben Kroeger */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 /* benchmark results for a simple loop 'if(iter.get<1>() != count++)'
37 
38  ********************
39  image size: 200^3, chunk size: 64^3, i.e. chunk count: 4^3
40  times in msec, excluding time to store file on disk
41 
42  win64/vs2012 (koethe-laptop): uint8 float double
43  plain array 18 18 18
44  chunked array (all in cache) 25 26 26
45  thread-safe chunked (all in cache) 27 28 29
46  thread-safe chunked (1 slice in cache) 29 33 39
47  thread-safe chunked (1 row in cache) 45 48 52
48  chunked (initial creation, all in cache) 33 43 57
49 
50  linux/gcc 4.7.3 (birdofprey): uint8 float double
51  plain array 16 20 21
52  chunked array (all in cache) 17 23 24
53  thread-safe chunked (all in cache) 19 24 25
54  thread-safe chunked (1 slice in cache) 20 29 34
55  thread-safe chunked (1 row in cache) 24 33 39
56  chunked (initial creation, all in cache) 22 34 48
57 
58  OS X 10.7: uint8 float double
59  plain array 11 22 24
60  chunked array (all in cache) -- -- --
61  thread-safe chunked (all in cache) 20 25 26
62  thread-safe chunked (1 slice in cache) 23 37 46
63  thread-safe chunked (1 row in cache) 32 50 56
64  chunked (initial creation, all in cache) 34 56 77
65  (These numbers refer to nested loop iteration. Scan-order iteration
66  is unfortunately 3.5 times slower on the Mac. On the other hand,
67  two-level indexing as faster on a Mac than on Linux and Windows --
68  the speed penalty is only a factor of 2 rather than 3.)
69 
70  **********************
71  image size: 400^3, chunk size: 127^3, i.e. chunk count: 4^3
72  times in msec, excluding time to store file on disk
73 
74  win64/vs2012 (koethe-laptop): uint8 float double
75  plain array 130 130 130
76  chunked array (all in cache) 190 190 200
77  thread-safe chunked (all in cache) 190 200 210
78  thread-safe chunked (1 slice in cache) 210 235 280
79  thread-safe chunked (1 row in cache) 240 270 300
80  chunked (initial creation, all in cache) 230 300 400
81 
82  linux/gcc 4.7.3 (birdofprey): uint8 float double
83  plain array 130 162 165
84  chunked array (all in cache) 131 180 184
85  thread-safe chunked (all in cache) 135 183 188
86  thread-safe chunked (1 slice in cache) 146 218 258
87  thread-safe chunked (1 row in cache) 154 229 270
88  chunked (initial creation, all in cache) 173 269 372
89 
90  ***********************
91  Compression:
92  * I tried ZLIB, LZO, SNAPPY, LZ4, LZFX and FASTLZ (with compression levels 1 -- faster
93  and level 2 -- higher compression). There are also QuickLZ and LZMAT which claim
94  to be fast, but they are under a GPL license.
95  * ZLIB compresses best, but is quite slow even at compression level 1
96  (times are in ms and include compression and decompression).
97  byte float double
98  ZLIB 121 3100 5800
99  ZLIB1 79 1110 1190
100  LZO 43 104 280
101  SNAPPY 46 71 305
102  LZ4 42 70 283
103  LZFX 76 278 330
104  FASTLZ1 52 280 309
105  FASTLZ1 53 286 339
106  * The fast compression algorithms are unable to compress the float array
107  and achieve ~50% for the double array, whereas ZLIB achieves 32% and 16%
108  respectively (at the fastest compression level 1, it is still 33% and 17%
109  respectively). LZFX cannot even compress the byte data (probably a bug?).
110  Average compression ratios for the byte array are
111  ZLIB: 2.3%
112  ZLIB1: 4.6%
113  LZO: 5.6%
114  SNAPPY: 9.8%
115  LZ4: 9.7%
116  FASTLZ1: 7.6%
117  FASTLZ2: 7.9%
118  * LZO is under GPL (but there is a Java implementation under Apache license at
119  http://svn.apache.org/repos/asf/hadoop/common/tags/release-0.19.2/src/core/org/apache/hadoop/io/compress/lzo/)
120  The others are BSD and MIT (FASTLZ).
121  * Snappy doesn't support Windows natively, but porting is simple (see my github repo)
122  * The source code for LZO, LZ4, LZFX, and FASTLZ can simply be copied to VIGRA,
123  but LZO's GPL license is unsuitable.
124  * HDF5 compression is already sufficient at level 1 (4-15%,
125  higher levels don't lead to big gains) and only a factor 3-10 slower
126  than without compression.
127 */
128 
129 #ifndef VIGRA_MULTI_ARRAY_CHUNKED_HXX
130 #define VIGRA_MULTI_ARRAY_CHUNKED_HXX
131 
132 #include <queue>
133 #include <string>
134 
135 #include "multi_fwd.hxx"
136 #include "multi_handle.hxx"
137 #include "multi_array.hxx"
138 #include "memory.hxx"
139 #include "metaprogramming.hxx"
140 #include "threading.hxx"
141 #include "compression.hxx"
142 
143 // // FIXME: why is this needed when compiling the Python bindng,
144 // // but not when compiling test_multiarray_chunked?
145 // #if defined(__GNUC__)
146 // # define memory_order_release memory_order_seq_cst
147 // # define memory_order_acquire memory_order_seq_cst
148 // #endif
149 
150 #ifdef _WIN32
151 # include "windows.h"
152 #else
153 # include <fcntl.h>
154 # include <stdlib.h>
155 # include <unistd.h>
156 # include <sys/stat.h>
157 # include <sys/mman.h>
158 #endif
159 
160 // Bounds checking Macro used if VIGRA_CHECK_BOUNDS is defined.
161 #ifdef VIGRA_CHECK_BOUNDS
162 #define VIGRA_ASSERT_INSIDE(diff) \
163  vigra_precondition(this->isInside(diff), "Index out of bounds")
164 #else
165 #define VIGRA_ASSERT_INSIDE(diff)
166 #endif
167 
168 
169 namespace vigra {
170 
171 #ifdef __APPLE__
172  #define VIGRA_NO_SPARSE_FILE
173 #endif
174 
175 #ifdef _WIN32
176 
177 inline
178 void winErrorToException(std::string message = "")
179 {
180  LPVOID lpMsgBuf;
181  DWORD dw = GetLastError();
182 
183  FormatMessage(
184  FORMAT_MESSAGE_ALLOCATE_BUFFER |
185  FORMAT_MESSAGE_FROM_SYSTEM |
186  FORMAT_MESSAGE_IGNORE_INSERTS,
187  NULL,
188  dw,
189  MAKELANGID(LANG_NEUTRAL, SUBLANG_DEFAULT),
190  (LPTSTR) &lpMsgBuf,
191  0, NULL );
192 
193  message += (char*)lpMsgBuf;
194  LocalFree(lpMsgBuf);
195 
196  throw std::runtime_error(message);
197 }
198 
199 inline
200 std::string winTempFileName(std::string path = "")
201 {
202  if(path == "")
203  {
204  TCHAR default_path[MAX_PATH];
205  if(!GetTempPath(MAX_PATH, default_path))
206  winErrorToException("winTempFileName(): ");
207  path = default_path;
208  }
209 
210  TCHAR name[MAX_PATH];
211  if(!GetTempFileName(path.c_str(), TEXT("vigra"), 0, name))
212  winErrorToException("winTempFileName(): ");
213 
214  return std::string(name);
215 }
216 
217 inline
218 std::size_t winClusterSize()
219 {
220  SYSTEM_INFO info;
221  ::GetSystemInfo(&info);
222  return info.dwAllocationGranularity;
223 }
224 
225 #endif
226 
227 namespace {
228 
229 #ifdef _WIN32
230 std::size_t mmap_alignment = winClusterSize();
231 #else
232 std::size_t mmap_alignment = sysconf(_SC_PAGE_SIZE);
233 #endif
234 
235 } // anonymous namespace
236 
237 template <unsigned int N, class T>
238 class IteratorChunkHandle;
239 
240 namespace detail {
241 
242 template <unsigned int N>
243 struct ChunkIndexing
244 {
245  template <class T, int M>
246  static void chunkIndex(TinyVector<T, M> const & p,
247  TinyVector<T, M> const & bits,
248  TinyVector<T, M> & index)
249  {
250  typedef std::size_t UI;
251  ChunkIndexing<N-1>::chunkIndex(p, bits, index);
252  index[N-1] = (UI)p[N-1] >> bits[N-1];
253  }
254 
255  template <class T, int M>
256  static std::size_t chunkOffset(TinyVector<T, M> const & p,
257  TinyVector<T, M> const & bits,
258  TinyVector<T, M> const & strides)
259  {
260  typedef std::size_t UI;
261  return ChunkIndexing<N-1>::chunkOffset(p, bits, strides) +
262  ((UI)p[N-1] >> bits[N-1]) * strides[N-1];
263  }
264 
265  template <class T, int M>
266  static std::size_t offsetInChunk(TinyVector<T, M> const & p,
267  TinyVector<T, M> const & mask,
268  TinyVector<T, M> const & strides)
269  {
270  typedef std::size_t UI;
271  return ChunkIndexing<N-1>::offsetInChunk(p, mask, strides) +
272  ((UI)p[N-1] & (UI)mask[N-1]) * strides[N-1];
273  }
274 };
275 
276 template <>
277 struct ChunkIndexing<1>
278 {
279  template <class T, int M>
280  static void chunkIndex(TinyVector<T, M> const & p,
281  TinyVector<T, M> const & bits,
282  TinyVector<T, M> & index)
283  {
284  typedef std::size_t UI;
285  index[0] = (UI)p[0] >> bits[0];
286  }
287 
288  template <class T, int M>
289  static std::size_t chunkOffset(TinyVector<T, M> const & p,
290  TinyVector<T, M> const & bits,
291  TinyVector<T, M> const & strides)
292  {
293  typedef std::size_t UI;
294  return ((UI)p[0] >> bits[0]) * strides[0];
295  }
296 
297  template <class T, int M>
298  static std::size_t offsetInChunk(TinyVector<T, M> const & p,
299  TinyVector<T, M> const & mask,
300  TinyVector<T, M> const & strides)
301  {
302  typedef std::size_t UI;
303  return ((UI)p[0] & (UI)mask[0]) * strides[0];
304  }
305 };
306 
307 template <class T, int M>
308 inline TinyVector<T, M>
309 computeChunkArrayShape(TinyVector<T, M> shape,
310  TinyVector<T, M> const & bits,
311  TinyVector<T, M> const & mask)
312 {
313  for(int k=0; k<M; ++k)
314  shape[k] = (shape[k] + mask[k]) >> bits[k];
315  return shape;
316 }
317 
318 template <class T, int M>
319 inline T
320 defaultCacheSize(TinyVector<T, M> const & shape)
321 {
322  T res = max(shape);
323  for(int k=0; k<M-1; ++k)
324  for(int j=k+1; j<M; ++j)
325  res = std::max(res, shape[k]*shape[j]);
326  return res + 1;
327 }
328 
329 } // namespace detail
330 
331 template <unsigned int N, class T>
332 class ChunkBase
333 {
334  public:
335  typedef typename MultiArrayShape<N>::type shape_type;
336  typedef T value_type;
337  typedef T* pointer;
338 
339  ChunkBase()
340  : strides_()
341  , pointer_()
342  {}
343 
344  ChunkBase(shape_type const & strides, pointer p = 0)
345  : strides_(strides)
346  , pointer_(p)
347  {}
348 
349  typename MultiArrayShape<N>::type strides_;
350  T * pointer_;
351 };
352 
353 template <unsigned int N, class T>
354 class SharedChunkHandle
355 {
356  public:
357  typedef typename MultiArrayShape<N>::type shape_type;
358 
359  static const long chunk_asleep = -2;
360  static const long chunk_uninitialized = -3;
361  static const long chunk_locked = -4;
362  static const long chunk_failed = -5;
363 
364  SharedChunkHandle()
365  : pointer_(0)
366  , chunk_state_()
367  {
368  chunk_state_ = chunk_uninitialized;
369  }
370 
371  SharedChunkHandle(SharedChunkHandle const & rhs)
372  : pointer_(rhs.pointer_)
373  , chunk_state_()
374  {
375  chunk_state_ = chunk_uninitialized;
376  }
377 
378  shape_type const & strides() const
379  {
380  return pointer_->strides_;
381  }
382 
383  ChunkBase<N, T> * pointer_;
384  mutable threading::atomic_long chunk_state_;
385 
386  private:
387  SharedChunkHandle & operator=(SharedChunkHandle const & rhs);
388 };
389 
390 template <unsigned int N, class T>
391 class ChunkedArrayBase
392 {
393  public:
394  enum ActualDimension{ actual_dimension = (N == 0) ? 1 : N };
395  typedef typename MultiArrayShape<N>::type shape_type;
396  typedef T value_type;
397  typedef value_type * pointer;
398  typedef value_type & reference;
399  typedef ChunkBase<N, T> Chunk;
400 
401  ChunkedArrayBase()
402  : shape_()
403  , chunk_shape_()
404  {}
405 
406  ChunkedArrayBase(shape_type const & shape, shape_type const & chunk_shape)
407  : shape_(shape)
408  , chunk_shape_(prod(chunk_shape) > 0 ? chunk_shape : detail::ChunkShape<N, T>::defaultShape())
409  {}
410 
411  virtual ~ChunkedArrayBase()
412  {}
413 
414  virtual void unrefChunk(IteratorChunkHandle<N, T> * h) const = 0;
415 
416  virtual pointer chunkForIterator(shape_type const & point,
417  shape_type & strides, shape_type & upper_bound,
418  IteratorChunkHandle<N, T> * h) = 0;
419 
420  virtual pointer chunkForIterator(shape_type const & point,
421  shape_type & strides, shape_type & upper_bound,
422  IteratorChunkHandle<N, T> * h) const = 0;
423 
424  virtual std::string backend() const = 0;
425 
426  virtual shape_type chunkArrayShape() const = 0;
427 
428  virtual bool isReadOnly() const
429  {
430  return false;
431  }
432 
433  MultiArrayIndex size() const
434  {
435  return prod(shape_);
436  }
437 
438  shape_type const & shape() const
439  {
440  return shape_;
441  }
442 
443  MultiArrayIndex shape(MultiArrayIndex d) const
444  {
445  return shape_[d];
446  }
447 
448  shape_type const & chunkShape() const
449  {
450  return chunk_shape_;
451  }
452 
453  MultiArrayIndex chunkShape(MultiArrayIndex d) const
454  {
455  return chunk_shape_[d];
456  }
457 
458  bool isInside(shape_type const & p) const
459  {
460  for(int d=0; d<N; ++d)
461  if(p[d] < 0 || p[d] >= shape_[d])
462  return false;
463  return true;
464  }
465 
466  shape_type shape_, chunk_shape_;
467 };
468 
469 template <unsigned int N, class T>
470 class ChunkedArray;
471 
472 struct ChunkUnrefProxyBase
473 {
474  virtual ~ChunkUnrefProxyBase() {}
475 };
476 
477 template <unsigned int N, class T_MaybeConst>
478 class MultiArrayView<N, T_MaybeConst, ChunkedArrayTag>
479 : public ChunkedArrayBase<N, typename UnqualifiedType<T_MaybeConst>::type>
480 {
481  public:
482  enum ActualDimension { actual_dimension = (N==0) ? 1 : N };
483  typedef typename UnqualifiedType<T_MaybeConst>::type T;
484  typedef T value_type; // FIXME: allow Multiband<T> ???
485  typedef T_MaybeConst & reference;
486  typedef const value_type &const_reference;
487  typedef T_MaybeConst * pointer;
488  typedef const value_type *const_pointer;
489  typedef typename MultiArrayShape<actual_dimension>::type difference_type;
490  typedef difference_type key_type;
491  typedef difference_type size_type;
492  typedef difference_type shape_type;
493  typedef MultiArrayIndex difference_type_1;
494  typedef ChunkIterator<actual_dimension, T_MaybeConst> chunk_iterator;
495  typedef ChunkIterator<actual_dimension, T const> chunk_const_iterator;
496  typedef StridedScanOrderIterator<actual_dimension, ChunkedMemory<T_MaybeConst>, T_MaybeConst&, T_MaybeConst*> iterator;
497  typedef StridedScanOrderIterator<actual_dimension, ChunkedMemory<T const>, T const &, T const *> const_iterator;
498  typedef MultiArrayView<N, T_MaybeConst, ChunkedArrayTag> view_type;
499  typedef MultiArrayView<N, T const, ChunkedArrayTag> const_view_type;
500  typedef ChunkedArrayTag StrideTag;
501  typedef ChunkBase<N, T> Chunk;
502 
503  typedef MultiArray<N, Chunk> ChunkHolder;
504 
505  struct UnrefProxy
506  : public ChunkUnrefProxyBase
507  {
508  UnrefProxy(int size, ChunkedArray<N, T> * array)
509  : chunks_(size)
510  , array_(array)
511  {}
512 
513  ~UnrefProxy()
514  {
515  if(array_)
516  array_->unrefChunks(chunks_);
517  }
518 
519  ArrayVector<SharedChunkHandle<N, T> *> chunks_;
520  ChunkedArray<N, T> * array_;
521  };
522 
523  virtual shape_type chunkArrayShape() const
524  {
525  return chunks_.shape();
526  }
527 
528  shape_type chunkStart(shape_type const & global_start) const
529  {
530  shape_type chunk_start(SkipInitialization);
531  detail::ChunkIndexing<N>::chunkIndex(global_start, bits_, chunk_start);
532  return chunk_start;
533  }
534 
535  shape_type chunkStop(shape_type global_stop) const
536  {
537  global_stop -= shape_type(1);
538  shape_type chunk_stop(SkipInitialization);
539  detail::ChunkIndexing<N>::chunkIndex(global_stop, bits_, chunk_stop);
540  chunk_stop += shape_type(1);
541  return chunk_stop;
542  }
543 
544  virtual void unrefChunk(IteratorChunkHandle<N, T> *) const {}
545 
546  virtual T* chunkForIterator(shape_type const & point,
547  shape_type & strides, shape_type & upper_bound,
548  IteratorChunkHandle<N, T> * h)
549  {
550  return const_cast<MultiArrayView const *>(this)->chunkForIterator(point, strides, upper_bound, h);
551  }
552 
553  virtual T* chunkForIterator(shape_type const & point,
554  shape_type & strides, shape_type & upper_bound,
555  IteratorChunkHandle<N, T> * h) const
556  {
557  shape_type global_point = point + h->offset_;
558 
559  if(!this->isInside(global_point))
560  {
561  upper_bound = point + this->chunk_shape_;
562  return 0;
563  }
564 
565  global_point += offset_;
566  shape_type coffset = offset_ + h->offset_;
567 
568  shape_type chunkIndex = chunkStart(global_point);
569  Chunk const * chunk = &chunks_[chunkIndex];
570  strides = chunk->strides_;
571  upper_bound = (chunkIndex + shape_type(1)) * this->chunk_shape_ - coffset;
572  std::size_t offset = detail::ChunkIndexing<N>::offsetInChunk(global_point, mask_, strides);
573  return const_cast<T*>(chunk->pointer_ + offset);
574  }
575 
576  virtual std::string backend() const
577  {
578  return "MultiArrayView<ChunkedArrayTag>";
579  }
580 
581  MultiArrayView()
582  : ChunkedArrayBase<N, T>()
583  {}
584 
585  MultiArrayView(shape_type const & shape, shape_type const & chunk_shape)
586  : ChunkedArrayBase<N, T>(shape, chunk_shape)
587  {}
588 
589  MultiArrayView & operator=(MultiArrayView const & rhs)
590  {
591  if(this != &rhs)
592  {
593  if(!hasData())
594  {
595  ChunkedArrayBase<N, T>::operator=(rhs);
596  chunks_ = rhs.chunks_;
597  offset_ = rhs.offset_;
598  bits_ = rhs.bits_;
599  mask_ = rhs.mask_;
600  unref_ = rhs.unref_;
601  }
602  else
603  {
604  vigra_precondition(this->shape() == rhs.shape(),
605  "MultiArrayView::operator=(): shape mismatch.");
606  iterator i = begin(), ie = end();
607  const_iterator j = rhs.begin();
608  for(; i != ie; ++i, ++j)
609  *i = *j;
610  }
611  }
612  return *this;
613  }
614 
615  #define VIGRA_CHUNKED_ARRAY_VIEW_ASSIGN(op) \
616  template<class U, class C1> \
617  MultiArrayView & operator op(MultiArrayView<N, U, C1> const & rhs) \
618  { \
619  vigra_precondition(this->shape() == rhs.shape(), \
620  "MultiArrayView::operator" #op "(): shape mismatch."); \
621  iterator i = begin(), ie = end(); \
622  typename MultiArrayView<N, U, C1>::const_iterator j = rhs.begin(); \
623  for(; i != ie; ++i, ++j) \
624  *i op detail::RequiresExplicitCast<value_type>::cast(*j); \
625  return *this; \
626  } \
627  \
628  MultiArrayView & operator op(value_type const & v) \
629  { \
630  if(hasData()) \
631  { \
632  iterator i = begin(), ie = end(); \
633  for(; i != ie; ++i) \
634  *i op v; \
635  } \
636  return *this; \
637  }
638 
639  VIGRA_CHUNKED_ARRAY_VIEW_ASSIGN(=)
640  VIGRA_CHUNKED_ARRAY_VIEW_ASSIGN(+=)
641  VIGRA_CHUNKED_ARRAY_VIEW_ASSIGN(-=)
642  VIGRA_CHUNKED_ARRAY_VIEW_ASSIGN(*=)
643  VIGRA_CHUNKED_ARRAY_VIEW_ASSIGN(/=)
644 
645  #undef VIGRA_CHUNKED_ARRAY_VIEW_ASSIGN
646 
647  // template<class Expression>
648  // MultiArrayView & operator=(multi_math::MultiMathOperand<Expression> const & rhs)
649  // {
650  // multi_math::math_detail::assign(*this, rhs);
651  // return *this;
652  // }
653 
654  // /** Add-assignment of an array expression. Fails with
655  // <tt>PreconditionViolation</tt> exception when the shapes do not match.
656  // */
657  // template<class Expression>
658  // MultiArrayView & operator+=(multi_math::MultiMathOperand<Expression> const & rhs)
659  // {
660  // multi_math::math_detail::plusAssign(*this, rhs);
661  // return *this;
662  // }
663 
664  // /** Subtract-assignment of an array expression. Fails with
665  // <tt>PreconditionViolation</tt> exception when the shapes do not match.
666  // */
667  // template<class Expression>
668  // MultiArrayView & operator-=(multi_math::MultiMathOperand<Expression> const & rhs)
669  // {
670  // multi_math::math_detail::minusAssign(*this, rhs);
671  // return *this;
672  // }
673 
674  // /** Multiply-assignment of an array expression. Fails with
675  // <tt>PreconditionViolation</tt> exception when the shapes do not match.
676  // */
677  // template<class Expression>
678  // MultiArrayView & operator*=(multi_math::MultiMathOperand<Expression> const & rhs)
679  // {
680  // multi_math::math_detail::multiplyAssign(*this, rhs);
681  // return *this;
682  // }
683 
684  // /** Divide-assignment of an array expression. Fails with
685  // <tt>PreconditionViolation</tt> exception when the shapes do not match.
686  // */
687  // template<class Expression>
688  // MultiArrayView & operator/=(multi_math::MultiMathOperand<Expression> const & rhs)
689  // {
690  // multi_math::math_detail::divideAssign(*this, rhs);
691  // return *this;
692  // }
693 
694  reference operator[](shape_type point)
695  {
696  VIGRA_ASSERT_INSIDE(point);
697  point += offset_;
698  Chunk * chunk = chunks_.data() +
699  detail::ChunkIndexing<N>::chunkOffset(point, bits_, chunks_.stride());
700  return *(chunk->pointer_ +
701  detail::ChunkIndexing<N>::offsetInChunk(point, mask_, chunk->strides_));
702  }
703 
704  const_reference operator[](shape_type const & point) const
705  {
706  return const_cast<MultiArrayView *>(this)->operator[](point);
707  }
708 
709  template <int M>
710  MultiArrayView <N-M, T, ChunkedArrayTag>
711  operator[](const TinyVector<MultiArrayIndex, M> &d) const
712  {
713  return bindInner(d);
714  }
715 
716  reference operator[](difference_type_1 d)
717  {
718  return operator[](scanOrderIndexToCoordinate(d));
719  }
720 
721  const_reference operator[](difference_type_1 d) const
722  {
723  return operator[](scanOrderIndexToCoordinate(d));
724  }
725 
726  difference_type scanOrderIndexToCoordinate(difference_type_1 d) const
727  {
728  difference_type coord(SkipInitialization);
729  detail::ScanOrderToCoordinate<actual_dimension>::exec(d, this->shape_, coord);
730  return coord;
731  }
732 
733  /** convert coordinate to scan-order index.
734  */
735  difference_type_1 coordinateToScanOrderIndex(const difference_type &d) const
736  {
737  return detail::CoordinateToScanOrder<actual_dimension>::exec(this->shape_, d);
738  }
739 
740  // /** 1D array access. Use only if N == 1.
741  // */
742  // reference operator() (difference_type_1 x)
743  // {
744  // VIGRA_ASSERT_INSIDE(difference_type(x));
745  // return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x)];
746  // }
747 
748  // /** 2D array access. Use only if N == 2.
749  // */
750  // reference operator() (difference_type_1 x, difference_type_1 y)
751  // {
752  // VIGRA_ASSERT_INSIDE(difference_type(x, y));
753  // return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x, y)];
754  // }
755 
756  // /** 3D array access. Use only if N == 3.
757  // */
758  // reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z)
759  // {
760  // VIGRA_ASSERT_INSIDE(difference_type(x, y, z));
761  // return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
762  // }
763 
764  // /** 4D array access. Use only if N == 4.
765  // */
766  // reference operator() (difference_type_1 x, difference_type_1 y,
767  // difference_type_1 z, difference_type_1 u)
768  // {
769  // VIGRA_ASSERT_INSIDE(difference_type(x, y, z, u));
770  // return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
771  // }
772 
773  // /** 5D array access. Use only if N == 5.
774  // */
775  // reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
776  // difference_type_1 u, difference_type_1 v)
777  // {
778  // VIGRA_ASSERT_INSIDE(difference_type(x, y,z, u,v));
779  // return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
780  // }
781 
782  // /** 1D const array access. Use only if N == 1.
783  // */
784  // const_reference operator() (difference_type_1 x) const
785  // {
786  // VIGRA_ASSERT_INSIDE(difference_type(x));
787  // return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x)];
788  // }
789 
790  // /** 2D const array access. Use only if N == 2.
791  // */
792  // const_reference operator() (difference_type_1 x, difference_type_1 y) const
793  // {
794  // VIGRA_ASSERT_INSIDE(difference_type(x, y));
795  // return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x, y)];
796  // }
797 
798  // /** 3D const array access. Use only if N == 3.
799  // */
800  // const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z) const
801  // {
802  // VIGRA_ASSERT_INSIDE(difference_type(x,y,z));
803  // return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
804  // }
805 
806  // /** 4D const array access. Use only if N == 4.
807  // */
808  // const_reference operator() (difference_type_1 x, difference_type_1 y,
809  // difference_type_1 z, difference_type_1 u) const
810  // {
811  // VIGRA_ASSERT_INSIDE(difference_type(x,y,z,u));
812  // return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
813  // }
814 
815  // /** 5D const array access. Use only if N == 5.
816  // */
817  // const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
818  // difference_type_1 u, difference_type_1 v) const
819  // {
820  // VIGRA_ASSERT_INSIDE(difference_type(x,y,z,u,v));
821  // return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
822  // }
823 
824  template <class U>
825  MultiArrayView & init(const U & init)
826  {
827  return operator=(init);
828  }
829 
830  template <class U, class CN>
831  void copy(const MultiArrayView <N, U, CN>& rhs)
832  {
833  operator=(rhs);
834  }
835 
836  template <class T2, class C2>
837  void swapData(MultiArrayView <N, T2, C2> rhs)
838  {
839  if(this == &rhs)
840  return;
841  vigra_precondition(this->shape() == rhs.shape(),
842  "MultiArrayView::swapData(): shape mismatch.");
843  iterator i = begin(), ie = end();
844  typename MultiArrayView<N, T2, C2>::iterator j = rhs.begin();
845  for(; i != ie; ++i, ++j)
846  std::swap(*i, *j);
847  }
848 
849  bool isUnstrided(unsigned int dimension = N-1) const
850  {
851  if(chunks_.size() > 1)
852  return false;
853  difference_type s = vigra::detail::defaultStride<actual_dimension>(this->shape());
854  for(unsigned int k = 0; k <= dimension; ++k)
855  if(chunks_.data()->strides_[k] != s[k])
856  return false;
857  return true;
858  }
859 
860  MultiArrayView<N-1, value_type, ChunkedArrayTag>
861  bindAt(MultiArrayIndex m, MultiArrayIndex d) const
862  {
863  MultiArrayView<N-1, value_type, ChunkedArrayTag> res(this->shape_.dropIndex(m), this->chunk_shape_.dropIndex(m));
864  res.offset_ = offset_.dropIndex(m);
865  res.bits_ = bits_.dropIndex(m);
866  res.mask_ = mask_.dropIndex(m);
867  res.chunks_.reshape(chunks_.shape().dropIndex(m));
868  res.unref_ = unref_;
869 
870  typedef std::size_t UI;
871  UI start = offset_[m] + d;
872  UI chunk_start = start >> bits_[m];
873  UI startInChunk = start - chunk_start * this->chunk_shape_[m];
874 
875  MultiArrayView<N-1, Chunk> view(chunks_.bindAt(m, chunk_start));
876  MultiCoordinateIterator<N-1> i(view.shape()),
877  end(i.getEndIterator());
878  for(; i != end; ++i)
879  {
880  res.chunks_[*i].pointer_ = view[*i].pointer_ + startInChunk*view[*i].strides_[m];
881  res.chunks_[*i].strides_ = view[*i].strides_.dropIndex(m);
882  }
883 
884  return res;
885  }
886 
887  template <unsigned int M>
888  MultiArrayView <N-1, value_type, ChunkedArrayTag>
889  bind (difference_type_1 d) const
890  {
891  return bindAt(M, d);
892  }
893 
894  MultiArrayView <N-1, value_type, ChunkedArrayTag>
895  bindOuter (difference_type_1 d) const
896  {
897  return bindAt(N-1, d);
898  }
899 
900  template <int M, class Index>
901  MultiArrayView <N-M, value_type, ChunkedArrayTag>
902  bindOuter(const TinyVector <Index, M> &d) const
903  {
904  return bindAt(N-1, d[M-1]).bindOuter(d.dropIndex(M-1));
905  }
906 
907  template <class Index>
908  MultiArrayView <N-1, value_type, ChunkedArrayTag>
909  bindOuter(const TinyVector <Index, 1> &d) const
910  {
911  return bindAt(N-1, d[0]);
912  }
913 
914  MultiArrayView <N-1, value_type, ChunkedArrayTag>
915  bindInner (difference_type_1 d) const
916  {
917  return bindAt(0, d);
918  }
919 
920  template <int M, class Index>
921  MultiArrayView <N-M, value_type, ChunkedArrayTag>
922  bindInner(const TinyVector <Index, M> &d) const
923  {
924  return bindAt(0, d[0]).bindInner(d.dropIndex(0));
925  }
926 
927  template <class Index>
928  MultiArrayView <N-1, value_type, ChunkedArrayTag>
929  bindInner(const TinyVector <Index, 1> &d) const
930  {
931  return bindAt(0, d[0]);
932  }
933 
934  // MultiArrayView <N, typename ExpandElementResult<T>::type, StridedArrayTag>
935  // bindElementChannel(difference_type_1 i) const
936  // {
937  // vigra_precondition(0 <= i && i < ExpandElementResult<T>::size,
938  // "MultiArrayView::bindElementChannel(i): 'i' out of range.");
939  // return expandElements(0).bindInner(i);
940  // }
941 
942  // MultiArrayView <N+1, typename ExpandElementResult<T>::type, StridedArrayTag>
943  // expandElements(difference_type_1 d) const;
944 
945  // MultiArrayView <N+1, T, StrideTag>
946  // insertSingletonDimension (difference_type_1 i) const;
947 
948  // MultiArrayView<N, Multiband<value_type>, StrideTag> multiband() const
949  // {
950  // return MultiArrayView<N, Multiband<value_type>, StrideTag>(*this);
951  // }
952 
953  // MultiArrayView<1, T, StridedArrayTag> diagonal() const
954  // {
955  // return MultiArrayView<1, T, StridedArrayTag>(Shape1(vigra::min(m_shape)),
956  // Shape1(vigra::sum(m_stride)), m_ptr);
957  // }
958 
959  inline void
960  checkSubarrayBounds(shape_type const & start, shape_type const & stop,
961  std::string message) const
962  {
963  message += ": subarray out of bounds.";
964  vigra_precondition(allLessEqual(shape_type(), start) &&
965  allLess(start, stop) &&
966  allLessEqual(stop, this->shape_),
967  message);
968  }
969 
970  MultiArrayView<N, value_type, ChunkedArrayTag>
971  subarray(shape_type start, shape_type stop)
972  {
973  checkSubarrayBounds(start, stop, "MultiArrayView<N-1, T, ChunkedArrayTag>::subarray()");
974  start += offset_;
975  stop += offset_;
976  shape_type chunk_start(chunkStart(start));
977 
978  MultiArrayView<N, value_type, ChunkedArrayTag> view(stop-start, this->chunk_shape_);
979  view.chunks_ = chunks_.subarray(chunk_start, chunkStop(stop));
980  view.offset_ = start - chunk_start * this->chunk_shape_;
981  view.bits_ = bits_;
982  view.mask_ = mask_;
983  view.unref_ = unref_;
984  return view;
985  }
986 
987  // /** apply an additional striding to the image, thereby reducing
988  // the shape of the array.
989  // for example, multiplying the stride of dimension one by three
990  // turns an appropriately laid out (interleaved) rgb image into
991  // a single band image.
992  // */
993  // MultiArrayView <N, T, StridedArrayTag>
994  // stridearray (const difference_type &s) const
995  // {
996  // difference_type shape = m_shape;
997  // for (unsigned int i = 0; i < actual_dimension; ++i)
998  // shape [i] /= s [i];
999  // return MultiArrayView <N, T, StridedArrayTag>(shape, m_stride * s, m_ptr);
1000  // }
1001 
1002  MultiArrayView <N, value_type, ChunkedArrayTag>
1003  transpose () const
1004  {
1005  return transpose(difference_type::linearSequence(N-1, -1));
1006  }
1007 
1008  MultiArrayView <N, value_type, ChunkedArrayTag>
1009  transpose(const difference_type &permutation) const
1010  {
1011  MultiArrayView<N, value_type, ChunkedArrayTag>
1012  view(vigra::transpose(this->shape_, permutation), vigra::transpose(this->chunk_shape_, permutation));
1013  view.chunks_ = chunks_.transpose(permutation); // also checks if permutation is valid
1014  view.offset_ = vigra::transpose(offset_, permutation);
1015  view.bits_ = vigra::transpose(bits_, permutation);
1016  view.mask_ = vigra::transpose(mask_, permutation);
1017  view.unref_ = unref_;
1018  typename MultiArray<N, Chunk>::iterator i = view.chunks_.begin(),
1019  iend = view.chunks_.end();
1020  for(; i != iend; ++i)
1021  i->strides_ = vigra::transpose(i->strides_, permutation);
1022  return view;
1023  }
1024 
1025  // MultiArrayView <N, T, StridedArrayTag>
1026  // permuteDimensions (const difference_type &s) const;
1027 
1028  // /** Permute the dimensions of the array so that the strides are in ascending order.
1029  // Determines the appropriate permutation and then calls permuteDimensions().
1030  // */
1031  // MultiArrayView <N, T, StridedArrayTag>
1032  // permuteStridesAscending() const;
1033 
1034  // /** Permute the dimensions of the array so that the strides are in descending order.
1035  // Determines the appropriate permutation and then calls permuteDimensions().
1036  // */
1037  // MultiArrayView <N, T, StridedArrayTag>
1038  // permuteStridesDescending() const;
1039 
1040  // /** Compute the ordering of the strides in this array.
1041  // The result is describes the current permutation of the axes relative
1042  // to the standard ascending stride order.
1043  // */
1044  // difference_type strideOrdering() const
1045  // {
1046  // return strideOrdering(m_stride);
1047  // }
1048 
1049  // /** Compute the ordering of the given strides.
1050  // The result is describes the current permutation of the axes relative
1051  // to the standard ascending stride order.
1052  // */
1053  // static difference_type strideOrdering(difference_type strides);
1054 
1055  template <class U, class C1>
1056  bool operator==(MultiArrayView<N, U, C1> const & rhs) const
1057  {
1058  if(this->shape() != rhs.shape())
1059  return false;
1060  const_iterator i = begin(), ie = end();
1061  typename MultiArrayView<N, U, C1>::const_iterator j = rhs.begin();
1062  for(; i != ie; ++i, ++j)
1063  if(*i != *j)
1064  return false;
1065  return true;
1066  }
1067 
1068  template <class U, class C1>
1069  bool operator!=(MultiArrayView<N, U, C1> const & rhs) const
1070  {
1071  return !operator==(rhs);
1072  }
1073 
1074  // bool all() const
1075  // {
1076  // bool res = true;
1077  // detail::reduceOverMultiArray(traverser_begin(), shape(),
1078  // res,
1079  // detail::AllTrueReduceFunctor(),
1080  // MetaInt<actual_dimension-1>());
1081  // return res;
1082  // }
1083 
1084  // bool any() const
1085  // {
1086  // bool res = false;
1087  // detail::reduceOverMultiArray(traverser_begin(), shape(),
1088  // res,
1089  // detail::AnyTrueReduceFunctor(),
1090  // MetaInt<actual_dimension-1>());
1091  // return res;
1092  // }
1093 
1094  // void minmax(T * minimum, T * maximum) const
1095  // {
1096  // std::pair<T, T> res(NumericTraits<T>::max(), NumericTraits<T>::min());
1097  // detail::reduceOverMultiArray(traverser_begin(), shape(),
1098  // res,
1099  // detail::MinmaxReduceFunctor(),
1100  // MetaInt<actual_dimension-1>());
1101  // *minimum = res.first;
1102  // *maximum = res.second;
1103  // }
1104 
1105  // template <class U>
1106  // void meanVariance(U * mean, U * variance) const
1107  // {
1108  // typedef typename NumericTraits<U>::RealPromote R;
1109  // R zero = R();
1110  // triple<double, R, R> res(0.0, zero, zero);
1111  // detail::reduceOverMultiArray(traverser_begin(), shape(),
1112  // res,
1113  // detail::MeanVarianceReduceFunctor(),
1114  // MetaInt<actual_dimension-1>());
1115  // *mean = res.second;
1116  // *variance = res.third / res.first;
1117  // }
1118 
1119  // template <class U>
1120  // U sum() const
1121  // {
1122  // U res = NumericTraits<U>::zero();
1123  // detail::reduceOverMultiArray(traverser_begin(), shape(),
1124  // res,
1125  // detail::SumReduceFunctor(),
1126  // MetaInt<actual_dimension-1>());
1127  // return res;
1128  // }
1129 
1130  // template <class U, class S>
1131  // void sum(MultiArrayView<N, U, S> sums) const
1132  // {
1133  // transformMultiArray(srcMultiArrayRange(*this),
1134  // destMultiArrayRange(sums),
1135  // FindSum<U>());
1136  // }
1137 
1138  // template <class U>
1139  // U product() const
1140  // {
1141  // U res = NumericTraits<U>::one();
1142  // detail::reduceOverMultiArray(traverser_begin(), shape(),
1143  // res,
1144  // detail::ProdReduceFunctor(),
1145  // MetaInt<actual_dimension-1>());
1146  // return res;
1147  // }
1148 
1149  // typename NormTraits<MultiArrayView>::SquaredNormType
1150  // squaredNorm() const
1151  // {
1152  // typedef typename NormTraits<MultiArrayView>::SquaredNormType SquaredNormType;
1153  // SquaredNormType res = NumericTraits<SquaredNormType>::zero();
1154  // detail::reduceOverMultiArray(traverser_begin(), shape(),
1155  // res,
1156  // detail::SquaredL2NormReduceFunctor(),
1157  // MetaInt<actual_dimension-1>());
1158  // return res;
1159  // }
1160 
1161  // typename NormTraits<MultiArrayView>::NormType
1162  // norm(int type = 2, bool useSquaredNorm = true) const;
1163 
1164  bool hasData () const
1165  {
1166  return chunks_.hasData();
1167  }
1168 
1169  iterator begin()
1170  {
1171  return createCoupledIterator(*this);
1172  }
1173 
1174  iterator end()
1175  {
1176  return begin().getEndIterator();
1177  }
1178 
1179  const_iterator cbegin() const
1180  {
1181  return createCoupledIterator(const_cast<MultiArrayView const &>(*this));
1182  }
1183 
1184  const_iterator cend() const
1185  {
1186  return cbegin().getEndIterator();
1187  }
1188 
1189  const_iterator begin() const
1190  {
1191  return createCoupledIterator(*this);
1192  }
1193 
1194  const_iterator end() const
1195  {
1196  return begin().getEndIterator();
1197  }
1198 
1199  chunk_iterator chunk_begin(shape_type const & start, shape_type const & stop)
1200  {
1201  checkSubarrayBounds(start, stop, "MultiArrayView<N-1, T, ChunkedArrayTag>::chunk_begin()");
1202  return chunk_iterator(this, start, stop, chunkStart(start), chunkStop(stop), this->chunk_shape_);
1203  }
1204 
1205  chunk_iterator chunk_end(shape_type const & start, shape_type const & stop)
1206  {
1207  return chunk_begin(start, stop).getEndIterator();
1208  }
1209 
1210  chunk_const_iterator chunk_begin(shape_type const & start, shape_type const & stop) const
1211  {
1212  checkSubarrayBounds(start, stop, "MultiArrayView<N-1, T, ChunkedArrayTag>::chunk_begin()");
1213  return chunk_const_iterator(this, start, stop, chunkStart(start), chunkStop(stop), this->chunk_shape_);
1214  }
1215 
1216  chunk_const_iterator chunk_end(shape_type const & start, shape_type const & stop) const
1217  {
1218  return chunk_begin(start, stop).getEndIterator();
1219  }
1220 
1221  chunk_const_iterator chunk_cbegin(shape_type const & start, shape_type const & stop) const
1222  {
1223  checkSubarrayBounds(start, stop, "MultiArrayView<N-1, T, ChunkedArrayTag>::chunk_cbegin()");
1224  return chunk_const_iterator(this, start, stop, chunkStart(start), chunkStop(stop), this->chunk_shape_);
1225  }
1226 
1227  chunk_const_iterator chunk_cend(shape_type const & start, shape_type const & stop) const
1228  {
1229  return chunk_cbegin(start, stop).getEndIterator();
1230  }
1231 
1232  view_type view ()
1233  {
1234  return *this;
1235  }
1236 
1237  MultiArray<N, Chunk> chunks_;
1238  shape_type offset_, bits_, mask_;
1239  VIGRA_SHARED_PTR<ChunkUnrefProxyBase> unref_;
1240 };
1241 
1242 template <unsigned int N, class T>
1243 typename MultiArrayView<N, T, ChunkedArrayTag>::iterator
1244 createCoupledIterator(MultiArrayView<N, T, ChunkedArrayTag> & m)
1245 {
1246  typedef typename MultiArrayView<N, T, ChunkedArrayTag>::iterator IteratorType;
1247  typedef typename IteratorType::handle_type P1;
1248  typedef typename P1::base_type P0;
1249 
1250  return IteratorType(P1(m,
1251  P0(m.shape())));
1252 }
1253 
1254 template <unsigned int N, class T>
1255 typename MultiArrayView<N, T, ChunkedArrayTag>::const_iterator
1256 createCoupledIterator(MultiArrayView<N, T, ChunkedArrayTag> const & m)
1257 {
1258  typedef typename MultiArrayView<N, T, ChunkedArrayTag>::const_iterator IteratorType;
1259  typedef typename IteratorType::handle_type P1;
1260  typedef typename P1::base_type P0;
1261 
1262  return IteratorType(P1(m,
1263  P0(m.shape())));
1264 }
1265 
1266 class ChunkedArrayOptions
1267 {
1268  public:
1269  ChunkedArrayOptions()
1270  : fill_value(0.0)
1271  , cache_max(-1)
1272  , compression_method(DEFAULT_COMPRESSION)
1273  {}
1274 
1275  ChunkedArrayOptions & fillValue(double v)
1276  {
1277  fill_value = v;
1278  return *this;
1279  }
1280 
1281  ChunkedArrayOptions fillValue(double v) const
1282  {
1283  return ChunkedArrayOptions(*this).fillValue(v);
1284  }
1285 
1286  ChunkedArrayOptions & cacheMax(int v)
1287  {
1288  cache_max = v;
1289  return *this;
1290  }
1291 
1292  ChunkedArrayOptions cacheMax(int v) const
1293  {
1294  return ChunkedArrayOptions(*this).cacheMax(v);
1295  }
1296 
1297  ChunkedArrayOptions & compression(CompressionMethod v)
1298  {
1299  compression_method = v;
1300  return *this;
1301  }
1302 
1303  ChunkedArrayOptions compression(CompressionMethod v) const
1304  {
1305  return ChunkedArrayOptions(*this).compression(v);
1306  }
1307 
1308  double fill_value;
1309  int cache_max;
1310  CompressionMethod compression_method;
1311 };
1312 
1313 /*
1314 The present implementation uses a memory-mapped sparse file to store the chunks.
1315 A sparse file is created on Linux using the O_TRUNC flag (this seems to be
1316 the default file behavior on Linux anyway), and on Windows by
1317 calling DeviceIoControl(file_handle, FSCTL_SET_SPARSE,...) after file creation.
1318 
1319 We can automatically delete the file upon closing. On Windows, this happens
1320 if the file was opened with FILE_FLAG_DELETE_ON_CLOSE. (It may be useful to
1321 combine this with the flag FILE_ATTRIBUTE_TEMPORARY, which tells the OS to
1322 avoid writing the file to disk if possible. However, judging from the timings,
1323 something is still written, or cleanup takes considerable time.)
1324 On Linux, you can call fileno(tmpfile()) for the same purpose.
1325 
1326 Alternatives are:
1327 * Keep the array in memory, but compress unused chunks.
1328 * Don't create a file explicitly, but use the swap file instead. This is
1329  achieved on Linux by mmap(..., MAP_PRIVATE | MAP_ANONYMOUS, -1, ...),
1330  on Windows by calling CreateFileMapping(INVALID_HANDLE_VALUE, ...).
1331  * On Linux, the memory must not be unmapped because this
1332  looses the data. In fact, anonymous mmap() is very similar to
1333  malloc(), and there is probably no good reason to use anonymous mmap().
1334  * On Windows, this is much faster, because the OS will never try to
1335  actually write something to disk (unless swapping is necessary).
1336 * Back chunks by HDF5 chunks, possibly using on-the-fly compression. This
1337  is in particular useful for existing HDF5 files.
1338 * Back chunks by HDF5 datasets. This can be combined with compression
1339  (both explicit and on-the-fly) or with memory mapping (using the
1340  function H5Dget_offset() to get the offset from the beginning of the file).
1341 
1342 FIXME:
1343 * backends:
1344  * allocators are not used
1345  * HDF5 only works for scalar types so far
1346  * HDF5 must support read-only and read/write mode
1347  * temp file arrays in swap (just an API addition to the constructor)
1348  * support TIFF chunked reading
1349 * the array implementations should go into cxx files in src/impex
1350  * this requires implementation of the low-level functions independently of dtype
1351  (use 'char *' and multiply shape and stride with sizeof(T))
1352  * don't forget to increment the soversion after the change
1353  * alternative: provide 'config_local.hxx' with flags for available packages
1354 * decide on chunk locking policies for array views (in particular, for index access)
1355  * array view has functions fetch()/release() (better names?) to lock/unlock
1356  _all_ chunks in the view
1357  * release() is automatically called in the destructor
1358  * it should be possible to call fetch in the constructor via a flag,
1359  but should the constructor fetch by default?
1360  * how should fetch() handle the case when the cache is too small
1361  * throw an exception?
1362  * silently enlarge the cache?
1363  * temporarily enlarge the cache?
1364  * provide an option to control the behavior?
1365  * also provide copySubarray() with ReadOnly and ReadWrite flags, where
1366  ReadWrite copies the subarray back in the destructor or on demand
1367  * locking is only required while each slice is copied
1368  * the copy functions can use normal array views and iterators
1369  * the ReadWrite version can store a checksum for each chunk (or part
1370  of a chunk) to detect collisions on write
1371  * use shared pointers to support memory management of the subarrays?
1372 * find efficient ways to support slicing and transposition in the indexing
1373  functions of a view.
1374  1. possibility: each view contains
1375  * an index object 'bound_index_' with original dimension whose values denote
1376  coordinates of bound axes and offsets for unbound coordinates
1377  * a permutation object 'permutation_' with dimension of the view that maps
1378  view coordinates to original coordinates
1379  * that is:
1380  operator[](index)
1381  {
1382  shape_type full_index(bound_index_);
1383  for(int k=0; k<N_view; ++k)
1384  full_index[permutation_[k]] += index[k];
1385  split full_index into chunk part and local part
1386  look up chunk
1387  return pixel
1388  }
1389  * maybe this is faster if it is combined with the stride computation?
1390  * an optimization for unsliced arrays is desirable
1391  2. possibility:
1392  * add the point offset to the low-dimensional index
1393  * split low-dimensional index into chunk part and local part
1394  * look up chunk
1395  * determine scalar pointer offset from local part and strides plus a
1396  chunk-specific correction that can be stored in a 3^N array
1397  - but can we efficiently determine where to look in that offset array?
1398  3. possibility:
1399  * don't care about speed - require copySubarray() if indexing should
1400  be fast
1401 * provide a ChunkIterator that iterates over all chunks in a given ROI and returns a
1402  MultiArrayView for the present chunk (which remains locked in cache until the
1403  iterator is advanced).
1404 * implement proper copy constructors and assignment for all backends
1405 * test HDF5 constructor from existing dataset
1406 * put HDF5 into header of its own
1407 * is the full chunkForIterator() function slow? Check this with a simplified one
1408  in a ChunkedArrayLazy where all chunlks are already implemented, so that
1409  we can simply can skip the check
1410 * add support for Multiband and TinyVector pixels
1411 
1412 */
1413 template <unsigned int N, class T>
1414 class ChunkedArray
1415 : public ChunkedArrayBase<N, T>
1416 {
1417  public:
1418  typedef ChunkedArrayBase<N, T> base_type;
1419  typedef typename MultiArrayShape<N>::type shape_type;
1420  typedef typename shape_type::value_type difference_type_1;
1421  typedef T value_type;
1422  typedef value_type * pointer;
1423  typedef value_type const * const_pointer;
1424  typedef value_type & reference;
1425  typedef value_type const & const_reference;
1426  typedef ChunkIterator<N, T> chunk_iterator;
1427  typedef ChunkIterator<N, T const> chunk_const_iterator;
1428  typedef StridedScanOrderIterator<N, ChunkedMemory<T>, reference, pointer> iterator;
1429  typedef StridedScanOrderIterator<N, ChunkedMemory<T const>, const_reference, const_pointer> const_iterator;
1430  typedef SharedChunkHandle<N, T> Handle;
1431  typedef ChunkBase<N, T> Chunk;
1432  typedef MultiArrayView<N, T, ChunkedArrayTag> view_type;
1433  typedef MultiArrayView<N, T const, ChunkedArrayTag> const_view_type;
1434  typedef std::queue<Handle*> CacheType;
1435 
1436  static const long chunk_asleep = Handle::chunk_asleep;
1437  static const long chunk_uninitialized = Handle::chunk_uninitialized;
1438  static const long chunk_locked = Handle::chunk_locked;
1439  static const long chunk_failed = Handle::chunk_failed;
1440 
1441  explicit ChunkedArray(shape_type const & shape,
1442  shape_type const & chunk_shape = shape_type(),
1443  ChunkedArrayOptions const & options = ChunkedArrayOptions())
1444  : ChunkedArrayBase<N, T>(shape, chunk_shape)
1445  , bits_(initBitMask(this->chunk_shape_))
1446  , mask_(this->chunk_shape_ -shape_type(1))
1447  , cache_max_size_(options.cache_max)
1448  , chunk_lock_(new threading::mutex())
1449  , fill_value_(T(options.fill_value))
1450  , fill_scalar_(options.fill_value)
1451  , handle_array_(detail::computeChunkArrayShape(shape, bits_, mask_))
1452  , data_bytes_()
1453  , overhead_bytes_(handle_array_.size()*sizeof(Handle))
1454  {
1455  fill_value_chunk_.pointer_ = &fill_value_;
1456  fill_value_handle_.pointer_ = &fill_value_chunk_;
1457  fill_value_handle_.chunk_state_.store(1);
1458  }
1459 
1460  static shape_type initBitMask(shape_type const & chunk_shape)
1461  {
1462  shape_type res;
1463  for(unsigned int k=0; k<N; ++k)
1464  {
1465  UInt32 bits = log2i(chunk_shape[k]);
1466  vigra_precondition(chunk_shape[k] == MultiArrayIndex(1 << bits),
1467  "ChunkedArray: chunk_shape elements must be powers of 2.");
1468  res[k] = bits;
1469  }
1470  return res;
1471  }
1472 
1473  virtual ~ChunkedArray()
1474  {
1475  // std::cerr << " final cache size: " << cacheSize() << " (max: " << cacheMaxSize() << ")\n";
1476  }
1477 
1478  int cacheSize() const
1479  {
1480  return cache_.size();
1481  }
1482 
1483  std::size_t dataBytes() const
1484  {
1485  return data_bytes_;
1486  }
1487 
1488  virtual shape_type chunkArrayShape() const
1489  {
1490  return handle_array_.shape();
1491  }
1492 
1493  virtual std::size_t dataBytes(Chunk * c) const = 0;
1494 
1495  std::size_t dataBytesPerChunk() const
1496  {
1497  return prod(this->chunk_shape_)*sizeof(T);
1498  }
1499 
1500  std::size_t overheadBytes() const
1501  {
1502  return overhead_bytes_;
1503  }
1504 
1505  virtual std::size_t overheadBytesPerChunk() const = 0;
1506 
1507  shape_type chunkStart(shape_type const & global_start) const
1508  {
1509  shape_type chunk_start(SkipInitialization);
1510  detail::ChunkIndexing<N>::chunkIndex(global_start, bits_, chunk_start);
1511  return chunk_start;
1512  }
1513 
1514  shape_type chunkStop(shape_type global_stop) const
1515  {
1516  global_stop -= shape_type(1);
1517  shape_type chunk_stop(SkipInitialization);
1518  detail::ChunkIndexing<N>::chunkIndex(global_stop, bits_, chunk_stop);
1519  chunk_stop += shape_type(1);
1520  return chunk_stop;
1521  }
1522 
1523  shape_type chunkShape(shape_type const & chunk_index) const
1524  {
1525  return min(this->chunk_shape_,
1526  this->shape_ - chunk_index*this->chunk_shape_);
1527  }
1528 
1529  using base_type::chunkShape;
1530 
1531  inline void
1532  checkSubarrayBounds(shape_type const & start, shape_type const & stop,
1533  std::string message) const
1534  {
1535  message += ": subarray out of bounds.";
1536  vigra_precondition(allLessEqual(shape_type(), start) &&
1537  allLess(start, stop) &&
1538  allLessEqual(stop, this->shape_),
1539  message);
1540  }
1541 
1542  template <class U, class C1>
1543  bool operator==(MultiArrayView<N, U, C1> const & rhs) const
1544  {
1545  if(this->shape() != rhs.shape())
1546  return false;
1547  const_iterator i = begin(), ie = end();
1548  typename MultiArrayView<N, U, C1>::const_iterator j = rhs.begin();
1549  for(; i != ie; ++i, ++j)
1550  if(*i != *j)
1551  return false;
1552  return true;
1553  }
1554 
1555  template <class U, class C1>
1556  bool operator!=(MultiArrayView<N, U, C1> const & rhs) const
1557  {
1558  return !operator==(rhs);
1559  }
1560 
1561  virtual pointer loadChunk(Chunk ** chunk, shape_type const & chunk_index) = 0;
1562 
1563  virtual bool unloadHandle(Handle * handle, bool destroy = false)
1564  {
1565  if(handle == &fill_value_handle_)
1566  return false;
1567  return unloadChunk(handle->pointer_, destroy);
1568  }
1569 
1570  virtual bool unloadChunk(Chunk * chunk, bool destroy = false) = 0;
1571 
1572  Handle * lookupHandle(shape_type const & index)
1573  {
1574  return &handle_array_[index];
1575  }
1576 
1577  virtual void unrefChunk(IteratorChunkHandle<N, T> * h) const
1578  {
1579  unrefChunk(h->chunk_);
1580  h->chunk_ = 0;
1581  }
1582 
1583  void unrefChunk(Handle * chunk) const
1584  {
1585  if(chunk)
1586  {
1587  long rc = chunk->chunk_state_.fetch_sub(1);
1588  #ifdef VIGRA_CHECK_BOUNDS
1589  vigra_invariant(rc >= 0,
1590  "ChunkedArray::unrefChunk(): chunk refcount got negative!");
1591  #endif
1592  }
1593  }
1594 
1595  void unrefChunks(ArrayVector<Handle*> const & chunks)
1596  {
1597  for(unsigned int k=0; k<chunks.size(); ++k)
1598  unrefChunk(chunks[k]);
1599 
1600  if(cacheMaxSize() > 0)
1601  {
1602  threading::lock_guard<threading::mutex> guard(*chunk_lock_);
1603  cleanCache(cache_.size());
1604  }
1605  }
1606 
1607  long acquireRef(Handle * handle) const
1608  {
1609  // Obtain a reference to the current chunk handle.
1610  // We use a simple spin-lock here because it is very fast in case of success
1611  // and failures (i.e. collisions with another thread) are presumably
1612  // very rare.
1613  //
1614  // the function returns the old value of chunk_state_
1615  long rc = handle->chunk_state_.load(threading::memory_order_acquire);
1616  while(true)
1617  {
1618  if(rc >= 0)
1619  {
1620  if(handle->chunk_state_.compare_exchange_weak(rc, rc+1, threading::memory_order_seq_cst))
1621  {
1622  return rc;
1623  }
1624  }
1625  else
1626  {
1627  if(rc == chunk_failed)
1628  {
1629  vigra_precondition(false,
1630  "ChunkedArray::acquireRef() attempt to access failed chunk.");
1631  }
1632  else if(rc == chunk_locked)
1633  {
1634  // cache management in progress => try again later
1635  threading::this_thread::yield();
1636  rc = handle->chunk_state_.load(threading::memory_order_acquire);
1637  }
1638  else if(handle->chunk_state_.compare_exchange_weak(rc, chunk_locked, threading::memory_order_seq_cst))
1639  {
1640  return rc;
1641  }
1642  }
1643  }
1644  }
1645 
1646  pointer getChunk(Handle * handle, bool isConst, bool insertInCache, shape_type const & chunk_index) const
1647  {
1648  ChunkedArray * self = const_cast<ChunkedArray *>(this);
1649 
1650  long rc = acquireRef(handle);
1651  if(rc >= 0)
1652  return handle->pointer_->pointer_;
1653 
1654  threading::lock_guard<threading::mutex> guard(*chunk_lock_);
1655  try
1656  {
1657  T * p = self->loadChunk(&handle->pointer_, chunk_index);
1658  Chunk * chunk = handle->pointer_;
1659  if(!isConst && rc == chunk_uninitialized)
1660  std::fill(p, p + prod(chunkShape(chunk_index)), this->fill_value_);
1661 
1662  self->data_bytes_ += dataBytes(chunk);
1663 
1664  if(cacheMaxSize() > 0 && insertInCache)
1665  {
1666  // insert in queue of mapped chunks
1667  self->cache_.push(handle);
1668 
1669  // do cache management if cache is full
1670  // (note that we still hold the chunk_lock_)
1671  self->cleanCache(2);
1672  }
1673  handle->chunk_state_.store(1, threading::memory_order_release);
1674  return p;
1675  }
1676  catch(...)
1677  {
1678  handle->chunk_state_.store(chunk_failed);
1679  throw;
1680  }
1681  }
1682 
1683  inline pointer
1684  chunkForIteratorImpl(shape_type const & point,
1685  shape_type & strides, shape_type & upper_bound,
1686  IteratorChunkHandle<N, T> * h,
1687  bool isConst) const
1688  {
1689  ChunkedArray * self = const_cast<ChunkedArray *>(this);
1690 
1691  unrefChunk(h->chunk_);
1692  h->chunk_ = 0;
1693 
1694  shape_type global_point = point + h->offset_;
1695 
1696  if(!this->isInside(global_point))
1697  {
1698  upper_bound = point + this->chunk_shape_;
1699  return 0;
1700  }
1701 
1702  shape_type chunkIndex(chunkStart(global_point));
1703 
1704  bool insertInCache = true;
1705  Handle * handle = self->lookupHandle(chunkIndex);
1706  if(isConst && handle->chunk_state_.load() == chunk_uninitialized)
1707  {
1708  handle = &self->fill_value_handle_;
1709  insertInCache = false;
1710  }
1711 
1712  pointer p = getChunk(handle, isConst, insertInCache, chunkIndex);
1713  strides = handle->strides();
1714  upper_bound = (chunkIndex + shape_type(1)) * this->chunk_shape_ - h->offset_;
1715  std::size_t offset = detail::ChunkIndexing<N>::offsetInChunk(global_point, mask_, strides);
1716  h->chunk_ = handle;
1717  return p + offset;
1718  }
1719 
1720  virtual pointer chunkForIterator(shape_type const & point,
1721  shape_type & strides, shape_type & upper_bound,
1722  IteratorChunkHandle<N, T> * h)
1723  {
1724  return chunkForIteratorImpl(point, strides, upper_bound, h, false);
1725  }
1726 
1727  virtual pointer chunkForIterator(shape_type const & point,
1728  shape_type & strides, shape_type & upper_bound,
1729  IteratorChunkHandle<N, T> * h) const
1730  {
1731  return chunkForIteratorImpl(point, strides, upper_bound, h, true);
1732  }
1733 
1734  // NOTE: This function must only be called while we hold the chunk_lock_.
1735  // This implies refcount != chunk_locked, so that race conditions are avoided.
1736  long releaseChunk(Handle * handle, bool destroy = false)
1737  {
1738  long rc = 0;
1739  bool mayUnload = handle->chunk_state_.compare_exchange_strong(rc, chunk_locked);
1740  if(!mayUnload && destroy)
1741  {
1742  rc = chunk_asleep;
1743  mayUnload = handle->chunk_state_.compare_exchange_strong(rc, chunk_locked);
1744  }
1745  if(mayUnload)
1746  {
1747  // refcount was zero or chunk_asleep => can unload
1748  try
1749  {
1750  vigra_invariant(handle != &fill_value_handle_,
1751  "ChunkedArray::releaseChunk(): attempt to release fill_value_handle_.");
1752  Chunk * chunk = handle->pointer_;
1753  this->data_bytes_ -= dataBytes(chunk);
1754  int didDestroy = unloadChunk(chunk, destroy);
1755  this->data_bytes_ += dataBytes(chunk);
1756  if(didDestroy)
1757  handle->chunk_state_.store(chunk_uninitialized);
1758  else
1759  handle->chunk_state_.store(chunk_asleep);
1760  }
1761  catch(...)
1762  {
1763  handle->chunk_state_.store(chunk_failed);
1764  throw;
1765  }
1766  }
1767  return rc;
1768  }
1769 
1770  // NOTE: this function must only be called while we hold the chunk_lock_
1771  void cleanCache(int how_many = -1)
1772  {
1773  if(how_many == -1)
1774  how_many = cache_.size();
1775  for(; cache_.size() > cacheMaxSize() && how_many > 0; --how_many)
1776  {
1777  Handle * handle = cache_.front();
1778  cache_.pop();
1779  long rc = releaseChunk(handle);
1780  if(rc > 0) // refcount was positive => chunk is still needed
1781  cache_.push(handle);
1782  }
1783  }
1784 
1785  // Sends all chunks asleep which are completely inside the given ROI.
1786  // If destroy == true and the backend supports destruction (currently:
1787  // ChunkedArrayLazy and ChunkedArrayCompressed), chunks will be deleted
1788  // entirely. The chunk's contents after releaseChunks() are undefined.
1789  // Currently, chunks retain their values when sent asleep, and assume the
1790  // array's fill_value when deleted, but applications should not rely on this
1791  // behavior.
1792  void releaseChunks(shape_type const & start, shape_type const & stop, bool destroy = false)
1793  {
1794  checkSubarrayBounds(start, stop, "ChunkedArray::releaseChunks()");
1795 
1796  MultiCoordinateIterator<N> i(chunkStart(start), chunkStop(stop)),
1797  end(i.getEndIterator());
1798  for(; i != end; ++i)
1799  {
1800  shape_type chunkOffset = *i * this->chunk_shape_;
1801  if(!allLessEqual(start, chunkOffset) ||
1802  !allLessEqual(min(chunkOffset+this->chunk_shape_, this->shape()), stop))
1803  {
1804  // chunk is only partially covered by the ROI
1805  continue;
1806  }
1807 
1808  Handle * handle = this->lookupHandle(*i);
1809  threading::lock_guard<threading::mutex> guard(*chunk_lock_);
1810  releaseChunk(handle, destroy);
1811  }
1812 
1813  // remove all chunks from the cache that are asleep or unitialized
1814  threading::lock_guard<threading::mutex> guard(*chunk_lock_);
1815  int cache_size = cache_.size();
1816  for(int k=0; k < cache_size; ++k)
1817  {
1818  Handle * handle = cache_.front();
1819  cache_.pop();
1820  if(handle->chunk_state_.load() >= 0)
1821  cache_.push(handle);
1822  }
1823  }
1824 
1825  template <class U, class Stride>
1826  void
1827  checkoutSubarray(shape_type const & start,
1828  MultiArrayView<N, U, Stride> & subarray) const
1829  {
1830  shape_type stop = start + subarray.shape();
1831 
1832  checkSubarrayBounds(start, stop, "ChunkedArray::checkoutSubarray()");
1833 
1834  chunk_const_iterator i = chunk_cbegin(start, stop);
1835  for(; i.isValid(); ++i)
1836  {
1837  subarray.subarray(i.chunkStart()-start, i.chunkStop()-start) = *i;
1838  }
1839  }
1840 
1841  template <class U, class Stride>
1842  void
1843  commitSubarray(shape_type const & start,
1844  MultiArrayView<N, U, Stride> const & subarray)
1845  {
1846  shape_type stop = start + subarray.shape();
1847 
1848  vigra_precondition(!this->isReadOnly(),
1849  "ChunkedArray::commitSubarray(): array is read-only.");
1850  checkSubarrayBounds(start, stop, "ChunkedArray::commitSubarray()");
1851 
1852  chunk_iterator i = chunk_begin(start, stop);
1853  for(; i.isValid(); ++i)
1854  {
1855  *i = subarray.subarray(i.chunkStart()-start, i.chunkStop()-start);
1856  }
1857  }
1858 
1859  template <class View>
1860  void subarrayImpl(shape_type const & start, shape_type const & stop,
1861  View & view,
1862  bool isConst) const
1863  {
1864  vigra_precondition(isConst || !this->isReadOnly(),
1865  "ChunkedArray::subarray(): array is read-only.");
1866  checkSubarrayBounds(start, stop, "ChunkedArray::subarray()");
1867  shape_type chunk_start(chunkStart(start)), chunk_stop(chunkStop(stop));
1868 
1869  view.shape_ = stop-start;
1870  view.chunk_shape_ = this->chunk_shape_;
1871  view.chunks_.reshape(chunk_stop-chunk_start);
1872  view.offset_ = start - chunk_start * this->chunk_shape_;
1873  view.bits_ = bits_;
1874  view.mask_ = mask_;
1875 
1876  typedef typename View::UnrefProxy Unref;
1877  ChunkedArray* self = const_cast<ChunkedArray*>(this);
1878  Unref * unref = new Unref(view.chunks_.size(), self);
1879  view.unref_ = VIGRA_SHARED_PTR<Unref>(unref);
1880 
1881  MultiCoordinateIterator<N> i(chunk_start, chunk_stop),
1882  end(i.getEndIterator());
1883  for(; i != end; ++i)
1884  {
1885  Handle * handle = self->lookupHandle(*i);
1886 
1887  if(isConst && handle->chunk_state_.load() == chunk_uninitialized)
1888  handle = &self->fill_value_handle_;
1889 
1890  // This potentially acquires the chunk_lock_ in each iteration.
1891  // Would it be better to acquire it once before the loop?
1892  pointer p = getChunk(handle, isConst, true, *i);
1893 
1894  ChunkBase<N, T> * mini_chunk = &view.chunks_[*i - chunk_start];
1895  mini_chunk->pointer_ = p;
1896  mini_chunk->strides_ = handle->strides();
1897  unref->chunks_[i.scanOrderIndex()] = handle;
1898  }
1899  }
1900 
1901  view_type
1902  subarray(shape_type const & start, shape_type const & stop)
1903  {
1904  view_type view;
1905  subarrayImpl(start, stop, view, false);
1906  return view;
1907  }
1908 
1909  const_view_type
1910  subarray(shape_type const & start, shape_type const & stop) const
1911  {
1912  const_view_type view;
1913  subarrayImpl(start, stop, view, true);
1914  return view;
1915  }
1916 
1917  const_view_type
1918  const_subarray(shape_type const & start, shape_type const & stop) const
1919  {
1920  const_view_type view;
1921  subarrayImpl(start, stop, view, true);
1922  return view;
1923  }
1924 
1925  value_type getItem(shape_type const & point) const
1926  {
1927  vigra_precondition(this->isInside(point),
1928  "ChunkedArray::getItem(): index out of bounds.");
1929 
1930  ChunkedArray * self = const_cast<ChunkedArray*>(this);
1931  shape_type chunk_index(chunkStart(point));
1932  Handle * handle = self->lookupHandle(chunk_index);
1933  if(handle->chunk_state_.load() == chunk_uninitialized)
1934  return fill_value_;
1935  pointer p = self->getChunk(handle, true, false, chunk_index);
1936  value_type res = *(p +
1937  detail::ChunkIndexing<N>::offsetInChunk(point, mask_, handle->strides()));
1938  self->unrefChunk(handle);
1939  return res;
1940  }
1941 
1942  void setItem(shape_type const & point, value_type const & v)
1943  {
1944  vigra_precondition(!this->isReadOnly(),
1945  "ChunkedArray::setItem(): array is read-only.");
1946  vigra_precondition(this->isInside(point),
1947  "ChunkedArray::setItem(): index out of bounds.");
1948 
1949  shape_type chunk_index(chunkStart(point));
1950  Handle * handle = lookupHandle(chunk_index);
1951  pointer p = getChunk(handle, false, false, chunk_index);
1952  *(p + detail::ChunkIndexing<N>::offsetInChunk(point, mask_, handle->strides())) = v;
1953  unrefChunk(handle);
1954  }
1955 
1956  MultiArrayView<N-1, T, ChunkedArrayTag>
1957  bindAt(MultiArrayIndex m, MultiArrayIndex d) const
1958  {
1959  shape_type start, stop(this->shape());
1960  start[m] = d;
1961  stop[m] = d+1;
1962  return subarray(start, stop).bindAt(m, 0);
1963  }
1964 
1965  template <unsigned int M>
1966  MultiArrayView <N-1, T, ChunkedArrayTag>
1967  bind (difference_type_1 d) const
1968  {
1969  return bindAt(M, d);
1970  }
1971 
1972  MultiArrayView <N-1, T, ChunkedArrayTag>
1973  bindOuter (difference_type_1 d) const
1974  {
1975  return bindAt(N-1, d);
1976  }
1977 
1978  template <int M, class Index>
1979  MultiArrayView <N-M, T, ChunkedArrayTag>
1980  bindOuter(const TinyVector <Index, M> &d) const
1981  {
1982  return bindAt(N-1, d[M-1]).bindOuter(d.dropIndex(M-1));
1983  }
1984 
1985  template <class Index>
1986  MultiArrayView <N-1, T, ChunkedArrayTag>
1987  bindOuter(const TinyVector <Index, 1> &d) const
1988  {
1989  return bindAt(N-1, d[0]);
1990  }
1991 
1992  MultiArrayView <N-1, T, ChunkedArrayTag>
1993  bindInner (difference_type_1 d) const
1994  {
1995  return bindAt(0, d);
1996  }
1997 
1998  template <int M, class Index>
1999  MultiArrayView <N-M, T, ChunkedArrayTag>
2000  bindInner(const TinyVector <Index, M> &d) const
2001  {
2002  return bindAt(0, d[0]).bindInner(d.dropIndex(0));
2003  }
2004 
2005  template <class Index>
2006  MultiArrayView <N-1, T, ChunkedArrayTag>
2007  bindInner(const TinyVector <Index, 1> &d) const
2008  {
2009  return bindAt(0, d[0]);
2010  }
2011 
2012  std::size_t cacheMaxSize() const
2013  {
2014  if(cache_max_size_ < 0)
2015  const_cast<int &>(cache_max_size_) = detail::defaultCacheSize(this->chunkArrayShape());
2016  return cache_max_size_;
2017  }
2018 
2019  void setCacheMaxSize(std::size_t c)
2020  {
2021  cache_max_size_ = c;
2022  if(c < cache_.size())
2023  {
2024  threading::lock_guard<threading::mutex> guard(*chunk_lock_);
2025  cleanCache();
2026  }
2027  }
2028 
2029  iterator begin()
2030  {
2031  return createCoupledIterator(*this);
2032  }
2033 
2034  iterator end()
2035  {
2036  return begin().getEndIterator();
2037  }
2038 
2039  const_iterator cbegin() const
2040  {
2041  return createCoupledIterator(const_cast<ChunkedArray const &>(*this));
2042  }
2043 
2044  const_iterator cend() const
2045  {
2046  return cbegin().getEndIterator();
2047  }
2048 
2049  const_iterator begin() const
2050  {
2051  return createCoupledIterator(*this);
2052  }
2053 
2054  const_iterator end() const
2055  {
2056  return begin().getEndIterator();
2057  }
2058 
2059  chunk_iterator chunk_begin(shape_type const & start, shape_type const & stop)
2060  {
2061  checkSubarrayBounds(start, stop, "ChunkedArray::chunk_begin()");
2062  return chunk_iterator(this, start, stop, chunkStart(start), chunkStop(stop), this->chunk_shape_);
2063  }
2064 
2065  chunk_iterator chunk_end(shape_type const & start, shape_type const & stop)
2066  {
2067  return chunk_begin(start, stop).getEndIterator();
2068  }
2069 
2070  chunk_const_iterator chunk_begin(shape_type const & start, shape_type const & stop) const
2071  {
2072  checkSubarrayBounds(start, stop, "ChunkedArray::chunk_begin()");
2073  return chunk_const_iterator(this, start, stop, chunkStart(start), chunkStop(stop), this->chunk_shape_);
2074  }
2075 
2076  chunk_const_iterator chunk_end(shape_type const & start, shape_type const & stop) const
2077  {
2078  return chunk_begin(start, stop).getEndIterator();
2079  }
2080 
2081  chunk_const_iterator chunk_cbegin(shape_type const & start, shape_type const & stop) const
2082  {
2083  checkSubarrayBounds(start, stop, "ChunkedArray::chunk_cbegin()");
2084  return chunk_const_iterator(this, start, stop, chunkStart(start), chunkStop(stop), this->chunk_shape_);
2085  }
2086 
2087  chunk_const_iterator chunk_cend(shape_type const & start, shape_type const & stop) const
2088  {
2089  return chunk_cbegin(start, stop).getEndIterator();
2090  }
2091 
2092  shape_type bits_, mask_;
2093  int cache_max_size_;
2094  VIGRA_SHARED_PTR<threading::mutex> chunk_lock_;
2095  CacheType cache_;
2096  Chunk fill_value_chunk_;
2097  Handle fill_value_handle_;
2098  value_type fill_value_;
2099  double fill_scalar_;
2100  MultiArray<N, Handle> handle_array_;
2101  std::size_t data_bytes_, overhead_bytes_;
2102 };
2103 
2104 /** Returns a CoupledScanOrderIterator to simultaneously iterate over image m1 and its coordinates.
2105  */
2106 template <unsigned int N, class T>
2107 typename ChunkedArray<N, T>::iterator
2108 createCoupledIterator(ChunkedArray<N, T> & m)
2109 {
2110  typedef typename ChunkedArray<N, T>::iterator IteratorType;
2111  typedef typename IteratorType::handle_type P1;
2112  typedef typename P1::base_type P0;
2113 
2114  return IteratorType(P1(m,
2115  P0(m.shape())));
2116 }
2117 
2118 template <unsigned int N, class T>
2119 typename ChunkedArray<N, T>::const_iterator
2120 createCoupledIterator(ChunkedArray<N, T> const & m)
2121 {
2122  typedef typename ChunkedArray<N, T>::const_iterator IteratorType;
2123  typedef typename IteratorType::handle_type P1;
2124  typedef typename P1::base_type P0;
2125 
2126  return IteratorType(P1(m,
2127  P0(m.shape())));
2128 }
2129 
2130 template <unsigned int N, class T, class Alloc = std::allocator<T> >
2131 class ChunkedArrayFull
2132 : public ChunkedArray<N, T>,
2133  public MultiArray<N, T, Alloc>
2134 {
2135  public:
2136 
2137  typedef MultiArray<N, T, Alloc> Storage;
2138  typedef typename Storage::value_type value_type;
2139  typedef typename Storage::pointer pointer;
2140  typedef typename Storage::const_pointer const_pointer;
2141  typedef typename Storage::reference reference;
2142  typedef typename Storage::const_reference const_reference;
2143  typedef typename Storage::difference_type difference_type;
2144  typedef typename Storage::difference_type shape_type;
2145  typedef typename Storage::key_type key_type;
2146  typedef typename Storage::size_type size_type;
2147  typedef typename Storage::difference_type_1 difference_type_1;
2148  typedef typename Storage::iterator iterator;
2149  typedef typename Storage::const_iterator const_iterator;
2150  typedef typename Storage::view_type view_type;
2151 
2152  typedef typename ChunkedArray<N, T>::Chunk Chunk;
2153 
2154  static shape_type computeChunkShape(shape_type s)
2155  {
2156  for(int k=0; k<N; ++k)
2157  s[k] = ceilPower2(s[k]);
2158  return s;
2159  }
2160 
2161  using Storage::subarray;
2162  using Storage::bindOuter;
2163  using Storage::bindInner;
2164  using Storage::bind;
2165  using Storage::bindAt;
2166  using Storage::isInside;
2167  using Storage::shape;
2168  using Storage::size;
2169  using Storage::begin;
2170  using Storage::end;
2171 
2172 #ifndef DOXYGEN // doxygen doesn't understand this
2173  using Storage::operator==;
2174  using Storage::operator!=;
2175 #endif
2176 
2177  explicit ChunkedArrayFull(shape_type const & shape,
2178  ChunkedArrayOptions const & options = ChunkedArrayOptions(),
2179  Alloc const & alloc = Alloc())
2180  : ChunkedArray<N, T>(shape, computeChunkShape(shape), options.cacheMax(0)),
2181  Storage(shape, this->fill_value_, alloc),
2182  upper_bound_(shape),
2183  chunk_(detail::defaultStride(shape), this->data())
2184  {
2185  this->handle_array_[0].pointer_ = &chunk_;
2186  this->handle_array_[0].chunk_state_.store(1);
2187  this->data_bytes_ = size()*sizeof(T);
2188  this->overhead_bytes_ = overheadBytesPerChunk();
2189  }
2190 
2191  ChunkedArrayFull(ChunkedArrayFull const & rhs)
2192  : ChunkedArray<N, T>(rhs),
2193  Storage(rhs),
2194  upper_bound_(rhs.upper_bound_),
2195  chunk_(detail::defaultStride(shape), this->data())
2196  {
2197  this->handle_array_[0].pointer_ = &chunk_;
2198  this->handle_array_[0].chunk_state_.store(1);
2199  }
2200 
2201  ChunkedArrayFull & operator=(ChunkedArrayFull const & rhs)
2202  {
2203  if(this != &rhs)
2204  {
2205  ChunkedArray<N, T>::operator=(rhs);
2206  Storage::operator=(rhs);
2207  upper_bound_ = rhs.upper_bound_;
2208  }
2209  return *this;
2210  }
2211 
2212  ~ChunkedArrayFull()
2213  {}
2214 
2215  virtual shape_type chunkArrayShape() const
2216  {
2217  return shape_type(1);
2218  }
2219 
2220  virtual pointer loadChunk(ChunkBase<N, T> **, shape_type const &)
2221  {
2222  return this->data();
2223  }
2224 
2225  virtual bool unloadChunk(ChunkBase<N, T> *, bool /* destroy */)
2226  {
2227  return false; // never destroys the data
2228  }
2229 
2230  virtual std::size_t dataBytes(Chunk * c) const
2231  {
2232  return prod(this->shape());
2233  }
2234 
2235  virtual std::size_t overheadBytesPerChunk() const
2236  {
2237  return sizeof(Chunk) + sizeof(SharedChunkHandle<N, T>);
2238  }
2239 
2240  virtual pointer chunkForIterator(shape_type const & point,
2241  shape_type & strides, shape_type & upper_bound,
2242  IteratorChunkHandle<N, T> * h) const
2243  {
2244  shape_type global_point = point + h->offset_;
2245 
2246  if(!this->isInside(global_point))
2247  {
2248  upper_bound = point + this->chunk_shape_;
2249  return 0;
2250  }
2251 
2252  strides = this->stride();
2253  upper_bound = upper_bound_;
2254  return const_cast<pointer>(&Storage::operator[](global_point));
2255  }
2256 
2257  virtual pointer chunkForIterator(shape_type const & point,
2258  shape_type & strides, shape_type & upper_bound,
2259  IteratorChunkHandle<N, T> * h)
2260  {
2261  shape_type global_point = point + h->offset_;
2262 
2263  if(!this->isInside(global_point))
2264  {
2265  upper_bound = point + this->chunk_shape_;
2266  return 0;
2267  }
2268 
2269  strides = this->stride();
2270  upper_bound = upper_bound_;
2271  return &Storage::operator[](global_point);
2272  }
2273 
2274  virtual std::string backend() const
2275  {
2276  return "ChunkedArrayFull";
2277  }
2278 
2279  shape_type upper_bound_;
2280  Chunk chunk_; // a dummy chunk to fulfill the API
2281 };
2282 
2283 template <unsigned int N, class T, class Alloc = std::allocator<T> >
2284 class ChunkedArrayLazy
2285 : public ChunkedArray<N, T>
2286 {
2287  public:
2288 
2289  class Chunk
2290  : public ChunkBase<N, T>
2291  {
2292  public:
2293  typedef typename MultiArrayShape<N>::type shape_type;
2294  typedef T value_type;
2295  typedef value_type * pointer;
2296  typedef value_type & reference;
2297 
2298  Chunk(shape_type const & shape, Alloc const & alloc = Alloc())
2299  : ChunkBase<N, T>(detail::defaultStride(shape))
2300  , size_(prod(shape))
2301  , alloc_(alloc)
2302  {}
2303 
2304  ~Chunk()
2305  {
2306  deallocate();
2307  }
2308 
2309  pointer allocate()
2310  {
2311  if(this->pointer_ == 0)
2312  this->pointer_ = detail::alloc_initialize_n<T>(size_, T(), alloc_);
2313  return this->pointer_;
2314  }
2315 
2316  void deallocate()
2317  {
2318  detail::destroy_dealloc_n(this->pointer_, size_, alloc_);
2319  this->pointer_ = 0;
2320  }
2321 
2322  MultiArrayIndex size_;
2323  Alloc alloc_;
2324 
2325  private:
2326  Chunk & operator=(Chunk const &);
2327  };
2328 
2329  typedef MultiArray<N, SharedChunkHandle<N, T> > ChunkStorage;
2330  typedef typename ChunkStorage::difference_type shape_type;
2331  typedef T value_type;
2332  typedef value_type * pointer;
2333  typedef value_type & reference;
2334 
2335  explicit ChunkedArrayLazy(shape_type const & shape,
2336  shape_type const & chunk_shape=shape_type(),
2337  ChunkedArrayOptions const & options = ChunkedArrayOptions(),
2338  Alloc const & alloc = Alloc())
2339  : ChunkedArray<N, T>(shape, chunk_shape, options.cacheMax(0))
2340  , alloc_(alloc)
2341  {}
2342 
2343  ~ChunkedArrayLazy()
2344  {
2345  typename ChunkStorage::iterator i = this->handle_array_.begin(),
2346  end = this->handle_array_.end();
2347  for(; i != end; ++i)
2348  {
2349  if(i->pointer_)
2350  delete static_cast<Chunk*>(i->pointer_);
2351  i->pointer_ = 0;
2352  }
2353  }
2354 
2355  virtual pointer loadChunk(ChunkBase<N, T> ** p, shape_type const & index)
2356  {
2357  if(*p == 0)
2358  {
2359  *p = new Chunk(this->chunkShape(index));
2360  this->overhead_bytes_ += sizeof(Chunk);
2361  }
2362  return static_cast<Chunk *>(*p)->allocate();
2363  }
2364 
2365  virtual bool unloadChunk(ChunkBase<N, T> * chunk, bool destroy)
2366  {
2367  if(destroy)
2368  static_cast<Chunk *>(chunk)->deallocate();
2369  return destroy;
2370  }
2371 
2372  virtual std::string backend() const
2373  {
2374  return "ChunkedArrayLazy";
2375  }
2376 
2377  virtual std::size_t dataBytes(ChunkBase<N,T> * c) const
2378  {
2379  return c->pointer_ == 0
2380  ? 0
2381  : static_cast<Chunk*>(c)->size_*sizeof(T);
2382  }
2383 
2384  virtual std::size_t overheadBytesPerChunk() const
2385  {
2386  return sizeof(Chunk) + sizeof(SharedChunkHandle<N, T>);
2387  }
2388 
2389  Alloc alloc_;
2390 };
2391 
2392 template <unsigned int N, class T, class Alloc = std::allocator<T> >
2393 class ChunkedArrayCompressed
2394 : public ChunkedArray<N, T>
2395 {
2396  public:
2397 
2398  class Chunk
2399  : public ChunkBase<N, T>
2400  {
2401  public:
2402  typedef typename MultiArrayShape<N>::type shape_type;
2403  typedef T value_type;
2404  typedef value_type * pointer;
2405  typedef value_type & reference;
2406 
2407  Chunk(shape_type const & shape)
2408  : ChunkBase<N, T>(detail::defaultStride(shape))
2409  , compressed_()
2410  , size_(prod(shape))
2411  {}
2412 
2413  ~Chunk()
2414  {
2415  deallocate();
2416  }
2417 
2418  pointer allocate()
2419  {
2420  if(this->pointer_ == 0)
2421  this->pointer_ = detail::alloc_initialize_n<T>(size_, T(), alloc_);
2422  return this->pointer_;
2423  }
2424 
2425  void deallocate()
2426  {
2427  detail::destroy_dealloc_n(this->pointer_, size_, alloc_);
2428  this->pointer_ = 0;
2429  compressed_.clear();
2430  }
2431 
2432  void compress(CompressionMethod method)
2433  {
2434  if(this->pointer_ != 0)
2435  {
2436  vigra_invariant(compressed_.size() == 0,
2437  "ChunkedArrayCompressed::Chunk::compress(): compressed and uncompressed pointer are both non-zero.");
2438 
2439  ::vigra::compress((char const *)this->pointer_, size_*sizeof(T), compressed_, method);
2440 
2441  // std::cerr << "compression ratio: " << double(compressed_.size())/(this->size()*sizeof(T)) << "\n";
2442  detail::destroy_dealloc_n(this->pointer_, size_, alloc_);
2443  this->pointer_ = 0;
2444  }
2445  }
2446 
2447  pointer uncompress(CompressionMethod method)
2448  {
2449  if(this->pointer_ == 0)
2450  {
2451  if(compressed_.size())
2452  {
2453  this->pointer_ = alloc_.allocate((typename Alloc::size_type)size_);
2454 
2455  ::vigra::uncompress(compressed_.data(), compressed_.size(),
2456  (char*)this->pointer_, size_*sizeof(T), method);
2457  compressed_.clear();
2458  }
2459  else
2460  {
2461  this->pointer_ = allocate();
2462  }
2463  }
2464  else
2465  {
2466  vigra_invariant(compressed_.size() == 0,
2467  "ChunkedArrayCompressed::Chunk::uncompress(): compressed and uncompressed pointer are both non-zero.");
2468  }
2469  return this->pointer_;
2470  }
2471 
2472  ArrayVector<char> compressed_;
2473  MultiArrayIndex size_;
2474  Alloc alloc_;
2475 
2476  private:
2477  Chunk & operator=(Chunk const &);
2478  };
2479 
2480  typedef MultiArray<N, SharedChunkHandle<N, T> > ChunkStorage;
2481  typedef typename ChunkStorage::difference_type shape_type;
2482  typedef T value_type;
2483  typedef value_type * pointer;
2484  typedef value_type & reference;
2485 
2486  explicit ChunkedArrayCompressed(shape_type const & shape,
2487  shape_type const & chunk_shape=shape_type(),
2488  ChunkedArrayOptions const & options = ChunkedArrayOptions())
2489  : ChunkedArray<N, T>(shape, chunk_shape, options),
2490  compression_method_(options.compression_method)
2491  {
2492  if(compression_method_ == DEFAULT_COMPRESSION)
2493  compression_method_ = LZ4;
2494  }
2495 
2496  ~ChunkedArrayCompressed()
2497  {
2498  typename ChunkStorage::iterator i = this->handle_array_.begin(),
2499  end = this->handle_array_.end();
2500  for(; i != end; ++i)
2501  {
2502  if(i->pointer_)
2503  delete static_cast<Chunk*>(i->pointer_);
2504  i->pointer_ = 0;
2505  }
2506  }
2507 
2508  virtual pointer loadChunk(ChunkBase<N, T> ** p, shape_type const & index)
2509  {
2510  if(*p == 0)
2511  {
2512  *p = new Chunk(this->chunkShape(index));
2513  this->overhead_bytes_ += sizeof(Chunk);
2514  }
2515  return static_cast<Chunk *>(*p)->uncompress(compression_method_);
2516  }
2517 
2518  virtual bool unloadChunk(ChunkBase<N, T> * chunk, bool destroy)
2519  {
2520  if(destroy)
2521  static_cast<Chunk *>(chunk)->deallocate();
2522  else
2523  static_cast<Chunk *>(chunk)->compress(compression_method_);
2524  return destroy;
2525  }
2526 
2527  virtual std::string backend() const
2528  {
2529  switch(compression_method_)
2530  {
2531  case ZLIB:
2532  return "ChunkedArrayCompressed<ZLIB>";
2533  case ZLIB_NONE:
2534  return "ChunkedArrayCompressed<ZLIB_NONE>";
2535  case ZLIB_FAST:
2536  return "ChunkedArrayCompressed<ZLIB_FAST>";
2537  case ZLIB_BEST:
2538  return "ChunkedArrayCompressed<ZLIB_BEST>";
2539  case LZ4:
2540  return "ChunkedArrayCompressed<LZ4>";
2541  default:
2542  return "unknown";
2543  }
2544  }
2545 
2546  virtual std::size_t dataBytes(ChunkBase<N,T> * c) const
2547  {
2548  return c->pointer_ == 0
2549  ? static_cast<Chunk*>(c)->compressed_.size()
2550  : static_cast<Chunk*>(c)->size_*sizeof(T);
2551  }
2552 
2553  virtual std::size_t overheadBytesPerChunk() const
2554  {
2555  return sizeof(Chunk) + sizeof(SharedChunkHandle<N, T>);
2556  }
2557 
2558  CompressionMethod compression_method_;
2559 };
2560 
2561 template <unsigned int N, class T>
2562 class ChunkedArrayTmpFile
2563 : public ChunkedArray<N, T>
2564 {
2565  public:
2566 #ifdef _WIN32
2567  typedef HANDLE FileHandle;
2568 #else
2569  typedef int FileHandle;
2570 #endif
2571 
2572  class Chunk
2573  : public ChunkBase<N, T>
2574  {
2575  public:
2576  typedef typename MultiArrayShape<N>::type shape_type;
2577  typedef T value_type;
2578  typedef value_type * pointer;
2579  typedef value_type & reference;
2580 
2581  Chunk(shape_type const & shape,
2582  std::size_t offset, size_t alloc_size,
2583  FileHandle file)
2584  : ChunkBase<N, T>(detail::defaultStride(shape))
2585  , offset_(offset)
2586  , alloc_size_(alloc_size)
2587  , file_(file)
2588  {}
2589 
2590  ~Chunk()
2591  {
2592  unmap();
2593  }
2594 
2595  pointer map()
2596  {
2597  if(this->pointer_ == 0)
2598  {
2599  #ifdef _WIN32
2600  static const std::size_t bits = sizeof(DWORD)*8,
2601  mask = (std::size_t(1) << bits) - 1;
2602  this->pointer_ = (pointer)MapViewOfFile(file_, FILE_MAP_ALL_ACCESS,
2603  std::size_t(offset_) >> bits, offset_ & mask, alloc_size_);
2604  if(this->pointer_ == 0)
2605  winErrorToException("ChunkedArrayChunk::map(): ");
2606  #else
2607  this->pointer_ = (pointer)mmap(0, alloc_size_, PROT_READ | PROT_WRITE, MAP_SHARED,
2608  file_, offset_);
2609  if(this->pointer_ == 0)
2610  throw std::runtime_error("ChunkedArrayChunk::map(): mmap() failed.");
2611  #endif
2612  }
2613  return this->pointer_;
2614  }
2615 
2616  void unmap()
2617  {
2618  if(this->pointer_ != 0)
2619  {
2620  #ifdef _WIN32
2621  ::UnmapViewOfFile(this->pointer_);
2622  #else
2623  munmap(this->pointer_, alloc_size_);
2624  #endif
2625  this->pointer_ = 0;
2626  }
2627  }
2628 
2629  std::size_t offset_, alloc_size_;
2630  FileHandle file_;
2631 
2632  private:
2633  Chunk & operator=(Chunk const &);
2634  };
2635 
2636  typedef MultiArray<N, SharedChunkHandle<N, T> > ChunkStorage;
2637  typedef MultiArray<N, std::size_t> OffsetStorage;
2638  typedef typename ChunkStorage::difference_type shape_type;
2639  typedef T value_type;
2640  typedef value_type * pointer;
2641  typedef value_type & reference;
2642 
2643  static std::size_t computeAllocSize(shape_type const & shape)
2644  {
2645  std::size_t size = prod(shape)*sizeof(T);
2646  std::size_t mask = mmap_alignment - 1;
2647  return (size + mask) & ~mask;
2648  }
2649 
2650  explicit ChunkedArrayTmpFile(shape_type const & shape,
2651  shape_type const & chunk_shape=shape_type(),
2652  ChunkedArrayOptions const & options = ChunkedArrayOptions(),
2653  std::string const & path = "")
2654  : ChunkedArray<N, T>(shape, chunk_shape, options)
2655  #ifndef VIGRA_NO_SPARSE_FILE
2656  , offset_array_(this->chunkArrayShape())
2657  #endif
2658  , file_size_()
2659  , file_capacity_()
2660  {
2661  #ifdef VIGRA_NO_SPARSE_FILE
2662  file_capacity_ = 4*prod(this->chunk_shape_)*sizeof(T);
2663  #else
2664  // compute offset in file
2665  typename OffsetStorage::iterator i = offset_array_.begin(),
2666  end = offset_array_.end();
2667  std::size_t size = 0;
2668  for(; i != end; ++i)
2669  {
2670  *i = size;
2671  size += computeAllocSize(this->chunkShape(i.point()));
2672  }
2673  file_capacity_ = size;
2674  this->overhead_bytes_ += offset_array_.size()*sizeof(std::size_t);
2675  // std::cerr << " file size: " << size << "\n";
2676  #endif
2677 
2678  #ifdef _WIN32
2679  // create a temp file
2680  file_ = ::CreateFile(winTempFileName(path).c_str(), GENERIC_READ | GENERIC_WRITE,
2681  0, NULL, CREATE_ALWAYS, FILE_ATTRIBUTE_TEMPORARY | FILE_FLAG_DELETE_ON_CLOSE, NULL);
2682  if (file_ == INVALID_HANDLE_VALUE)
2683  winErrorToException("ChunkedArrayTmpFile(): ");
2684 
2685  // make it a sparse file
2686  DWORD dwTemp;
2687  if(!::DeviceIoControl(file_, FSCTL_SET_SPARSE, NULL, 0, NULL, 0, &dwTemp, NULL))
2688  winErrorToException("ChunkedArrayTmpFile(): ");
2689 
2690  // place the data in the swap file
2691  // file_ = INVALID_HANDLE_VALUE;
2692 
2693  // resize and memory-map the file
2694  static const std::size_t bits = sizeof(LONG)*8, mask = (std::size_t(1) << bits) - 1;
2695  mappedFile_ = CreateFileMapping(file_, NULL, PAGE_READWRITE,
2696  file_capacity_ >> bits, file_capacity_ & mask, NULL);
2697  if(!mappedFile_)
2698  winErrorToException("ChunkedArrayTmpFile(): ");
2699  #else
2700  mappedFile_ = file_ = fileno(tmpfile());
2701  if(file_ == -1)
2702  throw std::runtime_error("ChunkedArrayTmpFile(): unable to open file.");
2703  lseek(file_, file_capacity_-1, SEEK_SET);
2704  if(write(file_, "0", 1) == -1)
2705  throw std::runtime_error("ChunkedArrayTmpFile(): unable to resize file.");
2706  #endif
2707  }
2708 
2709  ~ChunkedArrayTmpFile()
2710  {
2711  typename ChunkStorage::iterator i = this->handle_array_.begin(),
2712  end = this->handle_array_.end();
2713  for(; i != end; ++i)
2714  {
2715  if(i->pointer_)
2716  delete static_cast<Chunk*>(i->pointer_);
2717  i->pointer_ = 0;
2718  }
2719  #ifdef _WIN32
2720  ::CloseHandle(mappedFile_);
2721  ::CloseHandle(file_);
2722  #else
2723  ::close(file_);
2724  #endif
2725  }
2726 
2727  virtual pointer loadChunk(ChunkBase<N, T> ** p, shape_type const & index)
2728  {
2729  if(*p == 0)
2730  {
2731  shape_type shape = this->chunkShape(index);
2732  std::size_t chunk_size = computeAllocSize(shape);
2733  #ifdef VIGRA_NO_SPARSE_FILE
2734  std::size_t offset = file_size_;
2735  if(offset + chunk_size > file_capacity_)
2736  {
2737  file_capacity_ = max<std::size_t>(offset+chunk_size, file_capacity_ * 120 / 100); // extend file by 20%
2738  if(lseek(file_, file_capacity_-1, SEEK_SET) == -1)
2739  throw std::runtime_error("ChunkedArrayTmpFile(): unable to reset file size.");
2740  if(write(file_, "0", 1) == -1)
2741  throw std::runtime_error("ChunkedArrayTmpFile(): unable to resize file.");
2742  }
2743  file_size_ += chunk_size;
2744  #else
2745  std::size_t offset = offset_array_[index];
2746  #endif
2747  *p = new Chunk(shape, offset, chunk_size, mappedFile_);
2748  this->overhead_bytes_ += sizeof(Chunk);
2749  }
2750  return static_cast<Chunk*>(*p)->map();
2751  }
2752 
2753  virtual bool unloadChunk(ChunkBase<N, T> * chunk, bool /* destroy*/)
2754  {
2755  static_cast<Chunk *>(chunk)->unmap();
2756  return false; // never destroys the data
2757  }
2758 
2759  virtual std::string backend() const
2760  {
2761  return "ChunkedArrayTmpFile";
2762  }
2763 
2764  virtual std::size_t dataBytes(ChunkBase<N,T> * c) const
2765  {
2766  return c->pointer_ == 0
2767  ? 0
2768  : static_cast<Chunk*>(c)->alloc_size_;
2769  }
2770 
2771  virtual std::size_t overheadBytesPerChunk() const
2772  {
2773  #ifdef VIGRA_NO_SPARSE_FILE
2774  return sizeof(Chunk) + sizeof(SharedChunkHandle<N, T>);
2775  #else
2776  return sizeof(Chunk) + sizeof(SharedChunkHandle<N, T>) + sizeof(std::size_t);
2777  #endif
2778  }
2779 
2780  #ifndef VIGRA_NO_SPARSE_FILE
2781  OffsetStorage offset_array_; // the array of chunks
2782  #endif
2783  FileHandle file_, mappedFile_; // the file back-end
2784  std::size_t file_size_, file_capacity_;
2785 };
2786 
2787 template<unsigned int N, class U>
2788 class ChunkIterator
2789 : public MultiCoordinateIterator<N>
2790 , private MultiArrayView<N, typename UnqualifiedType<U>::type>
2791 {
2792  public:
2793  typedef typename UnqualifiedType<U>::type T;
2794  typedef MultiCoordinateIterator<N> base_type;
2795  typedef MultiArrayView<N, T> base_type2;
2796 
2797  typedef typename base_type::shape_type shape_type;
2798  typedef typename base_type::difference_type difference_type;
2799  typedef ChunkIterator iterator;
2800  typedef std::random_access_iterator_tag iterator_category;
2801 
2802  typedef MultiArrayView<N, T> value_type;
2803  typedef MultiArrayView<N, T> & reference;
2804  typedef MultiArrayView<N, T> const & const_reference;
2805  typedef MultiArrayView<N, T> * pointer;
2806  typedef MultiArrayView<N, T> const * const_pointer;
2807 
2808  typedef typename IfBool<UnqualifiedType<U>::isConst,
2809  ChunkedArrayBase<N, T> const,
2810  ChunkedArrayBase<N, T> >::type array_type;
2811  typedef IteratorChunkHandle<N, T> Chunk;
2812 
2813 
2814  ChunkIterator()
2815  : base_type()
2816  , base_type2()
2817  {}
2818 
2819  ChunkIterator(array_type * array,
2820  shape_type const & start, shape_type const & end,
2821  shape_type const & chunk_start, shape_type const & chunk_end,
2822  shape_type const & chunk_shape)
2823  : base_type(chunk_start, chunk_end)
2824  , array_(array)
2825  , chunk_(chunk_start * chunk_shape)
2826  , start_(start - chunk_.offset_)
2827  , stop_(end - chunk_.offset_)
2828  , chunk_shape_(chunk_shape)
2829  {
2830  getChunk();
2831  }
2832 
2833  ChunkIterator(ChunkIterator const & rhs)
2834  : base_type(rhs)
2835  , base_type2(rhs)
2836  , array_(rhs.array_)
2837  , chunk_(rhs.chunk_)
2838  , start_(rhs.start_)
2839  , stop_(rhs.stop_)
2840  , chunk_shape_(rhs.chunk_shape_)
2841  {
2842  getChunk();
2843  }
2844 
2845  ChunkIterator & operator=(ChunkIterator const & rhs)
2846  {
2847  if(this != &rhs)
2848  {
2849  base_type::operator=(rhs);
2850  array_ = rhs.array_;
2851  chunk_ = rhs.chunk_;
2852  start_ = rhs.start_;
2853  stop_ = rhs.stop_;
2854  chunk_shape_ = rhs.chunk_shape_;
2855  getChunk();
2856  }
2857  return *this;
2858  }
2859 
2860  reference operator*()
2861  {
2862  return *this;
2863  }
2864 
2865  const_reference operator*() const
2866  {
2867  return *this;
2868  }
2869 
2870  pointer operator->()
2871  {
2872  return this;
2873  }
2874 
2875  const_pointer operator->() const
2876  {
2877  return this;
2878  }
2879 
2880  value_type operator[](MultiArrayIndex i) const
2881  {
2882  return *(ChunkIterator(*this) += i);
2883  }
2884 
2885  value_type operator[](const shape_type &coordOffset) const
2886  {
2887  return *(ChunkIterator(*this) += coordOffset);
2888  }
2889 
2890  void getChunk()
2891  {
2892  if(array_)
2893  {
2894  shape_type array_point = max(start_, this->point()*chunk_shape_),
2895  upper_bound(SkipInitialization);
2896  this->m_ptr = array_->chunkForIterator(array_point, this->m_stride, upper_bound, &chunk_);
2897  this->m_shape = min(upper_bound, stop_) - array_point;
2898  }
2899  }
2900 
2901  shape_type chunkStart() const
2902  {
2903  return max(start_, this->point()*chunk_shape_) + chunk_.offset_;
2904  }
2905 
2906  shape_type chunkStop() const
2907  {
2908  return chunkStart() + this->m_shape;
2909  }
2910 
2911  ChunkIterator & operator++()
2912  {
2913  base_type::operator++();
2914  getChunk();
2915  return *this;
2916  }
2917 
2918  ChunkIterator operator++(int)
2919  {
2920  ChunkIterator res(*this);
2921  ++*this;
2922  return res;
2923  }
2924 
2925  ChunkIterator & operator+=(MultiArrayIndex i)
2926  {
2928  getChunk();
2929  return *this;
2930  }
2931 
2932  ChunkIterator & operator+=(const shape_type &coordOffset)
2933  {
2934  base_type::operator+=(coordOffset);
2935  getChunk();
2936  return *this;
2937  }
2938 
2939  ChunkIterator & operator--()
2940  {
2941  base_type::operator--();
2942  getChunk();
2943  return *this;
2944  }
2945 
2946  ChunkIterator operator--(int)
2947  {
2948  ChunkIterator res(*this);
2949  --*this;
2950  return res;
2951  }
2952 
2953  ChunkIterator & operator-=(MultiArrayIndex i)
2954  {
2955  return operator+=(-i);
2956  }
2957 
2958  ChunkIterator & operator-=(const shape_type &coordOffset)
2959  {
2960  return operator+=(-coordOffset);
2961  }
2962 
2963  ChunkIterator getEndIterator() const
2964  {
2965  ChunkIterator res(*this);
2966  static_cast<base_type &>(res) = base_type::getEndIterator();
2967  res.getChunk();
2968  return res;
2969  }
2970 
2971  ChunkIterator operator+(MultiArrayIndex d) const
2972  {
2973  return ChunkIterator(*this) += d;
2974  }
2975 
2976  ChunkIterator operator-(MultiArrayIndex d) const
2977  {
2978  return ChunkIterator(*this) -= d;
2979  }
2980 
2981  ChunkIterator operator+(const shape_type &coordOffset) const
2982  {
2983  return ChunkIterator(*this) += coordOffset;
2984  }
2985 
2986  ChunkIterator operator-(const shape_type &coordOffset) const
2987  {
2988  return ChunkIterator(*this) -= coordOffset;
2989  }
2990 
2991  MultiArrayIndex operator-(const ChunkIterator & other) const
2992  {
2993  return base_type::operator-(other);
2994  }
2995 
2996 #ifndef DOXYGEN // doxygen doesn't understand this
2997  using base_type::operator==;
2998  using base_type::operator!=;
2999 #endif
3000  using base_type::shape;
3001 
3002  array_type * array_;
3003  Chunk chunk_;
3004  shape_type start_, stop_, chunk_shape_, array_point_;
3005 };
3006 
3007 } // namespace vigra
3008 
3009 #undef VIGRA_ASSERT_INSIDE
3010 
3011 #endif /* VIGRA_MULTI_ARRAY_CHUNKED_HXX */

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.10.0 (Thu Jan 8 2015)