ecatfile.m

医学图像处理matlab工具箱

   [ 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 = [];