just sums up all the data values in all the pixels in the
image - in this case the answer is 17916010. If the image is linear
(which it is) and if it was calibrated (i.e. we knew the relation
between data numbers and brightness units) we could work out the
total brightness. Let's turn it round - we know that M51 has
a luminosity of about 1E36 Watts, so we can work out what
one data value corresponds to in physical units:
pdl> p 10**36/sum($gal)
5.58159992096455e+28
This is also about 200 solar luminosities, (Note we have switched to using C
as a shorthand for C - which only works in the C and C shells)
which gives 4 billion solar luminosities for the whole galaxy.
OK we do not need PDL for this simple arithmetic, let's get back to
computations that involve the whole image.
How can we get the sum of a piece of an image, e.g. near the centre?
Well in PDL there is more than one way to do it (Perl aficionados call
this phenomenon TIMTOWTDI). In this case, because we really want
the brightness in a circular aperture, we'll use the C
function:
pdl> $r = rvals $gal;
pdl> imag $r;
Remember C? It replaces all the pixels in an image with its distance
from the centre. We can turn this into a I with a simple
operation like:
pdl> $mask = $r<50;
pdl> imag $mask;
=for html
Figure: Using C to generate a mask image to isolate the galaxy bulge
and disk. Top row: radial gradient image C<$r>, and radial gradient
masked with less than operator C<$r < 50>. Bottom row: Bulge and disk of
the galaxy.
The Perl I is applied to all pixels in the image.
You can see the result is an image which is 0 on the outskirts and 1 in
the area of the nucleus. We can then simply use the mask image to
isolate in a simple way the bulge and disk components (lower row) and it
is then very easy to find the brightness of both pieces of the M51
galaxy:
pdl> $bulge = $mask * $gal
pdl> imag $bulge,0,300
Displaying 512 x 512 image from 0 to 300 ...
pdl> print sum $bulge;
3011125
pdl> $disk = $gal * (1-$mask)
pdl> imag $disk,0,300
Displaying 512 x 512 image from 0 to 300 ...
pdl> print sum $disk
14904884
You can see that the disk is about 5 times brighter than the bulge in
total, despite its more diffuse appearance. This is typical for
spiral galaxies. We might ask a different question: how does the average
I, the brightness per unit area on the sky,
compare between bulge and disk? This is again quite straight forward:
pdl> print sum($bulge)/sum($mask);
pdl> print sum($disk)/sum(1-$mask);
We work out the area by simply summing up the 0,1 pixels in the mask
image. The answer is the bulge has about 7 times the surface
brightness than the disk - something we might have guessed from
looking at the above figure, which tells astronomers its stellar density is
much higher.
Of course PDL being so powerful, we could have figured this out in one line:
pdl> print avg($gal->where(rvals($gal)<50)) / avg($gal->where(rvals($gal)>=50))
6.56590509414673
=head2 Twinkle, twinkle, little star
Let's look at something else, we'll zoom in on a small piece of the image:
pdl> $section = $gal(337:357,178:198);
pdl> imag $section; # the bright star
Here we are introducing something new - we can see that PDL supports
I to the Perl syntax. We can say C<$var(a:b,c:d...)> to specify
I. In this case we have produced a sub-image ranging
from pixel 337 to 357 along the first dimension, and 178 through 198 along the
second. Remember pdl data dimension indexes start from zero. We'll talk some
more about I later on. This sub-image happens to contain
a bright star.
At this point you will probably be able to work out for yourself the amount of
light coming from this star, compared to the whole galaxy. (Answer: about 2%)
But let's look at something more involved: the radial profile of the star.
Since stars are a long way away they are almost point sources, but our camera
will blur them out into little disks, and for our analysis we might want an
exact figure for this blurring.
We want to plot all the brightness of all the pixels in this section, against
the distance from the centre. (We've chosen the section to be conveniently
centred on the star, you could think if you want about how you might determine
the centroid automatically using the C and C functions). Well it
is simple enough to get the distance from the centre:
pdl> $r = rvals $section;
But to produce a one-dimensional plot of one against the other we need to
reduce the 2D data arrays to one dimension. (i.e our 21 by 21 image section
becomes a 441 element vector). This can be done using the PDL C
function, which 'clumps' together an arbitrary number of dimensions:
pdl> $rr = $r->clump(2); # Clump first two dimensions
pdl> $sec = $section->clump(2);
pdl> points $rr, $sec; # Radial plot
You should see a nice graph with points like those
in the figure below
showing the drop-off
from the bright centre of the star. The blurring is usually measured
by the 'Full Width Half Maximum' (FWHM) - or in plain terms how
fat the profile is across when it drops by half. Looking at Figure 1.6
it looks like this is about 2-3 pixels - pretty compact!
=for html
Figure: Radial light profile of the bright star with fitted curve.
Well we don't just want a guess - let's fit the profile with a function.
These blurring functions are usually represented by the C
function. PDL comes with a whole variety of general purpose and
special purpose fitting functions which people have written for
their own purposes (and so will you we hope!). Fitting Gaussians
is something that happens rather a lot and there is surpisingly
enough a special function for this very purpose. (One could use
more general fitting packages like C or
C but that would require more care).
pdl> use PDL::Fit::Gaussian;
This loads in the module to do this. PDL, like Perl, is modular. We
don't load all the available modules by default just a convenient
subset. How can we find useful PDL functions and modules? Well
C tells us more about what we already know, to find out
about what we don't know use C:
pdl> apropos gaussian
PDL::Fit::Gaussian ...
Module: routines for fitting gaussians
PDL::Gaussian Module: Gaussian distributions.
fitgauss1d Fit 1D Gassian to data piddle
fitgauss1dr Fit 1D Gassian to radial data piddle
gefa Factor a matrix using Gaussian elimination.
grandom Constructor which returns piddle of Gaussian random numbers
ndtri The value for which the area under the Gaussian probability density function (integrated from minus
infinity) is equal to the argument (cf erfi). Works inplace.
This tells us a whole lot about various functions and modules to do with
Gaussians. Note that we can abbreviate C and C
with 'C' and 'C' when using the C or C shells.
Let's fit a Gaussian:
pdl> use PDL::Fit::Gaussian;
pdl> ($peak, $fwhm, $background) = fitgauss1dr($rr, $sec);
pdl> p $peak, $fwhm, $background;
C is a function in the module L which fits
a Gaussian constrained to be radial (i.e. whose peak is at the origin).
You can see that, unlike C and FORTRAN, Perl functions can return
more than one result value. This is pretty convenient. You can see the
FWHM is more like 2.75 pixels. Let's generate a fitted curve with this
functional form.
pdl> $rrr = sequence(2000)/100; # Generate radial values 0,0.01,0,02..20
# Generate gaussian with given FWHM
pdl> $fit = $peak * exp(-2.772 * ($rrr/$fwhm)**2) + $background;
Note the use of a new function, C, which
generates a new piddle with N values ranging 0..(N-1).
We are simply using this to generate the horizontal axis values
for the plot. Now let's overlay it on the previous plot.
pdl> hold; # This command stops new plots starting new pages
pdl> line $rrr, $fit, {Colour=>2} ; # Line plot
The last C command shows the PDL syntax for optional function
arguments. This is based on the Perl's built in hash syntax. We'll say
more about this later in L. The result should look a
lot like the figure above. Not too bad. We could perhaps do a bit
better by exactly centroiding the image but it will do for now.
Let's make a I of the 2D stellar image. This is equally
easy:
pdl> $fit2d = $peak * exp(-2.772 * ($r/$fwhm)**2);
pdl> release; # Back to new page for new plots;
pdl> imag $fit2d;
pdl> wfits $fit2d, 'fake_star.fits'; # Save our work
But the figure below is a
boring. So far we have been using simple 2D graphics from the
C library. In fact PDL has more
than one graphics library (some see this as a flaw, some
as a feature!). Using the C library
which does OpenGL graphics we can look at our simulated
star in 3D (see the right hand panel);
=for html
Figure: Two different views of the 2D simulated Point Spread Function.
pdl> use PDL::Graphics::TriD; # Load the 3D graphics module
pdl> imag3d [$fit2d];
If you do this on your computer you should be able to look at the graphic from
different sides by simply dragging in the plot window with the mouse! You can
also zoom in and out with the right mouse button. Note that C has it's
a rather different syntax for processing it's arguments - for very good reasons
- we'll explore 3D graphics further in L.
Finally here's something interesting. Let's take our fake star and place it
elsewhere on the galaxy image.
pdl> $newsection = $gal(50:70,70:90);
pdl> $newsection += $fit2d;
pdl> imag $gal,0,300;
We have a bright new star where none existed before! The C-style C<+=>
increment operator is worth noting - it actually modifies the contents of
C<$newsection> in-place. And because C<$newsection> is a I of C<$gal>
the change also affects C<$gal>. This is an important property of slices - any
change to the slice affects the I. This kind of parent/child
relationship is a powerful property of many PDL functions, not just slicing.
What's more in many cases it leads to memory efficiency, when this kind of
linear slice is stored we only store the start/stop/step and not a new copy of
the actual data.
Of course sometimes we DO want a new copy of the actual data, for example if we
plan to do something evil to it. To do this we could use the alternative form:
pdl> $newsection = $newsection + $fit2d
Now a new version of C<$newsection> is created which has nothing to
do with the original C<$gal>. In fact there is more than one way to do
this as we will see in later chapters.
Just to amuse ourselves, lets write a short script to cover M51 with dozens of fake
stars of random brightnesses:
use PDL;
use PDL::Graphics::PGPLOT;
use PDL::NiceSlice; # must use in each program file
srand(42); # Set the random number seed
$gal = rfits "fixed_gal.fits";
$star = rfits "fake_star.fits";
sub addstar {
($x,$y) = @_;
$xx = $x+20; $yy = $y+20;
# Note use of slice on the LHS!
$gal($x:$xx,$y:$yy) += $star * rand(2);
}
for (1..100) {
$x1 = int(rand(470)+10);
$y1 = int(rand(470)+10);
addstar($x1,$y1);
}
imag $gal,0,1000;
This ought to give the casual reader some flavour of the Perl syntax - quite simple
and quite like C except that the entities being manipulated here are entire
arrays of data, not single numbers. The result is shown, for amusement,
in the figure below and takes virtually no time to compute.
=for html
Figure: M51 covered in fake stars.
=head2 Getting Complex with M51
To conclude this frantic whirl through the possibilities of PDL, let's look at
a moderately complex (sic) example. We'll take M51 and try to enhance it to reveal the
large-scale structure, and then subtract this to reveal small-scale structure.
Just to show off we'll use a method based on the Fourier transform - don't
worry if you don't know much about these, all you need to know is that the
Fourier transform turns the image into an 'inverse' image, with
complex numbers, where each pixel
represents the strength of wavelengths of different scales in the image.
Let's do it:
pdl> use PDL::FFT; # Load Fast Fourier Transform package
pdl> $gal = rfits "fixed_gal.fits";
Now C<$gal> contains real values, to do the Fourier transform it has to
have complex values. We create a variable C<$imag> to hold the imaginary
component and set to zero.(For reasons of efficiency complex numbers
are represented in PDL by seperate real and imaginary arrays - more about this
in Chapter 2.)
pdl> $imag = $gal * 0; # Create imaginary component, equal to zero
pdl> fftnd $gal, $imag; # Perform Fourier transform
C performs a Fast Fourier Transform, in-place, on arbitrary-dimensioned data (i.e.
it is 'N-dimensional'). You can display C<$gal> after the FFT but you won't see
much. If at this point we ran C to invert it we would get the original
C<$gal> back.
If we want to enhance the large-scale structure we want to make a filter to only
let through low-frequencies:
pdl> $tmp = rvals($gal)<10; # Radially-symmetric filter function
pdl> use PDL::ImageND; # provides kernctr()
pdl> $filter = kernctr $tmp, $tmp; # Shift origin to 0,0
pdl> imag $filter;
=for html
You can see from the image that C<$filter> is zero everywhere except near the origin
(0,0) (and the 3 reflected corners). As a result it only lets through
low-frequency wavelengths. So we multiply by the filter and FFT back to
see the result:
pdl> ($gal2, $imag2) = cmul $gal, $imag, $filter, 0; # Complex multiplication
pdl> ifftnd $gal2, $imag2;
pdl> imag $gal2,0,300;
=for html
Figure: Fourier filtered smoothed image and contrast enhanced image with
the smoothed image subtracted.
Well that looks quite a bit different! Just about all the
high-frequency information has vanished. To see the high-frequency
information we can just subtract our filtered image from the original to
form the right hand image.
pdl> $orig = rfits "fixed_gal.fits";
pdl> imag $orig-$gal2,0,100;
=head2 Roundoff
Well that is probably enough abuse of Messier 51. We have demonstrated the ease
of simple and complex data processing with PDL and how PDL fits neatly in to
the Perl syntax as well as extending it. You have come across basic
arithmetical operations and a scattering of useful functions - and learned how
to find more. You certainly ought now to have a good feel of what PDL is all
about. In the next chapter we'll take a more comprehensive look at the basic
parts of PDL that all keen PDL users should know.