FITSREAD reads a FITS image into a Matlab variable.
data = fitsread(filename);
data = fitsread(filename,range);
data = fitsread(filename,range,scaleQ);
The first form indicates that the full image should be read in. The
data is scaled by the BZERO and BSCALE keywords, which are assumed to
exist in the file.
The second form indicates that a subimage should be read in. The
subimage is specified by a string of the form 'xlow:xhigh,ylow:yhigh'.
The string 'full' can also be used if the full image is to be read in.
The data will be scaled using the BZERO and BSCALE keywords, which
are assumed to exist in the file.
The third form is like the second form, except a 1 or 0 is used to
indicate whether the BZERO and BSCALE keywords should be searched
for and used to scale the data (1=yes, 0=no). This should be
unnecessary unless you're reading in non-standard FITS files where
the BZERO and BSCALE keywords do not exist.
Notes:
1. If a filename with no extension is supplied, this routine looks
for a file with an extension of '.fits'.
2. If a sub-image is read in only the memory required to hold the
sub-image is used by the function (ie. the whole image is not read
in, which would be wasteful and possibly very slow).
3. To compare the image with one viewed with a standard image viewer
such as xv or SAOimage be sure to display the image with the
following settings:
axis xy;
axis image;
so that the origin is at the lower left and the correct aspect ratio
is shown.
Examples:
im=fitsread('n2715_J'); % Assumes .fits extension
im=fitsread('n2715_lJ.fits','10:313,15:175'); % Reads a subimage
im=fitsread('n2715_lJ.fits','full',0); % Full image, no scaling
Known deficiences:
1. Does not read anything more complicated than a simple 2D image
ie. no binary tables or data cubes.
2. Only does a very cursory inspection of possible architectures
in order to figure out if the data needs to be byte swapped.
This is because I only test for architectures I am familiar with.
(I believe this routine will work with PC-compatibles, Macs, Sun
workstations, and possibly DEC Alphas).
Version 2.5
R. Abraham, Institute of Astronomy, Cambridge University0001 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 = 'uint16=>uint16'; 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)=swapbytes(fread(file,nc,type)); 0273 fseek(file,byte_column_skip,'cof'); 0274 end 0275 0276 %Scale the data 0277 X = bscale*X + bzero; 0278 0279 %Clean up and output data 0280 fclose(file); 0281 0282 %Reshape to a 2D image 0283 output_image=reshape(X,nc,nr)'; 0284 0285 0286 0287 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0288 0289 function [keyword,value,comment] = parse_card(s) 0290 %Parses a FITS header card. 0291 %Reference: 0292 % NASA/Science Office of Standards and Technology 0293 % "Definition of the Flexible Image Transport System (FITS)" 0294 % NOST 100-1.0 June 19, 1993 0295 0296 %Set defaults 0297 keyword=[];value=[];comment=[]; 0298 0299 %Get keyword in bytes 1 - 8 0300 keyword=s(1:8); 0301 if nargout==1 0302 return; 0303 end 0304 0305 %If keyword is blank then the line is a comment 0306 if strmatch(keyword,' ') 0307 keyword=[]; 0308 value=[]; 0309 comment=deblank(s(11:80)); 0310 return; 0311 end; 0312 0313 0314 %Keyword is non-blank. Check if there is a corresponding value/comment. 0315 %If not then the only possibilities are that bytes 11:80 are a comment 0316 %or that they are blank 0317 if ~strmatch(s(9:10),'= ') 0318 keyword=deblank(keyword); 0319 value=[]; 0320 comment=deblank(s(11:80)); 0321 return; 0322 end; 0323 0324 %Card is a standard keyword/value/comment structure. Break the value/comment 0325 %string (bytes 11 - 80) into separate strings by tokenizing on "/" character. 0326 %Remove the leading and trailing blanks on the value and the trailing blanks 0327 %on the comment. 0328 0329 keyword=deblank(keyword); 0330 [value,comment]=strtok(s(11:80),'/'); 0331 comment=deblank(comment); 0332 value=fliplr(deblank(fliplr(deblank(value)))); 0333 0334 %Now figure out whether to output the value as a string or as a number. 0335 %The FITS standard requires all values that are strings to be in single 0336 %quotes like this: 'foo bar', so I can simply look for occurences of a 0337 %single quote to flag a string. However, logical variables can take the 0338 %values T or F without having any single quotes, so I'll have to look 0339 %out for those also. 0340 0341 %Test for logical. Return logical as a string. 0342 if strmatch(upper(value),'T') | strmatch(upper(value),'F') 0343 return; 0344 end; 0345 0346 %Test for string. Return string unconverted. 0347 if length(findstr('''',value)) ~= 0 0348 return; 0349 end; 0350 0351 %Only thing left is a number. Convert string to number. 0352 value=str2num(value); 0353 0354 0355 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0356 0357 function [lx,hx,ly,hy]=lowhigh(str) 0358 % LOWHIGH extracts range elements from a range string 0359 % of the form 'XL:XH,YL:YH'. The elements are returned 0360 % as numbers not as strings. 0361 0362 [dx,dy]=strtok(str,','); 0363 [lx,hx]=strtok(dx,':'); 0364 [ly,hy]=strtok(dy,':'); 0365 0366 % remove bogus characters 0367 hx = strrep(hx,':',''); 0368 ly = strrep(ly,',',''); 0369 hy = strrep(hy,':',''); 0370 0371 % return as a number instead of as a string 0372 lx = str2double(lx); 0373 hx = str2double(hx); 0374 ly = str2double(ly); 0375 hy = str2double(hy);