Arrays are a crucial component of any programming language, particularly for a data-oriented language like Julia. Arrays store values according to their location: in Julia, given a two-dimensional array `A`

, the expression `A[1,3]`

returns the value stored at a location known as `(1,3)`

. If, for example, `A`

stores `Float64`

numbers, the value returned by this expression will be a single `Float64`

number.

Julia's arrays conventionally start numbering their axes with 1, meaning that the first element of a one-dimensional array `a`

is `a[1]`

. The choice of 1 vs. 0 seems to generate a certain amount of discussion. A fairly recent addition to the Julia landscape is the ability to define arrays that start with an *arbitrary* index. The purpose of this blog post is to describe why this might be interesting. This is really a "user-oriented" blog post, hinting at some of the ways this new feature can make your life simpler. For developers who want to write code that supports arrays with arbitrary indices, see this documentation page.

Sometimes arrays are used as simple lists, in which case the indices may not matter to you. But in other cases, arrays are used as a discrete representation of a continuous quantity (e.g., values defined over space or time), and in such cases the array indices correspond to "location" in a way that may be meaningful.

As a simple example, consider the process of rotating an image:

img | img_rotated |
---|---|

Many languages provide functions for rotating an image; in Julia, you can do this with the `warp`

function defined in ImageTransformations.

Things get a little more "interesting" when you want to compare pixels in the rotated image to those of the original image. How, exactly, do these pixels match up? In other words, for a location `img[i,j]`

, what is the corresponding `i′,j′`

location in `img_rotated`

? In many languages, figuring out these types of geometric alignments may not be a simple task; it's no exaggeration to say that in complex situations (e.g., a three-dimensional image with a complex spatial deformation) that one can spend a day or more figuring out exactly how pixels/voxels in two images are supposed to be compared.

Why is this such a hard problem? The core problem is that, in most cases, the language is essentially "lying" to you about the location of pixels: if arrays always start indexing at 1 along any given axis, the array indices don't really correspond to an absolute spatial location. An index of 1 means "first index" rather than "spatial location 1."

So to fix this in Julia, starting with version 0.5 we support arrays with indices that don't start with 1. Let's illustrate this by specifying that we want the above rotation to be around a point in the head of the cameraman. Let's load the image:

```
julia> using Images, TestImages
julia> img = testimage("cameraman");
julia> summary(img)
"512×512 Array{Gray{N0f8},2}"
```

`summary`

shows that `img`

is a grayscale image indexed over the ranges `1:512×1:512`

. Using any of several approaches (e.g., ImageView and paying attention to the status bar to get the mouse pointer location), we can learn that the point `(y=126, x=251)`

is in the head of the cameraman. Consequently, let's define a rotation around this point:

```
julia> using Rotations, CoordinateTransformations
julia> tfm = Translation(125,250) ∘ LinearMap(RotMatrix(pi/6)) ∘ Translation(-125,-250)
AffineMap([0.866025 -0.5; 0.5 0.866025], [141.747, -29.0064])
```

This defines `tfm`

as the composition of a translation (shifting the head to the origin) followed by a rotation, and then translating back. (You can get the composition operator `∘`

by typing `\circ`

and then hitting TAB.) If we apply this transformation to the image, we get an interesting result:

```
julia> img_rotated = warp(img, tfm);
julia> summary(img_rotated)
"-107:592×-160:539 OffsetArray{Gray{N0f8},2}"
```

Perhaps surprisingly, `img_rotated`

is indexed over the range `-107:592×-160:539`

, meaning that we access the upper left corner by `img_rotated[-107,-160]`

and the lower right corner by `img_rotated[592,539]`

. It's not hard to see why these numbers arise, if we see how the corners of `img`

are transformed by `tfm`

:

```
julia> using StaticArrays
julia> itfm = inv(tfm)
AffineMap([0.866025 0.5; -0.5 0.866025], [-108.253, 95.9936])
julia> itfm(SVector(1,1))
2-element SVector{2,Float64}:
-106.887
96.3597
julia> itfm(SVector(512,1))
2-element SVector{2,Float64}:
335.652
-159.14
julia> itfm(SVector(1,512))
2-element SVector{2,Float64}:
148.613
538.899
julia> itfm(SVector(512,512))
2-element SVector{2,Float64}:
591.152
283.399
```

This makes it apparent that the output's indices span the region of the transformed coordinates.

The fact that the output preserves the coordinates makes it trivial to compare the images:

`julia> cv = colorview(RGB, paddedviews(0, img, img_rotated, img)...)`

`paddedviews(0, arrays...)`

pads input arrays with 0, as needed, to give them all the same indices, and `colorview(RGB, r, g, b)`

inserts the grayscale images `r`

, `g`

, and `b`

into the red, green, and blue channels respectively. If we visualize `cv`

, we see the following:

around image center | around head (cv) |
---|---|

The image on the left is for reference, showing what a rotation around the image center would look like when properly aligned. The image on the right corresponds to the steps taken above, and indeed confirms that the rotation is around the head. Alternatively, we can focus on the overlapping portions of these images like this:

```
julia> inds = map(intersect, indices(img), indices(img_rotated))
(1:512, 1:512)
julia> imgi = img[inds...];
julia> imgri = img_rotated[inds...];
```

so that `colorview(RGB, imgi, imgri, imgi)`

displays as

Since the indices of the pixels encode absolute spatial location, it's trivial to keep track of how different pixels align: pixel `i,j`

in one image corresponds to pixel `i,j`

in the other. This is true even if our coordinate transformation were far more complicated than a simple rotation.

Having motivated why this might be useful, let's take a step back and walk through array indices a bit more systematically.

In Julia, if we define an array

```
julia> A = collect(reshape(1:30, 5, 6))
5×6 Array{Int64,2}:
1 6 11 16 21 26
2 7 12 17 22 27
3 8 13 18 23 28
4 9 14 19 24 29
5 10 15 20 25 30
```

then we can refer to a rectangular region like this:

```
julia> B = A[1:3, 1:4]
3×4 Array{Int64,2}:
1 6 11 16
2 7 12 17
3 8 13 18
```

For certain applications, one negative to extracting blocks is that there is no record indicating where the new block originated from:

```
julia> B2 = A[2:4, 1:4]
3×4 Array{Int64,2}:
2 7 12 17
3 8 13 18
4 9 14 19
julia> B2[1,1]
2
```

So `B2[1,1]`

corresponds to `A[2,1]`

, despite the fact that, as measured by their indices, these are not the same location.

To maintain consistent "naming" of our indices, let's use the OffsetArrays package:

```
julia> using OffsetArrays
julia> B3 = OffsetArray(A[2:4, 1:4], 2:4, 1:4)
# wrap the snipped-out piece in an OffsetArray
OffsetArrays.OffsetArray{Int64,2,Array{Int64,2}} with indices 2:4×1:4:
2 7 12 17
3 8 13 18
4 9 14 19
julia> B3[3,4]
18
julia> A[3,4]
18
```

So the indices in `B3`

match those of `A`

. Indeed, `B3`

doesn't even have an element "named" `(1,1)`

:

```
julia> B3[1,1]
ERROR: BoundsError: attempt to access OffsetArrays.OffsetArray{Int64,2,Array{Int64,2}}
with indices 2:4×1:4 at index [1, 1]
Stacktrace:
[1] throw_boundserror(::OffsetArrays.OffsetArray{Int64,2,Array{Int64,2}}, ::Tuple{Int64,Int64}) at ./abstractarray.jl:426
[2] checkbounds at ./abstractarray.jl:355 [inlined]
[3] getindex(::OffsetArrays.OffsetArray{Int64,2,Array{Int64,2}}, ::Int64, ::Int64) at /home/tim/.julia/v0.6/OffsetArrays/src/OffsetArrays.jl:89
```

In this case we created `B3`

by explicitly "wrapping" the extracted array inside a type that allows you to supply custom indices. (You can retrieve just the extracted portion with `parent(B3)`

.) We could do the same thing by adjusting the *indices* instead:

```
julia> using IdentityRanges
julia> ind1, ind2 = IdentityRange(2:4), IdentityRange(1:4)
(IdentityRange(2:4), IdentityRange(1:4))
```

An `IdentityRange`

is a range with indices that match its values, `r[i] == i`

. (`ind1, ind2 = OffsetArray(2:4, 2:4), OffsetArray(1:4, 1:4)`

would be functionally equivalent.) Let's use `ind1`

and `ind2`

to snip out the region of the array:

```
julia> B4 = A[ind1, ind2]
OffsetArrays.OffsetArray{Int64,2,Array{Int64,2}} with indices 2:4×1:4:
2 7 12 17
3 8 13 18
4 9 14 19
julia> B4[3,4]
18
```

This implements a simple rule of composition:

`C = A[ind1, ind2]`

, then `C[i, j] == A[ind1[i], ind2[j]]`

Consequently, if your indices have their own unconventional indices, they will be propagated forward to the next stage.

This technique can also be used to create a "view":

```
julia> V = view(A, ind1, ind2)
SubArray{Int64,2,Array{Int64,2},Tuple{IdentityRanges.IdentityRange{Int64},IdentityRanges.IdentityRange{Int64}},false} with indices 2:4×1:4:
2 7 12 17
3 8 13 18
4 9 14 19
julia> V[3,4]
18
julia> V[1,1]
ERROR: BoundsError: attempt to access SubArray{Int64,2,Array{Int64,2},Tuple{IdentityRanges.IdentityRange{Int64},IdentityRanges.IdentityRange{Int64}},false} with indices 2:4×1:4 at index [1, 1]
Stacktrace:
[1] throw_boundserror(::SubArray{Int64,2,Array{Int64,2},Tuple{IdentityRanges.IdentityRange{Int64},IdentityRanges.IdentityRange{Int64}},false}, ::Tuple{Int64,Int64}) at ./abstractarray.jl:426
[2] checkbounds at ./abstractarray.jl:355 [inlined]
[3] getindex(::SubArray{Int64,2,Array{Int64,2},Tuple{IdentityRanges.IdentityRange{Int64},IdentityRanges.IdentityRange{Int64}},false}, ::Int64, ::Int64) at ./subarray.jl:184
```

Note that this object is a *conventional* `SubArray`

(it's not an `OffsetArray`

), but because it was passed `IdentityRange`

indices it preserves the indices of the indices.

As illustrated above for the image rotation example, a recent release (v0.6.0) of the Images package put both the power and responsibility for dealing with arrays with custom indices into the hands of users. One of the key functions in this package is `imfilter`

which can be used to smooth or otherwise "filter" arrays. The idea is that starting from an array `A`

, each local neighborhood is weighted by a "kernel" `kern`

, producing an output value according to the following formula:

This is the formula for correlation; the formula for another operation, convolution, is very similar.

Let's start with a trivial example: let's filter with a "delta function" kernel, meaning it has value `1`

at location 0 and is zero everywhere else. According to the correlation formula, because `kern[J]`

is 1 at `J==0`

, this should simply give us a copy of our original array:

```
julia> using Images
julia> imfilter(1:8, [1])
WARNING: assuming that the origin is at the center of the kernel;
to avoid this warning, call `centered(kernel)` or use an OffsetArray
Stacktrace:
[1] depwarn(::String, ::Symbol) at ./deprecated.jl:64
[2] _kernelshift at /home/tim/.julia/v0.6/ImageFiltering/src/imfilter.jl:1049 [inlined]
[3] kernelshift at /home/tim/.julia/v0.6/ImageFiltering/src/imfilter.jl:1046 [inlined]
[4] factorkernel(::Array{Int64,1}) at /home/tim/.julia/v0.6/ImageFiltering/src/imfilter.jl:1016
[5] imfilter at /home/tim/.julia/v0.6/ImageFiltering/src/imfilter.jl:10 [inlined]
[6] imfilter(::UnitRange{Int64}, ::Array{Int64,1}) at /home/tim/.julia/v0.6/ImageFiltering/src/imfilter.jl:5
[7] eval(::Module, ::Any) at ./boot.jl:235
[8] eval_user_input(::Any, ::Base.REPL.REPLBackend) at ./REPL.jl:66
[9] macro expansion at ./REPL.jl:97 [inlined]
[10] (::Base.REPL.##1#2{Base.REPL.REPLBackend})() at ./event.jl:73
while loading no file, in expression starting on line 0
8-element Array{Int64,1}:
1
2
3
4
5
6
7
8
```

The warning is telling you that Images decided to make a guess about your intention, that the kernel `[1]`

was intended to be centered around zero. You can suppress the warning by explicitly passing the following kernel instead:

```
julia> kern = centered([1])
OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}} with indices 0:0:
1
julia> kern[0]
1
julia> kern[1]
ERROR: BoundsError: attempt to access OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}} with indices 0:0 at index [1]
Stacktrace:
[1] throw_boundserror(::OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}}, ::Tuple{Int64}) at ./abstractarray.jl:426
[2] checkbounds at ./abstractarray.jl:355 [inlined]
[3] getindex(::OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}}, ::Int64) at /home/tim/.julia/v0.6/OffsetArrays/src/OffsetArrays.jl:94
```

By using an `OffsetArray`

you have clearly specified your intended indices for `kern`

.

This can be used to shift an image in the following way (by default, `imfilter`

returns its results over the same domain as the input):

```
julia> kern2 = OffsetArray([1], 2:2) # a delta function centered at 2
OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}} with indices 2:2:
1
julia> imfilter(1:8, kern2, Fill(0)) # pad the edges of the input with 0
8-element Array{Int64,1}:
3
4
5
6
7
8
0
0
julia> kern3 = OffsetArray([1], -5:-5) # a delta function centered at -5
OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}} with indices -5:-5:
1
julia> imfilter(1:8, kern3, Fill(0))
8-element Array{Int64,1}:
0
0
0
0
0
1
2
3
```

These are all illustrated in the following figure:

In this figure, we plotted the kernel as if it were at the location corresponding to convolution rather than correlation.

In other programming languages, when filtering with a kernel that has an even number of elements, it can be difficult to remember which of the two middle elements corresponds to the origin. In Julia, that's not an issue, because you can make that choice for yourself:

```
julia> kern = OffsetArray([0.5,0.5], 0:1)
OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}} with indices 0:1:
0.5
0.5
julia> imfilter(1:8, kern, Fill(0))
8-element Array{Float64,1}:
1.5
2.5
3.5
4.5
5.5
6.5
7.5
4.0
julia> kern = OffsetArray([0.5,0.5], -1:0)
OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}} with indices -1:0:
0.5
0.5
julia> imfilter(1:8, kern, Fill(0))
8-element Array{Float64,1}:
0.5
1.5
2.5
3.5
4.5
5.5
6.5
7.5
```

Likewise, sometimes we might have an application where we simply can't handle the edges properly, and we wish to discard them. For example, consider the following quadratic function:

```
julia> a = [(i-3)^2 for i = 1:9] # a quadratic function
9-element Array{Int64,1}:
4
1
0
1
4
9
16
25
36
julia> imfilter(a, Kernel.Laplacian((true,)))
9-element Array{Int64,1}:
-3
2
2
2
2
2
2
2
-11
```

Those weird values on the edges (for which there is no padding that will "extrapolate" the quadratic) might cause problems. Consequently, let's only extract the values that are well-defined, meaning that all inputs to the correlation formula have explicitly-assigned values:

```
julia> imfilter(a, Kernel.Laplacian((true,)), Inner())
OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}} with indices 2:8:
2
2
2
2
2
2
2
```

Notice that in this case, it returned an `OffsetArray`

so that the values in the result align properly with the original array.

There are many more things you can do with custom indices. As one last illustration, consider the Discrete Fourier Transform, which is defined on a periodic domain. Typically, it's rather difficult to emulate a periodic domain with arrays, because arrays have finite size. However, it's possible to define indexing objects which possess periodic behavior. Here we use the FFTViews package, demonstrating the technique on a simple sinusoid:

```
julia> using FFTViews
julia> a = [sin(2π*x)+0.1 for x in linspace(0,1,16)];
julia> afft = FFTView(fft(a))
FFTViews.FFTView{Complex{Float64},1,Array{Complex{Float64},1}} with indices FFTViews.URange(0,15):
1.6+0.0im
1.498-7.53098im
-0.288537+0.69659im
-0.236488+0.35393im
-0.222614+0.222614im
-0.216932+0.14495im
-0.214217+0.0887316im
-0.212937+0.0423558im
-0.212557+0.0im
-0.212937-0.0423558im
-0.214217-0.0887316im
-0.216932-0.14495im
-0.222614-0.222614im
-0.236488-0.35393im
-0.288537-0.69659im
1.498+7.53098im
```

Now, as every student of Fourier transforms learns, the 0-frequency bin holds the sum of the values in `a`

:

```
julia> afft[0]
1.6000000000000003 + 0.0im
```

Since the mean of a sinusoid is zero, this is (within roundoff error) 16*0.1 = 1.6.

We can also check the amplitude at the Fourier-peak, and explore the periodicity of the result:

```
julia> afft[1]
1.4980046017247872 - 7.53097769363728im
julia> afft[-1] # negative frequencies are OK
1.4980046017247872 + 7.53097769363728im
julia> afft[64+1] # look Ma, it's periodic!
1.4980046017247872 - 7.53097769363728im
julia> length(indices(afft,1)) # but we still know how big it is
16
```

Given the periodicity of `afft`

, the commonly-used `fftshift`

function (e.g., `fftshift(fft(a))`

) can be replaced by `afft[-8:7]`

. While very simple, these techniques make it surprisingly more pleasant to deal with what can sometimes become complex and error-prone index gymnastics.

This has only scratched the surface of what's possible with custom indices. In the opinion of the author, their main advantage is that they can increase the clarity of code by ensuring that "names" (indices) can be endowed with *absolute meaning*, rather than always being "referenced to whatever the first element of this particular array happens to encode."

There is quite a lot of code that hasn't yet properly accounted for the possibility of custom indices — surely, some of it written by the author of this post! So users should be prepared for the possibility that exploiting custom indices will trigger errors in base Julia or in packages. Rather than giving up, users are encouraged to report such errors as issues, as this is the only way that custom indices will, over the course of time, have solid support.

For some algorithms, there appears to be little reason to ever use arrays with anything but 1-based indices; in such cases, there may be no reason to modify existing code. But if your code has a "spatial" interpretation–where location has meaning–then you just might want to give the new facilities a try.

In transitioning existing code, the author of this post has observed the following tendencies:

algorithms that exploit custom indices are sometimes simpler to understand than their "1-locked" counterparts;

if you're porting old code to support custom indices, there's some bad news: if you had to think carefully about the indexing the first time you wrote it, it usually requires significant investment to re-think the indexing, even if the end result is somewhat simpler.

even when a specific algorithm might gain little advantage from supporting arbitrary indices, writing code that is "indices aware" from the beginning is often no harder than writing algorithms that implicitly assume indexing starts at 1.

Developers are referred to Julia's documentation for further guidance.