First Steps with PDL
"Maybe there are a few civilizations out there that have decided to stay home, piddle around and send out some radio waves once in a while."
- Annette Foglino, Space: Is Anyone Out There? Most astronomers say yes, Life, 1 Jul 1989.
It can be very frustrating to read an introductory book which takes a long time teaching you the very basics of a topic, in a "Janet and John" style. While you wish to learn, you are anxious to see something a bit more exciting and interesting to see what the language can do.
Fortunately our task in this book on PDL is made very much easier by the high-level of the language. We can take a tour through PDL, looking at the advanced features it offers without getting involved in complexity.
The aim of this section is to cover a breadth of PDL features rather than any in depth, to give the reader a flavour of what he or she can do using the language and a useful reference for getting started doing real work. Later sections will focus on looking at the features introduced here, in more depth.
Alright, let's do something
We'll assume PDL is correctly installed and set up on your computer system (see http://pdl.perl.org/ for details of obtaining and installing PDL).
For interactive use PDL comes with a program called perldl. This allows you to type raw PDL (and perl) commands and see the result right away. It also allows command line recall and editing (via the arrow keys) on most systems.
So we begin by running the perldl program from the system command line. On a Mac/UNIX/Linux system we would simply type perldl in a terminal window. On a Windows system we would type perldl in a command prompt window. If PDL is installed correctly this is all that is required to bring up perldl.
myhost% perldl
perlDL shell v1.357
PDL comes with ABSOLUTELY NO WARRANTY. For details, see the file
'COPYING' in the PDL distribution. This is free software and you
are welcome to redistribute it under certain conditions, see
the same file for details.
ReadLines, NiceSlice, MultiLines enabled
Reading PDL/default.perldlrc...
Found docs database /usr/lib/perl5/.../PDL/pdldoc.db
Type 'help' for online help
Type 'demo' for online demos
Loaded PDL v2.006 (supports bad values)
pdl>
We get a whole bunch of informational messages about what it is loading for startup and the help system. Note; the startup is completely configurable, an advanced user can completely customize which PDL modules are loaded. We are left with the pdl> prompt at which we can type commands. This kind of interactive program is called a 'shell'. There is also pdl2 which is a newer version of the PDL shell with additional features. It is still under development but completely usable.
Let's create something, and display it:
pdl> use PDL::Graphics::Simple
pdl> imag (sin(rvals(200,200)+1))
The result should look like the image below - a two dimensional sin function. rvals is a handy PDL function for creating an image whose pixel values are the radial distance from the central pixel of the image. With these arguments it creates a 200 by 200 'radial' image. (Try 'imag(rvals(200,200))' and you will see better what we mean!) sin() is the mathematical sine function, this already exists in perl but in the case of PDL is applied to all 40000 pixels at once, a topic we will come back to. The imag() function displays the image. You will see the syntax of perl/PDL is algebraic - by which we mean it is very similar to C and FORTRAN in how expressions are constructed. (In fact much more like C than FORTRAN). It is interesting to reflect on how much C code would be required to generate the same display, even given the existence of some convenient graphics library.
Figure of a two dimensional C<sin> function.
That's all fine. But what if we wanted to achieve the same results in a standalone perl script? Well it is pretty simple:
use PDL;
use PDL::Graphics::Simple;
imag (sin(rvals(200,200)+1));
That's it. This is a complete perl/PDL program. One could run it by typing perl filename. (In fact there are many ways of running it, most systems allows it to be setup so you can just type filename. See your local Perl documentation - then the perlrun manual page.)
Two comments:
1. The statements are all terminated by the ';' character. Perl is like C in this regard. When entering code at the pdl command line the final ';' may be omitted if you wish, note you can also use it to put multiple statements on one line. In our examples from now on we'll often omit the pdl prompt for clarity.
2. The directive use PDL; tells Perl to load the PDL module, which makes available all the standard PDL extensions. (Advanced users will be interested in knowing there are other ways of starting PDL which allows one to select which bits of it you want).
Whirling through the Whirlpool
Enough about the mechanics of using PDL, let's look at some real data! To work through these examples exactly you can download any needed input files from http://sourceforge.net/projects/pdl/files/PDL/PDL%20Book%20Example%20Data%20Set/ and we'll assume you are running any of these examples in the same directory as you have downloaded the input data files.
We'll be playing with an image of the famous spiral galaxy discovered by Charles Messier, known to astronomers as M51 and commonly as the Whirlpool Galaxy. This is a 'nearby' galaxy, a mere 25 million light years from Earth. The image file is stored in the 'FITS' format, a common astronomical format, which is one of the many formats standard PDL can read. (FITS stores more shades of gray than GIF or JPEG, but PDL can read these formats too).
pdl> $a = rfits("m51_raw.fits"); # m51_raw.fits is in current directory
Reading IMAGE data...
BITPIX = -32 size = 262144 pixels
Reading 1048576 bytes
BSCALE = && BZERO =
This looks pretty simple. As you can probably guess by now rfits is the PDL function to read a FITS file. This is stored in the perl variable $a.
This is an important PDL concept: PDL stores its data arrays in simple perl variables ($a, $x, $y, $MyData, etc.). PDL data arrays are special arrays which use a more efficient, compact storage than standard perl arrays (@a, @x, ...) and are much faster to access for numerical computations. To avoid confusion it is convenient to introduce a special name for them, we call them piddles (short for 'PDL variables') to distinguish them from ordinary Perl 'arrays', which are in fact really lists. We'll say more about this later.
Before we start seriously playing around with M51 it is worth noting that we can also say:
pdl> $a = rfits "m51_raw.fits";
Note we have now left off the brackets on the rfits function. Perl is rather simpler than C and allows one to omit the brackets on a function all together. It assumes all the items in a list are function arguments and can be pretty convenient. If you are calling more than one function it is however better to use some brackets so the meaning is clear. For the rules on this 'list operator' syntax see the Perl syntax documentation. From now on we'll mostly use the list operator syntax for conciseness
Let's look at M51:
pdl> imag $a;
Figure of the raw image C<m51_raw.fits> shown with progressively greater contrast using the C<imag> command.
A couple of bright spots can be seen, but where is the galaxy? It's the faint blob in the middle: by default the display range is autoscaled