Description of fitsread3

0001 function output_image=fitsread(filename,subset,scaleQ)
0002 % FITSREAD reads a FITS image into a Matlab variable.
0003 %
0004 %  data = fitsread(filename);
0005 %  data = fitsread(filename,range);
0006 %  data = fitsread(filename,range,scaleQ);
0007 %
0008 %  The first form indicates that the full image should be read in. The
0009 %  data is scaled by the BZERO and BSCALE keywords, which are assumed to
0010 %  exist in the file.
0011 %
0012 %  The second form indicates that a subimage should be read in. The
0013 %  subimage is specified by a string of the form 'xlow:xhigh,ylow:yhigh'.
0014 %  The string 'full' can also be used if the full image is to be read in.
0015 %  The data will be scaled using the BZERO and BSCALE keywords, which
0016 %  are assumed to exist in the file.
0017 %
0018 %  The third form is like the second form, except a 1 or 0 is used to
0019 %  indicate whether the BZERO and BSCALE keywords should be searched
0020 %  for and used to scale the data (1=yes, 0=no). This should be
0021 %  unnecessary unless you're reading in non-standard FITS files where
0022 %  the BZERO and BSCALE keywords do not exist.
0023 %
0024 %  Notes:
0025 %    1. If a filename with no extension is supplied, this routine looks
0026 %    for a file with an extension of '.fits'.
0027 %
0028 %    2. If a sub-image is read in only the memory required to hold the
0029 %    sub-image is used by the function (ie. the whole image is not read
0030 %    in, which would be wasteful and possibly very slow).
0031 %
0032 %    3. To compare the image with one viewed with a standard image viewer
0033 %    such as xv or SAOimage be sure to display the image with the
0034 %    following settings:
0035 %
0036 %              axis xy;
0037 %              axis image;
0038 %
0039 %    so that the origin is at the lower left and the correct aspect ratio
0040 %    is shown.
0041 %
0042 %  Examples:
0043 %    im=fitsread('n2715_J');                        % Assumes .fits extension
0044 %    im=fitsread('n2715_lJ.fits','10:313,15:175');  % Reads a subimage
0045 %    im=fitsread('n2715_lJ.fits','full',0);         % Full image, no scaling
0046 %
0047 %  Known deficiences:
0048 %  1. Does not read anything more complicated than a simple 2D image
0049 %     ie. no binary tables or data cubes.
0050 %
0051 %  2. Only does a very cursory inspection of possible architectures
0052 %     in order to figure out if the data needs to be byte swapped.
0053 %     This is because I only test for architectures I am familiar with.
0054 %     (I believe this routine will work with PC-compatibles, Macs, Sun
0055 %     workstations, and possibly DEC Alphas).
0056 %
0057 % Version 2.5
0058 % R. Abraham, Institute of Astronomy, Cambridge University
0059 
0060 
0061 % Changes
0062 % Mar. 99 -- added the capability of extracting a subimage from the
0063 %            full image without having to read the whole image into
0064 %            memory.
0065 % Dec. 98 -- removed option to specify number of 36 card header units.
0066 %            An option to scale the data has been substituted instead.
0067 % Oct. 98 -- added code to detect 80X86 and Alpha machines and to read
0068 %            data into these with the big endian switch.
0069 
0070 %
0071 %Useful FITS Definitions:
0072 %
0073 %  "card" = 80 byte line of ASCII data in a file. This
0074 %           contains keyword/value pairs and/or comments.
0075 %
0076 %  "Header area" = a set of 36 "cards". Note that
0077 %          FITS files must have an integer number of header
0078 %          areas (typically between 1 and 6, depending on how much
0079 %          header information is stored in the file.)
0080 %
0081 %
0082 %Reading FITS files:
0083 %
0084 %The first few cards of the header area must give information
0085 %in a pre-defined order (eg. the data format, number of axes,
0086 %size etc) but after that the header keywords can come in any
0087 %order. The end of the last card giving information is flagged by
0088 %the keyword END, after which blank lines are added to pad the
0089 %last header area so it also contains 36 cards. After the last card
0090 %the data begins, in a format specified by the BITPIX keyword.
0091 %The dimensions of the data are specified by the NAXIS1 and
0092 %NAXIS2 keywords.
0093 %
0094 %Reference:   NASA/Science Office of Standards and Technology
0095 %           "Definition of the Flexible Image Transport System"
0096 %                            NOST 100-1.0
0097 %
0098 %           This and other FITS documents are available on-line at:
0099 %           http://www.gsfc.nasa.gov/astro/fits/basics_info.html
0100 
0101 
0102 
0103 % First try to figure out if we need to swap bytes. This is
0104 % imperfect as I don't know the endian-ness of each
0105 % architecture, so I'm only looking for ones I know for
0106 % sure.
0107 friend = computer;
0108 if strmatch(friend,'PCWIN')
0109    bswap = 'b';
0110 elseif strmatch(friend,'LNX86')
0111    bswap = 'b';
0112 elseif strmatch(friend,'ALPHA')
0113    bswap = 'b';
0114 else
0115    bswap = 'l';
0116 end
0117 
0118 
0119 %Figure out how to scale the data by looking for the BSCALE and
0120 %BZERO keywords. I think these can be anywhere in the header (ie.
0121 %they don't have mandatory locations). We look for these using the
0122 %"fitsheader" routine which scans through the whole header block
0123 %looking for specific keywords (a slow process). We want to
0124 %do this here before opening the file, because the fitsheader function
0125 %also opens the file internally and we want to avoid mixing calls to
0126 %"fitsheader" with explicit reads from the data file within this
0127 %M-file, to avoid screwing up the  internal file indexing.
0128 if nargin<2 
0129    subset = 'full';   
0130    scaleQ = 1;
0131 elseif nargin<3
0132    scaleQ = 1;
0133 end
0134 
0135 if scaleQ==1
0136    bscale=fitsheader(filename,'BSCALE');
0137    bzero=fitsheader(filename,'BZERO');
0138    if isempty(bscale) | isempty(bzero)
0139       warning('BSCALE or BZERO keyword missing from the FITS file.')
0140       disp('Assuming BSCALE=1, BZERO=0');
0141       bscale=1.;
0142       bzero=0.;
0143    end
0144 else
0145    bscale=1;
0146    bzero=0;
0147 end
0148 
0149 
0150 %Now we open the file and grab the required keywords from the first
0151 %five cards, which MUST contain specific information according to
0152 %the FITS standard. We take advantage of the mandatory locations
0153 %to try to speed things up by reading all the required keywords in
0154 %a single pass.
0155 fid=-1;
0156 if ~isstr(filename)
0157    filename=setstr(filename);
0158 end;
0159 if (isempty(findstr(filename,'.'))==1)
0160    filename=[filename,'.fits'];
0161 end
0162 [file,message] = fopen(filename,'r',bswap);
0163 if file == -1
0164    error(message);
0165 end
0166 [d,simple,d]=parse_card(setstr(fread(file,80,'uchar')'));
0167 [d,bitpix,d]=parse_card(setstr(fread(file,80,'uchar')'));
0168 [d,naxis,d]=parse_card(setstr(fread(file,80,'uchar')'));
0169 [d,naxis1,d]=parse_card(setstr(fread(file,80,'uchar')'));
0170 [d,naxis2,d]=parse_card(setstr(fread(file,80,'uchar')'));
0171 n_card=5;
0172 
0173 %Keep reading cards until one turns up with the keyword 'END'.
0174 keyword='NIL';
0175 while(~strcmp(deblank(upper(keyword)),'END'))
0176    n_card=n_card+1;
0177    card=setstr(fread(file,80,'uchar')');
0178    [keyword]=parse_card(card);
0179 end;
0180 %Go past the blank lines of padding before the start of the data
0181 if (rem(n_card,36) ~= 0)
0182    n_blanks = 36 - rem(n_card,36);
0183    dummy=fread(file,n_blanks*80,'uchar');
0184 end;
0185 
0186 %work out the range of rows and columns to read
0187 if strcmp('FULL',upper(deblank(subset)))
0188    lx=1;
0189    hx=naxis1;
0190    ly=1;
0191    hy=naxis2;
0192    nr=naxis2;
0193    nc=naxis1;
0194 else
0195    [lx,hx,ly,hy] = lowhigh(subset);   %note order here
0196    nc = (hx-lx)+1;
0197    nr = (hy-ly)+1;
0198 end
0199 
0200 
0201 %
0202 %   Geometry of image reading: Matlab's fread begins reading
0203 %   at the lower left of the image and moves to the right, then
0204 %   up. The x-axis is naxis1, and the y-axis is naxis2.
0205 %
0206 %                ----------------------------------------------
0207 %                |                                            |
0208 %                |                                            |
0209 %                |                                            |
0210 %                |                                            |
0211 %                |                                            |
0212 %                |                                            |
0213 %       naxis2   |                                            |
0214 %                |                                            |
0215 %                | NB. Matlab starts reading at X (lower left |
0216 %                |     of frame) and moves to the right...    |
0217 %                |                                            |
0218 %                | X------->                                  |
0219 %                ----------------------------------------------
0220 %                                  naxis1
0221 %
0222 
0223 
0224 %work out how many pixels to skip to get to first byte of data
0225 pixel_start_skip = (ly - 1)*naxis1 + (lx - 1);
0226 
0227 %work out how many pixels to skip from the start of a row to the beginning
0228 %of the next row
0229 pixel_row_skip = naxis1 - nc; 
0230 
0231 % move to first byte
0232 if bitpix==-64
0233    byte_start_skip = pixel_start_skip*8;
0234    byte_column_skip = pixel_row_skip*8;
0235    type = 'float32';   
0236    fseek(file,byte_start_skip,'cof');
0237 elseif bitpix==-32
0238    byte_start_skip = pixel_start_skip*4;
0239    byte_column_skip = pixel_row_skip*4;   
0240    type = 'float';
0241    fseek(file,byte_start_skip,'cof');   
0242 elseif bitpix==8
0243    byte_start_skip = pixel_start_skip*1;
0244    byte_column_skip = pixel_row_skip*1;  
0245    type = 'uint8';   
0246    fseek(file,byte_start_skip,'cof');
0247 elseif bitpix==16
0248    byte_start_skip = pixel_start_skip*2;
0249    byte_column_skip = pixel_row_skip*2;
0250    type = 'short';   
0251    fseek(file,byte_start_skip,'cof');
0252 elseif bitpix==32
0253    byte_start_skip = pixel_start_skip*4;
0254    byte_column_skip = pixel_row_skip*4;
0255    type = 'long';   
0256    fseek(file,byte_start_skip,'cof');
0257 else
0258    error('data type specified by BITPIX keyword is not -64, -32, 8, 16, or 32');
0259 end;
0260 
0261 
0262 % now start reading in the data. We read a column at a time and then
0263 % skip to the next row.
0264 X=zeros(nr*nc,1);
0265 startIndex = 1;
0266 finIndex = nc;
0267 for i=1:nr  
0268    if i>1
0269       startIndex=finIndex+1;
0270       finIndex=startIndex+nc-1;
0271    end
0272    X(startIndex:finIndex)=fread(file,nc,type);
0273    fseek(file,byte_column_skip,'cof');
0274 end
0275 
0276 %Scale the data
0277 keyboard;
0278 X = bscale*X + bzero;
0279 
0280 %Clean up and output data
0281 fclose(file);
0282 
0283 %Reshape to a 2D image
0284 output_image=reshape(X,nc,nr)';
0285 
0286 
0287 
0288 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0289 
0290 function [keyword,value,comment] = parse_card(s)
0291 %Parses a FITS header card.
0292 %Reference:
0293 %                NASA/Science Office of Standards and Technology
0294 %           "Definition of the Flexible Image Transport System (FITS)"
0295 %                        NOST 100-1.0    June 19, 1993
0296 
0297 %Set defaults
0298 keyword=[];value=[];comment=[];
0299 
0300 %Get keyword in bytes 1 - 8
0301 keyword=s(1:8);
0302 if nargout==1
0303    return;
0304 end
0305 
0306 %If keyword is blank then the line is a comment
0307 if strmatch(keyword,'       ')
0308     keyword=[];
0309     value=[];
0310     comment=deblank(s(11:80));
0311     return;
0312 end;
0313 
0314 
0315 %Keyword is non-blank. Check if there is a corresponding value/comment.
0316 %If not then the only possibilities are that bytes 11:80 are a comment
0317 %or that they are blank
0318 if ~strmatch(s(9:10),'= ')
0319     keyword=deblank(keyword);
0320     value=[];
0321     comment=deblank(s(11:80));
0322     return;
0323 end;
0324 
0325 %Card is a standard keyword/value/comment structure. Break the value/comment
0326 %string (bytes 11 - 80) into separate strings by tokenizing on "/" character.
0327 %Remove the leading and trailing blanks on the value and the trailing blanks
0328 %on the comment.
0329 
0330 keyword=deblank(keyword);
0331 [value,comment]=strtok(s(11:80),'/');
0332 comment=deblank(comment);
0333 value=fliplr(deblank(fliplr(deblank(value))));
0334 
0335 %Now figure out whether to output the value as a string or as a number.
0336 %The FITS standard requires all values that are strings to be in single
0337 %quotes like this: 'foo bar', so I can simply look for occurences of a
0338 %single quote to flag a string. However, logical variables can take the
0339 %values T or F without having any single quotes, so I'll have to look
0340 %out for those also.
0341 
0342 %Test for logical. Return logical as a string.
0343 if strmatch(upper(value),'T') | strmatch(upper(value),'F')
0344     return;
0345 end;
0346 
0347 %Test for string. Return string unconverted.
0348 if length(findstr('''',value)) ~= 0
0349     return;
0350 end;
0351 
0352 %Only thing left is a number. Convert string to number.
0353 value=str2num(value);
0354 
0355 
0356 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0357 
0358 function [lx,hx,ly,hy]=lowhigh(str)
0359 % LOWHIGH extracts range elements from a range string
0360 % of the form 'XL:XH,YL:YH'. The elements are returned
0361 % as numbers not as strings.
0362 
0363 [dx,dy]=strtok(str,',');
0364 [lx,hx]=strtok(dx,':');
0365 [ly,hy]=strtok(dy,':');
0366 
0367 % remove bogus characters
0368 hx = strrep(hx,':','');
0369 ly = strrep(ly,',','');
0370 hy = strrep(hy,':','');
0371 
0372 % return as a number instead of as a string
0373 lx = str2double(lx);
0374 hx = str2double(hx);
0375 ly = str2double(ly);
0376 hy = str2double(hy);
fitsread3

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
