vigra - vigra/multi_array.hxx Source File

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

vigra/multi_array.hxx |

    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)