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

vigra/multi_array.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*     Copyright 2003-2008 by Gunnar Kedenburg and Ullrich Koethe       */
00004 /*                                                                      */
00005 /*    This file is part of the VIGRA computer vision library.           */
00006 /*    The VIGRA Website is                                              */
00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
00008 /*    Please direct questions, bug reports, and contributions to        */
00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00010 /*        vigra@informatik.uni-hamburg.de                               */
00011 /*                                                                      */
00012 /*    Permission is hereby granted, free of charge, to any person       */
00013 /*    obtaining a copy of this software and associated documentation    */
00014 /*    files (the "Software"), to deal in the Software without           */
00015 /*    restriction, including without limitation the rights to use,      */
00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00017 /*    sell copies of the Software, and to permit persons to whom the    */
00018 /*    Software is furnished to do so, subject to the following          */
00019 /*    conditions:                                                       */
00020 /*                                                                      */
00021 /*    The above copyright notice and this permission notice shall be    */
00022 /*    included in all copies or substantial portions of the             */
00023 /*    Software.                                                         */
00024 /*                                                                      */
00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00033 /*                                                                      */
00034 /************************************************************************/
00035 
00036 
00037 #ifndef VIGRA_MULTI_ARRAY_HXX
00038 #define VIGRA_MULTI_ARRAY_HXX
00039 
00040 #include <memory>
00041 #include <algorithm>
00042 #include "accessor.hxx"
00043 #include "tinyvector.hxx"
00044 #include "rgbvalue.hxx"
00045 #include "basicimageview.hxx"
00046 #include "imageiterator.hxx"
00047 #include "numerictraits.hxx"
00048 #include "multi_iterator.hxx"
00049 #include "metaprogramming.hxx"
00050 #include "mathutil.hxx"
00051 
00052 // Bounds checking Macro used if VIGRA_CHECK_BOUNDS is defined.
00053 #ifdef VIGRA_CHECK_BOUNDS
00054 #define VIGRA_ASSERT_INSIDE(diff) \
00055   vigra_precondition(this->isInside(diff), "Index out of bounds")
00056 #else
00057 #define VIGRA_ASSERT_INSIDE(diff)
00058 #endif
00059 
00060 namespace vigra
00061 {
00062 
00063 namespace detail
00064 {
00065 /********************************************************/
00066 /*                                                      */
00067 /*                    defaultStride                     */
00068 /*                                                      */
00069 /********************************************************/
00070 
00071 /* generates the stride for a gapless shape.
00072 
00073     Namespace: vigra::detail
00074 */
00075 template <unsigned int N>
00076 inline TinyVector <MultiArrayIndex, N>
00077 defaultStride(const TinyVector <MultiArrayIndex, N> &shape)
00078 {
00079     TinyVector <MultiArrayIndex, N> ret;
00080     ret [0] = 1;
00081     for (int i = 1; i < (int)N; ++i)
00082         ret [i] = ret [i-1] * shape [i-1];
00083     return ret;
00084 }
00085 
00086 /********************************************************/
00087 /*                                                      */
00088 /*                 ScanOrderToOffset                    */
00089 /*                                                      */
00090 /********************************************************/
00091 
00092 /* transforms an index in scan order sense to a pointer offset in a possibly
00093    strided, multi-dimensional array.
00094 
00095     Namespace: vigra::detail
00096 */
00097 
00098 template <int K>
00099 struct ScanOrderToOffset
00100 {
00101     template <int N>
00102     static MultiArrayIndex
00103     exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> &shape,
00104          const TinyVector <MultiArrayIndex, N> & stride)
00105     {
00106         return stride[N-K] * (d % shape[N-K]) +
00107                ScanOrderToOffset<K-1>::exec(d / shape[N-K], shape, stride);
00108     }
00109 };
00110 
00111 template <>
00112 struct ScanOrderToOffset<1>
00113 {
00114     template <int N>
00115     static MultiArrayIndex
00116     exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> & /*shape*/,
00117          const TinyVector <MultiArrayIndex, N> & stride)
00118     {
00119         return stride[N-1] * d;
00120     }
00121 };
00122 
00123 template <int K>
00124 struct ScanOrderToCoordinate
00125 {
00126     template <int N>
00127     static void
00128     exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> &shape,
00129          TinyVector <MultiArrayIndex, N> & result)
00130     {
00131         result[N-K] = (d % shape[N-K]);
00132         ScanOrderToCoordinate<K-1>::exec(d / shape[N-K], shape, result);
00133     }
00134 };
00135 
00136 template <>
00137 struct ScanOrderToCoordinate<1>
00138 {
00139     template <int N>
00140     static void
00141     exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> & /*shape*/,
00142          TinyVector <MultiArrayIndex, N> & result)
00143     {
00144         result[N-1] = d;
00145     }
00146 };
00147 
00148 
00149 template <class C>
00150 struct CoordinatesToOffest
00151 {
00152     template <int N>
00153     static MultiArrayIndex
00154     exec(const TinyVector <MultiArrayIndex, N> & stride, MultiArrayIndex x)
00155     {
00156         return stride[0] * x;
00157     }
00158     template <int N>
00159     static MultiArrayIndex
00160     exec(const TinyVector <MultiArrayIndex, N> & stride, MultiArrayIndex x, MultiArrayIndex y)
00161     {
00162         return stride[0] * x + stride[1] * y;
00163     }
00164 };
00165 
00166 template <>
00167 struct CoordinatesToOffest<UnstridedArrayTag>
00168 {
00169     template <int N>
00170     static MultiArrayIndex
00171     exec(const TinyVector <MultiArrayIndex, N> & /*stride*/, MultiArrayIndex x)
00172     {
00173         return x;
00174     }
00175     template <int N>
00176     static MultiArrayIndex
00177     exec(const TinyVector <MultiArrayIndex, N> & stride, MultiArrayIndex x, MultiArrayIndex y)
00178     {
00179         return x + stride[1] * y;
00180     }
00181 };
00182 
00183 /********************************************************/
00184 /*                                                      */
00185 /*                     MaybeStrided                     */
00186 /*                                                      */
00187 /********************************************************/
00188 
00189 /* metatag implementing a test for marking MultiArrays that were
00190     indexed at the zero'th dimension as strided, and all others as
00191     unstrided.
00192 
00193 <b>\#include</b>
00194 <vigra/multi_array.hxx>
00195 
00196 Namespace: vigra::detail
00197 */
00198 template <class StrideTag, unsigned int N>
00199 struct MaybeStrided
00200 {
00201     typedef StrideTag type;
00202 };
00203 
00204 template <class StrideTag>
00205 struct MaybeStrided <StrideTag, 0>
00206 {
00207     typedef StridedArrayTag type;
00208 };
00209 
00210 /********************************************************/
00211 /*                                                      */
00212 /*                MultiIteratorChooser                  */
00213 /*                                                      */
00214 /********************************************************/
00215 
00216 /* metatag implementing a test (by pattern matching) for marking
00217     MultiArrays that were indexed at the zero'th dimension as strided.
00218 
00219 <b>\#include</b>
00220 <vigra/multi_array.hxx>
00221 
00222 Namespace: vigra::detail
00223 */
00224 template <class O>
00225 struct MultiIteratorChooser
00226 {
00227     struct Nil {};
00228 
00229     template <unsigned int N, class T, class REFERENCE, class POINTER>
00230     struct Traverser
00231     {
00232         typedef Nil type;
00233     };
00234 };
00235 
00236 /********************************************************/
00237 /*                                                      */
00238 /*       MultiIteratorChooser <StridedArrayTag>         */
00239 /*                                                      */
00240 /********************************************************/
00241 
00242 /* specialization of the MultiIteratorChooser for strided arrays.
00243 
00244 <b>\#include</b>
00245 <vigra/multi_array.hxx>
00246 
00247 Namespace: vigra::detail
00248 */
00249 template <>
00250 struct MultiIteratorChooser <StridedArrayTag>
00251 {
00252     template <unsigned int N, class T, class REFERENCE, class POINTER>
00253     struct Traverser
00254     {
00255         typedef StridedMultiIterator <N, T, REFERENCE, POINTER> type;
00256     };
00257 };
00258 
00259 /********************************************************/
00260 /*                                                      */
00261 /*      MultiIteratorChooser <UnstridedArrayTag>        */
00262 /*                                                      */
00263 /********************************************************/
00264 
00265 /* specialization of the MultiIteratorChooser for unstrided arrays.
00266 
00267 <b>\#include</b>
00268 <vigra/multi_array.hxx>
00269 
00270 Namespace: vigra::detail
00271 */
00272 template <>
00273 struct MultiIteratorChooser <UnstridedArrayTag>
00274 {
00275     template <unsigned int N, class T, class REFERENCE, class POINTER>
00276     struct Traverser
00277     {
00278         typedef MultiIterator <N, T, REFERENCE, POINTER> type;
00279     };
00280 };
00281 
00282 /********************************************************/
00283 /*                                                      */
00284 /*                   helper functions                   */
00285 /*                                                      */
00286 /********************************************************/
00287 
00288 template <class DestIterator, class Shape, class T>
00289 inline void
00290 initMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<0>)
00291 {
00292     DestIterator dend = d + shape[0];
00293     for(; d < dend; ++d)
00294     {
00295         *d = init;
00296     }
00297 }
00298 
00299 template <class DestIterator, class Shape, class T, int N>
00300 void
00301 initMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<N>)
00302 {
00303     DestIterator dend = d + shape[N];
00304     for(; d < dend; ++d)
00305     {
00306         initMultiArrayData(d.begin(), shape, init, MetaInt<N-1>());
00307     }
00308 }
00309 
00310 // FIXME: the explicit overload for MultiIterator<1, UInt8, ... > works around a compiler crash in VisualStudio 2010
00311 #define VIGRA_COPY_MULTI_ARRAY_DATA(name, op) \
00312 template <class SrcIterator, class Shape, class DestIterator> \
00313 inline void \
00314 name##MultiArrayData(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>) \
00315 {     \
00316     SrcIterator send = s + shape[0]; \
00317     for(; s < send; ++s, ++d) \
00318     { \
00319         *d op detail::RequiresExplicitCast<typename DestIterator::value_type>::cast(*s); \
00320     } \
00321 } \
00322  \
00323 template <class Ref, class Ptr, class Shape, class DestIterator> \
00324 inline void \
00325 name##MultiArrayData(MultiIterator<1, UInt8, Ref, Ptr> si, Shape const & shape, DestIterator d, MetaInt<0>) \
00326 { \
00327     Ptr s = &(*si), send = s + shape[0]; \
00328     for(; s < send; ++s, ++d) \
00329     { \
00330         *d op detail::RequiresExplicitCast<typename DestIterator::value_type>::cast(*s); \
00331     } \
00332 } \
00333 \
00334 template <class SrcIterator, class Shape, class DestIterator, int N> \
00335 void \
00336 name##MultiArrayData(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>) \
00337 { \
00338     SrcIterator send = s + shape[N]; \
00339     for(; s < send; ++s, ++d) \
00340     { \
00341         name##MultiArrayData(s.begin(), shape, d.begin(), MetaInt<N-1>()); \
00342     } \
00343 } \
00344 \
00345 template <class DestIterator, class Shape, class T> \
00346 inline void \
00347 name##ScalarMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<0>) \
00348 {     \
00349     DestIterator dend = d + shape[0]; \
00350     for(; d < dend; ++d) \
00351     { \
00352         *d op detail::RequiresExplicitCast<typename DestIterator::value_type>::cast(init); \
00353     } \
00354 } \
00355  \
00356 template <class DestIterator, class Shape, class T, int N> \
00357 void \
00358 name##ScalarMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<N>) \
00359 {     \
00360     DestIterator dend = d + shape[N]; \
00361     for(; d < dend; ++d) \
00362     { \
00363         name##ScalarMultiArrayData(d.begin(), shape, init, MetaInt<N-1>()); \
00364     } \
00365 }
00366 
00367 VIGRA_COPY_MULTI_ARRAY_DATA(copy, =)
00368 VIGRA_COPY_MULTI_ARRAY_DATA(copyAdd, +=)
00369 VIGRA_COPY_MULTI_ARRAY_DATA(copySub, -=)
00370 VIGRA_COPY_MULTI_ARRAY_DATA(copyMul, *=)
00371 VIGRA_COPY_MULTI_ARRAY_DATA(copyDiv, /=)
00372 
00373 #undef VIGRA_COPY_MULTI_ARRAY_DATA
00374 
00375 template <class SrcIterator, class Shape, class T, class ALLOC>
00376 inline void
00377 uninitializedCopyMultiArrayData(SrcIterator s, Shape const & shape, T * & d, ALLOC & a, MetaInt<0>)
00378 {
00379     SrcIterator send = s + shape[0];
00380     for(; s < send; ++s, ++d)
00381     {
00382         a.construct(d, static_cast<T const &>(*s));
00383     }
00384 }
00385 
00386 // FIXME: this overload works around a compiler crash in VisualStudio 2010
00387 template <class Ref, class Ptr, class Shape, class T, class ALLOC>
00388 inline void
00389 uninitializedCopyMultiArrayData(MultiIterator<1, UInt8, Ref, Ptr> si, Shape const & shape, T * & d, ALLOC & a, MetaInt<0>)
00390 {
00391     Ptr s = &(*si), send = s + shape[0];
00392     for(; s < send; ++s, ++d)
00393     {
00394         a.construct(d, static_cast<T const &>(*s));
00395     }
00396 }
00397 
00398 template <class SrcIterator, class Shape, class T, class ALLOC, int N>
00399 void
00400 uninitializedCopyMultiArrayData(SrcIterator s, Shape const & shape, T * & d, ALLOC & a, MetaInt<N>)
00401 {
00402     SrcIterator send = s + shape[N];
00403     for(; s < send; ++s)
00404     {
00405         uninitializedCopyMultiArrayData(s.begin(), shape, d, a, MetaInt<N-1>());
00406     }
00407 }
00408 
00409 template <class SrcIterator, class Shape, class T, class Functor>
00410 inline void
00411 reduceOverMultiArray(SrcIterator s, Shape const & shape, T & result, Functor const & f, MetaInt<0>)
00412 {
00413     SrcIterator send = s + shape[0];
00414     for(; s < send; ++s)
00415     {
00416         f(result, *s);
00417     }
00418 }
00419 
00420 template <class SrcIterator, class Shape, class T, class Functor, int N>
00421 void
00422 reduceOverMultiArray(SrcIterator s, Shape const & shape, T & result, Functor const & f, MetaInt<N>)
00423 {
00424     SrcIterator send = s + shape[N];
00425     for(; s < send; ++s)
00426     {
00427         reduceOverMultiArray(s.begin(), shape, result, f, MetaInt<N-1>());
00428     }
00429 }
00430 
00431 struct MaxNormReduceFunctor
00432 {
00433     template <class T, class U>
00434     void operator()(T & result, U const & u) const
00435     {
00436         T v = norm(u);
00437         if(result < v)
00438             result = v;
00439     }
00440 };
00441 
00442 struct L1NormReduceFunctor
00443 {
00444     template <class T, class U>
00445     void operator()(T & result, U const & u) const
00446     {
00447         result += norm(u);
00448     }
00449 };
00450 
00451 struct SquaredL2NormReduceFunctor
00452 {
00453     template <class T, class U>
00454     void operator()(T & result, U const & u) const
00455     {
00456         result += squaredNorm(u);
00457     }
00458 };
00459 
00460 template <class T>
00461 struct WeightedL2NormReduceFunctor
00462 {
00463     T scale;
00464 
00465     WeightedL2NormReduceFunctor(T s)
00466     : scale(s)
00467     {}
00468 
00469     template <class U>
00470     void operator()(T & result, U const & u) const
00471     {
00472         result += squaredNorm(u * scale);
00473     }
00474 };
00475 
00476 struct SumReduceFunctor
00477 {
00478     template <class T, class U>
00479     void operator()(T & result, U const & u) const
00480     {
00481         result += u;
00482     }
00483 };
00484 
00485 struct ProdReduceFunctor
00486 {
00487     template <class T, class U>
00488     void operator()(T & result, U const & u) const
00489     {
00490         result *= u;
00491     }
00492 };
00493 
00494 struct MinmaxReduceFunctor
00495 {
00496     template <class T, class U>
00497     void operator()(T & result, U const & u) const
00498     {
00499         if(u < result.first)
00500             result.first = u;
00501         if(result.second < u)
00502             result.second = u;
00503     }
00504 };
00505 
00506 struct MeanVarianceReduceFunctor
00507 {
00508     template <class T, class U>
00509     void operator()(T & result, U const & u) const
00510     {
00511         ++result.first;
00512         typename T::second_type t1 = u - result.second;
00513         typename T::second_type t2 = t1 / result.first;
00514         result.second += t2;
00515         result.third += (result.first-1.0)*t1*t2;
00516     }
00517 };
00518 
00519 struct AllTrueReduceFunctor
00520 {
00521     template <class T, class U>
00522     void operator()(T & result, U const & u) const
00523     {
00524         result = result && (u != NumericTraits<U>::zero());
00525     }
00526 };
00527 
00528 struct AnyTrueReduceFunctor
00529 {
00530     template <class T, class U>
00531     void operator()(T & result, U const & u) const
00532     {
00533         result = result || (u != NumericTraits<U>::zero());
00534     }
00535 };
00536 
00537 template <class SrcIterator, class Shape, class DestIterator>
00538 inline bool
00539 equalityOfMultiArrays(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>)
00540 {
00541     SrcIterator send = s + shape[0];
00542     for(; s < send; ++s, ++d)
00543     {
00544         if(!(*s == *d))
00545             return false;
00546     }
00547     return true;
00548 }
00549 
00550 template <class SrcIterator, class Shape, class DestIterator, int N>
00551 bool
00552 equalityOfMultiArrays(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>)
00553 {
00554     SrcIterator send = s + shape[N];
00555     for(; s < send; ++s, ++d)
00556     {
00557         if(!equalityOfMultiArrays(s.begin(), shape, d.begin(), MetaInt<N-1>()))
00558             return false;
00559     }
00560     return true;
00561 }
00562 
00563 
00564 template <class SrcIterator, class Shape, class DestIterator>
00565 inline void
00566 swapDataImpl(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>)
00567 {
00568     SrcIterator send = s + shape[0];
00569     for(; s < send; ++s, ++d)
00570         std::swap(*s, *d);
00571 }
00572 
00573 template <class SrcIterator, class Shape, class DestIterator, int N>
00574 void
00575 swapDataImpl(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>)
00576 {
00577     SrcIterator send = s + shape[N];
00578     for(; s < send; ++s, ++d)
00579         swapDataImpl(s.begin(), shape, d.begin(), MetaInt<N-1>());
00580 }
00581 
00582 } // namespace detail
00583 
00584 /********************************************************/
00585 /*                                                      */
00586 /*                     MultiArrayView                   */
00587 /*                                                      */
00588 /********************************************************/
00589 
00590 // forward declarations
00591 
00592 template <unsigned int N, class T, class C = UnstridedArrayTag>
00593 class MultiArrayView;
00594 template <unsigned int N, class T, class A = std::allocator<T> >
00595 class MultiArray;
00596 
00597 namespace multi_math {
00598 
00599 template <class T>
00600 struct MultiMathOperand;
00601 
00602 namespace detail {
00603 
00604 template <unsigned int N, class T, class C, class E>
00605 void assign(MultiArrayView<N, T, C>, MultiMathOperand<E> const &);
00606 
00607 template <unsigned int N, class T, class C, class E>
00608 void plusAssign(MultiArrayView<N, T, C>, MultiMathOperand<E> const &);
00609 
00610 template <unsigned int N, class T, class C, class E>
00611 void minusAssign(MultiArrayView<N, T, C>, MultiMathOperand<E> const &);
00612 
00613 template <unsigned int N, class T, class C, class E>
00614 void multiplyAssign(MultiArrayView<N, T, C>, MultiMathOperand<E> const &);
00615 
00616 template <unsigned int N, class T, class C, class E>
00617 void divideAssign(MultiArrayView<N, T, C>, MultiMathOperand<E> const &);
00618 
00619 template <unsigned int N, class T, class A, class E>
00620 void assignOrResize(MultiArray<N, T, A> &, MultiMathOperand<E> const &);
00621 
00622 template <unsigned int N, class T, class A, class E>
00623 void plusAssignOrResize(MultiArray<N, T, A> &, MultiMathOperand<E> const &);
00624 
00625 template <unsigned int N, class T, class A, class E>
00626 void minusAssignOrResize(MultiArray<N, T, A> &, MultiMathOperand<E> const &);
00627 
00628 template <unsigned int N, class T, class A, class E>
00629 void multiplyAssignOrResize(MultiArray<N, T, A> &, MultiMathOperand<E> const &);
00630 
00631 template <unsigned int N, class T, class A, class E>
00632 void divideAssignOrResize(MultiArray<N, T, A> &, MultiMathOperand<E> const &);
00633 
00634 } // namespace detail
00635 
00636 } // namespace multi_math
00637 
00638 template <class T> class FindSum;
00639 
00640 struct UnsuitableTypeForExpandElements {};
00641 
00642 template <class T>
00643 struct ExpandElementResult
00644 {
00645     typedef UnsuitableTypeForExpandElements type;
00646 };
00647 
00648 template <class T>
00649 struct ExpandElementResult<std::complex<T> >
00650 {
00651     typedef T type;
00652     enum { size = 2 };
00653 };
00654 
00655 template <class T>
00656 class FFTWComplex;
00657 
00658 template <class T>
00659 struct ExpandElementResult<FFTWComplex<T> >
00660 {
00661     typedef T type;
00662     enum { size = 2 };
00663 };
00664 
00665 template <class T, int SIZE>
00666 struct ExpandElementResult<TinyVector<T, SIZE> >
00667 {
00668     typedef T type;
00669     enum { size = SIZE };
00670 };
00671 
00672 template <class T, unsigned int R, unsigned int G, unsigned int B>
00673 struct ExpandElementResult<RGBValue<T, R, G, B> >
00674 {
00675     typedef T type;
00676     enum { size = 3 };
00677 };
00678 
00679 #define VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(TYPE) \
00680 template <>  \
00681 struct ExpandElementResult<TYPE> \
00682 { \
00683     typedef TYPE type; \
00684     enum { size = 1 }; \
00685 }; \
00686 
00687 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(bool)
00688 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(char)
00689 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(signed char)
00690 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(signed short)
00691 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(signed int)
00692 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(signed long)
00693 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(signed long long)
00694 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(unsigned char)
00695 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(unsigned short)
00696 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(unsigned int)
00697 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(unsigned long)
00698 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(unsigned long long)
00699 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(float)
00700 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(double)
00701 VIGRA_DEFINE_EXPAND_ELEMENT_RESULT(long double)
00702 
00703 #undef VIGRA_DEFINE_EXPAND_ELEMENT_RESULT
00704 
00705 
00706 /********************************************************/
00707 /*                                                      */
00708 /*                       NormTraits                     */
00709 /*                                                      */
00710 /********************************************************/
00711 
00712 template <unsigned int N, class T, class C>
00713 struct NormTraits<MultiArrayView<N, T, C> >
00714 {
00715     typedef MultiArrayView<N, T, C>                                      Type;
00716     typedef typename NormTraits<T>::SquaredNormType SquaredNormType;
00717     typedef typename SquareRootTraits<SquaredNormType>::SquareRootResult NormType;
00718 };
00719 
00720 template <unsigned int N, class T, class A>
00721 struct NormTraits<MultiArray<N, T, A> >
00722 : public NormTraits<MultiArrayView<N, T, UnstridedArrayTag> >
00723 {
00724     typedef NormTraits<MultiArrayView<N, T, UnstridedArrayTag> > BaseType;
00725     typedef MultiArray<N, T, A>                                  Type;
00726     typedef typename BaseType::SquaredNormType                   SquaredNormType;
00727     typedef typename BaseType::NormType                          NormType;
00728 };
00729 
00730 /** \brief Base class for, and view to, \ref vigra::MultiArray.
00731 
00732 This class implements the interface of both MultiArray and
00733 MultiArrayView.  By default, MultiArrayViews are tagged as
00734 unstrided. If necessary, strided arrays are constructed automatically
00735 by calls to a variant of the bind...() function.
00736 
00737 In addition to the member functions described here, <tt>MultiArrayView</tt>
00738 and its subclasses support arithmetic and algebraic functions via the 
00739 module \ref MultiMathModule.
00740 
00741 If you want to apply an algorithm requiring an image to a
00742 <tt>MultiArrayView</tt> of appropriate (2-dimensional) shape, you can
00743 create a \ref vigra::BasicImageView that acts as a wrapper with the
00744 necessary interface -- see \ref MultiArrayToImage.
00745 
00746 The template parameter are as follows
00747 \code
00748     N: the array dimension
00749 
00750     T: the type of the array elements
00751 
00752     C: a tag determining whether the array's inner dimension is strided
00753        or not. An array is unstrided if the array elements occupy consecutive
00754        memory location, strided if there is an offset in between (e.g.
00755        when a view is created that skips every other array element).
00756        The compiler can generate faster code for unstrided arrays.
00757        Possible values: UnstridedArrayTag (default), StridedArrayTag
00758 \endcode
00759 
00760 <b>\#include</b>
00761 <vigra/multi_array.hxx>
00762 
00763 Namespace: vigra
00764 */
00765 template <unsigned int N, class T, class StrideTag>
00766 class MultiArrayView
00767 {
00768 public:
00769 
00770         /** the array's actual dimensionality.
00771             This ensures that MultiArrayView can also be used for
00772             scalars (that is, when <tt>N == 0</tt>). Calculated as:<br>
00773             \code
00774             actual_dimension = (N==0) ? 1 : N
00775             \endcode
00776          */
00777     enum ActualDimension { actual_dimension = (N==0) ? 1 : N };
00778 
00779         /** the array's value type
00780          */
00781     typedef T value_type;
00782 
00783         /** reference type (result of operator[])
00784          */
00785     typedef value_type &reference;
00786 
00787         /** const reference type (result of operator[] const)
00788          */
00789     typedef const value_type &const_reference;
00790 
00791         /** pointer type
00792          */
00793     typedef value_type *pointer;
00794 
00795         /** const pointer type
00796          */
00797     typedef const value_type *const_pointer;
00798 
00799         /** difference type (used for multi-dimensional offsets and indices)
00800          */
00801     typedef typename MultiArrayShape<actual_dimension>::type difference_type;
00802 
00803         /** size type
00804          */
00805     typedef difference_type size_type;
00806 
00807         /** difference and index type for a single dimension
00808          */
00809     typedef MultiArrayIndex difference_type_1;
00810 
00811         /** scan-order iterator (StridedScanOrderIterator) type
00812          */
00813     typedef StridedScanOrderIterator<actual_dimension, T, T &, T *> iterator;
00814 
00815         /** const scan-order iterator (StridedScanOrderIterator) type
00816          */
00817     typedef StridedScanOrderIterator<actual_dimension, T, T const &, T const *> const_iterator;
00818 
00819         /** traverser (MultiIterator) type
00820          */
00821     typedef typename vigra::detail::MultiIteratorChooser <
00822         StrideTag>::template Traverser <actual_dimension, T, T &, T *>::type traverser;
00823 
00824         /** const traverser (MultiIterator) type
00825          */
00826     typedef typename vigra::detail::MultiIteratorChooser <
00827         StrideTag>::template Traverser <actual_dimension, T, T const &, T const *>::type const_traverser;
00828 
00829         /** the view type associated with this array.
00830          */
00831     typedef MultiArrayView <N, T, StrideTag> view_type;
00832 
00833         /** the matrix type associated with this array.
00834          */
00835     typedef MultiArray <N, T> matrix_type;
00836 
00837 protected:
00838 
00839     typedef typename difference_type::value_type diff_zero_t;
00840 
00841         /** the shape of the image pointed to is stored here.
00842         */
00843     difference_type m_shape;
00844 
00845         /** the strides (offset of a sample to the next) for every dimension
00846             are stored here.
00847         */
00848     difference_type m_stride;
00849 
00850         /** pointer to the image.
00851          */
00852     pointer m_ptr;
00853 
00854     template <class U, class CN>
00855     void copyImpl(const MultiArrayView <N, U, CN>& rhs);
00856 
00857     template <class U, class CN>
00858     void swapDataImpl(MultiArrayView <N, U, CN> rhs);
00859 
00860     template <class CN>
00861     bool arraysOverlap(const MultiArrayView <N, T, CN>& rhs) const;
00862 
00863     template <class U, class CN>
00864     bool arraysOverlap(const MultiArrayView <N, U, CN>&) const
00865     {
00866         return false;
00867     }
00868     
00869     bool checkInnerStride(UnstridedArrayTag)
00870     {
00871         return m_stride[0] <= 1;
00872     }
00873     
00874     bool checkInnerStride(StridedArrayTag)
00875     {
00876         return true;
00877     }
00878 
00879 public:
00880 
00881         /** default constructor: create an invalid view,
00882          * i.e. hasData() returns false and size() is zero.
00883          */
00884     MultiArrayView ()
00885         : m_shape (diff_zero_t(0)), m_stride (diff_zero_t(0)), m_ptr (0)
00886     {}
00887 
00888         /** construct from shape and pointer
00889          */
00890     MultiArrayView (const difference_type &shape, pointer ptr)
00891     : m_shape (shape),
00892       m_stride (detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape)),
00893       m_ptr (ptr)
00894     {}
00895 
00896         /** Construct from shape, strides (offset of a sample to the
00897             next) for every dimension, and pointer.  (Note that
00898             strides are not given in bytes, but in offset steps of the
00899             respective pointer type.)
00900          */
00901     MultiArrayView (const difference_type &shape,
00902                     const difference_type &stride,
00903                     pointer ptr)
00904     : m_shape (shape),
00905       m_stride (stride),
00906       m_ptr (ptr)
00907     {
00908         vigra_precondition(checkInnerStride(StrideTag()),
00909             "MultiArrayView<..., UnstridedArrayTag>::MultiArrayView(): First dimension of given array is not unstrided.");
00910     }
00911     
00912         /** Conversion to a strided view.
00913          */
00914     operator MultiArrayView<N, T, StridedArrayTag>() const
00915     {
00916         return MultiArrayView<N, T, StridedArrayTag>(m_shape, m_stride, m_ptr);
00917     }
00918 
00919         /** Assignment. There are 3 cases:
00920 
00921             <ul>
00922             <li> When this <tt>MultiArrayView</tt> does not point to valid data
00923                  (e.g. after default construction), it becomes a copy of \a rhs.
00924             <li> When the shapes of the two arrays match, the array contents are copied.
00925             <li> Otherwise, a <tt>PreconditionViolation</tt> exception is thrown.
00926             </ul>
00927          */
00928     MultiArrayView & operator=(MultiArrayView const & rhs);
00929 
00930         /** Assignment of a differently typed MultiArrayView. Fails with
00931             <tt>PreconditionViolation</tt> exception when the shapes do not match.
00932          */
00933     template<class U, class C1>
00934     MultiArrayView & operator=(MultiArrayView<N, U, C1> const & rhs)
00935     {
00936         vigra_precondition(this->shape() == rhs.shape(),
00937             "MultiArrayView::operator=() size mismatch.");
00938         this->copyImpl(rhs);
00939         return *this;
00940     }
00941 
00942         /** Assignment of a scalar. Equivalent to MultiArrayView::init(v).
00943          */
00944     MultiArrayView & operator=(value_type const & v)
00945     {
00946         return init(v);
00947     }
00948 
00949         /** Add-assignment of a compatible MultiArrayView. Fails with
00950             <tt>PreconditionViolation</tt> exception when the shapes do not match.
00951          */
00952     template<class U, class C1>
00953     MultiArrayView & operator+=(MultiArrayView<N, U, C1> const & rhs);
00954 
00955         /** Subtract-assignment of a compatible MultiArrayView. Fails with
00956             <tt>PreconditionViolation</tt> exception when the shapes do not match.
00957          */
00958     template<class U, class C1>
00959     MultiArrayView & operator-=(MultiArrayView<N, U, C1> const & rhs);
00960 
00961         /** Multiply-assignment of a compatible MultiArrayView. Fails with
00962             <tt>PreconditionViolation</tt> exception when the shapes do not match.
00963          */
00964     template<class U, class C1>
00965     MultiArrayView & operator*=(MultiArrayView<N, U, C1> const & rhs);
00966 
00967         /** Divide-assignment of a compatible MultiArrayView. Fails with
00968             <tt>PreconditionViolation</tt> exception when the shapes do not match.
00969          */
00970     template<class U, class C1>
00971     MultiArrayView & operator/=(MultiArrayView<N, U, C1> const & rhs);
00972 
00973         /** Add-assignment of a scalar.
00974          */
00975     MultiArrayView & operator+=(T const & rhs)
00976     {
00977         detail::copyAddScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
00978         return *this;
00979     }
00980 
00981         /** Subtract-assignment of a scalar.
00982          */
00983     MultiArrayView & operator-=(T const & rhs)
00984     {
00985         detail::copySubScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
00986         return *this;
00987     }
00988 
00989         /** Multiply-assignment of a scalar.
00990          */
00991     MultiArrayView & operator*=(T const & rhs)
00992     {
00993         detail::copyMulScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
00994         return *this;
00995     }
00996 
00997         /** Divide-assignment of a scalar.
00998          */
00999     MultiArrayView & operator/=(T const & rhs)
01000     {
01001         detail::copyDivScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
01002         return *this;
01003     }
01004 
01005         /** Assignment of an array expression. Fails with
01006             <tt>PreconditionViolation</tt> exception when the shapes do not match.
01007          */
01008     template<class Expression>
01009     MultiArrayView & operator=(multi_math::MultiMathOperand<Expression> const & rhs)
01010     {
01011         multi_math::detail::assign(*this, rhs);
01012         return *this;
01013     }
01014 
01015         /** Add-assignment of an array expression. Fails with
01016             <tt>PreconditionViolation</tt> exception when the shapes do not match.
01017          */
01018     template<class Expression>
01019     MultiArrayView & operator+=(multi_math::MultiMathOperand<Expression> const & rhs)
01020     {
01021         multi_math::detail::plusAssign(*this, rhs);
01022         return *this;
01023     }
01024 
01025         /** Subtract-assignment of an array expression. Fails with
01026             <tt>PreconditionViolation</tt> exception when the shapes do not match.
01027          */
01028     template<class Expression>
01029     MultiArrayView & operator-=(multi_math::MultiMathOperand<Expression> const & rhs)
01030     {
01031         multi_math::detail::minusAssign(*this, rhs);
01032         return *this;
01033     }
01034 
01035         /** Multiply-assignment of an array expression. Fails with
01036             <tt>PreconditionViolation</tt> exception when the shapes do not match.
01037          */
01038     template<class Expression>
01039     MultiArrayView & operator*=(multi_math::MultiMathOperand<Expression> const & rhs)
01040     {
01041         multi_math::detail::multiplyAssign(*this, rhs);
01042         return *this;
01043     }
01044 
01045         /** Divide-assignment of an array expression. Fails with
01046             <tt>PreconditionViolation</tt> exception when the shapes do not match.
01047          */
01048     template<class Expression>
01049     MultiArrayView & operator/=(multi_math::MultiMathOperand<Expression> const & rhs)
01050     {
01051         multi_math::detail::divideAssign(*this, rhs);
01052         return *this;
01053     }
01054 
01055         /** array access.
01056          */
01057     reference operator[] (const difference_type &d)
01058     {
01059         VIGRA_ASSERT_INSIDE(d);
01060         return m_ptr [dot (d, m_stride)];
01061     }
01062 
01063         /** array access.
01064          */
01065     const_reference operator[] (const difference_type &d) const
01066     {
01067         VIGRA_ASSERT_INSIDE(d);
01068         return m_ptr [dot (d, m_stride)];
01069     }
01070 
01071         /** equivalent to bindInner(), when M < N.
01072          */
01073     template <int M>
01074     MultiArrayView <N-M, T, StridedArrayTag> operator[] (const TinyVector<MultiArrayIndex, M> &d) const
01075     {
01076         return bindInner(d);
01077     }
01078 
01079         /** Array access in scan-order sense.
01080             Mostly useful to support standard indexing for 1-dimensional multi-arrays,
01081             but works for any N. Use scanOrderIndexToCoordinate() and
01082             coordinateToScanOrderIndex() for conversion between indices and coordinates.
01083             
01084             <b>Note:</b> This function should not be used in the inner loop, because the
01085             conversion of the scan order index into a memory address is expensive
01086             (it must take into account that memory may not be consecutive for subarrays
01087             and/or strided arrays). Always prefer operator() if possible.            
01088          */
01089     reference operator[](difference_type_1 d)
01090     {
01091         VIGRA_ASSERT_INSIDE(scanOrderIndexToCoordinate(d));
01092         return m_ptr [detail::ScanOrderToOffset<actual_dimension>::exec(d, m_shape, m_stride)];
01093     }
01094 
01095         /** Array access in scan-order sense.
01096             Mostly useful to support standard indexing for 1-dimensional multi-arrays,
01097             but works for any N. Use scanOrderIndexToCoordinate() and
01098             coordinateToScanOrderIndex() for conversion between indices and coordinates.
01099              
01100             <b>Note:</b> This function should not be used in the inner loop, because the
01101             conversion of the scan order index into a memory address is expensive
01102             (it must take into account that memory may not be consecutive for subarrays
01103             and/or strided arrays). Always prefer operator() if possible.            
01104         */
01105     const_reference operator[](difference_type_1 d) const
01106     {
01107         VIGRA_ASSERT_INSIDE(scanOrderIndexToCoordinate(d));
01108         return m_ptr [detail::ScanOrderToOffset<actual_dimension>::exec(d, m_shape, m_stride)];
01109     }
01110 
01111         /** convert scan-order index to coordinate.
01112          */
01113     difference_type scanOrderIndexToCoordinate(difference_type_1 d) const
01114     {
01115         difference_type result;
01116         detail::ScanOrderToCoordinate<actual_dimension>::exec(d, m_shape, result);
01117         return result;
01118     }
01119 
01120         /** convert coordinate to scan-order index.
01121          */
01122     difference_type_1 coordinateToScanOrderIndex(const difference_type &d) const
01123     {
01124         return detail::CoordinateToScanOrder<actual_dimension>::exec(m_shape, d);
01125     }
01126 
01127         /** 1D array access. Use only if N == 1.
01128          */
01129     reference operator() (difference_type_1 x)
01130     {
01131         VIGRA_ASSERT_INSIDE(difference_type(x));
01132         return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x)];
01133     }
01134 
01135         /** 2D array access. Use only if N == 2.
01136          */
01137     reference operator() (difference_type_1 x, difference_type_1 y)
01138     {
01139         VIGRA_ASSERT_INSIDE(difference_type(x, y));
01140         return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x, y)];
01141     }
01142 
01143         /** 3D array access. Use only if N == 3.
01144          */
01145     reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z)
01146     {
01147         VIGRA_ASSERT_INSIDE(difference_type(x, y, z));
01148         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
01149     }
01150 
01151         /** 4D array access. Use only if N == 4.
01152          */
01153     reference operator() (difference_type_1 x, difference_type_1 y,
01154                           difference_type_1 z, difference_type_1 u)
01155     {
01156         VIGRA_ASSERT_INSIDE(difference_type(x, y, z, u));
01157         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
01158     }
01159 
01160         /** 5D array access. Use only if N == 5.
01161          */
01162     reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
01163                           difference_type_1 u, difference_type_1 v)
01164     {
01165         VIGRA_ASSERT_INSIDE(difference_type(x, y,z, u,v));
01166         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
01167     }
01168 
01169         /** 1D const array access. Use only if N == 1.
01170          */
01171     const_reference operator() (difference_type_1 x) const
01172     {
01173         VIGRA_ASSERT_INSIDE(difference_type(x));
01174         return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x)];
01175     }
01176 
01177         /** 2D const array access. Use only if N == 2.
01178          */
01179     const_reference operator() (difference_type_1 x, difference_type_1 y) const
01180     {
01181         VIGRA_ASSERT_INSIDE(difference_type(x, y));
01182         return m_ptr [detail::CoordinatesToOffest<StrideTag>::exec(m_stride, x, y)];
01183     }
01184 
01185         /** 3D const array access. Use only if N == 3.
01186          */
01187     const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z) const
01188     {
01189         VIGRA_ASSERT_INSIDE(difference_type(x,y,z));
01190         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
01191     }
01192 
01193         /** 4D const array access. Use only if N == 4.
01194          */
01195     const_reference operator() (difference_type_1 x, difference_type_1 y,
01196                                 difference_type_1 z, difference_type_1 u) const
01197     {
01198         VIGRA_ASSERT_INSIDE(difference_type(x,y,z,u));
01199         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
01200     }
01201 
01202         /** 5D const array access. Use only if N == 5.
01203          */
01204     const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
01205                                 difference_type_1 u, difference_type_1 v) const
01206     {
01207         VIGRA_ASSERT_INSIDE(difference_type(x,y,z,u,v));
01208         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
01209     }
01210 
01211         /** Init with a constant.
01212          */
01213     template <class U>
01214     MultiArrayView & init(const U & init)
01215     {
01216         if(hasData())
01217             detail::copyScalarMultiArrayData(traverser_begin(), shape(), init, MetaInt<actual_dimension-1>());
01218         return *this;
01219     }
01220 
01221 
01222         /** Copy the data of the right-hand array (array shapes must match).
01223          */
01224     void copy(const MultiArrayView & rhs)
01225     {
01226         if(this == &rhs)
01227             return;
01228         this->copyImpl(rhs);
01229     }
01230 
01231         /** Copy the data of the right-hand array (array shapes must match).
01232          */
01233     template <class U, class CN>
01234     void copy(const MultiArrayView <N, U, CN>& rhs)
01235     {
01236         this->copyImpl(rhs);
01237     }
01238 
01239         /** swap the data between two MultiArrayView objects.
01240 
01241             The shapes of the two array must match.
01242         */
01243     void swapData(MultiArrayView rhs)
01244     {
01245         if(this != &rhs)
01246             swapDataImpl(rhs);
01247     }
01248 
01249         /** swap the data between two MultiArrayView objects.
01250 
01251             The shapes of the two array must match.
01252         */
01253     template <class T2, class C2>
01254     void swapData(MultiArrayView <N, T2, C2> rhs)
01255     {
01256         swapDataImpl(rhs);
01257     }
01258     
01259         /** check whether the array is unstrided (i.e. has consecutive memory) up 
01260             to the given dimension.
01261 
01262             \a dimension can range from 0 ... N-1. If a certain dimension is unstrided, 
01263             all lower dimensions are also unstrided.
01264         */
01265     bool isUnstrided(unsigned int dimension = N-1) const
01266     {
01267         difference_type p = shape() - difference_type(1);
01268         for(unsigned int k = dimension+1; k < N; ++k)
01269             p[k] = 0;
01270         return (&operator[](p) - m_ptr) == coordinateToScanOrderIndex(p);
01271     }
01272 
01273         /** bind the M outmost dimensions to certain indices.
01274             this reduces the dimensionality of the image to
01275             max { 1, N-M }.
01276 
01277             <b>Usage:</b>
01278             \code
01279             // create a 3D array of size 40x30x20
01280             typedef MultiArray<3, double>::difference_type Shape;
01281             MultiArray<3, double> array3(Shape(40, 30, 20));
01282 
01283             // get a 1D array by fixing index 1 to 12, and index 2 to 10
01284             MultiArrayView <1, double> array1 = array3.bindOuter(TinyVector<MultiArrayIndex, 2>(12, 10));
01285             \endcode
01286         */
01287     template <int M, class Index>
01288     MultiArrayView <N-M, T, StrideTag> bindOuter(const TinyVector <Index, M> &d) const;
01289 
01290         /** bind the M innermost dimensions to certain indices.
01291             this reduces the dimensionality of the image to
01292             max { 1, N-M }.
01293 
01294             <b>Usage:</b>
01295             \code
01296             // create a 3D array of size 40x30x20
01297             typedef MultiArray<3, double>::difference_type Shape;
01298             MultiArray<3, double> array3(Shape(40, 30, 20));
01299 
01300             // get a 1D array by fixing index 0 to 12, and index 1 to 10
01301             MultiArrayView <1, double, StridedArrayTag> array1 = array3.bindInner(TinyVector<MultiArrayIndex, 2>(12, 10));
01302             \endcode
01303         */
01304     template <int M, class Index>
01305     MultiArrayView <N-M, T, StridedArrayTag> bindInner(const TinyVector <Index, M> &d) const;
01306 
01307         /** bind dimension M to index d.
01308             this reduces the dimensionality of the image to
01309             max { 1, N-1 }.
01310 
01311             <b>Usage:</b>
01312             \code
01313             // create a 3D array of size 40x30x20
01314             typedef MultiArray<3, double>::difference_type Shape;
01315             MultiArray<3, double> array3(Shape(40, 30, 20));
01316 
01317             // get a 2D array by fixing index 1 to 12
01318             MultiArrayView <2, double> array2 = array3.bind<1>(12);
01319 
01320             // get a 2D array by fixing index 0 to 23
01321             MultiArrayView <2, double, StridedArrayTag> array2a = array3.bind<0>(23);
01322             \endcode
01323          */
01324     template <unsigned int M>
01325     MultiArrayView <N-1, T, typename vigra::detail::MaybeStrided<StrideTag, M>::type >
01326     bind (difference_type_1 d) const;
01327 
01328         /** bind the outmost dimension to a certain index.
01329             this reduces the dimensionality of the image to
01330             max { 1, N-1 }.
01331 
01332             <b>Usage:</b>
01333             \code
01334             // create a 3D array of size 40x30x20
01335             typedef MultiArray<3, double>::difference_type Shape;
01336             MultiArray<3, double> array3(Shape(40, 30, 20));
01337 
01338             // get a 2D array by fixing the outermost index (i.e. index 2) to 12
01339             MultiArrayView <2, double> array2 = array3.bindOuter(12);
01340             \endcode
01341         */
01342     MultiArrayView <N-1, T, StrideTag> bindOuter (difference_type_1 d) const;
01343 
01344         /** bind the innermost dimension to a certain index.
01345             this reduces the dimensionality of the image to
01346             max { 1, N-1 }.
01347 
01348             <b>Usage:</b>
01349             \code
01350             // create a 3D array of size 40x30x20
01351             typedef MultiArray<3, double>::difference_type Shape;
01352             MultiArray<3, double> array3(Shape(40, 30, 20));
01353 
01354             // get a 2D array by fixing the innermost index (i.e. index 0) to 23
01355             MultiArrayView <2, double, StridedArrayTag> array2 = array3.bindInner(23);
01356             \endcode
01357         */
01358     MultiArrayView <N-1, T, StridedArrayTag> bindInner (difference_type_1 d) const;
01359 
01360         /** bind dimension m to index d.
01361             this reduces the dimensionality of the image to
01362             max { 1, N-1 }.
01363 
01364             <b>Usage:</b>
01365             \code
01366             // create a 3D array of size 40x30x20
01367             typedef MultiArray<3, double>::difference_type Shape;
01368             MultiArray<3, double> array3(Shape(40, 30, 20));
01369 
01370             // get a 2D array by fixing index 2 to 15
01371             MultiArrayView <2, double, StridedArrayTag> array2 = array3.bindAt(2, 15);
01372             \endcode
01373          */
01374     MultiArrayView <N-1, T, StridedArrayTag>
01375     bindAt (difference_type_1 m, difference_type_1 d) const;
01376     
01377     
01378         /** Create a view to channel 'i' of a vector-like value type. Possible value types
01379             (of the original array) are: \ref TinyVector, \ref RGBValue, \ref FFTWComplex, 
01380             and <tt>std::complex</tt>. The list can be extended to any type whose memory
01381             layout is equivalent to a fixed-size C array, by specializing 
01382             <tt>ExpandElementResult</tt>.
01383 
01384             <b>Usage:</b>
01385             \code
01386                 MultiArray<2, RGBValue<float> > rgb_image(Shape2(w, h));
01387                 
01388                 MultiArrayView<2, float, StridedArrayTag> red   = rgb_image.bindElementChannel(0);
01389                 MultiArrayView<2, float, StridedArrayTag> green = rgb_image.bindElementChannel(1);
01390                 MultiArrayView<2, float, StridedArrayTag> blue  = rgb_image.bindElementChannel(2);
01391             \endcode
01392         */
01393     MultiArrayView <N, typename ExpandElementResult<T>::type, StridedArrayTag> 
01394     bindElementChannel(difference_type_1 i) const
01395     {
01396         vigra_precondition(0 <= i && i < ExpandElementResult<T>::size,
01397               "MultiArrayView::bindElementChannel(i): 'i' out of range.");
01398         return expandElements(0).bindInner(i);
01399     }
01400 
01401         /** Create a view where a vector-like element type is expanded into a new 
01402             array dimension. The new dimension is inserted at index position 'd',
01403             which must be between 0 and N inclusive.
01404             
01405             Possible value types of the original array are: \ref TinyVector, \ref RGBValue, 
01406             \ref FFTWComplex, <tt>std::complex</tt>, and the built-in number types (in this 
01407             case, <tt>expandElements</tt> is equivalent to <tt>insertSingletonDimension</tt>). 
01408             The list of supported types can be extended to any type whose memory
01409             layout is equivalent to a fixed-size C array, by specializing 
01410             <tt>ExpandElementResult</tt>.
01411 
01412             <b>Usage:</b>
01413             \code
01414                 MultiArray<2, RGBValue<float> > rgb_image(Shape2(w, h));
01415                 
01416                 MultiArrayView<3, float, StridedArrayTag> multiband_image = rgb_image.expandElements(2);
01417             \endcode
01418         */
01419     MultiArrayView <N+1, typename ExpandElementResult<T>::type, StridedArrayTag> 
01420     expandElements(difference_type_1 d) const;
01421     
01422         /** Add a singleton dimension (dimension of length 1).
01423 
01424             Singleton dimensions don't change the size of the data, but introduce
01425             a new index that can only take the value 0. This is mainly useful for
01426             the 'reduce mode' of transformMultiArray() and combineTwoMultiArrays(),
01427             because these functions require the source and destination arrays to
01428             have the same number of dimensions.
01429 
01430             The range of \a i must be <tt>0 <= i <= N</tt>. The new dimension will become
01431             the i'th index, and the old indices from i upwards will shift one
01432             place to the right.
01433 
01434             <b>Usage:</b>
01435 
01436             Suppose we want have a 2D array and want to create a 1D array that contains
01437             the row average of the first array.
01438             \code
01439             typedef MultiArrayShape<2>::type Shape2;
01440             MultiArray<2, double> original(Shape2(40, 30));
01441 
01442             typedef MultiArrayShape<1>::type Shape1;
01443             MultiArray<1, double> rowAverages(Shape1(30));
01444 
01445             // temporarily add a singleton dimension to the destination array
01446             transformMultiArray(srcMultiArrayRange(original),
01447                                 destMultiArrayRange(rowAverages.insertSingletonDimension(0)),
01448                                 FindAverage<double>());
01449             \endcode
01450          */
01451     MultiArrayView <N+1, T, StrideTag>
01452     insertSingletonDimension (difference_type_1 i) const;
01453 
01454         /** Create a view to the diagonal elements of the array.
01455         
01456             This produces a 1D array view whose size equals the size
01457             of the shortest dimension of the original array.
01458 
01459             <b>Usage:</b>
01460             \code
01461             // create a 3D array of size 40x30x20
01462             typedef MultiArray<3, double>::difference_type Shape;
01463             MultiArray<3, double> array3(Shape(40, 30, 20));
01464 
01465             // get a view to the diagonal elements
01466             MultiArrayView <1, double, StridedArrayTag> diagonal = array3.diagonal();
01467             assert(diagonal.shape(0) == 20);
01468             \endcode
01469         */
01470     MultiArrayView<1, T, StridedArrayTag> diagonal() const
01471     {
01472         return MultiArrayView<1, T, StridedArrayTag>(Shape1(vigra::min(m_shape)), 
01473                                                      Shape1(vigra::sum(m_stride)), m_ptr);
01474     }
01475 
01476         /** create a rectangular subarray that spans between the
01477             points p and q, where p is in the subarray, q not.
01478 
01479             <b>Usage:</b>
01480             \code
01481             // create a 3D array of size 40x30x20
01482             typedef MultiArray<3, double>::difference_type Shape;
01483             MultiArray<3, double> array3(Shape(40, 30, 20));
01484 
01485             // get a subarray set is smaller by one element at all sides
01486             MultiArrayView <3, double> subarray = array3.subarray(Shape(1,1,1), Shape(39, 29, 19));
01487             \endcode
01488         */
01489     MultiArrayView subarray (const difference_type &p,
01490                              const difference_type &q) const
01491     {
01492         const difference_type_1 offset = dot (m_stride, p);
01493         return MultiArrayView (q - p, m_stride, m_ptr + offset);
01494     }
01495 
01496         /** apply an additional striding to the image, thereby reducing
01497             the shape of the array.
01498             for example, multiplying the stride of dimension one by three
01499             turns an appropriately laid out (interleaved) rgb image into
01500             a single band image.
01501         */
01502     MultiArrayView <N, T, StridedArrayTag>
01503     stridearray (const difference_type &s) const
01504     {
01505         difference_type shape = m_shape;
01506         for (unsigned int i = 0; i < actual_dimension; ++i)
01507             shape [i] /= s [i];
01508         return MultiArrayView <N, T, StridedArrayTag>(shape, m_stride * s, m_ptr);
01509     }
01510 
01511         /** Transpose an array. If N==2, this implements the usual matrix transposition.
01512             For N > 2, it reverses the order of the indices.
01513 
01514             <b>Usage:</b><br>
01515             \code
01516             typedef MultiArray<2, double>::difference_type Shape;
01517             MultiArray<2, double> array(10, 20);
01518 
01519             MultiArray<2, double, StridedArrayTag> transposed = array.transpose();
01520 
01521             for(int i=0; i<array.shape(0), ++i)
01522                 for(int j=0; j<array.shape(1); ++j)
01523                     assert(array(i, j) == transposed(j, i));
01524             \endcode
01525         */
01526     MultiArrayView <N, T, StridedArrayTag>
01527     transpose () const
01528     {
01529         difference_type shape(m_shape.begin(), difference_type::ReverseCopy),
01530                         stride(m_stride.begin(), difference_type::ReverseCopy);
01531         return MultiArrayView <N, T, StridedArrayTag>(shape, stride, m_ptr);
01532     }
01533 
01534         /** permute the dimensions of the array.
01535             The function exchanges the meaning of the dimensions without copying the data.
01536             In case of a 2-dimensional array, this is simply array transposition. In higher dimensions,
01537             there are more possibilities.
01538 
01539             <b>Usage:</b><br>
01540             \code
01541             typedef MultiArray<2, double>::difference_type Shape;
01542             MultiArray<2, double> array(10, 20);
01543 
01544             MultiArray<2, double, StridedArrayTag> transposed = array.permuteDimensions(Shape(1,0));
01545 
01546             for(int i=0; i<array.shape(0), ++i)
01547                 for(int j=0; j<array.shape(1); ++j)
01548                     assert(array(i, j) == transposed(j, i));
01549             \endcode
01550         */
01551     MultiArrayView <N, T, StridedArrayTag>
01552     permuteDimensions (const difference_type &s) const;
01553 
01554         /** Permute the dimensions of the array so that the strides are in ascending order.
01555             Determines the appropriate permutation and then calls permuteDimensions().
01556         */
01557     MultiArrayView <N, T, StridedArrayTag>
01558     permuteStridesAscending() const;
01559     
01560         /** Permute the dimensions of the array so that the strides are in descending order.
01561             Determines the appropriate permutation and then calls permuteDimensions().
01562         */
01563     MultiArrayView <N, T, StridedArrayTag>
01564     permuteStridesDescending() const;
01565     
01566         /** Compute the ordering of the strides in this array.
01567             The result is describes the current permutation of the axes relative 
01568             to the standard ascending stride order.
01569         */
01570     difference_type strideOrdering() const
01571     {
01572         return strideOrdering(m_stride);
01573     }
01574     
01575         /** Compute the ordering of the given strides.
01576             The result is describes the current permutation of the axes relative 
01577             to the standard ascending stride order.
01578         */
01579     static difference_type strideOrdering(difference_type strides);
01580 
01581         /** number of the elements in the array.
01582          */
01583     difference_type_1 elementCount () const
01584     {
01585         difference_type_1 ret = m_shape[0];
01586         for(int i = 1; i < actual_dimension; ++i)
01587             ret *= m_shape[i];
01588         return ret;
01589     }
01590 
01591         /** number of the elements in the array.
01592             Same as <tt>elementCount()</tt>. Mostly useful to support the std::vector interface.
01593          */
01594     difference_type_1 size () const
01595     {
01596         return elementCount();
01597     }
01598 
01599         /** return the array's shape.
01600          */
01601     const difference_type & shape () const
01602     {
01603         return m_shape;
01604     }
01605 
01606         /** return the array's size at a certain dimension.
01607          */
01608     difference_type_1 size (difference_type_1 n) const
01609     {
01610         return m_shape [n];
01611     }
01612 
01613         /** return the array's shape at a certain dimension
01614             (same as <tt>size(n)</tt>).
01615          */
01616     difference_type_1 shape (difference_type_1 n) const
01617     {
01618         return m_shape [n];
01619     }
01620 
01621         /** return the array's stride for every dimension.
01622          */
01623     const difference_type & stride () const
01624     {
01625         return m_stride;
01626     }
01627 
01628         /** return the array's stride at a certain dimension.
01629          */
01630     difference_type_1 stride (int n) const
01631     {
01632         return m_stride [n];
01633     }
01634 
01635         /** check whether two arrays are elementwise equal.
01636          */
01637     template <class U, class C1>
01638     bool operator==(MultiArrayView<N, U, C1> const & rhs) const
01639     {
01640         if(this->shape() != rhs.shape())
01641             return false;
01642         return detail::equalityOfMultiArrays(traverser_begin(), shape(), rhs.traverser_begin(), MetaInt<actual_dimension-1>());
01643     }
01644 
01645         /** check whether two arrays are not elementwise equal.
01646             Also true when the two arrays have different shapes.
01647          */
01648     template <class U, class C1>
01649     bool operator!=(MultiArrayView<N, U, C1> const & rhs) const
01650     {
01651         return !operator==(rhs);
01652     }
01653 
01654         /** check whether the given point is in the array range.
01655          */
01656     bool isInside (difference_type const & p) const
01657     {
01658         for(int d=0; d<actual_dimension; ++d)
01659             if(p[d] < 0 || p[d] >= shape(d))
01660                 return false;
01661         return true;
01662     }
01663 
01664         /** Check if the array contains only non-zero elements (or if all elements
01665             are 'true' if the value type is 'bool').
01666          */
01667     bool all() const
01668     {
01669         bool res = true;
01670         detail::reduceOverMultiArray(traverser_begin(), shape(),
01671                                      res, 
01672                                      detail::AllTrueReduceFunctor(),
01673                                      MetaInt<actual_dimension-1>());
01674         return res;
01675     }
01676 
01677         /** Check if the array contains a non-zero element (or an element
01678             that is 'true' if the value type is 'bool').
01679          */
01680     bool any() const
01681     {
01682         bool res = false;
01683         detail::reduceOverMultiArray(traverser_begin(), shape(),
01684                                      res, 
01685                                      detail::AnyTrueReduceFunctor(),
01686                                      MetaInt<actual_dimension-1>());
01687         return res;
01688     }
01689 
01690         /** Find the minimum and maximum element in this array. 
01691         See \ref FeatureAccumulators for a general feature 
01692         extraction framework.
01693          */
01694     void minmax(T * minimum, T * maximum) const
01695     {
01696         std::pair<T, T> res(NumericTraits<T>::max(), NumericTraits<T>::min());
01697         detail::reduceOverMultiArray(traverser_begin(), shape(),
01698                                      res, 
01699                                      detail::MinmaxReduceFunctor(),
01700                                      MetaInt<actual_dimension-1>());
01701         *minimum = res.first;
01702         *maximum = res.second;
01703     }
01704 
01705         /** Compute the mean and variance of the values in this array. 
01706         See \ref FeatureAccumulators for a general feature 
01707         extraction framework.
01708          */
01709     template <class U>
01710     void meanVariance(U * mean, U * variance) const
01711     {
01712         typedef typename NumericTraits<U>::RealPromote R;
01713     R zero;
01714         triple<double, R, R> res(0.0, zero, zero);
01715         detail::reduceOverMultiArray(traverser_begin(), shape(),
01716                                      res, 
01717                                      detail::MeanVarianceReduceFunctor(),
01718                                      MetaInt<actual_dimension-1>());
01719         *mean     = res.second;
01720         *variance = res.third / res.first;
01721     }
01722 
01723         /** Compute the sum of the array elements.
01724 
01725             You must provide the type of the result by an explicit template parameter:
01726             \code
01727             MultiArray<2, UInt8> A(width, height);
01728             
01729             double sum = A.sum<double>();
01730             \endcode
01731          */
01732     template <class U>
01733     U sum() const
01734     {
01735         U res = NumericTraits<U>::zero();
01736         detail::reduceOverMultiArray(traverser_begin(), shape(),
01737                                      res, 
01738                                      detail::SumReduceFunctor(),
01739                                      MetaInt<actual_dimension-1>());
01740         return res;
01741     }
01742 
01743         /** Compute the sum of the array elements over selected axes.
01744         
01745             \arg sums must have the same shape as this array, except for the
01746             axes along which the sum is to be accumulated. These axes must be 
01747             singletons. Note that you must include <tt>multi_pointoperators.hxx</tt>
01748             for this function to work.
01749 
01750             <b>Usage:</b>
01751             \code
01752             #include <vigra/multi_array.hxx>
01753             #include <vigra/multi_pointoperators.hxx>
01754             
01755             MultiArray<2, double> A(rows, cols);
01756             ... // fill A
01757             
01758             // make the first axis a singleton to sum over the first index
01759             MultiArray<2, double> rowSums(1, cols);
01760             A.sum(rowSums);
01761             
01762             // this is equivalent to
01763             transformMultiArray(srcMultiArrayRange(A),
01764                                 destMultiArrayRange(rowSums),
01765                                 FindSum<double>());
01766             \endcode
01767          */
01768     template <class U, class S>
01769     void sum(MultiArrayView<N, U, S> sums) const
01770     {
01771         transformMultiArray(srcMultiArrayRange(*this),
01772                             destMultiArrayRange(sums),
01773                             FindSum<U>());
01774     }
01775 
01776         /** Compute the product of the array elements.
01777 
01778             You must provide the type of the result by an explicit template parameter:
01779             \code
01780             MultiArray<2, UInt8> A(width, height);
01781             
01782             double prod = A.product<double>();
01783             \endcode
01784          */
01785     template <class U>
01786     U product() const
01787     {
01788         U res = NumericTraits<U>::one();
01789         detail::reduceOverMultiArray(traverser_begin(), shape(),
01790                                      res, 
01791                                      detail::ProdReduceFunctor(),
01792                                      MetaInt<actual_dimension-1>());
01793         return res;
01794     }
01795 
01796         /** Compute the squared Euclidean norm of the array (sum of squares of the array elements).
01797          */
01798     typename NormTraits<MultiArrayView>::SquaredNormType 
01799     squaredNorm() const
01800     {
01801         typedef typename NormTraits<MultiArrayView>::SquaredNormType SquaredNormType;
01802         SquaredNormType res = NumericTraits<SquaredNormType>::zero();
01803         detail::reduceOverMultiArray(traverser_begin(), shape(),
01804                                      res, 
01805                                      detail::SquaredL2NormReduceFunctor(),
01806                                      MetaInt<actual_dimension-1>());
01807         return res;
01808     }
01809 
01810         /** Compute various norms of the array.
01811             The norm is determined by parameter \a type:
01812 
01813             <ul>
01814             <li> type == 0: maximum norm (L-infinity): maximum of absolute values of the array elements
01815             <li> type == 1: Manhattan norm (L1): sum of absolute values of the array elements
01816             <li> type == 2: Euclidean norm (L2): square root of <tt>squaredNorm()</tt> when \a useSquaredNorm is <tt>true</tt>,<br>
01817                  or direct algorithm that avoids underflow/overflow otherwise.
01818             </ul>
01819 
01820             Parameter \a useSquaredNorm has no effect when \a type != 2. Defaults: compute L2 norm as square root of
01821             <tt>squaredNorm()</tt>.
01822          */
01823     typename NormTraits<MultiArrayView>::NormType 
01824     norm(int type = 2, bool useSquaredNorm = true) const;
01825 
01826         /** return the pointer to the image data
01827          */
01828     pointer data () const
01829     {
01830         return m_ptr;
01831     }
01832     
01833     pointer & unsafePtr()
01834     {
01835         return m_ptr;
01836     }
01837 
01838         /**
01839          * returns true iff this view refers to valid data,
01840          * i.e. data() is not a NULL pointer.  (this is false after
01841          * default construction.)
01842          */
01843     bool hasData () const
01844     {
01845         return m_ptr != 0;
01846     }
01847 
01848         /** returns a scan-order iterator pointing
01849             to the first array element.
01850         */
01851     iterator begin()
01852     {
01853         return iterator(m_ptr, m_shape, m_stride);
01854     }
01855 
01856         /** returns a const scan-order iterator pointing
01857             to the first array element.
01858         */
01859     const_iterator begin() const
01860     {
01861         return const_iterator(m_ptr, m_shape, m_stride);
01862     }
01863 
01864         /** returns a scan-order iterator pointing
01865             beyond the last array element.
01866         */
01867     iterator end()
01868     {
01869         return begin().getEndIterator();
01870     }
01871 
01872         /** returns a const scan-order iterator pointing
01873             beyond the last array element.
01874         */
01875     const_iterator end() const
01876     {
01877         return begin().getEndIterator();
01878     }
01879 
01880         /** returns the N-dimensional MultiIterator pointing
01881             to the first element in every dimension.
01882         */
01883     traverser traverser_begin ()
01884     {
01885         traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
01886         return ret;
01887     }
01888 
01889         /** returns the N-dimensional MultiIterator pointing
01890             to the const first element in every dimension.
01891         */
01892     const_traverser traverser_begin () const
01893     {
01894         const_traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
01895         return ret;
01896     }
01897 
01898         /** returns the N-dimensional MultiIterator pointing
01899             beyond the last element in dimension N, and to the
01900             first element in every other dimension.
01901         */
01902     traverser traverser_end ()
01903     {
01904         traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
01905         ret += m_shape [actual_dimension-1];
01906         return ret;
01907     }
01908 
01909         /** returns the N-dimensional const MultiIterator pointing
01910             beyond the last element in dimension N, and to the
01911             first element in every other dimension.
01912         */
01913     const_traverser traverser_end () const
01914     {
01915         const_traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
01916         ret += m_shape [actual_dimension-1];
01917         return ret;
01918     }
01919 
01920     view_type view ()
01921     {
01922         return *this;
01923     }
01924 };
01925 
01926 template <unsigned int N, class T, class StrideTag>
01927 MultiArrayView<N, T, StrideTag> &
01928 MultiArrayView <N, T, StrideTag>::operator=(MultiArrayView const & rhs)
01929 {
01930     if(this == &rhs)
01931         return *this;
01932     vigra_precondition(this->shape() == rhs.shape() || m_ptr == 0,
01933         "MultiArrayView::operator=(MultiArrayView const &) size mismatch.");
01934     if(m_ptr == 0)
01935     {
01936         m_shape  = rhs.m_shape;
01937         m_stride = rhs.m_stride;
01938         m_ptr    = rhs.m_ptr;
01939     }
01940     else
01941         this->copyImpl(rhs);
01942     return *this;
01943 }
01944 
01945 template <unsigned int N, class T, class StrideTag>
01946 template <class CN>
01947 bool
01948 MultiArrayView <N, T, StrideTag>::arraysOverlap(const MultiArrayView <N, T, CN>& rhs) const
01949 {
01950     vigra_precondition (shape () == rhs.shape (),
01951         "MultiArrayView::arraysOverlap(): shape mismatch.");
01952     const_pointer first_element = this->m_ptr,
01953                   last_element = first_element + dot(this->m_shape - difference_type(1), this->m_stride);
01954     typename MultiArrayView <N, T, CN>::const_pointer
01955            rhs_first_element = rhs.data(),
01956            rhs_last_element = rhs_first_element + dot(rhs.shape() - difference_type(1), rhs.stride());
01957     return !(last_element < rhs_first_element || rhs_last_element < first_element);
01958 }
01959 
01960 template <unsigned int N, class T, class StrideTag>
01961 template <class U, class CN>
01962 void
01963 MultiArrayView <N, T, StrideTag>::copyImpl(const MultiArrayView <N, U, CN>& rhs)
01964 {
01965     if(!arraysOverlap(rhs))
01966     {
01967         // no overlap -- can copy directly
01968         detail::copyMultiArrayData(rhs.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>());
01969     }
01970     else
01971     {
01972         // overlap: we got different views to the same data -- copy to intermediate memory in order to avoid
01973         // overwriting elements that are still needed on the rhs.
01974         MultiArray<N, T> tmp(rhs);
01975         detail::copyMultiArrayData(tmp.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>());
01976     }
01977 }
01978 
01979 #define VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(name, op) \
01980 template <unsigned int N, class T, class StrideTag> \
01981 template<class U, class C1> \
01982 MultiArrayView<N, T, StrideTag> &  \
01983 MultiArrayView <N, T, StrideTag>::operator op(MultiArrayView<N, U, C1> const & rhs) \
01984 { \
01985     vigra_precondition(this->shape() == rhs.shape(), "MultiArrayView::operator" #op "() size mismatch."); \
01986     if(!arraysOverlap(rhs)) \
01987     { \
01988         detail::name##MultiArrayData(rhs.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>()); \
01989     } \
01990     else \
01991     { \
01992         MultiArray<N, T> tmp(rhs); \
01993         detail::name##MultiArrayData(tmp.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>()); \
01994     } \
01995     return *this; \
01996 }
01997 
01998 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copyAdd, +=)
01999 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copySub, -=)
02000 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copyMul, *=)
02001 VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT(copyDiv, /=)
02002 
02003 #undef VIGRA_MULTI_ARRAY_COMPUTED_ASSIGNMENT
02004 
02005 template <unsigned int N, class T, class StrideTag>
02006 template <class U, class CN>
02007 void
02008 MultiArrayView <N, T, StrideTag>::swapDataImpl(MultiArrayView <N, U, CN> rhs)
02009 {
02010     vigra_precondition (shape () == rhs.shape (),
02011         "MultiArrayView::swapData(): shape mismatch.");
02012 
02013     // check for overlap of this and rhs
02014     const_pointer first_element = this->m_ptr,
02015                   last_element = first_element + dot(this->m_shape - difference_type(1), this->m_stride);
02016     typename MultiArrayView <N, U, CN>::const_pointer
02017            rhs_first_element = rhs.data(),
02018            rhs_last_element = rhs_first_element + dot(rhs.shape() - difference_type(1), rhs.stride());
02019     if(last_element < rhs_first_element || rhs_last_element < first_element)
02020     {
02021         // no overlap -- can swap directly
02022         detail::swapDataImpl(traverser_begin(), shape(), rhs.traverser_begin(), MetaInt<actual_dimension-1>());
02023     }
02024     else
02025     {
02026         // overlap: we got different views to the same data -- copy to intermediate memory in order to avoid
02027         // overwriting elements that are still needed.
02028         MultiArray<N, T> tmp(*this);
02029         copy(rhs);
02030         rhs.copy(tmp);
02031     }
02032 }
02033 
02034 template <unsigned int N, class T, class StrideTag>
02035 MultiArrayView <N, T, StridedArrayTag>
02036 MultiArrayView <N, T, StrideTag>::permuteDimensions (const difference_type &s) const
02037 {
02038     difference_type shape, stride, check((typename difference_type::value_type)0);
02039     for (unsigned int i = 0; i < actual_dimension; ++i)
02040     {
02041         shape[i]  = m_shape[s[i]];
02042         stride[i] = m_stride[s[i]];
02043         ++check[s[i]];
02044     }
02045     vigra_precondition(check == difference_type(1),
02046        "MultiArrayView::permuteDimensions(): every dimension must occur exactly once.");
02047     return MultiArrayView <N, T, StridedArrayTag>(shape, stride, m_ptr);
02048 }
02049 
02050 template <unsigned int N, class T, class StrideTag>
02051 typename MultiArrayView <N, T, StrideTag>::difference_type 
02052 MultiArrayView <N, T, StrideTag>::strideOrdering(difference_type stride)
02053 {
02054     difference_type permutation;
02055     for(int k=0; k<(int)N; ++k)
02056         permutation[k] = k;
02057     for(int k=0; k<(int)N-1; ++k)
02058     {
02059         int smallest = k;
02060         for(int j=k+1; j<(int)N; ++j)
02061         {
02062             if(stride[j] < stride[smallest])
02063                 smallest = j;
02064         }
02065         if(smallest != k)
02066         {
02067             std::swap(stride[k], stride[smallest]);
02068             std::swap(permutation[k], permutation[smallest]);
02069         }
02070     }
02071     difference_type ordering;
02072     for(unsigned int k=0; k<N; ++k)
02073         ordering[permutation[k]] = k;
02074     return ordering;
02075 }
02076 
02077 template <unsigned int N, class T, class StrideTag>
02078 MultiArrayView <N, T, StridedArrayTag>
02079 MultiArrayView <N, T, StrideTag>::permuteStridesAscending() const
02080 {
02081     difference_type ordering(strideOrdering(m_stride)), permutation;
02082     for(MultiArrayIndex k=0; k<N; ++k)
02083         permutation[ordering[k]] = k;
02084     return permuteDimensions(permutation);
02085 }
02086 
02087 template <unsigned int N, class T, class StrideTag>
02088 MultiArrayView <N, T, StridedArrayTag>
02089 MultiArrayView <N, T, StrideTag>::permuteStridesDescending() const
02090 {
02091     difference_type ordering(strideOrdering(m_stride)), permutation;
02092     for(MultiArrayIndex k=0; k<N; ++k)
02093         permutation[ordering[N-1-k]] = k;
02094     return permuteDimensions(permutation);
02095 }
02096 
02097 template <unsigned int N, class T, class StrideTag>
02098 template <int M, class Index>
02099 MultiArrayView <N-M, T, StrideTag>
02100 MultiArrayView <N, T, StrideTag>::bindOuter (const TinyVector <Index, M> &d) const
02101 {
02102     TinyVector <MultiArrayIndex, M> stride;
02103     stride.init (m_stride.begin () + N-M, m_stride.end ());
02104     pointer ptr = m_ptr + dot (d, stride);
02105     static const int NNew = (N-M == 0) ? 1 : N-M;
02106     TinyVector <MultiArrayIndex, NNew> inner_shape, inner_stride;
02107     if (N-M == 0)
02108     {
02109         inner_shape [0] = 1;
02110         inner_stride [0] = 0;
02111     }
02112     else
02113     {
02114         inner_shape.init (m_shape.begin (), m_shape.end () - M);
02115         inner_stride.init (m_stride.begin (), m_stride.end () - M);
02116     }
02117     return MultiArrayView <N-M, T, StrideTag> (inner_shape, inner_stride, ptr);
02118 }
02119 
02120 template <unsigned int N, class T, class StrideTag>
02121 template <int M, class Index>
02122 MultiArrayView <N - M, T, StridedArrayTag>
02123 MultiArrayView <N, T, StrideTag>::bindInner (const TinyVector <Index, M> &d) const
02124 {
02125     TinyVector <MultiArrayIndex, M> stride;
02126     stride.init (m_stride.begin (), m_stride.end () - N + M);
02127     pointer ptr = m_ptr + dot (d, stride);
02128     static const int NNew = (N-M == 0) ? 1 : N-M;
02129     TinyVector <MultiArrayIndex, NNew> outer_shape, outer_stride;
02130     if (N-M == 0)
02131     {
02132         outer_shape [0] = 1;
02133         outer_stride [0] = 0;
02134     }
02135     else
02136     {
02137         outer_shape.init (m_shape.begin () + M, m_shape.end ());
02138         outer_stride.init (m_stride.begin () + M, m_stride.end ());
02139     }
02140     return MultiArrayView <N-M, T, StridedArrayTag>
02141         (outer_shape, outer_stride, ptr);
02142 }
02143 
02144 template <unsigned int N, class T, class StrideTag>
02145 template <unsigned int M>
02146 MultiArrayView <N-1, T, typename detail::MaybeStrided<StrideTag, M>::type >
02147 MultiArrayView <N, T, StrideTag>::bind (difference_type_1 d) const
02148 {
02149     static const int NNew = (N-1 == 0) ? 1 : N-1;
02150     TinyVector <MultiArrayIndex, NNew> shape, stride;
02151     // the remaining dimensions are 0..n-1,n+1..N-1
02152     if (N-1 == 0)
02153     {
02154         shape[0] = 1;
02155         stride[0] = 0;
02156     }
02157     else
02158     {
02159         std::copy (m_shape.begin (), m_shape.begin () + M, shape.begin ());
02160         std::copy (m_shape.begin () + M+1, m_shape.end (),
02161                    shape.begin () + M);
02162         std::copy (m_stride.begin (), m_stride.begin () + M, stride.begin ());
02163         std::copy (m_stride.begin () + M+1, m_stride.end (),
02164                    stride.begin () + M);
02165     }
02166     return MultiArrayView <N-1, T, typename detail::MaybeStrided<StrideTag, M>::type>
02167         (shape, stride, m_ptr + d * m_stride[M]);
02168 }
02169 
02170 template <unsigned int N, class T, class StrideTag>
02171 MultiArrayView <N - 1, T, StrideTag>
02172 MultiArrayView <N, T, StrideTag>::bindOuter (difference_type_1 d) const
02173 {
02174     static const int NNew = (N-1 == 0) ? 1 : N-1;
02175     TinyVector <MultiArrayIndex, NNew> inner_shape, inner_stride;
02176     if (N-1 == 0)
02177     {
02178         inner_shape [0] = 1;
02179         inner_stride [0] = 0;
02180     }
02181     else
02182     {
02183         inner_shape.init (m_shape.begin (), m_shape.end () - 1);
02184         inner_stride.init (m_stride.begin (), m_stride.end () - 1);
02185     }
02186     return MultiArrayView <N-1, T, StrideTag> (inner_shape, inner_stride,
02187                                        m_ptr + d * m_stride [N-1]);
02188 }
02189 
02190 template <unsigned int N, class T, class StrideTag>
02191 MultiArrayView <N - 1, T, StridedArrayTag>
02192 MultiArrayView <N, T, StrideTag>::bindInner (difference_type_1 d) const
02193 {
02194     static const int NNew = (N-1 == 0) ? 1 : N-1;
02195     TinyVector <MultiArrayIndex, NNew> outer_shape, outer_stride;
02196     if (N-1 == 0)
02197     {
02198         outer_shape [0] = 1;
02199         outer_stride [0] = 0;
02200     }
02201     else
02202     {
02203         outer_shape.init (m_shape.begin () + 1, m_shape.end ());
02204         outer_stride.init (m_stride.begin () + 1, m_stride.end ());
02205     }
02206     return MultiArrayView <N-1, T, StridedArrayTag>
02207         (outer_shape, outer_stride, m_ptr + d * m_stride [0]);
02208 }
02209 
02210 template <unsigned int N, class T, class StrideTag>
02211 MultiArrayView <N - 1, T, StridedArrayTag>
02212 MultiArrayView <N, T, StrideTag>::bindAt (difference_type_1 n, difference_type_1 d) const
02213 {
02214     vigra_precondition (
02215         n < static_cast <int> (N),
02216         "MultiArrayView <N, T, StrideTag>::bindAt(): dimension out of range.");
02217     static const int NNew = (N-1 == 0) ? 1 : N-1;
02218     TinyVector <MultiArrayIndex, NNew> shape, stride;
02219     // the remaining dimensions are 0..n-1,n+1..N-1
02220     if (N-1 == 0)
02221     {
02222         shape [0] = 1;
02223         stride [0] = 0;
02224     }
02225     else
02226     {
02227         std::copy (m_shape.begin (), m_shape.begin () + n, shape.begin ());
02228         std::copy (m_shape.begin () + n+1, m_shape.end (),
02229                    shape.begin () + n);
02230         std::copy (m_stride.begin (), m_stride.begin () + n, stride.begin ());
02231         std::copy (m_stride.begin () + n+1, m_stride.end (),
02232                    stride.begin () + n);
02233     }
02234     return MultiArrayView <N-1, T, StridedArrayTag>
02235         (shape, stride, m_ptr + d * m_stride[n]);
02236 }
02237 
02238 
02239 template <unsigned int N, class T, class StrideTag>
02240 MultiArrayView <N+1, typename ExpandElementResult<T>::type, StridedArrayTag>
02241 MultiArrayView <N, T, StrideTag>::expandElements(difference_type_1 d) const
02242 {
02243     vigra_precondition(0 <= d && d <= static_cast <difference_type_1> (N),
02244           "MultiArrayView<N, ...>::expandElements(d): 0 <= 'd' <= N required.");
02245     
02246     int elementSize = ExpandElementResult<T>::size;
02247     typename MultiArrayShape<N+1>::type newShape, newStrides;
02248     for(int k=0; k<d; ++k)
02249     {
02250         newShape[k] = m_shape[k];
02251         newStrides[k] = m_stride[k]*elementSize;
02252     }   
02253     
02254     newShape[d] = elementSize;
02255     newStrides[d] = 1;
02256     
02257     for(int k=d; k<N; ++k)
02258     {
02259         newShape[k+1] = m_shape[k];
02260         newStrides[k+1] = m_stride[k]*elementSize;
02261     }   
02262     
02263     typedef typename ExpandElementResult<T>::type U;     
02264     return MultiArrayView<N+1, U, StridedArrayTag>(
02265                     newShape, newStrides, reinterpret_cast<U*>(m_ptr));
02266 }
02267 
02268 template <unsigned int N, class T, class StrideTag>
02269 MultiArrayView <N+1, T, StrideTag>
02270 MultiArrayView <N, T, StrideTag>::insertSingletonDimension (difference_type_1 i) const
02271 {
02272     vigra_precondition (
02273         0 <= i && i <= static_cast <difference_type_1> (N),
02274         "MultiArrayView <N, T, StrideTag>::insertSingletonDimension(): index out of range.");
02275     TinyVector <MultiArrayIndex, N+1> shape, stride;
02276     std::copy (m_shape.begin (), m_shape.begin () + i, shape.begin ());
02277     std::copy (m_shape.begin () + i, m_shape.end (), shape.begin () + i + 1);
02278     std::copy (m_stride.begin (), m_stride.begin () + i, stride.begin ());
02279     std::copy (m_stride.begin () + i, m_stride.end (), stride.begin () + i + 1);
02280     shape[i] = 1;
02281     stride[i] = 1;
02282 
02283     return MultiArrayView <N+1, T, StrideTag>(shape, stride, m_ptr);
02284 }
02285 
02286 template <unsigned int N, class T, class StrideTag>
02287 typename NormTraits<MultiArrayView <N, T, StrideTag> >::NormType
02288 MultiArrayView <N, T, StrideTag>::norm(int type, bool useSquaredNorm) const
02289 {
02290     typedef typename NormTraits<MultiArrayView>::NormType NormType;
02291 
02292     switch(type)
02293     {
02294       case 0:
02295       {
02296         NormType res = NumericTraits<NormType>::zero();
02297         detail::reduceOverMultiArray(traverser_begin(), shape(), 
02298                                      res, 
02299                                      detail::MaxNormReduceFunctor(),
02300                                      MetaInt<actual_dimension-1>());
02301         return res;
02302       }
02303       case 1:
02304       {
02305         NormType res = NumericTraits<NormType>::zero();
02306         detail::reduceOverMultiArray(traverser_begin(), shape(), 
02307                                      res, 
02308                                      detail::L1NormReduceFunctor(),
02309                                      MetaInt<actual_dimension-1>());
02310         return res;
02311       }
02312       case 2:
02313       {
02314         if(useSquaredNorm)
02315         {
02316             return sqrt((NormType)squaredNorm());
02317         }
02318         else
02319         {
02320             NormType normMax = NumericTraits<NormType>::zero();
02321             detail::reduceOverMultiArray(traverser_begin(), shape(), 
02322                                         normMax, 
02323                                         detail::MaxNormReduceFunctor(),
02324                                         MetaInt<actual_dimension-1>());
02325             if(normMax == NumericTraits<NormType>::zero())
02326                 return normMax;
02327             NormType res  = NumericTraits<NormType>::zero();
02328             detail::reduceOverMultiArray(traverser_begin(), shape(), 
02329                                          res, 
02330                                          detail::WeightedL2NormReduceFunctor<NormType>(1.0/normMax),
02331                                          MetaInt<actual_dimension-1>());
02332             return sqrt(res)*normMax;
02333         }
02334       }
02335       default:
02336         vigra_precondition(false, "MultiArrayView::norm(): Unknown norm type.");
02337         return NumericTraits<NormType>::zero(); // unreachable
02338     }
02339 }
02340 
02341 
02342 /********************************************************/
02343 /*                                                      */
02344 /*                          norm                        */
02345 /*                                                      */
02346 /********************************************************/
02347 
02348 template <unsigned int N, class T, class StrideTag>
02349 inline typename NormTraits<MultiArrayView <N, T, StrideTag> >::SquaredNormType
02350 squaredNorm(MultiArrayView <N, T, StrideTag> const & a)
02351 {
02352     return a.squaredNorm();
02353 }
02354 
02355 template <unsigned int N, class T, class StrideTag>
02356 inline typename NormTraits<MultiArrayView <N, T, StrideTag> >::NormType
02357 norm(MultiArrayView <N, T, StrideTag> const & a)
02358 {
02359     return a.norm();
02360 }
02361 
02362 /********************************************************/
02363 /*                                                      */
02364 /*                       MultiArray                     */
02365 /*                                                      */
02366 /********************************************************/
02367 
02368 /** \brief Main <TT>MultiArray</TT> class containing the memory
02369     management.
02370 
02371 This class inherits the interface of MultiArrayView, and implements
02372 the memory ownership.
02373 MultiArray's are always unstrided, striding them creates a MultiArrayView.
02374 
02375 
02376 The template parameters are as follows
02377 \code
02378     N: the array dimension
02379 
02380     T: the type of the array elements
02381 
02382     A: the allocator used for internal storage management
02383        (default: std::allocator<T>)
02384 \endcode
02385 
02386 <b>\#include</b>
02387 <vigra/multi_array.hxx>
02388 
02389 Namespace: vigra
02390 */
02391 template <unsigned int N, class T, class A /* default already declared above */>
02392 class MultiArray : public MultiArrayView <N, T>
02393 {
02394 
02395 public:
02396     using MultiArrayView <N, T>::actual_dimension;
02397 
02398         /** the allocator type used to allocate the memory
02399          */
02400     typedef A allocator_type;
02401 
02402         /** the view type associated with this array.
02403          */
02404     typedef MultiArrayView <N, T> view_type;
02405 
02406         /** the matrix type associated with this array.
02407          */
02408     typedef MultiArray <N, T, A> matrix_type;
02409 
02410         /** the array's value type
02411          */
02412     typedef typename view_type::value_type value_type;
02413 
02414         /** pointer type
02415          */
02416     typedef typename view_type::pointer pointer;
02417 
02418         /** const pointer type
02419          */
02420     typedef typename view_type::const_pointer const_pointer;
02421 
02422         /** reference type (result of operator[])
02423          */
02424     typedef typename view_type::reference reference;
02425 
02426         /** const reference type (result of operator[] const)
02427          */
02428     typedef typename view_type::const_reference const_reference;
02429 
02430         /** size type
02431          */
02432     typedef typename view_type::size_type size_type;
02433 
02434         /** difference type (used for multi-dimensional offsets and indices)
02435          */
02436     typedef typename view_type::difference_type difference_type;
02437 
02438         /** difference and index type for a single dimension
02439          */
02440     typedef typename view_type::difference_type_1 difference_type_1;
02441 
02442         /** traverser type
02443          */
02444     typedef typename vigra::detail::MultiIteratorChooser <
02445         UnstridedArrayTag>::template Traverser <N, T, T &, T *>::type
02446     traverser;
02447 
02448         /** traverser type to const data
02449          */
02450     typedef typename vigra::detail::MultiIteratorChooser <
02451         UnstridedArrayTag>::template Traverser <N, T, T const &, T const *>::type
02452     const_traverser;
02453 
02454         /** sequential (random access) iterator type
02455          */
02456     typedef T * iterator;
02457 
02458         /** sequential (random access) const iterator type
02459          */
02460     typedef T * const_iterator;
02461 
02462 protected:
02463 
02464     typedef typename difference_type::value_type diff_zero_t;
02465 
02466         /** the allocator used to allocate the memory
02467          */
02468     allocator_type m_alloc;
02469 
02470         /** allocate memory for s pixels, write its address into the given
02471             pointer and initialize the pixels with init.
02472         */
02473     void allocate (pointer &ptr, difference_type_1 s, const_reference init);
02474 
02475         /** allocate memory for s pixels, write its address into the given
02476             pointer and initialize the linearized pixels to the values of init.
02477         */
02478     template <class U>
02479     void allocate (pointer &ptr, difference_type_1 s, U const * init);
02480 
02481         /** allocate memory, write its address into the given
02482             pointer and initialize it by copying the data from the given MultiArrayView.
02483         */
02484     template <class U, class StrideTag>
02485     void allocate (pointer &ptr, MultiArrayView<N, U, StrideTag> const & init);
02486 
02487         /** deallocate the memory (of length s) starting at the given address.
02488          */
02489     void deallocate (pointer &ptr, difference_type_1 s);
02490 
02491     template <class U, class StrideTag>
02492     void copyOrReshape (const MultiArrayView<N, U, StrideTag> &rhs);
02493 public:
02494         /** default constructor
02495          */
02496     MultiArray ()
02497     : MultiArrayView <N, T> (difference_type (diff_zero_t(0)),
02498                              difference_type (diff_zero_t(0)), 0)
02499     {}
02500 
02501         /** construct with given allocator
02502          */
02503     MultiArray (allocator_type const & alloc)
02504     : MultiArrayView <N, T> (difference_type (diff_zero_t(0)),
02505                              difference_type (diff_zero_t(0)), 0),
02506       m_alloc(alloc)
02507     {}
02508 
02509         /** construct with given length
02510         
02511             Use only for 1-dimensional arrays (<tt>N==1</tt>).
02512          */
02513     explicit MultiArray (difference_type_1 length,
02514                          allocator_type const & alloc = allocator_type());
02515 
02516         /** construct with given shape
02517          */
02518     explicit MultiArray (const difference_type &shape,
02519                          allocator_type const & alloc = allocator_type());
02520 
02521         /** construct from shape with an initial value
02522          */
02523     MultiArray (const difference_type &shape, const_reference init,
02524                 allocator_type const & alloc = allocator_type());
02525 
02526         /** construct from shape and copy values from the given array
02527          */
02528     MultiArray (const difference_type &shape, const_pointer init,
02529                          allocator_type const & alloc = allocator_type());
02530 
02531         /** copy constructor
02532          */
02533     MultiArray (const MultiArray &rhs)
02534     : MultiArrayView <N, T> (rhs.m_shape, rhs.m_stride, 0),
02535       m_alloc (rhs.m_alloc)
02536     {
02537         allocate (this->m_ptr, this->elementCount (), rhs.data ());
02538     }
02539 
02540         /** constructor from an array expression
02541          */
02542     template<class Expression>
02543     MultiArray (multi_math::MultiMathOperand<Expression> const & rhs,
02544                 allocator_type const & alloc = allocator_type())
02545     : MultiArrayView <N, T> (difference_type (diff_zero_t(0)),
02546                              difference_type (diff_zero_t(0)), 0),
02547       m_alloc (alloc)
02548     {
02549         multi_math::detail::assignOrResize(*this, rhs);
02550     }
02551 
02552         /** construct by copying from a MultiArrayView
02553          */
02554     template <class U, class StrideTag>
02555     MultiArray (const MultiArrayView<N, U, StrideTag>  &rhs,
02556                 allocator_type const & alloc = allocator_type());
02557 
02558         /** assignment.<br>
02559             If the size of \a rhs is the same as the left-hand side arrays's old size, only
02560             the data are copied. Otherwise, new storage is allocated, which invalidates all
02561             objects (array views, iterators) depending on the lhs array.
02562          */
02563     MultiArray & operator= (const MultiArray &rhs)
02564     {
02565         if (this != &rhs)
02566             this->copyOrReshape(rhs);
02567         return *this;
02568     }
02569 
02570         /** assignment from arbitrary MultiArrayView.<br>
02571             If the size of \a rhs is the same as the left-hand side arrays's old size, only
02572             the data are copied. Otherwise, new storage is allocated, which invalidates all
02573             objects (array views, iterators) depending on the lhs array.
02574          */
02575     template <class U, class StrideTag>
02576     MultiArray &operator= (const MultiArrayView<N, U, StrideTag> &rhs)
02577     {
02578         this->copyOrReshape(rhs);
02579         return *this;
02580     }
02581 
02582         /** assignment from scalar.<br>
02583             Equivalent to MultiArray::init(v).
02584          */
02585     MultiArray & operator=(value_type const & v)
02586     {
02587         return this->init(v);
02588     }
02589 
02590         /** Add-assignment from arbitrary MultiArrayView. Fails with
02591             <tt>PreconditionViolation</tt> exception when the shapes do not match.
02592             If the left array has no data (hasData() is false), this function is 
02593             equivalent to a normal assignment (i.e. an empty
02594             array is interpreted as a zero-array of appropriate size).
02595          */
02596     template <class U, class StrideTag>
02597     MultiArray &operator+= (const MultiArrayView<N, U, StrideTag> &rhs)
02598     {
02599         if(this->hasData())
02600             view_type::operator+=(rhs);
02601         else
02602             *this = rhs;
02603         return *this;
02604     }
02605 
02606         /** Subtract-assignment from arbitrary MultiArrayView. Fails with
02607             <tt>PreconditionViolation</tt> exception when the shapes do not match.
02608             If the left array has no data (hasData() is false), this function is 
02609             equivalent to an assignment of the negated rhs (i.e. an empty
02610             array is interpreted as a zero-array of appropriate size).
02611          */
02612     template <class U, class StrideTag>
02613     MultiArray &operator-= (const MultiArrayView<N, U, StrideTag> &rhs)
02614     {
02615         if(!this->hasData())
02616             this->reshape(rhs.shape());
02617         view_type::operator-=(rhs);
02618         return *this;
02619     }
02620 
02621         /** Multiply-assignment from arbitrary MultiArrayView. Fails with
02622             <tt>PreconditionViolation</tt> exception when the shapes do not match.
02623             If the left array has no data (hasData() is false), this function is 
02624             equivalent to reshape(rhs.shape()) with zero initialisation (i.e. an empty
02625             array is interpreted as a zero-array of appropriate size).
02626          */
02627     template <class U, class StrideTag>
02628     MultiArray &operator*= (const MultiArrayView<N, U, StrideTag> &rhs)
02629     {
02630         if(this->hasData())
02631             view_type::operator*=(rhs);
02632         else
02633             this->reshape(rhs.shape());
02634         return *this;
02635     }
02636 
02637         /** Divide-assignment from arbitrary MultiArrayView. Fails with
02638             <tt>PreconditionViolation</tt> exception when the shapes do not match.
02639             If the left array has no data (hasData() is false), this function is 
02640             equivalent to reshape(rhs.shape()) with zero initialisation (i.e. an empty
02641             array is interpreted as a zero-array of appropriate size).
02642          */
02643     template <class U, class StrideTag>
02644     MultiArray &operator/= (const MultiArrayView<N, U, StrideTag> &rhs)
02645     {
02646         if(this->hasData())
02647             view_type::operator/=(rhs);
02648         else
02649             this->reshape(rhs.shape());
02650         return *this;
02651     }
02652 
02653         /** Add-assignment of a scalar.
02654          */
02655     MultiArray &operator+= (const T &rhs)
02656     {
02657         view_type::operator+=(rhs);
02658         return *this;
02659     }
02660 
02661         /** Subtract-assignment of a scalar.
02662          */
02663     MultiArray &operator-= (const T &rhs)
02664     {
02665         view_type::operator-=(rhs);
02666         return *this;
02667     }
02668 
02669         /** Multiply-assignment of a scalar.
02670          */
02671     MultiArray &operator*= (const T &rhs)
02672     {
02673         view_type::operator*=(rhs);
02674         return *this;
02675     }
02676 
02677         /** Divide-assignment of a scalar.
02678          */
02679     MultiArray &operator/= (const T &rhs)
02680     {
02681         view_type::operator/=(rhs);
02682         return *this;
02683     }
02684         /** Assignment of an array expression. Fails with
02685             <tt>PreconditionViolation</tt> exception when the shapes do not match.
02686          */
02687     template<class Expression>
02688     MultiArray & operator=(multi_math::MultiMathOperand<Expression> const & rhs)
02689     {
02690         multi_math::detail::assignOrResize(*this, rhs);
02691         return *this;
02692     }
02693 
02694         /** Add-assignment of an array expression. Fails with
02695             <tt>PreconditionViolation</tt> exception when the shapes do not match.
02696          */
02697     template<class Expression>
02698     MultiArray & operator+=(multi_math::MultiMathOperand<Expression> const & rhs)
02699     {
02700         multi_math::detail::plusAssignOrResize(*this, rhs);
02701         return *this;
02702     }
02703 
02704         /** Subtract-assignment of an array expression. Fails with
02705             <tt>PreconditionViolation</tt> exception when the shapes do not match.
02706          */
02707     template<class Expression>
02708     MultiArray & operator-=(multi_math::MultiMathOperand<Expression> const & rhs)
02709     {
02710         multi_math::detail::minusAssignOrResize(*this, rhs);
02711         return *this;
02712     }
02713 
02714         /** Multiply-assignment of an array expression. Fails with
02715             <tt>PreconditionViolation</tt> exception when the shapes do not match.
02716          */
02717     template<class Expression>
02718     MultiArray & operator*=(multi_math::MultiMathOperand<Expression> const & rhs)
02719     {
02720         multi_math::detail::multiplyAssignOrResize(*this, rhs);
02721         return *this;
02722     }
02723 
02724         /** Divide-assignment of an array expression. Fails with
02725             <tt>PreconditionViolation</tt> exception when the shapes do not match.
02726          */
02727     template<class Expression>
02728     MultiArray & operator/=(multi_math::MultiMathOperand<Expression> const & rhs)
02729     {
02730         multi_math::detail::divideAssignOrResize(*this, rhs);
02731         return *this;
02732     }
02733 
02734         /** destructor
02735          */
02736     ~MultiArray ()
02737     {
02738         deallocate (this->m_ptr, this->elementCount ());
02739     }
02740 
02741 
02742         /** init elements with a constant
02743          */
02744     template <class U>
02745     MultiArray & init(const U & init)
02746     {
02747         view_type::init(init);
02748         return *this;
02749     }
02750 
02751         /** Allocate new memory with the given shape and initialize with zeros.<br>
02752             <em>Note:</em> this operation invalidates all dependent objects
02753             (array views and iterators)
02754          */
02755     void reshape (const difference_type &shape)
02756     {
02757         reshape (shape, T());
02758     }
02759 
02760         /** Allocate new memory with the given shape and initialize it
02761             with the given value.<br>
02762             <em>Note:</em> this operation invalidates all dependent objects
02763             (array views and iterators)
02764          */
02765     void reshape (const difference_type &shape, const_reference init);
02766 
02767         /** Swap the contents with another MultiArray. This is fast,
02768             because no data are copied, but only pointers and shapes swapped.
02769             <em>Note:</em> this operation invalidates all dependent objects
02770             (array views and iterators)
02771          */
02772     void swap (MultiArray & other);
02773 
02774         /** sequential iterator pointing to the first array element.
02775          */
02776     iterator begin ()
02777     {
02778         return this->data();
02779     }
02780 
02781         /** sequential iterator pointing beyond the last array element.
02782          */
02783     iterator end ()
02784     {
02785         return this->data() + this->elementCount();
02786     }
02787 
02788         /** sequential const iterator pointing to the first array element.
02789          */
02790     const_iterator begin () const
02791     {
02792         return this->data();
02793     }
02794 
02795         /** sequential const iterator pointing beyond the last array element.
02796          */
02797     const_iterator end () const
02798     {
02799         return this->data() + this->elementCount();
02800     }
02801 
02802         /** get the allocator.
02803          */
02804     allocator_type const & allocator () const
02805     {
02806         return m_alloc;
02807     }
02808 };
02809 
02810 template <unsigned int N, class T, class A>
02811 MultiArray <N, T, A>::MultiArray (difference_type_1 length,
02812                                   allocator_type const & alloc)
02813 : MultiArrayView <N, T> (difference_type(length),
02814                          detail::defaultStride <1> (difference_type(length)),
02815                          0),
02816   m_alloc(alloc)
02817 {
02818     allocate (this->m_ptr, this->elementCount (), T());
02819 }
02820 
02821 template <unsigned int N, class T, class A>
02822 MultiArray <N, T, A>::MultiArray (const difference_type &shape,
02823                                   allocator_type const & alloc)
02824 : MultiArrayView <N, T> (shape,
02825                          detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
02826                          0),
02827   m_alloc(alloc)
02828 {
02829     if (N == 0)
02830     {
02831         this->m_shape [0] = 1;
02832         this->m_stride [0] = 0;
02833     }
02834     allocate (this->m_ptr, this->elementCount (), T());
02835 }
02836 
02837 template <unsigned int N, class T, class A>
02838 MultiArray <N, T, A>::MultiArray (const difference_type &shape, const_reference init,
02839                                   allocator_type const & alloc)
02840 : MultiArrayView <N, T> (shape,
02841                          detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
02842                          0),
02843   m_alloc(alloc)
02844 {
02845     if (N == 0)
02846     {
02847         this->m_shape [0] = 1;
02848         this->m_stride [0] = 0;
02849     }
02850     allocate (this->m_ptr, this->elementCount (), init);
02851 }
02852 
02853 template <unsigned int N, class T, class A>
02854 MultiArray <N, T, A>::MultiArray (const difference_type &shape, const_pointer init,
02855                                   allocator_type const & alloc)
02856 : MultiArrayView <N, T> (shape,
02857                          detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
02858                          0),
02859   m_alloc(alloc)
02860 {
02861     if (N == 0)
02862     {
02863         this->m_shape [0] = 1;
02864         this->m_stride [0] = 0;
02865     }
02866     allocate (this->m_ptr, this->elementCount (), init);
02867 }
02868 
02869 template <unsigned int N, class T, class A>
02870 template <class U, class StrideTag>
02871 MultiArray <N, T, A>::MultiArray(const MultiArrayView<N, U, StrideTag>  &rhs,
02872                                  allocator_type const & alloc)
02873 : MultiArrayView <N, T> (rhs.shape(),
02874                          detail::defaultStride <MultiArrayView<N,T>::actual_dimension>(rhs.shape()),
02875                          0),
02876   m_alloc (alloc)
02877 {
02878     allocate (this->m_ptr, rhs);
02879 }
02880 
02881 template <unsigned int N, class T, class A>
02882 template <class U, class StrideTag>
02883 void
02884 MultiArray <N, T, A>::copyOrReshape(const MultiArrayView<N, U, StrideTag> &rhs)
02885 {
02886     if (this->shape() == rhs.shape())
02887         this->copy(rhs);
02888     else
02889     {
02890         MultiArray t(rhs);
02891         this->swap(t);
02892     }
02893 }
02894 
02895 template <unsigned int N, class T, class A>
02896 void MultiArray <N, T, A>::reshape (const difference_type & new_shape,
02897                                     const_reference initial)
02898 {
02899     if (N== 0)
02900     {
02901         return;
02902     }
02903     else if(new_shape == this->shape())
02904     {
02905         this->init(initial);
02906     }
02907     else
02908     {
02909         difference_type new_stride = detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (new_shape);
02910         difference_type_1 new_size = new_shape [MultiArrayView<N,T>::actual_dimension-1] * new_stride [MultiArrayView<N,T>::actual_dimension-1];
02911         T *new_ptr;
02912         allocate (new_ptr, new_size, initial);
02913         deallocate (this->m_ptr, this->elementCount ());
02914         this->m_ptr = new_ptr;
02915         this->m_shape = new_shape;
02916         this->m_stride = new_stride;
02917     }
02918 }
02919 
02920 
02921 template <unsigned int N, class T, class A>
02922 inline void
02923 MultiArray <N, T, A>::swap (MultiArray & other)
02924 {
02925     if (this == &other)
02926         return;
02927     std::swap(this->m_shape,  other.m_shape);
02928     std::swap(this->m_stride, other.m_stride);
02929     std::swap(this->m_ptr,    other.m_ptr);
02930     std::swap(this->m_alloc,  other.m_alloc);
02931 }
02932 
02933 template <unsigned int N, class T, class A>
02934 void MultiArray <N, T, A>::allocate (pointer & ptr, difference_type_1 s,
02935                                      const_reference init)
02936 {
02937     ptr = m_alloc.allocate ((typename A::size_type)s);
02938     difference_type_1 i;
02939     try {
02940         for (i = 0; i < s; ++i)
02941             m_alloc.construct (ptr + i, init);
02942     }
02943     catch (...) {
02944         for (difference_type_1 j = 0; j < i; ++j)
02945             m_alloc.destroy (ptr + j);
02946         m_alloc.deallocate (ptr, (typename A::size_type)s);
02947         throw;
02948     }
02949 }
02950 
02951 template <unsigned int N, class T, class A>
02952 template <class U>
02953 void MultiArray <N, T, A>::allocate (pointer & ptr, difference_type_1 s,
02954                                      U const * init)
02955 {
02956     ptr = m_alloc.allocate ((typename A::size_type)s);
02957     difference_type_1 i;
02958     try {
02959         for (i = 0; i < s; ++i, ++init)
02960             m_alloc.construct (ptr + i, *init);
02961     }
02962     catch (...) {
02963         for (difference_type_1 j = 0; j < i; ++j)
02964             m_alloc.destroy (ptr + j);
02965         m_alloc.deallocate (ptr, (typename A::size_type)s);
02966         throw;
02967     }
02968 }
02969 
02970 template <unsigned int N, class T, class A>
02971 template <class U, class StrideTag>
02972 void MultiArray <N, T, A>::allocate (pointer & ptr, MultiArrayView<N, U, StrideTag> const & init)
02973 {
02974     difference_type_1 s = init.elementCount();
02975     ptr = m_alloc.allocate ((typename A::size_type)s);
02976     pointer p = ptr;
02977     try {
02978         detail::uninitializedCopyMultiArrayData(init.traverser_begin(), init.shape(),
02979                                                 p, m_alloc, MetaInt<actual_dimension-1>());
02980     }
02981     catch (...) {
02982         for (pointer pp = ptr; pp < p; ++pp)
02983             m_alloc.destroy (pp);
02984         m_alloc.deallocate (ptr, (typename A::size_type)s);
02985         throw;
02986     }
02987 }
02988 
02989 template <unsigned int N, class T, class A>
02990 inline void MultiArray <N, T, A>::deallocate (pointer & ptr, difference_type_1 s)
02991 {
02992     if (ptr == 0)
02993         return;
02994     for (difference_type_1 i = 0; i < s; ++i)
02995         m_alloc.destroy (ptr + i);
02996     m_alloc.deallocate (ptr, (typename A::size_type)s);
02997     ptr = 0;
02998 }
02999 
03000 /********************************************************/
03001 /*                                                      */
03002 /*              argument object factories               */
03003 /*                                                      */
03004 /********************************************************/
03005 
03006 template <unsigned int N, class T, class StrideTag>
03007 inline triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
03008               typename MultiArrayView<N,T,StrideTag>::difference_type,
03009               typename AccessorTraits<T>::default_const_accessor >
03010 srcMultiArrayRange( MultiArrayView<N,T,StrideTag> const & array )
03011 {
03012     return triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
03013                   typename MultiArrayView<N,T,StrideTag>::difference_type,
03014                   typename AccessorTraits<T>::default_const_accessor >
03015       ( array.traverser_begin(),
03016         array.shape(),
03017         typename AccessorTraits<T>::default_const_accessor() );
03018 }
03019 
03020 template <unsigned int N, class T, class StrideTag, class Accessor>
03021 inline triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
03022               typename MultiArrayView<N,T,StrideTag>::difference_type,
03023               Accessor >
03024 srcMultiArrayRange( MultiArrayView<N,T,StrideTag> const & array, Accessor a )
03025 {
03026     return triple<typename MultiArrayView<N,T,StrideTag>::const_traverser,
03027                   typename MultiArrayView<N,T,StrideTag>::difference_type,
03028                   Accessor >
03029       ( array.traverser_begin(),
03030         array.shape(),
03031         a);
03032 }
03033 
03034 template <unsigned int N, class T, class StrideTag>
03035 inline pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
03036             typename AccessorTraits<T>::default_const_accessor >
03037 srcMultiArray( MultiArrayView<N,T,StrideTag> const & array )
03038 {
03039     return pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
03040                 typename AccessorTraits<T>::default_const_accessor >
03041       ( array.traverser_begin(),
03042         typename AccessorTraits<T>::default_const_accessor() );
03043 }
03044 
03045 template <unsigned int N, class T, class StrideTag, class Accessor>
03046 inline pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
03047             Accessor >
03048 srcMultiArray( MultiArrayView<N,T,StrideTag> const & array, Accessor a )
03049 {
03050     return pair<typename MultiArrayView<N,T,StrideTag>::const_traverser,
03051                 Accessor >
03052       ( array.traverser_begin(), a );
03053 }
03054 
03055 template <unsigned int N, class T, class StrideTag>
03056 inline triple<typename MultiArrayView<N,T,StrideTag>::traverser,
03057               typename MultiArrayView<N,T,StrideTag>::difference_type,
03058               typename AccessorTraits<T>::default_accessor >
03059 destMultiArrayRange( MultiArrayView<N,T,StrideTag> & array )
03060 {
03061     return triple<typename MultiArrayView<N,T,StrideTag>::traverser,
03062                   typename MultiArrayView<N,T,StrideTag>::difference_type,
03063                   typename AccessorTraits<T>::default_accessor >
03064       ( array.traverser_begin(),
03065         array.shape(),
03066         typename AccessorTraits<T>::default_accessor() );
03067 }
03068 
03069 template <unsigned int N, class T, class StrideTag, class Accessor>
03070 inline triple<typename MultiArrayView<N,T,StrideTag>::traverser,
03071               typename MultiArrayView<N,T,StrideTag>::difference_type,
03072               Accessor >
03073 destMultiArrayRange( MultiArrayView<N,T,StrideTag> & array, Accessor a )
03074 {
03075     return triple<typename MultiArrayView<N,T,StrideTag>::traverser,
03076                   typename MultiArrayView<N,T,StrideTag>::difference_type,
03077                   Accessor >
03078       ( array.traverser_begin(),
03079         array.shape(),
03080         a );
03081 }
03082 
03083 template <unsigned int N, class T, class StrideTag>
03084 inline pair<typename MultiArrayView<N,T,StrideTag>::traverser,
03085             typename AccessorTraits<T>::default_accessor >
03086 destMultiArray( MultiArrayView<N,T,StrideTag> & array )
03087 {
03088     return pair<typename MultiArrayView<N,T,StrideTag>::traverser,
03089                 typename AccessorTraits<T>::default_accessor >
03090         ( array.traverser_begin(),
03091           typename AccessorTraits<T>::default_accessor() );
03092 }
03093 
03094 template <unsigned int N, class T, class StrideTag, class Accessor>
03095 inline pair<typename MultiArrayView<N,T,StrideTag>::traverser,
03096             Accessor >
03097 destMultiArray( MultiArrayView<N,T,StrideTag> & array, Accessor a )
03098 {
03099     return pair<typename MultiArrayView<N,T,StrideTag>::traverser,
03100                 Accessor >
03101         ( array.traverser_begin(), a );
03102 }
03103 
03104 /********************************************************************/
03105 
03106 template <class PixelType, class Accessor>
03107 inline triple<ConstStridedImageIterator<PixelType>,
03108               ConstStridedImageIterator<PixelType>, Accessor>
03109 srcImageRange(const MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
03110 {
03111     ConstStridedImageIterator<PixelType>
03112         ul(img.data(), 1, img.stride(0), img.stride(1));
03113     return triple<ConstStridedImageIterator<PixelType>,
03114                   ConstStridedImageIterator<PixelType>,
03115                   Accessor>(
03116                       ul, ul + Size2D(img.shape(0), img.shape(1)), a);
03117 }
03118 
03119 template <class PixelType, class Accessor>
03120 inline pair<ConstStridedImageIterator<PixelType>, Accessor>
03121 srcImage(const MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
03122 {
03123     ConstStridedImageIterator<PixelType>
03124         ul(img.data(), 1, img.stride(0), img.stride(1));
03125     return pair<ConstStridedImageIterator<PixelType>, Accessor>
03126         (ul, a);
03127 }
03128 
03129 template <class PixelType, class Accessor>
03130 inline triple<StridedImageIterator<PixelType>,
03131               StridedImageIterator<PixelType>, Accessor>
03132 destImageRange(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
03133 {
03134     StridedImageIterator<PixelType>
03135         ul(img.data(), 1, img.stride(0), img.stride(1));
03136     return triple<StridedImageIterator<PixelType>,
03137                   StridedImageIterator<PixelType>,
03138                   Accessor>(
03139                       ul, ul + Size2D(img.shape(0), img.shape(1)), a);
03140 }
03141 
03142 template <class PixelType, class Accessor>
03143 inline pair<StridedImageIterator<PixelType>, Accessor>
03144 destImage(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
03145 {
03146     StridedImageIterator<PixelType>
03147         ul(img.data(), 1, img.stride(0), img.stride(1));
03148     return pair<StridedImageIterator<PixelType>, Accessor>
03149         (ul, a);
03150 }
03151 
03152 template <class PixelType, class Accessor>
03153 inline pair<StridedImageIterator<PixelType>, Accessor>
03154 maskImage(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
03155 {
03156     StridedImageIterator<PixelType>
03157         ul(img.data(), 1, img.stride(0), img.stride(1));
03158     return pair<StridedImageIterator<PixelType>, Accessor>
03159         (ul, a);
03160 }
03161 
03162 // -------------------------------------------------------------------
03163 
03164 template <class PixelType>
03165 inline triple<ConstStridedImageIterator<PixelType>,
03166               ConstStridedImageIterator<PixelType>,
03167               typename AccessorTraits<PixelType>::default_const_accessor>
03168 srcImageRange(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
03169 {
03170     ConstStridedImageIterator<PixelType>
03171         ul(img.data(), 1, img.stride(0), img.stride(1));
03172     typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
03173     return triple<ConstStridedImageIterator<PixelType>,
03174                   ConstStridedImageIterator<PixelType>,
03175                   Accessor>
03176         (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
03177 }
03178 
03179 template <class PixelType>
03180 inline triple<ConstImageIterator<PixelType>,
03181               ConstImageIterator<PixelType>,
03182               typename AccessorTraits<PixelType>::default_const_accessor>
03183 srcImageRange(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
03184 {
03185     ConstImageIterator<PixelType>
03186         ul(img.data(), img.stride(1));
03187     typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
03188     return triple<ConstImageIterator<PixelType>,
03189                   ConstImageIterator<PixelType>,
03190                   Accessor>
03191         (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
03192 }
03193 
03194 template <class PixelType>
03195 inline pair< ConstStridedImageIterator<PixelType>,
03196              typename AccessorTraits<PixelType>::default_const_accessor>
03197 srcImage(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
03198 {
03199     ConstStridedImageIterator<PixelType>
03200         ul(img.data(), 1, img.stride(0), img.stride(1));
03201     typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
03202     return pair<ConstStridedImageIterator<PixelType>,
03203                 Accessor>
03204         (ul, Accessor());
03205 }
03206 
03207 template <class PixelType>
03208 inline pair< ConstImageIterator<PixelType>,
03209              typename AccessorTraits<PixelType>::default_const_accessor>
03210 srcImage(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
03211 {
03212     ConstImageIterator<PixelType>
03213         ul(img.data(), img.stride(1));
03214     typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
03215     return pair<ConstImageIterator<PixelType>,
03216                 Accessor>
03217         (ul, Accessor());
03218 }
03219 
03220 template <class PixelType>
03221 inline triple< StridedImageIterator<PixelType>,
03222                StridedImageIterator<PixelType>,
03223                typename AccessorTraits<PixelType>::default_accessor>
03224 destImageRange(MultiArrayView<2, PixelType, StridedArrayTag> & img)
03225 {
03226     StridedImageIterator<PixelType>
03227         ul(img.data(), 1, img.stride(0), img.stride(1));
03228     typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
03229     return triple<StridedImageIterator<PixelType>,
03230                   StridedImageIterator<PixelType>,
03231                   Accessor>
03232         (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
03233 }
03234 
03235 template <class PixelType>
03236 inline triple< ImageIterator<PixelType>,
03237                ImageIterator<PixelType>,
03238                typename AccessorTraits<PixelType>::default_accessor>
03239 destImageRange(MultiArrayView<2, PixelType, UnstridedArrayTag> & img)
03240 {
03241     ImageIterator<PixelType>
03242         ul(img.data(), img.stride(1));
03243     typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
03244     return triple<ImageIterator<PixelType>,
03245                   ImageIterator<PixelType>,
03246                   Accessor>
03247         (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
03248 }
03249 
03250 template <class PixelType>
03251 inline pair< StridedImageIterator<PixelType>,
03252              typename AccessorTraits<PixelType>::default_accessor>
03253 destImage(MultiArrayView<2, PixelType, StridedArrayTag> & img)
03254 {
03255     StridedImageIterator<PixelType>
03256         ul(img.data(), 1, img.stride(0), img.stride(1));
03257     typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
03258     return pair<StridedImageIterator<PixelType>, Accessor>
03259         (ul, Accessor());
03260 }
03261 
03262 template <class PixelType>
03263 inline pair< ImageIterator<PixelType>,
03264              typename AccessorTraits<PixelType>::default_accessor>
03265 destImage(MultiArrayView<2, PixelType, UnstridedArrayTag> & img)
03266 {
03267     ImageIterator<PixelType> ul(img.data(), img.stride(1));
03268     typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
03269     return pair<ImageIterator<PixelType>, Accessor>(ul, Accessor());
03270 }
03271 
03272 template <class PixelType>
03273 inline pair< ConstStridedImageIterator<PixelType>,
03274              typename AccessorTraits<PixelType>::default_accessor>
03275 maskImage(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
03276 {
03277     ConstStridedImageIterator<PixelType>
03278         ul(img.data(), 1, img.stride(0), img.stride(1));
03279     typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
03280     return pair<ConstStridedImageIterator<PixelType>, Accessor>
03281         (ul, Accessor());
03282 }
03283 
03284 template <class PixelType>
03285 inline pair< ConstImageIterator<PixelType>,
03286              typename AccessorTraits<PixelType>::default_accessor>
03287 maskImage(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
03288 {
03289     ConstImageIterator<PixelType>
03290         ul(img.data(), img.stride(1));
03291     typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
03292     return pair<ConstImageIterator<PixelType>, Accessor>
03293         (ul, Accessor());
03294 }
03295 
03296 /********************************************************/
03297 /*                                                      */
03298 /*                  makeBasicImageView                  */
03299 /*                                                      */
03300 /********************************************************/
03301 
03302 /** \addtogroup MultiArrayToImage Wrap a \ref vigra::MultiArrayView in
03303                                   a \ref vigra::BasicImageView
03304 */
03305 //@{
03306 /** Create a \ref vigra::BasicImageView from an unstrided 2-dimensional
03307     \ref vigra::MultiArrayView.
03308 
03309     The \ref vigra::BasicImageView will have the same <tt>value_type </tt>
03310     as the original \ref vigra::MultiArrayView.
03311 */
03312 template <class T>
03313 BasicImageView <T>
03314 makeBasicImageView (MultiArrayView <2, T, UnstridedArrayTag> const &array)
03315 {
03316     return BasicImageView <T> (array.data (), array.shape (0),
03317                                array.shape (1));
03318 }
03319 
03320 /** Create a \ref vigra::BasicImageView from a 3-dimensional
03321     \ref vigra::MultiArray.
03322 
03323     This wrapper flattens the two innermost dimensions of the array
03324     into single rows of the resulting image.
03325     The \ref vigra::BasicImageView will have the same <tt>value_type </tt>
03326     as the original \ref vigra::MultiArray.
03327 */
03328 template <class T>
03329 BasicImageView <T>
03330 makeBasicImageView (MultiArray <3, T> const &array)
03331 {
03332     return BasicImageView <T> (array.data (),
03333                                array.shape (0)*array.shape (1), array.shape (2));
03334 }
03335 
03336 /** Create a \ref vigra::BasicImageView from a 3-dimensional
03337     \ref vigra::MultiArray.
03338 
03339     This wrapper only works if <tt>T</tt> is a scalar type and the
03340     array's innermost dimension has size 3. It then re-interprets
03341     the data array as a 2-dimensional array with value_type
03342     <tt>RGBValue<T></tt>.
03343 */
03344 template <class T>
03345 BasicImageView <RGBValue<T> >
03346 makeRGBImageView (MultiArray<3, T> const &array)
03347 {
03348     vigra_precondition (
03349         array.shape (0) == 3, "makeRGBImageView(): array.shape(0) must be 3.");
03350     return BasicImageView <RGBValue<T> > (
03351         reinterpret_cast <RGBValue <T> *> (array.data ()),
03352         array.shape (1), array.shape (2));
03353 }
03354 
03355 //@}
03356 
03357 } // namespace vigra
03358 
03359 #undef VIGRA_ASSERT_INSIDE
03360 
03361 #endif // VIGRA_MULTI_ARRAY_HXX

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

html generated using doxygen and Python
vigra 1.9.0 (Tue Nov 6 2012)