search

imeka / ndarray-ndimage

Multidimensional image processing for ndarray
Apache License 2.0
14 stars 2 forks source link

Add sliding window optimisation for uniform_filter/1d (closes #21) #23

Closed matsjoyce-refeyn closed 1 year ago

matsjoyce-refeyn commented 1 year ago

As mentioned in #20 and #21, scipy optimises the uniform filter using a sliding window accumulator. This PR adds that optimisation to this repo's implementation. The result is a ~25% speed up, however this is still between 50% and 100% slower than scipy, the remaining time difference seems to be mainly in the padding and overall structure, not the actual maths (as noted in #11).

Same unit tests pass as before, under the same conditions. This is basically a copy of inner_correlate1d, with the weights code stripped out and the inner loop replaced by the sliding window code.

matsjoyce-refeyn commented 1 year ago

Thanks for the review, that looks neater now. I thought there must be a sum function, but my Haskell brain takes over with the "everything is a fold" mentality and I forgot to update it 😄.

matsjoyce-refeyn commented 1 year ago

LGTM. I checked your tests from the previous MR and I wonder if we have the exact same behaviour for integer data. Can you please add one test?

The uniform_filter in master has a type bound requiring floatyness (which is mostly required because the weights are usually less than 1). However, with this PR I can actually remove it, and so I'll add some tests for that as well. The behaviour matches scipy, but the scipy behaviour itself is slightly unintuitive:

>>> uniform_filter(x, 4)
array([[12, 12, 12, 12, 14, 16, 17],
       [15, 16, 15, 15, 17, 19, 20],
       [24, 24, 25, 26, 28, 30, 31],
       [36, 36, 38, 40, 42, 44, 45],
       [46, 46, 48, 50, 52, 54, 55]])
>>> uniform_filter(x.astype(float), 4) 
array([[12.25  , 12.75  , 12.125 , 12.    , 14.    , 16.    , 17.5   ],
       [15.75  , 16.25  , 15.625 , 15.5   , 17.5   , 19.5   , 21.    ],
       [24.125 , 24.625 , 25.0625, 26.    , 28.    , 30.    , 31.5   ],
       [36.    , 36.5   , 38.    , 40.    , 42.    , 44.    , 45.5   ],
       [46.5   , 47.    , 48.5   , 50.5   , 52.5   , 54.5   , 56.    ]])

I'm not really sure why the bottom-right corner switches from 55 to 56.

nilgoyette commented 1 year ago

The behaviour matches scipy, but the scipy behaviour itself is slightly unintuitive:

What do you mean by that?

matsjoyce-refeyn commented 1 year ago

What do you mean by that?

I would have expected the uniform_filter over the integer array to essentially return a rounded version of the same filter over the floating point array. But it does differ, even in cases when the floating point returned an integer (like the 56 vs 55 in the bottom right corner or the arrays in the comment above). I'm sure there's a good reason for it, it was just a bit unexpected.