MATLAB - Opening Astrometry Output (.FITS) Files

[Simon George]

Has anyone had any experience in opening Astrometry output (.FITS -new) files using MATLAB?

I have images collected in .fit format that open fine using MATLAB's default fitsread() function, but trying to feed them the same images after Astrometry has processed them and outut a "-new.fits" file causes an error in the CFITSIO library that MATLAB is using.

Before I'm forced to delve into coding something to read image data from these, has anyone got MATLAB script (/a solution) that they wouldn't mind sharing? It would be a great help!

(obviously, I'll post my results up here if there's not something existing that works) :)

Thanks!

Simon

[Dustin Lang]

Hi,

The -new.fits files is just the same image with a WCS header added... so what is it you're trying to do?

You could check with 'liststruc' or 'fitsverify' that the -new.fits files are indeed compliant with the FITS standard.

cheers,
--dustin

[Simon George]

I think the change in header confuses MATLAB. I'm sure it's nothing to do with the Astrometry files :)

I just read in the original files using MATLAB's built-in fitsread() function, which outputs a nice array of pixel values (16bit) for me that I can then manipulate and display as an image. Doing the same on the -new.fits files causes an error ["CFITSIO library error (108): error reading from FITS file"]

Error in matlab.io.fits.readImg (line 85)

imgdata = fitsiolib('read_subset',fptr,fpixel,lpixel,inc);

Error in fitsread>read_image_hdu (line 438)

    data = fits.readImg(fptr);

Error in fitsread (line 125)

        data = read_image_hdu(info,1,raw,pixelRegion);

Weird.

I've found an alternative fitsread() function [http://fits.gsfc.nasa.gov/software/fitsread.m] that, with some manipulation (there were bugs), reads the pixel data but for some reason doesn't give me the right numbers. In the original MATLAB fitsread() it gives me numbers between 0 and 65536 or thereabouts, but the alternative fitsread() gives values as -32000ish to +32000ish but not in the same scale, as if there's something strange going on with the BSCALE/BZERO stuff. I'll keep probing and trying to figure out what it's doing.

(FYI: the alternative fitsread() does the same value manipulation for my original .fit files as the Astrometry ones, so I know it's nothing wrong with Astrometry!)

I'm guessing that MATLAB gets confused by something in the WCS header, but I don't know for sure. If I get to the bottom of it, I'll report back.

Simon

[Simon George]

Okay, I've solved it.

I modified the code for fitsread() from http://uk.mathworks.com/matlabcentral/fileexchange/122-fits, and it needed some amendment to the way it was reading the binary data.

The lines:

% now start reading in the data. We read a column at a time and then

% skip to the next row.

X=zeros(nr*nc,1);

startIndex = 1;

finIndex = nc;

for i=1:nr  

   if i>1

      startIndex=finIndex+1;

      finIndex=startIndex+nc-1;

   end

  

   % X(startIndex:finIndex)=fread(file,nc,type);

   fseek(file,byte_column_skip,'cof');

end

can be replaced by:

X = fread(file,[nc nr], '*uint16', 'ieee-be');

to read in the image data in a way that matches MATLAB's own fitsread() function. The scaling and zeroing then needs to be applied, and the data converted to double format by

X = bscale*X - bzero;

X = im2double(X);

This should give you an output the same as the way MATLAB reads my original .fit images, and a workaround for where MATLAB doesn't like reading the Astrometry output files for whatever reason. 

I've attached the code here, in case anyone finds themselves with the same problem when reading in with MATLAB.

Simon