[Numpy-discussion] best way of speeding up a filtering-like algorithm

[Moroney, Catherine M (398E)]

Hello,

I have the following sample code (pretty simple algorithm that uses a rolling filter window) and am wondering what the best way is of speeding it up.  I tried rewriting it in Cython by pre-declaring the variables but that didn’t buy me a lot of time.  Then I rewrote it in Fortran (and compiled it with f2py) and now it’s lightning fast.  But I would still like to know if I could rewrite it in pure python/numpy/scipy or in Cython and get a similar speedup.

Here is the raw Python code:

def mixed_coastline_slow(nsidc, radius, count, mask=None):

    nsidc_copy = numpy.copy(nsidc)

    if (mask is None):
        idx_coastline = numpy.where(nsidc_copy == NSIDC_COASTLINE_MIXED)
    else:
        idx_coastline = numpy.where(mask & (nsidc_copy == NSIDC_COASTLINE_MIXED))

    for (irow0, icol0) in zip(idx_coastline[0], idx_coastline[1]):

        rows = ( max(irow0-radius, 0), min(irow0+radius+1, nsidc_copy.shape[0]) )
        cols = ( max(icol0-radius, 0), min(icol0+radius+1, nsidc_copy.shape[1]) )
        window = nsidc[rows[0]:rows[1], cols[0]:cols[1]]

        npoints = numpy.where(window != NSIDC_COASTLINE_MIXED, True, False).sum()
        nsnowice = numpy.where( (window >= NSIDC_SEAICE_LOW) & (window <= NSIDC_FRESHSNOW), \
                                True, False).sum()

        if (100.0*nsnowice/npoints >= count):
             nsidc_copy[irow0, icol0] = MISR_SEAICE_THRESHOLD

    return nsidc_copy

and here is my attempt at Cython-izing it:

import numpy
cimport numpy as cnumpy
cimport cython

cdef int NSIDC_SIZE  = 721
cdef int NSIDC_NO_SNOW = 0
cdef int NSIDC_ALL_SNOW = 100
cdef int NSIDC_FRESHSNOW = 103
cdef int NSIDC_PERMSNOW  = 101
cdef int NSIDC_SEAICE_LOW  = 1
cdef int NSIDC_SEAICE_HIGH = 100
cdef int NSIDC_COASTLINE_MIXED = 252
cdef int NSIDC_SUSPECT_ICE = 253

cdef int MISR_SEAICE_THRESHOLD = 6

def mixed_coastline(cnumpy.ndarray[cnumpy.uint8_t, ndim=2] nsidc, int radius, int count):

     cdef int irow, icol, irow1, irow2, icol1, icol2, npoints, nsnowice
     cdef cnumpy.ndarray[cnumpy.uint8_t, ndim=2] nsidc2 \
        = numpy.empty(shape=(NSIDC_SIZE, NSIDC_SIZE), dtype=numpy.uint8)
     cdef cnumpy.ndarray[cnumpy.uint8_t, ndim=2] window \
        = numpy.empty(shape=(2*radius+1, 2*radius+1), dtype=numpy.uint8)

     nsidc2 = numpy.copy(nsidc)

     idx_coastline = numpy.where(nsidc2 == NSIDC_COASTLINE_MIXED)

     for (irow, icol) in zip(idx_coastline[0], idx_coastline[1]):

          irow1 = max(irow-radius, 0)
          irow2 = min(irow+radius+1, NSIDC_SIZE)
          icol1 = max(icol-radius, 0)
          icol2 = min(icol+radius+1, NSIDC_SIZE)
          window = nsidc[irow1:irow2, icol1:icol2]

          npoints = numpy.where(window != NSIDC_COASTLINE_MIXED, True, False).sum()
          nsnowice = numpy.where( (window >= NSIDC_SEAICE_LOW) & (window <= NSIDC_FRESHSNOW), \
                                  True, False).sum()

          if (100.0*nsnowice/npoints >= count):
               nsidc2[irow, icol] = MISR_SEAICE_THRESHOLD

     return nsidc2

Thanks in advance for any advice!

Catherine

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20180329/3874cc1b/attachment.html>