[ 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 
149 template <class C>
150 struct CoordinatesToOffest
151 {
152  template <int N>
153  static MultiArrayIndex
154  exec(const TinyVector <MultiArrayIndex, N> & stride, MultiArrayIndex x)
155  {
156  return stride[0] * x;
157  }
158  template <int N>
159  static MultiArrayIndex
160  exec(const TinyVector <MultiArrayIndex, N> & stride, MultiArrayIndex x, MultiArrayIndex y)
161  {
162  return stride[0] * x + stride[1] * y;
163  }
164 };
165 
166 template <>
167 struct CoordinatesToOffest<UnstridedArrayTag>
168 {
169  template <int N>
170  static MultiArrayIndex
171  exec(const TinyVector <MultiArrayIndex, N> & /*stride*/, MultiArrayIndex x)
172  {
173  return x;
174  }
175  template <int N>
176  static MultiArrayIndex
177  exec(const TinyVector <MultiArrayIndex, N> & stride, MultiArrayIndex x, MultiArrayIndex y)
178  {
179  return x + stride[1] * y;
180  }
181 };
182 
183 /********************************************************/
184 /* */
185 /* MaybeStrided */
186 /* */
187 /********************************************************/
188 
189 /* metatag implementing a test for marking MultiArrays that were
190  indexed at the zero'th dimension as strided, and all others as
191  unstrided.
192 
193 <b>\#include</b>
194 <vigra/multi_array.hxx>
195 
196 Namespace: vigra::detail
197 */
198 template <class StrideTag, unsigned int N>
199 struct MaybeStrided
200 {
201  typedef StrideTag type;
202 };
203 
204 template <class StrideTag>
205 struct MaybeStrided <StrideTag, 0>
206 {
207  typedef StridedArrayTag type;
208 };
209 
210 /********************************************************/
211 /* */
212 /* MultiIteratorChooser */
213 /* */
214 /********************************************************/
215 
216 /* metatag implementing a test (by pattern matching) for marking
217  MultiArrays that were indexed at the zero'th dimension as strided.
218 
219 <b>\#include</b>
220 <vigra/multi_array.hxx>
221 
222 Namespace: vigra::detail
223 */
224 template <class O>
225 struct MultiIteratorChooser
226 {
227  struct Nil {};
228 
229  template <unsigned int N, class T, class REFERENCE, class POINTER>
230  struct Traverser
231  {
232  typedef Nil type;
233  };
234 };
235 
236 /********************************************************/
237 /* */
238 /* MultiIteratorChooser <StridedArrayTag> */
239 /* */
240 /********************************************************/
241 
242 /* specialization of the MultiIteratorChooser for strided arrays.
243 
244 <b>\#include</b>
245 <vigra/multi_array.hxx>
246 
247 Namespace: vigra::detail
248 */
249 template <>
250 struct MultiIteratorChooser <StridedArrayTag>
251 {
252  template <unsigned int N, class T, class REFERENCE, class POINTER>
253  struct Traverser
254  {
255  typedef StridedMultiIterator <N, T, REFERENCE, POINTER> type;
256  };
257 };
258 
259 /********************************************************/
260 /* */
261 /* MultiIteratorChooser <UnstridedArrayTag> */
262 /* */
263 /********************************************************/
264 
265 /* specialization of the MultiIteratorChooser for unstrided arrays.
266 
267 <b>\#include</b>
268 <vigra/multi_array.hxx>
269 
270 Namespace: vigra::detail
271 */
272 template <>
273 struct MultiIteratorChooser <UnstridedArrayTag>
274 {
275  template <unsigned int N, class T, class REFERENCE, class POINTER>
276  struct Traverser
277  {
278  typedef MultiIterator <N, T, REFERENCE, POINTER> type;
279  };
280 };
281 
282 /********************************************************/
283 /* */
284 /* helper functions */
285 /* */
286 /********************************************************/
287 
288 template <class DestIterator, class Shape, class T>
289 inline void
290 initMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<0>)
291 {
292  DestIterator dend = d + shape[0];
293  for(; d < dend; ++d)
294  {
295  *d = init;
296  }
297 }
298 
299 template <class DestIterator, class Shape, class T, int N>
300 void
301 initMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<N>)
302 {
303  DestIterator dend = d + shape[N];
304  for(; d < dend; ++d)
305  {
306  initMultiArrayData(d.begin(), shape, init, MetaInt<N-1>());
307  }
308 }
309 
310 // FIXME: the explicit overload for MultiIterator<1, UInt8, ... > works around a compiler crash in VisualStudio 2010
311 #define VIGRA_COPY_MULTI_ARRAY_DATA(name, op) \
312 template <class SrcIterator, class Shape, class DestIterator> \
313 inline void \
314 name##MultiArrayData(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>) \
315 { \
316  SrcIterator send = s + shape[0]; \
317  for(; s < send; ++s, ++d) \
318  { \
319  *d op detail::RequiresExplicitCast<typename DestIterator::value_type>::cast(*s); \
320  } \
321 } \
322  \
323 template <class Ref, class Ptr, class Shape, class DestIterator> \
324 inline void \
325 name##MultiArrayData(MultiIterator<1, UInt8, Ref, Ptr> si, Shape const & shape, DestIterator d, MetaInt<0>) \
326 { \
327  Ptr s = &(*si), send = s + shape[0]; \
328  for(; s < send; ++s, ++d) \
329  { \
330  *d op detail::RequiresExplicitCast<typename DestIterator::value_type>::cast(*s); \
331  } \
332 } \
333 \
334 template <class SrcIterator, class Shape, class DestIterator, int N> \
335 void \
336 name##MultiArrayData(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>) \
337 { \
338  SrcIterator send = s + shape[N]; \
339  for(; s < send; ++s, ++d) \
340  { \
341  name##MultiArrayData(s.begin(), shape, d.begin(), MetaInt<N-1>()); \
342  } \
343 } \
344 \
345 template <class DestIterator, class Shape, class T> \
346 inline void \
347 name##ScalarMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<0>) \
348 { \
349  DestIterator dend = d + shape[0]; \
350  for(; d < dend; ++d) \
351  { \
352  *d op detail::RequiresExplicitCast<typename DestIterator::value_type>::cast(init); \
353  } \
354 } \
355  \
356 template <class DestIterator, class Shape, class T, int N> \
357 void \
358 name##ScalarMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<N>) \
359 { \
360  DestIterator dend = d + shape[N]; \
361  for(; d < dend; ++d) \
362  { \
363  name##ScalarMultiArrayData(d.begin(), shape, init, MetaInt<N-1>()); \
364  } \
365 }
366 
367 VIGRA_COPY_MULTI_ARRAY_DATA(copy, =)
368 VIGRA_COPY_MULTI_ARRAY_DATA(copyAdd, +=)
369 VIGRA_COPY_MULTI_ARRAY_DATA(copySub, -=)
370 VIGRA_COPY_MULTI_ARRAY_DATA(copyMul, *=)
371 VIGRA_COPY_MULTI_ARRAY_DATA(copyDiv, /=)
372 
373 #undef VIGRA_COPY_MULTI_ARRAY_DATA
374 
375 template <class SrcIterator, class Shape, class T, class ALLOC>
376 inline void
377 uninitializedCopyMultiArrayData(SrcIterator s, Shape const & shape, T * & d, ALLOC & a, MetaInt<0>)
378 {
379  SrcIterator send = s + shape[0];
380  for(; s < send; ++s, ++d)
381  {
382  a.construct(d, static_cast<T const &>(*s));
383  }
384 }
385 
386 // FIXME: this overload works around a compiler crash in VisualStudio 2010
387 template <class Ref, class Ptr, class Shape, class T, class ALLOC>
388 inline void
389 uninitializedCopyMultiArrayData(MultiIterator<1, UInt8, Ref, Ptr> si, Shape const & shape, T * & d, ALLOC & a, MetaInt<0>)
390 {
391  Ptr s = &(*si), 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>
594 template <unsigned int N, class T, class A = std::allocator<T> >
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>
766 class MultiArrayView
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  /** Assignment of a scalar. Equivalent to MultiArrayView::init(v).
943  */
944  MultiArrayView & operator=(value_type const & v)
945  {
946  return init(v);
947  }
948 
949  /** Add-assignment of a compatible MultiArrayView. Fails with
950  <tt>PreconditionViolation</tt> exception when the shapes do not match.
951  */
952  template<class U, class C1>
954 
955  /** Subtract-assignment of a compatible MultiArrayView. Fails with
956  <tt>PreconditionViolation</tt> exception when the shapes do not match.
957  */
958  template<class U, class C1>
960 
961  /** Multiply-assignment of a compatible MultiArrayView. Fails with
962  <tt>PreconditionViolation</tt> exception when the shapes do not match.
963  */
964  template<class U, class C1>
966 
967  /** Divide-assignment of a compatible MultiArrayView. Fails with
968  <tt>PreconditionViolation</tt> exception when the shapes do not match.
969  */
970  template<class U, class C1>
972 
973  /** Add-assignment of a scalar.
974  */
975  MultiArrayView & operator+=(T const & rhs)
976  {
977  detail::copyAddScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
978  return *this;
979  }
980 
981  /** Subtract-assignment of a scalar.
982  */
983  MultiArrayView & operator-=(T const & rhs)
984  {
985  detail::copySubScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
986  return *this;
987  }
988 
989  /** Multiply-assignment of a scalar.
990  */
991  MultiArrayView & operator*=(T const & rhs)
992  {
993  detail::copyMulScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
994  return *this;
995  }
996 
997  /** Divide-assignment of a scalar.
998  */
999  MultiArrayView & operator/=(T const & rhs)
1000  {
1001  detail::copyDivScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
1002  return *this;
1003  }
1004 
1005  /** Assignment of an array expression. Fails with
1006  <tt>PreconditionViolation</tt> exception when the shapes do not match.
1007  */
1008  template<class Expression>
1009  MultiArrayView & operator=(multi_math::MultiMathOperand<Expression> const & rhs)
1010  {
1011  multi_math::detail::assign(*this, rhs);
1012  return *this;
1013  }
1014 
1015  /** Add-assignment of an array expression. Fails with
1016  <tt>PreconditionViolation</tt> exception when the shapes do not match.
1017  */
1018  template<class Expression>
1019  MultiArrayView & operator+=(multi_math::MultiMathOperand<Expression> const & rhs)
1020  {
1021  multi_math::detail::plusAssign(*this, rhs);
1022  return *this;
1023  }
1024 
1025  /** Subtract-assignment of an array expression. Fails with
1026  <tt>PreconditionViolation</tt> exception when the shapes do not match.
1027  */
1028  template<class Expression>
1029  MultiArrayView & operator-=(multi_math::MultiMathOperand<Expression> const & rhs)
1030  {
1031  multi_math::detail::minusAssign(*this, rhs);
1032  return *this;
1033  }
1034 
1035  /** Multiply-assignment of an array expression. Fails with
1036  <tt>PreconditionViolation</tt> exception when the shapes do not match.
1037  */
1038  template<class Expression>
1039  MultiArrayView & operator*=(multi_math::MultiMathOperand<Expression> const & rhs)
1040  {
1041  multi_math::detail::multiplyAssign(*this, rhs);
1042  return *this;
1043  }
1044 
1045  /** Divide-assignment of an array expression. Fails with
1046  <tt>PreconditionViolation</tt> exception when the shapes do not match.
1047  */
1048  template<class Expression>
1049  MultiArrayView & operator/=(multi_math::MultiMathOperand<Expression> const & rhs)
1050  {
1051  multi_math::detail::divideAssign(*this, rhs);
1052  return *this;
1053  }
1054 
1055  /** array access.
1056  */
1057  reference operator[] (const difference_type &d)
1058  {
1059  VIGRA_ASSERT_INSIDE(d);
1060  return m_ptr [dot (d, m_stride)];
1061  }
1062 
1063  /** array access.
1064  */
1065  const_reference operator[] (const difference_type &d) const
1066  {
1067  VIGRA_ASSERT_INSIDE(d);
1068  return m_ptr [dot (d, m_stride)];
1069  }
1070 
1071  /** equivalent to bindInner(), when M < N.
1072  */
1073  template <int M>
1075  {
1076  return bindInner(d);
1077  }
1078 
1079  /** Array access in scan-order sense.
1080  Mostly useful to support standard indexing for 1-dimensional multi-arrays,
1081  but works for any N. Use scanOrderIndexToCoordinate() and
1082  coordinateToScanOrderIndex() for conversion between indices and coordinates.
1083 
1084  <b>Note:</b> This function should not be used in the inner loop, because the
1085  conversion of the scan order index into a memory address is expensive
1086  (it must take into account that memory may not be consecutive for subarrays
1087  and/or strided arrays). Always prefer operator() if possible.
1088  */
1089  reference operator[](difference_type_1 d)
1090  {
1091  VIGRA_ASSERT_INSIDE(scanOrderIndexToCoordinate(d));
1092  return m_ptr [detail::ScanOrderToOffset<actual_dimension>::exec(d, m_shape, m_stride)];
1093  }
1094 
1095  /** Array access in scan-order sense.
1096  Mostly useful to support standard indexing for 1-dimensional multi-arrays,
1097  but works for any N. Use scanOrderIndexToCoordinate() and
1098  coordinateToScanOrderIndex() for conversion between indices and coordinates.
1099 
1100  <b>Note:</b> This function should not be used in the inner loop, because the
1101  conversion of the scan order index into a memory address is expensive
1102  (it must take into account that memory may not be consecutive for subarrays
1103  and/or strided arrays). Always prefer operator() if possible.
1104  */
1105  const_reference operator[](difference_type_1 d) const
1106  {
1107  VIGRA_ASSERT_INSIDE(scanOrderIndexToCoordinate(d));
1108  return m_ptr [detail::ScanOrderToOffset<actual_dimension>::exec(d, m_shape, m_stride)];
1109  }
1110 
1111  /** convert scan-order index to coordinate.
1112  */
1113  difference_type scanOrderIndexToCoordinate(difference_type_1 d) const
1114  {
1115  difference_type result;
1116  detail::ScanOrderToCoordinate<actual_dimension>::exec(d, m_shape, result);
1117  return result;
1118  }
1119 
1120  /** convert coordinate to scan-order index.
1121  */
1122  difference_type_1 coordinateToScanOrderIndex(const difference_type &d) const
1123  {
1124  return detail::CoordinateToScanOrder<actual_dimension>::exec(m_shape, d);
1125  }
1126 
1127  /** 1D array access. Use only if N == 1.
1128  */
1129  reference operator() (difference_type_1 x)
1130  {
1131  VIGRA_ASSERT_INSIDE(difference_type(x));
1132  return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x)];
1133  }
1134 
1135  /** 2D array access. Use only if N == 2.
1136  */
1137  reference operator() (difference_type_1 x, difference_type_1 y)
1138  {
1139  VIGRA_ASSERT_INSIDE(difference_type(x, y));
1140  return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x, y)];
1141  }
1142 
1143  /** 3D array access. Use only if N == 3.
1144  */
1145  reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z)
1146  {
1147  VIGRA_ASSERT_INSIDE(difference_type(x, y, z));
1148  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
1149  }
1150 
1151  /** 4D array access. Use only if N == 4.
1152  */
1153  reference operator() (difference_type_1 x, difference_type_1 y,
1154  difference_type_1 z, difference_type_1 u)
1155  {
1156  VIGRA_ASSERT_INSIDE(difference_type(x, y, z, u));
1157  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
1158  }
1159 
1160  /** 5D array access. Use only if N == 5.
1161  */
1162  reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
1163  difference_type_1 u, difference_type_1 v)
1164  {
1165  VIGRA_ASSERT_INSIDE(difference_type(x, y,z, u,v));
1166  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
1167  }
1168 
1169  /** 1D const array access. Use only if N == 1.
1170  */
1171  const_reference operator() (difference_type_1 x) const
1172  {
1173  VIGRA_ASSERT_INSIDE(difference_type(x));
1174  return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x)];
1175  }
1176 
1177  /** 2D const array access. Use only if N == 2.
1178  */
1179  const_reference operator() (difference_type_1 x, difference_type_1 y) const
1180  {
1181  VIGRA_ASSERT_INSIDE(difference_type(x, y));
1182  return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x, y)];
1183  }
1184 
1185  /** 3D const array access. Use only if N == 3.
1186  */
1187  const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z) const
1188  {
1189  VIGRA_ASSERT_INSIDE(difference_type(x,y,z));
1190  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
1191  }
1192 
1193  /** 4D const array access. Use only if N == 4.
1194  */
1195  const_reference operator() (difference_type_1 x, difference_type_1 y,
1196  difference_type_1 z, difference_type_1 u) const
1197  {
1198  VIGRA_ASSERT_INSIDE(difference_type(x,y,z,u));
1199  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
1200  }
1201 
1202  /** 5D const array access. Use only if N == 5.
1203  */
1204  const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
1205  difference_type_1 u, difference_type_1 v) const
1206  {
1207  VIGRA_ASSERT_INSIDE(difference_type(x,y,z,u,v));
1208  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
1209  }
1210 
1211  /** Init with a constant.
1212  */
1213  template <class U>
1214  MultiArrayView & init(const U & init)
1215  {
1216  if(hasData())
1217  detail::copyScalarMultiArrayData(traverser_begin(), shape(), init, MetaInt<actual_dimension-1>());
1218  return *this;
1219  }
1220 
1221 
1222  /** Copy the data of the right-hand array (array shapes must match).
1223  */
1224  void copy(const MultiArrayView & rhs)
1225  {
1226  if(this == &rhs)
1227  return;
1228  this->copyImpl(rhs);
1229  }
1230 
1231  /** Copy the data of the right-hand array (array shapes must match).
1232  */
1233  template <class U, class CN>
1235  {
1236  this->copyImpl(rhs);
1237  }
1238 
1239  /** swap the data between two MultiArrayView objects.
1240 
1241  The shapes of the two array must match.
1242  */
1244  {
1245  if(this != &rhs)
1246  swapDataImpl(rhs);
1247  }
1248 
1249  /** swap the data between two MultiArrayView objects.
1250 
1251  The shapes of the two array must match.
1252  */
1253  template <class T2, class C2>
1255  {
1256  swapDataImpl(rhs);
1257  }
1258 
1259  /** check whether the array is unstrided (i.e. has consecutive memory) up
1260  to the given dimension.
1261 
1262  \a dimension can range from 0 ... N-1. If a certain dimension is unstrided,
1263  all lower dimensions are also unstrided.
1264  */
1265  bool isUnstrided(unsigned int dimension = N-1) const
1266  {
1267  difference_type p = shape() - difference_type(1);
1268  for(unsigned int k = dimension+1; k < N; ++k)
1269  p[k] = 0;
1270  return (&operator[](p) - m_ptr) == coordinateToScanOrderIndex(p);
1271  }
1272 
1273  /** bind the M outmost dimensions to certain indices.
1274  this reduces the dimensionality of the image to
1275  max { 1, N-M }.
1276 
1277  <b>Usage:</b>
1278  \code
1279  // create a 3D array of size 40x30x20
1280  typedef MultiArray<3, double>::difference_type Shape;
1281  MultiArray<3, double> array3(Shape(40, 30, 20));
1282 
1283  // get a 1D array by fixing index 1 to 12, and index 2 to 10
1284  MultiArrayView <1, double> array1 = array3.bindOuter(TinyVector<MultiArrayIndex, 2>(12, 10));
1285  \endcode
1286  */
1287  template <int M, class Index>
1288  MultiArrayView <N-M, T, StrideTag> bindOuter(const TinyVector <Index, M> &d) const;
1289 
1290  /** bind the M innermost dimensions to certain indices.
1291  this reduces the dimensionality of the image to
1292  max { 1, N-M }.
1293 
1294  <b>Usage:</b>
1295  \code
1296  // create a 3D array of size 40x30x20
1297  typedef MultiArray<3, double>::difference_type Shape;
1298  MultiArray<3, double> array3(Shape(40, 30, 20));
1299 
1300  // get a 1D array by fixing index 0 to 12, and index 1 to 10
1301  MultiArrayView <1, double, StridedArrayTag> array1 = array3.bindInner(TinyVector<MultiArrayIndex, 2>(12, 10));
1302  \endcode
1303  */
1304  template <int M, class Index>
1306 
1307  /** bind dimension M to index d.
1308  this reduces the dimensionality of the image to
1309  max { 1, N-1 }.
1310 
1311  <b>Usage:</b>
1312  \code
1313  // create a 3D array of size 40x30x20
1314  typedef MultiArray<3, double>::difference_type Shape;
1315  MultiArray<3, double> array3(Shape(40, 30, 20));
1316 
1317  // get a 2D array by fixing index 1 to 12
1318  MultiArrayView <2, double> array2 = array3.bind<1>(12);
1319 
1320  // get a 2D array by fixing index 0 to 23
1321  MultiArrayView <2, double, StridedArrayTag> array2a = array3.bind<0>(23);
1322  \endcode
1323  */
1324  template <unsigned int M>
1325  MultiArrayView <N-1, T, typename vigra::detail::MaybeStrided<StrideTag, M>::type >
1326  bind (difference_type_1 d) const;
1327 
1328  /** bind the outmost dimension to a certain index.
1329  this reduces the dimensionality of the image to
1330  max { 1, N-1 }.
1331 
1332  <b>Usage:</b>
1333  \code
1334  // create a 3D array of size 40x30x20
1335  typedef MultiArray<3, double>::difference_type Shape;
1336  MultiArray<3, double> array3(Shape(40, 30, 20));
1337 
1338  // get a 2D array by fixing the outermost index (i.e. index 2) to 12
1339  MultiArrayView <2, double> array2 = array3.bindOuter(12);
1340  \endcode
1341  */
1342  MultiArrayView <N-1, T, StrideTag> bindOuter (difference_type_1 d) const;
1343 
1344  /** bind the innermost dimension to a certain index.
1345  this reduces the dimensionality of the image to
1346  max { 1, N-1 }.
1347 
1348  <b>Usage:</b>
1349  \code
1350  // create a 3D array of size 40x30x20
1351  typedef MultiArray<3, double>::difference_type Shape;
1352  MultiArray<3, double> array3(Shape(40, 30, 20));
1353 
1354  // get a 2D array by fixing the innermost index (i.e. index 0) to 23
1355  MultiArrayView <2, double, StridedArrayTag> array2 = array3.bindInner(23);
1356  \endcode
1357  */
1358  MultiArrayView <N-1, T, StridedArrayTag> bindInner (difference_type_1 d) const;
1359 
1360  /** bind dimension m to index d.
1361  this reduces the dimensionality of the image to
1362  max { 1, N-1 }.
1363 
1364  <b>Usage:</b>
1365  \code
1366  // create a 3D array of size 40x30x20
1367  typedef MultiArray<3, double>::difference_type Shape;
1368  MultiArray<3, double> array3(Shape(40, 30, 20));
1369 
1370  // get a 2D array by fixing index 2 to 15
1371  MultiArrayView <2, double, StridedArrayTag> array2 = array3.bindAt(2, 15);
1372  \endcode
1373  */
1375  bindAt (difference_type_1 m, difference_type_1 d) const;
1376 
1377 
1378  /** Create a view to channel 'i' of a vector-like value type. Possible value types
1379  (of the original array) are: \ref TinyVector, \ref RGBValue, \ref FFTWComplex,
1380  and <tt>std::complex</tt>. The list can be extended to any type whose memory
1381  layout is equivalent to a fixed-size C array, by specializing
1382  <tt>ExpandElementResult</tt>.
1383 
1384  <b>Usage:</b>
1385  \code
1386  MultiArray<2, RGBValue<float> > rgb_image(Shape2(w, h));
1387 
1388  MultiArrayView<2, float, StridedArrayTag> red = rgb_image.bindElementChannel(0);
1389  MultiArrayView<2, float, StridedArrayTag> green = rgb_image.bindElementChannel(1);
1390  MultiArrayView<2, float, StridedArrayTag> blue = rgb_image.bindElementChannel(2);
1391  \endcode
1392  */
1394  bindElementChannel(difference_type_1 i) const
1395  {
1396  vigra_precondition(0 <= i && i < ExpandElementResult<T>::size,
1397  "MultiArrayView::bindElementChannel(i): 'i' out of range.");
1398  return expandElements(0).bindInner(i);
1399  }
1400 
1401  /** Create a view where a vector-like element type is expanded into a new
1402  array dimension. The new dimension is inserted at index position 'd',
1403  which must be between 0 and N inclusive.
1404 
1405  Possible value types of the original array are: \ref TinyVector, \ref RGBValue,
1406  \ref FFTWComplex, <tt>std::complex</tt>, and the built-in number types (in this
1407  case, <tt>expandElements</tt> is equivalent to <tt>insertSingletonDimension</tt>).
1408  The list of supported types can be extended to any type whose memory
1409  layout is equivalent to a fixed-size C array, by specializing
1410  <tt>ExpandElementResult</tt>.
1411 
1412  <b>Usage:</b>
1413  \code
1414  MultiArray<2, RGBValue<float> > rgb_image(Shape2(w, h));
1415 
1416  MultiArrayView<3, float, StridedArrayTag> multiband_image = rgb_image.expandElements(2);
1417  \endcode
1418  */
1420  expandElements(difference_type_1 d) const;
1421 
1422  /** Add a singleton dimension (dimension of length 1).
1423 
1424  Singleton dimensions don't change the size of the data, but introduce
1425  a new index that can only take the value 0. This is mainly useful for
1426  the 'reduce mode' of transformMultiArray() and combineTwoMultiArrays(),
1427  because these functions require the source and destination arrays to
1428  have the same number of dimensions.
1429 
1430  The range of \a i must be <tt>0 <= i <= N</tt>. The new dimension will become
1431  the i'th index, and the old indices from i upwards will shift one
1432  place to the right.
1433 
1434  <b>Usage:</b>
1435 
1436  Suppose we want have a 2D array and want to create a 1D array that contains
1437  the row average of the first array.
1438  \code
1439  typedef MultiArrayShape<2>::type Shape2;
1440  MultiArray<2, double> original(Shape2(40, 30));
1441 
1442  typedef MultiArrayShape<1>::type Shape1;
1443  MultiArray<1, double> rowAverages(Shape1(30));
1444 
1445  // temporarily add a singleton dimension to the destination array
1446  transformMultiArray(srcMultiArrayRange(original),
1447  destMultiArrayRange(rowAverages.insertSingletonDimension(0)),
1448  FindAverage<double>());
1449  \endcode
1450  */
1452  insertSingletonDimension (difference_type_1 i) const;
1453 
1454  /** Create a view to the diagonal elements of the array.
1455 
1456  This produces a 1D array view whose size equals the size
1457  of the shortest dimension of the original array.
1458 
1459  <b>Usage:</b>
1460  \code
1461  // create a 3D array of size 40x30x20
1462  typedef MultiArray<3, double>::difference_type Shape;
1463  MultiArray<3, double> array3(Shape(40, 30, 20));
1464 
1465  // get a view to the diagonal elements
1466  MultiArrayView <1, double, StridedArrayTag> diagonal = array3.diagonal();
1467  assert(diagonal.shape(0) == 20);
1468  \endcode
1469  */
1471  {
1473  Shape1(vigra::sum(m_stride)), m_ptr);
1474  }
1475 
1476  /** create a rectangular subarray that spans between the
1477  points p and q, where p is in the subarray, q not.
1478 
1479  <b>Usage:</b>
1480  \code
1481  // create a 3D array of size 40x30x20
1482  typedef MultiArray<3, double>::difference_type Shape;
1483  MultiArray<3, double> array3(Shape(40, 30, 20));
1484 
1485  // get a subarray set is smaller by one element at all sides
1486  MultiArrayView <3, double> subarray = array3.subarray(Shape(1,1,1), Shape(39, 29, 19));
1487  \endcode
1488  */
1489  MultiArrayView subarray (const difference_type &p,
1490  const difference_type &q) const
1491  {
1492  const difference_type_1 offset = dot (m_stride, p);
1493  return MultiArrayView (q - p, m_stride, m_ptr + offset);
1494  }
1495 
1496  /** apply an additional striding to the image, thereby reducing
1497  the shape of the array.
1498  for example, multiplying the stride of dimension one by three
1499  turns an appropriately laid out (interleaved) rgb image into
1500  a single band image.
1501  */
1503  stridearray (const difference_type &s) const
1504  {
1505  difference_type shape = m_shape;
1506  for (unsigned int i = 0; i < actual_dimension; ++i)
1507  shape [i] /= s [i];
1508  return MultiArrayView <N, T, StridedArrayTag>(shape, m_stride * s, m_ptr);
1509  }
1510 
1511  /** Transpose an array. If N==2, this implements the usual matrix transposition.
1512  For N > 2, it reverses the order of the indices.
1513 
1514  <b>Usage:</b><br>
1515  \code
1516  typedef MultiArray<2, double>::difference_type Shape;
1517  MultiArray<2, double> array(10, 20);
1518 
1519  MultiArray<2, double, StridedArrayTag> transposed = array.transpose();
1520 
1521  for(int i=0; i<array.shape(0), ++i)
1522  for(int j=0; j<array.shape(1); ++j)
1523  assert(array(i, j) == transposed(j, i));
1524  \endcode
1525  */
1527  transpose () const
1528  {
1529  difference_type shape(m_shape.begin(), difference_type::ReverseCopy),
1530  stride(m_stride.begin(), difference_type::ReverseCopy);
1532  }
1533 
1534  /** permute the dimensions of the array.
1535  The function exchanges the meaning of the dimensions without copying the data.
1536  In case of a 2-dimensional array, this is simply array transposition. In higher dimensions,
1537  there are more possibilities.
1538 
1539  <b>Usage:</b><br>
1540  \code
1541  typedef MultiArray<2, double>::difference_type Shape;
1542  MultiArray<2, double> array(10, 20);
1543 
1544  MultiArray<2, double, StridedArrayTag> transposed = array.permuteDimensions(Shape(1,0));
1545 
1546  for(int i=0; i<array.shape(0), ++i)
1547  for(int j=0; j<array.shape(1); ++j)
1548  assert(array(i, j) == transposed(j, i));
1549  \endcode
1550  */
1552  permuteDimensions (const difference_type &s) const;
1553 
1554  /** Permute the dimensions of the array so that the strides are in ascending order.
1555  Determines the appropriate permutation and then calls permuteDimensions().
1556  */
1558  permuteStridesAscending() const;
1559 
1560  /** Permute the dimensions of the array so that the strides are in descending order.
1561  Determines the appropriate permutation and then calls permuteDimensions().
1562  */
1564  permuteStridesDescending() const;
1565 
1566  /** Compute the ordering of the strides in this array.
1567  The result is describes the current permutation of the axes relative
1568  to the standard ascending stride order.
1569  */
1570  difference_type strideOrdering() const
1571  {
1572  return strideOrdering(m_stride);
1573  }
1574 
1575  /** Compute the ordering of the given strides.
1576  The result is describes the current permutation of the axes relative
1577  to the standard ascending stride order.
1578  */
1579  static difference_type strideOrdering(difference_type strides);
1580 
1581  /** number of the elements in the array.
1582  */
1583  difference_type_1 elementCount () const
1584  {
1585  difference_type_1 ret = m_shape[0];
1586  for(int i = 1; i < actual_dimension; ++i)
1587  ret *= m_shape[i];
1588  return ret;
1589  }
1590 
1591  /** number of the elements in the array.
1592  Same as <tt>elementCount()</tt>. Mostly useful to support the std::vector interface.
1593  */
1594  difference_type_1 size () const
1595  {
1596  return elementCount();
1597  }
1598 
1599  /** return the array's shape.
1600  */
1601  const difference_type & shape () const
1602  {
1603  return m_shape;
1604  }
1605 
1606  /** return the array's size at a certain dimension.
1607  */
1608  difference_type_1 size (difference_type_1 n) const
1609  {
1610  return m_shape [n];
1611  }
1612 
1613  /** return the array's shape at a certain dimension
1614  (same as <tt>size(n)</tt>).
1615  */
1616  difference_type_1 shape (difference_type_1 n) const
1617  {
1618  return m_shape [n];
1619  }
1620 
1621  /** return the array's stride for every dimension.
1622  */
1623  const difference_type & stride () const
1624  {
1625  return m_stride;
1626  }
1627 
1628  /** return the array's stride at a certain dimension.
1629  */
1630  difference_type_1 stride (int n) const
1631  {
1632  return m_stride [n];
1633  }
1634 
1635  /** check whether two arrays are elementwise equal.
1636  */
1637  template <class U, class C1>
1638  bool operator==(MultiArrayView<N, U, C1> const & rhs) const
1639  {
1640  if(this->shape() != rhs.shape())
1641  return false;
1642  return detail::equalityOfMultiArrays(traverser_begin(), shape(), rhs.traverser_begin(), MetaInt<actual_dimension-1>());
1643  }
1644 
1645  /** check whether two arrays are not elementwise equal.
1646  Also true when the two arrays have different shapes.
1647  */
1648  template <class U, class C1>
1649  bool operator!=(MultiArrayView<N, U, C1> const & rhs) const
1650  {
1651  return !operator==(rhs);
1652  }
1653 
1654  /** check whether the given point is in the array range.
1655  */
1656  bool isInside (difference_type const & p) const
1657  {
1658  for(int d=0; d<actual_dimension; ++d)
1659  if(p[d] < 0 || p[d] >= shape(d))
1660  return false;
1661  return true;
1662  }
1663 
1664  /** Check if the array contains only non-zero elements (or if all elements
1665  are 'true' if the value type is 'bool').
1666  */
1667  bool all() const
1668  {
1669  bool res = true;
1670  detail::reduceOverMultiArray(traverser_begin(), shape(),
1671  res,
1672  detail::AllTrueReduceFunctor(),
1673  MetaInt<actual_dimension-1>());
1674  return res;
1675  }
1676 
1677  /** Check if the array contains a non-zero element (or an element
1678  that is 'true' if the value type is 'bool').
1679  */
1680  bool any() const
1681  {
1682  bool res = false;
1683  detail::reduceOverMultiArray(traverser_begin(), shape(),
1684  res,
1685  detail::AnyTrueReduceFunctor(),
1686  MetaInt<actual_dimension-1>());
1687  return res;
1688  }
1689 
1690  /** Find the minimum and maximum element in this array.
1691  See \ref FeatureAccumulators for a general feature
1692  extraction framework.
1693  */
1694  void minmax(T * minimum, T * maximum) const
1695  {
1696  std::pair<T, T> res(NumericTraits<T>::max(), NumericTraits<T>::min());
1697  detail::reduceOverMultiArray(traverser_begin(), shape(),
1698  res,
1699  detail::MinmaxReduceFunctor(),
1700  MetaInt<actual_dimension-1>());
1701  *minimum = res.first;
1702  *maximum = res.second;
1703  }
1704 
1705  /** Compute the mean and variance of the values in this array.
1706  See \ref FeatureAccumulators for a general feature
1707  extraction framework.
1708  */
1709  template <class U>
1710  void meanVariance(U * mean, U * variance) const
1711  {
1712  typedef typename NumericTraits<U>::RealPromote R;
1713  R zero;
1714  triple<double, R, R> res(0.0, zero, zero);
1715  detail::reduceOverMultiArray(traverser_begin(), shape(),
1716  res,
1717  detail::MeanVarianceReduceFunctor(),
1718  MetaInt<actual_dimension-1>());
1719  *mean = res.second;
1720  *variance = res.third / res.first;
1721  }
1722 
1723  /** Compute the sum of the array elements.
1724 
1725  You must provide the type of the result by an explicit template parameter:
1726  \code
1727  MultiArray<2, UInt8> A(width, height);
1728 
1729  double sum = A.sum<double>();
1730  \endcode
1731  */
1732  template <class U>
1733  U sum() const
1734  {
1735  U res = NumericTraits<U>::zero();
1736  detail::reduceOverMultiArray(traverser_begin(), shape(),
1737  res,
1738  detail::SumReduceFunctor(),
1739  MetaInt<actual_dimension-1>());
1740  return res;
1741  }
1742 
1743  /** Compute the sum of the array elements over selected axes.
1744 
1745  \arg sums must have the same shape as this array, except for the
1746  axes along which the sum is to be accumulated. These axes must be
1747  singletons. Note that you must include <tt>multi_pointoperators.hxx</tt>
1748  for this function to work.
1749 
1750  <b>Usage:</b>
1751  \code
1752  #include <vigra/multi_array.hxx>
1753  #include <vigra/multi_pointoperators.hxx>
1754 
1755  MultiArray<2, double> A(rows, cols);
1756  ... // fill A
1757 
1758  // make the first axis a singleton to sum over the first index
1759  MultiArray<2, double> rowSums(1, cols);
1760  A.sum(rowSums);
1761 
1762  // this is equivalent to
1763  transformMultiArray(srcMultiArrayRange(A),
1764  destMultiArrayRange(rowSums),
1765  FindSum<double>());
1766  \endcode
1767  */
1768  template <class U, class S>
1769  void sum(MultiArrayView<N, U, S> sums) const
1770  {
1771  transformMultiArray(srcMultiArrayRange(*this),
1772  destMultiArrayRange(sums),
1773  FindSum<U>());
1774  }
1775 
1776  /** Compute the product of the array elements.
1777 
1778  You must provide the type of the result by an explicit template parameter:
1779  \code
1780  MultiArray<2, UInt8> A(width, height);
1781 
1782  double prod = A.product<double>();
1783  \endcode
1784  */
1785  template <class U>
1786  U product() const
1787  {
1788  U res = NumericTraits<U>::one();
1789  detail::reduceOverMultiArray(traverser_begin(), shape(),
1790  res,
1791  detail::ProdReduceFunctor(),
1792  MetaInt<actual_dimension-1>());
1793  return res;
1794  }
1795 
1796  /** Compute the squared Euclidean norm of the array (sum of squares of the array elements).
1797  */
1798  typename NormTraits<MultiArrayView>::SquaredNormType
1799  squaredNorm() const
1800  {
1801  typedef typename NormTraits<MultiArrayView>::SquaredNormType SquaredNormType;
1802  SquaredNormType res = NumericTraits<SquaredNormType>::zero();
1803  detail::reduceOverMultiArray(traverser_begin(), shape(),
1804  res,
1805  detail::SquaredL2NormReduceFunctor(),
1806  MetaInt<actual_dimension-1>());
1807  return res;
1808  }
1809 
1810  /** Compute various norms of the array.
1811  The norm is determined by parameter \a type:
1812 
1813  <ul>
1814  <li> type == 0: maximum norm (L-infinity): maximum of absolute values of the array elements
1815  <li> type == 1: Manhattan norm (L1): sum of absolute values of the array elements
1816  <li> type == 2: Euclidean norm (L2): square root of <tt>squaredNorm()</tt> when \a useSquaredNorm is <tt>true</tt>,<br>
1817  or direct algorithm that avoids underflow/overflow otherwise.
1818  </ul>
1819 
1820  Parameter \a useSquaredNorm has no effect when \a type != 2. Defaults: compute L2 norm as square root of
1821  <tt>squaredNorm()</tt>.
1822  */
1823  typename NormTraits<MultiArrayView>::NormType
1824  norm(int type = 2, bool useSquaredNorm = true) const;
1825 
1826  /** return the pointer to the image data
1827  */
1828  pointer data () const
1829  {
1830  return m_ptr;
1831  }
1832 
1833  pointer & unsafePtr()
1834  {
1835  return m_ptr;
1836  }
1837 
1838  /**
1839  * returns true iff this view refers to valid data,
1840  * i.e. data() is not a NULL pointer. (this is false after
1841  * default construction.)
1842  */
1843  bool hasData () const
1844  {
1845  return m_ptr != 0;
1846  }
1847 
1848  /** returns a scan-order iterator pointing
1849  to the first array element.
1850  */
1851  iterator begin()
1852  {
1853  return iterator(m_ptr, m_shape, m_stride);
1854  }
1855 
1856  /** returns a const scan-order iterator pointing
1857  to the first array element.
1858  */
1859  const_iterator begin() const
1860  {
1861  return const_iterator(m_ptr, m_shape, m_stride);
1862  }
1863 
1864  /** returns a scan-order iterator pointing
1865  beyond the last array element.
1866  */
1867  iterator end()
1868  {
1869  return begin().getEndIterator();
1870  }
1871 
1872  /** returns a const scan-order iterator pointing
1873  beyond the last array element.
1874  */
1875  const_iterator end() const
1876  {
1877  return begin().getEndIterator();
1878  }
1879 
1880  /** returns the N-dimensional MultiIterator pointing
1881  to the first element in every dimension.
1882  */
1883  traverser traverser_begin ()
1884  {
1885  traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1886  return ret;
1887  }
1888 
1889  /** returns the N-dimensional MultiIterator pointing
1890  to the const first element in every dimension.
1891  */
1892  const_traverser traverser_begin () const
1893  {
1894  const_traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1895  return ret;
1896  }
1897 
1898  /** returns the N-dimensional MultiIterator pointing
1899  beyond the last element in dimension N, and to the
1900  first element in every other dimension.
1901  */
1902  traverser traverser_end ()
1903  {
1904  traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1905  ret += m_shape [actual_dimension-1];
1906  return ret;
1907  }
1908 
1909  /** returns the N-dimensional const MultiIterator pointing
1910  beyond the last element in dimension N, and to the
1911  first element in every other dimension.
1912  */
1913  const_traverser traverser_end () const
1914  {
1915  const_traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1916  ret += m_shape [actual_dimension-1];
1917  return ret;
1918  }
1919 
1920  view_type view ()
1921  {
1922  return *this;
1923  }
1924 };
1925 
1926 template <unsigned int N, class T, class StrideTag>
1927 MultiArrayView<N, T, StrideTag> &
1929 {
1930  if(this == &rhs)
1931  return *this;
1932  vigra_precondition(this->shape() == rhs.shape() || m_ptr == 0,
1933  "MultiArrayView::operator=(MultiArrayView const &) size mismatch.");
1934  if(m_ptr == 0)
1935  {
1936  m_shape = rhs.m_shape;
1937  m_stride = rhs.m_stride;
1938  m_ptr = rhs.m_ptr;
1939  }
1940  else
1941  this->copyImpl(rhs);
1942  return *this;
1943 }
1944 
1945 template <unsigned int N, class T, class StrideTag>
1946 template <class CN>
1947 bool
1949 {
1950  vigra_precondition (shape () == rhs.shape (),
1951  "MultiArrayView::arraysOverlap(): shape mismatch.");
1952  const_pointer first_element = this->m_ptr,
1953  last_element = first_element + dot(this->m_shape - difference_type(1), this->m_stride);
1955  rhs_first_element = rhs.data(),
1956  rhs_last_element = rhs_first_element + dot(rhs.shape() - difference_type(1), rhs.stride());
1957  return !(last_element < rhs_first_element || rhs_last_element < first_element);
1958 }
1959 
1960 template <unsigned int N, class T, class StrideTag>
1961 template <class U, class CN>
1962 void
1963 MultiArrayView <N, T, StrideTag>::copyImpl(const MultiArrayView <N, U, CN>& rhs)
1964 {
1965  if(!arraysOverlap(rhs))
1966  {
1967  // no overlap -- can copy directly
1968  detail::copyMultiArrayData(rhs.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>());
1969  }
1970  else
1971  {
1972  // overlap: we got different views to the same data -- copy to intermediate memory in order to avoid
1973  // overwriting elements that are still needed on the rhs.
1974  MultiArray<N, T> tmp(rhs);
1975  detail::copyMultiArrayData(tmp.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>());
1976  }
1977 }
1978 
1979 #define VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(name, op) \
1980 template <unsigned int N, class T, class StrideTag> \
1981 template<class U, class C1> \
1982 MultiArrayView<N, T, StrideTag> & \
1983 MultiArrayView <N, T, StrideTag>::operator op(MultiArrayView<N, U, C1> const & rhs) \
1984 { \
1985  vigra_precondition(this->shape() == rhs.shape(), "MultiArrayView::operator" #op "() size mismatch."); \
1986  if(!arraysOverlap(rhs)) \
1987  { \
1988  detail::name##MultiArrayData(rhs.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>()); \
1989  } \
1990  else \
1991  { \
1992  MultiArray<N, T> tmp(rhs); \
1993  detail::name##MultiArrayData(tmp.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>()); \
1994  } \
1995  return *this; \
1996 }
1997 
1998 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copyAdd, +=)
1999 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copySub, -=)
2000 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copyMul, *=)
2001 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copyDiv, /=)
2002 
2003 #undef VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT
2004 
2005 template <unsigned int N, class T, class StrideTag>
2006 template <class U, class CN>
2007 void
2008 MultiArrayView <N, T, StrideTag>::swapDataImpl(MultiArrayView <N, U, CN> rhs)
2009 {
2010  vigra_precondition (shape () == rhs.shape (),
2011  "MultiArrayView::swapData(): shape mismatch.");
2012 
2013  // check for overlap of this and rhs
2014  const_pointer first_element = this->m_ptr,
2015  last_element = first_element + dot(this->m_shape - difference_type(1), this->m_stride);
2016  typename MultiArrayView <N, U, CN>::const_pointer
2017  rhs_first_element = rhs.data(),
2018  rhs_last_element = rhs_first_element + dot(rhs.shape() - difference_type(1), rhs.stride());
2019  if(last_element < rhs_first_element || rhs_last_element < first_element)
2020  {
2021  // no overlap -- can swap directly
2022  detail::swapDataImpl(traverser_begin(), shape(), rhs.traverser_begin(), MetaInt<actual_dimension-1>());
2023  }
2024  else
2025  {
2026  // overlap: we got different views to the same data -- copy to intermediate memory in order to avoid
2027  // overwriting elements that are still needed.
2028  MultiArray<N, T> tmp(*this);
2029  copy(rhs);
2030  rhs.copy(tmp);
2031  }
2032 }
2033 
2034 template <unsigned int N, class T, class StrideTag>
2035 MultiArrayView <N, T, StridedArrayTag>
2037 {
2038  difference_type shape, stride, check((typename difference_type::value_type)0);
2039  for (unsigned int i = 0; i < actual_dimension; ++i)
2040  {
2041  shape[i] = m_shape[s[i]];
2042  stride[i] = m_stride[s[i]];
2043  ++check[s[i]];
2044  }
2045  vigra_precondition(check == difference_type(1),
2046  "MultiArrayView::permuteDimensions(): every dimension must occur exactly once.");
2048 }
2049 
2050 template <unsigned int N, class T, class StrideTag>
2053 {
2054  difference_type permutation;
2055  for(int k=0; k<(int)N; ++k)
2056  permutation[k] = k;
2057  for(int k=0; k<(int)N-1; ++k)
2058  {
2059  int smallest = k;
2060  for(int j=k+1; j<(int)N; ++j)
2061  {
2062  if(stride[j] < stride[smallest])
2063  smallest = j;
2064  }
2065  if(smallest != k)
2066  {
2067  std::swap(stride[k], stride[smallest]);
2068  std::swap(permutation[k], permutation[smallest]);
2069  }
2070  }
2071  difference_type ordering;
2072  for(unsigned int k=0; k<N; ++k)
2073  ordering[permutation[k]] = k;
2074  return ordering;
2075 }
2076 
2077 template <unsigned int N, class T, class StrideTag>
2080 {
2081  difference_type ordering(strideOrdering(m_stride)), permutation;
2082  for(MultiArrayIndex k=0; k<N; ++k)
2083  permutation[ordering[k]] = k;
2084  return permuteDimensions(permutation);
2085 }
2086 
2087 template <unsigned int N, class T, class StrideTag>
2090 {
2091  difference_type ordering(strideOrdering(m_stride)), permutation;
2092  for(MultiArrayIndex k=0; k<N; ++k)
2093  permutation[ordering[N-1-k]] = k;
2094  return permuteDimensions(permutation);
2095 }
2096 
2097 template <unsigned int N, class T, class StrideTag>
2098 template <int M, class Index>
2099 MultiArrayView <N-M, T, StrideTag>
2101 {
2103  stride.init (m_stride.begin () + N-M, m_stride.end ());
2104  pointer ptr = m_ptr + dot (d, stride);
2105  static const int NNew = (N-M == 0) ? 1 : N-M;
2106  TinyVector <MultiArrayIndex, NNew> inner_shape, inner_stride;
2107  if (N-M == 0)
2108  {
2109  inner_shape [0] = 1;
2110  inner_stride [0] = 0;
2111  }
2112  else
2113  {
2114  inner_shape.init (m_shape.begin (), m_shape.end () - M);
2115  inner_stride.init (m_stride.begin (), m_stride.end () - M);
2116  }
2117  return MultiArrayView <N-M, T, StrideTag> (inner_shape, inner_stride, ptr);
2118 }
2119 
2120 template <unsigned int N, class T, class StrideTag>
2121 template <int M, class Index>
2122 MultiArrayView <N - M, T, StridedArrayTag>
2124 {
2126  stride.init (m_stride.begin (), m_stride.end () - N + M);
2127  pointer ptr = m_ptr + dot (d, stride);
2128  static const int NNew = (N-M == 0) ? 1 : N-M;
2129  TinyVector <MultiArrayIndex, NNew> outer_shape, outer_stride;
2130  if (N-M == 0)
2131  {
2132  outer_shape [0] = 1;
2133  outer_stride [0] = 0;
2134  }
2135  else
2136  {
2137  outer_shape.init (m_shape.begin () + M, m_shape.end ());
2138  outer_stride.init (m_stride.begin () + M, m_stride.end ());
2139  }
2140  return MultiArrayView <N-M, T, StridedArrayTag>
2141  (outer_shape, outer_stride, ptr);
2142 }
2143 
2144 template <unsigned int N, class T, class StrideTag>
2145 template <unsigned int M>
2146 MultiArrayView <N-1, T, typename detail::MaybeStrided<StrideTag, M>::type >
2148 {
2149  static const int NNew = (N-1 == 0) ? 1 : N-1;
2151  // the remaining dimensions are 0..n-1,n+1..N-1
2152  if (N-1 == 0)
2153  {
2154  shape[0] = 1;
2155  stride[0] = 0;
2156  }
2157  else
2158  {
2159  std::copy (m_shape.begin (), m_shape.begin () + M, shape.begin ());
2160  std::copy (m_shape.begin () + M+1, m_shape.end (),
2161  shape.begin () + M);
2162  std::copy (m_stride.begin (), m_stride.begin () + M, stride.begin ());
2163  std::copy (m_stride.begin () + M+1, m_stride.end (),
2164  stride.begin () + M);
2165  }
2166  return MultiArrayView <N-1, T, typename detail::MaybeStrided<StrideTag, M>::type>
2167  (shape, stride, m_ptr + d * m_stride[M]);
2168 }
2169 
2170 template <unsigned int N, class T, class StrideTag>
2171 MultiArrayView <N - 1, T, StrideTag>
2173 {
2174  static const int NNew = (N-1 == 0) ? 1 : N-1;
2175  TinyVector <MultiArrayIndex, NNew> inner_shape, inner_stride;
2176  if (N-1 == 0)
2177  {
2178  inner_shape [0] = 1;
2179  inner_stride [0] = 0;
2180  }
2181  else
2182  {
2183  inner_shape.init (m_shape.begin (), m_shape.end () - 1);
2184  inner_stride.init (m_stride.begin (), m_stride.end () - 1);
2185  }
2186  return MultiArrayView <N-1, T, StrideTag> (inner_shape, inner_stride,
2187  m_ptr + d * m_stride [N-1]);
2188 }
2189 
2190 template <unsigned int N, class T, class StrideTag>
2191 MultiArrayView <N - 1, T, StridedArrayTag>
2193 {
2194  static const int NNew = (N-1 == 0) ? 1 : N-1;
2195  TinyVector <MultiArrayIndex, NNew> outer_shape, outer_stride;
2196  if (N-1 == 0)
2197  {
2198  outer_shape [0] = 1;
2199  outer_stride [0] = 0;
2200  }
2201  else
2202  {
2203  outer_shape.init (m_shape.begin () + 1, m_shape.end ());
2204  outer_stride.init (m_stride.begin () + 1, m_stride.end ());
2205  }
2206  return MultiArrayView <N-1, T, StridedArrayTag>
2207  (outer_shape, outer_stride, m_ptr + d * m_stride [0]);
2208 }
2209 
2210 template <unsigned int N, class T, class StrideTag>
2211 MultiArrayView <N - 1, T, StridedArrayTag>
2213 {
2214  vigra_precondition (
2215  n < static_cast <int> (N),
2216  "MultiArrayView <N, T, StrideTag>::bindAt(): dimension out of range.");
2217  static const int NNew = (N-1 == 0) ? 1 : N-1;
2219  // the remaining dimensions are 0..n-1,n+1..N-1
2220  if (N-1 == 0)
2221  {
2222  shape [0] = 1;
2223  stride [0] = 0;
2224  }
2225  else
2226  {
2227  std::copy (m_shape.begin (), m_shape.begin () + n, shape.begin ());
2228  std::copy (m_shape.begin () + n+1, m_shape.end (),
2229  shape.begin () + n);
2230  std::copy (m_stride.begin (), m_stride.begin () + n, stride.begin ());
2231  std::copy (m_stride.begin () + n+1, m_stride.end (),
2232  stride.begin () + n);
2233  }
2234  return MultiArrayView <N-1, T, StridedArrayTag>
2235  (shape, stride, m_ptr + d * m_stride[n]);
2236 }
2237 
2238 
2239 template <unsigned int N, class T, class StrideTag>
2242 {
2243  vigra_precondition(0 <= d && d <= static_cast <difference_type_1> (N),
2244  "MultiArrayView<N, ...>::expandElements(d): 0 <= 'd' <= N required.");
2245 
2246  int elementSize = ExpandElementResult<T>::size;
2247  typename MultiArrayShape<N+1>::type newShape, newStrides;
2248  for(int k=0; k<d; ++k)
2249  {
2250  newShape[k] = m_shape[k];
2251  newStrides[k] = m_stride[k]*elementSize;
2252  }
2253 
2254  newShape[d] = elementSize;
2255  newStrides[d] = 1;
2256 
2257  for(int k=d; k<N; ++k)
2258  {
2259  newShape[k+1] = m_shape[k];
2260  newStrides[k+1] = m_stride[k]*elementSize;
2261  }
2262 
2263  typedef typename ExpandElementResult<T>::type U;
2265  newShape, newStrides, reinterpret_cast<U*>(m_ptr));
2266 }
2267 
2268 template <unsigned int N, class T, class StrideTag>
2271 {
2272  vigra_precondition (
2273  0 <= i && i <= static_cast <difference_type_1> (N),
2274  "MultiArrayView <N, T, StrideTag>::insertSingletonDimension(): index out of range.");
2276  std::copy (m_shape.begin (), m_shape.begin () + i, shape.begin ());
2277  std::copy (m_shape.begin () + i, m_shape.end (), shape.begin () + i + 1);
2278  std::copy (m_stride.begin (), m_stride.begin () + i, stride.begin ());
2279  std::copy (m_stride.begin () + i, m_stride.end (), stride.begin () + i + 1);
2280  shape[i] = 1;
2281  stride[i] = 1;
2282 
2284 }
2285 
2286 template <unsigned int N, class T, class StrideTag>
2287 typename NormTraits<MultiArrayView <N, T, StrideTag> >::NormType
2288 MultiArrayView <N, T, StrideTag>::norm(int type, bool useSquaredNorm) const
2289 {
2290  typedef typename NormTraits<MultiArrayView>::NormType NormType;
2291 
2292  switch(type)
2293  {
2294  case 0:
2295  {
2296  NormType res = NumericTraits<NormType>::zero();
2297  detail::reduceOverMultiArray(traverser_begin(), shape(),
2298  res,
2299  detail::MaxNormReduceFunctor(),
2300  MetaInt<actual_dimension-1>());
2301  return res;
2302  }
2303  case 1:
2304  {
2305  NormType res = NumericTraits<NormType>::zero();
2306  detail::reduceOverMultiArray(traverser_begin(), shape(),
2307  res,
2308  detail::L1NormReduceFunctor(),
2309  MetaInt<actual_dimension-1>());
2310  return res;
2311  }
2312  case 2:
2313  {
2314  if(useSquaredNorm)
2315  {
2316  return sqrt((NormType)squaredNorm());
2317  }
2318  else
2319  {
2320  NormType normMax = NumericTraits<NormType>::zero();
2321  detail::reduceOverMultiArray(traverser_begin(), shape(),
2322  normMax,
2323  detail::MaxNormReduceFunctor(),
2324  MetaInt<actual_dimension-1>());
2325  if(normMax == NumericTraits<NormType>::zero())
2326  return normMax;
2327  NormType res = NumericTraits<NormType>::zero();
2328  detail::reduceOverMultiArray(traverser_begin(), shape(),
2329  res,
2330  detail::WeightedL2NormReduceFunctor<NormType>(1.0/normMax),
2331  MetaInt<actual_dimension-1>());
2332  return sqrt(res)*normMax;
2333  }
2334  }
2335  default:
2336  vigra_precondition(false, "MultiArrayView::norm(): Unknown norm type.");
2337  return NumericTraits<NormType>::zero(); // unreachable
2338  }
2339 }
2340 
2341 
2342 /********************************************************/
2343 /* */
2344 /* norm */
2345 /* */
2346 /********************************************************/
2347 
2348 template <unsigned int N, class T, class StrideTag>
2349 inline typename NormTraits<MultiArrayView <N, T, StrideTag> >::SquaredNormType
2351 {
2352  return a.squaredNorm();
2353 }
2354 
2355 template <unsigned int N, class T, class StrideTag>
2356 inline typename NormTraits<MultiArrayView <N, T, StrideTag> >::NormType
2357 norm(MultiArrayView <N, T, StrideTag> const & a)
2358 {
2359  return a.norm();
2360 }
2361 
2362 /********************************************************/
2363 /* */
2364 /* MultiArray */
2365 /* */
2366 /********************************************************/
2367 
2368 /** \brief Main <TT>MultiArray</TT> class containing the memory
2369  management.
2370 
2371 This class inherits the interface of MultiArrayView, and implements
2372 the memory ownership.
2373 MultiArray's are always unstrided, striding them creates a MultiArrayView.
2374 
2375 
2376 The template parameters are as follows
2377 \code
2378  N: the array dimension
2379 
2380  T: the type of the array elements
2381 
2382  A: the allocator used for internal storage management
2383  (default: std::allocator<T>)
2384 \endcode
2385 
2386 <b>\#include</b>
2387 <vigra/multi_array.hxx>
2388 
2389 Namespace: vigra
2390 */
2391 template <unsigned int N, class T, class A /* default already declared above */>
2392 class MultiArray : public MultiArrayView <N, T>
2393 {
2394 
2395 public:
2396  using MultiArrayView <N, T>::actual_dimension;
2397 
2398  /** the allocator type used to allocate the memory
2399  */
2400  typedef A allocator_type;
2401 
2402  /** the view type associated with this array.
2403  */
2405 
2406  /** the matrix type associated with this array.
2407  */
2409 
2410  /** the array's value type
2411  */
2413 
2414  /** pointer type
2415  */
2416  typedef typename view_type::pointer pointer;
2417 
2418  /** const pointer type
2419  */
2421 
2422  /** reference type (result of operator[])
2423  */
2425 
2426  /** const reference type (result of operator[] const)
2427  */
2429 
2430  /** size type
2431  */
2433 
2434  /** difference type (used for multi-dimensional offsets and indices)
2435  */
2437 
2438  /** difference and index type for a single dimension
2439  */
2441 
2442  /** traverser type
2443  */
2444  typedef typename vigra::detail::MultiIteratorChooser <
2445  UnstridedArrayTag>::template Traverser <N, T, T &, T *>::type
2447 
2448  /** traverser type to const data
2449  */
2450  typedef typename vigra::detail::MultiIteratorChooser <
2451  UnstridedArrayTag>::template Traverser <N, T, T const &, T const *>::type
2453 
2454  /** sequential (random access) iterator type
2455  */
2456  typedef T * iterator;
2457 
2458  /** sequential (random access) const iterator type
2459  */
2460  typedef T * const_iterator;
2461 
2462 protected:
2463 
2464  typedef typename difference_type::value_type diff_zero_t;
2465 
2466  /** the allocator used to allocate the memory
2467  */
2469 
2470  /** allocate memory for s pixels, write its address into the given
2471  pointer and initialize the pixels with init.
2472  */
2473  void allocate (pointer &ptr, difference_type_1 s, const_reference init);
2474 
2475  /** allocate memory for s pixels, write its address into the given
2476  pointer and initialize the linearized pixels to the values of init.
2477  */
2478  template <class U>
2479  void allocate (pointer &ptr, difference_type_1 s, U const * init);
2480 
2481  /** allocate memory, write its address into the given
2482  pointer and initialize it by copying the data from the given MultiArrayView.
2483  */
2484  template <class U, class StrideTag>
2485  void allocate (pointer &ptr, MultiArrayView<N, U, StrideTag> const & init);
2486 
2487  /** deallocate the memory (of length s) starting at the given address.
2488  */
2489  void deallocate (pointer &ptr, difference_type_1 s);
2490 
2491  template <class U, class StrideTag>
2492  void copyOrReshape (const MultiArrayView<N, U, StrideTag> &rhs);
2493 public:
2494  /** default constructor
2495  */
2497  : MultiArrayView <N, T> (difference_type (diff_zero_t(0)),
2498  difference_type (diff_zero_t(0)), 0)
2499  {}
2500 
2501  /** construct with given allocator
2502  */
2503  MultiArray (allocator_type const & alloc)
2504  : MultiArrayView <N, T> (difference_type (diff_zero_t(0)),
2505  difference_type (diff_zero_t(0)), 0),
2506  m_alloc(alloc)
2507  {}
2508 
2509  /** construct with given length
2510 
2511  Use only for 1-dimensional arrays (<tt>N==1</tt>).
2512  */
2513  explicit MultiArray (difference_type_1 length,
2514  allocator_type const & alloc = allocator_type());
2515 
2516  /** construct with given shape
2517  */
2518  explicit MultiArray (const difference_type &shape,
2519  allocator_type const & alloc = allocator_type());
2520 
2521  /** construct from shape with an initial value
2522  */
2523  MultiArray (const difference_type &shape, const_reference init,
2524  allocator_type const & alloc = allocator_type());
2525 
2526  /** construct from shape and copy values from the given array
2527  */
2528  MultiArray (const difference_type &shape, const_pointer init,
2529  allocator_type const & alloc = allocator_type());
2530 
2531  /** copy constructor
2532  */
2533  MultiArray (const MultiArray &rhs)
2534  : MultiArrayView <N, T> (rhs.m_shape, rhs.m_stride, 0),
2535  m_alloc (rhs.m_alloc)
2536  {
2537  allocate (this->m_ptr, this->elementCount (), rhs.data ());
2538  }
2539 
2540  /** constructor from an array expression
2541  */
2542  template<class Expression>
2543  MultiArray (multi_math::MultiMathOperand<Expression> const & rhs,
2544  allocator_type const & alloc = allocator_type())
2545  : MultiArrayView <N, T> (difference_type (diff_zero_t(0)),
2546  difference_type (diff_zero_t(0)), 0),
2547  m_alloc (alloc)
2548  {
2549  multi_math::detail::assignOrResize(*this, rhs);
2550  }
2551 
2552  /** construct by copying from a MultiArrayView
2553  */
2554  template <class U, class StrideTag>
2556  allocator_type const & alloc = allocator_type());
2557 
2558  /** assignment.<br>
2559  If the size of \a rhs is the same as the left-hand side arrays's old size, only
2560  the data are copied. Otherwise, new storage is allocated, which invalidates all
2561  objects (array views, iterators) depending on the lhs array.
2562  */
2564  {
2565  if (this != &rhs)
2566  this->copyOrReshape(rhs);
2567  return *this;
2568  }
2569 
2570  /** assignment from arbitrary MultiArrayView.<br>
2571  If the size of \a rhs is the same as the left-hand side arrays's old size, only
2572  the data are copied. Otherwise, new storage is allocated, which invalidates all
2573  objects (array views, iterators) depending on the lhs array.
2574  */
2575  template <class U, class StrideTag>
2577  {
2578  this->copyOrReshape(rhs);
2579  return *this;
2580  }
2581 
2582  /** assignment from scalar.<br>
2583  Equivalent to MultiArray::init(v).
2584  */
2586  {
2587  return this->init(v);
2588  }
2589 
2590  /** Add-assignment from arbitrary MultiArrayView. Fails with
2591  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2592  If the left array has no data (hasData() is false), this function is
2593  equivalent to a normal assignment (i.e. an empty
2594  array is interpreted as a zero-array of appropriate size).
2595  */
2596  template <class U, class StrideTag>
2598  {
2599  if(this->hasData())
2600  view_type::operator+=(rhs);
2601  else
2602  *this = rhs;
2603  return *this;
2604  }
2605 
2606  /** Subtract-assignment from arbitrary MultiArrayView. Fails with
2607  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2608  If the left array has no data (hasData() is false), this function is
2609  equivalent to an assignment of the negated rhs (i.e. an empty
2610  array is interpreted as a zero-array of appropriate size).
2611  */
2612  template <class U, class StrideTag>
2614  {
2615  if(!this->hasData())
2616  this->reshape(rhs.shape());
2617  view_type::operator-=(rhs);
2618  return *this;
2619  }
2620 
2621  /** Multiply-assignment from arbitrary MultiArrayView. Fails with
2622  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2623  If the left array has no data (hasData() is false), this function is
2624  equivalent to reshape(rhs.shape()) with zero initialisation (i.e. an empty
2625  array is interpreted as a zero-array of appropriate size).
2626  */
2627  template <class U, class StrideTag>
2629  {
2630  if(this->hasData())
2631  view_type::operator*=(rhs);
2632  else
2633  this->reshape(rhs.shape());
2634  return *this;
2635  }
2636 
2637  /** Divide-assignment from arbitrary MultiArrayView. Fails with
2638  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2639  If the left array has no data (hasData() is false), this function is
2640  equivalent to reshape(rhs.shape()) with zero initialisation (i.e. an empty
2641  array is interpreted as a zero-array of appropriate size).
2642  */
2643  template <class U, class StrideTag>
2645  {
2646  if(this->hasData())
2647  view_type::operator/=(rhs);
2648  else
2649  this->reshape(rhs.shape());
2650  return *this;
2651  }
2652 
2653  /** Add-assignment of a scalar.
2654  */
2655  MultiArray &operator+= (const T &rhs)
2656  {
2657  view_type::operator+=(rhs);
2658  return *this;
2659  }
2660 
2661  /** Subtract-assignment of a scalar.
2662  */
2663  MultiArray &operator-= (const T &rhs)
2664  {
2665  view_type::operator-=(rhs);
2666  return *this;
2667  }
2668 
2669  /** Multiply-assignment of a scalar.
2670  */
2671  MultiArray &operator*= (const T &rhs)
2672  {
2673  view_type::operator*=(rhs);
2674  return *this;
2675  }
2676 
2677  /** Divide-assignment of a scalar.
2678  */
2679  MultiArray &operator/= (const T &rhs)
2680  {
2681  view_type::operator/=(rhs);
2682  return *this;
2683  }
2684  /** Assignment of an array expression. Fails with
2685  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2686  */
2687  template<class Expression>
2688  MultiArray & operator=(multi_math::MultiMathOperand<Expression> const & rhs)
2689  {
2690  multi_math::detail::assignOrResize(*this, rhs);
2691  return *this;
2692  }
2693 
2694  /** Add-assignment of an array expression. Fails with
2695  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2696  */
2697  template<class Expression>
2698  MultiArray & operator+=(multi_math::MultiMathOperand<Expression> const & rhs)
2699  {
2700  multi_math::detail::plusAssignOrResize(*this, rhs);
2701  return *this;
2702  }
2703 
2704  /** Subtract-assignment of an array expression. Fails with
2705  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2706  */
2707  template<class Expression>
2708  MultiArray & operator-=(multi_math::MultiMathOperand<Expression> const & rhs)
2709  {
2710  multi_math::detail::minusAssignOrResize(*this, rhs);
2711  return *this;
2712  }
2713 
2714  /** Multiply-assignment of an array expression. Fails with
2715  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2716  */
2717  template<class Expression>
2718  MultiArray & operator*=(multi_math::MultiMathOperand<Expression> const & rhs)
2719  {
2720  multi_math::detail::multiplyAssignOrResize(*this, rhs);
2721  return *this;
2722  }
2723 
2724  /** Divide-assignment of an array expression. Fails with
2725  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2726  */
2727  template<class Expression>
2728  MultiArray & operator/=(multi_math::MultiMathOperand<Expression> const & rhs)
2729  {
2730  multi_math::detail::divideAssignOrResize(*this, rhs);
2731  return *this;
2732  }
2733 
2734  /** destructor
2735  */
2737  {
2738  deallocate (this->m_ptr, this->elementCount ());
2739  }
2740 
2741 
2742  /** init elements with a constant
2743  */
2744  template <class U>
2745  MultiArray & init(const U & init)
2746  {
2747  view_type::init(init);
2748  return *this;
2749  }
2750 
2751  /** Allocate new memory with the given shape and initialize with zeros.<br>
2752  <em>Note:</em> this operation invalidates all dependent objects
2753  (array views and iterators)
2754  */
2755  void reshape (const difference_type &shape)
2756  {
2757  reshape (shape, T());
2758  }
2759 
2760  /** Allocate new memory with the given shape and initialize it
2761  with the given value.<br>
2762  <em>Note:</em> this operation invalidates all dependent objects
2763  (array views and iterators)
2764  */
2765  void reshape (const difference_type &shape, const_reference init);
2766 
2767  /** Swap the contents with another MultiArray. This is fast,
2768  because no data are copied, but only pointers and shapes swapped.
2769  <em>Note:</em> this operation invalidates all dependent objects
2770  (array views and iterators)
2771  */
2772  void swap (MultiArray & other);
2773 
2774  /** sequential iterator pointing to the first array element.
2775  */
2777  {
2778  return this->data();
2779  }
2780 
2781  /** sequential iterator pointing beyond the last array element.
2782  */
2784  {
2785  return this->data() + this->elementCount();
2786  }
2787 
2788  /** sequential const iterator pointing to the first array element.
2789  */
2791  {
2792  return this->data();
2793  }
2794 
2795  /** sequential const iterator pointing beyond the last array element.
2796  */
2798  {
2799  return this->data() + this->elementCount();
2800  }
2801 
2802  /** get the allocator.
2803  */
2804  allocator_type const & allocator () const
2805  {
2806  return m_alloc;
2807  }
2808 };
2809 
2810 template <unsigned int N, class T, class A>
2812  allocator_type const & alloc)
2814  detail::defaultStride <1> (difference_type(length)),
2815  0),
2816  m_alloc(alloc)
2817 {
2818  allocate (this->m_ptr, this->elementCount (), T());
2819 }
2820 
2821 template <unsigned int N, class T, class A>
2823  allocator_type const & alloc)
2824 : MultiArrayView <N, T> (shape,
2825  detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
2826  0),
2827  m_alloc(alloc)
2828 {
2829  if (N == 0)
2830  {
2831  this->m_shape [0] = 1;
2832  this->m_stride [0] = 0;
2833  }
2834  allocate (this->m_ptr, this->elementCount (), T());
2835 }
2836 
2837 template <unsigned int N, class T, class A>
2839  allocator_type const & alloc)
2840 : MultiArrayView <N, T> (shape,
2841  detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
2842  0),
2843  m_alloc(alloc)
2844 {
2845  if (N == 0)
2846  {
2847  this->m_shape [0] = 1;
2848  this->m_stride [0] = 0;
2849  }
2850  allocate (this->m_ptr, this->elementCount (), init);
2851 }
2852 
2853 template <unsigned int N, class T, class A>
2855  allocator_type const & alloc)
2856 : MultiArrayView <N, T> (shape,
2857  detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
2858  0),
2859  m_alloc(alloc)
2860 {
2861  if (N == 0)
2862  {
2863  this->m_shape [0] = 1;
2864  this->m_stride [0] = 0;
2865  }
2866  allocate (this->m_ptr, this->elementCount (), init);
2867 }
2868 
2869 template <unsigned int N, class T, class A>
2870 template <class U, class StrideTag>
2872  allocator_type const & alloc)
2873 : MultiArrayView <N, T> (rhs.shape(),
2874  detail::defaultStride <MultiArrayView<N,T>::actual_dimension>(rhs.shape()),
2875  0),
2876  m_alloc (alloc)
2877 {
2878  allocate (this->m_ptr, rhs);
2879 }
2880 
2881 template <unsigned int N, class T, class A>
2882 template <class U, class StrideTag>
2883 void
2885 {
2886  if (this->shape() == rhs.shape())
2887  this->copy(rhs);
2888  else
2889  {
2890  MultiArray t(rhs);
2891  this->swap(t);
2892  }
2893 }
2894 
2895 template <unsigned int N, class T, class A>
2897  const_reference initial)
2898 {
2899  if (N== 0)
2900  {
2901  return;
2902  }
2903  else if(new_shape == this->shape())
2904  {
2905  this->init(initial);
2906  }
2907  else
2908  {
2909  difference_type new_stride = detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (new_shape);
2911  T *new_ptr;
2912  allocate (new_ptr, new_size, initial);
2913  deallocate (this->m_ptr, this->elementCount ());
2914  this->m_ptr = new_ptr;
2915  this->m_shape = new_shape;
2916  this->m_stride = new_stride;
2917  }
2918 }
2919 
2920 
2921 template <unsigned int N, class T, class A>
2922 inline void
2924 {
2925  if (this == &other)
2926  return;
2927  std::swap(this->m_shape, other.m_shape);
2928  std::swap(this->m_stride, other.m_stride);
2929  std::swap(this->m_ptr, other.m_ptr);
2930  std::swap(this->m_alloc, other.m_alloc);
2931 }
2932 
2933 template <unsigned int N, class T, class A>
2936 {
2937  ptr = m_alloc.allocate ((typename A::size_type)s);
2939  try {
2940  for (i = 0; i < s; ++i)
2941  m_alloc.construct (ptr + i, init);
2942  }
2943  catch (...) {
2944  for (difference_type_1 j = 0; j < i; ++j)
2945  m_alloc.destroy (ptr + j);
2946  m_alloc.deallocate (ptr, (typename A::size_type)s);
2947  throw;
2948  }
2949 }
2950 
2951 template <unsigned int N, class T, class A>
2952 template <class U>
2954  U const * init)
2955 {
2956  ptr = m_alloc.allocate ((typename A::size_type)s);
2958  try {
2959  for (i = 0; i < s; ++i, ++init)
2960  m_alloc.construct (ptr + i, *init);
2961  }
2962  catch (...) {
2963  for (difference_type_1 j = 0; j < i; ++j)
2964  m_alloc.destroy (ptr + j);
2965  m_alloc.deallocate (ptr, (typename A::size_type)s);
2966  throw;
2967  }
2968 }
2969 
2970 template <unsigned int N, class T, class A>
2971 template <class U, class StrideTag>
2973 {
2974  difference_type_1 s = init.elementCount();
2975  ptr = m_alloc.allocate ((typename A::size_type)s);
2976  pointer p = ptr;
2977  try {
2978  detail::uninitializedCopyMultiArrayData(init.traverser_begin(), init.shape(),
2979  p, m_alloc, MetaInt<actual_dimension-1>());
2980  }
2981  catch (...) {
2982  for (pointer pp = ptr; pp < p; ++pp)
2983  m_alloc.destroy (pp);
2984  m_alloc.deallocate (ptr, (typename A::size_type)s);
2985  throw;
2986  }
2987 }
2988 
2989 template <unsigned int N, class T, class A>
2991 {
2992  if (ptr == 0)
2993  return;
2994  for (difference_type_1 i = 0; i < s; ++i)
2995  m_alloc.destroy (ptr + i);
2996  m_alloc.deallocate (ptr, (typename A::size_type)s);
2997  ptr = 0;
2998 }
2999 
3000 /********************************************************/
3001 /* */
3002 /* argument object factories */
3003 /* */
3004 /********************************************************/
3005 
3006 template <unsigned int N, class T, class StrideTag>
3007 inline triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3010 srcMultiArrayRange( MultiArrayView<N,T,StrideTag> const & array )
3011 {
3012  return triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3015  ( array.traverser_begin(),
3016  array.shape(),
3018 }
3019 
3020 template <unsigned int N, class T, class StrideTag, class Accessor>
3021 inline triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3022  typename MultiArrayView<N,T,StrideTag>::difference_type,
3023  Accessor >
3024 srcMultiArrayRange( MultiArrayView<N,T,StrideTag> const & array, Accessor a )
3025 {
3026  return triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3027  typename MultiArrayView<N,T,StrideTag>::difference_type,
3028  Accessor >
3029  ( array.traverser_begin(),
3030  array.shape(),
3031  a);
3032 }
3033 
3034 template <unsigned int N, class T, class StrideTag>
3035 inline pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3036  typename AccessorTraits<T>::default_const_accessor >
3037 srcMultiArray( MultiArrayView<N,T,StrideTag> const & array )
3038 {
3039  return pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3040  typename AccessorTraits<T>::default_const_accessor >
3041  ( array.traverser_begin(),
3042  typename AccessorTraits<T>::default_const_accessor() );
3043 }
3044 
3045 template <unsigned int N, class T, class StrideTag, class Accessor>
3046 inline pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3047  Accessor >
3048 srcMultiArray( MultiArrayView<N,T,StrideTag> const & array, Accessor a )
3049 {
3050  return pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3051  Accessor >
3052  ( array.traverser_begin(), a );
3053 }
3054 
3055 template <unsigned int N, class T, class StrideTag>
3056 inline triple<typename MultiArrayView<N,T,StrideTag>::traverser,
3057  typename MultiArrayView<N,T,StrideTag>::difference_type,
3058  typename AccessorTraits<T>::default_accessor >
3059 destMultiArrayRange( MultiArrayView<N,T,StrideTag> & array )
3060 {
3061  return triple<typename MultiArrayView<N,T,StrideTag>::traverser,
3062  typename MultiArrayView<N,T,StrideTag>::difference_type,
3063  typename AccessorTraits<T>::default_accessor >
3064  ( array.traverser_begin(),
3065  array.shape(),
3066  typename AccessorTraits<T>::default_accessor() );
3067 }
3068 
3069 template <unsigned int N, class T, class StrideTag, class Accessor>
3070 inline triple<typename MultiArrayView<N,T,StrideTag>::traverser,
3071  typename MultiArrayView<N,T,StrideTag>::difference_type,
3072  Accessor >
3073 destMultiArrayRange( MultiArrayView<N,T,StrideTag> & array, Accessor a )
3074 {
3075  return triple<typename MultiArrayView<N,T,StrideTag>::traverser,
3076  typename MultiArrayView<N,T,StrideTag>::difference_type,
3077  Accessor >
3078  ( array.traverser_begin(),
3079  array.shape(),
3080  a );
3081 }
3082 
3083 template <unsigned int N, class T, class StrideTag>
3084 inline pair<typename MultiArrayView<N,T,StrideTag>::traverser,
3085  typename AccessorTraits<T>::default_accessor >
3086 destMultiArray( MultiArrayView<N,T,StrideTag> & array )
3087 {
3088  return pair<typename MultiArrayView<N,T,StrideTag>::traverser,
3089  typename AccessorTraits<T>::default_accessor >
3090  ( array.traverser_begin(),
3091  typename AccessorTraits<T>::default_accessor() );
3092 }
3093 
3094 template <unsigned int N, class T, class StrideTag, class Accessor>
3095 inline pair<typename MultiArrayView<N,T,StrideTag>::traverser,
3096  Accessor >
3097 destMultiArray( MultiArrayView<N,T,StrideTag> & array, Accessor a )
3098 {
3099  return pair<typename MultiArrayView<N,T,StrideTag>::traverser,
3100  Accessor >
3101  ( array.traverser_begin(), a );
3102 }
3103 
3104 /********************************************************************/
3105 
3106 template <class PixelType, class Accessor>
3107 inline triple<ConstStridedImageIterator<PixelType>,
3108  ConstStridedImageIterator<PixelType>, Accessor>
3109 srcImageRange(const MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
3110 {
3111  ConstStridedImageIterator<PixelType>
3112  ul(img.data(), 1, img.stride(0), img.stride(1));
3113  return triple<ConstStridedImageIterator<PixelType>,
3114  ConstStridedImageIterator<PixelType>,
3115  Accessor>(
3116  ul, ul + Size2D(img.shape(0), img.shape(1)), a);
3117 }
3118 
3119 template <class PixelType, class Accessor>
3120 inline pair<ConstStridedImageIterator<PixelType>, Accessor>
3121 srcImage(const MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
3122 {
3123  ConstStridedImageIterator<PixelType>
3124  ul(img.data(), 1, img.stride(0), img.stride(1));
3125  return pair<ConstStridedImageIterator<PixelType>, Accessor>
3126  (ul, a);
3127 }
3128 
3129 template <class PixelType, class Accessor>
3130 inline triple<StridedImageIterator<PixelType>,
3131  StridedImageIterator<PixelType>, Accessor>
3132 destImageRange(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
3133 {
3134  StridedImageIterator<PixelType>
3135  ul(img.data(), 1, img.stride(0), img.stride(1));
3136  return triple<StridedImageIterator<PixelType>,
3137  StridedImageIterator<PixelType>,
3138  Accessor>(
3139  ul, ul + Size2D(img.shape(0), img.shape(1)), a);
3140 }
3141 
3142 template <class PixelType, class Accessor>
3143 inline pair<StridedImageIterator<PixelType>, Accessor>
3144 destImage(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
3145 {
3146  StridedImageIterator<PixelType>
3147  ul(img.data(), 1, img.stride(0), img.stride(1));
3148  return pair<StridedImageIterator<PixelType>, Accessor>
3149  (ul, a);
3150 }
3151 
3152 template <class PixelType, class Accessor>
3153 inline pair<StridedImageIterator<PixelType>, Accessor>
3154 maskImage(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
3155 {
3156  StridedImageIterator<PixelType>
3157  ul(img.data(), 1, img.stride(0), img.stride(1));
3158  return pair<StridedImageIterator<PixelType>, Accessor>
3159  (ul, a);
3160 }
3161 
3162 // -------------------------------------------------------------------
3163 
3164 template <class PixelType>
3165 inline triple<ConstStridedImageIterator<PixelType>,
3166  ConstStridedImageIterator<PixelType>,
3167  typename AccessorTraits<PixelType>::default_const_accessor>
3168 srcImageRange(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
3169 {
3170  ConstStridedImageIterator<PixelType>
3171  ul(img.data(), 1, img.stride(0), img.stride(1));
3172  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
3173  return triple<ConstStridedImageIterator<PixelType>,
3174  ConstStridedImageIterator<PixelType>,
3175  Accessor>
3176  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
3177 }
3178 
3179 template <class PixelType>
3180 inline triple<ConstImageIterator<PixelType>,
3181  ConstImageIterator<PixelType>,
3182  typename AccessorTraits<PixelType>::default_const_accessor>
3183 srcImageRange(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
3184 {
3185  ConstImageIterator<PixelType>
3186  ul(img.data(), img.stride(1));
3187  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
3188  return triple<ConstImageIterator<PixelType>,
3189  ConstImageIterator<PixelType>,
3190  Accessor>
3191  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
3192 }
3193 
3194 template <class PixelType>
3195 inline pair< ConstStridedImageIterator<PixelType>,
3196  typename AccessorTraits<PixelType>::default_const_accessor>
3197 srcImage(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
3198 {
3199  ConstStridedImageIterator<PixelType>
3200  ul(img.data(), 1, img.stride(0), img.stride(1));
3201  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
3202  return pair<ConstStridedImageIterator<PixelType>,
3203  Accessor>
3204  (ul, Accessor());
3205 }
3206 
3207 template <class PixelType>
3208 inline pair< ConstImageIterator<PixelType>,
3209  typename AccessorTraits<PixelType>::default_const_accessor>
3210 srcImage(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
3211 {
3212  ConstImageIterator<PixelType>
3213  ul(img.data(), img.stride(1));
3214  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
3215  return pair<ConstImageIterator<PixelType>,
3216  Accessor>
3217  (ul, Accessor());
3218 }
3219 
3220 template <class PixelType>
3221 inline triple< StridedImageIterator<PixelType>,
3222  StridedImageIterator<PixelType>,
3223  typename AccessorTraits<PixelType>::default_accessor>
3224 destImageRange(MultiArrayView<2, PixelType, StridedArrayTag> & img)
3225 {
3226  StridedImageIterator<PixelType>
3227  ul(img.data(), 1, img.stride(0), img.stride(1));
3228  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3229  return triple<StridedImageIterator<PixelType>,
3230  StridedImageIterator<PixelType>,
3231  Accessor>
3232  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
3233 }
3234 
3235 template <class PixelType>
3236 inline triple< ImageIterator<PixelType>,
3237  ImageIterator<PixelType>,
3238  typename AccessorTraits<PixelType>::default_accessor>
3239 destImageRange(MultiArrayView<2, PixelType, UnstridedArrayTag> & img)
3240 {
3241  ImageIterator<PixelType>
3242  ul(img.data(), img.stride(1));
3243  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3244  return triple<ImageIterator<PixelType>,
3245  ImageIterator<PixelType>,
3246  Accessor>
3247  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
3248 }
3249 
3250 template <class PixelType>
3251 inline pair< StridedImageIterator<PixelType>,
3252  typename AccessorTraits<PixelType>::default_accessor>
3253 destImage(MultiArrayView<2, PixelType, StridedArrayTag> & img)
3254 {
3255  StridedImageIterator<PixelType>
3256  ul(img.data(), 1, img.stride(0), img.stride(1));
3257  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3258  return pair<StridedImageIterator<PixelType>, Accessor>
3259  (ul, Accessor());
3260 }
3261 
3262 template <class PixelType>
3263 inline pair< ImageIterator<PixelType>,
3264  typename AccessorTraits<PixelType>::default_accessor>
3265 destImage(MultiArrayView<2, PixelType, UnstridedArrayTag> & img)
3266 {
3267  ImageIterator<PixelType> ul(img.data(), img.stride(1));
3268  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3269  return pair<ImageIterator<PixelType>, Accessor>(ul, Accessor());
3270 }
3271 
3272 template <class PixelType>
3273 inline pair< ConstStridedImageIterator<PixelType>,
3274  typename AccessorTraits<PixelType>::default_accessor>
3275 maskImage(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
3276 {
3277  ConstStridedImageIterator<PixelType>
3278  ul(img.data(), 1, img.stride(0), img.stride(1));
3279  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3280  return pair<ConstStridedImageIterator<PixelType>, Accessor>
3281  (ul, Accessor());
3282 }
3283 
3284 template <class PixelType>
3285 inline pair< ConstImageIterator<PixelType>,
3286  typename AccessorTraits<PixelType>::default_accessor>
3287 maskImage(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
3288 {
3289  ConstImageIterator<PixelType>
3290  ul(img.data(), img.stride(1));
3291  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3292  return pair<ConstImageIterator<PixelType>, Accessor>
3293  (ul, Accessor());
3294 }
3295 
3296 /********************************************************/
3297 /* */
3298 /* makeBasicImageView */
3299 /* */
3300 /********************************************************/
3301 
3302 /** \addtogroup MultiArrayToImage Wrap a \ref vigra::MultiArrayView in
3303  a \ref vigra::BasicImageView
3304 */
3305 //@{
3306 /** Create a \ref vigra::BasicImageView from an unstrided 2-dimensional
3307  \ref vigra::MultiArrayView.
3308 
3309  The \ref vigra::BasicImageView will have the same <tt>value_type </tt>
3310  as the original \ref vigra::MultiArrayView.
3311 */
3312 template <class T>
3313 BasicImageView <T>
3315 {
3316  return BasicImageView <T> (array.data (), array.shape (0),
3317  array.shape (1));
3318 }
3319 
3320 /** Create a \ref vigra::BasicImageView from a 3-dimensional
3321  \ref vigra::MultiArray.
3322 
3323  This wrapper flattens the two innermost dimensions of the array
3324  into single rows of the resulting image.
3325  The \ref vigra::BasicImageView will have the same <tt>value_type </tt>
3326  as the original \ref vigra::MultiArray.
3327 */
3328 template <class T>
3329 BasicImageView <T>
3331 {
3332  return BasicImageView <T> (array.data (),
3333  array.shape (0)*array.shape (1), array.shape (2));
3334 }
3335 
3336 /** Create a \ref vigra::BasicImageView from a 3-dimensional
3337  \ref vigra::MultiArray.
3338 
3339  This wrapper only works if <tt>T</tt> is a scalar type and the
3340  array's innermost dimension has size 3. It then re-interprets
3341  the data array as a 2-dimensional array with value_type
3342  <tt>RGBValue<T></tt>.
3343 */
3344 template <class T>
3345 BasicImageView <RGBValue<T> >
3347 {
3348  vigra_precondition (
3349  array.shape (0) == 3, "makeRGBImageView(): array.shape(0) must be 3.");
3350  return BasicImageView <RGBValue<T> > (
3351  reinterpret_cast <RGBValue <T> *> (array.data ()),
3352  array.shape (1), array.shape (2));
3353 }
3354 
3355 //@}
3356 
3357 } // namespace vigra
3358 
3359 #undef VIGRA_ASSERT_INSIDE
3360 
3361 #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.9.0 (Tue Oct 22 2013)