Taking views of OffsetArrays

[juliohm]

My package is making use of the following call:

using Images

D = imfilter(ones(10,10), centered(ones(3,3)))
view(D,1:2,:) # ERROR: indices go from 2 to 9

I am not very familiar with the concept of offset arrays, but I am using views everywhere in my code to save memory allocations: https://github.com/juliohm/ImageQuilting.jl/blob/master/src/iqsim.jl#L258-L259

In the package, the D can either be a normal Array if I use "symmetric" padding or an OffsetArray if I use Inner() padding, both defined in ImageFiltering.jl. Can you help solve this case-by-case issue, or how to convert OffsetArrays to normal Arrays?

[timholy]

First, to clarify for anyone else who comes and reads this issue, the test script is misleading because as written it doesn't throw an error:

julia> using Images

julia> D = imfilter(ones(10,10), centered(ones(3,3)))
10×10 Array{Float64,2}:
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0

julia> view(D,1:2,:) # what error?
2×10 SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},Colon},false}:
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0

But the view fails if you add Inner() as padding:

julia> D = imfilter(ones(10,10), centered(ones(3,3)), Inner())
OffsetArrays.OffsetArray{Float64,2,Array{Float64,2}} with indices 2:9×2:9:
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0
 9.0  9.0  9.0  9.0  9.0  9.0  9.0  9.0

julia> view(D,1:2,:)
ERROR: BoundsError: attempt to access OffsetArrays.OffsetArray{Float64,2,Array{Float64,2}} with indices 2:9×2:9 at index [1:2,Colon()]
 in throw_boundserror(::OffsetArrays.OffsetArray{Float64,2,Array{Float64,2}}, ::Tuple{UnitRange{Int64},Colon}) at ./abstractarray.jl:363
 in checkbounds at ./abstractarray.jl:292 [inlined]
 in view(::OffsetArrays.OffsetArray{Float64,2,Array{Float64,2}}, ::UnitRange{Int64}, ::Colon) at ./subarray.jl:64

If you've trimmed off the edges by using Inner, then it's nice to know what region the output refers to. OffsetArrays allow you to specify that locations have meaning. Row "1" doesn't exist:

julia> indices(D)
(2:9,2:9)

Since the top row is called "2," just replace your line with view(D, 2:3, :).

[timholy]

There's a display bug, though, if you try to show the results of that view. Ref JuliaLang/julia#20431.

[juliohm]

Sorry for the typo in the example, I missed the Inner() option.

I argue that this manual offset handling is gonna cause more confusion and scary bugs than it will help users. I did read the concept of OffsetArrays, and it is super cool, never played with Fortran myself, but it can be a really useful trick. However, I don't think returning OffsetArrays sometimes and Arrays other times is productive. Indexing those 3D images with padding, centered vs. not centered kernels, etc, is already super trick to get right, even more if we decide to handle these offsets manually when taking views.

How can we avoid this manual check, is there a way to convert OffsetArray for plain simple Arrays? I don't want to play in this territory, the current interface is forcing me to do it though. I understand you are the master of arrays and N-dimensional indexing, but I am a normal human that spends hours debugging code to get indexing right 😊

[timholy]

parent(A).

But are you sure you're actually evaluating the proposal on its merits? Let's first think about the "traditional" scheme. If A is the input and B is the output, how do you match positions in A relative to positions in B? If A has the same size as B, maybe it's easy, but if B is smaller than A then you might be left wondering. B[1,1] always means "the first well-defined output," whatever that means, and B[2,3] is relative to that. Suppose your kernel has width 4, and you wonder where zero (i.e., the origin for the kernel's coordinates) is in there---at element 2? element 3? element 1?

So, to address this problem, there's a well-established tradition in computer science (and other fields) of reducing confusion by choosing consistent names for things. For example, every person in the US gets a social security number that stays invariant even if you change your first or last name---you can use the social security number to "link" what might otherwise be thought of as 2 different people. In the new scheme, every array that "descends" from a parent gets indexed by an absolute position: A[2,3] always refers to the value of A a particular position in space called (2,3), and it means exactly the same position as B[2,3]. If your kernel is of length 4, you specify where the origin (0) is; the centered function is just a simple utility to make common cases easy. Once you've specified the origin of the kernel, everything else is determined. For example:

julia> A = 1:8
1:8

julia> kern = OffsetArray([1], 0:0)   # a "delta function" centered at 0
OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}} with indices 0:0:
 1

julia> imfilter(A, kern, Fill(0))
8-element Array{Int128,1}:
 1
 2
 3
 4
 5
 6
 7
 8

julia> kern = OffsetArray([1], -2:-2)   # a "delta function" centered at -2
OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}} with indices -2:-2:
 1

julia> imfilter(A, kern, Fill(0))
8-element Array{Int128,1}:
 0
 0
 1
 2
 3
 4
 5
 6

Isn't that amazingly clear? Just like in the mathematics textbooks.

I see your concern about the views. But let me point this out:

julia> using Images

julia> D = imfilter(ones(10,10), centered(ones(3,3)), Inner());

julia> indices(view(D, 2:3, :))
(Base.OneTo(2),2:9)

julia> using OffsetArrays

julia> inds = OffsetArray(2:3, 2:3)
OffsetArrays.OffsetArray{Int64,1,UnitRange{Int64}} with indices 2:3:
 2
 3

julia> indices(view(D, inds, :))
(2:3,2:9)

This way of taking views keeps the names constant.

Let's separate the technical from the human issues: (1) have I done an adequate job of educating you/others about how this works? (2) from a technical standpoint, is it less confusing overall? (3) will users need to adapt their code and habits.

If there are good arguments to the contrary, I'm happy to consider changing it. But I don't want (1) and (3) contaminating the answer to (2). I suspect the answers to 1-3 are "no," "yes," and "certainly," but I'm happy to have you try to persuade me otherwise. There is some documentation at http://docs.julialang.org/en/latest/devdocs/offset-arrays.html#Arrays-with-custom-indices-1, but not a demo of why this is useful.

[juliohm]

Hi Tim, can you please elaborate on the last example you gave above? I am still trying to understand the message you tried to teach. I agree there is a lot of power in general indexing, but a lot of complexity too. Even for me, someone that is playing with scientific programming languages with different indexing schemes, I find it super tricky. I need more education to get it right.

[timholy]

If you're in a hurry, a WIP blog post is here. If you have further questions, please ask them as specific comments over at https://github.com/JuliaLang/julialang.github.com/pull/498/files.