1.67× fast_blur speedup + gaussian blur improvements

[art049]

Hello,

I've been profiling the blur operations and have a PR ready with three targeted optimizations: art049#2

Wall-time results (cargo bench)

Benchmark	Before	After	Speedup
fast blur σ=3.0	40.1 ms	24.0 ms	1.67×
fast blur σ=7.0	40.1 ms	24.1 ms	1.66×
fast blur σ=50.0	40.8 ms	24.6 ms	1.65×
gaussian blur σ=3.0	21.6 ms	18.5 ms	+17%
gaussian blur σ=7.0	36.3 ms	33.1 ms	+10%
gaussian blur σ=50.0	361.4 ms	348.7 ms	+4%

I also have a separate PR setting up CodSpeed: art049#1. That's what I used to find these hot spots and iterate.
Full disclaimer: I created CodSpeed, and we'd be happy to provide it for free to image-rs, but I prefer to be upfront about it.
Setting it up would require installing the GitHub App on the org (one-click, org admin needed), and you'd get PR comments like this: art049#2 (comment)

Both PRs are independent, happy to open the perf PR against upstream on its own. Let me know how you'd like to proceed!

[awxkee]

I can't speak on behalf on maintainers but to me it seems those PRs a bit controversial compared to the approach I tried to implement in those blur algorithms and to optimization itself:

• It introduces unsafe code, which we specifically wanted to avoid.
• It moves toward modular arithmetic but doesn't use actual integer reciprocals. Instead, it uses integer division which kills throughput since division requires ~20 cycles and is a blocking operation.
• It changes to_u8 and to_u16 under assumption they are faster, but it holds only in two cases: target platform is not aarch64 (or any other having direct rounding instruction ) and code is not compiled against sse4.1 on x86 which have rounding instructions. Also a lot of branching isn't optimization as well. Real fast rounding should move your float to a range where float can be the only integer and move it back as so:

 const SHIFTER: f32 = ((1u32 << 23) + (1u32 << 22)) as f32;
 let rounded_nearest = ((self + f32::EPSILON) + SHIFTER) - SHIFTER

• It goes against the goals to make less code, not more code.

[art049]

It introduces unsafe code, which we specifically wanted to avoid.

Obviously, it depends on the policy, but here this was the most significant hot spot (17%), and the pattern seemed simple enough since the whole buffers are overwritten before being read. Let me know what's your policy regarding those.

It moves toward modular arithmetic but doesn't use actual integer reciprocals.

Yes, definitely happy to implement it.

It goes against the goals to make less code, not more code.

Definitely agree on the u32 accumulator implementation, I first wanted to see if this was something of interest, but it's probably possible to create an abstraction around it to avoid duplicating the logic.

For to_u8/u16, I'll check it out, thanks for pointing it out.

[197g]

The assumption is also very strange, size_of::<P>() == 1 && P::DEFAULT_MIN_VALUE == P::zero() does not imply u8 or at least no robustly in a manner that we would ever notice if we accidentally broke this by i8 for instance. Note that the Subpixel (edit: Primitive, sorry for the confusion) trait is sealed on the main branch so please use that customization point—it's exactly so we can do this safely. Casting through byte slices instead of manually doing unsafe is highly preferred. Deviating from that should come with a good explanation of how it would not work otherwise. 'Hot spot' is not good enough; the cast itself is not in a tight loop.

[art049]

I implemented your suggestions. Here are the new results on x86 without sse4.1:

	Benchmark	Base	Head	Change
⚡	fast blur: sigma 3.0	40.1 ms	7.4 ms	5.4x
⚡	fast blur: sigma 7.0	40.1 ms	7.5 ms	5.4x
⚡	fast blur: sigma 50.0	40.8 ms	8.2 ms	5.0x
⚡	fast_blur (standalone)	40.9 ms	8.3 ms	5.0x
⚡	gaussian blur: sigma 3.0	22.2 ms	18.4 ms	+21%
⚡	gaussian blur: sigma 7.0	37.9 ms	34.1 ms	+11%
✅	gaussian blur: sigma 50.0	369.7 ms	367.1 ms	unchanged

(full report)

Integer reciprocal made the u32 accumulator way faster; that's really cool.

For the rounding in to_u8/to_u16, I implemented it as you mentioned when it's not natively supported (aarch64/sse4.1). Some other platforms have such instructions (at least wasm), happy to include them as well if you think it's worth it.

Thanks for those tricks @awxkee!

I refactored the accumulator logic to avoid duplication, as you mentioned. This gives flexibility on the accumulator and weight types.

Finally, after your comment @197g, I dropped the unsafe part that avoided initialization of buffers in the Gaussian blur. I'd still love to find a way to drop in a way that would match the safety requirements of the project, but I'm out of options right now.

If you find it better, happy to open the PR against the upstream for a proper review.

[art049]

After updating and taking most of the feedback it's actually a 5.5x speedup

Here is the full report with all the results

Walltime

Benchmark	Baseline	Optimized	Speedup
fast blur: sigma 3.0	41.8 ms	7.7 ms	×5.5
fast blur: sigma 7.0	42.1 ms	7.7 ms	×5.4
fast blur: sigma 50.0	42.5 ms	8.4 ms	×5.0
fast_blur	41.9 ms	8.5 ms	×4.9
gaussian blur: sigma 3.0	24.1 ms	20.6 ms	+17%
gaussian blur: sigma 7.0	39.6 ms	37.5 ms	+5% (noise)
gaussian blur: sigma 50.0	364.1 ms	379.1 ms	-4% (noise)

Simulation

Benchmark	Baseline	Optimized	Speedup
fast blur: sigma 3.0	254.3 ms	45.3 ms	×5.6
fast blur: sigma 7.0	254.8 ms	45.6 ms	×5.6
fast blur: sigma 50.0	260.6 ms	50.3 ms	×5.2
fast_blur	260.6 ms	50.3 ms	×5.2
gaussian blur: sigma 3.0	159.8 ms	137.1 ms	+17%
gaussian blur: sigma 7.0	270.1 ms	247.4 ms	+8%
gaussian blur: sigma 50.0	1.38 s	1.36 s	unchanged

[197g]

If you find it better, happy to open the PR against the upstream for a proper review.

Yes please. It looks quite nice already, and thank you for explaining the numbers. Then I have to say the simulation sounds quite compelling with that good of a match.