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:
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:
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:
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
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
