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

vigra/multi_array.hxx VIGRA

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

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

html generated using doxygen and Python
vigra 1.8.0 (20 Sep 2011)