`view` array indexes that vary as a function of another index of that subarray?


#1

In a view is there a way to define an index of an array that varies as a function of another index of that array?
i.e. the indexing of dim 1 varies along dim 2.

I’d love to avoid making a memory copy of the subarray. My application goal is to define a view to a 2D region of a video that changes position in time (the 3rd dimension), and read the full 3 dimensional output without creating a memory copy of the subarray.

An illustrative simplistic part-pseudocode example, where the resulting assignment of 1’s in data follows a sinusoidal shape:

f(x) = round.(Int,5 .+ (sin.((x)/2)*5))

data = zeros(Int,11,10)

dim2 = :                             #pseudocode
subarray = view(data,f(dim2),dim2)   #pseudocode

subarray .= 1

data
--
11×10 Array{Int64,2}:
 0  0  0  0  0  0  0  0  1  1
 0  0  0  0  0  0  0  1  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  1  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  1  0  0  0  0
 1  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  0  0  0  0  0
 0  1  0  0  0  0  0  0  0  0
 0  0  1  1  0  0  0  0  0  0

#2

Perhaps something like:

# Maps `getindex` by `f` and then indexes into `A`
struct MappedArray{F, T, N, A <: AbstractArray} <: AbstractArray{T, N}
    f::F
    a::A
end
MappedArray(f::F, A::AbstractArray{T, N}) where {F, T, N} = MappedArray{F, T, N, typeof(A)}(f, A)
Base.size(A::MappedArray) = size(A.a)
Base.getindex(A::MappedArray, i::Integer) = A.a[A.f(i)]

a = 1:10
z = length(a)
m = MappedArray(x -> rand(1:z), a)

Where we have an array that randomly picks an index for a:

julia> m[1]
9

julia> m[1]
1

julia> m[1]
8

Now we can create a view with it:

julia> v = rand(100);

julia> c = view(v, m)
10-element view(::Array{Float64,1}, [4, 6, 9, 5, 10, 3, 9, 10, 10, 8]) with eltype Float64:
 0.8159868211487065
 0.5112093449671173
 0.4704287700375531
 0.6906170863362935
 0.6906170863362935
 0.5112093449671173
 0.5112093449671173
 0.15705852167037837
 0.4704287700375531
 0.4589219238937361

julia> c
10-element view(::Array{Float64,1}, [4, 2, 1, 7, 1, 4, 9, 8, 7, 7]) with eltype Float64:
 0.4589219238937361
 0.6779349956175409
 0.8159868211487065
 0.8159868211487065
 0.6906170863362935
 0.15705852167037837
 0.8159868211487065
 0.48428708295030143
 0.4704287700375531
 0.48428708295030143

so in theory you should be able to just come up with your f.


#3

Wow, that’s a really nice approach!
I’m trying to figure out how to make it handle multidimensional source arrays & views.


#4

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)