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