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

vigra/numpy_array.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*       Copyright 2009 by Ullrich Koethe and Hans Meine                */
00004 /*                                                                      */
00005 /*    This file is part of the VIGRA computer vision library.           */
00006 /*    The VIGRA Website is                                              */
00007 /*        http://hci.iwr.uni-heidelberg.de/vigra/                       */
00008 /*    Please direct questions, bug reports, and contributions to        */
00009 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
00010 /*        vigra@informatik.uni-hamburg.de                               */
00011 /*                                                                      */
00012 /*    Permission is hereby granted, free of charge, to any person       */
00013 /*    obtaining a copy of this software and associated documentation    */
00014 /*    files (the "Software"), to deal in the Software without           */
00015 /*    restriction, including without limitation the rights to use,      */
00016 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00017 /*    sell copies of the Software, and to permit persons to whom the    */
00018 /*    Software is furnished to do so, subject to the following          */
00019 /*    conditions:                                                       */
00020 /*                                                                      */
00021 /*    The above copyright notice and this permission notice shall be    */
00022 /*    included in all copies or substantial portions of the             */
00023 /*    Software.                                                         */
00024 /*                                                                      */
00025 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00026 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00027 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00028 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00029 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00030 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00031 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00032 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00033 /*                                                                      */
00034 /************************************************************************/
00035 
00036 #ifndef VIGRA_NUMPY_ARRAY_HXX
00037 #define VIGRA_NUMPY_ARRAY_HXX
00038 
00039 #include <Python.h>
00040 #include <string>
00041 #include <iostream>
00042 #include <numpy/arrayobject.h>
00043 #include "multi_array.hxx"
00044 #include "array_vector.hxx"
00045 #include "python_utility.hxx"
00046 #include "numpy_array_traits.hxx"
00047 #include "numpy_array_taggedshape.hxx"
00048 
00049 // NumPy function called by NumPy’s import_array() macro (and our import_vigranumpy() below)
00050 int _import_array();
00051 
00052 namespace vigra {
00053 
00054 static inline void import_vigranumpy()
00055 {
00056     // roughly equivalent to import_array():
00057     if(_import_array() < 0)
00058         pythonToCppException(0);
00059 
00060     // in addition, import vigra.vigranumpycore:
00061     python_ptr module(PyImport_ImportModule("vigra.vigranumpycore"), python_ptr::keep_count);
00062     pythonToCppException(module);
00063 }
00064 
00065 /********************************************************/
00066 /*                                                      */
00067 /*               MultibandVectorAccessor                */
00068 /*                                                      */
00069 /********************************************************/
00070 
00071 template <class T>
00072 class MultibandVectorAccessor
00073 {
00074     MultiArrayIndex size_, stride_;
00075 
00076   public:
00077     MultibandVectorAccessor(MultiArrayIndex size, MultiArrayIndex stride)
00078     : size_(size),
00079       stride_(stride)
00080     {}
00081 
00082 
00083     typedef Multiband<T> value_type;
00084 
00085         /** the vector's value_type
00086         */
00087     typedef T component_type;
00088 
00089     typedef VectorElementAccessor<MultibandVectorAccessor<T> > ElementAccessor;
00090 
00091         /** Read the component data at given vector index
00092             at given iterator position
00093         */
00094     template <class ITERATOR>
00095     component_type const & getComponent(ITERATOR const & i, int idx) const
00096     {
00097         return *(&*i+idx*stride_);
00098     }
00099 
00100         /** Set the component data at given vector index
00101             at given iterator position. The type <TT>V</TT> of the passed
00102             in <TT>value</TT> is automatically converted to <TT>component_type</TT>.
00103             In case of a conversion floating point -> integral this includes rounding and clipping.
00104         */
00105     template <class V, class ITERATOR>
00106     void setComponent(V const & value, ITERATOR const & i, int idx) const
00107     {
00108         *(&*i+idx*stride_) = detail::RequiresExplicitCast<component_type>::cast(value);
00109     }
00110 
00111         /** Read the component data at given vector index
00112             at an offset of given iterator position
00113         */
00114     template <class ITERATOR, class DIFFERENCE>
00115     component_type const & getComponent(ITERATOR const & i, DIFFERENCE const & diff, int idx) const
00116     {
00117         return *(&i[diff]+idx*stride_);
00118     }
00119 
00120     /** Set the component data at given vector index
00121         at an offset of given iterator position. The type <TT>V</TT> of the passed
00122         in <TT>value</TT> is automatically converted to <TT>component_type</TT>.
00123             In case of a conversion floating point -> integral this includes rounding and clipping.
00124     */
00125     template <class V, class ITERATOR, class DIFFERENCE>
00126     void
00127     setComponent(V const & value, ITERATOR const & i, DIFFERENCE const & diff, int idx) const
00128     {
00129         *(&i[diff]+idx*stride_) = detail::RequiresExplicitCast<component_type>::cast(value);
00130     }
00131 
00132     template <class U>
00133     MultiArrayIndex size(U) const
00134     {
00135         return size_;
00136     }
00137 };
00138 
00139 /********************************************************/
00140 
00141 template <class TYPECODE> // pseudo-template to avoid inline expansion of the function
00142                           // will always be NPY_TYPES
00143 PyObject *
00144 constructArray(TaggedShape tagged_shape, TYPECODE typeCode, bool init,
00145                python_ptr arraytype = python_ptr());
00146 
00147 /********************************************************/
00148 /*                                                      */
00149 /*                    NumpyAnyArray                     */
00150 /*                                                      */
00151 /********************************************************/
00152 
00153 /** Wrapper class for a Python array.
00154 
00155     This class stores a reference-counted pointer to an Python numpy array object,
00156     i.e. an object where <tt>PyArray_Check(object)</tt> returns true (in Python, the
00157     object is then a subclass of <tt>numpy.ndarray</tt>). This class is mainly used
00158     as a smart pointer to these arrays, but some basic access and conversion functions
00159     are also provided.
00160 
00161     <b>\#include</b> <vigra/numpy_array.hxx><br>
00162     Namespace: vigra
00163 */
00164 class NumpyAnyArray
00165 {
00166   protected:
00167     python_ptr pyArray_;
00168 
00169   public:
00170 
00171         /// difference type
00172     typedef ArrayVector<npy_intp> difference_type;
00173 
00174     static python_ptr getArrayTypeObject()
00175     {
00176         return detail::getArrayTypeObject();
00177     }
00178 
00179     static std::string defaultOrder(std::string defaultValue = "C")
00180     {
00181         return detail::defaultOrder(defaultValue);
00182     }
00183 
00184     static python_ptr defaultAxistags(int ndim, std::string order = "")
00185     {
00186         return detail::defaultAxistags(ndim, order);
00187     }
00188 
00189     static python_ptr emptyAxistags(int ndim)
00190     {
00191         return detail::emptyAxistags(ndim);
00192     }
00193 
00194         /**
00195          Construct from a Python object. If \a obj is NULL, or is not a subclass
00196          of numpy.ndarray, the resulting NumpyAnyArray will have no data (i.e.
00197          hasData() returns false). Otherwise, it creates a new reference to the array
00198          \a obj, unless \a createCopy is true, where a new array is created by calling
00199          the C-equivalent of obj->copy().
00200          */
00201     explicit NumpyAnyArray(PyObject * obj = 0, bool createCopy = false, PyTypeObject * type = 0)
00202     {
00203         if(obj == 0)
00204             return;
00205         vigra_precondition(type == 0 || PyType_IsSubtype(type, &PyArray_Type),
00206              "NumpyAnyArray(obj, createCopy, type): type must be numpy.ndarray or a subclass thereof.");
00207         if(createCopy)
00208             makeCopy(obj, type);
00209         else
00210             vigra_precondition(makeReference(obj, type), "NumpyAnyArray(obj): obj isn't a numpy array.");
00211     }
00212 
00213         /**
00214          Copy constructor. By default, it creates a new reference to the array
00215          \a other. When \a createCopy is true, a new array is created by calling
00216          the C-equivalent of other.copy().
00217          */
00218     NumpyAnyArray(NumpyAnyArray const & other, bool createCopy = false, PyTypeObject * type = 0)
00219     {
00220         if(!other.hasData())
00221             return;
00222         vigra_precondition(type == 0 || PyType_IsSubtype(type, &PyArray_Type),
00223              "NumpyAnyArray(obj, createCopy, type): type must be numpy.ndarray or a subclass thereof.");
00224         if(createCopy)
00225             makeCopy(other.pyObject(), type);
00226         else
00227             makeReference(other.pyObject(), type);
00228     }
00229 
00230     // auto-generated destructor is ok
00231 
00232         /**
00233          * Assignment operator. If this is already a view with data
00234          * (i.e. hasData() is true) and the shapes match, the RHS
00235          * array contents are copied via the C-equivalent of
00236          * 'self[...] = other[...]'. If the shapes don't matched,
00237          * broadcasting is tried on the trailing (i.e. channel)
00238          * dimension.
00239          * If the LHS is an empty view, assignment is identical to
00240          * makeReference(other.pyObject()).
00241          */
00242     NumpyAnyArray & operator=(NumpyAnyArray const & other)
00243     {
00244         if(hasData())
00245         {
00246             vigra_precondition(other.hasData(),
00247                 "NumpyArray::operator=(): Cannot assign from empty array.");
00248 
00249             python_ptr arraytype = getArrayTypeObject();
00250             python_ptr f(PyString_FromString("_copyValuesImpl"), python_ptr::keep_count);
00251             if(PyObject_HasAttr(arraytype, f))
00252             {
00253                 python_ptr res(PyObject_CallMethodObjArgs(arraytype, f.get(),
00254                                                           pyArray_.get(), other.pyArray_.get(), NULL),
00255                                python_ptr::keep_count);
00256                 vigra_postcondition(res.get() != 0,
00257                        "NumpyArray::operator=(): VigraArray._copyValuesImpl() failed.");
00258             }
00259             else
00260             {
00261                 PyArrayObject * sarray = (PyArrayObject *)pyArray_.get();
00262                 PyArrayObject * tarray = (PyArrayObject *)other.pyArray_.get();
00263 
00264                 if(PyArray_CopyInto(tarray, sarray) == -1)
00265                     pythonToCppException(0);
00266             }
00267         }
00268         else
00269         {
00270             pyArray_ = other.pyArray_;
00271         }
00272         return *this;
00273     }
00274 
00275         /**
00276          Returns the number of dimensions of this array, or 0 if
00277          hasData() is false.
00278          */
00279     MultiArrayIndex ndim() const
00280     {
00281         if(hasData())
00282             return PyArray_NDIM(pyObject());
00283         return 0;
00284     }
00285 
00286         /**
00287          Returns the number of spatial dimensions of this array, or 0 if
00288          hasData() is false. If the enclosed Python array does not define
00289          the attribute spatialDimensions, ndim() is returned.
00290          */
00291     MultiArrayIndex spatialDimensions() const
00292     {
00293         if(!hasData())
00294             return 0;
00295         return pythonGetAttr(pyObject(), "spatialDimensions", ndim());
00296     }
00297 
00298     bool hasChannelAxis() const
00299     {
00300         if(!hasData())
00301             return false;
00302         return channelIndex() == ndim();
00303     }
00304 
00305     MultiArrayIndex channelIndex() const
00306     {
00307         if(!hasData())
00308             return 0;
00309         return pythonGetAttr(pyObject(), "channelIndex", ndim());
00310     }
00311 
00312     MultiArrayIndex innerNonchannelIndex() const
00313     {
00314         if(!hasData())
00315             return 0;
00316         return pythonGetAttr(pyObject(), "innerNonchannelIndex", ndim());
00317     }
00318 
00319         /**
00320          Returns the shape of this array. The size of
00321          the returned shape equals ndim().
00322          */
00323     difference_type shape() const
00324     {
00325         if(hasData())
00326             return difference_type(PyArray_DIMS(pyObject()), PyArray_DIMS(pyObject()) + ndim());
00327         return difference_type();
00328     }
00329 
00330         /** Compute the ordering of the strides of this array.
00331             The result is describes the current permutation of the axes relative
00332             to an ascending stride order.
00333         */
00334     difference_type strideOrdering() const
00335     {
00336         if(!hasData())
00337             return difference_type();
00338         MultiArrayIndex N = ndim();
00339         difference_type stride(PyArray_STRIDES(pyObject()), PyArray_STRIDES(pyObject()) + N),
00340                         permutation(N);
00341         for(MultiArrayIndex k=0; k<N; ++k)
00342             permutation[k] = k;
00343         for(MultiArrayIndex k=0; k<N-1; ++k)
00344         {
00345             MultiArrayIndex smallest = k;
00346             for(MultiArrayIndex j=k+1; j<N; ++j)
00347             {
00348                 if(stride[j] < stride[smallest])
00349                     smallest = j;
00350             }
00351             if(smallest != k)
00352             {
00353                 std::swap(stride[k], stride[smallest]);
00354                 std::swap(permutation[k], permutation[smallest]);
00355             }
00356         }
00357         difference_type ordering(N);
00358         for(MultiArrayIndex k=0; k<N; ++k)
00359             ordering[permutation[k]] = k;
00360         return ordering;
00361     }
00362 
00363         // /**
00364          // Returns the the permutation that will transpose this array into
00365          // canonical ordering (currently: F-order). The size of
00366          // the returned permutation equals ndim().
00367          // */
00368     // difference_type permutationToNormalOrder() const
00369     // {
00370         // if(!hasData())
00371             // return difference_type();
00372 
00373         // // difference_type res(detail::getAxisPermutationImpl(pyArray_,
00374                                                // // "permutationToNormalOrder", true));
00375         // difference_type res;
00376         // detail::getAxisPermutationImpl(res, pyArray_, "permutationToNormalOrder", true);
00377         // if(res.size() == 0)
00378         // {
00379             // res.resize(ndim());
00380             // linearSequence(res.begin(), res.end(), ndim()-1, MultiArrayIndex(-1));
00381         // }
00382         // return res;
00383     // }
00384 
00385         /**
00386          Returns the value type of the elements in this array, or -1
00387          when hasData() is false.
00388          */
00389     int dtype() const
00390     {
00391         if(hasData())
00392             return PyArray_DESCR(pyObject())->type_num;
00393         return -1;
00394     }
00395 
00396         /**
00397          * Return the AxisTags of this array or a NULL pointer when the attribute
00398            'axistags' is missing in the Python object or this array has no data.
00399          */
00400     python_ptr axistags() const
00401     {
00402         static python_ptr key(PyString_FromString("axistags"), python_ptr::keep_count);
00403 
00404         python_ptr axistags;
00405         if(pyObject())
00406         {
00407             axistags.reset(PyObject_GetAttr(pyObject(), key), python_ptr::keep_count);
00408             if(!axistags)
00409                 PyErr_Clear();
00410         }
00411         return axistags;
00412     }
00413 
00414         /**
00415          * Return a borrowed reference to the internal PyArrayObject.
00416          */
00417     PyArrayObject * pyArray() const
00418     {
00419         return (PyArrayObject *)pyArray_.get();
00420     }
00421 
00422         /**
00423          * Return a borrowed reference to the internal PyArrayObject
00424          * (see pyArray()), cast to PyObject for your convenience.
00425          */
00426     PyObject * pyObject() const
00427     {
00428         return pyArray_.get();
00429     }
00430 
00431         /**
00432            Reset the NumpyAnyArray to the given object. If \a obj is a numpy array object,
00433            a new reference to that array is created, and the function returns
00434            true. Otherwise, it returns false and the NumpyAnyArray remains unchanged.
00435            If \a type is given, the new reference will be a view with that type, provided
00436            that \a type is a numpy ndarray or a subclass thereof. Otherwise, an
00437            exception is thrown.
00438          */
00439     bool makeReference(PyObject * obj, PyTypeObject * type = 0)
00440     {
00441         if(obj == 0 || !PyArray_Check(obj))
00442             return false;
00443         if(type != 0)
00444         {
00445             vigra_precondition(PyType_IsSubtype(type, &PyArray_Type) != 0,
00446                 "NumpyAnyArray::makeReference(obj, type): type must be numpy.ndarray or a subclass thereof.");
00447             obj = PyArray_View((PyArrayObject*)obj, 0, type);
00448             pythonToCppException(obj);
00449         }
00450         pyArray_.reset(obj);
00451         return true;
00452     }
00453 
00454         /**
00455            Create a copy of the given array object. If \a obj is a numpy array object,
00456            a copy is created via the C-equivalent of 'obj->copy()'. If
00457            this call fails, or obj was not an array, an exception is thrown
00458            and the NumpyAnyArray remains unchanged.
00459          */
00460     void makeCopy(PyObject * obj, PyTypeObject * type = 0)
00461     {
00462         vigra_precondition(obj && PyArray_Check(obj),
00463              "NumpyAnyArray::makeCopy(obj): obj is not an array.");
00464         vigra_precondition(type == 0 || PyType_IsSubtype(type, &PyArray_Type),
00465              "NumpyAnyArray::makeCopy(obj, type): type must be numpy.ndarray or a subclass thereof.");
00466         python_ptr array(PyArray_NewCopy((PyArrayObject*)obj, NPY_ANYORDER), python_ptr::keep_count);
00467         pythonToCppException(array);
00468         makeReference(array, type);
00469     }
00470 
00471          /**
00472            Check whether this NumpyAnyArray actually points to a Python array.
00473          */
00474     bool hasData() const
00475     {
00476         return pyArray_ != 0;
00477     }
00478 };
00479 
00480 /********************************************************/
00481 /*                                                      */
00482 /*                    constructArray                    */
00483 /*                                                      */
00484 /********************************************************/
00485 
00486 namespace detail {
00487 
00488 inline bool
00489 nontrivialPermutation(ArrayVector<npy_intp> const & p)
00490 {
00491     for(unsigned int k=0; k<p.size(); ++k)
00492         if(p[k] != k)
00493             return true;
00494     return false;
00495 }
00496 
00497 } // namespace detail
00498 
00499 template <class TYPECODE> // pseudo-template to avoid inline expansion of the function
00500                           // will always be NPY_TYPES
00501 PyObject *
00502 constructArray(TaggedShape tagged_shape, TYPECODE typeCode, bool init, python_ptr arraytype)
00503 {
00504     ArrayVector<npy_intp> shape = finalizeTaggedShape(tagged_shape);
00505     PyAxisTags axistags(tagged_shape.axistags);
00506 
00507     int ndim = (int)shape.size();
00508     ArrayVector<npy_intp> inverse_permutation;
00509     int order = 1; // Fortran order
00510 
00511     if(axistags)
00512     {
00513         if(!arraytype)
00514             arraytype = NumpyAnyArray::getArrayTypeObject();
00515 
00516         inverse_permutation = axistags.permutationFromNormalOrder();
00517         vigra_precondition(ndim == (int)inverse_permutation.size(),
00518                      "axistags.permutationFromNormalOrder(): permutation has wrong size.");
00519     }
00520     else
00521     {
00522         arraytype = python_ptr((PyObject*)&PyArray_Type);
00523         order = 0; // C order
00524     }
00525 
00526 //    std::cerr << "constructArray: " << shape << "\n" << inverse_permutation << "\n";
00527 
00528     python_ptr array(PyArray_New((PyTypeObject *)arraytype.get(), ndim, shape.begin(),
00529                                   typeCode, 0, 0, 0, order, 0),
00530                      python_ptr::keep_count);
00531     pythonToCppException(array);
00532 
00533     if(detail::nontrivialPermutation(inverse_permutation))
00534     {
00535         PyArray_Dims permute = { inverse_permutation.begin(), ndim };
00536         array = python_ptr(PyArray_Transpose((PyArrayObject*)array.get(), &permute),
00537                            python_ptr::keep_count);
00538         pythonToCppException(array);
00539     }
00540 
00541     if(arraytype != (PyObject*)&PyArray_Type && axistags)
00542         pythonToCppException(PyObject_SetAttrString(array, "axistags", axistags.axistags) != -1);
00543 
00544     if(init)
00545         PyArray_FILLWBYTE((PyArrayObject *)array.get(), 0);
00546 
00547     return array.release();
00548 }
00549 
00550 // FIXME: reimplement in terms of TaggedShape?
00551 template <class TINY_VECTOR>
00552 inline
00553 python_ptr constructNumpyArrayFromData(
00554     TINY_VECTOR const & shape, npy_intp *strides,
00555     NPY_TYPES typeCode, void *data)
00556 {
00557     ArrayVector<npy_intp> pyShape(shape.begin(), shape.end());
00558 
00559     python_ptr array(PyArray_New(&PyArray_Type, shape.size(), pyShape.begin(),
00560                                  typeCode, strides, data, 0, NPY_WRITEABLE, 0),
00561                      python_ptr::keep_count);
00562     pythonToCppException(array);
00563 
00564     return array;
00565 }
00566 
00567 /********************************************************/
00568 /*                                                      */
00569 /*                     NumpyArray                       */
00570 /*                                                      */
00571 /********************************************************/
00572 
00573 /** Provide the MultiArrayView interface for a Python array.
00574 
00575     This class inherits from both \ref vigra::MultiArrayView and \ref vigra::NumpyAnyArray
00576     in order to support easy and save application of VIGRA functions to Python arrays.
00577 
00578     <b>\#include</b> <vigra/numpy_array.hxx><br>
00579     Namespace: vigra
00580 */
00581 template <unsigned int N, class T, class Stride = StridedArrayTag>
00582 class NumpyArray
00583 : public MultiArrayView<N, typename NumpyArrayTraits<N, T, Stride>::value_type, Stride>,
00584   public NumpyAnyArray
00585 {
00586   public:
00587     typedef NumpyArrayTraits<N, T, Stride> ArrayTraits;
00588     typedef typename ArrayTraits::dtype dtype;
00589     typedef T pseudo_value_type;
00590 
00591     static NPY_TYPES const typeCode = ArrayTraits::typeCode;
00592 
00593         /** the view type associated with this array.
00594          */
00595     typedef MultiArrayView<N, typename ArrayTraits::value_type, Stride> view_type;
00596 
00597     enum { actual_dimension = view_type::actual_dimension };
00598 
00599         /** the array's value type
00600          */
00601     typedef typename view_type::value_type value_type;
00602 
00603         /** pointer type
00604          */
00605     typedef typename view_type::pointer pointer;
00606 
00607         /** const pointer type
00608          */
00609     typedef typename view_type::const_pointer const_pointer;
00610 
00611         /** reference type (result of operator[])
00612          */
00613     typedef typename view_type::reference reference;
00614 
00615         /** const reference type (result of operator[] const)
00616          */
00617     typedef typename view_type::const_reference const_reference;
00618 
00619         /** size type
00620          */
00621     typedef typename view_type::size_type size_type;
00622 
00623         /** difference type (used for multi-dimensional offsets and indices)
00624          */
00625     typedef typename view_type::difference_type difference_type;
00626 
00627         /** difference and index type for a single dimension
00628          */
00629     typedef typename view_type::difference_type_1 difference_type_1;
00630 
00631         /** type of an array specifying an axis permutation
00632          */
00633     typedef typename NumpyAnyArray::difference_type permutation_type;
00634 
00635         /** traverser type
00636          */
00637     typedef typename view_type::traverser traverser;
00638 
00639         /** traverser type to const data
00640          */
00641     typedef typename view_type::const_traverser const_traverser;
00642 
00643         /** sequential (random access) iterator type
00644          */
00645     typedef typename view_type::iterator iterator;
00646 
00647         /** sequential (random access) const iterator type
00648          */
00649     typedef typename view_type::const_iterator const_iterator;
00650 
00651     using view_type::shape;   // resolve ambiguity of multiple inheritance
00652     using view_type::hasData; // resolve ambiguity of multiple inheritance
00653     using view_type::strideOrdering; // resolve ambiguity of multiple inheritance
00654 
00655   protected:
00656 
00657     // this function assumes that pyArray_ has already been set, and compatibility been checked
00658     void setupArrayView();
00659 
00660     static python_ptr init(difference_type const & shape, bool init = true,
00661                            std::string const & order = "")
00662     {
00663         vigra_precondition(order == "" || order == "C" || order == "F" ||
00664                            order == "V" || order == "A",
00665             "NumpyArray.init(): order must be in ['C', 'F', 'V', 'A', ''].");
00666         return python_ptr(constructArray(ArrayTraits::taggedShape(shape, order), typeCode, init),
00667                           python_ptr::keep_count);
00668     }
00669 
00670   public:
00671 
00672     using view_type::init;
00673 
00674         /**
00675          * Construct from a given PyObject pointer. When the given
00676          * python object is NULL, the internal python array will be
00677          * NULL and hasData() will return false.
00678          *
00679          * Otherwise, the function attempts to create a
00680          * new reference to the given Python object, unless
00681          * copying is forced by setting \a createCopy to true.
00682          * If either of this fails, the function throws an exception.
00683          * This will not happen if isReferenceCompatible(obj) (in case
00684          * of creating a new reference) or isCopyCompatible(obj)
00685          * (in case of copying) have returned true beforehand.
00686          */
00687     explicit NumpyArray(PyObject *obj = 0, bool createCopy = false)
00688     {
00689         if(obj == 0)
00690             return;
00691         if(createCopy)
00692             makeCopy(obj);
00693         else
00694             vigra_precondition(makeReference(obj),
00695                   "NumpyArray(obj): Cannot construct from incompatible array.");
00696     }
00697 
00698        /**
00699          * Copy constructor; does not copy the memory, but creates a
00700          * new reference to the same underlying python object, unless
00701          * a copy is forced by setting \a createCopy to true.
00702          * (If the source object has no data, this one will have
00703          * no data, too.)
00704          */
00705     NumpyArray(const NumpyArray &other, bool createCopy = false)
00706     : view_type(),
00707       NumpyAnyArray()
00708     {
00709         if(!other.hasData())
00710             return;
00711         if(createCopy)
00712             makeCopy(other.pyObject());
00713         else
00714             makeReferenceUnchecked(other.pyObject());
00715     }
00716 
00717        /**
00718          * Allocate new memory and copy data from a MultiArrayView.
00719          */
00720     template <class U, class S>
00721     explicit NumpyArray(const MultiArrayView<N, U, S> &other)
00722     {
00723         if(!other.hasData())
00724             return;
00725         vigra_postcondition(makeReference(init(other.shape(), false)),
00726                   "NumpyArray(MultiArrayView): Python constructor did not produce a compatible array.");
00727         view_type::operator=(other);
00728     }
00729 
00730         /**
00731          * Construct a new array object, allocating an internal python
00732          * ndarray of the given shape in the given order (default: VIGRA order), initialized
00733          * with zeros.
00734          *
00735          * An exception is thrown when construction fails.
00736          */
00737     explicit NumpyArray(difference_type const & shape, std::string const & order = "")
00738     {
00739         vigra_postcondition(makeReference(init(shape, true, order)),
00740                      "NumpyArray(shape): Python constructor did not produce a compatible array.");
00741     }
00742 
00743         /**
00744          * Construct a new array object, allocating an internal python
00745          * ndarray according to the given tagged shape, initialized with zeros.
00746          *
00747          * An exception is thrown when construction fails.
00748          */
00749     explicit NumpyArray(TaggedShape const & tagged_shape)
00750     {
00751         reshapeIfEmpty(tagged_shape,
00752            "NumpyArray(tagged_shape): Python constructor did not produce a compatible array.");
00753     }
00754 
00755         /**
00756          * Constructor from NumpyAnyArray.
00757          * Equivalent to NumpyArray(other.pyObject())
00758          */
00759     explicit NumpyArray(const NumpyAnyArray &other, bool createCopy = false)
00760     {
00761         if(!other.hasData())
00762             return;
00763         if(createCopy)
00764             makeCopy(other.pyObject());
00765         else
00766             vigra_precondition(makeReference(other.pyObject()), //, false),
00767                    "NumpyArray(NumpyAnyArray): Cannot construct from incompatible or empty array.");
00768     }
00769 
00770         /**
00771          * Assignment operator. If this is already a view with data
00772          * (i.e. hasData() is true) and the shapes match, the RHS
00773          * array contents are copied.  If this is an empty view,
00774          * assignment is identical to makeReferenceUnchecked(other.pyObject()).
00775          * See MultiArrayView::operator= for further information on
00776          * semantics.
00777          */
00778     NumpyArray &operator=(const NumpyArray &other)
00779     {
00780         if(hasData())
00781             view_type::operator=(other);
00782         else
00783             makeReferenceUnchecked(other.pyObject());
00784         return *this;
00785     }
00786 
00787         /**
00788          * Assignment operator. If this is already a view with data
00789          * (i.e. hasData() is true) and the shapes match, the RHS
00790          * array contents are copied.  If this is an empty view,
00791          * assignment is identical to makeReferenceUnchecked(other.pyObject()).
00792          * See MultiArrayView::operator= for further information on
00793          * semantics.
00794          */
00795     template <class U, class S>
00796     NumpyArray &operator=(const NumpyArray<N, U, S> &other)
00797     {
00798         if(hasData())
00799         {
00800             vigra_precondition(shape() == other.shape(),
00801                 "NumpyArray::operator=(): shape mismatch.");
00802             view_type::operator=(other);
00803         }
00804         else if(other.hasData())
00805         {
00806             NumpyArray copy;
00807             copy.reshapeIfEmpty(other.taggedShape(),
00808                 "NumpyArray::operator=(): reshape failed unexpectedly.");
00809             copy = other;
00810             makeReferenceUnchecked(copy.pyObject());
00811         }
00812         return *this;
00813     }
00814 
00815         /**
00816          * Assignment operator. If this is already a view with data
00817          * (i.e. hasData() is true) and the shapes match, the RHS
00818          * array contents are copied.  If this is an empty view,
00819          * a new buffer with the RHS shape is allocated before copying.
00820          */
00821     template <class U, class S>
00822     NumpyArray &operator=(const MultiArrayView<N, U, S> &other)
00823     {
00824         if(hasData())
00825         {
00826             vigra_precondition(shape() == other.shape(),
00827                 "NumpyArray::operator=(): shape mismatch.");
00828             view_type::operator=(other);
00829         }
00830         else if(other.hasData())
00831         {
00832             NumpyArray copy;
00833             copy.reshapeIfEmpty(other.shape(),
00834                 "NumpyArray::operator=(): reshape failed unexpectedly.");
00835             copy = other;
00836             makeReferenceUnchecked(copy.pyObject());
00837         }
00838         return *this;
00839     }
00840 
00841         /**
00842          * Assignment operator. If this is already a view with data
00843          * (i.e. hasData() is true) and the shapes match, the RHS
00844          * array contents are copied.
00845          * If this is an empty view, assignment is identical to
00846          * makeReference(other.pyObject()).
00847          * Otherwise, an exception is thrown.
00848          */
00849     NumpyArray &operator=(const NumpyAnyArray &other)
00850     {
00851         if(hasData())
00852         {
00853             NumpyAnyArray::operator=(other);
00854         }
00855         else if(isReferenceCompatible(other.pyObject()))
00856         {
00857             makeReferenceUnchecked(other.pyObject());
00858         }
00859         else
00860         {
00861             vigra_precondition(false,
00862                 "NumpyArray::operator=(): Cannot assign from incompatible array.");
00863         }
00864         return *this;
00865     }
00866 
00867         /**
00868          Permute the entries of the given array \a data exactly like the axes of this NumpyArray
00869          were permuted upon conversion from numpy.
00870          */
00871     template <class U>
00872     ArrayVector<U>
00873     permuteLikewise(ArrayVector<U> const & data) const
00874     {
00875         vigra_precondition(hasData(),
00876             "NumpyArray::permuteLikewise(): array has no data.");
00877 
00878         ArrayVector<U> res(data.size());
00879         ArrayTraits::permuteLikewise(this->pyArray_, data, res);
00880         return res;
00881     }
00882 
00883         /**
00884          Permute the entries of the given array \a data exactly like the axes of this NumpyArray
00885          were permuted upon conversion from numpy.
00886          */
00887     template <class U, int K>
00888     TinyVector<U, K>
00889     permuteLikewise(TinyVector<U, K> const & data) const
00890     {
00891         vigra_precondition(hasData(),
00892             "NumpyArray::permuteLikewise(): array has no data.");
00893 
00894         TinyVector<U, K> res;
00895         ArrayTraits::permuteLikewise(this->pyArray_, data, res);
00896         return res;
00897     }
00898 
00899         /**
00900          Get the permutation of the axes of this NumpyArray
00901          that was performed upon conversion from numpy.
00902          */
00903     template <int K>
00904     TinyVector<npy_intp, K>
00905     permuteLikewise() const
00906     {
00907         vigra_precondition(hasData(),
00908             "NumpyArray::permuteLikewise(): array has no data.");
00909 
00910         TinyVector<npy_intp, K> data, res;
00911         linearSequence(data.begin(), data.end());
00912         ArrayTraits::permuteLikewise(this->pyArray_, data, res);
00913         return res;
00914     }
00915 
00916         /**
00917          * Test whether a given python object is a numpy array that can be
00918          * converted (copied) into an array compatible to this NumpyArray type.
00919          * This means that the array's shape conforms to the requirements of
00920          * makeCopy().
00921          */
00922     static bool isCopyCompatible(PyObject *obj)
00923     {
00924 #if VIGRA_CONVERTER_DEBUG
00925         std::cerr << "class " << typeid(NumpyArray).name() << " got " << obj->ob_type->tp_name << "\n";
00926         std::cerr << "using traits " << typeid(ArrayTraits).name() << "\n";
00927         std::cerr<<"isArray: "<< ArrayTraits::isArray(obj)<<std::endl;
00928         std::cerr<<"isShapeCompatible: "<< ArrayTraits::isShapeCompatible((PyArrayObject *)obj)<<std::endl;
00929 #endif
00930 
00931         return ArrayTraits::isArray(obj) &&
00932                ArrayTraits::isShapeCompatible((PyArrayObject *)obj);
00933     }
00934 
00935         /**
00936          * Test whether a given python object is a numpy array with a
00937          * compatible dtype and the correct shape and strides, so that it
00938          * can be referenced as a view by this NumpyArray type (i.e.
00939          * it conforms to the requirements of makeReference()).
00940          */
00941     static bool isReferenceCompatible(PyObject *obj)
00942     {
00943         return ArrayTraits::isArray(obj) &&
00944                ArrayTraits::isPropertyCompatible((PyArrayObject *)obj);
00945     }
00946 
00947         /**
00948          * Deprecated, use isReferenceCompatible(obj) instead.
00949          */
00950     static bool isStrictlyCompatible(PyObject *obj)
00951     {
00952         return isReferenceCompatible(obj);
00953     }
00954 
00955         /**
00956          * Create a vector representing the standard stride ordering of a NumpyArray.
00957          * That is, we get a vector representing the range [0,...,N-1], which
00958          * denotes the stride ordering for Fortran order.
00959          */
00960     static difference_type standardStrideOrdering()
00961     {
00962         difference_type strideOrdering;
00963         for(unsigned int k=0; k<N; ++k)
00964             strideOrdering[k] = k;
00965         return strideOrdering;
00966     }
00967 
00968         /**
00969          * Set up a view to the given object without checking compatibility.
00970          * This function must not be used unless isReferenceCompatible(obj) returned
00971          * true on the given object (otherwise, a crash is likely).
00972          */
00973     void makeReferenceUnchecked(PyObject *obj)
00974     {
00975         NumpyAnyArray::makeReference(obj);
00976         setupArrayView();
00977     }
00978 
00979         /**
00980          * Try to set up a view referencing the given PyObject.
00981          * Returns false if the python object is not a compatible
00982          * numpy array (see isReferenceCompatible()).
00983          *
00984          * The second parameter ('strict') is deprecated and will be ignored.
00985          */
00986     bool makeReference(PyObject *obj, bool /* strict */ = false)
00987     {
00988         if(!isReferenceCompatible(obj))
00989             return false;
00990         makeReferenceUnchecked(obj);
00991         return true;
00992     }
00993 
00994         /**
00995          * Try to set up a view referencing the same data as the given
00996          * NumpyAnyArray.  This overloaded variant simply calls
00997          * makeReference() on array.pyObject(). The parameter \a strict
00998          * is deprecated and will be ignored.
00999          */
01000     bool makeReference(const NumpyAnyArray &array, bool strict = false)
01001     {
01002         return makeReference(array.pyObject(), strict);
01003     }
01004 
01005         /**
01006          * Set up an unsafe reference to the given MultiArrayView.
01007          * ATTENTION: This creates a numpy.ndarray that points to the
01008          * same data, but does not own it, so it must be ensured by
01009          * other means that the memory does not get freed before the
01010          * end of the ndarray's lifetime!  (One elegant way would be
01011          * to set the 'base' attribute of the resulting ndarray to a
01012          * python object which directly or indirectly holds the memory
01013          * of the given MultiArrayView.)
01014          */
01015     void makeUnsafeReference(const view_type &multiArrayView)
01016     {
01017         vigra_precondition(!hasData(),
01018             "makeUnsafeReference(): cannot replace existing view with given buffer");
01019 
01020         // construct an ndarray that points to our data (taking strides into account):
01021         python_ptr array(ArrayTraits::unsafeConstructorFromData(multiArrayView.shape(),
01022                                   multiArrayView.data(), multiArrayView.stride()));
01023 
01024         view_type::operator=(multiArrayView);
01025         pyArray_ = array;
01026     }
01027 
01028         /**
01029          Try to create a copy of the given PyObject.
01030          Raises an exception when obj is not a compatible array
01031          (see isCopyCompatible() or isReferenceCompatible(), according to the
01032          parameter \a strict) or the Python constructor call failed.
01033          */
01034     void makeCopy(PyObject *obj, bool strict = false)
01035     {
01036 #if VIGRA_CONVERTER_DEBUG
01037         int ndim = PyArray_NDIM((PyArrayObject *)obj);
01038         npy_intp * s = PyArray_DIMS((PyArrayObject *)obj);
01039         std::cerr << "makeCopy: " << ndim << " " <<  ArrayVectorView<npy_intp>(ndim, s) <<
01040                      ", strides " << ArrayVectorView<npy_intp>(ndim, PyArray_STRIDES((PyArrayObject *)obj)) << "\n";
01041         std::cerr << "for " << typeid(*this).name() << "\n";
01042 #endif
01043         vigra_precondition(strict ? isReferenceCompatible(obj) : isCopyCompatible(obj),
01044                      "NumpyArray::makeCopy(obj): Cannot copy an incompatible array.");
01045 
01046         NumpyAnyArray copy(obj, true);
01047         makeReferenceUnchecked(copy.pyObject());
01048     }
01049 
01050         /**
01051             Allocate new memory with the given shape and initialize with zeros.<br>
01052             If a stride ordering is given, the resulting array will have this stride
01053             ordering, when it is compatible with the array's memory layout (unstrided
01054             arrays only permit the standard ascending stride ordering).
01055 
01056             <em>Note:</em> this operation invalidates dependent objects
01057             (MultiArrayViews and iterators)
01058          */
01059     void reshape(difference_type const & shape)
01060     {
01061         vigra_postcondition(makeReference(init(shape)),
01062                 "NumpyArray.reshape(shape): Python constructor did not produce a compatible array.");
01063     }
01064 
01065         /**
01066             When this array has no data, allocate new memory with the given \a shape and
01067             initialize with zeros. Otherwise, check if the new shape matches the old shape
01068             and throw a precondition exception with the given \a message if not.
01069          */
01070     void reshapeIfEmpty(difference_type const & shape, std::string message = "")
01071     {
01072         // FIXME: is this really a good replacement?
01073         // reshapeIfEmpty(shape, standardStrideOrdering(), message);
01074         reshapeIfEmpty(TaggedShape(shape), message);
01075     }
01076 
01077         /**
01078             When this array has no data, allocate new memory with the given \a shape and
01079             initialize with zeros. Otherwise, check if the new shape matches the old shape
01080             and throw a precondition exception with the given \a message if not.
01081          */
01082     void reshapeIfEmpty(TaggedShape tagged_shape, std::string message = "")
01083     {
01084         ArrayTraits::finalizeTaggedShape(tagged_shape);
01085 
01086         if(hasData())
01087         {
01088             vigra_precondition(tagged_shape.compatible(taggedShape()), message.c_str());
01089         }
01090         else
01091         {
01092             python_ptr array(constructArray(tagged_shape, typeCode, true),
01093                              python_ptr::keep_count);
01094             vigra_postcondition(makeReference(NumpyAnyArray(array.get())),
01095                   "NumpyArray.reshapeIfEmpty(): Python constructor did not produce a compatible array.");
01096         }
01097     }
01098 
01099     TaggedShape taggedShape() const
01100     {
01101         return ArrayTraits::taggedShape(this->shape(), PyAxisTags(this->axistags(), true));
01102     }
01103 };
01104 
01105     // this function assumes that pyArray_ has already been set, and compatibility been checked
01106 template <unsigned int N, class T, class Stride>
01107 void NumpyArray<N, T, Stride>::setupArrayView()
01108 {
01109     if(NumpyAnyArray::hasData())
01110     {
01111         permutation_type permute;
01112         ArrayTraits::permutationToSetupOrder(this->pyArray_, permute);
01113 
01114         vigra_precondition(abs((int)permute.size() - actual_dimension) <= 1,
01115             "NumpyArray::setupArrayView(): got array of incompatible shape (should never happen).");
01116 
01117         applyPermutation(permute.begin(), permute.end(),
01118                          pyArray()->dimensions, this->m_shape.begin());
01119         applyPermutation(permute.begin(), permute.end(),
01120                          pyArray()->strides, this->m_stride.begin());
01121 
01122         if((int)permute.size() == actual_dimension - 1)
01123         {
01124             this->m_shape[actual_dimension-1] = 1;
01125             this->m_stride[actual_dimension-1] = sizeof(value_type);
01126         }
01127 
01128         this->m_stride /= sizeof(value_type);
01129         this->m_ptr = reinterpret_cast<pointer>(pyArray()->data);
01130         vigra_precondition(this->checkInnerStride(Stride()),
01131             "NumpyArray<..., UnstridedArrayTag>::setupArrayView(): First dimension of given array is not unstrided (should never happen).");
01132 
01133     }
01134     else
01135     {
01136         this->m_ptr = 0;
01137     }
01138 }
01139 
01140 
01141 typedef NumpyArray<2, float >  NumpyFArray2;
01142 typedef NumpyArray<3, float >  NumpyFArray3;
01143 typedef NumpyArray<4, float >  NumpyFArray4;
01144 typedef NumpyArray<2, Singleband<float> >  NumpyFImage;
01145 typedef NumpyArray<3, Singleband<float> >  NumpyFVolume;
01146 typedef NumpyArray<2, RGBValue<float> >  NumpyFRGBImage;
01147 typedef NumpyArray<3, RGBValue<float> >  NumpyFRGBVolume;
01148 typedef NumpyArray<3, Multiband<float> >  NumpyFMultibandImage;
01149 typedef NumpyArray<4, Multiband<float> >  NumpyFMultibandVolume;
01150 
01151 /********************************************************/
01152 /*                                                      */
01153 /*   NumpyArray Multiband Argument Object Factories     */
01154 /*                                                      */
01155 /********************************************************/
01156 
01157 template <class PixelType, class Stride>
01158 inline triple<ConstStridedImageIterator<PixelType>,
01159               ConstStridedImageIterator<PixelType>,
01160               MultibandVectorAccessor<PixelType> >
01161 srcImageRange(NumpyArray<3, Multiband<PixelType>, Stride> const & img)
01162 {
01163     ConstStridedImageIterator<PixelType>
01164         ul(img.data(), 1, img.stride(0), img.stride(1));
01165     return triple<ConstStridedImageIterator<PixelType>,
01166                   ConstStridedImageIterator<PixelType>,
01167                   MultibandVectorAccessor<PixelType> >
01168         (ul, ul + Size2D(img.shape(0), img.shape(1)), MultibandVectorAccessor<PixelType>(img.shape(2), img.stride(2)));
01169 }
01170 
01171 template <class PixelType, class Stride>
01172 inline pair< ConstStridedImageIterator<PixelType>,
01173              MultibandVectorAccessor<PixelType> >
01174 srcImage(NumpyArray<3, Multiband<PixelType>, Stride> const & img)
01175 {
01176     ConstStridedImageIterator<PixelType>
01177         ul(img.data(), 1, img.stride(0), img.stride(1));
01178     return pair<ConstStridedImageIterator<PixelType>, MultibandVectorAccessor<PixelType> >
01179         (ul, MultibandVectorAccessor<PixelType>(img.shape(2), img.stride(2)));
01180 }
01181 
01182 template <class PixelType, class Stride>
01183 inline triple< StridedImageIterator<PixelType>,
01184                StridedImageIterator<PixelType>,
01185                MultibandVectorAccessor<PixelType> >
01186 destImageRange(NumpyArray<3, Multiband<PixelType>, Stride> & img)
01187 {
01188     StridedImageIterator<PixelType>
01189         ul(img.data(), 1, img.stride(0), img.stride(1));
01190     typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
01191     return triple<StridedImageIterator<PixelType>,
01192                   StridedImageIterator<PixelType>,
01193                   MultibandVectorAccessor<PixelType> >
01194         (ul, ul + Size2D(img.shape(0), img.shape(1)),
01195         MultibandVectorAccessor<PixelType>(img.shape(2), img.stride(2)));
01196 }
01197 
01198 template <class PixelType, class Stride>
01199 inline pair< StridedImageIterator<PixelType>,
01200              MultibandVectorAccessor<PixelType> >
01201 destImage(NumpyArray<3, Multiband<PixelType>, Stride> & img)
01202 {
01203     StridedImageIterator<PixelType>
01204         ul(img.data(), 1, img.stride(0), img.stride(1));
01205     return pair<StridedImageIterator<PixelType>, MultibandVectorAccessor<PixelType> >
01206         (ul, MultibandVectorAccessor<PixelType>(img.shape(2), img.stride(2)));
01207 }
01208 
01209 template <class PixelType, class Stride>
01210 inline pair< ConstStridedImageIterator<PixelType>,
01211              MultibandVectorAccessor<PixelType> >
01212 maskImage(NumpyArray<3, Multiband<PixelType>, Stride> const & img)
01213 {
01214     ConstStridedImageIterator<PixelType>
01215         ul(img.data(), 1, img.stride(0), img.stride(1));
01216     typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
01217     return pair<ConstStridedImageIterator<PixelType>, MultibandVectorAccessor<PixelType> >
01218         (ul, MultibandVectorAccessor<PixelType>(img.shape(2), img.stride(2)));
01219 }
01220 
01221 } // namespace vigra
01222 
01223 #endif // VIGRA_NUMPY_ARRAY_HXX

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

html generated using doxygen and Python
vigra 1.9.0 (Tue Nov 6 2012)