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