[ 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 <vigra/multi_array.hxx>
219 
220 Namespace: vigra::detail
221 */
222 template <class StrideTag, unsigned int N>
223 struct MaybeStrided
224 {
225  typedef StrideTag type;
226 };
227 
228 template <class StrideTag>
229 struct MaybeStrided <StrideTag, 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 <vigra/multi_array.hxx>
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 <vigra/multi_array.hxx>
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 <vigra/multi_array.hxx>
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, class Functor>
410 inline void
411 reduceOverMultiArray(SrcIterator s, Shape const & shape, T & result, Functor const & f, MetaInt<0>)
412 {
413  SrcIterator send = s + shape[0];
414  for(; s < send; ++s)
415  {
416  f(result, *s);
417  }
418 }
419 
420 template <class SrcIterator, class Shape, class T, class Functor, int N>
421 void
422 reduceOverMultiArray(SrcIterator s, Shape const & shape, T & result, Functor const & f, MetaInt<N>)
423 {
424  SrcIterator send = s + shape[N];
425  for(; s < send; ++s)
426  {
427  reduceOverMultiArray(s.begin(), shape, result, f, MetaInt<N-1>());
428  }
429 }
430 
431 struct MaxNormReduceFunctor
432 {
433  template <class T, class U>
434  void operator()(T & result, U const & u) const
435  {
436  T v = norm(u);
437  if(result < v)
438  result = v;
439  }
440 };
441 
442 struct L1NormReduceFunctor
443 {
444  template <class T, class U>
445  void operator()(T & result, U const & u) const
446  {
447  result += norm(u);
448  }
449 };
450 
451 struct SquaredL2NormReduceFunctor
452 {
453  template <class T, class U>
454  void operator()(T & result, U const & u) const
455  {
456  result += squaredNorm(u);
457  }
458 };
459 
460 template <class T>
461 struct WeightedL2NormReduceFunctor
462 {
463  T scale;
464 
465  WeightedL2NormReduceFunctor(T s)
466  : scale(s)
467  {}
468 
469  template <class U>
470  void operator()(T & result, U const & u) const
471  {
472  result += squaredNorm(u * scale);
473  }
474 };
475 
476 struct SumReduceFunctor
477 {
478  template <class T, class U>
479  void operator()(T & result, U const & u) const
480  {
481  result += u;
482  }
483 };
484 
485 struct ProdReduceFunctor
486 {
487  template <class T, class U>
488  void operator()(T & result, U const & u) const
489  {
490  result *= u;
491  }
492 };
493 
494 struct MinmaxReduceFunctor
495 {
496  template <class T, class U>
497  void operator()(T & result, U const & u) const
498  {
499  if(u < result.first)
500  result.first = u;
501  if(result.second < u)
502  result.second = u;
503  }
504 };
505 
506 struct MeanVarianceReduceFunctor
507 {
508  template <class T, class U>
509  void operator()(T & result, U const & u) const
510  {
511  ++result.first;
512  typename T::second_type t1 = u - result.second;
513  typename T::second_type t2 = t1 / result.first;
514  result.second += t2;
515  result.third += (result.first-1.0)*t1*t2;
516  }
517 };
518 
519 struct AllTrueReduceFunctor
520 {
521  template <class T, class U>
522  void operator()(T & result, U const & u) const
523  {
524  result = result && (u != NumericTraits<U>::zero());
525  }
526 };
527 
528 struct AnyTrueReduceFunctor
529 {
530  template <class T, class U>
531  void operator()(T & result, U const & u) const
532  {
533  result = result || (u != NumericTraits<U>::zero());
534  }
535 };
536 
537 template <class SrcIterator, class Shape, class DestIterator>
538 inline bool
539 equalityOfMultiArrays(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>)
540 {
541  SrcIterator send = s + shape[0];
542  for(; s < send; ++s, ++d)
543  {
544  if(!(*s == *d))
545  return false;
546  }
547  return true;
548 }
549 
550 template <class SrcIterator, class Shape, class DestIterator, int N>
551 bool
552 equalityOfMultiArrays(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>)
553 {
554  SrcIterator send = s + shape[N];
555  for(; s < send; ++s, ++d)
556  {
557  if(!equalityOfMultiArrays(s.begin(), shape, d.begin(), MetaInt<N-1>()))
558  return false;
559  }
560  return true;
561 }
562 
563 
564 template <class SrcIterator, class Shape, class DestIterator>
565 inline void
566 swapDataImpl(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>)
567 {
568  SrcIterator send = s + shape[0];
569  for(; s < send; ++s, ++d)
570  std::swap(*s, *d);
571 }
572 
573 template <class SrcIterator, class Shape, class DestIterator, int N>
574 void
575 swapDataImpl(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>)
576 {
577  SrcIterator send = s + shape[N];
578  for(; s < send; ++s, ++d)
579  swapDataImpl(s.begin(), shape, d.begin(), MetaInt<N-1>());
580 }
581 
582 } // namespace detail
583 
584 /********************************************************/
585 /* */
586 /* MultiArrayView */
587 /* */
588 /********************************************************/
589 
590 // forward declarations
591 
592 template <unsigned int N, class T, class C = UnstridedArrayTag>
593 class MultiArrayView;
594 template <unsigned int N, class T, class A = std::allocator<T> >
595 class MultiArray;
596 
597 namespace multi_math {
598 
599 template <class T>
600 struct MultiMathOperand;
601 
602 namespace detail {
603 
604 template <unsigned int N, class T, class C, class E>
605 void assign(MultiArrayView<N, T, C>, MultiMathOperand<E> const &);
606 
607 template <unsigned int N, class T, class C, class E>
608 void plusAssign(MultiArrayView<N, T, C>, MultiMathOperand<E> const &);
609 
610 template <unsigned int N, class T, class C, class E>
611 void minusAssign(MultiArrayView<N, T, C>, MultiMathOperand<E> const &);
612 
613 template <unsigned int N, class T, class C, class E>
614 void multiplyAssign(MultiArrayView<N, T, C>, MultiMathOperand<E> const &);
615 
616 template <unsigned int N, class T, class C, class E>
617 void divideAssign(MultiArrayView<N, T, C>, MultiMathOperand<E> const &);
618 
619 template <unsigned int N, class T, class A, class E>
620 void assignOrResize(MultiArray<N, T, A> &, MultiMathOperand<E> const &);
621 
622 template <unsigned int N, class T, class A, class E>
623 void plusAssignOrResize(MultiArray<N, T, A> &, MultiMathOperand<E> const &);
624 
625 template <unsigned int N, class T, class A, class E>
626 void minusAssignOrResize(MultiArray<N, T, A> &, MultiMathOperand<E> const &);
627 
628 template <unsigned int N, class T, class A, class E>
629 void multiplyAssignOrResize(MultiArray<N, T, A> &, MultiMathOperand<E> const &);
630 
631 template <unsigned int N, class T, class A, class E>
632 void divideAssignOrResize(MultiArray<N, T, A> &, MultiMathOperand<E> const &);
633 
634 } // namespace detail
635 
636 } // namespace multi_math
637 
638 template <class T> class FindSum;
639 
640 struct UnsuitableTypeForExpandElements {};
641 
642 template <class T>
643 struct ExpandElementResult
644 {
645  typedef UnsuitableTypeForExpandElements type;
646 };
647 
648 template <class T>
649 struct ExpandElementResult<std::complex<T> >
650 {
651  typedef T type;
652  enum { size = 2 };
653 };
654 
655 template <class T>
656 class FFTWComplex;
657 
658 template <class T>
659 struct ExpandElementResult<FFTWComplex<T> >
660 {
661  typedef T type;
662  enum { size = 2 };
663 };
664 
665 template <class T, int SIZE>
666 struct ExpandElementResult<TinyVector<T, SIZE> >
667 {
668  typedef T type;
669  enum { size = SIZE };
670 };
671 
672 template <class T, unsigned int R, unsigned int G, unsigned int B>
673 struct ExpandElementResult<RGBValue<T, R, G, B> >
674 {
675  typedef T type;
676  enum { size = 3 };
677 };
678 
679 #define VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(TYPE) \
680 template <> \
681 struct ExpandElementResult<TYPE> \
682 { \
683  typedef TYPE type; \
684  enum { size = 1 }; \
685 }; \
686 
687 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(bool)
688 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(char)
689 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(signed char)
690 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(signed short)
691 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(signed int)
692 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(signed long)
693 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(signed long long)
694 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(unsigned char)
695 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(unsigned short)
696 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(unsigned int)
697 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(unsigned long)
698 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(unsigned long long)
699 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(float)
700 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(double)
701 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(long double)
702 
703 #undef VIGRA_DEFINE_EXPAND_ELEMENT_RESULT
704 
705 
706 /********************************************************/
707 /* */
708 /* NormTraits */
709 /* */
710 /********************************************************/
711 
712 template <unsigned int N, class T, class C>
713 struct NormTraits<MultiArrayView<N, T, C> >
714 {
715  typedef MultiArrayView<N, T, C> Type;
716  typedef typename NormTraits<T>::SquaredNormType SquaredNormType;
717  typedef typename SquareRootTraits<SquaredNormType>::SquareRootResult NormType;
718 };
719 
720 template <unsigned int N, class T, class A>
721 struct NormTraits<MultiArray<N, T, A> >
722 : public NormTraits<MultiArrayView<N, T, UnstridedArrayTag> >
723 {
724  typedef NormTraits<MultiArrayView<N, T, UnstridedArrayTag> > BaseType;
725  typedef MultiArray<N, T, A> Type;
726  typedef typename BaseType::SquaredNormType SquaredNormType;
727  typedef typename BaseType::NormType NormType;
728 };
729 
730 /** \brief Base class for, and view to, \ref vigra::MultiArray.
731 
732 This class implements the interface of both MultiArray and
733 MultiArrayView. By default, MultiArrayViews are tagged as
734 unstrided. If necessary, strided arrays are constructed automatically
735 by calls to a variant of the bind...() function.
736 
737 In addition to the member functions described here, <tt>MultiArrayView</tt>
738 and its subclasses support arithmetic and algebraic functions via the
739 module \ref MultiMathModule.
740 
741 If you want to apply an algorithm requiring an image to a
742 <tt>MultiArrayView</tt> of appropriate (2-dimensional) shape, you can
743 create a \ref vigra::BasicImageView that acts as a wrapper with the
744 necessary interface -- see \ref MultiArrayToImage.
745 
746 The template parameter are as follows
747 \code
748  N: the array dimension
749 
750  T: the type of the array elements
751 
752  C: a tag determining whether the array's inner dimension is strided
753  or not. An array is unstrided if the array elements occupy consecutive
754  memory location, strided if there is an offset in between (e.g.
755  when a view is created that skips every other array element).
756  The compiler can generate faster code for unstrided arrays.
757  Possible values: UnstridedArrayTag (default), StridedArrayTag
758 \endcode
759 
760 <b>\#include</b>
761 <vigra/multi_array.hxx>
762 
763 Namespace: vigra
764 */
765 template <unsigned int N, class T, class StrideTag>
767 {
768 public:
769 
770  /** the array's actual dimensionality.
771  This ensures that MultiArrayView can also be used for
772  scalars (that is, when <tt>N == 0</tt>). Calculated as:<br>
773  \code
774  actual_dimension = (N==0) ? 1 : N
775  \endcode
776  */
777  enum ActualDimension { actual_dimension = (N==0) ? 1 : N };
778 
779  /** the array's value type
780  */
781  typedef T value_type;
782 
783  /** reference type (result of operator[])
784  */
785  typedef value_type &reference;
786 
787  /** const reference type (result of operator[] const)
788  */
789  typedef const value_type &const_reference;
790 
791  /** pointer type
792  */
793  typedef value_type *pointer;
794 
795  /** const pointer type
796  */
797  typedef const value_type *const_pointer;
798 
799  /** difference type (used for multi-dimensional offsets and indices)
800  */
802 
803  /** size type
804  */
805  typedef difference_type size_type;
806 
807  /** difference and index type for a single dimension
808  */
810 
811  /** scan-order iterator (StridedScanOrderIterator) type
812  */
814 
815  /** const scan-order iterator (StridedScanOrderIterator) type
816  */
818 
819  /** traverser (MultiIterator) type
820  */
821  typedef typename vigra::detail::MultiIteratorChooser <
822  StrideTag>::template Traverser <actual_dimension, T, T &, T *>::type traverser;
823 
824  /** const traverser (MultiIterator) type
825  */
826  typedef typename vigra::detail::MultiIteratorChooser <
827  StrideTag>::template Traverser <actual_dimension, T, T const &, T const *>::type const_traverser;
828 
829  /** the view type associated with this array.
830  */
832 
833  /** the matrix type associated with this array.
834  */
836 
837 protected:
838 
839  typedef typename difference_type::value_type diff_zero_t;
840 
841  /** the shape of the image pointed to is stored here.
842  */
843  difference_type m_shape;
844 
845  /** the strides (offset of a sample to the next) for every dimension
846  are stored here.
847  */
848  difference_type m_stride;
849 
850  /** pointer to the image.
851  */
852  pointer m_ptr;
853 
854  template <class U, class CN>
855  void copyImpl(const MultiArrayView <N, U, CN>& rhs);
856 
857  template <class U, class CN>
858  void swapDataImpl(MultiArrayView <N, U, CN> rhs);
859 
860  template <class CN>
861  bool arraysOverlap(const MultiArrayView <N, T, CN>& rhs) const;
862 
863  template <class U, class CN>
864  bool arraysOverlap(const MultiArrayView <N, U, CN>&) const
865  {
866  return false;
867  }
868 
869  bool checkInnerStride(UnstridedArrayTag)
870  {
871  return m_stride[0] <= 1;
872  }
873 
874  bool checkInnerStride(StridedArrayTag)
875  {
876  return true;
877  }
878 
879 public:
880 
881  /** default constructor: create an invalid view,
882  * i.e. hasData() returns false and size() is zero.
883  */
885  : m_shape (diff_zero_t(0)), m_stride (diff_zero_t(0)), m_ptr (0)
886  {}
887 
888  /** construct from shape and pointer
889  */
890  MultiArrayView (const difference_type &shape, pointer ptr)
891  : m_shape (shape),
892  m_stride (detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape)),
893  m_ptr (ptr)
894  {}
895 
896  /** Construct from shape, strides (offset of a sample to the
897  next) for every dimension, and pointer. (Note that
898  strides are not given in bytes, but in offset steps of the
899  respective pointer type.)
900  */
901  MultiArrayView (const difference_type &shape,
902  const difference_type &stride,
903  pointer ptr)
904  : m_shape (shape),
905  m_stride (stride),
906  m_ptr (ptr)
907  {
908  vigra_precondition(checkInnerStride(StrideTag()),
909  "MultiArrayView<..., UnstridedArrayTag>::MultiArrayView(): First dimension of given array is not unstrided.");
910  }
911 
912  /** Conversion to a strided view.
913  */
915  {
917  }
918 
919  /** Assignment. There are 3 cases:
920 
921  <ul>
922  <li> When this <tt>MultiArrayView</tt> does not point to valid data
923  (e.g. after default construction), it becomes a copy of \a rhs.
924  <li> When the shapes of the two arrays match, the array contents are copied.
925  <li> Otherwise, a <tt>PreconditionViolation</tt> exception is thrown.
926  </ul>
927  */
928  MultiArrayView & operator=(MultiArrayView const & rhs);
929 
930  /** Assignment of a differently typed MultiArrayView. Fails with
931  <tt>PreconditionViolation</tt> exception when the shapes do not match.
932  */
933  template<class U, class C1>
935  {
936  vigra_precondition(this->shape() == rhs.shape(),
937  "MultiArrayView::operator=() size mismatch.");
938  this->copyImpl(rhs);
939  return *this;
940  }
941 
942  /** Add-assignment of a compatible MultiArrayView. Fails with
943  <tt>PreconditionViolation</tt> exception when the shapes do not match.
944  */
945  template<class U, class C1>
947 
948  /** Subtract-assignment of a compatible MultiArrayView. Fails with
949  <tt>PreconditionViolation</tt> exception when the shapes do not match.
950  */
951  template<class U, class C1>
953 
954  /** Multiply-assignment of a compatible MultiArrayView. Fails with
955  <tt>PreconditionViolation</tt> exception when the shapes do not match.
956  */
957  template<class U, class C1>
959 
960  /** Divide-assignment of a compatible MultiArrayView. Fails with
961  <tt>PreconditionViolation</tt> exception when the shapes do not match.
962  */
963  template<class U, class C1>
965 
966  /** Add-assignment of a scalar.
967  */
968  MultiArrayView & operator+=(T const & rhs)
969  {
970  detail::copyAddScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
971  return *this;
972  }
973 
974  /** Subtract-assignment of a scalar.
975  */
976  MultiArrayView & operator-=(T const & rhs)
977  {
978  detail::copySubScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
979  return *this;
980  }
981 
982  /** Multiply-assignment of a scalar.
983  */
984  MultiArrayView & operator*=(T const & rhs)
985  {
986  detail::copyMulScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
987  return *this;
988  }
989 
990  /** Divide-assignment of a scalar.
991  */
992  MultiArrayView & operator/=(T const & rhs)
993  {
994  detail::copyDivScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
995  return *this;
996  }
997 
998  /** Assignment of an array expression. Fails with
999  <tt>PreconditionViolation</tt> exception when the shapes do not match.
1000  */
1001  template<class Expression>
1002  MultiArrayView & operator=(multi_math::MultiMathOperand<Expression> const & rhs)
1003  {
1004  multi_math::detail::assign(*this, rhs);
1005  return *this;
1006  }
1007 
1008  /** Add-assignment of an array expression. Fails with
1009  <tt>PreconditionViolation</tt> exception when the shapes do not match.
1010  */
1011  template<class Expression>
1012  MultiArrayView & operator+=(multi_math::MultiMathOperand<Expression> const & rhs)
1013  {
1014  multi_math::detail::plusAssign(*this, rhs);
1015  return *this;
1016  }
1017 
1018  /** Subtract-assignment of an array expression. Fails with
1019  <tt>PreconditionViolation</tt> exception when the shapes do not match.
1020  */
1021  template<class Expression>
1022  MultiArrayView & operator-=(multi_math::MultiMathOperand<Expression> const & rhs)
1023  {
1024  multi_math::detail::minusAssign(*this, rhs);
1025  return *this;
1026  }
1027 
1028  /** Multiply-assignment of an array expression. Fails with
1029  <tt>PreconditionViolation</tt> exception when the shapes do not match.
1030  */
1031  template<class Expression>
1032  MultiArrayView & operator*=(multi_math::MultiMathOperand<Expression> const & rhs)
1033  {
1034  multi_math::detail::multiplyAssign(*this, rhs);
1035  return *this;
1036  }
1037 
1038  /** Divide-assignment of an array expression. Fails with
1039  <tt>PreconditionViolation</tt> exception when the shapes do not match.
1040  */
1041  template<class Expression>
1042  MultiArrayView & operator/=(multi_math::MultiMathOperand<Expression> const & rhs)
1043  {
1044  multi_math::detail::divideAssign(*this, rhs);
1045  return *this;
1046  }
1047 
1048  /** array access.
1049  */
1050  reference operator[] (const difference_type &d)
1051  {
1052  VIGRA_ASSERT_INSIDE(d);
1053  return m_ptr [dot (d, m_stride)];
1054  }
1055 
1056  /** array access.
1057  */
1058  const_reference operator[] (const difference_type &d) const
1059  {
1060  VIGRA_ASSERT_INSIDE(d);
1061  return m_ptr [dot (d, m_stride)];
1062  }
1063 
1064  /** equivalent to bindInner(), when M < N.
1065  */
1066  template <unsigned int M>
1068  {
1069  return bindInner(d);
1070  }
1071 
1072  /** Array access in scan-order sense.
1073  Mostly useful to support standard indexing for 1-dimensional multi-arrays,
1074  but works for any N. Use scanOrderIndexToCoordinate() and
1075  coordinateToScanOrderIndex() for conversion between indices and coordinates.
1076 
1077  <b>Note:</b> This function should not be used in the inner loop, because the
1078  conversion of the scan order index into a memory address is expensive
1079  (it must take into account that memory may not be consecutive for subarrays
1080  and/or strided arrays). Always prefer operator() if possible.
1081  */
1082  reference operator[](difference_type_1 d)
1083  {
1084  VIGRA_ASSERT_INSIDE(scanOrderIndexToCoordinate(d));
1085  return m_ptr [detail::ScanOrderToOffset<actual_dimension>::exec(d, m_shape, m_stride)];
1086  }
1087 
1088  /** Array access in scan-order sense.
1089  Mostly useful to support standard indexing for 1-dimensional multi-arrays,
1090  but works for any N. Use scanOrderIndexToCoordinate() and
1091  coordinateToScanOrderIndex() for conversion between indices and coordinates.
1092 
1093  <b>Note:</b> This function should not be used in the inner loop, because the
1094  conversion of the scan order index into a memory address is expensive
1095  (it must take into account that memory may not be consecutive for subarrays
1096  and/or strided arrays). Always prefer operator() if possible.
1097  */
1098  const_reference operator[](difference_type_1 d) const
1099  {
1100  VIGRA_ASSERT_INSIDE(scanOrderIndexToCoordinate(d));
1101  return m_ptr [detail::ScanOrderToOffset<actual_dimension>::exec(d, m_shape, m_stride)];
1102  }
1103 
1104  /** convert scan-order index to coordinate.
1105  */
1106  difference_type scanOrderIndexToCoordinate(difference_type_1 d) const
1107  {
1108  difference_type result;
1109  detail::ScanOrderToCoordinate<actual_dimension>::exec(d, m_shape, result);
1110  return result;
1111  }
1112 
1113  /** convert coordinate to scan-order index.
1114  */
1115  difference_type_1 coordinateToScanOrderIndex(const difference_type &d) const
1116  {
1117  return detail::CoordinateToScanOrder<actual_dimension>::exec(m_shape, d);
1118  }
1119 
1120  /** 1D array access. Use only if N == 1.
1121  */
1122  reference operator() (difference_type_1 x)
1123  {
1124  VIGRA_ASSERT_INSIDE(difference_type(x));
1125  return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x)];
1126  }
1127 
1128  /** 2D array access. Use only if N == 2.
1129  */
1130  reference operator() (difference_type_1 x, difference_type_1 y)
1131  {
1132  VIGRA_ASSERT_INSIDE(difference_type(x, y));
1133  return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x, y)];
1134  }
1135 
1136  /** 3D array access. Use only if N == 3.
1137  */
1138  reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z)
1139  {
1140  VIGRA_ASSERT_INSIDE(difference_type(x, y, z));
1141  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
1142  }
1143 
1144  /** 4D array access. Use only if N == 4.
1145  */
1146  reference operator() (difference_type_1 x, difference_type_1 y,
1147  difference_type_1 z, difference_type_1 u)
1148  {
1149  VIGRA_ASSERT_INSIDE(difference_type(x, y, z, u));
1150  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
1151  }
1152 
1153  /** 5D array access. Use only if N == 5.
1154  */
1155  reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
1156  difference_type_1 u, difference_type_1 v)
1157  {
1158  VIGRA_ASSERT_INSIDE(difference_type(x, y,z, u,v));
1159  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
1160  }
1161 
1162  /** 1D const array access. Use only if N == 1.
1163  */
1164  const_reference operator() (difference_type_1 x) const
1165  {
1166  VIGRA_ASSERT_INSIDE(difference_type(x));
1167  return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x)];
1168  }
1169 
1170  /** 2D const array access. Use only if N == 2.
1171  */
1172  const_reference operator() (difference_type_1 x, difference_type_1 y) const
1173  {
1174  VIGRA_ASSERT_INSIDE(difference_type(x, y));
1175  return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x, y)];
1176  }
1177 
1178  /** 3D const array access. Use only if N == 3.
1179  */
1180  const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z) const
1181  {
1182  VIGRA_ASSERT_INSIDE(difference_type(x,y,z));
1183  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
1184  }
1185 
1186  /** 4D const array access. Use only if N == 4.
1187  */
1188  const_reference operator() (difference_type_1 x, difference_type_1 y,
1189  difference_type_1 z, difference_type_1 u) const
1190  {
1191  VIGRA_ASSERT_INSIDE(difference_type(x,y,z,u));
1192  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
1193  }
1194 
1195  /** 5D const array access. Use only if N == 5.
1196  */
1197  const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
1198  difference_type_1 u, difference_type_1 v) const
1199  {
1200  VIGRA_ASSERT_INSIDE(difference_type(x,y,z,u,v));
1201  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
1202  }
1203 
1204  /** Init with a constant.
1205  */
1206  template <class U>
1207  MultiArrayView & init(const U & init)
1208  {
1209  detail::copyScalarMultiArrayData(traverser_begin(), shape(), init, MetaInt<actual_dimension-1>());
1210  return *this;
1211  }
1212 
1213 
1214  /** Copy the data of the right-hand array (array shapes must match).
1215  */
1216  void copy(const MultiArrayView & rhs)
1217  {
1218  if(this == &rhs)
1219  return;
1220  this->copyImpl(rhs);
1221  }
1222 
1223  /** Copy the data of the right-hand array (array shapes must match).
1224  */
1225  template <class U, class CN>
1227  {
1228  this->copyImpl(rhs);
1229  }
1230 
1231  /** swap the data between two MultiArrayView objects.
1232 
1233  The shapes of the two array must match.
1234  */
1236  {
1237  if(this != &rhs)
1238  swapDataImpl(rhs);
1239  }
1240 
1241  /** swap the data between two MultiArrayView objects.
1242 
1243  The shapes of the two array must match.
1244  */
1245  template <class T2, class C2>
1247  {
1248  swapDataImpl(rhs);
1249  }
1250 
1251  /** check whether the array is unstrided (i.e. has consecutive memory) up
1252  to the given dimension.
1253 
1254  \a dimension can range from 0 ... N-1. If a certain dimension is unstrided,
1255  all lower dimensions are also unstrided.
1256  */
1257  bool isUnstrided(unsigned int dimension = N-1) const
1258  {
1259  difference_type p = shape() - difference_type(1);
1260  for(unsigned int k = dimension+1; k < N; ++k)
1261  p[k] = 0;
1262  return (&operator[](p) - m_ptr) == coordinateToScanOrderIndex(p);
1263  }
1264 
1265  /** bind the M outmost dimensions to certain indices.
1266  this reduces the dimensionality of the image to
1267  max { 1, N-M }.
1268 
1269  <b>Usage:</b>
1270  \code
1271  // create a 3D array of size 40x30x20
1272  typedef MultiArray<3, double>::difference_type Shape;
1273  MultiArray<3, double> array3(Shape(40, 30, 20));
1274 
1275  // get a 1D array by fixing index 1 to 12, and index 2 to 10
1276  MultiArrayView <1, double> array1 = array3.bindOuter(TinyVector<MultiArrayIndex, 2>(12, 10));
1277  \endcode
1278  */
1279  template <unsigned int M>
1280  MultiArrayView <N-M, T, StrideTag> bindOuter (const TinyVector <MultiArrayIndex, M> &d) const;
1281 
1282  /** bind the M innermost dimensions to certain indices.
1283  this reduces the dimensionality of the image to
1284  max { 1, N-M }.
1285 
1286  <b>Usage:</b>
1287  \code
1288  // create a 3D array of size 40x30x20
1289  typedef MultiArray<3, double>::difference_type Shape;
1290  MultiArray<3, double> array3(Shape(40, 30, 20));
1291 
1292  // get a 1D array by fixing index 0 to 12, and index 1 to 10
1293  MultiArrayView <1, double, StridedArrayTag> array1 = array3.bindInner(TinyVector<MultiArrayIndex, 2>(12, 10));
1294  \endcode
1295  */
1296  template <unsigned int M>
1298  bindInner (const TinyVector <MultiArrayIndex, M> &d) const;
1299 
1300  /** bind dimension M to index d.
1301  this reduces the dimensionality of the image to
1302  max { 1, N-1 }.
1303 
1304  <b>Usage:</b>
1305  \code
1306  // create a 3D array of size 40x30x20
1307  typedef MultiArray<3, double>::difference_type Shape;
1308  MultiArray<3, double> array3(Shape(40, 30, 20));
1309 
1310  // get a 2D array by fixing index 1 to 12
1311  MultiArrayView <2, double> array2 = array3.bind<1>(12);
1312 
1313  // get a 2D array by fixing index 0 to 23
1314  MultiArrayView <2, double, StridedArrayTag> array2a = array3.bind<0>(23);
1315  \endcode
1316  */
1317  template <unsigned int M>
1318  MultiArrayView <N-1, T, typename vigra::detail::MaybeStrided<StrideTag, M>::type >
1319  bind (difference_type_1 d) const;
1320 
1321  /** bind the outmost dimension to a certain index.
1322  this reduces the dimensionality of the image to
1323  max { 1, N-1 }.
1324 
1325  <b>Usage:</b>
1326  \code
1327  // create a 3D array of size 40x30x20
1328  typedef MultiArray<3, double>::difference_type Shape;
1329  MultiArray<3, double> array3(Shape(40, 30, 20));
1330 
1331  // get a 2D array by fixing the outermost index (i.e. index 2) to 12
1332  MultiArrayView <2, double> array2 = array3.bindOuter(12);
1333  \endcode
1334  */
1335  MultiArrayView <N-1, T, StrideTag> bindOuter (difference_type_1 d) const;
1336 
1337  /** bind the innermost dimension to a certain index.
1338  this reduces the dimensionality of the image to
1339  max { 1, N-1 }.
1340 
1341  <b>Usage:</b>
1342  \code
1343  // create a 3D array of size 40x30x20
1344  typedef MultiArray<3, double>::difference_type Shape;
1345  MultiArray<3, double> array3(Shape(40, 30, 20));
1346 
1347  // get a 2D array by fixing the innermost index (i.e. index 0) to 23
1348  MultiArrayView <2, double, StridedArrayTag> array2 = array3.bindInner(23);
1349  \endcode
1350  */
1351  MultiArrayView <N-1, T, StridedArrayTag> bindInner (difference_type_1 d) const;
1352 
1353  /** bind dimension m to index d.
1354  this reduces the dimensionality of the image to
1355  max { 1, N-1 }.
1356 
1357  <b>Usage:</b>
1358  \code
1359  // create a 3D array of size 40x30x20
1360  typedef MultiArray<3, double>::difference_type Shape;
1361  MultiArray<3, double> array3(Shape(40, 30, 20));
1362 
1363  // get a 2D array by fixing index 2 to 15
1364  MultiArrayView <2, double, StridedArrayTag> array2 = array3.bindAt(2, 15);
1365  \endcode
1366  */
1368  bindAt (difference_type_1 m, difference_type_1 d) const;
1369 
1370 
1371  /** Create a view to channel 'i' of a vector-like value type. Possible value types
1372  (of the original array) are: \ref TinyVector, \ref RGBValue, \ref FFTWComplex,
1373  and <tt>std::complex</tt>. The list can be extended to any type whose memory
1374  layout is equivalent to a fixed-size C array, by specializing
1375  <tt>ExpandElementResult</tt>.
1376 
1377  <b>Usage:</b>
1378  \code
1379  MultiArray<2, RGBValue<float> > rgb_image(Shape2(w, h));
1380 
1381  MultiArrayView<2, float, StridedArrayTag> red = rgb_image.bindElementChannel(0);
1382  MultiArrayView<2, float, StridedArrayTag> green = rgb_image.bindElementChannel(1);
1383  MultiArrayView<2, float, StridedArrayTag> blue = rgb_image.bindElementChannel(2);
1384  \endcode
1385  */
1387  bindElementChannel(difference_type_1 i) const
1388  {
1389  vigra_precondition(0 <= i && i < ExpandElementResult<T>::size,
1390  "MultiArrayView::bindElementChannel(i): 'i' out of range.");
1391  return expandElements(0).bindInner(i);
1392  }
1393 
1394  /** Create a view where a vector-like element type is expanded into a new
1395  array dimension. The new dimension is inserted at index position 'd',
1396  which must be between 0 and N inclusive.
1397 
1398  Possible value types of the original array are: \ref TinyVector, \ref RGBValue,
1399  \ref FFTWComplex, <tt>std::complex</tt>, and the built-in number types (in this
1400  case, <tt>expandElements</tt> is equivalent to <tt>insertSingletonDimension</tt>).
1401  The list of supported types can be extended to any type whose memory
1402  layout is equivalent to a fixed-size C array, by specializing
1403  <tt>ExpandElementResult</tt>.
1404 
1405  <b>Usage:</b>
1406  \code
1407  MultiArray<2, RGBValue<float> > rgb_image(Shape2(w, h));
1408 
1409  MultiArrayView<3, float, StridedArrayTag> multiband_image = rgb_image.expandElements(2);
1410  \endcode
1411  */
1413  expandElements(difference_type_1 d) const;
1414 
1415  /** Add a singleton dimension (dimension of length 1).
1416 
1417  Singleton dimensions don't change the size of the data, but introduce
1418  a new index that can only take the value 0. This is mainly useful for
1419  the 'reduce mode' of transformMultiArray() and combineTwoMultiArrays(),
1420  because these functions require the source and destination arrays to
1421  have the same number of dimensions.
1422 
1423  The range of \a i must be <tt>0 <= i <= N</tt>. The new dimension will become
1424  the i'th index, and the old indices from i upwards will shift one
1425  place to the right.
1426 
1427  <b>Usage:</b>
1428 
1429  Suppose we want have a 2D array and want to create a 1D array that contains
1430  the row average of the first array.
1431  \code
1432  typedef MultiArrayShape<2>::type Shape2;
1433  MultiArray<2, double> original(Shape2(40, 30));
1434 
1435  typedef MultiArrayShape<1>::type Shape1;
1436  MultiArray<1, double> rowAverages(Shape1(30));
1437 
1438  // temporarily add a singleton dimension to the destination array
1439  transformMultiArray(srcMultiArrayRange(original),
1440  destMultiArrayRange(rowAverages.insertSingletonDimension(0)),
1441  FindAverage<double>());
1442  \endcode
1443  */
1445  insertSingletonDimension (difference_type_1 i) const;
1446 
1447  /** create a rectangular subarray that spans between the
1448  points p and q, where p is in the subarray, q not.
1449 
1450  <b>Usage:</b>
1451  \code
1452  // create a 3D array of size 40x30x20
1453  typedef MultiArray<3, double>::difference_type Shape;
1454  MultiArray<3, double> array3(Shape(40, 30, 20));
1455 
1456  // get a subarray set is smaller by one element at all sides
1457  MultiArrayView <3, double> subarray = array3.subarray(Shape(1,1,1), Shape(39, 29, 19));
1458  \endcode
1459  */
1460  MultiArrayView subarray (const difference_type &p,
1461  const difference_type &q) const
1462  {
1463  const difference_type_1 offset = dot (m_stride, p);
1464  return MultiArrayView (q - p, m_stride, m_ptr + offset);
1465  }
1466 
1467  /** apply an additional striding to the image, thereby reducing
1468  the shape of the array.
1469  for example, multiplying the stride of dimension one by three
1470  turns an appropriately laid out (interleaved) rgb image into
1471  a single band image.
1472  */
1474  stridearray (const difference_type &s) const
1475  {
1476  difference_type shape = m_shape;
1477  for (unsigned int i = 0; i < actual_dimension; ++i)
1478  shape [i] /= s [i];
1479  return MultiArrayView <N, T, StridedArrayTag>(shape, m_stride * s, m_ptr);
1480  }
1481 
1482  /** Transpose an array. If N==2, this implements the usual matrix transposition.
1483  For N > 2, it reverses the order of the indices.
1484 
1485  <b>Usage:</b><br>
1486  \code
1487  typedef MultiArray<2, double>::difference_type Shape;
1488  MultiArray<2, double> array(10, 20);
1489 
1490  MultiArray<2, double, StridedArrayTag> transposed = array.transpose();
1491 
1492  for(int i=0; i<array.shape(0), ++i)
1493  for(int j=0; j<array.shape(1); ++j)
1494  assert(array(i, j) == transposed(j, i));
1495  \endcode
1496  */
1498  transpose () const
1499  {
1500  difference_type shape(m_shape.begin(), difference_type::ReverseCopy),
1501  stride(m_stride.begin(), difference_type::ReverseCopy);
1503  }
1504 
1505  /** permute the dimensions of the array.
1506  The function exchanges the meaning of the dimensions without copying the data.
1507  In case of a 2-dimensional array, this is simply array transposition. In higher dimensions,
1508  there are more possibilities.
1509 
1510  <b>Usage:</b><br>
1511  \code
1512  typedef MultiArray<2, double>::difference_type Shape;
1513  MultiArray<2, double> array(10, 20);
1514 
1515  MultiArray<2, double, StridedArrayTag> transposed = array.permuteDimensions(Shape(1,0));
1516 
1517  for(int i=0; i<array.shape(0), ++i)
1518  for(int j=0; j<array.shape(1); ++j)
1519  assert(array(i, j) == transposed(j, i));
1520  \endcode
1521  */
1523  permuteDimensions (const difference_type &s) const;
1524 
1525  /** Permute the dimensions of the array so that the strides are in ascending order.
1526  Determines the appropriate permutation and then calls permuteDimensions().
1527  */
1529  permuteStridesAscending() const;
1530 
1531  /** Permute the dimensions of the array so that the strides are in descending order.
1532  Determines the appropriate permutation and then calls permuteDimensions().
1533  */
1535  permuteStridesDescending() const;
1536 
1537  /** Compute the ordering of the strides in this array.
1538  The result is describes the current permutation of the axes relative
1539  to the standard ascending stride order.
1540  */
1541  difference_type strideOrdering() const
1542  {
1543  return strideOrdering(m_stride);
1544  }
1545 
1546  /** Compute the ordering of the given strides.
1547  The result is describes the current permutation of the axes relative
1548  to the standard ascending stride order.
1549  */
1550  static difference_type strideOrdering(difference_type strides);
1551 
1552  /** number of the elements in the array.
1553  */
1554  difference_type_1 elementCount () const
1555  {
1556  difference_type_1 ret = m_shape[0];
1557  for(int i = 1; i < actual_dimension; ++i)
1558  ret *= m_shape[i];
1559  return ret;
1560  }
1561 
1562  /** number of the elements in the array.
1563  Same as <tt>elementCount()</tt>. Mostly useful to support the std::vector interface.
1564  */
1565  difference_type_1 size () const
1566  {
1567  return elementCount();
1568  }
1569 
1570  /** return the array's shape.
1571  */
1572  const difference_type & shape () const
1573  {
1574  return m_shape;
1575  }
1576 
1577  /** return the array's size at a certain dimension.
1578  */
1579  difference_type_1 size (difference_type_1 n) const
1580  {
1581  return m_shape [n];
1582  }
1583 
1584  /** return the array's shape at a certain dimension
1585  (same as <tt>size(n)</tt>).
1586  */
1587  difference_type_1 shape (difference_type_1 n) const
1588  {
1589  return m_shape [n];
1590  }
1591 
1592  /** return the array's stride for every dimension.
1593  */
1594  const difference_type & stride () const
1595  {
1596  return m_stride;
1597  }
1598 
1599  /** return the array's stride at a certain dimension.
1600  */
1601  difference_type_1 stride (int n) const
1602  {
1603  return m_stride [n];
1604  }
1605 
1606  /** check whether two arrays are elementwise equal.
1607  */
1608  template <class U, class C1>
1609  bool operator==(MultiArrayView<N, U, C1> const & rhs) const
1610  {
1611  if(this->shape() != rhs.shape())
1612  return false;
1613  return detail::equalityOfMultiArrays(traverser_begin(), shape(), rhs.traverser_begin(), MetaInt<actual_dimension-1>());
1614  }
1615 
1616  /** check whether two arrays are not elementwise equal.
1617  Also true when the two arrays have different shapes.
1618  */
1619  template <class U, class C1>
1620  bool operator!=(MultiArrayView<N, U, C1> const & rhs) const
1621  {
1622  return !operator==(rhs);
1623  }
1624 
1625  /** check whether the given point is in the array range.
1626  */
1627  bool isInside (difference_type const & p) const
1628  {
1629  for(int d=0; d<actual_dimension; ++d)
1630  if(p[d] < 0 || p[d] >= shape(d))
1631  return false;
1632  return true;
1633  }
1634 
1635  /** Check if the array contains only non-zero elements (or if all elements
1636  are 'true' if the value type is 'bool').
1637  */
1638  bool all() const
1639  {
1640  bool res = true;
1641  detail::reduceOverMultiArray(traverser_begin(), shape(),
1642  res,
1643  detail::AllTrueReduceFunctor(),
1644  MetaInt<actual_dimension-1>());
1645  return res;
1646  }
1647 
1648  /** Check if the array contains a non-zero element (or an element
1649  that is 'true' if the value type is 'bool').
1650  */
1651  bool any() const
1652  {
1653  bool res = false;
1654  detail::reduceOverMultiArray(traverser_begin(), shape(),
1655  res,
1656  detail::AnyTrueReduceFunctor(),
1657  MetaInt<actual_dimension-1>());
1658  return res;
1659  }
1660 
1661  /** Find the minimum and maximum element in this array.
1662  */
1663  void minmax(T * minimum, T * maximum) const
1664  {
1665  std::pair<T, T> res(NumericTraits<T>::max(), NumericTraits<T>::min());
1666  detail::reduceOverMultiArray(traverser_begin(), shape(),
1667  res,
1668  detail::MinmaxReduceFunctor(),
1669  MetaInt<actual_dimension-1>());
1670  *minimum = res.first;
1671  *maximum = res.second;
1672  }
1673 
1674  /** Compute the mean and variance of the values in this array.
1675  */
1676  template <class U>
1677  void meanVariance(U * mean, U * variance) const
1678  {
1679  typedef typename NumericTraits<U>::RealPromote R;
1680  triple<R, R, R> res(0.0, 0.0, 0.0);
1681  detail::reduceOverMultiArray(traverser_begin(), shape(),
1682  res,
1683  detail::MeanVarianceReduceFunctor(),
1684  MetaInt<actual_dimension-1>());
1685  *mean = res.second;
1686  *variance = res.third / res.first;
1687  }
1688 
1689  /** Compute the sum of the array elements.
1690 
1691  You must provide the type of the result by an explicit template parameter:
1692  \code
1693  MultiArray<2, UInt8> A(width, height);
1694 
1695  double sum = A.sum<double>();
1696  \endcode
1697  */
1698  template <class U>
1699  U sum() const
1700  {
1701  U res = NumericTraits<U>::zero();
1702  detail::reduceOverMultiArray(traverser_begin(), shape(),
1703  res,
1704  detail::SumReduceFunctor(),
1705  MetaInt<actual_dimension-1>());
1706  return res;
1707  }
1708 
1709  /** Compute the sum of the array elements over selected axes.
1710 
1711  \arg sums must have the same shape as this array, except for the
1712  axes along which the sum is to be accumulated. These axes must be
1713  singletons. Note that you must include <tt>multi_pointoperators.hxx</tt>
1714  for this function to work.
1715 
1716  <b>Usage:</b>
1717  \code
1718  #include <vigra/multi_array.hxx>
1719  #include <vigra/multi_pointoperators.hxx>
1720 
1721  MultiArray<2, double> A(rows, cols);
1722  ... // fill A
1723 
1724  // make the first axis a singleton to sum over the first index
1725  MultiArray<2, double> rowSums(1, cols);
1726  A.sum(rowSums);
1727 
1728  // this is equivalent to
1729  transformMultiArray(srcMultiArrayRange(A),
1730  destMultiArrayRange(rowSums),
1731  FindSum<double>());
1732  \endcode
1733  */
1734  template <class U, class S>
1735  void sum(MultiArrayView<N, U, S> sums) const
1736  {
1737  transformMultiArray(srcMultiArrayRange(*this),
1738  destMultiArrayRange(sums),
1739  FindSum<U>());
1740  }
1741 
1742  /** Compute the product of the array elements.
1743 
1744  You must provide the type of the result by an explicit template parameter:
1745  \code
1746  MultiArray<2, UInt8> A(width, height);
1747 
1748  double prod = A.product<double>();
1749  \endcode
1750  */
1751  template <class U>
1752  U product() const
1753  {
1754  U res = NumericTraits<U>::one();
1755  detail::reduceOverMultiArray(traverser_begin(), shape(),
1756  res,
1757  detail::ProdReduceFunctor(),
1758  MetaInt<actual_dimension-1>());
1759  return res;
1760  }
1761 
1762  /** Compute the squared Euclidean norm of the array (sum of squares of the array elements).
1763  */
1764  typename NormTraits<MultiArrayView>::SquaredNormType
1765  squaredNorm() const
1766  {
1767  typedef typename NormTraits<MultiArrayView>::SquaredNormType SquaredNormType;
1768  SquaredNormType res = NumericTraits<SquaredNormType>::zero();
1769  detail::reduceOverMultiArray(traverser_begin(), shape(),
1770  res,
1771  detail::SquaredL2NormReduceFunctor(),
1772  MetaInt<actual_dimension-1>());
1773  return res;
1774  }
1775 
1776  /** Compute various norms of the array.
1777  The norm is determined by parameter \a type:
1778 
1779  <ul>
1780  <li> type == 0: maximum norm (L-infinity): maximum of absolute values of the array elements
1781  <li> type == 1: Manhattan norm (L1): sum of absolute values of the array elements
1782  <li> type == 2: Euclidean norm (L2): square root of <tt>squaredNorm()</tt> when \a useSquaredNorm is <tt>true</tt>,<br>
1783  or direct algorithm that avoids underflow/overflow otherwise.
1784  </ul>
1785 
1786  Parameter \a useSquaredNorm has no effect when \a type != 2. Defaults: compute L2 norm as square root of
1787  <tt>squaredNorm()</tt>.
1788  */
1789  typename NormTraits<MultiArrayView>::NormType
1790  norm(int type = 2, bool useSquaredNorm = true) const;
1791 
1792  /** return the pointer to the image data
1793  */
1794  pointer data () const
1795  {
1796  return m_ptr;
1797  }
1798 
1799  /**
1800  * returns true iff this view refers to valid data,
1801  * i.e. data() is not a NULL pointer. (this is false after
1802  * default construction.)
1803  */
1804  bool hasData () const
1805  {
1806  return m_ptr != 0;
1807  }
1808 
1809  /** returns a scan-order iterator pointing
1810  to the first array element.
1811  */
1812  iterator begin()
1813  {
1814  return iterator(m_ptr, m_shape, m_stride);
1815  }
1816 
1817  /** returns a const scan-order iterator pointing
1818  to the first array element.
1819  */
1820  const_iterator begin() const
1821  {
1822  return const_iterator(m_ptr, m_shape, m_stride);
1823  }
1824 
1825  /** returns a scan-order iterator pointing
1826  beyond the last array element.
1827  */
1828  iterator end()
1829  {
1830  return begin().getEndIterator();
1831  }
1832 
1833  /** returns a const scan-order iterator pointing
1834  beyond the last array element.
1835  */
1836  const_iterator end() const
1837  {
1838  return begin().getEndIterator();
1839  }
1840 
1841  /** returns the N-dimensional MultiIterator pointing
1842  to the first element in every dimension.
1843  */
1844  traverser traverser_begin ()
1845  {
1846  traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1847  return ret;
1848  }
1849 
1850  /** returns the N-dimensional MultiIterator pointing
1851  to the const first element in every dimension.
1852  */
1853  const_traverser traverser_begin () const
1854  {
1855  const_traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1856  return ret;
1857  }
1858 
1859  /** returns the N-dimensional MultiIterator pointing
1860  beyond the last element in dimension N, and to the
1861  first element in every other dimension.
1862  */
1863  traverser traverser_end ()
1864  {
1865  traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1866  ret += m_shape [actual_dimension-1];
1867  return ret;
1868  }
1869 
1870  /** returns the N-dimensional const MultiIterator pointing
1871  beyond the last element in dimension N, and to the
1872  first element in every other dimension.
1873  */
1874  const_traverser traverser_end () const
1875  {
1876  const_traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1877  ret += m_shape [actual_dimension-1];
1878  return ret;
1879  }
1880 
1881  view_type view ()
1882  {
1883  return *this;
1884  }
1885 };
1886 
1887 template <unsigned int N, class T, class StrideTag>
1888 MultiArrayView<N, T, StrideTag> &
1890 {
1891  if(this == &rhs)
1892  return *this;
1893  vigra_precondition(this->shape() == rhs.shape() || m_ptr == 0,
1894  "MultiArrayView::operator=(MultiArrayView const &) size mismatch.");
1895  if(m_ptr == 0)
1896  {
1897  m_shape = rhs.m_shape;
1898  m_stride = rhs.m_stride;
1899  m_ptr = rhs.m_ptr;
1900  }
1901  else
1902  this->copyImpl(rhs);
1903  return *this;
1904 }
1905 
1906 template <unsigned int N, class T, class StrideTag>
1907 template <class CN>
1908 bool
1910 {
1911  vigra_precondition (shape () == rhs.shape (),
1912  "MultiArrayView::arraysOverlap(): shape mismatch.");
1913  const_pointer first_element = this->m_ptr,
1914  last_element = first_element + dot(this->m_shape - difference_type(1), this->m_stride);
1916  rhs_first_element = rhs.data(),
1917  rhs_last_element = rhs_first_element + dot(rhs.shape() - difference_type(1), rhs.stride());
1918  return !(last_element < rhs_first_element || rhs_last_element < first_element);
1919 }
1920 
1921 template <unsigned int N, class T, class StrideTag>
1922 template <class U, class CN>
1923 void
1924 MultiArrayView <N, T, StrideTag>::copyImpl(const MultiArrayView <N, U, CN>& rhs)
1925 {
1926  if(!arraysOverlap(rhs))
1927  {
1928  // no overlap -- can copy directly
1929  detail::copyMultiArrayData(rhs.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>());
1930  }
1931  else
1932  {
1933  // overlap: we got different views to the same data -- copy to intermediate memory in order to avoid
1934  // overwriting elements that are still needed on the rhs.
1935  MultiArray<N, T> tmp(rhs);
1936  detail::copyMultiArrayData(tmp.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>());
1937  }
1938 }
1939 
1940 #define VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(name, op) \
1941 template <unsigned int N, class T, class StrideTag> \
1942 template<class U, class C1> \
1943 MultiArrayView<N, T, StrideTag> & \
1944 MultiArrayView <N, T, StrideTag>::operator op(MultiArrayView<N, U, C1> const & rhs) \
1945 { \
1946  vigra_precondition(this->shape() == rhs.shape(), "MultiArrayView::operator" #op "() size mismatch."); \
1947  if(!arraysOverlap(rhs)) \
1948  { \
1949  detail::name##MultiArrayData(rhs.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>()); \
1950  } \
1951  else \
1952  { \
1953  MultiArray<N, T> tmp(rhs); \
1954  detail::name##MultiArrayData(tmp.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>()); \
1955  } \
1956  return *this; \
1957 }
1958 
1959 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copyAdd, +=)
1960 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copySub, -=)
1961 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copyMul, *=)
1962 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copyDiv, /=)
1963 
1964 #undef VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT
1965 
1966 template <unsigned int N, class T, class StrideTag>
1967 template <class U, class CN>
1968 void
1969 MultiArrayView <N, T, StrideTag>::swapDataImpl(MultiArrayView <N, U, CN> rhs)
1970 {
1971  vigra_precondition (shape () == rhs.shape (),
1972  "MultiArrayView::swapData(): shape mismatch.");
1973 
1974  // check for overlap of this and rhs
1975  const_pointer first_element = this->m_ptr,
1976  last_element = first_element + dot(this->m_shape - difference_type(1), this->m_stride);
1977  typename MultiArrayView <N, U, CN>::const_pointer
1978  rhs_first_element = rhs.data(),
1979  rhs_last_element = rhs_first_element + dot(rhs.shape() - difference_type(1), rhs.stride());
1980  if(last_element < rhs_first_element || rhs_last_element < first_element)
1981  {
1982  // no overlap -- can swap directly
1983  detail::swapDataImpl(traverser_begin(), shape(), rhs.traverser_begin(), MetaInt<actual_dimension-1>());
1984  }
1985  else
1986  {
1987  // overlap: we got different views to the same data -- copy to intermediate memory in order to avoid
1988  // overwriting elements that are still needed.
1989  MultiArray<N, T> tmp(*this);
1990  copy(rhs);
1991  rhs.copy(tmp);
1992  }
1993 }
1994 
1995 template <unsigned int N, class T, class StrideTag>
1996 MultiArrayView <N, T, StridedArrayTag>
1998 {
1999  difference_type shape, stride, check((typename difference_type::value_type)0);
2000  for (unsigned int i = 0; i < actual_dimension; ++i)
2001  {
2002  shape[i] = m_shape[s[i]];
2003  stride[i] = m_stride[s[i]];
2004  ++check[s[i]];
2005  }
2006  vigra_precondition(check == difference_type(1),
2007  "MultiArrayView::permuteDimensions(): every dimension must occur exactly once.");
2009 }
2010 
2011 template <unsigned int N, class T, class StrideTag>
2014 {
2015  difference_type permutation;
2016  for(int k=0; k<(int)N; ++k)
2017  permutation[k] = k;
2018  for(int k=0; k<(int)N-1; ++k)
2019  {
2020  int smallest = k;
2021  for(int j=k+1; j<(int)N; ++j)
2022  {
2023  if(stride[j] < stride[smallest])
2024  smallest = j;
2025  }
2026  if(smallest != k)
2027  {
2028  std::swap(stride[k], stride[smallest]);
2029  std::swap(permutation[k], permutation[smallest]);
2030  }
2031  }
2032  difference_type ordering;
2033  for(unsigned int k=0; k<N; ++k)
2034  ordering[permutation[k]] = k;
2035  return ordering;
2036 }
2037 
2038 template <unsigned int N, class T, class StrideTag>
2041 {
2042  difference_type ordering(strideOrdering(m_stride)), permutation;
2043  for(MultiArrayIndex k=0; k<N; ++k)
2044  permutation[ordering[k]] = k;
2045  return permuteDimensions(permutation);
2046 }
2047 
2048 template <unsigned int N, class T, class StrideTag>
2051 {
2052  difference_type ordering(strideOrdering(m_stride)), permutation;
2053  for(MultiArrayIndex k=0; k<N; ++k)
2054  permutation[ordering[N-1-k]] = k;
2055  return permuteDimensions(permutation);
2056 }
2057 
2058 template <unsigned int N, class T, class StrideTag>
2059 template <unsigned int M>
2060 MultiArrayView <N-M, T, StrideTag>
2062 {
2064  stride.init (m_stride.begin () + N-M, m_stride.end ());
2065  pointer ptr = m_ptr + dot (d, stride);
2066  static const int NNew = (N-M == 0) ? 1 : N-M;
2067  TinyVector <MultiArrayIndex, NNew> inner_shape, inner_stride;
2068  if (N-M == 0)
2069  {
2070  inner_shape [0] = 1;
2071  inner_stride [0] = 0;
2072  }
2073  else
2074  {
2075  inner_shape.init (m_shape.begin (), m_shape.end () - M);
2076  inner_stride.init (m_stride.begin (), m_stride.end () - M);
2077  }
2078  return MultiArrayView <N-M, T, StrideTag> (inner_shape, inner_stride, ptr);
2079 }
2080 
2081 template <unsigned int N, class T, class StrideTag>
2082 template <unsigned int M>
2083 MultiArrayView <N - M, T, StridedArrayTag>
2085 {
2087  stride.init (m_stride.begin (), m_stride.end () - N + M);
2088  pointer ptr = m_ptr + dot (d, stride);
2089  static const int NNew = (N-M == 0) ? 1 : N-M;
2090  TinyVector <MultiArrayIndex, NNew> outer_shape, outer_stride;
2091  if (N-M == 0)
2092  {
2093  outer_shape [0] = 1;
2094  outer_stride [0] = 0;
2095  }
2096  else
2097  {
2098  outer_shape.init (m_shape.begin () + M, m_shape.end ());
2099  outer_stride.init (m_stride.begin () + M, m_stride.end ());
2100  }
2101  return MultiArrayView <N-M, T, StridedArrayTag>
2102  (outer_shape, outer_stride, ptr);
2103 }
2104 
2105 template <unsigned int N, class T, class StrideTag>
2106 template <unsigned int M>
2107 MultiArrayView <N-1, T, typename detail::MaybeStrided<StrideTag, M>::type >
2109 {
2110  static const int NNew = (N-1 == 0) ? 1 : N-1;
2112  // the remaining dimensions are 0..n-1,n+1..N-1
2113  if (N-1 == 0)
2114  {
2115  shape[0] = 1;
2116  stride[0] = 0;
2117  }
2118  else
2119  {
2120  std::copy (m_shape.begin (), m_shape.begin () + M, shape.begin ());
2121  std::copy (m_shape.begin () + M+1, m_shape.end (),
2122  shape.begin () + M);
2123  std::copy (m_stride.begin (), m_stride.begin () + M, stride.begin ());
2124  std::copy (m_stride.begin () + M+1, m_stride.end (),
2125  stride.begin () + M);
2126  }
2127  return MultiArrayView <N-1, T, typename detail::MaybeStrided<StrideTag, M>::type>
2128  (shape, stride, m_ptr + d * m_stride[M]);
2129 }
2130 
2131 template <unsigned int N, class T, class StrideTag>
2132 MultiArrayView <N - 1, T, StrideTag>
2134 {
2135  static const int NNew = (N-1 == 0) ? 1 : N-1;
2136  TinyVector <MultiArrayIndex, NNew> inner_shape, inner_stride;
2137  if (N-1 == 0)
2138  {
2139  inner_shape [0] = 1;
2140  inner_stride [0] = 0;
2141  }
2142  else
2143  {
2144  inner_shape.init (m_shape.begin (), m_shape.end () - 1);
2145  inner_stride.init (m_stride.begin (), m_stride.end () - 1);
2146  }
2147  return MultiArrayView <N-1, T, StrideTag> (inner_shape, inner_stride,
2148  m_ptr + d * m_stride [N-1]);
2149 }
2150 
2151 template <unsigned int N, class T, class StrideTag>
2152 MultiArrayView <N - 1, T, StridedArrayTag>
2154 {
2155  static const int NNew = (N-1 == 0) ? 1 : N-1;
2156  TinyVector <MultiArrayIndex, NNew> outer_shape, outer_stride;
2157  if (N-1 == 0)
2158  {
2159  outer_shape [0] = 1;
2160  outer_stride [0] = 0;
2161  }
2162  else
2163  {
2164  outer_shape.init (m_shape.begin () + 1, m_shape.end ());
2165  outer_stride.init (m_stride.begin () + 1, m_stride.end ());
2166  }
2167  return MultiArrayView <N-1, T, StridedArrayTag>
2168  (outer_shape, outer_stride, m_ptr + d * m_stride [0]);
2169 }
2170 
2171 template <unsigned int N, class T, class StrideTag>
2172 MultiArrayView <N - 1, T, StridedArrayTag>
2174 {
2175  vigra_precondition (
2176  n < static_cast <int> (N),
2177  "MultiArrayView <N, T, StrideTag>::bindAt(): dimension out of range.");
2178  static const int NNew = (N-1 == 0) ? 1 : N-1;
2180  // the remaining dimensions are 0..n-1,n+1..N-1
2181  if (N-1 == 0)
2182  {
2183  shape [0] = 1;
2184  stride [0] = 0;
2185  }
2186  else
2187  {
2188  std::copy (m_shape.begin (), m_shape.begin () + n, shape.begin ());
2189  std::copy (m_shape.begin () + n+1, m_shape.end (),
2190  shape.begin () + n);
2191  std::copy (m_stride.begin (), m_stride.begin () + n, stride.begin ());
2192  std::copy (m_stride.begin () + n+1, m_stride.end (),
2193  stride.begin () + n);
2194  }
2195  return MultiArrayView <N-1, T, StridedArrayTag>
2196  (shape, stride, m_ptr + d * m_stride[n]);
2197 }
2198 
2199 
2200 template <unsigned int N, class T, class StrideTag>
2203 {
2204  vigra_precondition(0 <= d && d <= static_cast <difference_type_1> (N),
2205  "MultiArrayView<N, ...>::expandElements(d): 0 <= 'd' <= N required.");
2206 
2207  int elementSize = ExpandElementResult<T>::size;
2208  typename MultiArrayShape<N+1>::type newShape, newStrides;
2209  for(int k=0; k<d; ++k)
2210  {
2211  newShape[k] = m_shape[k];
2212  newStrides[k] = m_stride[k]*elementSize;
2213  }
2214 
2215  newShape[d] = elementSize;
2216  newStrides[d] = 1;
2217 
2218  for(int k=d; k<N; ++k)
2219  {
2220  newShape[k+1] = m_shape[k];
2221  newStrides[k+1] = m_stride[k]*elementSize;
2222  }
2223 
2224  typedef typename ExpandElementResult<T>::type U;
2226  newShape, newStrides, reinterpret_cast<U*>(m_ptr));
2227 }
2228 
2229 template <unsigned int N, class T, class StrideTag>
2232 {
2233  vigra_precondition (
2234  0 <= i && i <= static_cast <difference_type_1> (N),
2235  "MultiArrayView <N, T, StrideTag>::insertSingletonDimension(): index out of range.");
2237  std::copy (m_shape.begin (), m_shape.begin () + i, shape.begin ());
2238  std::copy (m_shape.begin () + i, m_shape.end (), shape.begin () + i + 1);
2239  std::copy (m_stride.begin (), m_stride.begin () + i, stride.begin ());
2240  std::copy (m_stride.begin () + i, m_stride.end (), stride.begin () + i + 1);
2241  shape[i] = 1;
2242  stride[i] = 1;
2243 
2245 }
2246 
2247 template <unsigned int N, class T, class StrideTag>
2248 typename NormTraits<MultiArrayView <N, T, StrideTag> >::NormType
2249 MultiArrayView <N, T, StrideTag>::norm(int type, bool useSquaredNorm) const
2250 {
2251  typedef typename NormTraits<MultiArrayView>::NormType NormType;
2252 
2253  switch(type)
2254  {
2255  case 0:
2256  {
2257  NormType res = NumericTraits<NormType>::zero();
2258  detail::reduceOverMultiArray(traverser_begin(), shape(),
2259  res,
2260  detail::MaxNormReduceFunctor(),
2261  MetaInt<actual_dimension-1>());
2262  return res;
2263  }
2264  case 1:
2265  {
2266  NormType res = NumericTraits<NormType>::zero();
2267  detail::reduceOverMultiArray(traverser_begin(), shape(),
2268  res,
2269  detail::L1NormReduceFunctor(),
2270  MetaInt<actual_dimension-1>());
2271  return res;
2272  }
2273  case 2:
2274  {
2275  if(useSquaredNorm)
2276  {
2277  return sqrt((NormType)squaredNorm());
2278  }
2279  else
2280  {
2281  NormType normMax = NumericTraits<NormType>::zero();
2282  detail::reduceOverMultiArray(traverser_begin(), shape(),
2283  normMax,
2284  detail::MaxNormReduceFunctor(),
2285  MetaInt<actual_dimension-1>());
2286  if(normMax == NumericTraits<NormType>::zero())
2287  return normMax;
2288  NormType res = NumericTraits<NormType>::zero();
2289  detail::reduceOverMultiArray(traverser_begin(), shape(),
2290  res,
2291  detail::WeightedL2NormReduceFunctor<NormType>(1.0/normMax),
2292  MetaInt<actual_dimension-1>());
2293  return sqrt(res)*normMax;
2294  }
2295  }
2296  default:
2297  vigra_precondition(false, "MultiArrayView::norm(): Unknown norm type.");
2298  return NumericTraits<NormType>::zero(); // unreachable
2299  }
2300 }
2301 
2302 
2303 /********************************************************/
2304 /* */
2305 /* norm */
2306 /* */
2307 /********************************************************/
2308 
2309 template <unsigned int N, class T, class StrideTag>
2310 inline typename NormTraits<MultiArrayView <N, T, StrideTag> >::SquaredNormType
2312 {
2313  return a.squaredNorm();
2314 }
2315 
2316 template <unsigned int N, class T, class StrideTag>
2317 inline typename NormTraits<MultiArrayView <N, T, StrideTag> >::NormType
2318 norm(MultiArrayView <N, T, StrideTag> const & a)
2319 {
2320  return a.norm();
2321 }
2322 
2323 /********************************************************/
2324 /* */
2325 /* MultiArray */
2326 /* */
2327 /********************************************************/
2328 
2329 /** \brief Main <TT>MultiArray</TT> class containing the memory
2330  management.
2331 
2332 This class inherits the interface of MultiArrayView, and implements
2333 the memory ownership.
2334 MultiArray's are always unstrided, striding them creates a MultiArrayView.
2335 
2336 
2337 The template parameters are as follows
2338 \code
2339  N: the array dimension
2340 
2341  T: the type of the array elements
2342 
2343  A: the allocator used for internal storage management
2344  (default: std::allocator<T>)
2345 \endcode
2346 
2347 <b>\#include</b>
2348 <vigra/multi_array.hxx>
2349 
2350 Namespace: vigra
2351 */
2352 template <unsigned int N, class T, class A /* default already declared above */>
2353 class MultiArray : public MultiArrayView <N, T>
2354 {
2355 
2356 public:
2358 
2359  /** the allocator type used to allocate the memory
2360  */
2361  typedef A allocator_type;
2362 
2363  /** the view type associated with this array.
2364  */
2366 
2367  /** the matrix type associated with this array.
2368  */
2370 
2371  /** the array's value type
2372  */
2374 
2375  /** pointer type
2376  */
2377  typedef typename view_type::pointer pointer;
2378 
2379  /** const pointer type
2380  */
2382 
2383  /** reference type (result of operator[])
2384  */
2386 
2387  /** const reference type (result of operator[] const)
2388  */
2390 
2391  /** size type
2392  */
2394 
2395  /** difference type (used for multi-dimensional offsets and indices)
2396  */
2398 
2399  /** difference and index type for a single dimension
2400  */
2402 
2403  /** traverser type
2404  */
2405  typedef typename vigra::detail::MultiIteratorChooser <
2406  UnstridedArrayTag>::template Traverser <N, T, T &, T *>::type
2408 
2409  /** traverser type to const data
2410  */
2411  typedef typename vigra::detail::MultiIteratorChooser <
2412  UnstridedArrayTag>::template Traverser <N, T, T const &, T const *>::type
2414 
2415  /** sequential (random access) iterator type
2416  */
2417  typedef T * iterator;
2418 
2419  /** sequential (random access) const iterator type
2420  */
2421  typedef T * const_iterator;
2422 
2423 protected:
2424 
2425  typedef typename difference_type::value_type diff_zero_t;
2426 
2427  /** the allocator used to allocate the memory
2428  */
2430 
2431  /** allocate memory for s pixels, write its address into the given
2432  pointer and initialize the pixels with init.
2433  */
2434  void allocate (pointer &ptr, difference_type_1 s, const_reference init);
2435 
2436  /** allocate memory for s pixels, write its address into the given
2437  pointer and initialize the linearized pixels to the values of init.
2438  */
2439  template <class U>
2440  void allocate (pointer &ptr, difference_type_1 s, U const * init);
2441 
2442  /** allocate memory, write its address into the given
2443  pointer and initialize it by copying the data from the given MultiArrayView.
2444  */
2445  template <class U, class StrideTag>
2446  void allocate (pointer &ptr, MultiArrayView<N, U, StrideTag> const & init);
2447 
2448  /** deallocate the memory (of length s) starting at the given address.
2449  */
2450  void deallocate (pointer &ptr, difference_type_1 s);
2451 
2452  template <class U, class StrideTag>
2453  void copyOrReshape (const MultiArrayView<N, U, StrideTag> &rhs);
2454 public:
2455  /** default constructor
2456  */
2458  : MultiArrayView <N, T> (difference_type (diff_zero_t(0)),
2459  difference_type (diff_zero_t(0)), 0)
2460  {}
2461 
2462  /** construct with given allocator
2463  */
2464  MultiArray (allocator_type const & alloc)
2465  : MultiArrayView <N, T> (difference_type (diff_zero_t(0)),
2466  difference_type (diff_zero_t(0)), 0),
2467  m_alloc(alloc)
2468  {}
2469 
2470  /** construct with given shape
2471  */
2472  explicit MultiArray (const difference_type &shape,
2473  allocator_type const & alloc = allocator_type());
2474 
2475  /** construct from shape with an initial value
2476  */
2477  MultiArray (const difference_type &shape, const_reference init,
2478  allocator_type const & alloc = allocator_type());
2479 
2480  /** construct from shape and copy values from the given array
2481  */
2482  MultiArray (const difference_type &shape, const_pointer init,
2483  allocator_type const & alloc = allocator_type());
2484 
2485  /** copy constructor
2486  */
2487  MultiArray (const MultiArray &rhs)
2488  : MultiArrayView <N, T> (rhs.m_shape, rhs.m_stride, 0),
2489  m_alloc (rhs.m_alloc)
2490  {
2491  allocate (this->m_ptr, this->elementCount (), rhs.data ());
2492  }
2493 
2494  /** constructor from an array expression
2495  */
2496  template<class Expression>
2497  MultiArray (multi_math::MultiMathOperand<Expression> const & rhs,
2498  allocator_type const & alloc = allocator_type())
2499  : MultiArrayView <N, T> (difference_type (diff_zero_t(0)),
2500  difference_type (diff_zero_t(0)), 0),
2501  m_alloc (alloc)
2502  {
2503  multi_math::detail::assignOrResize(*this, rhs);
2504  }
2505 
2506  /** construct by copying from a MultiArrayView
2507  */
2508  template <class U, class StrideTag>
2510  allocator_type const & alloc = allocator_type());
2511 
2512  /** assignment.<br>
2513  If the size of \a rhs is the same as the left-hand side arrays's old size, only
2514  the data are copied. Otherwise, new storage is allocated, which invalidates all
2515  objects (array views, iterators) depending on the lhs array.
2516  */
2518  {
2519  if (this != &rhs)
2520  this->copyOrReshape(rhs);
2521  return *this;
2522  }
2523 
2524  /** assignment from arbitrary MultiArrayView.<br>
2525  If the size of \a rhs is the same as the left-hand side arrays's old size, only
2526  the data are copied. Otherwise, new storage is allocated, which invalidates all
2527  objects (array views, iterators) depending on the lhs array.
2528  */
2529  template <class U, class StrideTag>
2531  {
2532  this->copyOrReshape(rhs);
2533  return *this;
2534  }
2535 
2536  /** Add-assignment from arbitrary MultiArrayView. Fails with
2537  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2538  */
2539  template <class U, class StrideTag>
2541  {
2542  view_type::operator+=(rhs);
2543  return *this;
2544  }
2545 
2546  /** Subtract-assignment from arbitrary MultiArrayView. Fails with
2547  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2548  */
2549  template <class U, class StrideTag>
2551  {
2552  view_type::operator-=(rhs);
2553  return *this;
2554  }
2555 
2556  /** Multiply-assignment from arbitrary MultiArrayView. Fails with
2557  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2558  */
2559  template <class U, class StrideTag>
2561  {
2562  view_type::operator*=(rhs);
2563  return *this;
2564  }
2565 
2566  /** Divide-assignment from arbitrary MultiArrayView. Fails with
2567  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2568  */
2569  template <class U, class StrideTag>
2571  {
2572  view_type::operator/=(rhs);
2573  return *this;
2574  }
2575 
2576  /** Add-assignment of a scalar.
2577  */
2578  MultiArray &operator+= (const T &rhs)
2579  {
2580  view_type::operator+=(rhs);
2581  return *this;
2582  }
2583 
2584  /** Subtract-assignment of a scalar.
2585  */
2586  MultiArray &operator-= (const T &rhs)
2587  {
2588  view_type::operator-=(rhs);
2589  return *this;
2590  }
2591 
2592  /** Multiply-assignment of a scalar.
2593  */
2594  MultiArray &operator*= (const T &rhs)
2595  {
2596  view_type::operator*=(rhs);
2597  return *this;
2598  }
2599 
2600  /** Divide-assignment of a scalar.
2601  */
2602  MultiArray &operator/= (const T &rhs)
2603  {
2604  view_type::operator/=(rhs);
2605  return *this;
2606  }
2607  /** Assignment of an array expression. Fails with
2608  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2609  */
2610  template<class Expression>
2611  MultiArray & operator=(multi_math::MultiMathOperand<Expression> const & rhs)
2612  {
2613  multi_math::detail::assignOrResize(*this, rhs);
2614  return *this;
2615  }
2616 
2617  /** Add-assignment of an array expression. Fails with
2618  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2619  */
2620  template<class Expression>
2621  MultiArray & operator+=(multi_math::MultiMathOperand<Expression> const & rhs)
2622  {
2623  multi_math::detail::plusAssignOrResize(*this, rhs);
2624  return *this;
2625  }
2626 
2627  /** Subtract-assignment of an array expression. Fails with
2628  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2629  */
2630  template<class Expression>
2631  MultiArray & operator-=(multi_math::MultiMathOperand<Expression> const & rhs)
2632  {
2633  multi_math::detail::minusAssignOrResize(*this, rhs);
2634  return *this;
2635  }
2636 
2637  /** Multiply-assignment of an array expression. Fails with
2638  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2639  */
2640  template<class Expression>
2641  MultiArray & operator*=(multi_math::MultiMathOperand<Expression> const & rhs)
2642  {
2643  multi_math::detail::multiplyAssignOrResize(*this, rhs);
2644  return *this;
2645  }
2646 
2647  /** Divide-assignment of an array expression. Fails with
2648  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2649  */
2650  template<class Expression>
2651  MultiArray & operator/=(multi_math::MultiMathOperand<Expression> const & rhs)
2652  {
2653  multi_math::detail::divideAssignOrResize(*this, rhs);
2654  return *this;
2655  }
2656 
2657  /** destructor
2658  */
2660  {
2661  deallocate (this->m_ptr, this->elementCount ());
2662  }
2663 
2664 
2665  /** init elements with a constant
2666  */
2667  template <class U>
2668  MultiArray & init(const U & init)
2669  {
2670  view_type::init(init);
2671  return *this;
2672  }
2673 
2674  /** Allocate new memory with the given shape and initialize with zeros.<br>
2675  <em>Note:</em> this operation invalidates all dependent objects
2676  (array views and iterators)
2677  */
2678  void reshape (const difference_type &shape)
2679  {
2680  reshape (shape, T());
2681  }
2682 
2683  /** Allocate new memory with the given shape and initialize it
2684  with the given value.<br>
2685  <em>Note:</em> this operation invalidates all dependent objects
2686  (array views and iterators)
2687  */
2688  void reshape (const difference_type &shape, const_reference init);
2689 
2690  /** Swap the contents with another MultiArray. This is fast,
2691  because no data are copied, but only pointers and shapes swapped.
2692  <em>Note:</em> this operation invalidates all dependent objects
2693  (array views and iterators)
2694  */
2695  void swap (MultiArray & other);
2696 
2697  /** sequential iterator pointing to the first array element.
2698  */
2700  {
2701  return this->data();
2702  }
2703 
2704  /** sequential iterator pointing beyond the last array element.
2705  */
2707  {
2708  return this->data() + this->elementCount();
2709  }
2710 
2711  /** sequential const iterator pointing to the first array element.
2712  */
2714  {
2715  return this->data();
2716  }
2717 
2718  /** sequential const iterator pointing beyond the last array element.
2719  */
2721  {
2722  return this->data() + this->elementCount();
2723  }
2724 
2725  /** get the allocator.
2726  */
2727  allocator_type const & allocator () const
2728  {
2729  return m_alloc;
2730  }
2731 };
2732 
2733 template <unsigned int N, class T, class A>
2735  allocator_type const & alloc)
2736 : MultiArrayView <N, T> (shape,
2737  detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
2738  0),
2739  m_alloc(alloc)
2740 {
2741  if (N == 0)
2742  {
2743  this->m_shape [0] = 1;
2744  this->m_stride [0] = 0;
2745  }
2746  allocate (this->m_ptr, this->elementCount (), T());
2747 }
2748 
2749 template <unsigned int N, class T, class A>
2751  allocator_type const & alloc)
2752 : MultiArrayView <N, T> (shape,
2753  detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
2754  0),
2755  m_alloc(alloc)
2756 {
2757  if (N == 0)
2758  {
2759  this->m_shape [0] = 1;
2760  this->m_stride [0] = 0;
2761  }
2762  allocate (this->m_ptr, this->elementCount (), init);
2763 }
2764 
2765 template <unsigned int N, class T, class A>
2767  allocator_type const & alloc)
2768 : MultiArrayView <N, T> (shape,
2769  detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
2770  0),
2771  m_alloc(alloc)
2772 {
2773  if (N == 0)
2774  {
2775  this->m_shape [0] = 1;
2776  this->m_stride [0] = 0;
2777  }
2778  allocate (this->m_ptr, this->elementCount (), init);
2779 }
2780 
2781 template <unsigned int N, class T, class A>
2782 template <class U, class StrideTag>
2784  allocator_type const & alloc)
2785 : MultiArrayView <N, T> (rhs.shape(),
2786  detail::defaultStride <MultiArrayView<N,T>::actual_dimension>(rhs.shape()),
2787  0),
2788  m_alloc (alloc)
2789 {
2790  allocate (this->m_ptr, rhs);
2791 }
2792 
2793 template <unsigned int N, class T, class A>
2794 template <class U, class StrideTag>
2795 void
2797 {
2798  if (this->shape() == rhs.shape())
2799  this->copy(rhs);
2800  else
2801  {
2802  MultiArray t(rhs);
2803  this->swap(t);
2804  }
2805 }
2806 
2807 template <unsigned int N, class T, class A>
2809  const_reference initial)
2810 {
2811  if (N== 0)
2812  {
2813  return;
2814  }
2815  else if(new_shape == this->shape())
2816  {
2817  this->init(initial);
2818  }
2819  else
2820  {
2821  difference_type new_stride = detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (new_shape);
2823  T *new_ptr;
2824  allocate (new_ptr, new_size, initial);
2825  deallocate (this->m_ptr, this->elementCount ());
2826  this->m_ptr = new_ptr;
2827  this->m_shape = new_shape;
2828  this->m_stride = new_stride;
2829  }
2830 }
2831 
2832 
2833 template <unsigned int N, class T, class A>
2834 inline void
2836 {
2837  if (this == &other)
2838  return;
2839  std::swap(this->m_shape, other.m_shape);
2840  std::swap(this->m_stride, other.m_stride);
2841  std::swap(this->m_ptr, other.m_ptr);
2842  std::swap(this->m_alloc, other.m_alloc);
2843 }
2844 
2845 template <unsigned int N, class T, class A>
2848 {
2849  ptr = m_alloc.allocate ((typename A::size_type)s);
2851  try {
2852  for (i = 0; i < s; ++i)
2853  m_alloc.construct (ptr + i, init);
2854  }
2855  catch (...) {
2856  for (difference_type_1 j = 0; j < i; ++j)
2857  m_alloc.destroy (ptr + j);
2858  m_alloc.deallocate (ptr, (typename A::size_type)s);
2859  throw;
2860  }
2861 }
2862 
2863 template <unsigned int N, class T, class A>
2864 template <class U>
2866  U const * init)
2867 {
2868  ptr = m_alloc.allocate ((typename A::size_type)s);
2870  try {
2871  for (i = 0; i < s; ++i, ++init)
2872  m_alloc.construct (ptr + i, *init);
2873  }
2874  catch (...) {
2875  for (difference_type_1 j = 0; j < i; ++j)
2876  m_alloc.destroy (ptr + j);
2877  m_alloc.deallocate (ptr, (typename A::size_type)s);
2878  throw;
2879  }
2880 }
2881 
2882 template <unsigned int N, class T, class A>
2883 template <class U, class StrideTag>
2885 {
2886  difference_type_1 s = init.elementCount();
2887  ptr = m_alloc.allocate ((typename A::size_type)s);
2888  pointer p = ptr;
2889  try {
2890  detail::uninitializedCopyMultiArrayData(init.traverser_begin(), init.shape(),
2891  p, m_alloc, MetaInt<actual_dimension-1>());
2892  }
2893  catch (...) {
2894  for (pointer pp = ptr; pp < p; ++pp)
2895  m_alloc.destroy (pp);
2896  m_alloc.deallocate (ptr, (typename A::size_type)s);
2897  throw;
2898  }
2899 }
2900 
2901 template <unsigned int N, class T, class A>
2903 {
2904  if (ptr == 0)
2905  return;
2906  for (difference_type_1 i = 0; i < s; ++i)
2907  m_alloc.destroy (ptr + i);
2908  m_alloc.deallocate (ptr, (typename A::size_type)s);
2909  ptr = 0;
2910 }
2911 
2912 /********************************************************/
2913 /* */
2914 /* argument object factories */
2915 /* */
2916 /********************************************************/
2917 
2918 template <unsigned int N, class T, class StrideTag>
2919 inline triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
2922 srcMultiArrayRange( MultiArrayView<N,T,StrideTag> const & array )
2923 {
2924  return triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
2927  ( array.traverser_begin(),
2928  array.shape(),
2930 }
2931 
2932 template <unsigned int N, class T, class StrideTag, class Accessor>
2933 inline triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
2934  typename MultiArrayView<N,T,StrideTag>::difference_type,
2935  Accessor >
2936 srcMultiArrayRange( MultiArrayView<N,T,StrideTag> const & array, Accessor a )
2937 {
2938  return triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
2939  typename MultiArrayView<N,T,StrideTag>::difference_type,
2940  Accessor >
2941  ( array.traverser_begin(),
2942  array.shape(),
2943  a);
2944 }
2945 
2946 template <unsigned int N, class T, class StrideTag>
2947 inline pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
2948  typename AccessorTraits<T>::default_const_accessor >
2949 srcMultiArray( MultiArrayView<N,T,StrideTag> const & array )
2950 {
2951  return pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
2952  typename AccessorTraits<T>::default_const_accessor >
2953  ( array.traverser_begin(),
2954  typename AccessorTraits<T>::default_const_accessor() );
2955 }
2956 
2957 template <unsigned int N, class T, class StrideTag, class Accessor>
2958 inline pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
2959  Accessor >
2960 srcMultiArray( MultiArrayView<N,T,StrideTag> const & array, Accessor a )
2961 {
2962  return pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
2963  Accessor >
2964  ( array.traverser_begin(), a );
2965 }
2966 
2967 template <unsigned int N, class T, class StrideTag>
2968 inline triple<typename MultiArrayView<N,T,StrideTag>::traverser,
2969  typename MultiArrayView<N,T,StrideTag>::difference_type,
2970  typename AccessorTraits<T>::default_accessor >
2971 destMultiArrayRange( MultiArrayView<N,T,StrideTag> & array )
2972 {
2973  return triple<typename MultiArrayView<N,T,StrideTag>::traverser,
2974  typename MultiArrayView<N,T,StrideTag>::difference_type,
2975  typename AccessorTraits<T>::default_accessor >
2976  ( array.traverser_begin(),
2977  array.shape(),
2978  typename AccessorTraits<T>::default_accessor() );
2979 }
2980 
2981 template <unsigned int N, class T, class StrideTag, class Accessor>
2982 inline triple<typename MultiArrayView<N,T,StrideTag>::traverser,
2983  typename MultiArrayView<N,T,StrideTag>::difference_type,
2984  Accessor >
2985 destMultiArrayRange( MultiArrayView<N,T,StrideTag> & array, Accessor a )
2986 {
2987  return triple<typename MultiArrayView<N,T,StrideTag>::traverser,
2988  typename MultiArrayView<N,T,StrideTag>::difference_type,
2989  Accessor >
2990  ( array.traverser_begin(),
2991  array.shape(),
2992  a );
2993 }
2994 
2995 template <unsigned int N, class T, class StrideTag>
2996 inline pair<typename MultiArrayView<N,T,StrideTag>::traverser,
2997  typename AccessorTraits<T>::default_accessor >
2998 destMultiArray( MultiArrayView<N,T,StrideTag> & array )
2999 {
3000  return pair<typename MultiArrayView<N,T,StrideTag>::traverser,
3001  typename AccessorTraits<T>::default_accessor >
3002  ( array.traverser_begin(),
3003  typename AccessorTraits<T>::default_accessor() );
3004 }
3005 
3006 template <unsigned int N, class T, class StrideTag, class Accessor>
3007 inline pair<typename MultiArrayView<N,T,StrideTag>::traverser,
3008  Accessor >
3009 destMultiArray( MultiArrayView<N,T,StrideTag> & array, Accessor a )
3010 {
3011  return pair<typename MultiArrayView<N,T,StrideTag>::traverser,
3012  Accessor >
3013  ( array.traverser_begin(), a );
3014 }
3015 
3016 /********************************************************************/
3017 
3018 template <class PixelType, class Accessor>
3019 inline triple<ConstStridedImageIterator<PixelType>,
3020  ConstStridedImageIterator<PixelType>, Accessor>
3021 srcImageRange(const MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
3022 {
3023  ConstStridedImageIterator<PixelType>
3024  ul(img.data(), 1, img.stride(0), img.stride(1));
3025  return triple<ConstStridedImageIterator<PixelType>,
3026  ConstStridedImageIterator<PixelType>,
3027  Accessor>(
3028  ul, ul + Size2D(img.shape(0), img.shape(1)), a);
3029 }
3030 
3031 template <class PixelType, class Accessor>
3032 inline pair<ConstStridedImageIterator<PixelType>, Accessor>
3033 srcImage(const MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
3034 {
3035  ConstStridedImageIterator<PixelType>
3036  ul(img.data(), 1, img.stride(0), img.stride(1));
3037  return pair<ConstStridedImageIterator<PixelType>, Accessor>
3038  (ul, a);
3039 }
3040 
3041 template <class PixelType, class Accessor>
3042 inline triple<StridedImageIterator<PixelType>,
3043  StridedImageIterator<PixelType>, Accessor>
3044 destImageRange(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
3045 {
3046  StridedImageIterator<PixelType>
3047  ul(img.data(), 1, img.stride(0), img.stride(1));
3048  return triple<StridedImageIterator<PixelType>,
3049  StridedImageIterator<PixelType>,
3050  Accessor>(
3051  ul, ul + Size2D(img.shape(0), img.shape(1)), a);
3052 }
3053 
3054 template <class PixelType, class Accessor>
3055 inline pair<StridedImageIterator<PixelType>, Accessor>
3056 destImage(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
3057 {
3058  StridedImageIterator<PixelType>
3059  ul(img.data(), 1, img.stride(0), img.stride(1));
3060  return pair<StridedImageIterator<PixelType>, Accessor>
3061  (ul, a);
3062 }
3063 
3064 template <class PixelType, class Accessor>
3065 inline pair<StridedImageIterator<PixelType>, Accessor>
3066 maskImage(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
3067 {
3068  StridedImageIterator<PixelType>
3069  ul(img.data(), 1, img.stride(0), img.stride(1));
3070  return pair<StridedImageIterator<PixelType>, Accessor>
3071  (ul, a);
3072 }
3073 
3074 // -------------------------------------------------------------------
3075 
3076 template <class PixelType>
3077 inline triple<ConstStridedImageIterator<PixelType>,
3078  ConstStridedImageIterator<PixelType>,
3079  typename AccessorTraits<PixelType>::default_const_accessor>
3080 srcImageRange(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
3081 {
3082  ConstStridedImageIterator<PixelType>
3083  ul(img.data(), 1, img.stride(0), img.stride(1));
3084  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
3085  return triple<ConstStridedImageIterator<PixelType>,
3086  ConstStridedImageIterator<PixelType>,
3087  Accessor>
3088  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
3089 }
3090 
3091 template <class PixelType>
3092 inline triple<ConstImageIterator<PixelType>,
3093  ConstImageIterator<PixelType>,
3094  typename AccessorTraits<PixelType>::default_const_accessor>
3095 srcImageRange(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
3096 {
3097  ConstImageIterator<PixelType>
3098  ul(img.data(), img.stride(1));
3099  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
3100  return triple<ConstImageIterator<PixelType>,
3101  ConstImageIterator<PixelType>,
3102  Accessor>
3103  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
3104 }
3105 
3106 template <class PixelType>
3107 inline pair< ConstStridedImageIterator<PixelType>,
3108  typename AccessorTraits<PixelType>::default_const_accessor>
3109 srcImage(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
3110 {
3111  ConstStridedImageIterator<PixelType>
3112  ul(img.data(), 1, img.stride(0), img.stride(1));
3113  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
3114  return pair<ConstStridedImageIterator<PixelType>,
3115  Accessor>
3116  (ul, Accessor());
3117 }
3118 
3119 template <class PixelType>
3120 inline pair< ConstImageIterator<PixelType>,
3121  typename AccessorTraits<PixelType>::default_const_accessor>
3122 srcImage(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
3123 {
3124  ConstImageIterator<PixelType>
3125  ul(img.data(), img.stride(1));
3126  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
3127  return pair<ConstImageIterator<PixelType>,
3128  Accessor>
3129  (ul, Accessor());
3130 }
3131 
3132 template <class PixelType>
3133 inline triple< StridedImageIterator<PixelType>,
3134  StridedImageIterator<PixelType>,
3135  typename AccessorTraits<PixelType>::default_accessor>
3136 destImageRange(MultiArrayView<2, PixelType, StridedArrayTag> & img)
3137 {
3138  StridedImageIterator<PixelType>
3139  ul(img.data(), 1, img.stride(0), img.stride(1));
3140  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3141  return triple<StridedImageIterator<PixelType>,
3142  StridedImageIterator<PixelType>,
3143  Accessor>
3144  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
3145 }
3146 
3147 template <class PixelType>
3148 inline triple< ImageIterator<PixelType>,
3149  ImageIterator<PixelType>,
3150  typename AccessorTraits<PixelType>::default_accessor>
3151 destImageRange(MultiArrayView<2, PixelType, UnstridedArrayTag> & img)
3152 {
3153  ImageIterator<PixelType>
3154  ul(img.data(), img.stride(1));
3155  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3156  return triple<ImageIterator<PixelType>,
3157  ImageIterator<PixelType>,
3158  Accessor>
3159  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
3160 }
3161 
3162 template <class PixelType>
3163 inline pair< StridedImageIterator<PixelType>,
3164  typename AccessorTraits<PixelType>::default_accessor>
3165 destImage(MultiArrayView<2, PixelType, StridedArrayTag> & img)
3166 {
3167  StridedImageIterator<PixelType>
3168  ul(img.data(), 1, img.stride(0), img.stride(1));
3169  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3170  return pair<StridedImageIterator<PixelType>, Accessor>
3171  (ul, Accessor());
3172 }
3173 
3174 template <class PixelType>
3175 inline pair< ImageIterator<PixelType>,
3176  typename AccessorTraits<PixelType>::default_accessor>
3177 destImage(MultiArrayView<2, PixelType, UnstridedArrayTag> & img)
3178 {
3179  ImageIterator<PixelType> ul(img.data(), img.stride(1));
3180  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3181  return pair<ImageIterator<PixelType>, Accessor>(ul, Accessor());
3182 }
3183 
3184 template <class PixelType>
3185 inline pair< ConstStridedImageIterator<PixelType>,
3186  typename AccessorTraits<PixelType>::default_accessor>
3187 maskImage(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
3188 {
3189  ConstStridedImageIterator<PixelType>
3190  ul(img.data(), 1, img.stride(0), img.stride(1));
3191  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3192  return pair<ConstStridedImageIterator<PixelType>, Accessor>
3193  (ul, Accessor());
3194 }
3195 
3196 template <class PixelType>
3197 inline pair< ConstImageIterator<PixelType>,
3198  typename AccessorTraits<PixelType>::default_accessor>
3199 maskImage(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
3200 {
3201  ConstImageIterator<PixelType>
3202  ul(img.data(), img.stride(1));
3203  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3204  return pair<ConstImageIterator<PixelType>, Accessor>
3205  (ul, Accessor());
3206 }
3207 
3208 /********************************************************/
3209 /* */
3210 /* makeBasicImageView */
3211 /* */
3212 /********************************************************/
3213 
3214 /** \addtogroup MultiArrayToImage Wrap a \ref vigra::MultiArrayView in
3215  a \ref vigra::BasicImageView
3216 */
3217 //@{
3218 /** Create a \ref vigra::BasicImageView from an unstrided 2-dimensional
3219  \ref vigra::MultiArrayView.
3220 
3221  The \ref vigra::BasicImageView will have the same <tt>value_type </tt>
3222  as the original \ref vigra::MultiArrayView.
3223 */
3224 template <class T>
3225 BasicImageView <T>
3227 {
3228  return BasicImageView <T> (array.data (), array.shape (0),
3229  array.shape (1));
3230 }
3231 
3232 /** Create a \ref vigra::BasicImageView from a 3-dimensional
3233  \ref vigra::MultiArray.
3234 
3235  This wrapper flattens the two innermost dimensions of the array
3236  into single rows of the resulting image.
3237  The \ref vigra::BasicImageView will have the same <tt>value_type </tt>
3238  as the original \ref vigra::MultiArray.
3239 */
3240 template <class T>
3241 BasicImageView <T>
3243 {
3244  return BasicImageView <T> (array.data (),
3245  array.shape (0)*array.shape (1), array.shape (2));
3246 }
3247 
3248 /** Create a \ref vigra::BasicImageView from a 3-dimensional
3249  \ref vigra::MultiArray.
3250 
3251  This wrapper only works if <tt>T</tt> is a scalar type and the
3252  array's innermost dimension has size 3. It then re-interprets
3253  the data array as a 2-dimensional array with value_type
3254  <tt>RGBValue<T></tt>.
3255 */
3256 template <class T>
3257 BasicImageView <RGBValue<T> >
3259 {
3260  vigra_precondition (
3261  array.shape (0) == 3, "makeRGBImageView(): array.shape(0) must be 3.");
3262  return BasicImageView <RGBValue<T> > (
3263  reinterpret_cast <RGBValue <T> *> (array.data ()),
3264  array.shape (1), array.shape (2));
3265 }
3266 
3267 //@}
3268 
3269 } // namespace vigra
3270 #undef VIGRA_ASSERT_INSIDE
3271 #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.8.0 (Wed Sep 26 2012)