[17de4c]: content / FirstSteps.html Maximize Restore History

Download this file

FirstSteps.html    611 lines (570 with data), 30.9 kB

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
<head>
<title>The PDL Book</title>
<meta http-equiv="content-type" content="text/html; charset=utf-8" />
</head>
<body style="background-color: white">
<p>
</p>
<hr />
<h1><a name="first_steps_with_pdl">First Steps with PDL</a></h1>
<p><em>&quot;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.&quot;</em></p>
<p><em>- Annette Foglino, Space: Is Anyone Out There? Most astronomers say yes, Life, 1 Jul 1989.</em></p>
<p>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 &quot;Janet and John&quot;
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.</p>
<p>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.</p>
<p>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.</p>
<p>
</p>
<h2><a name="alright__let_s_do_something">Alright, let's <em>do</em> something</a></h2>
<p>We'll assume PDL is correctly installed and set up on
your computer system (see <a href="http://pdl.perl.org/">http://pdl.perl.org/</a> for details
of obtaining and installing PDL).</p>
<p>For interactive use PDL comes with a program called <code>perldl</code>. 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.</p>
<p>So we begin by running the <code>perldl</code> program from the system
command line. On a Mac/UNIX/Linux system we would simply type <code>perldl</code>
in a terminal window. On a Windows system we would type <code>perldl</code>
in a command prompt window. If PDL is installed correctly this is
all that is required to bring up <code>perldl</code>.</p>
<pre>
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&gt;</pre>
<p>We get a whole bunch of informational messages about what it is loading for
startup and the help system. Note; the startup is <em>completely</em> configurable,
an advanced user can completely customize which PDL modules
are loaded. We are left with the <code>pdl</code>&gt; prompt at which we can type commands. This kind
of interactive program is called a 'shell'. There is also <code>pdl2</code>
which is a newer version of the PDL shell with additional features.
It is still under development but completely usable.</p>
<p>Let's create something, and display it:</p>
<pre>
pdl&gt; use PDL::Graphics::Simple
pdl&gt; imag (sin(rvals(200,200)+1))</pre>
<p>The result should look like the image below - a two dimensional <code>sin</code> function.
<code>rvals</code> 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 '<code>imag(rvals(200,200))</code>' and you
will see better what we mean!) <code>sin()</code> 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 <code>imag()</code> 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.</p>
<table>
<caption align="bottom">Figure of a two dimensional <code>sin</code> function.</caption>
<tr><td><img WIDTH=300 src="images/sepia-small/whirl-sync.png"></td></tr>
</table>
<p>That's all fine. But what if we wanted to achieve the same results in a standalone
perl script? Well it is pretty simple:</p>
<pre>
use PDL;
use PDL::Graphics::Simple;
imag (sin(rvals(200,200)+1));</pre>
<p>That's it. This is a complete perl/PDL program. One could run it by typing
<code>perl filename</code>. (In fact there are many ways of running it, most systems
allows it to be setup so you can just type <em>filename</em>. See your local
Perl documentation - then the <code>perlrun</code> manual page.)</p>
<p>Two comments:</p>
<ol>
<li>
<p>The statements are all terminated by the '<code>;</code>' character. Perl is like C
in this regard. When entering code at the <code>pdl</code> command line the final
'<code>;</code>' 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
<code>pdl</code> prompt for clarity.</p>
</li>
<li>
<p>The directive <code>use PDL;</code> 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).</p>
</li>
</ol>
<p>
</p>
<h2><a name="whirling_through_the_whirlpool">Whirling through the Whirlpool</a></h2>
<p>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
<a href="http://sourceforge.net/projects/pdl/files/PDL/PDL%20Book%20Example%20Data%20Set/">http://sourceforge.net/projects/pdl/files/PDL/PDL%20Book%20Example%20Data%20Set/</a>
and we'll assume you are running any of these examples in the same
directory as you have downloaded the input data files.</p>
<p>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).</p>
<pre>
pdl&gt; $a = rfits(&quot;m51_raw.fits&quot;); # m51_raw.fits is in current directory
Reading IMAGE data...
BITPIX = -32 size = 262144 pixels
Reading 1048576 bytes
BSCALE = &amp;&amp; BZERO =</pre>
<p>This looks pretty simple. As you can probably guess by now <code>rfits</code> is the PDL
function to read a FITS file. This is stored in the perl variable <code>$a</code>.</p>
<p><strong>This is an important PDL concept: PDL stores its data arrays in simple perl
variables</strong> (<code>$a, $x, $y, $MyData</code>, etc.). PDL data arrays are special arrays
which use a more efficient, compact storage than standard perl arrays (<code>@a,
@x, ...</code>) 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
<em>piddles</em> (short for 'PDL variables') to distinguish them from ordinary Perl
'arrays', which are in fact really lists. We'll say more about this later.</p>
<p>Before we start seriously playing around with M51 it is worth noting that we
can also say:</p>
<pre>
pdl&gt; $a = rfits &quot;m51_raw.fits&quot;;</pre>
<p>Note we have now left off the brackets on the <code>rfits</code> 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</p>
<p>Let's look at M51:</p>
<pre>
pdl&gt; imag $a;
pdl&gt imag $a,0,1000; # More contrast
pdl&gt imag $a,0,300; # Even more contrast
</pre>
<p>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 linearly from
the faintest to the brightest pixel, and only the bright star slightly to the
bottom right of the center can be seen without contrast enhancement. We can
easily change that by specifying the black/white data values (Note: <code>#</code> starts
a Perl comment and can be ignored - i.e. no need to type the stuff after it!):</p>
<table>
<caption align="bottom">Figure of the raw image <code>m51_raw.fits</code> shown
with progressively greater contrast using the <code>imag</code> command.</caption>
<tr>
<td><img WIDTH=300 src="images/sepia-small/whirl-m51-default.png"></td>
<td><img WIDTH=300 src="images/sepia-small/whirl-m51-1000.png"></td>
<td><img WIDTH=300 src="images/sepia-small/whirl-m51-300.png"></td>
</tr>
</table>
<p>You can see that <code>imag</code> takes additional arguments to specify the display
range. In fact <code>imag</code> takes quite a few arguments, many of them optional. By
typing '<code>help imag</code>' at the <code>pdl</code> prompt we can find out all about the
function.</p>
<p>It is certainly a spiral galaxy with a few foreground stars thrown in for good
measure. But what is that horrible stripey pattern running from bottom right to
top left? That certainly is not part of the galaxy? Well no. What we have here
is the uneven sensitivity of the detector used to record the image, a common
artifact in digital imaging. We can correct for this using an image of a
uniformly illuminated screen, what is commonly known as a 'flatfield'.</p>
<pre>
pdl&gt; $flat = rfits &quot;m51_flatfield.fits&quot;;
pdl&gt; imag $flat;</pre>
<p>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
image, we merely have to divide the image by the flatfield:</p>
<table>
<caption align="bottom">Figure: The 'flatfield' image showing the detector sensitivity of the raw data.</caption>
<tr><td> <img WIDTH=300 src="images/sepia-small/whirl-flat.png"> </td></tr>
</table>
<pre>
pdl&gt; $gal = $a / $flat;
pdl&gt; imag $gal,0,300;
pdl&gt; wfits $gal, 'fixed_gal.fits'; # Save our work as a FITS file</pre>
<p>Well that's a lot better. But think what we have
just done. Both <code>$a</code> and <code>$flat</code> are <em>images</em>, with 512 pixels by
512 pixels. <strong>The divide operator '<code>/</code>' has been applied over all
262144 data values in the piddles <code>$a</code> and <code>$flat</code>.</strong> And it was
pretty fast too - these are what are known as <em>vectorized</em>
operations. In PDL each of these is implemented by heavily optimized
C code, which is what makes PDL very efficient for procession of
large chunks of data. If you did the same operation using normal
perl arrays rather than piddles it would be about ten to twenty times slower
(and use ten times more memory). In fact we can
do whatever arithmetic operations we like on image piddles:</p>
<table>
<caption align="bottom">Figure: The M51 image corrected for the flatfield.</caption>
<tr><td> <img WIDTH=300 src="images/sepia-small/whirl-flattened.png"> </td></tr>
</table>
<pre>
pdl&gt; $funny = log(($gal/300)**2 - $gal/100 + 4);
pdl&gt; imag $funny; # Surprise!</pre>
<p>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).</p>
<p><strong>This the key to PDL: the ability to process large chunks of data at once.</strong></p>
<p>
</p>
<h2><a name="measuring_the_brightness_of_m51">Measuring the brightness of M51</a></h2>
<p>How might we extract some useful
scientific information out of this image? A simple
quantity an astronomer might want to know is how the brightness of the
the 'disk' of the galaxy (the outer region which contains the spiral
arms) compares with the 'bulge' (the compact inner nucleus). Well
let's find out the total sum of all the light in the image:</p>
<pre>
pdl&gt; print sum($gal);
17916010</pre>
<p><code>sum</code> 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:</p>
<pre>
pdl&gt; p 10**36/sum($gal)
5.58159992096455e+28</pre>
<p>This is also about 200 solar luminosities, (Note we have switched to using <code>p</code>
as a shorthand for <code>print</code> - which only works in the <code>pdl</code> and <code>pdl2</code> shells)
which gives 4 billion solar luminosities for the whole galaxy.</p>
<p>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 <code>rvals</code>
function:</p>
<pre>
pdl&gt; $r = rvals $gal;
pdl&gt; imag $r;
...</pre>
<p>Remember <code>rvals</code>? It replaces all the pixels in an image with its distance
from the centre. We can turn this into a <em>mask</em> with a simple
operation like:</p>
<pre>
pdl&gt; $mask = $r&lt;50;
pdl&gt; imag $mask;
...</pre>
<table>
<caption align="bottom">Figure: Using <code>rvals</code> to generate a mask image
to isolate the galaxy bulge and disk. Top row: radial gradient
image <code>$r</code>, and radial gradient masked with less than operator
<code>$r < 50</code>. Bottom row: Bulge and disk of the galaxy.
</caption>
<tr>
<td> <img WIDTH=300 src="images/sepia-small/whirl-maska.png"> </td>
<td> <img WIDTH=300 src="images/sepia-small/whirl-maskb.png"> </td>
</tr>
<tr>
<td> <img WIDTH=300 src="images/sepia-small/whirl-maskc.png"> </td>
<td> <img WIDTH=300 src="images/sepia-small/whirl-maskd.png"> </td>
</tr>
</table>
<p>The Perl <em>less than operator</em> 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:</p>
<pre>
pdl&gt; $bulge = $mask * $gal
pdl&gt; imag $bulge,0,300
...
pdl&gt; print sum $bulge;
3011125
pdl&gt; $disk = $gal * (1-$mask)
pdl&gt; imag $disk,0,300
...
pdl&gt; print sum $disk
14904884</pre>
<p>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
<em>surface brightness</em>, the brightness per unit area on the sky,
compare between bulge and disk? This is again quite straight forward:</p>
<pre>
pdl&gt; print sum($bulge)/sum($mask);
pdl&gt; print sum($disk)/sum(1-$mask);</pre>
<p>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.</p>
<p>Of course PDL being so powerful, we could have figured this out in one line:</p>
<pre>
pdl&gt; print ( avg($gal-&gt;where(rvals($gal)&lt;50)) / avg($gal-&gt;where(rvals($gal)&gt;=50)) )
6.56590509414673</pre>
<p>
</p>
<h2><a name="twinkle__twinkle__little_star">Twinkle, twinkle, little star</a></h2>
<p>Let's look at something else, we'll zoom in on a small piece of the image:</p>
<pre>
pdl&gt; $section = $gal(337:357,178:198);
pdl&gt; imag $section; # the bright star</pre>
<p>Here we are introducing something new - we can see that PDL supports
<em>extensions</em> to the Perl syntax. We can say <code>$var(a:b,c:d...)</code> to specify
<em>multidimensional slices</em>. 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 <em>slicing and dicing</em> later on. This sub-image happens to contain
a bright star.</p>
<p>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.</p>
<p>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
centered on the star, you could think if you want about how you might determine
the centroid automatically using the <code>xvals</code> and <code>yvals</code> functions). Well it
is simple enough to get the distance from the centre:</p>
<pre>
pdl&gt; $r = rvals $section;</pre>
<p>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 <code>clump</code>
function, which 'clumps' together an arbitrary number of dimensions:</p>
<pre>
pdl&gt; $rr = $r-&gt;clump(2); # Clump first two dimensions
pdl&gt; $sec = $section-&gt;clump(2);
pdl&gt; points $rr, $sec; # Radial plot</pre>
<p>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 the plot
it looks like this is about 2-3 pixels - pretty compact!</p>
<table>
<caption align="bottom">Figure: Radial light profile of the bright star with fitted curve. </caption>
<tr> <td> <img WIDTH=300 src="images/sepia-small/whirl-star-radial.png"> </td> </tr>
</table>
<p>Well we don't just want a guess - let's fit the profile with a function.
These blurring functions are usually represented by the <code>Gaussian</code>
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 surprisingly
enough a special function for this very purpose. (One could use
more general fitting packages like <code>PDL::Fit::LM</code> or
<code>PDL::Opt::Simplex</code> but that would require more care).</p>
<pre>
pdl&gt; use PDL::Fit::Gaussian;</pre>
<p>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
<code>help</code> tells us more about what we already know, to find out
about what we don't know use <code>apropos</code>:</p>
<pre>
pdl&gt; 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.</pre>
<p>This tells us a whole lot about various functions and modules to do with
Gaussians. Note that we can abbreviate <code>help</code> and <code>apropos</code>
with '<code>?</code>' and '<code>??</code>' when using the <code>pdl</code> or <code>pdl2</code> shells.</p>
<p>Let's fit a Gaussian:</p>
<pre>
pdl&gt; use PDL::Fit::Gaussian;
pdl&gt; ($peak, $fwhm, $background) = fitgauss1dr($rr, $sec);
pdl&gt; p $peak, $fwhm, $background;</pre>
<p><code>fitgauss1dr</code> is a function in the module <a href="/PDL/Fit/Gaussian.html">the PDL::Fit::Gaussian manpage</a> 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.</p>
<pre>
pdl&gt; $rrr = sequence(2000)/100; # Generate radial values 0,0.01,0,02..20
# Generate Gaussian with given FWHM</pre>
<pre>
pdl&gt; $fit = $peak * exp(-2.772 * ($rrr/$fwhm)**2) + $background;</pre>
<p>Note the use of a new function, <code>sequence(N)</code>, 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.</p>
<pre>
pdl&gt; hold; # This command stops new plots starting new pages
pdl&gt; line $rrr, $fit, {Colour=&gt;2} ; # Line plot</pre>
<p>The last <code>line</code> 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 <a href="/PDL/Book/PGPLOT.html">the PDL::Book::PGPLOT manpage</a>. 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.</p>
<p>Let's make a <em>simulation</em> of the 2D stellar image. This is equally
easy:</p>
<pre>
pdl&gt; $fit2d = $peak * exp(-2.772 * ($r/$fwhm)**2);
pdl&gt; release; # Back to new page for new plots;
pdl&gt; imag $fit2d;
...
pdl&gt; wfits $fit2d, 'fake_star.fits'; # Save our work</pre>
<p>But the figure below is a
boring. So far we have been using simple 2D graphics from the
<code>PDL::Graphics::Simple</code> library. In fact PDL has more
than one graphics library (some see this as a flaw, some
as a feature!). Using the <code>PDL::Graphics::TriD</code> library
which does OpenGL graphics we can look at our simulated
star in 3D (see the right hand panel);</p>
<table>
<caption align="bottom">Figure: Two different views of the 2D simulated Point Spread Function.</caption>
<tr>
<td> <img WIDTH=300 src="images/sepia-small/whirl-starsima.png"> </td>
<td> <img WIDTH=300 src="images/sepia-small/whirl-starsimb.png"> </td>
</tr>
</table>
<pre>
pdl&gt; use PDL::Graphics::TriD; # Load the 3D graphics module
pdl&gt; imag3d [$fit2d];</pre>
<p>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 <code>imag3d</code> has it's
a rather different syntax for processing it's arguments - for very good reasons
- we'll explore 3D graphics further in <a href="/PDL/Book/TriD.html">the PDL::Book::TriD manpage</a>.</p>
<pre>
To continue: Select the TriD window and type q</pre>
<p>Finally here's something interesting. Let's take our fake star and place it
elsewhere on the galaxy image.</p>
<pre>
pdl&gt; $newsection = $gal(50:70,70:90);
pdl&gt; $newsection += $fit2d;
pdl&gt; imag $gal,0,300;</pre>
<p>We have a bright new star where none existed before! The C-style <code>+=</code>
increment operator is worth noting - it actually modifies the contents of
<code>$newsection</code> in-place. And because <code>$newsection</code> is a <em>slice</em> of <code>$gal</code>
the change also affects <code>$gal</code>. This is an important property of slices - any
change to the slice affects the <em>parent</em>. 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.</p>
<p>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:</p>
<pre>
pdl&gt; $newsection = $newsection + $fit2d</pre>
<p>Now a new version of <code>$newsection</code> is created which has nothing to
do with the original <code>$gal</code>. In fact there is more than one way to do
this as we will see in later chapters.</p>
<p>Just to amuse ourselves, lets write a short script to cover M51 with dozens of fake
stars of random brightnesses:</p>
<pre>
use PDL;
use PDL::Graphics::Simple;
use PDL::NiceSlice; # must use in each program file</pre>
<pre>
srand(42); # Set the random number seed
$gal = rfits &quot;fixed_gal.fits&quot;;
$star = rfits &quot;fake_star.fits&quot;;
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);
}</pre>
<pre>
for (1..100) {
$x1 = int(rand(470)+10);
$y1 = int(rand(470)+10);
addstar($x1,$y1);
}
imag $gal,0,1000;</pre>
<p>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.</p>
<table>
<caption align="bottom">Figure: M51 covered in fake stars.</caption>
<tr> <td> <img WIDTH=300 src="images/sepia-small/whirl-fakestars.png"> </td> </tr>
</table>
<p>
</p>
<h2><a name="getting_complex_with_m51">Getting Complex with M51</a></h2>
<p>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.</p>
<p>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:</p>
<pre>
pdl&gt; use PDL::FFT; # Load Fast Fourier Transform package
pdl&gt; $gal = rfits &quot;fixed_gal.fits&quot;;</pre>
<p>Now <code>$gal</code> contains real values, to do the Fourier transform it has to
have complex values. We create a variable <code>$imag</code> to hold the imaginary
component and set to zero.(For reasons of efficiency complex numbers
are represented in PDL by separate real and imaginary arrays - more about this
in Chapter 2.)</p>
<pre>
pdl&gt; $imag = $gal * 0; # Create imaginary component, equal to zero
pdl&gt; fftnd $gal, $imag; # Perform Fourier transform</pre>
<p><code>fftnd</code> performs a Fast Fourier Transform, in-place, on arbitrary-dimensioned data (i.e.
it is 'N-dimensional'). You can display <code>$gal</code> after the FFT but you won't see
much. If at this point we ran <code>ifftnd</code> to invert it we would get the original
<code>$gal</code> back.</p>
<p>If we want to enhance the large-scale structure we want to make a filter to only
let through low-frequencies:</p>
<pre>
pdl&gt; $tmp = rvals($gal)&lt;10; # Radially-symmetric filter function
pdl&gt; use PDL::ImageND; # provides kernctr()
pdl&gt; $filter = kernctr $tmp, $tmp; # Shift origin to 0,0
pdl&gt; imag $filter;</pre>
<table>
<caption align="bottom"> Figure: The result of <code>kernctr()</code> </caption>
<tr> <td> <img WIDTH=300 src="images/sepia-small/gal-filter.png"> </td> </tr>
</table>
<p>You can see from the image that <code>$filter</code> is zero everywhere except near the origin <code>(0,0)</code> (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 (<code>cmul</code> is complex multiplication):</p>
<pre>
pdl&gt; ($gal2, $imag2) = cmul $gal, $imag, $filter, 0;
pdl&gt; ifftnd $gal2, $imag2;
pdl&gt; imag $gal2,0,300;</pre>
<table>
<caption align="bottom">Figure: Fourier filtered smoothed image
and contrast enhanced image with the smoothed image subtracted. </caption>
<tr>
<td> <img WIDTH=300 src="images/sepia-small/whirl-ffta.png"> </td>
<td> <img WIDTH=300 src="images/sepia-small/whirl-fftb.png"> </td>
</tr>
</table>
<p>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.</p>
<pre>
pdl&gt; $orig = rfits &quot;fixed_gal.fits&quot;;
pdl&gt; imag $orig-$gal2,0,100;</pre>
<p>
</p>
<h2><a name="roundoff">Roundoff</a></h2>
<p>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.</p>
<h3><a name="pdlbook">Where to Read More</a></h3>
<p> We hope this excerpt from the
<a href="http://sourceforge.net/projects/pdl/files/PDL_2013/PDL-Book/">PDL Book</a>
has given you a taste of what PDL is capable of. Please download the full text to
find out more.</p>
<p>As always your one-stop-shop for all things PDL is at
<a href="http://pdl.perl.org">http://pdl.perl.org</a>.</p>
</body>
</html>