Q: array of vectors?

General Usage
1

Hello,

What is the optimum representation to use to store vectors over a discrete field.
For example, the discrete field being the 2D pixels of an image or the 3D voxels of an image volume?

I’ll use the 2D pixel field as an example, as the generalization to 3D should be straightforward.

A 2D array with can be instantiated as

an_array::Array{Float64, 2} = zeros(Float64, (n_rows, n_cols))

wheree n_rows and n_cols are the number of rows and columns of, say, an image, respectively.

A 3D vector can be instantiated as

a_vector::Vector{Float64} = [0.0, 0.0, 0.0]

How can I combine the two data structures such that a vector is associated with each pixel?

A concrete example would be the gradient vector at each pixel of an image.

2 
a = fill(Vector{Float32}, nrows, ncols)

But I think you’d want a static type (StaticVector?).

Edit:

using StaticArrays
a = fill(SVector{3,Float32}, nrows, ncols)
3

Unless something special is happening, a dense multi-axis array (tensor) is the most efficient representation.

The axes can be [n_rows, n_cols, n_vector] which would represent an n_vector-dimensional vector over n_rows-by-n_vector grid of points (which covers your example of a gradient).

Then you probably can abstract that away a bit and provide a pretty API instead of indexing it manually.

Depending on the type of computations you will be performing, the order of the axes could also be crucially important for performance.

1 Like
4

This would not be the “optimum” representation as it would be quite non-local in memory. But it is convenient and pretty if speed is not of upmost importance.

5

Definitely seconding the Matrix{SVector{3, Float64}} suggestion. Having the 3-element value at each point represented as an actual 3-element array creates a clear distinction between the 2D coordinates of the pixel and its stored values: You index a pixel with M[i, j] and you access a particular element of that pixel with M[i, j][2]. Otherwise it’s easy to accidentally forget whether the 2nd element of some pixel is M[i, j, 2] or M[2, i, j] in an unlabeled 3D array.

Furthermore, an Array{SVector{3, Float64}, 2} is stored in exactly the same layout as a contiguous array of Float64s, and you can reinterpret one as the other for free. That means you can always use the 3D array representation if you need to:

julia> M = rand(SVector{3, Float64}, 2, 4)
2×4 Matrix{SVector{3, Float64}}:
 [0.365512, 0.245407, 0.960372]  [0.670765, 0.445472, 0.670487]  [0.386643, 0.217216, 0.364761]   [0.0798776, 0.318391, 0.297401]
 [0.301753, 0.307553, 0.800906]  [0.899388, 0.033749, 0.285046]  [0.362781, 0.590045, 0.0470405]  [0.916995, 0.729333, 0.066537]

julia> reshape(reinterpret(Float64, M), :, 2, 4)
3×2×4 reshape(reinterpret(Float64, ::Matrix{SVector{3, Float64}}), 3, 2, 4) with eltype Float64:
[:, :, 1] =
 0.365512  0.301753
 0.245407  0.307553
 0.960372  0.800906
3 Likes
6

BTW, this is probably the wrong order since Julia is column-major (or first-axis-major in general). You probably want to store the n_vector elements adjacent to one another, which suggests [n_vector, n_rows, n_cols] as the order. This is all the more reason to use an Array{SVector{...}} instead, since it ensures that the individual pixel elements are contiguous without you having to remember which order the axes should go in.

2 Likes
7

It might be worth it to look at Images.jl and the whole ecosystem around it. Images are basically generic because of the many pixel types available. Leveraging the tools from the ecosystem is bound to payoff later.

1 Like
8

Agreed. No point in reinventing an advanced, well-developed, wheel.
However, my understanding is that the
Images.jl
package is for 2D images, whereas my eventual application involves 3D image volumes.

9

My thanks to all who provided comments and solutions.
Much appreciated.

1 Like
10

Actually, my code fragment was not right. To initialize, use

nrows, ncols = 2, 3
a = fill(Vector{Float32}([0,0,0]), nrows, ncols)
11

Not at all! Images.jl actually works very well for images of any dimension (a big advantage over tools in other languages). The core image type in Images.jl is simply a standard Julia Array of pixels. For example, a standard RGB image using float values for the channels is just an Array{RGB{Float32}, 2} i.e. a 2D array whose element type is RGB. As with the SVector example posted above, this is exactly equivalent to a 3D array of Float32, and you can use the channelview() function to map between the two. Likewise, you can create a 3D image volume as an Array{RGB{Float32}, 3}.

Here’s a random 3D image volume:

julia> img = rand(RGB{Float32}, 2, 2, 2)
2×2×2 Array{RGB{Float32},3} with eltype RGB{Float32}:

and here it is reinterpreted as a 4D array of numbers:

julia> channelview(img)
3×2×2×2 reinterpret(reshape, Float32, ::Array{RGB{Float32},3}) with eltype Float32:
3 Likes
12

Thanks for correcting my misunderstanding and for your examples.
I’m too used to image analysis packages being limited to only two dimensions in other languages and jumped to conclusions after a cursory look.

Will investigate Images.jl further. Looks like there are a lot of useful tools.