ecatfile.m

医学图像处理matlab工具箱

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'