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