Reading FITS images into Mathematica

Roberto G. Abraham

Practical image processing now seems possible with Version 4 of Mathematica (which now includes packed arrays, and has vastly improved numerical performance and storage characteristics). At least one commercial quality software add-on package for image processing in Mathematica will soon become available.  Mathematica may well prove useful for analysing astronomical CCD data.  However the vast majority of professional astronomical images are distributed in FITS (Flexible Image Transport System) format, as opposed to more common image formats such as GIF or JPEG.  Mathematica cannot directly handle FITS images with the Import[] function.    I have therefore written a package which allows FITS data to be read into Mathematica.   This package uses MathLink to call the publicly available NASA CFITSIO library.

FITS is a data format based upon the notion of a "header data unit" (HDU).  Each set of data (which can be a one-dimensional list, a two dimensional image, a multi-dimensional data cube, or a table of data in binary or ascii format) is preceeded by a number of ASCII header lines giving useful information (such as the dimensionality of the data, the origin of the data, format for the data, etc.). A single FITS file may contain several HDUs. The present version of the Mathematica FITS reader is designed to work only with two-dimensional images.  More information on the FITS standard is given in Appendix II.

Setup

I am using Mathematica on Windows 98/NT systems, and have thus pre-compiled the FITS reader for use with Windows. If you are using a non-Windows  operating system you will need to compile the MathLink source files and link these against the freely available CFITSIO library.  The CFITSIO library exists for Windows, many flavours of Unix, VMS, and Macintosh.  I expect building the MathLink code and linking it against CFITSIO is straightforward for almost any architecture. The source code for the MathLink and glue routines is given in Appendix I.  Information on downloading the CFITSIO library is given in Appendix II.

Assuming you want to use the pre-compiled Windows binary, all you need to do is modify one line in the "Fits.m" package file to indicate the location of the executible file on your system.  Open tbe "Fits.m" package file, and change the following line to correspond to the correct location of the fitsread.exe binary on your machine:

Install["C:\\Mathematica FITS Reader\\fitsread.exe"];

Usage

Functions

The Fits.m package exports 3 functions:

FitsRead[] - Reads a FITS file into a Mathematica Real array.
FitsHeader[] - Extracts a FITS header as a newline-delimited string.
FitsHeaderToRules[] - Converts a FITS header string into a set of rules for convenient manipulation.

Reading images

Fits files are read as follows:

image = FitsRead[filename]

An interesting thing to note is that the filename specification can include a template specifying a subset of the image. This is very useful for quickly reading small subsets of enormous images (the FitsRead[] routine jumps to the subset in the data file and does not waste memory). The template format is as follows:

<filename>[xlow:xhigh,ylow:yhigh]

In other words, to read a 100x100 "postage stamp" subset from the files "foo.fits" beginning at pixel x=50,y=100 you would use the following syntax:

im = FitsRead["foo.fits[50:150,100:200"];

Multiple header data units can also be handled. The syntax for this is:

<filename>[hdu_number][xlow:xhigh,ylow:yhigh]

For more information on this capability (and other features) see the documentation for the CFITSIO library (linked to in Appendix II).

Examining and manipulating image headers

The header from a FITS file can be examined as follows:

header = FitsHeader[filename]

This extracts the header as a long string, broken into lines known as FITS "cards". To extract individual header entries or to manipulate the header information, it's more convenient to convert this into a set of rules:

rules = FitsHeader[filename] // FitsHeaderToRules

Note that this routine will change some of FITS keywords from "illegal" Mathematica symbols to ones that are acceptable to Mathematica. For example, keywords with underscores and minus signs will be transformed into corresponding symbols with the illegal characters omitted.  In very pathological cases this may lead to problems. For example, if your FITS file has keywords "FOO_BAR", "FOO-BAR" and "FOOBAR", then you will end up with multiple rules for the FOOBAR symbol. 

This routine also attempts to convert keyword values from the Fortran-esque notation often used in FITS files to a valid Mathematica number. For example, the routine will convert from -1.2E11 to  -1.2 ×(10^11).  I have tested this conversion carefully and I believe the conversion works for all valid Fortran-formatted numbers, but please check this on your own data as I would not be surprised if very pathological cases do not convert properly.

An example

<< "C:\\Mathematica FITS Reader\\Fits.m"

Read in and display an entire FITS image. Reading a 313x313 image takes around 7 seconds on a rather slow P166 machine.

im = FitsRead["C:\\Windows\\Desktop\\N2715_lr.fit"];

ListDensityPlot[im, Mesh -> False, 
      PlotRange -> {1000, 3000}]

Read in and display a 100x100 subset of this FITS image. Small subsets of images are read in very quickly... this example takes under a second on my slowest machine.

\(\(im = 
        FitsRead["C:\Windows\Desktop\N2715_lr.fit[115:215,105:205]"];\)\)

ListDensityPlot[im, Mesh -> False, 
      PlotRange -> {1000, 3000}]

Read in the FITS header. Note that only FITS cards with keywords are converted (eg. comments are omitted).

\(\(\(h\)\(=\)\(FitsHeader["C:\Windows\Desktop\N2715_lr.fit"]\)\(\ \)\)\)

Let's manipulate some header information by converting the header to a list of rules. Note that FITS keywords with underscores and minus signs are convted to valid Mathematica symbol names by dropping the offending characters.

hr =  FitsHeaderToRules[h]

Extract the right ascension, declination, time, airmass, and datamax keywords:

{RA, DEC, DATE,  TIME, AIRMASS, DATAMAX}  /.  hr

Appendix I: MathLink source code

The following routines need to linked against the freely available CFITSIO library. See Appendix II for information on downloading this library (which is, incidentally, extremely well designed and trivial to build on every machine I've ever tried it on).

Routines calling MathLink

All explicit MathLink calls are isolated to the following routines.

fitsread.c

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include "mathlink.h"
#include "fitsread.h"

/*-----------------*/
/* return FITS 2D image */

void fitsread( char *name)
{
	double *arr;
	long nelem;
	long dims[2];
	long depth = 2;
	int nr,nc;
	int status;             /* non-zero indicates a CFITSIO error */ 
	long mlerr;             /* 0 indicates an error for MathLink */
	char *better_name;      /* "\\" converted to "/" in directory name */

	status = 0;
	mlerr = 1;

	/* Deal with Mathematica's "\\" directory separator for Windows */
	better_name = convert_to_single_backslash(name);

	arr = readimage(better_name, &nr, &nc, &status); 
	if (check_for_errors(status,mlerr)){
		free(better_name);
		return;
	}

	/* Send the array back to Mathematica */
	dims[0] = (long) nc;
	dims[1] = (long) nr;
	nelem = (long) nr*nc;
	mlerr = MLPutRealArray(stdlink, arr, dims, NULL, depth);
	free(arr);
	free(better_name);
	check_for_errors(status,mlerr);
	return;
}


/*--------------------------------------------------------------------------*/
/* return FITS file header */

void fitsheader( char *name)
{
	char *header;
	int status;             /* non-zero indicates a CFITSIO error */ 
	long mlerr;             /* 0 indicates an error for MathLink */
	char *better_name;      /* "\\" converted to "/" in directory name */

	status = 0;
	mlerr = 1;

	/* Deal with Mathematica's "\\" directory separator for Windows */
	better_name = convert_to_single_backslash(name);

	header = readheader(better_name, &status); 
	if (check_for_errors(status,mlerr)){
		free(better_name);
		return;
	}

	/* Send the string back to Mathematica */
	mlerr = MLPutString(stdlink, header);
	free(header);
	free(better_name);
	check_for_errors(status,mlerr);
	return;
}



/*--------------------------------------------------------------------------*/
/* extract filename from a string with region subset suffix attatched */

char *get_filename(char s[])
{
	int i;
	char *new_s;

	new_s = malloc(MAX_STRING);
	i=0;
	while (s[i] != '\0' && s[i] != '[' && i < MAX_STRING){
			new_s[i]=s[i];
		++i;
	}
	new_s[i]='\0';
	return(new_s);
};


/*--------------------------------------------------------------------------*/
/* convert "\\" to "\" in a string. */

char* convert_to_single_backslash(char s[])
{
	int i,j;
	char *new_s;

	new_s = malloc(MAX_STRING);
	i=0;
	j=0;
	while (s[i] != '\0'){
		if ((s[i] == 92) && (s[i+1] == 92) && i < MAX_STRING)
			++i;
		new_s[j]=s[i];
		++i;
		++j;
	}
	new_s[j]='\0';
	return(new_s);
};



/*--------------------------------------------------------------------------*/
/* General error routine that reports both CFITSIO and MathLink errors and
tries to recover from the latter. If an error was found return 1. */

int check_for_errors(int status, long mlerr)
{
	char err_msg[100];
	int found;
	found = 0;
	if (status>0){
		sprintf(err_msg, "%s\"%.76s\"%s",
			"Message[fits::cfitsio,",errormessage(status),"]");
		found = 1;
	}
	if (status==115){
		sprintf(err_msg, "%s",
			"Message[fits::nofile]");
	}
	if(mlerr==0){
		sprintf(err_msg, "%s\"%.76s\"%s",
			"Message[fits::mlink,",MLErrorMessage(stdlink),"]");
		found = 1;
	}
	if (found) {
		MLClearError(stdlink);
		MLNewPacket(stdlink);
		MLEvaluate(stdlink, err_msg);         
		MLNextPacket(stdlink);
		MLNewPacket(stdlink);
		MLPutSymbol(stdlink, "$Failed");
	}
	return(found);
}



/*--------------------------------------------------------------------------*/
#if MACINTOSH_MATHLINK

int main( int argc, char* argv[])
{
	/* Due to a bug in some standard C libraries that have shipped with
	 * MPW, zero is passed to MLMain below.  (If you build this program
	 * as an MPW tool, you can change the zero to argc.)
	 */
	argc = argc; /* suppress warning */
	return MLMain( 0, argv);
}

#elif WINDOWS_MATHLINK

#if __BORLANDC__
#pragma argsused
#endif

int PASCAL WinMain( HINSTANCE hinstCurrent, HINSTANCE hinstPrevious, LPSTR lpszCmdLine, int nCmdShow)
{
	char  buff[512];
	char FAR * buff_start = buff;
	char FAR * argv[32];
	char FAR * FAR * argv_end = argv + 32;

	hinstPrevious = hinstPrevious; /* suppress warning */

	if( !MLInitializeIcon( hinstCurrent, nCmdShow)) return 1;
	MLScanString( argv, &argv_end, &lpszCmdLine, &buff_start);
	return MLMain( argv_end - argv, argv);
}

#else

int main(argc, argv)
	int argc; char* argv[];
{
	return MLMain(argc, argv);
}

#endif

fitsread.h

double *readimage(char *filename, int *nx, int *ny, int *status);
char *readheader(char *filename, int *error_status);
char *errormessage(int status);
char *convert_to_single_backslash(char *s);
char *get_filename(char *s);
int check_for_errors(int status, long mlerr);


#define MAX_STRING 1000
#define CFITSIO_ERROR 100
#define MATHLINK_ERROR 200

fitsread.tm

void fitsread P(( char *));

:Begin:
:Function:       fitsread
:Pattern:        FitsRead[file_String]
:Arguments:      { file }
:ArgumentTypes:  { String }
:ReturnType:     Manual
:End:


void fitsheader P(( char *));

:Begin:
:Function:       fitsheader
:Pattern:        FitsHeader[file_String]
:Arguments:      { file }
:ArgumentTypes:  { String }
:ReturnType:     Manual
:End:

Routines calling the CFITSIO library

For convenience, all CFITSIO calls are isolated in these two routines

readimage.c

#include <string.h>
#include <stdlib.h>
#include <malloc.h>

#include "fitsio.h"

double *readimage(char *filename, int *nx, int *ny, int *error_status)

{
	double *arr;
    fitsfile *fptr;       /* pointer to the FITS file, defined in fitsio.h */
    long naxes[2], npixels, firstelem;
    int nfound, anynull, status;
    double nullval;
	char filestring[FLEN_FILENAME];
	strncpy(filestring,filename,FLEN_FILENAME);

    status = 0;
	*nx = 0;
	*ny = 0;

    if(fits_open_file(&fptr, filestring, READONLY, &status)){
		fits_close_file(fptr, &status);
		*error_status = status;		
		return (arr);
	}
    if(fits_read_keys_lng(fptr, "NAXIS", 1, 2, naxes, &nfound, &status)){
		fits_close_file(fptr, &status);
		*error_status = status;
		return(arr);
	}

	firstelem = 1;
    npixels  = naxes[0] * naxes[1];     /* number of pixels in the image */
    nullval  = 0;         /* don't check for null values in the image */
	arr = (double *) malloc(sizeof(double)*npixels);
	if(fits_read_img(fptr,TDOUBLE,firstelem,npixels,&nullval,arr,&anynull,&status)){
		fits_close_file(fptr, &status);
		*error_status = status;
		return(arr);
	}

    if(fits_close_file(fptr, &status)){
		*error_status = status;
		return(arr);
	}

	*nx = naxes[0];
	*ny = naxes[1];
	*error_status = status;

    return(arr);
};


char *errormessage(int status)
{
	char *message;
	message=malloc(FLEN_ERRMSG);
	fits_get_errstatus(status,message);
	return(message);
}

readheader.c

#include <string.h>
#include <stdlib.h>
#include <malloc.h>

#include "fitsio.h"

char *readheader(char *filename, int *error_status)
{
    fitsfile *fptr;       /* pointer to the FITS file, defined in fitsio.h */
    int status;
	int nkey,nblank;
	char filestring[FLEN_FILENAME];
	char *header;
	char card[80];
	int i,j;
	int index;
	
	strncpy(filestring,filename,FLEN_FILENAME);
    status = 0;
	
    if(fits_open_file(&fptr, filestring, READONLY, &status)){
		fits_close_file(fptr, &status);
		*error_status = status;		
		return (header);
	}
	
    if(fits_get_hdrspace(fptr, &nkey, &nblank, &status)){
		fits_close_file(fptr, &status);
		*error_status = status;
		return(header);
	}
	
	header = malloc(nkey*81+1);   /* A FITS line is 80 characters by definition. Add 1 for newlines 
	                                 and the final '\0' on the last line. */
	index = 0;
	for(i=0;i<nkey;i++){
		if(fits_read_record(fptr, i+1, card, &status)){
			fits_close_file(fptr, &status);
			*error_status = status;		
			return (header);
		}
		
		for(j=0;j<80;j++){
			if (card[j]!='\0'){
		 		header[index++] = card[j];
			}
			else{
				header[index++] = ' ';
			}
		}
		header[index++] = '\n';
	}
	header[81*nkey-1] = '\0';
	
    if(fits_close_file(fptr, &status)){
		*error_status = status;
		return(header);
	}
	
	*error_status = status;
	
    return(header);
};

Appendix II: Useful FITS links

An excellent site for general FITS information is:  http://fits.cv.nrao.edu

A general introduction to FITS can be found at: http://fits.gsfc.nasa.gov/fits_home.html

The CFITSIO library exists for various flavours of Unix, VMS, Windows, and MacOS, and is freely available at http://legacy.gsfc.nasa.gov/docs/software/fitsio/fitsio.html

Appendix III: Tips for image processing in Mathematica

Users of other software packages (such as IDL or MATLAB) should remember that Mathematica keeps a precise record of the results from each step in your calculations in the Out[] symbols. This means that even if you remove a variable holding an image from memory, you're still storing the results you obtained earlier from your calculations with that image.  If, as is often the case, the intermediate steps in a calculation are also images, they could be taking up a collosal amount of storage.  For example, if the variable "a" holds an image and you evaluate it a few times by hitting "shift-return" in succession,  even if you delete the image from memory you can always get it back because it's stored in the Out[] variables corresponding to your previous steps. 

Having access to the results from intermediate steps in long calculations is often a very useful feature, but if you're running into memory limitations you may not appreciate it much.  If you're willing to forego this feature of Mathematica and would prefer to recover the memory being used to store the results from previous calculations, then you can use the CleanSlate package.  This is shipped with Mathematica version 4 (on my system it lives in "C:\\Program Files\\Wolfram Research\\Mathematica\\4.0\\AddOns\\ExtraPackages\\Utilities\\CleanSlate.m"), although it is not documented in the on-line help.  This package has a ClearInOut[] function that deletes the intermediate values from your calculations stored in the Out[] symbols,  while preserving the current values of your variables.