function [hd, img, bytes_read] = geread(FILENAME, FIELDS)
% [hd, img, bytes_read] = geread(FILENAME, FIELDS)
%
% Read in a GE 5.x MRI image file.
%
% hd is a structure containing all of the header information.
%  hd.ih is the overall file header
%  
%  Other fields in the header are only read if the string fields contains
%  the appropriate one-letter flag.  The most useful fields are:
%  
%   'e' reads hd.exam (Exam number and patient information)
%   's' reads hd.series (Series number and MR spin echo information)
%   'm' reads hd.mri  (MR image number and image location information)
%
%  Other less useful fields include:
%
%   'i' reads hd.id (unique file id)
%   'p' reads hd.pack (image packing information, if defined)
%   'c' reads hd.compress (image decompression information, if defined)
%   'h' reads hd.histo (pixel histogram information)
%   't' reads hd.text (text annotation of the image)
%   'u' reads hd.user (user-defined data)
%   'S' reads hd.suite (ID of MRI machine that took the pictures)
%   'd' reads hd.debug, which contains info for debugging this function
%
%  EXAMPLE 
%  Read hd.ih, hd.exam, hd.series, hd.mri, and the image from file 'foo.mr':
%
%  [hd, img] = geread('foo.mr','esm'); 
%
%  HEADER VERSION INFORMATION
%
%  This script is based on the file headers.h obtained from an older GE 5.x 
%  machine with header version 3 (that is, hd.ih.header_version == 3),
%  and with Geneses version 7.0 (h.suite.verscre=='07', h.series.delta_cnt==0).
%
%  I have also loaded files from a Genesis 8.1, header version 2 machine.
%  Since I don't have access to a headers.h file, I don't know how accurate the
%  results are, but most of hd.suite, hd.exam, hd.series, and hd.mri look correct.
%  In particular, the all-important field hd.mri.corners loads correctly.
%  hd.histo looks wrong, and I think the pointers in hd.ih are wrong.
%

% Here is the list of possible fields, in the order they usually appear
POSSIBLE_FIELDS = { 'id', 'pack', 'compress', 'histo', 'text', 'graphics', ...
    'user', 'suite', 'exam', 'series', 'mri' };

% Default: don't read any fields other than hd.ih
if nargin<2, FIELDS=''; end;

% First, for safety, we find out if the file exists, and if so, what is its length
FILEINFO = dir(FILENAME);
if length(FILEINFO) < 1,
  disp(sprintf('Unable to read from file %s', FILENAME));
end

% Open the file, and read in the 'fixed-length' initial file header
[fid, ih, bytes_read] = read_file_header(FILENAME);
ih_length = bytes_read;
hd = struct('ih', ih);

% Add 'pack' to the FIELDS list if hd.ih.compress is 2 or 4
if (hd.ih.compress == 2 | hd.ih.compress == 4),
  FIELDS = strcat(FIELDS,'p');
end;

% Go through each of the possible fields, and read it if asked to do so
for fi=1:length(POSSIBLE_FIELDS),
  field = POSSIBLE_FIELDS{fi};
  
  % Check to see whether this field was requested by first letter
  tag = field(1);
  if strcmp(field,'suite'), tag='S'; end
  if findstr(FIELDS, tag),
  
    % Check to see whether the field pointer and length were given in IH
    p_string = sprintf('p_%s', field);
    l_string = sprintf('l_%s', field);
    if isfield(hd.ih, p_string) & isfield(hd.ih, l_string),
      
      % Check to make sure that we will not over-run the length of the file
      p_head = getfield(hd.ih, p_string);
      l_head = getfield(hd.ih, l_string);
      if (p_head >= 0) & (p_head+l_head <= FILEINFO.bytes),
	
	% Depending on the field, switch to appropriate header-read function
	switch field,
	  case 'id', [hd.id, bytes]=read_id(fid, p_head, l_head, hd.ih);
	  case 'pack', [hd.pack,bytes]=read_pack(fid, p_head, l_head, hd.ih);
	  case 'compress', [hd.compress,bytes]=read_compress(fid, p_head, l_head, hd.ih);
	  case 'histo', [hd.histo,bytes]=read_histo(fid, p_head, l_head, hd.ih);
	  case 'text', [hd.text,bytes]=read_text(fid,p_head, l_head, hd.ih); 
	  case 'user', [hd.user,bytes]=read_user(fid, p_head, l_head, hd.ih);
	  case 'suite', [hd.suite,bytes]=read_suite(fid, p_head, l_head, hd.ih); 
	  case 'exam', [hd.exam,bytes]=read_exam(fid, p_head, l_head, hd.ih); 
	  case 'series', [hd.series,bytes]=read_series(fid, p_head, l_head, hd.ih); 
	  case 'mri', [hd.mri,bytes]=read_mri(fid, p_head, l_head, hd.ih); 
	  otherwise, error(strcat('Unknown header field requested: ', field));
	end

	bytes_read = bytes_read + bytes;
	
      end % if p_head+l_head < FILEINFO.bytes
    end % if isfield(hd.ih, p_string) & isfield(hd.ih, l_string)
  end % if findstr(FIELDS, tag)
end; % for fi=1:length(POSSIBLE_FIELDS)

% Read in the image data, if an image is requested
%
if nargout>1,
  % First, move to the end of the header, assuming the image starts after header
  p_image = hd.ih.hdr_length;
  fseek(fid, p_image, 'bof');
    
  % Read off the rest of the file
  [tail, Ntail] = fread(fid, Inf, 'uint8');
  tail = uint8(tail);

  % If packed or compacked, call geread_i with left and wide, else call without
  if (hd.ih.compress == 2 | hd.ih.compress == 4),
    [ img, bytes ] = geread_i(tail,Ntail,hd.ih.width,hd.ih.height,hd.ih.depth,...
	hd.ih.compress,hd.pack.left,hd.pack.wide);
  else,
    [ img, bytes ] = geread_i(tail,Ntail,hd.ih.width,hd.ih.height,hd.ih.depth,hd.ih.compress);
  end;
  bytes_read = bytes_read + bytes;
end;

% Put tail into the header if requested
if findstr(FIELDS, 'd'),
	hd.debug.tail = tail;
end;

fclose(fid);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Private function file_header, visible only within this file
function [ fid, ih, bytes ] = read_file_header(FILENAME)
% [fid, ih, bytes] = read_file_header(FILENAME)
%
% Read in a 'fixed-length' file header from the file given by fid.
% The header must be at least 60 bytes in length, because bytes 57-60 give a pointer to
% the first variable header component, and this first pointer is assumed to also
% specify the length of the 'fixed-length' header.
%

% Open the specified file in read-only mode with big-endian data
fid=fopen(FILENAME,'r','b');
if fid<0,
  disp(sprintf('geread: Unable to open file %s for reading', FILENAME));
  return;
end

% Read in the 'magic number', and check to make sure it is 'IMGF'
[ magic, n ] = fread(fid, 4, 'char'); bytes = n;
magic = char(magic');
if ~strcmp(magic, 'IMGF'),
  disp(sprintf('geread: %s does not look like a GE file, but I will try to read it anyway', ...
      FILENAME));
end
ih = struct('magic', magic); bytes = 4;

% Read in basic file header information, starting from 5th byte 
fseek(fid, 4, 'bof');
[ ih.hdr_length, n ] = fread(fid, 1, 'int32'); bytes = bytes + 4*n;
[ ih.width, n ] = fread(fid, 1, 'int32'); bytes = bytes +4*n;
[ ih.height, n ] = fread(fid, 1, 'int32'); bytes = bytes + 4*n;
[ ih.depth, n ] = fread(fid, 1, 'int32'); bytes = bytes + 4*n;
[ ih.compress, n ] = fread(fid, 1, 'int32'); bytes = bytes + 4*n;
[ ih.default_window, n ] = fread(fid, 1, 'int32'); bytes = bytes + 4*n;
[ ih.default_level, n ] = fread(fid, 1, 'int32'); bytes = bytes + 4*n;
[ ih.background, n ] = fread(fid, 1, 'int32'); bytes = bytes + 4*n;
[ ih.ovrflow, n ] = fread(fid, 1, 'int32'); bytes = bytes + 4*n;
[ ih.undflow, n ] = fread(fid, 1, 'int32'); bytes = bytes + 4*n;
[ ih.top_offset, n ] = fread(fid, 1, 'int32'); bytes = bytes + 4*n;
[ ih.bot_offset, n ] = fread(fid, 1, 'int32'); bytes = bytes + 4*n;
[ ih.header_version, n ] = fread(fid, 1, 'short'); bytes = bytes + 2*n;
[ ih.checksum, n ] = fread(fid, 1, 'ushort'); bytes = bytes + 2*n;

% Set IH header length according to the header version
if ih.header_version==2,
  ih_length = 124;
else
  ih_length = 156;
end

% Read in the first pointer, and clip IH header length accordingly
[ ih.p_id, n ] = fread(fid, 1, 'int32'); bytes = bytes + 4*n;
[ ih.l_id, n ] = fread(fid, 1, 'int32'); bytes = bytes + 4*n;
ih_length = min(ih_length, ih.p_id);

% Keep track of the highest header byte
first_nonheader_byte=0;

% Read the first set of pointers 
fields = {'pack', 'compress', 'histo', 'text', 'graphics', 'dbHdr' };
for i=1:length(fields),
  if bytes < ih_length,
    [ dat, n ] = fread(fid, 2, 'int32'); bytes = bytes + 4*n;
    ih = setfield(ih, sprintf('p_%s', fields{i}), dat(1));
    ih = setfield(ih, sprintf('l_%s', fields{i}), dat(2));
    first_nonheader_byte = max(first_nonheader_byte, getfield(ih, sprintf('p_%s', fields{i})) + ...
	getfield(ih, sprintf('l_%s', fields{i})));
  end
end

% Read the levelOffset field
if bytes < ih_length,
  [ ih.levelOffset, n ] = fread(fid, 1, 'int32'); bytes = bytes + 4*n;
end

% In the Genesis 8.1 header ("version 2"), pointers to suite, exam, series, and mri seem to
% be undefined.  The correct pointers seem to be the following 'default values', but I don't
% know what is the relationship between these values and the pointers given in hd.ih.
default_pointers = [ 0, 124, 240, 1280, 2308 ];
default_lengths = [ 0, 114, 1024, 1020, 1022 ];

% Read the second set of pointers
fields = {'user', 'suite', 'exam', 'series', 'mri'};
for i=1:length(fields),
  if bytes < ih_length,
    % If there are still bytes in the specified header, read them
    [ dat, n ] = fread(fid, 2, 'int32'); bytes = bytes + 4*n;
    ih = setfield(ih, sprintf('p_%s', fields{i}), dat(1));
    ih = setfield(ih, sprintf('l_%s', fields{i}), dat(2));  
  else
    % If there are not still bytes in the specified header, we guess the next field's length
    ih = setfield(ih, sprintf('p_%s', fields{i}), default_pointers(i));
    ih = setfield(ih, sprintf('l_%s', fields{i}), default_lengths(i));
  end
  % Update the pointer to first nonheader byte
  first_nonheader_byte = max(first_nonheader_byte, getfield(ih, sprintf('p_%s', fields{i})) + ...
      getfield(ih, sprintf('l_%s', fields{i})));
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Private function read_id, visible only within this file
function [s,bytes]=read_id(fid, p_head, length, ih)
% [s,bytes]=read_id(fid, p_head, length, ih)
%
% Read the unique identifier from a GE 5.x file.
%
% Field names are invented, to fit the apparent use.
%

% Move to start of the ID header
fseek(fid, p_head, 'bof');

% Read in the apparent fields
[d,n]=fread(fid,[1,11],'uchar'); bytes=n;
s=struct('type',sscanf(char(d),'%s'));
[s.no,n]=fread(fid,[1,2],'uint8'); bytes=bytes+n;
[d,n]=fread(fid,[1,11],'uchar'); bytes=bytes+n;
s.image=printable_part(d);

% Move back to field beginning, and read in byte data
fseek(fid, p_head, 'bof');
[s.data, n]=fread(fid, [1, length], 'uint8');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Private function read_pack, visible only within this file
function [s,bytes]=read_pack(fid, p_head, length, ih)
% [s,bytes]=read_pack(fid, p_head, length, ih)
%
% Read the pack information from a GE 5.x file.
%
% According to the Medical Image FAQ, this is a table of
% shorts, two for every row in the image: pixels_to_left_of_line, 
% and pixels_stored_in_line.
%

% Move to start of the pack header
fseek(fid, p_head, 'bof');

% How many rows specified?
n_rows = floor(length/4);
if n_rows < ih.height,
  sprintf('pack header contains only %g rows instead of expected %d', n/2, ih.height),
end
if 4*n_rows < length,
  sprintf('pack header length is odd (%d bytes, %g rows); throwing away extra byte',length,n_rows),
end

% Read in the data
[d,n] = fread(fid,[2, n_rows],'short'); bytes = 2*n;

% Move to left, wide
s = struct('left', d(1,:)', 'wide', d(2,:)');

% Move back to field beginning, and read in byte data
fseek(fid, p_head, 'bof');
[s.data, n]=fread(fid, [1, length], 'uint8');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Private function read_compress, visible only within this file
function [s,bytes]=read_compress(fid, p_head, length, ih)
% [s,bytes]=read_compress(fid, p_head, length, ih)
%
% Read the compress information from a GE 5.x file.
%

% Move to start of the header section
fseek(fid, p_head, 'bof');

% According to the image format FAQ, this section gives start values 
% at different points in the compressed data.  I think these values should
% probably be short integers, but I'm not sure
%
[d,n] = fread(fid,[2,floor(length/4)],'short');
bytes = 2*n;
s = struct('index',d(1,:),'seed',d(2,:));

% Move back to field beginning, and read in byte data
fseek(fid, p_head, 'bof');
[s.data, n]=fread(fid, [1, length], 'uint8');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Private function read_histo, visible only within this file
function [s,bytes]=read_histo(fid, p_head, length, ih)
% [s,bytes]=read_histo(fid, p_head, length, ih)
%
% Read the histogram information from a GE 5.x file.
%
% s is a structure containing data using from the offsets 
% given in if_offsets.h.
%

% Move to start of the header section
fseek(fid, p_head, 'bof');

% Length of fixed part of the structure
s_length = 20;

% Read data in predefined fields
%
bytes = 0;
[d, n] = fread(fid,1,'int'); bytes = 4*n;
s = struct('version',d);
[s.sd, n] = fread(fid,1,'float'); bytes = bytes+4*n;
[s.mean, n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.min, n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.max, n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.first, n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.region, n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.length, n] = fread(fid,1,'short'); bytes = bytes+2*n;

% Expand s_length by the number of bins
s_length = s_length + 2*s.length + 4;

% Read the histogram data
%
[s.bins, n] = fread(fid, s.length, 'short'); bytes = bytes + 2*n;
[ s.hist_left, n ] = fread(fid, 1, 'short'); bytes = bytes + 2*n;
[ s.hist_right, n ] = fread(fid, 1, 'short'); bytes = bytes + 2*n;

% Move back to field beginning, and read in byte data
fseek(fid, p_head, 'bof');
[s.data, n]=fread(fid, [1, length], 'uint8');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Private function read_text, visible only within this file
function [s,bytes]=read_text(fid, p_head, length, ih)
% [s,bytes]=read_text(fid, p_head, length, ih)
%
% Read the text information from a GE 5.x file.
%
% Assume that this is just a string of the given length.
%

% Move to start of the header section
fseek(fid, p_head, 'bof');

% Read characters from the file
[dat, bytes] = fread(fid,[1,length],'uchar');

% Create a structure with both dat and text
s = struct('text', printable_part(dat), 'dat', dat);

% Move back to field beginning, and read in byte data
fseek(fid, p_head, 'bof');
[s.data, n]=fread(fid, [1, length], 'uint8');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Private function read_user, visible only within this file
function [s,bytes]=read_user(fid, p_head, length, ih)
% [s,bytes]=read_user(fid, p_head, length, ih)
%
% Read the user information from a GE 5.x file.
%
% I don't know how to read this information, but it seems
% to come in groups of four, so for now I assume that 
% it is an array of type int.
%

% Move to start of the header section
fseek(fid, p_head, 'bof');

[d,n] = fread(fid,[1,length],'uint8');
s = struct('text',printable_part(d),'data',d);
bytes = n;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Private function read_suite, visible only within this file
function [s,bytes]=read_suite(fid, p_head, length, ih)
% [s,bytes]=read_suite(fid, p_head, length, ih)
%
% Read the suite header from a GE 5.x file.
%
% s is a structure containing data using from the offsets 
% given in if_offsets.h.
%

% Move to start of the header section
fseek(fid, p_head, 'bof');

% Declared structure length
s_length = 114;

% Read data in predefined fields
%
[d,n]=fread(fid,[1,4],'char'); bytes=n;
s = struct('id',printable_part(d));
[s.uniq, n] = fread(fid,1,'short'); bytes=bytes+2*n;
[d,n]=fread(fid,[1,1],'char'); bytes=bytes+n;
s.diskid = printable_part(d);

[d,n]=fread(fid,[1,13],'char'); bytes=bytes+n;
s.prodid = printable_part(d);
[d,n]=fread(fid,[1,2],'char'); bytes=bytes+n;
s.verscre = printable_part(d);
[d,n]=fread(fid,[1,2],'char'); bytes=bytes+n;
s.verscur = printable_part(d);
[s.checksum, n] = fread(fid,1,'uint32'); bytes=bytes+4*n;

% Move back to field beginning, and read in byte data
fseek(fid, p_head, 'bof');
[s.data, n]=fread(fid, [1, length], 'uint8');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Private function read_exam, visible only within this file
function [s, bytes]=read_exam(fid, p_head, length, ih)
% [s, bytes]=read_exam(fid, p_head, length, ih)
%
% Read the exam header from a GE 5.x file.
%
% s is a structure containing data using from the offsets 
% given in if_offsets.h.
%

% Move to start of the header section
fseek(fid, p_head, 'bof');

% Declared structure length
s_length = 1024;

% Warn if length is too small
%
if length < s_length,
 sprintf('Warning: hd.exam length %d bytes, less than the expected %d; returning one-field struct.', length, s_length);
 s = struct('dat',fread(fid,[1,length],'uint8'));
 return;
end

% Read data in predefined fields
%
[d,n]=fread(fid,[1,4],'uchar'); bytes=n;
s = struct('suid',printable_part(d));
[s.uniq,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[d,n]=fread(fid,[1,2],'uchar'); bytes=bytes+n;
s.diskid = printable_part(d);
[s.no,n] = fread(fid,1,'ushort'); bytes=bytes+2*n;
[d,n]=fread(fid,[1,34],'uchar'); bytes=bytes+n;
s.hospital = printable_part(d);

[s.detect,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.numcells,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.zerocell,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.cellspace,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.srctodet,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.srctoiso,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.tubetyp,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.dastyp,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.num_dcnk,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.dcn_len,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.dcn_density,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.dcn_stepsize,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.dcn_shiftcnt,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.magstrength,n] = fread(fid,1,'int'); bytes=bytes+4*n;

% Check header version; if header is version #2, read off another 2 bytes here
if ih.header_version == 2,  trash = fread(fid, 1, 'short'); end

[d,n]=fread(fid,[1,13],'uchar'); bytes=bytes+n;
s.patid = printable_part(d);

% Check header version; if header is version #2, read off another 2 bytes here
if ih.header_version == 2,  trash = fread(fid, 1, 'short'); end

[d,n]=fread(fid,[1,25],'uchar'); bytes=bytes+n;
s.patname = printable_part(d);
[s.patage,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.patian,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.patsex,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.patweight,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.trauma,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[d,n]=fread(fid,[1,61],'uchar'); bytes=bytes+n;
s.hist = printable_part(d);
[d,n]=fread(fid,[1,13],'uchar'); bytes=bytes+n;
s.reqnum = printable_part(d);
[s.datetime,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[d,n]=fread(fid,[1,33],'uchar'); bytes=bytes+n;
s.refphy = printable_part(d);
[d,n]=fread(fid,[1,33],'uchar'); bytes=bytes+n;
s.diagrad = printable_part(d);
[d,n]=fread(fid,[1,4],'uchar'); bytes=bytes+n;
s.op = printable_part(d);
[d,n]=fread(fid,[1,23],'uchar'); bytes=bytes+n;
s.desc = printable_part(d);
[d,n]=fread(fid,[1,3],'uchar'); bytes=bytes+n;
s.typ = printable_part(d);
[s.format,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.firstaxtime,n] = fread(fid,1,'double'); bytes=bytes+8*n;
[d,n]=fread(fid,[1,10],'uchar'); bytes=bytes+n;
s.sysid = printable_part(d);
[s.lastmod,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.protocolflag,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[d,n]=fread(fid,[1,14],'uchar'); bytes=bytes+n;
s.alloc_key = printable_part(d);

% Check header version; if header is version #2, read off another 8 bytes here
if ih.header_version == 2,  trash = fread(fid, 8, 'uint8'); end

[s.delta_cnt,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[d,n]=fread(fid,[1,2],'uchar'); bytes=bytes+n;
s.verscre = printable_part(d);
[d,n]=fread(fid,[1,2],'uchar'); bytes=bytes+n;
s.verscur = printable_part(d);
[s.checksum,n] = fread(fid,1,'uint'); bytes=bytes+4*n;
[s.complete,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.seriesct,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.numarch,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.numseries,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.series,n] = fread(fid,[1,8],'uint8'); bytes=bytes+n;
[s.numunser,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.unseries,n] = fread(fid,[1,8],'uint8'); bytes=bytes+n;
[s.toarchcnt,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.toarchive,n] = fread(fid,[1,8],'uint8'); bytes=bytes+n;
[s.prospcnt,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.prosp,n] = fread(fid,[1,8],'uint8'); bytes=bytes+n;
[s.modelnum,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.modelcnt,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.models,n] = fread(fid,[1,8],'uint8'); bytes=bytes+n;
[s.stat,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[d,n]=fread(fid,[1,16],'uchar'); bytes=bytes+n;
s.uniq_sys_id = printable_part(d);
[d,n]=fread(fid,[1,16],'uchar'); bytes=bytes+n;
s.service_id = printable_part(d);
[d,n]=fread(fid,[1,4],'uchar'); bytes=bytes+n;
s.mobile_loc = printable_part(d);
[d,n]=fread(fid,[1,32],'uchar'); bytes=bytes+n;
s.study_uid = printable_part(d);
[s.study_status,n] = fread(fid,1,'short'); bytes=bytes+2*n;

% Move back to field beginning, and read in byte data
fseek(fid, p_head, 'bof');
[s.data, n]=fread(fid, [1, length], 'uint8');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Private function read_series, visible only within this file
function [s, bytes]=read_series(fid, p_head, length, ih)
% [s, bytes]=read_series(fid, p_head, length, ih)
%
% Read the series header from a GE 5.x file.
%
% s is a structure containing data using from the offsets 
% given in if_offsets.h.
%

% Move to start of the header section
fseek(fid, p_head, 'bof');

% Declared structure length
s_length = 1020;

% Warn if length is too small
%
if length < s_length,
 sprintf('Warning: hd.series length %d bytes, less than the expected %d; returning one-field struct.', length, s_length);
 s = struct('dat',fread(fid,[1,length],'uint8'));
 return;
end

% Read data in predefined fields
%
[d,n]=fread(fid,[1,4],'uchar'); bytes = n;
s = struct('suid',printable_part(d));
[s.uniq,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[d,n]=fread(fid,[1,2],'uchar'); bytes=bytes+n;
s.diskid = printable_part(d);
[s.exno,n] = fread(fid,1,'ushort'); bytes=bytes+2*n;
[s.no,n] = fread(fid,1,'short'); bytes=bytes+2*n;

[s.datetime,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.actual_dt,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[d,n]=fread(fid,[1,30],'uchar'); bytes=bytes+n;
s.desc = printable_part(d);
[d,n]=fread(fid,[1,9],'uchar'); bytes=bytes+n;
s.pr_sysid = printable_part(d);
[d,n]=fread(fid,[1,9],'uchar'); bytes=bytes+n;
s.pansysid = printable_part(d);
[s.typ,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.source,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.plane,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.scan_type,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.position,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.entry,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[d,n]=fread(fid,[1,4],'uchar'); bytes=bytes+n;
s.anref = printable_part(d);
[s.lmhor,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[d,n]=fread(fid,[1,26],'uchar'); bytes=bytes+n;
s.prtcl = printable_part(d);
[s.contrast,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[d,n]=fread(fid,[1,2],'uchar'); bytes=bytes+n;
s.start_ras = printable_part(d);

% Check header version; if header is version #2, read off another 2 bytes here
if ih.header_version == 2,  trash = fread(fid, 1, 'short'); end

[s.start_loc,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[d,n]=fread(fid,[1,2],'uchar'); bytes=bytes+n;
s.end_ras = printable_part(d);

% Check header version; if header is version #2, read off another 2 bytes here
if ih.header_version == 2,  trash = fread(fid, 1, 'short'); end

[s.end_loc,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.pseq,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.sortorder,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.lndmrkcnt,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.nacq,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.xbasest,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.xbaseend,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.xenhst,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.xenhend,n] = fread(fid,1,'short'); bytes=bytes+2*n;

% Check header version; if header is version #2, read off another 2 bytes here
if ih.header_version == 2,  trash = fread(fid, 1, 'short'); end

[s.lastmod,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[d,n]=fread(fid,[1,14],'uchar'); bytes=bytes+n;
s.alloc_key = printable_part(d);

% Check header version; if header is version #2, read off another 2 bytes here
if ih.header_version == 2,  trash = fread(fid, 1, 'short'); end

[s.delta_cnt,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[d,n]=fread(fid,[1,2],'uchar'); bytes=bytes+n;
s.verscre = printable_part(d);
[d,n]=fread(fid,[1,2],'uchar'); bytes=bytes+n;
s.verscur = printable_part(d);
[s.pds_a,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.pds_c,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.pds_u,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.checksum,n] = fread(fid,1,'uint'); bytes=bytes+4*n;
[s.complete,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.numarch,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.imagect,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.numimages,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.images,n] = fread(fid,[1,8],'uint8'); bytes=bytes+n;
[s.numunimg,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.unimages,n] = fread(fid,[1,8],'uint8'); bytes=bytes+n;
[s.toarchcnt,n] = fread(fid,1,'int'); bytes=bytes+4*n;
[s.toarchive,n] = fread(fid,[1,8],'uint8'); bytes=bytes+n;
[s.echo1_alpha,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo1_beta,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo1_window,n] = fread(fid,1,'ushort'); bytes=bytes+2*n;
[s.echo1_level,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.echo2_alpha,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo2_beta,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo2_window,n] = fread(fid,1,'ushort'); bytes=bytes+2*n;
[s.echo2_level,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.echo3_alpha,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo3_beta,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo3_window,n] = fread(fid,1,'ushort'); bytes=bytes+2*n;
[s.echo3_level,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.echo4_alpha,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo4_beta,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo4_window,n] = fread(fid,1,'ushort'); bytes=bytes+2*n;
[s.echo4_level,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.echo5_alpha,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo5_beta,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo5_window,n] = fread(fid,1,'ushort'); bytes=bytes+2*n;
[s.echo5_level,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.echo6_alpha,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo6_beta,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo6_window,n] = fread(fid,1,'ushort'); bytes=bytes+2*n;
[s.echo6_level,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.echo7_alpha,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo7_beta,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo7_window,n] = fread(fid,1,'ushort'); bytes=bytes+2*n;
[s.echo7_level,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[s.echo8_alpha,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo8_beta,n] = fread(fid,1,'float'); bytes=bytes+4*n;
[s.echo8_window,n] = fread(fid,1,'ushort'); bytes=bytes+2*n;
[s.echo8_level,n] = fread(fid,1,'short'); bytes=bytes+2*n;
[d,n]=fread(fid,[1,32],'uchar'); bytes=bytes+n;
s.series_uid = printable_part(d);
[d,n]=fread(fid,[1,32],'uchar'); bytes=bytes+n;
s.landmark_uid = printable_part(d);
[d,n]=fread(fid,[1,32],'uchar'); bytes=bytes+n;
s.equipmnt_uid = printable_part(d);

% Move back to field beginning, and read in byte data
fseek(fid, p_head, 'bof');
s.data = fread(fid, [1, length], 'uint8');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Private function read_mri, visible only within this file
function [s, bytes]=read_mri(fid, p_head, length, ih)
% [s, bytes]=read_mri(fid, p_head, length, ih)
%
% Read the mri header from a GE 5.x file.
%
% s is a structure containing data using from the offsets 
% given in if_offsets.h.
%
% bytes tells how many bytes were read from the file.
%

% Move to start of the header section
fseek(fid, p_head, 'bof');

% Declared structure length
s_length = 1022;

% Warn if length is too small
%
if length < s_length,
 sprintf('Warning: hd.mri length %d bytes, less than the expected %d; returning one-field struct.', length, s_length);
 s = struct('dat',fread(fid,[1,length],'uint8'));
 return;
end

% Read data in predefined fields
%
[d,n]=fread(fid,[1,4],'uchar'); bytes = n;
s = struct('suid',printable_part(d));
[s.uniq,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[d,n]=fread(fid,[1,2],'uchar'); bytes = bytes+n;
s.diskid = printable_part(d);
[s.exno,n] = fread(fid,1,'ushort'); bytes = bytes+2*n;
[s.seno,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.no,n] = fread(fid,1,'short'); bytes = bytes+2*n;

% Check header version; if header is version #2, read off another 2 bytes here
if ih.header_version == 2,  trash = fread(fid, 1, 'short'); end

[s.datetime,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[s.actual_dt,n] = fread(fid,1,'int'); bytes = bytes+4*n;

[s.scan_duration,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[s.slice_thickness,n] = fread(fid,1,'float'); bytes = bytes+4*n;

% s.imatrix -- image matrix size
[imatrix_X,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[imatrix_Y,n] = fread(fid,1,'short'); bytes = bytes+2*n;
s.imatrix = [imatrix_X imatrix_Y];

[s.dfov,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[s.dfov_rect,n] = fread(fid,1,'float'); bytes = bytes+4*n;

% s.dim -- dimension of the FFT in X and Y directions
[dim_X,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[dim_Y,n] = fread(fid,1,'float'); bytes = bytes+4*n;
s.dim = [dim_X dim_Y];

% s.pixsize -- size of each pixel in the X and Y directions
[pixsize_X,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[pixsize_Y,n] = fread(fid,1,'float'); bytes = bytes+4*n;
s.pixsize = [pixsize_X pixsize_Y];

[pdid,n]=fread(fid,[1,14],'uchar'); bytes = bytes+n;
[d,n]=fread(fid,[1,17],'uchar'); bytes = bytes+n;
s.contrastIV = printable_part(d);
[d,n]=fread(fid,[1,17],'uchar'); bytes = bytes+n;
s.contrastOral = printable_part(d);
[s.contmode,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.serrx,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.imgrx,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.screenformat,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.plane,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.scanspacing,n] = fread(fid,1,'float'); bytes = bytes+4*n;

% Check header version; if header is version #2, read off another 2 bytes here
if ih.header_version == 2,  trash = fread(fid, 1, 'short'); end

[s.compress,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.scouttype,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[d,n]=fread(fid,[1,2],'uchar'); bytes = bytes+n;
s.loc_ras = printable_part(d);

% Check header version; if header is version #2, read off another 2 bytes here
if ih.header_version == 2,  trash = fread(fid, 1, 'short'); end

[s.loc,n] = fread(fid,1,'float'); bytes = bytes+4*n;

% s.center = [center_R, center_A, center_S]
[ctr_R,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[ctr_A,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[ctr_S,n] = fread(fid,1,'float'); bytes = bytes+4*n;
s.center = [ctr_R, ctr_A, ctr_S];

% s.normalvector = [norm_R, norm_A, norm_S]
[norm_R,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[norm_A,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[norm_S,n] = fread(fid,1,'float'); bytes = bytes+4*n;
s.normalvector = [norm_R, norm_A, norm_S];

% s.corners is [top_left; top_right; bottom_left; bottom_right]
[tlhc_R,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[tlhc_A,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[tlhc_S,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[trhc_R,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[trhc_A,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[trhc_S,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[brhc_R,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[brhc_A,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[brhc_S,n] = fread(fid,1,'float'); bytes = bytes+4*n;
s.corners = [tlhc_R,tlhc_A,tlhc_S; trhc_R,trhc_A,trhc_S; ...
  0,0,0; brhc_R,brhc_A,brhc_S];
topdown = (s.corners(2,:) == s.corners(4,:));
s.corners(3,topdown) = s.corners(1,topdown);
s.corners(3,~topdown) = s.corners(4,~topdown);

% Continue with other fields
[d,n]=fread(fid,[1,4],'uchar'); bytes = bytes+n;
s.forimgrev = printable_part(d);

% Repetition, imaging, and echo times
[s.tr,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[s.ti,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[s.te,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[te2,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[s.numecho,n] = fread(fid,1,'short'); bytes = bytes+2*n;
if s.numecho==2, s.te = [s.te te2]; end;
[s.echonum,n] = fread(fid,1,'short'); bytes = bytes+2*n;

[s.tbldlta,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[s.nex,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[s.contig,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.hrtrate,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.tdel,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[s.saravg,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[s.sarpeak,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[s.monsar,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.trgwindow,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.reptime,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[s.imgpcyc,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.xmtgain,n] = fread(fid,1,'short'); bytes = bytes+2*n;

% Receive gains
[rcvgain1,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[rcvgain2,n] = fread(fid,1,'short'); bytes = bytes+2*n;
s.rcvgains = [rcvgain1 rcvgain2];

[s.mr_flip,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.mindat,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[s.cphase,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.swappf,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.pauseint,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.pausetime,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[s.obplane,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[s.slocfov,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[s.xmtfreq,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[s.autoxmtfreq,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[s.autoxmtgain,n] = fread(fid,1,'short'); bytes = bytes+2*n;

% Prescan repetitions?
[prescan_r1,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[prescan_r2,n] = fread(fid,1,'short'); bytes = bytes+2*n;
s.prescan = [prescan_r1 prescan_r2];

[s.user_bitmap,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[s.cenfreq,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.imode,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.iopt,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[s.pseq,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.pseqmode,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[d,n]=fread(fid,[1,34],'uchar'); bytes = bytes+n;
s.psdname = printable_part(d);
[s.psd_datetime,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[d,n]=fread(fid,[1,14],'uchar'); bytes = bytes+n;
s.psd_iname = printable_part(d);
[s.ctyp,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[d,n]=fread(fid,[1,18],'uchar'); bytes = bytes+n;
s.cname = printable_part(d);
[s.surfctyp,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.surfcext,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.rawrunnum,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[s.cal_fldstr,n] = fread(fid,1,'uint'); bytes = bytes+4*n;
[s.supp_tech,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.vbw,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[s.slquant,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.gpre,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.intr_del,n] = fread(fid,1,'int'); bytes = bytes+4*n;

% User data
[s.user,n] = fread(fid,[1 25],'float'); bytes = bytes+4*n;

% Check header version; if header is version #2, read off another 12 bytes here
if ih.header_version == 2,  trash = fread(fid, 12, 'uchar'); end

% Continue
[d,n]=fread(fid,[1,14],'uchar'); bytes = bytes+n;
s.alloc_key = printable_part(d);

% Check header version; if header is version #2, read off another 2 bytes here
if ih.header_version == 2,  trash = fread(fid, 2, 'uchar'); end

[s.lastmod,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[d,n]=fread(fid,[1,2],'uchar'); bytes = bytes+n;
s.verscre = printable_part(d);
[d,n]=fread(fid,[1,2],'uchar'); bytes = bytes+n;
s.verscur = printable_part(d);

% More possibly (?) useful fields
[pds_a,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[pds_c,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[pds_u,n] = fread(fid,1,'int'); bytes = bytes+4*n;
s.pds = [pds_a pds_c pds_u];

[s.checksum,n] = fread(fid,1,'uint'); bytes = bytes+4*n;
[s.archived,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[s.complete,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[s.satbits,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.scic,n] = fread(fid,1,'short'); bytes = bytes+2*n;

% Sagittal (?) x, y, and z locations 1 and 2
[satxloc1,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[satxloc2,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[satyloc1,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[satyloc2,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[satzloc1,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[satzloc2,n] = fread(fid,1,'short'); bytes = bytes+2*n;
s.satloc = [satxloc1, satyloc1, satzloc1; ...
  satxloc2, satyloc2, satzloc2];

% Image thickness (?) in x, y, and z directions
[satxthick,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[satythick,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[satzthick,n] = fread(fid,1,'short'); bytes = bytes+2*n;
s.satthick = [satxthick, satythick, satzthick];

[s.flax,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.venc,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.thk_disclmr,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.ps_flag,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.ps_status,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.image_type,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.vas_collapse,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.usern,n] = fread(fid,[1 2],'float'); bytes = bytes+4*n;
[s.proj_alg,n] = fread(fid,1,'short'); bytes = bytes+2*n;

% More data
[d,n]=fread(fid,[1,14],'uchar'); bytes = bytes+n;
s.proj_name = printable_part(d);

% s.rotation = [x_axis_rot, y_axis_rot, z_axis_rot]
[x_axis_rot,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[y_axis_rot,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[z_axis_rot,n] = fread(fid,1,'float'); bytes = bytes+4*n;
s.rotation = [x_axis_rot, y_axis_rot, z_axis_rot];

% s.thresh_min and s.thresh_max
[thresh_min1,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[thresh_max1,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[thresh_min2,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[thresh_max2,n] = fread(fid,1,'int'); bytes = bytes+4*n;
s.thresh_min_max = [thresh_min1, thresh_min2; thresh_max1, thresh_max2];

% More data
[s.echo_trn_len,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.frac_echo,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.prep_pulse,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.cphasenum,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.var_echo,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[d,n]=fread(fid,[1,1],'uchar'); bytes = bytes+n;
s.ref_img = printable_part(d);
[d,n]=fread(fid,[1,1],'uchar'); bytes = bytes+n;
s.sum_img = printable_part(d);
[s.img_window,n] = fread(fid,1,'ushort'); bytes = bytes+2*n;
[s.img_level,n] = fread(fid,1,'short'); bytes = bytes+2*n;

% Integer slopes
[slop_int_1,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[slop_int_2,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[slop_int_3,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[slop_int_4,n] = fread(fid,1,'int'); bytes = bytes+4*n;
[slop_int_5,n] = fread(fid,1,'int'); bytes = bytes+4*n;
s.slop_int = [slop_int_1, slop_int_2, slop_int_3, slop_int_4, slop_int_5];

% Floating point slopes
[slop_float_1,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[slop_float_2,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[slop_float_3,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[slop_float_4,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[slop_float_5,n] = fread(fid,1,'float'); bytes = bytes+4*n;
s.slop_float=[slop_float_1,slop_float_2,slop_float_3,slop_float_4,slop_float_5];

[d,n]=fread(fid,[1,16],'uchar'); bytes = bytes+n;
s.slop_str_1 = printable_part(d);
[d,n]=fread(fid,[1,16],'uchar'); bytes = bytes+n;
s.slop_str_2 = printable_part(d);
[s.scanactno,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.vasflags,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.vencscale,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[s.integrity,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.fphase,n] = fread(fid,1,'int'); bytes = bytes+4*n;

% Check header version; if header is version #2, read off another 2 bytes here
if ih.header_version == 2,  trash = fread(fid, 2, 'uchar'); end

[s.freq_dir,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.vas_mode,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.image_uid,n]=fread(fid,[1,32],'uchar'); bytes = bytes+n;
[s.sop_uid,n]=fread(fid,[1,32],'uchar'); bytes = bytes+n;

% "Don't use" fields
[dont_use_1,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[dont_use_2,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[dont_use_3,n] = fread(fid,1,'short'); bytes = bytes+2*n;

% More fields
[s.pscopts,n] = fread(fid,1,'short'); bytes = bytes+2*n;

% s.asoffset
[asoffsetx,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[asoffsety,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[asoffsetz,n] = fread(fid,1,'short'); bytes = bytes+2*n;
s.asoffset = [asoffsetx asoffsety asoffsetz];

[s.unoriginal,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.interleaves,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.effechospace,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.viewsperseg,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.rbpm,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.rtpoint,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.rcvrtype,n] = fread(fid,1,'short'); bytes = bytes+2*n;
[s.dbdt,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[s.dbdtper,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[s.estdbdtper,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[s.estdbdtts,n] = fread(fid,1,'float'); bytes = bytes+4*n;
[s.saravghead,n] = fread(fid,1,'float'); bytes = bytes+4*n;

% Move back to field beginning, and read in byte data
fseek(fid, p_head, 'bof');
s.data =fread(fid, [1, length], 'uint8');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Utility function printable_part -- remove non-printable characters
%
function text=printable_part(dat)

% Which characters do I replace?
replace = find(dat < 9 | (dat > 13 & dat < 32) | dat > 126);

% Build the text string
text = char(max(0,min(127,dat(:)')));
text(replace) = repmat(' ', size(replace));

