Diff of /content/FirstSteps.html [000000] .. [17de4c]  Maximize  Restore

Switch to unified view

a b/content/FirstSteps.html
1
<head>
2
<title>The PDL Book</title>
3
<meta http-equiv="content-type" content="text/html; charset=utf-8" />
4
</head>
5
6
<body style="background-color: white">
7
8
<p>
9
</p>
10
<hr />
11
<h1><a name="first_steps_with_pdl">First Steps with PDL</a></h1>
12
<p><em>&quot;Maybe there are a few civilizations out there that have decided to stay
13
home, piddle around and send out some radio waves once in a while.&quot;</em></p>
14
<p><em>- Annette Foglino, Space: Is Anyone Out There? Most astronomers say yes, Life, 1 Jul 1989.</em></p>
15
<p>It can be very frustrating to read an introductory book which takes a
16
long time teaching you the very basics of a topic, in a &quot;Janet and John&quot;
17
style.  While you wish to learn, you are anxious to see something a bit
18
more exciting and interesting to see what the language can do.</p>
19
<p>Fortunately our task in this book on PDL is made very much easier by the
20
high-level of the language. We can take a tour through PDL, looking at
21
the advanced features it offers without getting involved in complexity.</p>
22
<p>The aim of this section is to cover a breadth of PDL features rather
23
than any in depth, to give the reader a flavour of what he or she can do
24
using the language and a useful reference for getting started doing real
25
work. Later sections will focus on looking at the features introduced
26
here, in more depth.</p>
27
<p>
28
</p>
29
<h2><a name="alright__let_s_do_something">Alright, let's <em>do</em> something</a></h2>
30
<p>We'll assume PDL is correctly installed and set up on
31
your computer system (see <a href="http://pdl.perl.org/">http://pdl.perl.org/</a> for details
32
of obtaining and installing PDL).</p>
33
<p>For interactive use PDL comes with a program called <code>perldl</code>. This allows
34
you to type raw PDL (and perl) commands and see the result right away. It
35
also allows command line recall and editing (via the arrow keys) on most
36
systems.</p>
37
<p>So we begin by running the <code>perldl</code> program from the system
38
command line. On a Mac/UNIX/Linux system we would simply type <code>perldl</code>
39
in a terminal window. On a Windows system we would type <code>perldl</code>
40
in a command prompt window. If PDL is installed correctly this is
41
all that is required to bring up <code>perldl</code>.</p>
42
<pre>
43
  myhost% perldl
44
  perlDL shell v1.357
45
   PDL comes with ABSOLUTELY NO WARRANTY. For details, see the file
46
   'COPYING' in the PDL distribution. This is free software and you
47
   are welcome to redistribute it under certain conditions, see
48
   the same file for details.
49
  ReadLines, NiceSlice, MultiLines  enabled
50
  Reading PDL/default.perldlrc...
51
  Found docs database /usr/lib/perl5/.../PDL/pdldoc.db
52
  Type 'help' for online help
53
  Type 'demo' for online demos
54
  Loaded PDL v2.006 (supports bad values)
55
  pdl&gt;</pre>
56
<p>We get a whole bunch of informational messages about what it is loading for
57
startup and the help system. Note; the startup is <em>completely</em> configurable,
58
an advanced user can completely customize which PDL modules
59
are loaded. We are left with the <code>pdl</code>&gt; prompt at which we can type commands. This kind
60
of interactive program is called a 'shell'.  There is also <code>pdl2</code>
61
which is a newer version of the PDL shell with additional features.
62
It is still under development but completely usable.</p>
63
<p>Let's create something, and display it:</p>
64
<pre>
65
  pdl&gt; use PDL::Graphics::Simple
66
  pdl&gt; imag (sin(rvals(200,200)+1))</pre>
67
<p>The result should look like the image below - a two dimensional <code>sin</code> function.
68
<code>rvals</code> is a handy PDL function for creating an image whose pixel values are
69
the radial distance from the central pixel of the image.  With these arguments
70
it creates a 200 by 200 'radial' image.  (Try '<code>imag(rvals(200,200))</code>' and you
71
will see better what we mean!) <code>sin()</code> is the mathematical sine function, this
72
already exists in perl but in the case of PDL is applied to all 40000 pixels at
73
once, a topic we will come back to. The <code>imag()</code> function displays the image.
74
You will see the syntax of perl/PDL is algebraic - by which we mean it is very
75
similar to C and FORTRAN in how expressions are constructed.  (In fact much
76
more like C than FORTRAN). It is interesting to reflect on how much C code
77
would be required to generate the same display, even given the existence of
78
some convenient graphics library.</p>
79
<table>
80
  <caption align="bottom">Figure of a two dimensional <code>sin</code> function.</caption>
81
  <tr><td><img WIDTH=300 src="images/sepia-small/whirl-sync.png"></td></tr>
82
</table>
83
  <p>That's all fine. But what if we wanted to achieve the same results in a standalone
84
perl script? Well it is pretty simple:</p>
85
<pre>
86
  use PDL;
87
  use PDL::Graphics::Simple;
88
  imag (sin(rvals(200,200)+1));</pre>
89
<p>That's it. This is a complete perl/PDL program. One could run it by typing 
90
<code>perl filename</code>. (In fact there are many ways of running it, most systems
91
allows it to be setup so you can just type <em>filename</em>. See your local
92
Perl documentation - then the <code>perlrun</code> manual page.)</p>
93
<p>Two comments:</p>
94
<ol>
95
<li>
96
<p>The statements are all terminated by the '<code>;</code>' character. Perl is like C
97
in this regard. When entering code at the <code>pdl</code> command line the final
98
'<code>;</code>' may be omitted if you wish, note you can also use it to put multiple
99
statements on one line. In our examples from now on we'll often omit the 
100
<code>pdl</code> prompt for clarity.</p>
101
</li>
102
<li>
103
<p>The directive <code>use PDL;</code> tells Perl to load the PDL module, which makes
104
available all the standard PDL extensions. (Advanced users 
105
will be interested in knowing
106
there are other ways of starting PDL which allows one to select which bits
107
of it you want).</p>
108
</li>
109
</ol>
110
<p>
111
</p>
112
<h2><a name="whirling_through_the_whirlpool">Whirling through the Whirlpool</a></h2>
113
<p>Enough about the mechanics of using PDL, let's look at some real data! To work
114
through these examples exactly you can download any needed input files from
115
<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>
116
and we'll assume you are running any of these examples in the same
117
directory as you have downloaded the input data files.</p>
118
<p>We'll be playing with an image of the famous spiral galaxy discovered by
119
Charles Messier, known to astronomers as M51 and commonly as the Whirlpool
120
Galaxy. This is a 'nearby' galaxy, a mere 25 million light years from Earth.
121
The image file is stored in the 'FITS' format, a common astronomical format,
122
which is one of the many formats standard PDL can read. (FITS stores more
123
shades of gray than GIF or JPEG, but PDL can read these formats too).</p>
124
<pre>
125
  pdl&gt; $a = rfits(&quot;m51_raw.fits&quot;);   # m51_raw.fits is in current directory
126
  Reading IMAGE data...
127
  BITPIX =  -32  size = 262144 pixels 
128
  Reading  1048576  bytes
129
  BSCALE =  &amp;&amp;  BZERO =</pre>
130
<p>This looks pretty simple. As you can probably guess by now <code>rfits</code> is the PDL
131
function to read a FITS file. This is stored in the perl variable <code>$a</code>.</p>
132
<p><strong>This is an important PDL concept: PDL stores its data arrays in simple perl
133
variables</strong> (<code>$a, $x, $y, $MyData</code>, etc.).  PDL data arrays are special arrays
134
which use a more efficient, compact storage than standard perl arrays (<code>@a,
135
@x, ...</code>) and are much faster to access for numerical computations. To avoid
136
confusion it is convenient to introduce a special name for them, we call them
137
<em>piddles</em> (short for 'PDL variables') to distinguish them from ordinary Perl
138
'arrays', which are in fact really lists.  We'll say more about this later.</p>
139
<p>Before we start seriously playing around with M51 it is worth noting that we
140
can also say:</p>
141
<pre>
142
  pdl&gt; $a = rfits &quot;m51_raw.fits&quot;;</pre>
143
<p>Note we have now left off the brackets on the <code>rfits</code> function. Perl is rather
144
simpler than C and allows one to omit the brackets on a function all together.
145
It assumes all the items in a list are function arguments and can be pretty
146
convenient. If you are calling more than one function it is however better to
147
use some brackets so the meaning is clear. For the rules on this 'list
148
operator' syntax see the Perl syntax documentation.  From now on we'll mostly
149
use the list operator syntax for conciseness</p>
150
<p>Let's look at M51:</p>
151
<pre>
152
  pdl&gt; imag $a;
153
  pdl&gt imag $a,0,1000; # More contrast
154
  pdl&gt imag $a,0,300;  # Even more contrast
155
</pre>
156
157
<p>A couple of bright spots can be seen, but where is the galaxy? It's the faint
158
blob in the middle: by default the display range is autoscaled linearly from
159
the faintest to the brightest pixel, and only the bright star slightly to the
160
bottom right of the center can be seen without contrast enhancement. We can
161
easily change that by specifying the black/white data values (Note: <code>#</code> starts
162
a Perl comment and can be ignored - i.e. no need to type the stuff after it!):</p>
163
164
<table>
165
  <caption align="bottom">Figure of the raw image <code>m51_raw.fits</code> shown
166
      with progressively greater contrast using the <code>imag</code> command.</caption>
167
  <tr>
168
      <td><img WIDTH=300 src="images/sepia-small/whirl-m51-default.png"></td>
169
      <td><img WIDTH=300 src="images/sepia-small/whirl-m51-1000.png"></td>
170
      <td><img WIDTH=300 src="images/sepia-small/whirl-m51-300.png"></td>
171
  </tr>
172
</table>
173
174
175
176
<p>You can see that <code>imag</code> takes additional arguments to specify the display
177
range. In fact <code>imag</code> takes quite a few arguments, many of them optional. By
178
typing '<code>help imag</code>' at the <code>pdl</code> prompt we can find out all about the
179
function.</p>
180
<p>It is certainly a spiral galaxy with a few foreground stars thrown in for good
181
measure. But what is that horrible stripey pattern running from bottom right to
182
top left? That certainly is not part of the galaxy? Well no. What we have here
183
is the uneven sensitivity of the detector used to record the image, a common
184
artifact in digital imaging. We can correct for this using an image of a
185
uniformly illuminated screen, what is commonly known as a 'flatfield'.</p>
186
<pre>
187
  pdl&gt; $flat = rfits &quot;m51_flatfield.fits&quot;;
188
  pdl&gt; imag $flat;</pre>
189
<p>This is shown in the next figure. Because the image is of a uniform field,
190
the actual image reflects the detector sensitivity. To correct our M51
191
image, we merely have to divide the image by the flatfield:</p>
192
<table>
193
  <caption align="bottom">Figure: The 'flatfield' image showing the detector sensitivity of the raw data.</caption>
194
  <tr><td> <img WIDTH=300 src="images/sepia-small/whirl-flat.png"> </td></tr>
195
</table>
196
197
<pre>
198
  pdl&gt; $gal = $a / $flat;
199
  pdl&gt; imag $gal,0,300;  
200
  pdl&gt; wfits $gal, 'fixed_gal.fits'; # Save our work as a FITS file</pre>
201
202
<p>Well that's a lot better.  But think what we have
203
just done. Both <code>$a</code>  and <code>$flat</code> are <em>images</em>, with 512 pixels by
204
512 pixels. <strong>The divide operator '<code>/</code>' has been applied over all
205
262144 data values in the piddles <code>$a</code>  and <code>$flat</code>.</strong> And it was
206
pretty fast too - these are what are known as <em>vectorized</em>
207
operations. In PDL each of these is implemented by heavily optimized
208
C code, which is what makes PDL very efficient for procession of
209
large chunks of data. If you did the same operation using normal
210
perl arrays rather than piddles it would be about ten to twenty times slower
211
(and use ten times more memory). In fact we can
212
do whatever arithmetic operations we like on image piddles:</p>
213
<table>
214
  <caption align="bottom">Figure: The M51 image corrected for the flatfield.</caption>
215
  <tr><td> <img WIDTH=300 src="images/sepia-small/whirl-flattened.png"> </td></tr>
216
</table>
217
218
<pre>
219
  pdl&gt; $funny = log(($gal/300)**2 - $gal/100  + 4); 
220
  pdl&gt; imag $funny; # Surprise!</pre>
221
222
<p>Or on 1-D line piddles. On on 3-D cubic piddles. In fact piddles can support an infinite
223
number of dimensions (though your computers memory won't).</p>
224
<p><strong>This the key to PDL: the ability to process large chunks of data at once.</strong></p>
225
<p>
226
</p>
227
<h2><a name="measuring_the_brightness_of_m51">Measuring the brightness of M51</a></h2>
228
<p>How might we extract some useful
229
scientific information out of this image? A simple
230
quantity an astronomer might want to know is how the brightness of the
231
the 'disk' of the galaxy (the outer region which contains the spiral
232
arms) compares with the 'bulge' (the compact inner nucleus). Well
233
let's find out the total sum of all the light in the image:</p>
234
<pre>
235
  pdl&gt; print sum($gal);
236
  17916010</pre>
237
<p><code>sum</code> just sums up all the data values in all the pixels in the 
238
image - in this case the answer is 17916010. If the image is linear
239
(which it is) and if it was calibrated (i.e. we knew the relation
240
between data numbers and brightness units) we could work out the
241
total brightness. Let's turn it round - we know that M51 has
242
a luminosity of about 1E36 Watts, so we can work out what
243
one data value corresponds to in physical units:</p>
244
<pre>
245
  pdl&gt; p 10**36/sum($gal)
246
  5.58159992096455e+28</pre>
247
<p>This is also about 200 solar luminosities, (Note we have switched to using <code>p</code>
248
as a shorthand for <code>print</code> - which only works in the <code>pdl</code> and <code>pdl2</code> shells)
249
which gives 4 billion solar luminosities for the whole galaxy.</p>
250
<p>OK we do not need PDL for this simple arithmetic, let's get back to
251
computations that involve the whole image.
252
How can we get the sum of a piece of an image, e.g. near the centre? 
253
Well in PDL there is more than one way to do it (Perl aficionados call
254
this phenomenon TIMTOWTDI). In this case, because we really want
255
the brightness in a circular aperture, we'll use the <code>rvals</code>
256
function:</p>
257
<pre>
258
  pdl&gt; $r = rvals $gal;
259
  pdl&gt; imag $r;
260
  ...</pre>
261
<p>Remember <code>rvals</code>? It replaces all the pixels in an image with its distance
262
from the centre. We can turn this into a <em>mask</em> with a simple
263
operation like:</p>
264
<pre>
265
  pdl&gt; $mask = $r&lt;50;
266
  pdl&gt; imag $mask;
267
  ...</pre>
268
269
270
<table>
271
  <caption align="bottom">Figure: Using  <code>rvals</code> to generate a mask image
272
      to isolate the galaxy bulge and disk.  Top row: radial gradient
273
      image <code>$r</code>, and radial gradient masked with less than operator
274
      <code>$r < 50</code>.  Bottom row: Bulge and disk of the galaxy.
275
  </caption>
276
  <tr>
277
      <td> <img WIDTH=300 src="images/sepia-small/whirl-maska.png"> </td>
278
      <td> <img WIDTH=300 src="images/sepia-small/whirl-maskb.png"> </td>
279
  </tr>
280
  <tr>
281
      <td> <img WIDTH=300 src="images/sepia-small/whirl-maskc.png"> </td>
282
      <td> <img WIDTH=300 src="images/sepia-small/whirl-maskd.png"> </td>
283
  </tr>
284
</table>
285
286
287
<p>The Perl <em>less than operator</em> is applied to all pixels in the image.
288
You can see the result is an image which is 0 on the outskirts and 1 in
289
the area of the nucleus. We can then simply use the mask image to
290
isolate in a simple way the bulge and disk components (lower row) and it
291
is then very easy to find the brightness of both pieces of the M51
292
galaxy:</p>
293
<pre>
294
 pdl&gt; $bulge = $mask * $gal
295
 pdl&gt; imag $bulge,0,300
296
 ...
297
 pdl&gt; print sum $bulge;
298
 3011125
299
 
300
 pdl&gt; $disk = $gal * (1-$mask)
301
 pdl&gt; imag $disk,0,300
302
 ...
303
 pdl&gt; print sum $disk
304
 14904884</pre>
305
<p>You can see that the disk is about 5 times brighter than the bulge in
306
total, despite its more diffuse appearance. This is typical for
307
spiral galaxies. We might ask a different question: how does the average
308
<em>surface brightness</em>, the brightness per unit area on the sky,
309
compare between bulge and disk? This is again quite straight forward:</p>
310
<pre>
311
  pdl&gt; print sum($bulge)/sum($mask);
312
  pdl&gt; print sum($disk)/sum(1-$mask);</pre>
313
<p>We work out the area by simply summing up the 0,1 pixels in the mask
314
image. The answer is the bulge has about 7 times the surface
315
brightness than the disk - something we might have guessed from
316
looking at the above figure, which tells astronomers its stellar density is
317
much higher.</p>
318
<p>Of course PDL being so powerful, we could have figured this out in one line:</p>
319
<pre>
320
  pdl&gt; print ( avg($gal-&gt;where(rvals($gal)&lt;50)) / avg($gal-&gt;where(rvals($gal)&gt;=50)) )
321
  6.56590509414673</pre>
322
<p>
323
</p>
324
<h2><a name="twinkle__twinkle__little_star">Twinkle, twinkle, little star</a></h2>
325
<p>Let's look at something else, we'll zoom in on a small piece of the image:</p>
326
<pre>
327
  pdl&gt; $section = $gal(337:357,178:198);
328
  pdl&gt; imag $section; # the bright star</pre>
329
<p>Here we are introducing something new - we can see that PDL supports
330
<em>extensions</em> to the Perl syntax. We can say <code>$var(a:b,c:d...)</code> to specify
331
<em>multidimensional slices</em>. In this case we have produced a sub-image ranging
332
from pixel 337 to 357 along the first dimension, and 178 through 198 along the
333
second. Remember pdl data dimension indexes start from zero.  We'll talk some
334
more about <em>slicing and dicing</em> later on. This sub-image happens to contain
335
a bright star.</p>
336
<p>At this point you will probably be able to work out for yourself the amount of
337
light coming from this star, compared to the whole galaxy. (Answer: about 2%)
338
But let's look at something more involved: the radial profile of the star.
339
Since stars are a long way away they are almost point sources, but our camera
340
will blur them out into little disks, and for our analysis we might want an
341
exact figure for this blurring.</p>
342
<p>We want to plot all the brightness of all the pixels in this section, against
343
the distance from the centre. (We've chosen the section to be conveniently
344
centered on the star, you could think if you want about how you might determine
345
the centroid automatically using the <code>xvals</code> and <code>yvals</code> functions). Well it
346
is simple enough to get the distance from the centre:</p>
347
<pre>
348
  pdl&gt; $r = rvals $section;</pre>
349
<p>But to produce a one-dimensional plot of one against the other we need to
350
reduce the 2D data arrays to one dimension. (i.e our 21 by 21 image section
351
becomes a 441 element vector). This can be done using the PDL <code>clump</code>
352
function, which 'clumps' together an arbitrary number of dimensions:</p>
353
<pre>
354
  pdl&gt; $rr  = $r-&gt;clump(2); # Clump first two dimensions
355
  pdl&gt; $sec = $section-&gt;clump(2);
356
  
357
  pdl&gt; points $rr, $sec;  # Radial plot</pre>
358
<p>You should see a nice graph with points like those
359
in the figure below showing the drop-off from the bright centre of the star.
360
The blurring is usually measured
361
by the 'Full Width Half Maximum' (FWHM) - or in plain terms how
362
fat the profile is across when it drops by half. Looking at the plot
363
it looks like this is about 2-3 pixels - pretty compact!</p>
364
365
<table>
366
  <caption align="bottom">Figure: Radial light profile of the bright star with fitted curve. </caption>
367
  <tr> <td> <img WIDTH=300 src="images/sepia-small/whirl-star-radial.png"> </td> </tr>
368
</table>
369
370
371
<p>Well we don't just want a guess - let's fit the profile with a function.
372
These blurring functions are usually represented by the <code>Gaussian</code>
373
function. PDL comes with a whole variety of general purpose and
374
special purpose fitting functions which people have written for
375
their own purposes (and so will you we hope!). Fitting Gaussians
376
is something that happens rather a lot and there is surprisingly
377
enough a special function for this very purpose. (One could use
378
more general fitting packages like <code>PDL::Fit::LM</code> or
379
<code>PDL::Opt::Simplex</code> but that would require more care).</p>
380
<pre>
381
  pdl&gt; use PDL::Fit::Gaussian;</pre>
382
<p>This loads in the module to do this. PDL, like Perl, is modular. We
383
don't load all the available modules by default just a convenient
384
subset. How can we find useful PDL functions and modules? Well
385
<code>help</code> tells us more about what we already know, to find out
386
about what we don't know use <code>apropos</code>:</p>
387
<pre>
388
  pdl&gt; apropos gaussian
389
  PDL::Fit::Gaussian ...
390
                                Module: routines for fitting gaussians
391
  PDL::Gaussian   Module: Gaussian distributions.
392
  fitgauss1d      Fit 1D Gassian to data piddle
393
  fitgauss1dr     Fit 1D Gassian to radial data piddle
394
  gefa            Factor a matrix using Gaussian elimination.
395
  grandom         Constructor which returns piddle of Gaussian random numbers
396
  ndtri           The value for which the area under the Gaussian probability density function (integrated from minus
397
                                infinity) is equal to the argument (cf erfi). Works inplace.</pre>
398
<p>This tells us a whole lot about various functions and modules to do with
399
Gaussians. Note that we can abbreviate <code>help</code> and <code>apropos</code>
400
with '<code>?</code>' and '<code>??</code>' when using the <code>pdl</code> or <code>pdl2</code> shells.</p>
401
<p>Let's fit a Gaussian:</p>
402
<pre>
403
  pdl&gt; use PDL::Fit::Gaussian;
404
  pdl&gt; ($peak, $fwhm, $background) = fitgauss1dr($rr, $sec); 
405
  pdl&gt; p $peak, $fwhm, $background;</pre>
406
<p><code>fitgauss1dr</code> is a function in the module <a href="/PDL/Fit/Gaussian.html">the PDL::Fit::Gaussian manpage</a> which fits
407
a Gaussian constrained to be radial (i.e. whose peak is at the origin).
408
You can see that, unlike C and FORTRAN, Perl functions can return
409
more than one result value. This is pretty convenient. You can see the
410
FWHM is more like 2.75 pixels. Let's generate a fitted curve with this
411
functional form.</p>
412
<pre>
413
  pdl&gt; $rrr = sequence(2000)/100;  # Generate radial values 0,0.01,0,02..20
414
       
415
  # Generate Gaussian with given FWHM</pre>
416
<pre>
417
  pdl&gt; $fit = $peak * exp(-2.772 * ($rrr/$fwhm)**2) + $background;</pre>
418
<p>Note the use of a new function, <code>sequence(N)</code>, which 
419
generates a new piddle with N values ranging 0..(N-1).
420
We are simply using this to generate the horizontal axis values
421
for the plot.  Now let's overlay it on the previous plot.</p>
422
<pre>
423
  pdl&gt; hold; # This command stops new plots starting new pages
424
  pdl&gt; line $rrr, $fit, {Colour=&gt;2} ; # Line plot</pre>
425
<p>The last <code>line</code> command shows the PDL syntax for optional function
426
arguments. This is based on the Perl's built in hash syntax. We'll say
427
more about this later in <a href="/PDL/Book/PGPLOT.html">the PDL::Book::PGPLOT manpage</a>.  The result should look a
428
lot like the figure above.  Not too bad. We could perhaps do a bit
429
better by exactly centroiding the image but it will do for now.</p>
430
<p>Let's make a <em>simulation</em> of the 2D stellar image. This is equally
431
easy:</p>
432
<pre>
433
  pdl&gt; $fit2d = $peak * exp(-2.772 * ($r/$fwhm)**2);
434
  pdl&gt; release; # Back to new page for new plots;
435
  pdl&gt; imag $fit2d;
436
  ...
437
  pdl&gt; wfits $fit2d, 'fake_star.fits'; # Save our work</pre>
438
<p>But the figure below is a
439
boring. So far we have been using simple 2D graphics from the
440
<code>PDL::Graphics::Simple</code> library. In fact PDL has more
441
than one graphics library (some see this as a flaw, some
442
as a feature!). Using the <code>PDL::Graphics::TriD</code> library
443
which does OpenGL graphics we can look at our simulated
444
star in 3D (see the right hand panel);</p>
445
446
<table>
447
  <caption align="bottom">Figure: Two different views of the 2D simulated Point Spread Function.</caption>
448
  <tr>
449
      <td> <img WIDTH=300 src="images/sepia-small/whirl-starsima.png"> </td>
450
      <td> <img WIDTH=300 src="images/sepia-small/whirl-starsimb.png"> </td>
451
  </tr>
452
</table>
453
454
455
<pre>
456
   pdl&gt; use PDL::Graphics::TriD; # Load the 3D graphics module
457
   pdl&gt; imag3d [$fit2d];</pre>
458
<p>If you do this on your computer you should be able to look at the graphic from
459
different sides by simply dragging in the plot window with the mouse! You can
460
also zoom in and out with the right mouse button.  Note that <code>imag3d</code> has it's
461
a rather different syntax for processing it's arguments - for very good reasons
462
- we'll explore 3D graphics further in <a href="/PDL/Book/TriD.html">the PDL::Book::TriD manpage</a>.</p>
463
<pre>
464
   To continue: Select the TriD window and type q</pre>
465
<p>Finally here's something interesting. Let's take our fake star and place it
466
elsewhere on the galaxy image.</p>
467
<pre>
468
  pdl&gt; $newsection = $gal(50:70,70:90);
469
  pdl&gt; $newsection +=  $fit2d;
470
  pdl&gt; imag $gal,0,300;</pre>
471
<p>We have a bright new star where none existed before! The C-style <code>+=</code>
472
increment operator is worth noting - it actually modifies the contents of
473
<code>$newsection</code> in-place. And because <code>$newsection</code> is a <em>slice</em> of <code>$gal</code>
474
the change also affects <code>$gal</code>. This is an important property of slices - any
475
change to the slice affects the <em>parent</em>. This kind of parent/child
476
relationship is a powerful property of many PDL functions, not just slicing.
477
What's more in many cases it leads to memory efficiency, when this kind of
478
linear slice is stored we only store the start/stop/step and not a new copy of
479
the actual data.</p>
480
<p>Of course sometimes we DO want a new copy of the actual data, for example if we
481
plan to do something evil to it. To do this we could use the alternative form:</p>
482
<pre>
483
  pdl&gt; $newsection = $newsection +  $fit2d</pre>
484
<p>Now a new version of <code>$newsection</code>  is created which has nothing to 
485
do with the original <code>$gal</code>. In fact there is more than one way to do
486
this as we will see in later chapters.</p>
487
<p>Just to amuse ourselves, lets write a short script to cover M51 with dozens of fake
488
stars of random brightnesses:</p>
489
<pre>
490
   use PDL;
491
   use PDL::Graphics::Simple;
492
   use PDL::NiceSlice;  # must use in each program file</pre>
493
<pre>
494
   srand(42); # Set the random number seed
495
   $gal  = rfits &quot;fixed_gal.fits&quot;;
496
   $star = rfits &quot;fake_star.fits&quot;;
497
   
498
   sub addstar {
499
      ($x,$y) = @_;
500
      $xx = $x+20; $yy = $y+20;
501
      # Note use of slice on the LHS!
502
      $gal($x:$xx,$y:$yy) += $star * rand(2);
503
   }</pre>
504
<pre>
505
   for (1..100) {
506
      $x1 = int(rand(470)+10);
507
      $y1 = int(rand(470)+10);
508
      addstar($x1,$y1);
509
   }
510
   imag $gal,0,1000;</pre>
511
<p>This ought to give the casual reader some flavour of the Perl syntax - quite simple
512
and quite like C except that the entities being manipulated here are entire
513
arrays of data, not single numbers. The result is shown, for amusement,
514
in the figure below and takes virtually no time to compute.</p>
515
516
<table>
517
  <caption align="bottom">Figure: M51 covered in fake stars.</caption>
518
  <tr> <td> <img WIDTH=300 src="images/sepia-small/whirl-fakestars.png"> </td> </tr>
519
</table>
520
521
522
<p>
523
</p>
524
<h2><a name="getting_complex_with_m51">Getting Complex with M51</a></h2>
525
<p>To conclude this frantic whirl through the possibilities of PDL, let's look at
526
a moderately complex (sic) example. We'll take M51 and try to enhance it to reveal the
527
large-scale structure, and then subtract this to reveal small-scale structure.</p>
528
<p>Just to show off we'll use a method based on the Fourier transform - don't
529
worry if you don't know much about these, all you need to know is that the
530
Fourier transform turns the image into an 'inverse' image, with
531
complex numbers, where each pixel
532
represents the strength of wavelengths of different scales in the image. 
533
Let's do it:</p>
534
<pre>
535
  pdl&gt; use PDL::FFT; # Load Fast Fourier Transform package
536
  pdl&gt; $gal = rfits &quot;fixed_gal.fits&quot;;</pre>
537
<p>Now <code>$gal</code> contains real values, to do the Fourier transform it has to
538
have complex values. We create a variable <code>$imag</code> to hold the imaginary
539
component and set to zero.(For reasons of efficiency complex numbers
540
are represented in PDL by separate real and imaginary arrays - more about this
541
in Chapter 2.)</p>
542
<pre>
543
  pdl&gt; $imag = $gal * 0;       # Create imaginary component, equal to zero
544
  pdl&gt; fftnd $gal, $imag;      # Perform Fourier transform</pre>
545
<p><code>fftnd</code> performs a Fast Fourier Transform, in-place, on arbitrary-dimensioned data (i.e.
546
it is 'N-dimensional').  You can display <code>$gal</code>  after the FFT but you won't see
547
much. If at this point we ran <code>ifftnd</code> to invert it we would get the original
548
<code>$gal</code>  back.</p>
549
<p>If we want to enhance the large-scale structure we want to make a filter to only
550
let through low-frequencies:</p>
551
<pre>
552
  pdl&gt; $tmp = rvals($gal)&lt;10;        # Radially-symmetric filter function
553
  pdl&gt; use PDL::ImageND;             # provides kernctr()
554
  pdl&gt; $filter = kernctr $tmp, $tmp; # Shift origin to 0,0
555
  pdl&gt; imag $filter;</pre>
556
557
<table>
558
  <caption align="bottom"> Figure: The result of <code>kernctr()</code> </caption>
559
  <tr> <td> <img WIDTH=300 src="images/sepia-small/gal-filter.png"> </td> </tr>
560
</table>
561
562
<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>
563
564
<pre>
565
  pdl&gt; ($gal2, $imag2) = cmul $gal, $imag, $filter, 0;
566
  pdl&gt; ifftnd $gal2, $imag2;
567
  pdl&gt; imag $gal2,0,300;</pre>
568
569
570
<table>
571
  <caption align="bottom">Figure: Fourier filtered smoothed image
572
      and contrast enhanced image with the smoothed image subtracted. </caption>
573
  <tr>
574
      <td> <img WIDTH=300 src="images/sepia-small/whirl-ffta.png"> </td>
575
      <td> <img WIDTH=300 src="images/sepia-small/whirl-fftb.png"> </td>
576
  </tr>
577
</table>
578
579
580
<p>Well that looks quite a bit different!  Just about all the
581
high-frequency information has vanished. To see the high-frequency
582
information we can just subtract our filtered image from the original to
583
form the right hand image.</p>
584
585
<pre>
586
  pdl&gt; $orig = rfits &quot;fixed_gal.fits&quot;;
587
  pdl&gt; imag $orig-$gal2,0,100;</pre>
588
<p>
589
</p>
590
<h2><a name="roundoff">Roundoff</a></h2>
591
<p>Well that is probably enough abuse of Messier 51. We have demonstrated the ease
592
of simple and complex data processing with PDL and how PDL fits neatly in to
593
the Perl syntax as well as extending it.  You have come across basic
594
arithmetical operations and a scattering of useful functions - and learned how
595
to find more.  You certainly ought now to have a good feel of what PDL is all
596
about.  In the next chapter we'll take a more comprehensive look at the basic
597
parts of PDL that all keen PDL users should know.</p>
598
599
<h3><a name="pdlbook">Where to Read More</a></h3>
600
<p> We hope this excerpt from the
601
<a href="http://sourceforge.net/projects/pdl/files/PDL_2013/PDL-Book/">PDL Book</a>
602
has given you a taste of what PDL is capable of.  Please download the full text to
603
find out more.</p>
604
605
<p>As always your one-stop-shop for all things PDL is at
606
<a href="http://pdl.perl.org">http://pdl.perl.org</a>.</p>
607
608
</body>
609
610
</html>

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks