The more I think about the nature of indexing required for this, I’m realising that for my particular application of subregion video registration without creating memory copies, the below may actually be an acceptable approach.
I’ll need to explore the overhead of interacting with the array of 2D views vs. a 3D array, but it allows the 2D subregions to be shifted for registration.
One question though is what the best way to preallocate the subdata_arrayofviews
object as I’m only defining the number of array SubArray
elements, without preallocating the SubArray
itself. Perhaps thats insignificant though.
using BenchmarkTools
data = rand(UInt8,1000,1000,1000)
x_sub = 100:200
y_sub = 200:400
function method1(data,x_sub,y_sub)
subdata_singleview = view(data,x_sub,y_sub,:) #Standard non-translating view method.
end
function method2(data,x_sub,y_sub)
subdata_arrayofviews = Array{SubArray}(undef,1000)
for i = 1:size(data,3)
x_sub_jitter = x_sub .+ rand(-10:10)
y_sub_jitter = y_sub .+ rand(-10:10)
subdata_arrayofviews[i] = view(data,x_sub_jitter,y_sub_jitter,i)
end
return subdata_arrayofviews
end
@btime m1 = method1(data,x_sub,y_sub) #688.592 ns (6 allocations: 256 bytes)
@btime m2 = method2(data,x_sub,y_sub) #369.149 μs (7490 allocations: 281.20 KiB)
–
I thought to try building an indexing vector based on collating the translated subregion indices for each 2D slice, then reshaping the view based on this discussion: https://github.com/JuliaLang/julia/issues/4211
However (and obviously… in hindsight) it seems that the indexing vector is sufficiently large to do away with any of the benefits of using a view, given each pixel index occupies an Int64.
If there was a way to represent this vector through a series of UnitRange
perhaps it might save a decent amount of memory/allocations.
function method3(data,x_sub,y_sub)
w = size(data,1)
h = size(data,2)
npix::UInt64 = w*h
indexplane = reshape(1:npix,w,h)
indexvec = Vector{Int64}(undef,0)
for i = 1:size(data,3)
x_sub_jitter = x_sub .+ rand(-10:10)
y_sub_jitter = y_sub .+ rand(-10:10)
append!(indexvec,vec(indexplane[x_sub_jitter,y_sub_jitter] .+ ((i-1)*npix)))
end
subdata_reshaped = reshape(view(data,indexvec),length(x_sub),length(y_sub),size(data,3))
end
@btime m3 = method3(data,x_sub,y_sub) #490.177 ms (6018 allocations: 350.78 MiB)