Here is my attempt, with G’MIC.
The command foo below does the same kind of box blurring with box sizes that increases horizontally and vertically along the X and Y axes.
# $1 = min radius, $2 = max radius.
foo : skip ${1=1},${2=max(w,h)/10}
repeat $! l[$>]
100%,100%,1,2,"round([ lerp($1,$2,(x/w)^3), lerp($1,$2,(y/h)^3) ])"
+cumulate[0] xy
f[0] "*
const W = w - 1; const H = h - 1;
px(X,Y) = I(#2,X,Y,0,c);
T(X,Y)=Y>1?px(X,Y-1):Y*px(X,0);
B(X,Y)=Y<H?px(X,Y):(1+Y-H)*px(X,H)-(Y-H)*px(X,H-1);
G(X,Y)=X<W?B(X,Y):(1+X-W)*B(W,Y)-(X-W)*B(W-1,Y);
D(X,Y)=X>1?T(X-1,Y):X*T(0,Y);
E(X,Y)=X<W?T(X,Y):(1+X-W)*T(W,Y)-(X-W)*T(W-1,Y);
F(X,Y)=X>1?B(X-1,Y):X*B(0,Y);
bw = i(#1,x,y,0,0); # Horizontal size of blur window
bh = i(#1,x,y,0,1); # Vertical size of blur window
bw1 = int(bw/2); bh1 = int(bh/2); # Left/Top half-sizes
bw0 = bw - bw1 - 1; bh0 = bh - bh1 - 1; # Right/Bottom half-sizes
x0 = x - bw0;
y0 = y - bh0;
x1 = x + bw1;
y1 = y + bh1;
(G(x1,y1) + D(x0,y0) - E(x1,y0) - F(x0,y1))/(bw*bh)"
k[0]
endl done
It can be used like this:
$ gmic sp colorful tic +foo , toc
[gmic]-0./ Start G'MIC interpreter.
[gmic]-1./ Input sample image 'colorful' (1 image 800x800x1x3).
[gmic]-1./ Initialize timer.
[gmic]-2./ Elapsed time: 0.093 s.
For a large image (8000x8000) :
$ gmic sp colorful,8000 tic +foo , toc
[gmic]-0./ Start G'MIC interpreter.
[gmic]-1./ Input sample image 'colorful' (1 image 8000x8000x1x3).
[gmic]-1./ Initialize timer.
[gmic]-2./ Elapsed time: 8.538 s.
[gmic]-2./ Display images [0,1] = 'colorful, colorful_c1'.
Timings look quite good to me.
EDIT : edited to better manage boundary conditions (by extending the image with Neumann conditions before processing).
EDIT2: Used @garagecoder trick to manage boundary conditions of SAT. Very efficient stuff !
. Someday but never on my decrepit laptop.)
Anyway, that is the cost of participating in my threads. Mine, including PlayRaws, typically see tumbleweeds.