A “low-tech” but relatively fast and simple approach may be to consider explicit diffusion as smoothing mechanism. Let’s say you want to smooth a field A, then following would do the job in 2D (and can be easily extended to n-D, GPU parallel processing, …):
using Plots
@views function do_smoothing()
nsmooth = 50
A = rand(100,100)
CFL = 0.5/2.1/length(size(A))
for ismooth = 1:nsmooth
A[2:end-1,2:end-1] .= A[2:end-1,2:end-1] + CFL*(diff(diff(A[:,2:end-1],dims=1),dims=1) + diff(diff(A[2:end-1,:],dims=2),dims=2))
A[1,:] = A[2,:]; A[end,:] = A[end-1,:]; A[:,1] = A[:,2]; A[:,end] = A[:,end-1]
end
display(heatmap(A'))
end
@time do_smoothing()
Then you could as well decouple the diffusion in different dimensions and apply various number of smoothing steps individually for different dimensions.
EDIT: To put this in context with the moving mean filter, it would e.g. take about 10 diffusion smoothing steps nsmooth=10 to achieve a ~10-cell window smoothing:
modifying the above code as such
using Plots
@views function do_smoothing()
nsmooth = 10
A0 = 0.0*rand(100,100)
A0[40:60,40:60] .= 1.0
A = copy(A0)
CFL = 0.5/2.1/length(size(A))
for ismooth = 1:nsmooth
A[2:end-1,2:end-1] .= A[2:end-1,2:end-1] + CFL*(diff(diff(A[:,2:end-1],dims=1),dims=1) + diff(diff(A[2:end-1,:],dims=2),dims=2))
A[1,:] = A[2,:]; A[end,:] = A[end-1,:]; A[:,1] = A[:,2]; A[:,end] = A[:,end-1]
end
p=plot(A0[:,50],label="initial", linewidth=3); plot!(A[:,50], linewidth=3, markersize=5, markershape=:circle, label="smooth", framestyle=:box)
display(p)
end
@time do_smoothing()
