2023_Blaschke_Stroke_tDCS_rsfMRI/scripts/toolbox/spm12/spm_ecat2nifti.m

function N = spm_ecat2nifti(fname,opts)
% Import ECAT 7 images from CTI PET scanners
% FORMAT N = spm_ecat2nifti(fname)
% fname    - name of ECAT file
% opts     - options structure
%
% N        - NIfTI object (written in current directory)
%__________________________________________________________________________
% Copyright (C) 2005-2015 Wellcome Trust Centre for Neuroimaging

% John Ashburner & Roger Gunn
% $Id: spm_ecat2nifti.m 6315 2015-01-23 17:09:06Z guillaume $


if nargin==1
    opts = struct('ext',spm_file_ext);
else
    if opts.ext(1) ~= '.', opts.ext = ['.' opts.ext]; end
end

fp = fopen(fname,'r','ieee-be');
if fp == -1,
    error(['Can''t open "' fname '".']);
end

mh = ECAT7_mheader(fp);
if ~strcmp(mh.MAGIC_NUMBER,'MATRIX70v') &&...
        ~strcmp(mh.MAGIC_NUMBER,'MATRIX71v') &&...
        ~strcmp(mh.MAGIC_NUMBER,'MATRIX72v')
    fclose(fp);
    error(['"' fname '" does not appear to be ECAT 7 format.']);
end

if mh.FILE_TYPE ~= 7
    fclose(fp);
    error(['"' fname '" does not appear to be an image file.']);
end

list    = s7_matlist(fp);
matches = find((list(:,4) == 1) | (list(:,4) == 2));
llist   = list(matches,:);

for i=1:size(llist,1)
    sh(i) = ECAT7_sheader(fp,llist(i,2));
end
fclose(fp);

for i=1:size(llist,1)
    dim      = [sh(i).X_DIMENSION sh(i).Y_DIMENSION sh(i).Z_DIMENSION];
    dtype    = [4 1];
    off      = 512*llist(i,2);
    scale    = sh(i).SCALE_FACTOR*mh.ECAT_CALIBRATION_FACTOR;
    inter    = 0;
    dati     = file_array(fname,dim,dtype,off,scale,inter);

    dircos   = diag([-1 -1 -1]);
    step     = ([sh(i).X_PIXEL_SIZE sh(i).Y_PIXEL_SIZE sh(i).Z_PIXEL_SIZE]*10);
    start    = -(dim(1:3)'/2).*step';
    mat      = [[dircos*diag(step) dircos*start] ; [0 0 0 1]];

    matnum   = sprintf('%.8x',list(i,1));
    nam      = spm_file(fname,'basename');
    fnameo   = fullfile(pwd,[nam '_' matnum opts.ext]);
    dato     = file_array(fnameo,dim,[4 spm_platform('bigend')],0,scale,inter);

    N        = nifti;
    N.dat    = dato;
    N.mat    = mat;
    N.mat0   = mat;
    N.mat_intent  = 'aligned';
    N.mat0_intent = 'scanner';
    N.descrip     = sh(i).ANNOTATION;
    N.timing = struct('toffset',sh(i).FRAME_START_TIME/1000,'tspace',sh(i).FRAME_DURATION/1000);
    create(N);
    for j=1:dim(3)
        N.dat(:,:,j) = dati(:,:,j);
    end
    N.extras     = struct('mh',mh,'sh',sh(i));
end


%==========================================================================
% function list = s7_matlist(fp)
%==========================================================================
function list = s7_matlist(fp)
%S7_MATLIST List the available matrixes in an ECAT 7 file.
%   LIST = S7_MATLIST(FP) lists the available matrixes
%   in the file specified by FP.
%
%   Columns in LIST:
%       1 - Matrix identifier.
%       2 - Matrix subheader record number
%       3 - Last record number of matrix data block.
%       4 - Matrix status:
%           1 - exists - rw
%           2 - exists - ro
%           3 - matrix deleted
%

% I believe fp should be opened with:
%     fp = fopen(filename,'r','ieee-be');

fseek(fp,512,'bof');
block = fread(fp,128,'int');
if size(block,1) ~= 128
    list = [];
    return;
end
block = reshape(block,4,32);
list  = [];
while block(2,1) ~= 2
    if block(1,1)+block(4,1) ~= 31
        list = []; return;
    end
    list = [list block(:,2:32)];

    fseek(fp,512*(block(2,1)-1),'bof');
    block = fread(fp,128,'int');
    if size(block,1) ~= 128, list = []; return; end
    block = reshape(block,4,32);
end
list = [list block(:,2:(block(4,1)+1))];
list = list';

%==========================================================================
% function SHEADER=ECAT7_sheader(fid,record)
%==========================================================================
function SHEADER=ECAT7_sheader(fid,record)
%
% Sub header read routine for ECAT 7 image files
%
% Roger Gunn, 260298
off                     = (record-1)*512;
status                  = fseek(fid, off,'bof');
data_type               = fread(fid,1,'uint16',0);
num_dimensions          = fread(fid,1,'uint16',0);
x_dimension             = fread(fid,1,'uint16',0);
y_dimension             = fread(fid,1,'uint16',0);
z_dimension             = fread(fid,1,'uint16',0);
x_offset                = fread(fid,1,'float32',0);
y_offset                = fread(fid,1,'float32',0);
z_offset                = fread(fid,1,'float32',0);
recon_zoom              = fread(fid,1,'float32',0);
scale_factor            = fread(fid,1,'float32',0);
image_min               = fread(fid,1,'int16',0);
image_max               = fread(fid,1,'int16',0);
x_pixel_size            = fread(fid,1,'float32',0);
y_pixel_size            = fread(fid,1,'float32',0);
z_pixel_size            = fread(fid,1,'float32',0);
frame_duration          = fread(fid,1,'uint32',0);
frame_start_time        = fread(fid,1,'uint32',0);
filter_code             = fread(fid,1,'uint16',0);
x_resolution            = fread(fid,1,'float32',0);
y_resolution            = fread(fid,1,'float32',0);
z_resolution            = fread(fid,1,'float32',0);
num_r_elements          = fread(fid,1,'float32',0);
num_angles              = fread(fid,1,'float32',0);
z_rotation_angle        = fread(fid,1,'float32',0);
decay_corr_fctr         = fread(fid,1,'float32',0);
corrections_applied     = fread(fid,1,'uint32',0);
gate_duration           = fread(fid,1,'uint32',0);
r_wave_offset           = fread(fid,1,'uint32',0);
num_accepted_beats      = fread(fid,1,'uint32',0);
filter_cutoff_frequency = fread(fid,1,'float32',0);
filter_resolution       = fread(fid,1,'float32',0);
filter_ramp_slope       = fread(fid,1,'float32',0);
filter_order            = fread(fid,1,'uint16',0);
filter_scatter_fraction = fread(fid,1,'float32',0);
filter_scatter_slope    = fread(fid,1,'float32',0);
annotation              = fread(fid,40,'char',0);
mt_1_1                  = fread(fid,1,'float32',0);
mt_1_2                  = fread(fid,1,'float32',0);
mt_1_3                  = fread(fid,1,'float32',0);
mt_2_1                  = fread(fid,1,'float32',0);
mt_2_2                  = fread(fid,1,'float32',0);
mt_2_3                  = fread(fid,1,'float32',0);
mt_3_1                  = fread(fid,1,'float32',0);
mt_3_2                  = fread(fid,1,'float32',0);
mt_3_3                  = fread(fid,1,'float32',0);
rfilter_cutoff          = fread(fid,1,'float32',0);
rfilter_resolution      = fread(fid,1,'float32',0);
rfilter_code            = fread(fid,1,'uint16',0);
rfilter_order           = fread(fid,1,'uint16',0);
zfilter_cutoff          = fread(fid,1,'float32',0);
zfilter_resolution      = fread(fid,1,'float32',0);
zfilter_code            = fread(fid,1,'uint16',0);
zfilter_order           = fread(fid,1,'uint16',0);
mt_4_1                  = fread(fid,1,'float32',0);
mt_4_2                  = fread(fid,1,'float32',0);
mt_4_3                  = fread(fid,1,'float32',0);
scatter_type            = fread(fid,1,'uint16',0);
recon_type              = fread(fid,1,'uint16',0);
recon_views             = fread(fid,1,'uint16',0);
fill                    = fread(fid,1,'uint16',0);
annotation              = deblank(char(annotation.*(annotation>0))');

SHEADER = struct('DATA_TYPE', data_type, ...
    'NUM_DIMENSIONS', num_dimensions, ...
    'X_DIMENSION', x_dimension, ...
    'Y_DIMENSION', y_dimension, ...
    'Z_DIMENSION', z_dimension, ...
    'X_OFFSET', x_offset, ...
    'Y_OFFSET', y_offset, ...
    'Z_OFFSET', z_offset, ...
    'RECON_ZOOM', recon_zoom, ...
    'SCALE_FACTOR', scale_factor, ...
    'IMAGE_MIN', image_min, ...
    'IMAGE_MAX', image_max, ...
    'X_PIXEL_SIZE', x_pixel_size, ...
    'Y_PIXEL_SIZE', y_pixel_size, ...
    'Z_PIXEL_SIZE', z_pixel_size, ...
    'FRAME_DURATION', frame_duration, ...
    'FRAME_START_TIME', frame_start_time, ...
    'FILTER_CODE', filter_code, ...
    'X_RESOLUTION', x_resolution, ...
    'Y_RESOLUTION', y_resolution, ...
    'Z_RESOLUTION', z_resolution, ...
    'NUM_R_ELEMENTS', num_r_elements, ...
    'NUM_ANGLES', num_angles, ...
    'Z_ROTATION_ANGLE', z_rotation_angle, ...
    'DECAY_CORR_FCTR', decay_corr_fctr, ...
    'CORRECTIONS_APPLIED', corrections_applied, ...
    'GATE_DURATION', gate_duration, ...
    'R_WAVE_OFFSET', r_wave_offset, ...
    'NUM_ACCEPTED_BEATS', num_accepted_beats, ...
    'FILTER_CUTOFF_FREQUENCY', filter_cutoff_frequency, ...
    'FILTER_RESOLUTION', filter_resolution, ...
    'FILTER_RAMP_SLOPE', filter_ramp_slope, ...
    'FILTER_ORDER', filter_order, ...
    'FILTER_SCATTER_CORRECTION', filter_scatter_fraction, ...
    'FILTER_SCATTER_SLOPE', filter_scatter_slope, ...
    'ANNOTATION', annotation, ...
    'MT_1_1', mt_1_1, ...
    'MT_1_2', mt_1_2, ...
    'MT_1_3', mt_1_3, ...
    'MT_2_1', mt_2_1, ...
    'MT_2_2', mt_2_2, ...
    'MT_2_3', mt_2_3, ...
    'MT_3_1', mt_3_1, ...
    'MT_3_2', mt_3_2, ...
    'MT_3_3', mt_3_3, ...
    'RFILTER_CUTOFF', rfilter_cutoff, ...
    'RFILTER_RESOLUTION', rfilter_resolution, ...
    'RFILTER_CODE', rfilter_code, ...
    'RFILTER_ORDER', rfilter_order, ...
    'ZFILTER_CUTOFF', zfilter_cutoff, ...
    'ZFILTER_RESOLUTION', zfilter_resolution, ...
    'ZFILTER_CODE', zfilter_code, ...
    'ZFILTER_ORDER', zfilter_order, ...
    'MT_4_1', mt_4_1, ...
    'MT_4_2', mt_4_2, ...
    'MT_4_3', mt_4_3, ...
    'SCATTER_TYPE', scatter_type, ...
    'RECON_TYPE', recon_type, ...
    'RECON_VIEWS', recon_views, ...
    'FILL', fill);

%==========================================================================
% function [MHEADER]=ECAT7_mheader(fid)
%==========================================================================
function [MHEADER]=ECAT7_mheader(fid)
%
% Main header read routine for ECAT 7 image files
%
% Roger Gunn, 260298

status                  = fseek(fid, 0,'bof');
magic_number            = fread(fid,14,'char',0);
original_file_name      = fread(fid,32,'char',0);
sw_version              = fread(fid,1,'uint16',0);
system_type             = fread(fid,1,'uint16',0);
file_type               = fread(fid,1,'uint16',0);
serial_number           = fread(fid,10,'char',0);
scan_start_time         = fread(fid,1,'uint32',0);
isotope_name            = fread(fid,8,'char',0);
isotope_halflife        = fread(fid,1,'float32',0);
radiopharmaceutical     = fread(fid,32,'char',0);
gantry_tilt             = fread(fid,1,'float32',0);
gantry_rotation         = fread(fid,1,'float32',0);
bed_elevation           = fread(fid,1,'float32',0);
intrinsic_tilt           = fread(fid,1,'float32',0);
wobble_speed             = fread(fid,1,'uint16',0);
transm_source_type       = fread(fid,1,'uint16',0);
distance_scanned         = fread(fid,1,'float32',0);
transaxial_fov           = fread(fid,1,'float32',0);
angular_compression      = fread(fid,1,'uint16',0);
coin_samp_mode           = fread(fid,1,'uint16',0);
axial_samp_mode          = fread(fid,1,'uint16',0);
ecat_calibration_factor  = fread(fid,1,'float32',0);
calibration_units        = fread(fid,1,'uint16',0);
calibration_units_type   = fread(fid,1,'uint16',0);
compression_code         = fread(fid,1,'uint16',0);
study_type               = fread(fid,12,'char',0);
patient_id               = fread(fid,16,'char',0);
patient_name             = fread(fid,32,'char',0);
patient_sex              = fread(fid,1,'char',0);
patient_dexterity        = fread(fid,1,'char',0);
patient_age              = fread(fid,1,'float32',0);
patient_height           = fread(fid,1,'float32',0);
patient_weight           = fread(fid,1,'float32',0);
patient_birth_date       = fread(fid,1,'uint32',0);
physician_name           = fread(fid,32,'char',0);
operator_name            = fread(fid,32,'char',0);
study_description        = fread(fid,32,'char',0);
acquisition_type         = fread(fid,1,'uint16',0);
patient_orientation      = fread(fid,1,'uint16',0);
facility_name            = fread(fid,20,'char',0);
num_planes               = fread(fid,1,'uint16',0);
num_frames               = fread(fid,1,'uint16',0);
num_gates                = fread(fid,1,'uint16',0);
num_bed_pos              = fread(fid,1,'uint16',0);
init_bed_position        = fread(fid,1,'float32',0);
bed_position = zeros(15,1);
for bed=1:15
    tmp = fread(fid,1,'float32',0);
    if ~isempty(tmp), bed_position(bed) = tmp; end
end
plane_separation         = fread(fid,1,'float32',0);
lwr_sctr_thres           = fread(fid,1,'uint16',0);
lwr_true_thres           = fread(fid,1,'uint16',0);
upr_true_thres           = fread(fid,1,'uint16',0);
user_process_code        = fread(fid,10,'char',0);
acquisition_mode         = fread(fid,1,'uint16',0);
bin_size                 = fread(fid,1,'float32',0);
branching_fraction       = fread(fid,1,'float32',0);
dose_start_time          = fread(fid,1,'uint32',0);
dosage                   = fread(fid,1,'float32',0);
well_counter_corr_factor = fread(fid,1,'float32',0);
data_units               = fread(fid,32,'char',0);
septa_state              = fread(fid,1,'uint16',0);
fill                     = fread(fid,1,'uint16',0);
    
magic_number             = deblank(char(magic_number.*(magic_number>32))');
original_file_name       = deblank(char(original_file_name.*(original_file_name>0))');
serial_number            = deblank(char(serial_number.*(serial_number>0))');
isotope_name             = deblank(char(isotope_name.*(isotope_name>0))');
radiopharmaceutical      = deblank(char(radiopharmaceutical.*(radiopharmaceutical>0))');
study_type               = deblank(char(study_type.*(study_type>0))');
patient_id               = deblank(char(patient_id.*(patient_id>0))');
patient_name             = deblank(char(patient_name.*(patient_name>0))');
patient_sex              = deblank(char(patient_sex.*(patient_sex>0))');
patient_dexterity        = deblank(char(patient_dexterity.*(patient_dexterity>0))');
physician_name           = deblank(char(physician_name.*(physician_name>0))');
operator_name            = deblank(char(operator_name.*(operator_name>0))');
study_description        = deblank(char(study_description.*(study_description>0))');
facility_name            = deblank(char(facility_name.*(facility_name>0))');
user_process_code        = deblank(char(user_process_code.*(user_process_code>0))');
data_units               = deblank(char(data_units.*(data_units>0))');

MHEADER = struct('MAGIC_NUMBER', magic_number, ...
    'ORIGINAL_FILE_NAME', original_file_name, ...
    'SW_VERSION', sw_version, ...
    'SYSTEM_TYPE', system_type, ...
    'FILE_TYPE', file_type, ...
    'SERIAL_NUMBER', serial_number, ...
    'SCAN_START_TIME', scan_start_time, ...
    'ISOTOPE_NAME', isotope_name, ...
    'ISOTOPE_HALFLIFE', isotope_halflife, ...
    'RADIOPHARMACEUTICAL', radiopharmaceutical, ... 
    'GANTRY_TILT', gantry_tilt, ...
    'GANTRY_ROTATION', gantry_rotation, ... 
    'BED_ELEVATION', bed_elevation, ...
    'INTRINSIC_TILT', intrinsic_tilt, ... 
    'WOBBLE_SPEED', wobble_speed, ...
    'TRANSM_SOURCE_TYPE', transm_source_type, ...
    'DISTANCE_SCANNED', distance_scanned, ...
    'TRANSAXIAL_FOV', transaxial_fov, ...
    'ANGULAR_COMPRESSION', angular_compression, ... 
    'COIN_SAMP_MODE', coin_samp_mode, ...
    'AXIAL_SAMP_MODE', axial_samp_mode, ...
    'ECAT_CALIBRATION_FACTOR', ecat_calibration_factor, ...
    'CALIBRATION_UNITS', calibration_units, ...
    'CALIBRATION_UNITS_TYPE', calibration_units_type, ...
    'COMPRESSION_CODE', compression_code, ...
    'STUDY_TYPE', study_type, ...
    'PATIENT_ID', patient_id, ...
    'PATIENT_NAME', patient_name, ...
    'PATIENT_SEX', patient_sex, ...
    'PATIENT_DEXTERITY', patient_dexterity, ...
    'PATIENT_AGE', patient_age, ...
    'PATIENT_HEIGHT', patient_height, ...
    'PATIENT_WEIGHT', patient_weight, ...
    'PATIENT_BIRTH_DATE', patient_birth_date, ...
    'PHYSICIAN_NAME', physician_name, ...
    'OPERATOR_NAME', operator_name, ...
    'STUDY_DESCRIPTION', study_description, ...
    'ACQUISITION_TYPE', acquisition_type, ...
    'PATIENT_ORIENTATION', patient_orientation, ...
    'FACILITY_NAME', facility_name, ...
    'NUM_PLANES', num_planes, ...
    'NUM_FRAMES', num_frames, ...
    'NUM_GATES', num_gates, ...
    'NUM_BED_POS', num_bed_pos, ...
    'INIT_BED_POSITION', init_bed_position, ...
    'BED_POSITION', bed_position, ...
    'PLANE_SEPARATION', plane_separation, ...
    'LWR_SCTR_THRES', lwr_sctr_thres, ...
    'LWR_TRUE_THRES', lwr_true_thres, ...
    'UPR_TRUE_THRES', upr_true_thres, ...
    'USER_PROCESS_CODE', user_process_code, ...
    'ACQUISITION_MODE', acquisition_mode, ...
    'BIN_SIZE', bin_size, ...
    'BRANCHING_FRACTION', branching_fraction, ...
    'DOSE_START_TIME', dose_start_time, ...
    'DOSAGE', dosage, ...
    'WELL_COUNTER_CORR_FACTOR', well_counter_corr_factor, ...
    'DATA_UNITS', data_units, ...
    'SEPTA_STATE', septa_state, ...
    'FILL', fill);