[Numpy-discussion] C++ Example
Luis Pedro Coelho
lpc@cmu....
Sat Mar 3 10:07:16 CST 2012
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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 198 bytes
Desc: This is a digitally signed message part.
Url : http://mail.scipy.org/pipermail/numpy-discussion/attachments/20120303/6aa473e5/attachment-0001.bin
More information about the NumPy-Discussion
mailing list