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

numpy_array.hxx
1 /************************************************************************/
2 /* */
3 /* Copyright 2009 by Ullrich Koethe and Hans Meine */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_NUMPY_ARRAY_HXX
37 #define VIGRA_NUMPY_ARRAY_HXX
38 
39 #include <Python.h>
40 #include <string>
41 #include <iostream>
42 #include <numpy/arrayobject.h>
43 #include "multi_array.hxx"
44 #include "array_vector.hxx"
45 #include "python_utility.hxx"
46 #include "numpy_array_traits.hxx"
47 #include "numpy_array_taggedshape.hxx"
48 
49 // NumPy function called by NumPy’s import_array() macro (and our import_vigranumpy() below)
50 int _import_array();
51 
52 namespace vigra {
53 
54 static inline void import_vigranumpy()
55 {
56  // roughly equivalent to import_array():
57  if(_import_array() < 0)
58  pythonToCppException(0);
59 
60  // in addition, import vigra.vigranumpycore:
61  python_ptr module(PyImport_ImportModule("vigra.vigranumpycore"), python_ptr::keep_count);
62  pythonToCppException(module);
63 }
64 
65 /********************************************************/
66 /* */
67 /* MultibandVectorAccessor */
68 /* */
69 /********************************************************/
70 
71 template <class T>
72 class MultibandVectorAccessor
73 {
74  MultiArrayIndex size_, stride_;
75 
76  public:
77  MultibandVectorAccessor(MultiArrayIndex size, MultiArrayIndex stride)
78  : size_(size),
79  stride_(stride)
80  {}
81 
82 
83  typedef Multiband<T> value_type;
84 
85  /** the vector's value_type
86  */
87  typedef T component_type;
88 
89  typedef VectorElementAccessor<MultibandVectorAccessor<T> > ElementAccessor;
90 
91  /** Read the component data at given vector index
92  at given iterator position
93  */
94  template <class ITERATOR>
95  component_type const & getComponent(ITERATOR const & i, int idx) const
96  {
97  return *(&*i+idx*stride_);
98  }
99 
100  /** Set the component data at given vector index
101  at given iterator position. The type <TT>V</TT> of the passed
102  in <TT>value</TT> is automatically converted to <TT>component_type</TT>.
103  In case of a conversion floating point -> integral this includes rounding and clipping.
104  */
105  template <class V, class ITERATOR>
106  void setComponent(V const & value, ITERATOR const & i, int idx) const
107  {
108  *(&*i+idx*stride_) = detail::RequiresExplicitCast<component_type>::cast(value);
109  }
110 
111  /** Read the component data at given vector index
112  at an offset of given iterator position
113  */
114  template <class ITERATOR, class DIFFERENCE>
115  component_type const & getComponent(ITERATOR const & i, DIFFERENCE const & diff, int idx) const
116  {
117  return *(&i[diff]+idx*stride_);
118  }
119 
120  /** Set the component data at given vector index
121  at an offset of given iterator position. The type <TT>V</TT> of the passed
122  in <TT>value</TT> is automatically converted to <TT>component_type</TT>.
123  In case of a conversion floating point -> integral this includes rounding and clipping.
124  */
125  template <class V, class ITERATOR, class DIFFERENCE>
126  void
127  setComponent(V const & value, ITERATOR const & i, DIFFERENCE const & diff, int idx) const
128  {
129  *(&i[diff]+idx*stride_) = detail::RequiresExplicitCast<component_type>::cast(value);
130  }
131 
132  template <class U>
133  MultiArrayIndex size(U) const
134  {
135  return size_;
136  }
137 };
138 
139 /********************************************************/
140 
141 template <class TYPECODE> // pseudo-template to avoid inline expansion of the function
142  // will always be NPY_TYPES
143 PyObject *
144 constructArray(TaggedShape tagged_shape, TYPECODE typeCode, bool init,
145  python_ptr arraytype = python_ptr());
146 
147 /********************************************************/
148 /* */
149 /* NumpyAnyArray */
150 /* */
151 /********************************************************/
152 
153 /** Wrapper class for a Python array.
154 
155  This class stores a reference-counted pointer to an Python numpy array object,
156  i.e. an object where <tt>PyArray_Check(object)</tt> returns true (in Python, the
157  object is then a subclass of <tt>numpy.ndarray</tt>). This class is mainly used
158  as a smart pointer to these arrays, but some basic access and conversion functions
159  are also provided.
160 
161  <b>\#include</b> <vigra/numpy_array.hxx><br>
162  Namespace: vigra
163 */
165 {
166  protected:
167  python_ptr pyArray_;
168 
169  public:
170 
171  /// difference type
173 
174  static python_ptr getArrayTypeObject()
175  {
176  return detail::getArrayTypeObject();
177  }
178 
179  static std::string defaultOrder(std::string defaultValue = "C")
180  {
181  return detail::defaultOrder(defaultValue);
182  }
183 
184  static python_ptr defaultAxistags(int ndim, std::string order = "")
185  {
186  return detail::defaultAxistags(ndim, order);
187  }
188 
189  static python_ptr emptyAxistags(int ndim)
190  {
191  return detail::emptyAxistags(ndim);
192  }
193 
194  /**
195  Construct from a Python object. If \a obj is NULL, or is not a subclass
196  of numpy.ndarray, the resulting NumpyAnyArray will have no data (i.e.
197  hasData() returns false). Otherwise, it creates a new reference to the array
198  \a obj, unless \a createCopy is true, where a new array is created by calling
199  the C-equivalent of obj->copy().
200  */
201  explicit NumpyAnyArray(PyObject * obj = 0, bool createCopy = false, PyTypeObject * type = 0)
202  {
203  if(obj == 0)
204  return;
205  vigra_precondition(type == 0 || PyType_IsSubtype(type, &PyArray_Type),
206  "NumpyAnyArray(obj, createCopy, type): type must be numpy.ndarray or a subclass thereof.");
207  if(createCopy)
208  makeCopy(obj, type);
209  else
210  vigra_precondition(makeReference(obj, type), "NumpyAnyArray(obj): obj isn't a numpy array.");
211  }
212 
213  /**
214  Copy constructor. By default, it creates a new reference to the array
215  \a other. When \a createCopy is true, a new array is created by calling
216  the C-equivalent of other.copy().
217  */
218  NumpyAnyArray(NumpyAnyArray const & other, bool createCopy = false, PyTypeObject * type = 0)
219  {
220  if(!other.hasData())
221  return;
222  vigra_precondition(type == 0 || PyType_IsSubtype(type, &PyArray_Type),
223  "NumpyAnyArray(obj, createCopy, type): type must be numpy.ndarray or a subclass thereof.");
224  if(createCopy)
225  makeCopy(other.pyObject(), type);
226  else
227  makeReference(other.pyObject(), type);
228  }
229 
230  // auto-generated destructor is ok
231 
232  /**
233  * Assignment operator. If this is already a view with data
234  * (i.e. hasData() is true) and the shapes match, the RHS
235  * array contents are copied via the C-equivalent of
236  * 'self[...] = other[...]'. If the shapes don't matched,
237  * broadcasting is tried on the trailing (i.e. channel)
238  * dimension.
239  * If the LHS is an empty view, assignment is identical to
240  * makeReference(other.pyObject()).
241  */
243  {
244  if(hasData())
245  {
246  vigra_precondition(other.hasData(),
247  "NumpyArray::operator=(): Cannot assign from empty array.");
248 
249  python_ptr arraytype = getArrayTypeObject();
250  python_ptr f(PyString_FromString("_copyValuesImpl"), python_ptr::keep_count);
251  if(PyObject_HasAttr(arraytype, f))
252  {
253  python_ptr res(PyObject_CallMethodObjArgs(arraytype, f.get(),
254  pyArray_.get(), other.pyArray_.get(), NULL),
255  python_ptr::keep_count);
256  vigra_postcondition(res.get() != 0,
257  "NumpyArray::operator=(): VigraArray._copyValuesImpl() failed.");
258  }
259  else
260  {
261  PyArrayObject * sarray = (PyArrayObject *)pyArray_.get();
262  PyArrayObject * tarray = (PyArrayObject *)other.pyArray_.get();
263 
264  if(PyArray_CopyInto(tarray, sarray) == -1)
265  pythonToCppException(0);
266  }
267  }
268  else
269  {
270  pyArray_ = other.pyArray_;
271  }
272  return *this;
273  }
274 
275  /**
276  Returns the number of dimensions of this array, or 0 if
277  hasData() is false.
278  */
280  {
281  if(hasData())
282  return PyArray_NDIM(pyObject());
283  return 0;
284  }
285 
286  /**
287  Returns the number of spatial dimensions of this array, or 0 if
288  hasData() is false. If the enclosed Python array does not define
289  the attribute spatialDimensions, ndim() is returned.
290  */
292  {
293  if(!hasData())
294  return 0;
295  return pythonGetAttr(pyObject(), "spatialDimensions", ndim());
296  }
297 
298  bool hasChannelAxis() const
299  {
300  if(!hasData())
301  return false;
302  return channelIndex() == ndim();
303  }
304 
305  MultiArrayIndex channelIndex() const
306  {
307  if(!hasData())
308  return 0;
309  return pythonGetAttr(pyObject(), "channelIndex", ndim());
310  }
311 
312  MultiArrayIndex innerNonchannelIndex() const
313  {
314  if(!hasData())
315  return 0;
316  return pythonGetAttr(pyObject(), "innerNonchannelIndex", ndim());
317  }
318 
319  /**
320  Returns the shape of this array. The size of
321  the returned shape equals ndim().
322  */
324  {
325  if(hasData())
326  return difference_type(PyArray_DIMS(pyObject()), PyArray_DIMS(pyObject()) + ndim());
327  return difference_type();
328  }
329 
330  /** Compute the ordering of the strides of this array.
331  The result is describes the current permutation of the axes relative
332  to an ascending stride order.
333  */
335  {
336  if(!hasData())
337  return difference_type();
338  MultiArrayIndex N = ndim();
339  difference_type stride(PyArray_STRIDES(pyObject()), PyArray_STRIDES(pyObject()) + N),
340  permutation(N);
341  for(MultiArrayIndex k=0; k<N; ++k)
342  permutation[k] = k;
343  for(MultiArrayIndex k=0; k<N-1; ++k)
344  {
345  MultiArrayIndex smallest = k;
346  for(MultiArrayIndex j=k+1; j<N; ++j)
347  {
348  if(stride[j] < stride[smallest])
349  smallest = j;
350  }
351  if(smallest != k)
352  {
353  std::swap(stride[k], stride[smallest]);
354  std::swap(permutation[k], permutation[smallest]);
355  }
356  }
357  difference_type ordering(N);
358  for(MultiArrayIndex k=0; k<N; ++k)
359  ordering[permutation[k]] = k;
360  return ordering;
361  }
362 
363  // /**
364  // Returns the the permutation that will transpose this array into
365  // canonical ordering (currently: F-order). The size of
366  // the returned permutation equals ndim().
367  // */
368  // difference_type permutationToNormalOrder() const
369  // {
370  // if(!hasData())
371  // return difference_type();
372 
373  // // difference_type res(detail::getAxisPermutationImpl(pyArray_,
374  // // "permutationToNormalOrder", true));
375  // difference_type res;
376  // detail::getAxisPermutationImpl(res, pyArray_, "permutationToNormalOrder", true);
377  // if(res.size() == 0)
378  // {
379  // res.resize(ndim());
380  // linearSequence(res.begin(), res.end(), ndim()-1, MultiArrayIndex(-1));
381  // }
382  // return res;
383  // }
384 
385  /**
386  Returns the value type of the elements in this array, or -1
387  when hasData() is false.
388  */
389  int dtype() const
390  {
391  if(hasData())
392  return PyArray_DESCR(pyObject())->type_num;
393  return -1;
394  }
395 
396  /**
397  * Return the AxisTags of this array or a NULL pointer when the attribute
398  'axistags' is missing in the Python object or this array has no data.
399  */
400  python_ptr axistags() const
401  {
402  static python_ptr key(PyString_FromString("axistags"), python_ptr::keep_count);
403 
404  python_ptr axistags;
405  if(pyObject())
406  {
407  axistags.reset(PyObject_GetAttr(pyObject(), key), python_ptr::keep_count);
408  if(!axistags)
409  PyErr_Clear();
410  }
411  return axistags;
412  }
413 
414  /**
415  * Return a borrowed reference to the internal PyArrayObject.
416  */
417  PyArrayObject * pyArray() const
418  {
419  return (PyArrayObject *)pyArray_.get();
420  }
421 
422  /**
423  * Return a borrowed reference to the internal PyArrayObject
424  * (see pyArray()), cast to PyObject for your convenience.
425  */
426  PyObject * pyObject() const
427  {
428  return pyArray_.get();
429  }
430 
431  /**
432  Reset the NumpyAnyArray to the given object. If \a obj is a numpy array object,
433  a new reference to that array is created, and the function returns
434  true. Otherwise, it returns false and the NumpyAnyArray remains unchanged.
435  If \a type is given, the new reference will be a view with that type, provided
436  that \a type is a numpy ndarray or a subclass thereof. Otherwise, an
437  exception is thrown.
438  */
439  bool makeReference(PyObject * obj, PyTypeObject * type = 0)
440  {
441  if(obj == 0 || !PyArray_Check(obj))
442  return false;
443  if(type != 0)
444  {
445  vigra_precondition(PyType_IsSubtype(type, &PyArray_Type) != 0,
446  "NumpyAnyArray::makeReference(obj, type): type must be numpy.ndarray or a subclass thereof.");
447  obj = PyArray_View((PyArrayObject*)obj, 0, type);
448  pythonToCppException(obj);
449  }
450  pyArray_.reset(obj);
451  return true;
452  }
453 
454  /**
455  Create a copy of the given array object. If \a obj is a numpy array object,
456  a copy is created via the C-equivalent of 'obj->copy()'. If
457  this call fails, or obj was not an array, an exception is thrown
458  and the NumpyAnyArray remains unchanged.
459  */
460  void makeCopy(PyObject * obj, PyTypeObject * type = 0)
461  {
462  vigra_precondition(obj && PyArray_Check(obj),
463  "NumpyAnyArray::makeCopy(obj): obj is not an array.");
464  vigra_precondition(type == 0 || PyType_IsSubtype(type, &PyArray_Type),
465  "NumpyAnyArray::makeCopy(obj, type): type must be numpy.ndarray or a subclass thereof.");
466  python_ptr array(PyArray_NewCopy((PyArrayObject*)obj, NPY_ANYORDER), python_ptr::keep_count);
467  pythonToCppException(array);
468  makeReference(array, type);
469  }
470 
471  /**
472  Check whether this NumpyAnyArray actually points to a Python array.
473  */
474  bool hasData() const
475  {
476  return pyArray_ != 0;
477  }
478 };
479 
480 /********************************************************/
481 /* */
482 /* constructArray */
483 /* */
484 /********************************************************/
485 
486 namespace detail {
487 
488 inline bool
489 nontrivialPermutation(ArrayVector<npy_intp> const & p)
490 {
491  for(unsigned int k=0; k<p.size(); ++k)
492  if(p[k] != k)
493  return true;
494  return false;
495 }
496 
497 } // namespace detail
498 
499 template <class TYPECODE> // pseudo-template to avoid inline expansion of the function
500  // will always be NPY_TYPES
501 PyObject *
502 constructArray(TaggedShape tagged_shape, TYPECODE typeCode, bool init, python_ptr arraytype)
503 {
504  ArrayVector<npy_intp> shape = finalizeTaggedShape(tagged_shape);
505  PyAxisTags axistags(tagged_shape.axistags);
506 
507  int ndim = (int)shape.size();
508  ArrayVector<npy_intp> inverse_permutation;
509  int order = 1; // Fortran order
510 
511  if(axistags)
512  {
513  if(!arraytype)
514  arraytype = NumpyAnyArray::getArrayTypeObject();
515 
516  inverse_permutation = axistags.permutationFromNormalOrder();
517  vigra_precondition(ndim == (int)inverse_permutation.size(),
518  "axistags.permutationFromNormalOrder(): permutation has wrong size.");
519  }
520  else
521  {
522  arraytype = python_ptr((PyObject*)&PyArray_Type);
523  order = 0; // C order
524  }
525 
526 // std::cerr << "constructArray: " << shape << "\n" << inverse_permutation << "\n";
527 
528  python_ptr array(PyArray_New((PyTypeObject *)arraytype.get(), ndim, shape.begin(),
529  typeCode, 0, 0, 0, order, 0),
530  python_ptr::keep_count);
531  pythonToCppException(array);
532 
533  if(detail::nontrivialPermutation(inverse_permutation))
534  {
535  PyArray_Dims permute = { inverse_permutation.begin(), ndim };
536  array = python_ptr(PyArray_Transpose((PyArrayObject*)array.get(), &permute),
537  python_ptr::keep_count);
538  pythonToCppException(array);
539  }
540 
541  if(arraytype != (PyObject*)&PyArray_Type && axistags)
542  pythonToCppException(PyObject_SetAttrString(array, "axistags", axistags.axistags) != -1);
543 
544  if(init)
545  PyArray_FILLWBYTE((PyArrayObject *)array.get(), 0);
546 
547  return array.release();
548 }
549 
550 // FIXME: reimplement in terms of TaggedShape?
551 template <class TINY_VECTOR>
552 inline
553 python_ptr constructNumpyArrayFromData(
554  TINY_VECTOR const & shape, npy_intp *strides,
555  NPY_TYPES typeCode, void *data)
556 {
557  ArrayVector<npy_intp> pyShape(shape.begin(), shape.end());
558 
559  python_ptr array(PyArray_New(&PyArray_Type, shape.size(), pyShape.begin(),
560  typeCode, strides, data, 0, NPY_WRITEABLE, 0),
561  python_ptr::keep_count);
562  pythonToCppException(array);
563 
564  return array;
565 }
566 
567 /********************************************************/
568 /* */
569 /* NumpyArray */
570 /* */
571 /********************************************************/
572 
573 /** Provide the MultiArrayView interface for a Python array.
574 
575  This class inherits from both \ref vigra::MultiArrayView and \ref vigra::NumpyAnyArray
576  in order to support easy and save application of VIGRA functions to Python arrays.
577 
578  <b>\#include</b> <vigra/numpy_array.hxx><br>
579  Namespace: vigra
580 */
581 template <unsigned int N, class T, class Stride = StridedArrayTag>
583 : public MultiArrayView<N, typename NumpyArrayTraits<N, T, Stride>::value_type, Stride>,
584  public NumpyAnyArray
585 {
586  public:
587  typedef NumpyArrayTraits<N, T, Stride> ArrayTraits;
588  typedef typename ArrayTraits::dtype dtype;
589  typedef T pseudo_value_type;
590 
591  static NPY_TYPES const typeCode = ArrayTraits::typeCode;
592 
593  /** the view type associated with this array.
594  */
596 
597  enum { actual_dimension = view_type::actual_dimension };
598 
599  /** the array's value type
600  */
602 
603  /** pointer type
604  */
605  typedef typename view_type::pointer pointer;
606 
607  /** const pointer type
608  */
610 
611  /** reference type (result of operator[])
612  */
613  typedef typename view_type::reference reference;
614 
615  /** const reference type (result of operator[] const)
616  */
618 
619  /** size type
620  */
621  typedef typename view_type::size_type size_type;
622 
623  /** difference type (used for multi-dimensional offsets and indices)
624  */
626 
627  /** difference and index type for a single dimension
628  */
630 
631  /** type of an array specifying an axis permutation
632  */
634 
635  /** traverser type
636  */
637  typedef typename view_type::traverser traverser;
638 
639  /** traverser type to const data
640  */
642 
643  /** sequential (random access) iterator type
644  */
645  typedef typename view_type::iterator iterator;
646 
647  /** sequential (random access) const iterator type
648  */
650 
651  using view_type::shape; // resolve ambiguity of multiple inheritance
652  using view_type::hasData; // resolve ambiguity of multiple inheritance
653  using view_type::strideOrdering; // resolve ambiguity of multiple inheritance
654 
655  protected:
656 
657  // this function assumes that pyArray_ has already been set, and compatibility been checked
658  void setupArrayView();
659 
660  static python_ptr init(difference_type const & shape, bool init = true,
661  std::string const & order = "")
662  {
663  vigra_precondition(order == "" || order == "C" || order == "F" ||
664  order == "V" || order == "A",
665  "NumpyArray.init(): order must be in ['C', 'F', 'V', 'A', ''].");
666  return python_ptr(constructArray(ArrayTraits::taggedShape(shape, order), typeCode, init),
667  python_ptr::keep_count);
668  }
669 
670  public:
671 
672  using view_type::init;
673 
674  /**
675  * Construct from a given PyObject pointer. When the given
676  * python object is NULL, the internal python array will be
677  * NULL and hasData() will return false.
678  *
679  * Otherwise, the function attempts to create a
680  * new reference to the given Python object, unless
681  * copying is forced by setting \a createCopy to true.
682  * If either of this fails, the function throws an exception.
683  * This will not happen if isReferenceCompatible(obj) (in case
684  * of creating a new reference) or isCopyCompatible(obj)
685  * (in case of copying) have returned true beforehand.
686  */
687  explicit NumpyArray(PyObject *obj = 0, bool createCopy = false)
688  {
689  if(obj == 0)
690  return;
691  if(createCopy)
692  makeCopy(obj);
693  else
694  vigra_precondition(makeReference(obj),
695  "NumpyArray(obj): Cannot construct from incompatible array.");
696  }
697 
698  /**
699  * Copy constructor; does not copy the memory, but creates a
700  * new reference to the same underlying python object, unless
701  * a copy is forced by setting \a createCopy to true.
702  * (If the source object has no data, this one will have
703  * no data, too.)
704  */
705  NumpyArray(const NumpyArray &other, bool createCopy = false)
706  : view_type(),
707  NumpyAnyArray()
708  {
709  if(!other.hasData())
710  return;
711  if(createCopy)
712  makeCopy(other.pyObject());
713  else
715  }
716 
717  /**
718  * Allocate new memory and copy data from a MultiArrayView.
719  */
720  template <class U, class S>
721  explicit NumpyArray(const MultiArrayView<N, U, S> &other)
722  {
723  if(!other.hasData())
724  return;
725  vigra_postcondition(makeReference(init(other.shape(), false)),
726  "NumpyArray(MultiArrayView): Python constructor did not produce a compatible array.");
727  view_type::operator=(other);
728  }
729 
730  /**
731  * Construct a new array object, allocating an internal python
732  * ndarray of the given shape in the given order (default: VIGRA order), initialized
733  * with zeros.
734  *
735  * An exception is thrown when construction fails.
736  */
737  explicit NumpyArray(difference_type const & shape, std::string const & order = "")
738  {
739  vigra_postcondition(makeReference(init(shape, true, order)),
740  "NumpyArray(shape): Python constructor did not produce a compatible array.");
741  }
742 
743  /**
744  * Construct a new array object, allocating an internal python
745  * ndarray according to the given tagged shape, initialized with zeros.
746  *
747  * An exception is thrown when construction fails.
748  */
749  explicit NumpyArray(TaggedShape const & tagged_shape)
750  {
751  reshapeIfEmpty(tagged_shape,
752  "NumpyArray(tagged_shape): Python constructor did not produce a compatible array.");
753  }
754 
755  /**
756  * Constructor from NumpyAnyArray.
757  * Equivalent to NumpyArray(other.pyObject())
758  */
759  explicit NumpyArray(const NumpyAnyArray &other, bool createCopy = false)
760  {
761  if(!other.hasData())
762  return;
763  if(createCopy)
764  makeCopy(other.pyObject());
765  else
766  vigra_precondition(makeReference(other.pyObject()), //, false),
767  "NumpyArray(NumpyAnyArray): Cannot construct from incompatible or empty array.");
768  }
769 
770  /**
771  * Assignment operator. If this is already a view with data
772  * (i.e. hasData() is true) and the shapes match, the RHS
773  * array contents are copied. If this is an empty view,
774  * assignment is identical to makeReferenceUnchecked(other.pyObject()).
775  * See MultiArrayView::operator= for further information on
776  * semantics.
777  */
779  {
780  if(hasData())
781  view_type::operator=(other);
782  else
784  return *this;
785  }
786 
787  /**
788  * Assignment operator. If this is already a view with data
789  * (i.e. hasData() is true) and the shapes match, the RHS
790  * array contents are copied. If this is an empty view,
791  * assignment is identical to makeReferenceUnchecked(other.pyObject()).
792  * See MultiArrayView::operator= for further information on
793  * semantics.
794  */
795  template <class U, class S>
797  {
798  if(hasData())
799  {
800  vigra_precondition(shape() == other.shape(),
801  "NumpyArray::operator=(): shape mismatch.");
802  view_type::operator=(other);
803  }
804  else if(other.hasData())
805  {
807  copy.reshapeIfEmpty(other.taggedShape(),
808  "NumpyArray::operator=(): reshape failed unexpectedly.");
809  copy = other;
811  }
812  return *this;
813  }
814 
815  /**
816  * Assignment operator. If this is already a view with data
817  * (i.e. hasData() is true) and the shapes match, the RHS
818  * array contents are copied. If this is an empty view,
819  * a new buffer with the RHS shape is allocated before copying.
820  */
821  template <class U, class S>
823  {
824  if(hasData())
825  {
826  vigra_precondition(shape() == other.shape(),
827  "NumpyArray::operator=(): shape mismatch.");
828  view_type::operator=(other);
829  }
830  else if(other.hasData())
831  {
833  copy.reshapeIfEmpty(other.shape(),
834  "NumpyArray::operator=(): reshape failed unexpectedly.");
835  copy = other;
837  }
838  return *this;
839  }
840 
841  /**
842  * Assignment operator. If this is already a view with data
843  * (i.e. hasData() is true) and the shapes match, the RHS
844  * array contents are copied.
845  * If this is an empty view, assignment is identical to
846  * makeReference(other.pyObject()).
847  * Otherwise, an exception is thrown.
848  */
850  {
851  if(hasData())
852  {
854  }
855  else if(isReferenceCompatible(other.pyObject()))
856  {
858  }
859  else
860  {
861  vigra_precondition(false,
862  "NumpyArray::operator=(): Cannot assign from incompatible array.");
863  }
864  return *this;
865  }
866 
867  /**
868  Permute the entries of the given array \a data exactly like the axes of this NumpyArray
869  were permuted upon conversion from numpy.
870  */
871  template <class U>
873  permuteLikewise(ArrayVector<U> const & data) const
874  {
875  vigra_precondition(hasData(),
876  "NumpyArray::permuteLikewise(): array has no data.");
877 
878  ArrayVector<U> res(data.size());
879  ArrayTraits::permuteLikewise(this->pyArray_, data, res);
880  return res;
881  }
882 
883  /**
884  Permute the entries of the given array \a data exactly like the axes of this NumpyArray
885  were permuted upon conversion from numpy.
886  */
887  template <class U, int K>
889  permuteLikewise(TinyVector<U, K> const & data) const
890  {
891  vigra_precondition(hasData(),
892  "NumpyArray::permuteLikewise(): array has no data.");
893 
894  TinyVector<U, K> res;
895  ArrayTraits::permuteLikewise(this->pyArray_, data, res);
896  return res;
897  }
898 
899  /**
900  Get the permutation of the axes of this NumpyArray
901  that was performed upon conversion from numpy.
902  */
903  template <int K>
906  {
907  vigra_precondition(hasData(),
908  "NumpyArray::permuteLikewise(): array has no data.");
909 
911  linearSequence(data.begin(), data.end());
912  ArrayTraits::permuteLikewise(this->pyArray_, data, res);
913  return res;
914  }
915 
916  /**
917  * Test whether a given python object is a numpy array that can be
918  * converted (copied) into an array compatible to this NumpyArray type.
919  * This means that the array's shape conforms to the requirements of
920  * makeCopy().
921  */
922  static bool isCopyCompatible(PyObject *obj)
923  {
924 #if VIGRA_CONVERTER_DEBUG
925  std::cerr << "class " << typeid(NumpyArray).name() << " got " << obj->ob_type->tp_name << "\n";
926  std::cerr << "using traits " << typeid(ArrayTraits).name() << "\n";
927  std::cerr<<"isArray: "<< ArrayTraits::isArray(obj)<<std::endl;
928  std::cerr<<"isShapeCompatible: "<< ArrayTraits::isShapeCompatible((PyArrayObject *)obj)<<std::endl;
929 #endif
930 
931  return ArrayTraits::isArray(obj) &&
932  ArrayTraits::isShapeCompatible((PyArrayObject *)obj);
933  }
934 
935  /**
936  * Test whether a given python object is a numpy array with a
937  * compatible dtype and the correct shape and strides, so that it
938  * can be referenced as a view by this NumpyArray type (i.e.
939  * it conforms to the requirements of makeReference()).
940  */
941  static bool isReferenceCompatible(PyObject *obj)
942  {
943  return ArrayTraits::isArray(obj) &&
944  ArrayTraits::isPropertyCompatible((PyArrayObject *)obj);
945  }
946 
947  /**
948  * Deprecated, use isReferenceCompatible(obj) instead.
949  */
950  static bool isStrictlyCompatible(PyObject *obj)
951  {
952  return isReferenceCompatible(obj);
953  }
954 
955  /**
956  * Create a vector representing the standard stride ordering of a NumpyArray.
957  * That is, we get a vector representing the range [0,...,N-1], which
958  * denotes the stride ordering for Fortran order.
959  */
961  {
963  for(unsigned int k=0; k<N; ++k)
964  strideOrdering[k] = k;
965  return strideOrdering;
966  }
967 
968  /**
969  * Set up a view to the given object without checking compatibility.
970  * This function must not be used unless isReferenceCompatible(obj) returned
971  * true on the given object (otherwise, a crash is likely).
972  */
973  void makeReferenceUnchecked(PyObject *obj)
974  {
976  setupArrayView();
977  }
978 
979  /**
980  * Try to set up a view referencing the given PyObject.
981  * Returns false if the python object is not a compatible
982  * numpy array (see isReferenceCompatible()).
983  *
984  * The second parameter ('strict') is deprecated and will be ignored.
985  */
986  bool makeReference(PyObject *obj, bool /* strict */ = false)
987  {
988  if(!isReferenceCompatible(obj))
989  return false;
991  return true;
992  }
993 
994  /**
995  * Try to set up a view referencing the same data as the given
996  * NumpyAnyArray. This overloaded variant simply calls
997  * makeReference() on array.pyObject(). The parameter \a strict
998  * is deprecated and will be ignored.
999  */
1000  bool makeReference(const NumpyAnyArray &array, bool strict = false)
1001  {
1002  return makeReference(array.pyObject(), strict);
1003  }
1004 
1005  /**
1006  * Set up an unsafe reference to the given MultiArrayView.
1007  * ATTENTION: This creates a numpy.ndarray that points to the
1008  * same data, but does not own it, so it must be ensured by
1009  * other means that the memory does not get freed before the
1010  * end of the ndarray's lifetime! (One elegant way would be
1011  * to set the 'base' attribute of the resulting ndarray to a
1012  * python object which directly or indirectly holds the memory
1013  * of the given MultiArrayView.)
1014  */
1015  void makeUnsafeReference(const view_type &multiArrayView)
1016  {
1017  vigra_precondition(!hasData(),
1018  "makeUnsafeReference(): cannot replace existing view with given buffer");
1019 
1020  // construct an ndarray that points to our data (taking strides into account):
1021  python_ptr array(ArrayTraits::unsafeConstructorFromData(multiArrayView.shape(),
1022  multiArrayView.data(), multiArrayView.stride()));
1023 
1024  view_type::operator=(multiArrayView);
1025  pyArray_ = array;
1026  }
1027 
1028  /**
1029  Try to create a copy of the given PyObject.
1030  Raises an exception when obj is not a compatible array
1031  (see isCopyCompatible() or isReferenceCompatible(), according to the
1032  parameter \a strict) or the Python constructor call failed.
1033  */
1034  void makeCopy(PyObject *obj, bool strict = false)
1035  {
1036 #if VIGRA_CONVERTER_DEBUG
1037  int ndim = PyArray_NDIM((PyArrayObject *)obj);
1038  npy_intp * s = PyArray_DIMS((PyArrayObject *)obj);
1039  std::cerr << "makeCopy: " << ndim << " " << ArrayVectorView<npy_intp>(ndim, s) <<
1040  ", strides " << ArrayVectorView<npy_intp>(ndim, PyArray_STRIDES((PyArrayObject *)obj)) << "\n";
1041  std::cerr << "for " << typeid(*this).name() << "\n";
1042 #endif
1043  vigra_precondition(strict ? isReferenceCompatible(obj) : isCopyCompatible(obj),
1044  "NumpyArray::makeCopy(obj): Cannot copy an incompatible array.");
1045 
1046  NumpyAnyArray copy(obj, true);
1048  }
1049 
1050  /**
1051  Allocate new memory with the given shape and initialize with zeros.<br>
1052  If a stride ordering is given, the resulting array will have this stride
1053  ordering, when it is compatible with the array's memory layout (unstrided
1054  arrays only permit the standard ascending stride ordering).
1055 
1056  <em>Note:</em> this operation invalidates dependent objects
1057  (MultiArrayViews and iterators)
1058  */
1059  void reshape(difference_type const & shape)
1060  {
1061  vigra_postcondition(makeReference(init(shape)),
1062  "NumpyArray.reshape(shape): Python constructor did not produce a compatible array.");
1063  }
1064 
1065  /**
1066  When this array has no data, allocate new memory with the given \a shape and
1067  initialize with zeros. Otherwise, check if the new shape matches the old shape
1068  and throw a precondition exception with the given \a message if not.
1069  */
1070  void reshapeIfEmpty(difference_type const & shape, std::string message = "")
1071  {
1072  // FIXME: is this really a good replacement?
1073  // reshapeIfEmpty(shape, standardStrideOrdering(), message);
1074  reshapeIfEmpty(TaggedShape(shape), message);
1075  }
1076 
1077  /**
1078  When this array has no data, allocate new memory with the given \a shape and
1079  initialize with zeros. Otherwise, check if the new shape matches the old shape
1080  and throw a precondition exception with the given \a message if not.
1081  */
1082  void reshapeIfEmpty(TaggedShape tagged_shape, std::string message = "")
1083  {
1084  ArrayTraits::finalizeTaggedShape(tagged_shape);
1085 
1086  if(hasData())
1087  {
1088  vigra_precondition(tagged_shape.compatible(taggedShape()), message.c_str());
1089  }
1090  else
1091  {
1092  python_ptr array(constructArray(tagged_shape, typeCode, true),
1093  python_ptr::keep_count);
1094  vigra_postcondition(makeReference(NumpyAnyArray(array.get())),
1095  "NumpyArray.reshapeIfEmpty(): Python constructor did not produce a compatible array.");
1096  }
1097  }
1098 
1099  TaggedShape taggedShape() const
1100  {
1101  return ArrayTraits::taggedShape(this->shape(), PyAxisTags(this->axistags(), true));
1102  }
1103 };
1104 
1105  // this function assumes that pyArray_ has already been set, and compatibility been checked
1106 template <unsigned int N, class T, class Stride>
1107 void NumpyArray<N, T, Stride>::setupArrayView()
1108 {
1110  {
1111  permutation_type permute;
1112  ArrayTraits::permutationToSetupOrder(this->pyArray_, permute);
1113 
1114  vigra_precondition(abs((int)permute.size() - actual_dimension) <= 1,
1115  "NumpyArray::setupArrayView(): got array of incompatible shape (should never happen).");
1116 
1117  applyPermutation(permute.begin(), permute.end(),
1118  pyArray()->dimensions, this->m_shape.begin());
1119  applyPermutation(permute.begin(), permute.end(),
1120  pyArray()->strides, this->m_stride.begin());
1121 
1122  if((int)permute.size() == actual_dimension - 1)
1123  {
1124  this->m_shape[actual_dimension-1] = 1;
1125  this->m_stride[actual_dimension-1] = sizeof(value_type);
1126  }
1127 
1128  this->m_stride /= sizeof(value_type);
1129  this->m_ptr = reinterpret_cast<pointer>(pyArray()->data);
1130  vigra_precondition(this->checkInnerStride(Stride()),
1131  "NumpyArray<..., UnstridedArrayTag>::setupArrayView(): First dimension of given array is not unstrided (should never happen).");
1132 
1133  }
1134  else
1135  {
1136  this->m_ptr = 0;
1137  }
1138 }
1139 
1140 
1141 typedef NumpyArray<2, float > NumpyFArray2;
1142 typedef NumpyArray<3, float > NumpyFArray3;
1143 typedef NumpyArray<4, float > NumpyFArray4;
1144 typedef NumpyArray<2, Singleband<float> > NumpyFImage;
1145 typedef NumpyArray<3, Singleband<float> > NumpyFVolume;
1146 typedef NumpyArray<2, RGBValue<float> > NumpyFRGBImage;
1147 typedef NumpyArray<3, RGBValue<float> > NumpyFRGBVolume;
1148 typedef NumpyArray<3, Multiband<float> > NumpyFMultibandImage;
1149 typedef NumpyArray<4, Multiband<float> > NumpyFMultibandVolume;
1150 
1151 /********************************************************/
1152 /* */
1153 /* NumpyArray Multiband Argument Object Factories */
1154 /* */
1155 /********************************************************/
1156 
1157 template <class PixelType, class Stride>
1158 inline triple<ConstStridedImageIterator<PixelType>,
1159  ConstStridedImageIterator<PixelType>,
1160  MultibandVectorAccessor<PixelType> >
1161 srcImageRange(NumpyArray<3, Multiband<PixelType>, Stride> const & img)
1162 {
1163  ConstStridedImageIterator<PixelType>
1164  ul(img.data(), 1, img.stride(0), img.stride(1));
1165  return triple<ConstStridedImageIterator<PixelType>,
1166  ConstStridedImageIterator<PixelType>,
1167  MultibandVectorAccessor<PixelType> >
1168  (ul, ul + Size2D(img.shape(0), img.shape(1)), MultibandVectorAccessor<PixelType>(img.shape(2), img.stride(2)));
1169 }
1170 
1171 template <class PixelType, class Stride>
1172 inline pair< ConstStridedImageIterator<PixelType>,
1173  MultibandVectorAccessor<PixelType> >
1174 srcImage(NumpyArray<3, Multiband<PixelType>, Stride> const & img)
1175 {
1176  ConstStridedImageIterator<PixelType>
1177  ul(img.data(), 1, img.stride(0), img.stride(1));
1178  return pair<ConstStridedImageIterator<PixelType>, MultibandVectorAccessor<PixelType> >
1179  (ul, MultibandVectorAccessor<PixelType>(img.shape(2), img.stride(2)));
1180 }
1181 
1182 template <class PixelType, class Stride>
1183 inline triple< StridedImageIterator<PixelType>,
1184  StridedImageIterator<PixelType>,
1185  MultibandVectorAccessor<PixelType> >
1186 destImageRange(NumpyArray<3, Multiband<PixelType>, Stride> & img)
1187 {
1188  StridedImageIterator<PixelType>
1189  ul(img.data(), 1, img.stride(0), img.stride(1));
1190  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
1191  return triple<StridedImageIterator<PixelType>,
1192  StridedImageIterator<PixelType>,
1193  MultibandVectorAccessor<PixelType> >
1194  (ul, ul + Size2D(img.shape(0), img.shape(1)),
1195  MultibandVectorAccessor<PixelType>(img.shape(2), img.stride(2)));
1196 }
1197 
1198 template <class PixelType, class Stride>
1199 inline pair< StridedImageIterator<PixelType>,
1200  MultibandVectorAccessor<PixelType> >
1201 destImage(NumpyArray<3, Multiband<PixelType>, Stride> & img)
1202 {
1203  StridedImageIterator<PixelType>
1204  ul(img.data(), 1, img.stride(0), img.stride(1));
1205  return pair<StridedImageIterator<PixelType>, MultibandVectorAccessor<PixelType> >
1206  (ul, MultibandVectorAccessor<PixelType>(img.shape(2), img.stride(2)));
1207 }
1208 
1209 template <class PixelType, class Stride>
1210 inline pair< ConstStridedImageIterator<PixelType>,
1211  MultibandVectorAccessor<PixelType> >
1212 maskImage(NumpyArray<3, Multiband<PixelType>, Stride> const & img)
1213 {
1214  ConstStridedImageIterator<PixelType>
1215  ul(img.data(), 1, img.stride(0), img.stride(1));
1216  typedef typename AccessorTraits<PixelType>::default_accessor Accessor;
1217  return pair<ConstStridedImageIterator<PixelType>, MultibandVectorAccessor<PixelType> >
1218  (ul, MultibandVectorAccessor<PixelType>(img.shape(2), img.stride(2)));
1219 }
1220 
1221 } // namespace vigra
1222 
1223 #endif // VIGRA_NUMPY_ARRAY_HXX

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

html generated using doxygen and Python
vigra 1.9.0 (Wed Feb 27 2013)