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

vigra/multi_array.hxx
00001 /************************************************************************/
00002 /*                                                                      */
00003 /*     Copyright 2003-2008 by Gunnar Kedenburg and Ullrich Koethe       */
00004 /*                                                                      */
00005 /*    This file is part of the VIGRA computer vision library.           */
00006 /*    The VIGRA Website is                                              */
00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
00008 /*    Please direct questions, bug reports, and contributions to        */
00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00010 /*        vigra@informatik.uni-hamburg.de                               */
00011 /*                                                                      */
00012 /*    Permission is hereby granted, free of charge, to any person       */
00013 /*    obtaining a copy of this software and associated documentation    */
00014 /*    files (the "Software"), to deal in the Software without           */
00015 /*    restriction, including without limitation the rights to use,      */
00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00017 /*    sell copies of the Software, and to permit persons to whom the    */
00018 /*    Software is furnished to do so, subject to the following          */
00019 /*    conditions:                                                       */
00020 /*                                                                      */
00021 /*    The above copyright notice and this permission notice shall be    */
00022 /*    included in all copies or substantial portions of the             */
00023 /*    Software.                                                         */
00024 /*                                                                      */
00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00033 /*                                                                      */
00034 /************************************************************************/
00035 
00036 
00037 #ifndef VIGRA_MULTI_ARRAY_HXX
00038 #define VIGRA_MULTI_ARRAY_HXX
00039 
00040 #include <memory>
00041 #include <algorithm>
00042 #include "accessor.hxx"
00043 #include "tinyvector.hxx"
00044 #include "rgbvalue.hxx"
00045 #include "basicimageview.hxx"
00046 #include "imageiterator.hxx"
00047 #include "numerictraits.hxx"
00048 #include "multi_iterator.hxx"
00049 #include "metaprogramming.hxx"
00050 #include "mathutil.hxx"
00051 
00052 // Bounds checking Macro used if VIGRA_CHECK_BOUNDS is defined.
00053 #ifdef VIGRA_CHECK_BOUNDS
00054 #define VIGRA_ASSERT_INSIDE(diff) \
00055   vigra_precondition(this->isInside(diff), "Index out of bounds")
00056 #else
00057 #define VIGRA_ASSERT_INSIDE(diff)
00058 #endif
00059 
00060 namespace vigra
00061 {
00062 
00063 namespace detail
00064 {
00065 /********************************************************/
00066 /*                                                      */
00067 /*                    defaultStride                     */
00068 /*                                                      */
00069 /********************************************************/
00070 
00071 /* generates the stride for a gapless shape.
00072 
00073     Namespace: vigra::detail
00074 */
00075 template <unsigned int N>
00076 inline TinyVector <MultiArrayIndex, N>
00077 defaultStride(const TinyVector <MultiArrayIndex, N> &shape)
00078 {
00079     TinyVector <MultiArrayIndex, N> ret;
00080     ret [0] = 1;
00081     for (int i = 1; i < (int)N; ++i)
00082         ret [i] = ret [i-1] * shape [i-1];
00083     return ret;
00084 }
00085 
00086 /********************************************************/
00087 /*                                                      */
00088 /*                 ScanOrderToOffset                    */
00089 /*                                                      */
00090 /********************************************************/
00091 
00092 /* transforms an index in scan order sense to a pointer offset in a possibly
00093    strided, multi-dimensional array.
00094 
00095     Namespace: vigra::detail
00096 */
00097 
00098 template <int K>
00099 struct ScanOrderToOffset
00100 {
00101     template <int N>
00102     static MultiArrayIndex
00103     exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> &shape,
00104          const TinyVector <MultiArrayIndex, N> & stride)
00105     {
00106         return stride[N-K] * (d % shape[N-K]) +
00107                ScanOrderToOffset<K-1>::exec(d / shape[N-K], shape, stride);
00108     }
00109 };
00110 
00111 template <>
00112 struct ScanOrderToOffset<1>
00113 {
00114     template <int N>
00115     static MultiArrayIndex
00116     exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> & /*shape*/,
00117          const TinyVector <MultiArrayIndex, N> & stride)
00118     {
00119         return stride[N-1] * d;
00120     }
00121 };
00122 
00123 template <int K>
00124 struct ScanOrderToCoordinate
00125 {
00126     template <int N>
00127     static void
00128     exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> &shape,
00129          TinyVector <MultiArrayIndex, N> & result)
00130     {
00131         result[N-K] = (d % shape[N-K]);
00132         ScanOrderToCoordinate<K-1>::exec(d / shape[N-K], shape, result);
00133     }
00134 };
00135 
00136 template <>
00137 struct ScanOrderToCoordinate<1>
00138 {
00139     template <int N>
00140     static void
00141     exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> & /*shape*/,
00142          TinyVector <MultiArrayIndex, N> & result)
00143     {
00144         result[N-1] = d;
00145     }
00146 };
00147 
00148 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 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
00219 
00220 Namespace: vigra::detail
00221 */
00222 template <unsigned int N>
00223 struct MaybeStrided
00224 {
00225     typedef UnstridedArrayTag type;
00226 };
00227 
00228 template <>
00229 struct MaybeStrided <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 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
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 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
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 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
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>
00410 inline void
00411 normMaxOfMultiArray(SrcIterator s, Shape const & shape, T & result, MetaInt<0>)
00412 {
00413     SrcIterator send = s + shape[0];
00414     for(; s < send; ++s)
00415     {
00416         T v = norm(*s);
00417         if(result < v)
00418             result = v;
00419     }
00420 }
00421 
00422 template <class SrcIterator, class Shape, class T, int N>
00423 void
00424 normMaxOfMultiArray(SrcIterator s, Shape const & shape, T & result, MetaInt<N>)
00425 {
00426     SrcIterator send = s + shape[N];
00427     for(; s < send; ++s)
00428     {
00429         normMaxOfMultiArray(s.begin(), shape, result, MetaInt<N-1>());
00430     }
00431 }
00432 
00433 template <class T>
00434 struct MultiArrayL1Functor
00435 {
00436     template <class U>
00437     T operator()(U t) const
00438     { return norm(t); }
00439 };
00440 
00441 template <class T>
00442 struct MultiArrayL2Functor
00443 {
00444     template <class U>
00445     T operator()(U t) const
00446     { return squaredNorm(t); }
00447 };
00448 
00449 template <class T>
00450 struct MultiArrayScaledL2Functor
00451 {
00452     T scale;
00453 
00454     MultiArrayScaledL2Functor(T s)
00455     : scale(s)
00456     {}
00457 
00458     template <class U>
00459     T operator()(U t) const
00460     { return squaredNorm(T(t) / scale); }
00461 };
00462 
00463 template <class SrcIterator, class Shape, class Functor, class T>
00464 inline void
00465 sumOverMultiArray(SrcIterator s, Shape const & shape, Functor f, T & result, MetaInt<0>)
00466 {
00467     SrcIterator send = s + shape[0];
00468     for(; s < send; ++s)
00469     {
00470         result += f(*s);
00471     }
00472 }
00473 
00474 template <class SrcIterator, class Shape, class Functor, class T, int N>
00475 void
00476 sumOverMultiArray(SrcIterator s, Shape const & shape, Functor f, T & result, MetaInt<N>)
00477 {
00478     SrcIterator send = s + shape[N];
00479     for(; s < send; ++s)
00480     {
00481         sumOverMultiArray(s.begin(), shape, f, result, MetaInt<N-1>());
00482     }
00483 }
00484 
00485 template <class SrcIterator, class Shape, class DestIterator>
00486 inline bool
00487 equalityOfMultiArrays(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>)
00488 {
00489     SrcIterator send = s + shape[0];
00490     for(; s < send; ++s, ++d)
00491     {
00492         if(!(*s == *d))
00493             return false;
00494     }
00495     return true;
00496 }
00497 
00498 template <class SrcIterator, class Shape, class DestIterator, int N>
00499 bool
00500 equalityOfMultiArrays(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>)
00501 {
00502     SrcIterator send = s + shape[N];
00503     for(; s < send; ++s, ++d)
00504     {
00505         if(!equalityOfMultiArrays(s.begin(), shape, d.begin(), MetaInt<N-1>()))
00506             return false;
00507     }
00508     return true;
00509 }
00510 
00511 
00512 template <class SrcIterator, class Shape, class DestIterator>
00513 inline void
00514 swapDataImpl(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>)
00515 {
00516     SrcIterator send = s + shape[0];
00517     for(; s < send; ++s, ++d)
00518         std::swap(*s, *d);
00519 }
00520 
00521 template <class SrcIterator, class Shape, class DestIterator, int N>
00522 void
00523 swapDataImpl(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>)
00524 {
00525     SrcIterator send = s + shape[N];
00526     for(; s < send; ++s, ++d)
00527         swapDataImpl(s.begin(), shape, d.begin(), MetaInt<N-1>());
00528 }
00529 
00530 } // namespace detail
00531 
00532 /********************************************************/
00533 /*                                                      */
00534 /*                     MultiArrayView                   */
00535 /*                                                      */
00536 /********************************************************/
00537 
00538 // forward declaration
00539 template <unsigned int N, class T, class C = UnstridedArrayTag>
00540 class MultiArrayView;
00541 template <unsigned int N, class T, class A = std::allocator<T> >
00542 class MultiArray;
00543 
00544 /********************************************************/
00545 /*                                                      */
00546 /*                       NormTraits                     */
00547 /*                                                      */
00548 /********************************************************/
00549 
00550 template <unsigned int N, class T, class C>
00551 struct NormTraits<MultiArrayView<N, T, C> >
00552 {
00553     typedef MultiArrayView<N, T, C>                                      Type;
00554     typedef typename NormTraits<T>::SquaredNormType SquaredNormType;
00555     typedef typename SquareRootTraits<SquaredNormType>::SquareRootResult NormType;
00556 };
00557 
00558 template <unsigned int N, class T, class A>
00559 struct NormTraits<MultiArray<N, T, A> >
00560 : public NormTraits<MultiArrayView<N, T, UnstridedArrayTag> >
00561 {
00562     typedef NormTraits<MultiArrayView<N, T, UnstridedArrayTag> > BaseType;
00563     typedef MultiArray<N, T, A>                                  Type;
00564     typedef typename BaseType::SquaredNormType                   SquaredNormType;
00565     typedef typename BaseType::NormType                          NormType;
00566 };
00567 
00568 /** \brief Base class for, and view to, \ref vigra::MultiArray.
00569 
00570 This class implements the interface of both MultiArray and
00571 MultiArrayView.  By default, MultiArrayViews are tagged as
00572 unstrided. If necessary, strided arrays are constructed automatically
00573 by calls to a variant of the bind...() function.
00574 
00575 If you want to apply an algorithm requiring an image to a
00576 <tt>MultiArrayView</tt> of appropriate (2-dimensional) shape, you can
00577 create a \ref vigra::BasicImageView that acts as a wrapper with the
00578 necessary interface -- see \ref MultiArrayToImage.
00579 
00580 The template parameter are as follows
00581 \code
00582     N: the array dimension
00583 
00584     T: the type of the array elements
00585 
00586     C: a tag determining whether the array's inner dimension is strided
00587        or not. An array is unstrided if the array elements occupy consecutive
00588        memory location, strided if there is an offset in between (e.g.
00589        when a view is created that skips every other array element).
00590        The compiler can generate faster code for unstrided arrays.
00591        Possible values: UnstridedArrayTag (default), StridedArrayTag
00592 \endcode
00593 
00594 <b>\#include</b>
00595 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
00596 
00597 Namespace: vigra
00598 */
00599 template <unsigned int N, class T, class C>
00600 class MultiArrayView
00601 {
00602 public:
00603 
00604         /** the array's actual dimensionality.
00605             This ensures that MultiArrayView can also be used for
00606             scalars (that is, when <tt>N == 0</tt>). Calculated as:<br>
00607             \code
00608             actual_dimension = (N==0) ? 1 : N
00609             \endcode
00610          */
00611     enum ActualDimension { actual_dimension = (N==0) ? 1 : N };
00612 
00613         /** the array's value type
00614          */
00615     typedef T value_type;
00616 
00617         /** reference type (result of operator[])
00618          */
00619     typedef value_type &reference;
00620 
00621         /** const reference type (result of operator[] const)
00622          */
00623     typedef const value_type &const_reference;
00624 
00625         /** pointer type
00626          */
00627     typedef value_type *pointer;
00628 
00629         /** const pointer type
00630          */
00631     typedef const value_type *const_pointer;
00632 
00633         /** difference type (used for multi-dimensional offsets and indices)
00634          */
00635     typedef typename MultiArrayShape<actual_dimension>::type difference_type;
00636 
00637         /** size type
00638          */
00639     typedef difference_type size_type;
00640 
00641         /** difference and index type for a single dimension
00642          */
00643     typedef MultiArrayIndex difference_type_1;
00644 
00645         /** traverser (MultiIterator) type
00646          */
00647     typedef typename vigra::detail::MultiIteratorChooser <
00648         C>::template Traverser <actual_dimension, T, T &, T *>::type traverser;
00649 
00650         /** const traverser (MultiIterator) type
00651          */
00652     typedef typename vigra::detail::MultiIteratorChooser <
00653         C>::template Traverser <actual_dimension, T, T const &, T const *>::type const_traverser;
00654 
00655         /** the view type associated with this array.
00656          */
00657     typedef MultiArrayView <N, T, C> view_type;
00658 
00659         /** the matrix type associated with this array.
00660          */
00661     typedef MultiArray <N, T> matrix_type;
00662 
00663 protected:
00664 
00665     typedef typename difference_type::value_type diff_zero_t;
00666 
00667         /** the shape of the image pointed to is stored here.
00668         */
00669     difference_type m_shape;
00670 
00671         /** the strides (offset of a sample to the next) for every dimension
00672             are stored here.
00673         */
00674     difference_type m_stride;
00675 
00676         /** pointer to the image.
00677          */
00678     pointer m_ptr;
00679 
00680     template <class U, class CN>
00681     void copyImpl(const MultiArrayView <N, U, CN>& rhs);
00682 
00683     template <class U, class CN>
00684     void swapDataImpl(MultiArrayView <N, U, CN> rhs);
00685 
00686     template <class CN>
00687     bool arraysOverlap(const MultiArrayView <N, T, CN>& rhs) const;
00688 
00689     template <class U, class CN>
00690     bool arraysOverlap(const MultiArrayView <N, U, CN>&) const
00691     {
00692         return false;
00693     }
00694 
00695 public:
00696 
00697         /** default constructor: create an invalid view,
00698          * i.e. hasData() returns false and size() is zero.
00699          */
00700     MultiArrayView ()
00701         : m_shape (diff_zero_t(0)), m_stride (diff_zero_t(0)), m_ptr (0)
00702     {}
00703 
00704         /** construct from shape and pointer
00705          */
00706     MultiArrayView (const difference_type &shape, pointer ptr)
00707     : m_shape (shape),
00708       m_stride (detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape)),
00709       m_ptr (ptr)
00710     {}
00711 
00712         /** Construct from shape, strides (offset of a sample to the
00713             next) for every dimension, and pointer.  (Note that
00714             strides are not given in bytes, but in offsets of the
00715             respective pointer type.)
00716          */
00717     MultiArrayView (const difference_type &shape,
00718                     const difference_type &stride,
00719                     pointer ptr)
00720     : m_shape (shape),
00721       m_stride (stride),
00722       m_ptr (ptr)
00723     {}
00724 
00725 
00726         /** Assignment. There are 3 cases:
00727 
00728             <ul>
00729             <li> When this <tt>MultiArrayView</tt> does not point to valid data
00730                  (e.g. after default construction), it becomes a copy of \a rhs.
00731             <li> When the shapes of the two arrays match, the array contents are copied.
00732             <li> Otherwise, a <tt>PreconditionViolation</tt> exception is thrown.
00733             </ul>
00734          */
00735     MultiArrayView & operator=(MultiArrayView const & rhs);
00736 
00737         /** Assignment of a differently typed MultiArrayView. Fails with
00738             <tt>PreconditionViolation</tt> exception when the shapes do not match.
00739          */
00740     template<class U, class C1>
00741     MultiArrayView & operator=(MultiArrayView<N, U, C1> const & rhs)
00742     {
00743         vigra_precondition(this->shape() == rhs.shape(),
00744             "MultiArrayView::operator=() size mismatch.");
00745         this->copyImpl(rhs);
00746         return *this;
00747     }
00748 
00749         /** Add-assignment of a compatible MultiArrayView. Fails with
00750             <tt>PreconditionViolation</tt> exception when the shapes do not match.
00751          */
00752     template<class U, class C1>
00753     MultiArrayView & operator+=(MultiArrayView<N, U, C1> const & rhs);
00754 
00755         /** Subtract-assignment of a compatible MultiArrayView. Fails with
00756             <tt>PreconditionViolation</tt> exception when the shapes do not match.
00757          */
00758     template<class U, class C1>
00759     MultiArrayView & operator-=(MultiArrayView<N, U, C1> const & rhs);
00760 
00761         /** Multiply-assignment of a compatible MultiArrayView. Fails with
00762             <tt>PreconditionViolation</tt> exception when the shapes do not match.
00763          */
00764     template<class U, class C1>
00765     MultiArrayView & operator*=(MultiArrayView<N, U, C1> const & rhs);
00766 
00767         /** Divide-assignment of a compatible MultiArrayView. Fails with
00768             <tt>PreconditionViolation</tt> exception when the shapes do not match.
00769          */
00770     template<class U, class C1>
00771     MultiArrayView & operator/=(MultiArrayView<N, U, C1> const & rhs);
00772 
00773         /** Add-assignment of a scalar.
00774          */
00775     MultiArrayView & operator+=(T const & rhs)
00776     {
00777         detail::copyAddScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
00778         return *this;
00779     }
00780 
00781         /** Subtract-assignment of a scalar.
00782          */
00783     MultiArrayView & operator-=(T const & rhs)
00784     {
00785         detail::copySubScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
00786         return *this;
00787     }
00788 
00789         /** Multiply-assignment of a scalar.
00790          */
00791     MultiArrayView & operator*=(T const & rhs)
00792     {
00793         detail::copyMulScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
00794         return *this;
00795     }
00796 
00797         /** Divide-assignment of a scalar.
00798          */
00799     MultiArrayView & operator/=(T const & rhs)
00800     {
00801         detail::copyDivScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
00802         return *this;
00803     }
00804 
00805         /** array access.
00806          */
00807     reference operator[] (const difference_type &d)
00808     {
00809         VIGRA_ASSERT_INSIDE(d);
00810         return m_ptr [dot (d, m_stride)];
00811     }
00812 
00813         /** array access.
00814          */
00815     const_reference operator[] (const difference_type &d) const
00816     {
00817         VIGRA_ASSERT_INSIDE(d);
00818         return m_ptr [dot (d, m_stride)];
00819     }
00820 
00821         /** equvalent to bindInner(), when M < N.
00822          */
00823     template <unsigned int M>
00824     MultiArrayView <N-M, T, StridedArrayTag> operator[] (const TinyVector<MultiArrayIndex, M> &d) const
00825     {
00826         return bindInner(d);
00827     }
00828 
00829         /** Array access in scan-order sense.
00830             Mostly useful to support standard indexing for 1-dimensional multi-arrays,
00831             but works for any N. Use scanOrderIndexToCoordinate() and
00832             coordinateToScanOrderIndex() for conversion between indices and coordinates.
00833             
00834             <b>Note:</b> This function should not be used in the inner loop, because the
00835             conversion of the scan order index into a memory address is expensive
00836             (it must take into account that memory may not be consecutive for subarrays
00837             and/or strided arrays). Always prefer operator() if possible.            
00838          */
00839     reference operator[](difference_type_1 d)
00840     {
00841         VIGRA_ASSERT_INSIDE(scanOrderIndexToCoordinate(d));
00842         return m_ptr [detail::ScanOrderToOffset<actual_dimension>::exec(d, m_shape, m_stride)];
00843     }
00844 
00845         /** Array access in scan-order sense.
00846             Mostly useful to support standard indexing for 1-dimensional multi-arrays,
00847             but works for any N. Use scanOrderIndexToCoordinate() and
00848             coordinateToScanOrderIndex() for conversion between indices and coordinates.
00849              
00850             <b>Note:</b> This function should not be used in the inner loop, because the
00851             conversion of the scan order index into a memory address is expensive
00852             (it must take into account that memory may not be consecutive for subarrays
00853             and/or strided arrays). Always prefer operator() if possible.            
00854         */
00855     const_reference operator[](difference_type_1 d) const
00856     {
00857         VIGRA_ASSERT_INSIDE(scanOrderIndexToCoordinate(d));
00858         return m_ptr [detail::ScanOrderToOffset<actual_dimension>::exec(d, m_shape, m_stride)];
00859     }
00860 
00861         /** convert scan-order index to coordinate.
00862          */
00863     difference_type scanOrderIndexToCoordinate(difference_type_1 d) const
00864     {
00865         difference_type result;
00866         detail::ScanOrderToCoordinate<actual_dimension>::exec(d, m_shape, result);
00867         return result;
00868     }
00869 
00870         /** convert coordinate to scan-order index.
00871          */
00872     difference_type_1 coordinateToScanOrderIndex(const difference_type &d) const
00873     {
00874         return detail::CoordinateToScanOrder<actual_dimension>::exec(m_shape, d);
00875     }
00876 
00877         /** 1D array access. Use only if N == 1.
00878          */
00879     reference operator() (difference_type_1 x)
00880     {
00881         VIGRA_ASSERT_INSIDE(difference_type(x));
00882         return m_ptr [detail::CoordinatesToOffest<C>::exec(m_stride, x)];
00883     }
00884 
00885         /** 2D array access. Use only if N == 2.
00886          */
00887     reference operator() (difference_type_1 x, difference_type_1 y)
00888     {
00889         VIGRA_ASSERT_INSIDE(difference_type(x, y));
00890         return m_ptr [detail::CoordinatesToOffest<C>::exec(m_stride, x, y)];
00891     }
00892 
00893         /** 3D array access. Use only if N == 3.
00894          */
00895     reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z)
00896     {
00897         VIGRA_ASSERT_INSIDE(difference_type(x, y, z));
00898         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
00899     }
00900 
00901         /** 4D array access. Use only if N == 4.
00902          */
00903     reference operator() (difference_type_1 x, difference_type_1 y,
00904                           difference_type_1 z, difference_type_1 u)
00905     {
00906         VIGRA_ASSERT_INSIDE(difference_type(x, y, z, u));
00907         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
00908     }
00909 
00910         /** 5D array access. Use only if N == 5.
00911          */
00912     reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
00913                           difference_type_1 u, difference_type_1 v)
00914     {
00915         VIGRA_ASSERT_INSIDE(difference_type(x, y,z, u,v));
00916         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
00917     }
00918 
00919         /** 1D const array access. Use only if N == 1.
00920          */
00921     const_reference operator() (difference_type_1 x) const
00922     {
00923         VIGRA_ASSERT_INSIDE(difference_type(x));
00924         return m_ptr [detail::CoordinatesToOffest<C>::exec(m_stride, x)];
00925     }
00926 
00927         /** 2D const array access. Use only if N == 2.
00928          */
00929     const_reference operator() (difference_type_1 x, difference_type_1 y) const
00930     {
00931         VIGRA_ASSERT_INSIDE(difference_type(x, y));
00932         return m_ptr [detail::CoordinatesToOffest<C>::exec(m_stride, x, y)];
00933     }
00934 
00935         /** 3D const array access. Use only if N == 3.
00936          */
00937     const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z) const
00938     {
00939         VIGRA_ASSERT_INSIDE(difference_type(x,y,z));
00940         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
00941     }
00942 
00943         /** 4D const array access. Use only if N == 4.
00944          */
00945     const_reference operator() (difference_type_1 x, difference_type_1 y,
00946                                 difference_type_1 z, difference_type_1 u) const
00947     {
00948         VIGRA_ASSERT_INSIDE(difference_type(x,y,z,u));
00949         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
00950     }
00951 
00952         /** 5D const array access. Use only if N == 5.
00953          */
00954     const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
00955                                 difference_type_1 u, difference_type_1 v) const
00956     {
00957         VIGRA_ASSERT_INSIDE(difference_type(x,y,z,u,v));
00958         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
00959     }
00960 
00961         /** Init with a constant.
00962          */
00963     template <class U>
00964     MultiArrayView & init(const U & init)
00965     {
00966         detail::copyScalarMultiArrayData(traverser_begin(), shape(), init, MetaInt<actual_dimension-1>());
00967         return *this;
00968     }
00969 
00970 
00971         /** Copy the data of the right-hand array (array shapes must match).
00972          */
00973     void copy(const MultiArrayView & rhs)
00974     {
00975         if(this == &rhs)
00976             return;
00977         this->copyImpl(rhs);
00978     }
00979 
00980         /** Copy the data of the right-hand array (array shapes must match).
00981          */
00982     template <class U, class CN>
00983     void copy(const MultiArrayView <N, U, CN>& rhs)
00984     {
00985         this->copyImpl(rhs);
00986     }
00987 
00988         /** swap the data between two MultiArrayView objects.
00989 
00990             The shapes of the two array must match.
00991         */
00992     void swapData(MultiArrayView rhs)
00993     {
00994         if(this != &rhs)
00995             swapDataImpl(rhs);
00996     }
00997 
00998         /** swap the data between two MultiArrayView objects.
00999 
01000             The shapes of the two array must match.
01001         */
01002     template <class T2, class C2>
01003     void swapData(MultiArrayView <N, T2, C2> rhs)
01004     {
01005         swapDataImpl(rhs);
01006     }
01007 
01008         /** bind the M outmost dimensions to certain indices.
01009             this reduces the dimensionality of the image to
01010             max { 1, N-M }.
01011 
01012             <b>Usage:</b>
01013             \code
01014             // create a 3D array of size 40x30x20
01015             typedef MultiArray<3, double>::difference_type Shape;
01016             MultiArray<3, double> array3(Shape(40, 30, 20));
01017 
01018             // get a 1D array by fixing index 1 to 12, and index 2 to 10
01019             MultiArrayView <1, double> array1 = array3.bindOuter(TinyVector<MultiArrayIndex, 2>(12, 10));
01020             \endcode
01021         */
01022     template <unsigned int M>
01023     MultiArrayView <N-M, T, C> bindOuter (const TinyVector <MultiArrayIndex, M> &d) const;
01024 
01025         /** bind the M innermost dimensions to certain indices.
01026             this reduces the dimensionality of the image to
01027             max { 1, N-M }.
01028 
01029             <b>Usage:</b>
01030             \code
01031             // create a 3D array of size 40x30x20
01032             typedef MultiArray<3, double>::difference_type Shape;
01033             MultiArray<3, double> array3(Shape(40, 30, 20));
01034 
01035             // get a 1D array by fixing index 0 to 12, and index 1 to 10
01036             MultiArrayView <1, double, StridedArrayTag> array1 = array3.bindInner(TinyVector<MultiArrayIndex, 2>(12, 10));
01037             \endcode
01038         */
01039     template <unsigned int M>
01040     MultiArrayView <N-M, T, StridedArrayTag>
01041     bindInner (const TinyVector <MultiArrayIndex, M> &d) const;
01042 
01043         /** bind dimension M to index d.
01044             this reduces the dimensionality of the image to
01045             max { 1, N-1 }.
01046 
01047             <b>Usage:</b>
01048             \code
01049             // create a 3D array of size 40x30x20
01050             typedef MultiArray<3, double>::difference_type Shape;
01051             MultiArray<3, double> array3(Shape(40, 30, 20));
01052 
01053             // get a 2D array by fixing index 1 to 12
01054             MultiArrayView <2, double> array2 = array3.bind<1>(12);
01055 
01056             // get a 2D array by fixing index 0 to 23
01057             MultiArrayView <2, double, StridedArrayTag> array2a = array3.bind<0>(23);
01058             \endcode
01059          */
01060     template <unsigned int M>
01061     MultiArrayView <N-1, T, typename vigra::detail::MaybeStrided <M>::type >
01062     bind (difference_type_1 d) const;
01063 
01064         /** bind the outmost dimension to a certain index.
01065             this reduces the dimensionality of the image to
01066             max { 1, N-1 }.
01067 
01068             <b>Usage:</b>
01069             \code
01070             // create a 3D array of size 40x30x20
01071             typedef MultiArray<3, double>::difference_type Shape;
01072             MultiArray<3, double> array3(Shape(40, 30, 20));
01073 
01074             // get a 2D array by fixing the outermost index (i.e. index 2) to 12
01075             MultiArrayView <2, double> array2 = array3.bindOuter(12);
01076             \endcode
01077         */
01078     MultiArrayView <N-1, T, C> bindOuter (difference_type_1 d) const;
01079 
01080         /** bind the innermost dimension to a certain index.
01081             this reduces the dimensionality of the image to
01082             max { 1, N-1 }.
01083 
01084             <b>Usage:</b>
01085             \code
01086             // create a 3D array of size 40x30x20
01087             typedef MultiArray<3, double>::difference_type Shape;
01088             MultiArray<3, double> array3(Shape(40, 30, 20));
01089 
01090             // get a 2D array by fixing the innermost index (i.e. index 0) to 23
01091             MultiArrayView <2, double, StridedArrayTag> array2 = array3.bindInner(23);
01092             \endcode
01093         */
01094     MultiArrayView <N-1, T, StridedArrayTag> bindInner (difference_type_1 d) const;
01095 
01096         /** bind dimension m to index d.
01097             this reduces the dimensionality of the image to
01098             max { 1, N-1 }.
01099 
01100             <b>Usage:</b>
01101             \code
01102             // create a 3D array of size 40x30x20
01103             typedef MultiArray<3, double>::difference_type Shape;
01104             MultiArray<3, double> array3(Shape(40, 30, 20));
01105 
01106             // get a 2D array by fixing index 2 to 15
01107             MultiArrayView <2, double, StridedArrayTag> array2 = array3.bindAt(2, 15);
01108             \endcode
01109          */
01110     MultiArrayView <N-1, T, StridedArrayTag>
01111     bindAt (difference_type_1 m, difference_type_1 d) const;
01112 
01113         /** Add a singleton dimension (dimension of legth 1).
01114 
01115             Singleton dimensions don't change the size of the data, but introduce
01116             a new index that can only take the value 0. This is mainly useful for
01117             the 'reduce mode' of transformMultiArray() and combineTwoMultiArrays(),
01118             because these functions require the source and destination arrays to
01119             have the same number of dimensions.
01120 
01121             The range of \a i must be <tt>0 <= i <= N</tt>. The new dimension will become
01122             the i'th index, and the old indices from i upwards will shift one
01123             place to the right.
01124 
01125             <b>Usage:</b>
01126 
01127             Suppose we want have a 2D array and want to create a 1D array that contains
01128             the row average of the first array.
01129             \code
01130             typedef MultiArrayShape<2>::type Shape2;
01131             MultiArray<2, double> original(Shape2(40, 30));
01132 
01133             typedef MultiArrayShape<1>::type Shape1;
01134             MultiArray<1, double> rowAverages(Shape1(30));
01135 
01136             // temporarily add a singleton dimension to the destination array
01137             transformMultiArray(srcMultiArrayRange(original),
01138                                 destMultiArrayRange(rowAverages.insertSingletonDimension(0)),
01139                                 FindAverage<double>());
01140             \endcode
01141          */
01142     MultiArrayView <N+1, T, C>
01143     insertSingletonDimension (difference_type_1 i) const;
01144 
01145         /** create a rectangular subarray that spans between the
01146             points p and q, where p is in the subarray, q not.
01147 
01148             <b>Usage:</b>
01149             \code
01150             // create a 3D array of size 40x30x20
01151             typedef MultiArray<3, double>::difference_type Shape;
01152             MultiArray<3, double> array3(Shape(40, 30, 20));
01153 
01154             // get a subarray set is smaller by one element at all sides
01155             MultiArrayView <3, double> subarray = array3.subarray(Shape(1,1,1), Shape(39, 29, 19));
01156             \endcode
01157         */
01158     MultiArrayView subarray (const difference_type &p,
01159                              const difference_type &q) const
01160     {
01161         const difference_type_1 offset = dot (m_stride, p);
01162         return MultiArrayView (q - p, m_stride, m_ptr + offset);
01163     }
01164 
01165         /** apply an additional striding to the image, thereby reducing
01166             the shape of the array.
01167             for example, multiplying the stride of dimension one by three
01168             turns an appropriately layed out (interleaved) rgb image into
01169             a single band image.
01170         */
01171     MultiArrayView <N, T, StridedArrayTag>
01172     stridearray (const difference_type &s) const
01173     {
01174         difference_type shape = m_shape;
01175         for (unsigned int i = 0; i < actual_dimension; ++i)
01176             shape [i] /= s [i];
01177         return MultiArrayView <N, T, StridedArrayTag>(shape, m_stride * s, m_ptr);
01178     }
01179 
01180         /** Transpose an array. If N==2, this implements the usual matrix transposition.
01181             For N > 2, it reverses the order of the indices.
01182 
01183             <b>Usage:</b><br>
01184             \code
01185             typedef MultiArray<2, double>::difference_type Shape;
01186             MultiArray<2, double> array(10, 20);
01187 
01188             MultiArray<2, double, StridedArrayTag> transposed = array.transpose();
01189 
01190             for(int i=0; i<array.shape(0), ++i)
01191                 for(int j=0; j<array.shape(1); ++j)
01192                     assert(array(i, j) == transposed(j, i));
01193             \endcode
01194         */
01195     MultiArrayView <N, T, StridedArrayTag>
01196     transpose () const
01197     {
01198         difference_type shape(m_shape.begin(), difference_type::ReverseCopy),
01199                         stride(m_stride.begin(), difference_type::ReverseCopy);
01200         return MultiArrayView <N, T, StridedArrayTag>(shape, stride, m_ptr);
01201     }
01202 
01203         /** permute the dimensions of the array.
01204             The function exchanges the meaning of the dimensions without copying the data.
01205             In case of a 2-dimensional array, this is simply array transposition. In higher dimensions,
01206             there are more possibilities.
01207 
01208             <b>Usage:</b><br>
01209             \code
01210             typedef MultiArray<2, double>::difference_type Shape;
01211             MultiArray<2, double> array(10, 20);
01212 
01213             MultiArray<2, double, StridedArrayTag> transposed = array.permuteDimensions(Shape(1,0));
01214 
01215             for(int i=0; i<array.shape(0), ++i)
01216                 for(int j=0; j<array.shape(1); ++j)
01217                     assert(array(i, j) == transposed(j, i));
01218             \endcode
01219         */
01220     MultiArrayView <N, T, StridedArrayTag>
01221     permuteDimensions (const difference_type &s) const;
01222 
01223         /** Permute the dimensions of the array so that the strides are in ascending order.
01224             Determines the appropriate permutation and then calls permuteDimensions().
01225         */
01226     MultiArrayView <N, T, StridedArrayTag>
01227     permuteStridesAscending() const;
01228     
01229         /** Permute the dimensions of the array so that the strides are in descending order.
01230             Determines the appropriate permutation and then calls permuteDimensions().
01231         */
01232     MultiArrayView <N, T, StridedArrayTag>
01233     permuteStridesDescending() const;
01234     
01235         /** Compute the ordering of the strides in this array.
01236             The result is describes the current permutation of the axes relative 
01237             to the standard ascending stride order.
01238         */
01239     difference_type strideOrdering() const
01240     {
01241         return strideOrdering(m_stride);
01242     }
01243     
01244         /** Compute the ordering of the given strides.
01245             The result is describes the current permutation of the axes relative 
01246             to the standard ascending stride order.
01247         */
01248     static difference_type strideOrdering(difference_type strides);
01249 
01250         /** number of the elements in the array.
01251          */
01252     difference_type_1 elementCount () const
01253     {
01254         difference_type_1 ret = m_shape[0];
01255         for(int i = 1; i < actual_dimension; ++i)
01256             ret *= m_shape[i];
01257         return ret;
01258     }
01259 
01260         /** number of the elements in the array.
01261             Same as <tt>elementCount()</tt>. Mostly useful to support the std::vector interface.
01262          */
01263     difference_type_1 size () const
01264     {
01265         return elementCount();
01266     }
01267 
01268         /** return the array's shape.
01269          */
01270     const difference_type & shape () const
01271     {
01272         return m_shape;
01273     }
01274 
01275         /** return the array's size at a certain dimension.
01276          */
01277     difference_type_1 size (difference_type_1 n) const
01278     {
01279         return m_shape [n];
01280     }
01281 
01282         /** return the array's shape at a certain dimension
01283             (same as <tt>size(n)</tt>).
01284          */
01285     difference_type_1 shape (difference_type_1 n) const
01286     {
01287         return m_shape [n];
01288     }
01289 
01290         /** return the array's stride for every dimension.
01291          */
01292     const difference_type & stride () const
01293     {
01294         return m_stride;
01295     }
01296 
01297         /** return the array's stride at a certain dimension.
01298          */
01299     difference_type_1 stride (int n) const
01300     {
01301         return m_stride [n];
01302     }
01303 
01304         /** check whether two arrays are elementwise equal.
01305          */
01306     template <class U, class C1>
01307     bool operator==(MultiArrayView<N, U, C1> const & rhs) const
01308     {
01309         if(this->shape() != rhs.shape())
01310             return false;
01311         return detail::equalityOfMultiArrays(traverser_begin(), shape(), rhs.traverser_begin(), MetaInt<actual_dimension-1>());
01312     }
01313 
01314         /** check whether two arrays are not elementwise equal.
01315             Also true when the two arrays have different shapes.
01316          */
01317     template <class U, class C1>
01318     bool operator!=(MultiArrayView<N, U, C1> const & rhs) const
01319     {
01320         return !operator==(rhs);
01321     }
01322 
01323         /** check whether the given point is in the array range.
01324          */
01325     bool isInside (difference_type const & p) const
01326     {
01327         for(int d=0; d<actual_dimension; ++d)
01328             if(p[d] < 0 || p[d] >= shape(d))
01329                 return false;
01330         return true;
01331     }
01332 
01333         /** Compute the squared Euclidean norm of the array (sum of squares of the array elements).
01334          */
01335     typename NormTraits<MultiArrayView>::SquaredNormType squaredNorm() const
01336     {
01337         typedef typename NormTraits<MultiArrayView>::SquaredNormType SquaredNormType;
01338         SquaredNormType res = NumericTraits<SquaredNormType>::zero();
01339         detail::sumOverMultiArray(traverser_begin(), shape(), detail::MultiArrayL2Functor<SquaredNormType>(),
01340                                   res, MetaInt<actual_dimension-1>());
01341         return res;
01342     }
01343 
01344         /** Compute various norms of the array.
01345             The norm is determined by parameter \a type:
01346 
01347             <ul>
01348             <li> type == 0: maximum norm (L-infinity): maximum of absolute values of the array elements
01349             <li> type == 1: Manhattan norm (L1): sum of absolute values of the array elements
01350             <li> type == 2: Euclidean norm (L2): square root of <tt>squaredNorm()</tt> when \a useSquaredNorm is <tt>true</tt>,<br>
01351                  or direct algorithm that avoids underflow/overflow otherwise.
01352             </ul>
01353 
01354             Parameter \a useSquaredNorm has no effect when \a type != 2. Defaults: compute L2 norm as square root of
01355             <tt>squaredNorm()</tt>.
01356          */
01357     typename NormTraits<MultiArrayView>::NormType norm(int type = 2, bool useSquaredNorm = true) const;
01358 
01359         /** return the pointer to the image data
01360          */
01361     pointer data () const
01362     {
01363         return m_ptr;
01364     }
01365 
01366         /**
01367          * returns true iff this view refers to valid data,
01368          * i.e. data() is not a NULL pointer.  (this is false after
01369          * default construction.)
01370          */
01371     bool hasData () const
01372     {
01373         return m_ptr != 0;
01374     }
01375 
01376         /** returns the N-dimensional MultiIterator pointing
01377             to the first element in every dimension.
01378         */
01379     traverser traverser_begin ()
01380     {
01381         traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
01382         return ret;
01383     }
01384 
01385         /** returns the N-dimensional MultiIterator pointing
01386             to the const first element in every dimension.
01387         */
01388     const_traverser traverser_begin () const
01389     {
01390         const_traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
01391         return ret;
01392     }
01393 
01394         /** returns the N-dimensional MultiIterator pointing
01395             beyond the last element in dimension N, and to the
01396             first element in every other dimension.
01397         */
01398     traverser traverser_end ()
01399     {
01400         traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
01401         ret += m_shape [actual_dimension-1];
01402         return ret;
01403     }
01404 
01405         /** returns the N-dimensional const MultiIterator pointing
01406             beyond the last element in dimension N, and to the
01407             first element in every other dimension.
01408         */
01409     const_traverser traverser_end () const
01410     {
01411         const_traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
01412         ret += m_shape [actual_dimension-1];
01413         return ret;
01414     }
01415 
01416     view_type view ()
01417     {
01418         return *this;
01419     }
01420 };
01421 
01422 template <unsigned int N, class T, class C>
01423 MultiArrayView<N, T, C> &
01424 MultiArrayView <N, T, C>::operator=(MultiArrayView<N, T, C> const & rhs)
01425 {
01426     if(this == &rhs)
01427         return *this;
01428     vigra_precondition(this->shape() == rhs.shape() || m_ptr == 0,
01429         "MultiArrayView::operator=(MultiArrayView const &) size mismatch - use MultiArrayView::reset().");
01430     if(m_ptr == 0)
01431     {
01432         m_shape  = rhs.m_shape;
01433         m_stride = rhs.m_stride;
01434         m_ptr    = rhs.m_ptr;
01435     }
01436     else
01437         this->copyImpl(rhs);
01438     return *this;
01439 }
01440 
01441 template <unsigned int N, class T, class C>
01442 template <class CN>
01443 bool
01444 MultiArrayView <N, T, C>::arraysOverlap(const MultiArrayView <N, T, CN>& rhs) const
01445 {
01446     vigra_precondition (shape () == rhs.shape (),
01447         "MultiArrayView::arraysOverlap(): shape mismatch.");
01448     const_pointer first_element = this->m_ptr,
01449                   last_element = first_element + dot(this->m_shape - difference_type(1), this->m_stride);
01450     typename MultiArrayView <N, T, CN>::const_pointer
01451            rhs_first_element = rhs.data(),
01452            rhs_last_element = rhs_first_element + dot(rhs.shape() - difference_type(1), rhs.stride());
01453     return !(last_element < rhs_first_element || rhs_last_element < first_element);
01454 }
01455 
01456 template <unsigned int N, class T, class C>
01457 template <class U, class CN>
01458 void
01459 MultiArrayView <N, T, C>::copyImpl(const MultiArrayView <N, U, CN>& rhs)
01460 {
01461     if(!arraysOverlap(rhs))
01462     {
01463         // no overlap -- can copy directly
01464         detail::copyMultiArrayData(rhs.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>());
01465     }
01466     else
01467     {
01468         // overlap: we got different views to the same data -- copy to intermediate memory in order to avoid
01469         // overwriting elements that are still needed on the rhs.
01470         MultiArray<N, T> tmp(rhs);
01471         detail::copyMultiArrayData(tmp.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>());
01472     }
01473 }
01474 
01475 #define VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT(name, op) \
01476 template <unsigned int N, class T, class C> \
01477 template<class U, class C1> \
01478 MultiArrayView<N, T, C> &  \
01479 MultiArrayView <N, T, C>::operator op(MultiArrayView<N, U, C1> const & rhs) \
01480 { \
01481     vigra_precondition(this->shape() == rhs.shape(), "MultiArrayView::operator" #op "() size mismatch."); \
01482     if(!arraysOverlap(rhs)) \
01483     { \
01484         detail::name##MultiArrayData(rhs.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>()); \
01485     } \
01486     else \
01487     { \
01488         MultiArray<N, T> tmp(rhs); \
01489         detail::name##MultiArrayData(tmp.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>()); \
01490     } \
01491     return *this; \
01492 }
01493 
01494 VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT(copyAdd, +=)
01495 VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT(copySub, -=)
01496 VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT(copyMul, *=)
01497 VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT(copyDiv, /=)
01498 
01499 #undef VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT
01500 
01501 template <unsigned int N, class T, class C>
01502 template <class U, class CN>
01503 void
01504 MultiArrayView <N, T, C>::swapDataImpl(MultiArrayView <N, U, CN> rhs)
01505 {
01506     vigra_precondition (shape () == rhs.shape (),
01507         "MultiArrayView::swapData(): shape mismatch.");
01508 
01509     // check for overlap of this and rhs
01510     const_pointer first_element = this->m_ptr,
01511                   last_element = first_element + dot(this->m_shape - difference_type(1), this->m_stride);
01512     typename MultiArrayView <N, U, CN>::const_pointer
01513            rhs_first_element = rhs.data(),
01514            rhs_last_element = rhs_first_element + dot(rhs.shape() - difference_type(1), rhs.stride());
01515     if(last_element < rhs_first_element || rhs_last_element < first_element)
01516     {
01517         // no overlap -- can swap directly
01518         detail::swapDataImpl(traverser_begin(), shape(), rhs.traverser_begin(), MetaInt<actual_dimension-1>());
01519     }
01520     else
01521     {
01522         // overlap: we got different views to the same data -- copy to intermediate memory in order to avoid
01523         // overwriting elements that are still needed.
01524         MultiArray<N, T> tmp(*this);
01525         copy(rhs);
01526         rhs.copy(tmp);
01527     }
01528 }
01529 
01530 template <unsigned int N, class T, class C>
01531 MultiArrayView <N, T, StridedArrayTag>
01532 MultiArrayView <N, T, C>::permuteDimensions (const difference_type &s) const
01533 {
01534     difference_type shape, stride, check((typename difference_type::value_type)0);
01535     for (unsigned int i = 0; i < actual_dimension; ++i)
01536     {
01537         shape[i]  = m_shape[s[i]];
01538         stride[i] = m_stride[s[i]];
01539         ++check[s[i]];
01540     }
01541     vigra_precondition(check == difference_type(1),
01542        "MultiArrayView::permuteDimensions(): every dimension must occur exactly once.");
01543     return MultiArrayView <N, T, StridedArrayTag>(shape, stride, m_ptr);
01544 }
01545 
01546 template <unsigned int N, class T, class C>
01547 typename MultiArrayView <N, T, C>::difference_type 
01548 MultiArrayView <N, T, C>::strideOrdering(difference_type stride)
01549 {
01550     difference_type permutation;
01551     for(unsigned int k=0; k<N; ++k)
01552         permutation[k] = k;
01553     for(unsigned int k=0; k<N-1; ++k)
01554     {
01555         unsigned int smallest = k;
01556         for(unsigned int j=k+1; j<N; ++j)
01557         {
01558             if(stride[j] < stride[smallest])
01559                 smallest = j;
01560         }
01561         if(smallest != k)
01562         {
01563             std::swap(stride[k], stride[smallest]);
01564             std::swap(permutation[k], permutation[smallest]);
01565         }
01566     }
01567     difference_type ordering;
01568     for(unsigned int k=0; k<N; ++k)
01569         ordering[permutation[k]] = k;
01570     return ordering;
01571 }
01572 
01573 template <unsigned int N, class T, class C>
01574 MultiArrayView <N, T, StridedArrayTag>
01575 MultiArrayView <N, T, C>::permuteStridesAscending() const
01576 {
01577     difference_type ordering(strideOrdering(m_stride)), permutation;
01578     for(MultiArrayIndex k=0; k<N; ++k)
01579         permutation[ordering[k]] = k;
01580     return permuteDimensions(permutation);
01581 }
01582 
01583 template <unsigned int N, class T, class C>
01584 MultiArrayView <N, T, StridedArrayTag>
01585 MultiArrayView <N, T, C>::permuteStridesDescending() const
01586 {
01587     difference_type ordering(strideOrdering(m_stride)), permutation;
01588     for(MultiArrayIndex k=0; k<N; ++k)
01589         permutation[ordering[N-1-k]] = k;
01590     return permuteDimensions(permutation);
01591 }
01592 
01593 template <unsigned int N, class T, class C>
01594 template <unsigned int M>
01595 MultiArrayView <N-M, T, C>
01596 MultiArrayView <N, T, C>::bindOuter (const TinyVector <MultiArrayIndex, M> &d) const
01597 {
01598     TinyVector <MultiArrayIndex, M> stride;
01599     stride.init (m_stride.begin () + N-M, m_stride.end ());
01600     pointer ptr = m_ptr + dot (d, stride);
01601     static const int NNew = (N-M == 0) ? 1 : N-M;
01602     TinyVector <MultiArrayIndex, NNew> inner_shape, inner_stride;
01603     if (N-M == 0)
01604     {
01605         inner_shape [0] = 1;
01606         inner_stride [0] = 0;
01607     }
01608     else
01609     {
01610         inner_shape.init (m_shape.begin (), m_shape.end () - M);
01611         inner_stride.init (m_stride.begin (), m_stride.end () - M);
01612     }
01613     return MultiArrayView <N-M, T, C> (inner_shape, inner_stride, ptr);
01614 }
01615 
01616 template <unsigned int N, class T, class C>
01617 template <unsigned int M>
01618 MultiArrayView <N - M, T, StridedArrayTag>
01619 MultiArrayView <N, T, C>::bindInner (const TinyVector <MultiArrayIndex, M> &d) const
01620 {
01621     TinyVector <MultiArrayIndex, M> stride;
01622     stride.init (m_stride.begin (), m_stride.end () - N + M);
01623     pointer ptr = m_ptr + dot (d, stride);
01624     static const int NNew = (N-M == 0) ? 1 : N-M;
01625     TinyVector <MultiArrayIndex, NNew> outer_shape, outer_stride;
01626     if (N-M == 0)
01627     {
01628         outer_shape [0] = 1;
01629         outer_stride [0] = 0;
01630     }
01631     else
01632     {
01633         outer_shape.init (m_shape.begin () + M, m_shape.end ());
01634         outer_stride.init (m_stride.begin () + M, m_stride.end ());
01635     }
01636     return MultiArrayView <N-M, T, StridedArrayTag>
01637         (outer_shape, outer_stride, ptr);
01638 }
01639 
01640 template <unsigned int N, class T, class C>
01641 template <unsigned int M>
01642 MultiArrayView <N-1, T, typename detail::MaybeStrided <M>::type >
01643 MultiArrayView <N, T, C>::bind (difference_type_1 d) const
01644 {
01645     static const int NNew = (N-1 == 0) ? 1 : N-1;
01646     TinyVector <MultiArrayIndex, NNew> shape, stride;
01647     // the remaining dimensions are 0..n-1,n+1..N-1
01648     if (N-1 == 0)
01649     {
01650         shape[0] = 1;
01651         stride[0] = 0;
01652     }
01653     else
01654     {
01655         std::copy (m_shape.begin (), m_shape.begin () + M, shape.begin ());
01656         std::copy (m_shape.begin () + M+1, m_shape.end (),
01657                    shape.begin () + M);
01658         std::copy (m_stride.begin (), m_stride.begin () + M, stride.begin ());
01659         std::copy (m_stride.begin () + M+1, m_stride.end (),
01660                    stride.begin () + M);
01661     }
01662     return MultiArrayView <N-1, T, typename detail::MaybeStrided <M>::type>
01663         (shape, stride, m_ptr + d * m_stride[M]);
01664 }
01665 
01666 template <unsigned int N, class T, class C>
01667 MultiArrayView <N - 1, T, C>
01668 MultiArrayView <N, T, C>::bindOuter (difference_type_1 d) const
01669 {
01670     static const int NNew = (N-1 == 0) ? 1 : N-1;
01671     TinyVector <MultiArrayIndex, NNew> inner_shape, inner_stride;
01672     if (N-1 == 0)
01673     {
01674         inner_shape [0] = 1;
01675         inner_stride [0] = 0;
01676     }
01677     else
01678     {
01679         inner_shape.init (m_shape.begin (), m_shape.end () - 1);
01680         inner_stride.init (m_stride.begin (), m_stride.end () - 1);
01681     }
01682     return MultiArrayView <N-1, T, C> (inner_shape, inner_stride,
01683                                        m_ptr + d * m_stride [N-1]);
01684 }
01685 
01686 template <unsigned int N, class T, class C>
01687 MultiArrayView <N - 1, T, StridedArrayTag>
01688 MultiArrayView <N, T, C>::bindInner (difference_type_1 d) const
01689 {
01690     static const int NNew = (N-1 == 0) ? 1 : N-1;
01691     TinyVector <MultiArrayIndex, NNew> outer_shape, outer_stride;
01692     if (N-1 == 0)
01693     {
01694         outer_shape [0] = 1;
01695         outer_stride [0] = 0;
01696     }
01697     else
01698     {
01699         outer_shape.init (m_shape.begin () + 1, m_shape.end ());
01700         outer_stride.init (m_stride.begin () + 1, m_stride.end ());
01701     }
01702     return MultiArrayView <N-1, T, StridedArrayTag>
01703         (outer_shape, outer_stride, m_ptr + d * m_stride [0]);
01704 }
01705 
01706 template <unsigned int N, class T, class C>
01707 MultiArrayView <N - 1, T, StridedArrayTag>
01708 MultiArrayView <N, T, C>::bindAt (difference_type_1 n, difference_type_1 d) const
01709 {
01710     vigra_precondition (
01711         n < static_cast <int> (N),
01712         "MultiArrayView <N, T, C>::bindAt(): dimension out of range.");
01713     static const int NNew = (N-1 == 0) ? 1 : N-1;
01714     TinyVector <MultiArrayIndex, NNew> shape, stride;
01715     // the remaining dimensions are 0..n-1,n+1..N-1
01716     if (N-1 == 0)
01717     {
01718         shape [0] = 1;
01719         stride [0] = 0;
01720     }
01721     else
01722     {
01723         std::copy (m_shape.begin (), m_shape.begin () + n, shape.begin ());
01724         std::copy (m_shape.begin () + n+1, m_shape.end (),
01725                    shape.begin () + n);
01726         std::copy (m_stride.begin (), m_stride.begin () + n, stride.begin ());
01727         std::copy (m_stride.begin () + n+1, m_stride.end (),
01728                    stride.begin () + n);
01729     }
01730     return MultiArrayView <N-1, T, StridedArrayTag>
01731         (shape, stride, m_ptr + d * m_stride[n]);
01732 }
01733 
01734 template <unsigned int N, class T, class C>
01735 MultiArrayView <N + 1, T, C>
01736 MultiArrayView <N, T, C>::insertSingletonDimension (difference_type_1 i) const
01737 {
01738     vigra_precondition (
01739         0 <= i && i <= static_cast <difference_type_1> (N),
01740         "MultiArrayView <N, T, C>::insertSingletonDimension(): index out of range.");
01741     TinyVector <MultiArrayIndex, N+1> shape, stride;
01742     std::copy (m_shape.begin (), m_shape.begin () + i, shape.begin ());
01743     std::copy (m_shape.begin () + i, m_shape.end (), shape.begin () + i + 1);
01744     std::copy (m_stride.begin (), m_stride.begin () + i, stride.begin ());
01745     std::copy (m_stride.begin () + i, m_stride.end (), stride.begin () + i + 1);
01746     shape[i] = 1;
01747     stride[i] = 1;
01748 
01749     return MultiArrayView <N+1, T, C>(shape, stride, m_ptr);
01750 }
01751 
01752 template <unsigned int N, class T, class C>
01753 typename NormTraits<MultiArrayView <N, T, C> >::NormType
01754 MultiArrayView <N, T, C>::norm(int type, bool useSquaredNorm) const
01755 {
01756     typedef typename NormTraits<MultiArrayView>::NormType NormType;
01757 
01758     switch(type)
01759     {
01760       case 0:
01761       {
01762         NormType res = NumericTraits<NormType>::zero();
01763         detail::normMaxOfMultiArray(traverser_begin(), shape(), res, MetaInt<actual_dimension-1>());
01764         return res;
01765       }
01766       case 1:
01767       {
01768         NormType res = NumericTraits<NormType>::zero();
01769         detail::sumOverMultiArray(traverser_begin(), shape(), detail::MultiArrayL1Functor<NormType>(),
01770                                 res, MetaInt<actual_dimension-1>());
01771         return res;
01772       }
01773       case 2:
01774       {
01775         if(useSquaredNorm)
01776         {
01777             return sqrt((NormType)squaredNorm());
01778         }
01779         else
01780         {
01781             NormType normMax = NumericTraits<NormType>::zero();
01782             detail::normMaxOfMultiArray(traverser_begin(), shape(), normMax, MetaInt<actual_dimension-1>());
01783             if(normMax == NumericTraits<NormType>::zero())
01784                 return normMax;
01785             NormType res  = NumericTraits<NormType>::zero();
01786             detail::sumOverMultiArray(traverser_begin(), shape(), detail::MultiArrayScaledL2Functor<NormType>(normMax),
01787                                     res, MetaInt<actual_dimension-1>());
01788             return sqrt(res)*normMax;
01789         }
01790       }
01791       default:
01792         vigra_precondition(false, "MultiArrayView::norm(): Unknown norm type.");
01793         return NumericTraits<NormType>::zero(); // unreachable
01794     }
01795 }
01796 
01797 
01798 /********************************************************/
01799 /*                                                      */
01800 /*                          norm                        */
01801 /*                                                      */
01802 /********************************************************/
01803 
01804 template <unsigned int N, class T, class C>
01805 inline typename NormTraits<MultiArrayView <N, T, C> >::SquaredNormType
01806 squaredNorm(MultiArrayView <N, T, C> const & a)
01807 {
01808     return a.squaredNorm();
01809 }
01810 
01811 template <unsigned int N, class T, class C>
01812 inline typename NormTraits<MultiArrayView <N, T, C> >::NormType
01813 norm(MultiArrayView <N, T, C> const & a)
01814 {
01815     return a.norm();
01816 }
01817 
01818 /********************************************************/
01819 /*                                                      */
01820 /*                       MultiArray                     */
01821 /*                                                      */
01822 /********************************************************/
01823 
01824 /** \brief Main <TT>MultiArray</TT> class containing the memory
01825     management.
01826 
01827 This class inherits the interface of MultiArrayView, and implements
01828 the memory ownership.
01829 MultiArray's are always unstrided, striding them creates a MultiArrayView.
01830 
01831 
01832 The template parameters are as follows
01833 \code
01834     N: the array dimension
01835 
01836     T: the type of the array elements
01837 
01838     A: the allocator used for internal storage management
01839        (default: std::allocator<T>)
01840 \endcode
01841 
01842 <b>\#include</b>
01843 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
01844 
01845 Namespace: vigra
01846 */
01847 template <unsigned int N, class T, class A /* default already declared above */>
01848 class MultiArray : public MultiArrayView <N, T>
01849 {
01850 
01851 public:
01852     using MultiArrayView <N, T>::actual_dimension;
01853 
01854         /** the allocator type used to allocate the memory
01855          */
01856     typedef A allocator_type;
01857 
01858         /** the view type associated with this array.
01859          */
01860     typedef MultiArrayView <N, T> view_type;
01861 
01862         /** the matrix type associated with this array.
01863          */
01864     typedef MultiArray <N, T, A> matrix_type;
01865 
01866         /** the array's value type
01867          */
01868     typedef typename view_type::value_type value_type;
01869 
01870         /** pointer type
01871          */
01872     typedef typename view_type::pointer pointer;
01873 
01874         /** const pointer type
01875          */
01876     typedef typename view_type::const_pointer const_pointer;
01877 
01878         /** reference type (result of operator[])
01879          */
01880     typedef typename view_type::reference reference;
01881 
01882         /** const reference type (result of operator[] const)
01883          */
01884     typedef typename view_type::const_reference const_reference;
01885 
01886         /** size type
01887          */
01888     typedef typename view_type::size_type size_type;
01889 
01890         /** difference type (used for multi-dimensional offsets and indices)
01891          */
01892     typedef typename view_type::difference_type difference_type;
01893 
01894         /** difference and index type for a single dimension
01895          */
01896     typedef typename view_type::difference_type_1 difference_type_1;
01897 
01898         /** traverser type
01899          */
01900     typedef typename vigra::detail::MultiIteratorChooser <
01901         UnstridedArrayTag>::template Traverser <N, T, T &, T *>::type
01902     traverser;
01903 
01904         /** traverser type to const data
01905          */
01906     typedef typename vigra::detail::MultiIteratorChooser <
01907         UnstridedArrayTag>::template Traverser <N, T, T const &, T const *>::type
01908     const_traverser;
01909 
01910         /** sequential (random access) iterator type
01911          */
01912     typedef T * iterator;
01913 
01914         /** sequential (random access) const iterator type
01915          */
01916     typedef T * const_iterator;
01917 
01918 protected:
01919 
01920     typedef typename difference_type::value_type diff_zero_t;
01921 
01922         /** the allocator used to allocate the memory
01923          */
01924     allocator_type m_alloc;
01925 
01926         /** allocate memory for s pixels, write its address into the given
01927             pointer and initialize the pixels with init.
01928         */
01929     void allocate (pointer &ptr, difference_type_1 s, const_reference init);
01930 
01931         /** allocate memory for s pixels, write its address into the given
01932             pointer and initialize the linearized pixels to the values of init.
01933         */
01934     template <class U>
01935     void allocate (pointer &ptr, difference_type_1 s, U const * init);
01936 
01937         /** allocate memory, write its address into the given
01938             pointer and initialize it by copying the data from the given MultiArrayView.
01939         */
01940     template <class U, class C>
01941     void allocate (pointer &ptr, MultiArrayView<N, U, C> const & init);
01942 
01943         /** deallocate the memory (of length s) starting at the given address.
01944          */
01945     void deallocate (pointer &ptr, difference_type_1 s);
01946 
01947     template <class U, class C>
01948     void copyOrReshape (const MultiArrayView<N, U, C> &rhs);
01949 public:
01950         /** default constructor
01951          */
01952     MultiArray ()
01953     : MultiArrayView <N, T> (difference_type (diff_zero_t(0)),
01954                              difference_type (diff_zero_t(0)), 0)
01955     {}
01956 
01957         /** construct with given allocator
01958          */
01959     MultiArray (allocator_type const & alloc)
01960     : MultiArrayView <N, T> (difference_type (diff_zero_t(0)),
01961                              difference_type (diff_zero_t(0)), 0),
01962       m_alloc(alloc)
01963     {}
01964 
01965         /** construct with given shape
01966          */
01967     explicit MultiArray (const difference_type &shape,
01968                          allocator_type const & alloc = allocator_type());
01969 
01970         /** construct from shape with an initial value
01971          */
01972     MultiArray (const difference_type &shape, const_reference init,
01973                 allocator_type const & alloc = allocator_type());
01974 
01975         /** construct from shape and copy values from the given array
01976          */
01977     MultiArray (const difference_type &shape, const_pointer init,
01978                          allocator_type const & alloc = allocator_type());
01979 
01980         /** copy constructor
01981          */
01982     MultiArray (const MultiArray &rhs)
01983     : MultiArrayView <N, T> (rhs.m_shape, rhs.m_stride, 0),
01984       m_alloc (rhs.m_alloc)
01985     {
01986         allocate (this->m_ptr, this->elementCount (), rhs.data ());
01987     }
01988 
01989         /** construct by copying from a MultiArrayView
01990          */
01991     template <class U, class C>
01992     MultiArray (const MultiArrayView<N, U, C>  &rhs,
01993                 allocator_type const & alloc = allocator_type());
01994 
01995         /** assignment.<br>
01996             If the size of \a rhs is the same as the left-hand side arrays's old size, only
01997             the data are copied. Otherwise, new storage is allocated, which invalidates all
01998             objects (array views, iterators) depending on the lhs array.
01999          */
02000     MultiArray & operator= (const MultiArray &rhs)
02001     {
02002         if (this != &rhs)
02003             this->copyOrReshape(rhs);
02004         return *this;
02005     }
02006 
02007         /** assignment from arbitrary MultiArrayView.<br>
02008             If the size of \a rhs is the same as the left-hand side arrays's old size, only
02009             the data are copied. Otherwise, new storage is allocated, which invalidates all
02010             objects (array views, iterators) depending on the lhs array.
02011          */
02012     template <class U, class C>
02013     MultiArray &operator= (const MultiArrayView<N, U, C> &rhs)
02014     {
02015         this->copyOrReshape(rhs);
02016         return *this;
02017     }
02018 
02019         /** Add-assignment from arbitrary MultiArrayView. Fails with
02020             <tt>PreconditionViolation</tt> exception when the shapes do not match.
02021          */
02022     template <class U, class C>
02023     MultiArray &operator+= (const MultiArrayView<N, U, C> &rhs)
02024     {
02025         view_type::operator+=(rhs);
02026         return *this;
02027     }
02028 
02029         /** Subtract-assignment from arbitrary MultiArrayView. Fails with
02030             <tt>PreconditionViolation</tt> exception when the shapes do not match.
02031          */
02032     template <class U, class C>
02033     MultiArray &operator-= (const MultiArrayView<N, U, C> &rhs)
02034     {
02035         view_type::operator-=(rhs);
02036         return *this;
02037     }
02038 
02039         /** Multiply-assignment from arbitrary MultiArrayView. Fails with
02040             <tt>PreconditionViolation</tt> exception when the shapes do not match.
02041          */
02042     template <class U, class C>
02043     MultiArray &operator*= (const MultiArrayView<N, U, C> &rhs)
02044     {
02045         view_type::operator*=(rhs);
02046         return *this;
02047     }
02048 
02049         /** Divide-assignment from arbitrary MultiArrayView. Fails with
02050             <tt>PreconditionViolation</tt> exception when the shapes do not match.
02051          */
02052     template <class U, class C>
02053     MultiArray &operator/= (const MultiArrayView<N, U, C> &rhs)
02054     {
02055         view_type::operator/=(rhs);
02056         return *this;
02057     }
02058 
02059         /** Add-assignment of a scalar.
02060          */
02061     MultiArray &operator+= (const T &rhs)
02062     {
02063         view_type::operator+=(rhs);
02064         return *this;
02065     }
02066 
02067         /** Subtract-assignment of a scalar.
02068          */
02069     MultiArray &operator-= (const T &rhs)
02070     {
02071         view_type::operator-=(rhs);
02072         return *this;
02073     }
02074 
02075         /** Multiply-assignment of a scalar.
02076          */
02077     MultiArray &operator*= (const T &rhs)
02078     {
02079         view_type::operator*=(rhs);
02080         return *this;
02081     }
02082 
02083         /** Divide-assignment of a scalar.
02084          */
02085     MultiArray &operator/= (const T &rhs)
02086     {
02087         view_type::operator/=(rhs);
02088         return *this;
02089     }
02090 
02091         /** destructor
02092          */
02093     ~MultiArray ()
02094     {
02095         deallocate (this->m_ptr, this->elementCount ());
02096     }
02097 
02098 
02099         /** init elements with a constant
02100          */
02101     template <class U>
02102     MultiArray & init(const U & init)
02103     {
02104         view_type::init(init);
02105         return *this;
02106     }
02107 
02108         /** Allocate new memory with the given shape and initialize with zeros.<br>
02109             <em>Note:</em> this operation invalidates all dependent objects
02110             (array views and iterators)
02111          */
02112     void reshape (const difference_type &shape)
02113     {
02114         reshape (shape, NumericTraits <T>::zero ());
02115     }
02116 
02117         /** Allocate new memory with the given shape and initialize it
02118             with the given value.<br>
02119             <em>Note:</em> this operation invalidates all dependent objects
02120             (array views and iterators)
02121          */
02122     void reshape (const difference_type &shape, const_reference init);
02123 
02124         /** Swap the contents with another MultiArray. This is fast,
02125             because no data are copied, but only pointers and shapes swapped.
02126             <em>Note:</em> this operation invalidates all dependent objects
02127             (array views and iterators)
02128          */
02129     void swap (MultiArray & other);
02130 
02131         /** sequential iterator pointing to the first array element.
02132          */
02133     iterator begin ()
02134     {
02135         return this->data();
02136     }
02137 
02138         /** sequential iterator pointing beyond the last array element.
02139          */
02140     iterator end ()
02141     {
02142         return this->data() + this->elementCount();
02143     }
02144 
02145         /** sequential const iterator pointing to the first array element.
02146          */
02147     const_iterator begin () const
02148     {
02149         return this->data();
02150     }
02151 
02152         /** sequential const iterator pointing beyond the last array element.
02153          */
02154     const_iterator end () const
02155     {
02156         return this->data() + this->elementCount();
02157     }
02158 
02159         /** get the allocator.
02160          */
02161     allocator_type const & allocator () const
02162     {
02163         return m_alloc;
02164     }
02165 };
02166 
02167 template <unsigned int N, class T, class A>
02168 MultiArray <N, T, A>::MultiArray (const difference_type &shape,
02169                                   allocator_type const & alloc)
02170 : MultiArrayView <N, T> (shape,
02171                          detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
02172                          0),
02173   m_alloc(alloc)
02174 {
02175     if (N == 0)
02176     {
02177         this->m_shape [0] = 1;
02178         this->m_stride [0] = 0;
02179     }
02180     allocate (this->m_ptr, this->elementCount (), NumericTraits<T>::zero ());
02181 }
02182 
02183 template <unsigned int N, class T, class A>
02184 MultiArray <N, T, A>::MultiArray (const difference_type &shape, const_reference init,
02185                                   allocator_type const & alloc)
02186 : MultiArrayView <N, T> (shape,
02187                          detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
02188                          0),
02189   m_alloc(alloc)
02190 {
02191     if (N == 0)
02192     {
02193         this->m_shape [0] = 1;
02194         this->m_stride [0] = 0;
02195     }
02196     allocate (this->m_ptr, this->elementCount (), init);
02197 }
02198 
02199 template <unsigned int N, class T, class A>
02200 MultiArray <N, T, A>::MultiArray (const difference_type &shape, const_pointer init,
02201                                   allocator_type const & alloc)
02202 : MultiArrayView <N, T> (shape,
02203                          detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
02204                          0),
02205   m_alloc(alloc)
02206 {
02207     if (N == 0)
02208     {
02209         this->m_shape [0] = 1;
02210         this->m_stride [0] = 0;
02211     }
02212     allocate (this->m_ptr, this->elementCount (), init);
02213 }
02214 
02215 template <unsigned int N, class T, class A>
02216 template <class U, class C>
02217 MultiArray <N, T, A>::MultiArray(const MultiArrayView<N, U, C>  &rhs,
02218                                  allocator_type const & alloc)
02219 : MultiArrayView <N, T> (rhs.shape(),
02220                          detail::defaultStride <MultiArrayView<N,T>::actual_dimension>(rhs.shape()),
02221                          0),
02222   m_alloc (alloc)
02223 {
02224     allocate (this->m_ptr, rhs);
02225 }
02226 
02227 template <unsigned int N, class T, class A>
02228 template <class U, class C>
02229 void
02230 MultiArray <N, T, A>::copyOrReshape(const MultiArrayView<N, U, C> &rhs)
02231 {
02232     if (this->shape() == rhs.shape())
02233         this->copy(rhs);
02234     else
02235     {
02236         MultiArray t(rhs);
02237         this->swap(t);
02238     }
02239 }
02240 
02241 template <unsigned int N, class T, class A>
02242 void MultiArray <N, T, A>::reshape (const difference_type & new_shape,
02243                                     const_reference initial)
02244 {
02245     if (N== 0)
02246     {
02247         return;
02248     }
02249     else if(new_shape == this->shape())
02250     {
02251         this->init(initial);
02252     }
02253     else
02254     {
02255         difference_type new_stride = detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (new_shape);
02256         difference_type_1 new_size = new_shape [MultiArrayView<N,T>::actual_dimension-1] * new_stride [MultiArrayView<N,T>::actual_dimension-1];
02257         T *new_ptr;
02258         allocate (new_ptr, new_size, initial);
02259         deallocate (this->m_ptr, this->elementCount ());
02260         this->m_ptr = new_ptr;
02261         this->m_shape = new_shape;
02262         this->m_stride = new_stride;
02263     }
02264 }
02265 
02266 
02267 template <unsigned int N, class T, class A>
02268 inline void
02269 MultiArray <N, T, A>::swap (MultiArray <N, T, A> & other)
02270 {
02271     if (this == &other)
02272         return;
02273     std::swap(this->m_shape,  other.m_shape);
02274     std::swap(this->m_stride, other.m_stride);
02275     std::swap(this->m_ptr,    other.m_ptr);
02276     std::swap(this->m_alloc,  other.m_alloc);
02277 }
02278 
02279 template <unsigned int N, class T, class A>
02280 void MultiArray <N, T, A>::allocate (pointer & ptr, difference_type_1 s,
02281                                      const_reference init)
02282 {
02283     ptr = m_alloc.allocate ((typename A::size_type)s);
02284     difference_type_1 i;
02285     try {
02286         for (i = 0; i < s; ++i)
02287             m_alloc.construct (ptr + i, init);
02288     }
02289     catch (...) {
02290         for (difference_type_1 j = 0; j < i; ++j)
02291             m_alloc.destroy (ptr + j);
02292         m_alloc.deallocate (ptr, (typename A::size_type)s);
02293         throw;
02294     }
02295 }
02296 
02297 template <unsigned int N, class T, class A>
02298 template <class U>
02299 void MultiArray <N, T, A>::allocate (pointer & ptr, difference_type_1 s,
02300                                      U const * init)
02301 {
02302     ptr = m_alloc.allocate ((typename A::size_type)s);
02303     difference_type_1 i;
02304     try {
02305         for (i = 0; i < s; ++i, ++init)
02306             m_alloc.construct (ptr + i, *init);
02307     }
02308     catch (...) {
02309         for (difference_type_1 j = 0; j < i; ++j)
02310             m_alloc.destroy (ptr + j);
02311         m_alloc.deallocate (ptr, (typename A::size_type)s);
02312         throw;
02313     }
02314 }
02315 
02316 template <unsigned int N, class T, class A>
02317 template <class U, class C>
02318 void MultiArray <N, T, A>::allocate (pointer & ptr, MultiArrayView<N, U, C> const & init)
02319 {
02320     difference_type_1 s = init.elementCount();
02321     ptr = m_alloc.allocate ((typename A::size_type)s);
02322     pointer p = ptr;
02323     try {
02324         detail::uninitializedCopyMultiArrayData(init.traverser_begin(), init.shape(),
02325                                                 p, m_alloc, MetaInt<actual_dimension-1>());
02326     }
02327     catch (...) {
02328         for (pointer pp = ptr; pp < p; ++pp)
02329             m_alloc.destroy (pp);
02330         m_alloc.deallocate (ptr, (typename A::size_type)s);
02331         throw;
02332     }
02333 }
02334 
02335 template <unsigned int N, class T, class A>
02336 inline void MultiArray <N, T, A>::deallocate (pointer & ptr, difference_type_1 s)
02337 {
02338     if (ptr == 0)
02339         return;
02340     for (difference_type_1 i = 0; i < s; ++i)
02341         m_alloc.destroy (ptr + i);
02342     m_alloc.deallocate (ptr, (typename A::size_type)s);
02343     ptr = 0;
02344 }
02345 
02346 /********************************************************/
02347 /*                                                      */
02348 /*              argument object factories               */
02349 /*                                                      */
02350 /********************************************************/
02351 
02352 template <unsigned int N, class T, class C>
02353 inline triple<typename MultiArrayView<N,T,C>::const_traverser,
02354               typename MultiArrayView<N,T,C>::difference_type,
02355               typename AccessorTraits<T>::default_const_accessor >
02356 srcMultiArrayRange( MultiArrayView<N,T,C> const & array )
02357 {
02358     return triple<typename MultiArrayView<N,T,C>::const_traverser,
02359                   typename MultiArrayView<N,T,C>::difference_type,
02360                   typename AccessorTraits<T>::default_const_accessor >
02361       ( array.traverser_begin(),
02362         array.shape(),
02363         typename AccessorTraits<T>::default_const_accessor() );
02364 }
02365 
02366 template <unsigned int N, class T, class C, class Accessor>
02367 inline triple<typename MultiArrayView<N,T,C>::const_traverser,
02368               typename MultiArrayView<N,T,C>::difference_type,
02369               Accessor >
02370 srcMultiArrayRange( MultiArrayView<N,T,C> const & array, Accessor a )
02371 {
02372     return triple<typename MultiArrayView<N,T,C>::const_traverser,
02373                   typename MultiArrayView<N,T,C>::difference_type,
02374                   Accessor >
02375       ( array.traverser_begin(),
02376         array.shape(),
02377         a);
02378 }
02379 
02380 template <unsigned int N, class T, class C>
02381 inline pair<typename MultiArrayView<N,T,C>::const_traverser,
02382             typename AccessorTraits<T>::default_const_accessor >
02383 srcMultiArray( MultiArrayView<N,T,C> const & array )
02384 {
02385     return pair<typename MultiArrayView<N,T,C>::const_traverser,
02386                 typename AccessorTraits<T>::default_const_accessor >
02387       ( array.traverser_begin(),
02388         typename AccessorTraits<T>::default_const_accessor() );
02389 }
02390 
02391 template <unsigned int N, class T, class C, class Accessor>
02392 inline pair<typename MultiArrayView<N,T,C>::const_traverser,
02393             Accessor >
02394 srcMultiArray( MultiArrayView<N,T,C> const & array, Accessor a )
02395 {
02396     return pair<typename MultiArrayView<N,T,C>::const_traverser,
02397                 Accessor >
02398       ( array.traverser_begin(), a );
02399 }
02400 
02401 template <unsigned int N, class T, class C>
02402 inline triple<typename MultiArrayView<N,T,C>::traverser,
02403               typename MultiArrayView<N,T,C>::difference_type,
02404               typename AccessorTraits<T>::default_accessor >
02405 destMultiArrayRange( MultiArrayView<N,T,C> & array )
02406 {
02407     return triple<typename MultiArrayView<N,T,C>::traverser,
02408                   typename MultiArrayView<N,T,C>::difference_type,
02409                   typename AccessorTraits<T>::default_accessor >
02410       ( array.traverser_begin(),
02411         array.shape(),
02412         typename AccessorTraits<T>::default_accessor() );
02413 }
02414 
02415 template <unsigned int N, class T, class C, class Accessor>
02416 inline triple<typename MultiArrayView<N,T,C>::traverser,
02417               typename MultiArrayView<N,T,C>::difference_type,
02418               Accessor >
02419 destMultiArrayRange( MultiArrayView<N,T,C> & array, Accessor a )
02420 {
02421     return triple<typename MultiArrayView<N,T,C>::traverser,
02422                   typename MultiArrayView<N,T,C>::difference_type,
02423                   Accessor >
02424       ( array.traverser_begin(),
02425         array.shape(),
02426         a );
02427 }
02428 
02429 template <unsigned int N, class T, class C>
02430 inline pair<typename MultiArrayView<N,T,C>::traverser,
02431             typename AccessorTraits<T>::default_accessor >
02432 destMultiArray( MultiArrayView<N,T,C> & array )
02433 {
02434     return pair<typename MultiArrayView<N,T,C>::traverser,
02435                 typename AccessorTraits<T>::default_accessor >
02436         ( array.traverser_begin(),
02437           typename AccessorTraits<T>::default_accessor() );
02438 }
02439 
02440 template <unsigned int N, class T, class C, class Accessor>
02441 inline pair<typename MultiArrayView<N,T,C>::traverser,
02442             Accessor >
02443 destMultiArray( MultiArrayView<N,T,C> & array, Accessor a )
02444 {
02445     return pair<typename MultiArrayView<N,T,C>::traverser,
02446                 Accessor >
02447         ( array.traverser_begin(), a );
02448 }
02449 
02450 /********************************************************************/
02451 
02452 template <class PixelType, class Accessor>
02453 inline triple<ConstStridedImageIterator<PixelType>,
02454               ConstStridedImageIterator<PixelType>, Accessor>
02455 srcImageRange(const MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
02456 {
02457     ConstStridedImageIterator<PixelType>
02458         ul(img.data(), 1, img.stride(0), img.stride(1));
02459     return triple<ConstStridedImageIterator<PixelType>,
02460                   ConstStridedImageIterator<PixelType>,
02461                   Accessor>(
02462                       ul, ul + Size2D(img.shape(0), img.shape(1)), a);
02463 }
02464 
02465 template <class PixelType, class Accessor>
02466 inline pair<ConstStridedImageIterator<PixelType>, Accessor>
02467 srcImage(const MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
02468 {
02469     ConstStridedImageIterator<PixelType>
02470         ul(img.data(), 1, img.stride(0), img.stride(1));
02471     return pair<ConstStridedImageIterator<PixelType>, Accessor>
02472         (ul, a);
02473 }
02474 
02475 template <class PixelType, class Accessor>
02476 inline triple<StridedImageIterator<PixelType>,
02477               StridedImageIterator<PixelType>, Accessor>
02478 destImageRange(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
02479 {
02480     StridedImageIterator<PixelType>
02481         ul(img.data(), 1, img.stride(0), img.stride(1));
02482     return triple<StridedImageIterator<PixelType>,
02483                   StridedImageIterator<PixelType>,
02484                   Accessor>(
02485                       ul, ul + Size2D(img.shape(0), img.shape(1)), a);
02486 }
02487 
02488 template <class PixelType, class Accessor>
02489 inline pair<StridedImageIterator<PixelType>, Accessor>
02490 destImage(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
02491 {
02492     StridedImageIterator<PixelType>
02493         ul(img.data(), 1, img.stride(0), img.stride(1));
02494     return pair<StridedImageIterator<PixelType>, Accessor>
02495         (ul, a);
02496 }
02497 
02498 template <class PixelType, class Accessor>
02499 inline pair<StridedImageIterator<PixelType>, Accessor>
02500 maskImage(MultiArrayView<2, PixelType, StridedArrayTag> & img, Accessor a)
02501 {
02502     StridedImageIterator<PixelType>
02503         ul(img.data(), 1, img.stride(0), img.stride(1));
02504     return pair<StridedImageIterator<PixelType>, Accessor>
02505         (ul, a);
02506 }
02507 
02508 // -------------------------------------------------------------------
02509 
02510 template <class PixelType>
02511 inline triple<ConstStridedImageIterator<PixelType>,
02512               ConstStridedImageIterator<PixelType>,
02513               typename AccessorTraits<PixelType>::default_const_accessor>
02514 srcImageRange(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
02515 {
02516     ConstStridedImageIterator<PixelType>
02517         ul(img.data(), 1, img.stride(0), img.stride(1));
02518     typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
02519     return triple<ConstStridedImageIterator<PixelType>,
02520                   ConstStridedImageIterator<PixelType>,
02521                   Accessor>
02522         (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
02523 }
02524 
02525 template <class PixelType>
02526 inline triple<ConstImageIterator<PixelType>,
02527               ConstImageIterator<PixelType>,
02528               typename AccessorTraits<PixelType>::default_const_accessor>
02529 srcImageRange(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
02530 {
02531     ConstImageIterator<PixelType>
02532         ul(img.data(), img.stride(1));
02533     typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
02534     return triple<ConstImageIterator<PixelType>,
02535                   ConstImageIterator<PixelType>,
02536                   Accessor>
02537         (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
02538 }
02539 
02540 template <class PixelType>
02541 inline pair< ConstStridedImageIterator<PixelType>,
02542              typename AccessorTraits<PixelType>::default_const_accessor>
02543 srcImage(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
02544 {
02545     ConstStridedImageIterator<PixelType>
02546         ul(img.data(), 1, img.stride(0), img.stride(1));
02547     typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
02548     return pair<ConstStridedImageIterator<PixelType>,
02549                 Accessor>
02550         (ul, Accessor());
02551 }
02552 
02553 template <class PixelType>
02554 inline pair< ConstImageIterator<PixelType>,
02555              typename AccessorTraits<PixelType>::default_const_accessor>
02556 srcImage(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
02557 {
02558     ConstImageIterator<PixelType>
02559         ul(img.data(), img.stride(1));
02560     typedef typename AccessorTraits<PixelType>::default_const_accessor Accessor;
02561     return pair<ConstImageIterator<PixelType>,
02562                 Accessor>
02563         (ul, Accessor());
02564 }
02565 
02566 template <class PixelType>
02567 inline triple< StridedImageIterator<PixelType>,
02568                StridedImageIterator<PixelType>,
02569                typename AccessorTraits<PixelType>::default_accessor>
02570 destImageRange(MultiArrayView<2, PixelType, StridedArrayTag> & img)
02571 {
02572     StridedImageIterator<PixelType>
02573         ul(img.data(), 1, img.stride(0), img.stride(1));
02574     typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
02575     return triple<StridedImageIterator<PixelType>,
02576                   StridedImageIterator<PixelType>,
02577                   Accessor>
02578         (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
02579 }
02580 
02581 template <class PixelType>
02582 inline triple< ImageIterator<PixelType>,
02583                ImageIterator<PixelType>,
02584                typename AccessorTraits<PixelType>::default_accessor>
02585 destImageRange(MultiArrayView<2, PixelType, UnstridedArrayTag> & img)
02586 {
02587     ImageIterator<PixelType>
02588         ul(img.data(), img.stride(1));
02589     typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
02590     return triple<ImageIterator<PixelType>,
02591                   ImageIterator<PixelType>,
02592                   Accessor>
02593         (ul, ul + Size2D(img.shape(0), img.shape(1)), Accessor());
02594 }
02595 
02596 template <class PixelType>
02597 inline pair< StridedImageIterator<PixelType>,
02598              typename AccessorTraits<PixelType>::default_accessor>
02599 destImage(MultiArrayView<2, PixelType, StridedArrayTag> & img)
02600 {
02601     StridedImageIterator<PixelType>
02602         ul(img.data(), 1, img.stride(0), img.stride(1));
02603     typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
02604     return pair<StridedImageIterator<PixelType>, Accessor>
02605         (ul, Accessor());
02606 }
02607 
02608 template <class PixelType>
02609 inline pair< ImageIterator<PixelType>,
02610              typename AccessorTraits<PixelType>::default_accessor>
02611 destImage(MultiArrayView<2, PixelType, UnstridedArrayTag> & img)
02612 {
02613     ImageIterator<PixelType> ul(img.data(), img.stride(1));
02614     typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
02615     return pair<ImageIterator<PixelType>, Accessor>(ul, Accessor());
02616 }
02617 
02618 template <class PixelType>
02619 inline pair< ConstStridedImageIterator<PixelType>,
02620              typename AccessorTraits<PixelType>::default_accessor>
02621 maskImage(MultiArrayView<2, PixelType, StridedArrayTag> const & img)
02622 {
02623     ConstStridedImageIterator<PixelType>
02624         ul(img.data(), 1, img.stride(0), img.stride(1));
02625     typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
02626     return pair<ConstStridedImageIterator<PixelType>, Accessor>
02627         (ul, Accessor());
02628 }
02629 
02630 template <class PixelType>
02631 inline pair< ConstImageIterator<PixelType>,
02632              typename AccessorTraits<PixelType>::default_accessor>
02633 maskImage(MultiArrayView<2, PixelType, UnstridedArrayTag> const & img)
02634 {
02635     ConstImageIterator<PixelType>
02636         ul(img.data(), img.stride(1));
02637     typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
02638     return pair<ConstImageIterator<PixelType>, Accessor>
02639         (ul, Accessor());
02640 }
02641 
02642 /********************************************************/
02643 /*                                                      */
02644 /*                  makeBasicImageView                  */
02645 /*                                                      */
02646 /********************************************************/
02647 
02648 /** \addtogroup MultiArrayToImage Wrap a \ref vigra::MultiArrayView in
02649                                   a \ref vigra::BasicImageView
02650 */
02651 //@{
02652 /** Create a \ref vigra::BasicImageView from an unstrided 2-dimensional
02653     \ref vigra::MultiArrayView.
02654 
02655     The \ref vigra::BasicImageView will have the same <tt>value_type </tt>
02656     as the original \ref vigra::MultiArrayView.
02657 */
02658 template <class T>
02659 BasicImageView <T>
02660 makeBasicImageView (MultiArrayView <2, T, UnstridedArrayTag> const &array)
02661 {
02662     return BasicImageView <T> (array.data (), array.shape (0),
02663                                array.shape (1));
02664 }
02665 
02666 /** Create a \ref vigra::BasicImageView from a 3-dimensional
02667     \ref vigra::MultiArray.
02668 
02669     This wrapper flattens the two innermost dimensions of the array
02670     into single rows of the resulting image.
02671     The \ref vigra::BasicImageView will have the same <tt>value_type </tt>
02672     as the original \ref vigra::MultiArray.
02673 */
02674 template <class T>
02675 BasicImageView <T>
02676 makeBasicImageView (MultiArray <3, T> const &array)
02677 {
02678     return BasicImageView <T> (array.data (),
02679                                array.shape (0)*array.shape (1), array.shape (2));
02680 }
02681 
02682 /** Create a \ref vigra::BasicImageView from a 3-dimensional
02683     \ref vigra::MultiArray.
02684 
02685     This wrapper only works if <tt>T</tt> is a scalar type and the
02686     array's innermost dimension has size 3. It then re-interprets
02687     the data array as a 2-dimensional array with value_type
02688     <tt>RGBValue<T></tt>.
02689 */
02690 template <class T>
02691 BasicImageView <RGBValue<T> >
02692 makeRGBImageView (MultiArray<3, T> const &array)
02693 {
02694     vigra_precondition (
02695         array.shape (0) == 3, "makeRGBImageView(): array.shape(0) must be 3.");
02696     return BasicImageView <RGBValue<T> > (
02697         reinterpret_cast <RGBValue <T> *> (array.data ()),
02698         array.shape (1), array.shape (2));
02699 }
02700 
02701 //@}
02702 
02703 } // namespace vigra
02704 #undef VIGRA_ASSERT_INSIDE
02705 #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.7.1 (Sun Feb 19 2012)