7

Is there a filter similar to ndimage's generic_filter that supports vector output? I did not manage to make scipy.ndimage.filters.generic_filter return more than a scalar. Uncomment the line in the code below to get the error: TypeError: only length-1 arrays can be converted to Python scalars.

I'm looking for a generic filter that process 2D or 3D arrays and returns a vector at each point. Thus the output would have one added dimension. For the example below I'd expect something like this:

m.shape    # (10,10)
res.shape  # (10,10,2)

Example Code

import numpy as np
from scipy import ndimage

a = np.ones((10, 10)) * np.arange(10)

footprint = np.array([[1,1,1],
                    [1,0,1],
                    [1,1,1]])

def myfunc(x):
    r = sum(x)
    #r = np.array([1,1])  # uncomment this
    return r

res = ndimage.generic_filter(a, myfunc, footprint=footprint)
user12719
  • 211
  • 1
  • 6

2 Answers2

6

The generic_filter expects myfunc to return a scalar, never a vector. However, there is nothing that precludes myfunc from also adding information to, say, a list which is passed to myfunc as an extra argument.

Instead of using the array returned by generic_filter, we can generate our vector-valued array by reshaping this list.

For example,

import numpy as np
from scipy import ndimage

a = np.ones((10, 10)) * np.arange(10)

footprint = np.array([[1,1,1],
                      [1,0,1],
                      [1,1,1]])

ndim = 2
def myfunc(x, out):
    r = np.arange(ndim, dtype='float64')
    out.extend(r)
    return 0

result = []
ndimage.generic_filter(
    a, myfunc, footprint=footprint, extra_arguments=(result,))
result = np.array(result).reshape(a.shape+(ndim,))
unutbu
  • 842,883
  • 184
  • 1,785
  • 1,677
  • Using a "mutable default value" is indeed a clever trick, which I hadn't seen yet. I think your solution is elegant and concise and answers my question. Could you please comment on the performance. It looks like there is no overhead except for the manipulation of the `result` variable. Also, is `r[1:][::-1]` equivalent to `r[:0:-1]`? – user12719 Mar 01 '15 at 02:27
  • The performance of `generic_filter` is never excellent because it requires calling a Python function once for each pixel of the image. I would try to use it only as a last resort -- after searching for a way to express the computation with faster built-in NumPy/SciPy functions applied to whole arrays which off-load the majority of the work to underlying functions written in C/C++ or Fortran. – unutbu Mar 01 '15 at 02:36
  • You are right `r[:0:-1]` is equivalent to `r[1:][::-1]`. Due to my ability to easily get confused, I went with the (to my mind) simpler `r[1:][::-1]`. `r[:0:-1]` is a bit faster, but watch out for preoptimization. The real bottleneck in `myfunc` is likely to be in the computation of `r`, not in reversing it. – unutbu Mar 01 '15 at 02:39
  • If you are looking for speed and are willing to sacrifice accuracy, there may be an alternative trick you could use: `myfunc` must return a scalar, such as a `np.float128`. But NumPy also supports `np.float16`, so you can pack 8 `np.float16`s into one `np.float128`. So if you wanted to return the vector `np.arange(8, dtype=np.float16)` you could instead return the single scalar `np.arange(8, dtype=np.float16).view(np.float128)` and then use `res = res.view(np.float16)` to recover the 8 `np.float16` values. There is also a `np.complex256` dtype, which may allow you to pack 16 `np.float16`s. – unutbu Mar 01 '15 at 02:56
  • 1
    I've found a simpler way to achieve the desired result. Simply append all the vectors to a list. After `generic_filter` completes, convert the list to an array and then reshape it. This is also slightly faster than my first answer. – unutbu Mar 01 '15 at 11:00
  • I'm using unutbu's updated solution that is based on a list, `out.extend(r)`, and it works fine. As a warning to everybody applying this to 3D arrays. It's slow! If you don't need a flexible kernel, it is worth thinking about faster solutions. If your kernel is rectangular or cubic, combining 'shifted' copies of your array might be faster. – user12719 Mar 06 '15 at 17:39
0

I think I get what you're asking, but I'm not completely sure how does the ndimage.generic_filter work (how abstruse is the source!).

Here's just a simple wrapper function. This function will take in an array, all the parameters ndimage.generic_filter needs. Function returns an array where each element of the former array is now represented by an array with shape (2,), result of the function is stored as the second element of that array.

def generic_expand_filter(inarr, func, **kwargs):
    shape = inarr.shape
    res = np.empty((  shape+(2,) ))
    temp = ndimage.generic_filter(inarr, func, **kwargs)
    for row in range(shape[0]):
        for val in range(shape[1]):
            res[row][val][0] = inarr[row][val]
            res[row][val][1] = temp[row][val]
    return res

Output, where res denotes just the generic_filter and res2 denotes generic_expand_filter, of this function is:

>>> a.shape #same as res.shape
(10, 10)
>>> res2.shape
(10, 10, 2)

>>> a[0]
array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9.])
>>> res[0]
array([  3.,   8.,  16.,  24.,  32.,  40.,  48.,  56.,  64.,  69.])
>>> print(*res2[0], sep=", ") #this is just to avoid the vertical default output
[ 0.  3.], [ 1.  8.], [  2.  16.], [  3.  24.], [  4.  32.], [  5.  40.], [  6.  48.], [  7.  56.], [  8.  64.], [  9.  69.]

>>> a[0][0]
0.0
>>> res[0][0]
3.0
>>> res2[0][0]
array([ 0.,  3.])

Of course you probably don't want to save the old array, but instead have both fields as new results. Except I don't know what exactly you had in mind, if the two values you want stored are unrelated, just add a temp2 and func2 and call another generic_filter with the same **kwargs and store that as the first value.

However if you want an actual vector quantity that is calculated using multiple inarr elements, meaning that the two new created fields aren't independent, you are just going to have to write that kind of a function, one that takes in an array, idx, idy indices and returns a tuple\list\array value which you can then unpack and assign to the result.

ljetibo
  • 3,048
  • 19
  • 25