The PDL Book by PDL 2.006 Various Contributors - HTML preview

PLEASE NOTE: This is an HTML preview only and some elements such as links or page numbers may be incorrect.
Download the book in PDF, ePub, Kindle for a complete version.

The other important information supplied by the signature is the dimensionality of each of these arguments in an elementary operation. Each formal parameter carries this information in a list of formal dimension identifiers enclosed in parentheses. So indeed a(n) marks a as a one-dimensional argument. Additionally, each dimension has a named size in a signature, in this example n . c() has an empty list of dimension sizes: it is declared to be zero-dimensional (a scalar).

If piddles that are supplied as runtime arguments to a function have more dimensions than specified for their respective formal arguments in the signature then these dimensions are treated by PDL as extra dimensions and lead to the operation being threaded over the appropriate subslices, just what we have seen in the simple examples above.

As mentioned before a higher dimensional piddle can be viewed as an array (again not in the perl array sense) of lower dimensional subslices. Anybody who has ever worked with matrix algebra will be familiar with the concept. For some of the following examples it will be useful to illustrate this concept in somewhat more detail. Let's make a piddle first, a simple 3D piddle:

$pdl = sequence(3,4,5);

A boring piddle, you say? Yes, boring, but simple enough to clearly see what is going on in the following. First we look at it as a 3D array of 0D subslices. Since we know the syntax of the slice method already we can write down all 0D subslices, no problem:

$pdl(($i),($j),($k));

Well, obviously we have not written down all 3*4*5 = 60 subslices literally but rather in a more concise way. It is understood that $i can have any value between 0 and 2, $j between 0 and 3 and $k between 0 and 4. To emphasize this we sometimes write

$pdl(($i),($j),($k)) $i=0..2; $j=0..3; $k=0..4

With the meaning as above (and '..' not meaning the perl list operator). In that way we enumerate all the subslices. Quite analogously, when dealing with an elementary operation that consumes 1D slices we want to view $pdl as an [4,5] array of 1D subslices:

$pdl(:,($i),($j)) $i=0..3; $j=0..4

And similarly, as a [5] array of 2D subslices:

$pdl(:,:,($i)) $i=0..4

You see how we just insert a ":" for each complete dimension we include in the subslice. In fig. XXX the situation is illustrated graphically for a 2D piddle. Depending on the dimensions involved in an elementary operation we therefore often group the dimensions (what we call the shape) of a piddle in a form that suggests the interpretation as an array of subslices. For example, given our 3D piddle above that has a shape [3,4,5] we have the following different interpretations:

    ()[3,4,5]  a shape [3,4,5] array of 0D slices

    (3)[4,5]   a shape [4,5]   array of 1D slices (of shape [3])

    (3,4)[5]   a shape [5]     array of 2D slices (of shape [3,4])

    (3,4,5)[]  a        0D     array of 3D slices (of shape [3,4,5])

The dimensions in parentheses suggest that these are used in the elementary operation (mimicking the signature syntax); in the context of threading we call these the elementary dimensions . The following group of dimensions in rectangular brackets are the extra dimensions. Conversely, given the elementary/extra dims notation we can easily obtain the shape of the underlying piddle by appending the extradims to the elementary dims. For example, a [3,6,1] array of 2D subslices of shape [3,4]:

(3,4)[3,6,1]

identifies our piddle's shape as [3,4,3,6,1]

Alright, the principles are simple. But nothing is better than a few examples. Again a typical imaging processing task is our starting point. We want to convert a colour image to grayscale. The input image is represented as a two-dimensional array of triples of RGB colour coordinates, or in other words, a piddle of shape [3,n,m] . Without delving too deeply into the details of digital colour representation it suffices to note that commonly a gray value i corresponding to a colour represented by a triple of red, green and blue intensities (r,g,b) is obtained as a weighted sum:

       77        150       29

   i = --- r  +  --- g  +  --- b

       256       256       256

A straight forward way to compute this weighted sum in PDL uses the inner function. This function implements the well-known inner product between two vectors. In a elementary operation inner computes the sum of the element-by-element product of two one-dimensional subslices (vectors) of equal length:

        __ n-1

   c = \          a   b

       /__ i=0     i   i

Now you should already be able to guess inner's signature:

pdl> sig inner

Signature: inner(a(n); b(n); [o]c())

a(n); b(n); [o]c(); : two one-dimensional input parameters a(n) and b(n) and a scalar output parameter c() . Since a and b both have the same named dimension size n the corresponding dimension sizes of the actual arguments will have to match at runtime (which will be checked by PDL!). We demonstrate the computation starting with a colour triple that produces a sort of yellow/orange on an RGB display:

    $yel = byte [255, 214, 0]; # a yellowish pixel

    $conv = float([77,150,29])/256; # conversion factor

    $i = inner($yel,$conv)->byte; # compute and convert to byte

    print "$i \backslash n";

    202

Now threading makes extending this example to a whole RGB image very straightforward:

    use PDL::IO::Pic; # IO for popular image formats

    $im = rpic 'pdllogo.jpg'; # a colour image from the book dataset

    $gray = inner($im->pdim('COLOR'),$conv);

       # threaded inner product over all pixels

    $gb = $gray->byte; # back to byte type

    COLOR Byte D [3,500,300]

The code needs no modification! Let us analyze what is going on. We know that $conv has just the required number of dimensions (namely one of size 3). So this argument doesn't require PDL to perform threading. However, the first argument $im has two extra dimensions (shape [500,300]). In this case threading works (as you would probably expect) by iterating the inner product over the combination of all 1D subslices of $im with the one and only subslice of $conv creating a resulting piddle (the grayscale image) that is made up of all results of these elementary operations: a 500x300 array of scalars, or in other words, a 2D piddle of shape [500,300].

We can more concisely express what we have said in words above in our new way to split piddle arguments in elementary dimensions and extra dimensions. At the top we write inner's signature and at the bottom the slice expressions that show the subslices involved in each elementary operation:

    Piddles $im $conv $gray

    Signature a(n); b(n); [o]c()

    Dims (3)[500,300] (3)[] ()[500,300]

    Slices ":,($i),($j)" ":" "($i)($j)"

Remember that the slice notation at the bottom does not mean that you have to generate all these slices yourself. It rather tells you which subslices are used in a elementary operation. It is a way to keep track what is going on behind the scenes when PDL threading is at work. Threading makes it possible that we can call the grayscale conversion with piddles representing just one RGB pixel (shape [3]), a line of RGB pixels (shape [3,n]), RGB images (shape [3,m,n]), volumes of RGB data (shape [3,m,n,o]), etc. All we have to do is wrap the code above into a small subroutine that also does some type conversion to complete it:

    sub rgbtogr {

       my ($im) = @_;

       my $conv = float([77,150,29])/256; # conversion factor

       my $gray = inner $im, $conv;

       return $gray->convert($im->type); # convert to original type

    }

You can write your own threading routines

Did you notice? By writing this little routine we have created a new function with its own signature that will thread as appropriate. It has inherited the ability to thread from inner. So what is the signature of rgbtogr? It is nowhere written explicitly and we can't use the sig function to find out about it sig will only know about functions that were created using PDL::PP or if we explicitly specified the signature in the PDL documentation but from the properties of inner and the definition of rgbtogr we can work it out. As input it takes piddles with a size of the first dimension of 3 and returns for each of the 1D subslices a 0D result (the gray value). In other words, the signature is

a(tri = 3); [o] b()

There is some new syntax in this signature that we haven't seen before: writing tri=3 signifies that in a elementary operation rgbtogr will work on 1D subslices (we have encountered this before); additionally, the size of the first dimension (named