--- a/PDL/Book/FirstSteps.pod +++ b/PDL/Book/FirstSteps.pod @@ -37,26 +37,26 @@ For interactive use PDL comes with a program called C<pdl>. This allows you to type raw PDL (and perl) commands and see the result right away. It -also allows command line recall and editting (via the arrow keys) on most +also allows command line recall and editing (via the arrow keys) on most systems. So we begin by running the C<pdl> program from the system command line. On a Mac/UNIX/Linux system we would simply type C<pdl> in a 'terminal window'. On a Windows system we would type C<pdl> in a command prompt window. If PDL is installed correctly this is all that is required to bring up C<pdl>. - myhost% pdl - perlDL shell v1.354_001 - 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 /Users/kenworthy/.perldlrc... - Found docs database /Library/Perl/5.12/darwin-thread-multi-2level/PDL/pdldoc.db - Type 'help' for online help - Type 'demo' for online demos - Loaded PDL v2.4.9_009 (supports bad values) - pdl> + myhost% pdl + perlDL shell v1.354_001 + 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 /Users/kenworthy/.perldlrc... + Found docs database /Library/Perl/5.12/darwin-thread-multi-2level/PDL/pdldoc.db + Type 'help' for online help + Type 'demo' for online demos + Loaded PDL v2.4.9_009 (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 I<completely> configurable, @@ -67,9 +67,9 @@ Let's create something, and display it: - pdl> use PDL::Graphics::PGPLOT - pdl> imag (sin(rvals(200,200)+1)) - Displaying 200 x 200 image from -1 to 0.999999940395355 ... + pdl> use PDL::Graphics::PGPLOT + pdl> imag (sin(rvals(200,200)+1)) + Displaying 200 x 200 image from -1 to 0.999999940395355 ... The result should look like the image below a two dimensional C<sin> function. C<rvals> is a handy PDL function for creating an image whose pixel values are @@ -92,6 +92,7 @@ perl script? Well it is pretty simple: use PDL; + use PDL::Graphics::PGPLOT; imag (sin(rvals(200,200)+1)); That's it. This is a complete perl/PDL program. One could run it by typing @@ -134,12 +135,11 @@ which is one of the many formats standard PDL can read. (FITS stores more shades of grey than GIF or JPEG, PDL can also read these formats). - $a = rfits("PDL/Book/m51_raw.fits"); - - Reading IMAGE data... - BITPIX = -32 size = 262144 pixels - Reading 1048576 bytes - BSCALE = && BZERO = + pdl> $a = rfits("PDL/Book/m51_raw.fits"); + Reading IMAGE data... + BITPIX = -32 size = 262144 pixels + Reading 1048576 bytes + BSCALE = && BZERO = This looks pretty simple. As you can probably guess by now C<rfits> is the PDL function to read a FITS file. This is stored in the perl variable C<$a>. @@ -155,7 +155,7 @@ Before we start seriously playing around with M51 it is worth noting that we can also say: - $a = rfits "m51_raw.fits"; + pdl> $a = rfits "m51_raw.fits"; Note we have now left off the brackets on the C<rfits> function. Perl is rather simpler than C and allows one to omit the brackets on a function all together. @@ -167,7 +167,7 @@ Let's look at M51: - imag $a; + pdl> imag $a; =for html <img WIDTH=400 src="firststeps/whirl-m51a.png"> @@ -181,8 +181,8 @@ easily change that by specifying the black/white data values (Note: C<#> starts a Perl comment and can be ignored - i.e. no need to type the stuff after it!): - imag $a,0,1000; # More contrast - imag $a,0,300; # Even more contrast + pdl> imag $a,0,1000; # More contrast + pdl> imag $a,0,300; # Even more contrast You can see that C<imag> takes additional arguments to specify the display range. In fact C<imag> takes quite a few arguments, many of them optional. By @@ -196,8 +196,8 @@ artifact in digital imaging. We can correct for this using an image of a uniformly illuminated screen, what is commonly known as a 'flatfield'. - $flat = rfits "m51_flatfield.fits"; - imag $flat; + pdl> $flat = rfits "m51_flatfield.fits"; + pdl> imag $flat; This is shown in the next figure. Because the image is of a uniform field, the actual image reflects the detector sensitivity. To correct our M51 @@ -207,9 +207,9 @@ Figure: The 'flatfield' image showing the detector sensitivity of the raw data. - $gal = $a / $flat; - imag $gal,0,300; - wfits $gal, 'fixed_gal.fits'; # Save our work as a FITS file + pdl> $gal = $a / $flat; + pdl> imag $gal,0,300; + pdl> wfits $gal, 'fixed_gal.fits'; # Save our work as a FITS file Well that's a lot better. But think what we have just done. Both C<$a> and C<$flat> are I<images>, with 512 pixels by @@ -227,8 +227,8 @@ Figure: The M51 image corrected for the flatfield. - $funny = log(($gal/300)**2 - $gal/100 + 4); - imag $funny; # Surprise! + pdl> $funny = log(($gal/300)**2 - $gal/100 + 4); + pdl> imag $funny; # Surprise! Or on 1-D line piddles. On on 3-D cubic piddles. In fact piddles can support an infinite number of dimensions (though your computers memory won't). @@ -256,8 +256,8 @@ one data value corresponds to in physical units: - pdl> p 10**36/sum($gal) - 5.58159992096455e+28 + pdl> p 10**36/sum($gal) + 5.58159992096455e+28 This is also about 200 solar luminosities, (Note we have switched to using 'C<p>' @@ -272,15 +272,15 @@ the brightness in a circular aperture, we'll use the C<rvals> function: - $r = rvals $gal; - imag $r; + pdl> $r = rvals $gal; + pdl> imag $r; Remember C<rvals>? It replaces all the pixels in an image with its distance from the centre. We can turn this into a I<mask> with a simple operation like: - $mask = $r<50; - imag $mask; + pdl> $mask = $r<50; + pdl> imag $mask; =for html <img WIDTH=400 src="firststeps/whirl-maska.png"> =for html <img WIDTH=400 src="firststeps/whirl-maskb.png"> @@ -318,8 +318,8 @@ I<surface brightness>, the brightness per unit area on the sky, compare between bulge and disk? This is again quite straight forward: - print sum($bulge)/sum($mask); - print sum($disk)/sum(1-$mask); + 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 @@ -329,24 +329,24 @@ Of course PDL being so powerful, we could have figured this out in one line: - pdl> print average($gal->where(rvals($gal)<50)) - / average($gal->where(rvals($gal)>=50)) - 6.56590509414673 + 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: - $section = $gal(337:357,178:198); - imag $section; # the bright star + 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<extensions> to the Perl syntax. We can say C<$var(a:b,c:d...)> to specify I<multidimensional slices>. In this case we have produced a sub-image ranging -from pixel 327 to 367 along the first dimension, and 168 through 208 along the -second. We'll talk some more about I<slicing and dicing> later on. This -sub-image happens to contain a bright star. +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<slicing and dicing> 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%) @@ -361,17 +361,17 @@ the centroid automatically using the C<xvals> and C<yvals> functions). Well it is simple enough to get the distance from the centre: - $r = rvals $section; + 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<clump> function, which 'clumps' together an arbitrary number of dimensions: - $rr = $r->clump(2); # Clump first two dimensions - $sec = $section->clump(2); + pdl> $rr = $r->clump(2); # Clump first two dimensions + pdl> $sec = $section->clump(2); - points $rr, $sec; # Radial plot + pdl> points $rr, $sec; # Radial plot You should see a nice graph with points like those in Figure 1.6 @@ -395,7 +395,7 @@ more general fitting packages like C<PDL::Fit::LM> or C<PDL::Opt::Simplex> but that would require more care). - use PDL::Fit::Gaussian; + 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 @@ -403,16 +403,16 @@ C<help> tells us more about what we already know, to find out about what we don't know use C<apropos>: - perldl> 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. + 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<help> and C<apropos> @@ -421,9 +421,9 @@ Let's fit a Gaussian: - use PDL::Fit::Gaussian - ($peak, $fwhm, $background) = fitgauss1dr($rr, $sec); - print $peak, $fwhm, $background; + pdl> use PDL::Fit::Gaussian; + pdl> ($peak, $fwhm, $background) = fitgauss1dr($rr, $sec); + pdl> p $peak, $fwhm, $background; C<fitgauss1dr> is a function in the module L<PDL::Fit::Gaussian> which fits @@ -433,11 +433,11 @@ FWHM is more like 2.75 pixels. Let's generate a fitted curve with this functional form. - $rrr = sequence(2000)/100 # Generate radial values 0,0.01,0,02..20 + pdl> $rrr = sequence(2000)/100; # Generate radial values 0,0.01,0,02..20 # Generate gaussian with given FWHM - - $fit = $peak * exp(-2.772 * ($rrr/$fwhm)**2) + $background; + + pdl> $fit = $peak * exp(-2.772 * ($rrr/$fwhm)**2) + $background; Note the use of a new function, C<sequence(N)>, which @@ -446,8 +446,8 @@ Now let's overlay it on the previous plot. - hold; # This command stops new plots starting new pages - line $rrr, $fit, {Colour=>2} ; # Line plot + pdl> hold; # This command stops new plots starting new pages + pdl> line $rrr, $fit, {Colour=>2} ; # Line plot The last C<line> command shows @@ -462,10 +462,10 @@ Let's make a I<simulation> of the 2D stellar image. This is equally easy: - $fit2d = $peak * exp(-2.772 * ($r/$fwhm)**2); - release; # Back to new page for new plots; - imag $fit2d; - wfits $fit2d, 'fake_star.fits'; # Save our work + 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 that plot (Figure 1.7(a)) is a bit @@ -481,8 +481,8 @@ Figure 1.7: Two different views of the 2D simulated Point Spread Function. B<XXX need to mosaic the panels> - use PDL::Graphics::TriD; # Load the 3D graphics module - imag3d [$fit2d]; + 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 @@ -493,9 +493,9 @@ Finally here's something interesting. Let's take our fake star and place it elsewhere on the galaxy image. - $newsection = $gal(50:70,70:90); - $newsection += $fit2d; - imag $gal,0,300; + 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 @@ -510,7 +510,7 @@ 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: - $newsection = $newsection + $fit2d + 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 @@ -521,6 +521,7 @@ 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"; @@ -562,9 +563,9 @@ represents the strength of wavelengths of different scales in the image. Let's do it: - use PDL::FFT; # Load Fast Fourier Transform package - - $gal = rfits "fixed_gal.fits"; + 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 @@ -572,9 +573,9 @@ are represented in PDL by seperate real and imaginary arrays - more about this in Chapter 2.) - $imag = $gal * 0; # Create imaginary component, equal to zero - - fftnd $gal, $imag; # Perform Fourier transform + pdl> $imag = $gal * 0; # Create imaginary component, equal to zero + + pdl> fftnd $gal, $imag; # Perform Fourier transform C<fftnd> 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 @@ -584,18 +585,18 @@ If we want to enhance the large-scale structure we want to make a filter to only let through low-frequencies: - $tmp = rvals($gal)<10; # Radially-symmetric filter function - $filter = kernctr $tmp, $tmp; # Shift origin to 0,0 - imag $filter; + pdl> $tmp = rvals($gal)<10; # Radially-symmetric filter function + pdl> $filter = kernctr $tmp, $tmp; # Shift origin to 0,0 + pdl> imag $filter; 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 let's through low-frequency wavelengths. So we multiply by the filter and FFT back to see the result: - ($gal2, $imag2) = cmul $gal, $imag, $filter, 0; # Complex multiplication - ifftnd $gal2, $imag2; - imag $gal2,0,300; + pdl> ($gal2, $imag2) = cmul $gal, $imag, $filter, 0; # Complex multiplication + pdl> ifftnd $gal2, $imag2; + pdl> imag $gal2,0,300; =for html <img WIDTH=400 src="firststeps/whirl-ffta.png"> =for html <img WIDTH=400 src="firststeps/whirl-fftb.png"> @@ -609,8 +610,8 @@ subtract our filtered image from the original (Figure 1.9(b) - double yuk!): - $orig = rfits "fixed_gal.fits"; - imag $orig-$gal2,0,100; + pdl> $orig = rfits "fixed_gal.fits"; + pdl> imag $orig-$gal2,0,100; =head2 Roundoff