ecatfile.m

医学图像处理matlab工具箱

               '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