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

multi_array.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2003-2008 by Gunnar Kedenburg and Ullrich Koethe */
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 
37 #ifndef VIGRA_MULTI_ARRAY_HXX
38 #define VIGRA_MULTI_ARRAY_HXX
39 
40 #include <memory>
41 #include <algorithm>
42 #include "accessor.hxx"
43 #include "tinyvector.hxx"
44 #include "rgbvalue.hxx"
45 #include "basicimageview.hxx"
46 #include "imageiterator.hxx"
47 #include "numerictraits.hxx"
48 #include "multi_iterator.hxx"
49 #include "metaprogramming.hxx"
50 #include "mathutil.hxx"
51 
52 // Bounds checking Macro used if VIGRA_CHECK_BOUNDS is defined.
53 #ifdef VIGRA_CHECK_BOUNDS
54 #define VIGRA_ASSERT_INSIDE(diff) \
55  vigra_precondition(this->isInside(diff), "Index out of bounds")
56 #else
57 #define VIGRA_ASSERT_INSIDE(diff)
58 #endif
59 
60 namespace vigra
61 {
62 
63 namespace detail
64 {
65 /********************************************************/
66 /* */
67 /* defaultStride */
68 /* */
69 /********************************************************/
70 
71 /* generates the stride for a gapless shape.
72 
73  Namespace: vigra::detail
74 */
75 template <unsigned int N>
76 inline TinyVector <MultiArrayIndex, N>
77 defaultStride(const TinyVector <MultiArrayIndex, N> &shape)
78 {
79  TinyVector <MultiArrayIndex, N> ret;
80  ret [0] = 1;
81  for (int i = 1; i < (int)N; ++i)
82  ret [i] = ret [i-1] * shape [i-1];
83  return ret;
84 }
85 
86 /********************************************************/
87 /* */
88 /* ScanOrderToOffset */
89 /* */
90 /********************************************************/
91 
92 /* transforms an index in scan order sense to a pointer offset in a possibly
93  strided, multi-dimensional array.
94 
95  Namespace: vigra::detail
96 */
97 
98 template <int K>
99 struct ScanOrderToOffset
100 {
101  template <int N>
102  static MultiArrayIndex
103  exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> &shape,
104  const TinyVector <MultiArrayIndex, N> & stride)
105  {
106  return stride[N-K] * (d % shape[N-K]) +
107  ScanOrderToOffset<K-1>::exec(d / shape[N-K], shape, stride);
108  }
109 };
110 
111 template <>
112 struct ScanOrderToOffset<1>
113 {
114  template <int N>
115  static MultiArrayIndex
116  exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> & /*shape*/,
117  const TinyVector <MultiArrayIndex, N> & stride)
118  {
119  return stride[N-1] * d;
120  }
121 };
122 
123 template <int K>
124 struct ScanOrderToCoordinate
125 {
126  template <int N>
127  static void
128  exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> &shape,
129  TinyVector <MultiArrayIndex, N> & result)
130  {
131  result[N-K] = (d % shape[N-K]);
132  ScanOrderToCoordinate<K-1>::exec(d / shape[N-K], shape, result);
133  }
134 };
135 
136 template <>
137 struct ScanOrderToCoordinate<1>
138 {
139  template <int N>
140  static void
141  exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> & /*shape*/,
142  TinyVector <MultiArrayIndex, N> & result)
143  {
144  result[N-1] = d;
145  }
146 };
147 
148 template <int K>
149 struct CoordinateToScanOrder
150 {
151  template <int N>
152  static MultiArrayIndex
153  exec(const TinyVector <MultiArrayIndex, N> &shape,
154  const TinyVector <MultiArrayIndex, N> & coordinate)
155  {
156  return coordinate[N-K] + shape[N-K] * CoordinateToScanOrder<K-1>::exec(shape, coordinate);
157  }
158 };
159 
160 template <>
161 struct CoordinateToScanOrder<1>
162 {
163  template <int N>
164  static MultiArrayIndex
165  exec(const TinyVector <MultiArrayIndex, N> & /*shape*/,
166  const TinyVector <MultiArrayIndex, N> & coordinate)
167  {
168  return coordinate[N-1];
169  }
170 };
171 
172 
173 template <class C>
174 struct CoordinatesToOffest
175 {
176  template <int N>
177  static MultiArrayIndex
178  exec(const TinyVector <MultiArrayIndex, N> & stride, MultiArrayIndex x)
179  {
180  return stride[0] * x;
181  }
182  template <int N>
183  static MultiArrayIndex
184  exec(const TinyVector <MultiArrayIndex, N> & stride, MultiArrayIndex x, MultiArrayIndex y)
185  {
186  return stride[0] * x + stride[1] * y;
187  }
188 };
189 
190 template <>
191 struct CoordinatesToOffest<UnstridedArrayTag>
192 {
193  template <int N>
194  static MultiArrayIndex
195  exec(const TinyVector <MultiArrayIndex, N> & /*stride*/, MultiArrayIndex x)
196  {
197  return x;
198  }
199  template <int N>
200  static MultiArrayIndex
201  exec(const TinyVector <MultiArrayIndex, N> & stride, MultiArrayIndex x, MultiArrayIndex y)
202  {
203  return x + stride[1] * y;
204  }
205 };
206 
207 /********************************************************/
208 /* */
209 /* MaybeStrided */
210 /* */
211 /********************************************************/
212 
213 /* metatag implementing a test for marking MultiArrays that were
214  indexed at the zero'th dimension as strided, and all others as
215  unstrided.
216 
217 <b>\#include</b>
218 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
219 
220 Namespace: vigra::detail
221 */
222 template <unsigned int N>
223 struct MaybeStrided
224 {
225  typedef UnstridedArrayTag type;
226 };
227 
228 template <>
229 struct MaybeStrided <0>
230 {
231  typedef StridedArrayTag type;
232 };
233 
234 /********************************************************/
235 /* */
236 /* MultiIteratorChooser */
237 /* */
238 /********************************************************/
239 
240 /* metatag implementing a test (by pattern matching) for marking
241  MultiArrays that were indexed at the zero'th dimension as strided.
242 
243 <b>\#include</b>
244 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
245 
246 Namespace: vigra::detail
247 */
248 template <class O>
249 struct MultiIteratorChooser
250 {
251  struct Nil {};
252 
253  template <unsigned int N, class T, class REFERENCE, class POINTER>
254  struct Traverser
255  {
256  typedef Nil type;
257  };
258 };
259 
260 /********************************************************/
261 /* */
262 /* MultiIteratorChooser <StridedArrayTag> */
263 /* */
264 /********************************************************/
265 
266 /* specialization of the MultiIteratorChooser for strided arrays.
267 
268 <b>\#include</b>
269 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
270 
271 Namespace: vigra::detail
272 */
273 template <>
274 struct MultiIteratorChooser <StridedArrayTag>
275 {
276  template <unsigned int N, class T, class REFERENCE, class POINTER>
277  struct Traverser
278  {
279  typedef StridedMultiIterator <N, T, REFERENCE, POINTER> type;
280  };
281 };
282 
283 /********************************************************/
284 /* */
285 /* MultiIteratorChooser <UnstridedArrayTag> */
286 /* */
287 /********************************************************/
288 
289 /* specialization of the MultiIteratorChooser for unstrided arrays.
290 
291 <b>\#include</b>
292 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
293 
294 Namespace: vigra::detail
295 */
296 template <>
297 struct MultiIteratorChooser <UnstridedArrayTag>
298 {
299  template <unsigned int N, class T, class REFERENCE, class POINTER>
300  struct Traverser
301  {
302  typedef MultiIterator <N, T, REFERENCE, POINTER> type;
303  };
304 };
305 
306 /********************************************************/
307 /* */
308 /* helper functions */
309 /* */
310 /********************************************************/
311 
312 template <class DestIterator, class Shape, class T>
313 inline void
314 initMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<0>)
315 {
316  DestIterator dend = d + shape[0];
317  for(; d < dend; ++d)
318  {
319  *d = init;
320  }
321 }
322 
323 template <class DestIterator, class Shape, class T, int N>
324 void
325 initMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<N>)
326 {
327  DestIterator dend = d + shape[N];
328  for(; d < dend; ++d)
329  {
330  initMultiArrayData(d.begin(), shape, init, MetaInt<N-1>());
331  }
332 }
333 
334 #define VIGRA_COPY_MULTI_ARRAY_DATA(name, op) \
335 template <class SrcIterator, class Shape, class DestIterator> \
336 inline void \
337 name##MultiArrayData(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>) \
338 { \
339  SrcIterator send = s + shape[0]; \
340  for(; s < send; ++s, ++d) \
341  { \
342  *d op detail::RequiresExplicitCast<typename DestIterator::value_type>::cast(*s); \
343  } \
344 } \
345  \
346 template <class SrcIterator, class Shape, class DestIterator, int N> \
347 void \
348 name##MultiArrayData(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>) \
349 { \
350  SrcIterator send = s + shape[N]; \
351  for(; s < send; ++s, ++d) \
352  { \
353  name##MultiArrayData(s.begin(), shape, d.begin(), MetaInt<N-1>()); \
354  } \
355 } \
356 \
357 template <class DestIterator, class Shape, class T> \
358 inline void \
359 name##ScalarMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<0>) \
360 { \
361  DestIterator dend = d + shape[0]; \
362  for(; d < dend; ++d) \
363  { \
364  *d op detail::RequiresExplicitCast<typename DestIterator::value_type>::cast(init); \
365  } \
366 } \
367  \
368 template <class DestIterator, class Shape, class T, int N> \
369 void \
370 name##ScalarMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<N>) \
371 { \
372  DestIterator dend = d + shape[N]; \
373  for(; d < dend; ++d) \
374  { \
375  name##ScalarMultiArrayData(d.begin(), shape, init, MetaInt<N-1>()); \
376  } \
377 }
378 
379 VIGRA_COPY_MULTI_ARRAY_DATA(copy, =)
380 VIGRA_COPY_MULTI_ARRAY_DATA(copyAdd, +=)
381 VIGRA_COPY_MULTI_ARRAY_DATA(copySub, -=)
382 VIGRA_COPY_MULTI_ARRAY_DATA(copyMul, *=)
383 VIGRA_COPY_MULTI_ARRAY_DATA(copyDiv, /=)
384 
385 #undef VIGRA_COPY_MULTI_ARRAY_DATA
386 
387 template <class SrcIterator, class Shape, class T, class ALLOC>
388 inline void
389 uninitializedCopyMultiArrayData(SrcIterator s, Shape const & shape, T * & d, ALLOC & a, MetaInt<0>)
390 {
391  SrcIterator send = s + shape[0];
392  for(; s < send; ++s, ++d)
393  {
394  a.construct(d, static_cast<T const &>(*s));
395  }
396 }
397 
398 template <class SrcIterator, class Shape, class T, class ALLOC, int N>
399 void
400 uninitializedCopyMultiArrayData(SrcIterator s, Shape const & shape, T * & d, ALLOC & a, MetaInt<N>)
401 {
402  SrcIterator send = s + shape[N];
403  for(; s < send; ++s)
404  {
405  uninitializedCopyMultiArrayData(s.begin(), shape, d, a, MetaInt<N-1>());
406  }
407 }
408 
409 template <class SrcIterator, class Shape, class T>
410 inline void
411 normMaxOfMultiArray(SrcIterator s, Shape const & shape, T & result, MetaInt<0>)
412 {
413  SrcIterator send = s + shape[0];
414  for(; s < send; ++s)
415  {
416  T v = norm(*s);
417  if(result < v)
418  result = v;
419  }
420 }
421 
422 template <class SrcIterator, class Shape, class T, int N>
423 void
424 normMaxOfMultiArray(SrcIterator s, Shape const & shape, T & result, MetaInt<N>)
425 {
426  SrcIterator send = s + shape[N];
427  for(; s < send; ++s)
428  {
429  normMaxOfMultiArray(s.begin(), shape, result, MetaInt<N-1>());
430  }
431 }
432 
433 template <class T>
434 struct MultiArrayL1Functor
435 {
436  template <class U>
437  T operator()(U t) const
438  { return norm(t); }
439 };
440 
441 template <class T>
442 struct MultiArrayL2Functor
443 {
444  template <class U>
445  T operator()(U t) const
446  { return squaredNorm(t); }
447 };
448 
449 template <class T>
450 struct MultiArrayScaledL2Functor
451 {
452  T scale;
453 
454  MultiArrayScaledL2Functor(T s)
455  : scale(s)
456  {}
457 
458  template <class U>
459  T operator()(U t) const
460  { return squaredNorm(T(t) / scale); }
461 };
462 
463 template <class SrcIterator, class Shape, class Functor, class T>
464 inline void
465 sumOverMultiArray(SrcIterator s, Shape const & shape, Functor f, T & result, MetaInt<0>)
466 {
467  SrcIterator send = s + shape[0];
468  for(; s < send; ++s)
469  {
470  result += f(*s);
471  }
472 }
473 
474 template <class SrcIterator, class Shape, class Functor, class T, int N>
475 void
476 sumOverMultiArray(SrcIterator s, Shape const & shape, Functor f, T & result, MetaInt<N>)
477 {
478  SrcIterator send = s + shape[N];
479  for(; s < send; ++s)
480  {
481  sumOverMultiArray(s.begin(), shape, f, result, MetaInt<N-1>());
482  }
483 }
484 
485 template <class SrcIterator, class Shape, class DestIterator>
486 inline bool
487 equalityOfMultiArrays(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>)
488 {
489  SrcIterator send = s + shape[0];
490  for(; s < send; ++s, ++d)
491  {
492  if(!(*s == *d))
493  return false;
494  }
495  return true;
496 }
497 
498 template <class SrcIterator, class Shape, class DestIterator, int N>
499 bool
500 equalityOfMultiArrays(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>)
501 {
502  SrcIterator send = s + shape[N];
503  for(; s < send; ++s, ++d)
504  {
505  if(!equalityOfMultiArrays(s.begin(), shape, d.begin(), MetaInt<N-1>()))
506  return false;
507  }
508  return true;
509 }
510 
511 
512 template <class SrcIterator, class Shape, class DestIterator>
513 inline void
514 swapDataImpl(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>)
515 {
516  SrcIterator send = s + shape[0];
517  for(; s < send; ++s, ++d)
518  std::swap(*s, *d);
519 }
520 
521 template <class SrcIterator, class Shape, class DestIterator, int N>
522 void
523 swapDataImpl(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>)
524 {
525  SrcIterator send = s + shape[N];
526  for(; s < send; ++s, ++d)
527  swapDataImpl(s.begin(), shape, d.begin(), MetaInt<N-1>());
528 }
529 
530 } // namespace detail
531 
532 /********************************************************/
533 /* */
534 /* MultiArrayView */
535 /* */
536 /********************************************************/
537 
538 // forward declaration
539 template <unsigned int N, class T, class C = UnstridedArrayTag>
540 class MultiArrayView;
541 template <unsigned int N, class T, class A = std::allocator<T> >
542 class MultiArray;
543 
544 /********************************************************/
545 /* */
546 /* NormTraits */
547 /* */
548 /********************************************************/
549 
550 template <unsigned int N, class T, class C>
551 struct NormTraits<MultiArrayView<N, T, C> >
552 {
553  typedef MultiArrayView<N, T, C> Type;
554  typedef typename NormTraits<T>::SquaredNormType SquaredNormType;
555  typedef typename SquareRootTraits<SquaredNormType>::SquareRootResult NormType;
556 };
557 
558 template <unsigned int N, class T, class A>
559 struct NormTraits<MultiArray<N, T, A> >
560 : public NormTraits<MultiArrayView<N, T, UnstridedArrayTag> >
561 {
562  typedef NormTraits<MultiArrayView<N, T, UnstridedArrayTag> > BaseType;
563  typedef MultiArray<N, T, A> Type;
564  typedef typename BaseType::SquaredNormType SquaredNormType;
565  typedef typename BaseType::NormType NormType;
566 };
567 
568 /** \brief Base class for, and view to, \ref vigra::MultiArray.
569 
570 This class implements the interface of both MultiArray and
571 MultiArrayView. By default, MultiArrayViews are tagged as
572 unstrided. If necessary, strided arrays are constructed automatically
573 by calls to a variant of the bind...() function.
574 
575 If you want to apply an algorithm requiring an image to a
576 <tt>MultiArrayView</tt> of appropriate (2-dimensional) shape, you can
577 create a \ref vigra::BasicImageView that acts as a wrapper with the
578 necessary interface -- see \ref MultiArrayToImage.
579 
580 The template parameter are as follows
581 \code
582  N: the array dimension
583 
584  T: the type of the array elements
585 
586  C: a tag determining whether the array's inner dimension is strided
587  or not. An array is unstrided if the array elements occupy consecutive
588  memory location, strided if there is an offset in between (e.g.
589  when a view is created that skips every other array element).
590  The compiler can generate faster code for unstrided arrays.
591  Possible values: UnstridedArrayTag (default), StridedArrayTag
592 \endcode
593 
594 <b>\#include</b>
595 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
596 
597 Namespace: vigra
598 */
599 template <unsigned int N, class T, class C>
601 {
602 public:
603 
604  /** the array's actual dimensionality.
605  This ensures that MultiArrayView can also be used for
606  scalars (that is, when <tt>N == 0</tt>). Calculated as:<br>
607  \code
608  actual_dimension = (N==0) ? 1 : N
609  \endcode
610  */
611  enum ActualDimension { actual_dimension = (N==0) ? 1 : N };
612 
613  /** the array's value type
614  */
615  typedef T value_type;
616 
617  /** reference type (result of operator[])
618  */
619  typedef value_type &reference;
620 
621  /** const reference type (result of operator[] const)
622  */
623  typedef const value_type &const_reference;
624 
625  /** pointer type
626  */
627  typedef value_type *pointer;
628 
629  /** const pointer type
630  */
631  typedef const value_type *const_pointer;
632 
633  /** difference type (used for multi-dimensional offsets and indices)
634  */
636 
637  /** size type
638  */
639  typedef difference_type size_type;
640 
641  /** difference and index type for a single dimension
642  */
644 
645  /** traverser (MultiIterator) type
646  */
647  typedef typename vigra::detail::MultiIteratorChooser <
648  C>::template Traverser <actual_dimension, T, T &, T *>::type traverser;
649 
650  /** const traverser (MultiIterator) type
651  */
652  typedef typename vigra::detail::MultiIteratorChooser <
653  C>::template Traverser <actual_dimension, T, T const &, T const *>::type const_traverser;
654 
655  /** the view type associated with this array.
656  */
658 
659  /** the matrix type associated with this array.
660  */
662 
663 protected:
664 
665  typedef typename difference_type::value_type diff_zero_t;
666 
667  /** the shape of the image pointed to is stored here.
668  */
669  difference_type m_shape;
670 
671  /** the strides (offset of a sample to the next) for every dimension
672  are stored here.
673  */
674  difference_type m_stride;
675 
676  /** pointer to the image.
677  */
678  pointer m_ptr;
679 
680  template <class U, class CN>
681  void copyImpl(const MultiArrayView <N, U, CN>& rhs);
682 
683  template <class U, class CN>
684  void swapDataImpl(MultiArrayView <N, U, CN> rhs);
685 
686  template <class CN>
687  bool arraysOverlap(const MultiArrayView <N, T, CN>& rhs) const;
688 
689  template <class U, class CN>
690  bool arraysOverlap(const MultiArrayView <N, U, CN>&) const
691  {
692  return false;
693  }
694 
695 public:
696 
697  /** default constructor: create an invalid view,
698  * i.e. hasData() returns false and size() is zero.
699  */
701  : m_shape (diff_zero_t(0)), m_stride (diff_zero_t(0)), m_ptr (0)
702  {}
703 
704  /** construct from shape and pointer
705  */
706  MultiArrayView (const difference_type &shape, pointer ptr)
707  : m_shape (shape),
708  m_stride (detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape)),
709  m_ptr (ptr)
710  {}
711 
712  /** Construct from shape, strides (offset of a sample to the
713  next) for every dimension, and pointer. (Note that
714  strides are not given in bytes, but in offsets of the
715  respective pointer type.)
716  */
717  MultiArrayView (const difference_type &shape,
718  const difference_type &stride,
719  pointer ptr)
720  : m_shape (shape),
721  m_stride (stride),
722  m_ptr (ptr)
723  {}
724 
725 
726  /** Assignment. There are 3 cases:
727 
728  <ul>
729  <li> When this <tt>MultiArrayView</tt> does not point to valid data
730  (e.g. after default construction), it becomes a copy of \a rhs.
731  <li> When the shapes of the two arrays match, the array contents are copied.
732  <li> Otherwise, a <tt>PreconditionViolation</tt> exception is thrown.
733  </ul>
734  */
735  MultiArrayView & operator=(MultiArrayView const & rhs);
736 
737  /** Assignment of a differently typed MultiArrayView. Fails with
738  <tt>PreconditionViolation</tt> exception when the shapes do not match.
739  */
740  template<class U, class C1>
742  {
743  vigra_precondition(this->shape() == rhs.shape(),
744  "MultiArrayView::operator=() size mismatch.");
745  this->copyImpl(rhs);
746  return *this;
747  }
748 
749  /** Add-assignment of a compatible MultiArrayView. Fails with
750  <tt>PreconditionViolation</tt> exception when the shapes do not match.
751  */
752  template<class U, class C1>
754 
755  /** Subtract-assignment of a compatible MultiArrayView. Fails with
756  <tt>PreconditionViolation</tt> exception when the shapes do not match.
757  */
758  template<class U, class C1>
760 
761  /** Multiply-assignment of a compatible MultiArrayView. Fails with
762  <tt>PreconditionViolation</tt> exception when the shapes do not match.
763  */
764  template<class U, class C1>
766 
767  /** Divide-assignment of a compatible MultiArrayView. Fails with
768  <tt>PreconditionViolation</tt> exception when the shapes do not match.
769  */
770  template<class U, class C1>
772 
773  /** Add-assignment of a scalar.
774  */
775  MultiArrayView & operator+=(T const & rhs)
776  {
777  detail::copyAddScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
778  return *this;
779  }
780 
781  /** Subtract-assignment of a scalar.
782  */
783  MultiArrayView & operator-=(T const & rhs)
784  {
785  detail::copySubScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
786  return *this;
787  }
788 
789  /** Multiply-assignment of a scalar.
790  */
791  MultiArrayView & operator*=(T const & rhs)
792  {
793  detail::copyMulScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
794  return *this;
795  }
796 
797  /** Divide-assignment of a scalar.
798  */
799  MultiArrayView & operator/=(T const & rhs)
800  {
801  detail::copyDivScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
802  return *this;
803  }
804 
805  /** array access.
806  */
807  reference operator[] (const difference_type &d)
808  {
809  VIGRA_ASSERT_INSIDE(d);
810  return m_ptr [dot (d, m_stride)];
811  }
812 
813  /** array access.
814  */
815  const_reference operator[] (const difference_type &d) const
816  {
817  VIGRA_ASSERT_INSIDE(d);
818  return m_ptr [dot (d, m_stride)];
819  }
820 
821  /** equvalent to bindInner(), when M < N.
822  */
823  template <unsigned int M>
825  {
826  return bindInner(d);
827  }
828 
829  /** Array access in scan-order sense.
830  Mostly useful to support standard indexing for 1-dimensional multi-arrays,
831  but works for any N. Use scanOrderIndexToCoordinate() and
832  coordinateToScanOrderIndex() for conversion between indices and coordinates.
833 
834  <b>Note:</b> This function should not be used in the inner loop, because the
835  conversion of the scan order index into a memory address is expensive
836  (it must take into account that memory may not be consecutive for subarrays
837  and/or strided arrays). Always prefer operator() if possible.
838  */
839  reference operator[](difference_type_1 d)
840  {
841  VIGRA_ASSERT_INSIDE(scanOrderIndexToCoordinate(d));
842  return m_ptr [detail::ScanOrderToOffset<actual_dimension>::exec(d, m_shape, m_stride)];
843  }
844 
845  /** Array access in scan-order sense.
846  Mostly useful to support standard indexing for 1-dimensional multi-arrays,
847  but works for any N. Use scanOrderIndexToCoordinate() and
848  coordinateToScanOrderIndex() for conversion between indices and coordinates.
849 
850  <b>Note:</b> This function should not be used in the inner loop, because the
851  conversion of the scan order index into a memory address is expensive
852  (it must take into account that memory may not be consecutive for subarrays
853  and/or strided arrays). Always prefer operator() if possible.
854  */
855  const_reference operator[](difference_type_1 d) const
856  {
857  VIGRA_ASSERT_INSIDE(scanOrderIndexToCoordinate(d));
858  return m_ptr [detail::ScanOrderToOffset<actual_dimension>::exec(d, m_shape, m_stride)];
859  }
860 
861  /** convert scan-order index to coordinate.
862  */
863  difference_type scanOrderIndexToCoordinate(difference_type_1 d) const
864  {
865  difference_type result;
866  detail::ScanOrderToCoordinate<actual_dimension>::exec(d, m_shape, result);
867  return result;
868  }
869 
870  /** convert coordinate to scan-order index.
871  */
872  difference_type_1 coordinateToScanOrderIndex(const difference_type &d) const
873  {
874  return detail::CoordinateToScanOrder<actual_dimension>::exec(m_shape, d);
875  }
876 
877  /** 1D array access. Use only if N == 1.
878  */
879  reference operator() (difference_type_1 x)
880  {
881  VIGRA_ASSERT_INSIDE(difference_type(x));
882  return m_ptr [detail::CoordinatesToOffest<C>::exec(m_stride, x)];
883  }
884 
885  /** 2D array access. Use only if N == 2.
886  */
887  reference operator() (difference_type_1 x, difference_type_1 y)
888  {
889  VIGRA_ASSERT_INSIDE(difference_type(x, y));
890  return m_ptr [detail::CoordinatesToOffest<C>::exec(m_stride, x, y)];
891  }
892 
893  /** 3D array access. Use only if N == 3.
894  */
895  reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z)
896  {
897  VIGRA_ASSERT_INSIDE(difference_type(x, y, z));
898  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
899  }
900 
901  /** 4D array access. Use only if N == 4.
902  */
903  reference operator() (difference_type_1 x, difference_type_1 y,
904  difference_type_1 z, difference_type_1 u)
905  {
906  VIGRA_ASSERT_INSIDE(difference_type(x, y, z, u));
907  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
908  }
909 
910  /** 5D array access. Use only if N == 5.
911  */
912  reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
913  difference_type_1 u, difference_type_1 v)
914  {
915  VIGRA_ASSERT_INSIDE(difference_type(x, y,z, u,v));
916  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
917  }
918 
919  /** 1D const array access. Use only if N == 1.
920  */
921  const_reference operator() (difference_type_1 x) const
922  {
923  VIGRA_ASSERT_INSIDE(difference_type(x));
924  return m_ptr [detail::CoordinatesToOffest<C>::exec(m_stride, x)];
925  }
926 
927  /** 2D const array access. Use only if N == 2.
928  */
929  const_reference operator() (difference_type_1 x, difference_type_1 y) const
930  {
931  VIGRA_ASSERT_INSIDE(difference_type(x, y));
932  return m_ptr [detail::CoordinatesToOffest<C>::exec(m_stride, x, y)];
933  }
934 
935  /** 3D const array access. Use only if N == 3.
936  */
937  const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z) const
938  {
939  VIGRA_ASSERT_INSIDE(difference_type(x,y,z));
940  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
941  }
942 
943  /** 4D const array access. Use only if N == 4.
944  */
945  const_reference operator() (difference_type_1 x, difference_type_1 y,
946  difference_type_1 z, difference_type_1 u) const
947  {
948  VIGRA_ASSERT_INSIDE(difference_type(x,y,z,u));
949  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
950  }
951 
952  /** 5D const array access. Use only if N == 5.
953  */
954  const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
955  difference_type_1 u, difference_type_1 v) const
956  {
957  VIGRA_ASSERT_INSIDE(difference_type(x,y,z,u,v));
958  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
959  }
960 
961  /** Init with a constant.
962  */
963  template <class U>
964  MultiArrayView & init(const U & init)
965  {
966  detail::copyScalarMultiArrayData(traverser_begin(), shape(), init, MetaInt<actual_dimension-1>());
967  return *this;
968  }
969 
970 
971  /** Copy the data of the right-hand array (array shapes must match).
972  */
973  void copy(const MultiArrayView & rhs)
974  {
975  if(this == &rhs)
976  return;
977  this->copyImpl(rhs);
978  }
979 
980  /** Copy the data of the right-hand array (array shapes must match).
981  */
982  template <class U, class CN>
984  {
985  this->copyImpl(rhs);
986  }
987 
988  /** swap the data between two MultiArrayView objects.
989 
990  The shapes of the two array must match.
991  */
993  {
994  if(this != &rhs)
995  swapDataImpl(rhs);
996  }
997 
998  /** swap the data between two MultiArrayView objects.
999 
1000  The shapes of the two array must match.
1001  */
1002  template <class T2, class C2>
1004  {
1005  swapDataImpl(rhs);
1006  }
1007 
1008  /** bind the M outmost dimensions to certain indices.
1009  this reduces the dimensionality of the image to
1010  max { 1, N-M }.
1011 
1012  <b>Usage:</b>
1013  \code
1014  // create a 3D array of size 40x30x20
1015  typedef MultiArray<3, double>::difference_type Shape;
1016  MultiArray<3, double> array3(Shape(40, 30, 20));
1017 
1018  // get a 1D array by fixing index 1 to 12, and index 2 to 10
1019  MultiArrayView <1, double> array1 = array3.bindOuter(TinyVector<MultiArrayIndex, 2>(12, 10));
1020  \endcode
1021  */
1022  template <unsigned int M>
1023  MultiArrayView <N-M, T, C> bindOuter (const TinyVector <MultiArrayIndex, M> &d) const;
1024 
1025  /** bind the M innermost dimensions to certain indices.
1026  this reduces the dimensionality of the image to
1027  max { 1, N-M }.
1028 
1029  <b>Usage:</b>
1030  \code
1031  // create a 3D array of size 40x30x20
1032  typedef MultiArray<3, double>::difference_type Shape;
1033  MultiArray<3, double> array3(Shape(40, 30, 20));
1034 
1035  // get a 1D array by fixing index 0 to 12, and index 1 to 10
1036  MultiArrayView <1, double, StridedArrayTag> array1 = array3.bindInner(TinyVector<MultiArrayIndex, 2>(12, 10));
1037  \endcode
1038  */
1039  template <unsigned int M>
1041  bindInner (const TinyVector <MultiArrayIndex, M> &d) const;
1042 
1043  /** bind dimension M to index d.
1044  this reduces the dimensionality of the image to
1045  max { 1, N-1 }.
1046 
1047  <b>Usage:</b>
1048  \code
1049  // create a 3D array of size 40x30x20
1050  typedef MultiArray<3, double>::difference_type Shape;
1051  MultiArray<3, double> array3(Shape(40, 30, 20));
1052 
1053  // get a 2D array by fixing index 1 to 12
1054  MultiArrayView <2, double> array2 = array3.bind<1>(12);
1055 
1056  // get a 2D array by fixing index 0 to 23
1057  MultiArrayView <2, double, StridedArrayTag> array2a = array3.bind<0>(23);
1058  \endcode
1059  */
1060  template <unsigned int M>
1062  bind (difference_type_1 d) const;
1063 
1064  /** bind the outmost dimension to a certain index.
1065  this reduces the dimensionality of the image to
1066  max { 1, N-1 }.
1067 
1068  <b>Usage:</b>
1069  \code
1070  // create a 3D array of size 40x30x20
1071  typedef MultiArray<3, double>::difference_type Shape;
1072  MultiArray<3, double> array3(Shape(40, 30, 20));
1073 
1074  // get a 2D array by fixing the outermost index (i.e. index 2) to 12
1075  MultiArrayView <2, double> array2 = array3.bindOuter(12);
1076  \endcode
1077  */
1078  MultiArrayView <N-1, T, C> bindOuter (difference_type_1 d) const;
1079 
1080  /** bind the innermost dimension to a certain index.
1081  this reduces the dimensionality of the image to
1082  max { 1, N-1 }.
1083 
1084  <b>Usage:</b>
1085  \code
1086  // create a 3D array of size 40x30x20
1087  typedef MultiArray<3, double>::difference_type Shape;
1088  MultiArray<3, double> array3(Shape(40, 30, 20));
1089 
1090  // get a 2D array by fixing the innermost index (i.e. index 0) to 23
1091  MultiArrayView <2, double, StridedArrayTag> array2 = array3.bindInner(23);
1092  \endcode
1093  */
1094  MultiArrayView <N-1, T, StridedArrayTag> bindInner (difference_type_1 d) const;
1095 
1096  /** bind dimension m to index d.
1097  this reduces the dimensionality of the image to
1098  max { 1, N-1 }.
1099 
1100  <b>Usage:</b>
1101  \code
1102  // create a 3D array of size 40x30x20
1103  typedef MultiArray<3, double>::difference_type Shape;
1104  MultiArray<3, double> array3(Shape(40, 30, 20));
1105 
1106  // get a 2D array by fixing index 2 to 15
1107  MultiArrayView <2, double, StridedArrayTag> array2 = array3.bindAt(2, 15);
1108  \endcode
1109  */
1111  bindAt (difference_type_1 m, difference_type_1 d) const;
1112 
1113  /** Add a singleton dimension (dimension of legth 1).
1114 
1115  Singleton dimensions don't change the size of the data, but introduce
1116  a new index that can only take the value 0. This is mainly useful for
1117  the 'reduce mode' of transformMultiArray() and combineTwoMultiArrays(),
1118  because these functions require the source and destination arrays to
1119  have the same number of dimensions.
1120 
1121  The range of \a i must be <tt>0 <= i <= N</tt>. The new dimension will become
1122  the i'th index, and the old indices from i upwards will shift one
1123  place to the right.
1124 
1125  <b>Usage:</b>
1126 
1127  Suppose we want have a 2D array and want to create a 1D array that contains
1128  the row average of the first array.
1129  \code
1130  typedef MultiArrayShape<2>::type Shape2;
1131  MultiArray<2, double> original(Shape2(40, 30));
1132 
1133  typedef MultiArrayShape<1>::type Shape1;
1134  MultiArray<1, double> rowAverages(Shape1(30));
1135 
1136  // temporarily add a singleton dimension to the destination array
1137  transformMultiArray(srcMultiArrayRange(original),
1138  destMultiArrayRange(rowAverages.insertSingletonDimension(0)),
1139  FindAverage<double>());
1140  \endcode
1141  */
1143  insertSingletonDimension (difference_type_1 i) const;
1144 
1145  /** create a rectangular subarray that spans between the
1146  points p and q, where p is in the subarray, q not.
1147 
1148  <b>Usage:</b>
1149  \code
1150  // create a 3D array of size 40x30x20
1151  typedef MultiArray<3, double>::difference_type Shape;
1152  MultiArray<3, double> array3(Shape(40, 30, 20));
1153 
1154  // get a subarray set is smaller by one element at all sides
1155  MultiArrayView <3, double> subarray = array3.subarray(Shape(1,1,1), Shape(39, 29, 19));
1156  \endcode
1157  */
1158  MultiArrayView subarray (const difference_type &p,
1159  const difference_type &q) const
1160  {
1161  const difference_type_1 offset = dot (m_stride, p);
1162  return MultiArrayView (q - p, m_stride, m_ptr + offset);
1163  }
1164 
1165  /** apply an additional striding to the image, thereby reducing
1166  the shape of the array.
1167  for example, multiplying the stride of dimension one by three
1168  turns an appropriately layed out (interleaved) rgb image into
1169  a single band image.
1170  */
1172  stridearray (const difference_type &s) const
1173  {
1174  difference_type shape = m_shape;
1175  for (unsigned int i = 0; i < actual_dimension; ++i)
1176  shape [i] /= s [i];
1177  return MultiArrayView <N, T, StridedArrayTag>(shape, m_stride * s, m_ptr);
1178  }
1179 
1180  /** Transpose an array. If N==2, this implements the usual matrix transposition.
1181  For N > 2, it reverses the order of the indices.
1182 
1183  <b>Usage:</b><br>
1184  \code
1185  typedef MultiArray<2, double>::difference_type Shape;
1186  MultiArray<2, double> array(10, 20);
1187 
1188  MultiArray<2, double, StridedArrayTag> transposed = array.transpose();
1189 
1190  for(int i=0; i<array.shape(0), ++i)
1191  for(int j=0; j<array.shape(1); ++j)
1192  assert(array(i, j) == transposed(j, i));
1193  \endcode
1194  */
1196  transpose () const
1197  {
1198  difference_type shape(m_shape.begin(), difference_type::ReverseCopy),
1199  stride(m_stride.begin(), difference_type::ReverseCopy);
1201  }
1202 
1203  /** permute the dimensions of the array.
1204  The function exchanges the meaning of the dimensions without copying the data.
1205  In case of a 2-dimensional array, this is simply array transposition. In higher dimensions,
1206  there are more possibilities.
1207 
1208  <b>Usage:</b><br>
1209  \code
1210  typedef MultiArray<2, double>::difference_type Shape;
1211  MultiArray<2, double> array(10, 20);
1212 
1213  MultiArray<2, double, StridedArrayTag> transposed = array.permuteDimensions(Shape(1,0));
1214 
1215  for(int i=0; i<array.shape(0), ++i)
1216  for(int j=0; j<array.shape(1); ++j)
1217  assert(array(i, j) == transposed(j, i));
1218  \endcode
1219  */
1221  permuteDimensions (const difference_type &s) const;
1222 
1223  /** Permute the dimensions of the array so that the strides are in ascending order.
1224  Determines the appropriate permutation and then calls permuteDimensions().
1225  */
1227  permuteStridesAscending() const;
1228 
1229  /** Permute the dimensions of the array so that the strides are in descending order.
1230  Determines the appropriate permutation and then calls permuteDimensions().
1231  */
1233  permuteStridesDescending() const;
1234 
1235  /** Compute the ordering of the strides in this array.
1236  The result is describes the current permutation of the axes relative
1237  to the standard ascending stride order.
1238  */
1239  difference_type strideOrdering() const
1240  {
1241  return strideOrdering(m_stride);
1242  }
1243 
1244  /** Compute the ordering of the given strides.
1245  The result is describes the current permutation of the axes relative
1246  to the standard ascending stride order.
1247  */
1248  static difference_type strideOrdering(difference_type strides);
1249 
1250  /** number of the elements in the array.
1251  */
1252  difference_type_1 elementCount () const
1253  {
1254  difference_type_1 ret = m_shape[0];
1255  for(int i = 1; i < actual_dimension; ++i)
1256  ret *= m_shape[i];
1257  return ret;
1258  }
1259 
1260  /** number of the elements in the array.
1261  Same as <tt>elementCount()</tt>. Mostly useful to support the std::vector interface.
1262  */
1263  difference_type_1 size () const
1264  {
1265  return elementCount();
1266  }
1267 
1268  /** return the array's shape.
1269  */
1270  const difference_type & shape () const
1271  {
1272  return m_shape;
1273  }
1274 
1275  /** return the array's size at a certain dimension.
1276  */
1277  difference_type_1 size (difference_type_1 n) const
1278  {
1279  return m_shape [n];
1280  }
1281 
1282  /** return the array's shape at a certain dimension
1283  (same as <tt>size(n)</tt>).
1284  */
1285  difference_type_1 shape (difference_type_1 n) const
1286  {
1287  return m_shape [n];
1288  }
1289 
1290  /** return the array's stride for every dimension.
1291  */
1292  const difference_type & stride () const
1293  {
1294  return m_stride;
1295  }
1296 
1297  /** return the array's stride at a certain dimension.
1298  */
1299  difference_type_1 stride (int n) const
1300  {
1301  return m_stride [n];
1302  }
1303 
1304  /** check whether two arrays are elementwise equal.
1305  */
1306  template <class U, class C1>
1307  bool operator==(MultiArrayView<N, U, C1> const & rhs) const
1308  {
1309  if(this->shape() != rhs.shape())
1310  return false;
1311  return detail::equalityOfMultiArrays(traverser_begin(), shape(), rhs.traverser_begin(), MetaInt<actual_dimension-1>());
1312  }
1313 
1314  /** check whether two arrays are not elementwise equal.
1315  Also true when the two arrays have different shapes.
1316  */
1317  template <class U, class C1>
1318  bool operator!=(MultiArrayView<N, U, C1> const & rhs) const
1319  {
1320  return !operator==(rhs);
1321  }
1322 
1323  /** check whether the given point is in the array range.
1324  */
1325  bool isInside (difference_type const & p) const
1326  {
1327  for(int d=0; d<actual_dimension; ++d)
1328  if(p[d] < 0 || p[d] >= shape(d))
1329  return false;
1330  return true;
1331  }
1332 
1333  /** Compute the squared Euclidean norm of the array (sum of squares of the array elements).
1334  */
1335  typename NormTraits<MultiArrayView>::SquaredNormType squaredNorm() const
1336  {
1337  typedef typename NormTraits<MultiArrayView>::SquaredNormType SquaredNormType;
1338  SquaredNormType res = NumericTraits<SquaredNormType>::zero();
1339  detail::sumOverMultiArray(traverser_begin(), shape(), detail::MultiArrayL2Functor<SquaredNormType>(),
1340  res, MetaInt<actual_dimension-1>());
1341  return res;
1342  }
1343 
1344  /** Compute various norms of the array.
1345  The norm is determined by parameter \a type:
1346 
1347  <ul>
1348  <li> type == 0: maximum norm (L-infinity): maximum of absolute values of the array elements
1349  <li> type == 1: Manhattan norm (L1): sum of absolute values of the array elements
1350  <li> type == 2: Euclidean norm (L2): square root of <tt>squaredNorm()</tt> when \a useSquaredNorm is <tt>true</tt>,<br>
1351  or direct algorithm that avoids underflow/overflow otherwise.
1352  </ul>
1353 
1354  Parameter \a useSquaredNorm has no effect when \a type != 2. Defaults: compute L2 norm as square root of
1355  <tt>squaredNorm()</tt>.
1356  */
1357  typename NormTraits<MultiArrayView>::NormType norm(int type = 2, bool useSquaredNorm = true) const;
1358 
1359  /** return the pointer to the image data
1360  */
1361  pointer data () const
1362  {
1363  return m_ptr;
1364  }
1365 
1366  /**
1367  * returns true iff this view refers to valid data,
1368  * i.e. data() is not a NULL pointer. (this is false after
1369  * default construction.)
1370  */
1371  bool hasData () const
1372  {
1373  return m_ptr != 0;
1374  }
1375 
1376  /** returns the N-dimensional MultiIterator pointing
1377  to the first element in every dimension.
1378  */
1379  traverser traverser_begin ()
1380  {
1381  traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1382  return ret;
1383  }
1384 
1385  /** returns the N-dimensional MultiIterator pointing
1386  to the const first element in every dimension.
1387  */
1388  const_traverser traverser_begin () const
1389  {
1390  const_traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1391  return ret;
1392  }
1393 
1394  /** returns the N-dimensional MultiIterator pointing
1395  beyond the last element in dimension N, and to the
1396  first element in every other dimension.
1397  */
1398  traverser traverser_end ()
1399  {
1400  traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1401  ret += m_shape [actual_dimension-1];
1402  return ret;
1403  }
1404 
1405  /** returns the N-dimensional const MultiIterator pointing
1406  beyond the last element in dimension N, and to the
1407  first element in every other dimension.
1408  */
1409  const_traverser traverser_end () const
1410  {
1411  const_traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1412  ret += m_shape [actual_dimension-1];
1413  return ret;
1414  }
1415 
1416  view_type view ()
1417  {
1418  return *this;
1419  }
1420 };
1421 
1422 template <unsigned int N, class T, class C>
1423 MultiArrayView<N, T, C> &
1424 MultiArrayView <N, T, C>::operator=(MultiArrayView<N, T, C> const & rhs)
1425 {
1426  if(this == &rhs)
1427  return *this;
1428  vigra_precondition(this->shape() == rhs.shape() || m_ptr == 0,
1429  "MultiArrayView::operator=(MultiArrayView const &) size mismatch - use MultiArrayView::reset().");
1430  if(m_ptr == 0)
1431  {
1432  m_shape = rhs.m_shape;
1433  m_stride = rhs.m_stride;
1434  m_ptr = rhs.m_ptr;
1435  }
1436  else
1437  this->copyImpl(rhs);
1438  return *this;
1439 }
1440 
1441 template <unsigned int N, class T, class C>
1442 template <class CN>
1443 bool
1444 MultiArrayView <N, T, C>::arraysOverlap(const MultiArrayView <N, T, CN>& rhs) const
1445 {
1446  vigra_precondition (shape () == rhs.shape (),
1447  "MultiArrayView::arraysOverlap(): shape mismatch.");
1448  const_pointer first_element = this->m_ptr,
1449  last_element = first_element + dot(this->m_shape - difference_type(1), this->m_stride);
1450  typename MultiArrayView <N, T, CN>::const_pointer
1451  rhs_first_element = rhs.data(),
1452  rhs_last_element = rhs_first_element + dot(rhs.shape() - difference_type(1), rhs.stride());
1453  return !(last_element < rhs_first_element || rhs_last_element < first_element);
1454 }
1455 
1456 template <unsigned int N, class T, class C>
1457 template <class U, class CN>
1458 void
1459 MultiArrayView <N, T, C>::copyImpl(const MultiArrayView <N, U, CN>& rhs)
1460 {
1461  if(!arraysOverlap(rhs))
1462  {
1463  // no overlap -- can copy directly
1464  detail::copyMultiArrayData(rhs.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>());
1465  }
1466  else
1467  {
1468  // overlap: we got different views to the same data -- copy to intermediate memory in order to avoid
1469  // overwriting elements that are still needed on the rhs.
1470  MultiArray<N, T> tmp(rhs);
1471  detail::copyMultiArrayData(tmp.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>());
1472  }
1473 }
1474 
1475 #define VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT(name, op) \
1476 template <unsigned int N, class T, class C> \
1477 template<class U, class C1> \
1478 MultiArrayView<N, T, C> & \
1479 MultiArrayView <N, T, C>::operator op(MultiArrayView<N, U, C1> const & rhs) \
1480 { \
1481  vigra_precondition(this->shape() == rhs.shape(), "MultiArrayView::operator" #op "() size mismatch."); \
1482  if(!arraysOverlap(rhs)) \
1483  { \
1484  detail::name##MultiArrayData(rhs.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>()); \
1485  } \
1486  else \
1487  { \
1488  MultiArray<N, T> tmp(rhs); \
1489  detail::name##MultiArrayData(tmp.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>()); \
1490  } \
1491  return *this; \
1492 }
1493 
1494 VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT(copyAdd, +=)
1495 VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT(copySub, -=)
1496 VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT(copyMul, *=)
1497 VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT(copyDiv, /=)
1498 
1499 #undef VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT
1500 
1501 template <unsigned int N, class T, class C>
1502 template <class U, class CN>
1503 void
1504 MultiArrayView <N, T, C>::swapDataImpl(MultiArrayView <N, U, CN> rhs)
1505 {
1506  vigra_precondition (shape () == rhs.shape (),
1507  "MultiArrayView::swapData(): shape mismatch.");
1508 
1509  // check for overlap of this and rhs
1510  const_pointer first_element = this->m_ptr,
1511  last_element = first_element + dot(this->m_shape - difference_type(1), this->m_stride);
1512  typename MultiArrayView <N, U, CN>::const_pointer
1513  rhs_first_element = rhs.data(),
1514  rhs_last_element = rhs_first_element + dot(rhs.shape() - difference_type(1), rhs.stride());
1515  if(last_element < rhs_first_element || rhs_last_element < first_element)
1516  {
1517  // no overlap -- can swap directly
1518  detail::swapDataImpl(traverser_begin(), shape(), rhs.traverser_begin(), MetaInt<actual_dimension-1>());
1519  }
1520  else
1521  {
1522  // overlap: we got different views to the same data -- copy to intermediate memory in order to avoid
1523  // overwriting elements that are still needed.
1524  MultiArray<N, T> tmp(*this);
1525  copy(rhs);
1526  rhs.copy(tmp);
1527  }
1528 }
1529 
1530 template <unsigned int N, class T, class C>
1531 MultiArrayView <N, T, StridedArrayTag>
1533 {
1534  difference_type shape, stride, check((typename difference_type::value_type)0);
1535  for (unsigned int i = 0; i < actual_dimension; ++i)
1536  {
1537  shape[i] = m_shape[s[i]];
1538  stride[i] = m_stride[s[i]];
1539  ++check[s[i]];
1540  }
1541  vigra_precondition(check == difference_type(1),
1542  "MultiArrayView::permuteDimensions(): every dimension must occur exactly once.");
1544 }
1545 
1546 template <unsigned int N, class T, class C>
1549 {
1550  difference_type permutation;
1551  for(unsigned int k=0; k<N; ++k)
1552  permutation[k] = k;
1553  for(unsigned int k=0; k<N-1; ++k)
1554  {
1555  unsigned int smallest = k;
1556  for(unsigned int j=k+1; j<N; ++j)
1557  {
1558  if(stride[j] < stride[smallest])
1559  smallest = j;
1560  }
1561  if(smallest != k)
1562  {
1563  std::swap(stride[k], stride[smallest]);
1564  std::swap(permutation[k], permutation[smallest]);
1565  }
1566  }
1567  difference_type ordering;
1568  for(unsigned int k=0; k<N; ++k)
1569  ordering[permutation[k]] = k;
1570  return ordering;
1571 }
1572 
1573 template <unsigned int N, class T, class C>
1576 {
1577  difference_type ordering(strideOrdering(m_stride)), permutation;
1578  for(MultiArrayIndex k=0; k<N; ++k)
1579  permutation[ordering[k]] = k;
1580  return permuteDimensions(permutation);
1581 }
1582 
1583 template <unsigned int N, class T, class C>
1586 {
1587  difference_type ordering(strideOrdering(m_stride)), permutation;
1588  for(MultiArrayIndex k=0; k<N; ++k)
1589  permutation[ordering[N-1-k]] = k;
1590  return permuteDimensions(permutation);
1591 }
1592 
1593 template <unsigned int N, class T, class C>
1594 template <unsigned int M>
1595 MultiArrayView <N-M, T, C>
1597 {
1599  stride.init (m_stride.begin () + N-M, m_stride.end ());
1600  pointer ptr = m_ptr + dot (d, stride);
1601  static const int NNew = (N-M == 0) ? 1 : N-M;
1602  TinyVector <MultiArrayIndex, NNew> inner_shape, inner_stride;
1603  if (N-M == 0)
1604  {
1605  inner_shape [0] = 1;
1606  inner_stride [0] = 0;
1607  }
1608  else
1609  {
1610  inner_shape.init (m_shape.begin (), m_shape.end () - M);
1611  inner_stride.init (m_stride.begin (), m_stride.end () - M);
1612  }
1613  return MultiArrayView <N-M, T, C> (inner_shape, inner_stride, ptr);
1614 }
1615 
1616 template <unsigned int N, class T, class C>
1617 template <unsigned int M>
1618 MultiArrayView <N - M, T, StridedArrayTag>
1620 {
1622  stride.init (m_stride.begin (), m_stride.end () - N + M);
1623  pointer ptr = m_ptr + dot (d, stride);
1624  static const int NNew = (N-M == 0) ? 1 : N-M;
1625  TinyVector <MultiArrayIndex, NNew> outer_shape, outer_stride;
1626  if (N-M == 0)
1627  {
1628  outer_shape [0] = 1;
1629  outer_stride [0] = 0;
1630  }
1631  else
1632  {
1633  outer_shape.init (m_shape.begin () + M, m_shape.end ());
1634  outer_stride.init (m_stride.begin () + M, m_stride.end ());
1635  }
1636  return MultiArrayView <N-M, T, StridedArrayTag>
1637  (outer_shape, outer_stride, ptr);
1638 }
1639 
1640 template <unsigned int N, class T, class C>
1641 template <unsigned int M>
1644 {
1645  static const int NNew = (N-1 == 0) ? 1 : N-1;
1647  // the remaining dimensions are 0..n-1,n+1..N-1
1648  if (N-1 == 0)
1649  {
1650  shape[0] = 1;
1651  stride[0] = 0;
1652  }
1653  else
1654  {
1655  std::copy (m_shape.begin (), m_shape.begin () + M, shape.begin ());
1656  std::copy (m_shape.begin () + M+1, m_shape.end (),
1657  shape.begin () + M);
1658  std::copy (m_stride.begin (), m_stride.begin () + M, stride.begin ());
1659  std::copy (m_stride.begin () + M+1, m_stride.end (),
1660  stride.begin () + M);
1661  }
1662  return MultiArrayView <N-1, T, typename detail::MaybeStrided <M>::type>
1663  (shape, stride, m_ptr + d * m_stride[M]);
1664 }
1665 
1666 template <unsigned int N, class T, class C>
1667 MultiArrayView <N - 1, T, C>
1669 {
1670  static const int NNew = (N-1 == 0) ? 1 : N-1;
1671  TinyVector <MultiArrayIndex, NNew> inner_shape, inner_stride;
1672  if (N-1 == 0)
1673  {
1674  inner_shape [0] = 1;
1675  inner_stride [0] = 0;
1676  }
1677  else
1678  {
1679  inner_shape.init (m_shape.begin (), m_shape.end () - 1);
1680  inner_stride.init (m_stride.begin (), m_stride.end () - 1);
1681  }
1682  return MultiArrayView <N-1, T, C> (inner_shape, inner_stride,
1683  m_ptr + d * m_stride [N-1]);
1684 }
1685 
1686 template <unsigned int N, class T, class C>
1687 MultiArrayView <N - 1, T, StridedArrayTag>
1689 {
1690  static const int NNew = (N-1 == 0) ? 1 : N-1;
1691  TinyVector <MultiArrayIndex, NNew> outer_shape, outer_stride;
1692  if (N-1 == 0)
1693  {
1694  outer_shape [0] = 1;
1695  outer_stride [0] = 0;
1696  }
1697  else
1698  {
1699  outer_shape.init (m_shape.begin () + 1, m_shape.end ());
1700  outer_stride.init (m_stride.begin () + 1, m_stride.end ());
1701  }
1702  return MultiArrayView <N-1, T, StridedArrayTag>
1703  (outer_shape, outer_stride, m_ptr + d * m_stride [0]);
1704 }
1705 
1706 template <unsigned int N, class T, class C>
1707 MultiArrayView <N - 1, T, StridedArrayTag>
1709 {
1710  vigra_precondition (
1711  n < static_cast <int> (N),
1712  "MultiArrayView <N, T, C>::bindAt(): dimension out of range.");
1713  static const int NNew = (N-1 == 0) ? 1 : N-1;
1715  // the remaining dimensions are 0..n-1,n+1..N-1
1716  if (N-1 == 0)
1717  {
1718  shape [0] = 1;
1719  stride [0] = 0;
1720  }
1721  else
1722  {
1723  std::copy (m_shape.begin (), m_shape.begin () + n, shape.begin ());
1724  std::copy (m_shape.begin () + n+1, m_shape.end (),
1725  shape.begin () + n);
1726  std::copy (m_stride.begin (), m_stride.begin () + n, stride.begin ());
1727  std::copy (m_stride.begin () + n+1, m_stride.end (),
1728  stride.begin () + n);
1729  }
1730  return MultiArrayView <N-1, T, StridedArrayTag>
1731  (shape, stride, m_ptr + d * m_stride[n]);
1732 }
1733 
1734 template <unsigned int N, class T, class C>
1737 {
1738  vigra_precondition (
1739  0 <= i && i <= static_cast <difference_type_1> (N),
1740  "MultiArrayView <N, T, C>::insertSingletonDimension(): index out of range.");
1742  std::copy (m_shape.begin (), m_shape.begin () + i, shape.begin ());
1743  std::copy (m_shape.begin () + i, m_shape.end (), shape.begin () + i + 1);
1744  std::copy (m_stride.begin (), m_stride.begin () + i, stride.begin ());
1745  std::copy (m_stride.begin () + i, m_stride.end (), stride.begin () + i + 1);
1746  shape[i] = 1;
1747  stride[i] = 1;
1748 
1750 }
1751 
1752 template <unsigned int N, class T, class C>
1753 typename NormTraits<MultiArrayView <N, T, C> >::NormType
1754 MultiArrayView <N, T, C>::norm(int type, bool useSquaredNorm) const
1755 {
1756  typedef typename NormTraits<MultiArrayView>::NormType NormType;
1757 
1758  switch(type)
1759  {
1760  case 0:
1761  {
1762  NormType res = NumericTraits<NormType>::zero();
1763  detail::normMaxOfMultiArray(traverser_begin(), shape(), res, MetaInt<actual_dimension-1>());
1764  return res;
1765  }
1766  case 1:
1767  {
1768  NormType res = NumericTraits<NormType>::zero();
1769  detail::sumOverMultiArray(traverser_begin(), shape(), detail::MultiArrayL1Functor<NormType>(),
1770  res, MetaInt<actual_dimension-1>());
1771  return res;
1772  }
1773  case 2:
1774  {
1775  if(useSquaredNorm)
1776  {
1777  return sqrt((NormType)squaredNorm());
1778  }
1779  else
1780  {
1781  NormType normMax = NumericTraits<NormType>::zero();
1782  detail::normMaxOfMultiArray(traverser_begin(), shape(), normMax, MetaInt<actual_dimension-1>());
1783  if(normMax == NumericTraits<NormType>::zero())
1784  return normMax;
1785  NormType res = NumericTraits<NormType>::zero();
1786  detail::sumOverMultiArray(traverser_begin(), shape(), detail::MultiArrayScaledL2Functor<NormType>(normMax),
1787  res, MetaInt<actual_dimension-1>());
1788  return sqrt(res)*normMax;
1789  }
1790  }
1791  default:
1792  vigra_precondition(false, "MultiArrayView::norm(): Unknown norm type.");
1793  return NumericTraits<NormType>::zero(); // unreachable
1794  }
1795 }
1796 
1797 
1798 /********************************************************/
1799 /* */
1800 /* norm */
1801 /* */
1802 /********************************************************/
1803 
1804 template <unsigned int N, class T, class C>
1805 inline typename NormTraits<MultiArrayView <N, T, C> >::SquaredNormType
1807 {
1808  return a.squaredNorm();
1809 }
1810 
1811 template <unsigned int N, class T, class C>
1812 inline typename NormTraits<MultiArrayView <N, T, C> >::NormType
1813 norm(MultiArrayView <N, T, C> const & a)
1814 {
1815  return a.norm();
1816 }
1817 
1818 /********************************************************/
1819 /* */
1820 /* MultiArray */
1821 /* */
1822 /********************************************************/
1823 
1824 /** \brief Main <TT>MultiArray</TT> class containing the memory
1825  management.
1826 
1827 This class inherits the interface of MultiArrayView, and implements
1828 the memory ownership.
1829 MultiArray's are always unstrided, striding them creates a MultiArrayView.
1830 
1831 
1832 The template parameters are as follows
1833 \code
1834  N: the array dimension
1835 
1836  T: the type of the array elements
1837 
1838  A: the allocator used for internal storage management
1839  (default: std::allocator<T>)
1840 \endcode
1841 
1842 <b>\#include</b>
1843 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
1844 
1845 Namespace: vigra
1846 */
1847 template <unsigned int N, class T, class A /* default already declared above */>
1848 class MultiArray : public MultiArrayView <N, T>
1849 {
1850 
1851 public:
1853 
1854  /** the allocator type used to allocate the memory
1855  */
1856  typedef A allocator_type;
1857 
1858  /** the view type associated with this array.
1859  */
1861 
1862  /** the matrix type associated with this array.
1863  */
1865 
1866  /** the array's value type
1867  */
1869 
1870  /** pointer type
1871  */
1872  typedef typename view_type::pointer pointer;
1873 
1874  /** const pointer type
1875  */
1877 
1878  /** reference type (result of operator[])
1879  */
1881 
1882  /** const reference type (result of operator[] const)
1883  */
1885 
1886  /** size type
1887  */
1889 
1890  /** difference type (used for multi-dimensional offsets and indices)
1891  */
1893 
1894  /** difference and index type for a single dimension
1895  */
1897 
1898  /** traverser type
1899  */
1900  typedef typename vigra::detail::MultiIteratorChooser <
1901  UnstridedArrayTag>::template Traverser <N, T, T &, T *>::type
1903 
1904  /** traverser type to const data
1905  */
1906  typedef typename vigra::detail::MultiIteratorChooser <
1907  UnstridedArrayTag>::template Traverser <N, T, T const &, T const *>::type
1909 
1910  /** sequential (random access) iterator type
1911  */
1912  typedef T * iterator;
1913 
1914  /** sequential (random access) const iterator type
1915  */
1916  typedef T * const_iterator;
1917 
1918 protected:
1919 
1920  typedef typename difference_type::value_type diff_zero_t;
1921 
1922  /** the allocator used to allocate the memory
1923  */
1925 
1926  /** allocate memory for s pixels, write its address into the given
1927  pointer and initialize the pixels with init.
1928  */
1929  void allocate (pointer &ptr, difference_type_1 s, const_reference init);
1930 
1931  /** allocate memory for s pixels, write its address into the given
1932  pointer and initialize the linearized pixels to the values of init.
1933  */
1934  template <class U>
1935  void allocate (pointer &ptr, difference_type_1 s, U const * init);
1936 
1937  /** allocate memory, write its address into the given
1938  pointer and initialize it by copying the data from the given MultiArrayView.
1939  */
1940  template <class U, class C>
1941  void allocate (pointer &ptr, MultiArrayView<N, U, C> const & init);
1942 
1943  /** deallocate the memory (of length s) starting at the given address.
1944  */
1945  void deallocate (pointer &ptr, difference_type_1 s);
1946 
1947  template <class U, class C>
1948  void copyOrReshape (const MultiArrayView<N, U, C> &rhs);
1949 public:
1950  /** default constructor
1951  */
1953  : MultiArrayView <N, T> (difference_type (diff_zero_t(0)),
1954  difference_type (diff_zero_t(0)), 0)
1955  {}
1956 
1957  /** construct with given allocator
1958  */
1959  MultiArray (allocator_type const & alloc)
1960  : MultiArrayView <N, T> (difference_type (diff_zero_t(0)),
1961  difference_type (diff_zero_t(0)), 0),
1962  m_alloc(alloc)
1963  {}
1964 
1965  /** construct with given shape
1966  */
1967  explicit MultiArray (const difference_type &shape,
1968  allocator_type const & alloc = allocator_type());
1969 
1970  /** construct from shape with an initial value
1971  */
1972  MultiArray (const difference_type &shape, const_reference init,
1973  allocator_type const & alloc = allocator_type());
1974 
1975  /** construct from shape and copy values from the given array
1976  */
1977  MultiArray (const difference_type &shape, const_pointer init,
1978  allocator_type const & alloc = allocator_type());
1979 
1980  /** copy constructor
1981  */
1982  MultiArray (const MultiArray &rhs)
1983  : MultiArrayView <N, T> (rhs.m_shape, rhs.m_stride, 0),
1984  m_alloc (rhs.m_alloc)
1985  {
1986  allocate (this->m_ptr, this->elementCount (), rhs.data ());
1987  }
1988 
1989  /** construct by copying from a MultiArrayView
1990  */
1991  template <class U, class C>
1992  MultiArray (const MultiArrayView<N, U, C> &rhs,
1993  allocator_type const & alloc = allocator_type());
1994 
1995  /** assignment.<br>
1996  If the size of \a rhs is the same as the left-hand side arrays's old size, only
1997  the data are copied. Otherwise, new storage is allocated, which invalidates all
1998  objects (array views, iterators) depending on the lhs array.
1999  */
2001  {
2002  if (this != &rhs)
2003  this->copyOrReshape(rhs);
2004  return *this;
2005  }
2006 
2007  /** assignment from arbitrary MultiArrayView.<br>
2008  If the size of \a rhs is the same as the left-hand side arrays's old size, only
2009  the data are copied. Otherwise, new storage is allocated, which invalidates all
2010  objects (array views, iterators) depending on the lhs array.
2011  */
2012  template <class U, class C>
2014  {
2015  this->copyOrReshape(rhs);
2016  return *this;
2017  }
2018 
2019  /** Add-assignment from arbitrary MultiArrayView. Fails with
2020  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2021  */
2022  template <class U, class C>
2024  {
2025  view_type::operator+=(rhs);
2026  return *this;
2027  }
2028 
2029  /** Subtract-assignment from arbitrary MultiArrayView. Fails with
2030  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2031  */
2032  template <class U, class C>
2034  {
2035  view_type::operator-=(rhs);
2036  return *this;
2037  }
2038 
2039  /** Multiply-assignment from arbitrary MultiArrayView. Fails with
2040  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2041  */
2042  template <class U, class C>
2044  {
2045  view_type::operator*=(rhs);
2046  return *this;
2047  }
2048 
2049  /** Divide-assignment from arbitrary MultiArrayView. Fails with
2050  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2051  */
2052  template <class U, class C>
2054  {
2055  view_type::operator/=(rhs);
2056  return *this;
2057  }
2058 
2059  /** Add-assignment of a scalar.
2060  */
2061  MultiArray &operator+= (const T &rhs)
2062  {
2063  view_type::operator+=(rhs);
2064  return *this;
2065  }
2066 
2067  /** Subtract-assignment of a scalar.
2068  */
2069  MultiArray &operator-= (const T &rhs)
2070  {
2071  view_type::operator-=(rhs);
2072  return *this;
2073  }
2074 
2075  /** Multiply-assignment of a scalar.
2076  */
2077  MultiArray &operator*= (const T &rhs)
2078  {
2079  view_type::operator*=(rhs);
2080  return *this;
2081  }
2082 
2083  /** Divide-assignment of a scalar.
2084  */
2085  MultiArray &operator/= (const T &rhs)
2086  {
2087  view_type::operator/=(rhs);
2088  return *this;
2089  }
2090 
2091  /** destructor
2092  */
2094  {
2095  deallocate (this->m_ptr, this->elementCount ());
2096  }
2097 
2098 
2099  /** init elements with a constant
2100  */
2101  template <class U>
2102  MultiArray & init(const U & init)
2103  {
2104  view_type::init(init);
2105  return *this;
2106  }
2107 
2108  /** Allocate new memory with the given shape and initialize with zeros.<br>
2109  <em>Note:</em> this operation invalidates all dependent objects
2110  (array views and iterators)
2111  */
2112  void reshape (const difference_type &shape)
2113  {
2114  reshape (shape, NumericTraits <T>::zero ());
2115  }
2116 
2117  /** Allocate new memory with the given shape and initialize it
2118  with the given value.<br>
2119  <em>Note:</em> this operation invalidates all dependent objects
2120  (array views and iterators)
2121  */
2122  void reshape (const difference_type &shape, const_reference init);
2123 
2124  /** Swap the contents with another MultiArray. This is fast,
2125  because no data are copied, but only pointers and shapes swapped.
2126  <em>Note:</em> this operation invalidates all dependent objects
2127  (array views and iterators)
2128  */
2129  void swap (MultiArray & other);
2130 
2131  /** sequential iterator pointing to the first array element.
2132  */
2133  iterator begin ()
2134  {
2135  return this->data();
2136  }
2137 
2138  /** sequential iterator pointing beyond the last array element.
2139  */
2140  iterator end ()
2141  {
2142  return this->data() + this->elementCount();
2143  }
2144 
2145  /** sequential const iterator pointing to the first array element.
2146  */
2147  const_iterator begin () const
2148  {
2149  return this->data();
2150  }
2151 
2152  /** sequential const iterator pointing beyond the last array element.
2153  */
2154  const_iterator end () const
2155  {
2156  return this->data() + this->elementCount();
2157  }
2158 
2159  /** get the allocator.
2160  */
2161  allocator_type const & allocator () const
2162  {
2163  return m_alloc;
2164  }
2165 };
2166 
2167 template <unsigned int N, class T, class A>
2169  allocator_type const & alloc)
2170 : MultiArrayView <N, T> (shape,
2171  detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
2172  0),
2173  m_alloc(alloc)
2174 {
2175  if (N == 0)
2176  {
2177  this->m_shape [0] = 1;
2178  this->m_stride [0] = 0;
2179  }
2180  allocate (this->m_ptr, this->elementCount (), NumericTraits<T>::zero ());
2181 }
2182 
2183 template <unsigned int N, class T, class A>
2185  allocator_type const & alloc)
2186 : MultiArrayView <N, T> (shape,
2187  detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
2188  0),
2189  m_alloc(alloc)
2190 {
2191  if (N == 0)
2192  {
2193  this->m_shape [0] = 1;
2194  this->m_stride [0] = 0;
2195  }
2196  allocate (this->m_ptr, this->elementCount (), init);
2197 }
2198 
2199 template <unsigned int N, class T, class A>
2201  allocator_type const & alloc)
2202 : MultiArrayView <N, T> (shape,
2203  detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
2204  0),
2205  m_alloc(alloc)
2206 {
2207  if (N == 0)
2208  {
2209  this->m_shape [0] = 1;
2210  this->m_stride [0] = 0;
2211  }
2212  allocate (this->m_ptr, this->elementCount (), init);
2213 }
2214 
2215 template <unsigned int N, class T, class A>
2216 template <class U, class C>
2218  allocator_type const & alloc)
2219 : MultiArrayView <N, T> (rhs.shape(),
2220  detail::defaultStride <MultiArrayView<N,T>::actual_dimension>(rhs.shape()),
2221  0),
2222  m_alloc (alloc)
2223 {
2224  allocate (this->m_ptr, rhs);
2225 }
2226 
2227 template <unsigned int N, class T, class A>
2228 template <class U, class C>
2229 void
2231 {
2232  if (this->shape() == rhs.shape())
2233  this->copy(rhs);
2234  else
2235  {
2236  MultiArray t(rhs);
2237  this->swap(t);
2238  }
2239 }
2240 
2241 template <unsigned int N, class T, class A>
2243  const_reference initial)
2244 {
2245  if (N== 0)
2246  {
2247  return;
2248  }
2249  else if(new_shape == this->shape())
2250  {
2251  this->init(initial);
2252  }
2253  else
2254  {
2255  difference_type new_stride = detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (new_shape);
2257  T *new_ptr;
2258  allocate (new_ptr, new_size, initial);
2259  deallocate (this->m_ptr, this->elementCount ());
2260  this->m_ptr = new_ptr;
2261  this->m_shape = new_shape;
2262  this->m_stride = new_stride;
2263  }
2264 }
2265 
2266 
2267 template <unsigned int N, class T, class A>
2268 inline void
2270 {
2271  if (this == &other)
2272  return;
2273  std::swap(this->m_shape, other.m_shape);
2274  std::swap(this->m_stride, other.m_stride);
2275  std::swap(this->m_ptr, other.m_ptr);
2276  std::swap(this->m_alloc, other.m_alloc);
2277 }
2278 
2279 template <unsigned int N, class T, class A>
2282 {
2283  ptr = m_alloc.allocate ((typename A::size_type)s);
2285  try {
2286  for (i = 0; i < s; ++i)
2287  m_alloc.construct (ptr + i, init);
2288  }
2289  catch (...) {
2290  for (difference_type_1 j = 0; j < i; ++j)
2291  m_alloc.destroy (ptr + j);
2292  m_alloc.deallocate (ptr, (typename A::size_type)s);
2293  throw;
2294  }
2295 }
2296 
2297 template <unsigned int N, class T, class A>
2298 template <class U>
2300  U const * init)
2301 {
2302  ptr = m_alloc.allocate ((typename A::size_type)s);
2304  try {
2305  for (i = 0; i < s; ++i, ++init)
2306  m_alloc.construct (ptr + i, *init);
2307  }
2308  catch (...) {
2309  for (difference_type_1 j = 0; j < i; ++j)
2310  m_alloc.destroy (ptr + j);
2311  m_alloc.deallocate (ptr, (typename A::size_type)s);
2312  throw;
2313  }
2314 }
2315 
2316 template <unsigned int N, class T, class A>
2317 template <class U, class C>
2319 {
2320  difference_type_1 s = init.elementCount();
2321  ptr = m_alloc.allocate ((typename A::size_type)s);
2322  pointer p = ptr;
2323  try {
2324  detail::uninitializedCopyMultiArrayData(init.traverser_begin(), init.shape(),
2325  p, m_alloc, MetaInt<actual_dimension-1>());
2326  }
2327  catch (...) {
2328  for (pointer pp = ptr; pp < p; ++pp)
2329  m_alloc.destroy (pp);
2330  m_alloc.deallocate (ptr, (typename A::size_type)s);
2331  throw;
2332  }
2333 }
2334 
2335 template <unsigned int N, class T, class A>
2337 {
2338  if (ptr == 0)
2339  return;
2340  for (difference_type_1 i = 0; i < s; ++i)
2341  m_alloc.destroy (ptr + i);
2342  m_alloc.deallocate (ptr, (typename A::size_type)s);
2343  ptr = 0;
2344 }
2345 
2346 /********************************************************/
2347 /* */
2348 /* argument object factories */
2349 /* */
2350 /********************************************************/
2351 
2352 template <unsigned int N, class T, class C>
2353 inline triple<typename MultiArrayView<N,T,C>::const_traverser,
2356 srcMultiArrayRange( MultiArrayView<N,T,C> const & array )
2357 {
2358  return triple<typename MultiArrayView<N,T,C>::const_traverser,
2361  ( array.traverser_begin(),
2362  array.shape(),
2364 }
2365 
2366 template <unsigned int N, class T, class C, class Accessor>
2367 inline triple<typename MultiArrayView<N,T,C>::const_traverser,
2368  typename MultiArrayView<N,T,C>::difference_type,
2369  Accessor >
2370 srcMultiArrayRange( MultiArrayView<N,T,C> const & array, Accessor a )
2371 {
2372  return triple<typename MultiArrayView<N,T,C>::const_traverser,
2373  typename MultiArrayView<N,T,C>::difference_type,
2374  Accessor >
2375  ( array.traverser_begin(),
2376  array.shape(),
2377  a);
2378 }
2379 
2380 template <unsigned int N, class T, class C>
2381 inline pair<typename MultiArrayView<N,T,C>::const_traverser,
2382  typename AccessorTraits<T>::default_const_accessor >
2383 srcMultiArray( MultiArrayView<N,T,C> const & array )
2384 {
2385  return pair<typename MultiArrayView<N,T,C>::const_traverser,
2386  typename AccessorTraits<T>::default_const_accessor >
2387  ( array.traverser_begin(),
2388  typename AccessorTraits<T>::default_const_accessor() );
2389 }
2390 
2391 template <unsigned int N, class T, class C, class Accessor>
2392 inline pair<typename MultiArrayView<N,T,C>::const_traverser,
2393  Accessor >
2394 srcMultiArray( MultiArrayView<N,T,C> const & array, Accessor a )
2395 {
2396  return pair<typename MultiArrayView<N,T,C>::const_traverser,
2397  Accessor >
2398  ( array.traverser_begin(), a );
2399 }
2400 
2401 template <unsigned int N, class T, class C>
2402 inline triple<typename MultiArrayView<N,T,C>::traverser,
2403  typename MultiArrayView<N,T,C>::difference_type,
2404  typename AccessorTraits<T>::default_accessor >
2405 destMultiArrayRange( MultiArrayView<N,T,C> & array )
2406 {
2407  return triple<typename MultiArrayView<N,T,C>::traverser,
2408  typename MultiArrayView<N,T,C>::difference_type,
2409  typename AccessorTraits<T>::default_accessor >
2410  ( array.traverser_begin(),
2411  array.shape(),
2412  typename AccessorTraits<T>::default_accessor() );
2413 }
2414 
2415 template <unsigned int N, class T, class C, class Accessor>
2416 inline triple<typename MultiArrayView<N,T,C>::traverser,
2417  typename MultiArrayView<N,T,C>::difference_type,
2418  Accessor >
2419 destMultiArrayRange( MultiArrayView<N,T,C> & array, Accessor a )
2420 {
2421  return triple<typename MultiArrayView<N,T,C>::traverser,
2422  typename MultiArrayView<N,T,C>::difference_type,
2423  Accessor >
2424  ( array.traverser_begin(),
2425  array.shape(),
2426  a );
2427 }
2428 
2429 template <unsigned int N, class T, class C>
2430 inline pair<typename MultiArrayView<N,T,C>::traverser,
2431  typename AccessorTraits<T>::default_accessor >
2432 destMultiArray( MultiArrayView<N,T,C> & array )
2433 {
2434  return pair<typename MultiArrayView<N,T,C>::traverser,
2435  typename AccessorTraits<T>::default_accessor >
2436  ( array.traverser_begin(),
2437  typename AccessorTraits<T>::default_accessor() );
2438 }
2439 
2440 template <unsigned int N, class T, class C, class Accessor>
2441 inline pair<typename MultiArrayView<N,T,C>::traverser,
2442  Accessor >
2443 destMultiArray( MultiArrayView<N,T,C> & array, Accessor a )
2444 {
2445  return pair<typename MultiArrayView<N,T,C>::traverser,
2446  Accessor >
2447  ( array.traverser_begin(), a );
2448 }
2449 
2450 /********************************************************************/
2451 
2452 template <class PixelType, class Accessor>
2453 inline triple<ConstStridedImageIterator<PixelType>,
2454  ConstStridedImageIterator<PixelType>, Accessor>
2455 srcImageRange(const MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
2456 {
2457  ConstStridedImageIterator<PixelType>
2458  ul(img.data(), 1, img.stride(0), img.stride(1));
2459  return triple<ConstStridedImageIterator<PixelType>,
2460  ConstStridedImageIterator<PixelType>,
2461  Accessor>(
2462  ul, ul + Size2D(img.shape(0), img.shape(1)), a);
2463 }
2464 
2465 template <class PixelType, class Accessor>
2466 inline pair<ConstStridedImageIterator<PixelType>, Accessor>
2467 srcImage(const MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
2468 {
2469  ConstStridedImageIterator<PixelType>
2470  ul(img.data(), 1, img.stride(0), img.stride(1));
2471  return pair<ConstStridedImageIterator<PixelType>, Accessor>
2472  (ul, a);
2473 }
2474 
2475 template <class PixelType, class Accessor>
2476 inline triple<StridedImageIterator<PixelType>,
2477  StridedImageIterator<PixelType>, Accessor>
2478 destImageRange(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
2479 {
2480  StridedImageIterator<PixelType>
2481  ul(img.data(), 1, img.stride(0), img.stride(1));
2482  return triple<StridedImageIterator<PixelType>,
2483  StridedImageIterator<PixelType>,
2484  Accessor>(
2485  ul, ul + Size2D(img.shape(0), img.shape(1)), a);
2486 }
2487 
2488 template <class PixelType, class Accessor>
2489 inline pair<StridedImageIterator<PixelType>, Accessor>
2490 destImage(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
2491 {
2492  StridedImageIterator<PixelType>
2493  ul(img.data(), 1, img.stride(0), img.stride(1));
2494  return pair<StridedImageIterator<PixelType>, Accessor>
2495  (ul, a);
2496 }
2497 
2498 template <class PixelType, class Accessor>
2499 inline pair<StridedImageIterator<PixelType>, Accessor>
2500 maskImage(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
2501 {
2502  StridedImageIterator<PixelType>
2503  ul(img.data(), 1, img.stride(0), img.stride(1));
2504  return pair<StridedImageIterator<PixelType>, Accessor>
2505  (ul, a);
2506 }
2507 
2508 // -------------------------------------------------------------------
2509 
2510 template <class PixelType>
2511 inline triple<ConstStridedImageIterator<PixelType>,
2512  ConstStridedImageIterator<PixelType>,
2513  typename AccessorTraits<PixelType>::default_const_accessor>
2514 srcImageRange(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
2515 {
2516  ConstStridedImageIterator<PixelType>
2517  ul(img.data(), 1, img.stride(0), img.stride(1));
2518  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
2519  return triple<ConstStridedImageIterator<PixelType>,
2520  ConstStridedImageIterator<PixelType>,
2521  Accessor>
2522  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
2523 }
2524 
2525 template <class PixelType>
2526 inline triple<ConstImageIterator<PixelType>,
2527  ConstImageIterator<PixelType>,
2528  typename AccessorTraits<PixelType>::default_const_accessor>
2529 srcImageRange(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
2530 {
2531  ConstImageIterator<PixelType>
2532  ul(img.data(), img.stride(1));
2533  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
2534  return triple<ConstImageIterator<PixelType>,
2535  ConstImageIterator<PixelType>,
2536  Accessor>
2537  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
2538 }
2539 
2540 template <class PixelType>
2541 inline pair< ConstStridedImageIterator<PixelType>,
2542  typename AccessorTraits<PixelType>::default_const_accessor>
2543 srcImage(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
2544 {
2545  ConstStridedImageIterator<PixelType>
2546  ul(img.data(), 1, img.stride(0), img.stride(1));
2547  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
2548  return pair<ConstStridedImageIterator<PixelType>,
2549  Accessor>
2550  (ul, Accessor());
2551 }
2552 
2553 template <class PixelType>
2554 inline pair< ConstImageIterator<PixelType>,
2555  typename AccessorTraits<PixelType>::default_const_accessor>
2556 srcImage(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
2557 {
2558  ConstImageIterator<PixelType>
2559  ul(img.data(), img.stride(1));
2560  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
2561  return pair<ConstImageIterator<PixelType>,
2562  Accessor>
2563  (ul, Accessor());
2564 }
2565 
2566 template <class PixelType>
2567 inline triple< StridedImageIterator<PixelType>,
2568  StridedImageIterator<PixelType>,
2569  typename AccessorTraits<PixelType>::default_accessor>
2570 destImageRange(MultiArrayView<2, PixelType, StridedArrayTag> & img)
2571 {
2572  StridedImageIterator<PixelType>
2573  ul(img.data(), 1, img.stride(0), img.stride(1));
2574  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
2575  return triple<StridedImageIterator<PixelType>,
2576  StridedImageIterator<PixelType>,
2577  Accessor>
2578  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
2579 }
2580 
2581 template <class PixelType>
2582 inline triple< ImageIterator<PixelType>,
2583  ImageIterator<PixelType>,
2584  typename AccessorTraits<PixelType>::default_accessor>
2585 destImageRange(MultiArrayView<2, PixelType, UnstridedArrayTag> & img)
2586 {
2587  ImageIterator<PixelType>
2588  ul(img.data(), img.stride(1));
2589  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
2590  return triple<ImageIterator<PixelType>,
2591  ImageIterator<PixelType>,
2592  Accessor>
2593  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
2594 }
2595 
2596 template <class PixelType>
2597 inline pair< StridedImageIterator<PixelType>,
2598  typename AccessorTraits<PixelType>::default_accessor>
2599 destImage(MultiArrayView<2, PixelType, StridedArrayTag> & img)
2600 {
2601  StridedImageIterator<PixelType>
2602  ul(img.data(), 1, img.stride(0), img.stride(1));
2603  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
2604  return pair<StridedImageIterator<PixelType>, Accessor>
2605  (ul, Accessor());
2606 }
2607 
2608 template <class PixelType>
2609 inline pair< ImageIterator<PixelType>,
2610  typename AccessorTraits<PixelType>::default_accessor>
2611 destImage(MultiArrayView<2, PixelType, UnstridedArrayTag> & img)
2612 {
2613  ImageIterator<PixelType> ul(img.data(), img.stride(1));
2614  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
2615  return pair<ImageIterator<PixelType>, Accessor>(ul, Accessor());
2616 }
2617 
2618 template <class PixelType>
2619 inline pair< ConstStridedImageIterator<PixelType>,
2620  typename AccessorTraits<PixelType>::default_accessor>
2621 maskImage(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
2622 {
2623  ConstStridedImageIterator<PixelType>
2624  ul(img.data(), 1, img.stride(0), img.stride(1));
2625  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
2626  return pair<ConstStridedImageIterator<PixelType>, Accessor>
2627  (ul, Accessor());
2628 }
2629 
2630 template <class PixelType>
2631 inline pair< ConstImageIterator<PixelType>,
2632  typename AccessorTraits<PixelType>::default_accessor>
2633 maskImage(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
2634 {
2635  ConstImageIterator<PixelType>
2636  ul(img.data(), img.stride(1));
2637  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
2638  return pair<ConstImageIterator<PixelType>, Accessor>
2639  (ul, Accessor());
2640 }
2641 
2642 /********************************************************/
2643 /* */
2644 /* makeBasicImageView */
2645 /* */
2646 /********************************************************/
2647 
2648 /** \addtogroup MultiArrayToImage Wrap a \ref vigra::MultiArrayView in
2649  a \ref vigra::BasicImageView
2650 */
2651 //@{
2652 /** Create a \ref vigra::BasicImageView from an unstrided 2-dimensional
2653  \ref vigra::MultiArrayView.
2654 
2655  The \ref vigra::BasicImageView will have the same <tt>value_type </tt>
2656  as the original \ref vigra::MultiArrayView.
2657 */
2658 template <class T>
2659 BasicImageView <T>
2661 {
2662  return BasicImageView <T> (array.data (), array.shape (0),
2663  array.shape (1));
2664 }
2665 
2666 /** Create a \ref vigra::BasicImageView from a 3-dimensional
2667  \ref vigra::MultiArray.
2668 
2669  This wrapper flattens the two innermost dimensions of the array
2670  into single rows of the resulting image.
2671  The \ref vigra::BasicImageView will have the same <tt>value_type </tt>
2672  as the original \ref vigra::MultiArray.
2673 */
2674 template <class T>
2675 BasicImageView <T>
2677 {
2678  return BasicImageView <T> (array.data (),
2679  array.shape (0)*array.shape (1), array.shape (2));
2680 }
2681 
2682 /** Create a \ref vigra::BasicImageView from a 3-dimensional
2683  \ref vigra::MultiArray.
2684 
2685  This wrapper only works if <tt>T</tt> is a scalar type and the
2686  array's innermost dimension has size 3. It then re-interprets
2687  the data array as a 2-dimensional array with value_type
2688  <tt>RGBValue<T></tt>.
2689 */
2690 template <class T>
2691 BasicImageView <RGBValue<T> >
2693 {
2694  vigra_precondition (
2695  array.shape (0) == 3, "makeRGBImageView(): array.shape(0) must be 3.");
2696  return BasicImageView <RGBValue<T> > (
2697  reinterpret_cast <RGBValue <T> *> (array.data ()),
2698  array.shape (1), array.shape (2));
2699 }
2700 
2701 //@}
2702 
2703 } // namespace vigra
2704 #undef VIGRA_ASSERT_INSIDE
2705 #endif // VIGRA_MULTI_ARRAY_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.7.1 (Tue Jul 10 2012)