Calculate the average of neighbor pixels for raster edge

Question

I have a small raster (or matrix) and I would like to know what is the difference between each pixel and the average of its' 8 neighbors in the case where the kernel window is on the edges and no values exist. When the 8 neighbor pixels have "real" values I can calculate the average of the 8 for each pixel, by using scipy.ndimage.convolve with kernel window as np.array([[0.125,0.125,0.125],[0.125,0,0.125],[0.125,0.125,0.125]]).

The problem is on the edges.

For example, if this is my matrix:

array([[13, 21, 13,  8],
       [ 5, 10, 22, 14],
       [21, 33,  9,  0],
       [ 0,  0,  0,  0]])

I would like to calculate the average of the neighbors of the upper left cell (value 13). The desired calculation will be

(21+10+5)/3 = 12

For the pixel in the first column and second row (value 5) the desired calculation will be -

(13+21+10+5)/4 = 12.25

Note that I need to ignored the central pixel (value 13 or 5) and all the cells that are outside of the 3x3 kernel window and also that the number of neighboring pixels with "real" values is different from pixel to pixel.

I read about the mode that deals with edges here and in other places but could not find a solution.

Answer 1

The easiest way is to combine numpy's NaN functions (in this case nanmean) and ndimage.generic_filter, like so:

import numpy as np
from scipy import ndimage

m = np.array([[13, 21, 13,  8],
              [ 5, 10, 22, 14],
              [21, 33,  9,  0],
              [ 0,  0,  0,  0]], dtype=np.float)

result = ndimage.generic_filter(a, np.nanmean, size=3, mode='constant', cval=np.NaN)

At this point result comes out as:

array([[ 12.25      ,  14.        ,  14.66666667,  14.25      ],
       [ 17.16666667,  16.33333333,  14.44444444,  11.        ],
       [ 11.5       ,  11.11111111,   9.77777778,   7.5       ],
       [ 13.5       ,  10.5       ,   7.        ,   2.25      ]])

If you want to ignore the centre cell as well you can do so by supplying a footprint:

mask = np.ones((3, 3))
mask[1, 1] = 0

result = ndimage.generic_filter(a, np.nanmean, footprint=mask, mode='constant', cval=np.NaN)

This returns the result

array([[ 12.        ,  12.6       ,  15.        ,  16.33333333],
       [ 19.6       ,  17.125     ,  13.5       ,  10.4       ],
       [  9.6       ,   8.375     ,   9.875     ,   9.        ],
       [ 18.        ,  12.6       ,   8.4       ,   3.        ]])

Which we can check for the top left cell as (21 + 5 + 10) / 3 == 12 and the second row, first column (value 5) becomes (13 + 21 + 10 + 21 + 33) / 5 == 19.6, and the bottom right (value 0) becomres (9 + 0 + 0) / 3 == 3.

Answer 2

@om_henners gives a generic_filter method that works well for small arrays, which is the intended use case from the original question; however, this method can be slow for medium and large arrays. A similar approach using convolve2d will produce identical results and can provide substantial speed improvements, as demonstrated below. With a (2048, 512) array, I see a speedup of ~300 when using the convolve2d approach; even with a (4, 4) array, I see a speedup of ~20 when using the convolve2d approach. Cheers!

import numpy as np
from scipy.signal import convolve2d
from scipy.ndimage import generic_filter
import time


def eight_neighbor_average_generic_filter(x):
    footprint = np.ones((3, 3))
    footprint[1, 1] = 0

    result = generic_filter(
        x.astype('float'), np.nanmean, footprint=footprint,
        mode='constant', cval=np.nan)

    return result


def eight_neighbor_average_convolve2d(x):
    kernel = np.ones((3, 3))
    kernel[1, 1] = 0

    neighbor_sum = convolve2d(
        x, kernel, mode='same',
        boundary='fill', fillvalue=0)

    num_neighbor = convolve2d(
        np.ones(x.shape), kernel, mode='same',
        boundary='fill', fillvalue=0)

    return neighbor_sum / num_neighbor


if __name__ == '__main__':
    # Test 8-neighbor averaging on small example array
    x = np.array([
        [13, 21, 13,  8],
        [ 5, 10, 22, 14],
        [21, 33,  9,  0],
        [ 0,  0,  0,  0]])

    y1 = eight_neighbor_average_generic_filter(x)
    y2 = eight_neighbor_average_convolve2d(x)
    print('\nmethods equivalent: %r' % np.alltrue(y1 == y2))

    # Test 8-neighbor averaging on "large" array
    x_large = np.random.randn(2048, 512)

    t0 = time.time()
    y1_large = eight_neighbor_average_generic_filter(x_large)
    t_generic_filter = time.time() - t0

    t0 = time.time()
    y2_large = eight_neighbor_average_convolve2d(x_large)
    t_convolve2d = time.time() - t0

    print('\nExecution times for "large" array:')
    print('    generic filter: %f s' % t_generic_filter)
    print('    convolve2d: %f s' % t_convolve2d)