[ 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 /*    ( Version 1.6.0, Aug 13 2008 )                                    */
00007 /*    The VIGRA Website is                                              */
00008 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
00009 /*    Please direct questions, bug reports, and contributions to        */
00010 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00011 /*        vigra@informatik.uni-hamburg.de                               */
00012 /*                                                                      */
00013 /*    Permission is hereby granted, free of charge, to any person       */
00014 /*    obtaining a copy of this software and associated documentation    */
00015 /*    files (the "Software"), to deal in the Software without           */
00016 /*    restriction, including without limitation the rights to use,      */
00017 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00018 /*    sell copies of the Software, and to permit persons to whom the    */
00019 /*    Software is furnished to do so, subject to the following          */
00020 /*    conditions:                                                       */
00021 /*                                                                      */
00022 /*    The above copyright notice and this permission notice shall be    */
00023 /*    included in all copies or substantial portions of the             */
00024 /*    Software.                                                         */
00025 /*                                                                      */
00026 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00027 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00028 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00029 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00030 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00031 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00032 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00033 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00034 /*                                                                      */
00035 /************************************************************************/
00036 
00037 
00038 #ifndef VIGRA_MULTI_ARRAY_HXX
00039 #define VIGRA_MULTI_ARRAY_HXX
00040 
00041 #include <memory>
00042 #include <algorithm>
00043 #include "accessor.hxx"
00044 #include "tinyvector.hxx"
00045 #include "rgbvalue.hxx"
00046 #include "basicimageview.hxx"
00047 #include "imageiterator.hxx"
00048 #include "numerictraits.hxx"
00049 #include "multi_iterator.hxx"
00050 #include "metaprogramming.hxx"
00051 #include "mathutil.hxx"
00052 
00053 namespace vigra
00054 {
00055 
00056 namespace detail
00057 {
00058 /********************************************************/
00059 /*                                                      */
00060 /*                    defaultStride                     */
00061 /*                                                      */
00062 /********************************************************/
00063 
00064 /* generates the stride for a gapless shape.
00065 
00066     Namespace: vigra::detail
00067 */
00068 template <unsigned int N>
00069 inline TinyVector <MultiArrayIndex, N>
00070 defaultStride(const TinyVector <MultiArrayIndex, N> &shape)
00071 {
00072     TinyVector <MultiArrayIndex, N> ret;
00073     ret [0] = 1;
00074     for (int i = 1; i < (int)N; ++i)
00075         ret [i] = ret [i-1] * shape [i-1];
00076     return ret;
00077 }
00078 
00079 /********************************************************/
00080 /*                                                      */
00081 /*                 ScanOrderToOffset                    */
00082 /*                                                      */
00083 /********************************************************/
00084 
00085 /* transforms an index in scan order sense to a pointer offset in a possibly
00086    strided, multi-dimensional array.
00087 
00088     Namespace: vigra::detail
00089 */
00090 
00091 template <int K>
00092 struct ScanOrderToOffset
00093 {
00094     template <int N>
00095     static MultiArrayIndex
00096     exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> &shape,
00097          const TinyVector <MultiArrayIndex, N> & stride)
00098     {
00099         return stride[N-K] * (d % shape[N-K]) +
00100                ScanOrderToOffset<K-1>::exec(d / shape[N-K], shape, stride);
00101     }
00102 };
00103 
00104 template <>
00105 struct ScanOrderToOffset<1>
00106 {
00107     template <int N>
00108     static MultiArrayIndex
00109     exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> & /*shape*/,
00110          const TinyVector <MultiArrayIndex, N> & stride)
00111     {
00112         return stride[N-1] * d;
00113     }
00114 };
00115 
00116 template <int K>
00117 struct ScanOrderToCoordinate
00118 {
00119     template <int N>
00120     static void
00121     exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> &shape,
00122          TinyVector <MultiArrayIndex, N> & result)
00123     {
00124         result[N-K] = (d % shape[N-K]);
00125         ScanOrderToCoordinate<K-1>::exec(d / shape[N-K], shape, result);
00126     }
00127 };
00128 
00129 template <>
00130 struct ScanOrderToCoordinate<1>
00131 {
00132     template <int N>
00133     static void
00134     exec(MultiArrayIndex d, const TinyVector <MultiArrayIndex, N> & /*shape*/,
00135          TinyVector <MultiArrayIndex, N> & result)
00136     {
00137         result[N-1] = d;
00138     }
00139 };
00140 
00141 template <int K>
00142 struct CoordinateToScanOrder
00143 {
00144     template <int N>
00145     static MultiArrayIndex
00146     exec(const TinyVector <MultiArrayIndex, N> &shape,
00147          const TinyVector <MultiArrayIndex, N> & coordinate)
00148     {
00149         return coordinate[N-K] + shape[N-K] * CoordinateToScanOrder<K-1>::exec(shape, coordinate);
00150     }
00151 };
00152 
00153 template <>
00154 struct CoordinateToScanOrder<1>
00155 {
00156     template <int N>
00157     static MultiArrayIndex
00158     exec(const TinyVector <MultiArrayIndex, N> & /*shape*/,
00159          const TinyVector <MultiArrayIndex, N> & coordinate)
00160     {
00161         return coordinate[N-1];
00162     }
00163 };
00164 
00165 
00166 template <class C>
00167 struct CoordinatesToOffest
00168 {
00169     template <int N>
00170     static MultiArrayIndex
00171     exec(const TinyVector <MultiArrayIndex, N> & stride, MultiArrayIndex x)
00172     {
00173         return stride[0] * x;
00174     }
00175     template <int N>
00176     static MultiArrayIndex
00177     exec(const TinyVector <MultiArrayIndex, N> & stride, MultiArrayIndex x, MultiArrayIndex y)
00178     {
00179         return stride[0] * x + stride[1] * y;
00180     }
00181 };
00182 
00183 template <>
00184 struct CoordinatesToOffest<UnstridedArrayTag>
00185 {
00186     template <int N>
00187     static MultiArrayIndex
00188     exec(const TinyVector <MultiArrayIndex, N> & /*stride*/, MultiArrayIndex x)
00189     {
00190         return x;
00191     }
00192     template <int N>
00193     static MultiArrayIndex
00194     exec(const TinyVector <MultiArrayIndex, N> & stride, MultiArrayIndex x, MultiArrayIndex y)
00195     {
00196         return x + stride[1] * y;
00197     }
00198 };
00199 
00200 /********************************************************/
00201 /*                                                      */
00202 /*                     MaybeStrided                     */
00203 /*                                                      */
00204 /********************************************************/
00205 
00206 /* metatag implementing a test for marking MultiArrays that were
00207     indexed at the zero'th dimension as strided, and all others as
00208     unstrided.
00209 
00210 <b>\#include</b>
00211 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
00212 
00213 Namespace: vigra::detail
00214 */
00215 template <unsigned int N>
00216 struct MaybeStrided
00217 {
00218     typedef UnstridedArrayTag type;
00219 };
00220 
00221 template <>
00222 struct MaybeStrided <0>
00223 {
00224     typedef StridedArrayTag type;
00225 };
00226 
00227 /********************************************************/
00228 /*                                                      */
00229 /*                MultiIteratorChooser                  */
00230 /*                                                      */
00231 /********************************************************/
00232 
00233 /* metatag implementing a test (by pattern matching) for marking
00234     MultiArrays that were indexed at the zero'th dimension as strided.
00235 
00236 <b>\#include</b>
00237 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
00238 
00239 Namespace: vigra::detail
00240 */
00241 template <class O>
00242 struct MultiIteratorChooser
00243 {
00244     struct Nil {};
00245 
00246     template <unsigned int N, class T, class REFERENCE, class POINTER>
00247     struct Traverser
00248     {
00249         typedef Nil type;
00250     };
00251 };
00252 
00253 /********************************************************/
00254 /*                                                      */
00255 /*       MultiIteratorChooser <StridedArrayTag>         */
00256 /*                                                      */
00257 /********************************************************/
00258 
00259 /* specialization of the MultiIteratorChooser for strided arrays.
00260 
00261 <b>\#include</b>
00262 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
00263 
00264 Namespace: vigra::detail
00265 */
00266 template <>
00267 struct MultiIteratorChooser <StridedArrayTag>
00268 {
00269     template <unsigned int N, class T, class REFERENCE, class POINTER>
00270     struct Traverser
00271     {
00272         typedef StridedMultiIterator <N, T, REFERENCE, POINTER> type;
00273     };
00274 };
00275 
00276 /********************************************************/
00277 /*                                                      */
00278 /*      MultiIteratorChooser <UnstridedArrayTag>        */
00279 /*                                                      */
00280 /********************************************************/
00281 
00282 /* specialization of the MultiIteratorChooser for unstrided arrays.
00283 
00284 <b>\#include</b>
00285 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
00286 
00287 Namespace: vigra::detail
00288 */
00289 template <>
00290 struct MultiIteratorChooser <UnstridedArrayTag>
00291 {
00292     template <unsigned int N, class T, class REFERENCE, class POINTER>
00293     struct Traverser
00294     {
00295         typedef MultiIterator <N, T, REFERENCE, POINTER> type;
00296     };
00297 };
00298 
00299 /********************************************************/
00300 /*                                                      */
00301 /*                   helper functions                   */
00302 /*                                                      */
00303 /********************************************************/
00304 
00305 template <class DestIterator, class Shape, class T>
00306 inline void
00307 initMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<0>)
00308 {
00309     DestIterator dend = d + shape[0];
00310     for(; d < dend; ++d)
00311     {
00312         *d = init;
00313     }
00314 }
00315 
00316 template <class DestIterator, class Shape, class T, int N>
00317 void
00318 initMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<N>)
00319 {
00320     DestIterator dend = d + shape[N];
00321     for(; d < dend; ++d)
00322     {
00323         initMultiArrayData(d.begin(), shape, init, MetaInt<N-1>());
00324     }
00325 }
00326 
00327 #define VIGRA_COPY_MULTI_ARRAY_DATA(name, op) \
00328 template <class SrcIterator, class Shape, class DestIterator> \
00329 inline void \
00330 name##MultiArrayData(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>) \
00331 {     \
00332     SrcIterator send = s + shape[0]; \
00333     for(; s < send; ++s, ++d) \
00334     { \
00335         *d op *s; \
00336     } \
00337 } \
00338  \
00339 template <class SrcIterator, class Shape, class DestIterator, int N> \
00340 void \
00341 name##MultiArrayData(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>) \
00342 { \
00343     SrcIterator send = s + shape[N]; \
00344     for(; s < send; ++s, ++d) \
00345     { \
00346         name##MultiArrayData(s.begin(), shape, d.begin(), MetaInt<N-1>()); \
00347     } \
00348 } \
00349 \
00350 template <class DestIterator, class Shape, class T> \
00351 inline void \
00352 name##ScalarMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<0>) \
00353 {     \
00354     DestIterator dend = d + shape[0]; \
00355     for(; d < dend; ++d) \
00356     { \
00357         *d op init; \
00358     } \
00359 } \
00360  \
00361 template <class DestIterator, class Shape, class T, int N> \
00362 void \
00363 name##ScalarMultiArrayData(DestIterator d, Shape const & shape, T const & init, MetaInt<N>) \
00364 {     \
00365     DestIterator dend = d + shape[N]; \
00366     for(; d < dend; ++d) \
00367     { \
00368         name##ScalarMultiArrayData(d.begin(), shape, init, MetaInt<N-1>()); \
00369     } \
00370 }
00371 
00372 VIGRA_COPY_MULTI_ARRAY_DATA(copy, =)
00373 VIGRA_COPY_MULTI_ARRAY_DATA(copyAdd, +=)
00374 VIGRA_COPY_MULTI_ARRAY_DATA(copySub, -=)
00375 VIGRA_COPY_MULTI_ARRAY_DATA(copyMul, *=)
00376 VIGRA_COPY_MULTI_ARRAY_DATA(copyDiv, /=)
00377 
00378 #undef VIGRA_COPY_MULTI_ARRAY_DATA
00379 
00380 template <class SrcIterator, class Shape, class T, class ALLOC>
00381 inline void
00382 uninitializedCopyMultiArrayData(SrcIterator s, Shape const & shape, T * & d, ALLOC & a, MetaInt<0>)
00383 {
00384     SrcIterator send = s + shape[0];
00385     for(; s < send; ++s, ++d)
00386     {
00387         a.construct(d, static_cast<T const &>(*s));
00388     }
00389 }
00390 
00391 template <class SrcIterator, class Shape, class T, class ALLOC, int N>
00392 void
00393 uninitializedCopyMultiArrayData(SrcIterator s, Shape const & shape, T * & d, ALLOC & a, MetaInt<N>)
00394 {
00395     SrcIterator send = s + shape[N];
00396     for(; s < send; ++s)
00397     {
00398         uninitializedCopyMultiArrayData(s.begin(), shape, d, a, MetaInt<N-1>());
00399     }
00400 }
00401 
00402 template <class SrcIterator, class Shape, class T>
00403 inline void
00404 normMaxOfMultiArray(SrcIterator s, Shape const & shape, T & result, MetaInt<0>)
00405 {
00406     SrcIterator send = s + shape[0];
00407     for(; s < send; ++s)
00408     {
00409         T v = norm(*s);
00410         if(result < v)
00411             result = v;
00412     }
00413 }
00414 
00415 template <class SrcIterator, class Shape, class T, int N>
00416 void
00417 normMaxOfMultiArray(SrcIterator s, Shape const & shape, T & result, MetaInt<N>)
00418 {
00419     SrcIterator send = s + shape[N];
00420     for(; s < send; ++s)
00421     {
00422         normMaxOfMultiArray(s.begin(), shape, result, MetaInt<N-1>());
00423     }
00424 }
00425 
00426 template <class T>
00427 struct MultiArrayL1Functor
00428 {
00429     template <class U>
00430     T operator()(U t) const
00431     { return norm(t); }
00432 };
00433 
00434 template <class T>
00435 struct MultiArrayL2Functor
00436 {
00437     template <class U>
00438     T operator()(U t) const
00439     { return squaredNorm(t); }
00440 };
00441 
00442 template <class T>
00443 struct MultiArrayScaledL2Functor
00444 {
00445     T scale;
00446 
00447     MultiArrayScaledL2Functor(T s)
00448     : scale(s)
00449     {}
00450 
00451     template <class U>
00452     T operator()(U t) const
00453     { return squaredNorm(T(t) / scale); }
00454 };
00455 
00456 template <class SrcIterator, class Shape, class Functor, class T>
00457 inline void
00458 sumOverMultiArray(SrcIterator s, Shape const & shape, Functor f, T & result, MetaInt<0>)
00459 {
00460     SrcIterator send = s + shape[0];
00461     for(; s < send; ++s)
00462     {
00463         result += f(*s);
00464     }
00465 }
00466 
00467 template <class SrcIterator, class Shape, class Functor, class T, int N>
00468 void
00469 sumOverMultiArray(SrcIterator s, Shape const & shape, Functor f, T & result, MetaInt<N>)
00470 {
00471     SrcIterator send = s + shape[N];
00472     for(; s < send; ++s)
00473     {
00474         sumOverMultiArray(s.begin(), shape, f, result, MetaInt<N-1>());
00475     }
00476 }
00477 
00478 template <class SrcIterator, class Shape, class DestIterator>
00479 inline bool
00480 equalityOfMultiArrays(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>)
00481 {
00482     SrcIterator send = s + shape[0];
00483     for(; s < send; ++s, ++d)
00484     {
00485         if(!(*s == *d))
00486             return false;
00487     }
00488     return true;
00489 }
00490 
00491 template <class SrcIterator, class Shape, class DestIterator, int N>
00492 bool
00493 equalityOfMultiArrays(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>)
00494 {
00495     SrcIterator send = s + shape[N];
00496     for(; s < send; ++s, ++d)
00497     {
00498         if(!equalityOfMultiArrays(s.begin(), shape, d.begin(), MetaInt<N-1>()))
00499             return false;
00500     }
00501     return true;
00502 }
00503 
00504 
00505 template <class SrcIterator, class Shape, class DestIterator>
00506 inline void
00507 swapDataImpl(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<0>)
00508 {
00509     SrcIterator send = s + shape[0];
00510     for(; s < send; ++s, ++d)
00511         std::swap(*s, *d);
00512 }
00513 
00514 template <class SrcIterator, class Shape, class DestIterator, int N>
00515 void
00516 swapDataImpl(SrcIterator s, Shape const & shape, DestIterator d, MetaInt<N>)
00517 {
00518     SrcIterator send = s + shape[N];
00519     for(; s < send; ++s, ++d)
00520         swapDataImpl(s.begin(), shape, d.begin(), MetaInt<N-1>());
00521 }
00522 
00523 } // namespace detail
00524 
00525 /********************************************************/
00526 /*                                                      */
00527 /*                     MultiArrayView                   */
00528 /*                                                      */
00529 /********************************************************/
00530 
00531 // forward declaration
00532 template <unsigned int N, class T, class C = UnstridedArrayTag>
00533 class MultiArrayView;
00534 template <unsigned int N, class T, class A = std::allocator<T> >
00535 class MultiArray;
00536 
00537 /********************************************************/
00538 /*                                                      */
00539 /*                       NormTraits                     */
00540 /*                                                      */
00541 /********************************************************/
00542 
00543 template <unsigned int N, class T, class C>
00544 struct NormTraits<MultiArrayView<N, T, C> >
00545 {
00546     typedef MultiArrayView<N, T, C>                                      Type;
00547     typedef typename NormTraits<T>::SquaredNormType SquaredNormType;
00548     typedef typename SquareRootTraits<SquaredNormType>::SquareRootResult NormType;
00549 };
00550 
00551 template <unsigned int N, class T, class A>
00552 struct NormTraits<MultiArray<N, T, A> >
00553 : public NormTraits<MultiArrayView<N, T, UnstridedArrayTag> >
00554 {
00555     typedef NormTraits<MultiArrayView<N, T, UnstridedArrayTag> > BaseType;
00556     typedef MultiArray<N, T, A>                                  Type;
00557     typedef typename BaseType::SquaredNormType                   SquaredNormType;
00558     typedef typename BaseType::NormType                          NormType;
00559 };
00560 
00561 /** \brief Base class for, and view to, \ref vigra::MultiArray.
00562 
00563 This class implements the interface of both MultiArray and
00564 MultiArrayView.  By default, MultiArrayViews are tagged as
00565 unstrided. If necessary, strided arrays are constructed automatically
00566 by calls to a variant of the bind...() function.
00567 
00568 If you want to apply an algorithm requiring an image to a
00569 <tt>MultiArrayView</tt> of appropriate (2-dimensional) shape, you can
00570 create a \ref vigra::BasicImageView that acts as a wrapper with the
00571 necessary interface -- see \ref MultiArrayToImage.
00572 
00573 The template parameter are as follows
00574 \code
00575     N: the array dimension
00576 
00577     T: the type of the array elements
00578 
00579     C: a tag determining whether the array's inner dimension is strided
00580        or not. An array is unstrided if the array elements occupy consecutive
00581        memory location, strided if there is an offset in between (e.g.
00582        when a view is created that skips every other array element).
00583        The compiler can generate faster code for unstrided arrays.
00584        Possible values: UnstridedArrayTag (default), StridedArrayTag
00585 \endcode
00586 
00587 <b>\#include</b>
00588 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
00589 
00590 Namespace: vigra
00591 */
00592 template <unsigned int N, class T, class C>
00593 class MultiArrayView
00594 {
00595 public:
00596 
00597         /** the array's actual dimensionality.
00598             This ensures that MultiArrayView can also be used for
00599             scalars (that is, when <tt>N == 0</tt>). Calculated as:<br>
00600             \code
00601             actual_dimension = (N==0) ? 1 : N
00602             \endcode
00603          */
00604     enum ActualDimension { actual_dimension = (N==0) ? 1 : N };
00605 
00606         /** the array's value type
00607          */
00608     typedef T value_type;
00609 
00610         /** reference type (result of operator[])
00611          */
00612     typedef value_type &reference;
00613 
00614         /** const reference type (result of operator[] const)
00615          */
00616     typedef const value_type &const_reference;
00617 
00618         /** pointer type
00619          */
00620     typedef value_type *pointer;
00621 
00622         /** const pointer type
00623          */
00624     typedef const value_type *const_pointer;
00625 
00626         /** difference type (used for multi-dimensional offsets and indices)
00627          */
00628     typedef typename MultiArrayShape<actual_dimension>::type difference_type;
00629 
00630         /** size type
00631          */
00632     typedef difference_type size_type;
00633 
00634         /** difference and index type for a single dimension
00635          */
00636     typedef MultiArrayIndex difference_type_1;
00637 
00638         /** traverser (MultiIterator) type
00639          */
00640     typedef typename vigra::detail::MultiIteratorChooser <
00641         C>::template Traverser <actual_dimension, T, T &, T *>::type traverser;
00642 
00643         /** const traverser (MultiIterator) type
00644          */
00645     typedef typename vigra::detail::MultiIteratorChooser <
00646         C>::template Traverser <actual_dimension, T, T const &, T const *>::type const_traverser;
00647 
00648         /** the view type associated with this array.
00649          */
00650     typedef MultiArrayView <N, T, C> view_type;
00651 
00652         /** the matrix type associated with this array.
00653          */
00654     typedef MultiArray <N, T> matrix_type;
00655 
00656 protected:
00657 
00658     typedef typename difference_type::value_type diff_zero_t;
00659 
00660         /** the shape of the image pointed to is stored here.
00661         */
00662     difference_type m_shape;
00663 
00664         /** the strides (offset of a sample to the next) for every dimension
00665             are stored here.
00666         */
00667     difference_type m_stride;
00668 
00669         /** pointer to the image.
00670          */
00671     pointer m_ptr;
00672 
00673     template <class U, class CN>
00674     void copyImpl(const MultiArrayView <N, U, CN>& rhs);
00675 
00676     template <class U, class CN>
00677     void swapDataImpl(MultiArrayView <N, U, CN> rhs);
00678 
00679     template <class U, class CN>
00680     bool arraysOverlap(const MultiArrayView <N, U, CN>& rhs) const;
00681 
00682 public:
00683 
00684         /** default constructor: create an empty image of size 0.
00685          */
00686     MultiArrayView ()
00687         : m_shape (diff_zero_t(0)), m_stride (diff_zero_t(0)), m_ptr (0)
00688     {}
00689 
00690         /** construct from shape and pointer
00691          */
00692     MultiArrayView (const difference_type &shape, pointer ptr)
00693     : m_shape (shape),
00694       m_stride (detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape)),
00695       m_ptr (ptr)
00696     {}
00697 
00698         /** construct from shape, strides (offset of a sample to the next)
00699             for every dimension) and pointer
00700          */
00701     MultiArrayView (const difference_type &shape,
00702                     const difference_type &stride,
00703                     pointer ptr)
00704     : m_shape (shape),
00705       m_stride (stride),
00706       m_ptr (ptr)
00707     {}
00708 
00709 
00710         /** Assignment. There are 3 cases:
00711 
00712             <ul>
00713             <li> When this <tt>MultiArrayView</tt> does not point to valid data
00714                  (e.g. after default construction), it becomes a copy of \a rhs.
00715             <li> When the shapes of the two arrays match, the array contents are copied.
00716             <li> Otherwise, a <tt>PreconditionViolation</tt> exception is thrown.
00717             </ul>
00718          */
00719     MultiArrayView & operator=(MultiArrayView const & rhs);
00720 
00721         /** Assignment of a differently typed MultiArrayView. Fails with
00722             <tt>PreconditionViolation</tt> exception when the shapes do not match.
00723          */
00724     template<class U, class C1>
00725     MultiArrayView & operator=(MultiArrayView<N, U, C1> const & rhs)
00726     {
00727         vigra_precondition(this->shape() == rhs.shape(),
00728             "MultiArrayView::operator=() size mismatch.");
00729         this->copyImpl(rhs);
00730         return *this;
00731     }
00732 
00733         /** Add-assignment of a compatible MultiArrayView. Fails with
00734             <tt>PreconditionViolation</tt> exception when the shapes do not match.
00735          */
00736     template<class U, class C1>
00737     MultiArrayView & operator+=(MultiArrayView<N, U, C1> const & rhs);
00738 
00739         /** Subtract-assignment of a compatible MultiArrayView. Fails with
00740             <tt>PreconditionViolation</tt> exception when the shapes do not match.
00741          */
00742     template<class U, class C1>
00743     MultiArrayView & operator-=(MultiArrayView<N, U, C1> const & rhs);
00744 
00745         /** Multiply-assignment of a compatible MultiArrayView. Fails with
00746             <tt>PreconditionViolation</tt> exception when the shapes do not match.
00747          */
00748     template<class U, class C1>
00749     MultiArrayView & operator*=(MultiArrayView<N, U, C1> const & rhs);
00750 
00751         /** Divide-assignment of a compatible MultiArrayView. Fails with
00752             <tt>PreconditionViolation</tt> exception when the shapes do not match.
00753          */
00754     template<class U, class C1>
00755     MultiArrayView & operator/=(MultiArrayView<N, U, C1> const & rhs);
00756 
00757         /** Add-assignment of a scalar.
00758          */
00759     MultiArrayView & operator+=(T const & rhs)
00760     {
00761         detail::copyAddScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
00762         return *this;
00763     }
00764 
00765         /** Subtract-assignment of a scalar.
00766          */
00767     MultiArrayView & operator-=(T const & rhs)
00768     {
00769         detail::copySubScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
00770         return *this;
00771     }
00772 
00773         /** Multiply-assignment of a scalar.
00774          */
00775     MultiArrayView & operator*=(T const & rhs)
00776     {
00777         detail::copyMulScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
00778         return *this;
00779     }
00780 
00781         /** Divide-assignment of a scalar.
00782          */
00783     MultiArrayView & operator/=(T const & rhs)
00784     {
00785         detail::copyDivScalarMultiArrayData(traverser_begin(), shape(), rhs, MetaInt<actual_dimension-1>());
00786         return *this;
00787     }
00788 
00789         /** array access.
00790          */
00791     reference operator[] (const difference_type &d)
00792     {
00793         return m_ptr [dot (d, m_stride)];
00794     }
00795 
00796         /** array access.
00797          */
00798     const_reference operator[] (const difference_type &d) const
00799     {
00800         return m_ptr [dot (d, m_stride)];
00801     }
00802 
00803         /** array access in scan-order sense.
00804             Mostly useful to support standard indexing for 1-dimensional multi-arrays,
00805             but works for any N. Use scanOrderIndexToCoordinate() and
00806             coordinateToScanOrderIndex() for conversion between indices and coordinates.
00807          */
00808     reference operator[](difference_type_1 d)
00809     {
00810         return m_ptr [detail::ScanOrderToOffset<actual_dimension>::exec(d, m_shape, m_stride)];
00811     }
00812 
00813         /** array access in scan-order sense.
00814             Mostly useful to support standard indexing for 1-dimensional multi-arrays,
00815             but works for any N. Use scanOrderIndexToCoordinate() and
00816             coordinateToScanOrderIndex() for conversion between indices and coordinates.
00817          */
00818     const_reference operator[](difference_type_1 d) const
00819     {
00820         return m_ptr [detail::ScanOrderToOffset<actual_dimension>::exec(d, m_shape, m_stride)];
00821     }
00822 
00823         /** convert scan-order index to coordinate.
00824          */
00825     difference_type scanOrderIndexToCoordinate(difference_type_1 d) const
00826     {
00827         difference_type result;
00828         detail::ScanOrderToCoordinate<actual_dimension>::exec(d, m_shape, result);
00829         return result;
00830     }
00831 
00832         /** convert coordinate to scan-order index.
00833          */
00834     difference_type_1 coordinateToScanOrderIndex(const difference_type &d) const
00835     {
00836         return detail::CoordinateToScanOrder<actual_dimension>::exec(m_shape, d);
00837     }
00838 
00839         /** 1D array access. Use only if N == 1.
00840          */
00841     reference operator() (difference_type_1 x)
00842     {
00843         return m_ptr [detail::CoordinatesToOffest<C>::exec(m_stride, x)];
00844     }
00845 
00846         /** 2D array access. Use only if N == 2.
00847          */
00848     reference operator() (difference_type_1 x, difference_type_1 y)
00849     {
00850         return m_ptr [detail::CoordinatesToOffest<C>::exec(m_stride, x, y)];
00851     }
00852 
00853         /** 3D array access. Use only if N == 3.
00854          */
00855     reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z)
00856     {
00857         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
00858     }
00859 
00860         /** 4D array access. Use only if N == 4.
00861          */
00862     reference operator() (difference_type_1 x, difference_type_1 y,
00863                           difference_type_1 z, difference_type_1 u)
00864     {
00865         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
00866     }
00867 
00868         /** 5D array access. Use only if N == 5.
00869          */
00870     reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
00871                           difference_type_1 u, difference_type_1 v)
00872     {
00873         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
00874     }
00875 
00876         /** 1D const array access. Use only if N == 1.
00877          */
00878     const_reference operator() (difference_type_1 x) const
00879     {
00880         return m_ptr [detail::CoordinatesToOffest<C>::exec(m_stride, x)];
00881     }
00882 
00883         /** 2D const array access. Use only if N == 2.
00884          */
00885     const_reference operator() (difference_type_1 x, difference_type_1 y) const
00886     {
00887         return m_ptr [detail::CoordinatesToOffest<C>::exec(m_stride, x, y)];
00888     }
00889 
00890         /** 3D const array access. Use only if N == 3.
00891          */
00892     const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z) const
00893     {
00894         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z];
00895     }
00896 
00897         /** 4D const array access. Use only if N == 4.
00898          */
00899     const_reference operator() (difference_type_1 x, difference_type_1 y,
00900                                 difference_type_1 z, difference_type_1 u) const
00901     {
00902         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u];
00903     }
00904 
00905         /** 5D const array access. Use only if N == 5.
00906          */
00907     const_reference operator() (difference_type_1 x, difference_type_1 y, difference_type_1 z,
00908                                 difference_type_1 u, difference_type_1 v) const
00909     {
00910         return m_ptr [m_stride[0]*x + m_stride[1]*y + m_stride[2]*z + m_stride[3]*u + m_stride[4]*v];
00911     }
00912 
00913         /** Init with a constant.
00914          */
00915     template <class U>
00916     MultiArrayView & init(const U & init)
00917     {
00918         detail::copyScalarMultiArrayData(traverser_begin(), shape(), init, MetaInt<actual_dimension-1>());
00919         return *this;
00920     }
00921 
00922 
00923         /** Copy the data of the right-hand array (array shapes must match).
00924          */
00925     void copy(const MultiArrayView & rhs)
00926     {
00927         if(this == &rhs)
00928             return;
00929         this->copyImpl(rhs);
00930     }
00931 
00932         /** Copy the data of the right-hand array (array shapes must match).
00933          */
00934     template <class U, class CN>
00935     void copy(const MultiArrayView <N, U, CN>& rhs)
00936     {
00937         this->copyImpl(rhs);
00938     }
00939 
00940         /** swap the data between two MultiArrayView objects.
00941 
00942             The shapes of the two array must match.
00943         */
00944     void swapData(MultiArrayView rhs)
00945     {
00946         if(this != &rhs)
00947             swapDataImpl(rhs);
00948     }
00949 
00950         /** swap the data between two MultiArrayView objects.
00951 
00952             The shapes of the two array must match.
00953         */
00954     template <class T2, class C2>
00955     void swapData(MultiArrayView <N, T2, C2> rhs)
00956     {
00957         swapDataImpl(rhs);
00958     }
00959 
00960         /** bind the M outmost dimensions to certain indices.
00961             this reduces the dimensionality of the image to
00962             max { 1, N-M }.
00963 
00964             <b>Usage:</b>
00965             \code
00966             // create a 3D array of size 40x30x20
00967             typedef MultiArray<3, double>::difference_type Shape;
00968             MultiArray<3, double> array3(Shape(40, 30, 20));
00969 
00970             // get a 1D array by fixing index 1 to 12, and index 2 to 10
00971             MultiArrayView <1, double> array1 = array3.bindOuter(TinyVector<int, 2>(12, 10));
00972             \endcode
00973         */
00974     template <unsigned int M>
00975     MultiArrayView <N-M, T, C> bindOuter (const TinyVector <MultiArrayIndex, M> &d) const;
00976 
00977         /** bind the M innermost dimensions to certain indices.
00978             this reduces the dimensionality of the image to
00979             max { 1, N-M }.
00980 
00981             <b>Usage:</b>
00982             \code
00983             // create a 3D array of size 40x30x20
00984             typedef MultiArray<3, double>::difference_type Shape;
00985             MultiArray<3, double> array3(Shape(40, 30, 20));
00986 
00987             // get a 1D array by fixing index 0 to 12, and index 1 to 10
00988             MultiArrayView <1, double, StridedArrayTag> array1 = array3.bindInner(TinyVector<int, 2>(12, 10));
00989             \endcode
00990         */
00991     template <unsigned int M>
00992     MultiArrayView <N-M, T, StridedArrayTag>
00993     bindInner (const TinyVector <MultiArrayIndex, M> &d) const;
00994 
00995         /** bind dimension M to index d.
00996             this reduces the dimensionality of the image to
00997             max { 1, N-1 }.
00998 
00999             <b>Usage:</b>
01000             \code
01001             // create a 3D array of size 40x30x20
01002             typedef MultiArray<3, double>::difference_type Shape;
01003             MultiArray<3, double> array3(Shape(40, 30, 20));
01004 
01005             // get a 2D array by fixing index 1 to 12
01006             MultiArrayView <2, double> array2 = array3.bind<1>(12);
01007 
01008             // get a 2D array by fixing index 0 to 23
01009             MultiArrayView <2, double, StridedArrayTag> array2a = array3.bind<0>(23);
01010             \endcode
01011          */
01012     template <unsigned int M>
01013     MultiArrayView <N-1, T, typename vigra::detail::MaybeStrided <M>::type >
01014     bind (difference_type_1 d) const;
01015 
01016         /** bind the outmost dimension to a certain index.
01017             this reduces the dimensionality of the image to
01018             max { 1, N-1 }.
01019 
01020             <b>Usage:</b>
01021             \code
01022             // create a 3D array of size 40x30x20
01023             typedef MultiArray<3, double>::difference_type Shape;
01024             MultiArray<3, double> array3(Shape(40, 30, 20));
01025 
01026             // get a 2D array by fixing the outermost index (i.e. index 2) to 12
01027             MultiArrayView <2, double> array2 = array3.bindOuter(12);
01028             \endcode
01029         */
01030     MultiArrayView <N-1, T, C> bindOuter (difference_type_1 d) const;
01031 
01032         /** bind the innermost dimension to a certain index.
01033             this reduces the dimensionality of the image to
01034             max { 1, N-1 }.
01035 
01036             <b>Usage:</b>
01037             \code
01038             // create a 3D array of size 40x30x20
01039             typedef MultiArray<3, double>::difference_type Shape;
01040             MultiArray<3, double> array3(Shape(40, 30, 20));
01041 
01042             // get a 2D array by fixing the innermost index (i.e. index 0) to 23
01043             MultiArrayView <2, double, StridedArrayTag> array2 = array3.bindInner(23);
01044             \endcode
01045         */
01046     MultiArrayView <N-1, T, StridedArrayTag> bindInner (difference_type_1 d) const;
01047 
01048         /** bind dimension m to index d.
01049             this reduces the dimensionality of the image to
01050             max { 1, N-1 }.
01051 
01052             <b>Usage:</b>
01053             \code
01054             // create a 3D array of size 40x30x20
01055             typedef MultiArray<3, double>::difference_type Shape;
01056             MultiArray<3, double> array3(Shape(40, 30, 20));
01057 
01058             // get a 2D array by fixing index 2 to 15
01059             MultiArrayView <2, double, StridedArrayTag> array2 = array3.bindAt(2, 15);
01060             \endcode
01061          */
01062     MultiArrayView <N-1, T, StridedArrayTag>
01063     bindAt (difference_type_1 m, difference_type_1 d) const;
01064 
01065         /** create a rectangular subarray that spans between the
01066             points p and q, where p is in the subarray, q not.
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 subarray set is smaller by one element at all sides
01075             MultiArrayView <3, double> subarray = array3.subarray(Shape(1,1,1), Shape(39, 29, 19));
01076             \endcode
01077         */
01078     MultiArrayView subarray (const difference_type &p,
01079                              const difference_type &q) const
01080     {
01081         const difference_type_1 offset = dot (m_stride, p);
01082         return MultiArrayView (q - p, m_stride, m_ptr + offset);
01083     }
01084 
01085         /** apply an additional striding to the image, thereby reducing
01086             the shape of the array.
01087             for example, multiplying the stride of dimension one by three
01088             turns an appropriately layed out (interleaved) rgb image into
01089             a single band image.
01090         */
01091     MultiArrayView <N, T, StridedArrayTag>
01092     stridearray (const difference_type &s) const
01093     {
01094         difference_type shape = m_shape;
01095         for (unsigned int i = 0; i < actual_dimension; ++i)
01096             shape [i] /= s [i];
01097         return MultiArrayView <N, T, StridedArrayTag>
01098             (shape, m_stride * s, m_ptr);
01099     }
01100 
01101         /** permute the dimensions of the array.
01102             The function exchanges the meaning of the dimensions without copying the data.
01103             In case of a 2-dimensional array, this is simply array transposition. In higher dimensions,
01104             there are more possibilities.
01105 
01106             <b>Usage:</b><br>
01107             \code
01108             typedef MultiArray<2, double>::difference_type Shape;
01109             MultiArray<2, double> array(10, 20);
01110 
01111             MultiArray<2, double, StridedArrayTag> transposed = array.permuteDimensions(Shape(1,0));
01112 
01113             for(int i=0; i<array.shape(0), ++i)
01114                 for(int j=0; j<array.shape(1); ++j)
01115                     assert(array(i, j) == transposed(j, i));
01116             \endcode
01117         */
01118     MultiArrayView <N, T, StridedArrayTag>
01119     permuteDimensions (const difference_type &s) const
01120     {
01121         difference_type shape, stride, check((typename difference_type::value_type)0);
01122         for (unsigned int i = 0; i < actual_dimension; ++i)
01123         {
01124             shape[i]  = m_shape[s[i]];
01125             stride[i] = m_stride[s[i]];
01126             ++check[s[i]];
01127         }
01128         vigra_precondition(check == difference_type(1),
01129            "MultiArrayView::permuteDimensions(): every dimension must occur exactly once.");
01130         return MultiArrayView <N, T, StridedArrayTag>(shape, stride, m_ptr);
01131     }
01132 
01133         /** transpose a 2-dimensional array. Use only if N==2.
01134 
01135             <b>Usage:</b><br>
01136             \code
01137             typedef MultiArray<2, double>::difference_type Shape;
01138             MultiArray<2, double> array(10, 20);
01139 
01140             MultiArray<2, double, StridedArrayTag> transposed = array.transpose();
01141 
01142             for(int i=0; i<array.shape(0), ++i)
01143                 for(int j=0; j<array.shape(1); ++j)
01144                     assert(array(i, j) == transposed(j, i));
01145             \endcode
01146         */
01147     MultiArrayView <2, T, StridedArrayTag>
01148     transpose () const
01149     {
01150         difference_type shape(m_shape[1], m_shape[0]),
01151                         stride(m_stride[1], m_stride[0]);
01152         return MultiArrayView <2, T, StridedArrayTag>(shape, stride, m_ptr);
01153     }
01154 
01155         /** number of the elements in the array.
01156          */
01157     difference_type_1 elementCount () const
01158     {
01159         difference_type_1 ret = m_shape[0];
01160         for(int i = 1; i < actual_dimension; ++i)
01161             ret *= m_shape[i];
01162         return ret;
01163     }
01164 
01165         /** number of the elements in the array.
01166             Same as <tt>elementCount()</tt>. Mostly useful to support the std::vector interface.
01167          */
01168     difference_type_1 size () const
01169     {
01170         return elementCount();
01171     }
01172 
01173         /** return the array's shape.
01174          */
01175     const difference_type & shape () const
01176     {
01177         return m_shape;
01178     }
01179 
01180         /** return the array's size at a certain dimension.
01181          */
01182     difference_type_1 size (difference_type_1 n) const
01183     {
01184         return m_shape [n];
01185     }
01186 
01187         /** return the array's shape at a certain dimension
01188             (same as <tt>size(n)</tt>).
01189          */
01190     difference_type_1 shape (difference_type_1 n) const
01191     {
01192         return m_shape [n];
01193     }
01194 
01195         /** return the array's stride for every dimension.
01196          */
01197     const difference_type & stride () const
01198     {
01199         return m_stride;
01200     }
01201 
01202         /** return the array's stride at a certain dimension.
01203          */
01204     difference_type_1 stride (int n) const
01205     {
01206         return m_stride [n];
01207     }
01208 
01209         /** check whether two arrays are elementwise equal.
01210          */
01211     template <class U, class C1>
01212     bool operator==(MultiArrayView<N, U, C1> const & rhs) const
01213     {
01214         if(this->shape() != rhs.shape())
01215             return false;
01216         return detail::equalityOfMultiArrays(traverser_begin(), shape(), rhs.traverser_begin(), MetaInt<actual_dimension-1>());
01217     }
01218 
01219         /** check whether two arrays are not elementwise equal.
01220             Also true when the two arrays have different shapes.
01221          */
01222     template <class U, class C1>
01223     bool operator!=(MultiArrayView<N, U, C1> const & rhs) const
01224     {
01225         return !operator==(rhs);
01226     }
01227 
01228         /** check whether the given point is in the array range.
01229          */
01230     bool isInside (difference_type const & p) const
01231     {
01232         for(int d=0; d<actual_dimension; ++d)
01233             if(p[d] < 0 || p[d] >= shape(d))
01234                 return false;
01235         return true;
01236     }
01237 
01238         /** Compute the squared Euclidean norm of the array (sum of squares of the array elements).
01239          */
01240     typename NormTraits<MultiArrayView>::SquaredNormType squaredNorm() const
01241     {
01242         typedef typename NormTraits<MultiArrayView>::SquaredNormType SquaredNormType;
01243         SquaredNormType res = NumericTraits<SquaredNormType>::zero();
01244         detail::sumOverMultiArray(traverser_begin(), shape(), detail::MultiArrayL2Functor<SquaredNormType>(),
01245                                   res, MetaInt<actual_dimension-1>());
01246         return res;
01247     }
01248 
01249         /** Compute various norms of the array.
01250             The norm is determined by parameter \a type:
01251 
01252             <ul>
01253             <li> type == 0: maximum norm (L-infinity): maximum of absolute values of the array elements
01254             <li> type == 1: Manhattan norm (L1): sum of absolute values of the array elements
01255             <li> type == 2: Euclidean norm (L2): square root of <tt>squaredNorm()</tt> when \a useSquaredNorm is <tt>true</tt>,<br>
01256                  or direct algorithm that avoids underflow/overflow otherwise.
01257             </ul>
01258 
01259             Parameter \a useSquaredNorm has no effect when \a type != 2. Defaults: compute L2 norm as square root of
01260             <tt>squaredNorm()</tt>.
01261          */
01262     typename NormTraits<MultiArrayView>::NormType norm(int type = 2, bool useSquaredNorm = true) const;
01263 
01264         /** return the pointer to the image data
01265          */
01266     pointer data () const
01267     {
01268         return m_ptr;
01269     }
01270 
01271         /** returns the N-dimensional MultiIterator pointing
01272             to the first element in every dimension.
01273         */
01274     traverser traverser_begin ()
01275     {
01276         traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
01277         return ret;
01278     }
01279 
01280         /** returns the N-dimensional MultiIterator pointing
01281             to the const first element in every dimension.
01282         */
01283     const_traverser traverser_begin () const
01284     {
01285         const_traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
01286         return ret;
01287     }
01288 
01289         /** returns the N-dimensional MultiIterator pointing
01290             beyond the last element in dimension N, and to the
01291             first element in every other dimension.
01292         */
01293     traverser traverser_end ()
01294     {
01295         traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
01296         ret += m_shape [actual_dimension-1];
01297         return ret;
01298     }
01299 
01300         /** returns the N-dimensional const MultiIterator pointing
01301             beyond the last element in dimension N, and to the
01302             first element in every other dimension.
01303         */
01304     const_traverser traverser_end () const
01305     {
01306         const_traverser ret (m_ptr, m_stride.begin (), m_shape.begin ());
01307         ret += m_shape [actual_dimension-1];
01308         return ret;
01309     }
01310 
01311     view_type view ()
01312     {
01313         return *this;
01314     }
01315 };
01316 
01317 template <unsigned int N, class T, class C>
01318 MultiArrayView<N, T, C> &
01319 MultiArrayView <N, T, C>::operator=(MultiArrayView<N, T, C> const & rhs)
01320 {
01321     if(this == &rhs)
01322         return *this;
01323     vigra_precondition(this->shape() == rhs.shape() || m_ptr == 0,
01324         "MultiArrayView::operator=(MultiArrayView const &) size mismatch - use MultiArrayView::reset().");
01325     if(m_ptr == 0)
01326     {
01327         m_shape  = rhs.m_shape;
01328         m_stride = rhs.m_stride;
01329         m_ptr    = rhs.m_ptr;
01330     }
01331     else
01332         this->copyImpl(rhs);
01333     return *this;
01334 }
01335 
01336 template <unsigned int N, class T, class C>
01337 template <class U, class CN>
01338 bool
01339 MultiArrayView <N, T, C>::arraysOverlap(const MultiArrayView <N, U, CN>& rhs) const
01340 {
01341     vigra_precondition (shape () == rhs.shape (),
01342         "MultiArrayView::arraysOverlap(): shape mismatch.");
01343     const_pointer first_element = this->m_ptr,
01344                   last_element = first_element + dot(this->m_shape - difference_type(1), this->m_stride);
01345     typename MultiArrayView <N, U, CN>::const_pointer
01346            rhs_first_element = rhs.data(),
01347            rhs_last_element = rhs_first_element + dot(rhs.shape() - difference_type(1), rhs.stride());
01348     return !(last_element < rhs_first_element || rhs_last_element < first_element);
01349 }
01350 
01351 template <unsigned int N, class T, class C>
01352 template <class U, class CN>
01353 void
01354 MultiArrayView <N, T, C>::copyImpl(const MultiArrayView <N, U, CN>& rhs)
01355 {
01356     if(!arraysOverlap(rhs))
01357     {
01358         // no overlap -- can copy directly
01359         detail::copyMultiArrayData(rhs.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>());
01360     }
01361     else
01362     {
01363         // overlap: we got different views to the same data -- copy to intermediate memory in order to avoid
01364         // overwriting elements that are still needed on the rhs.
01365         MultiArray<N, T> tmp(rhs);
01366         detail::copyMultiArrayData(tmp.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>());
01367     }
01368 }
01369 
01370 #define VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT(name, op) \
01371 template <unsigned int N, class T, class C> \
01372 template<class U, class C1> \
01373 MultiArrayView<N, T, C> &  \
01374 MultiArrayView <N, T, C>::operator op(MultiArrayView<N, U, C1> const & rhs) \
01375 { \
01376     vigra_precondition(this->shape() == rhs.shape(), "MultiArrayView::operator" #op "() size mismatch."); \
01377     if(!arraysOverlap(rhs)) \
01378     { \
01379         detail::name##MultiArrayData(rhs.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>()); \
01380     } \
01381     else \
01382     { \
01383         MultiArray<N, T> tmp(rhs); \
01384         detail::name##MultiArrayData(tmp.traverser_begin(), shape(), traverser_begin(), MetaInt<actual_dimension-1>()); \
01385     } \
01386     return *this; \
01387 }
01388 
01389 VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT(copyAdd, +=)
01390 VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT(copySub, -=)
01391 VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT(copyMul, *=)
01392 VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT(copyDiv, /=)
01393 
01394 #undef VIGRA_MULTI_ARRY_COMPUTED_ASSIGNMENT
01395 
01396 template <unsigned int N, class T, class C>
01397 template <class U, class CN>
01398 void
01399 MultiArrayView <N, T, C>::swapDataImpl(MultiArrayView <N, U, CN> rhs)
01400 {
01401     vigra_precondition (shape () == rhs.shape (),
01402         "MultiArrayView::swapData(): shape mismatch.");
01403 
01404     // check for overlap of this and rhs
01405     const_pointer first_element = this->m_ptr,
01406                   last_element = first_element + dot(this->m_shape - difference_type(1), this->m_stride);
01407     typename MultiArrayView <N, U, CN>::const_pointer
01408            rhs_first_element = rhs.data(),
01409            rhs_last_element = rhs_first_element + dot(rhs.shape() - difference_type(1), rhs.stride());
01410     if(last_element < rhs_first_element || rhs_last_element < first_element)
01411     {
01412         // no overlap -- can swap directly
01413         detail::swapDataImpl(traverser_begin(), shape(), rhs.traverser_begin(), MetaInt<actual_dimension-1>());
01414     }
01415     else
01416     {
01417         // overlap: we got different views to the same data -- copy to intermediate memory in order to avoid
01418         // overwriting elements that are still needed.
01419         MultiArray<N, T> tmp(*this);
01420         copy(rhs);
01421         rhs.copy(tmp);
01422     }
01423 }
01424 
01425 template <unsigned int N, class T, class C>
01426 template <unsigned int M>
01427 MultiArrayView <N-M, T, C>
01428 MultiArrayView <N, T, C>::bindOuter (const TinyVector <MultiArrayIndex, M> &d) const
01429 {
01430     TinyVector <MultiArrayIndex, M> stride;
01431     stride.init (m_stride.begin () + N-M, m_stride.end ());
01432     pointer ptr = m_ptr + dot (d, stride);
01433     static const int NNew = (N-M == 0) ? 1 : N-M;
01434     TinyVector <MultiArrayIndex, NNew> inner_shape, inner_stride;
01435     if (N-M == 0)
01436     {
01437         inner_shape [0] = 1;
01438         inner_stride [0] = 0;
01439     }
01440     else
01441     {
01442         inner_shape.init (m_shape.begin (), m_shape.end () - M);
01443         inner_stride.init (m_stride.begin (), m_stride.end () - M);
01444     }
01445     return MultiArrayView <N-M, T, C> (inner_shape, inner_stride, ptr);
01446 }
01447 
01448 template <unsigned int N, class T, class C>
01449 template <unsigned int M>
01450 MultiArrayView <N - M, T, StridedArrayTag>
01451 MultiArrayView <N, T, C>::bindInner (const TinyVector <MultiArrayIndex, M> &d) const
01452 {
01453     TinyVector <MultiArrayIndex, M> stride;
01454     stride.init (m_stride.begin (), m_stride.end () - N + M);
01455     pointer ptr = m_ptr + dot (d, stride);
01456     static const int NNew = (N-M == 0) ? 1 : N-M;
01457     TinyVector <MultiArrayIndex, NNew> outer_shape, outer_stride;
01458     if (N-M == 0)
01459     {
01460         outer_shape [0] = 1;
01461         outer_stride [0] = 0;
01462     }
01463     else
01464     {
01465         outer_shape.init (m_shape.begin () + M, m_shape.end ());
01466         outer_stride.init (m_stride.begin () + M, m_stride.end ());
01467     }
01468     return MultiArrayView <N-M, T, StridedArrayTag>
01469         (outer_shape, outer_stride, ptr);
01470 }
01471 
01472 template <unsigned int N, class T, class C>
01473 template <unsigned int M>
01474 MultiArrayView <N-1, T, typename detail::MaybeStrided <M>::type >
01475 MultiArrayView <N, T, C>::bind (difference_type_1 d) const
01476 {
01477     static const int NNew = (N-1 == 0) ? 1 : N-1;
01478     TinyVector <MultiArrayIndex, NNew> shape, stride;
01479     // the remaining dimensions are 0..n-1,n+1..N-1
01480     if (N-1 == 0)
01481     {
01482         shape[0] = 1;
01483         stride[0] = 0;
01484     }
01485     else
01486     {
01487         std::copy (m_shape.begin (), m_shape.begin () + M, shape.begin ());
01488         std::copy (m_shape.begin () + M+1, m_shape.end (),
01489                    shape.begin () + M);
01490         std::copy (m_stride.begin (), m_stride.begin () + M, stride.begin ());
01491         std::copy (m_stride.begin () + M+1, m_stride.end (),
01492                    stride.begin () + M);
01493     }
01494     return MultiArrayView <N-1, T, typename detail::MaybeStrided <M>::type>
01495         (shape, stride, m_ptr + d * m_stride[M]);
01496 }
01497 
01498 template <unsigned int N, class T, class C>
01499 MultiArrayView <N - 1, T, C>
01500 MultiArrayView <N, T, C>::bindOuter (difference_type_1 d) const
01501 {
01502     static const int NNew = (N-1 == 0) ? 1 : N-1;
01503     TinyVector <MultiArrayIndex, NNew> inner_shape, inner_stride;
01504     if (N-1 == 0)
01505     {
01506         inner_shape [0] = 1;
01507         inner_stride [0] = 0;
01508     }
01509     else
01510     {
01511         inner_shape.init (m_shape.begin (), m_shape.end () - 1);
01512         inner_stride.init (m_stride.begin (), m_stride.end () - 1);
01513     }
01514     return MultiArrayView <N-1, T, C> (inner_shape, inner_stride,
01515                                        m_ptr + d * m_stride [N-1]);
01516 }
01517 
01518 template <unsigned int N, class T, class C>
01519 MultiArrayView <N - 1, T, StridedArrayTag>
01520 MultiArrayView <N, T, C>::bindInner (difference_type_1 d) const
01521 {
01522     static const int NNew = (N-1 == 0) ? 1 : N-1;
01523     TinyVector <MultiArrayIndex, NNew> outer_shape, outer_stride;
01524     if (N-1 == 0)
01525     {
01526         outer_shape [0] = 1;
01527         outer_stride [0] = 0;
01528     }
01529     else
01530     {
01531         outer_shape.init (m_shape.begin () + 1, m_shape.end ());
01532         outer_stride.init (m_stride.begin () + 1, m_stride.end ());
01533     }
01534     return MultiArrayView <N-1, T, StridedArrayTag>
01535         (outer_shape, outer_stride, m_ptr + d * m_stride [0]);
01536 }
01537 
01538 template <unsigned int N, class T, class C>
01539 MultiArrayView <N - 1, T, StridedArrayTag>
01540 MultiArrayView <N, T, C>::bindAt (difference_type_1 n, difference_type_1 d) const
01541 {
01542     vigra_precondition (
01543         n < static_cast <int> (N),
01544         "MultiArrayView <N, T, C>::bindAt(): dimension out of range.");
01545     static const int NNew = (N-1 == 0) ? 1 : N-1;
01546     TinyVector <MultiArrayIndex, NNew> shape, stride;
01547     // the remaining dimensions are 0..n-1,n+1..N-1
01548     if (N-1 == 0)
01549     {
01550         shape [0] = 1;
01551         stride [0] = 0;
01552     }
01553     else
01554     {
01555         std::copy (m_shape.begin (), m_shape.begin () + n, shape.begin ());
01556         std::copy (m_shape.begin () + n+1, m_shape.end (),
01557                    shape.begin () + n);
01558         std::copy (m_stride.begin (), m_stride.begin () + n, stride.begin ());
01559         std::copy (m_stride.begin () + n+1, m_stride.end (),
01560                    stride.begin () + n);
01561     }
01562     return MultiArrayView <N-1, T, StridedArrayTag>
01563         (shape, stride, m_ptr + d * m_stride[n]);
01564 }
01565 
01566 template <unsigned int N, class T, class C>
01567 typename NormTraits<MultiArrayView <N, T, C> >::NormType
01568 MultiArrayView <N, T, C>::norm(int type, bool useSquaredNorm) const
01569 {
01570     typedef typename NormTraits<MultiArrayView>::NormType NormType;
01571 
01572     switch(type)
01573     {
01574       case 0:
01575       {
01576         NormType res = NumericTraits<NormType>::zero();
01577         detail::normMaxOfMultiArray(traverser_begin(), shape(), res, MetaInt<actual_dimension-1>());
01578         return res;
01579       }
01580       case 1:
01581       {
01582         NormType res = NumericTraits<NormType>::zero();
01583         detail::sumOverMultiArray(traverser_begin(), shape(), detail::MultiArrayL1Functor<NormType>(),
01584                                 res, MetaInt<actual_dimension-1>());
01585         return res;
01586       }
01587       case 2:
01588       {
01589         if(useSquaredNorm)
01590         {
01591             return sqrt((NormType)squaredNorm());
01592         }
01593         else
01594         {
01595             NormType normMax = NumericTraits<NormType>::zero();
01596             detail::normMaxOfMultiArray(traverser_begin(), shape(), normMax, MetaInt<actual_dimension-1>());
01597             if(normMax == NumericTraits<NormType>::zero())
01598                 return normMax;
01599             NormType res  = NumericTraits<NormType>::zero();
01600             detail::sumOverMultiArray(traverser_begin(), shape(), detail::MultiArrayScaledL2Functor<NormType>(normMax),
01601                                     res, MetaInt<actual_dimension-1>());
01602             return sqrt(res)*normMax;
01603         }
01604       }
01605       default:
01606         vigra_precondition(false, "MultiArrayView::norm(): Unknown norm type.");
01607         return NumericTraits<NormType>::zero(); // unreachable
01608     }
01609 }
01610 
01611 
01612 /********************************************************/
01613 /*                                                      */
01614 /*                          norm                        */
01615 /*                                                      */
01616 /********************************************************/
01617 
01618 template <unsigned int N, class T, class C>
01619 inline typename NormTraits<MultiArrayView <N, T, C> >::SquaredNormType
01620 squaredNorm(MultiArrayView <N, T, C> const & a)
01621 {
01622     return a.squaredNorm();
01623 }
01624 
01625 template <unsigned int N, class T, class C>
01626 inline typename NormTraits<MultiArrayView <N, T, C> >::NormType
01627 norm(MultiArrayView <N, T, C> const & a)
01628 {
01629     return a.norm();
01630 }
01631 
01632 /********************************************************/
01633 /*                                                      */
01634 /*                       MultiArray                     */
01635 /*                                                      */
01636 /********************************************************/
01637 
01638 /** \brief Main <TT>MultiArray</TT> class containing the memory
01639     management.
01640 
01641 This class inherits the interface of MultiArrayView, and implements
01642 the memory ownership.
01643 MultiArray's are always unstrided, striding them creates a MultiArrayView.
01644 
01645 
01646 The template parameters are as follows
01647 \code
01648     N: the array dimension
01649 
01650     T: the type of the array elements
01651 
01652     A: the allocator used for internal storage management
01653        (default: std::allocator<T>)
01654 \endcode
01655 
01656 <b>\#include</b>
01657 <<a href="multi__array_8hxx-source.html">vigra/multi_array.hxx</a>>
01658 
01659 Namespace: vigra
01660 */
01661 template <unsigned int N, class T, class A /* default already declared above */>
01662 class MultiArray : public MultiArrayView <N, T>
01663 {
01664 
01665 public:
01666     using MultiArrayView <N, T>::actual_dimension;
01667 
01668         /** the allocator type used to allocate the memory
01669          */
01670     typedef A allocator_type;
01671 
01672         /** the view type associated with this array.
01673          */
01674     typedef MultiArrayView <N, T> view_type;
01675 
01676         /** the matrix type associated with this array.
01677          */
01678     typedef MultiArray <N, T, A> matrix_type;
01679 
01680         /** the array's value type
01681          */
01682     typedef typename view_type::value_type value_type;
01683 
01684         /** pointer type
01685          */
01686     typedef typename view_type::pointer pointer;
01687 
01688         /** const pointer type
01689          */
01690     typedef typename view_type::const_pointer const_pointer;
01691 
01692         /** reference type (result of operator[])
01693          */
01694     typedef typename view_type::reference reference;
01695 
01696         /** const reference type (result of operator[] const)
01697          */
01698     typedef typename view_type::const_reference const_reference;
01699 
01700         /** size type
01701          */
01702     typedef typename view_type::size_type size_type;
01703 
01704         /** difference type (used for multi-dimensional offsets and indices)
01705          */
01706     typedef typename view_type::difference_type difference_type;
01707 
01708         /** difference and index type for a single dimension
01709          */
01710     typedef typename view_type::difference_type_1 difference_type_1;
01711 
01712         /** traverser type
01713          */
01714     typedef typename vigra::detail::MultiIteratorChooser <
01715         UnstridedArrayTag>::template Traverser <N, T, T &, T *>::type
01716     traverser;
01717 
01718         /** traverser type to const data
01719          */
01720     typedef typename vigra::detail::MultiIteratorChooser <
01721         UnstridedArrayTag>::template Traverser <N, T, T const &, T const *>::type
01722     const_traverser;
01723 
01724         /** sequential (random access) iterator type
01725          */
01726     typedef T * iterator;
01727 
01728         /** sequential (random access) const iterator type
01729          */
01730     typedef T * const_iterator;
01731 
01732 protected:
01733 
01734     typedef typename difference_type::value_type diff_zero_t;
01735 
01736         /** the allocator used to allocate the memory
01737          */
01738     allocator_type m_alloc;
01739 
01740         /** allocate memory for s pixels, write its address into the given
01741             pointer and initialize the pixels with init.
01742         */
01743     void allocate (pointer &ptr, difference_type_1 s, const_reference init);
01744 
01745         /** allocate memory for s pixels, write its address into the given
01746             pointer and initialize the linearized pixels to the values of init.
01747         */
01748     template <class U>
01749     void allocate (pointer &ptr, difference_type_1 s, U const * init);
01750 
01751         /** allocate memory, write its address into the given
01752             pointer and initialize it by copying the data from the given MultiArrayView.
01753         */
01754     template <class U, class C>
01755     void allocate (pointer &ptr, MultiArrayView<N, U, C> const & init);
01756 
01757         /** deallocate the memory (of length s) starting at the given address.
01758          */
01759     void deallocate (pointer &ptr, difference_type_1 s);
01760 
01761     template <class U, class C>
01762     void copyOrReshape (const MultiArrayView<N, U, C> &rhs);
01763 public:
01764         /** default constructor
01765          */
01766     MultiArray ()
01767     : MultiArrayView <N, T> (difference_type (diff_zero_t(0)),
01768                              difference_type (diff_zero_t(0)), 0)
01769     {}
01770 
01771         /** construct with given allocator
01772          */
01773     MultiArray (allocator_type const & alloc)
01774     : MultiArrayView <N, T> (difference_type (diff_zero_t(0)),
01775                              difference_type (diff_zero_t(0)), 0),
01776       m_alloc(alloc)
01777     {}
01778 
01779         /** construct with given shape
01780          */
01781     explicit MultiArray (const difference_type &shape,
01782                          allocator_type const & alloc = allocator_type());
01783 
01784         /** construct from shape with an initial value
01785          */
01786     MultiArray (const difference_type &shape, const_reference init,
01787                 allocator_type const & alloc = allocator_type());
01788 
01789         /** construct from shape and copy values from the given array
01790          */
01791     MultiArray (const difference_type &shape, const_pointer init,
01792                          allocator_type const & alloc = allocator_type());
01793 
01794         /** copy constructor
01795          */
01796     MultiArray (const MultiArray &rhs)
01797     : MultiArrayView <N, T> (rhs.m_shape, rhs.m_stride, 0),
01798       m_alloc (rhs.m_alloc)
01799     {
01800         allocate (this->m_ptr, this->elementCount (), rhs.data ());
01801     }
01802 
01803         /** construct by copying from a MultiArrayView
01804          */
01805     template <class U, class C>
01806     MultiArray (const MultiArrayView<N, U, C>  &rhs,
01807                 allocator_type const & alloc = allocator_type());
01808 
01809         /** assignment.<br>
01810             If the size of \a rhs is the same as the left-hand side arrays's old size, only
01811             the data are copied. Otherwise, new storage is allocated, which invalidates all
01812             objects (array views, iterators) depending on the lhs array.
01813          */
01814     MultiArray & operator= (const MultiArray &rhs)
01815     {
01816         if (this != &rhs)
01817             this->copyOrReshape(rhs);
01818         return *this;
01819     }
01820 
01821         /** assignment from arbitrary MultiArrayView.<br>
01822             If the size of \a rhs is the same as the left-hand side arrays's old size, only
01823             the data are copied. Otherwise, new storage is allocated, which invalidates all
01824             objects (array views, iterators) depending on the lhs array.
01825          */
01826     template <class U, class C>
01827     MultiArray &operator= (const MultiArrayView<N, U, C> &rhs)
01828     {
01829         this->copyOrReshape(rhs);
01830         return *this;
01831     }
01832 
01833         /** Add-assignment from arbitrary MultiArrayView. Fails with
01834             <tt>PreconditionViolation</tt> exception when the shapes do not match.
01835          */
01836     template <class U, class C>
01837     MultiArray &operator+= (const MultiArrayView<N, U, C> &rhs)
01838     {
01839         view_type::operator+=(rhs);
01840         return *this;
01841     }
01842 
01843         /** Subtract-assignment from arbitrary MultiArrayView. Fails with
01844             <tt>PreconditionViolation</tt> exception when the shapes do not match.
01845          */
01846     template <class U, class C>
01847     MultiArray &operator-= (const MultiArrayView<N, U, C> &rhs)
01848     {
01849         view_type::operator-=(rhs);
01850         return *this;
01851     }
01852 
01853         /** Multiply-assignment from arbitrary MultiArrayView. Fails with
01854             <tt>PreconditionViolation</tt> exception when the shapes do not match.
01855          */
01856     template <class U, class C>
01857     MultiArray &operator*= (const MultiArrayView<N, U, C> &rhs)
01858     {
01859         view_type::operator*=(rhs);
01860         return *this;
01861     }
01862 
01863         /** Divide-assignment from arbitrary MultiArrayView. Fails with
01864             <tt>PreconditionViolation</tt> exception when the shapes do not match.
01865          */
01866     template <class U, class C>
01867     MultiArray &operator/= (const MultiArrayView<N, U, C> &rhs)
01868     {
01869         view_type::operator/=(rhs);
01870         return *this;
01871     }
01872 
01873         /** Add-assignment of a scalar.
01874          */
01875     MultiArray &operator+= (const T &rhs)
01876     {
01877         view_type::operator+=(rhs);
01878         return *this;
01879     }
01880 
01881         /** Subtract-assignment of a scalar.
01882          */
01883     MultiArray &operator-= (const T &rhs)
01884     {
01885         view_type::operator-=(rhs);
01886         return *this;
01887     }
01888 
01889         /** Multiply-assignment of a scalar.
01890          */
01891     MultiArray &operator*= (const T &rhs)
01892     {
01893         view_type::operator*=(rhs);
01894         return *this;
01895     }
01896 
01897         /** Divide-assignment of a scalar.
01898          */
01899     MultiArray &operator/= (const T &rhs)
01900     {
01901         view_type::operator/=(rhs);
01902         return *this;
01903     }
01904 
01905         /** destructor
01906          */
01907     ~MultiArray ()
01908     {
01909         deallocate (this->m_ptr, this->elementCount ());
01910     }
01911 
01912 
01913         /** init elements with a constant
01914          */
01915     template <class U>
01916     MultiArray & init(const U & init)
01917     {
01918         view_type::init(init);
01919         return *this;
01920     }
01921 
01922         /** change the shape and allocate new memory.<br>
01923             <em>Note:</em> this operation invalidates all dependent objects
01924             (array views and iterators)
01925          */
01926     void reshape (const difference_type &shape)
01927     {
01928         reshape (shape, NumericTraits <T>::zero ());
01929     }
01930 
01931         /** change the shape, allocate new memory and initialize it
01932             with the given value.<br>
01933             <em>Note:</em> this operation invalidates all dependent objects
01934             (array views and iterators)
01935          */
01936     void reshape (const difference_type &shape, const_reference init);
01937 
01938         /** Swap the contents with another MultiArray. This is fast,
01939             because no data are copied, but only pointers and shapes swapped.
01940             <em>Note:</em> this operation invalidates all dependent objects
01941             (array views and iterators)
01942          */
01943     void swap (MultiArray & other);
01944 
01945         /** sequential iterator pointing to the first array element.
01946          */
01947     iterator begin ()
01948     {
01949         return this->data();
01950     }
01951 
01952         /** sequential iterator pointing beyond the last array element.
01953          */
01954     iterator end ()
01955     {
01956         return this->data() + this->elementCount();
01957     }
01958 
01959         /** sequential const iterator pointing to the first array element.
01960          */
01961     const_iterator begin () const
01962     {
01963         return this->data();
01964     }
01965 
01966         /** sequential const iterator pointing beyond the last array element.
01967          */
01968     const_iterator end () const
01969     {
01970         return this->data() + this->elementCount();
01971     }
01972 
01973         /** get the allocator.
01974          */
01975     allocator_type const & allocator () const
01976     {
01977         return m_alloc;
01978     }
01979 };
01980 
01981 template <unsigned int N, class T, class A>
01982 MultiArray <N, T, A>::MultiArray (const difference_type &shape,
01983                                   allocator_type const & alloc)
01984 : MultiArrayView <N, T> (shape,
01985                          detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
01986                          0),
01987   m_alloc(alloc)
01988 {
01989     if (N == 0)
01990     {
01991         this->m_shape [0] = 1;
01992         this->m_stride [0] = 0;
01993     }
01994     allocate (this->m_ptr, this->elementCount (), NumericTraits<T>::zero ());
01995 }
01996 
01997 template <unsigned int N, class T, class A>
01998 MultiArray <N, T, A>::MultiArray (const difference_type &shape, const_reference init,
01999                                   allocator_type const & alloc)
02000 : MultiArrayView <N, T> (shape,
02001                          detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
02002                          0),
02003   m_alloc(alloc)
02004 {
02005     if (N == 0)
02006     {
02007         this->m_shape [0] = 1;
02008         this->m_stride [0] = 0;
02009     }
02010     allocate (this->m_ptr, this->elementCount (), init);
02011 }
02012 
02013 template <unsigned int N, class T, class A>
02014 MultiArray <N, T, A>::MultiArray (const difference_type &shape, const_pointer init,
02015                                   allocator_type const & alloc)
02016 : MultiArrayView <N, T> (shape,
02017                          detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (shape),
02018                          0),
02019   m_alloc(alloc)
02020 {
02021     if (N == 0)
02022     {
02023         this->m_shape [0] = 1;
02024         this->m_stride [0] = 0;
02025     }
02026     allocate (this->m_ptr, this->elementCount (), init);
02027 }
02028 
02029 template <unsigned int N, class T, class A>
02030 template <class U, class C>
02031 MultiArray <N, T, A>::MultiArray(const MultiArrayView<N, U, C>  &rhs,
02032                                  allocator_type const & alloc)
02033 : MultiArrayView <N, T> (rhs.shape(),
02034                          detail::defaultStride <MultiArrayView<N,T>::actual_dimension>(rhs.shape()),
02035                          0),
02036   m_alloc (alloc)
02037 {
02038     allocate (this->m_ptr, rhs);
02039 }
02040 
02041 template <unsigned int N, class T, class A>
02042 template <class U, class C>
02043 void
02044 MultiArray <N, T, A>::copyOrReshape(const MultiArrayView<N, U, C> &rhs)
02045 {
02046     if (this->shape() == rhs.shape())
02047         this->copy(rhs);
02048     else
02049     {
02050         MultiArray t(rhs);
02051         this->swap(t);
02052     }
02053 }
02054 
02055 template <unsigned int N, class T, class A>
02056 void MultiArray <N, T, A>::reshape (const difference_type & new_shape,
02057                                     const_reference initial)
02058 {
02059     if (N== 0)
02060     {
02061         return;
02062     }
02063     else if(new_shape == this->shape())
02064     {
02065         this->init(initial);
02066     }
02067     else
02068     {
02069         difference_type new_stride = detail::defaultStride <MultiArrayView<N,T>::actual_dimension> (new_shape);
02070         difference_type_1 new_size = new_shape [MultiArrayView<N,T>::actual_dimension-1] * new_stride [MultiArrayView<N,T>::actual_dimension-1];
02071         T *new_ptr;
02072         allocate (new_ptr, new_size, initial);
02073         deallocate (this->m_ptr, this->elementCount ());
02074         this->m_ptr = new_ptr;
02075         this->m_shape = new_shape;
02076         this->m_stride = new_stride;
02077     }
02078 }
02079 
02080 
02081 template <unsigned int N, class T, class A>
02082 inline void
02083 MultiArray <N, T, A>::swap (MultiArray <N, T, A> & other)
02084 {
02085     if (this == &other)
02086         return;
02087     std::swap(this->m_shape,  other.m_shape);
02088     std::swap(this->m_stride, other.m_stride);
02089     std::swap(this->m_ptr,    other.m_ptr);
02090     std::swap(this->m_alloc,  other.m_alloc);
02091 }
02092 
02093 template <unsigned int N, class T, class A>
02094 void MultiArray <N, T, A>::allocate (pointer & ptr, difference_type_1 s,
02095                                      const_reference init)
02096 {
02097     ptr = m_alloc.allocate ((typename A::size_type)s);
02098     difference_type_1 i;
02099     try {
02100         for (i = 0; i < s; ++i)
02101             m_alloc.construct (ptr + i, init);
02102     }
02103     catch (...) {
02104         for (difference_type_1 j = 0; j < i; ++j)
02105             m_alloc.destroy (ptr + j);
02106         m_alloc.deallocate (ptr, (typename A::size_type)s);
02107         throw;
02108     }
02109 }
02110 
02111 template <unsigned int N, class T, class A>
02112 template <class U>
02113 void MultiArray <N, T, A>::allocate (pointer & ptr, difference_type_1 s,
02114                                      U const * init)
02115 {
02116     ptr = m_alloc.allocate ((typename A::size_type)s);
02117     difference_type_1 i;
02118     try {
02119         for (i = 0; i < s; ++i, ++init)
02120             m_alloc.construct (ptr + i, *init);
02121     }
02122     catch (...) {
02123         for (difference_type_1 j = 0; j < i; ++j)
02124             m_alloc.destroy (ptr + j);
02125         m_alloc.deallocate (ptr, (typename A::size_type)s);
02126         throw;
02127     }
02128 }
02129 
02130 template <unsigned int N, class T, class A>
02131 template <class U, class C>
02132 void MultiArray <N, T, A>::allocate (pointer & ptr, MultiArrayView<N, U, C> const & init)
02133 {
02134     difference_type_1 s = init.elementCount();
02135     ptr = m_alloc.allocate ((typename A::size_type)s);
02136     pointer p = ptr;
02137     try {
02138         detail::uninitializedCopyMultiArrayData(init.traverser_begin(), init.shape(),
02139                                                 p, m_alloc, MetaInt<actual_dimension-1>());
02140     }
02141     catch (...) {
02142         for (pointer pp = ptr; pp < p; ++pp)
02143             m_alloc.destroy (pp);
02144         m_alloc.deallocate (ptr, (typename A::size_type)s);
02145         throw;
02146     }
02147 }
02148 
02149 template <unsigned int N, class T, class A>
02150 inline void MultiArray <N, T, A>::deallocate (pointer & ptr, difference_type_1 s)
02151 {
02152     if (ptr == 0)
02153         return;
02154     for (difference_type_1 i = 0; i < s; ++i)
02155         m_alloc.destroy (ptr + i);
02156     m_alloc.deallocate (ptr, (typename A::size_type)s);
02157     ptr = 0;
02158 }
02159 
02160 /********************************************************/
02161 /*                                                      */
02162 /*              argument object factories               */
02163 /*                                                      */
02164 /********************************************************/
02165 
02166 template <unsigned int N, class T, class C>
02167 inline triple<typename MultiArrayView<N,T,C>::const_traverser,
02168               typename MultiArrayView<N,T,C>::difference_type,
02169               typename AccessorTraits<T>::default_const_accessor >
02170 srcMultiArrayRange( MultiArrayView<N,T,C> const & array )
02171 {
02172     return triple<typename MultiArrayView<N,T,C>::const_traverser,
02173                   typename MultiArrayView<N,T,C>::difference_type,
02174                   typename AccessorTraits<T>::default_const_accessor >
02175       ( array.traverser_begin(),
02176         array.shape(),
02177         typename AccessorTraits<T>::default_const_accessor() );
02178 }
02179 
02180 template <unsigned int N, class T, class C, class Accessor>
02181 inline triple<typename MultiArrayView<N,T,C>::const_traverser,
02182               typename MultiArrayView<N,T,C>::difference_type,
02183               Accessor >
02184 srcMultiArrayRange( MultiArrayView<N,T,C> const & array, Accessor a )
02185 {
02186     return triple<typename MultiArrayView<N,T,C>::const_traverser,
02187                   typename MultiArrayView<N,T,C>::difference_type,
02188                   Accessor >
02189       ( array.traverser_begin(),
02190         array.shape(),
02191         a);
02192 }
02193 
02194 template <unsigned int N, class T, class C>
02195 inline pair<typename MultiArrayView<N,T,C>::const_traverser,
02196             typename AccessorTraits<T>::default_const_accessor >
02197 srcMultiArray( MultiArrayView<N,T,C> const & array )
02198 {
02199     return pair<typename MultiArrayView<N,T,C>::const_traverser,
02200                 typename AccessorTraits<T>::default_const_accessor >
02201       ( array.traverser_begin(),
02202         typename AccessorTraits<T>::default_const_accessor() );
02203 }
02204 
02205 template <unsigned int N, class T, class C, class Accessor>
02206 inline pair<typename MultiArrayView<N,T,C>::const_traverser,
02207             Accessor >
02208 srcMultiArray( MultiArrayView<N,T,C> const & array, Accessor a )
02209 {
02210     return pair<typename MultiArrayView<N,T,C>::const_traverser,
02211                 Accessor >
02212       ( array.traverser_begin(), a );
02213 }
02214 
02215 template <unsigned int N, class T, class C>
02216 inline triple<typename MultiArrayView<N,T,C>::traverser,
02217               typename MultiArrayView<N,T,C>::difference_type,
02218               typename AccessorTraits<T>::default_accessor >
02219 destMultiArrayRange( MultiArrayView<N,T,C> & array )
02220 {
02221     return triple<typename MultiArrayView<N,T,C>::traverser,
02222                   typename MultiArrayView<N,T,C>::difference_type,
02223                   typename AccessorTraits<T>::default_accessor >
02224       ( array.traverser_begin(),
02225         array.shape(),
02226         typename AccessorTraits<T>::default_accessor() );
02227 }
02228 
02229 template <unsigned int N, class T, class C, class Accessor>
02230 inline triple<typename MultiArrayView<N,T,C>::traverser,
02231               typename MultiArrayView<N,T,C>::difference_type,
02232               Accessor >
02233 destMultiArrayRange( MultiArrayView<N,T,C> & array, Accessor a )
02234 {
02235     return triple<typename MultiArrayView<N,T,C>::traverser,
02236                   typename MultiArrayView<N,T,C>::difference_type,
02237                   Accessor >
02238       ( array.traverser_begin(),
02239         array.shape(),
02240         a );
02241 }
02242 
02243 template <unsigned int N, class T, class C>
02244 inline pair<typename MultiArrayView<N,T,C>::traverser,
02245             typename AccessorTraits<T>::default_accessor >
02246 destMultiArray( MultiArrayView<N,T,C> & array )
02247 {
02248     return pair<typename MultiArrayView<N,T,C>::traverser,
02249                 typename AccessorTraits<T>::default_accessor >
02250         ( array.traverser_begin(),
02251           typename AccessorTraits<T>::default_accessor() );
02252 }
02253 
02254 template <unsigned int N, class T, class C, class Accessor>
02255 inline pair<typename MultiArrayView<N,T,C>::traverser,
02256             Accessor >
02257 destMultiArray( MultiArrayView<N,T,C> & array, Accessor a )
02258 {
02259     return pair<typename MultiArrayView<N,T,C>::traverser,
02260                 Accessor >
02261         ( array.traverser_begin(), a );
02262 }
02263 
02264 /********************************************************/
02265 /*                                                      */
02266 /*                  makeBasicImageView                  */
02267 /*                                                      */
02268 /********************************************************/
02269 
02270 /** \addtogroup MultiArrayToImage Wrap a \ref vigra::MultiArrayView in
02271                                   a \ref vigra::BasicImageView
02272 */
02273 //@{
02274 /** Create a \ref vigra::BasicImageView from an unstrided 2-dimensional
02275     \ref vigra::MultiArrayView.
02276 
02277     The \ref vigra::BasicImageView will have the same <tt>value_type </tt>
02278     as the original \ref vigra::MultiArrayView.
02279 */
02280 template <class T>
02281 BasicImageView <T>
02282 makeBasicImageView (MultiArrayView <2, T, UnstridedArrayTag> const &array)
02283 {
02284     return BasicImageView <T> (array.data (), array.shape (0),
02285                                array.shape (1));
02286 }
02287 
02288 /** Create a \ref vigra::BasicImageView from a 3-dimensional
02289     \ref vigra::MultiArray.
02290 
02291     This wrapper flattens the two innermost dimensions of the array
02292     into single rows of the resulting image.
02293     The \ref vigra::BasicImageView will have the same <tt>value_type </tt>
02294     as the original \ref vigra::MultiArray.
02295 */
02296 template <class T>
02297 BasicImageView <T>
02298 makeBasicImageView (MultiArray <3, T> const &array)
02299 {
02300     return BasicImageView <T> (array.data (),
02301                                array.shape (0)*array.shape (1), array.shape (2));
02302 }
02303 
02304 /** Create a \ref vigra::BasicImageView from a 3-dimensional
02305     \ref vigra::MultiArray.
02306 
02307     This wrapper only works if <tt>T</tt> is a scalar type and the
02308     array's innermost dimension has size 3. It then re-interprets
02309     the data array as a 2-dimensional array with value_type
02310     <tt>RGBValue<T></tt>.
02311 */
02312 template <class T>
02313 BasicImageView <RGBValue<T> >
02314 makeRGBImageView (MultiArray<3, T> const &array)
02315 {
02316     vigra_precondition (
02317         array.shape (0) == 3, "makeRGBImageView(): array.shape(0) must be 3.");
02318     return BasicImageView <RGBValue<T> > (
02319         reinterpret_cast <RGBValue <T> *> (array.data ()),
02320         array.shape (1), array.shape (2));
02321 }
02322 
02323 //@}
02324 
02325 } // namespace vigra
02326 
02327 #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.6.0 (13 Aug 2008)