File Coverage

blib/lib/PDLA/IO/FITS.pm
Criterion Covered Total %
statement 324 969 33.4
branch 155 578 26.8
condition 63 233 27.0
subroutine 28 53 52.8
pod 3 10 30.0
total 573 1843 31.0


line stmt bran cond sub pod time code
1             =head1 NAME
2              
3             PDLA::IO::FITS -- Simple FITS support for PDLA
4              
5             =head1 SYNOPSIS
6              
7             use PDLA;
8             use PDLA::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 PDLA, 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 PDLA distribution. If this file is separated from the
44             PDLA distribution, the copyright notice should be pasted into in this file.
45              
46             =head1 FUNCTIONS
47              
48             =cut
49              
50 15     15   70634 use strict;
  15         34  
  15         467  
51 15     15   87 use warnings;
  15         49  
  15         1465  
52              
53             BEGIN {
54              
55             package PDLA::IO::FITS;
56              
57 15     15   60 $PDLA::IO::FITS::VERSION = 0.92; # Will be 1.0 when ascii table read/write works.
58              
59 15         45 our @EXPORT_OK = qw( rfits rfitshdr wfits );
60 15         52 our %EXPORT_TAGS = (Func=>[@EXPORT_OK]);
61 15         318 our @ISA = ('PDLA::Exporter');
62              
63 15     15   106 use PDLA::Core;
  15         43  
  15         116  
64 15     15   4874 use PDLA::Config;
  15         44  
  15         505  
65 15     15   1998 use PDLA::IO::Misc;
  15         50  
  15         1292  
66 15     15   1111 use PDLA::Exporter;
  15         75  
  15         103  
67 15     15   112 use PDLA::Primitive;
  15         36  
  15         85  
68 15     15   425 use PDLA::Types;
  15         231  
  15         1823  
69 15     15   105 use PDLA::Options;
  15         31  
  15         820  
70 15     15   106 use PDLA::Bad;
  15         51  
  15         93  
71             # use PDLA::NiceSlice;
72 15     15   103 use Carp;
  15         28  
  15         838  
73 15     15   94 use strict;
  15         52  
  15         2663  
74              
75             ##############################
76             #
77             # Check if there's Astro::FITS::Header support, and set flag.
78             # Kludgy but it only has to run once, on first load. --CED
79             #
80 15     15   1289 eval "use Astro::FITS::Header;";
  15         2589  
  0         0  
  0         0  
81 15         67 $PDLA::Astro_FITS_Header = (defined $Astro::FITS::Header::VERSION);
82 15 50       72 if($PDLA::Astro_FITS_Header) {
83 0         0 my($x) = $Astro::FITS::Header::VERSION;
84 0         0 $x =~ s/[^0-9\.].*//;
85 0 0       0 $PDLA::Astro_FITS_Header = 0 if($x < 1.12);
86             }
87              
88 15 50       52 unless($PDLA::Astro_FITS_Header) {
89 15 50 33     14372 unless($ENV{"PDLA_FITS_LEGACY"} || $PDLA::Config{FITS_LEGACY}) {
90 0         0 print(STDERR "\n\nWARNING: Can't find the Astro::FITS::Header module, limiting FITS support.\n\n PDLA 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 'PDLA_FITS_LEGACY' or the global PDLA config value (in perldl.conf)\n \$PDLA::Config{FITS_LEGACY} to 1.\n\n");
91             }
92             }
93             }
94              
95             package PDLA::IO::FITS;
96              
97             ## declare subroutines
98              
99             sub _wfits_nullhdu ($);
100             sub _wfits_table ($$$);
101              
102             =head2 rfits()
103              
104             =for ref
105              
106             Simple piddle FITS reader.
107              
108             =for example
109              
110             $pdl = rfits('file.fits'); # Read a simple FITS image
111              
112             Suffix magic:
113              
114             $pdl = rfits('file.fits.gz'); # Read a file with gunzip(1)
115             $pdl = rfits('file.fits.Z'); # Read a file with uncompress(1)
116              
117             $pdl = rfits('file.fits[2]'); # Read 2nd extension
118             $pdl = rfits('file.fits.gz[3]'); # Read 3rd extension
119             @pdls = rfits('file.fits'); # Read primary data and extensions
120              
121             Tilde expansion:
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 PDLA;
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             PDLA. Setting the flag will cause an explicit deep copy of the header whenever
159             you use the returned PDLA in an arithmetic or slicing operation. That is useful
160             in many circumstances but also causes a hit in speed. When two or more PDLAs
161             with hdrcpy set are used in an expression, the result gets the header of the
162             first PDLA 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 PDLA and can be retrieved
187             with L or L. The
188             L flag of the PDLA 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 PDLA 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 PDLA
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 PDLA (for numerical values) or a
236             perl list (for string values). The PDLA'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 PDLA 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-PDLA'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 PDLAs 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 PDLA::Options( { bscale=>1, data=>1, hdrcpy=>0, expand=>1, afh=>1 } );
322              
323             sub PDLA::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 PDLA::rfits($class,$file,$u_opt);
329             }
330              
331             sub PDLA::rfits {
332            
333 7     7 0 18 my $class = shift;
334            
335 7 50 33     48 barf 'Usage: $x = rfits($file) -or- $x = PDLA->rfits($file)' if (@_ < 1 || @_ > 2);
336            
337 7         23 my $file = shift;
338            
339 7         34 my $u_opt = ifhref(shift);
340 7         34 my $opt = $rfits_options->options($u_opt);
341            
342 7         18 my($nbytes, $line, $name, $rest, $size, $i, $bscale, $bzero, $extnum);
343              
344 7         12 $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 7 50       25 my $explicit_extension = ($file =~ m/\[\d+\]$/ ? 1 : 0);
349 7 50       25 $extnum = ( ($file =~ s/\[(\d+)\]$//) ? $1 : 0 );
350              
351 7         15 $file =~ s/^(~)/glob($1)/e; #tilde expansion.
  0         0  
352              
353 7 50       24 $file = "gunzip -c $file |" if $file =~ /\.gz$/; # Handle compression
354 7 50       20 $file = "uncompress -c $file |" if $file =~ /\.Z$/;
355            
356 7 50       40 my $fh = IO::File->new( $file )
357             or barf "FITS file $file not found";
358 7         560 binmode $fh;
359              
360 7         17 my @extensions; # This accumulates the list in list context...
361 7         12 my $currentext=0;
362 7         17 my $pdl;
363              
364 7   100     16 hdu:{do { # Runs over extensions, in list context
  7         12  
365 9         30 my $ext_type = 'IMAGE'; # Gets the type of XTENSION if one is detected.
366 9         17 my $foo={}; # To go in pdl
367 9         21 my @history=();
368 9         17 my @cards = ();
369            
370 9         49 $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 9 50 66     54 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 9         51 my $ct = $fh->read($line,80);
387 9 50 66     258 barf "file $file is not in FITS-format:\n$line\n"
388             if( $nbytes==0 && ($line !~ /^SIMPLE = +T/));
389 9 50 33     49 last hdu if($fh->eof() || !$ct);
390             }
391              
392 9         89 $nbytes = 80; # Number of bytes read from this extension (1 line so far)
393            
394 9 100       42 if($line =~ /^XTENSION= \'(\w+)\s*\'/) {
    50          
395 2         6 $ext_type = $1;
396             } elsif( @extensions ) {
397              
398 0 0       0 print "Warning: expected XTENSION, found '$line'. Exiting.\n"
399             if($PDLA::verbose);
400 0         0 last hdu;
401             }
402              
403 9 50       23 push(@cards,$line) if($PDLA::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 9 50 66     37 if(!wantarray and $currentext != $extnum) {
415              
416 0         0 skipper: while(1) {
417             # Move to next record
418 0         0 $nbytes += $fh->read($line,2880-80);
419 0 0       0 barf "Unexpected end of FITS file\n" if $fh->eof();
420             # Read start of next record
421 0         0 $nbytes += $fh->read($line,80);
422 0 0       0 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 0 0       0 $currentext++ if $line =~ /^XTENSION\= \'(\w+)\s*\'/;
427 0 0       0 if ($currentext == $extnum) {
428 0         0 $ext_type = $1;
429 0         0 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 9 50 33     42 if($PDLA::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 0   0     0 do {
445 0         0 $nbytes += $fh->read($line, 80);
446 0         0 push(@cards,$line);
447             } while(!$fh->eof() && $line !~ m/^END(\s|\000)/);
448              
449 0         0 $nbytes += $fh->read(my $dummy, 2879 - ($nbytes-1)%2880);
450              
451 0         0 my($hdr) = Astro::FITS::Header->new(Cards => \@cards);
452 0         0 my(%hdrhash);
453 0         0 tie %hdrhash,"Astro::FITS::Header",$hdr;
454 0         0 $foo = \%hdrhash;
455            
456             } else {
457            
458             ## Legacy (straight header-to-hash-ref) parsing.
459             ## Cheesy but fast.
460            
461 9         15 hdr_legacy: { do {
  9         12  
462 15     15   132 no strict 'refs';
  15         39  
  15         110561  
463             # skip if the first eight characters are ' '
464             # - as seen in headers from the DSS at STScI
465 64 50       578 if (substr($line,0,8) ne " " x 8) { # If non-blank
466            
467 64         152 $name = (split(' ',substr($line,0,8)))[0];
468              
469 64         119 $rest = substr($line,8);
470            
471 64 50       115 if ($name =~ m/^HISTORY/) {
472 0         0 push @history, $rest;
473             } else {
474 64         125 $$foo{$name} = "";
475            
476 64 100       274 $$foo{$name}=$1 if $rest =~ m|^= +([^\/\' ][^\/ ]*) *( +/(.*))?$| ;
477 64 100       165 $$foo{$name}=$1 if $rest =~ m|^= \'(.*)\' *( +/(.*))?$| ;
478 64 100       142 $$foo{COMMENT}{$name} = $3 if defined($3);
479             }
480             } # non-blank
481 64 100 66     201 last hdr_legacy if ((defined $name) && $name eq "END");
482 55         131 $nbytes += $fh->read($line, 80);
483             } while(!$fh->eof()); }
484              
485             # Clean up HISTORY card
486 9 50       38 $$foo{"HISTORY"} = \@history if $#history >= 0;
487            
488             # Step to end of header block in file
489 9         22 my $skip = 2879 - ($nbytes-1)%2880;
490 9 50       45 $fh->read(my $dummy, $skip) if $skip;
491 9         79 $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 9 50 100     58 if( !(defined $foo->{XTENSION}) # Primary header
      66        
      33        
503             and $foo->{NAXIS} == 0 # No data
504             and !wantarray # Scalar context
505             and !$explicit_extension # No HDU specifier
506             ) {
507 0 0       0 print "rfits: Skipping null primary HDU (use [0] to force read of primary)...\n"
508             if($PDLA::verbose);
509 0         0 return PDLA::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 9 100 66     53 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         8 $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 8 50       31 if (ref $PDLA::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 8         13 $pdl = &{$PDLA::IO::FITS::_Extension->{$ext_type}}($fh,$foo,$opt,$pdl);
  8         31  
534            
535             } else {
536              
537 0 0 0     0 print STDERR "rfits: Ignoring unknown extension '$ext_type'...\n"
538             if($PDLA::verbose || $PDLA::debug);
539            
540 0         0 $pdl = undef;
541              
542             }
543             }
544              
545             #
546             # Note -- $pdl isn't necessarily a PDLA. It's only a $pdl if
547             # the extension was an IMAGE.
548             #
549 9 100       222 push(@extensions,$pdl) if(wantarray);
550 9         41 $currentext++;
551            
552             } while( wantarray && !$fh->eof() );} # Repeat if we are in list context
553            
554 7         46 $fh->close;
555            
556 7 100       264 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         3 shift @extensions;
564             }
565             # Return all the extensions
566 1         12 return @extensions;
567             }
568 6         67 return $pdl;
569             }
570              
571              
572 7     7 1 846 sub rfits { PDLA->rfits(@_); }
573              
574             sub rfitshdr {
575 0     0 1 0 my($file,$opt) = shift;
576 0         0 $opt->{data} =0;
577 0         0 PDLA->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             $PDLA::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=>$PDLA_B,
600             16=>$PDLA_S,
601             32=>$PDLA_L,
602             64=>$PDLA_LL,
603             -32=>$PDLA_F,
604             -64=>$PDLA_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 8 50   8   21 print "Reading IMAGE data...\n" if($PDLA::verbose);
618 8         15 my $fh = shift; # file handle to read from
619 8         12 my $foo = shift; # $foo contains the pre-read header
620 8         11 my $opt = shift; # $opt contains the option hash
621 8         24 my $pdl = shift; # $pdl contains a pre-blessed virgin PDLA
622              
623             # Setup piddle structure
624            
625 8 50       38 if( defined($type_table->{0 + $foo->{"BITPIX"}}) ) {
626 8         47 $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 8         12 my @dims; # Store the dimenions 1..N, compute total number of pixels
632 8         15 my $i = 1;
633 8         14 my $size = 1;
634              
635             ##second part of the conditional guards against a poorly-written hdr.
636 8   66     48 while(defined( $$foo{"NAXIS$i"} ) && $i <= $$foo{"NAXIS"}) {
637 12         28 $size *= $$foo{"NAXIS$i"};
638 12         30 push @dims, $$foo{"NAXIS$i"} ; $i++;
  12         41  
639             }
640 8         57 $pdl->setdims([@dims]);
641            
642 8         36 my $dref = $pdl->get_dataref();
643            
644 8 50       29 print "BITPIX = ",$$foo{"BITPIX"}," size = $size pixels \n"
645             if $PDLA::verbose;
646            
647             # Slurp the FITS binary data
648            
649 8 50       25 print "Reading ",$size*PDLA::Core::howbig($pdl->get_datatype) , " bytes\n"
650             if $PDLA::verbose;
651            
652             # Read the data and pad to the next HDU
653 8         64 my $rdct = $size * PDLA::Core::howbig($pdl->get_datatype);
654 8         32 $fh->read( $$dref, $rdct );
655 8         61 $fh->read( my $dummy, 2880 - (($rdct-1) % 2880) - 1 );
656 8         68 $pdl->upd_data();
657              
658 8 50       31 if (!isbigendian() ) {
659             # Need to byte swap on little endian machines
660 8 50       38 bswap2($pdl) if $pdl->get_datatype == $PDLA_S;
661 8 100 66     76 bswap4($pdl) if $pdl->get_datatype == $PDLA_L || $pdl->get_datatype == $PDLA_F;
662 8 100 100     115 bswap8($pdl) if $pdl->get_datatype == $PDLA_D || $pdl->get_datatype==$PDLA_LL;
663             }
664            
665 8 50       30 if(exists $opt->{bscale}) {
666 8         22 $pdl = treat_bscale($pdl, $foo);
667             }
668            
669             # Header
670            
671 8         36 $pdl->sethdr($foo);
672              
673 8         32 $pdl->hdrcpy($opt->{hdrcpy});
674              
675 8         22 return $pdl;
676             }
677              
678             sub treat_bscale($$){
679 8     8 0 15 my $pdl = shift;
680 8         14 my $foo = shift;
681              
682 8 50       22 print "treating bscale...\n" if($PDLA::debug);
683              
684 8 50       21 if ( $PDLA::Bad::Status ) {
685             # do we have bad values? - needs to be done before BSCALE/BZERO
686             # (at least for integers)
687             #
688 8 50 66     47 if ( $$foo{BITPIX} > 0 and exists $$foo{BLANK} ) {
    100          
689             # integer, so bad value == BLANK keyword
690 0         0 my $blank = $foo->{BLANK};
691             # do we have to do any conversion?
692 0 0       0 if ( $blank == $pdl->badvalue() ) {
693 0         0 $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 6         26 $pdl->inplace->setnantobad();
703             }
704 8 50 66     70 print "FITS file may contain bad values.\n"
705             if $pdl->badflag() and $PDLA::verbose;
706             } # if: PDLA::Bad::Status
707            
708 8         19 my ($bscale, $bzero);
709 8         16 $bscale = $$foo{"BSCALE"}; $bzero = $$foo{"BZERO"};
  8         15  
710 8 50       27 print "BSCALE = $bscale && BZERO = $bzero\n" if $PDLA::verbose;
711 8 50 33     37 $bscale = 1 if (!defined($bscale) || $bscale eq "");
712 8 50 33     33 $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 8         12 my $tmp;
720              
721 8 100       47 if ( $pdl->badflag() == 0 ) {
    50          
722 4         13 $tmp = $pdl->flat()->slice("0:0");
723             } elsif ( $pdl->ngood > 0 ) {
724 4         15 my $index = which( $pdl->flat()->isbad() == 0 )->at(0);
725 4         40 $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 $PDLA::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 8 50       41 $tmp = $tmp*$bscale if $bscale != 1; # Dummy run on one element
738 8 50       24 $tmp = $tmp+$bzero if $bzero != 0;
739            
740 8 50       35 $pdl = $pdl->convert($tmp->type) if $tmp->get_datatype != $pdl->get_datatype;
741            
742 8 50       40 $pdl *= $bscale if $bscale != 1;
743 8 50       21 $pdl += $bzero if $bzero != 0;
744            
745 8         31 delete $$foo{"BSCALE"}; delete $$foo{"BZERO"};
  8         15  
746 8         46 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 PDLA 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             $PDLA::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 += PDLA::Core::howbig($type) * $nrows * 2;
876 0         0 return PDLA->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 * PDLA::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 * PDLA::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'.(PDLA::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 PDLA. 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 = PDLA->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 * PDLA::Core::howbig($type); # rpt should be unity, otherwise we'd have to multiply it in.
935 0         0 my $readlen = $oflen->at(0) * PDLA::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($PDLA::debug);
942              
943             # Copy the data into the output PDLA.
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 PDLA::wfits.\n";
951             }
952              
953             sub _fnP {
954 0     0   0 my( $type, $pdl, $n, $hdr, $opt ) = @_;
955 0         0 my $post = PDLA::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 PDLA or list ref according to the
979             # header.
980             #
981              
982             sub _rfits_bintable ($$$$) {
983 0     0   0 my $fh = shift;
984 0         0 my $hdr = shift;
985 0         0 my $opt = shift;
986             ##shift; ### (ignore $pdl argument)
987              
988 0 0       0 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 0         0 my $tbl = {}; # Table is indexed by name
992 0         0 $tbl->{hdr} = $hdr;
993 0         0 $tbl->{tbl} = 'binary';
994              
995 0         0 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 0 0       0 barf "Binary extension has no fields (TFIELDS=0)" unless($hdr->{TFIELDS});
1002 0         0 my $rowlen = 0;
1003            
1004 0         0 for my $i(1..$hdr->{TFIELDS}) {
1005 0         0 my $iter;
1006 0   0     0 my $name = $tmp->[$i]->{name} = $hdr->{"TTYPE$i"} || "COL";
1007            
1008             ### Allocate some temp space for dealing with this column
1009 0         0 my $tmpcol = $tmp->[$i] = {};
1010            
1011             ### Check for duplicate name and change accordingly...
1012 0   0     0 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 0 0       0 $hdr->{"TTYPE$i"} = $name unless($hdr->{"TTYPE$i"} eq $name);
1019 0         0 $tmpcol->{name} = $name;
1020              
1021 0 0       0 if( ($hdr->{"TFORM$i"}) =~ m/(\d*)(P?.)(.*)/ ) {
1022 0         0 ($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 0   0     0 $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 0 0       0 $tmpcol->{rpt} = PDLA::ceil($tmpcol->{rpt}/8) if ($tmpcol->{type} eq 'X');
1036              
1037             $tmpcol->{handler} = # sic - assignment
1038             $PDLA::IO::FITS_bintable_handlers->{ $tmpcol->{type} }
1039             or
1040             barf "Unknown type ".$hdr->{"TFORM$i"}." in BINTABLE column $i "."("
1041 0 0       0 . $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 0         0 my $foo = $tmpcol->{handler}->[0];
1048 0 0       0 if( ref ($foo) eq 'CODE' ) {
1049             $tmpcol->{data} = $tbl->{$name} =
1050 0         0 &{$foo}( $tmpcol->{rpt}
1051             , $tmpcol->{extra}
1052             , $hdr->{NAXIS2}
1053 0         0 , \$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             PDLA->new_from_specification(
1061             $foo
1062             , $tmpcol->{rpt},
1063 0   0     0 , $hdr->{NAXIS2} || 1
1064             );
1065              
1066 0         0 $rowlen += PDLA::Core::howbig($foo) * $tmpcol->{rpt};
1067             }
1068            
1069 0 0       0 print "Prefrobnicated col. $i "."(".$hdr->{"TTYPE$i"}.")\ttype is ".$hdr->{"TFORM$i"}."\t length is now $rowlen\n" if($PDLA::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 \$PDLA::debug for more detailed info)\n"
1077 0 0       0 if($rowlen != $hdr->{NAXIS1});
1078            
1079             ### Snarf up the whole extension, and pad to 2880 bytes...
1080 0         0 my ($rawtable, $heap, $n1, $n2);
1081              
1082             # n1 gets number of bytes in table plus gap
1083 0         0 $n1 = $hdr->{NAXIS1} * $hdr->{NAXIS2};
1084 0 0       0 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 0 0       0 $n2 = $hdr->{PCOUNT} + ($hdr->{THEAP} ? ($hdr->{NAXIS1}*$hdr->{NAXIS2} - $hdr->{THEAP}) : 0);
1094 0         0 $n2 = ($n1+$n2-1)+2880 - (($n1+$n2-1) % 2880) - $n1;
1095              
1096 0 0       0 print "Reading $n1 bytes of table data and $n2 bytes of heap data....\n"
1097             if($PDLA::verbose);
1098 0         0 $fh->read($rawtable, $n1);
1099              
1100 0 0       0 if($n2) {
1101 0         0 $fh->read($heap, $n2)
1102             } else {
1103 0         0 $heap = which(pdl(0)); # empty PDLA
1104             }
1105              
1106             ### Frobnicate the rows, one at a time.
1107 0         0 for my $row(0..$hdr->{NAXIS2}-1) {
1108 0         0 my $prelen = length($rawtable);
1109 0         0 for my $i(1..$hdr->{TFIELDS}) {
1110 0         0 my $tmpcol = $tmp->[$i];
1111 0         0 my $reader = $tmpcol->{handler}->[1];
1112 0 0       0 if(ref $reader eq 'CODE') {
    0          
1113 0         0 &{$reader}( $tmpcol->{data}
1114             , $row
1115             , \$rawtable
1116             , $tmpcol->{rpt}
1117             , $tmpcol->{extra}
1118 0         0 , \$heap
1119             , $tbl
1120             , $i
1121             );
1122             } elsif(ref $tmpcol->{data} eq 'PDLA') {
1123 0         0 my $rlen = $reader * $tmpcol->{rpt};
1124              
1125 0         0 substr( ${$tmpcol->{data}->get_dataref()}, $rlen * $row, $rlen ) =
  0         0  
1126             substr( $rawtable, 0, $rlen, '');
1127 0         0 $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 0 0       0 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 0         0 for my $i(1..$hdr->{TFIELDS}) { # Postfrobnication loop
1148 0         0 my $tmpcol = $tmp->[$i];
1149 0         0 my $post = $tmpcol->{handler}->[3];
1150            
1151 0 0       0 if(ref $post eq 'CODE') {
    0          
    0          
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 'PDLA' ) {
1157             # Do standard PDLA-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 0 0       0 unless( isbigendian() ) {
1166 0 0       0 if( $post == 2 ) { bswap2($tmpcol->{data}); }
  0 0       0  
    0          
    0          
1167 0         0 elsif( $post == 4 ) { bswap4($tmpcol->{data}); }
1168 0         0 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 0 0       0 if($opt->{bscale}) {
1180 0 0       0 my $tzero = defined($hdr->{"TZERO$i"}) ? $hdr->{"TZERO$i"} : 0.0;
1181 0 0       0 my $tscal = defined($hdr->{"TSCAL$i"}) ? $hdr->{"TSCAL$i"} : 1.0;
1182            
1183             # The $valid_ flags let us avoid unnecessary arithmetic.
1184 0         0 my $valid_tzero = ($tzero != 0.0);
1185 0         0 my $valid_tscal = ($tscal != 1.0);
1186            
1187 0 0 0     0 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 PDLA'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 PDLA 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 0         0 delete $hdr->{"TZERO$i"};
1229 0         0 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 0         0 my @tdims = ();
1244 0         0 $tmpcol->{data}->hdrcpy(1);
1245              
1246 0 0       0 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 PDLA out to the table itself.
1266 0 0 0     0 if($hdr->{NAXIS2} > 0 && $tmpcol->{rpt}>0) {
1267             $tbl->{$tmpcol->{name}} =
1268             ( ( $tmpcol->{data}->dim(0) == 1 )
1269             ? $tmpcol->{data}->slice("(0)")
1270 0 0       0 : $tmpcol->{data}->xchg(0,1)
1271             );
1272             }
1273             }
1274              
1275             # End of PDLA 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 0 0 0     0 if($hdr->{ZIMAGE} && $hdr->{ZCMPTYPE} && $opt->{expand}) {
      0        
1286 0         0 eval 'use PDLA::Compression;';
1287 0 0       0 if($@) {
1288 0         0 die "rfits: error while loading PDLA::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 0         0 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} = PDLA::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 PDLA 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( PDLA::howbig($compressed->get_datatype) != $bytepix ) {
1352             my @dims = $compressed->dims;
1353             my $newdim0;
1354             my $scaledim0;
1355              
1356             $scaledim0 = $dims[0] * PDLA::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 / PDLA::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 = PDLA->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 = PDLA->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 = PDLA::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 if(exists $hdr->{BSCALE}) {
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 PDLA 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             #expand leading ~ to home directory (using glob())
1574             wfits $pdl, '~/filename.fits';
1575              
1576             =over 3
1577              
1578             =item * Ordinary (PDLA) data handling:
1579              
1580             If the first argument is a PDLA, then the PDLA is written out as an
1581             ordinary FITS file with a single Header/Data Unit of data.
1582              
1583             $BITPIX is then optional and coerces the output data type according to
1584             the standard FITS convention for the BITPIX field (with positive
1585             values representing integer types and negative values representing
1586             floating-point types).
1587              
1588             If C<$pdl> has a FITS header attached to it (actually, any hash that
1589             contains a C<< SIMPLE=>T >> keyword), then that FITS header is written
1590             out to the file. The image dimension tags are adjusted to the actual
1591             dataset. If there's a mismatch between the dimensions of the data and
1592             the dimensions in the FITS header, then the header gets corrected and
1593             a warning is printed.
1594              
1595             If C<$pdl> is a slice of another PDLA with a FITS header already
1596             present (and header copying enabled), then you must be careful.
1597             C will remove any extraneous C keywords (per the FITS
1598             standard), and also remove the other keywords associated with that
1599             axis: C, C, C, C, and C. This
1600             may cause confusion if the slice is NOT out of the last dimension:
1601             C and you would be best off adjusting
1602             the header yourself before calling C.
1603              
1604             You can tile-compress images according to the CFITSIO extension to the
1605             FITS standard, by adding an option hash to the arguments:
1606              
1607             =over 3
1608              
1609             =item compress
1610              
1611             This can be either unity, in which case Rice compression is used,
1612             or a (case-insensitive) string matching the CFITSIO compression
1613             type names. Currently supported compression algorithms are:
1614              
1615             =over 3
1616              
1617             =item * RICE_1 - linear Rice compression
1618              
1619             This uses limited-symbol-length Rice compression, which works well on
1620             low entropy image data (where most pixels differ from their neighbors
1621             by much less than the dynamic range of the image).
1622              
1623             =back
1624              
1625             =item tilesize (default C<[-1,1]>)
1626              
1627             This specifies the dimension of the compression tiles, in pixels. You
1628             can hand in a PDLA, a scalar, or an array ref. If you specify fewer
1629             dimensions than exist in the image, the last dim is repeated - so "32"
1630             yields 32x32 pixel tiles in a 2-D image. A dim of -1 in any dimension
1631             duplicates the image size, so the default C<[-1,1]> causes compression
1632             along individual rows.
1633              
1634             =item tilesize (RICE_1 only; default C<32>)
1635              
1636             For RICE_1, BLOCKSIZE indicates the number of pixel samples to use
1637             for each compression block within the compression algorithm. The
1638             blocksize is independent of the tile dimensions. For RICE
1639             compression the pixels from each tile are arranged in normal pixel
1640             order (early dims fastest) and compressed as a linear stream.
1641              
1642             =back
1643              
1644             =item * Table handling:
1645              
1646             If you feed in a hash ref instead of a PDLA, then the hash ref is
1647             written out as a binary table extension. The hash ref keys are
1648             treated as column names, and their values are treated as the data to
1649             be put in each column.
1650              
1651             For numeric information, the hash values should contain PDLAs. The 0th
1652             dim of the PDLA runs across rows, and higher dims are written as
1653             multi-value entries in the table (e.g. a 7x5 PDLA will yield a single
1654             named column with 7 rows and 5 numerical entries per row, in a binary
1655             table). Note that this is slightly different from the usual concept
1656             of threading, in which dimension 1 runs across rows.
1657              
1658             ASCII tables only allow one entry per column in each row, so
1659             if you plan to write an ASCII table then all of the values of C<$hash>
1660             should have at most one dim.
1661              
1662             All of the columns' 0 dims must agree in the threading sense. That is to
1663             say, the 0th dimension of all of the values of C<$hash> should be the
1664             same (indicating that all columns have the same number of rows). As
1665             an exception, if the 0th dim of any of the values is 1, or if that
1666             value is a PDLA scalar (with 0 dims), then that value is "threaded"
1667             over -- copied into all rows.
1668              
1669             Data dimensions higher than 2 are preserved in binary tables,
1670             via the TDIMn field (e.g. a 7x5x3 PDLA is stored internally as
1671             seven rows with 15 numerical entries per row, and reconstituted
1672             as a 7x5x3 PDLA on read).
1673              
1674             Non-PDLA Perl scalars are treated as strings, even if they contain
1675             numerical values. For example, a list ref containing 7 values is
1676             treated as 7 rows containing one string each. There is no such thing
1677             as a multi-string column in FITS tables, so any nonscalar values in
1678             the list are stringified before being written. For example, if you
1679             pass in a perl list of 7 PDLAs, each PDLA will be stringified before
1680             being written, just as if you printed it to the screen. This is
1681             probably not what you want -- you should use L to connect
1682             the separate PDLAs into a single one. (e.g. C<$x-Eglue(1,$y,$c)-Emv(1,0)>)
1683              
1684             The column names are case-insensitive, but by convention the keys of
1685             C<$hash> should normally be ALL CAPS, containing only digits, capital
1686             letters, hyphens, and underscores. If you include other characters,
1687             then case is smashed to ALL CAPS, whitespace is converted to
1688             underscores, and unrecognized characters are ignored -- so if you
1689             include the key "Au Purity (%)", it will be written to the file as a
1690             column that is named "AU_PURITY". Since this is not guaranteed to
1691             produce unique column names, subsequent columns by the same name are
1692             disambiguated by the addition of numbers.
1693              
1694             You can specify the use of variable-length rows in the output, saving
1695             space in the file. To specify variable length rows for a column named
1696             "FOO", you can include a separate key "len_FOO" in the hash to be
1697             written. The key's value should be a PDLA containing the number of
1698             actual samples in each row. The result is a FITS P-type variable
1699             length column that, upon read with C, will restore to a field
1700             named FOO and a corresponding field named "len_FOO". Invalid data in
1701             the final PDLA consist of a padding value (which defaults to 0 but
1702             which you may set by including a TNULL field in the hdr specificaion).
1703             Variable length arrays must be 2-D PDLAs, with the variable length in
1704             the 1 dimension.
1705              
1706             Two further special keys, 'hdr' and 'tbl', can contain
1707             meta-information about the type of table you want to write. You may
1708             override them by including an C<$OPTIONS> hash with a 'hdr' and/or
1709             'tbl' key.
1710              
1711             The 'tbl' key, if it exists, must contain either 'ASCII' or 'binary'
1712             (case-insensitive), indicating whether to write an ascii or binary
1713             table. The default is binary. [ASCII table writing is planned but
1714             does not yet exist].
1715              
1716             You can specify the format of the table quite specifically with the
1717             'hdr' key or option field. If it exists, then the 'hdr' key should
1718             contain fields appropriate to the table extension being used. Any
1719             field information that you don't specify will be filled in
1720             automatically, so (for example) you can specify that a particular
1721             column name goes in a particular position, but allow C to
1722             arrange the other columns in the usual alphabetical order into any
1723             unused slots that you leave behind. The C, C,
1724             C, C, C, and C keywords are ignored:
1725             their values are calculated based on the hash that you supply. Any
1726             other fields are passed into the final FITS header verbatim.
1727              
1728             As an example, the following
1729              
1730             $x = long(1,2,4);
1731             $y = double(1,2,4);
1732             wfits { 'COLA'=>$x, 'COLB'=>$y }, "table1.fits";
1733              
1734             will create a binary FITS table called F which
1735             contains two columns called C and C. The order
1736             of the columns is controlled by setting the C
1737             keywords in the header array, so
1738              
1739             $h = { 'TTYPE1'=>'Y', 'TTYPE2'=>'X' };
1740             wfits { 'X'=>$x, 'Y'=>$y, hdr=>$h }, "table2.fits";
1741              
1742             creates F where the first column is
1743             called C and the second column is C.
1744              
1745             =item * multi-value handling
1746              
1747             If you feed in a perl list rather than a PDLA or a hash, then
1748             each element is written out as a separate HDU in the FITS file.
1749             Each element of the list must be a PDLA or a hash. [This is not implemented
1750             yet but should be soon!]
1751              
1752             =item * DEVEL NOTES
1753              
1754             ASCII tables are not yet handled but should be.
1755              
1756             Binary tables currently only handle one vector (up to 1-D array)
1757             per table entry; the standard allows more, and should be fully implemented.
1758             This means that PDLA::Complex piddles currently can not be written to
1759             disk.
1760              
1761             Handling multidim arrays implies that perl multidim lists should also be
1762             handled.
1763              
1764             =back
1765              
1766             =for bad
1767              
1768             For integer types (ie C 0>), the C keyword is set
1769             to the bad value. For floating-point types, the bad value is
1770             converted to NaN (if necessary) before writing.
1771              
1772             =cut
1773              
1774             *wfits = \&PDLA::wfits;
1775              
1776             BEGIN {
1777 15     15   118 @PDLA::IO::FITS::wfits_keyword_order =
1778             ('SIMPLE','BITPIX','NAXIS','NAXIS1','BUNIT','BSCALE','BZERO');
1779 15         100732 @PDLA::IO::FITS::wfits_numbered_keywords =
1780             ('CTYPE','CRPIX','CRVAL','CDELT','CROTA');
1781             }
1782              
1783             # Until we do a rewrite these have to be file global since they
1784             # are used by the wheader routine
1785             my (%hdr, $nbytes);
1786              
1787             # Local utility routine of wfits()
1788             sub wheader ($$) {
1789 48     48 0 70 my $fh = shift;
1790 48         68 my $k = shift;
1791            
1792 48 100       135 if ($k =~ m/(HISTORY|COMMENT)/) {
1793 8         27 my $hc = $1;
1794 8 50       32 return unless ref($hdr{$k}) eq 'ARRAY';
1795 0         0 foreach my $line (@{$hdr{$k}}) {
  0         0  
1796 0         0 $fh->printf( "$hc %-72s", substr($line,0,72) );
1797 0         0 $nbytes += 80;
1798             }
1799 0         0 delete $hdr{$k};
1800             } else {
1801             # Check that we are dealing with a scalar value in the header
1802             # Need to make sure that the header does not include PDLAs or
1803             # other structures. Return unless $hdr{$k} is a scalar.
1804 40         78 my($hdrk) = $hdr{$k};
1805            
1806 40 100       84 if(ref $hdrk eq 'ARRAY') {
1807 1         6 $hdrk = join("\n",@$hdrk);
1808             }
1809            
1810 40 50       74 return unless not ref($hdrk);
1811            
1812 40 50       80 if ($hdrk eq "") {
1813 0         0 $fh->printf( "%-80s", substr($k,0,8) );
1814             } else {
1815 40         130 $fh->printf( "%-8s= ", substr($k,0,8) );
1816            
1817             my $com = ( ref $hdr{COMMENT} eq 'HASH' ) ?
1818 40 50       341 $hdr{COMMENT}{$k} : undef;
1819            
1820 40 100 33     255 if ($hdrk =~ /^ *([+-]?)(?=\d|\.\d)\d*(\.\d*)?([Ee]([+-]?\d+))? *$/) { # Number?
    50          
1821 31 50       70 my $cl=60-($com ? 2 : 0);
1822 31         61 my $end=' ' x $cl;
1823 31 50       61 $end =' /'. $com if($com);
1824 31         103 $fh->printf( "%20s%-50s", substr($hdrk,0,20),
1825             substr($end, 0, 50) );
1826             } elsif ($hdrk eq 'F' or $hdrk eq 'T') {
1827             # Logical flags ?
1828 0         0 $fh->printf( "%20s", $hdrk );
1829 0         0 my $end=' ' x 50;
1830 0 0       0 $end =' /'.$com if($com);
1831 0         0 $fh->printf( "%-50s", $end );
1832             } else {
1833             # Handle strings, truncating as necessary
1834             # (doesn't do multicard strings like Astro::FITS::Header does)
1835            
1836             # Convert single quotes to doubled single quotes
1837             # (per FITS standard)
1838 9         25 my($st) = $hdrk;
1839 9         29 $st =~ s/\'/\'\'/g;
1840            
1841 9         18 my $sl=length($st)+2;
1842 9 50       41 my $cl=70-$sl-($com ? 2 : 0);
1843 9         56 $fh->print( "'$st'" );
1844            
1845 9 50       76 if (defined $com) {
1846 0         0 $fh->printf( " /%-$ {cl}s", substr($com, 0, $cl) );
1847             } else {
1848 9         47 $fh->printf( "%-$ {cl}s", ' ' x $cl );
1849             }
1850             }
1851             }
1852 40         327 $nbytes += 80; delete $hdr{$k};
  40         79  
1853             }
1854 40 50       92 delete $hdr{COMMENT}{$k} if ref $hdr{COMMENT} eq 'HASH';
1855 40         68 1;
1856             }
1857              
1858             # Write a PDLA to a FITS format file
1859             #
1860             sub PDLA::wfits {
1861 7 50 33 7 0 173 barf 'Usage: wfits($pdl,$file,[$BITPIX],[{options}])' if $#_<1 || $#_>3;
1862 7         26 my ($pdl,$file,$x,$y) = @_;
1863 7         15 my ($opt, $BITPIX);
1864              
1865 7         31 local $\ = undef; # fix sf.net bug #3394327
1866              
1867 7 50       31 if(ref $x eq 'HASH') {
    50          
1868 0         0 $x = $opt;
1869 0         0 $BITPIX = $y;
1870             } elsif(ref $y eq 'HASH') {
1871 0         0 $BITPIX = $x;
1872 0         0 $opt = $y;
1873             }
1874              
1875 7         16 my ($k, $buff, $off, $ndims, $sz);
1876              
1877 7         20 $file =~ s/^(~)/glob($1)/e; # tilde expansion
  0         0  
1878              
1879 7         172 local $SIG{PIPE};
1880              
1881 7 50       50 if ($file =~ /\.gz$/) { # Handle suffix-style compression
    50          
1882 0     0   0 $SIG{PIPE}= sub {}; # Prevent crashing if gzip dies
1883 0         0 $file = "|gzip -9 > $file";
1884             }
1885             elsif ($file =~ /\.Z$/) {
1886 0     0   0 $SIG{PIPE}= sub {}; # Prevent crashing if compress dies
1887 0         0 $file = "|compress > $file";
1888             }
1889             else{
1890 7         27 $file = ">$file";
1891             }
1892            
1893             #### Figure output type
1894              
1895 7         18 my @outputs = ();
1896 7         11 my $issue_nullhdu;
1897              
1898 7 100       38 if( UNIVERSAL::isa($pdl,'PDLA') ) {
    50          
    50          
    0          
1899 6         15 $issue_nullhdu = 0;
1900 6         12 push(@outputs, $pdl);
1901             } elsif( ref($pdl) eq 'HASH' ) {
1902 0         0 $issue_nullhdu = 1;
1903 0         0 push(@outputs, $pdl);
1904             } elsif( ref($pdl) eq 'ARRAY' ) {
1905 1         4 $issue_nullhdu = 1;
1906 1         12 @outputs = @$pdl;
1907             } elsif( length(ref($pdl))==0 ) {
1908 0         0 barf "wfits: needs a HASH or PDLA argument to write out\n";
1909             } else {
1910 0         0 barf "wfits: unknown ref type ".ref($pdl)."\n";
1911             }
1912              
1913             ## Open file & prepare to write binary info
1914 7 50       61 my $fh = IO::File->new( $file )
1915             or barf "Unable to create FITS file $file\n";
1916 7         1123 binmode $fh;
1917              
1918 7 100       22 if($issue_nullhdu) {
1919 1         6 _wfits_nullhdu($fh);
1920             }
1921              
1922 7         41 for $pdl(@outputs) {
1923              
1924 8 50       139 if(ref($pdl) eq 'HASH') {
    50          
1925             my $table_type = ( exists($pdl->{tbl}) ?
1926 0 0       0 ($pdl->{tbl} =~ m/^a/i ? 'ascii' : 'binary') :
    0          
1927             "binary"
1928             );
1929 0         0 _wfits_table($fh,$pdl,$table_type);
1930             } elsif( UNIVERSAL::isa($pdl,'PDLA') ) {
1931              
1932             ### Regular image writing.
1933            
1934 8 100       30 $BITPIX = "" unless defined $BITPIX;
1935 8 100       25 if ($BITPIX eq "") {
1936 7 50       40 $BITPIX = 8 if $pdl->get_datatype == $PDLA_B;
1937 7 50 33     59 $BITPIX = 16 if $pdl->get_datatype == $PDLA_S || $pdl->get_datatype == $PDLA_US;
1938 7 100       34 $BITPIX = 32 if $pdl->get_datatype == $PDLA_L;
1939 7 100       28 $BITPIX = 64 if $pdl->get_datatype == $PDLA_LL;
1940 7 50       24 $BITPIX = -32 if $pdl->get_datatype == $PDLA_F;
1941 7 100       25 $BITPIX = -64 if $pdl->get_datatype == $PDLA_D;
1942 7 50       25 $BITPIX = 8 * PDLA::Core::howbig($PDLA_IND) if($pdl->get_datatype==$PDLA_IND);
1943             }
1944 8 50       23 if ($BITPIX eq "") {
1945 0         0 $BITPIX = -64;
1946 0         0 warn "wfits: PDLA has an unsupported datatype -- defaulting to 64-bit float.\n";
1947             }
1948              
1949 8     0   44 my $convert = sub { return $_[0] }; # Default - do nothing
  0         0  
1950 8 50   0   24 $convert = sub {byte($_[0])} if $BITPIX == 8;
  0         0  
1951 8 50   0   27 $convert = sub {short($_[0])} if $BITPIX == 16;
  0         0  
1952 8 100   2   34 $convert = sub {long($_[0])} if $BITPIX == 32;
  2         9  
1953 8 100   2   25 $convert = sub {longlong($_[0])} if $BITPIX == 64;
  2         9  
1954 8 50   0   19 $convert = sub {float($_[0])} if $BITPIX == -32;
  0         0  
1955 8 100   12   62 $convert = sub {double($_[0])} if $BITPIX == -64;
  12         46  
1956            
1957             # Automatically figure output scaling
1958            
1959 8         14 my $bzero = 0; my $bscale = 1;
  8         17  
1960 8 100       19 if ($BITPIX>0) {
1961 2         39 my $min = $pdl->min;
1962 2         9 my $max = $pdl->max;
1963 2 50       7 my ($dmin,$dmax) = (0, 2**8-1) if $BITPIX == 8;
1964 2 50       7 ($dmin,$dmax) = (-2**15, 2**15-1) if $BITPIX == 16;
1965 2 100       7 ($dmin,$dmax) = (-2**31, 2**31-1) if $BITPIX == 32;
1966 2 100       9 ($dmin,$dmax) = (-(pdl(longlong,1)<<63), (pdl(longlong,1)<<63)-1) if $BITPIX==64;
1967              
1968 2 50 33     76 if ($min<$dmin || $max>$dmax) {
1969 0         0 $bzero = $min - $dmin;
1970 0         0 $max -= $bzero;
1971 0 0       0 $bscale = $max/$dmax if $max>$dmax;
1972             }
1973 2 50       22 print "BSCALE = $bscale && BZERO = $bzero\n" if $PDLA::verbose;
1974             }
1975            
1976             # Check for tile-compression format for the image, and handle it.
1977             # We add the image-compression format tags and reprocess the whole
1978             # shebang as a binary table.
1979 8 50       26 if($opt->{compress}) {
1980 0         0 croak "Placeholder -- tile compression not yet supported\n";
1981             }
1982            
1983            
1984             ##############################
1985             ## Check header and prepare to write it out
1986            
1987 8         37 my($h) = $pdl->gethdr();
1988            
1989             # Extra logic: if we got handed a vanilla hash that that is *not* an Astro::FITS::Header, but
1990             # looks like it's a FITS header encoded in a hash, then attempt to process it with
1991             # Astro::FITS::Header before writing it out -- this helps with cleanup of tags.
1992 8 0 33     29 if($PDLA::Astro_FITS_Header and
      33        
      0        
1993             defined($h) and
1994             ref($h) eq 'HASH' and
1995             !defined( tied %$h )
1996             ) {
1997            
1998 0         0 my $all_valid_fits = 1;
1999 0         0 for my $k(keys %$h) {
2000 0 0 0     0 if(length($k) > 8 or
2001             $k !~ m/^[A-Z_][A-Z\d\_]*$/i
2002             ) {
2003 0         0 $all_valid_fits = 0;
2004 0         0 last;
2005             }
2006             }
2007            
2008 0 0       0 if($all_valid_fits) {
2009             # All the keys look like valid FITS header keywords -- so
2010             # create a tied FITS header object and use that instead.
2011 0         0 my $afh = new Astro::FITS::Header( );
2012 0         0 my %hh;
2013 0         0 tie %hh, "Astro::FITS::Header", $afh;
2014 0         0 for (keys %$h) {
2015 0         0 $hh{$_} = $h->{$_};
2016             }
2017 0         0 $h = \%hh;
2018             }
2019             }
2020            
2021             # Now decide whether to emit a hash or an AFH object
2022 8 50 33     57 if(defined($h) &&
      66        
2023             ( (defined (tied %$h)) &&
2024             (UNIVERSAL::isa(tied %$h,"Astro::FITS::Header")))
2025             ){
2026 0         0 my $k;
2027            
2028             ##n############################
2029             ## Tied-hash code -- I'm too lazy to incorporate this into KGB's
2030             ## direct hash handler below, so I've more or less just copied and
2031             ## pasted with some translation. --CED
2032             ##
2033 0         0 my $hdr = tied %$h;
2034            
2035             #
2036             # Put advertising comment in the SIMPLE field
2037             #n
2038 0 0       0 if($issue_nullhdu) {
2039 0         0 $h->{XTENSION} = "IMAGE";
2040             } else {
2041 0         0 $h->{SIMPLE} = 'T';
2042 0         0 my(@a) = $hdr->itembyname('SIMPLE');
2043 0         0 $a[0]->comment('Created with PDLA (http://pdl.perl.org)');
2044             # and register it as a LOGICAL rather than a string
2045 0         0 $a[0]->type('LOGICAL');
2046             }
2047            
2048             #
2049             # Use tied interface to set all the keywords. Note that this
2050             # preserves existing per-line comments, only changing the values.
2051             #
2052 0         0 $h->{BITPIX} = $BITPIX;
2053 0         0 $h->{NAXIS} = $pdl->getndims;
2054 0         0 my $correction = 0;
2055 0         0 for $k(1..$h->{NAXIS}) {
2056             $correction |= (exists $h->{"NAXIS$k"} and
2057 0   0     0 $h->{"NAXIS$k"} != $pdl->dim($k-1)
2058             );
2059 0         0 $h->{"NAXIS$k"} = $pdl->getdim($k-1);
2060             }
2061 0 0       0 carp("Warning: wfits corrected dimensions of FITS header")
2062             if($correction);
2063            
2064 0 0       0 $h->{BUNIT} = "Data Value" unless exists $h->{BUNIT};
2065 0 0       0 $h->{BSCALE} = $bscale if($bscale != 1);
2066 0 0       0 $h->{BZERO} = $bzero if($bzero != 0);
2067            
2068 0 0       0 if ( $pdl->badflag() ) {
2069 0 0       0 if ( $BITPIX > 0 ) { my $x = &$convert(pdl(0.0));
  0         0  
2070 0         0 $h->{BLANK} = $x->badvalue(); }
2071 0         0 else { delete $h->{BLANK}; }
2072             }
2073            
2074             # Use object interface to sort the lines. This is complicated by
2075             # the need for an arbitrary number of NAXIS lines in the middle
2076             # of the sorting. Keywords with a trailing '1' in the sorted-order
2077             # list get looped over.
2078 0         0 my($kk) = 0;
2079 0         0 my(@removed_naxis) = ();
2080 0         0 for $k(0..$#PDLA::IO::FITS::wfits_keyword_order) {
2081 0         0 my($kn) = 0;
2082            
2083 0         0 my @index;
2084 0   0     0 do { # Loop over numericised keywords (e.g. NAXIS1)
2085            
2086 0         0 my $kw = $PDLA::IO::FITS::wfits_keyword_order[$k]; # $kw get keyword
2087 0 0       0 $kw .= (++$kn) if( $kw =~ s/\d$//); # NAXIS1 -> NAXIS
2088 0         0 @index = $hdr->index($kw);
2089            
2090 0 0       0 if(defined $index[0]) {
2091 0 0       0 if($kn <= $pdl->getndims){
2092 0 0       0 $hdr->insert($kk, $hdr->remove($index[0]))
2093             unless ($index[0] == $kk) ;
2094 0         0 $kk++;
2095             }
2096             else{ #remove e.g. NAXIS3 from hdr if NAXIS==2
2097 0         0 $hdr->removebyname($kw);
2098 0         0 push(@removed_naxis,$kw);
2099             }
2100             }
2101             } while((defined $index[0]) && $kn);
2102             }
2103            
2104 0         0 foreach my $naxis(@removed_naxis){
2105 0         0 $naxis =~ m/(\d)$/;
2106 0         0 my $n = $1;
2107 0         0 foreach my $kw(@PDLA::IO::FITS::wfits_numbered_keywords){
2108 0         0 $hdr->removebyname($kw . $n);
2109             }
2110             }
2111             #
2112             # Delete the END card if necessary (for later addition at the end)
2113             #
2114 0         0 $hdr->removebyname('END');
2115            
2116             #
2117             # Make sure that the HISTORY lines all come at the end
2118             #
2119 0         0 my @hindex = $hdr->index('HISTORY');
2120 0         0 for $k(0..$#hindex) {
2121 0         0 $hdr->insert(-1-$k, $hdr->remove($hindex[-1-$k]));
2122             }
2123            
2124             #
2125             # Make sure the last card is an END
2126             #
2127 0         0 $hdr->insert(scalar($hdr->cards),
2128             Astro::FITS::Header::Item->new(Keyword=>'END'));
2129            
2130             #
2131             # Write out all the cards, and note how many bytes for later padding.
2132             #
2133 0         0 my $s = join("",$hdr->cards);
2134            
2135 0         0 $fh->print( $s );
2136 0         0 $nbytes = length $s;
2137             } else {
2138             ##
2139             ## Legacy emitter (note different advertisement in the SIMPLE
2140             ## comment, for debugging!)
2141             ##
2142              
2143 8 100       22 if($issue_nullhdu) {
2144 2         9 $fh->printf( "%-80s", "XTENSION= 'IMAGE'" );
2145             } else {
2146 6         42 $fh->printf( "%-80s", "SIMPLE = T / PDLA::IO::FITS::wfits (http://pdl.perl.org)" );
2147             }
2148            
2149 8         200 $nbytes = 80; # Number of bytes written so far
2150            
2151             # Write FITS header
2152            
2153 8         18 %hdr = ();
2154 8 100       66 if (defined($h)) {
2155 1         5 for (keys %$h) { $hdr{uc $_} = $$h{$_} } # Copy (ensuring keynames are uppercase)
  4         12  
2156             }
2157            
2158 8         19 delete $hdr{SIMPLE}; delete $hdr{'END'};
  8         15  
2159            
2160 8         48 $hdr{BITPIX} = $BITPIX;
2161 8 50       33 $hdr{BUNIT} = "Data Value" unless exists $hdr{BUNIT};
2162 8         28 wheader($fh, 'BITPIX');
2163            
2164 8         33 $ndims = $pdl->getndims; # Dimensions of data array
2165 8         20 $hdr{NAXIS} = $ndims;
2166 8         24 wheader($fh, 'NAXIS');
2167 8         26 for $k (1..$ndims) { $hdr{"NAXIS$k"} = $pdl->getdim($k-1) }
  12         66  
2168 8         22 for $k (1..$ndims) { wheader($fh,"NAXIS$k") }
  12         33  
2169            
2170 8 50 33     44 if ($bscale != 1 || $bzero != 0) {
2171 0         0 $hdr{BSCALE} = $bscale;
2172 0         0 $hdr{BZERO} = $bzero;
2173 0         0 wheader($fh,'BSCALE');
2174 0         0 wheader($fh,'BZERO');
2175             }
2176 8         21 wheader($fh,'BUNIT');
2177            
2178             # IF badflag is set
2179             # and BITPIX > 0 - ensure the header contains the BLANK keyword
2180             # (make sure it's for the correct type)
2181             # otherwise - make sure the BLANK keyword is removed
2182 8 100       62 if ( $pdl->badflag() ) {
2183 4 50       12 if ( $BITPIX > 0 ) { my $x = &$convert(pdl(0.0)); $hdr{BLANK} = $x->badvalue(); }
  0         0  
  0         0  
2184 4         8 else { delete $hdr{BLANK}; }
2185             }
2186            
2187 8         60 for $k (sort fits_field_cmp keys %hdr) {
2188 4 50       13 wheader($fh,$k) unless $k =~ m/HISTORY/;
2189             }
2190 8         28 wheader($fh, 'HISTORY'); # Make sure that HISTORY entries come last.
2191 8         40 $fh->printf( "%-80s", "END" );
2192 8         62 $nbytes += 80;
2193             }
2194            
2195             #
2196             # Pad the header to a legal value and write the rest of the FITS file.
2197             #
2198 8         22 $nbytes %= 2880;
2199 8 50       103 $fh->print( " "x(2880-$nbytes) )
2200             if $nbytes != 0; # Fill up HDU
2201            
2202             # Decide how to byte swap - note does not quite work yet. Needs hack
2203             # to IO.xs
2204            
2205 8     0   76 my $bswap = sub {}; # Null routine
2206 8 50       36 if ( !isbigendian() ) { # Need to set a byte swap routine
2207 8 50       23 $bswap = \&bswap2 if $BITPIX==16;
2208 8 100 66     42 $bswap = \&bswap4 if $BITPIX==32 || $BITPIX==-32;
2209 8 100 100     43 $bswap = \&bswap8 if $BITPIX==-64 || $BITPIX==64;
2210             }
2211            
2212             # Write FITS data
2213            
2214 8         34 my $p1d = $pdl->copy->reshape($pdl->nelem); # Data as 1D stream;
2215            
2216 8         16 $off = 0;
2217 8         33 $sz = PDLA::Core::howbig(&$convert($p1d->slice('0:0'))->get_datatype);
2218            
2219 8         67 $nbytes = $p1d->getdim(0) * $sz;
2220            
2221             # Transfer data in blocks (because might need to byte swap)
2222             # Buffer is also type converted on the fly
2223            
2224 8         16 my $BUFFSZ = 360*2880; # = ~1Mb - must be multiple of 2880
2225 8         14 my $tmp;
2226            
2227 8 50 66     52 if ( $pdl->badflag() and $BITPIX < 0 and $PDLA::Bad::UseNaN == 0 ) {
      66        
2228             # just print up a message - conversion is actually done in the loop
2229 4 50       14 print "Converting PDLA bad value to NaN\n" if $PDLA::verbose;
2230             }
2231            
2232 8         35 while ($nbytes - $off > $BUFFSZ) {
2233            
2234             # Data to be transferred
2235            
2236 0         0 $buff = &$convert( ($p1d->slice( ($off/$sz).":". (($off+$BUFFSZ)/$sz-1))
2237             -$bzero)/$bscale );
2238            
2239             # if there are bad values present, and output type is floating-point,
2240             # convert the bad values to NaN's. We can ignore integer types, since
2241             # we have set the BLANK keyword
2242             #
2243 0 0 0     0 if ( $pdl->badflag() and $BITPIX < 0 and $PDLA::Bad::UseNaN == 0 ) {
      0        
2244 0         0 $buff->inplace->setbadtonan();
2245             }
2246            
2247 0         0 &$bswap($buff);
2248 0         0 $fh->print( ${$buff->get_dataref} );
  0         0  
2249 0         0 $off += $BUFFSZ;
2250             }
2251 8         45 $buff = &$convert( ($p1d->slice($off/$sz.":-1") - $bzero)/$bscale );
2252            
2253 8 50 66     166 if ( $pdl->badflag() and $BITPIX < 0 and $PDLA::Bad::UseNaN == 0 ) {
      66        
2254 4         17 $buff->inplace->setbadtonan();
2255             }
2256            
2257 8         103 &$bswap($buff);
2258 8         26 $fh->print( ${$buff->get_dataref} );
  8         65  
2259             # Fill HDU and close
2260             # note that for the data space the fill character is \0 not " "
2261             #
2262 8         151 $fh->print( "\0"x(($BUFFSZ - $buff->getdim(0) * $sz)%2880) );
2263             } # end of image writing block
2264             else {
2265             # Not a PDLA and not a hash ref
2266 0         0 barf("wfits: unknown data type - quitting");
2267             }
2268             } # end of output loop
2269 7         138 $fh->close();
2270 7         673 1;
2271             }
2272              
2273             ######################################################################
2274             ######################################################################
2275              
2276             # Compare FITS headers in a sensible manner.
2277              
2278             =head2 fits_field_cmp
2279              
2280             =for ref
2281              
2282             fits_field_cmp
2283              
2284             Sorting comparison routine that makes proper sense of the digits at the end
2285             of some FITS header fields. Sort your hash keys using "fits_field_cmp" and
2286             you will get (e.g.) your "TTYPE" fields in the correct order even if there
2287             are 140 of them.
2288              
2289             This is a standard kludgey perl comparison sub -- it uses the magical
2290             $a and $b variables, rather than normal argument passing.
2291              
2292             =cut
2293              
2294             sub fits_field_cmp {
2295 5 50   5 1 18 if( $a=~m/^(.*[^\d])(\d+)$/ ) {
2296 0         0 my ($an,$ad) = ($1,$2);
2297 0 0       0 if( $b=~m/^(.*[^\d])(\d+)$/ ) {
2298 0         0 my($bn,$bd) = ($1,$2);
2299            
2300 0 0       0 if($an eq $bn) {
2301 0         0 return $ad<=>$bd;
2302             }
2303             }
2304             }
2305 5         12 $a cmp $b;
2306             }
2307              
2308             =head2 _rows()
2309              
2310             =for ref
2311              
2312             Return the number of rows in a variable for table entry
2313              
2314             You feed in a PDLA or a list ref, and you get back the 0th dimension.
2315              
2316             =cut
2317              
2318             sub _rows {
2319 0     0   0 my $var = shift;
2320              
2321 0 0       0 return $var->dim(0) if( UNIVERSAL::isa($var,'PDLA') );
2322 0 0       0 return 1+$#$var if(ref $var eq 'ARRAY');
2323 0 0       0 return 1 unless(ref $var);
2324            
2325 0 0       0 print STDERR "Warning: _rows found an unacceptable ref. ".ref $var.". Ignoring...\n"
2326             if($PDLA::verbose);
2327            
2328 0         0 return undef;
2329             }
2330              
2331              
2332              
2333             =head2 _prep_table()
2334              
2335             =for ref
2336              
2337             Accept a hash ref containing a table, and return a header describing the table
2338             and a string to be written out as the table, or barf.
2339              
2340             You can indicate whether the table should be binary or ascii. The default
2341             is binary; it can be overridden by the "tbl" field of the hash (if present)
2342             or by parameter.
2343              
2344             =cut
2345              
2346             our %bintable_types = (
2347             'byte'=>['B',1],
2348             'short'=>['I',2],
2349             'ushort'=>['J',4, sub {return long shift;}],
2350             'long'=>['J',4],
2351             'longlong'=>['D', 8, sub {return double shift;}],
2352             'float'=>['E',4],
2353             'double'=>['D',8],
2354             # 'complex'=>['M',8] # Complex doubles are supported (actually, they aren't at the moment)
2355             );
2356              
2357              
2358             sub _prep_table {
2359 0     0   0 my ($hash,$tbl,$nosquish) = @_;
2360            
2361 0         0 my $ohash;
2362              
2363 0         0 my $hdr = $hash->{hdr};
2364            
2365 0         0 my $heap = "";
2366              
2367             # Make a local copy of the header.
2368 0         0 my $h = {};
2369 0 0       0 if(defined $hdr) {
2370 0         0 local $_;
2371 0         0 for (keys %$hdr) {$h->{$_} = $hdr->{$_}};
  0         0  
2372             }
2373 0         0 $hdr = $h;
2374              
2375 0 0       0 $tbl = $hash->{tbl} unless defined($tbl);
2376              
2377 0 0       0 barf "_prep_table called without a HASH reference as the first argument"
2378             unless ref $hash eq 'HASH';
2379              
2380             #####
2381             # Figure out how many columns are in the table
2382 0   0     0 my @colkeys = grep( ( !m/^(hdr|tbl)$/ and !m/^len_/ and defined $hash->{$_}),
2383             sort fits_field_cmp keys %$hash
2384             );
2385 0         0 my $cols = @colkeys;
2386              
2387 0 0       0 print "Table seems to have $cols columns...\n"
2388             if($PDLA::verbose);
2389            
2390             #####
2391             # Figure out how many rows are in the table, and store counts...
2392             #
2393 0         0 my $rows;
2394             my $rkey;
2395 0         0 for my $key(@colkeys) {
2396 0         0 my $r = _rows($hash->{$key});
2397 0 0 0     0 ($rows,$rkey) = ($r,$key) unless(defined($rows) && $rows != 1);
2398 0 0 0     0 if($r != $rows && $r != 1) {
2399 0         0 barf "_prep_table: inconsistent number of rows ($rkey: $rows vs. $key: $r)\n";
2400             }
2401             }
2402            
2403 0 0       0 print "Table seems to have $rows rows...\n"
2404             if($PDLA::verbose);
2405              
2406             #####
2407             # Squish and disambiguate column names
2408             #
2409 0         0 my %keysbyname;
2410             my %namesbykey;
2411              
2412 0 0       0 print "Renaming hash keys...\n"
2413             if($PDLA::debug);
2414              
2415 0         0 for my $key(@colkeys) {
2416 0         0 my $name = $key;
2417              
2418 0         0 $name =~ tr/[a-z]/[A-Z]/; # Uppercaseify (required by FITS standard)
2419 0         0 $name =~ s/\s+/_/g; # Remove whitespace (required by FITS standard)
2420            
2421 0 0       0 unless($nosquish) {
2422 0         0 $name =~ s/[^A-Z0-9_-]//g; # Squish (recommended by FITS standard)
2423             }
2424            
2425             ### Disambiguate...
2426 0 0       0 if(defined $ohash->{$name}) {
2427 0         0 my $iter = 1;
2428 0         0 my $name2;
2429 0         0 do { $name2 = $name."_".($iter++); }
2430 0         0 while(defined $ohash->{$name2});
2431 0         0 $name = $name2;
2432             }
2433              
2434 0         0 $ohash->{$name} = $hash->{$key};
2435 0         0 $keysbyname{$name} = $key;
2436 0         0 $namesbykey{$key} = $name;
2437              
2438 0 0 0     0 print "\tkey '$key'\t-->\tname '$name'\n"
      0        
2439             if($PDLA::debug || (($name ne $key) and $PDLA::verbose));
2440             }
2441              
2442              
2443             # The first element of colnames is ignored (since FITS starts the
2444             # count at 1)
2445             #
2446 0         0 my @colnames; # Names by number
2447             my %colnums; # Numbers by name
2448              
2449             ### Allocate any table columns that are already in the header...
2450 0         0 local $_;
2451 0         0 map { for my $x(1) { # [Shenanigans to make next work right]
  0         0  
2452 0 0       0 next unless m/^TTYPE(\d*)$/;
2453              
2454 0         0 my $num = $1;
2455            
2456 0 0 0     0 if($num > $cols || $num < 1) {
2457 0 0       0 print "Ignoring illegal column number $num ( should be in range 1..$cols )\n"
2458             if($PDLA::verbose);
2459 0         0 delete $hdr->{$_};
2460 0         0 next;
2461             }
2462              
2463 0         0 my $key = $hdr->{$_};
2464              
2465 0         0 my $name;
2466 0 0       0 unless( $name = $namesbykey{$key}) { # assignment
2467 0         0 $name = $key;
2468 0 0       0 unless( $key = $keysbyname{$key}) {
2469 0 0       0 print "Ignoring column $num in existing header (unknown name '$key')\n"
2470             if($PDLA::verbose);
2471 0         0 next;
2472             }
2473             }
2474              
2475 0         0 $colnames[$num] = $name;
2476 0         0 $colnums{$name} = $num;
2477             } } sort fits_field_cmp keys %$hdr;
2478              
2479             ### Allocate all other table columns in alphabetical order...
2480 0         0 my $i = 0;
2481 0         0 for my $k (@colkeys) {
2482 0         0 my $name = $namesbykey{$k};
2483              
2484 0 0       0 unless($colnums{$name}) {
2485 0         0 while($colnames[++$i]) { }
2486 0         0 $colnames[$i] = $name;
2487 0         0 $colnums{$name} = $i;
2488 0         0 } else { $i++; }
2489             }
2490            
2491 0 0 0     0 print "Assertion failed: i ($i) != colnums ($cols)\n"
2492             if($PDLA::debug && $i != $cols);
2493              
2494             print "colnames: " .
2495 0 0       0 join(",", map { $colnames[$_]; } (1..$cols) ) ."\n"
  0         0  
2496             if($PDLA::debug);
2497              
2498             ########
2499             # OK, now the columns are allocated -- spew out a header.
2500              
2501 0         0 my @converters = (); # Will fill up with conversion routines for each column
2502 0         0 my @field_len = (); # Will fill up with field lengths for each column
2503 0         0 my @internaltype = (); # Gets flag for PDLAhood
2504 0         0 my @fieldvars = (); # Gets refs to all the fields of the hash.
2505              
2506 0 0       0 if($tbl eq 'binary') {
    0          
2507 0         0 $hdr->{XTENSION} = 'BINTABLE';
2508 0         0 $hdr->{BITPIX} = 8;
2509 0         0 $hdr->{NAXIS} = 2;
2510             #$hdr->{NAXIS1} = undef; # Calculated below; inserted here as placeholder.
2511 0         0 $hdr->{NAXIS2} = $rows;
2512 0         0 $hdr->{PCOUNT} = 0; # Change this is variable-arrays are adopted
2513 0         0 $hdr->{GCOUNT} = 1;
2514 0         0 $hdr->{TFIELDS} = $cols;
2515              
2516             # Figure out data types, and accumulate a row length at the same time.
2517              
2518 0         0 my $rowlen = 0;
2519              
2520             # NOTE:
2521             # the conversion from ushort to long below is a hack to work
2522             # around the issue that otherwise perl treats it as a 2-byte
2523             # NOT 4-byte string on writing out, which leads to data corruption
2524             # Really ushort arrays should be written out using SCALE/ZERO
2525             # so that it can be written as an Int2 rather than Int4
2526             #
2527 0         0 for my $i(1..$cols) {
2528 0         0 $fieldvars[$i] = $hash->{$keysbyname{$colnames[$i]}};
2529 0         0 my $var = $fieldvars[$i];
2530              
2531 0         0 $hdr->{"TTYPE$i"} = $colnames[$i];
2532 0         0 my $tform;
2533            
2534             my $tstr;
2535 0         0 my $rpt;
2536 0         0 my $bytes;
2537            
2538 0 0       0 if( UNIVERSAL::isa($var,'PDLA') ) {
    0          
    0          
2539              
2540 0         0 $internaltype[$i] = 'P';
2541              
2542 0         0 my $t;
2543              
2544 0         0 my $dims = pdl($var->dims);
2545 0         0 ($t = $dims->slice("(0)")) .= 1;
2546 0         0 $rpt = $dims->prod;
2547              
2548             =pod
2549              
2550             =begin WHENCOMPLEXVALUESWORK
2551            
2552             if( UNIVERSAL::isa($var,'PDLA::Complex') ) {
2553             $rpt = $var->dim(1);
2554             $t = 'complex'
2555             } else {
2556             $t = type $var;
2557             }
2558              
2559             =end WHENCOMPLEXVALUESWORK
2560              
2561             =cut
2562              
2563 0 0       0 barf "Error: wfits() currently can not handle PDLA::Complex arrays (column $colnames[$i])\n"
2564             if UNIVERSAL::isa($var,'PDLA::Complex');
2565 0         0 $t = $var->type;
2566              
2567 0         0 $t = $bintable_types{$t};
2568            
2569 0 0       0 unless(defined($t)) {
2570 0 0       0 print "Warning: converting unknown type $t (column $colnames[$i]) to double...\n"
2571             if($PDLA::verbose);
2572 0         0 $t = $bintable_types{'double'};
2573             }
2574              
2575 0         0 ($tstr, $bytes, $converters[$i]) = @$t;
2576              
2577             } elsif( ref $var eq 'ARRAY' ) {
2578              
2579 0         0 $internaltype[$i] = 'A';
2580 0         0 $bytes = 1;
2581            
2582             # Got an array (of strings) -- find the longest element
2583 0         0 $rpt = 0;
2584 0         0 for(@$var) {
2585 0         0 my $l = length($_);
2586 0 0       0 $rpt = $l if($l>$rpt);
2587             }
2588 0         0 ($tstr, $bytes, $converters[$i]) = ('A',1,undef);
2589              
2590             } elsif( ref $var ) {
2591 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";
2592            
2593             } else { # Scalar
2594 0         0 $internaltype[$i] = 'A';
2595 0         0 ($tstr, $bytes, $converters[$i]) = ('A',1,undef);
2596 0         0 $rpt = length($var);
2597             }
2598              
2599              
2600             # Now check if it's a variable-length array and, if so, insert an
2601             # extra converter
2602 0         0 my $lname = "len_".$keysbyname{$colnames[$i]};
2603 0 0       0 if(exists $hash->{$lname}) {
2604 0         0 my $lengths = $hash->{$lname};
2605              
2606             # Variable length array - add extra handling logic.
2607              
2608             # First, check we're legit
2609 0 0 0     0 if( !UNIVERSAL::isa($var, 'PDLA') ||
      0        
      0        
      0        
2610             $var->ndims != 2 ||
2611             !UNIVERSAL::isa($lengths,'PDLA') ||
2612             $lengths->ndims != 1 ||
2613             $lengths->dim(0) != $var->dim(0)
2614             ) {
2615 0         0 die <<'FOO';
2616             wfits(): you specified a 'len_$keysbyname{$colnames[$i]}' field in
2617             your binary table output hash, indicating a variable-length array for
2618             each row of the output table, but I'm having trouble interpreting it.
2619             Either your source column isn't a 2-D PDLA, or your length column isn't
2620             a 1-D PDLA, or the two lengths don't match. I give up.
2621             FOO
2622             }
2623              
2624             # The definition below wraps around the existing converter,
2625             # dumping the variable to the heap and returning the length
2626             # and index of the data for the current row as a PDLA LONG.
2627             # This does the Right Thing below in the write loop, with
2628             # the side effect of putting the data into the heap.
2629             #
2630             # The main downside here is that the heap gets copied
2631             # multiple times as we accumulate it, since we are using
2632             # string concatenation to add onto it. It might be better
2633             # to preallocate a large heap, but I'm too lazy to figure
2634             # out how to do that.
2635 0         0 my $csub = $converters[$i];
2636             $converters[$i] = sub {
2637 0     0   0 my $var = shift;
2638 0         0 my $row = shift;
2639 0         0 my $col = shift;
2640            
2641 0         0 my $len = $hash->{"len_".$keysbyname{$colnames[$i]}};
2642 0         0 my $l;
2643 0 0       0 if(ref $len eq 'ARRAY') {
    0          
    0          
2644 0         0 $l = $len->[$row];
2645             } elsif( UNIVERSAL::isa($len,'PDLA') ) {
2646 0         0 $l = $len->dice_axis(0,$row);
2647             } elsif( ref $len ) {
2648 0         0 die "wfits: Couldn't understand length spec 'len_".$keysbyname{$colnames[$i]}."' in bintable output (length spec must be a PDLA or array ref).\n";
2649             } else {
2650 0         0 $l = $len;
2651             }
2652            
2653             # The standard says we should give a zero-offset
2654             # pointer if the current row is zero-length; hence
2655             # the ternary operator.
2656 0 0       0 my $ret = pdl( $l, $l ? length($heap) : 0)->long;
2657              
2658              
2659 0 0       0 if($l) {
2660             # This echoes the normal-table swap and accumulation
2661             # stuff below, except we're accumulating into the heap.
2662 0 0       0 my $tmp = $csub ? &$csub($var, $row, $col) : $var;
2663 0         0 $tmp = $tmp->slice("0:".($l-1))->sever;
2664            
2665 0 0       0 if(!isbigendian()) {
2666 0 0       0 bswap2($tmp) if($tmp->get_datatype == $PDLA_S);
2667 0 0 0     0 bswap4($tmp) if($tmp->get_datatype == $PDLA_L ||
2668             $tmp->get_datatype == $PDLA_F);
2669 0 0       0 bswap8($tmp) if($tmp->get_datatype == $PDLA_D);
2670             }
2671 0         0 my $t = $tmp->get_dataref;
2672 0         0 $heap .= $$t;
2673             }
2674              
2675 0         0 return $ret;
2676 0         0 };
2677              
2678             # Having defined the conversion routine, now modify tstr to make this a heap-array
2679             # reference.
2680 0         0 $tstr = sprintf("P%s(%d)",$tstr, $hash->{"len_".$keysbyname{$colnames[$i]}}->max );
2681 0         0 $rpt = 1;
2682 0         0 $bytes = 8; # two longints per row in the main table.
2683             }
2684              
2685            
2686 0         0 $hdr->{"TFORM$i"} = "$rpt$tstr";
2687              
2688 0 0 0     0 if(UNIVERSAL::isa($var, 'PDLA') and $var->ndims > 1) {
2689 0         0 $hdr->{"TDIM$i"} = "(".join(",",$var->slice("(0)")->dims).")";
2690             }
2691              
2692 0         0 $rowlen += ($field_len[$i] = $rpt * $bytes);
2693             }
2694            
2695 0         0 $hdr->{NAXIS1} = $rowlen;
2696            
2697             ## Now accumulate the binary table
2698              
2699 0         0 my $table = "";
2700            
2701 0         0 for my $r(0..$rows-1) {
2702 0         0 my $row = "";
2703 0         0 for my $c(1..$cols) {
2704 0         0 my $tmp;
2705 0         0 my $x = $fieldvars[$c];
2706            
2707 0 0       0 if($internaltype[$c] eq 'P') { # PDLA handling
2708             $tmp = $converters[$c]
2709 0 0       0 ? &{$converters[$c]}($x->slice("$r")->flat->sever, $r, $c)
  0         0  
2710             : $x->slice("$r")->flat->sever ;
2711              
2712             ## This would go faster if moved outside the loop but I'm too
2713             ## lazy to do it Right just now. Perhaps after it actually works.
2714             ##
2715 0 0       0 if(!isbigendian()) {
2716 0 0       0 bswap2($tmp) if($tmp->get_datatype == $PDLA_S);
2717 0 0 0     0 bswap4($tmp) if($tmp->get_datatype == $PDLA_L ||
2718             $tmp->get_datatype == $PDLA_F);
2719 0 0       0 bswap8($tmp) if($tmp->get_datatype == $PDLA_D);
2720             }
2721              
2722 0         0 my $t = $tmp->get_dataref;
2723 0         0 $tmp = $$t;
2724             } else { # Only other case is ASCII just now...
2725 0 0       0 $tmp = ( ref $x eq 'ARRAY' ) ? # Switch on array or string
    0          
2726             ( $#$x == 0 ? $x->[0] : $x->[$r] ) # Thread arrays as needed
2727             : $x;
2728            
2729 0         0 $tmp .= " " x ($field_len[$c] - length($tmp));
2730             }
2731            
2732             # Now $tmp contains the bytes to be written out...
2733             #
2734 0         0 $row .= substr($tmp,0,$field_len[$c]);
2735             } # for: $c
2736 0         0 $table .= $row;
2737             } # for: $r
2738              
2739 0         0 my $table_size = $rowlen * $rows;
2740 0 0       0 if( (length $table) != $table_size ) {
2741 0         0 print "Warning: Table length is ".(length $table)."; expected $table_size\n";
2742             }
2743              
2744 0         0 return ($hdr,$table, $heap);
2745            
2746             } elsif($tbl eq 'ascii') {
2747 0         0 barf "ASCII tables not yet supported...\n";
2748             } else {
2749 0         0 barf "unknown table type '$tbl' -- giving up.";
2750             }
2751             }
2752              
2753             # the header fill value can be blanks, but the data fill value must
2754             # be zeroes in non-ASCII tables
2755             #
2756             sub _print_to_fits ($$$) {
2757 1     1   3 my $fh = shift;
2758 1         2 my $data = shift;
2759 1         2 my $blank = shift;
2760              
2761 1         4 my $len = ((length $data) - 1) % 2880 + 1;
2762 1         19 $fh->print( $data . ($blank x (2880-$len)) );
2763             }
2764            
2765             {
2766             my $ctr = 0;
2767              
2768 0     0 0 0 sub reset_hdr_ctr() { $ctr = 0; }
2769              
2770             sub add_hdr_item ($$$$;$) {
2771 0     0 0 0 my ( $hdr, $key, $value, $type, $comment ) = @_;
2772 0 0       0 $type = uc($type) if defined $type;
2773 0         0 my $item = Astro::FITS::Header::Item->new( Keyword=>$key,
2774             Value=>$value,
2775             Type=>$type );
2776 0 0       0 $item->comment( $comment ) if defined $comment;
2777 0         0 $hdr->replace( $ctr++, $item );
2778             };
2779             }
2780              
2781             ##############################
2782             #
2783             # _wfits_table -- given a hash ref, try to write it out as a
2784             # table extension. The file FITS should be open when you call it.
2785             # Most of the work of creating the extension header, and all of
2786             # the work of creating the table, is handled by _prep_table().
2787             #
2788             # NOTE:
2789             # can not think of a sensible name for the extension so calling
2790             # it TABLE for now
2791             #
2792             sub _wfits_table ($$$) {
2793 0     0   0 my $fh = shift;
2794 0         0 my $hash = shift;
2795 0         0 my $tbl = shift;
2796            
2797 0 0       0 barf "FITS BINTABLES are not supported without the Astro::FITS::Header module.\nGet it from www.cpan.org.\n"
2798             unless($PDLA::Astro_FITS_Header);
2799              
2800 0         0 my ($hdr,$table, $heap) = _prep_table($hash,$tbl,0);
2801 0 0       0 $heap="" unless defined($heap);
2802              
2803             # Copy the prepared fields into the extension header.
2804 0         0 tie my %newhdr,'Astro::FITS::Header',my $h = Astro::FITS::Header->new;
2805            
2806 0         0 reset_hdr_ctr();
2807 0 0       0 add_hdr_item $h, "XTENSION", ($tbl eq 'ascii'?"TABLE":"BINTABLE"), 'string', "from perl hash";
2808 0         0 add_hdr_item $h, "BITPIX", $hdr->{BITPIX}, 'int';
2809 0         0 add_hdr_item $h, "NAXIS", 2, 'int';
2810 0         0 add_hdr_item $h, "NAXIS1", $hdr->{NAXIS1}, 'int', 'Bytes per row';
2811 0         0 add_hdr_item $h, "NAXIS2", $hdr->{NAXIS2}, 'int', 'Number of rows';
2812 0 0       0 add_hdr_item $h, "PCOUNT", length($heap), 'int', ($tbl eq 'ascii' ? undef : "No heap") ;
2813 0 0       0 add_hdr_item $h, "THEAP", "0", "(No gap before heap)" if(length($heap));
2814 0         0 add_hdr_item $h, "GCOUNT", 1, 'int';
2815 0         0 add_hdr_item $h, "TFIELDS", $hdr->{TFIELDS},'int';
2816 0         0 add_hdr_item $h, "HDUNAME", "TABLE", 'string';
2817              
2818 0         0 for my $field( sort fits_field_cmp keys %$hdr ) {
2819 0 0 0     0 next if( defined $newhdr{$field} or $field =~ m/^end|simple|xtension$/i);
2820             my $type = (UNIVERSAL::isa($hdr->{field},'PDLA') ?
2821             $hdr->{$field}->type :
2822 0 0       0 ((($hdr->{$field})=~ m/^[tf]$/i) ?
    0          
2823             'logical' :
2824             undef ## 'string' seems to have a bug - 'undef' works OK
2825             ));
2826              
2827 0         0 add_hdr_item $h, $field, $hdr->{$field}, $type, $hdr->{"${field}_COMMENT"};
2828             }
2829              
2830 0         0 add_hdr_item $h, "END", undef, 'undef';
2831              
2832 0         0 $hdr = join("",$h->cards);
2833 0         0 _print_to_fits( $fh, $hdr, " " );
2834 0         0 _print_to_fits( $fh, $table.$heap, "\0" ); # use " " if it is an ASCII table
2835             # Add heap dump
2836             }
2837              
2838             sub _wfits_nullhdu ($) {
2839 1     1   3 my $fh = shift;
2840 1 50       4 if($Astro::FITS::Header) {
2841 0         0 my $h = Astro::FITS::Header->new();
2842            
2843 0         0 reset_hdr_ctr();
2844 0         0 add_hdr_item $h, "SIMPLE", "T", 'logical', "Null HDU (no data, only extensions)";
2845 0         0 add_hdr_item $h, "BITPIX", -32, 'int', "Needed to make fverify happy";
2846 0         0 add_hdr_item $h, "NAXIS", 0, 'int';
2847 0         0 add_hdr_item $h, "EXTEND", "T", 'logical', "File contains extensions";
2848 0         0 add_hdr_item $h, "COMMENT", "", "comment",
2849             " File written by perl (PDLA::IO::FITS::wfits)";
2850             #
2851             # The following seems to cause a problem so removing for now (I don't
2852             # believe it is required, but may be useful for people who aren't
2853             # FITS connoisseurs). It could also be down to a version issue in
2854             # Astro::FITS::Header since it worked on linux with a newer version
2855             # than on Solaris with an older version of the header module)
2856             #
2857             ## add_hdr_item $h, "COMMENT", "", "comment",
2858             ## " FITS (Flexible Image Transport System) format is defined in 'Astronomy";
2859             ## add_hdr_item $h, "COMMENT", "", "comment",
2860             ## " and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H";
2861 0         0 add_hdr_item $h, "HDUNAME", "PRIMARY", 'string';
2862 0         0 add_hdr_item $h, "END", undef, 'undef';
2863            
2864 0         0 my $hdr = join("",$h->cards);
2865 0         0 _print_to_fits( $fh, $hdr, " " );
2866             } else {
2867 1         6 _print_to_fits( $fh,
2868             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 (PDLA::IO::FITS::wfits) legacy code. COMMENT For best results, install Astro::FITS::Header. HDUNAME = 'PRIMARY ' END +,
2869             " ");
2870             }
2871             }
2872              
2873            
2874             1;
2875