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

multi_array.hxx VIGRA

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 "basicimage.hxx"
46 #include "imageiterator.hxx"
47 #include "numerictraits.hxx"
48 #include "multi_iterator.hxx"
49 #include "multi_pointoperators.hxx"
50 #include "metaprogramming.hxx"
51 #include "mathutil.hxx"
52 #include "algorithm.hxx"
53 
54 // Bounds checking Macro used if VIGRA_CHECK_BOUNDS is defined.
55 #ifdef VIGRA_CHECK_BOUNDS
56 #define VIGRA_ASSERT_INSIDE(diff) \
57  vigra_precondition(this->isInside(diff), "Index out of bounds")
58 #else
59 #define VIGRA_ASSERT_INSIDE(diff)
60 #endif
61 
62 namespace vigra
63 {
64 
65 namespace detail
66 {
67 
68 /********************************************************/
69 /* */
70 /* MaybeStrided */
71 /* */
72 /********************************************************/
73 
74 /* metatag implementing a test for marking MultiArrays that were
75  indexed at the zero'th dimension as strided, and all others as
76  unstrided.
77 
78  <b>\#include</b> <vigra/multi_array.hxx> <br/>
79  Namespace: vigra::detail
80 */
81 template <class StrideTag, unsigned int N>
82 struct MaybeStrided
83 {
84  typedef StrideTag type;
85 };
86 
87 template <class StrideTag>
88 struct MaybeStrided <StrideTag, 0>
89 {
90  typedef StridedArrayTag type;
91 };
92 
93 /********************************************************/
94 /* */
95 /* MultiIteratorChooser */
96 /* */
97 /********************************************************/
98 
99 /* metatag implementing a test (by pattern matching) for marking
100  MultiArrays that were indexed at the zero'th dimension as strided.
101 
102  <b>\#include</b> <vigra/multi_array.hxx> <br/>
103  Namespace: vigra::detail
104 */
105 template <class O>
106 struct MultiIteratorChooser;
107 
108 /********************************************************/
109 /* */
110 /* MultiIteratorChooser <StridedArrayTag> */
111 /* */
112 /********************************************************/
113 
114 /* specialization of the MultiIteratorChooser for strided arrays.
115 
116  <b>\#include</b> <vigra/multi_array.hxx> <br/>
117  Namespace: vigra::detail
118 */
119 template <>
120 struct MultiIteratorChooser <StridedArrayTag>
121 {
122  template <unsigned int N, class T, class REFERENCE, class POINTER>
123  struct Traverser
124  {
125  typedef StridedMultiIterator <N, T, REFERENCE, POINTER> type;
126  };
127 
128  template <unsigned int N, class T, class REFERENCE, class POINTER>
129  struct Iterator
130  {
131  typedef StridedScanOrderIterator <N, T, REFERENCE, POINTER> type;
132  };
133 
134  template <class Iter, class View>
135  static Iter constructIterator(View * v)
136  {
137  return v->begin();
138  }
139 };
140 
141 /********************************************************/
142 /* */
143 /* MultiIteratorChooser <UnstridedArrayTag> */
144 /* */
145 /********************************************************/
146 
147 /* specialization of the MultiIteratorChooser for unstrided arrays.
148 
149  <b>\#include</b> <vigra/multi_array.hxx> <br/>
150  Namespace: vigra::detail
151 */
152 template <>
153 struct MultiIteratorChooser <UnstridedArrayTag>
154 {
155  template <unsigned int N, class T, class REFERENCE, class POINTER>
156  struct Traverser
157  {
158  typedef MultiIterator <N, T, REFERENCE, POINTER> type;
159  };
160 
161  template <unsigned int N, class T, class REFERENCE, class POINTER>
162  struct Iterator
163  {
164  typedef POINTER type;
165  };
166 
167  template <class Iter, class View>
168  static Iter constructIterator(View * v)
169  {
170  return v->data();
171  }
172 };
173 
174 /********************************************************/
175 /* */
176 /* helper functions */
177 /* */
178 /********************************************************/
179 
180 template <class DestIterator, class Shape, class T>
181 inline void
182 initMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<0>)
183 {
184  DestIterator dend = d + shape[0];
185  for(; d < dend; ++d)
186  {
187  *d = init;
188  }
189 }
190 
191 template <class DestIterator, class Shape, class T, int N>
192 void
193 initMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<N>)
194 {
195  DestIterator dend = d + shape[N];
196  for(; d < dend; ++d)
197  {
198  initMultiArrayData(d.begin(), shape, init, MetaInt<N-1>());
199  }
200 }
201 
202 // FIXME: the explicit overload for MultiIterator<1, UInt8, ... > works around a compiler crash in VisualStudio 2010
203 #define VIGRA_COPY_MULTI_ARRAY_DATA(name, op) \
204 template <class SrcIterator, class Shape, class DestIterator> \
205 inline void \
206 name##MultiArrayData(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>) \
207 { \
208  for(MultiArrayIndex i=0; i < shape[0]; ++i, ++s, ++d) \
209  { \
210  *d op detail::RequiresExplicitCast<typename DestIterator::value_type>::cast(*s); \
211  } \
212 } \
213  \
214 template <class Ref, class Ptr, class Shape, class DestIterator> \
215 inline void \
216 name##MultiArrayData(MultiIterator<1, UInt8, Ref, Ptr> si, Shape const & shape, DestIterator d, MetaInt<0>) \
217 { \
218  Ptr s = &(*si); \
219  for(MultiArrayIndex i=0; i < shape[0]; ++i, ++s, ++d) \
220  { \
221  *d op detail::RequiresExplicitCast<typename DestIterator::value_type>::cast(*s); \
222  } \
223 } \
224 \
225 template <class SrcIterator, class Shape, class DestIterator, int N> \
226 void \
227 name##MultiArrayData(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>) \
228 { \
229  for(MultiArrayIndex i=0; i < shape[N]; ++i, ++s, ++d) \
230  { \
231  name##MultiArrayData(s.begin(), shape, d.begin(), MetaInt<N-1>()); \
232  } \
233 } \
234 \
235 template <class DestIterator, class Shape, class T> \
236 inline void \
237 name##ScalarMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<0>) \
238 { \
239  for(MultiArrayIndex i=0; i < shape[0]; ++i, ++d) \
240  { \
241  *d op detail::RequiresExplicitCast<typename DestIterator::value_type>::cast(init); \
242  } \
243 } \
244  \
245 template <class DestIterator, class Shape, class T, int N> \
246 void \
247 name##ScalarMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<N>) \
248 { \
249  for(MultiArrayIndex i=0; i < shape[N]; ++i, ++d) \
250  { \
251  name##ScalarMultiArrayData(d.begin(), shape, init, MetaInt<N-1>()); \
252  } \
253 }
254 
255 VIGRA_COPY_MULTI_ARRAY_DATA(copy, =)
256 VIGRA_COPY_MULTI_ARRAY_DATA(copyAdd, +=)
257 VIGRA_COPY_MULTI_ARRAY_DATA(copySub, -=)
258 VIGRA_COPY_MULTI_ARRAY_DATA(copyMul, *=)
259 VIGRA_COPY_MULTI_ARRAY_DATA(copyDiv, /=)
260 
261 #undef VIGRA_COPY_MULTI_ARRAY_DATA
262 
263 template <class SrcIterator, class Shape, class T, class ALLOC>
264 inline void
265 uninitializedCopyMultiArrayData(SrcIterator s, Shape const & shape, T * & d, ALLOC & a, MetaInt<0>)
266 {
267  SrcIterator send = s + shape[0];
268  for(; s < send; ++s, ++d)
269  {
270  a.construct(d, static_cast<T const &>(*s));
271  }
272 }
273 
274 // FIXME: this overload works around a compiler crash in VisualStudio 2010
275 template <class Ref, class Ptr, class Shape, class T, class ALLOC>
276 inline void
277 uninitializedCopyMultiArrayData(MultiIterator<1, UInt8, Ref, Ptr> si, Shape const & shape, T * & d, ALLOC & a, MetaInt<0>)
278 {
279  Ptr s = &(*si), send = s + shape[0];
280  for(; s < send; ++s, ++d)
281  {
282  a.construct(d, static_cast<T const &>(*s));
283  }
284 }
285 
286 template <class SrcIterator, class Shape, class T, class ALLOC, int N>
287 void
288 uninitializedCopyMultiArrayData(SrcIterator s, Shape const & shape, T * & d, ALLOC & a, MetaInt<N>)
289 {
290  SrcIterator send = s + shape[N];
291  for(; s < send; ++s)
292  {
293  uninitializedCopyMultiArrayData(s.begin(), shape, d, a, MetaInt<N-1>());
294  }
295 }
296 
297 template <class SrcIterator, class Shape, class T, class Functor>
298 inline void
299 reduceOverMultiArray(SrcIterator s, Shape const & shape, T & result, Functor const & f, MetaInt<0>)
300 {
301  SrcIterator send = s + shape[0];
302  for(; s < send; ++s)
303  {
304  f(result, *s);
305  }
306 }
307 
308 template <class SrcIterator, class Shape, class T, class Functor, int N>
309 void
310 reduceOverMultiArray(SrcIterator s, Shape const & shape, T & result, Functor const & f, MetaInt<N>)
311 {
312  SrcIterator send = s + shape[N];
313  for(; s < send; ++s)
314  {
315  reduceOverMultiArray(s.begin(), shape, result, f, MetaInt<N-1>());
316  }
317 }
318 
319 struct MaxNormReduceFunctor
320 {
321  template <class T, class U>
322  void operator()(T & result, U const & u) const
323  {
324  T v = norm(u);
325  if(result < v)
326  result = v;
327  }
328 };
329 
330 struct L1NormReduceFunctor
331 {
332  template <class T, class U>
333  void operator()(T & result, U const & u) const
334  {
335  result += norm(u);
336  }
337 };
338 
339 struct SquaredL2NormReduceFunctor
340 {
341  template <class T, class U>
342  void operator()(T & result, U const & u) const
343  {
344  result += squaredNorm(u);
345  }
346 };
347 
348 template <class T>
349 struct WeightedL2NormReduceFunctor
350 {
351  T scale;
352 
353  WeightedL2NormReduceFunctor(T s)
354  : scale(s)
355  {}
356 
357  template <class U>
358  void operator()(T & result, U const & u) const
359  {
360  result += squaredNorm(u * scale);
361  }
362 };
363 
364 struct SumReduceFunctor
365 {
366  template <class T, class U>
367  void operator()(T & result, U const & u) const
368  {
369  result += u;
370  }
371 };
372 
373 struct ProdReduceFunctor
374 {
375  template <class T, class U>
376  void operator()(T & result, U const & u) const
377  {
378  result *= u;
379  }
380 };
381 
382 struct MinmaxReduceFunctor
383 {
384  template <class T, class U>
385  void operator()(T & result, U const & u) const
386  {
387  if(u < result.first)
388  result.first = u;
389  if(result.second < u)
390  result.second = u;
391  }
392 };
393 
394 struct MeanVarianceReduceFunctor
395 {
396  template <class T, class U>
397  void operator()(T & result, U const & u) const
398  {
399  ++result.first;
400  typename T::second_type t1 = u - result.second;
401  typename T::second_type t2 = t1 / result.first;
402  result.second += t2;
403  result.third += (result.first-1.0)*t1*t2;
404  }
405 };
406 
407 struct AllTrueReduceFunctor
408 {
409  template <class T, class U>
410  void operator()(T & result, U const & u) const
411  {
412  result = result && (u != NumericTraits<U>::zero());
413  }
414 };
415 
416 struct AnyTrueReduceFunctor
417 {
418  template <class T, class U>
419  void operator()(T & result, U const & u) const
420  {
421  result = result || (u != NumericTraits<U>::zero());
422  }
423 };
424 
425 template <class SrcIterator, class Shape, class DestIterator>
426 inline bool
427 equalityOfMultiArrays(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>)
428 {
429  SrcIterator send = s + shape[0];
430  for(; s < send; ++s, ++d)
431  {
432  if(!(*s == *d))
433  return false;
434  }
435  return true;
436 }
437 
438 template <class SrcIterator, class Shape, class DestIterator, int N>
439 bool
440 equalityOfMultiArrays(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>)
441 {
442  SrcIterator send = s + shape[N];
443  for(; s < send; ++s, ++d)
444  {
445  if(!equalityOfMultiArrays(s.begin(), shape, d.begin(), MetaInt<N-1>()))
446  return false;
447  }
448  return true;
449 }
450 
451 
452 template <class SrcIterator, class Shape, class DestIterator>
453 inline void
454 swapDataImpl(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>)
455 {
456  SrcIterator send = s + shape[0];
457  for(; s < send; ++s, ++d)
458  std::swap(*s, *d);
459 }
460 
461 template <class SrcIterator, class Shape, class DestIterator, int N>
462 void
463 swapDataImpl(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>)
464 {
465  SrcIterator send = s + shape[N];
466  for(; s < send; ++s, ++d)
467  swapDataImpl(s.begin(), shape, d.begin(), MetaInt<N-1>());
468 }
469 
470 } // namespace detail
471 
472 /********************************************************/
473 /* */
474 /* MultiArrayView */
475 /* */
476 /********************************************************/
477 
478 // forward declarations
479 
480 namespace multi_math {
481 
482 template <class T>
483 struct MultiMathOperand;
484 
485 namespace math_detail {
486 
487 template <unsigned int N, class T, class C, class E>
488 void assign(MultiArrayView<N, T, C>, MultiMathOperand<E> const &);
489 
490 template <unsigned int N, class T, class C, class E>
491 void plusAssign(MultiArrayView<N, T, C>, MultiMathOperand<E> const &);
492 
493 template <unsigned int N, class T, class C, class E>
494 void minusAssign(MultiArrayView<N, T, C>, MultiMathOperand<E> const &);
495 
496 template <unsigned int N, class T, class C, class E>
497 void multiplyAssign(MultiArrayView<N, T, C>, MultiMathOperand<E> const &);
498 
499 template <unsigned int N, class T, class C, class E>
500 void divideAssign(MultiArrayView<N, T, C>, MultiMathOperand<E> const &);
501 
502 template <unsigned int N, class T, class A, class E>
503 void assignOrResize(MultiArray<N, T, A> &, MultiMathOperand<E> const &);
504 
505 template <unsigned int N, class T, class A, class E>
506 void plusAssignOrResize(MultiArray<N, T, A> &, MultiMathOperand<E> const &);
507 
508 template <unsigned int N, class T, class A, class E>
509 void minusAssignOrResize(MultiArray<N, T, A> &, MultiMathOperand<E> const &);
510 
511 template <unsigned int N, class T, class A, class E>
512 void multiplyAssignOrResize(MultiArray<N, T, A> &, MultiMathOperand<E> const &);
513 
514 template <unsigned int N, class T, class A, class E>
515 void divideAssignOrResize(MultiArray<N, T, A> &, MultiMathOperand<E> const &);
516 
517 } // namespace math_detail
518 
519 } // namespace multi_math
520 
521 template <class T> class FindSum;
522 
523 struct UnsuitableTypeForExpandElements {};
524 
525 template <class T>
526 struct ExpandElementResult
527 {
528  typedef UnsuitableTypeForExpandElements type;
529 };
530 
531 template <class T>
532 struct ExpandElementResult<std::complex<T> >
533 {
534  typedef T type;
535  enum { size = 2 };
536 };
537 
538 template <class T>
539 class FFTWComplex;
540 
541 template <class T>
542 struct ExpandElementResult<FFTWComplex<T> >
543 {
544  typedef T type;
545  enum { size = 2 };
546 };
547 
548 template <class T, int SIZE>
549 struct ExpandElementResult<TinyVector<T, SIZE> >
550 {
551  typedef T type;
552  enum { size = SIZE };
553 };
554 
555 template <class T, unsigned int R, unsigned int G, unsigned int B>
556 struct ExpandElementResult<RGBValue<T, R, G, B> >
557 {
558  typedef T type;
559  enum { size = 3 };
560 };
561 
562 #define VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(TYPE) \
563 template <> \
564 struct ExpandElementResult<TYPE> \
565 { \
566  typedef TYPE type; \
567  enum { size = 1 }; \
568 }; \
569 
570 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(bool)
571 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(char)
572 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(signed char)
573 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(signed short)
574 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(signed int)
575 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(signed long)
576 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(signed long long)
577 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(unsigned char)
578 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(unsigned short)
579 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(unsigned int)
580 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(unsigned long)
581 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(unsigned long long)
582 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(float)
583 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(double)
584 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(long double)
585 
586 #undef VIGRA_DEFINE_EXPAND_ELEMENT_RESULT
587 
588 
589 /********************************************************/
590 /* */
591 /* NormTraits */
592 /* */
593 /********************************************************/
594 
595 template <unsigned int N, class T, class C>
596 struct NormTraits<MultiArrayView<N, T, C> >
597 {
598  typedef MultiArrayView<N, T, C> Type;
599  typedef typename NormTraits<T>::SquaredNormType SquaredNormType;
600  typedef typename SquareRootTraits<SquaredNormType>::SquareRootResult NormType;
601 };
602 
603 template <unsigned int N, class T, class A>
604 struct NormTraits<MultiArray<N, T, A> >
605 : public NormTraits<typename MultiArray<N, T, A>::view_type>
606 {
607  typedef NormTraits<typename MultiArray<N, T, A>::view_type> BaseType;
608  typedef MultiArray<N, T, A> Type;
609  typedef typename BaseType::SquaredNormType SquaredNormType;
610  typedef typename BaseType::NormType NormType;
611 };
612 
613 /** \brief Base class for, and view to, \ref vigra::MultiArray.
614 
615 This class implements the interface of both MultiArray and
616 MultiArrayView. By default, MultiArrayViews are tagged as
617 strided (using <tt>StridedArrayTag</tt> as third template parameter).
618 This means that the array elements need not be consecutive in memory,
619 making the view flexible to represent all kinds of subarrays and slices.
620 In certain cases (which have become rare due to improvements of
621 optimizer and processor technology), an array may be tagged with
622 <tt>UnstridedArrayTag</tt> which indicates that the first array dimension
623 is guaranteed to be unstrided, i.e. has consecutive elements in memory.
624 
625 In addition to the member functions described here, <tt>MultiArrayView</tt>
626 and its subclasses support arithmetic and algebraic functions via the
627 module \ref MultiMathModule.
628 
629 If you want to apply an algorithm requiring an image to a
630 <tt>MultiArrayView</tt> of appropriate (2-dimensional) shape, you can
631 create a \ref vigra::BasicImageView that acts as a wrapper with the
632 necessary interface -- see \ref MultiArrayToImage.
633 
634 The template parameter are as follows
635 \code
636  N: the array dimension
637 
638  T: the type of the array elements
639 
640  C: a tag determining whether the array's inner dimension is strided
641  or not. An array is unstrided if the array elements occupy consecutive
642  memory location, strided if there is an offset in between (e.g.
643  when a view is created that skips every other array element).
644  The compiler can generate faster code for unstrided arrays.
645  Possible values: StridedArrayTag (default), UnstridedArrayTag
646 \endcode
647 
648 <b>\#include</b> <vigra/multi_array.hxx> <br/>
649 Namespace: vigra
650 */
651 template <unsigned int N, class T, class StrideTag>
653 {
654 public:
655 
656  /** the array's actual dimensionality.
657  This ensures that MultiArrayView can also be used for
658  scalars (that is, when <tt>N == 0</tt>). Calculated as:<br>
659  \code
660  actual_dimension = (N==0) ? 1 : N
661  \endcode
662  */
663  enum ActualDimension { actual_dimension = (N==0) ? 1 : N };
664 
665  /** the array's value type
666  */
667  typedef T value_type;
668 
669  /** reference type (result of operator[])
670  */
671  typedef value_type &reference;
672 
673  /** const reference type (result of operator[] const)
674  */
675  typedef const value_type &const_reference;
676 
677  /** pointer type
678  */
679  typedef value_type *pointer;
680 
681  /** const pointer type
682  */
683  typedef const value_type *const_pointer;
684 
685  /** difference type (used for multi-dimensional offsets and indices)
686  */
688 
689  /** key type (argument of index operator array[i] -- same as difference_type)
690  */
691  typedef difference_type key_type;
692 
693  /** size type
694  */
695  typedef difference_type size_type;
696 
697  /** difference and index type for a single dimension
698  */
700 
701  /** scan-order iterator (StridedScanOrderIterator) type
702  */
704 
705  /** const scan-order iterator (StridedScanOrderIterator) type
706  */
708 
709  /** traverser (MultiIterator) type
710  */
711  typedef typename vigra::detail::MultiIteratorChooser <
712  StrideTag>::template Traverser <actual_dimension, T, T &, T *>::type traverser;
713 
714  /** const traverser (MultiIterator) type
715  */
716  typedef typename vigra::detail::MultiIteratorChooser <
717  StrideTag>::template Traverser <actual_dimension, T, T const &, T const *>::type const_traverser;
718 
719  /** the view type associated with this array.
720  */
722 
723  /** the matrix type associated with this array.
724  */
726 
727  bool checkInnerStride(UnstridedArrayTag) const
728  {
729  return m_stride[0] <= 1;
730  }
731 
732  bool checkInnerStride(StridedArrayTag) const
733  {
734  return true;
735  }
736 
737  protected:
738 
739  typedef typename difference_type::value_type diff_zero_t;
740 
741  /** the shape of the image pointed to is stored here.
742  */
743  difference_type m_shape;
744 
745  /** the strides (offset of a sample to the next) for every dimension
746  are stored here.
747  */
748  difference_type m_stride;
749 
750  /** pointer to the image.
751  */
752  pointer m_ptr;
753 
754  template <class CN>
755  void assignImpl(const MultiArrayView <N, T, CN>& rhs);
756 
757  template <class U, class CN>
758  void copyImpl(const MultiArrayView <N, U, CN>& rhs);
759 
760  template <class U, class CN>
761  void swapDataImpl(MultiArrayView <N, U, CN> rhs);
762 
763  template <class CN>
764  bool arraysOverlap(const MultiArrayView <N, T, CN>& rhs) const;
765 
766  template <class U, class CN>
767  bool arraysOverlap(const MultiArrayView <N, U, CN>&) const
768  {
769  return false;
770  }
771 
772 public:
773 
774  /** default constructor: create an invalid view,
775  * i.e. hasData() returns false and size() is zero.
776  */
778  : m_shape (diff_zero_t(0)), m_stride (diff_zero_t(0)), m_ptr (0)
779  {}
780 
781  /** construct from another array view.
782  Throws a precondition error if this array has UnstridedArrayTag, but the
783  innermost dimension of \a other is strided.
784  */
785  template <class Stride>
787  : m_shape (other.shape()),
788  m_stride (other.stride()),
789  m_ptr (other.data())
790  {
791  vigra_precondition(other.checkInnerStride(StrideTag()),
792  "MultiArrayView<..., UnstridedArrayTag>(MultiArrayView const &): cannot create unstrided view from strided array.");
793  }
794 
795  /** construct from shape and pointer
796  */
797  MultiArrayView (const difference_type &shape, const_pointer ptr)
798  : m_shape (shape),
799  m_stride (detail::defaultStride<actual_dimension>(shape)),
800  m_ptr (const_cast<pointer>(ptr))
801  {}
802 
803  /** Construct from shape, strides (offset of a sample to the
804  next) for every dimension, and pointer. (Note that
805  strides are not given in bytes, but in offset steps of the
806  respective pointer type.)
807  */
808  MultiArrayView (const difference_type &shape,
809  const difference_type &stride,
810  const_pointer ptr)
811  : m_shape (shape),
812  m_stride (stride),
813  m_ptr (const_cast<pointer>(ptr))
814  {
815  vigra_precondition(checkInnerStride(StrideTag()),
816  "MultiArrayView<..., UnstridedArrayTag>::MultiArrayView(): First dimension of given array is not unstrided.");
817  }
818 
819  /** Construct from an old-style BasicImage.
820  */
821  template <class ALLOC>
823  : m_shape (Shape2(image.width(), image.height())),
824  m_stride (detail::defaultStride<actual_dimension>(m_shape)),
825  m_ptr (const_cast<pointer>(image.data()))
826  {}
827 
828  /** Conversion to a strided view.
829  */
831  {
833  }
834 
835  /** Assignment. There are 3 cases:
836 
837  <ul>
838  <li> When this <tt>MultiArrayView</tt> does not point to valid data
839  (e.g. after default construction), it becomes a new view of \a rhs.
840  <li> Otherwise, when the shapes of the two arrays match, the contents
841  (i.e. the elements) of \a rhs are copied.
842  <li> Otherwise, a <tt>PreconditionViolation</tt> exception is thrown.
843  </ul>
844  */
846  {
847  if(this != &rhs)
848  assignImpl(rhs);
849  return *this;
850  }
851 
852  template<class Stride2>
854  {
855  assignImpl(rhs);
856  return *this;
857  }
858 
859  /** Assignment of a differently typed MultiArrayView. It copies the elements
860  of\a rhs or fails with <tt>PreconditionViolation</tt> exception when
861  the shapes do not match.
862  */
863  template<class U, class C1>
865  {
866  vigra_precondition(this->shape() == rhs.shape(),
867  "MultiArrayView::operator=(): shape mismatch.");
868  this->copyImpl(rhs);
869  return *this;
870  }
871 
872  /** Assignment of a scalar. Equivalent to MultiArrayView::init(v).
873  */
874  MultiArrayView & operator=(value_type const & v)
875  {
876  return init(v);
877  }
878 
879  /** Add-assignment of a compatible MultiArrayView. Fails with
880  <tt>PreconditionViolation</tt> exception when the shapes do not match.
881  */
882  template<class U, class C1>
884 
885  /** Subtract-assignment of a compatible MultiArrayView. Fails with
886  <tt>PreconditionViolation</tt> exception when the shapes do not match.
887  */
888  template<class U, class C1>
890 
891  /** Multiply-assignment of a compatible MultiArrayView. Fails with
892  <tt>PreconditionViolation</tt> exception when the shapes do not match.
893  */
894  template<class U, class C1>
896 
897  /** Divide-assignment of a compatible MultiArrayView. Fails with
898  <tt>PreconditionViolation</tt> exception when the shapes do not match.
899  */
900  template<class U, class C1>
902 
903  /** Add-assignment of a scalar.
904  */
905  MultiArrayView & operator+=(T const & rhs)
906  {
907  detail::copyAddScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
908  return *this;
909  }
910 
911  /** Subtract-assignment of a scalar.
912  */
913  MultiArrayView & operator-=(T const & rhs)
914  {
915  detail::copySubScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
916  return *this;
917  }
918 
919  /** Multiply-assignment of a scalar.
920  */
921  MultiArrayView & operator*=(T const & rhs)
922  {
923  detail::copyMulScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
924  return *this;
925  }
926 
927  /** Divide-assignment of a scalar.
928  */
929  MultiArrayView & operator/=(T const & rhs)
930  {
931  detail::copyDivScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
932  return *this;
933  }
934 
935  /** Assignment of an array expression. Fails with
936  <tt>PreconditionViolation</tt> exception when the shapes do not match.
937  */
938  template<class Expression>
939  MultiArrayView & operator=(multi_math::MultiMathOperand<Expression> const & rhs)
940  {
941  multi_math::math_detail::assign(*this, rhs);
942  return *this;
943  }
944 
945  /** Add-assignment of an array expression. Fails with
946  <tt>PreconditionViolation</tt> exception when the shapes do not match.
947  */
948  template<class Expression>
949  MultiArrayView & operator+=(multi_math::MultiMathOperand<Expression> const & rhs)
950  {
951  multi_math::math_detail::plusAssign(*this, rhs);
952  return *this;
953  }
954 
955  /** Subtract-assignment of an array expression. Fails with
956  <tt>PreconditionViolation</tt> exception when the shapes do not match.
957  */
958  template<class Expression>
959  MultiArrayView & operator-=(multi_math::MultiMathOperand<Expression> const & rhs)
960  {
961  multi_math::math_detail::minusAssign(*this, rhs);
962  return *this;
963  }
964 
965  /** Multiply-assignment of an array expression. Fails with
966  <tt>PreconditionViolation</tt> exception when the shapes do not match.
967  */
968  template<class Expression>
969  MultiArrayView & operator*=(multi_math::MultiMathOperand<Expression> const & rhs)
970  {
971  multi_math::math_detail::multiplyAssign(*this, rhs);
972  return *this;
973  }
974 
975  /** Divide-assignment of an array expression. Fails with
976  <tt>PreconditionViolation</tt> exception when the shapes do not match.
977  */
978  template<class Expression>
979  MultiArrayView & operator/=(multi_math::MultiMathOperand<Expression> const & rhs)
980  {
981  multi_math::math_detail::divideAssign(*this, rhs);
982  return *this;
983  }
984 
985  /** array access.
986  */
987  reference operator[] (const difference_type &d)
988  {
989  VIGRA_ASSERT_INSIDE(d);
990  return m_ptr [dot (d, m_stride)];
991  }
992 
993  /** array access.
994  */
995  const_reference operator[] (const difference_type &d) const
996  {
997  VIGRA_ASSERT_INSIDE(d);
998  return m_ptr [dot (d, m_stride)];
999  }
1000 
1001  /** equivalent to bindInner(), when M < N.
1002  */
1003  template <int M>
1005  {
1006  return bindInner(d);
1007  }
1008 
1009  /** Array access in scan-order sense.
1010  Mostly useful to support standard indexing for 1-dimensional multi-arrays,
1011  but works for any N. Use scanOrderIndexToCoordinate() and
1012  coordinateToScanOrderIndex() for conversion between indices and coordinates.
1013 
1014  <b>Note:</b> This function should not be used in the inner loop, because the
1015  conversion of the scan order index into a memory address is expensive
1016  (it must take into account that memory may not be consecutive for subarrays
1017  and/or strided arrays). Always prefer operator() if possible.
1018  */
1019  reference operator[](difference_type_1 d)
1020  {
1021  VIGRA_ASSERT_INSIDE(scanOrderIndexToCoordinate(d));
1022  return m_ptr [detail::ScanOrderToOffset<actual_dimension>::exec(d, m_shape, m_stride)];
1023  }
1024 
1025  /** Array access in scan-order sense.
1026  Mostly useful to support standard indexing for 1-dimensional multi-arrays,
1027  but works for any N. Use scanOrderIndexToCoordinate() and
1028  coordinateToScanOrderIndex() for conversion between indices and coordinates.
1029 
1030  <b>Note:</b> This function should not be used in the inner loop, because the
1031  conversion of the scan order index into a memory address is expensive
1032  (it must take into account that memory may not be consecutive for subarrays
1033  and/or strided arrays). Always prefer operator() if possible.
1034  */
1035  const_reference operator[](difference_type_1 d) const
1036  {
1037  VIGRA_ASSERT_INSIDE(scanOrderIndexToCoordinate(d));
1038  return m_ptr [detail::ScanOrderToOffset<actual_dimension>::exec(d, m_shape, m_stride)];
1039  }
1040 
1041  /** convert scan-order index to coordinate.
1042  */
1043  difference_type scanOrderIndexToCoordinate(difference_type_1 d) const
1044  {
1045  difference_type result;
1046  detail::ScanOrderToCoordinate<actual_dimension>::exec(d, m_shape, result);
1047  return result;
1048  }
1049 
1050  /** convert coordinate to scan-order index.
1051  */
1052  difference_type_1 coordinateToScanOrderIndex(const difference_type &d) const
1053  {
1054  return detail::CoordinateToScanOrder<actual_dimension>::exec(m_shape, d);
1055  }
1056 
1057  /** 1D array access. Use only if N == 1.
1058  */
1059  reference operator() (difference_type_1 x)
1060  {
1061  VIGRA_ASSERT_INSIDE(difference_type(x));
1062  return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x)];
1063  }
1064 
1065  /** 2D array access. Use only if N == 2.
1066  */
1067  reference operator() (difference_type_1 x, difference_type_1 y)
1068  {
1069  VIGRA_ASSERT_INSIDE(difference_type(x, y));
1070  return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x, y)];
1071  }
1072 
1073  /** 3D array access. Use only if N == 3.
1074  */
1075  reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z)
1076  {
1077  VIGRA_ASSERT_INSIDE(difference_type(x, y, z));
1078  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
1079  }
1080 
1081  /** 4D array access. Use only if N == 4.
1082  */
1083  reference operator() (difference_type_1 x, difference_type_1 y,
1084  difference_type_1 z, difference_type_1 u)
1085  {
1086  VIGRA_ASSERT_INSIDE(difference_type(x, y, z, u));
1087  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
1088  }
1089 
1090  /** 5D array access. Use only if N == 5.
1091  */
1092  reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
1093  difference_type_1 u, difference_type_1 v)
1094  {
1095  VIGRA_ASSERT_INSIDE(difference_type(x, y,z, u,v));
1096  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
1097  }
1098 
1099  /** 1D const array access. Use only if N == 1.
1100  */
1101  const_reference operator() (difference_type_1 x) const
1102  {
1103  VIGRA_ASSERT_INSIDE(difference_type(x));
1104  return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x)];
1105  }
1106 
1107  /** 2D const array access. Use only if N == 2.
1108  */
1109  const_reference operator() (difference_type_1 x, difference_type_1 y) const
1110  {
1111  VIGRA_ASSERT_INSIDE(difference_type(x, y));
1112  return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x, y)];
1113  }
1114 
1115  /** 3D const array access. Use only if N == 3.
1116  */
1117  const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z) const
1118  {
1119  VIGRA_ASSERT_INSIDE(difference_type(x,y,z));
1120  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
1121  }
1122 
1123  /** 4D const array access. Use only if N == 4.
1124  */
1125  const_reference operator() (difference_type_1 x, difference_type_1 y,
1126  difference_type_1 z, difference_type_1 u) const
1127  {
1128  VIGRA_ASSERT_INSIDE(difference_type(x,y,z,u));
1129  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
1130  }
1131 
1132  /** 5D const array access. Use only if N == 5.
1133  */
1134  const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
1135  difference_type_1 u, difference_type_1 v) const
1136  {
1137  VIGRA_ASSERT_INSIDE(difference_type(x,y,z,u,v));
1138  return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
1139  }
1140 
1141  /** Init with a constant.
1142  */
1143  template <class U>
1144  MultiArrayView & init(const U & init)
1145  {
1146  if(hasData())
1147  detail::copyScalarMultiArrayData(traverser_begin(), shape(), init, MetaInt<actual_dimension-1>());
1148  return *this;
1149  }
1150 
1151 
1152  /** Copy the data of the right-hand array (array shapes must match).
1153  */
1154  void copy(const MultiArrayView & rhs)
1155  {
1156  if(this == &rhs)
1157  return;
1158  this->copyImpl(rhs);
1159  }
1160 
1161  /** Copy the data of the right-hand array (array shapes must match).
1162  */
1163  template <class U, class CN>
1165  {
1166  this->copyImpl(rhs);
1167  }
1168 
1169  /** swap the data between two MultiArrayView objects.
1170 
1171  The shapes of the two array must match.
1172  */
1174  {
1175  if(this != &rhs)
1176  swapDataImpl(rhs);
1177  }
1178 
1179  /** swap the data between two MultiArrayView objects.
1180 
1181  The shapes of the two array must match.
1182  */
1183  template <class T2, class C2>
1185  {
1186  swapDataImpl(rhs);
1187  }
1188 
1189  /** check whether the array is unstrided (i.e. has consecutive memory) up
1190  to the given dimension.
1191 
1192  \a dimension can range from 0 ... N-1. If a certain dimension is unstrided,
1193  all lower dimensions are also unstrided.
1194  */
1195  bool isUnstrided(unsigned int dimension = N-1) const
1196  {
1197  difference_type s = vigra::detail::defaultStride<actual_dimension>(shape());
1198  for(unsigned int k = 0; k <= dimension; ++k)
1199  if(stride(k) != s[k])
1200  return false;
1201  return true;
1202  }
1203 
1204  /** bind the M outmost dimensions to certain indices.
1205  this reduces the dimensionality of the image to
1206  max { 1, N-M }.
1207 
1208  <b>Usage:</b>
1209  \code
1210  // create a 3D array of size 40x30x20
1211  typedef MultiArray<3, double>::difference_type Shape;
1212  MultiArray<3, double> array3(Shape(40, 30, 20));
1213 
1214  // get a 1D array by fixing index 1 to 12, and index 2 to 10
1215  MultiArrayView <1, double> array1 = array3.bindOuter(TinyVector<MultiArrayIndex, 2>(12, 10));
1216  \endcode
1217  */
1218  template <int M, class Index>
1219  MultiArrayView <N-M, T, StrideTag> bindOuter(const TinyVector <Index, M> &d) const;
1220 
1221  /** bind the M innermost dimensions to certain indices.
1222  this reduces the dimensionality of the image to
1223  max { 1, N-M }.
1224 
1225  <b>Usage:</b>
1226  \code
1227  // create a 3D array of size 40x30x20
1228  typedef MultiArray<3, double>::difference_type Shape;
1229  MultiArray<3, double> array3(Shape(40, 30, 20));
1230 
1231  // get a 1D array by fixing index 0 to 12, and index 1 to 10
1232  MultiArrayView <1, double, StridedArrayTag> array1 = array3.bindInner(TinyVector<MultiArrayIndex, 2>(12, 10));
1233  \endcode
1234  */
1235  template <int M, class Index>
1237 
1238  /** bind dimension M to index d.
1239  this reduces the dimensionality of the image to
1240  max { 1, N-1 }.
1241 
1242  <b>Usage:</b>
1243  \code
1244  // create a 3D array of size 40x30x20
1245  typedef MultiArray<3, double>::difference_type Shape;
1246  MultiArray<3, double> array3(Shape(40, 30, 20));
1247 
1248  // get a 2D array by fixing index 1 to 12
1249  MultiArrayView <2, double> array2 = array3.bind<1>(12);
1250 
1251  // get a 2D array by fixing index 0 to 23
1252  MultiArrayView <2, double, StridedArrayTag> array2a = array3.bind<0>(23);
1253  \endcode
1254  */
1255  template <unsigned int M>
1256  MultiArrayView <N-1, T, typename vigra::detail::MaybeStrided<StrideTag, M>::type >
1257  bind (difference_type_1 d) const;
1258 
1259  /** bind the outmost dimension to a certain index.
1260  this reduces the dimensionality of the image to
1261  max { 1, N-1 }.
1262 
1263  <b>Usage:</b>
1264  \code
1265  // create a 3D array of size 40x30x20
1266  typedef MultiArray<3, double>::difference_type Shape;
1267  MultiArray<3, double> array3(Shape(40, 30, 20));
1268 
1269  // get a 2D array by fixing the outermost index (i.e. index 2) to 12
1270  MultiArrayView <2, double> array2 = array3.bindOuter(12);
1271  \endcode
1272  */
1273  MultiArrayView <N-1, T, StrideTag> bindOuter (difference_type_1 d) const;
1274 
1275  /** bind the innermost dimension to a certain index.
1276  this reduces the dimensionality of the image to
1277  max { 1, N-1 }.
1278 
1279  <b>Usage:</b>
1280  \code
1281  // create a 3D array of size 40x30x20
1282  typedef MultiArray<3, double>::difference_type Shape;
1283  MultiArray<3, double> array3(Shape(40, 30, 20));
1284 
1285  // get a 2D array by fixing the innermost index (i.e. index 0) to 23
1286  MultiArrayView <2, double, StridedArrayTag> array2 = array3.bindInner(23);
1287  \endcode
1288  */
1289  MultiArrayView <N-1, T, StridedArrayTag> bindInner (difference_type_1 d) const;
1290 
1291  /** bind dimension m to index d.
1292  this reduces the dimensionality of the image to
1293  max { 1, N-1 }.
1294 
1295  <b>Usage:</b>
1296  \code
1297  // create a 3D array of size 40x30x20
1298  typedef MultiArray<3, double>::difference_type Shape;
1299  MultiArray<3, double> array3(Shape(40, 30, 20));
1300 
1301  // get a 2D array by fixing index 2 to 15
1302  MultiArrayView <2, double, StridedArrayTag> array2 = array3.bindAt(2, 15);
1303  \endcode
1304  */
1306  bindAt (difference_type_1 m, difference_type_1 d) const;
1307 
1308  /** Create a view to channel 'i' of a vector-like value type. Possible value types
1309  (of the original array) are: \ref TinyVector, \ref RGBValue, \ref FFTWComplex,
1310  and <tt>std::complex</tt>. The list can be extended to any type whose memory
1311  layout is equivalent to a fixed-size C array, by specializing
1312  <tt>ExpandElementResult</tt>.
1313 
1314  <b>Usage:</b>
1315  \code
1316  MultiArray<2, RGBValue<float> > rgb_image(Shape2(w, h));
1317 
1318  MultiArrayView<2, float, StridedArrayTag> red = rgb_image.bindElementChannel(0);
1319  MultiArrayView<2, float, StridedArrayTag> green = rgb_image.bindElementChannel(1);
1320  MultiArrayView<2, float, StridedArrayTag> blue = rgb_image.bindElementChannel(2);
1321  \endcode
1322  */
1324  bindElementChannel(difference_type_1 i) const
1325  {
1326  vigra_precondition(0 <= i && i < ExpandElementResult<T>::size,
1327  "MultiArrayView::bindElementChannel(i): 'i' out of range.");
1328  return expandElements(0).bindInner(i);
1329  }
1330 
1331  /** Create a view where a vector-like element type is expanded into a new
1332  array dimension. The new dimension is inserted at index position 'd',
1333  which must be between 0 and N inclusive.
1334 
1335  Possible value types of the original array are: \ref TinyVector, \ref RGBValue,
1336  \ref FFTWComplex, <tt>std::complex</tt>, and the built-in number types (in this
1337  case, <tt>expandElements</tt> is equivalent to <tt>insertSingletonDimension</tt>).
1338  The list of supported types can be extended to any type whose memory
1339  layout is equivalent to a fixed-size C array, by specializing
1340  <tt>ExpandElementResult</tt>.
1341 
1342  <b>Usage:</b>
1343  \code
1344  MultiArray<2, RGBValue<float> > rgb_image(Shape2(w, h));
1345 
1346  MultiArrayView<3, float, StridedArrayTag> multiband_image = rgb_image.expandElements(2);
1347  \endcode
1348  */
1350  expandElements(difference_type_1 d) const;
1351 
1352  /** Add a singleton dimension (dimension of length 1).
1353 
1354  Singleton dimensions don't change the size of the data, but introduce
1355  a new index that can only take the value 0. This is mainly useful for
1356  the 'reduce mode' of transformMultiArray() and combineTwoMultiArrays(),
1357  because these functions require the source and destination arrays to
1358  have the same number of dimensions.
1359 
1360  The range of \a i must be <tt>0 <= i <= N</tt>. The new dimension will become
1361  the i'th index, and the old indices from i upwards will shift one
1362  place to the right.
1363 
1364  <b>Usage:</b>
1365 
1366  Suppose we want have a 2D array and want to create a 1D array that contains
1367  the row average of the first array.
1368  \code
1369  typedef MultiArrayShape<2>::type Shape2;
1370  MultiArray<2, double> original(Shape2(40, 30));
1371 
1372  typedef MultiArrayShape<1>::type Shape1;
1373  MultiArray<1, double> rowAverages(Shape1(30));
1374 
1375  // temporarily add a singleton dimension to the destination array
1376  transformMultiArray(srcMultiArrayRange(original),
1377  destMultiArrayRange(rowAverages.insertSingletonDimension(0)),
1378  FindAverage<double>());
1379  \endcode
1380  */
1382  insertSingletonDimension (difference_type_1 i) const;
1383 
1384  /** create a multiband view for this array.
1385 
1386  The type <tt>MultiArrayView<N, Multiband<T> ></tt> tells VIGRA
1387  algorithms which recognize the <tt>Multiband</tt> modifier to
1388  interpret the outermost (last) dimension as a channel dimension.
1389  In effect, these algorithms will treat the data as a set of
1390  (N-1)-dimensional arrays instead of a single N-dimensional array.
1391  */
1393  {
1394  return MultiArrayView<N, Multiband<value_type>, StrideTag>(*this);
1395  }
1396 
1397  /** Create a view to the diagonal elements of the array.
1398 
1399  This produces a 1D array view whose size equals the size
1400  of the shortest dimension of the original array.
1401 
1402  <b>Usage:</b>
1403  \code
1404  // create a 3D array of size 40x30x20
1405  typedef MultiArray<3, double>::difference_type Shape;
1406  MultiArray<3, double> array3(Shape(40, 30, 20));
1407 
1408  // get a view to the diagonal elements
1409  MultiArrayView <1, double, StridedArrayTag> diagonal = array3.diagonal();
1410  assert(diagonal.shape(0) == 20);
1411  \endcode
1412  */
1414  {
1416  Shape1(vigra::sum(m_stride)), m_ptr);
1417  }
1418 
1419  /** create a rectangular subarray that spans between the
1420  points p and q, where p is in the subarray, q not.
1421  If an element of p or q is negative, it is subtracted
1422  from the correspongng shape.
1423 
1424  <b>Usage:</b>
1425  \code
1426  // create a 3D array of size 40x30x20
1427  typedef MultiArray<3, double>::difference_type Shape;
1428  MultiArray<3, double> array3(Shape(40, 30, 20));
1429 
1430  // get a subarray set is smaller by one element at all sides
1431  MultiArrayView <3, double> subarray = array3.subarray(Shape(1,1,1), Shape(39, 29, 19));
1432 
1433  // specifying the end point with a vector of '-1' is equivalent
1434  MultiArrayView <3, double> subarray2 = array3.subarray(Shape(1,1,1), Shape(-1, -1, -1));
1435  \endcode
1436  */
1437  MultiArrayView subarray (difference_type p, difference_type q) const
1438  {
1439  detail::RelativeToAbsoluteCoordinate<actual_dimension-1>::exec(shape(), p);
1440  detail::RelativeToAbsoluteCoordinate<actual_dimension-1>::exec(shape(), q);
1441  const difference_type_1 offset = dot (m_stride, p);
1442  return MultiArrayView (q - p, m_stride, m_ptr + offset);
1443  }
1444 
1445  /** apply an additional striding to the image, thereby reducing
1446  the shape of the array.
1447  for example, multiplying the stride of dimension one by three
1448  turns an appropriately laid out (interleaved) rgb image into
1449  a single band image.
1450  */
1452  stridearray (const difference_type &s) const
1453  {
1454  difference_type shape = m_shape;
1455  for (unsigned int i = 0; i < actual_dimension; ++i)
1456  shape [i] /= s [i];
1457  return MultiArrayView <N, T, StridedArrayTag>(shape, m_stride * s, m_ptr);
1458  }
1459 
1460  /** Transpose an array. If N==2, this implements the usual matrix transposition.
1461  For N > 2, it reverses the order of the indices.
1462 
1463  <b>Usage:</b><br>
1464  \code
1465  typedef MultiArray<2, double>::difference_type Shape;
1466  MultiArray<2, double> array(10, 20);
1467 
1468  MultiArrayView<2, double, StridedArrayTag> transposed = array.transpose();
1469 
1470  for(int i=0; i<array.shape(0), ++i)
1471  for(int j=0; j<array.shape(1); ++j)
1472  assert(array(i, j) == transposed(j, i));
1473  \endcode
1474  */
1476  transpose () const
1477  {
1478  difference_type shape(m_shape.begin(), difference_type::ReverseCopy),
1479  stride(m_stride.begin(), difference_type::ReverseCopy);
1481  }
1482 
1483  /** Permute the dimensions of the array.
1484  The function exchanges the orer of the array's axes without copying the data.
1485  Argument\a permutation specifies the desired order such that
1486  <tt>permutation[k] = j</tt> means that axis <tt>j</tt> in the original array
1487  becomes axis <tt>k</tt> in the transposed array.
1488 
1489  <b>Usage:</b><br>
1490  \code
1491  typedef MultiArray<2, double>::difference_type Shape;
1492  MultiArray<2, double> array(10, 20);
1493 
1494  MultiArrayView<2, double, StridedArrayTag> transposed = array.transpose(Shape(1,0));
1495 
1496  for(int i=0; i<array.shape(0), ++i)
1497  for(int j=0; j<array.shape(1); ++j)
1498  assert(array(i, j) == transposed(j, i));
1499  \endcode
1500  */
1502  transpose(const difference_type &permutation) const
1503  {
1504  return permuteDimensions(permutation);
1505  }
1506 
1508  permuteDimensions (const difference_type &s) const;
1509 
1510  /** Permute the dimensions of the array so that the strides are in ascending order.
1511  Determines the appropriate permutation and then calls permuteDimensions().
1512  */
1514  permuteStridesAscending() const;
1515 
1516  /** Permute the dimensions of the array so that the strides are in descending order.
1517  Determines the appropriate permutation and then calls permuteDimensions().
1518  */
1520  permuteStridesDescending() const;
1521 
1522  /** Compute the ordering of the strides in this array.
1523  The result is describes the current permutation of the axes relative
1524  to the standard ascending stride order.
1525  */
1526  difference_type strideOrdering() const
1527  {
1528  return strideOrdering(m_stride);
1529  }
1530 
1531  /** Compute the ordering of the given strides.
1532  The result is describes the current permutation of the axes relative
1533  to the standard ascending stride order.
1534  */
1535  static difference_type strideOrdering(difference_type strides);
1536 
1537  /** number of the elements in the array.
1538  */
1539  difference_type_1 elementCount () const
1540  {
1541  difference_type_1 ret = m_shape[0];
1542  for(int i = 1; i < actual_dimension; ++i)
1543  ret *= m_shape[i];
1544  return ret;
1545  }
1546 
1547  /** number of the elements in the array.
1548  Same as <tt>elementCount()</tt>. Mostly useful to support the std::vector interface.
1549  */
1550  difference_type_1 size () const
1551  {
1552  return elementCount();
1553  }
1554 
1555  /** return the array's shape.
1556  */
1557  const difference_type & shape () const
1558  {
1559  return m_shape;
1560  }
1561 
1562  /** return the array's size at a certain dimension.
1563  */
1564  difference_type_1 size (difference_type_1 n) const
1565  {
1566  return m_shape [n];
1567  }
1568 
1569  /** return the array's shape at a certain dimension
1570  (same as <tt>size(n)</tt>).
1571  */
1572  difference_type_1 shape (difference_type_1 n) const
1573  {
1574  return m_shape [n];
1575  }
1576 
1577  /** return the array's width (same as <tt>shape(0)</tt>).
1578  */
1579  difference_type_1 width() const
1580  {
1581  return m_shape [0];
1582  }
1583 
1584  /** return the array's height (same as <tt>shape(1)</tt>).
1585  */
1586  difference_type_1 height() const
1587  {
1588  return m_shape [1];
1589  }
1590 
1591  /** return the array's stride for every dimension.
1592  */
1593  const difference_type & stride () const
1594  {
1595  return m_stride;
1596  }
1597 
1598  /** return the array's stride at a certain dimension.
1599  */
1600  difference_type_1 stride (int n) const
1601  {
1602  return m_stride [n];
1603  }
1604 
1605  /** check whether two arrays are elementwise equal.
1606  */
1607  template <class U, class C1>
1608  bool operator==(MultiArrayView<N, U, C1> const & rhs) const
1609  {
1610  if(this->shape() != rhs.shape())
1611  return false;
1612  return detail::equalityOfMultiArrays(traverser_begin(), shape(), rhs.traverser_begin(), MetaInt<actual_dimension-1>());
1613  }
1614 
1615  /** check whether two arrays are not elementwise equal.
1616  Also true when the two arrays have different shapes.
1617  */
1618  template <class U, class C1>
1619  bool operator!=(MultiArrayView<N, U, C1> const & rhs) const
1620  {
1621  return !operator==(rhs);
1622  }
1623 
1624  /** check whether the given point is in the array range.
1625  */
1626  bool isInside (difference_type const & p) const
1627  {
1628  for(int d=0; d<actual_dimension; ++d)
1629  if(p[d] < 0 || p[d] >= shape(d))
1630  return false;
1631  return true;
1632  }
1633  /** check whether the given point is not in the array range.
1634  */
1635  bool isOutside (difference_type const & p) const
1636  {
1637  for(int d=0; d<actual_dimension; ++d)
1638  if(p[d] < 0 || p[d] >= shape(d))
1639  return true;
1640  return false;
1641  }
1642 
1643  /** Check if the array contains only non-zero elements (or if all elements
1644  are 'true' if the value type is 'bool').
1645  */
1646  bool all() const
1647  {
1648  bool res = true;
1649  detail::reduceOverMultiArray(traverser_begin(), shape(),
1650  res,
1651  detail::AllTrueReduceFunctor(),
1652  MetaInt<actual_dimension-1>());
1653  return res;
1654  }
1655 
1656  /** Check if the array contains a non-zero element (or an element
1657  that is 'true' if the value type is 'bool').
1658  */
1659  bool any() const
1660  {
1661  bool res = false;
1662  detail::reduceOverMultiArray(traverser_begin(), shape(),
1663  res,
1664  detail::AnyTrueReduceFunctor(),
1665  MetaInt<actual_dimension-1>());
1666  return res;
1667  }
1668 
1669  /** Find the minimum and maximum element in this array.
1670  See \ref FeatureAccumulators for a general feature
1671  extraction framework.
1672  */
1673  void minmax(T * minimum, T * maximum) const
1674  {
1675  std::pair<T, T> res(NumericTraits<T>::max(), NumericTraits<T>::min());
1676  detail::reduceOverMultiArray(traverser_begin(), shape(),
1677  res,
1678  detail::MinmaxReduceFunctor(),
1679  MetaInt<actual_dimension-1>());
1680  *minimum = res.first;
1681  *maximum = res.second;
1682  }
1683 
1684  /** Compute the mean and variance of the values in this array.
1685  See \ref FeatureAccumulators for a general feature
1686  extraction framework.
1687  */
1688  template <class U>
1689  void meanVariance(U * mean, U * variance) const
1690  {
1691  typedef typename NumericTraits<U>::RealPromote R;
1692  R zero = R();
1693  triple<double, R, R> res(0.0, zero, zero);
1694  detail::reduceOverMultiArray(traverser_begin(), shape(),
1695  res,
1696  detail::MeanVarianceReduceFunctor(),
1697  MetaInt<actual_dimension-1>());
1698  *mean = res.second;
1699  *variance = res.third / res.first;
1700  }
1701 
1702  /** Compute the sum of the array elements.
1703 
1704  You must provide the type of the result by an explicit template parameter:
1705  \code
1706  MultiArray<2, UInt8> A(width, height);
1707 
1708  double sum = A.sum<double>();
1709  \endcode
1710  */
1711  template <class U>
1712  U sum() const
1713  {
1714  U res = NumericTraits<U>::zero();
1715  detail::reduceOverMultiArray(traverser_begin(), shape(),
1716  res,
1717  detail::SumReduceFunctor(),
1718  MetaInt<actual_dimension-1>());
1719  return res;
1720  }
1721 
1722  /** Compute the sum of the array elements over selected axes.
1723 
1724  \arg sums must have the same shape as this array, except for the
1725  axes along which the sum is to be accumulated. These axes must be
1726  singletons. Note that you must include <tt>multi_pointoperators.hxx</tt>
1727  for this function to work.
1728 
1729  <b>Usage:</b>
1730  \code
1731  #include <vigra/multi_array.hxx>
1732  #include <vigra/multi_pointoperators.hxx>
1733 
1734  MultiArray<2, double> A(Shape2(rows, cols));
1735  ... // fill A
1736 
1737  // make the first axis a singleton to sum over the first index
1738  MultiArray<2, double> rowSums(Shape2(1, cols));
1739  A.sum(rowSums);
1740 
1741  // this is equivalent to
1742  transformMultiArray(srcMultiArrayRange(A),
1743  destMultiArrayRange(rowSums),
1744  FindSum<double>());
1745  \endcode
1746  */
1747  template <class U, class S>
1748  void sum(MultiArrayView<N, U, S> sums) const
1749  {
1750  transformMultiArray(srcMultiArrayRange(*this),
1751  destMultiArrayRange(sums),
1752  FindSum<U>());
1753  }
1754 
1755  /** Compute the product of the array elements.
1756 
1757  You must provide the type of the result by an explicit template parameter:
1758  \code
1759  MultiArray<2, UInt8> A(width, height);
1760 
1761  double prod = A.product<double>();
1762  \endcode
1763  */
1764  template <class U>
1765  U product() const
1766  {
1767  U res = NumericTraits<U>::one();
1768  detail::reduceOverMultiArray(traverser_begin(), shape(),
1769  res,
1770  detail::ProdReduceFunctor(),
1771  MetaInt<actual_dimension-1>());
1772  return res;
1773  }
1774 
1775  /** Compute the squared Euclidean norm of the array (sum of squares of the array elements).
1776  */
1777  typename NormTraits<MultiArrayView>::SquaredNormType
1778  squaredNorm() const
1779  {
1780  typedef typename NormTraits<MultiArrayView>::SquaredNormType SquaredNormType;
1781  SquaredNormType res = NumericTraits<SquaredNormType>::zero();
1782  detail::reduceOverMultiArray(traverser_begin(), shape(),
1783  res,
1784  detail::SquaredL2NormReduceFunctor(),
1785  MetaInt<actual_dimension-1>());
1786  return res;
1787  }
1788 
1789  /** Compute various norms of the array.
1790  The norm is determined by parameter \a type:
1791 
1792  <ul>
1793  <li> type == 0: maximum norm (L-infinity): maximum of absolute values of the array elements
1794  <li> type == 1: Manhattan norm (L1): sum of absolute values of the array elements
1795  <li> type == 2: Euclidean norm (L2): square root of <tt>squaredNorm()</tt> when \a useSquaredNorm is <tt>true</tt>,<br>
1796  or direct algorithm that avoids underflow/overflow otherwise.
1797  </ul>
1798 
1799  Parameter \a useSquaredNorm has no effect when \a type != 2. Defaults: compute L2 norm as square root of
1800  <tt>squaredNorm()</tt>.
1801  */
1802  typename NormTraits<MultiArrayView>::NormType
1803  norm(int type = 2, bool useSquaredNorm = true) const;
1804 
1805  /** return the pointer to the image data
1806  */
1807  pointer data () const
1808  {
1809  return m_ptr;
1810  }
1811 
1812  pointer & unsafePtr()
1813  {
1814  return m_ptr;
1815  }
1816 
1817  /**
1818  * returns true iff this view refers to valid data,
1819  * i.e. data() is not a NULL pointer. (this is false after
1820  * default construction.)
1821  */
1822  bool hasData () const
1823  {
1824  return m_ptr != 0;
1825  }
1826 
1827  /** returns a scan-order iterator pointing
1828  to the first array element.
1829  */
1830  iterator begin()
1831  {
1832  return iterator(*this);
1833  }
1834 
1835  /** returns a const scan-order iterator pointing
1836  to the first array element.
1837  */
1838  const_iterator begin() const
1839  {
1840  return const_iterator(*this);
1841  }
1842 
1843  /** returns a scan-order iterator pointing
1844  beyond the last array element.
1845  */
1846  iterator end()
1847  {
1848  return begin().getEndIterator();
1849  }
1850 
1851  /** returns a const scan-order iterator pointing
1852  beyond the last array element.
1853  */
1854  const_iterator end() const
1855  {
1856  return begin().getEndIterator();
1857  }
1858 
1859  /** returns the N-dimensional MultiIterator pointing
1860  to the first element in every dimension.
1861  */
1862  traverser traverser_begin ()
1863  {
1864  traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1865  return ret;
1866  }
1867 
1868  /** returns the N-dimensional MultiIterator pointing
1869  to the const first element in every dimension.
1870  */
1871  const_traverser traverser_begin () const
1872  {
1873  const_traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1874  return ret;
1875  }
1876 
1877  /** returns the N-dimensional MultiIterator pointing
1878  beyond the last element in dimension N, and to the
1879  first element in every other dimension.
1880  */
1881  traverser traverser_end ()
1882  {
1883  traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1884  ret += m_shape [actual_dimension-1];
1885  return ret;
1886  }
1887 
1888  /** returns the N-dimensional const MultiIterator pointing
1889  beyond the last element in dimension N, and to the
1890  first element in every other dimension.
1891  */
1892  const_traverser traverser_end () const
1893  {
1894  const_traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
1895  ret += m_shape [actual_dimension-1];
1896  return ret;
1897  }
1898 
1899  view_type view () const
1900  {
1901  return *this;
1902  }
1903 };
1904 
1905 template <unsigned int N, class T, class StrideTag>
1906 class MultiArrayView<N, Multiband<T>, StrideTag>
1907 : public MultiArrayView<N, T, StrideTag>
1908 {
1909  public:
1910  MultiArrayView(MultiArrayView<N, T, StrideTag> const & v)
1911  : MultiArrayView<N, T, StrideTag>(v)
1912  {}
1913 };
1914 
1915 
1916 template <unsigned int N, class T, class Stride1>
1917 template <class Stride2>
1918 void
1919 MultiArrayView <N, T, Stride1>::assignImpl(MultiArrayView<N, T, Stride2> const & rhs)
1920 {
1921  if(m_ptr == 0)
1922  {
1923  vigra_precondition(rhs.checkInnerStride(Stride1()),
1924  "MultiArrayView<..., UnstridedArrayTag>::operator=(MultiArrayView const &): cannot create unstrided view from strided array.");
1925 
1926  m_shape = rhs.shape();
1927  m_stride = rhs.stride();
1928  m_ptr = rhs.data();
1929  }
1930  else
1931  {
1932  vigra_precondition(this->shape() == rhs.shape(),
1933  "MultiArrayView::operator=(MultiArrayView const &): shape mismatch.");
1934  this->copyImpl(rhs);
1935  }
1936 }
1937 
1938 template <unsigned int N, class T, class StrideTag>
1939 template <class CN>
1940 bool
1941 MultiArrayView <N, T, StrideTag>::arraysOverlap(const MultiArrayView <N, T, CN>& rhs) const
1942 {
1943  vigra_precondition (shape () == rhs.shape (),
1944  "MultiArrayView::arraysOverlap(): shape mismatch.");
1945  const_pointer first_element = this->m_ptr,
1946  last_element = first_element + dot(this->m_shape - difference_type(1), this->m_stride);
1947  typename MultiArrayView <N, T, CN>::const_pointer
1948  rhs_first_element = rhs.data(),
1949  rhs_last_element = rhs_first_element + dot(rhs.shape() - difference_type(1), rhs.stride());
1950  return !(last_element < rhs_first_element || rhs_last_element < first_element);
1951 }
1952 
1953 template <unsigned int N, class T, class StrideTag>
1954 template <class U, class CN>
1955 void
1956 MultiArrayView <N, T, StrideTag>::copyImpl(const MultiArrayView <N, U, CN>& rhs)
1957 {
1958  if(!arraysOverlap(rhs))
1959  {
1960  // no overlap -- can copy directly
1961  detail::copyMultiArrayData(rhs.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>());
1962  }
1963  else
1964  {
1965  // overlap: we got different views to the same data -- copy to intermediate memory in order to avoid
1966  // overwriting elements that are still needed on the rhs.
1967  MultiArray<N, T> tmp(rhs);
1968  detail::copyMultiArrayData(tmp.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>());
1969  }
1970 }
1971 
1972 #define VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(name, op) \
1973 template <unsigned int N, class T, class StrideTag> \
1974 template<class U, class C1> \
1975 MultiArrayView<N, T, StrideTag> & \
1976 MultiArrayView <N, T, StrideTag>::operator op(MultiArrayView<N, U, C1> const & rhs) \
1977 { \
1978  vigra_precondition(this->shape() == rhs.shape(), "MultiArrayView::operator" #op "() size mismatch."); \
1979  if(!arraysOverlap(rhs)) \
1980  { \
1981  detail::name##MultiArrayData(rhs.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>()); \
1982  } \
1983  else \
1984  { \
1985  MultiArray<N, T> tmp(rhs); \
1986  detail::name##MultiArrayData(tmp.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>()); \
1987  } \
1988  return *this; \
1989 }
1990 
1991 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copyAdd, +=)
1992 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copySub, -=)
1993 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copyMul, *=)
1994 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copyDiv, /=)
1995 
1996 #undef VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT
1997 
1998 template <unsigned int N, class T, class StrideTag>
1999 template <class U, class CN>
2000 void
2001 MultiArrayView <N, T, StrideTag>::swapDataImpl(MultiArrayView <N, U, CN> rhs)
2002 {
2003  vigra_precondition (shape () == rhs.shape (),
2004  "MultiArrayView::swapData(): shape mismatch.");
2005 
2006  // check for overlap of this and rhs
2007  const_pointer first_element = this->m_ptr,
2008  last_element = first_element + dot(this->m_shape - difference_type(1), this->m_stride);
2009  typename MultiArrayView <N, U, CN>::const_pointer
2010  rhs_first_element = rhs.data(),
2011  rhs_last_element = rhs_first_element + dot(rhs.shape() - difference_type(1), rhs.stride());
2012  if(last_element < rhs_first_element || rhs_last_element < first_element)
2013  {
2014  // no overlap -- can swap directly
2015  detail::swapDataImpl(traverser_begin(), shape(), rhs.traverser_begin(), MetaInt<actual_dimension-1>());
2016  }
2017  else
2018  {
2019  // overlap: we got different views to the same data -- copy to intermediate memory in order to avoid
2020  // overwriting elements that are still needed.
2021  MultiArray<N, T> tmp(*this);
2022  copy(rhs);
2023  rhs.copy(tmp);
2024  }
2025 }
2026 
2027 template <unsigned int N, class T, class StrideTag>
2028 MultiArrayView <N, T, StridedArrayTag>
2029 MultiArrayView <N, T, StrideTag>::permuteDimensions (const difference_type &s) const
2030 {
2031  difference_type shape, stride, check((typename difference_type::value_type)0);
2032  for (unsigned int i = 0; i < actual_dimension; ++i)
2033  {
2034  shape[i] = m_shape[s[i]];
2035  stride[i] = m_stride[s[i]];
2036  ++check[s[i]];
2037  }
2038  vigra_precondition(check == difference_type(1),
2039  "MultiArrayView::transpose(): every dimension must occur exactly once.");
2040  return MultiArrayView <N, T, StridedArrayTag>(shape, stride, m_ptr);
2041 }
2042 
2043 template <unsigned int N, class T, class StrideTag>
2044 typename MultiArrayView <N, T, StrideTag>::difference_type
2046 {
2047  difference_type permutation;
2048  for(int k=0; k<(int)N; ++k)
2049  permutation[k] = k;
2050  for(int k=0; k<(int)N-1; ++k)
2051  {
2052  int smallest = k;
2053  for(int j=k+1; j<(int)N; ++j)
2054  {
2055  if(stride[j] < stride[smallest])
2056  smallest = j;
2057  }
2058  if(smallest != k)
2059  {
2060  std::swap(stride[k], stride[smallest]);
2061  std::swap(permutation[k], permutation[smallest]);
2062  }
2063  }
2064  difference_type ordering;
2065  for(unsigned int k=0; k<N; ++k)
2066  ordering[permutation[k]] = k;
2067  return ordering;
2068 }
2069 
2070 template <unsigned int N, class T, class StrideTag>
2073 {
2074  difference_type ordering(strideOrdering(m_stride)), permutation;
2075  for(MultiArrayIndex k=0; k<N; ++k)
2076  permutation[ordering[k]] = k;
2077  return permuteDimensions(permutation);
2078 }
2079 
2080 template <unsigned int N, class T, class StrideTag>
2083 {
2084  difference_type ordering(strideOrdering(m_stride)), permutation;
2085  for(MultiArrayIndex k=0; k<N; ++k)
2086  permutation[N-1-ordering[k]] = k;
2087  return permuteDimensions(permutation);
2088 }
2089 
2090 template <unsigned int N, class T, class StrideTag>
2091 template <int M, class Index>
2092 MultiArrayView <N-M, T, StrideTag>
2094 {
2096  stride.init (m_stride.begin () + N-M, m_stride.end ());
2097  pointer ptr = m_ptr + dot (d, stride);
2098  static const int NNew = (N-M == 0) ? 1 : N-M;
2099  TinyVector <MultiArrayIndex, NNew> inner_shape, inner_stride;
2100  if (N-M == 0)
2101  {
2102  inner_shape [0] = 1;
2103  inner_stride [0] = 1;
2104  }
2105  else
2106  {
2107  inner_shape.init (m_shape.begin (), m_shape.end () - M);
2108  inner_stride.init (m_stride.begin (), m_stride.end () - M);
2109  }
2110  return MultiArrayView <N-M, T, StrideTag> (inner_shape, inner_stride, ptr);
2111 }
2112 
2113 template <unsigned int N, class T, class StrideTag>
2114 template <int M, class Index>
2115 MultiArrayView <N - M, T, StridedArrayTag>
2117 {
2119  stride.init (m_stride.begin (), m_stride.end () - N + M);
2120  pointer ptr = m_ptr + dot (d, stride);
2121  static const int NNew = (N-M == 0) ? 1 : N-M;
2122  TinyVector <MultiArrayIndex, NNew> outer_shape, outer_stride;
2123  if (N-M == 0)
2124  {
2125  outer_shape [0] = 1;
2126  outer_stride [0] = 1;
2127  }
2128  else
2129  {
2130  outer_shape.init (m_shape.begin () + M, m_shape.end ());
2131  outer_stride.init (m_stride.begin () + M, m_stride.end ());
2132  }
2133  return MultiArrayView <N-M, T, StridedArrayTag>
2134  (outer_shape, outer_stride, ptr);
2135 }
2136 
2137 template <unsigned int N, class T, class StrideTag>
2138 template <unsigned int M>
2139 MultiArrayView <N-1, T, typename detail::MaybeStrided<StrideTag, M>::type >
2141 {
2142  static const int NNew = (N-1 == 0) ? 1 : N-1;
2144  // the remaining dimensions are 0..n-1,n+1..N-1
2145  if (N-1 == 0)
2146  {
2147  shape[0] = 1;
2148  stride[0] = 1;
2149  }
2150  else
2151  {
2152  std::copy (m_shape.begin (), m_shape.begin () + M, shape.begin ());
2153  std::copy (m_shape.begin () + M+1, m_shape.end (),
2154  shape.begin () + M);
2155  std::copy (m_stride.begin (), m_stride.begin () + M, stride.begin ());
2156  std::copy (m_stride.begin () + M+1, m_stride.end (),
2157  stride.begin () + M);
2158  }
2159  return MultiArrayView <N-1, T, typename detail::MaybeStrided<StrideTag, M>::type>
2160  (shape, stride, m_ptr + d * m_stride[M]);
2161 }
2162 
2163 template <unsigned int N, class T, class StrideTag>
2164 MultiArrayView <N - 1, T, StrideTag>
2166 {
2167  static const int NNew = (N-1 == 0) ? 1 : N-1;
2168  TinyVector <MultiArrayIndex, NNew> inner_shape, inner_stride;
2169  if (N-1 == 0)
2170  {
2171  inner_shape [0] = 1;
2172  inner_stride [0] = 1;
2173  }
2174  else
2175  {
2176  inner_shape.init (m_shape.begin (), m_shape.end () - 1);
2177  inner_stride.init (m_stride.begin (), m_stride.end () - 1);
2178  }
2179  return MultiArrayView <N-1, T, StrideTag> (inner_shape, inner_stride,
2180  m_ptr + d * m_stride [N-1]);
2181 }
2182 
2183 template <unsigned int N, class T, class StrideTag>
2184 MultiArrayView <N - 1, T, StridedArrayTag>
2186 {
2187  static const int NNew = (N-1 == 0) ? 1 : N-1;
2188  TinyVector <MultiArrayIndex, NNew> outer_shape, outer_stride;
2189  if (N-1 == 0)
2190  {
2191  outer_shape [0] = 1;
2192  outer_stride [0] = 1;
2193  }
2194  else
2195  {
2196  outer_shape.init (m_shape.begin () + 1, m_shape.end ());
2197  outer_stride.init (m_stride.begin () + 1, m_stride.end ());
2198  }
2199  return MultiArrayView <N-1, T, StridedArrayTag>
2200  (outer_shape, outer_stride, m_ptr + d * m_stride [0]);
2201 }
2202 
2203 template <unsigned int N, class T, class StrideTag>
2204 MultiArrayView <N - 1, T, StridedArrayTag>
2206 {
2207  vigra_precondition (
2208  n < static_cast <int> (N),
2209  "MultiArrayView <N, T, StrideTag>::bindAt(): dimension out of range.");
2210  static const int NNew = (N-1 == 0) ? 1 : N-1;
2212  // the remaining dimensions are 0..n-1,n+1..N-1
2213  if (N-1 == 0)
2214  {
2215  shape [0] = 1;
2216  stride [0] = 1;
2217  }
2218  else
2219  {
2220  std::copy (m_shape.begin (), m_shape.begin () + n, shape.begin ());
2221  std::copy (m_shape.begin () + n+1, m_shape.end (),
2222  shape.begin () + n);
2223  std::copy (m_stride.begin (), m_stride.begin () + n, stride.begin ());
2224  std::copy (m_stride.begin () + n+1, m_stride.end (),
2225  stride.begin () + n);
2226  }
2227  return MultiArrayView <N-1, T, StridedArrayTag>
2228  (shape, stride, m_ptr + d * m_stride[n]);
2229 }
2230 
2231 
2232 template <unsigned int N, class T, class StrideTag>
2235 {
2236  vigra_precondition(0 <= d && d <= static_cast <difference_type_1> (N),
2237  "MultiArrayView<N, ...>::expandElements(d): 0 <= 'd' <= N required.");
2238 
2239  int elementSize = ExpandElementResult<T>::size;
2240  typename MultiArrayShape<N+1>::type newShape, newStrides;
2241  for(int k=0; k<d; ++k)
2242  {
2243  newShape[k] = m_shape[k];
2244  newStrides[k] = m_stride[k]*elementSize;
2245  }
2246 
2247  newShape[d] = elementSize;
2248  newStrides[d] = 1;
2249 
2250  for(int k=d; k<N; ++k)
2251  {
2252  newShape[k+1] = m_shape[k];
2253  newStrides[k+1] = m_stride[k]*elementSize;
2254  }
2255 
2256  typedef typename ExpandElementResult<T>::type U;
2258  newShape, newStrides, reinterpret_cast<U*>(m_ptr));
2259 }
2260 
2261 template <unsigned int N, class T, class StrideTag>
2264 {
2265  vigra_precondition (
2266  0 <= i && i <= static_cast <difference_type_1> (N),
2267  "MultiArrayView <N, T, StrideTag>::insertSingletonDimension(): index out of range.");
2269  std::copy (m_shape.begin (), m_shape.begin () + i, shape.begin ());
2270  std::copy (m_shape.begin () + i, m_shape.end (), shape.begin () + i + 1);
2271  std::copy (m_stride.begin (), m_stride.begin () + i, stride.begin ());
2272  std::copy (m_stride.begin () + i, m_stride.end (), stride.begin () + i + 1);
2273  shape[i] = 1;
2274  stride[i] = 1;
2275 
2277 }
2278 
2279 template <unsigned int N, class T, class StrideTag>
2280 typename NormTraits<MultiArrayView <N, T, StrideTag> >::NormType
2281 MultiArrayView <N, T, StrideTag>::norm(int type, bool useSquaredNorm) const
2282 {
2283  typedef typename NormTraits<MultiArrayView>::NormType NormType;
2284 
2285  switch(type)
2286  {
2287  case 0:
2288  {
2289  NormType res = NumericTraits<NormType>::zero();
2290  detail::reduceOverMultiArray(traverser_begin(), shape(),
2291  res,
2292  detail::MaxNormReduceFunctor(),
2293  MetaInt<actual_dimension-1>());
2294  return res;
2295  }
2296  case 1:
2297  {
2298  NormType res = NumericTraits<NormType>::zero();
2299  detail::reduceOverMultiArray(traverser_begin(), shape(),
2300  res,
2301  detail::L1NormReduceFunctor(),
2302  MetaInt<actual_dimension-1>());
2303  return res;
2304  }
2305  case 2:
2306  {
2307  if(useSquaredNorm)
2308  {
2309  return sqrt((NormType)squaredNorm());
2310  }
2311  else
2312  {
2313  NormType normMax = NumericTraits<NormType>::zero();
2314  detail::reduceOverMultiArray(traverser_begin(), shape(),
2315  normMax,
2316  detail::MaxNormReduceFunctor(),
2317  MetaInt<actual_dimension-1>());
2318  if(normMax == NumericTraits<NormType>::zero())
2319  return normMax;
2320  NormType res = NumericTraits<NormType>::zero();
2321  detail::reduceOverMultiArray(traverser_begin(), shape(),
2322  res,
2323  detail::WeightedL2NormReduceFunctor<NormType>(1.0/normMax),
2324  MetaInt<actual_dimension-1>());
2325  return sqrt(res)*normMax;
2326  }
2327  }
2328  default:
2329  vigra_precondition(false, "MultiArrayView::norm(): Unknown norm type.");
2330  return NumericTraits<NormType>::zero(); // unreachable
2331  }
2332 }
2333 
2334 
2335 /********************************************************/
2336 /* */
2337 /* norm */
2338 /* */
2339 /********************************************************/
2340 
2341 template <unsigned int N, class T, class StrideTag>
2342 inline typename NormTraits<MultiArrayView <N, T, StrideTag> >::SquaredNormType
2344 {
2345  return a.squaredNorm();
2346 }
2347 
2348 template <unsigned int N, class T, class StrideTag>
2349 inline typename NormTraits<MultiArrayView <N, T, StrideTag> >::NormType
2350 norm(MultiArrayView <N, T, StrideTag> const & a)
2351 {
2352  return a.norm();
2353 }
2354 
2355 /********************************************************/
2356 /* */
2357 /* MultiArray */
2358 /* */
2359 /********************************************************/
2360 
2361 /** \brief Main <TT>MultiArray</TT> class containing the memory
2362  management.
2363 
2364 This class inherits the interface of MultiArrayView, and implements
2365 the memory ownership.
2366 MultiArray's are always unstrided, striding them creates a MultiArrayView.
2367 
2368 
2369 The template parameters are as follows
2370 \code
2371  N: the array dimension
2372 
2373  T: the type of the array elements
2374 
2375  A: the allocator used for internal storage management
2376  (default: std::allocator<T>)
2377 \endcode
2378 
2379 <b>\#include</b> <vigra/multi_array.hxx> <br/>
2380 Namespace: vigra
2381 */
2382 template <unsigned int N, class T, class A /* default already declared above */>
2384 : public MultiArrayView <N, typename vigra::detail::ResolveMultiband<T>::type,
2385  typename vigra::detail::ResolveMultiband<T>::Stride>
2386 {
2387  public:
2389 
2390  /** the view type associated with this array.
2391  */
2394 
2395  using view_type::actual_dimension;
2396 
2397  /** the allocator type used to allocate the memory
2398  */
2399  typedef A allocator_type;
2400 
2401  /** the matrix type associated with this array.
2402  */
2404 
2405  /** the array's value type
2406  */
2408 
2409  /** pointer type
2410  */
2411  typedef typename view_type::pointer pointer;
2412 
2413  /** const pointer type
2414  */
2416 
2417  /** reference type (result of operator[])
2418  */
2420 
2421  /** const reference type (result of operator[] const)
2422  */
2424 
2425  /** size type
2426  */
2428 
2429  /** difference type (used for multi-dimensional offsets and indices)
2430  */
2432 
2433  /** difference and index type for a single dimension
2434  */
2436 
2437  /** traverser type
2438  */
2440 
2441  /** traverser type to const data
2442  */
2444 
2445  // /** sequential (random access) iterator type
2446  // */
2447  // typedef typename vigra::detail::MultiIteratorChooser<actual_stride>::template Iterator<N, value_type, reference, pointer>::type
2448  // iterator;
2449 
2450  // /** sequential (random access) const iterator type
2451  // */
2452  // typedef typename vigra::detail::MultiIteratorChooser<actual_stride>::template Iterator<N, value_type, const_reference, const_pointer>::type
2453  // const_iterator;
2454 
2455  /** sequential (random access) iterator type
2456  */
2457  typedef typename view_type::iterator iterator;
2458 
2459  /** sequential (random access) const iterator type
2460  */
2462 
2463 protected:
2464 
2465  typedef typename difference_type::value_type diff_zero_t;
2466 
2467  /** the allocator used to allocate the memory
2468  */
2470 
2471  /** allocate memory for s pixels, write its address into the given
2472  pointer and initialize the pixels with init.
2473  */
2474  void allocate (pointer &ptr, difference_type_1 s, const_reference init);
2475 
2476  /** allocate memory for s pixels, write its address into the given
2477  pointer and initialize the linearized pixels to the values of init.
2478  */
2479  template <class U>
2480  void allocate (pointer &ptr, difference_type_1 s, U const * init);
2481 
2482  /** allocate memory, write its address into the given
2483  pointer and initialize it by copying the data from the given MultiArrayView.
2484  */
2485  template <class U, class StrideTag>
2486  void allocate (pointer &ptr, MultiArrayView<N, U, StrideTag> const & init);
2487 
2488  /** deallocate the memory (of length s) starting at the given address.
2489  */
2490  void deallocate (pointer &ptr, difference_type_1 s);
2491 
2492  template <class U, class StrideTag>
2493  void copyOrReshape (const MultiArrayView<N, U, StrideTag> &rhs);
2494 public:
2495  /** default constructor
2496  */
2498  : view_type (difference_type (diff_zero_t(0)),
2499  difference_type (diff_zero_t(0)), 0)
2500  {}
2501 
2502  /** construct with given allocator
2503  */
2504  MultiArray (allocator_type const & alloc)
2505  : view_type(difference_type (diff_zero_t(0)),
2506  difference_type (diff_zero_t(0)), 0),
2507  m_alloc(alloc)
2508  {}
2509 
2510  /** construct with given length
2511 
2512  Use only for 1-dimensional arrays (<tt>N==1</tt>).
2513  */
2514  explicit MultiArray (difference_type_1 length,
2515  allocator_type const & alloc = allocator_type());
2516 
2517 
2518  /** construct with given width and height
2519 
2520  Use only for 2-dimensional arrays (<tt>N==2</tt>).
2521  */
2523  allocator_type const & alloc = allocator_type());
2524 
2525  /** construct with given shape
2526  */
2527  explicit MultiArray (const difference_type &shape,
2528  allocator_type const & alloc = allocator_type());
2529 
2530  /** construct from shape with an initial value
2531  */
2532  MultiArray (const difference_type &shape, const_reference init,
2533  allocator_type const & alloc = allocator_type());
2534 
2535  /** construct from shape and initialize with a linear sequence in scan order
2536  (i.e. first pixel gets value 0, second on gets value 1 and so on).
2537  */
2539  allocator_type const & alloc = allocator_type());
2540 
2541  /** construct from shape and copy values from the given array
2542  */
2543  MultiArray (const difference_type &shape, const_pointer init,
2544  allocator_type const & alloc = allocator_type());
2545 
2546  /** copy constructor
2547  */
2548  MultiArray (const MultiArray &rhs)
2549  : view_type(rhs.m_shape, rhs.m_stride, 0),
2550  m_alloc (rhs.m_alloc)
2551  {
2552  allocate (this->m_ptr, this->elementCount (), rhs.data ());
2553  }
2554 
2555  /** constructor from an array expression
2556  */
2557  template<class Expression>
2558  MultiArray (multi_math::MultiMathOperand<Expression> const & rhs,
2559  allocator_type const & alloc = allocator_type())
2560  : view_type(difference_type (diff_zero_t(0)),
2561  difference_type (diff_zero_t(0)), 0),
2562  m_alloc (alloc)
2563  {
2564  multi_math::math_detail::assignOrResize(*this, rhs);
2565  }
2566 
2567  /** construct by copying from a MultiArrayView
2568  */
2569  template <class U, class StrideTag>
2571  allocator_type const & alloc = allocator_type());
2572 
2573  /** assignment.<br>
2574  If the size of \a rhs is the same as the left-hand side arrays's old size, only
2575  the data are copied. Otherwise, new storage is allocated, which invalidates all
2576  objects (array views, iterators) depending on the lhs array.
2577  */
2579  {
2580  if (this != &rhs)
2581  this->copyOrReshape(rhs);
2582  return *this;
2583  }
2584 
2585  /** assignment from arbitrary MultiArrayView.<br>
2586  If the size of \a rhs is the same as the left-hand side arrays's old size, only
2587  the data are copied. Otherwise, new storage is allocated, which invalidates all
2588  objects (array views, iterators) depending on the lhs array.
2589  */
2590  template <class U, class StrideTag>
2592  {
2593  this->copyOrReshape(rhs);
2594  return *this;
2595  }
2596 
2597  /** assignment from scalar.<br>
2598  Equivalent to MultiArray::init(v).
2599  */
2601  {
2602  return this->init(v);
2603  }
2604 
2605  /** Add-assignment from arbitrary MultiArrayView. Fails with
2606  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2607  If the left array has no data (hasData() is false), this function is
2608  equivalent to a normal assignment (i.e. an empty
2609  array is interpreted as a zero-array of appropriate size).
2610  */
2611  template <class U, class StrideTag>
2613  {
2614  if(this->hasData())
2615  view_type::operator+=(rhs);
2616  else
2617  *this = rhs;
2618  return *this;
2619  }
2620 
2621  /** Subtract-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 an assignment of the negated rhs (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  this->reshape(rhs.shape());
2632  view_type::operator-=(rhs);
2633  return *this;
2634  }
2635 
2636  /** Multiply-assignment from arbitrary MultiArrayView. Fails with
2637  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2638  If the left array has no data (hasData() is false), this function is
2639  equivalent to reshape(rhs.shape()) with zero initialisation (i.e. an empty
2640  array is interpreted as a zero-array of appropriate size).
2641  */
2642  template <class U, class StrideTag>
2644  {
2645  if(this->hasData())
2646  view_type::operator*=(rhs);
2647  else
2648  this->reshape(rhs.shape());
2649  return *this;
2650  }
2651 
2652  /** Divide-assignment from arbitrary MultiArrayView. Fails with
2653  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2654  If the left array has no data (hasData() is false), this function is
2655  equivalent to reshape(rhs.shape()) with zero initialisation (i.e. an empty
2656  array is interpreted as a zero-array of appropriate size).
2657  */
2658  template <class U, class StrideTag>
2660  {
2661  if(this->hasData())
2662  view_type::operator/=(rhs);
2663  else
2664  this->reshape(rhs.shape());
2665  return *this;
2666  }
2667 
2668  /** Add-assignment of a scalar.
2669  */
2670  MultiArray &operator+= (const T &rhs)
2671  {
2672  view_type::operator+=(rhs);
2673  return *this;
2674  }
2675 
2676  /** Subtract-assignment of a scalar.
2677  */
2678  MultiArray &operator-= (const T &rhs)
2679  {
2680  view_type::operator-=(rhs);
2681  return *this;
2682  }
2683 
2684  /** Multiply-assignment of a scalar.
2685  */
2686  MultiArray &operator*= (const T &rhs)
2687  {
2688  view_type::operator*=(rhs);
2689  return *this;
2690  }
2691 
2692  /** Divide-assignment of a scalar.
2693  */
2694  MultiArray &operator/= (const T &rhs)
2695  {
2696  view_type::operator/=(rhs);
2697  return *this;
2698  }
2699  /** Assignment of an array expression. Fails with
2700  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2701  */
2702  template<class Expression>
2703  MultiArray & operator=(multi_math::MultiMathOperand<Expression> const & rhs)
2704  {
2705  multi_math::math_detail::assignOrResize(*this, rhs);
2706  return *this;
2707  }
2708 
2709  /** Add-assignment of an array expression. Fails with
2710  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2711  */
2712  template<class Expression>
2713  MultiArray & operator+=(multi_math::MultiMathOperand<Expression> const & rhs)
2714  {
2715  multi_math::math_detail::plusAssignOrResize(*this, rhs);
2716  return *this;
2717  }
2718 
2719  /** Subtract-assignment of an array expression. Fails with
2720  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2721  */
2722  template<class Expression>
2723  MultiArray & operator-=(multi_math::MultiMathOperand<Expression> const & rhs)
2724  {
2725  multi_math::math_detail::minusAssignOrResize(*this, rhs);
2726  return *this;
2727  }
2728 
2729  /** Multiply-assignment of an array expression. Fails with
2730  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2731  */
2732  template<class Expression>
2733  MultiArray & operator*=(multi_math::MultiMathOperand<Expression> const & rhs)
2734  {
2735  multi_math::math_detail::multiplyAssignOrResize(*this, rhs);
2736  return *this;
2737  }
2738 
2739  /** Divide-assignment of an array expression. Fails with
2740  <tt>PreconditionViolation</tt> exception when the shapes do not match.
2741  */
2742  template<class Expression>
2743  MultiArray & operator/=(multi_math::MultiMathOperand<Expression> const & rhs)
2744  {
2745  multi_math::math_detail::divideAssignOrResize(*this, rhs);
2746  return *this;
2747  }
2748 
2749  /** destructor
2750  */
2752  {
2753  deallocate (this->m_ptr, this->elementCount ());
2754  }
2755 
2756 
2757  /** init elements with a constant
2758  */
2759  template <class U>
2760  MultiArray & init(const U & init)
2761  {
2762  view_type::init(init);
2763  return *this;
2764  }
2765 
2766  /** Allocate new memory with the given shape and initialize with zeros.<br>
2767  <em>Note:</em> this operation invalidates all dependent objects
2768  (array views and iterators)
2769  */
2770  void reshape (const difference_type &shape)
2771  {
2772  reshape (shape, value_type());
2773  }
2774 
2775  /** Allocate new memory with the given shape and initialize it
2776  with the given value.<br>
2777  <em>Note:</em> this operation invalidates all dependent objects
2778  (array views and iterators)
2779  */
2780  void reshape (const difference_type &shape, const_reference init);
2781 
2782  /** Swap the contents with another MultiArray. This is fast,
2783  because no data are copied, but only pointers and shapes swapped.
2784  <em>Note:</em> this operation invalidates all dependent objects
2785  (array views and iterators)
2786  */
2787  void swap (MultiArray & other);
2788 
2789  // /** sequential iterator pointing to the first array element.
2790  // */
2791  // iterator begin ()
2792  // {
2793  // return vigra::detail::MultiIteratorChooser<actual_stride>::template constructIterator<iterator>((view_type *)this);
2794  // }
2795 
2796  // /** sequential iterator pointing beyond the last array element.
2797  // */
2798  // iterator end ()
2799  // {
2800  // return begin() + this->elementCount();
2801  // }
2802 
2803  // /** sequential const iterator pointing to the first array element.
2804  // */
2805  // const_iterator begin () const
2806  // {
2807  // return vigra::detail::MultiIteratorChooser<actual_stride>::template constructIterator<iterator>((view_type const *)this);
2808  // }
2809 
2810  // /** sequential const iterator pointing beyond the last array element.
2811  // */
2812  // const_iterator end () const
2813  // {
2814  // return begin() + this->elementCount();
2815  // }
2816 
2817  /** get the allocator.
2818  */
2819  allocator_type const & allocator () const
2820  {
2821  return m_alloc;
2822  }
2823 
2824  static difference_type defaultStride(difference_type const & shape)
2825  {
2826  return vigra::detail::ResolveMultiband<T>::defaultStride(shape);
2827  }
2828 };
2829 
2830 template <unsigned int N, class T, class A>
2832  allocator_type const & alloc)
2833 : view_type(difference_type(length),
2834  defaultStride(difference_type(length)),
2835  0),
2836  m_alloc(alloc)
2837 {
2838  allocate (this->m_ptr, this->elementCount (), value_type());
2839 }
2840 
2841 template <unsigned int N, class T, class A>
2843  allocator_type const & alloc)
2844 : view_type(difference_type(width, height),
2845  defaultStride(difference_type(width, height)),
2846  0),
2847  m_alloc(alloc)
2848 {
2849  allocate (this->m_ptr, this->elementCount (), value_type());
2850 }
2851 
2852 template <unsigned int N, class T, class A>
2854  allocator_type const & alloc)
2855 : view_type(shape,
2856  defaultStride(shape),
2857  0),
2858  m_alloc(alloc)
2859 {
2860  if (N == 0)
2861  {
2862  this->m_shape [0] = 1;
2863  this->m_stride [0] = 1;
2864  }
2865  allocate (this->m_ptr, this->elementCount (), value_type());
2866 }
2867 
2868 template <unsigned int N, class T, class A>
2870  allocator_type const & alloc)
2871 : view_type(shape,
2872  defaultStride(shape),
2873  0),
2874  m_alloc(alloc)
2875 {
2876  if (N == 0)
2877  {
2878  this->m_shape [0] = 1;
2879  this->m_stride [0] = 1;
2880  }
2881  allocate (this->m_ptr, this->elementCount (), init);
2882 }
2883 
2884 template <unsigned int N, class T, class A>
2886  allocator_type const & alloc)
2887 : view_type(shape,
2888  defaultStride(shape),
2889  0),
2890  m_alloc(alloc)
2891 {
2892  if (N == 0)
2893  {
2894  this->m_shape [0] = 1;
2895  this->m_stride [0] = 1;
2896  }
2897  allocate (this->m_ptr, this->elementCount (), value_type());
2898  switch(init)
2899  {
2900  case LinearSequence:
2901  linearSequence(this->begin(), this->end());
2902  break;
2903  default:
2904  vigra_precondition(false,
2905  "MultiArray(): invalid MultiArrayInitializationTag.");
2906  }
2907 }
2908 
2909 template <unsigned int N, class T, class A>
2911  allocator_type const & alloc)
2912 : view_type(shape,
2913  defaultStride(shape),
2914  0),
2915  m_alloc(alloc)
2916 {
2917  if (N == 0)
2918  {
2919  this->m_shape [0] = 1;
2920  this->m_stride [0] = 1;
2921  }
2922  allocate (this->m_ptr, this->elementCount (), init);
2923 }
2924 
2925 template <unsigned int N, class T, class A>
2926 template <class U, class StrideTag>
2928  allocator_type const & alloc)
2929 : view_type(rhs.shape(),
2930  defaultStride(rhs.shape()),
2931  0),
2932  m_alloc (alloc)
2933 {
2934  allocate (this->m_ptr, rhs);
2935 }
2936 
2937 template <unsigned int N, class T, class A>
2938 template <class U, class StrideTag>
2939 void
2941 {
2942  if (this->shape() == rhs.shape())
2943  this->copy(rhs);
2944  else
2945  {
2946  MultiArray t(rhs);
2947  this->swap(t);
2948  }
2949 }
2950 
2951 template <unsigned int N, class T, class A>
2953  const_reference initial)
2954 {
2955  if (N == 0)
2956  {
2957  return;
2958  }
2959  else if(new_shape == this->shape())
2960  {
2961  this->init(initial);
2962  }
2963  else
2964  {
2965  difference_type new_stride = defaultStride(new_shape);
2966  difference_type_1 new_size = prod(new_shape);
2967  pointer new_ptr = pointer();
2968  allocate (new_ptr, new_size, initial);
2969  deallocate (this->m_ptr, this->elementCount ());
2970  this->m_ptr = new_ptr;
2971  this->m_shape = new_shape;
2972  this->m_stride = new_stride;
2973  }
2974 }
2975 
2976 
2977 template <unsigned int N, class T, class A>
2978 inline void
2980 {
2981  if (this == &other)
2982  return;
2983  std::swap(this->m_shape, other.m_shape);
2984  std::swap(this->m_stride, other.m_stride);
2985  std::swap(this->m_ptr, other.m_ptr);
2986  std::swap(this->m_alloc, other.m_alloc);
2987 }
2988 
2989 template <unsigned int N, class T, class A>
2992 {
2993  if(s == 0)
2994  {
2995  ptr = 0;
2996  return;
2997  }
2998  ptr = m_alloc.allocate ((typename A::size_type)s);
2999  difference_type_1 i = 0;
3000  try {
3001  for (; i < s; ++i)
3002  m_alloc.construct (ptr + i, init);
3003  }
3004  catch (...) {
3005  for (difference_type_1 j = 0; j < i; ++j)
3006  m_alloc.destroy (ptr + j);
3007  m_alloc.deallocate (ptr, (typename A::size_type)s);
3008  throw;
3009  }
3010 }
3011 
3012 template <unsigned int N, class T, class A>
3013 template <class U>
3015  U const * init)
3016 {
3017  if(s == 0)
3018  {
3019  ptr = 0;
3020  return;
3021  }
3022  ptr = m_alloc.allocate ((typename A::size_type)s);
3023  difference_type_1 i = 0;
3024  try {
3025  for (; i < s; ++i, ++init)
3026  m_alloc.construct (ptr + i, *init);
3027  }
3028  catch (...) {
3029  for (difference_type_1 j = 0; j < i; ++j)
3030  m_alloc.destroy (ptr + j);
3031  m_alloc.deallocate (ptr, (typename A::size_type)s);
3032  throw;
3033  }
3034 }
3035 
3036 template <unsigned int N, class T, class A>
3037 template <class U, class StrideTag>
3039 {
3040  difference_type_1 s = init.elementCount();
3041  if(s == 0)
3042  {
3043  ptr = 0;
3044  return;
3045  }
3046  ptr = m_alloc.allocate ((typename A::size_type)s);
3047  pointer p = ptr;
3048  try {
3049  detail::uninitializedCopyMultiArrayData(init.traverser_begin(), init.shape(),
3050  p, m_alloc, MetaInt<actual_dimension-1>());
3051  }
3052  catch (...) {
3053  for (pointer pp = ptr; pp < p; ++pp)
3054  m_alloc.destroy (pp);
3055  m_alloc.deallocate (ptr, (typename A::size_type)s);
3056  throw;
3057  }
3058 }
3059 
3060 template <unsigned int N, class T, class A>
3062 {
3063  if (ptr == 0)
3064  return;
3065  for (difference_type_1 i = 0; i < s; ++i)
3066  m_alloc.destroy (ptr + i);
3067  m_alloc.deallocate (ptr, (typename A::size_type)s);
3068  ptr = 0;
3069 }
3070 
3071 /********************************************************/
3072 /* */
3073 /* argument object factories */
3074 /* */
3075 /********************************************************/
3076 
3077 template <unsigned int N, class T, class StrideTag>
3078 inline triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3081 srcMultiArrayRange( MultiArrayView<N,T,StrideTag> const & array )
3082 {
3083  return triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3086  ( array.traverser_begin(),
3087  array.shape(),
3089 }
3090 
3091 template <unsigned int N, class T, class StrideTag, class Accessor>
3092 inline triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3093  typename MultiArrayView<N,T,StrideTag>::difference_type,
3094  Accessor >
3095 srcMultiArrayRange( MultiArrayView<N,T,StrideTag> const & array, Accessor a )
3096 {
3097  return triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3098  typename MultiArrayView<N,T,StrideTag>::difference_type,
3099  Accessor >
3100  ( array.traverser_begin(),
3101  array.shape(),
3102  a);
3103 }
3104 
3105 template <unsigned int N, class T, class StrideTag>
3106 inline pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3107  typename AccessorTraits<T>::default_const_accessor >
3108 srcMultiArray( MultiArrayView<N,T,StrideTag> const & array )
3109 {
3110  return pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3111  typename AccessorTraits<T>::default_const_accessor >
3112  ( array.traverser_begin(),
3113  typename AccessorTraits<T>::default_const_accessor() );
3114 }
3115 
3116 template <unsigned int N, class T, class StrideTag, class Accessor>
3117 inline pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3118  Accessor >
3119 srcMultiArray( MultiArrayView<N,T,StrideTag> const & array, Accessor a )
3120 {
3121  return pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
3122  Accessor >
3123  ( array.traverser_begin(), a );
3124 }
3125 
3126 template <unsigned int N, class T, class StrideTag>
3127 inline triple<typename MultiArrayView<N,T,StrideTag>::traverser,
3128  typename MultiArrayView<N,T,StrideTag>::difference_type,
3129  typename AccessorTraits<T>::default_accessor >
3130 destMultiArrayRange( MultiArrayView<N,T,StrideTag> & array )
3131 {
3132  return triple<typename MultiArrayView<N,T,StrideTag>::traverser,
3133  typename MultiArrayView<N,T,StrideTag>::difference_type,
3134  typename AccessorTraits<T>::default_accessor >
3135  ( array.traverser_begin(),
3136  array.shape(),
3137  typename AccessorTraits<T>::default_accessor() );
3138 }
3139 
3140 template <unsigned int N, class T, class StrideTag, class Accessor>
3141 inline triple<typename MultiArrayView<N,T,StrideTag>::traverser,
3142  typename MultiArrayView<N,T,StrideTag>::difference_type,
3143  Accessor >
3144 destMultiArrayRange( MultiArrayView<N,T,StrideTag> & array, Accessor a )
3145 {
3146  return triple<typename MultiArrayView<N,T,StrideTag>::traverser,
3147  typename MultiArrayView<N,T,StrideTag>::difference_type,
3148  Accessor >
3149  ( array.traverser_begin(),
3150  array.shape(),
3151  a );
3152 }
3153 
3154 template <unsigned int N, class T, class StrideTag>
3155 inline pair<typename MultiArrayView<N,T,StrideTag>::traverser,
3156  typename AccessorTraits<T>::default_accessor >
3157 destMultiArray( MultiArrayView<N,T,StrideTag> & array )
3158 {
3159  return pair<typename MultiArrayView<N,T,StrideTag>::traverser,
3160  typename AccessorTraits<T>::default_accessor >
3161  ( array.traverser_begin(),
3162  typename AccessorTraits<T>::default_accessor() );
3163 }
3164 
3165 template <unsigned int N, class T, class StrideTag, class Accessor>
3166 inline pair<typename MultiArrayView<N,T,StrideTag>::traverser,
3167  Accessor >
3168 destMultiArray( MultiArrayView<N,T,StrideTag> & array, Accessor a )
3169 {
3170  return pair<typename MultiArrayView<N,T,StrideTag>::traverser,
3171  Accessor >
3172  ( array.traverser_begin(), a );
3173 }
3174 
3175 /********************************************************************/
3176 
3177 template <class PixelType, class Accessor>
3178 inline triple<ConstStridedImageIterator<PixelType>,
3179  ConstStridedImageIterator<PixelType>, Accessor>
3180 srcImageRange(const MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
3181 {
3182  ConstStridedImageIterator<PixelType>
3183  ul(img.data(), 1, img.stride(0), img.stride(1));
3184  return triple<ConstStridedImageIterator<PixelType>,
3185  ConstStridedImageIterator<PixelType>,
3186  Accessor>(
3187  ul, ul + Size2D(img.shape(0), img.shape(1)), a);
3188 }
3189 
3190 template <class PixelType, class Accessor>
3191 inline pair<ConstStridedImageIterator<PixelType>, Accessor>
3192 srcImage(const MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
3193 {
3194  ConstStridedImageIterator<PixelType>
3195  ul(img.data(), 1, img.stride(0), img.stride(1));
3196  return pair<ConstStridedImageIterator<PixelType>, Accessor>
3197  (ul, a);
3198 }
3199 
3200 template <class PixelType, class Accessor>
3201 inline triple<StridedImageIterator<PixelType>,
3202  StridedImageIterator<PixelType>, Accessor>
3203 destImageRange(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
3204 {
3205  StridedImageIterator<PixelType>
3206  ul(img.data(), 1, img.stride(0), img.stride(1));
3207  return triple<StridedImageIterator<PixelType>,
3208  StridedImageIterator<PixelType>,
3209  Accessor>(
3210  ul, ul + Size2D(img.shape(0), img.shape(1)), a);
3211 }
3212 
3213 template <class PixelType, class Accessor>
3214 inline pair<StridedImageIterator<PixelType>, Accessor>
3215 destImage(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
3216 {
3217  StridedImageIterator<PixelType>
3218  ul(img.data(), 1, img.stride(0), img.stride(1));
3219  return pair<StridedImageIterator<PixelType>, Accessor>
3220  (ul, a);
3221 }
3222 
3223 template <class PixelType, class Accessor>
3224 inline pair<StridedImageIterator<PixelType>, Accessor>
3225 maskImage(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
3226 {
3227  StridedImageIterator<PixelType>
3228  ul(img.data(), 1, img.stride(0), img.stride(1));
3229  return pair<StridedImageIterator<PixelType>, Accessor>
3230  (ul, a);
3231 }
3232 
3233 // -------------------------------------------------------------------
3234 
3235 template <class PixelType>
3236 inline triple<ConstStridedImageIterator<PixelType>,
3237  ConstStridedImageIterator<PixelType>,
3238  typename AccessorTraits<PixelType>::default_const_accessor>
3239 srcImageRange(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
3240 {
3241  ConstStridedImageIterator<PixelType>
3242  ul(img.data(), 1, img.stride(0), img.stride(1));
3243  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
3244  return triple<ConstStridedImageIterator<PixelType>,
3245  ConstStridedImageIterator<PixelType>,
3246  Accessor>
3247  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
3248 }
3249 
3250 template <class PixelType>
3251 inline triple<ConstImageIterator<PixelType>,
3252  ConstImageIterator<PixelType>,
3253  typename AccessorTraits<PixelType>::default_const_accessor>
3254 srcImageRange(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
3255 {
3256  ConstImageIterator<PixelType>
3257  ul(img.data(), img.stride(1));
3258  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
3259  return triple<ConstImageIterator<PixelType>,
3260  ConstImageIterator<PixelType>,
3261  Accessor>
3262  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
3263 }
3264 
3265 template <class PixelType>
3266 inline pair< ConstStridedImageIterator<PixelType>,
3267  typename AccessorTraits<PixelType>::default_const_accessor>
3268 srcImage(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
3269 {
3270  ConstStridedImageIterator<PixelType>
3271  ul(img.data(), 1, img.stride(0), img.stride(1));
3272  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
3273  return pair<ConstStridedImageIterator<PixelType>,
3274  Accessor>
3275  (ul, Accessor());
3276 }
3277 
3278 template <class PixelType>
3279 inline pair< ConstImageIterator<PixelType>,
3280  typename AccessorTraits<PixelType>::default_const_accessor>
3281 srcImage(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
3282 {
3283  ConstImageIterator<PixelType>
3284  ul(img.data(), img.stride(1));
3285  typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
3286  return pair<ConstImageIterator<PixelType>,
3287  Accessor>
3288  (ul, Accessor());
3289 }
3290 
3291 template <class PixelType>
3292 inline triple< StridedImageIterator<PixelType>,
3293  StridedImageIterator<PixelType>,
3294  typename AccessorTraits<PixelType>::default_accessor>
3295 destImageRange(MultiArrayView<2, PixelType, StridedArrayTag> & img)
3296 {
3297  StridedImageIterator<PixelType>
3298  ul(img.data(), 1, img.stride(0), img.stride(1));
3299  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3300  return triple<StridedImageIterator<PixelType>,
3301  StridedImageIterator<PixelType>,
3302  Accessor>
3303  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
3304 }
3305 
3306 template <class PixelType>
3307 inline triple< ImageIterator<PixelType>,
3308  ImageIterator<PixelType>,
3309  typename AccessorTraits<PixelType>::default_accessor>
3310 destImageRange(MultiArrayView<2, PixelType, UnstridedArrayTag> & img)
3311 {
3312  ImageIterator<PixelType>
3313  ul(img.data(), img.stride(1));
3314  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3315  return triple<ImageIterator<PixelType>,
3316  ImageIterator<PixelType>,
3317  Accessor>
3318  (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
3319 }
3320 
3321 template <class PixelType>
3322 inline pair< StridedImageIterator<PixelType>,
3323  typename AccessorTraits<PixelType>::default_accessor>
3324 destImage(MultiArrayView<2, PixelType, StridedArrayTag> & img)
3325 {
3326  StridedImageIterator<PixelType>
3327  ul(img.data(), 1, img.stride(0), img.stride(1));
3328  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3329  return pair<StridedImageIterator<PixelType>, Accessor>
3330  (ul, Accessor());
3331 }
3332 
3333 template <class PixelType>
3334 inline pair< ImageIterator<PixelType>,
3335  typename AccessorTraits<PixelType>::default_accessor>
3336 destImage(MultiArrayView<2, PixelType, UnstridedArrayTag> & img)
3337 {
3338  ImageIterator<PixelType> ul(img.data(), img.stride(1));
3339  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3340  return pair<ImageIterator<PixelType>, Accessor>(ul, Accessor());
3341 }
3342 
3343 template <class PixelType>
3344 inline pair< ConstStridedImageIterator<PixelType>,
3345  typename AccessorTraits<PixelType>::default_accessor>
3346 maskImage(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
3347 {
3348  ConstStridedImageIterator<PixelType>
3349  ul(img.data(), 1, img.stride(0), img.stride(1));
3350  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3351  return pair<ConstStridedImageIterator<PixelType>, Accessor>
3352  (ul, Accessor());
3353 }
3354 
3355 template <class PixelType>
3356 inline pair< ConstImageIterator<PixelType>,
3357  typename AccessorTraits<PixelType>::default_accessor>
3358 maskImage(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
3359 {
3360  ConstImageIterator<PixelType>
3361  ul(img.data(), img.stride(1));
3362  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
3363  return pair<ConstImageIterator<PixelType>, Accessor>
3364  (ul, Accessor());
3365 }
3366 
3367 /********************************************************/
3368 /* */
3369 /* makeBasicImageView */
3370 /* */
3371 /********************************************************/
3372 
3373 /** \addtogroup MultiArrayToImage Create BasicImageView from MultiArrayViews
3374 
3375  Some convenience functions for wrapping a \ref vigra::MultiArrayView's
3376  data in a \ref vigra::BasicImageView.
3377 */
3378 //@{
3379 /** Create a \ref vigra::BasicImageView from an unstrided 2-dimensional
3380  \ref vigra::MultiArrayView.
3381 
3382  The \ref vigra::BasicImageView will have the same <tt>value_type </tt>
3383  as the original \ref vigra::MultiArrayView.
3384 */
3385 template <class T, class Stride>
3386 BasicImageView <T>
3388 {
3389  vigra_precondition(array.isUnstrided(),
3390  "makeBasicImageView(array): array must be unstrided (i.e. array.isUnstrided() == true).");
3391  return BasicImageView <T> (array.data (), array.shape (0),
3392  array.shape (1), array.stride(1));
3393 }
3394 
3395 /** Create a \ref vigra::BasicImageView from a 3-dimensional
3396  \ref vigra::MultiArray.
3397 
3398  This wrapper flattens the two innermost dimensions of the array
3399  into single rows of the resulting image.
3400  The \ref vigra::BasicImageView will have the same <tt>value_type </tt>
3401  as the original \ref vigra::MultiArray.
3402 */
3403 template <class T>
3404 BasicImageView <T>
3406 {
3407  vigra_precondition(array.stride(1) == array.shape(0),
3408  "makeBasicImageView(): cannot join strided dimensions");
3409  return BasicImageView <T> (array.data (),
3410  array.shape (0)*array.shape (1), array.shape (2), array.stride(2));
3411 }
3412 
3413 /** Create a \ref vigra::BasicImageView from a 3-dimensional
3414  \ref vigra::MultiArray.
3415 
3416  This wrapper only works if <tt>T</tt> is a scalar type and the
3417  array's innermost dimension has size 3. It then re-interprets
3418  the data array as a 2-dimensional array with value_type
3419  <tt>RGBValue<T></tt>.
3420 */
3421 template <class T, class Stride>
3422 BasicImageView <RGBValue<T> >
3424 {
3425  vigra_precondition(array.shape (0) == 3,
3426  "makeRGBImageView(): array.shape(0) must be 3.");
3427  vigra_precondition(array.isUnstrided(),
3428  "makeRGBImageView(array): array must be unstrided (i.e. array.isUnstrided() == true).");
3429  return BasicImageView <RGBValue<T> > (
3430  reinterpret_cast <RGBValue <T> *> (array.data ()),
3431  array.shape (1), array.shape (2));
3432 }
3433 
3434 //@}
3435 
3436 } // namespace vigra
3437 
3438 #undef VIGRA_ASSERT_INSIDE
3439 
3440 #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.10.0 (Thu Jan 8 2015)