| Multi-Dimensional Arrays |  |
Section Contents
- Basic MultiArray Usage
- Array Indexing via Coordinates
- One-dimensional Indexing and Scan-Order Iterator
- Important MultiArray Methods
Basic MultiArray Usage
vigra::MultiArray is the most fundamental data structure in VIGRA. It holds a rectangular block of values in arbitrary many dimensions. Most VIGRA functions operate on top of MultiArray or the associated class MultiArrayView (see The MultiArrayView Interface).
A 2D image can be interpreted as a matrix, i.e. a 2-dimensional array, where each element holds the information of a specific pixel. Internally, the data are stored in a single 1-dimensional piece of memory, and MultiArray encapsulates the entire mapping between our familiar 2-dimensional notation and the raw memory. Pixels in an image are identified by a coordinate pair (x,y), where indexing starts at 0. That is, the pixels in a 800x600 image are indexed by x = 0,...,799 and y = 0,...,599. The principle analoguously extends to higher dimensions.
The structure of a multidimensional array is given by its shape vector, and the length of the shape vector is the array's dimension. The dimension must be fixed as a template parameter at compule time, while the shape is passed to the array's constructor. The second important template parameter is the pixel's value_type, as you know it form std::vector.
To represent the data of a gray scale image, we just need to store one value per pixel, so we choose a 2-dimensional array, where each element has the unsigned char type (in VIGRA, this type is also available as vigra::UInt8). We instantiate a gray scale image object like this:
By default, VIGRA arrays are always zero-initialized. Another initial value can be provided in the constructor, or later via the init() function or the assignment operator:
The Shape2 typedef also exists for higher dimensions up to five as Shape3 etc. If you need even more dimensions, use MultiArrayShape<N>::type instead, were N is the number of dimensions:
When storing RGB images we obviously can't simply use the unsigned char type because every pixel contains 3 numbers: values for red, green and blue. Mathematically, you want to store a data vector for each pixel. To this end, VIGRA provides the vigra::RGBValue<ValueType> class. So for RGB-images just use:
vigra::RGBValue<ValueType> is a specialised 3-dimensional vector containing ValueType elements. Arbitrary short vectors can be stored in the TinyVector<ValueType, SIZE> class, which is the base class of RGBValue. It's length must be specified at compile time in the template parameter SIZE. Vectors whose length is known at compile time are very useful for the compiler to produce highly optimized code. Therefore, Shape2 and it's higher-dimensional cousins are instances of TinyVector as well.
Alternatively you can use a 3-dimensional array vigra::MultiArray<3, unsigned char> to represent a color image. The third dimension has size 3 and contains the information for the red, green and blue channel.
However, we are not going to use this form in the tutorial because vigra::RGBValue<ValueType> provides many helpful methods that are not available when color is just another array dimension.
Array Indexing via Coordinates
The easiest way to access the values of an array is via the coordinates. A tuple of coordinates can again be specified by the appropriate shape object, which must be passed to the array's indexing operator:
Important Remark: Notice that VIGRA follows the mathematical convention of the index order: dimension 0 corresponds to the x (horizontal) coordinate, dimension 1 to the y (vertical) coordinate, and so on. Accordingly, dimension 0 is changing fastest in memory: when we increase x by one, we get to the next memory location. In matrix jargon, this is also known as Fortran order. Many image processing libraries (e.g. Image Magick, OpenCV, and Qt) use the same convention. However, others like Matlab and numpy, use the reverse order (so called C order). Don't be confused!
Internally, shape objects are implemented in terms of the powerful vigra::TinyVector class. This means that shape objects support the usual mathematical operations like addition, multiplication and scalar products. Coordinate computations can thus be performed on entire coordinate objects at once - there is no need to manipulate the individual coordinates separately.
Nonetheless, in some circumstances it is more convenient to provide coordinates individually rather than in a shape object. This is possible with round brackets (x,y):
This kind of indexing is supported for dimensions up to five, and only if the array's dimension is known (this is not always the case: in a generic function where the dimension is a template parameter, you must use shape objects). In combination with the method shape(n), that returns the length of the n-th dimension, we can use the coordinates to set the element of an entire row or column:
On first sight, individual coordinates may seem to be necessary for iterating over the image or parts of it. But the following example shows that the same effect can be achieved with a shape object that is allocated outside the loop: 3rd column of a 8x4-matrix (initialized with 5) to 7.
We will discuss more powerful methods to access certain parts of an array in section Important MultiArray Methods.
One-dimensional Indexing and Scan-Order Iterator
Regardless of the array's dimension, it is always possible to access elements with 1-dimensional index, its scan-order index, via the normal indexing operator. For example, array[1] refers to the index of the second array element. Defining a scan order is often called flattening of an array, because a high-dimensional data structure is accessed like a 1-dimensional vector. Notive that scan-order access in VIGRA does not require the data to be copied.
VIGRA defines scan-order by considering the dimensions from front to back. Thus, items are accessed such that only the x coordinate is incremented, while y (and possibly further coordinates) are held fixed at 0. When x is exhausted, y is incremented by one and the iteration starts again at x=0. To control iteration, the function array.size() returns the total number of elements:
Alternatively, scan-order access can be achieved with an STL-compatible iterator pair obtained by calling array.begin() and array.end(). Continuing with the example above, we can write:
The iterator is implemented by class StridedScanOrderIterator which encapsulates all the bookkeeping necessary to get the elements in the correct order, even when the array was transposed (see below).
Scan-order access is useful to implement pointwise operations, e.g. the addition of two matrices. The following code adds two matrices and stores the result in the first one:
This is convenient because the actual high-dimensional shapes of the arrays are of no significance for point-wise operations as long as the shapes match. Be careful: the arrays themselves have no way of checking this condition. So thefollowing code using two transposed shapes is perfectly valid C++, but has probably not the intended effect:
By the way: VIGRA provides the += operator (and is cousins) to write this more concisely, and this operator throws an exception if the shapes don't match:
For more information on mathematical operations on arrays see the multi_math module.
As mentioned, VIGRA's scan order is similar to the NumPy-method array.flatten(). You use it, to copy a multi-dimensional array into an one-dimensional array, or to access elements in flattened order. The only difference is that NumPy uses "C-order" , i.e. the rightmost dimension takes priority, whereas VIGRA uses Fortran-order, i.e. the leftmost dimension takes priority. A method like flatten can be implemented in VIGRA as:
To show the difference between VIGRA and NumPy we'll add the NumPy output, i.e. the result when we had used C-order in the code above:
flatArray - index 0 1 2 3 4 5 ----------------------------------------------------------------- VIGRA-output: 1 2 3 4 5 6 intArray - index [0,0] [1,0] [2,0] [0,1] [1,1] [2,1] ----------------------------------------------------------------- NumPy-output: 1 4 2 5 3 6 intArray - index [0,0] [0,1] [1,0] [1,1] [2,0] [2,1]
To change the axis priorities of the StridedScanOrderIterator, look at the transpose-function in the next section.
Important MultiArray Methods
This part of the tutorial explains important methods of MultiArray. However, before we proceed, we need to introduce the class vigra::MultiArrayView.
The MultiArrayView Interface
A vigra::MultiArrayView has the same interface as a MultiArray (with the exception of reshape() and swap()), but it doesn't own its data. Instead, it provides a view onto the data of some other array. In contrast, a MultiArray owns its data and is responsible for freeing it in the destructor. MultiArrays are automatically converted into MultiArrayViews when needed.
The point of this distinction is that MultiArrayViews can be used to access and manipulate the same data in many different ways without any need for creating copies. For example, we can work with a 2-dimensional slice of a volume dataset (i.e. a lower dimensional part of a 3D array) without first copying the slice into a 2D image. This is possible whenever the desired view can be realized by just manipulating the internal mapping from indices and shapes to memory locations, and not the memory layout itself.
This possibility – which is similarly implemented in other packages like Matlab and numpy – is a key ingredient for efficient computations with multi-dimensional arrays. Thus, most VIGRA functions actually receive MultiArrayViews to maximize flexibility. This section describes the most important ways to create new MultiArrayViews from an existing array or view. The complete documentation is available in the vigra::MultiArrayView reference.
subarray(p,q)
This method creates a rectangular subarray of your array between the points p and q, where p (the starting point of the subregion) is included, q (the ending point) is not. subarray does not change the dimension of the array (this is the task of the various bind-methods).
To give an example, we create a 4x4 array that consitst of a checkerboard with 2x2 squares:
The positions p and q are specified with the familiar Shape objects. In this example we simply overwrite parts of the array. The following larger example uses subarray to output a half-sized subimage around the center of the original image: subimage_tutorial.cxx
The relevant part of this code is shown here (the functions importImage and exportImage are described in section Image Input and Output):
After reading the (here: gray scale) image data to an array we need to calculate the coordinates of our subimage. In this case we want to cut out the middle part of the image. Afterwards we write the subimage into a new array. Look at the result:
input file | 
subimage output file |
bind<M>(i) and bindAt(M, i)
These methods bind axis M to the index i and thus reduce the dimension of the array by one. The only difference between the two forms is that the axis to be bound must be known at compile time in the first form, whereas it can be specified at runtime in the second.
Binding is useful when we want to access and/or manipulate a particular row or column of an image, or a single slice of a volume. In principle, the same can also be achieved by explicit loops, but use of bind often leads to more elegant and more generic code. Consider the following code to initialize the third column of an image with the constant 5:
NumPy-Users are familiar with the bind mechanism as "slicing". The example above written in numpy syntax becomes:
array2d[2, :] = 5 // NumPy-equivalent of array2d.bind<0>(2) = 5
You can also initialize a lower-dimensional array with the bind-method:
The array array1d contains the elements the 3rd column of array2d. This bahavior nicely illustrates the difference between a copy and a view: array1d contains a copy of the 3rd column, whereas the bind function only creates a new view to the existing data in array2d.
At this point we have to distinguish between the classes MultiArray and MultiArrayView . MultiArray inherits from MultiArrayView and contains the memory management of the array. With MultiArrayView we can view the data stored in a MultiArray. The code above produces a copy of the 3rd column of intArray. If we change the elements of lowArray nothing happens to intArray .
Moving on to image processing we'll give an example how you can flip an image by using bind. We read a gray scale image into a 2-dimensional array called imageArray . Then we initalize a new array newImageArray of the same dimension and size and set the first row of newImageArray to the values of the last row of imageArray , the second row to the values of the second last row and so on. Hence, we flip the image top to bottom.
However, you don't need to implement a method like this yourself because VIGRA already provides the function reflectImage(). We use this function to flip the image left-to-right:
The complete example can be found in mirror_tutorial.cxx. (This program needs an infile and an outfile as command-line arguments and contains additional I/O code which will be explained in section Image Input and Output.) Here you can see what happens to an input file:
|
input file | 
mirrored top to bottom | 
mirrored left to right |
For completeness, there are five additional versions of the bind()-method:
- bindInner(i) with scalar or multi-dimensional index i:
- if i is an
integer, the innermost dimension (axis 0) is fixed to i,
if i isMultiArrayShape<M>::type(a shape of size M), then the M innermost dimensions (axes 0...M-1) are fixed to the values in the shape vector - bindOuter(i) with scalar or multi-dimensional index i:
- if i is an
integer, the outmost dimension (axis N-1) is fixed to i,
if i isMultiArrayShape<M>::type(a shape of size M), then the M outmost dimensions (axes N-M ... N-1) are fixed to the values in the shape vector - diagonal() :
- Create a 1-dimensional view to the diagonal elements of the original array (i.e.
view[i] == array(i,i,i)for a 3D original array).
The opposite of binding - inserting a new axis - is also possible. However, since we cannot alter the internal memory layout and thus cannot insert additional data elements, a new axis must be singleton axis, i.e. an axis with length 1. The argument of insertSingletonDimension(k) determines the position of the new axis, with 0 <= k <= N when the original array has N dimensions:
expandElements(k) and bindElementChannel(i)
When the array elements are vectors (i.e. vigra::TinyVector or vigra::RGBValue), we can expand these elements into an addtional array dimension:
The argument k of expandElements(k) determines the desired position of the channel axis, i.e. the index that refers to the vector elements. When the original vector array has N dimensions (not counting the channel axis), it is required that 0 <= k <= N.
Often, we are only interested in a single channel of a vector-valued array. This can be achieved with the function bindElementChannel(i). For example, we can extract the green channel (i.e. channel 1) from an RGB image like this:
This is simply an abbreviation for rgb_array.expandElements(0).bindInner(1).
transpose()
Everyone is familiar with the transpose() function of a matrix (i.e. a 2-dimensional array). Once again, this operation is possible without copying the data by just manipulating the internal access functions. The following example demonstrates the difference between a transposed copy and view:
The output is:
base_array: 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 transarray: 0 0 0 0 1 1 1 1 2 2 2 2 3 3 3 3 transArrayView: 0 0 0 0 1 1 1 1 2 2 2 2 3 3 3 3 base_array after setting transarray to 5 (no change, since transarray is a copy): 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 base_array after setting transarrayView to 5 (base_array changes because transarrayView is a view): 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
The function MultiArrayView::transpose() generalizes transposition to arrays of arbitrary dimensions. Here, it just reverses the order of the axes: axis 0 becomes axis N-1, axis 1 becomes axis N-2 and so on. In the following example we transpose a 5D array and print out its shape.
The output is:
Shape of array5D: (1, 2, 3, 4, 5) Shape of array5D view after default transpose(): (5, 4, 3, 2, 1)
Finally, MultiArrayView::transpose() can also be called with a shape object that specifies the desired permutation of the axes: When permutation[k] = j, axis j of the original array becomes axis k of the transposed array (remember, that VIGRA counts the axes from 0):
The permutation in the example is 2,1,3,4,0. Thus, original dimension 0 appears in the last position of the new view, original dimension 2 appears in the first position, and so on as demonstrated by the output of the example:
Shape of array5D view after user-defined transpose(): (3, 2, 4, 5, 1)
(applied permutation 2 => 0, 1 => 1, 3 => 2, 4 => 3, 0 => 4 to the axes)
When we transpose an image about the major diagonal, we can simply use the view created by MultiArrayView::transpose(). However, transposition about the minor diagonal requires a new image, which can be filled by transposeImage() like this:
The result is:
|
input file | 
transpose about major diagonal | 
transpose about minor diagonal |
In VIGRA, image transposition is also implemented in function vigra::transposeImage(...). The difference is that transposeImage() copies the image data, whereas MultiArrayView::transpose() just changes the internal mapping from indices to memory locations.
Important note: Transposing an array also changes the direction of the StridedScanOrderIterator. Imagine a 3x4- matrix. Scan-order means that we iterate from left to right, row by row. Now, let's transpose the matrix to a 4x3 view. Than, scan-order in the new view is again left to right, row by row. However, in the original matrix this now corresponds to a transposed scan: from top to bottom, column by column. The same applies to the array's index operator with integer argument.
isUnstrided(k)
A MultiArray always accesses its elements in consecutive memory order, i.e. &array[i] == &array.data()[i] for all i in the range [0, array.size()). However, this does in general not hold for MultiArrayViews, because changing array access is the whole point of view creation. Sometimes, it is necessary to find out if a view still has consecutive, unstrided memory access, for example when you want to pass on the view's data to an external library that only accepts plain C arrays: When the view happens to be unstrided, you can avoid to create a copy of the data. You can determine this with the function isUnstrided(k) which returns true when the array is unstrided up to dimension k (k defaults to N-1, i.e. the entire array must be unstrided):