Code covered by the BSD License  

Highlights from
Mia 2.4

image thumbnail

Mia 2.4

by

 

02 Apr 2004 (Updated )

Medical Image Analysis GUI.

ecatfile( action, p1, p2, p3, p4, p5 )
function [ res1, res2, res3 ] = ecatfile( action, p1, p2, p3, p4, p5 )
% USAGE:
% [ fid, message ]       = ecatfile( 'open', filename, file_system );
% [ fid, message ]       = ecatfile( 'open', filename, 'ecat7' );
% [ fid, message ]       = ecatfile( 'open', filename, 'ecat6.4' );
% [ fid, message ]       = ecatfile( 'open', filename, '' );
% [ fid, message ]       = ecatfile( 'create', filename, mh );
% [ mh, message ]        = ecatfile( 'mh', fid );
% [ vol, hd, message ]   = ecatfile( 'read', fid, selmatrix );
% [ vol, hd, message ]   = ecatfile( 'read', fid, selmatrix, segment ); % 3D sinogram
% [ descr, hd, message ] = ecatfile( 'read_offset', fid, selmatrix );
% [ descr, hd, message ] = ecatfile( 'read_offset', fid, selmatrix, segment ); % 3D sinogram
% [ hd, message ]        = ecatfile( 'hd', fid, selmatrix );
% message                = ecatfile( 'writemh', fid, mh );
% message                = ecatfile( 'write', fid, vol, hd, selmatrix );
% message                = ecatfile( 'write', fid, vol, hd, selmatrix, segment ); % 3D sinogram
% message                = ecatfile( 'close', fid );
% message                = ecatfile( 'close', 'all' );
% [ matranges, message ] = ecatfile( 'matranges', fid );
% [ matlist, matstatus, message ] = ecatfile( 'matlist', fid );
% [ hd, hds, message ]   = ecatfile( 'blankhd', file_system );
% [ mh, hds, message ]   = ecatfile( 'blankhd', 'ecat7' );
% [ mh, hds, message ]   = ecatfile( 'blankhd', 'ecat6.4' );
% [ sh, hds, message ]   = ecatfile( 'blankhd', 'ecat7', file_type );
% [ sh, hds, message ]   = ecatfile( 'blankhd', 'ecat6.4', file_type );
%
% The user has the responsibility that the written headers are compatible
%   with the file and data structure
% Revision: April 24, 2001.
% ecatfile has been tested under MATLAB5.3. It has not been tested thoroughly under MATLAB6.
% 1999, 2000, 2001 by Flemming Hermansen.

% Known problem with opened files:
% The ecatfile.m reads all headers and all the index blocks when the ECAT file is opened. 
% Further calls to open will not reread the headers and index blocks if the file is 
% already open. Each call to open or create will generate different fid values for the 
% same file. The headers and the index blocks will only be discarded, when all fid's have 
% been closed (or when 'close all' is issued). This will give problems if the file is 
% modified outside the Matlab process, because the headers and index blocks in the Maltab 
% process then don't correspond to those in the physical file. This may also happen, if a 
% file is renamed or copied. In this situation the user may call 
%    message = ecatfile( 'close', 'all' )
% to make sure that ecatfile.m will reread the headers and index blocks. Usually, each 
% program should close a file after using it. Still, files may be left opened, if a 
% program crashes before closing the file. A subroutine that is used to read portions of a 
% file will thus open and close the file many times. This will impose a high overhead 
% because the headers and index blocks will be read repeatedly. This overhead may be 
% prevented if the calling program, itself, opens the file before calling the 
% subroutine, and closes the file when it has performed the last call to the subroutine.

% UNRESOLVED PROBLEM: I have not seen a definite report from CTI stating all possible
% formats. There may, therefore, be errors due to misunderstandings and ignorance of 
% certain formats. 
% Especially, I do not know the machine format of the file. I have assumed that 6.4
% files can be read by using machineformat 'vaxg' in the fopen routine. Likewise, I
% assume that ecat7 files can be opened using 'ieee-be'. I assume that the entire file
% be read by the same machineformat.

% The routines have been tested, and they can read and write the
% following ecat6.4 files types: 1:*.scn, 2:*.img, 3:*.atn, 4:*.nrm
% and the following ecat7 files types: 3:*.a, 7:*.v, 11:*.S, 14:*.S
% It can read the header from the following file type: 13:*.N
% It should also be able to read other headers (see the tables below), but they have not been tested.
%
% The written files have been tested binarily, and there are some differences:
% The CTI generated files do not contain a truly doubly linked list of index blocks, as index block 2
% points back to block 0.
% Floating point numbers are not reproduced binarily exact, but the read value is the same, 
% except for values very close to zero. (assumedly a problem with the representation of 
% underflow, or illegal values)
% The last block in CTI generated ecat7 3D sinograms are filled with non-zeros.
% The last block in CTI generated ecat7 3D attenuation is truncated.

%NOT IMPLEMENTED and known errors:
% Problem not yet solved (requires a call to the c-routine stat in stat.h):
% A file may be opened twice or more times at the same time. The filename must, 
% however, be spelled exactly the same way. Otherwise, the ecatfile-system 
% can not keep the headers and directoryblocks updated for all fids.
% 
% check that seek succeeds.
% seek during write should write zeros in front of the matrix, if the file is too short.
%
% The system updates the index blocks after the data have been written, so that
% the file always will contain valid data. The only exception is when 3D files are written.
% Then the index block is updated after the first 3D-segment has been written. The file
% will thus be too short, until the last 3D-segment has been written.

% ------------------------------------------------------

% pers.handlelist
% pers.filelist

% pers.filelist{}.filename
% pers.filelist{}.mainheader
% pers.filelist{}.dirblocks.data( 1:n )
% pers.filelist{}.dirblocks.filepos( 1:n )
% pers.filelist{}.pos
% pers.filelist{}.pos.dirblock()
% pers.filelist{}.pos.fpgdb()
% pers.filelist{}.pos.blockstart()
% pers.filelist{}.pos.blockend()
% pers.filelist{}.pos.status()
% pers.filelist{}.pos.header{} % may be empty

   global pers
%   persistent pers

if isempty( pers )
   pers = [];
   pers.handlelist = [];
   pers.filelist = [];
end

switch action
   
case 'blankhd'
% [ hd, hds, message ]   = ecatfile( 'blankhd', file_system );
% [ mh, hds, message ]   = ecatfile( 'blankhd', 'ecat7' );
% [ mh, hds, message ]   = ecatfile( 'blankhd', 'ecat6.4' );
   file_system = p1;
   if nargin == 2
      [ hds , message ] = getmhs( file_system );
      if isempty( message )
         [ hd , message ] = blankhd( hds );
      end
      res1 = hd;
      res2 = hds;
      res3 = message;
      return
   end
   file_type = p2;
   [ hds , message ] = getshs( file_system, file_type );
   if isempty( message )
      [ hd , message ] = blankhd( hds );
   end
   res1 = hd;
   res2 = hds;
   res3 = message;
   
case 'matlist'
   % [ matlist, matstatus, message ] = ecatfile( 'matlist', fid );
   res1 = [];
   res2 = [];
   fid = p1;
   
   res3 = 'Illegal handle';
   if fid > length( pers.handlelist ) | fid < 1
      res3 = 'Illegal handle';
      return;
   elseif pers.handlelist( fid ) == 0
      res3 = 'Illegal handle';
      return;
   end
   file = pers.filelist{ pers.handlelist( fid )  };
   
   res1 = file.pos.fpgdb;
   res2 = file.pos.status;
   res3 = '';
   return
   
case 'matranges'
   % [ matranges, message ] = ecatfile( 'matranges', fid );
   res1 = [];
   res2 = '';
   fid = p1;
   
   if fid > length( pers.handlelist ) | fid < 1
      res2 = 'Illegal handle';
      return;
   elseif pers.handlelist( fid ) == 0
      res2 = 'Illegal handle';
      return;
   end
   file = pers.filelist{ pers.handlelist( fid )  };
   
   [ ranges, message ] = getranges( file );
   res1 = ranges;
   res2 = message;
   return
   
case 'close'
   % message                 = ecatfile( 'close', fid );
   % message                 = ecatfile( 'close', 'all' );
   
   res1 = '';
   fid = p1;
   if strcmp( fid,'all' )
      pers = [];
      pers.handlelist = [];
      pers.filelist = [];
      res1 = '';
      return
   end
   
   if fid > length( pers.handlelist ) | fid < 1
      res1 = 'Illegal handle';
      return;
   elseif pers.handlelist( fid ) == 0
      res1 = 'Illegal handle';
      return;
   end
   
   fileno = pers.handlelist( fid );
   pers.handlelist( fid ) = 0;
   
   if any( pers.handlelist == fileno ) | fid < 1
      return
   end
   
   pers.filelist{ fileno } = [];
   return
   
case 'write'
   % message = ecatfile( 'write', fid, vol, hd, selmatrix );
   % message = ecatfile( 'write', fid, vol, hd, selmatrix, segment ); % 3D sinogram
   fid = p1;
   vol = p2;
   hd = p3;
   selmatrix = p4;
   % 3D sinogram segment:
   if nargin < 6
      selsegment = 1;
   else
      selsegment = p5;
   end

   res1 = '';
   
   postponeindexupdate = 0;
   
   xyz_dimension = size( vol );
   while length( xyz_dimension ) < 3
      xyz_dimension = [ xyz_dimension, 1 ];
   end
   if prod( xyz_dimension ) ~= prod( xyz_dimension(1:3 ) )
      error( 'dimension of vol is higher than 3' )
      return
   end
   
   if fid > length( pers.handlelist )
      res1 = 'Illegal handle';
      return;
   elseif pers.handlelist( fid ) == 0
      res1 = 'Illegal handle';
      return;
   end
   
   file = pers.filelist{ pers.handlelist( fid ) };
   
   [ datadescr, message ] = getdatadescr( file.mainheader.file_system, ...
      file.mainheader.file_type, hd.sh.data_type );
   
   if isempty( message )
      [ hd, offsetfromheader, writefill, numdatablocks, message ] = ...
         set_xyz_dimension( hd, datadescr, xyz_dimension, selsegment  );
   end
   
   if isempty( message )
      [ hds, message ] = getshs( file.mainheader.file_system, file.mainheader.file_type );
   end
   
   if ~isempty( message )
      return
   end
   
   % total number of blocks == numdatablocks + hds.header_size / 512
   
   fidphys = fopen( file.filename, 'r+', file.mhmachineformat );
   
   if isempty( file.pos.fpgdb )
      matpos = 0;
   else
      matpos = firstno( ...
         file.pos.fpgdb( 1, : ) == selmatrix( 1 ) & ...
         file.pos.fpgdb( 2, : ) == selmatrix( 2 ) & ...
         file.pos.fpgdb( 3, : ) == selmatrix( 3 ) & ...
         file.pos.fpgdb( 4, : ) == selmatrix( 4 ) & ...
         file.pos.fpgdb( 5, : ) == selmatrix( 5 ) );
   end
   
   % calculate number of blocks
   %   data_type.ftype = 'int16';
   %   data_type.size = 2;
   %   data_type.fopentype = 'ieee-be';
   
   if matpos == 0
      % entry does not exist, create entry
      % allocate space in the index block, 
      % do not write the updated index block yet
      % function matnum = fpgdb2matnum( fpgdb )
      
      % find an index block that is not full
      indexblockno = firstno( file.dirblocks.data( 1, : ) ~= 0 );
      
      if indexblockno == 0
         % all index blocks are full, create a new index block
         newindexblock = max( max( file.dirblocks.data( 1, : ) ), ...
            max( file.pos.blockend ) ) + 1;
         lastindexblock = file.dirblocks.filepos( end );
         file.dirblocks.data( 2, end ) = newindexblock;
         file.dirblocks.data( 3, 1 ) = newindexblock;
         file.dirblocks.filepos = [ file.dirblocks.filepos, newindexblock ];
         file.dirblocks.data( :, end + 1 ) = ...
            [ 31, 2, lastindexblock, 0, zeros( 1, 124 ) ];
         
         % the following writing sequence ensures, that next will be consistent in 
         % the links, even in case of interruption during writing
         % write new index block;
         status = fseek( fidphys, ( newindexblock - 1 ) * 512, -1 );
         if status ~= 0
            fclose( fidphys );
            res1 = 'Error writing the index block';
            return
         end
         count = fwrite( fidphys, file.dirblocks.data( :, end ), 'uint32' );
         if count ~= 128
            fclose( fidphys );
            res1 = 'Error writing the index block';
            return
         end
         % write first index block
         status = fseek( fidphys, ( 2 - 1 ) * 512, -1 );
         if status ~= 0
            fclose( fidphys );
            res1 = 'Error writing the index block';
            return
         end
         count = fwrite( fidphys, file.dirblocks.data( :, 1 ), 'uint32' );
         if count ~= 128
            res1 = 'Error writing the index block';
            return
         end
         % write last index block
         if length( file.dirblocks.filepos ) > 2
            status = fseek( fidphys, ( file.dirblocks.filepos( end - 1 ) - 1 ) * 512, -1 );
            if status ~= 0
               fclose( fidphys );
               res1 = 'Error writing the index block';
               return
            end
            count = fwrite( fidphys, file.dirblocks.data( :, end - 1 ), ...
               'uint32' );
            if count ~= 128
               fclose( fidphys );
               res1 = 'Error writing the index block';
               return
            end
         end
         pers.filelist{ pers.handlelist( fid )  } = file;
         
         indexblockno = length( file.dirblocks.filepos );
      end % create a new index block
      
      % add item into indexblock
      % 0: matnum
      % 1: first data block (subheader)
      % 2: last data block
      % 3: status: 1:exists; 0: allocated, no data yet written, -1: matrix deleted
      startblock = max( max( file.dirblocks.data( 2, : ) ), ...
         max( [ file.pos.blockend, 0 ] ) ) + 1;
      p = file.dirblocks.data( 4, indexblockno ) * 4 + 4;
      file.dirblocks.data( p+1:p+4, indexblockno ) = [ ...
            fpgdb2matnum( selmatrix ), ...
            startblock, ...
            startblock + numdatablocks + hds.header_size / 512 - 1, ...
            1 ];
      file.dirblocks.data( 1, indexblockno ) = file.dirblocks.data( 1, indexblockno ) - 1;
      file.dirblocks.data( 4, indexblockno ) = file.dirblocks.data( 4, indexblockno ) + 1;
      
      matpos = length( file.pos.dirblock ) + 1;
      file.pos.dirblock( matpos ) = indexblockno;
      file.pos.fpgdb( :, matpos ) = selmatrix( : );
      file.pos.blockstart( matpos ) = startblock;
      file.pos.blockend( matpos ) = ...
         startblock + numdatablocks + hds.header_size / 512 - 1;
      file.pos.status( matpos ) = 1;
      file.pos.header{ matpos } = hd.sh;
      file.pos.relpos( matpos ) = p / 4;
      
      postponeindexupdate = 1;
      % postpone the updating of file and writing the above changes to the hard disk
      % until the matrix has been succesfully written
   else
      % entry exists, check size of existing entry
      if file.pos.blockend( matpos ) - file.pos.blockstart( matpos ) + 1 ~= ...
            numdatablocks + hds.header_size / 512
         error( 'Trying to change the size of a matrix' )
         return
      end      
      file.pos.header{ matpos } = hd.sh;
      if file.pos.status ~= 1
         % change status to 1
         p = file.pos.relpos( matpos );
         file.pos.status( matpos ) = 1;
         file.dirblocks.data( p * 4 + 4, indexblockno ) = 1;
         postponeindexupdate = 1;
      end
   end
   
   % pers.filelist{}.pos.dirblock()
   % pers.filelist{}.pos.fpgdb()
   % pers.filelist{}.pos.blockstart()
   % pers.filelist{}.pos.blockend()
   % pers.filelist{}.pos.status()
   % pers.filelist{}.pos.header{} % may be empty
   
   status = fseek( fidphys, ( file.pos.blockstart( matpos ) - 1 ) * 512, -1 );
   if status ~= 0
      fclose( fidphys );
      res1 = 'Error writing matrix';
      return
   end
   
   message = writehd( fidphys, hd.sh, hds );
   if ~isempty( message )
      fclose( fidphys );
      res1 = message;
      return
   end
   
   % 'a', ftell( fidphys )
   
   % skip to the right 3D segment
   if offsetfromheader ~= 0
      fseek( fidphys, offsetfromheader, 0 );
   end
   
   for pl = 1:xyz_dimension( 3 )
      count = fwrite( fidphys, double( vol( :, :, pl ) ), datadescr.ftype );
      if count ~= prod( xyz_dimension( 1:2 ) )
         fclose( fidphys );
         res1 = 'Error writing matrix';
         return
      end
   end
   
   % fill last block with zeros
   count = fwrite( fidphys, zeros( 1, writefill / datadescr.size ), datadescr.ftype );
   if count ~= writefill / datadescr.size
      fclose( fidphys );
      res1 = 'Error writing matrix';
      return
   end
   
   if hds.header_size + ...
         ( file.pos.blockstart( matpos ) - 1 + numdatablocks ) * 512 ...
         ~= ftell( fidphys )
      % warning( 'The length of the written matrix is wrong.' );
   end
   
   % write the updated indexblock
   if postponeindexupdate
      status = fseek( fidphys, ( file.dirblocks.filepos( indexblockno ) - 1 ) * 512, -1 );
      if status ~= 0
         fclose( fidphys );
         res1 = 'Error writing the index block';
         return
      end
      count = fwrite( fidphys, file.dirblocks.data( :, indexblockno ), 'uint32' );
      if count ~= 128
         fclose( fidphys );
         res1 = 'Error writing the index block';
         return
      end
   end
   
   pers.filelist{ pers.handlelist( fid )  } = file;
   
   fclose( fidphys );
   return
case { 'read', 'hd', 'read_offset' }
% [ vol, hd, message ]   = ecatfile( 'read', fid, selmatrix );
% [ vol, hd, message ]   = ecatfile( 'read', fid, selmatrix, segment ); % 3D sinogram
% [ descr, hd, message ] = ecatfile( 'read_offset', fid, selmatrix );
% [ descr, hd, message ] = ecatfile( 'read_offset', fid, selmatrix, segment ); % 3D sinogram
% [ hd, message ]        = ecatfile( 'hd', fid, selmatrix );
   fid = p1;
   selmatrix = p2;
   % 3D sinogram segment:
   if nargin < 4
      selsegment = 1;
   else
      selsegment = p3;
   end
   vol = [ ];
   hd = [ ];
   descr = [];
   message = '';
   if fid > length( pers.handlelist ) | fid <= 0
      message = 'Illegal handle';
   elseif pers.handlelist( fid ) == 0
      message = 'Illegal handle';
   else
      file = pers.filelist{ pers.handlelist( fid )  };
      if isempty( file.pos.fpgdb )
         message = 'Matrix does not exist';
      else
         matpos = firstno( ...
            file.pos.fpgdb( 1, : ) == selmatrix( 1 ) & ...
            file.pos.fpgdb( 2, : ) == selmatrix( 2 ) & ...
            file.pos.fpgdb( 3, : ) == selmatrix( 3 ) & ...
            file.pos.fpgdb( 4, : ) == selmatrix( 4 ) & ...
            file.pos.fpgdb( 5, : ) == selmatrix( 5 ) );
         if matpos == 0
            message = 'Matrix does not exist';
         end
      end
   end
   
   % pers.filelist{}.pos.dirblock()
   % pers.filelist{}.pos.fpgdb()
   % pers.filelist{}.pos.blockstart()
   % pers.filelist{}.pos.blockend()
   % pers.filelist{}.pos.status()
   % pers.filelist{}.pos.header{} % may be empty
   
   if isempty( message )
      fidphys = fopen( file.filename, 'r', file.mhmachineformat );
      fseek( fidphys, ( file.pos.blockstart( matpos ) - 1 ) * 512, -1 );
      [ hds, message ] = getshs( file.mainheader.file_system, file.mainheader.file_type );
   else
      fidphys = -1;
   end
   
   if isempty( message )
      [ sh, message ] = readhd( fidphys, hds );
      hd.sh =sh;
      hd.mh = file.mainheader;
   end
   
   if strcmp( action, 'hd' )
      res1 = hd;
      res2 = message;
      if fidphys > 0
         fclose( fidphys );
      end
      return
   end
   
   % 'a', ftell( fidphys )
   
   if isempty( message )
      [ descr, message ] = getdatadescr( file.mainheader.file_system, file.mainheader.file_type, sh.data_type );
   end
   
   if isempty( message )
      [ xyz_dimension, offsetfromheader ] = get_xyz_dimension( hd, descr, selsegment );
      descr.xyz_dimension = xyz_dimension;
      descr.headeroffset = ( file.pos.blockstart( matpos ) - 1 ) * 512;
      descr.header_size = hds.header_size;
      descr.dataoffset = descr.headeroffset + hds.header_size + offsetfromheader;
   end
   
   if strcmp( action, 'read' ) & isempty( message )
      
      % skip to the right 3D segment
      if offsetfromheader ~= 0
         fseek( fidphys, offsetfromheader, 0 );
      end
      
      % preallocate vol:
      switch descr.ftype
      case 'float32'
         vol1 = zeros( xyz_dimension( 1 ), xyz_dimension( 2 ) );
      case { 'int8', 'int16' }
         vol1 = int16( zeros( xyz_dimension( 1 ), xyz_dimension( 2 ) ) );
      end
      vol = vol1( :, :, ones( xyz_dimension( 3 ), 1 ) );
      for pl = 1:xyz_dimension( 3 )
         [ vol1, count ] = fread( fidphys, prod( xyz_dimension( 1:2 ) ), ...
            descr.ftype );
         if count ~= prod( xyz_dimension( 1:2 ) ) 
            message = 'file too short';
            break
         end
         switch descr.ftype
         case 'float32'
            vol( :, :, pl ) = ...
               reshape( vol1, xyz_dimension( 1:2 ) );
         case { 'int8', 'int16' }
            vol( :, :, pl ) = ...
               reshape( int16( vol1 ), xyz_dimension( 1:2 ) );
         end
      end
   end 
   
   if strcmp( action, 'read' )
      if ~isempty( message )
         vol = [];
      end
      res1 = vol;
   else
      res1 = descr;
   end
   
   res2 = hd;
   res3 = message;
   if fidphys > 0
      fclose( fidphys );
   end
   
   return
   
case 'writemh'
   % message = ecatfile( 'writemh', fidphys, mh );
   fid = p1;
   mh = p2;
   res1 = '';
   if fid > length( pers.handlelist ) | fid <= 0
      res1 = 'Illegal handle';
   elseif pers.handlelist( fid ) == 0
      res1 = 'Illegal handle';
   end
   file = pers.filelist{ pers.handlelist( fid ) };
   fidphys = fopen( file.filename, ...
      'r+', file.mhmachineformat );
   if fidphys <= 0
      res1 = [ 'Can''t open file: ', file.filename ];
      return
   end
   
   [ hds, message ] = getmhs( mh.file_system );
   if isempty( message )
      message = writehd( fidphys, mh, hds )
   end
   fclose( fidphys );
   if ~isempty( message )
      res1 = message;
      return
   end
   pers.filelist{ pers.handlelist( fid )  }.mainheader = mh;
   return
   
case 'mh'
   % [ mh, message ] = ecatfile( 'mh', fid );
   res1 = [];
   res2 = '';
   fid = p1;
   if fid <= 0 | fid > length( pers.handlelist ) | fid <= 0
      res2 = 'Illegal handle';
      return;
   elseif pers.handlelist( fid ) == 0
      res2 = 'Illegal handle';
      return;
   end
   res1 = pers.filelist{ pers.handlelist( fid ) }.mainheader;
   return
   
case { 'open', 'create' }
% [ fid, message ]       = ecatfile( 'open', filename, file_system );
% [ fid, message ]       = ecatfile( 'open', filename, 'ecat7' );
% [ fid, message ]       = ecatfile( 'open', filename, 'ecat6.4' );
% [ fid, message ]       = ecatfile( 'open', filename, '' );
% [ fid, message ]       = ecatfile( 'create', filename, mh );
   filename = p1;
   if strcmp( filename, '' )
      fid = -1;
      message = 'Illegal filename';
      res1 = fid;
      res2 = message;
      return
   end
   
   if nargin < 3
      p2 = '';
   end
   
   for i = 1:length( pers.filelist )
      if ~isempty( pers.filelist{ i } )
         if strcmp( pers.filelist{ i }.filename, filename )
            if strcmp( action, 'create' )
               % file already open. Truncate the file.
               [ file, message ] = openecat( filename, 'create', p2 );
               if ~isempty( message )
                  fid = -1;
                  res1 = fid;
                  res2 = message;
                  return
               end
               pers.filelist{ i } = file;
            end
            % file already open, return a duplicate fid
            fid = firstno( [ pers.handlelist, 0 ] == 0 );
            pers.handlelist( fid ) = i;
            message = '';
            res1 = fid;
            res2 = message;
            return
         end
      end
   end
   
   [ file, message ] = openecat( filename, action, p2 );
   if ~isempty( message )
      fid = -1;
      res1 = fid;
      res2 = message;
      return
   end
   
   pos = length( pers.filelist ) + 1;
   for i = 1:length( pers.filelist )
      if isempty( pers.filelist{ i } ) 
         pos = i;
         break
      end
   end
   
   pers.filelist{ pos } = file;
   fid = firstno( [ pers.handlelist, 0 ] == 0 );
   pers.handlelist( fid ) = pos;
   message = '';
   res1 = fid;
   res2 = message;
   return
   
otherwise
   error( 'Illegal action' )
   
end % switch action

% ------------------------------------------------------------
function matnum = fpgdb2matnum( fpgdb )
frame = fpgdb( 1 );
plane = fpgdb( 2 );
gate = fpgdb( 3 );
data = fpgdb( 4 );
bed = fpgdb( 5 );
matnum = bitand( frame, 511 ) + ...
   bitshift( bitand( plane, 3 * 256 ), 9 - 8 ) + ...
   bitshift( bitand( data, 4 ), 11 - 2 ) + ...
   bitshift( bitand( bed, 15 ), 12 ) + ...
   bitshift( bitand( plane, 255 ), 16 ) + ...
   bitshift( bitand( gate, 63 ), 24 ) + ...
   bitshift( bitand( data, 3 ), 30 );

% ------------------------------------------------------------

function [ file, message ] = openecat( filename, create, file_system );
% function [ file, message ] = openecat( filename, 'open', file_system );
% function [ file, message ] = openecat( filename, 'create', mh );
if nargin < 2
   create = 'open';
end

message = '';
file = [];

% read mh, matdir, fpgdb

if strcmp( create, 'create' )
   % truncate or create
   mh = file_system;
   file_system = mh.file_system;
   switch mh.file_system
   case 'ecat7'
      fidphys = fopen( filename, 'w', 'ieee-be'  );
   case 'ecat6.4'
      fidphys = fopen( filename, 'w', 'vaxg'  );
   otherwise
      error( 'Wrong header' )
   end
   if fidphys == -1
      message = [ 'Can''t open file: ', filename ];
      return
   end
   [ hds, message ] = getmhs( mh.file_system );
   if isempty( message )
      message = writehd( fidphys, mh, hds );
   end
   if ~isempty( message )
      fclose( fidphys );
      return
   end
   
   fwrite( fidphys, 31, 'uint32' ); % 31 unused entries
   fwrite( fidphys,  2, 'uint32' ); % 1: next index block
   fwrite( fidphys,  2, 'uint32' ); % 2: previous index block ? 0 
   fwrite( fidphys,  0, 'uint32' ); % 3: number of used entries in the block?
   fwrite( fidphys, zeros( 4 * 31, 1 ), 'uint32' );
   if ftell( fidphys ) ~= 1024
      message = 'File write error';
   end
   fclose( fidphys );
end % create or truncate

[ file, message ] = readmh( filename, file_system );

if ~isempty( message )
   file = [];
   return
end

if strcmp( create, 'create' )
   fidphys = fopen( filename, 'r+', file.mhmachineformat  ); % 'w+' ?
else
   fidphys = fopen( filename, 'r', file.mhmachineformat  ); % 'w+' ?
end

if fidphys == -1
   file = [];
   message = [ 'Can''t open file: ', filename ];
   return
end

% read indexblocks

% first entry of an index block:
% 0: number of unused entries?
% 1: next index block
% 2: previous index block
% 3: number of used entries in the block?

% succeeding entries in an index block:
% 0: matnum
% 1: first data block (subheader)
% 2: last data block
% 3: status: 1:exists; 0: allocated, no data yet written, -1: matrix deleted

fseek( fidphys, 512, -1 );

blockno = 2;
matdir = [ ];
ok = 1;

file.dirblocks.data = [ ];
file.dirblocks.filepos = [ ];

while ok
   fseek( fidphys, 512 * ( blockno - 1), -1 );
   [ a1, count]= fread( fidphys, 128, 'uint32');
   if count < 128
      file = [];
      message = [ 'Can''t read index block: ', filename ];
      return
   end
   
   file.dirblocks.data = [ file.dirblocks.data, a1 ];
   file.dirblocks.filepos = [ file.dirblocks.filepos, blockno ];
   
   if a1(2) == 2
      ok = 0;
   elseif any( a1(2) == file.dirblocks.filepos( 2:end ) ) | a1( 2 ) <= 0
      % error if the directory block numbers are circular or not positive
      message = 'Illegal index block numbers';
      file = [];
      return;
   end
   
   blockno = a1(2);
   
end

% matdir(1:8, :)

%  fclose(fidphys);

sz = size( file.dirblocks.data );

direntries = [ ];
file.pos.dirblock = [ ];
file.pos.relpos = [ ];
for i = 1:sz( 2 )
   jj = file.dirblocks.data( 4, i );
   if jj < 0 | jj > 31
      message = 'Wrong index block';
      file = [];
      return
   end
   direntries = [ direntries, reshape( file.dirblocks.data( 5:4 + jj * 4, i ), 4, jj ) ];
   file.pos.dirblock = [ file.pos.dirblock, i + zeros( 1, jj ) ];
   file.pos.relpos = [ file.pos.relpos, 1:jj ];
end

fpgdb = file.dirblocks.data( 5:128, : );

fpgdb = [ zeros(7, sz(2)*31); reshape( fpgdb, 4, sz(2)*31 ) ];

m = direntries( 1, : );
frame   = rem( m, 512 ); m = floor( m / 512 );
hiplane = rem( m, 4 ); m = floor( m / 4 );
hidata  = rem( m, 2 ); m = floor( m / 2 );
bed     = rem( m, 16 ); m = floor( m / 16 );
plane   = rem( m, 256 ) + hiplane * 256; m = floor( m / 256 );
gate    = rem( m, 64 ); m = floor( m / 64 );
data    = rem( m, 4 ) + hidata * 4; m = floor( m / 4 );

file.pos.fpgdb = [ frame; plane; gate; data; bed ];
file.pos.blockstart = direntries( 2, : );
file.pos.blockend = direntries( 3, : );
file.pos.status = direntries( 4, : );

fclose( fidphys );

% ------------------------------------------------------
function [ ranges, message ] = getranges( file )

ranges = [ ];
message = '';

% check homogeneity
if length( file.pos.fpgdb( 1, : ) ) == 1
   m1 = file.pos.fpgdb;
   m2 = m1;
else
   m1 = min( file.pos.fpgdb' )';
   m2 = max( file.pos.fpgdb' )';
end
tot = m2 - m1 + 1;
f = 1;
for i = 1:5
   f( i, 1 ) = f( end ) * tot( i );
end

if f(5) ~= length( file.pos.status )
   message = 'Inhomogeneous fpgdb\n';
   return
end
testexists = zeros( f(5), 1 );
p = 1 + [ 1, f(1:4 )'] * ...
   ( file.pos.fpgdb - m1( :, ones( 1, size( file.pos.fpgdb, 2 ) ) ) );

testexists( p ) = 1;

if sum(testexists) ~= length( file.pos.status )
   message = 'Inhomogeneous fpgdb\n';
   return
end

ranges = [ m1, m2 ]';
if any( file.pos.status ~= 1 )
   message = 'Matrices with wrong status\n';
end

return

% ------------------------------------------------------

function [ hd, message ] = blankhd( hds )
message = '';
hd = [];
n = 0;
pos = [];

for i =1:length( hds.s )
   pos( end + 1 ) = n;
   switch hds.s(i).type
   case 'int16'
      n = n + 2 * hds.s(i).noelements;
      hd = setfield( hd, hds.s(i).fieldname, zeros( 1, hds.s(i).noelements ) );
   case { 'float32', 'mat1st', 'mat2nd', 'int32' }
      n = n + 4 * hds.s(i).noelements;
      hd = setfield( hd, hds.s(i).fieldname, zeros( 1, hds.s(i).noelements ) );
   case 'char'
      n = n + 1 * hds.s(i).noelements;
      hd = setfield( hd, hds.s(i).fieldname, '' );
   end % switch hds.s(i).type
end
if n ~= hds.header_size
   error( 'Blank header, wrong internal description' )
end

% make sure that ecat 7 mh contains a non-zero magic number
if strcmp( hds.file_system, 'ecat7' )
   hd.magic_number = 'MATRIX7';   
end

% ---------------------------------------------

function message = writehd( fidphys, hd, hds )
message = '';
initialpos = ftell( fidphys );

% make sure that ecat 7 mh contains a non-zero magic number
if strcmp( hd.file_system, 'ecat7' )
   if isempty( hd.magic_number )
      hd.magic_number = 'MATRIX7'; 
   elseif all( hd.magic_number == 0 )
      hd.magic_number = 'MATRIX7'; 
   end
end

ok = 1;
for i =1:length( hds.s )
   val = getfield( hd, hds.s(i).fieldname );
   if strcmp( hds.s(i).type, 'char' )
      if length( val(:) ) < hds.s(i).noelements
         val = [ val(:); char( ...
               zeros( hds.s(i).noelements - length( val(:) ), 1 ) ) ];
      end
   end
   if strcmp( hds.s(i).type, 'mat1st' )
      if length( val(:) ) ~= 12
         error( 'illegal header length\n' )
      end
   elseif strcmp( hds.s(i).type, 'mat2nd' )
      % the size has already been tested
   elseif length( val(:) ) ~= hds.s(i).noelements
      error( 'illegal header length\n' )
   end
   
   switch hds.s(i).type 
   case 'int16'
      ok = ok & length( val(:) ) == fwrite( fidphys, val, 'int16' );
   case 'int32'
      ok = ok & length( val(:) ) == fwrite( fidphys, val, 'int32' );
   case 'float32'
      ok = ok & length( val(:) ) == fwrite( fidphys, val, 'float32' );
   case 'mat1st'
      val = val( 1:3, 1:3 )';
      ok = ok & length( val(:) ) == fwrite( fidphys, val(:), 'float32' );
   case 'mat2nd'
      val = val(:,4);
      ok = ok & length( val(:) ) == fwrite( fidphys, val, 'float32' );
   case 'char'
      ok = ok & length( val(:) ) == fwrite( fidphys, val, 'uchar' );
   end % switch hds.s(i).type
end
if ~ok
   message = 'Can''t write header';
   return
end

if ftell( fidphys ) - initialpos ~= hds.header_size
   error( 'Writing header, wrong internal description' )
end

% ---------------------------------------------------
function [ file, message ] = readmh( filename, file_system )
% this function determines the header type and the machineformat 
% that is used for reading the main header and the index blocks.
% if the file is too short, it returns file = [], otherwise, it returns
%   file.mainheader = mh;
%   file.mhmachineformat = mhmachineformat;
%   file.filename = filename;
% the routine uses ad hoc testing to determine the format. This testing may
% need to be changed if the routine can not determine the format correctly   

message = '';
file = [];
mhmachineformat = 'ieee-be';
mh = [];

fidphys = fopen( filename, 'r', mhmachineformat );
if fidphys == -1
   message = [ 'Can''t open file: ', filename ];
   return
end

magic_number = fread( fidphys, 14, 'char' );
% check length of main header
[ mh, count ]= fread( fidphys, 256-7, 'uint16');
if count < 256-7
   message = 'Header read error';
   return
end

if ( all( magic_number == 0 ) & strcmp( file_system, '' ) ) | ...
      strcmp( file_system, 'ecat6.4' )
   fclose( fidphys );
   % ecat 6.4
   mhmachineformat = 'vaxg';
   fidphys = fopen( filename, 'r', mhmachineformat );
   if fidphys == -1
      return
   end
   hds = getmhs( 'ecat6.4' );
elseif strcmp( file_system, '' ) | strcmp( file_system, 'ecat7' ) 
   % ecat7
   fseek( fidphys, 0, -1 );
   hds = getmhs( 'ecat7' );
else
   message = 'Illegal file system';
   fclose( fidphys );
   return
end

[ mh, message ] = readhd( fidphys, hds );
if ~isempty( mh )
   file.mainheader = mh;
   file.mhmachineformat = mhmachineformat;
   file.filename = filename;
end

fclose( fidphys );

% -----------------------------------------------------

function [ hd, message ] = readhd( fidphys, hds )
message = '';
ok = 1;
hd = [];
initialpos = ftell( fidphys );

hd.file_system = hds.file_system;

for i =1:length( hds.s )
   switch hds.s(i).type
   case 'int16'
      [ val, count ] = fread( fidphys, hds.s(i).noelements, 'int16' );
      hd = setfield( hd, hds.s(i).fieldname, val' );
   case 'int32'
      [ val, count ] = fread( fidphys, hds.s(i).noelements, 'int32' );
      hd = setfield( hd, hds.s(i).fieldname, val' );
   case { 'float32', 'mat1st' }
      [ val, count ] = fread( fidphys, hds.s(i).noelements, 'float32' );
      hd = setfield( hd, hds.s(i).fieldname, val' );
   case 'mat2nd'
      [ val, count ] = fread( fidphys, hds.s(i).noelements, 'float32' );
      mt = getfield( hd, hds.s(i).fieldname );
      mt = [ reshape( mt, 3, 3 )', val ];
      hd = setfield( hd, hds.s(i).fieldname, mt );
   case 'char'
      [ val, count ] = fread( fidphys, hds.s(i).noelements, 'uchar' );
      val = char( val( val ~= 0 ) );
      hd = setfield( hd, hds.s(i).fieldname, val' );
   otherwise
      count = hds.s(i).noelements;
   end % switch hds.s(i).type
   ok = ok & count == hds.s(i).noelements;
   
end
if ~ok
   message = 'Header read error';
end
if ftell( fidphys ) - initialpos ~= hds.header_size
   error( 'Reading header, wrong internal description' )
end

function [ hds, message ] = getmhs( file_system )
message = '';

switch file_system
   
case 'ecat7'
   a = { 'magic_number',                   14, 'char'
      'original_file_name',             32, 'char'
      'sw_version',                      1, 'int16'
      'system_type',                     1, 'int16'
      'file_type',                       1, 'int16'
      'serial_number',                  10, 'char'
      'scan_start_time',                 1, 'int32'
      'isotope_name',                    8, 'char'
      'isotope_halflife',                1, 'float32'
      'radiopharmaceutical',            32, 'char'
      'gantry_tilt',                     1, 'float32'
      'gantry_rotation',                 1, 'float32'
      'bed_elevation',                   1, 'float32'
      'intrinsic_tilt',                  1, 'float32'
      'wobble_speed',                    1, 'int16'
      'transm_source_type',              1, 'int16'
      'distance_scanned',                1, 'float32'
      'transaxial_fov',                  1, 'float32'
      'angular_compression',             1, 'int16'
      'coin_samp_mode',                  1, 'int16'
      'axial_samp_mode',                 1, 'int16'
      'ecat_calibration_factor',         1, 'float32'
      'calibration_units',               1, 'int16'
      'calibration_units_label',         1, 'int16'
      'compression_code',                1, 'int16'
      'study_type',                     12, 'char'
      'patient_id',                     16, 'char'
      'patient_name',                   32, 'char'
      'patient_sex',                     1, 'char'
      'patient_dexterity',               1, 'char'
      'patient_age',                     1, 'float32'
      'patient_height',                  1, 'float32'
      'patient_weight',                  1, 'float32'
      'patient_birth_date',              1, 'int32'
      'physician_name',                 32, 'char'
      'operator_name',                  32, 'char'
      'study_description',              32, 'char'
      'acquisition_type',                1, 'int16'
      'patient_orientation',             1, 'int16'
      'facility_name',                  20, 'char'
      'num_planes',                      1, 'int16'
      'num_frames',                      1, 'int16'
      'num_gates',                       1, 'int16'
      'num_bed_pos',                     1, 'int16'
      'bed_position',                   16, 'float32'
      'plane_separation',                1, 'float32'
      'lwr_sctr_thres',                  1, 'int16'
      'lwr_true_thres',                  1, 'int16'
      'upr_true_thres',                  1, 'int16'
      'user_process_code',              10, 'char'
      'acquisition_mode',                1, 'int16'
      'bin_size',                        1, 'float32'
      'branching_fraction',              1, 'float32'
      'dose_start_time',                 1, 'int32'
      'dosage',                          1, 'float32'
      'well_counter_corr_factor',        1, 'float32'
      'data_units',                     32, 'char'
      'septa_state',                     1, 'int16'
      'fill',                            6, 'int16' };
   
   hds.s = struct( 'fieldname', a(:,1), 'noelements', a(:,2), ...
      'type', a(:,3) );
   hds.header_size = 512;
   hds.file_system = 'ecat7';
   
case 'ecat6.4'
   a = { 'fill1',                          14, 'int16' % changed type
      'original_file_name',             20, 'char' % changed length
      'sw_version',                      1, 'int16'
      'data_type',                       1, 'int16' % ecat6.4
      'system_type',                     1, 'int16'
      'file_type',                       1, 'int16'
      'node_id',                        10, 'char' % serial_number
      ... %'scan_start_time', 1, 'int32'
      'scan_start_day',                  1, 'int16' % ecat6.4
      'scan_start_month',                1, 'int16' % ecat6.4
      'scan_start_year',                 1, 'int16' % ecat6.4
      'scan_start_hour',                 1, 'int16' % ecat6.4
      'scan_start_minute',               1, 'int16' % ecat6.4
      'scan_start_second',               1, 'int16' % ecat6.4
      'isotope_code',                    8, 'char' % isotope_name
      'isotope_halflife',                1, 'float32'
      'radiopharmaceutical',            32, 'char'
      'gantry_tilt',                     1, 'float32'
      'gantry_rotation',                 1, 'float32'
      'bed_elevation',                   1, 'float32'
      ... %'intrinsic_tilt', 1, 'float32'
      'rot_source_speed',                1, 'int16' % ecat6.4
      'wobble_speed',                    1, 'int16'
      'transm_source_type',              1, 'int16'
      ... %'distance_scanned', 1, 'float32'
      'axial_fov',                       1, 'float32'
      'transaxial_fov',                  1, 'float32'
      'transaxial_samp_mode',            1, 'int16' % ecat6.4
      ... %'angular_compression', 1, 'int16'
      'coin_samp_mode',                  1, 'int16'
      'axial_samp_mode',                 1, 'int16' % subtracted 1 in ecat7
      'calibration_factor',              1, 'float32' % ecat_calibration_factor
      'calibration_units',               1, 'int16'
      ... %'calibration_units_label', 1, 'int16'
      'compression_code',                1, 'int16'
      'study_name',                     12, 'char' % study_type
      'patient_id',                     16, 'char'
      'patient_name',                   32, 'char' % changed length
      'patient_sex',                     1, 'char'
      'patient_age',                    10, 'char' % changed type
      'patient_height',                 10, 'char' % changed type
      'patient_weight',                 10, 'char' % changed type
      'patient_dexterity',               1, 'char' % changed position
      ... %'patient_birth_date', 1, 'int32'
      'physician_name',                 32, 'char'
      'operator_name',                  32, 'char'
      'study_description',              32, 'char'
      'acquisition_type',                1, 'int16'
      'bed_type',                        1, 'int16' % ecat6.4
      'septa_type',                      1, 'int16'
      ... %'patient_orientation', 1, 'int16'
      'facility_name',                  20, 'char'
      'num_planes',                      1, 'int16'
      'num_frames',                      1, 'int16'
      'num_gates',                       1, 'int16'
      'num_bed_pos',                     1, 'int16'
      'bed_position',                   16, 'float32'
      'plane_separation',                1, 'float32'
      'lwr_sctr_thres',                  1, 'int16'
      'lwr_true_thres',                  1, 'int16'
      'upr_true_thres',                  1, 'int16'
      'collimator',                      1, 'float32' % ecat6.4
      'user_process_code',              10, 'char'
      ... %'acquisition_mode', 1, 'int16'
      'fill2',                          20, 'int16' };
      %'bin_size', 1, 'float32'
      %'branching_fraction', 1, 'float32'
      %'dose_start_time', 1, 'int32'
      %'dosage', 1, 'float32'
      %'well_counter_corr_factor', 1, 'float32'
      %'data_units', 32, 'char'
      %'septa_state', 1, 'int16'
      %'fill', 6, 'int16' };
   
   hds.s = struct( 'fieldname', a(:,1), 'noelements', a(:,2), ...
      'type', a(:,3) );
   hds.header_size = 512;
   hds.file_system = 'ecat6.4';
otherwise
   message = 'Illegal file system';
   hds = [];

end % switch file_system

function [ hds, message ] = getshs( file_system, file_type );
message = '';
hds = [];

switch file_system
case 'ecat7'
   
   switch file_type
      
      % implemented:
      % 01=Sinogram, 03=Attenuation Correction, 04=Normalization, 
      % 05=Polar Map, 07=Volume 16, 11=3D Sinogram 16, 13=3D Normalization
      
      % not implemented:
      % 00=unknown, 02=Image-16, 06=Volume 8, 08=Projection 8, 
      % 09=Projection 16, 10=Image 8, 12=3D Sinogram 8, 14=3D Sinogram Fit
      
   case 1
      % sinogram
      a = { ...
            'data_type',                       1, 'int16'
         'num_dimensions',                  1, 'int16'
         'num_r_elements',                  1, 'int16'
         'num_angles',                      1, 'int16'
         'corrections_applied',             1, 'int16'
         'num_z_elements',                  1, 'int16'
         'ring_difference',                 1, 'int16'
         'xyz_resolution',                  3, 'float32'
         'w_resolution',                    1, 'float32'
         'fill1',                           6, 'int16'
         'gate_duration',                   1, 'int32'
         'r_wave_offset',                   1, 'int32'
         'num_accepted_beats',              1, 'int32'
         'scale_factor',                    1, 'float32'
         'scan_min',                        1, 'int16'
         'scan_max',                        1, 'int16'
         'prompts',                         1, 'int32'
         'delayed',                         1, 'int32'
         'multiples',                       1, 'int32'
         'net_trues',                       1, 'int32'
         'cor_singles',                    16, 'float32'
         'uncor_singles',                  16, 'float32'
         'tot_avg_cor',                     1, 'float32'
         'tot_avg_uncor',                   1, 'float32'
         'total_coin_rate',                 1, 'int32'
         'frame_start_time',                1, 'int32'
         'frame_duration',                  1, 'int32'
         'deadtime_correction_factor',      1, 'float32'
         'physical_planes',                 8, 'int16'
         'fill2',                          83, 'int16'
         'fill3',                          50, 'int16' };
      
      hds.s = struct( 'fieldname', a(:,1), 'noelements', a(:,2), ...
         'type', a(:,3) );
      hds.header_size = 512;
      
   case 3
      % attenuation
      a = { ...
            'data_type',                       1, 'int16'
         'num_dimensions',                  1, 'int16'
         'attenuation_type',                1, 'int16'
         'num_r_elements',                  1, 'int16'
         'num_angles',                      1, 'int16'
         'num_z_elements',                  1, 'int16'
         'ring_difference',                 1, 'int16'
         'x_resolution',                    1, 'float32'
         'y_resolution',                    1, 'float32'
         'z_resolution',                    1, 'float32'
         'w_resolution',                    1, 'float32'
         'scale_factor',                    1, 'float32'
         'x_offset',                        1, 'float32'
         'y_offset',                        1, 'float32'
         'x_radius',                        1, 'float32'
         'y_radius',                        1, 'float32'
         'tilt_angle',                      1, 'float32'
         'attenuation_coeff',               1, 'float32'
         'attenuation_min',                 1, 'float32'
         'attenuation_max',                 1, 'float32'
         'skull_thickness',                 1, 'float32'
         'num_additional_atten_coeff',      1, 'int16'
         'additional_atten_coeff',          8, 'float32'
         'edge_finding_threshold',          1, 'float32'
         'storage_order',                   1, 'int16'
         'span',                            1, 'int16'
         'z_elements',                     64, 'int16'
         'fill1',                          86, 'int16'
         'fill2',                          50, 'int16' };
      
      hds.s = struct( 'fieldname', a(:,1), 'noelements', a(:,2), ...
         'type', a(:,3) );
      hds.header_size = 512;
      
   case 4
      % 6.5 normalization
      a = { ...
            'data_type',                       1, 'int16'
         'num_dimensions',                  1, 'int16'
         'num_r_elements',                  1, 'int16'
         'num_angles',                      1, 'int16'
         'num_z_elements',                  1, 'int16'
         'ring_difference',                 1, 'int16'
         'scale_factor',                    1, 'float32'
         'norm_min',                        1, 'float32'
         'norm_max',                        1, 'float32'
         'fov_source_width',                1, 'float32'
         'norm_quality_factor',             1, 'float32'
         'norm_quality_factor_code',        1, 'int16'
         'storage_order',                   1, 'int16'
         'span',                            1, 'int16'
         'z_elements',                     64, 'int16'
         'fill1',                         123, 'int16'
         'fill2',                          50, 'int16' };
      
      hds.s = struct( 'fieldname', a(:,1), 'noelements', a(:,2), ...
         'type', a(:,3) );
      hds.header_size = 512;
      
      % polar map
   case 5
      a = { ...
            'data_type',                       1, 'int16'
         'polar_map_type',                  1, 'int16'
         'num_rings',                       1, 'int16'
         'sectors_per_ring',               32, 'int16'
         'ring_position',                  32, 'float32'
         'ring_angle',                     32, 'int16'
         'start_angle',                     1, 'int16'
         'long_axis_left',                  3, 'int16'
         'long_axis_right',                 3, 'int16'
         'position_data',                   1, 'int16'
         'image_min',                       1, 'int16'
         'image_max',                       1, 'int16'
         'scale_factor',                    1, 'float32'
         'pixel_size',                      1, 'float32'
         'frame_duration',                  1, 'int32'
         'frame_start_time',                1, 'int32'
         'processing_code',                 1, 'int16'
         'quant_units',                     1, 'int16'
         'annotation',                     40, 'char'
         'gate_duration',                   1, 'int32'
         'r_wave_offset',                   1, 'int32'
         'num_accepted_beats',              1, 'int32'
         'polar_map_protocol',             20, 'char'
         'database_name',                  30, 'char'
         'fill1',                          27, 'int16'
         'fill2',                          27, 'int16' };
      
      hds.s = struct( 'fieldname', a(:,1), 'noelements', a(:,2), ...
         'type', a(:,3) );
      hds.header_size = 512;
      
   case 7
      % image
      a = { ...
            'data_type',                       1, 'int16'
         'num_dimensions',                  1, 'int16'
         'xyz_dimension',                   3, 'int16'
         'xyz_offset',                      3, 'float32'
         'recon_zoom',                      1, 'float32'
         'scale_factor',                    1, 'float32'
         'image_min',                       1, 'int16'
         'image_max',                       1, 'int16'
         'xyz_pixel_size',                  3, 'float32'
         'frame_duration',                  1, 'int32'
         'frame_start_time',                1, 'int32'
         'filter_code',                     1, 'int16'
         'xyz_resolution',                  3, 'float32'
         'num_r_elements',                  1, 'float32'
         'num_angles',                      1, 'float32'
         'z_rotation_angle',                1, 'float32'
         'decay_corr_fctr',                 1, 'float32'
         'processing_code',                 1, 'int32'
         'gate_duration',                   1, 'int32'
         'r_wave_offset',                   1, 'int32'
         'num_accepted_beats',              1, 'int32'
         'filter_cutoff_frequency',         1, 'float32'
         'filter_resolution',               1, 'float32'
         'filter_ramp_slope',               1, 'float32'
         'filter_order',                    1, 'int16'
         'filter_scatter_fraction',         1, 'float32'
         'filter_scatter_slope',            1, 'float32'
         'annotation',                     40, 'char'
         'mt',                              9, 'mat1st'
         'rfilter_cutoff',                  1, 'float32'
         'rfilter_resolution',              1, 'float32'
         'rfilter_code',                    1, 'int16'
         'rfilter_order',                   1, 'int16'
         'zfilter_cutoff',                  1, 'float32'
         'zfilter_resolution',              1, 'float32'
         'zfilter_code',                    1, 'int16'
         'zfilter_order',                   1, 'int16'
         'mt',                              3, 'mat2nd'
         'scatter_type',                    1, 'int16'
         'recon_type',                      1, 'int16'
         'recon_views',                     1, 'int16'
         'fill',                          136, 'int16' };
      
      hds.s = struct( 'fieldname', a(:,1), 'noelements', a(:,2), ...
         'type', a(:,3) );
      hds.header_size = 512;
      
   case { 11, 14 }
      % 3D sinogram
      a = { ...
            'data_type',                       1, 'int16'
         'num_dimensions',                  1, 'int16'
         'num_r_elements',                  1, 'int16'
         'num_angles',                      1, 'int16'
         'corrections_applied',             1, 'int16'
         'num_z_elements',                 64, 'int16'
         'ring_difference',                 1, 'int16'
         'storage_order',                   1, 'int16'
         'axial_compression',               1, 'int16'
         'x_resolution',                    1, 'float32'
         'v_resolution',                    1, 'float32'
         'z_resolution',                    1, 'float32'
         'w_resolution',                    1, 'float32'
         'fill1',                           6, 'int16'
         'gate_duration',                   1, 'int32'
         'r_wave_offset',                   1, 'int32'
         'num_accepted_beats',              1, 'int32'
         'scale_factor',                    1, 'float32'
         'scan_min',                        1, 'int16'
         'scan_max',                        1, 'int16'
         'prompts',                         1, 'int32'
         'delayed',                         1, 'int32'
         'multiples',                       1, 'int32'
         'net_trues',                       1, 'int32'
         'tot_avg_cor',                     1, 'float32'
         'tot_avg_uncor',                   1, 'float32'
         'total_coin_rate',                 1, 'int32'
         'frame_start_time',                1, 'int32'
         'frame_duration',                  1, 'int32'
         'deadtime_correction_factor',      1, 'float32'
         'fill2',                          90, 'int16'
         'fill3',                          50, 'int16'
         'uncor_singles',                 128, 'float32' };
      
      hds.s = struct( 'fieldname', a(:,1), 'noelements', a(:,2), ...
         'type', a(:,3) );
      hds.header_size = 1024;
      
   case 13
      % 3D normalization
      a = { ...
            'data_type',                       1, 'int16'
         'num_r_elements',                  1, 'int16'
         'num_transaxial_crystals',         1, 'int16'
         'num_crystal_rings',               1, 'int16'
         'crystals_per_ring',               1, 'int16'
         'num_geo_corr_planes',             1, 'int16'
         'uld',                             1, 'int16'
         'lld',                             1, 'int16'
         'scatter_energy',                  1, 'int16'
         'norm_quality_factor',             1, 'float32'
         'norm_quality_factor_code',        1, 'int16'
         'ring_dtcor1',                    32, 'float32'
         'ring_dtcor2',                    32, 'float32'
         'crystal_dtcor',                   8, 'float32'
         'span',                            1, 'int16'
         'max_ring_diff',                   1, 'int16'
         'fill1',                          48, 'int16'
         'fill2',                          50, 'int16' };
      
      hds.s = struct( 'fieldname', a(:,1), 'noelements', a(:,2), ...
         'type', a(:,3) );
      hds.header_size = 512;
   otherwise
      message = 'Illegal file type';
      
   end % switch file_type
   
   hds.file_system = 'ecat7sh';
   
case 'ecat6.4'
   
   switch file_type
      % implemented:
      % 01=.scn, 02=.img, 03=.atn, 04=.nrm 
      
   case 1 % .scn
      
      a = { ...
            'fill1',                          63, 'int16'
         'data_type',                       1, 'int16'
         'fill2',                           2, 'int16'
         'dimension_1',                     1, 'int16'
         'dimension_2',                     1, 'int16'
         'smoothing',                       1, 'int16'
         'processing_code',                 1, 'int16'
         'fill3',                           3, 'int16'
         'sample_distance',                 1, 'float32'
         'fill4',                           8, 'int16'
         'isotope_halflife',                1, 'float32'
         'frame_duration_sec',              1, 'int16'
         'gate_duration',                   1, 'int32'
         'r_wave_offset',                   1, 'int32'
         'fill5',                           1, 'int16'
         'scale_factor',                    1, 'float32'
         'fill6',                           3, 'int16'
         'scan_min',                        1, 'int16'
         'scan_max',                        1, 'int16'
         'prompts',                         1, 'int32'
         'delayed',                         1, 'int32'
         'multiples',                       1, 'int32'
         'net_trues',                       1, 'int32'
         'fill7',                          52, 'int16'
         'cor_singles',                    16, 'float32'
         'uncor_singles',                  16, 'float32'
         'tot_avg_cor',                     1, 'float32'
         'tot_avg_uncor',                   1, 'float32'
         'total_coin_rate',                 1, 'int32'
         'frame_start_time',                1, 'int32'
         'frame_duration',                  1, 'int32'
         'loss_correction_fctr'             1, 'float32'
         'fill8',                          22, 'int16' };
      
      hds.s = struct( 'fieldname', a(:,1), 'noelements', a(:,2), ...
         'type', a(:,3) );
      hds.header_size = 512;
      
   case 2 % .img
      
      a = { ...
            'fill1',                          63, 'int16'
         'data_type',                       1, 'int16'
         'num_dimensions',                  1, 'int16'
         'fill2',                           1, 'int16'
         'dimension_1',                     1, 'int16'
         'dimension_2',                     1, 'int16'
         'fill3',                          12, 'int16'
         'x_origin',                        1, 'float32'
         'y_origin',                        1, 'float32'
         'recon_scale',                     1, 'float32'
         'quant_scale',                     1, 'float32'
         'image_min',                       1, 'int16'
         'image_max',                       1, 'int16'
         'fill4',                           2, 'int16'
         'pixel_size',                      1, 'float32'
         'slice_width',                     1, 'float32'
         'frame_duration',                  1, 'int32'
         'frame_start_time',                1, 'int32'
         'slice_location',                  1, 'int16'
         'recon_start_hour',                1, 'int16'
         'recon_start_min',                 1, 'int16'
         'recon_start_sec',                 1, 'int16'
         'recon_duration',                  1, 'int32'
         'fill5',                          12, 'int16'
         'filter_code',                     1, 'int16'
         'scan_matrix_num',                 1, 'int32'
         'norm_matrix_num',                 1, 'int32'
         'atten_cor_mat_num',               1, 'int32'
         'fill6',                          23, 'int16'
         'image_rotation',                  1, 'float32'
         'plane_eff_corr_fctr',             1, 'float32'
         'decay_corr_fctr',                 1, 'float32'
         'loss_corr_fctr',                  1, 'float32'
         'fill7',                          32, 'int16'
         'processing_code',                 1, 'int16'
         'fill8',                           1, 'int16'
         'quant_units',                     1, 'int16'
         'recon_start_day',                 1, 'int16'
         'recon_start_month',               1, 'int16'
         'recon_start_year',                1, 'int16'
         'ecat_calibration_fctr',           1, 'float32'
         'well_counter_cal_fctr',           1, 'float32'
         'filter_params',                   6, 'float32'
         'annotation',                     40, 'char'
         'fill9',                          26, 'int16' };
         
         hds.s = struct( 'fieldname', a(:,1), 'noelements', a(:,2), ...
            'type', a(:,3) );
         hds.header_size = 512;
      
   case 3 % .atn
      
      a = { ...
            'fill1',                          63, 'int16'
         'data_type',                       1, 'int16'
         'attenuation_type',                1, 'int16'
         'fill2',                           1, 'int16'
         'dimension_1',                     1, 'int16'
         'dimension_2',                     1, 'int16'
         'fill3',                          23, 'int16'
         'scale_factor',                    1, 'float32'
         'x_origin',                        1, 'float32'
         'y_origin',                        1, 'float32'
         'x_radius',                        1, 'float32'
         'y_radius',                        1, 'float32'
         'tilt_angle',                      1, 'float32'
         'attenuation_coeff',               1, 'float32'
         'sample_distance',                 1, 'float32'
         'fill4',                         149, 'int16' };
         
         hds.s = struct( 'fieldname', a(:,1), 'noelements', a(:,2), ...
            'type', a(:,3) );
         hds.header_size = 512;
      
   case 4 % .nrm
      
      a = { ...
            'fill1',                          63, 'int16'
         'data_type',                       1, 'int16'
         'fill2',                           2, 'int16'
         'dimension_1',                     1, 'int16'
         'dimension_2',                     1, 'int16'
         'fill3',                          23, 'int16'
         'scale_factor',                    1, 'float32'
         'norm_hour',                       1, 'int16'
         'norm_minute',                     1, 'int16'
         'norm_second',                     1, 'int16'
         'norm_day',                        1, 'int16'
         'norm_month',                      1, 'int16'
         'norm_year',                       1, 'int16'
         'fov_source_width',                1, 'float32'
         'fill4',                         155, 'int16' };
      
      hds.s = struct( 'fieldname', a(:,1), 'noelements', a(:,2), ...
         'type', a(:,3) );
      hds.header_size = 512;
   otherwise
      message = 'Illegal file type';
      
   end % switch file_type
   
   hds.file_system = 'ecat6.4sh';
   
otherwise
   message = 'Illegal file system';
   
end % switch file_system

function [ datadescr, message ] = getdatadescr( file_system, file_type, data_type );

datadescr = [];
message = '';

switch file_system
   
case 'ecat6.4'
   switch file_type
      
   case { 1, 2, 3, 4 }
      % 1 = .scn, 2 = .img, 3 = .atn, 4 = .nrm
      switch data_type
      case 2 % VAX_Ix2
         datadescr.ftype = 'int16';
         datadescr.size = 2;
      case 4 % vax float
         datadescr.ftype = 'float32';
         datadescr.size = 4;
      case 5 % ieee float
         datadescr.ftype = 'float32';
         datadescr.size = 4;
      otherwise
         message = 'Data type not implemented';
      end
   otherwise 
      message = 'File type not implemented';
   end %switch file_type
   
case 'ecat7'
   
   % implemented:
   % 01=Sinogram, 03=Attenuation Correction, 04=Normalization, 
   % 05=Polar Map, 07=Volume 16, 11=3D Sinogram 16, 13=3D Normalization
   
   % not implemented:
   % 00=unknown, 02=Image-16, 06=Volume 8, 08=Projection 8, 
   % 09=Projection 16, 10=Image 8, 12=3D Sinogram 8, 14=3D Sinogram Fit
   switch file_type
   case { 11, 14 }
      switch data_type
      case 1 % byte
         datadescr.ftype = 'int8';
         datadescr.size = 1;
      case { 2, 6 } % Sun short, sun_in
         datadescr.ftype = 'int16';
         datadescr.size = 2;
      case 5 % ieee float
         datadescr.ftype = 'float32';
         datadescr.size = 4;
      otherwise
         message = 'Data type not implemented';
      end   
   case 13
      switch data_type
      case 1 % ieee-float
         datadescr.ftype = 'float32';
         datadescr.size = 4;
      case 5 % sun-float
         datadescr.ftype = 'float32';
         datadescr.size = 4;
      otherwise
         message = 'Data type not implemented';
      end   
   case { 1, 3, 4, 5, 7 }
      % (0=Unkonwn Matrix Data Type, 1=Byte Data, 2=VAX_Ix2, 3=VAX_Ix4, 
      %  4=VAX_Rx4, 5=IEEE Float, 6=Sun short, 7=Sun long) 
      % DTYPE_BYTES, _I2, _I4, _VAXR4, _SUNFL, _SUNIN 
      switch data_type    
      case 1 % byte
         datadescr.ftype = 'int8';
         datadescr.size = 1;
      case 6 % Sun short
         datadescr.ftype = 'int16';
         datadescr.size = 2;
      case 5 % ieee float
         datadescr.ftype = 'float32';
         datadescr.size = 4;
      otherwise
         message = 'Data type not implemented';
      end
   otherwise 
      error( 'File type not implemented' );
   end %switch file_type
   
otherwise
   message = 'File system not implemented';
   
end % switch file_system
% ------------------------------------------------------------

function [ hd, offsetfromheader, writefill, numdatablocks, message ] = ...
      set_xyz_dimension( hd, datadescr, xyz_dimension, selsegment  );
% implemented:
% 01=Sinogram, 03=Attenuation Correction, 04=Normalization,
% 05=Polar Map, 07=Volume 16, 11=3D Sinogram 16, 13=3D Normalization

message = '';
offsetfromheader = 0;
writefill = 0;
numdatablocks = 0;

switch hd.mh.file_system
case 'ecat6.4'
   
   switch hd.mh.file_type
   case { 01, 02, 03, 04 } % .scn, .img, .atn, .nrm
      hd.sh.dimension_1 = xyz_dimension( 1 );
      hd.sh.dimension_2 = xyz_dimension( 2 );
      numdatablocks = floor( ...
         ( prod( xyz_dimension ) * datadescr.size + 511 ) / 512 );
      writefill = ( 511 - rem( prod( xyz_dimension ) * datadescr.size + 511, 512 ) );
      
   otherwise
      message = 'Writing of the file type not implemented';
      hd = [];
      
   end
   
case 'ecat7'
   
   switch hd.mh.file_type
   case 07 % Volume 16
      hd.sh.xyz_dimension = xyz_dimension;
      numdatablocks = floor( ...
         ( prod( xyz_dimension ) * datadescr.size + 511 ) / 512 );
      writefill = ( 511 - rem( prod( xyz_dimension ) * datadescr.size + 511, 512 ) );
      
   case { 11, 14 } % 3D Sinogram 16, FORE sinogram float
      if hd.sh.storage_order == 0
         hd.sh.num_r_elements = xyz_dimension( 1 );
         hd.sh.num_z_elements( selsegment ) = xyz_dimension( 2 );
         hd.sh.num_angles = xyz_dimension( 3 );
      else
         hd.sh.num_r_elements = xyz_dimension( 1 );
         hd.sh.num_angles = xyz_dimension( 2 );
         hd.sh.num_z_elements( selsegment ) = xyz_dimension( 3 );
      end
      
      nobytes = hd.sh.num_r_elements * hd.sh.num_angles * datadescr.size * ...
         hd.sh.num_z_elements;
      %   numdatablocksarr = floor( ( nobytes + 511 ) / 512 );
      %   offsetfromheader = sum( numdatablocksarr( 1:selsegment-1 ) ) * 512;
      offsetfromheader = sum( nobytes( 1:selsegment-1 ) );
      
      %   numdatablocks = sum( numdatablocksarr );
      numdatablocks = floor( ( sum( nobytes ) + 511 ) / 512 );
      
      if sum( hd.sh.num_z_elements( selsegment + 1 : end ) ) == 0
         writefill = ( 511 - rem( prod( xyz_dimension ) * datadescr.size + offsetfromheader + 511, 512 ) );
      else
         writefill = 0;
      end
      
   case 3 % Attenuation
      if hd.sh.storage_order == 0
         hd.sh.num_r_elements = xyz_dimension( 1 );
         hd.sh.z_elements( selsegment ) = xyz_dimension( 2 );
         hd.sh.num_angles = xyz_dimension( 3 );
      else
         hd.sh.num_r_elements = xyz_dimension( 1 );
         hd.sh.num_angles = xyz_dimension( 2 );
         hd.sh.z_elements( selsegment ) = xyz_dimension( 3 );
      end
      if selsegment == 1
         hd.sh.num_z_elements = hd.sh.z_elements( 1 );
      end
      
      nobytes = hd.sh.num_r_elements * hd.sh.num_angles * datadescr.size * ...
         hd.sh.z_elements;
      %   numdatablocksarr = floor( ( nobytes + 511 ) / 512 );
      %   offsetfromheader = sum( numdatablocksarr( 1:selsegment-1 ) ) * 512;
      offsetfromheader = sum( nobytes( 1:selsegment-1 ) );
      
      %   numdatablocks = sum( numdatablocksarr );
      numdatablocks = floor( ( sum( nobytes ) + 511 ) / 512 );
      
      if sum( hd.sh.num_z_elements( selsegment + 1 : end ) ) == 0
         writefill = ( 511 - rem( prod( xyz_dimension ) * datadescr.size + offsetfromheader + 511, 512 ) );
      else
         writefill = 0;
      end
      
   otherwise
      message = 'Writing of the file type not implemented';
      hd = [];
      
   end
   
otherwise
   message = 'File system not implemented';
   hd = [];
   
end % switch hd.mh.file_system

% ------------------------------------------------------------
function [ xyz_dimension, offsetfromheader, message ] = get_xyz_dimension( hd, datadescr, selsegment )
% implemented:
% 01=Sinogram, 03=Attenuation Correction, 04=Normalization, 
% 05=Polar Map, 07=Volume 16, 11=3D Sinogram 16, 13=3D Normalization

message = '';
xyz_dimension = [];
offsetfromheader = 0;

switch hd.mh.file_system
case 'ecat7'
   
   switch hd.mh.file_type
   case 07 % Volume 16
      xyz_dimension = hd.sh.xyz_dimension;
      
   case { 11, 14 } % 3D Sinogram 16, FORE sinogram float
      if hd.sh.storage_order == 0
         xyz_dimension = [ hd.sh.num_r_elements, ...
               hd.sh.num_z_elements( selsegment ), hd.sh.num_angles ];
      else
         xyz_dimension = [ hd.sh.num_r_elements, hd.sh.num_angles, ...
               hd.sh.num_z_elements( selsegment ) ];
      end
      nobytes = hd.sh.num_r_elements * hd.sh.num_angles * datadescr.size * ...
         hd.sh.num_z_elements;
      %   numdatablocksarr = floor( ( nobytes + 511 ) / 512 );
      %   offsetfromheader = sum( numdatablocksarr( 1:selsegment-1 ) ) * 512;
      offsetfromheader = sum( nobytes( 1:selsegment-1 ) );
      
   case 3 % Attenuation
      if hd.sh.storage_order == 0
         xyz_dimension = [ hd.sh.num_r_elements, ...
               hd.sh.z_elements( selsegment ), hd.sh.num_angles ];
      else
         xyz_dimension = [ hd.sh.num_r_elements, hd.sh.num_angles, ...
               hd.sh.z_elements( selsegment ) ];
      end
      nobytes = hd.sh.num_r_elements * hd.sh.num_angles * datadescr.size * ...
         hd.sh.z_elements;
      %   numdatablocksarr = floor( ( nobytes + 511 ) / 512 );
      %   offsetfromheader = sum( numdatablocksarr( 1:selsegment-1 ) ) * 512;
      offsetfromheader = sum( nobytes( 1:selsegment-1 ) );
      
   otherwise
      message = 'Reading of the file type not implemented';
      
   end
   
case 'ecat6.4'
   
   switch hd.mh.file_type
   case { 01, 02, 03, 04 } % .scn, .img, .atn, .nrm
      xyz_dimension = [ hd.sh.dimension_1, hd.sh.dimension_2, 1 ];
      
   otherwise
      message = 'Reading of the file type not implemented';
      
   end
   
otherwise
   message = 'File system not implemented';
   
end % switch hd.mh.file_system

% ------------------------------------------------------------

function n = firstno( sel )
n = [ find( sel(:) ); 0 ];
n = n( 1 );

% ------------------------------------------------------------



Contact us