5

I have spent lots of time on this, and I know how to manually do it by slicing and indexing the boundary rows/cols, but there has to be a simpler way with SciPy.

I need to set the CVAL (values to fill past the edges when mode=constant) to NaN, however, this will return NaN.

I will explain it with code and figures:

import numpy as np
from scipy import ndimage
m = np.reshape(np.arange(0,100),(10,10)).astype(np.float)

picture of the input array Use SciPy ndimage uniform filter to calculate the mean using a 3x3 kernel:

filter = ndimage.uniform_filter(m, size=3, mode='constant')
print(filter[1][1]) # equal to 11
print(filter[9][9]) # I need 93.5, however it gets 41.55 due to zeros

As you can see, the first value comes out as 11, which is as expected, however, for any cell along the border, it will fill the values with zero (I have also tried all the other modes).

Here is what I need to achieve (left) vs mode=constant and CVAL=0 (default 0)

scipy avg and what I need

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
dfresh22
  • 961
  • 1
  • 15
  • 23

2 Answers2

4

One simple approach is to use Normalized Convolution:

import numpy as np
from scipy import ndimage
m = np.reshape(np.arange(0,100),(10,10)).astype(np.float)

filter = ndimage.uniform_filter(m, size=3, mode='constant')    # normal filter result

weights = ndimage.uniform_filter(np.ones(m.shape), size=3, mode='constant')
filter = filter / weights    # normalized convolution result

print(filter[1][1]) # equal to 11
print(filter[9][9]) # equal to 93.49999999999994 -- rounding error! :)

We computed the result of the filter if all data points were 1 (weights). This shows how many data elements there are in each filter window, and returns a value of 1 everywhere except near the boundary, where this value decreases proportionally. By dividing the filtering result with these weights, we correct for the averaging taking zeros into account that were outside the data domain.

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
1

This suggestion is less than ideal, because it will be slow compared to uniform_filter, but it will do what you want.

Using your idea of using nan for the constant value, you can implement the uniform filter by using ndimage.generic_filter instead of uniform_filter, with numpy.nanmean as the generic filter function.

For example, here's your sample array m:

In [102]: import numpy as np

In [103]: m = np.reshape(np.arange(0,100),(10,10)).astype(np.float)

Apply generic_filter, with numpy.nanmean as the function to be applied:

In [104]: from scipy.ndimage import generic_filter

In [105]: generic_filter(m, np.nanmean, mode='constant', cval=np.nan, size=3)
Out[105]: 
array([[ 5.5,  6. ,  7. ,  8. ,  9. , 10. , 11. , 12. , 13. , 13.5],
       [10.5, 11. , 12. , 13. , 14. , 15. , 16. , 17. , 18. , 18.5],
       [20.5, 21. , 22. , 23. , 24. , 25. , 26. , 27. , 28. , 28.5],
       [30.5, 31. , 32. , 33. , 34. , 35. , 36. , 37. , 38. , 38.5],
       [40.5, 41. , 42. , 43. , 44. , 45. , 46. , 47. , 48. , 48.5],
       [50.5, 51. , 52. , 53. , 54. , 55. , 56. , 57. , 58. , 58.5],
       [60.5, 61. , 62. , 63. , 64. , 65. , 66. , 67. , 68. , 68.5],
       [70.5, 71. , 72. , 73. , 74. , 75. , 76. , 77. , 78. , 78.5],
       [80.5, 81. , 82. , 83. , 84. , 85. , 86. , 87. , 88. , 88.5],
       [85.5, 86. , 87. , 88. , 89. , 90. , 91. , 92. , 93. , 93.5]])
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
  • ahh yes the Generic Filter, thank you! I will check the performance. typically I will be doing on arrays of shape (1e6, 1e6). – dfresh22 Dec 17 '18 at 21:17
  • I suspect you won't be happy. – Warren Weckesser Dec 17 '18 at 21:19
  • 1
    yup you were right, that is awfully slow. m = np.random.rand(1000,1000 the uniform filter is about 0.01 seconds versus the generic filter which is about 22 seconds on my computer. – dfresh22 Dec 17 '18 at 21:25
  • Probably your best bet is to use `uniform_filter` on the full image, and then add a few more lines to fix the boundary. From what you wrote in the question, I think that's what you were already doing. – Warren Weckesser Dec 17 '18 at 21:31