Diff of /PDL/Book/FirstSteps.pod [aa2d63] .. [845b10] Maximize Restore

  Switch to unified view

a/PDL/Book/FirstSteps.pod b/PDL/Book/FirstSteps.pod
...
...
85
similar to C and FORTRAN in how expressions are constructed.  (In fact much
85
similar to C and FORTRAN in how expressions are constructed.  (In fact much
86
more like C than FORTRAN). It is interesting to reflect on how much C code
86
more like C than FORTRAN). It is interesting to reflect on how much C code
87
would be required to generate the same display, even given the existence of
87
would be required to generate the same display, even given the existence of
88
some convenient graphics library.
88
some convenient graphics library.
89
89
90
=for html <img WIDTH=400 src="firststeps/whirl-sync.png">
90
=for html <img WIDTH=400 src="firststeps/sepia/crop/whirl-sync.png">
91
91
92
  Figure of a two dimensional C<sin>  function.
92
  Figure of a two dimensional C<sin>  function.
93
93
94
That's all fine. But what if we wanted to achieve the same results in a standalone
94
That's all fine. But what if we wanted to achieve the same results in a standalone
95
perl script? Well it is pretty simple:
95
perl script? Well it is pretty simple:
...
...
172
172
173
Let's look at M51:
173
Let's look at M51:
174
174
175
  pdl> imag $a;
175
  pdl> imag $a;
176
176
177
=for html <img WIDTH=600 src="firststeps/whirl-m51.png">
177
=for html <img WIDTH=600 src="firststeps/sepia/crop/whirl-m51-default.png">
178
178
179
  Figure of the raw image C<m51_raw.fits> shown with
179
  Figure of the raw image C<m51_raw.fits> shown with
180
  progressively greater contrast using the C<imag> command.
180
  progressively greater contrast using the C<imag> command.
181
181
182
A couple of bright spots can be seen, but where is the galaxy? It's the faint
182
A couple of bright spots can be seen, but where is the galaxy? It's the faint
...
...
184
the faintest to the brightest pixel, and only the bright star slightly to the
184
the faintest to the brightest pixel, and only the bright star slightly to the
185
bottom right of the center can be seen without contrast enhancement. We can
185
bottom right of the center can be seen without contrast enhancement. We can
186
easily change that by specifying the black/white data values (Note: C<#> starts
186
easily change that by specifying the black/white data values (Note: C<#> starts
187
a Perl comment and can be ignored - i.e. no need to type the stuff after it!):
187
a Perl comment and can be ignored - i.e. no need to type the stuff after it!):
188
188
189
=for html <img WIDTH=600 src="firststeps/sepia/crop/whirl-m51-1000.png">
190
189
  pdl> imag $a,0,1000; # More contrast
191
  pdl> imag $a,0,1000; # More contrast
192
193
194
=for html <img WIDTH=600 src="firststeps/sepia/crop/whirl-m51-300.png">
195
190
  pdl> imag $a,0,300;  # Even more contrast
196
  pdl> imag $a,0,300;  # Even more contrast
197
191
198
192
You can see that C<imag> takes additional arguments to specify the display
199
You can see that C<imag> takes additional arguments to specify the display
193
range. In fact C<imag> takes quite a few arguments, many of them optional. By
200
range. In fact C<imag> takes quite a few arguments, many of them optional. By
194
typing 'C<help imag>' at the C<pdl> prompt we can find out all about the
201
typing 'C<help imag>' at the C<pdl> prompt we can find out all about the
195
function.
202
function.
...
...
206
213
207
This is shown in the next figure. Because the image is of a uniform field,
214
This is shown in the next figure. Because the image is of a uniform field,
208
the actual image reflects the detector sensitivity. To correct our M51
215
the actual image reflects the detector sensitivity. To correct our M51
209
image, we merely have to divide the image by the flatfield:
216
image, we merely have to divide the image by the flatfield:
210
217
211
=for html <img WIDTH=400 src="firststeps/whirl-flat.png">
218
=for html <img WIDTH=400 src="firststeps/sepia/crop/whirl-flat.png">
212
219
213
  Figure: The 'flatfield' image showing the detector sensitivity of the raw data.
220
Figure: The 'flatfield' image showing the detector sensitivity of the raw data.
214
221
215
222
216
  pdl> $gal = $a / $flat;
223
  pdl> $gal = $a / $flat;
217
  pdl> imag $gal,0,300;  
224
  pdl> imag $gal,0,300;  
218
  pdl> wfits $gal, 'fixed_gal.fits'; # Save our work as a FITS file
225
  pdl> wfits $gal, 'fixed_gal.fits'; # Save our work as a FITS file
...
...
227
large chunks of data. If you did the same operation using normal
234
large chunks of data. If you did the same operation using normal
228
perl arrays rather than piddles it would be about ten to twenty times slower
235
perl arrays rather than piddles it would be about ten to twenty times slower
229
(and use ten times more memory). In fact we can
236
(and use ten times more memory). In fact we can
230
do whatever arithmetic operations we like on image piddles:
237
do whatever arithmetic operations we like on image piddles:
231
238
232
=for html <img WIDTH=400 src="firststeps/whirl-flattened.png">
239
=for html <img WIDTH=400 src="firststeps/sepia/crop/whirl-flattened.png">
233
240
234
  Figure: The M51 image corrected for the flatfield.
241
Figure: The M51 image corrected for the flatfield.
235
242
236
243
237
  pdl> $funny = log(($gal/300)**2 - $gal/100  + 4); 
244
  pdl> $funny = log(($gal/300)**2 - $gal/100  + 4); 
238
  pdl> imag $funny; # Surprise!
245
  pdl> imag $funny; # Surprise!
239
246
...
...
289
296
290
  pdl> $mask = $r<50;
297
  pdl> $mask = $r<50;
291
  pdl> imag $mask;
298
  pdl> imag $mask;
292
  ...
299
  ...
293
300
294
=for html <img WIDTH=400 src="firststeps/whirl-mask.png">
301
=for html <img WIDTH=400 src="firststeps/sepia/crop/whirl-maska.png">
295
302
296
 Figure: Using  C<rvals> to generate a mask image
303
=for html <img WIDTH=400 src="firststeps/sepia/crop/whirl-maskb.png">
297
 to isolate the galaxy bulge and disk.
304
298
 
305
=for html <img WIDTH=400 src="firststeps/sepia/crop/whirl-maskc.png">
299
 Top row: radial gradient image C<$r>, and radial gradient
306
300
 masked with less than operator C<$r < 50>.
307
=for html <img WIDTH=400 src="firststeps/sepia/crop/whirl-maskd.png">
301
 
308
309
Figure: Using  C<rvals> to generate a mask image to isolate the galaxy bulge and disk.
310
Top row: radial gradient image C<$r>, and radial gradient masked with less than operator C<$r < 50>.
302
 Bottom row: Bulge and disk of the galaxy.
311
Bottom row: Bulge and disk of the galaxy.
303
312
304
The Perl I<less than operator> is applied to all pixels in the image.
313
The Perl I<less than operator> is applied to all pixels in the image.
305
You can see the result is an image which is 0 on the outskirts and 1 in
314
You can see the result is an image which is 0 on the outskirts and 1 in
306
the area of the nucleus. We can then simply use the mask image to
315
the area of the nucleus. We can then simply use the mask image to
307
isolate in a simple way the bulge and disk components (lower row) and it
316
isolate in a simple way the bulge and disk components (lower row) and it
...
...
386
The blurring is usually measured
395
The blurring is usually measured
387
by the 'Full Width Half Maximum' (FWHM) - or in plain terms how
396
by the 'Full Width Half Maximum' (FWHM) - or in plain terms how
388
fat the profile is across when it drops by half. Looking at the plot
397
fat the profile is across when it drops by half. Looking at the plot
389
it looks like this is about 2-3 pixels - pretty compact!
398
it looks like this is about 2-3 pixels - pretty compact!
390
399
391
=for html <img WIDTH=400 src="firststeps/whirl-starradial.png">
400
=for html <img WIDTH=400 src="firststeps/sepia/crop/whirl-star-radial.png">
392
401
393
  Figure: Radial light profile of the bright star with fitted curve.
402
Figure: Radial light profile of the bright star with fitted curve.
394
403
395
Well we don't just want a guess - let's fit the profile with a function.
404
Well we don't just want a guess - let's fit the profile with a function.
396
These blurring functions are usually represented by the C<Gaussian>
405
These blurring functions are usually represented by the C<Gaussian>
397
function. PDL comes with a whole variety of general purpose and
406
function. PDL comes with a whole variety of general purpose and
398
special purpose fitting functions which people have written for
407
special purpose fitting functions which people have written for
...
...
479
than one graphics library (some see this as a flaw, some
488
than one graphics library (some see this as a flaw, some
480
as a feature!). Using the C<PDL::Graphics::TriD> library
489
as a feature!). Using the C<PDL::Graphics::TriD> library
481
which does OpenGL graphics we can look at our simulated
490
which does OpenGL graphics we can look at our simulated
482
star in 3D (see the right hand panel);
491
star in 3D (see the right hand panel);
483
492
493
=for html <img WIDTH=300 src="firststeps/sepia/crop/whirl-starsima.png">
494
484
=for html <img WIDTH=600 src="firststeps/whirl-starsim.png">
495
=for html <img WIDTH=300 src="firststeps/grey/work/whirl-starsimb.png">
485
496
486
  Figure: Two different views of the 2D simulated Point Spread Function.
497
Figure: Two different views of the 2D simulated Point Spread Function.
487
498
488
   pdl> use PDL::Graphics::TriD; # Load the 3D graphics module
499
   pdl> use PDL::Graphics::TriD; # Load the 3D graphics module
489
   pdl> imag3d [$fit2d];
500
   pdl> imag3d [$fit2d];
490
501
491
If you do this on your computer you should be able to look at the graphic from
502
If you do this on your computer you should be able to look at the graphic from
...
...
550
This ought to give the casual reader some flavour of the Perl syntax - quite simple
561
This ought to give the casual reader some flavour of the Perl syntax - quite simple
551
and quite like C except that the entities being manipulated here are entire
562
and quite like C except that the entities being manipulated here are entire
552
arrays of data, not single numbers. The result is shown, for amusement,
563
arrays of data, not single numbers. The result is shown, for amusement,
553
in the figure below and takes virtually no time to compute.
564
in the figure below and takes virtually no time to compute.
554
565
555
=for html <img WIDTH=400 src="firststeps/whirl-fakestars.png">
566
=for html <img WIDTH=400 src="firststeps/sepia/crop/whirl-fakestars.png">
556
567
557
  Figure: M51 covered in fake stars.
568
Figure: M51 covered in fake stars.
558
569
559
=head2 Getting Complex with M51
570
=head2 Getting Complex with M51
560
571
561
To conclude this frantic whirl through the possibilities of PDL, let's look at
572
To conclude this frantic whirl through the possibilities of PDL, let's look at
562
a moderately complex (sic) example. We'll take M51 and try to enhance it to reveal the
573
a moderately complex (sic) example. We'll take M51 and try to enhance it to reveal the
...
...
592
  pdl> $tmp = rvals($gal)<10;        # Radially-symmetric filter function
603
  pdl> $tmp = rvals($gal)<10;        # Radially-symmetric filter function
593
  pdl> use PDL::ImageND;             # provides kernctr()
604
  pdl> use PDL::ImageND;             # provides kernctr()
594
  pdl> $filter = kernctr $tmp, $tmp; # Shift origin to 0,0
605
  pdl> $filter = kernctr $tmp, $tmp; # Shift origin to 0,0
595
  pdl> imag $filter;
606
  pdl> imag $filter;
596
607
597
=for html <img WIDTH=400 src="firststeps/gal-filter.png">
608
=for html <img WIDTH=400 src="firststeps/sepia/crop/gal-filter.png">
598
609
599
You can see from the image that C<$filter> is zero everywhere except near the origin
610
You can see from the image that C<$filter> is zero everywhere except near the origin
600
(0,0) (and the 3 reflected corners). As a result it only lets through
611
(0,0) (and the 3 reflected corners). As a result it only lets through
601
low-frequency wavelengths. So we multiply by the filter and  FFT back to
612
low-frequency wavelengths. So we multiply by the filter and  FFT back to
602
see the result (C<cmul> is complex multiplication):
613
see the result (C<cmul> is complex multiplication):
603
614
604
  pdl> ($gal2, $imag2) = cmul $gal, $imag, $filter, 0;
615
  pdl> ($gal2, $imag2) = cmul $gal, $imag, $filter, 0;
605
  pdl> ifftnd $gal2, $imag2;
616
  pdl> ifftnd $gal2, $imag2;
606
  pdl> imag $gal2,0,300;
617
  pdl> imag $gal2,0,300;
607
618
608
=for html <img WIDTH=400 src="firststeps/whirl-fft.png">
619
=for html <img WIDTH=400 src="firststeps/sepia/crop/whirl-ffta.png">
609
620
610
  Figure: Fourier filtered smoothed image and contrast enhanced
621
=for html <img WIDTH=400 src="firststeps/sepia/crop/whirl-fftb.png">
611
  image with the smoothed image subtracted.
622
623
Figure: Fourier filtered smoothed image and contrast enhanced image with the smoothed image subtracted.
612
624
613
Well that looks quite a bit different!  Just about all the
625
Well that looks quite a bit different!  Just about all the
614
high-frequency information has vanished. To see the high-frequency
626
high-frequency information has vanished. To see the high-frequency
615
information we can just subtract our filtered image from the original to
627
information we can just subtract our filtered image from the original to
616
form the right hand image.
628
form the right hand image.