Re: FITS format files

[jfc]

Laura Richter <llrichter <at> gmail.com> writes:
> I would like to be able to read FITS (Flexible Image Transport System)
> format files into octave.

Hi Laura,

you will find public domain code lying around the web to do just this if
you google for FITS+Matlab.  I have found that these m.files need serious
editing to work in octave, though.  I gave up up trying to use them (also
because they implement only the most basic features of FITS) and went
another route: write oct files which wrap around the cfitsio library.

If you want to go this way, you need to install the cfitsio
library available from NASA at http://heasarc.nasa.gov/fitsio/fitsio.html
and then "mkoctfile" an oct file.  I append below the code
to read images from a fits file. 

By the way: this is naive code, mostly resulting from trial-and-error.
Comments from C++ litterati would be welcome.

Cheers, JF.



========= beginning of code  ==============
#include <octave/oct.h>

extern "C" { 
#include "fitsio.h"
}

/*
  read a "standard" (?) image from hdu 1 of a fits file
  Author: JF Cardoso, nov 2006.


  TODO:  We are duplicating the data, aren't we ? 
         We should clarify this octave_value thing...

*/


DEFUN_DLD (read_fits_image, args, , \
" matrix = read_fits_image (fits_file_name) reads an image from hdu1.\n.")

{ octave_value fitsimage ;
  
  if ( args.length() !=  1 ) { 
     error("Needs a single argument: a file name\n");
     return fitsimage = -1 ; }

  // args(0) = file name
  if ( !args(0).is_string() ) 
     { error("Error: argument is not string type");
       return fitsimage = -1 ;  }
  std::string infile = args(0).string_value ();


  /* Initializations */
  int status = 0;
  fitsfile *fptr;
  

  /* Open FITS file and position to first HDU containing an image*/
  if ( fits_open_image(&fptr, infile.c_str(), READONLY, &status) ) 
    {  
      fprintf(stderr, "%s (line %d): File %s not found or image not found!\n",
              __FILE__, __LINE__, infile.c_str() );     
      fits_report_error(stderr, status);
      return fitsimage = -1 ;  
    }  
 
  int  naxis ; // Number of dimensions, should be 2
  long naxes [2] ;
  int maxdim = 2 ; // ????
  int bitpix ;
  if ( fits_get_img_param(fptr, maxdim, &bitpix, &naxis, naxes, &status) )
    {
      fprintf(stderr,
 "%s (line %d): Problem reading information about current HDU from file %s\n",
              __FILE__, __LINE__, infile.c_str() );     
      fits_report_error(stderr, status);
      return fitsimage = -1 ;  
      }  
  
  if ( naxis != 2) {
      fprintf(stderr, "Are you insane?  We expect an image here.\n");
      return fitsimage = -1 ;  
  }

  long xsize = naxes[0] ;
  long ysize = naxes[1] ;
  long npix = xsize * ysize ;

  /* Read the array */    
  Matrix mapvec(xsize, ysize ); 
  double *map = mapvec.fortran_vec() ;

  long fpixel[2] ;
  fpixel[0] = 1 ;  // To start at
  fpixel[1] = 1 ;  // first pixel

  int  anynul;
  if ( fits_read_pix
       (fptr, 
        TDOUBLE,  // This is the type of the C array to be filled,
                   // NOT the type of the FITS table
        fpixel,   // indices of the first pixels to read in each dimension
        npix,     // number of pixels
        NULL,     // should be the POINTER to the null value. 
                  // If NULL, no checks are made (am I understanding myself?)
        map, 
        &anynul,
        &status ) ) 
    { 
      fprintf(stderr, "%s (line %d): Cannot read image %i.\n", 
     __FILE__,
     __LINE__ ); 
      fits_report_error(stderr, status);
      return fitsimage = -1 ; 
    }
     
  fprintf(stderr, "Done reading \n");
  
  if ( fits_close_file(fptr, &status) ) {  /* Close the file */
    fprintf(stderr,
"%s (line %d): Cannot close the FITS file!  This is really unexpected...\n",
  __FILE__, __LINE__ );     
    fits_report_error(stderr, status);
    return fitsimage = -1 ;
  }
  
  return fitsimage =  mapvec ; // What does this do exactly ???
}
========= end of code