36 #ifndef VIGRA_NUMPY_ARRAY_HXX
37 #define VIGRA_NUMPY_ARRAY_HXX
39 #ifndef NPY_NO_DEPRECATED_API
40 # define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
46 #include <numpy/arrayobject.h>
47 #include "multi_array.hxx"
48 #include "array_vector.hxx"
49 #include "python_utility.hxx"
50 #include "numpy_array_traits.hxx"
51 #include "numpy_array_taggedshape.hxx"
58 static inline void import_vigranumpy()
61 if(_import_array() < 0)
62 pythonToCppException(0);
66 char const * load_vigra =
68 "if not sys.modules.has_key('vigra.vigranumpycore'):\n"
70 pythonToCppException(PyRun_SimpleString(load_vigra) == 0);
80 class MultibandVectorAccessor
91 typedef Multiband<T> value_type;
95 typedef T component_type;
97 typedef VectorElementAccessor<MultibandVectorAccessor<T> > ElementAccessor;
102 template <
class ITERATOR>
103 component_type
const & getComponent(ITERATOR
const & i,
int idx)
const
105 return *(&*i+idx*stride_);
113 template <
class V,
class ITERATOR>
114 void setComponent(V
const & value, ITERATOR
const & i,
int idx)
const
116 *(&*i+idx*stride_) = detail::RequiresExplicitCast<component_type>::cast(value);
122 template <
class ITERATOR,
class DIFFERENCE>
123 component_type
const & getComponent(ITERATOR
const & i, DIFFERENCE
const & diff,
int idx)
const
125 return *(&i[diff]+idx*stride_);
133 template <
class V,
class ITERATOR,
class DIFFERENCE>
135 setComponent(V
const & value, ITERATOR
const & i, DIFFERENCE
const & diff,
int idx)
const
137 *(&i[diff]+idx*stride_) = detail::RequiresExplicitCast<component_type>::cast(value);
149 template <
class TYPECODE>
152 constructArray(TaggedShape tagged_shape, TYPECODE typeCode,
bool init,
153 python_ptr arraytype = python_ptr());
157 template <
class Shape>
158 void numpyParseSlicing(Shape
const & shape, PyObject * idx, Shape & start, Shape & stop)
160 int N = shape.size();
161 for(
int k=0; k<N; ++k)
167 python_ptr index(idx);
168 if(!PySequence_Check(index))
170 index = python_ptr(PyTuple_Pack(1, index.ptr()), python_ptr::new_nonzero_reference);
172 int lindex = PyTuple_Size(index);
174 for(; kindex<lindex; ++kindex)
176 if(PyTuple_GET_ITEM((PyTupleObject *)index.ptr(), kindex) == Py_Ellipsis)
179 if(kindex == lindex && lindex < N)
181 python_ptr ellipsis = python_ptr(PyTuple_Pack(1, Py_Ellipsis), python_ptr::new_nonzero_reference);
182 index = python_ptr(PySequence_Concat(index, ellipsis), python_ptr::new_nonzero_reference);
186 for(
int k=0; k < N; ++k)
188 PyObject * item = PyTuple_GET_ITEM((PyTupleObject *)index.ptr(), kindex);
189 if(PyInt_Check(item))
194 start[k] += shape[k];
198 else if(PySlice_Check(item))
200 Py_ssize_t sstart, sstop, step;
201 if(PySlice_GetIndices((PySliceObject *)item, shape[k], &sstart, &sstop, &step) != 0)
202 pythonToCppException(0);
203 vigra_precondition(step == 1,
204 "numpyParseSlicing(): only unit steps are supported.");
209 else if(item == Py_Ellipsis)
218 vigra_precondition(
false,
219 "numpyParseSlicing(): unsupported index object.");
252 static python_ptr getArrayTypeObject()
254 return detail::getArrayTypeObject();
257 static std::string defaultOrder(std::string defaultValue =
"C")
259 return detail::defaultOrder(defaultValue);
262 static python_ptr defaultAxistags(
int ndim, std::string order =
"")
264 return detail::defaultAxistags(ndim, order);
267 static python_ptr emptyAxistags(
int ndim)
269 return detail::emptyAxistags(ndim);
279 explicit NumpyAnyArray(PyObject * obj = 0,
bool createCopy =
false, PyTypeObject * type = 0)
283 vigra_precondition(type == 0 || PyType_IsSubtype(type, &PyArray_Type),
284 "NumpyAnyArray(obj, createCopy, type): type must be numpy.ndarray or a subclass thereof.");
288 vigra_precondition(
makeReference(obj, type),
"NumpyAnyArray(obj): obj isn't a numpy array.");
300 vigra_precondition(type == 0 || PyType_IsSubtype(type, &PyArray_Type),
301 "NumpyAnyArray(obj, createCopy, type): type must be numpy.ndarray or a subclass thereof.");
324 vigra_precondition(other.
hasData(),
325 "NumpyArray::operator=(): Cannot assign from empty array.");
327 python_ptr arraytype = getArrayTypeObject();
328 python_ptr f(PyString_FromString(
"_copyValuesImpl"), python_ptr::keep_count);
329 if(PyObject_HasAttr(arraytype, f))
331 python_ptr res(PyObject_CallMethodObjArgs(arraytype, f.get(),
332 pyArray_.get(), other.pyArray_.get(), NULL),
333 python_ptr::keep_count);
334 vigra_postcondition(res.get() != 0,
335 "NumpyArray::operator=(): VigraArray._copyValuesImpl() failed.");
339 PyArrayObject * sarray = (PyArrayObject *)pyArray_.get();
340 PyArrayObject * tarray = (PyArrayObject *)other.pyArray_.get();
342 if(PyArray_CopyInto(tarray, sarray) == -1)
343 pythonToCppException(0);
348 pyArray_ = other.pyArray_;
360 return PyArray_NDIM(
pyArray());
373 return pythonGetAttr(
pyObject(),
"spatialDimensions",
ndim());
376 bool hasChannelAxis()
const
380 return channelIndex() ==
ndim();
387 return pythonGetAttr(
pyObject(),
"channelIndex",
ndim());
394 return pythonGetAttr(
pyObject(),
"innerNonchannelIndex",
ndim());
426 if(stride[j] < stride[smallest])
431 std::swap(stride[k], stride[smallest]);
432 std::swap(permutation[k], permutation[smallest]);
437 ordering[permutation[k]] = k;
470 return PyArray_DESCR(
pyArray())->type_num;
477 template <
class Shape>
481 unsigned int size =
ndim();
482 vigra_precondition(start.size() == size && stop.size() == size,
483 "NumpyAnyArray::getitem(): shape has wrong dimension.");
487 python_ptr index(PyTuple_New(size), python_ptr::new_nonzero_reference);
488 for(
unsigned int k=0; k<size; ++k)
494 vigra_precondition(0 <= start[k] && start[k] <= stop[k] && stop[k] <= s[k],
495 "NumpyAnyArray::getitem(): slice out of bounds.");
497 if(start[k] == stop[k])
499 item = PyInt_FromLong(start[k]);
503 python_ptr s0(PyInt_FromLong(start[k]), python_ptr::new_nonzero_reference);
504 python_ptr s1(PyInt_FromLong(stop[k]), python_ptr::new_nonzero_reference);
505 item = PySlice_New(s0, s1, 0);
507 pythonToCppException(item);
508 PyTuple_SET_ITEM((PyTupleObject *)index.ptr(), k, item);
511 python_ptr func(PyString_FromString(
"__getitem__"), python_ptr::new_nonzero_reference);
512 python_ptr res(PyObject_CallMethodObjArgs(
pyObject(), func.ptr(), index.ptr(), NULL),
513 python_ptr::new_nonzero_reference);
527 python_ptr key(PyString_FromString(
"axistags"), python_ptr::keep_count);
528 axistags.reset(PyObject_GetAttr(
pyObject(), key), python_ptr::keep_count);
540 return (PyArrayObject *)pyArray_.get();
549 return pyArray_.get();
562 if(obj == 0 || !PyArray_Check(obj))
566 vigra_precondition(PyType_IsSubtype(type, &PyArray_Type) != 0,
567 "NumpyAnyArray::makeReference(obj, type): type must be numpy.ndarray or a subclass thereof.");
568 obj = PyArray_View((PyArrayObject*)obj, 0, type);
569 pythonToCppException(obj);
581 void makeCopy(PyObject * obj, PyTypeObject * type = 0)
583 vigra_precondition(obj && PyArray_Check(obj),
584 "NumpyAnyArray::makeCopy(obj): obj is not an array.");
585 vigra_precondition(type == 0 || PyType_IsSubtype(type, &PyArray_Type),
586 "NumpyAnyArray::makeCopy(obj, type): type must be numpy.ndarray or a subclass thereof.");
587 python_ptr array(PyArray_NewCopy((PyArrayObject*)obj, NPY_ANYORDER), python_ptr::keep_count);
588 pythonToCppException(array);
597 return pyArray_ != 0;
610 nontrivialPermutation(ArrayVector<npy_intp>
const & p)
612 for(
unsigned int k=0; k<p.size(); ++k)
620 template <
class TYPECODE>
623 constructArray(TaggedShape tagged_shape, TYPECODE typeCode,
bool init, python_ptr arraytype)
625 ArrayVector<npy_intp> shape = finalizeTaggedShape(tagged_shape);
626 PyAxisTags axistags(tagged_shape.axistags);
628 int ndim = (int)shape.size();
629 ArrayVector<npy_intp> inverse_permutation;
635 arraytype = NumpyAnyArray::getArrayTypeObject();
637 inverse_permutation = axistags.permutationFromNormalOrder();
638 vigra_precondition(ndim == (
int)inverse_permutation.size(),
639 "axistags.permutationFromNormalOrder(): permutation has wrong size.");
643 arraytype = python_ptr((PyObject*)&PyArray_Type);
649 python_ptr array(PyArray_New((PyTypeObject *)arraytype.get(), ndim, shape.begin(),
650 typeCode, 0, 0, 0, order, 0),
651 python_ptr::keep_count);
652 pythonToCppException(array);
654 if(detail::nontrivialPermutation(inverse_permutation))
656 PyArray_Dims permute = { inverse_permutation.begin(), ndim };
657 array = python_ptr(PyArray_Transpose((PyArrayObject*)array.get(), &permute),
658 python_ptr::keep_count);
659 pythonToCppException(array);
662 if(arraytype != (PyObject*)&PyArray_Type && axistags)
663 pythonToCppException(PyObject_SetAttrString(array,
"axistags", axistags.axistags) != -1);
666 PyArray_FILLWBYTE((PyArrayObject *)array.get(), 0);
668 return array.release();
672 template <
class TINY_VECTOR>
674 python_ptr constructNumpyArrayFromData(
675 TINY_VECTOR
const & shape, npy_intp *strides,
676 NPY_TYPES typeCode,
void *data)
678 ArrayVector<npy_intp> pyShape(shape.begin(), shape.end());
680 #ifndef NPY_ARRAY_WRITEABLE
681 # define NPY_ARRAY_WRITEABLE NPY_WRITEABLE // old API compatibility
684 python_ptr array(PyArray_New(&PyArray_Type, shape.size(), pyShape.begin(),
685 typeCode, strides, data, 0, NPY_ARRAY_WRITEABLE, 0),
686 python_ptr::keep_count);
687 pythonToCppException(array);
706 template <
unsigned int N,
class T,
class Str
ide = Str
idedArrayTag>
708 :
public MultiArrayView<N, typename NumpyArrayTraits<N, T, Stride>::value_type, Stride>,
712 typedef NumpyArrayTraits<N, T, Stride> ArrayTraits;
713 typedef typename ArrayTraits::dtype
dtype;
714 typedef T pseudo_value_type;
716 static NPY_TYPES
const typeCode = ArrayTraits::typeCode;
722 enum { actual_dimension = view_type::actual_dimension };
783 void setupArrayView();
785 static python_ptr init(
difference_type const & shape,
bool init =
true,
786 std::string
const & order =
"")
788 vigra_precondition(order ==
"" || order ==
"C" || order ==
"F" ||
789 order ==
"V" || order ==
"A",
790 "NumpyArray.init(): order must be in ['C', 'F', 'V', 'A', ''].");
791 return python_ptr(constructArray(ArrayTraits::taggedShape(shape, order), typeCode, init),
792 python_ptr::keep_count);
812 explicit NumpyArray(PyObject *obj = 0,
bool createCopy =
false)
820 "NumpyArray(obj): Cannot construct from incompatible array.");
845 template <
class U,
class S>
851 "NumpyArray(MultiArrayView): Python constructor did not produce a compatible array.");
865 "NumpyArray(shape): Python constructor did not produce a compatible array.");
877 "NumpyArray(tagged_shape): Python constructor did not produce a compatible array.");
892 "NumpyArray(NumpyAnyArray): Cannot construct from incompatible or empty array.");
920 template <
class U,
class S>
926 "NumpyArray::operator=(): shape mismatch.");
933 "NumpyArray::operator=(): reshape failed unexpectedly.");
946 template <
class U,
class S>
952 "NumpyArray::operator=(): shape mismatch.");
959 "NumpyArray::operator=(): reshape failed unexpectedly.");
986 vigra_precondition(
false,
987 "NumpyArray::operator=(): Cannot assign from incompatible array.");
1001 "NumpyArray::permuteLikewise(): array has no data.");
1004 ArrayTraits::permuteLikewise(this->pyArray_, data, res);
1012 template <
class U,
int K>
1017 "NumpyArray::permuteLikewise(): array has no data.");
1020 ArrayTraits::permuteLikewise(this->pyArray_, data, res);
1033 "NumpyArray::permuteLikewise(): array has no data.");
1037 ArrayTraits::permuteLikewise(this->pyArray_, data, res);
1049 #if VIGRA_CONVERTER_DEBUG
1050 std::cerr <<
"class " <<
typeid(
NumpyArray).name() <<
" got " << obj->ob_type->tp_name <<
"\n";
1051 std::cerr <<
"using traits " <<
typeid(ArrayTraits).name() <<
"\n";
1052 std::cerr<<
"isArray: "<< ArrayTraits::isArray(obj)<<std::endl;
1053 std::cerr<<
"isShapeCompatible: "<< ArrayTraits::isShapeCompatible((PyArrayObject *)obj)<<std::endl;
1056 return ArrayTraits::isArray(obj) &&
1057 ArrayTraits::isShapeCompatible((PyArrayObject *)obj);
1068 return ArrayTraits::isArray(obj) &&
1069 ArrayTraits::isPropertyCompatible((PyArrayObject *)obj);
1088 for(
unsigned int k=0; k<N; ++k)
1089 strideOrdering[k] = k;
1142 vigra_precondition(!
hasData(),
1143 "makeUnsafeReference(): cannot replace existing view with given buffer");
1146 python_ptr array(ArrayTraits::unsafeConstructorFromData(multiArrayView.
shape(),
1147 multiArrayView.
data(), multiArrayView.
stride()));
1161 #if VIGRA_CONVERTER_DEBUG
1162 int ndim = PyArray_NDIM((PyArrayObject *)obj);
1163 npy_intp * s = PyArray_DIMS((PyArrayObject *)obj);
1166 std::cerr <<
"for " <<
typeid(*this).name() <<
"\n";
1169 "NumpyArray::makeCopy(obj): Cannot copy an incompatible array.");
1187 "NumpyArray.reshape(shape): Python constructor did not produce a compatible array.");
1209 ArrayTraits::finalizeTaggedShape(tagged_shape);
1213 vigra_precondition(tagged_shape.compatible(taggedShape()), message.c_str());
1217 python_ptr array(constructArray(tagged_shape, typeCode,
true),
1218 python_ptr::keep_count);
1220 "NumpyArray.reshapeIfEmpty(): Python constructor did not produce a compatible array.");
1224 TaggedShape taggedShape()
const
1226 return ArrayTraits::taggedShape(this->
shape(), PyAxisTags(this->
axistags(),
true));
1231 template <
unsigned int N,
class T,
class Str
ide>
1232 void NumpyArray<N, T, Stride>::setupArrayView()
1236 permutation_type permute;
1237 ArrayTraits::permutationToSetupOrder(this->pyArray_, permute);
1239 vigra_precondition(
abs((
int)permute.size() - actual_dimension) <= 1,
1240 "NumpyArray::setupArrayView(): got array of incompatible shape (should never happen).");
1243 PyArray_DIMS(pyArray()), this->m_shape.begin());
1245 PyArray_STRIDES(pyArray()), this->m_stride.begin());
1247 if((
int)permute.size() == actual_dimension - 1)
1249 this->m_shape[actual_dimension-1] = 1;
1250 this->m_stride[actual_dimension-1] =
sizeof(value_type);
1253 this->m_stride /=
sizeof(value_type);
1254 this->m_ptr =
reinterpret_cast<pointer
>(PyArray_DATA(pyArray()));
1255 vigra_precondition(this->checkInnerStride(Stride()),
1256 "NumpyArray<..., UnstridedArrayTag>::setupArrayView(): First dimension of given array is not unstrided (should never happen).");
1266 typedef NumpyArray<2, float > NumpyFArray2;
1267 typedef NumpyArray<3, float > NumpyFArray3;
1268 typedef NumpyArray<4, float > NumpyFArray4;
1269 typedef NumpyArray<2, Singleband<float> > NumpyFImage;
1270 typedef NumpyArray<3, Singleband<float> > NumpyFVolume;
1271 typedef NumpyArray<2, RGBValue<float> > NumpyFRGBImage;
1272 typedef NumpyArray<3, RGBValue<float> > NumpyFRGBVolume;
1273 typedef NumpyArray<3, Multiband<float> > NumpyFMultibandImage;
1274 typedef NumpyArray<4, Multiband<float> > NumpyFMultibandVolume;
1282 template <
class PixelType,
class Str
ide>
1283 inline triple<ConstStridedImageIterator<PixelType>,
1284 ConstStridedImageIterator<PixelType>,
1285 MultibandVectorAccessor<PixelType> >
1286 srcImageRange(NumpyArray<3, Multiband<PixelType>, Stride>
const & img)
1288 ConstStridedImageIterator<PixelType>
1289 ul(img.data(), 1, img.stride(0), img.stride(1));
1290 return triple<ConstStridedImageIterator<PixelType>,
1291 ConstStridedImageIterator<PixelType>,
1292 MultibandVectorAccessor<PixelType> >
1293 (ul, ul + Size2D(img.shape(0), img.shape(1)), MultibandVectorAccessor<PixelType>(img.shape(2), img.stride(2)));
1296 template <
class PixelType,
class Str
ide>
1297 inline pair< ConstStridedImageIterator<PixelType>,
1298 MultibandVectorAccessor<PixelType> >
1299 srcImage(NumpyArray<3, Multiband<PixelType>, Stride>
const & img)
1301 ConstStridedImageIterator<PixelType>
1302 ul(img.data(), 1, img.stride(0), img.stride(1));
1303 return pair<ConstStridedImageIterator<PixelType>, MultibandVectorAccessor<PixelType> >
1304 (ul, MultibandVectorAccessor<PixelType>(img.shape(2), img.stride(2)));
1307 template <
class PixelType,
class Str
ide>
1308 inline triple< StridedImageIterator<PixelType>,
1309 StridedImageIterator<PixelType>,
1310 MultibandVectorAccessor<PixelType> >
1311 destImageRange(NumpyArray<3, Multiband<PixelType>, Stride> & img)
1313 StridedImageIterator<PixelType>
1314 ul(img.data(), 1, img.stride(0), img.stride(1));
1315 return triple<StridedImageIterator<PixelType>,
1316 StridedImageIterator<PixelType>,
1317 MultibandVectorAccessor<PixelType> >
1318 (ul, ul + Size2D(img.shape(0), img.shape(1)),
1319 MultibandVectorAccessor<PixelType>(img.shape(2), img.stride(2)));
1322 template <
class PixelType,
class Str
ide>
1323 inline pair< StridedImageIterator<PixelType>,
1324 MultibandVectorAccessor<PixelType> >
1325 destImage(NumpyArray<3, Multiband<PixelType>, Stride> & img)
1327 StridedImageIterator<PixelType>
1328 ul(img.data(), 1, img.stride(0), img.stride(1));
1329 return pair<StridedImageIterator<PixelType>, MultibandVectorAccessor<PixelType> >
1330 (ul, MultibandVectorAccessor<PixelType>(img.shape(2), img.stride(2)));
1333 template <
class PixelType,
class Str
ide>
1334 inline pair< ConstStridedImageIterator<PixelType>,
1335 MultibandVectorAccessor<PixelType> >
1336 maskImage(NumpyArray<3, Multiband<PixelType>, Stride>
const & img)
1338 ConstStridedImageIterator<PixelType>
1339 ul(img.data(), 1, img.stride(0), img.stride(1));
1340 return pair<ConstStridedImageIterator<PixelType>, MultibandVectorAccessor<PixelType> >
1341 (ul, MultibandVectorAccessor<PixelType>(img.shape(2), img.stride(2)));
1346 #endif // VIGRA_NUMPY_ARRAY_HXX