ecatfile.m

医学图像处理matlab工具箱

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.

% Copyright (C) 1999, 2000, 2001 by Flemming Hermansen.

% Disclaimer, copyright & licencing
% ecatfile is free but copyright software, distributed under the
% terms of the GNU General Public Licence as published by the Free
% Software Foundation (either version 2, or
% at your option, any later version). Further details on "copyleft" can
% be found at http://www.gnu.org/copyleft/. In particular, ecatfile is
% supplied as is. No formal support or maintenance is provided or
% implied. The author will have no responsibility of the use and possible errors.

% 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 ), ...