Hi everyone,
I’m new to neurodata and I am learn to preprocess 18F-FDG Pet images for deep learning models. Raw data (left one) n output (right one) look like this:
But after preprocessing (spatial norm using MNI template + intensity normalization + gaussian smoothing), they are all look wrong, can you tell me help me identify what is the problem with my matlab code, and should i use spm matlab for preprocessing PET data or is there any robust python toolbox for preprocessing PET data, here is my code. Thank you so much for your help.
function preprocess_pet_batch(root_dir, output_dir)
spm(‘defaults’, ‘PET’);spm_jobman(‘initcfg’);
template = fullfile(spm(‘Dir’), ‘toolbox’, ‘OldNorm’, ‘PET.nii’);if ~exist(template, ‘file’)error(‘SPM PET template not found at:\n %s\nCheck your SPM installation.’, template);end
if ~exist(output_dir, ‘dir’)mkdir(output_dir);end
subj_entries = dir(root_dir);subj_dirs = subj_entries([subj_entries.isdir] & …~ismember({subj_entries.name}, {‘.’,‘..’}));
if isempty(subj_dirs)warning(‘No subject folders found in: %s’, root_dir);returnend
fprintf(‘\n=== PET Preprocessing Pipeline ===\n’);fprintf(‘Root dir : %s\n’, root_dir);fprintf(‘Output dir: %s\n’, output_dir);fprintf(‘Subjects : %d\n\n’, numel(subj_dirs));
for s = 1:numel(subj_dirs)subj_name = subj_dirs(s).name;subj_path = fullfile(root_dir, subj_name);fprintf('[%d/%d] Subject: %s\n', s, numel(subj_dirs), subj_name);
% Find all .nii files recursively under this subject folder
nii_files = find_nii_files(subj_path);
if isempty(nii_files)
fprintf(' >> No .nii files found – skipping.\n');
continue
end
for f = 1:numel(nii_files)
src_nii = nii_files{f};
fprintf(' [%d/%d] %s\n', f, numel(nii_files), src_nii);
try
% ---- mirror folder structure in output_dir ---- %
rel_path = strrep(src_nii, root_dir, ''); % relative path
rel_path = strtrim(rel_path);
if rel_path(1) == filesep
rel_path(1) = [];
end
dst_nii = fullfile(output_dir, rel_path);
dst_folder = fileparts(dst_nii);
if ~exist(dst_folder, 'dir')
mkdir(dst_folder);
end
% Copy source .nii (and companion .mat if present) to output
copyfile(src_nii, dst_nii);
src_mat = strrep(src_nii, '.nii', '.mat');
if exist(src_mat, 'file')
copyfile(src_mat, strrep(dst_nii, '.nii', '.mat'));
end
% From here on, work on the copy so the source is untouched
pet_file = dst_nii;
reset_origin(pet_file);
[pet_dir, pet_stem, pet_ext] = fileparts(pet_file);
warped_file = fullfile(pet_dir, ['w' pet_stem pet_ext]);
run_normalise(pet_file, template);
[~, w_stem, w_ext] = fileparts(warped_file);
global_normed_file = fullfile(pet_dir, ['g' w_stem w_ext]);
run_intensity_norm(warped_file, global_normed_file);
[~, gn_stem, gn_ext] = fileparts(global_normed_file);
smoothed_file = fullfile(pet_dir, ['s' gn_stem gn_ext]);
spm_smooth(global_normed_file, smoothed_file, [8 8 8]);
delete_if_exists(warped_file);
delete_if_exists(global_normed_file);
sn_mat = fullfile(pet_dir, [pet_stem '_sn.mat']);
delete_if_exists(sn_mat);
% delete_if_exists(pet_file);
fprintf(' >> Done: %s\n', smoothed_file);
catch ME
fprintf(' !! ERROR on %s\n %s\n', src_nii, ME.message);
end
end
end
fprintf('\n=== All subjects processed. ===\n');
end
function nii_list = find_nii_files(folder)
% Return cell array of absolute paths to all .nii files under folder.
items = dir(fullfile(folder, '**', '*.nii'));
nii_list = fullfile({items.folder}, {items.name})';
end
function reset_origin(nii_file)
% Reset the origin of a NIfTI volume to its centre of mass.
V = spm_vol(nii_file);
dim = V.dim;
M = V.mat;
M(1:3, 4) = -M(1:3, 1:3) * (dim(:) / 2);
spm_get_space(nii_file, M);
end
function run_normalise(pet_file, template)
% Estimate + write spatial normalisation to MNI space using SPM OldNorm.
clear matlabbatch
matlabbatch{1}.spm.spatial.normalise.estwrite.subj.source = {pet_file};
matlabbatch{1}.spm.spatial.normalise.estwrite.subj.wtsrc = '';
matlabbatch{1}.spm.spatial.normalise.estwrite.subj.resample = {pet_file};
matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.template = {template};
matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.weight = '';
matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.smosrc = 8;
matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.smoref = 0;
matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.regtype = 'mni';
matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.cutoff = 25;
matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.nits = 16;
matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.reg = 1;
matlabbatch{1}.spm.spatial.normalise.estwrite.roptions.preserve = 0;
matlabbatch{1}.spm.spatial.normalise.estwrite.roptions.bb = ...
[-90 -126 -72; 90 90 108];
matlabbatch{1}.spm.spatial.normalise.estwrite.roptions.vox = [2 2 2];
matlabbatch{1}.spm.spatial.normalise.estwrite.roptions.interp = 4;
matlabbatch{1}.spm.spatial.normalise.estwrite.roptions.wrap = [0 0 0];
matlabbatch{1}.spm.spatial.normalise.estwrite.roptions.prefix = 'w';
spm_jobman('run', matlabbatch);
end
function run_intensity_norm(in_file, out_file)
V = spm_vol(in_file);
Y = spm_read_vols(V);
mask = Y > 0;
if ~any(mask(:))
error('Intensity normalisation: no positive voxels found in %s', in_file);
end
global_mean = mean(Y(mask));
fprintf(' Global mean = %.4f\n', global_mean);
Y_norm = Y / global_mean;
Vout = V;
Vout.fname = out_file;
spm_write_vol(Vout, Y_norm);
end
function delete_if_exists(fpath)
if exist(fpath, 'file')
delete(fpath);
end
end
root_dir = '...';
output_dir = '...';
preprocess_pet_batch(root_dir, output_dir);