[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/numpy_array.hxx | ![]() |
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) |
html generated using doxygen and Python
|