File Coverage

blib/lib/PDL/IO/FITS.pm
Criterion Covered Total %
statement 627 966 64.9
branch 288 578 49.8
condition 104 239 43.5
subroutine 37 52 71.1
pod 3 10 30.0
total 1059 1845 57.4


line stmt bran cond sub pod time code
1             =head1 NAME
2              
3             PDL::IO::FITS -- Simple FITS support for PDL
4              
5             =head1 SYNOPSIS
6              
7             use PDL;
8             use PDL::IO::FITS;
9              
10             $x = rfits('foo.fits'); # read a FITS file
11             $x->wfits('bar.fits'); # write a FITS file
12              
13             =head1 DESCRIPTION
14              
15             This module provides basic FITS support for PDL, in the sense of
16             reading and writing whole FITS files. (For more complex operations,
17             such as prefiltering rows out of tables or performing operations on
18             the FITS file in-place on disk), you can use the Astro::FITS::CFITSIO
19             module that is available on CPAN.
20              
21             Basic FITS image files are supported, along with BINTABLE and IMAGE extensions.
22             ASCII Table support is planned, as are the HEASARC bintable extensions that
23             are recommended in the 1999 FITS standard.
24              
25             Table support is based on hashes and named columns, rather than the
26             less convenient (but slightly more congruent) technique of perl lists
27             of numbered columns.
28              
29             The principle interface routines are C and C, for
30             reading and writing respectively. FITS headers are returned as perl
31             hashes or (if the module is present) Astro::FITS::Header objects that
32             are tied to perl hashes. Astro::FITS::Header objects provide
33             convenient access through the tied hash interface, but also allow you
34             to control the card structure in more detail using a separate method
35             interface; see the L
36             documentation for details.
37              
38             =head1 AUTHOR
39              
40             Copyright (C) Karl Glazebrook, Craig DeForest, and Doug Burke, 1997-2010.
41             There is no warranty. You are allowed to redistribute and/or modify
42             this software under certain conditions. For details, see the file
43             COPYING in the PDL distribution. If this file is separated from the
44             PDL distribution, the copyright notice should be pasted into in this file.
45              
46             =head1 FUNCTIONS
47              
48             =cut
49              
50 49     49   77436 use strict;
  49         116  
  49         4765  
51              
52             BEGIN {
53              
54             package PDL::IO::FITS;
55              
56 49     49   236 $PDL::IO::FITS::VERSION = 0.92; # Will be 1.0 when ascii table read/write works.
57              
58 49         184 our @EXPORT_OK = qw( rfits rfitshdr wfits );
59 49         216 our %EXPORT_TAGS = (Func=>[@EXPORT_OK]);
60 49         978 our @ISA = ('PDL::Exporter');
61              
62 49     49   332 use PDL::Core;
  49         115  
  49         338  
63 49     49   19905 use PDL::Config;
  49         2885  
  49         1965  
64 49     49   958 use PDL::IO::Misc;
  49         239  
  49         358  
65 49     49   338 use PDL::Exporter;
  49         105  
  49         281  
66 49     49   257 use PDL::Primitive;
  49         103  
  49         283  
67 49     49   386 use PDL::Types;
  49         100  
  49         6735  
68 49     49   347 use PDL::Options;
  49         109  
  49         3428  
69 49     49   328 use PDL::Bad;
  49         101  
  49         356  
70             # use PDL::NiceSlice;
71 49     49   439 use Carp;
  49         100  
  49         2831  
72 49     49   375 use strict;
  49         126  
  49         9672  
73              
74             ##############################
75             #
76             # Check if there's Astro::FITS::Header support, and set flag.
77             # Kludgy but it only has to run once, on first load. --CED
78             #
79 49     49   4098 eval "use Astro::FITS::Header;";
  49         33389  
  49         487670  
  49         1108  
80 49         207 $PDL::Astro_FITS_Header = (defined $Astro::FITS::Header::VERSION);
81 49 50       230 if($PDL::Astro_FITS_Header) {
82 49         130 my($x) = $Astro::FITS::Header::VERSION;
83 49         453 $x =~ s/[^0-9\.].*//;
84 49 50       310 $PDL::Astro_FITS_Header = 0 if($x < 1.12);
85             }
86              
87 49 50       44108 unless($PDL::Astro_FITS_Header) {
88 0 0 0     0 unless($ENV{"PDL_FITS_LEGACY"} || $PDL::Config{FITS_LEGACY}) {
89 0         0 print(STDERR "\n\nWARNING: Can't find the Astro::FITS::Header module, limiting FITS support.\n\n PDL will use the deprecated legacy perl hash handling code but will not\n properly support tables, FITS extensions, or COMMENT cards. You really\n ought to install the Astro::FITS::Header module, available from\n 'http://www.cpan.org'. (You can also get rid of this message by setting\n the environment variable 'PDL_FITS_LEGACY' or the global PDL config value (in perldl.conf)\n \$PDL::Config{FITS_LEGACY} to 1.\n\n");
90             }
91             }
92             }
93              
94             package PDL::IO::FITS;
95              
96             ## declare subroutines
97              
98             sub _wfits_nullhdu ($);
99             sub _wfits_table ($$$);
100              
101             =head2 rfits()
102              
103             =for ref
104              
105             Simple piddle FITS reader.
106              
107             =for example
108              
109             $pdl = rfits('file.fits'); # Read a simple FITS image
110              
111             Suffix magic:
112              
113             $pdl = rfits('file.fits.gz'); # Read a file with gunzip(1)
114             $pdl = rfits('file.fits.Z'); # Read a file with uncompress(1)
115              
116             $pdl = rfits('file.fits[2]'); # Read 2nd extension
117             $pdl = rfits('file.fits.gz[3]'); # Read 3rd extension
118             @pdls = rfits('file.fits'); # Read primary data and extensions
119              
120             Tilde expansion:
121              
122             #expand leading ~ to home directory (using glob())
123             $pdl = rfits '~/filename.fits';
124              
125             $hdr = rfits('file.fits',{data=>0}); # Options hash changes behavior
126              
127             In list context, C reads the primary image and all possible
128             extensions, returning them in the same order that they occurred in the
129             file -- except that, by default, the primary HDU is skipped if it
130             contains no data. In scalar context, the default is to read the first
131             HDU that contains data. One can read other HDU's by using the [n]
132             syntax. Using the [0] syntax forces a read of the first HDU,
133             regardless of whether it contains data or no. Currently recognized
134             extensions are IMAGE and BINTABLE. (See the addendum on EXTENSIONS
135             for details).
136              
137             C accepts several options that may be passed in as a hash ref
138             if desired:
139              
140             =over 3
141              
142             =item bscale (default=1)
143              
144             Determines whether the data are linearly scaled using the BSCALE/BZERO keywords
145             in the FITS header. To read in the exact data values in the file, set this
146             to 0.
147              
148             =item data (default=1)
149              
150             Determines whether to read the data, or just the header. If you set this to
151             0, you will get back the FITS header rather than the data themselves. (Note
152             that the header is normally returned as the C field of the returned PDL;
153             this causes it to be returned as a hash ref directly.)
154              
155             =item hdrcpy (default=0)
156              
157             Determines whether the L flag is set in the returned
158             PDL. Setting the flag will cause an explicit deep copy of the header whenever
159             you use the returned PDL in an arithmetic or slicing operation. That is useful
160             in many circumstances but also causes a hit in speed. When two or more PDLs
161             with hdrcpy set are used in an expression, the result gets the header of the
162             first PDL in the expression. See L for an example.
163              
164             =item expand (default=1)
165              
166             Determines whether auto-expansion of tile-compressed images should happen.
167             Tile-compressed images are transmitted as binary tables with particular
168             fields ("ZIMAGE") set. Leaving this alone does what you want most of the
169             time, unpacking such images transparently and returning the data and header
170             as if they were part of a normal IMAGE extension. Setting "expand" to 0
171             delivers the binary table, rather than unpacking it into an image.
172              
173             =item afh (default=1)
174              
175             By default rfits uses Astro::FITS::Header tied-hash objects to contain
176             the FITS header information. This permits explicit control over FITS
177             card information, and conforms well with the FITS specification. But
178             Astro::FITS::Header objects are about 40-60x more memory intensive
179             than comparable perl hashes, and also use ~10x more CPU to manage.
180             For jobs where header processing performance is important (e.g. reading
181             just the headers of 1,000 FITS files), set afh to 0 to use the legacy parser
182             and get a large boost in speed.
183              
184             =back
185              
186             FITS image headers are stored in the output PDL and can be retrieved
187             with L or L. The
188             L flag of the PDL is set so that the header
189             is copied to derived piddles by default. (This is inefficient if you
190             are planning to do lots of small operations on the data; clear
191             the flag with "->hcpy(0)" or via the options hash if that's the case.)
192              
193             The header is a hash whose keys are the keywords in the FITS header.
194             If you have the "Astro::FITS::Header" module installed, the header is
195             actually a tied hash to a FITS header object, which can give you
196             more control over card order, comment fields, and variable types.
197             (see L for details).
198              
199             The header keywords are converted to I per the FITS
200             standard. Access is case-insensitive on the perl side, provided that
201             Astro::FITS::Header is installed.
202              
203             If Astro::FITS::Header is not installed, then a built-in legacy parser
204             is used to generate the header hash. Keyword-associated comments in
205             the headers are stored under the hash key
206             C<< _COMMENT> >>. All HISTORY cards in the header are
207             collected into a single multiline string stored in the C key.
208             All COMMENT cards are similarly collected under the C key.
209              
210             =head3 BSCALE/BZERO
211              
212             If the BSCALE and/or BZERO keywords are set, they are applied to the
213             image before it is returned. The returned PDL is promoted as
214             necessary to contain the multiplied values, and the BSCALE and BZERO
215             keywords are deleted from the header for clarity. If you don't want
216             this type of processing, set 'bscale=>0' in the options hash.
217              
218             =head3 EXTENSIONS
219              
220             Sometimes a FITS file contains only extensions and a stub header in
221             the first header/data unit ("primary HDU"). In scalar context, you
222             normally only get back the primary HDU -- but in this special case,
223             you get back the first extension HDU. You can force a read of the
224             primary HDU by adding a '[0]' suffix to the file name.
225              
226             =head3 BINTABLE EXTENSIONS
227              
228             Binary tables are handled. Currently only the following PDL
229             datatypes are supported: byte, short, ushort, long, float, and
230             double. At present ushort() data is written as a long rather than
231             as a short with TSCAL/ZERO; this may change.
232              
233             The return value for a binary table is a hash ref containing the names
234             of the columns in the table (in UPPER CASE as per the FITS standard).
235             Each element of the hash contains a PDL (for numerical values) or a
236             perl list (for string values). The PDL's 0th dimension runs across
237             rows; the 1st dimension runs across the repeat index within the row
238             (for rows with more than one value). (Note that this is different from
239             standard threading order - but it allows Least Surprise to work when
240             adding more complicated objects such as collections of numbers (via
241             the repeat count) or variable length arrays.)
242              
243             Thus, if your table contains a column named C with type C<5D>,
244             the expression
245              
246             $x->{FOO}->((2))
247              
248             returns a 5-element double-precision PDL containing the values of FOO
249             from the third row of the table.
250              
251             The header of the table itself is parsed as with a normal FITS HDU,
252             and is returned in the element 'hdr' of the returned hash. You can
253             use that to preserve the original column order or access the table at a low
254             level, if you like.
255              
256             Scaling and zero-point adjustment are performed as with BSCALE/BZERO:
257             the appropriate keywords are deleted from the as-returned header. To avoid
258             this behavior, set 'bscale=>0' in the options hash.
259              
260             As appropriate, TSCAL/ZERO and TUNIT are copied into each column-PDL's
261             header as BSCALE/BZERO and BUNIT.
262              
263             The main hash also contains the element 'tbl', which is set
264             to 'binary' to distinguish it from an ASCII table.
265              
266             Because different columns in the table might have identical names in a
267             FITS file, the binary table reader practices collision avoidance. If
268             you have multiple columns named "FOO", then the first one encountered
269             (numerically) gets the name "FOO", the next one gets "FOO_1", and the
270             next "FOO_2", etc. The appropriate TTYPEn fields in the header are
271             changed to match the renamed column fields.
272              
273             Columns with no name are assigned the name "COL_", where starts
274             at 1 and increments for each no-name column found.
275              
276             Variable-length arrays are supported for reading. They are unpacked
277             into PDLs that appear exactly the same as the output for fixed-length
278             rows, except that each row is padded to the maximum length given
279             in the extra characters -- e.g. a row with TFORM of 1PB(300) will
280             yield an NAXIS2x300 output field in the final hash. The padding
281             uses the TNULL keyword for the column, or 0 if TNULL is not
282             present. The output hash also gets an additional field, "len_",
283             that contains the number of elements in each table row.
284              
285             =head3 TILE-COMPRESSED IMAGES
286              
287             CFITSIO and several large projects (including NASA's Solar Dynamics
288             Observatory) now support an unofficial extension to FITS that stores
289             images as a collection of individually compressed tiles within a
290             BINTABLE extension. These images are automagically uncompressed by
291             default, and delivered as if they were normal image files. You can
292             override this behavior by supplying the "expand" key in the options hash.
293              
294             Currently, only Rice compression is supported, though there is a framework
295             in place for adding other compression schemes.
296              
297             =for bad
298              
299             =head3 BAD VALUE HANDLING
300              
301             If a FITS file contains the C keyword (and has C 0>),
302             the piddle will have its bad flag set, and those elements which equal the
303             C value will be set bad. For C 0>, any NaN's are
304             converted to bad (if necessary).
305              
306              
307             =head2 rfitshdr()
308              
309             =for ref
310              
311             Read only the header of a FITS file or an extension within it.
312              
313             This is syntactic sugar for the C0> option to L.
314              
315             See L for details on header handling. rfitshdr() runs
316             the same code to read the header, but returns it rather than
317             reading in a data structure as well.
318              
319             =cut
320              
321             our $rfits_options = new PDL::Options( { bscale=>1, data=>1, hdrcpy=>0, expand=>1, afh=>1 } );
322              
323             sub PDL::rfitshdr {
324 0     0 0 0 my $class = shift;
325 0         0 my $file = shift;
326 0         0 my $u_opt = ifhref(shift);
327 0         0 $u_opt->{data} = 0;
328 0         0 PDL::rfits($class,$file,$u_opt);
329             }
330              
331             sub PDL::rfits {
332            
333 46     46 0 274 my $class = shift;
334            
335 46 50 33     365 barf 'Usage: $x = rfits($file) -or- $x = PDL->rfits($file)' if (@_ < 1 || @_ > 2);
336            
337 46         114 my $file = shift;
338            
339 46         297 my $u_opt = ifhref(shift);
340 46         299 my $opt = $rfits_options->options($u_opt);
341            
342 46         139 my($nbytes, $line, $name, $rest, $size, $i, $bscale, $bzero, $extnum);
343              
344 46         102 $nbytes = 0;
345              
346             # Modification 02/04/2005 - JB. Earlier version stripped the extension
347             # indicator which cancelled the check for empty primary data array at the end.
348 46 100       234 my $explicit_extension = ($file =~ m/\[\d+\]$/ ? 1 : 0);
349 46 100       216 $extnum = ( ($file =~ s/\[(\d+)\]$//) ? $1 : 0 );
350              
351 46         115 $file =~ s/^(~)/glob($1)/e; #tilde expansion.
  0         0  
352              
353 46 50       191 $file = "gunzip -c $file |" if $file =~ /\.gz$/; # Handle compression
354 46 50       195 $file = "uncompress -c $file |" if $file =~ /\.Z$/;
355            
356 46 50       348 my $fh = IO::File->new( $file )
357             or barf "FITS file $file not found";
358 46         5381 binmode $fh;
359              
360 46         113 my @extensions; # This accumulates the list in list context...
361 46         80 my $currentext=0;
362 46         132 my $pdl;
363              
364 46   100     126 hdu:{do { # Runs over extensions, in list context
  46         90  
365 48         123 my $ext_type = 'IMAGE'; # Gets the type of XTENSION if one is detected.
366 48         166 my $foo={}; # To go in pdl
367 48         126 my @history=();
368 48         114 my @cards = ();
369            
370 48         417 $pdl = $class->new;
371            
372              
373             # If $opt->{data} is false, then the reading routines leave the
374             # file alone, so the file pointer is left at the end of the last
375             # header. Skip over the unread data to the next extension...
376            
377 48 50 66     302 if( wantarray and !$opt->{data} and @extensions) {
      33        
378 0   0     0 while( $fh->read($line,80) && ($line !~ /^XTENSION=/) && !$fh->eof() ) {
      0        
379 0         0 $fh->read($line,2880-80);
380             };
381            
382             return @extensions
383 0 0       0 if($fh->eof());
384              
385             } else {
386 48         307 my $ct = $fh->read($line,80);
387 48 50 66     2380 barf "file $file is not in FITS-format:\n$line\n"
388             if( $nbytes==0 && ($line !~ /^SIMPLE = +T/));
389 48 50 33     278 last hdu if($fh->eof() || !$ct);
390             }
391              
392 48         568 $nbytes = 80; # Number of bytes read from this extension (1 line so far)
393            
394 48 100       244 if($line =~ /^XTENSION= \'(\w+)\s*\'/) {
    50          
395 2         10 $ext_type = $1;
396             } elsif( @extensions ) {
397              
398 0 0       0 print "Warning: expected XTENSION, found '$line'. Exiting.\n"
399             if($PDL::verbose);
400 0         0 last hdu;
401             }
402              
403 48 50       206 push(@cards,$line) if($PDL::Astro_FITS_Header);
404            
405             #
406             # If we are in scalar context, skip to the desired extension
407             # number. [This implementation is really slow since we have
408             # to read the whole file. Someone Really Ought To rework it to
409             # read individual headers and skip forward an extension at a
410             # a time with seek() calls. ]
411             # --CD
412             #
413            
414 48 100 100     296 if(!wantarray and $currentext != $extnum) {
415              
416 3         7 skipper: while(1) {
417             # Move to next record
418 3         12 $nbytes += $fh->read($line,2880-80);
419 3 50       37 barf "Unexpected end of FITS file\n" if $fh->eof();
420             # Read start of next record
421 3         24 $nbytes += $fh->read($line,80);
422 3 50       20 barf "Unexpected end of FITS file\n" if $fh->eof();
423             # Check if we have found the new extension
424             # if not move on
425              
426 3 50       118 $currentext++ if $line =~ /^XTENSION\= \'(\w+)\s*\'/;
427 3 50       10 if ($currentext == $extnum) {
428 3         11 $ext_type = $1;
429 3         11 last skipper;
430             }
431             }
432             } # End of skipping to desired extension
433            
434             #
435             # Snarf up the found header, and parse it if Astro::FITS::Header
436             # does not exist.
437             #
438              
439 48 50 33     264 if($PDL::Astro_FITS_Header and $opt->{afh}) {
440              
441             ## Astro::FITS::Header parsing. Snarf lines to the END card,
442             ## and pass them to Astro::FITS::Header.
443              
444 48   66     104 do {
445 437         2566 $nbytes += $fh->read($line, 80);
446 437         2535 push(@cards,$line);
447             } while(!$fh->eof() && $line !~ m/^END(\s|\000)/);
448              
449 48         533 $nbytes += $fh->read(my $dummy, 2879 - ($nbytes-1)%2880);
450              
451 48         841 my($hdr) = Astro::FITS::Header->new(Cards => \@cards);
452 48         69038 my(%hdrhash);
453 48         333 tie %hdrhash,"Astro::FITS::Header",$hdr;
454 48         704 $foo = \%hdrhash;
455            
456             } else {
457            
458             ## Legacy (straight header-to-hash-ref) parsing.
459             ## Cheesy but fast.
460            
461 0         0 hdr_legacy: { do {
  0         0  
462 49     49   544 no strict 'refs';
  49         132  
  49         360351  
463             # skip if the first eight characters are ' '
464             # - as seen in headers from the DSS at STScI
465 0 0       0 if (substr($line,0,8) ne " " x 8) { # If non-blank
466            
467 0         0 $name = (split(' ',substr($line,0,8)))[0];
468              
469 0         0 $rest = substr($line,8);
470            
471 0 0       0 if ($name =~ m/^HISTORY/) {
472 0         0 push @history, $rest;
473             } else {
474 0         0 $$foo{$name} = "";
475            
476 0 0       0 $$foo{$name}=$1 if $rest =~ m|^= +([^\/\' ][^\/ ]*) *( +/(.*))?$| ;
477 0 0       0 $$foo{$name}=$1 if $rest =~ m|^= \'(.*)\' *( +/(.*))?$| ;
478 0 0       0 $$foo{COMMENT}{$name} = $3 if defined($3);
479             }
480             } # non-blank
481 0 0 0     0 last hdr_legacy if ((defined $name) && $name eq "END");
482 0         0 $nbytes += $fh->read($line, 80);
483             } while(!$fh->eof()); }
484              
485             # Clean up HISTORY card
486 0 0       0 $$foo{"HISTORY"} = \@history if $#history >= 0;
487            
488             # Step to end of header block in file
489 0         0 my $skip = 2879 - ($nbytes-1)%2880;
490 0 0       0 $fh->read(my $dummy, $skip) if $skip;
491 0         0 $nbytes += $skip;
492              
493             } # End of legacy header parsing
494            
495              
496            
497             ##############################
498             # Special case: if the file only contains
499             # extensions then read the first extension in scalar context,
500             # instead of the null zeroth extension.
501             #
502 48 100 100     314 if( !(defined $foo->{XTENSION}) # Primary header
      100        
      66        
503             and $foo->{NAXIS} == 0 # No data
504             and !wantarray # Scalar context
505             and !$explicit_extension # No HDU specifier
506             ) {
507 3 50       545 print "rfits: Skipping null primary HDU (use [0] to force read of primary)...\n"
508             if($PDL::verbose);
509 3         31 return PDL::rfits($class,$file.'[1]',$opt);
510             }
511            
512              
513             ##########
514             # If we don't want data, return the header from the HDU. Likewise,
515             # if NAXIS is 0 then there are no data, so return the header instead.
516 45 100 66     7486 if( ! $opt->{data} || $foo->{NAXIS}==0 ) {
517             # If we're NOT reading data, then return the header instead of the
518             # image.
519              
520 1         77 $pdl = $foo;
521              
522             } else {
523            
524             ##########
525             # Switch based on extension type to do the dirty work of reading
526             # the data. Handlers are listed in the _Extension patch-panel.
527            
528 44 50       3071 if (ref $PDL::IO::FITS::_Extension->{$ext_type} ) {
529            
530             # Pass $pdl into the extension reader for easier use -- but
531             # it just gets overwritten (and disappears) if ignored.
532            
533 44         130 $pdl = &{$PDL::IO::FITS::_Extension->{$ext_type}}($fh,$foo,$opt,$pdl);
  44         172  
534            
535             } else {
536              
537 0 0 0     0 print STDERR "rfits: Ignoring unknown extension '$ext_type'...\n"
538             if($PDL::verbose || $PDL::debug);
539            
540 0         0 $pdl = undef;
541              
542             }
543             }
544              
545             #
546             # Note -- $pdl isn't necessarily a PDL. It's only a $pdl if
547             # the extension was an IMAGE.
548             #
549 45 100       235 push(@extensions,$pdl) if(wantarray);
550 45         359 $currentext++;
551            
552             } while( wantarray && !$fh->eof() );} # Repeat if we are in list context
553            
554 43         348 $fh->close;
555            
556 43 100       1367 if(wantarray) {
557             ## By default, ditch primary HDU placeholder
558 1 50 33     15 if( ref($extensions[0]) eq 'HASH' and
      33        
      33        
559             $extensions[0]->{SIMPLE} and
560             exists($extensions[0]->{NAXIS}) and
561             $extensions[0]->{NAXIS} == 0
562             ) {
563 1         193 shift @extensions;
564             }
565             # Return all the extensions
566 1         40 return @extensions;
567             }
568 42         1104 return $pdl;
569             }
570              
571              
572 23     23 1 2465 sub rfits { PDL->rfits(@_); }
573              
574             sub rfitshdr {
575 0     0 1 0 my($file,$opt) = shift;
576 0         0 $opt->{data} =0;
577 0         0 PDL->rfitshdr($file,$opt);
578             }
579              
580              
581             ##############################
582             #
583             # FITS extensions patch-table links extension name to the supported reader.
584             # IMAGE extensions are a special case that gets read just like a normal
585             # FITS file.
586             #
587              
588             $PDL::IO::FITS::_Extension = {
589             IMAGE => \&_rfits_image
590             , BINTABLE => \&_rfits_bintable
591             };
592              
593              
594            
595             ##############################
596             #
597             # IMAGE extension -- this is also the default reader.
598             our $type_table = {
599             8=>$PDL_B,
600             16=>$PDL_S,
601             32=>$PDL_L,
602             64=>$PDL_LL,
603             -32=>$PDL_F,
604             -64=>$PDL_D
605             };
606              
607             our $type_table_2 = {
608             8=>byte,
609             16=>short,
610             32=>long,
611             64=>longlong,
612             -32=>float,
613             -64=>double
614             };
615              
616             sub _rfits_image($$$$) {
617 41 50   41   151 print "Reading IMAGE data...\n" if($PDL::verbose);
618 41         206 my $fh = shift; # file handle to read from
619 41         91 my $foo = shift; # $foo contains the pre-read header
620 41         76 my $opt = shift; # $opt contains the option hash
621 41         79 my $pdl = shift; # $pdl contains a pre-blessed virgin PDL
622              
623             # Setup piddle structure
624            
625 41 50       212 if( defined($type_table->{0 + $foo->{"BITPIX"}}) ) {
626 41         2987 $pdl->set_datatype( $type_table->{$foo->{"BITPIX"}} );
627             } else {
628 0         0 die("rfits: strange BITPIX value ".$foo->{"BITPIX"}." in header - I give up!\n");
629             }
630            
631 41         2938 my @dims; # Store the dimenions 1..N, compute total number of pixels
632 41         91 my $i = 1;
633 41         84 my $size = 1;
634              
635             ##second part of the conditional guards against a poorly-written hdr.
636 41   66     238 while(defined( $$foo{"NAXIS$i"} ) && $i <= $$foo{"NAXIS"}) {
637 68         9036 $size *= $$foo{"NAXIS$i"};
638 68         4472 push @dims, $$foo{"NAXIS$i"} ; $i++;
  68         4413  
639             }
640 41         2260 $pdl->setdims([@dims]);
641            
642 41         345 my $dref = $pdl->get_dataref();
643            
644 41 50       153 print "BITPIX = ",$$foo{"BITPIX"}," size = $size pixels \n"
645             if $PDL::verbose;
646            
647             # Slurp the FITS binary data
648            
649 41 50       118 print "Reading ",$size*PDL::Core::howbig($pdl->get_datatype) , " bytes\n"
650             if $PDL::verbose;
651            
652             # Read the data and pad to the next HDU
653 41         354 my $rdct = $size * PDL::Core::howbig($pdl->get_datatype);
654 41         245 $fh->read( $$dref, $rdct );
655 41         3036 $fh->read( my $dummy, 2880 - (($rdct-1) % 2880) - 1 );
656 41         498 $pdl->upd_data();
657              
658 41 50       328 if (!isbigendian() ) {
659             # Need to byte swap on little endian machines
660 41 100       278 bswap2($pdl) if $pdl->get_datatype == $PDL_S;
661 41 100 100     6558 bswap4($pdl) if $pdl->get_datatype == $PDL_L || $pdl->get_datatype == $PDL_F;
662 41 100 100     719 bswap8($pdl) if $pdl->get_datatype == $PDL_D || $pdl->get_datatype==$PDL_LL;
663             }
664            
665 41 50       170 if(exists $opt->{bscale}) {
666 41         150 $pdl = treat_bscale($pdl, $foo);
667             }
668            
669             # Header
670            
671 41         289 $pdl->sethdr($foo);
672              
673 41         199 $pdl->hdrcpy($opt->{hdrcpy});
674              
675 41         153 return $pdl;
676             }
677              
678             sub treat_bscale($$){
679 41     41 0 86 my $pdl = shift;
680 41         79 my $foo = shift;
681              
682 41 50       172 print "treating bscale...\n" if($PDL::debug);
683              
684 41 50       129 if ( $PDL::Bad::Status ) {
685             # do we have bad values? - needs to be done before BSCALE/BZERO
686             # (at least for integers)
687             #
688 41 100 100     260 if ( $$foo{BITPIX} > 0 and exists $$foo{BLANK} ) {
    100          
689             # integer, so bad value == BLANK keyword
690 1         135 my $blank = $foo->{BLANK};
691             # do we have to do any conversion?
692 1 50       69 if ( $blank == $pdl->badvalue() ) {
693 1         8 $pdl->badflag(1);
694             } else {
695             # we change all BLANK values to the current bad value
696             # (would not be needed with a per-piddle bad value)
697 0         0 $pdl->inplace->setvaltobad( $blank );
698             }
699             } elsif ( $foo->{BITPIX} < 0 ) {
700             # bad values are stored as NaN's in FITS
701             # let setnanbad decide if we need to change anything
702 23         3403 $pdl->inplace->setnantobad();
703             }
704 41 50 66     3365 print "FITS file may contain bad values.\n"
705             if $pdl->badflag() and $PDL::verbose;
706             } # if: PDL::Bad::Status
707            
708 41         108 my ($bscale, $bzero);
709 41         202 $bscale = $$foo{"BSCALE"}; $bzero = $$foo{"BZERO"};
  41         2162  
710 41 50       1987 print "BSCALE = $bscale && BZERO = $bzero\n" if $PDL::verbose;
711 41 50 33     180 $bscale = 1 if (!defined($bscale) || $bscale eq "");
712 41 50 33     227 $bzero = 0 if (!defined($bzero) || $bzero eq "");
713            
714             # Be clever and work out the final datatype before eating
715             # memory
716             #
717             # ensure we pick an element that is not equal to the bad value
718             # (is this OTT?)
719 41         97 my $tmp;
720              
721 41 100       227 if ( $pdl->badflag() == 0 ) {
    50          
722 39         197 $tmp = $pdl->flat()->slice("0:0");
723             } elsif ( $pdl->ngood > 0 ) {
724 2         9 my $index = which( $pdl->flat()->isbad() == 0 )->at(0);
725 2         22 $tmp = $pdl->flat()->slice("${index}:${index}");
726             } else {
727             # all bad, so ignore the type conversion and return
728             # -- too lazy to include this check in the code below,
729             # so just copy the header clean up stuff
730 0 0       0 print "All elements are bad.\n" if $PDL::verbose;
731            
732 0         0 delete $$foo{"BSCALE"}; delete $$foo{"BZERO"};
  0         0  
733 0         0 $tmp = $pdl;
734             } #end of BSCALE section (whew!)
735            
736            
737 41 50       258 $tmp = $tmp*$bscale if $bscale != 1; # Dummy run on one element
738 41 50       131 $tmp = $tmp+$bzero if $bzero != 0;
739            
740 41 50       204 $pdl = $pdl->convert($tmp->type) if $tmp->get_datatype != $pdl->get_datatype;
741            
742 41 50       130 $pdl *= $bscale if $bscale != 1;
743 41 50       157 $pdl += $bzero if $bzero != 0;
744            
745 41         274 delete $$foo{"BSCALE"}; delete $$foo{"BZERO"};
  41         7143  
746 41         6442 return $pdl;
747             }
748              
749              
750             ##########
751             #
752             # bintable_handlers -- helper table for bintable_row, below.
753             #
754             # Each element of the table is named by the appropriate type letter
755             # from the FITS specification. The value is a list containing the
756             # reading and writing methods.
757             #
758             # This probably ought to be a separate class, but instead it's a tawdry
759             # imitation. Too bad -- except that the ersatz really does run faster than
760             # genuine.
761             #
762             # 0: either a data type or a constructor.
763             # 1: either a length per element or a read method.
764             # 2: either a length per element or a write method.
765             # 3: 'finish' contains finishing-up code or a byte-count to swap.
766             #
767             # Main bintable type handler table.
768             # Elements: (constructor or type, reader or nbytes, writer or nbytes,
769             # finisher or nbytes). The finisher should convert the internal reading
770             # format into the final output format, e.g. by swapping (which is done
771             # automatically in the basic case). Output has row in the 0th dim.
772             #
773             # If present, the constructor should
774             # accept ($rowlen $extra, $nrows, $$size), where $rowlen is the repeat
775             # specifier in the TFORM field, $extra is the extra characters if any
776             # (for added flavor later, if desired), and $nrows is the number of rows in
777             # the table. $$size points to a scalar value that should be incremented by the
778             # size (in bytes) of a single row of the data, for accounting purposes.
779             #
780             # If a read method is specified, it should accept:
781             # ($thing, $rownum, $strptr, $rpt, $extra)
782             # where $rpt is the repeat count and $extra is the extra characters in the
783             # specifier; and it should cut the used characters off the front of the string.
784             #
785             # If a writer is specified it should accept:
786             # ($thing, $row, $rpt, $extra)
787             # and return the generated binary string.
788             #
789             # The finisher just takes the data itself. It should:
790             # * Byteswap
791             # * Condition the data to final dimensional form (if necessary)
792             # * Apply TSCAL/TZERO keywords (as necessary)
793             #
794             # The magic numbers in the table (1,2,4,8, etc.) are kludgey -- they break
795             # the isolation of PDL size and local code -- but they are a part of the FITS
796             # standard. The code will break anyway (damn) on machines that have other
797             # sizes for these datatypes.
798             #
799              
800             $PDL::IO::FITS_bintable_handlers = {
801             'X' => [ byte # Packed bit field
802             , sub {
803             my( $pdl, $row, $strptr ) = @_; # (ignore repeat and extra)
804             my $n = $pdl->dim(0);
805             my $s = unpack( "B".$n, substr(${$strptr}, 0, int(($n+7)/8),''));
806             $s =~ tr/[01]/[\000\001]/;
807             substr( ${$pdl->get_dataref}, $n * $row, length($s)) = $s;
808             }
809             , sub {
810             my( $pdl, $row ) = @_; # Ignore extra and rpt
811             my $n = $pdl->dim(0);
812             my $p2 = byte(($pdl->slice("($row)") != 0));
813             my $s = ${$p2->get_dataref};
814             $s =~ tr/[\000\001]/[01]/;
815             pack( "B".$pdl->dim(0), $s );
816             }
817             , 1
818             ]
819             ,'A' => [ sub { # constructor # String - handle as perl list
820             my($rowlen, $extra, $nrows, $szptr) = @_;
821             my($i,@a);
822             $$szptr += $rowlen;
823             for $i(1..$nrows) { push(@a,' 'x$rowlen); }
824             \@a;
825             }
826             , sub { # reader
827             my( $list, $row, $strptr, $rpt ) = @_;
828             $list->[$row] = substr(${$strptr},0,$rpt,'');
829             }
830             , sub { # writer
831             my($strs, $row, $rpt ) = @_;
832             my $s = substr($strs->[$row],0,$rpt);
833             $s . ' 'x($rpt - length $s);
834             }
835             , undef # no finisher needed
836             ]
837             ,'B' => [ byte, 1, 1, 1 ] # byte
838             ,'L' => [ byte, 1, 1, 1 ] # logical - treat as byte
839             ,'I' => [ short, 2, 2, 2 ] # short (no unsigned shorts?)
840             ,'J' => [ long, 4, 4, 4 ] # long
841             ,'K' => [ longlong,8, 8, 8 ] # longlong
842             ,'E' => [ float, 4, 4, 4 ] # single-precision
843             ,'D' => [ double, 8, 8, 8 ] # double-precision
844             ,'C' => [ sub { _nucomplx(float, eval '@_') }, sub { _rdcomplx(float, eval '@_') },
845             sub { _wrcomplx(float, eval '@_') }, sub { _fncomplx(float, eval '@_') }
846             ]
847             ,'M' => [ sub { _nucomplx(double, eval '@_') }, sub { _rdcomplx(double, eval '@_') },
848             sub { _wrcomplx(double, eval '@_') }, sub { _fncomplx(double, eval '@_') }
849             ]
850             ,'PB' => [ sub { _nuP(byte, eval '@_') }, sub { _rdP(byte, eval '@_') },
851             sub { _wrP(byte, eval '@_') }, sub { _fnP(byte, eval '@_') }
852             ]
853             ,'PL' => [ sub { _nuP(byte, eval '@_') }, sub { _rdP(byte, eval '@_') },
854             sub { _wrP(byte, eval '@_') }, sub { _fnP(byte, eval '@_') }
855             ]
856             ,'PI' => [ sub { _nuP(short, eval '@_') }, sub { _rdP(short, eval '@_') },
857             sub { _wrP(short, eval '@_') }, sub { _fnP(short, eval '@_') }
858             ]
859             ,'PJ' => [ sub { _nuP(long, eval '@_') }, sub { _rdP(long, eval '@_') },
860             sub { _wrP(long, eval '@_') }, sub { _fnP(long, eval '@_') }
861             ]
862             ,'PE' => [ sub { _nuP(float, eval '@_') }, sub { _rdP(float, eval '@_') },
863             sub { _wrP(float, eval '@_') }, sub { _fnP(float, eval '@_') }
864             ]
865             ,'PD' => [ sub { _nuP(double, eval '@_') }, sub { _rdP(double, eval '@_') },
866             sub { _wrP(double, eval '@_') }, sub { _fnP(double, eval '@_') }
867             ]
868             };
869              
870              
871             ##############################
872             # Helpers for complex numbers (construct/read/write/finish)
873             sub _nucomplx { # complex-number constructor
874 0     0   0 my($type, $rowlen, $extra, $nrows, $szptr) = @_;
875 0         0 $szptr += PDL::Core::howbig($type) * $nrows * 2;
876 0         0 return PDL->new_from_specification($type,2,$rowlen,$nrows);
877             }
878             sub _rdcomplx { # complex-number reader
879 0     0   0 my( $type, $pdl, $row, $strptr, $rpt ) = @_; # ignore extra
880 0         0 my $s = $pdl->get_dataref;
881 0         0 my $rlen = 2 * PDL::Core::howbig($type) * $rpt;
882 0         0 substr($$s, $row*$rlen, $rlen) = substr($strptr, 0, $rlen, '');
883             }
884             sub _wrcomplx { # complex-number writer
885 0     0   0 my( $type, $pdl, $row, $rpt ) = @_; # ignore extra
886 0         0 my $rlen = 2 * PDL::Core::howbig($type) * $rpt;
887 0         0 substr( ${$pdl->get_dataref}, $rlen * $row, $rlen );
  0         0  
888             }
889             sub _fncomplx { # complex-number finisher-upper
890 0     0   0 my( $type, $pdl, $n, $hdr, $opt) = shift;
891 0         0 eval 'bswap'.(PDL::Core::howbig($type)).'($pdl)';
892             print STDERR "Ignoring poorly-defined TSCAL/TZERO for complex data in col. $n (".$hdr->{"TTYPE$n"}.").\n"
893 0 0 0     0 if( length($hdr->{"TSCAL$n"}) or length($hdr->{"TZERO$n"}) );
894 0         0 return $pdl->reorder(2,1,0);
895             }
896              
897             ##############################
898             # Helpers for variable-length array types (construct/read/write/finish)
899             # These look just like the complex-number case, except that they use $extra to determine the
900             # size of the 0 dimension.
901             sub _nuP {
902 0     0   0 my( $type, $rowlen, $extra, $nrows, $szptr, $hdr, $i, $tbl ) = @_;
903 0         0 $extra =~ s/\((.*)\)/$1/; # strip parens from $extra in-place
904 0         0 $$szptr += 8;
905              
906 0 0       0 if($rowlen != 1) {
907 0         0 die("rfits: variable-length record has a repeat count that isn't unity! (got $rowlen); I give up.");
908             }
909              
910             # declare the PDL. Fill it with the blank value or (failing that) 0.
911             # Since P repeat count is required to be 0 or 1, we don't need an additional dimension for the
912             # repeat count -- the variable-length rows take that role.
913 0         0 my $pdl = PDL->new_from_specification($type, $extra, $nrows);
914 0   0     0 $pdl .= ($hdr->{"TNULL$i"} || 0);
915              
916 0         0 my $lenpdl = zeroes(long, $nrows);
917 0         0 $tbl->{"len_".$hdr->{"TTYPE$i"}} = $lenpdl;
918              
919 0         0 return $pdl;
920             }
921             sub _rdP {
922 0     0   0 my( $type, $pdl, $row, $strptr, $rpt, $extra, $heap_ptr, $tbl, $i ) = @_;
923 0         0 $extra =~ s/\((.*)\)/$1/;
924 0         0 my $s = $pdl->get_dataref;
925              
926             # Read current offset and length
927 0         0 my $oflen = pdl(long,0,0);
928 0         0 my $ofs = $oflen->get_dataref;
929 0         0 substr($$ofs,0,8) = substr($$strptr, 0, 8, '');
930 0         0 $oflen->upd_data;
931 0         0 bswap4($oflen);
932            
933             # Now get 'em
934 0         0 my $rlen = $extra * PDL::Core::howbig($type); # rpt should be unity, otherwise we'd have to multiply it in.
935 0         0 my $readlen = $oflen->at(0) * PDL::Core::howbig($type);
936              
937             # Store the length of this row in the header field.
938 0         0 $tbl->{"len_".$tbl->{hdr}->{"TTYPE$i"}}->dice_axis(0,$row) .= $oflen->at(0);
939              
940 0 0       0 print "_rdP: pdl is ",join("x",$pdl->dims),"; reading row $row - readlen is $readlen\n"
941             if($PDL::debug);
942              
943             # Copy the data into the output PDL.
944 0         0 my $of = $oflen->at(1);
945 0         0 substr($$s, $row*$rlen, $readlen) = substr($$heap_ptr, $of, $readlen);
946 0         0 $pdl->upd_data;
947             }
948              
949             sub _wrP {
950 0     0   0 die "This code path should never execute - you are trying to write a variable-length array via direct handler, which is wrong. Check the code path in PDL::wfits.\n";
951             }
952              
953             sub _fnP {
954 0     0   0 my( $type, $pdl, $n, $hdr, $opt ) = @_;
955 0         0 my $post = PDL::Core::howbig($type);
956 0 0       0 unless( isbigendian() ) {
957 0 0       0 if( $post == 2 ) { bswap2($pdl); }
  0 0       0  
    0          
    0          
958 0         0 elsif( $post == 4 ) { bswap4($pdl); }
959 0         0 elsif( $post == 8 ) { bswap8($pdl); }
960             elsif( $post != 1 ) {
961 0         0 print STDERR "Unknown swapsize $post! This is a bug. You (may) lose..\n";
962             }
963             }
964              
965 0 0       0 my $tzero = defined($hdr->{"TZERO$n"}) ? $hdr->{"TZERO$n"} : 0.0;
966 0 0       0 my $tscal = defined($hdr->{"TSCAL$n"}) ? $hdr->{"TSCAL$n"} : 1.0;
967 0         0 my $valid_tzero = ($tzero != 0.0);
968 0         0 my $valid_tscal = ($tscal != 1.0);
969 0 0 0     0 if( length($hdr->{"TZERO$n"}) or length($hdr->{"TSCAL$n"})) {
970 0         0 print STDERR "Ignoring TSCAL/TZERO keywords for binary table array column - sorry, my mind is blown!\n";
971             }
972 0         0 return $pdl->mv(-1,0);
973             }
974              
975             ##############################
976             #
977             # _rfits_bintable -- snarf up a binary table, returning the named columns
978             # in a hash ref, each element of which is a PDL or list ref according to the
979             # header.
980             #
981              
982             sub _rfits_bintable ($$$$) {
983 3     3   7 my $fh = shift;
984 3         22 my $hdr = shift;
985 3         7 my $opt = shift;
986             ##shift; ### (ignore $pdl argument)
987              
988 3 50       15 print STDERR "Warning: BINTABLE extension should have BITPIX=8, found ".$hdr->{BITPIX}.". Winging it...\n" unless($hdr->{BITPIX} == 8);
989            
990             ### Allocate the main table hash
991 3         218 my $tbl = {}; # Table is indexed by name
992 3         15 $tbl->{hdr} = $hdr;
993 3         9 $tbl->{tbl} = 'binary';
994              
995 3         8 my $tmp = []; # Temporary space is indexed by col. no.
996            
997            
998             ### Allocate all the columns of the table, checking for consistency
999             ### and name duplication.
1000            
1001 3 50       12 barf "Binary extension has no fields (TFIELDS=0)" unless($hdr->{TFIELDS});
1002 3         224 my $rowlen = 0;
1003            
1004 3         11 for my $i(1..$hdr->{TFIELDS}) {
1005 9         236 my $iter;
1006 9   50     44 my $name = $tmp->[$i]->{name} = $hdr->{"TTYPE$i"} || "COL";
1007            
1008             ### Allocate some temp space for dealing with this column
1009 9         741 my $tmpcol = $tmp->[$i] = {};
1010            
1011             ### Check for duplicate name and change accordingly...
1012 9   33     46 while( defined( $tbl->{ $name } ) || ($name eq "COL") ) {
1013 0         0 $iter++;
1014 0         0 $name = ($hdr->{"TTYPE$i"} )."_$iter";
1015             }
1016            
1017             # (Check avoids scrozzling comment fields unnecessarily)
1018 9 50       38 $hdr->{"TTYPE$i"} = $name unless($hdr->{"TTYPE$i"} eq $name);
1019 9         645 $tmpcol->{name} = $name;
1020              
1021 9 50       34 if( ($hdr->{"TFORM$i"}) =~ m/(\d*)(P?.)(.*)/ ) {
1022 9         685 ($tmpcol->{rpt}, $tmpcol->{type}, $tmpcol->{extra}) = ($1,$2,$3);
1023             # added by DJB 03.18/04 - works for my data file but is it correct?
1024 9   50     25 $tmpcol->{rpt} ||= 1;
1025             } else {
1026             barf "Couldn't parse BINTABLE form '"
1027             . $hdr->{"TFORM$i"}
1028             . "' for column $i ("
1029             . $hdr->{"TTYPE$i"}
1030 0 0       0 . ")\n" if($hdr->{"TFORM$i"});
1031 0         0 barf "BINTABLE header is missing a crucial field, TFORM$i. I give up.\n";
1032             }
1033              
1034             # "A bit array consists of an integral number of bytes with trailing bits zero"
1035 9 50       30 $tmpcol->{rpt} = PDL::ceil($tmpcol->{rpt}/8) if ($tmpcol->{type} eq 'X');
1036              
1037             $tmpcol->{handler} = # sic - assignment
1038             $PDL::IO::FITS_bintable_handlers->{ $tmpcol->{type} }
1039             or
1040             barf "Unknown type ".$hdr->{"TFORM$i"}." in BINTABLE column $i "."("
1041 9 50       43 . $hdr->{"TTYPE$i"}
1042             . ")\n That invalidates the byte count, so I give up.\n" ;
1043            
1044            
1045             ### Allocate the actual data space and increment the row length
1046            
1047 9         21 my $foo = $tmpcol->{handler}->[0];
1048 9 100       32 if( ref ($foo) eq 'CODE' ) {
1049             $tmpcol->{data} = $tbl->{$name} =
1050 1         5 &{$foo}( $tmpcol->{rpt}
1051             , $tmpcol->{extra}
1052             , $hdr->{NAXIS2}
1053 1         5 , \$rowlen
1054             , $hdr # hdr and column number are passed in, in case extra info needs to be gleaned.
1055             , $i
1056             , $tbl
1057             );
1058             } else {
1059             $tmpcol->{data} = $tbl->{$name} =
1060             PDL->new_from_specification(
1061             $foo
1062             , $tmpcol->{rpt},
1063 8   50     29 , $hdr->{NAXIS2} || 1
1064             );
1065              
1066 8         50 $rowlen += PDL::Core::howbig($foo) * $tmpcol->{rpt};
1067             }
1068            
1069 9 50       31 print "Prefrobnicated col. $i "."(".$hdr->{"TTYPE$i"}.")\ttype is ".$hdr->{"TFORM$i"}."\t length is now $rowlen\n" if($PDL::debug);
1070            
1071            
1072             } ### End of prefrobnication loop...
1073            
1074              
1075             barf "Calculated row length is $rowlen, hdr claims ".$hdr->{NAXIS1}
1076             . ". Giving up. (Set \$PDL::debug for more detailed info)\n"
1077 3 50       16 if($rowlen != $hdr->{NAXIS1});
1078            
1079             ### Snarf up the whole extension, and pad to 2880 bytes...
1080 3         221 my ($rawtable, $heap, $n1, $n2);
1081              
1082             # n1 gets number of bytes in table plus gap
1083 3         12 $n1 = $hdr->{NAXIS1} * $hdr->{NAXIS2};
1084 3 50       380 if($hdr->{THEAP}) {
1085 0 0       0 if($hdr->{THEAP} < $n1) {
1086 0         0 die("Inconsistent THEAP keyword in binary table\n");
1087             } else {
1088 0         0 $n1 = $hdr->{THEAP};
1089             }
1090             }
1091              
1092             # n2 gets number of bytes in heap (PCOUNT - gap).
1093 3 50       155 $n2 = $hdr->{PCOUNT} + ($hdr->{THEAP} ? ($hdr->{NAXIS1}*$hdr->{NAXIS2} - $hdr->{THEAP}) : 0);
1094 3         426 $n2 = ($n1+$n2-1)+2880 - (($n1+$n2-1) % 2880) - $n1;
1095              
1096 3 50       13 print "Reading $n1 bytes of table data and $n2 bytes of heap data....\n"
1097             if($PDL::verbose);
1098 3         30 $fh->read($rawtable, $n1);
1099              
1100 3 50       46 if($n2) {
1101 3         10 $fh->read($heap, $n2)
1102             } else {
1103 0         0 $heap = which(pdl(0)); # empty PDL
1104             }
1105              
1106             ### Frobnicate the rows, one at a time.
1107 3         118 for my $row(0..$hdr->{NAXIS2}-1) {
1108 12         963 my $prelen = length($rawtable);
1109 12         39 for my $i(1..$hdr->{TFIELDS}) {
1110 36         847 my $tmpcol = $tmp->[$i];
1111 36         53 my $reader = $tmpcol->{handler}->[1];
1112 36 100       99 if(ref $reader eq 'CODE') {
    50          
1113 4         13 &{$reader}( $tmpcol->{data}
1114             , $row
1115             , \$rawtable
1116             , $tmpcol->{rpt}
1117             , $tmpcol->{extra}
1118 4         12 , \$heap
1119             , $tbl
1120             , $i
1121             );
1122             } elsif(ref $tmpcol->{data} eq 'PDL') {
1123 32         46 my $rlen = $reader * $tmpcol->{rpt};
1124              
1125 32         59 substr( ${$tmpcol->{data}->get_dataref()}, $rlen * $row, $rlen ) =
  32         103  
1126             substr( $rawtable, 0, $rlen, '');
1127 32         95 $tmpcol->{data}->upd_data;
1128              
1129             } else {
1130 0         0 die ("rfits: Bug detected: inconsistent types in BINTABLE reader\n");
1131             }
1132              
1133             } # End of TFIELDS loop
1134              
1135 12 50       53 if(length($rawtable) ne $prelen - $hdr->{NAXIS1}) {
1136 0         0 die "rfits BINTABLE: Something got screwed up -- expected a length of $prelen - $hdr->{NAXIS1}, got ".length($rawtable).". Giving up.\n";
1137             }
1138             } # End of NAXIS2 loop
1139              
1140             #
1141             # Note: the above code tickles a bug in most versions of the emacs
1142             # prettyprinter. The following "for my $i..." should be indented
1143             # two spaces.
1144             #
1145              
1146             ### Postfrobnicate the columns.
1147 3         223 for my $i(1..$hdr->{TFIELDS}) { # Postfrobnication loop
1148 9         210 my $tmpcol = $tmp->[$i];
1149 9         22 my $post = $tmpcol->{handler}->[3];
1150            
1151 9 50       62 if(ref $post eq 'CODE') {
    100          
    50          
1152             # Do postprocessing on all special types
1153            
1154 0         0 $tbl->{$tmpcol->{name}} = &$post($tmpcol->{data}, $i, $hdr, $opt);
1155            
1156             } elsif( (ref ($tmpcol->{data})) eq 'PDL' ) {
1157             # Do standard PDL-type postprocessing
1158            
1159             ## Is this call to upd_data necessary?
1160             ## I think not. (reinstate if there are bugs)
1161             # $tmpcol->{data}->upd_data;
1162            
1163            
1164             # Do swapping as necessary
1165 8 50       40 unless( isbigendian() ) {
1166 8 100       57 if( $post == 2 ) { bswap2($tmpcol->{data}); }
  1 100       26  
    100          
    50          
1167 4         78 elsif( $post == 4 ) { bswap4($tmpcol->{data}); }
1168 2         38 elsif( $post == 8 ) { bswap8($tmpcol->{data}); }
1169             elsif( $post != 1 ) {
1170             print STDERR "Unknown swapsize $post for column $i ("
1171 0         0 . $tmpcol->{name} . ")! This is a bug. Winging it.\n";
1172             }
1173             }
1174            
1175             # Apply scaling and badval keys, which are illegal for A, L, and X
1176             # types but legal for anyone else. (A shouldn't be here, L and X
1177             # might be)
1178            
1179 8 50       28 if($opt->{bscale}) {
1180 8 50       51 my $tzero = defined($hdr->{"TZERO$i"}) ? $hdr->{"TZERO$i"} : 0.0;
1181 8 50       462 my $tscal = defined($hdr->{"TSCAL$i"}) ? $hdr->{"TSCAL$i"} : 1.0;
1182            
1183             # The $valid_ flags let us avoid unnecessary arithmetic.
1184 8         409 my $valid_tzero = ($tzero != 0.0);
1185 8         18 my $valid_tscal = ($tscal != 1.0);
1186            
1187 8 50 33     44 if ( $valid_tzero or $valid_tscal ) {
1188 0 0       0 if ( $tmpcol->{type} =~ m/[ALX]/i ) {
1189            
1190             print STDERR "Ignoring illegal TSCAL/TZERO keywords for col $i (" .
1191 0         0 $tmpcol->{name} . "); type is $tmpcol->{type})\n";
1192            
1193             } else { # Not an illegal type -- do the scaling
1194            
1195             # (Normal execution path)
1196             # Use PDL's cleverness to work out the final datatype...
1197            
1198 0         0 my $tmp;
1199 0         0 my $pdl = $tmpcol->{data};
1200            
1201 0 0       0 if($pdl->badflag() == 0) {
    0          
1202            
1203 0         0 $tmp = $pdl->flat()->slice("0:0");
1204            
1205             } elsif($pdl->ngood > 0) {
1206            
1207 0         0 my $index = which( $pdl->flat()->isbad()==0 )->at(0);
1208 0         0 $tmp = $pdl->flat()->slice("${index}:${index}");
1209            
1210             } else { # Do nothing if it's all bad....
1211 0         0 $tmp = $pdl;
1212             }
1213            
1214             # Figure out the type by scaling the single element.
1215 0         0 $tmp = ($tmp - $tzero) * $tscal;
1216            
1217             # Convert the whole PDL as necessary for the scaling.
1218 0 0       0 $tmpcol->{data} = $pdl->convert($tmp->type)
1219             if($tmp->get_datatype != $pdl->get_datatype);
1220            
1221             # Do the scaling.
1222 0         0 $tmpcol->{data} -= $tzero;
1223 0         0 $tmpcol->{data} *= $tscal;
1224              
1225             } # End of legal-type conditional
1226             } # End of valid_ conditional
1227              
1228 8         49 delete $hdr->{"TZERO$i"};
1229 8         2094 delete $hdr->{"TSCAL$i"};
1230            
1231             } else { # $opt->{bscale} is zero; don't scale.
1232             # Instead, copy factors into individual column headers.
1233 0         0 my %foo = ("TZERO$i"=>"BZERO",
1234             "TSCAL$i"=>"BSCALE",
1235             "TUNIT$i"=>"BUNIT");
1236 0         0 for my $x(keys %foo) {
1237             $tmpcol->{data}->hdr->{$foo{$x}} = $hdr->{$x}
1238 0 0       0 if( defined($hdr->{$x}) );
1239             }
1240             } # End of bscale checking...
1241              
1242             # Try to grab a TDIM dimension list...
1243 8         1831 my @tdims = ();
1244 8         42 $tmpcol->{data}->hdrcpy(1);
1245              
1246 8 50       37 if(exists($hdr->{"TDIM$i"})) {
1247 0 0       0 if($hdr->{"TDIM$i"} =~ m/\((\s*\d+(\s*\,\s*\d+)*\s*)\)/) {
1248 0         0 my $x = $1;
1249 0         0 @tdims = map { $_+0 } split(/\,/,$x);
  0         0  
1250 0         0 my $tdims = pdl(@tdims);
1251 0         0 my $tds = $tdims->prodover;
1252 0 0       0 if($tds > $tmpcol->{data}->dim(0)) {
    0          
1253 0         0 die("rfits: TDIM$i is too big in binary table. I give up.\n");
1254             } elsif($tds < $tmpcol->{data}->dim(0)) {
1255 0         0 print STDERR "rfits: WARNING: TDIM$i is too small in binary table. Carrying on...\n";
1256             }
1257              
1258 0         0 $tmpcol->{data}->hdrcpy(1);
1259 0         0 my $td = $tmpcol->{data}->xchg(0,1);
1260 0         0 $tbl->{$tmpcol->{name}} = $td->reshape($td->dim(0),@tdims);
1261             } else {
1262 0         0 print STDERR "rfits: WARNING: invalid TDIM$i field in binary table. Ignoring.\n";
1263             }
1264             } else {
1265             # Copy the PDL out to the table itself.
1266 8 50 33     138 if($hdr->{NAXIS2} > 0 && $tmpcol->{rpt}>0) {
1267             $tbl->{$tmpcol->{name}} =
1268             ( ( $tmpcol->{data}->dim(0) == 1 )
1269             ? $tmpcol->{data}->slice("(0)")
1270 8 50       684 : $tmpcol->{data}->xchg(0,1)
1271             );
1272             }
1273             }
1274              
1275             # End of PDL postfrobnication case
1276             } elsif(defined $post) {
1277            
1278             print STDERR "Postfrobnication bug detected in column $i ("
1279 0         0 . $tmpcol->{name}. "). Winging it.\n";
1280            
1281             }
1282             } # End of postfrobnication loop over columns
1283              
1284             ### Check whether this is actually a compressed image, in which case we hand it off to the image decompressor
1285 3 0 33     15 if($hdr->{ZIMAGE} && $hdr->{ZCMPTYPE} && $opt->{expand}) {
      0        
1286 0         0 eval 'use PDL::Compression;';
1287 0 0       0 if($@) {
1288 0         0 die "rfits: error while loading PDL::Compression to unpack tile-compressed image.\n\t$@\n\tUse option expand=>0 to get the binary table.\n";
1289             }
1290 0         0 return _rfits_unpack_zimage($tbl,$opt);
1291             }
1292              
1293             ### Done!
1294 3         211 return $tbl;
1295             }
1296              
1297              
1298             ##############################
1299             ##############################
1300             #
1301             # _rfits_unpack_zimage - unpack a binary table that actually contains a compressed image
1302             #
1303             # This is implemented to support the partial spec by White, Greenfield, Pence, & Tody dated Oct 21, 1999,
1304             # with reverse-engineered bits from the CFITSIO3240 library where the spec comes up short.
1305             #
1306              
1307             ## keyword is a compression algorithm name; value is an array ref containing tile compressor/uncompressor.
1308             ## The compressor/uncompressor takes (nx, ny, data) and returns the compressed/uncompressed data.
1309             ## The four currently (2010) supported-by-CFITSIO compressors are listed. Not all have been ported, hence
1310             ## the "undef"s in the table. --CED.
1311              
1312             ## Master jump table for compressors/uncompressors.
1313             ## 0 element of each array ref is the compressor; 1 element is the uncompressor.
1314             ## Uncompressed tiles are reshaped to rows of a tile table handed in (to the compressor)
1315             ## or out (of the uncompressor); actual tile shape is fed in as $params->{tiledims}, so
1316             ## higher-than-1D compression algorithms can be used.
1317             our $tile_compressors = {
1318             'GZIP_1' => undef
1319             , 'RICE_1' => [ ### RICE_1 compressor
1320             sub { my ($tiles, $tbl, $params) = @_;
1321             my ($compressed,$len) = $tiles->rice_compress($params->{BLOCKSIZE} || 32);
1322             $tbl->{ZNAME1} = "BLOCKSIZE";
1323             $tbl->{ZVAL1} = $params->{BLOCKSIZE};
1324             $tbl->{ZNAME2} = "BYTEPIX";
1325             $tbl->{ZVAL2} = PDL::howbig($tiles->get_datatype);
1326             # Convert the compressed data to a byte array...
1327             if($tbl->{ZVAL2} != 1) {
1328             my @dims = $compressed->dims;
1329             $dims[0] *= $tbl->{ZVAL2};
1330             my $cd2 = zeroes( byte, @dims );
1331             my $cdr = $compressed->get_dataref;
1332             my $cd2r = $cd2->get_dataref;
1333             $$cd2r = $$cdr;
1334             $cd2->upd_data;
1335             $compressed = $cd2;
1336             }
1337             $tbl->{COMPRESSED_DATA} = $compressed->mv(0,-1);
1338             $tbl->{len_COMPRESSED_DATA} = $len;
1339             },
1340             ### RICE_1 expander
1341             sub { my ($tilesize, $tbl, $params) = @_;
1342             my $compressed = $tbl->{COMPRESSED_DATA} -> mv(-1,0);
1343             my $bytepix = $params->{BYTEPIX} || 4;
1344              
1345             # Put the compressed tile bitstream into a variable of appropriate type.
1346             # This works by direct copying of the PDL data, which sidesteps local
1347             # byteswap issues in the usual case that the compressed stream is type
1348             # byte. But it does add the extra complication that we have to pad the
1349             # compressed array out to a factor-of-n elements in certain cases.
1350              
1351             if( PDL::howbig($compressed->get_datatype) != $bytepix ) {
1352             my @dims = $compressed->dims;
1353             my $newdim0;
1354             my $scaledim0;
1355              
1356             $scaledim0 = $dims[0] * PDL::howbig($compressed->get_datatype) / $bytepix;
1357             $newdim0 = pdl($scaledim0)->ceil;
1358              
1359             if($scaledim0 != $newdim0) {
1360             my $padding = zeroes($compressed->type,
1361             ($newdim0-$scaledim0) * $bytepix / PDL::howbig($compressed->get_datatype),
1362             @dims[1..$#dims]
1363             );
1364             $compressed = $compressed->append($padding);
1365             }
1366            
1367             my $c2 = zeroes( $type_table_2->{$bytepix * 8}, $newdim0, @dims[1..$#dims] );
1368              
1369             my $c2dr = $c2->get_dataref;
1370             my $cdr = $compressed->get_dataref;
1371             substr($$c2dr,0,length($$cdr)) = $$cdr;
1372             $c2->upd_data;
1373             $compressed = $c2;
1374             }
1375              
1376             return $compressed->rice_expand( $tilesize, $params->{BLOCKSIZE} || 32);
1377             }
1378             ]
1379             , 'PLIO_1' => undef
1380             , 'HCOMPRESS_1' => undef
1381             };
1382              
1383             ## List of the eight mandatory keywords and their ZIMAGE preservation pigeonholes, for copying after we
1384             ## expand an image.
1385             our $hdrconv = {
1386             "ZSIMPLE" => "SIMPLE",
1387             "ZTENSION" => "XTENSION",
1388             "ZEXTEND" => "EXTEND",
1389             "ZBLOCKED" => "BLOCKED",
1390             "ZPCOUNT" => "PCOUNT",
1391             "ZGCOUNT" => "GCOUNT",
1392             "ZHECKSUM" => "CHECKSUM",
1393             "ZDATASUM" => "DATASUM"
1394             };
1395              
1396              
1397             sub _rfits_unpack_zimage($$$) {
1398 0     0   0 my $tbl = shift;
1399 0         0 my $opt = shift;
1400              
1401 0         0 my $hdr = $tbl->{hdr};
1402              
1403 0         0 my $tc = $tile_compressors->{$hdr->{ZCMPTYPE}};
1404 0 0       0 unless(defined $tc) {
1405 0         0 print STDERR "WARNING: rfits: Compressed image has unsupported comp. type ('$hdr->{ZCMPTYPE}').\n";
1406 0         0 return $tbl;
1407             }
1408              
1409             #############
1410             # Declare the output image
1411 0         0 my $type;
1412 0 0       0 unless($type_table->{$hdr->{ZBITPIX}}) {
1413 0         0 print STDERR "WARNING: rfits: unrecognized ZBITPIX value $hdr->{ZBITPIX} in compressed image. Assuming -64.\n";
1414 0         0 $type = $type_table_2->{-64};
1415             } else {
1416 0         0 $type = $type_table_2->{$hdr->{ZBITPIX}};
1417             }
1418 0         0 my @dims;
1419 0         0 for my $i(1..$hdr->{ZNAXIS}) {
1420 0         0 push(@dims,$hdr->{"ZNAXIS$i"});
1421             }
1422              
1423 0         0 my $pdl = PDL->new_from_specification( $type, @dims );
1424              
1425             ############
1426             # Calculate tile size and allocate a working tile.
1427 0         0 my @tiledims;
1428 0         0 for my $i(1..$hdr->{ZNAXIS}) {
1429 0 0       0 if($hdr->{"ZTILE$i"}) {
1430 0         0 push(@tiledims, $hdr->{"ZTILE$i"});
1431             } else {
1432 0 0       0 push(@tiledims, (($i==1) ? $hdr->{ZNAXIS1} : 1) );
1433             }
1434             }
1435              
1436              
1437             # my $tile = PDL->new_from_specification( $type, @tiledims );
1438 0         0 my $tiledims = pdl(@tiledims);
1439 0         0 my $tilesize = $tiledims->prodover;
1440             ###########
1441             # Calculate tile counts and compare to the number of stored tiles
1442 0         0 my $ntiles = ( pdl(@dims) / pdl(@tiledims) )->ceil;
1443 0         0 my $tilecount = $ntiles->prodover;
1444              
1445 0 0       0 if($tilecount != $tbl->{COMPRESSED_DATA}->dim(0)) {
1446 0         0 printf STDERR "WARNING: rfits: compressed data has $hdr->{NAXIS2} rows; we expected $tilecount (",join("x",list $ntiles),"). Winging it...\n";
1447             }
1448              
1449             ##########
1450             # Quantization - ignore for now
1451 0 0       0 if($hdr->{ZQUANTIZ}) {
1452 0         0 printf STDERR "WARNING: rfits: ignoring quantization/dithering (ZQUANTIZ=$hdr->{ZQUANTIZ})\n";
1453             }
1454            
1455             ##########
1456             # Snarf up compression parameters
1457 0         0 my $params = {};
1458 0         0 my $i = 1;
1459 0         0 while( $hdr->{"ZNAME$i"} ) {
1460 0         0 $params->{ $hdr->{"ZNAME$i"} } = $hdr->{"ZVAL$i"};
1461 0         0 $i++;
1462             }
1463              
1464             ##########
1465             # Enumerate tile coordinates for looping, and the corresponding row number
1466 0         0 my ($step, @steps, $steps);
1467 0         0 $step = 1;
1468 0         0 for my $i(0..$ntiles->nelem-1) {
1469 0         0 push(@steps, $step);
1470 0         0 $step *= $ntiles->at($i);
1471             }
1472 0         0 $step = pdl(@steps);
1473              
1474             # $tiledex is 2-D (coordinate-index, list-index) and enumerates all tiles by image
1475             # location; $tilerow is 1-D (list-index) and enumerates all tiles by row in the bintable
1476 0         0 my $tiledex = PDL::ndcoords($ntiles->list)->mv(0,-1)->clump($ntiles->dim(0))->mv(-1,0);
1477 0         0 $TMP::tiledex = $tiledex;
1478 0         0 my $tilerow = ($tiledex * $step)->sumover;
1479            
1480             ##########
1481             # Restore all the tiles at once
1482 0         0 my $tiles = &{$tc->[1]}( $tilesize, $tbl, $params ); # gets a (tilesize x ntiles) output
  0         0  
1483 0         0 my $patchup = which($tbl->{len_COMPRESSED_DATA} <= 0);
1484 0 0       0 if($patchup->nelem) {
1485 0 0       0 unless(defined $tbl->{UNCOMPRESSED_DATA}) {
1486 0         0 die "rfits: need some uncompressed data for missing compressed rows, but none were found!\n";
1487             }
1488 0 0       0 if($tbl->{UNCOMPRESSED_DATA}->dim(1) != $tilesize) {
1489 0         0 die "rfits: tile size is $tilesize, but uncompressed data rows have size ".$tbl->{UNCOMPRESSED_DATA}->dim(1)."\n";
1490             }
1491 0         0 $tiles->dice_axis(1,$patchup) .= $tbl->{UNCOMPRESSED_DATA}->dice_axis(0,$patchup)->xchg(0,1);
1492             }
1493              
1494             ##########
1495             # Slice up the output image plane into tiles, and use the threading engine
1496             # to assign everything to them.
1497 0         0 my $cutup = $pdl->range( $tiledex, [@tiledims], 't') # < ntiles, tilesize0..tilesizen >
1498             ->mv(0,-1) # < tilesize0..tilesizen, ntiles >
1499             ->clump($tiledims->nelem); # < tilesize, ntiles >
1500              
1501 0         0 $cutup .= $tiles; # dump all the tiles at once into the image - they flow back to $pdl.
1502 0         0 undef $cutup; # sever connection to prevent expensive future dataflow.
1503              
1504             ##########
1505             # Perform scaling if necessary ( Just the ZIMAGE quantization step )
1506             # bscaling is handled farther down with treat_bscale.
1507              
1508 0 0       0 $pdl *= $hdr->{ZSCALE} if defined($hdr->{ZSCALE});
1509 0 0       0 $pdl += $hdr->{ZZERO} if defined($hdr->{ZZERO});
1510              
1511             ##########
1512             # Put the FITS header into the newly reconstructed image.
1513 0         0 delete $hdr->{PCOUNT};
1514 0         0 delete $hdr->{GCOUNT};
1515              
1516             # Copy the mandated name-conversions
1517 0         0 for my $k(keys %$hdrconv) {
1518 0 0       0 if($hdr->{$k}) {
1519 0         0 $hdr->{$hdrconv->{$k}} = $hdr->{$k};
1520 0         0 delete $hdr->{$k};
1521             }
1522             }
1523              
1524             # Copy the ZNAXIS* keywords to NAXIS*
1525 0         0 foreach (grep /^NAXIS/,keys %$hdr){
1526 0 0 0     0 if (exists($hdr->{'Z'.$_}) && defined($hdr->{'Z'.$_})){
1527 0         0 $hdr->{$_} = $hdr->{'Z'.$_};
1528             }
1529             }
1530              
1531             # Clean up the ZFOO extensions and table cruft
1532 0         0 for my $k(keys %{$hdr}) {
  0         0  
1533 0 0 0     0 delete $hdr->{$k} if(
      0        
      0        
1534             $k=~ m/^Z/ ||
1535             $k eq "TFIELDS" ||
1536             $k =~ m/^TTYPE/ ||
1537             $k =~ m/^TFORM/
1538             );
1539             }
1540              
1541 0 0 0     0 if(exists $hdr->{BSCALE} || exists $hdr->{BLANK}) {
1542 0         0 $pdl = treat_bscale($pdl, $hdr);
1543             }
1544 0         0 $pdl->sethdr($hdr);
1545 0         0 $pdl->hdrcpy($opt->{hdrcpy});
1546              
1547 0         0 return $pdl;
1548             }
1549              
1550              
1551            
1552              
1553             =head2 wfits()
1554              
1555             =for ref
1556              
1557             Simple PDL FITS writer
1558              
1559             =for example
1560              
1561             wfits $pdl, 'filename.fits', [$BITPIX], [$COMPRESSION_OPTIONS];
1562             wfits $hash, 'filename.fits', [$OPTIONS];
1563             $pdl->wfits('foo.fits',-32);
1564              
1565             Suffix magic:
1566              
1567             # Automatically compress through pipe to gzip
1568             wfits $pdl, 'filename.fits.gz';
1569             # Automatically compress through pipe to compress
1570             wfits $pdl, 'filename.fits.Z';
1571              
1572             Tilde expansion:
1573              
1574             #expand leading ~ to home directory (using glob())
1575             wfits $pdl, '~/filename.fits';
1576              
1577             =over 3
1578              
1579             =item * Ordinary (PDL) data handling:
1580              
1581             If the first argument is a PDL, then the PDL is written out as an
1582             ordinary FITS file with a single Header/Data Unit of data.
1583              
1584             $BITPIX is then optional and coerces the output data type according to
1585             the standard FITS convention for the BITPIX field (with positive
1586             values representing integer types and negative values representing
1587             floating-point types).
1588              
1589             If C<$pdl> has a FITS header attached to it (actually, any hash that
1590             contains a C<< SIMPLE=>T >> keyword), then that FITS header is written
1591             out to the file. The image dimension tags are adjusted to the actual
1592             dataset. If there's a mismatch between the dimensions of the data and
1593             the dimensions in the FITS header, then the header gets corrected and
1594             a warning is printed.
1595              
1596             If C<$pdl> is a slice of another PDL with a FITS header already
1597             present (and header copying enabled), then you must be careful.
1598             C will remove any extraneous C keywords (per the FITS
1599             standard), and also remove the other keywords associated with that
1600             axis: C, C, C, C, and C. This
1601             may cause confusion if the slice is NOT out of the last dimension:
1602             C and you would be best off adjusting
1603             the header yourself before calling C.
1604              
1605             You can tile-compress images according to the CFITSIO extension to the
1606             FITS standard, by adding an option hash to the arguments:
1607              
1608             =over 3
1609              
1610             =item compress
1611              
1612             CURRENTLY UNIMPLEMENTED. Below describes the envisioned usage.
1613              
1614             This can be either unity, in which case Rice compression is used,
1615             or a (case-insensitive) string matching the CFITSIO compression
1616             type names. Currently supported compression algorithms are:
1617              
1618             =over 3
1619              
1620             =item * RICE_1 - linear Rice compression
1621              
1622             This uses limited-symbol-length Rice compression, which works well on
1623             low entropy image data (where most pixels differ from their neighbors
1624             by much less than the dynamic range of the image).
1625              
1626             =back
1627              
1628             =item tilesize (default C<[-1,1]>)
1629              
1630             This specifies the dimension of the compression tiles, in pixels. You
1631             can hand in a PDL, a scalar, or an array ref. If you specify fewer
1632             dimensions than exist in the image, the last dim is repeated - so "32"
1633             yields 32x32 pixel tiles in a 2-D image. A dim of -1 in any dimension
1634             duplicates the image size, so the default C<[-1,1]> causes compression
1635             along individual rows.
1636              
1637             =item tilesize (RICE_1 only; default C<32>)
1638              
1639             For RICE_1, BLOCKSIZE indicates the number of pixel samples to use
1640             for each compression block within the compression algorithm. The
1641             blocksize is independent of the tile dimensions. For RICE
1642             compression the pixels from each tile are arranged in normal pixel
1643             order (early dims fastest) and compressed as a linear stream.
1644              
1645             =back
1646              
1647             =item * Table handling:
1648              
1649             If you feed in a hash ref instead of a PDL, then the hash ref is
1650             written out as a binary table extension. The hash ref keys are
1651             treated as column names, and their values are treated as the data to
1652             be put in each column.
1653              
1654             For numeric information, the hash values should contain PDLs. The 0th
1655             dim of the PDL runs across rows, and higher dims are written as
1656             multi-value entries in the table (e.g. a 7x5 PDL will yield a single
1657             named column with 7 rows and 5 numerical entries per row, in a binary
1658             table). Note that this is slightly different from the usual concept
1659             of threading, in which dimension 1 runs across rows.
1660              
1661             ASCII tables only allow one entry per column in each row, so
1662             if you plan to write an ASCII table then all of the values of C<$hash>
1663             should have at most one dim.
1664              
1665             All of the columns' 0 dims must agree in the threading sense. That is to
1666             say, the 0th dimension of all of the values of C<$hash> should be the
1667             same (indicating that all columns have the same number of rows). As
1668             an exception, if the 0th dim of any of the values is 1, or if that
1669             value is a PDL scalar (with 0 dims), then that value is "threaded"
1670             over -- copied into all rows.
1671              
1672             Data dimensions higher than 2 are preserved in binary tables,
1673             via the TDIMn field (e.g. a 7x5x3 PDL is stored internally as
1674             seven rows with 15 numerical entries per row, and reconstituted
1675             as a 7x5x3 PDL on read).
1676              
1677             Non-PDL Perl scalars are treated as strings, even if they contain
1678             numerical values. For example, a list ref containing 7 values is
1679             treated as 7 rows containing one string each. There is no such thing
1680             as a multi-string column in FITS tables, so any nonscalar values in
1681             the list are stringified before being written. For example, if you
1682             pass in a perl list of 7 PDLs, each PDL will be stringified before
1683             being written, just as if you printed it to the screen. This is
1684             probably not what you want -- you should use L to connect
1685             the separate PDLs into a single one. (e.g. C<$x-Eglue(1,$y,$c)-Emv(1,0)>)
1686              
1687             The column names are case-insensitive, but by convention the keys of
1688             C<$hash> should normally be ALL CAPS, containing only digits, capital
1689             letters, hyphens, and underscores. If you include other characters,
1690             then case is smashed to ALL CAPS, whitespace is converted to
1691             underscores, and unrecognized characters are ignored -- so if you
1692             include the key "Au Purity (%)", it will be written to the file as a
1693             column that is named "AU_PURITY". Since this is not guaranteed to
1694             produce unique column names, subsequent columns by the same name are
1695             disambiguated by the addition of numbers.
1696              
1697             You can specify the use of variable-length rows in the output, saving
1698             space in the file. To specify variable length rows for a column named
1699             "FOO", you can include a separate key "len_FOO" in the hash to be
1700             written. The key's value should be a PDL containing the number of
1701             actual samples in each row. The result is a FITS P-type variable
1702             length column that, upon read with C, will restore to a field
1703             named FOO and a corresponding field named "len_FOO". Invalid data in
1704             the final PDL consist of a padding value (which defaults to 0 but
1705             which you may set by including a TNULL field in the hdr specificaion).
1706             Variable length arrays must be 2-D PDLs, with the variable length in
1707             the 1 dimension.
1708              
1709             Two further special keys, 'hdr' and 'tbl', can contain
1710             meta-information about the type of table you want to write. You may
1711             override them by including an C<$OPTIONS> hash with a 'hdr' and/or
1712             'tbl' key.
1713              
1714             The 'tbl' key, if it exists, must contain either 'ASCII' or 'binary'
1715             (case-insensitive), indicating whether to write an ascii or binary
1716             table. The default is binary. [ASCII table writing is planned but
1717             does not yet exist].
1718              
1719             You can specify the format of the table quite specifically with the
1720             'hdr' key or option field. If it exists, then the 'hdr' key should
1721             contain fields appropriate to the table extension being used. Any
1722             field information that you don't specify will be filled in
1723             automatically, so (for example) you can specify that a particular
1724             column name goes in a particular position, but allow C to
1725             arrange the other columns in the usual alphabetical order into any
1726             unused slots that you leave behind. The C, C,
1727             C, C, C, and C keywords are ignored:
1728             their values are calculated based on the hash that you supply. Any
1729             other fields are passed into the final FITS header verbatim.
1730              
1731             As an example, the following
1732              
1733             $x = long(1,2,4);
1734             $y = double(1,2,4);
1735             wfits { 'COLA'=>$x, 'COLB'=>$y }, "table1.fits";
1736              
1737             will create a binary FITS table called F which
1738             contains two columns called C and C. The order
1739             of the columns is controlled by setting the C
1740             keywords in the header array, so
1741              
1742             $h = { 'TTYPE1'=>'Y', 'TTYPE2'=>'X' };
1743             wfits { 'X'=>$x, 'Y'=>$y, hdr=>$h }, "table2.fits";
1744              
1745             creates F where the first column is
1746             called C and the second column is C.
1747              
1748             =item * multi-value handling
1749              
1750             If you feed in a perl list rather than a PDL or a hash, then
1751             each element is written out as a separate HDU in the FITS file.
1752             Each element of the list must be a PDL or a hash. [This is not implemented
1753             yet but should be soon!]
1754              
1755             =item * DEVEL NOTES
1756              
1757             ASCII tables are not yet handled but should be.
1758              
1759             Binary tables currently only handle one vector (up to 1-D array)
1760             per table entry; the standard allows more, and should be fully implemented.
1761             This means that PDL::Complex piddles currently can not be written to
1762             disk.
1763              
1764             Handling multidim arrays implies that perl multidim lists should also be
1765             handled.
1766              
1767             =back
1768              
1769             =for bad
1770              
1771             For integer types (ie C 0>), the C keyword is set
1772             to the bad value. For floating-point types, the bad value is
1773             converted to NaN (if necessary) before writing.
1774              
1775             =cut
1776              
1777             *wfits = \&PDL::wfits;
1778              
1779             BEGIN {
1780 49     49   458 @PDL::IO::FITS::wfits_keyword_order =
1781             ('SIMPLE','BITPIX','NAXIS','NAXIS1','BUNIT','BSCALE','BZERO');
1782 49         334368 @PDL::IO::FITS::wfits_numbered_keywords =
1783             ('CTYPE','CRPIX','CRVAL','CDELT','CROTA');
1784             }
1785              
1786             # Until we do a rewrite these have to be file global since they
1787             # are used by the wheader routine
1788             my (%hdr, $nbytes);
1789              
1790             # Local utility routine of wfits()
1791             sub wheader ($$) {
1792 176     176 0 271 my $fh = shift;
1793 176         292 my $k = shift;
1794            
1795 176 100       553 if ($k =~ m/(HISTORY|COMMENT)/) {
1796 31         105 my $hc = $1;
1797 31 50       122 return unless ref($hdr{$k}) eq 'ARRAY';
1798 0         0 foreach my $line (@{$hdr{$k}}) {
  0         0  
1799 0         0 $fh->printf( "$hc %-72s", substr($line,0,72) );
1800 0         0 $nbytes += 80;
1801             }
1802 0         0 delete $hdr{$k};
1803             } else {
1804             # Check that we are dealing with a scalar value in the header
1805             # Need to make sure that the header does not include PDLs or
1806             # other structures. Return unless $hdr{$k} is a scalar.
1807 145         274 my($hdrk) = $hdr{$k};
1808            
1809 145 50       342 if(ref $hdrk eq 'ARRAY') {
1810 0         0 $hdrk = join("\n",@$hdrk);
1811             }
1812            
1813 145 50       269 return unless not ref($hdrk);
1814            
1815 145 50       288 if ($hdrk eq "") {
1816 0         0 $fh->printf( "%-80s", substr($k,0,8) );
1817             } else {
1818 145         525 $fh->printf( "%-8s= ", substr($k,0,8) );
1819            
1820             my $com = ( ref $hdr{COMMENT} eq 'HASH' ) ?
1821 145 50       1324 $hdr{COMMENT}{$k} : undef;
1822            
1823 145 100 33     927 if ($hdrk =~ /^ *([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))? *$/) { # Number?
    50          
1824 114 50       270 my $cl=60-($com ? 2 : 0);
1825 114         237 my $end=' ' x $cl;
1826 114 50       257 $end =' /'. $com if($com);
1827 114         357 $fh->printf( "%20s%-50s", substr($hdrk,0,20),
1828             substr($end, 0, 50) );
1829             } elsif ($hdrk eq 'F' or $hdrk eq 'T') {
1830             # Logical flags ?
1831 0         0 $fh->printf( "%20s", $hdrk );
1832 0         0 my $end=' ' x 50;
1833 0 0       0 $end =' /'.$com if($com);
1834 0         0 $fh->printf( "%-50s", $end );
1835             } else {
1836             # Handle strings, truncating as necessary
1837             # (doesn't do multicard strings like Astro::FITS::Header does)
1838            
1839             # Convert single quotes to doubled single quotes
1840             # (per FITS standard)
1841 31         100 my($st) = $hdrk;
1842 31         113 $st =~ s/\'/\'\'/g;
1843            
1844 31         85 my $sl=length($st)+2;
1845 31 50       101 my $cl=70-$sl-($com ? 2 : 0);
1846 31         205 $fh->print( "'$st'" );
1847            
1848 31 50       286 if (defined $com) {
1849 0         0 $fh->printf( " /%-$ {cl}s", substr($com, 0, $cl) );
1850             } else {
1851 31         153 $fh->printf( "%-$ {cl}s", ' ' x $cl );
1852             }
1853             }
1854             }
1855 145         1180 $nbytes += 80; delete $hdr{$k};
  145         277  
1856             }
1857 145 50       325 delete $hdr{COMMENT}{$k} if ref $hdr{COMMENT} eq 'HASH';
1858 145         282 1;
1859             }
1860              
1861             # Write a PDL to a FITS format file
1862             #
1863             sub PDL::wfits {
1864 38 50 33 38 0 12674 barf 'Usage: wfits($pdl,$file,[$BITPIX],[{options}])' if $#_<1 || $#_>3;
1865 38         142 my ($pdl,$file,$x,$y) = @_;
1866 38         123 my ($opt, $BITPIX);
1867              
1868 38         196 local $\ = undef; # fix sf.net bug #3394327
1869              
1870 38 50 33     322 if(ref $x eq 'HASH') {
    50          
1871 0         0 $opt = $x;
1872 0         0 $BITPIX = $y;
1873             } elsif(ref $y eq 'HASH' || ! defined $y) {
1874 38         79 $BITPIX = $x;
1875 38         78 $opt = $y;
1876             }
1877              
1878 38         84 my ($k, $buff, $off, $ndims, $sz);
1879              
1880 38         114 $file =~ s/^(~)/glob($1)/e; # tilde expansion
  0         0  
1881              
1882 38         880 local $SIG{PIPE};
1883              
1884 38 50       280 if ($file =~ /\.gz$/) { # Handle suffix-style compression
    50          
1885 0     0   0 $SIG{PIPE}= sub {}; # Prevent crashing if gzip dies
1886 0         0 $file = "|gzip -9 > $file";
1887             }
1888             elsif ($file =~ /\.Z$/) {
1889 0     0   0 $SIG{PIPE}= sub {}; # Prevent crashing if compress dies
1890 0         0 $file = "|compress > $file";
1891             }
1892             else{
1893 38         144 $file = ">$file";
1894             }
1895            
1896             #### Figure output type
1897              
1898 38         106 my @outputs = ();
1899 38         84 my $issue_nullhdu;
1900              
1901 38 100       235 if( UNIVERSAL::isa($pdl,'PDL') ) {
    100          
    100          
    50          
1902 31         77 $issue_nullhdu = 0;
1903 31         82 push(@outputs, $pdl);
1904             } elsif( ref($pdl) eq 'HASH' ) {
1905 3         8 $issue_nullhdu = 1;
1906 3         6 push(@outputs, $pdl);
1907             } elsif( ref($pdl) eq 'ARRAY' ) {
1908 1         4 $issue_nullhdu = 1;
1909 1         4 @outputs = @$pdl;
1910             } elsif( length(ref($pdl))==0 ) {
1911 3         12 barf "wfits: needs a HASH or PDL argument to write out\n";
1912             } else {
1913 0         0 barf "wfits: unknown ref type ".ref($pdl)."\n";
1914             }
1915              
1916             ## Open file & prepare to write binary info
1917 35 50       375 my $fh = IO::File->new( $file )
1918             or barf "Unable to create FITS file $file\n";
1919 35         7934 binmode $fh;
1920              
1921 35 100       128 if($issue_nullhdu) {
1922 4         25 _wfits_nullhdu($fh);
1923             }
1924              
1925 35         227 for $pdl(@outputs) {
1926              
1927 36 100       340 if(ref($pdl) eq 'HASH') {
    50          
1928             my $table_type = ( exists($pdl->{tbl}) ?
1929 3 0       16 ($pdl->{tbl} =~ m/^a/i ? 'ascii' : 'binary') :
    50          
1930             "binary"
1931             );
1932 3         16 _wfits_table($fh,$pdl,$table_type);
1933             } elsif( UNIVERSAL::isa($pdl,'PDL') ) {
1934              
1935             ### Regular image writing.
1936            
1937 33 100       158 $BITPIX = "" unless defined $BITPIX;
1938 33 100       124 if ($BITPIX eq "") {
1939 21 100       164 $BITPIX = 8 if $pdl->get_datatype == $PDL_B;
1940 21 100 66     208 $BITPIX = 16 if $pdl->get_datatype == $PDL_S || $pdl->get_datatype == $PDL_US;
1941 21 100       100 $BITPIX = 32 if $pdl->get_datatype == $PDL_L;
1942 21 100       114 $BITPIX = 64 if $pdl->get_datatype == $PDL_LL;
1943 21 100       103 $BITPIX = -32 if $pdl->get_datatype == $PDL_F;
1944 21 100       99 $BITPIX = -64 if $pdl->get_datatype == $PDL_D;
1945 21 50       112 $BITPIX = 8 * PDL::Core::howbig($PDL_IND) if($pdl->get_datatype==$PDL_IND);
1946             }
1947 33 50       120 if ($BITPIX eq "") {
1948 0         0 $BITPIX = -64;
1949 0         0 warn "wfits: PDL has an unsupported datatype -- defaulting to 64-bit float.\n";
1950             }
1951              
1952 33     0   245 my $convert = sub { return $_[0] }; # Default - do nothing
  0         0  
1953 33 100   8   144 $convert = sub {byte($_[0])} if $BITPIX == 8;
  8         36  
1954 33 100   11   115 $convert = sub {short($_[0])} if $BITPIX == 16;
  11         50  
1955 33 100   10   122 $convert = sub {long($_[0])} if $BITPIX == 32;
  10         60  
1956 33 100   2   125 $convert = sub {longlong($_[0])} if $BITPIX == 64;
  2         9  
1957 33 100   8   115 $convert = sub {float($_[0])} if $BITPIX == -32;
  8         36  
1958 33 100   28   164 $convert = sub {double($_[0])} if $BITPIX == -64;
  28         126  
1959            
1960             # Automatically figure output scaling
1961            
1962 33         77 my $bzero = 0; my $bscale = 1;
  33         74  
1963 33 100       107 if ($BITPIX>0) {
1964 15         158 my $min = $pdl->min;
1965 15         85 my $max = $pdl->max;
1966 15 100       67 my ($dmin,$dmax) = (0, 2**8-1) if $BITPIX == 8;
1967 15 100       60 ($dmin,$dmax) = (-2**15, 2**15-1) if $BITPIX == 16;
1968 15 100       56 ($dmin,$dmax) = (-2**31, 2**31-1) if $BITPIX == 32;
1969 15 100       58 ($dmin,$dmax) = (-(pdl(longlong,1)<<63), (pdl(longlong,1)<<63)-1) if $BITPIX==64;
1970              
1971 15 50 33     143 if ($min<$dmin || $max>$dmax) {
1972 0         0 $bzero = $min - $dmin;
1973 0         0 $max -= $bzero;
1974 0 0       0 $bscale = $max/$dmax if $max>$dmax;
1975             }
1976 15 50       64 print "BSCALE = $bscale && BZERO = $bzero\n" if $PDL::verbose;
1977             }
1978            
1979             # Check for tile-compression format for the image, and handle it.
1980             # We add the image-compression format tags and reprocess the whole
1981             # shebang as a binary table.
1982 33 50       122 if($opt->{compress}) {
1983 0         0 croak "Placeholder -- tile compression not yet supported\n";
1984             }
1985            
1986            
1987             ##############################
1988             ## Check header and prepare to write it out
1989            
1990 33         219 my($h) = $pdl->gethdr();
1991            
1992             # Extra logic: if we got handed a vanilla hash that that is *not* an Astro::FITS::Header, but
1993             # looks like it's a FITS header encoded in a hash, then attempt to process it with
1994             # Astro::FITS::Header before writing it out -- this helps with cleanup of tags.
1995 33 50 66     299 if($PDL::Astro_FITS_Header and
      66        
      66        
1996             defined($h) and
1997             ref($h) eq 'HASH' and
1998             !defined( tied %$h )
1999             ) {
2000            
2001 2         5 my $all_valid_fits = 1;
2002 2         10 for my $k(keys %$h) {
2003 5 50 33     41 if(length($k) > 8 or
2004             $k !~ m/^[A-Z_][A-Z\d\_]*$/i
2005             ) {
2006 0         0 $all_valid_fits = 0;
2007 0         0 last;
2008             }
2009             }
2010            
2011 2 50       10 if($all_valid_fits) {
2012             # All the keys look like valid FITS header keywords -- so
2013             # create a tied FITS header object and use that instead.
2014 2         30 my $afh = new Astro::FITS::Header( );
2015 2         1230 my %hh;
2016 2         15 tie %hh, "Astro::FITS::Header", $afh;
2017 2         33 for (keys %$h) {
2018 5         930 $hh{$_} = $h->{$_};
2019             }
2020 2         587 $h = \%hh;
2021             }
2022             }
2023            
2024             # Now decide whether to emit a hash or an AFH object
2025 33 100 33     200 if(defined($h) &&
      66        
2026             ( (defined (tied %$h)) &&
2027             (UNIVERSAL::isa(tied %$h,"Astro::FITS::Header")))
2028             ){
2029 2         17 my $k;
2030            
2031             ##n############################
2032             ## Tied-hash code -- I'm too lazy to incorporate this into KGB's
2033             ## direct hash handler below, so I've more or less just copied and
2034             ## pasted with some translation. --CED
2035             ##
2036 2         30 my $hdr = tied %$h;
2037            
2038             #
2039             # Put advertising comment in the SIMPLE field
2040             #n
2041 2 50       8 if($issue_nullhdu) {
2042 0         0 $h->{XTENSION} = "IMAGE";
2043             } else {
2044 2         15 $h->{SIMPLE} = 'T';
2045 2         104 my(@a) = $hdr->itembyname('SIMPLE');
2046 2         48 $a[0]->comment('Created with PDL (http://pdl.perl.org)');
2047             # and register it as a LOGICAL rather than a string
2048 2         21 $a[0]->type('LOGICAL');
2049             }
2050            
2051             #
2052             # Use tied interface to set all the keywords. Note that this
2053             # preserves existing per-line comments, only changing the values.
2054             #
2055 2         27 $h->{BITPIX} = $BITPIX;
2056 2         554 $h->{NAXIS} = $pdl->getndims;
2057 2         556 my $correction = 0;
2058 2         14 for $k(1..$h->{NAXIS}) {
2059             $correction |= (exists $h->{"NAXIS$k"} and
2060 4   33     848 $h->{"NAXIS$k"} != $pdl->dim($k-1)
2061             );
2062 4         94 $h->{"NAXIS$k"} = $pdl->getdim($k-1);
2063             }
2064 2 50       696 carp("Warning: wfits corrected dimensions of FITS header")
2065             if($correction);
2066            
2067 2 50       12 $h->{BUNIT} = "Data Value" unless exists $h->{BUNIT};
2068 2 50       744 $h->{BSCALE} = $bscale if($bscale != 1);
2069 2 50       9 $h->{BZERO} = $bzero if($bzero != 0);
2070            
2071 2 50       54 if ( $pdl->badflag() ) {
2072 0 0       0 if ( $BITPIX > 0 ) { my $x = &$convert(pdl(0.0));
  0         0  
2073 0         0 $h->{BLANK} = $x->badvalue(); }
2074 0         0 else { delete $h->{BLANK}; }
2075             }
2076            
2077             # Use object interface to sort the lines. This is complicated by
2078             # the need for an arbitrary number of NAXIS lines in the middle
2079             # of the sorting. Keywords with a trailing '1' in the sorted-order
2080             # list get looped over.
2081 2         8 my($kk) = 0;
2082 2         280 my(@removed_naxis) = ();
2083 2         14 for $k(0..$#PDL::IO::FITS::wfits_keyword_order) {
2084 14         30 my($kn) = 0;
2085            
2086 14         163 my @index;
2087 14   100     21 do { # Loop over numericised keywords (e.g. NAXIS1)
2088            
2089 18         36 my $kw = $PDL::IO::FITS::wfits_keyword_order[$k]; # $kw get keyword
2090 18 100       78 $kw .= (++$kn) if( $kw =~ s/\d$//); # NAXIS1 -> NAXIS
2091 18         49 @index = $hdr->index($kw);
2092            
2093 18 100       258 if(defined $index[0]) {
2094 12 50       49 if($kn <= $pdl->getndims){
2095 12 100       38 $hdr->insert($kk, $hdr->remove($index[0]))
2096             unless ($index[0] == $kk) ;
2097 12         2801 $kk++;
2098             }
2099             else{ #remove e.g. NAXIS3 from hdr if NAXIS==2
2100 0         0 $hdr->removebyname($kw);
2101 0         0 push(@removed_naxis,$kw);
2102             }
2103             }
2104             } while((defined $index[0]) && $kn);
2105             }
2106            
2107 2         7 foreach my $naxis(@removed_naxis){
2108 0         0 $naxis =~ m/(\d)$/;
2109 0         0 my $n = $1;
2110 0         0 foreach my $kw(@PDL::IO::FITS::wfits_numbered_keywords){
2111 0         0 $hdr->removebyname($kw . $n);
2112             }
2113             }
2114             #
2115             # Delete the END card if necessary (for later addition at the end)
2116             #
2117 2         12 $hdr->removebyname('END');
2118            
2119             #
2120             # Make sure that the HISTORY lines all come at the end
2121             #
2122 2         303 my @hindex = $hdr->index('HISTORY');
2123 2         36 for $k(0..$#hindex) {
2124 0         0 $hdr->insert(-1-$k, $hdr->remove($hindex[-1-$k]));
2125             }
2126            
2127             #
2128             # Make sure the last card is an END
2129             #
2130 2         10 $hdr->insert(scalar($hdr->cards),
2131             Astro::FITS::Header::Item->new(Keyword=>'END'));
2132            
2133             #
2134             # Write out all the cards, and note how many bytes for later padding.
2135             #
2136 2         1744 my $s = join("",$hdr->cards);
2137            
2138 2         308 $fh->print( $s );
2139 2         76 $nbytes = length $s;
2140             } else {
2141             ##
2142             ## Legacy emitter (note different advertisement in the SIMPLE
2143             ## comment, for debugging!)
2144             ##
2145              
2146 31 100       107 if($issue_nullhdu) {
2147 2         12 $fh->printf( "%-80s", "XTENSION= 'IMAGE'" );
2148             } else {
2149 29         235 $fh->printf( "%-80s", "SIMPLE = T / PDL::IO::FITS::wfits (http://pdl.perl.org)" );
2150             }
2151            
2152 31         738 $nbytes = 80; # Number of bytes written so far
2153            
2154             # Write FITS header
2155            
2156 31         101 %hdr = ();
2157 31 50       116 if (defined($h)) {
2158 0         0 for (keys %$h) { $hdr{uc $_} = $$h{$_} } # Copy (ensuring keynames are uppercase)
  0         0  
2159             }
2160            
2161 31         81 delete $hdr{SIMPLE}; delete $hdr{'END'};
  31         85  
2162            
2163 31         125 $hdr{BITPIX} = $BITPIX;
2164 31 50       136 $hdr{BUNIT} = "Data Value" unless exists $hdr{BUNIT};
2165 31         146 wheader($fh, 'BITPIX');
2166            
2167 31         159 $ndims = $pdl->getndims; # Dimensions of data array
2168 31         93 $hdr{NAXIS} = $ndims;
2169 31         93 wheader($fh, 'NAXIS');
2170 31         146 for $k (1..$ndims) { $hdr{"NAXIS$k"} = $pdl->getdim($k-1) }
  51         329  
2171 31         103 for $k (1..$ndims) { wheader($fh,"NAXIS$k") }
  51         149  
2172            
2173 31 50 33     213 if ($bscale != 1 || $bzero != 0) {
2174 0         0 $hdr{BSCALE} = $bscale;
2175 0         0 $hdr{BZERO} = $bzero;
2176 0         0 wheader($fh,'BSCALE');
2177 0         0 wheader($fh,'BZERO');
2178             }
2179 31         104 wheader($fh,'BUNIT');
2180            
2181             # IF badflag is set
2182             # and BITPIX > 0 - ensure the header contains the BLANK keyword
2183             # (make sure it's for the correct type)
2184             # otherwise - make sure the BLANK keyword is removed
2185 31 100       262 if ( $pdl->badflag() ) {
2186 2 100       8 if ( $BITPIX > 0 ) { my $x = &$convert(pdl(0.0)); $hdr{BLANK} = $x->badvalue(); }
  1         4  
  1         8  
2187 1         3 else { delete $hdr{BLANK}; }
2188             }
2189            
2190 31         229 for $k (sort fits_field_cmp keys %hdr) {
2191 1 50       6 wheader($fh,$k) unless $k =~ m/HISTORY/;
2192             }
2193 31         181 wheader($fh, 'HISTORY'); # Make sure that HISTORY entries come last.
2194 31         130 $fh->printf( "%-80s", "END" );
2195 31         235 $nbytes += 80;
2196             }
2197            
2198             #
2199             # Pad the header to a legal value and write the rest of the FITS file.
2200             #
2201 33         127 $nbytes %= 2880;
2202 33 50       406 $fh->print( " "x(2880-$nbytes) )
2203             if $nbytes != 0; # Fill up HDU
2204            
2205             # Decide how to byte swap - note does not quite work yet. Needs hack
2206             # to IO.xs
2207            
2208 33     4   337 my $bswap = sub {}; # Null routine
2209 33 50       185 if ( !isbigendian() ) { # Need to set a byte swap routine
2210 33 100       116 $bswap = \&bswap2 if $BITPIX==16;
2211 33 100 100     187 $bswap = \&bswap4 if $BITPIX==32 || $BITPIX==-32;
2212 33 100 100     188 $bswap = \&bswap8 if $BITPIX==-64 || $BITPIX==64;
2213             }
2214            
2215             # Write FITS data
2216            
2217 33         182 my $p1d = $pdl->copy->reshape($pdl->nelem); # Data as 1D stream;
2218            
2219 33         104 $off = 0;
2220 33         233 $sz = PDL::Core::howbig(&$convert($p1d->slice('0:0'))->get_datatype);
2221            
2222 33         332 $nbytes = $p1d->getdim(0) * $sz;
2223            
2224             # Transfer data in blocks (because might need to byte swap)
2225             # Buffer is also type converted on the fly
2226            
2227 33         70 my $BUFFSZ = 360*2880; # = ~1Mb - must be multiple of 2880
2228 33         75 my $tmp;
2229            
2230 33 100 100     204 if ( $pdl->badflag() and $BITPIX < 0 and $PDL::Bad::UseNaN == 0 ) {
      66        
2231             # just print up a message - conversion is actually done in the loop
2232 1 50       5 print "Converting PDL bad value to NaN\n" if $PDL::verbose;
2233             }
2234            
2235 33         146 while ($nbytes - $off > $BUFFSZ) {
2236            
2237             # Data to be transferred
2238            
2239 0         0 $buff = &$convert( ($p1d->slice( ($off/$sz).":". (($off+$BUFFSZ)/$sz-1))
2240             -$bzero)/$bscale );
2241            
2242             # if there are bad values present, and output type is floating-point,
2243             # convert the bad values to NaN's. We can ignore integer types, since
2244             # we have set the BLANK keyword
2245             #
2246 0 0 0     0 if ( $pdl->badflag() and $BITPIX < 0 and $PDL::Bad::UseNaN == 0 ) {
      0        
2247 0         0 $buff->inplace->setbadtonan();
2248             }
2249            
2250 0         0 &$bswap($buff);
2251 0         0 $fh->print( ${$buff->get_dataref} );
  0         0  
2252 0         0 $off += $BUFFSZ;
2253             }
2254 33         202 $buff = &$convert( ($p1d->slice($off/$sz.":-1") - $bzero)/$bscale );
2255            
2256 33 100 100     532 if ( $pdl->badflag() and $BITPIX < 0 and $PDL::Bad::UseNaN == 0 ) {
      66        
2257 1         4 $buff->inplace->setbadtonan();
2258             }
2259            
2260 33         520 &$bswap($buff);
2261 33         78 $fh->print( ${$buff->get_dataref} );
  33         243  
2262             # Fill HDU and close
2263             # note that for the data space the fill character is \0 not " "
2264             #
2265 33         737 $fh->print( "\0"x(($BUFFSZ - $buff->getdim(0) * $sz)%2880) );
2266             } # end of image writing block
2267             else {
2268             # Not a PDL and not a hash ref
2269 0         0 barf("wfits: unknown data type - quitting");
2270             }
2271             } # end of output loop
2272 35         1351 $fh->close();
2273 35         5250 1;
2274             }
2275              
2276             ######################################################################
2277             ######################################################################
2278              
2279             # Compare FITS headers in a sensible manner.
2280              
2281             =head2 fits_field_cmp
2282              
2283             =for ref
2284              
2285             fits_field_cmp
2286              
2287             Sorting comparison routine that makes proper sense of the digits at the end
2288             of some FITS header fields. Sort your hash keys using "fits_field_cmp" and
2289             you will get (e.g.) your "TTYPE" fields in the correct order even if there
2290             are 140 of them.
2291              
2292             This is a standard kludgey perl comparison sub -- it uses the magical
2293             $a and $b variables, rather than normal argument passing.
2294              
2295             =cut
2296              
2297             sub fits_field_cmp {
2298 126 100   126 1 385 if( $a=~m/^(.*[^\d])(\d+)$/ ) {
2299 64         137 my ($an,$ad) = ($1,$2);
2300 64 100       183 if( $b=~m/^(.*[^\d])(\d+)$/ ) {
2301 42         76 my($bn,$bd) = ($1,$2);
2302            
2303 42 100       87 if($an eq $bn) {
2304 23         57 return $ad<=>$bd;
2305             }
2306             }
2307             }
2308 103         223 $a cmp $b;
2309             }
2310              
2311             =head2 _rows()
2312              
2313             =for ref
2314              
2315             Return the number of rows in a variable for table entry
2316              
2317             You feed in a PDL or a list ref, and you get back the 0th dimension.
2318              
2319             =cut
2320              
2321             sub _rows {
2322 9     9   74 my $var = shift;
2323              
2324 9 100       63 return $var->dim(0) if( UNIVERSAL::isa($var,'PDL') );
2325 1 50       9 return 1+$#$var if(ref $var eq 'ARRAY');
2326 0 0       0 return 1 unless(ref $var);
2327            
2328 0 0       0 print STDERR "Warning: _rows found an unacceptable ref. ".ref $var.". Ignoring...\n"
2329             if($PDL::verbose);
2330            
2331 0         0 return undef;
2332             }
2333              
2334              
2335              
2336             =head2 _prep_table()
2337              
2338             =for ref
2339              
2340             Accept a hash ref containing a table, and return a header describing the table
2341             and a string to be written out as the table, or barf.
2342              
2343             You can indicate whether the table should be binary or ascii. The default
2344             is binary; it can be overridden by the "tbl" field of the hash (if present)
2345             or by parameter.
2346              
2347             =cut
2348              
2349             our %bintable_types = (
2350             'byte'=>['B',1],
2351             'short'=>['I',2],
2352             'ushort'=>['J',4, sub {return long shift;}],
2353             'long'=>['J',4],
2354             'longlong'=>['D', 8, sub {return double shift;}],
2355             'float'=>['E',4],
2356             'double'=>['D',8],
2357             # 'complex'=>['M',8] # Complex doubles are supported (actually, they aren't at the moment)
2358             );
2359              
2360              
2361             sub _prep_table {
2362 3     3   10 my ($hash,$tbl,$nosquish) = @_;
2363            
2364 3         6 my $ohash;
2365              
2366 3         8 my $hdr = $hash->{hdr};
2367            
2368 3         9 my $heap = "";
2369              
2370             # Make a local copy of the header.
2371 3         7 my $h = {};
2372 3 100       11 if(defined $hdr) {
2373 1         5 local $_;
2374 1         7 for (keys %$hdr) {$h->{$_} = $hdr->{$_}};
  2         6  
2375             }
2376 3         9 $hdr = $h;
2377              
2378 3 50       11 $tbl = $hash->{tbl} unless defined($tbl);
2379              
2380 3 50       10 barf "_prep_table called without a HASH reference as the first argument"
2381             unless ref $hash eq 'HASH';
2382              
2383             #####
2384             # Figure out how many columns are in the table
2385 3   66     44 my @colkeys = grep( ( !m/^(hdr|tbl)$/ and !m/^len_/ and defined $hash->{$_}),
2386             sort fits_field_cmp keys %$hash
2387             );
2388 3         9 my $cols = @colkeys;
2389              
2390 3 50       12 print "Table seems to have $cols columns...\n"
2391             if($PDL::verbose);
2392            
2393             #####
2394             # Figure out how many rows are in the table, and store counts...
2395             #
2396 3         10 my $rows;
2397             my $rkey;
2398 3         9 for my $key(@colkeys) {
2399 9         25 my $r = _rows($hash->{$key});
2400 9 100 66     40 ($rows,$rkey) = ($r,$key) unless(defined($rows) && $rows != 1);
2401 9 50 33     28 if($r != $rows && $r != 1) {
2402 0         0 barf "_prep_table: inconsistent number of rows ($rkey: $rows vs. $key: $r)\n";
2403             }
2404             }
2405            
2406 3 50       11 print "Table seems to have $rows rows...\n"
2407             if($PDL::verbose);
2408              
2409             #####
2410             # Squish and disambiguate column names
2411             #
2412 3         5 my %keysbyname;
2413             my %namesbykey;
2414              
2415 3 50       10 print "Renaming hash keys...\n"
2416             if($PDL::debug);
2417              
2418 3         10 for my $key(@colkeys) {
2419 9         20 my $name = $key;
2420              
2421 9         24 $name =~ tr/[a-z]/[A-Z]/; # Uppercaseify (required by FITS standard)
2422 9         22 $name =~ s/\s+/_/g; # Remove whitespace (required by FITS standard)
2423            
2424 9 50       43 unless($nosquish) {
2425 9         20 $name =~ s/[^A-Z0-9_-]//g; # Squish (recommended by FITS standard)
2426             }
2427            
2428             ### Disambiguate...
2429 9 50       23 if(defined $ohash->{$name}) {
2430 0         0 my $iter = 1;
2431 0         0 my $name2;
2432 0         0 do { $name2 = $name."_".($iter++); }
2433 0         0 while(defined $ohash->{$name2});
2434 0         0 $name = $name2;
2435             }
2436              
2437 9         22 $ohash->{$name} = $hash->{$key};
2438 9         16 $keysbyname{$name} = $key;
2439 9         15 $namesbykey{$key} = $name;
2440              
2441 9 0 33     85 print "\tkey '$key'\t-->\tname '$name'\n"
      33        
2442             if($PDL::debug || (($name ne $key) and $PDL::verbose));
2443             }
2444              
2445              
2446             # The first element of colnames is ignored (since FITS starts the
2447             # count at 1)
2448             #
2449 3         8 my @colnames; # Names by number
2450             my %colnums; # Numbers by name
2451              
2452             ### Allocate any table columns that are already in the header...
2453 3         7 local $_;
2454 3         20 map { for my $x(1) { # [Shenanigans to make next work right]
  2         7  
2455 2 50       11 next unless m/^TTYPE(\d*)$/;
2456              
2457 2         6 my $num = $1;
2458            
2459 2 50 33     12 if($num > $cols || $num < 1) {
2460 0 0       0 print "Ignoring illegal column number $num ( should be in range 1..$cols )\n"
2461             if($PDL::verbose);
2462 0         0 delete $hdr->{$_};
2463 0         0 next;
2464             }
2465              
2466 2         4 my $key = $hdr->{$_};
2467              
2468 2         3 my $name;
2469 2 50       8 unless( $name = $namesbykey{$key}) { # assignment
2470 0         0 $name = $key;
2471 0 0       0 unless( $key = $keysbyname{$key}) {
2472 0 0       0 print "Ignoring column $num in existing header (unknown name '$key')\n"
2473             if($PDL::verbose);
2474 0         0 next;
2475             }
2476             }
2477              
2478 2         7 $colnames[$num] = $name;
2479 2         8 $colnums{$name} = $num;
2480             } } sort fits_field_cmp keys %$hdr;
2481              
2482             ### Allocate all other table columns in alphabetical order...
2483 3         9 my $i = 0;
2484 3         10 for my $k (@colkeys) {
2485 9         15 my $name = $namesbykey{$k};
2486              
2487 9 100       19 unless($colnums{$name}) {
2488 7         19 while($colnames[++$i]) { }
2489 7         11 $colnames[$i] = $name;
2490 7         24 $colnums{$name} = $i;
2491 2         6 } else { $i++; }
2492             }
2493            
2494 3 50 33     15 print "Assertion failed: i ($i) != colnums ($cols)\n"
2495             if($PDL::debug && $i != $cols);
2496              
2497             print "colnames: " .
2498 3 50       11 join(",", map { $colnames[$_]; } (1..$cols) ) ."\n"
  0         0  
2499             if($PDL::debug);
2500              
2501             ########
2502             # OK, now the columns are allocated -- spew out a header.
2503              
2504 3         8 my @converters = (); # Will fill up with conversion routines for each column
2505 3         7 my @field_len = (); # Will fill up with field lengths for each column
2506 3         7 my @internaltype = (); # Gets flag for PDLhood
2507 3         5 my @fieldvars = (); # Gets refs to all the fields of the hash.
2508              
2509 3 50       10 if($tbl eq 'binary') {
    0          
2510 3         9 $hdr->{XTENSION} = 'BINTABLE';
2511 3         8 $hdr->{BITPIX} = 8;
2512 3         8 $hdr->{NAXIS} = 2;
2513             #$hdr->{NAXIS1} = undef; # Calculated below; inserted here as placeholder.
2514 3         8 $hdr->{NAXIS2} = $rows;
2515 3         8 $hdr->{PCOUNT} = 0; # Change this is variable-arrays are adopted
2516 3         8 $hdr->{GCOUNT} = 1;
2517 3         9 $hdr->{TFIELDS} = $cols;
2518              
2519             # Figure out data types, and accumulate a row length at the same time.
2520              
2521 3         7 my $rowlen = 0;
2522              
2523             # NOTE:
2524             # the conversion from ushort to long below is a hack to work
2525             # around the issue that otherwise perl treats it as a 2-byte
2526             # NOT 4-byte string on writing out, which leads to data corruption
2527             # Really ushort arrays should be written out using SCALE/ZERO
2528             # so that it can be written as an Int2 rather than Int4
2529             #
2530 3         18 for my $i(1..$cols) {
2531 9         27 $fieldvars[$i] = $hash->{$keysbyname{$colnames[$i]}};
2532 9         14 my $var = $fieldvars[$i];
2533              
2534 9         31 $hdr->{"TTYPE$i"} = $colnames[$i];
2535 9         31 my $tform;
2536            
2537             my $tstr;
2538 9         0 my $rpt;
2539 9         0 my $bytes;
2540            
2541 9 100       80 if( UNIVERSAL::isa($var,'PDL') ) {
    50          
    0          
2542              
2543 8         17 $internaltype[$i] = 'P';
2544              
2545 8         13 my $t;
2546              
2547 8         32 my $dims = pdl($var->dims);
2548 8         45 ($t = $dims->slice("(0)")) .= 1;
2549 8         40 $rpt = $dims->prod;
2550              
2551             =pod
2552              
2553             =begin WHENCOMPLEXVALUESWORK
2554            
2555             if( UNIVERSAL::isa($var,'PDL::Complex') ) {
2556             $rpt = $var->dim(1);
2557             $t = 'complex'
2558             } else {
2559             $t = type $var;
2560             }
2561              
2562             =end WHENCOMPLEXVALUESWORK
2563              
2564             =cut
2565              
2566 8 50       57 barf "Error: wfits() currently can not handle PDL::Complex arrays (column $colnames[$i])\n"
2567             if UNIVERSAL::isa($var,'PDL::Complex');
2568 8         27 $t = $var->type;
2569              
2570 8         30 $t = $bintable_types{$t};
2571            
2572 8 50       28 unless(defined($t)) {
2573 0 0       0 print "Warning: converting unknown type $t (column $colnames[$i]) to double...\n"
2574             if($PDL::verbose);
2575 0         0 $t = $bintable_types{'double'};
2576             }
2577              
2578 8         45 ($tstr, $bytes, $converters[$i]) = @$t;
2579              
2580             } elsif( ref $var eq 'ARRAY' ) {
2581              
2582 1         3 $internaltype[$i] = 'A';
2583 1         4 $bytes = 1;
2584            
2585             # Got an array (of strings) -- find the longest element
2586 1         3 $rpt = 0;
2587 1         3 for(@$var) {
2588 4         6 my $l = length($_);
2589 4 100       11 $rpt = $l if($l>$rpt);
2590             }
2591 1         5 ($tstr, $bytes, $converters[$i]) = ('A',1,undef);
2592              
2593             } elsif( ref $var ) {
2594 0         0 barf "You seem to be writing out a ".(ref $var)." as a table column. I\ndon't know how to do that (yet).\n";
2595            
2596             } else { # Scalar
2597 0         0 $internaltype[$i] = 'A';
2598 0         0 ($tstr, $bytes, $converters[$i]) = ('A',1,undef);
2599 0         0 $rpt = length($var);
2600             }
2601              
2602              
2603             # Now check if it's a variable-length array and, if so, insert an
2604             # extra converter
2605 9         39 my $lname = "len_".$keysbyname{$colnames[$i]};
2606 9 50       33 if(exists $hash->{$lname}) {
2607 0         0 my $lengths = $hash->{$lname};
2608              
2609             # Variable length array - add extra handling logic.
2610              
2611             # First, check we're legit
2612 0 0 0     0 if( !UNIVERSAL::isa($var, 'PDL') ||
      0        
      0        
      0        
2613             $var->ndims != 2 ||
2614             !UNIVERSAL::isa($lengths,'PDL') ||
2615             $lengths->ndims != 1 ||
2616             $lengths->dim(0) != $var->dim(0)
2617             ) {
2618 0         0 die <<'FOO';
2619             wfits(): you specified a 'len_$keysbyname{$colnames[$i]}' field in
2620             your binary table output hash, indicating a variable-length array for
2621             each row of the output table, but I'm having trouble interpreting it.
2622             Either your source column isn't a 2-D PDL, or your length column isn't
2623             a 1-D PDL, or the two lengths don't match. I give up.
2624             FOO
2625             }
2626              
2627             # The definition below wraps around the existing converter,
2628             # dumping the variable to the heap and returning the length
2629             # and index of the data for the current row as a PDL LONG.
2630             # This does the Right Thing below in the write loop, with
2631             # the side effect of putting the data into the heap.
2632             #
2633             # The main downside here is that the heap gets copied
2634             # multiple times as we accumulate it, since we are using
2635             # string concatenation to add onto it. It might be better
2636             # to preallocate a large heap, but I'm too lazy to figure
2637             # out how to do that.
2638 0         0 my $csub = $converters[$i];
2639             $converters[$i] = sub {
2640 0     0   0 my $var = shift;
2641 0         0 my $row = shift;
2642 0         0 my $col = shift;
2643            
2644 0         0 my $len = $hash->{"len_".$keysbyname{$colnames[$i]}};
2645 0         0 my $l;
2646 0 0       0 if(ref $len eq 'ARRAY') {
    0          
    0          
2647 0         0 $l = $len->[$row];
2648             } elsif( UNIVERSAL::isa($len,'PDL') ) {
2649 0         0 $l = $len->dice_axis(0,$row);
2650             } elsif( ref $len ) {
2651 0         0 die "wfits: Couldn't understand length spec 'len_".$keysbyname{$colnames[$i]}."' in bintable output (length spec must be a PDL or array ref).\n";
2652             } else {
2653 0         0 $l = $len;
2654             }
2655            
2656             # The standard says we should give a zero-offset
2657             # pointer if the current row is zero-length; hence
2658             # the ternary operator.
2659 0 0       0 my $ret = pdl( $l, $l ? length($heap) : 0)->long;
2660              
2661              
2662 0 0       0 if($l) {
2663             # This echoes the normal-table swap and accumulation
2664             # stuff below, except we're accumulating into the heap.
2665 0 0       0 my $tmp = $csub ? &$csub($var, $row, $col) : $var;
2666 0         0 $tmp = $tmp->slice("0:".($l-1))->sever;
2667            
2668 0 0       0 if(!isbigendian()) {
2669 0 0       0 bswap2($tmp) if($tmp->get_datatype == $PDL_S);
2670 0 0 0     0 bswap4($tmp) if($tmp->get_datatype == $PDL_L ||
2671             $tmp->get_datatype == $PDL_F);
2672 0 0       0 bswap8($tmp) if($tmp->get_datatype == $PDL_D);
2673             }
2674 0         0 my $t = $tmp->get_dataref;
2675 0         0 $heap .= $$t;
2676             }
2677              
2678 0         0 return $ret;
2679 0         0 };
2680              
2681             # Having defined the conversion routine, now modify tstr to make this a heap-array
2682             # reference.
2683 0         0 $tstr = sprintf("P%s(%d)",$tstr, $hash->{"len_".$keysbyname{$colnames[$i]}}->max );
2684 0         0 $rpt = 1;
2685 0         0 $bytes = 8; # two longints per row in the main table.
2686             }
2687              
2688            
2689 9         75 $hdr->{"TFORM$i"} = "$rpt$tstr";
2690              
2691 9 50 66     59 if(UNIVERSAL::isa($var, 'PDL') and $var->ndims > 1) {
2692 0         0 $hdr->{"TDIM$i"} = "(".join(",",$var->slice("(0)")->dims).")";
2693             }
2694              
2695 9         33 $rowlen += ($field_len[$i] = $rpt * $bytes);
2696             }
2697            
2698 3         9 $hdr->{NAXIS1} = $rowlen;
2699            
2700             ## Now accumulate the binary table
2701              
2702 3         7 my $table = "";
2703            
2704 3         13 for my $r(0..$rows-1) {
2705 12         19 my $row = "";
2706 12         24 for my $c(1..$cols) {
2707 36         54 my $tmp;
2708 36         47 my $x = $fieldvars[$c];
2709            
2710 36 100       70 if($internaltype[$c] eq 'P') { # PDL handling
2711             $tmp = $converters[$c]
2712 32 100       109 ? &{$converters[$c]}($x->slice("$r")->flat->sever, $r, $c)
  4         10  
2713             : $x->slice("$r")->flat->sever ;
2714              
2715             ## This would go faster if moved outside the loop but I'm too
2716             ## lazy to do it Right just now. Perhaps after it actually works.
2717             ##
2718 32 50       116 if(!isbigendian()) {
2719 32 100       143 bswap2($tmp) if($tmp->get_datatype == $PDL_S);
2720 32 100 100     241 bswap4($tmp) if($tmp->get_datatype == $PDL_L ||
2721             $tmp->get_datatype == $PDL_F);
2722 32 100       150 bswap8($tmp) if($tmp->get_datatype == $PDL_D);
2723             }
2724              
2725 32         74 my $t = $tmp->get_dataref;
2726 32         71 $tmp = $$t;
2727             } else { # Only other case is ASCII just now...
2728 4 50       18 $tmp = ( ref $x eq 'ARRAY' ) ? # Switch on array or string
    50          
2729             ( $#$x == 0 ? $x->[0] : $x->[$r] ) # Thread arrays as needed
2730             : $x;
2731            
2732 4         14 $tmp .= " " x ($field_len[$c] - length($tmp));
2733             }
2734            
2735             # Now $tmp contains the bytes to be written out...
2736             #
2737 36         161 $row .= substr($tmp,0,$field_len[$c]);
2738             } # for: $c
2739 12         34 $table .= $row;
2740             } # for: $r
2741              
2742 3         9 my $table_size = $rowlen * $rows;
2743 3 50       9 if( (length $table) != $table_size ) {
2744 0         0 print "Warning: Table length is ".(length $table)."; expected $table_size\n";
2745             }
2746              
2747 3         35 return ($hdr,$table, $heap);
2748            
2749             } elsif($tbl eq 'ascii') {
2750 0         0 barf "ASCII tables not yet supported...\n";
2751             } else {
2752 0         0 barf "unknown table type '$tbl' -- giving up.";
2753             }
2754             }
2755              
2756             # the header fill value can be blanks, but the data fill value must
2757             # be zeroes in non-ASCII tables
2758             #
2759             sub _print_to_fits ($$$) {
2760 10     10   22 my $fh = shift;
2761 10         22 my $data = shift;
2762 10         22 my $blank = shift;
2763              
2764 10         40 my $len = ((length $data) - 1) % 2880 + 1;
2765 10         167 $fh->print( $data . ($blank x (2880-$len)) );
2766             }
2767            
2768             {
2769             my $ctr = 0;
2770              
2771 3     3 0 7 sub reset_hdr_ctr() { $ctr = 0; }
2772              
2773             sub add_hdr_item ($$$$;$) {
2774 48     48 0 125 my ( $hdr, $key, $value, $type, $comment ) = @_;
2775 48 100       104 $type = uc($type) if defined $type;
2776 48         124 my $item = Astro::FITS::Header::Item->new( Keyword=>$key,
2777             Value=>$value,
2778             Type=>$type );
2779 48 100       4839 $item->comment( $comment ) if defined $comment;
2780 48         197 $hdr->replace( $ctr++, $item );
2781             };
2782             }
2783              
2784             ##############################
2785             #
2786             # _wfits_table -- given a hash ref, try to write it out as a
2787             # table extension. The file FITS should be open when you call it.
2788             # Most of the work of creating the extension header, and all of
2789             # the work of creating the table, is handled by _prep_table().
2790             #
2791             # NOTE:
2792             # can not think of a sensible name for the extension so calling
2793             # it TABLE for now
2794             #
2795             sub _wfits_table ($$$) {
2796 3     3   7 my $fh = shift;
2797 3         7 my $hash = shift;
2798 3         6 my $tbl = shift;
2799            
2800 3 50       24 barf "FITS BINTABLES are not supported without the Astro::FITS::Header module.\nGet it from www.cpan.org.\n"
2801             unless($PDL::Astro_FITS_Header);
2802              
2803 3         15 my ($hdr,$table, $heap) = _prep_table($hash,$tbl,0);
2804 3 50       10 $heap="" unless defined($heap);
2805              
2806             # Copy the prepared fields into the extension header.
2807 3         42 tie my %newhdr,'Astro::FITS::Header',my $h = Astro::FITS::Header->new;
2808            
2809 3         1216 reset_hdr_ctr();
2810 3 50       20 add_hdr_item $h, "XTENSION", ($tbl eq 'ascii'?"TABLE":"BINTABLE"), 'string', "from perl hash";
2811 3         230 add_hdr_item $h, "BITPIX", $hdr->{BITPIX}, 'int';
2812 3         183 add_hdr_item $h, "NAXIS", 2, 'int';
2813 3         207 add_hdr_item $h, "NAXIS1", $hdr->{NAXIS1}, 'int', 'Bytes per row';
2814 3         339 add_hdr_item $h, "NAXIS2", $hdr->{NAXIS2}, 'int', 'Number of rows';
2815 3 50       285 add_hdr_item $h, "PCOUNT", length($heap), 'int', ($tbl eq 'ascii' ? undef : "No heap") ;
2816 3 50       307 add_hdr_item $h, "THEAP", "0", "(No gap before heap)" if(length($heap));
2817 3         12 add_hdr_item $h, "GCOUNT", 1, 'int';
2818 3         342 add_hdr_item $h, "TFIELDS", $hdr->{TFIELDS},'int';
2819 3         380 add_hdr_item $h, "HDUNAME", "TABLE", 'string';
2820              
2821 3         464 for my $field( sort fits_field_cmp keys %$hdr ) {
2822 42 100 66     5061 next if( defined $newhdr{$field} or $field =~ m/^end|simple|xtension$/i);
2823             my $type = (UNIVERSAL::isa($hdr->{field},'PDL') ?
2824             $hdr->{$field}->type :
2825 18 50       1094 ((($hdr->{$field})=~ m/^[tf]$/i) ?
    50          
2826             'logical' :
2827             undef ## 'string' seems to have a bug - 'undef' works OK
2828             ));
2829              
2830 18         76 add_hdr_item $h, $field, $hdr->{$field}, $type, $hdr->{"${field}_COMMENT"};
2831             }
2832              
2833 3         246 add_hdr_item $h, "END", undef, 'undef';
2834              
2835 3         669 $hdr = join("",$h->cards);
2836 3         3039 _print_to_fits( $fh, $hdr, " " );
2837 3         51 _print_to_fits( $fh, $table.$heap, "\0" ); # use " " if it is an ASCII table
2838             # Add heap dump
2839             }
2840              
2841             sub _wfits_nullhdu ($) {
2842 4     4   12 my $fh = shift;
2843 4 50       16 if($Astro::FITS::Header) {
2844 0         0 my $h = Astro::FITS::Header->new();
2845            
2846 0         0 reset_hdr_ctr();
2847 0         0 add_hdr_item $h, "SIMPLE", "T", 'logical', "Null HDU (no data, only extensions)";
2848 0         0 add_hdr_item $h, "BITPIX", -32, 'int', "Needed to make fverify happy";
2849 0         0 add_hdr_item $h, "NAXIS", 0, 'int';
2850 0         0 add_hdr_item $h, "EXTEND", "T", 'logical', "File contains extensions";
2851 0         0 add_hdr_item $h, "COMMENT", "", "comment",
2852             " File written by perl (PDL::IO::FITS::wfits)";
2853             #
2854             # The following seems to cause a problem so removing for now (I don't
2855             # believe it is required, but may be useful for people who aren't
2856             # FITS connoisseurs). It could also be down to a version issue in
2857             # Astro::FITS::Header since it worked on linux with a newer version
2858             # than on Solaris with an older version of the header module)
2859             #
2860             ## add_hdr_item $h, "COMMENT", "", "comment",
2861             ## " FITS (Flexible Image Transport System) format is defined in 'Astronomy";
2862             ## add_hdr_item $h, "COMMENT", "", "comment",
2863             ## " and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H";
2864 0         0 add_hdr_item $h, "HDUNAME", "PRIMARY", 'string';
2865 0         0 add_hdr_item $h, "END", undef, 'undef';
2866            
2867 0         0 my $hdr = join("",$h->cards);
2868 0         0 _print_to_fits( $fh, $hdr, " " );
2869             } else {
2870 4         21 _print_to_fits( $fh,
2871             q+SIMPLE = T / Null HDU (no data, only extensions) BITPIX = -32 / Needed to make fverify happy NAXIS = 0 EXTEND = T / File contains extensions COMMENT Written by perl (PDL::IO::FITS::wfits) legacy code. COMMENT For best results, install Astro::FITS::Header. HDUNAME = 'PRIMARY ' END +,
2872             " ");
2873             }
2874             }
2875              
2876            
2877             1;
2878