[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