Re: Flexible Image Transport System (FITS) file format support

[Thomas D. Dean]

This octave code reads all the fits files I have.  This code is most
likely not the best design, but, ...

The Image Magick, linux application display will display fits images.

## Copyright (C) 2010 Thomas Dean <address@hidden>
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.
##
###############################################################
##
## several functions to do with reading fits files
##
## load this by
## octave:1> cd to the correct directory
## octave:2> readfits
## octave:3> fits=readfitsfile('rosat_pspc_rdf2_3_im1.fits');
##
## the structure fits contains
## fits.idx
## fits.hdr{1:fits.idx}
## fits.bin{1:fits.idx}
##
## for the rosat image file above,
## image(fits.bin{1}/64)
##
##
1; ## make sure octave does not think this is a function file
##
###############################################################
## read fits header
function [hdr,numread]=readfitshdr(fid)
  ## [hdr,numread]=readfitshdr(fid)
  ## read a fits header from the file
  hdr=[];
  numread=0;
  while (sum(strcmp(hdr,'END'))==0)
    [tmp_hdr,cnt]=fread(fid,[80,36],"uchar",0);
    numread = numread + cnt;
    hdr=[hdr;cellstr(tmp_hdr')];
  end;
endfunction;
##
###############################################################
## what are the units of the binary part of the file
function [prec,readmult]=whatprec(bitpix)
  switch bitpix
    case -64
      prec="float32";
      readmult = 4;
    case -32
      prec="float16";
      readmult=2;
    case 8
      prec="uint8";
      readmult=1;
    case 16
      prec="uint16";
      readmult=2;
    case 32
      prec="uint32";
      readmult=4;
    otherwise
      error("bitpix not in { -64 -32 8 16 32 64 }");
  endswitch
endfunction;
##
###############################################################
## read fits file
function [fits]=readfitsfile(fname)
  ## read a fits file
  ##   [hdr,hcnt,bin,bcnt]=readfits(fname)
  ##     fname   - a file name 'rosat_pspc_rdf2_3_bk1.fits'
  ##     binsize - size of the binary array [80,200]
  ##     ndr     - an array containing the header of the fits file
  ##     hcnt    - number of bytes in the fits file
  ##     bin     - binary data from the fits file header
  ##     bcnt    - number of elements in the binary array
  ##
  ## work on header
  ## [bitpix,c]=sscanf(hdr{strncmp(hdr,'BITPIX  =',9)},'BITPIX  =%d');
  ## [naxis,c] =sscanf(hdr{strncmp(hdr,'NAXIS ',6)},'NAXIS   =%d');
  ## [naxis1,c] =sscanf(hdr{strncmp(hdr,'NAXIS1',6)},'NAXIS1  =%d');
  ## [naxis2,c] =sscanf(hdr{strncmp(hdr,'NAXIS2',6)},'NAXIS2  =%d');
  ##
  ## use image(bin) or image(bin/32), etc.
  numread=0;
  fitsidx=0;
  [fid, msg]=fopen(fname,"r");
  rc = fseek(fid,0,SEEK_END);
  filesize = ftell(fid);
  rc = fseek(fid,0,SEEK_SET);
  printf("The file has %d bytes\n",filesize);
  while (numread < filesize)
    ## read a header
    [hdr,cnt]=readfitshdr(fid);
    fitsidx = fitsidx+1;
    fits.hdr{fitsidx}=hdr;
    numread = numread + cnt;
    ## check number of bits per pixel
    [bitpix,c]=sscanf(hdr{strncmp(hdr,"BITPIX  =",9)},"BITPIX  =%d");
    [prec,readmult]=whatprec(bitpix);
    ## check the number of axis, 0, 1, 2, 3, ...
    ## for now, only support 3?????
    [naxis,c] =sscanf(hdr{strncmp(hdr,"NAXIS ",6)},"NAXIS   =%d");
    ## check for 1, 2, or 3 axis
    if (naxis >= 1)
      [naxis1,c] =sscanf(hdr{strncmp(hdr,"NAXIS1",6)},"NAXIS1  =%d");
      if (naxis >= 2)
        [naxis2,c] =sscanf(hdr{strncmp(hdr,"NAXIS2",6)},"NAXIS2  =%d");
      endif;
      if (naxis == 3)
        [naxis3,c] =sscanf(hdr{strncmp(hdr,"NAXIS3",6)},"NAXIS3  =%d");
      endif;
      ## read the binary part of the file
      [bin,cnt]=fread(fid,[naxis1,naxis2],prec,0);
      ## read any remaining data to fill out a integral number of
blocks,
      ## 2880 bytes per block
      numread = numread + readmult*cnt;
      printf("So far, have read %d bytes\n",numread);
      numtoread = ceil(numread/2880)*2880-numread;
      printf("Need to dummy read %d bytes\n",numtoread);
      [dummy,cnt]=fread(fid,numtoread);
      #### copy results to output
      fits.bin{fitsidx}=bin;
      numread = numread + cnt;
    endif;
    printf("Now, have read %d bytes or %d blocks
\n",numread,numread/2880);
  end
  ## close the input file
  fclose(fid);
  fits.idx = fitsidx;
endfunction;