JuliaImages / ImageFiltering.jl

Julia implementations of multidimensional array convolution and nonlinear stencil operations
Other
99 stars 51 forks source link

Question: mapwindow with multi-threads? #239

Open mdhe1248 opened 2 years ago

mdhe1248 commented 2 years ago

Hi, I was wondering if mapwindow can work with multi-threads. Best,

rafaqz commented 1 year ago

If people here are interested (@johnnychen94 ?), I'm extracting the neighborhoods from DynamicGrids.jl into a KernelAbstractions.jl based neighborhoods/windowing package. https://github.com/rafaqz/Neighborhoods.jl

Its threaded using KernelAbstractions, and runs on GPU as well for #209.

There may be some issues of differing window implementations, and would be some things to add to support everything here. But it would be nice to have a shared back-end for these things.

Benchmarks seem pretty good, especially for smaller windows (on a six core machine, so mostly actually SVector vs SubArray window performance):

julia> @benchmark mapwindow(f2, r, (3, 3))
BenchmarkTools.Trial: 67 samples with 1 evaluation.
 Range (min … max):  67.756 ms … 82.823 ms  ┊ GC (min … max): 12.14% … 25.00%
 Time  (median):     75.261 ms              ┊ GC (median):    10.94%
 Time  (mean ± σ):   75.295 ms ±  4.377 ms  ┊ GC (mean ± σ):  14.06% ±  4.96%

  █ ▄      ▁              ▁ ▄█▄▄█▄▄       █         █    ▁  ▁
  █▆█▆▁▆▁▆▆█▁▆▁▁▁▆▁▁▁▁▁▁▆▁█▁███████▆▆▆▁▁▁▆█▁▆▁▁▆▁▁▆▁█▆▁▆▆█▆▆█ ▁
  67.8 ms         Histogram: frequency by time        82.4 ms <

 Memory estimate: 132.63 MiB, allocs estimate: 1055963.

julia> @benchmark broadcast_neighborhood(f, r; neighborhood=Window{1}()) # Radius 1 means diameter 3
BenchmarkTools.Trial: 1759 samples with 1 evaluation.
 Range (min … max):  2.281 ms …   4.462 ms  ┊ GC (min … max): 0.00% … 37.36%
 Time  (median):     2.757 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.836 ms ± 321.669 μs  ┊ GC (mean ± σ):  2.88% ±  7.64%

             ▂▅▇▇█▆▅▃▂                               ▁ ▁▁     ▁
  ▄▄▅▅▄▆▇▇▇▅▆██████████▇▇▆▅▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▅▅█████▆▄▆ █
  2.28 ms      Histogram: log(frequency) by time      4.13 ms <

 Memory estimate: 7.64 MiB, allocs estimate: 126.
mdhe1248 commented 1 year ago

Because I couldn't figure out how to apply Neighborhood.jl, I instead tried TiledIteration.jl (https://github.com/JuliaArrays/TiledIteration.jl). Somehow, it did not much improve the performance. Could anyone advise me?


using TiledIteration, Images, Statistics

function TileIterator_overlap(ax, tilesize, overlap)
  tiles = collect(TileIterator(ax, tilesize))
  tiles1 = similar(tiles)
  for (i,tile) in enumerate(tiles)
    tiles1[i] = Tuple(first(tile[j])-overlap:last(tile[j])+overlap for j in eachindex(ax))
  end
  tiles1
end

function tiled_mapwindow(ax, tilesize, overlap, f, a, ws)
  tileids = collect(TileIterator(ax, tilesize))
  tileids_overlap = TileIterator_overlap(ax, tilesize, overlap)
  ## mapwindow setting
  padded = padarray(a, Pad(:replicate, ntuple(x -> overlap, length(ws)))) # paddedview?
  b = zeros(eltype(padded), size(a)) #MMAP?
  ## run mapwnidow
  Threads.@threads for i = 1:length(tileids)
    tileaxes = tileids[i]
    tileaxes_overlap = tileids_overlap[i]
    mapwindow!(f, b[tileaxes...], padded[tileaxes_overlap...], ws, Inner())
  end
  b
end

## Run test
a = zeros(300,300,300)
a[3,3,3] = 1
ax = axes(a)
tilesize = (100,100,100)
overlap = 1
ws = (3,3,3)
@time tiled_mapwindow(ax, tilesize, overlap, mean, a, ws);
@time mapwindow(mean, a, ws);
rafaqz commented 1 year ago

Its Neighborhoods.jl, Neighborhood.jl is an unrelated package ;)

I need to think of a new name to be able to register it now.

Neighborhoods.jl will definately improve the performance. There is example code above.

rafaqz commented 1 year ago

Now registered: https://github.com/rafaqz/Stencils.jl