Hi, I sort of missed the big C++ discussion, but I'd like to give some examples of how writing code can become much simpler if you are based on C++. This is from my mahotas package, which has a thin C++ wrapper around numpy's C API
https://github.com/luispedro/mahotas/blob/master/mahotas/_morph.cpp and it implements multi-type greyscale erosion. // numpy::aligned_array wraps PyArrayObject* template<typename T> void erode(numpy::aligned_array<T> res, numpy::aligned_array<T> array, numpy::aligned_array<T> Bc) { // Release the GIL using RAII gil_release nogil; const int N = res.size(); typename numpy::aligned_array<T>::iterator iter = array.begin(); // this is adapted from scipy.ndimage. // it implements the convolution-like filtering. filter_iterator<T> filter(res.raw_array(), Bc.raw_array(), EXTEND_NEAREST, is_bool(T())); const int N2 = filter.size(); T* rpos = res.data(); for (int i = 0; i != N; ++i, ++rpos, filter.iterate_both(iter)) { T value = std::numeric_limits<T>::max(); for (int j = 0; j != N2; ++j) { T arr_val = T(); filter.retrieve(iter, j, arr_val); value = std::min<T>(value, erode_sub(arr_val, filter[j])); } *rpos = value; } } If you compare this with the equivalent scipy.ndimage function, which is very good C code (but mostly write-only—in fact, ndimage has not been maintainable because it is so hard [at least for me, I've tried]): int NI_BinaryErosion(PyArrayObject* input, PyArrayObject* strct, PyArrayObject* mask, PyArrayObject* output, int bdr_value, npy_intp *origins, int invert, int center_is_true, int* changed, NI_CoordinateList **coordinate_list) { npy_intp struct_size = 0, *offsets = NULL, size, *oo, jj; npy_intp ssize, block_size = 0, *current = NULL, border_flag_value; int kk, true, false, msk_value; NI_Iterator ii, io, mi; NI_FilterIterator fi; Bool *ps, out = 0; char *pi, *po, *pm = NULL; NI_CoordinateBlock *block = NULL; ps = (Bool*)PyArray_DATA(strct); ssize = 1; for(kk = 0; kk < strct->nd; kk++) ssize *= strct->dimensions[kk]; for(jj = 0; jj < ssize; jj++) if (ps[jj]) ++struct_size; if (mask) { if (!NI_InitPointIterator(mask, &mi)) return 0; pm = (void *)PyArray_DATA(mask); } /* calculate the filter offsets: */ if (!NI_InitFilterOffsets(input, ps, strct->dimensions, origins, NI_EXTEND_CONSTANT, &offsets, &border_flag_value, NULL)) goto exit; /* initialize input element iterator: */ if (!NI_InitPointIterator(input, &ii)) goto exit; /* initialize output element iterator: */ if (!NI_InitPointIterator(output, &io)) goto exit; /* initialize filter iterator: */ if (!NI_InitFilterIterator(input->nd, strct->dimensions, struct_size, input->dimensions, origins, &fi)) goto exit; /* get data pointers an size: */ pi = (void *)PyArray_DATA(input); po = (void *)PyArray_DATA(output); size = 1; for(kk = 0; kk < input->nd; kk++) size *= input->dimensions[kk]; if (invert) { bdr_value = bdr_value ? 0 : 1; true = 0; false = 1; } else { bdr_value = bdr_value ? 1 : 0; true = 1; false = 0; } if (coordinate_list) { block_size = LIST_SIZE / input->nd / sizeof(int); if (block_size < 1) block_size = 1; if (block_size > size) block_size = size; *coordinate_list = NI_InitCoordinateList(block_size, input->nd); if (!*coordinate_list) goto exit; } /* iterator over the elements: */ oo = offsets; *changed = 0; msk_value = 1; for(jj = 0; jj < size; jj++) { int pchange = 0; if (mask) { switch(mask->descr->type_num) { CASE_GET_MASK(msk_value, pm, Bool); CASE_GET_MASK(msk_value, pm, UInt8); CASE_GET_MASK(msk_value, pm, UInt16); CASE_GET_MASK(msk_value, pm, UInt32); #if HAS_UINT64 CASE_GET_MASK(msk_value, pm, UInt64); #endif CASE_GET_MASK(msk_value, pm, Int8); CASE_GET_MASK(msk_value, pm, Int16); CASE_GET_MASK(msk_value, pm, Int32); CASE_GET_MASK(msk_value, pm, Int64); CASE_GET_MASK(msk_value, pm, Float32); CASE_GET_MASK(msk_value, pm, Float64); default: PyErr_SetString(PyExc_RuntimeError, "data type not supported"); return 0; } } switch (input->descr->type_num) { CASE_NI_ERODE_POINT(pi, out, oo, struct_size, Bool, msk_value, bdr_value, border_flag_value, center_is_true, true, false, pchange); CASE_NI_ERODE_POINT(pi, out, oo, struct_size, UInt8, msk_value, bdr_value, border_flag_value, center_is_true, true, false, pchange); CASE_NI_ERODE_POINT(pi, out, oo, struct_size, UInt16, msk_value, bdr_value, border_flag_value, center_is_true, true, false, pchange); CASE_NI_ERODE_POINT(pi, out, oo, struct_size, UInt32, msk_value, bdr_value, border_flag_value, center_is_true, true, false, pchange); #if HAS_UINT64 CASE_NI_ERODE_POINT(pi, out, oo, struct_size, UInt64, msk_value, bdr_value, border_flag_value, center_is_true, true, false, pchange); #endif CASE_NI_ERODE_POINT(pi, out, oo, struct_size, Int8, msk_value, bdr_value, border_flag_value, center_is_true, true, false, pchange); CASE_NI_ERODE_POINT(pi, out, oo, struct_size, Int16, msk_value, bdr_value, border_flag_value, center_is_true, true, false, pchange); CASE_NI_ERODE_POINT(pi, out, oo, struct_size, Int32, msk_value, bdr_value, border_flag_value, center_is_true, true, false, pchange); CASE_NI_ERODE_POINT(pi, out, oo, struct_size, Int64, msk_value, bdr_value, border_flag_value, center_is_true, true, false, pchange); CASE_NI_ERODE_POINT(pi, out, oo, struct_size, Float32, msk_value, bdr_value, border_flag_value, center_is_true, true, false, pchange); CASE_NI_ERODE_POINT(pi, out, oo, struct_size, Float64, msk_value, bdr_value, border_flag_value, center_is_true, true, false, pchange); default: PyErr_SetString(PyExc_RuntimeError, "data type not supported"); goto exit; } switch (output->descr->type_num) { CASE_OUTPUT(po, out, Bool); CASE_OUTPUT(po, out, UInt8); CASE_OUTPUT(po, out, UInt16); CASE_OUTPUT(po, out, UInt32); #if HAS_UINT64 CASE_OUTPUT(po, out, UInt64); #endif CASE_OUTPUT(po, out, Int8); CASE_OUTPUT(po, out, Int16); CASE_OUTPUT(po, out, Int32); CASE_OUTPUT(po, out, Int64); CASE_OUTPUT(po, out, Float32); CASE_OUTPUT(po, out, Float64); default: PyErr_SetString(PyExc_RuntimeError, "data type not supported"); goto exit; } if (pchange) { *changed = 1; if (coordinate_list) { if (block == NULL || block->size == block_size) { block = NI_CoordinateListAddBlock(*coordinate_list); current = block->coordinates; } for(kk = 0; kk < input->nd; kk++) *current++ = ii.coordinates[kk]; block->size++; } } if (mask) { NI_FILTER_NEXT3(fi, ii, io, mi, oo, pi, po, pm); } else { NI_FILTER_NEXT2(fi, ii, io, oo, pi, po); } } exit: if (offsets) free(offsets); if (PyErr_Occurred()) { if (coordinate_list) { NI_FreeCoordinateList(*coordinate_list); *coordinate_list = NULL; } return 0; } else { return 1; } return PyErr_Occurred() ? 0 : 1; } HTH -- Luis Pedro Coelho | Institute for Molecular Medicine | http://luispedro.org
signature.asc
Description: This is a digitally signed message part.
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
