# Best approach to smooth array? None of my approaches are very useful
**URL:** https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633
**Category:** Performance
**Tags:** images, statistics, smoothing
**Created:** [November 17, 2018, 5:33am UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633 "2018-11-17T05:33:30Z")
**Posts on this page:** 19
**Page:** 1
### Author:  [@aramirezreyes](https://discourse.julialang.org/u/aramirezreyes)
#### Post date: [November 17, 2018, 5:33am UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/1 "2018-11-17T05:33:30Z")
Hi!
I´ve been trying this week to get something accomplished but I cannot get it done.
Problem: I have group of variables coming from a NetCDF. When imported using NCDatasets, I get arrays of 64 bits with a size (128,128,80,2400).
A code I was using (in matlab) has a line of the style `movemean(moveman(movmean(array,window_i),window_j),window_l)`
Which performs a smoothing in the 1st, 2nd, and 4th dimensions.
To do this in Julia, I have tried different approaches:
## First approach: using Images.ImageFiltering with the kernel:
```julia-auto
function kernel4d(window_h,window_t)
kern = ones(window_h,window_h,1,window_t)/(window_t*window_h*window_h);
end
```
And passing it to imfilter like this:
```julia-auto
function filter_array(array,smooth_x,smooth_time)
filtered = similar(array)
imfilter!(filtered,array, kernel4d(smooth_x,smooth_time))
return filtered
end
```
## Second approach: using Images.ImageFiltering with mapwindow:
```julia-auto
function filter_array2(array,smooth_x,smooth_time)
filtered = mapwindow(mean,mapwindow(mean,mapwindow(mean,array, (smooth_x,1,1,1)),(1,smooth_x,1,1)),(1,1,1,smooth_time))
end
```
## Third approach: using Satistics.mean
```julia-auto
function filter_array3(array,windowh,windowt)
filtered = similar(array)
fach = div((windowh + 1),2)
fact = div((windowt + 1),2)
xend,yend,zend,tend = size(array);
for I in CartesianIndices(array)
filtered[I] = mean(array[max(1,I[1]-fach):min(xend,I[1]+fach),I[2],I[3],I[4]]);
end
for I in CartesianIndices(array)
filtered[I] = mean(filtered[[1],max(1,I[2]-fach):min(yend,I[2]+fach),I[3],I[4]]);
end
for I in CartesianIndices(array)
filtered[I] = mean(filtered[I[1],I[2],I[3],max(1,I[4]-fact):min(tend,I[4]+fact)]);
end
return filtered
end
```
Although these approaches are not strictly equivalent (some take into account some padding in the borders), they are significantly slower than the matlab one line code (10x). My code has two major bottlenecks file reading (I don’t really know how to improve this) and this smoothing function. The window sizes are small in the first and second dimensions (around 5 elements) but for the fourth it is around 30 elements.
The first approach also uses impressive amounts of memory (it breaks with an out of memory error in a system with 50Gb).
The second approach has a problem: the second and the third loops are increasingly expensive due to accesing elements very far in memory.
I must say I call these process 8 times in a row for different variables. And I am using Julia 0.7.0
Is there something you could tell me about this?
Thanks
---
### Author:  [@mauro3](https://discourse.julialang.org/u/mauro3)
#### Post date: [November 17, 2018, 9:04am UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/2 "2018-11-17T09:04:29Z")
Have you tried factored kernels: [https://juliaimages.github.io/latest/imagefiltering.html#Factored-kernels-1](https://juliaimages.github.io/latest/imagefiltering.html#Factored-kernels-1)? Also, potentially it could be worth moving the last dim to be the first as that would lead to better memory locality.
---
### Author:  [@spaceLem](https://discourse.julialang.org/u/spaceLem)
#### Post date: [November 17, 2018, 9:42am UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/3 "2018-11-17T09:42:08Z")
Is approach 3 doing the right thing? You’ve got 3 `for` loops, each overwriting the previous one. Also in the second loop you have `[1]`, instead of `I[1]`, which might be doing something weird (not sure, don’t have Julia in front of me to test).
---
### Author:  [@aramirezreyes](https://discourse.julialang.org/u/aramirezreyes)
#### Post date: [November 18, 2018, 1:27am UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/4 "2018-11-18T01:27:25Z")
Thanks for noticing that. Indeed the result was wrong but fixing it didn’t help much. The thing about the three loops is what I want, as I want to filter the same array in three of its dimensions.
---
### Author:  [@aramirezreyes](https://discourse.julialang.org/u/aramirezreyes)
#### Post date: [November 18, 2018, 1:56am UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/5 "2018-11-18T01:56:24Z")
Hey mauro3, thanks for your answer:
I’ve tried Factored Kernels in this way, I am not so sure if this is the correct one, though.
```julia
function kernel4d_2(window_h,window_t)
kern1 = centered(ones(window_h,1,1,1)/(window_h));
kern2 = centered(ones(1,window_h,1,1)/(window_h));
kern3 = centered(ones(1,1,1,window_t)/(window_t));
kernz = centered(ones(1,1,1,1))
return kernelfactors((kern1,kern2,kernz,kern3))
end
```
Is this the correct way to use it, and should I expect it to be better than using the 4-d Kernel directly?
---
### Author:  [@y4lu](https://discourse.julialang.org/u/y4lu)
#### Post date: [November 18, 2018, 2:14am UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/6 "2018-11-18T02:14:37Z")
If it does any decimation to the output afterwards, it would probably be worth trying to build that in to the filter. The x2400 by 30 dimension seems like a good candidate for a stride of 2 or more, but it depends on how you want to use the output
Add: there maybe something in BlasPack or the like that can help?
---
### Author:  [@tim.holy](https://discourse.julialang.org/u/tim.holy)
#### Post date: [November 18, 2018, 11:33am UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/7 "2018-11-18T11:33:08Z")
You can use it that way but this also works:
```julia
julia> function kernel4d(window_h, window_t)
kernel_h = ones(window_h)/window_h
kernel_t = ones(window_t)/window_t
return kernelfactors((kernel_h, kernel_h, [1.0], kernel_t))
end
```
(see `?kernelfactors` and the part about being able to pass vectors). To explain what you’re seeing, try `?KernelFactors.ReshapedOneD` and look at the individual elements of the return value.
> and should I expect it to be better than using the 4-d Kernel directly?
That’s explained in the docs that @mauro3 linked:
> If the kernel is of size `m×n` , then the upper version line requires `mn` operations for each point of `filtered` , whereas the lower version requires `m+n` operations. Especially when `m` and `n` are larger, this can result in a substantial savings.
In your case it’s the difference between 5 \* 5 \* 30 = 750 and 5 + 5 + 30 = 40. This therefore explains the factor of 10 difference with Matlab: notice that your Matlab code is using a factored kernel: `window_i`, `window_j`, and `window_l` are separate factors of the kernel. So you were asking Julia to use a fundamentally worse (more computationally demanding) approach than you were asking Matlab.
For even greater speed see [https://juliaimages.github.io/latest/imagefiltering.html#Multithreading-1](https://juliaimages.github.io/latest/imagefiltering.html#Multithreading-1).
---
### Author:  [@aramirezreyes](https://discourse.julialang.org/u/aramirezreyes)
#### Post date: [November 18, 2018, 6:13pm UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/8 "2018-11-18T18:13:22Z")
Thanks for the suggestion! I’m gonna try it making that if everything else fails. Initially all the data is kept.
---
### Author:  [@aramirezreyes](https://discourse.julialang.org/u/aramirezreyes)
#### Post date: [November 21, 2018, 5:34am UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/9 "2018-11-21T05:34:08Z")
Thanks a lot Tim!
The suggested new form of the call to kernel factors causes a surprisingly long error message. I will stay with the old one until I can figure out what is happening.
I let you the message in case is interesting/useful. I am using julia 0.7, would it be useful to start an issue? This only happens apparently when using “export JULIA\_NUM\_THREADS=N” where N is greater than 1
The error line message too long, so I pasted it in pastebin:
> **[using ImageFiltering function kernel4d\_2(window\_h,window\_t) kernel\_h = o...](https://pastebin.com/E3V6CvJH)**
>
> Pastebin.com is the number one paste tool since 2002. Pastebin is a website where you can store text online for a set period of time.
---
### Author:  [@tim.holy](https://discourse.julialang.org/u/tim.holy)
#### Post date: [November 21, 2018, 8:25am UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/10 "2018-11-21T08:25:12Z")
As I think you realized, that’s a Julia bug, not an ImageFiltering bug. For me it works on Julia 1.0.2 (and I can reproduce the error on Julia 0.7), so perhaps it might be worth trying on a newer version of Julia. If you can reproduce it on the latest release, it’s definitely worth reporting as a bug.
---
### Author:  [@aramirezreyes](https://discourse.julialang.org/u/aramirezreyes)
#### Post date: [November 21, 2018, 8:31am UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/11 "2018-11-21T08:31:20Z")
I just tried Julia 1.0.2 and works flawlessly.
Thanks so much for all the help!
---
### Author:  [@aramirezreyes](https://discourse.julialang.org/u/aramirezreyes)
#### Post date: [January 30, 2019, 12:39am UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/12 "2019-01-30T00:39:45Z")
Hi Tim:
Are this two things supposed to give the same output?
```julia
a = rand(10,10,10)
b = similar(a)
for i in 1:10
b[:,:,i] = imfilter(a[:,:,i], kernelfactors(( ones(3)/3, ones(3)/3)))
end
c = imfilter(a,kernelfactors(( ones(3)/3 , ones(3)/3) , [1.0])) )
```
The question is,
`shouldn't b==c?`
It doesn’t, but I don’t understand why.
---
### Author:  [@tim.holy](https://discourse.julialang.org/u/tim.holy)
#### Post date: [January 30, 2019, 2:44pm UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/13 "2019-01-30T14:44:02Z")
It’s a bug, I filed an issue for you: [Bug in 3d cascade · Issue #92 · JuliaImages/ImageFiltering.jl · GitHub](https://github.com/JuliaImages/ImageFiltering.jl/issues/92)
---
### Author:  [@aramirezreyes](https://discourse.julialang.org/u/aramirezreyes)
#### Post date: [February 1, 2019, 10:15pm UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/14 "2019-02-01T22:15:39Z")
Just updating this in case someone finds this post and uses it as a reference.
> [@](#):
>
> You can use it that way but this also works:
>
> ```julia
> julia> function kernel4d(window_h, window_t)
> kernel_h = ones(window_h)/window_h
> kernel_t = ones(window_t)/window_t
> return kernelfactors((kernel_h, kernel_h, [1.0], kernel_t))
> end
>
> ```
The correct form of smoothing in the three directions without displacing the array is:
```julia
function kernel4d(window_h, window_t)
kernel_h = ones(window_h)/window_h
kernel_t = ones(window_t)/window_t
return kernelfactors((centered(kernel_h), centered(kernel_h), centered([1.0]), centered(kernel_t)))
end
```
---
### Author:  [@tim.holy](https://discourse.julialang.org/u/tim.holy)
#### Post date: [February 1, 2019, 10:18pm UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/15 "2019-02-01T22:18:38Z")
> The correct form of smoothing in the three directions without displacing the array is
It’s not just three directions, you want to use `centered` even in 1d or 2d if you are using a kernel that is symmetric around 0. The original form shifted the array along all 3 axes, not just the 3rd axis.
(The behavior above turned out not to be a bug, it was a feature.)
---
### Author:  [@aramirezreyes](https://discourse.julialang.org/u/aramirezreyes)
#### Post date: [February 1, 2019, 10:50pm UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/16 "2019-02-01T22:50:11Z")
Right, I was just trying to make it a final answer for the initial question of the post 😃
Thanks!
---
### Author:  [@aramirezreyes](https://discourse.julialang.org/u/aramirezreyes)
#### Post date: [May 5, 2021, 7:16pm UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/17 "2021-05-05T19:16:30Z")
I am coming back to this issue two year later:
Smoothing is a part of one calculation I do in the framework I am using for my PhD research.
The thing here is that the arrays I work with are of size (512,512,80,1200) of Float32s.
In one part of the process I smooth my data out and this is without question the part that takes the longest in the orders of hours while the rest of the process is done in less than one hour.
I am thinking an important part of this issue is the data non-locality. Especifically, to smooth in the time direction, julia needs to fetch 50 values that are 512_512_80 floats away from each other, and it needs to repeat this process for 512_512_80.
I know I can probably live with this and let it sink in the sands of time after my PhD, but I am obsessed with it so I come back every once in a while.
I ask you, dear Julia community. What would be the best approach to smooth a 512_512_80\*1200 array using moving means of 20 points in the first dimension, 20 in the second and 50 in the fourth?
I know intelligent people must have come across this in some moment in time.
---
### Author:  [@luraess](https://discourse.julialang.org/u/luraess)
#### Post date: [May 6, 2021, 12:54pm UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/18 "2021-05-06T12:54:08Z")
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, …):
```julia
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
```julia
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()
```
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [May 6, 2021, 10:23pm UTC](https://discourse.julialang.org/t/best-approach-to-smooth-array-none-of-my-approaches-are-very-useful/17633/19 "2021-05-06T22:23:23Z")
@aramirezreyes, have you looked at [FastConv.jl](https://github.com/aamini/FastConv.jl)?
The package seems to be broken but looking at the code and description they seem to implement the same idea of performing direct convolution with small kernels as in the code posted by @sijo [here](https://discourse.julialang.org/t/seven-lines-of-julia-examples-sought/50416/48).
From your posts above it seems that you had tried something similar but if you could please confirm, it would be helpful.