Skip to content
Snippets Groups Projects
Commit 6027f14d authored by Vandewalle Gilles's avatar Vandewalle Gilles
Browse files

Upload New File

parent c4be0003
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
MRI_processing/Struct_Preprocessing_Synstrip.m 0 → 100644
+ 290
− 0
View file @ 6027f14d
%% Script for the processing of 7T structural images
%
% Cyclotron Research Centre, University of Liege, Belgium
%
% - Use of Synthstrip method (FreeSurfer) for brain extraction.
%
% Preprocessing of UNI images in different steps:
% 1- Remove background
% 2- Auto-reorient
% 3- Segment
% 4- Normalise
% 5- Brain extraction
%
% Arg:
% - overwrite: 0 or 1, checks if an ouptut already exists and re-run the
% preprocessing of structural data (1) or not (0).
%
% Output:
% - Created folder "ses-struct/output_preprocessing"
%
% --> BE CAREFUL : this code only reads "prepared" data. You need to have
% the data downloaded on your computer in the prepared format:
%
%..../prepared
% --> sub-xxx
% --> ses-evening
% --> MRI
% --> 005-b1map_sag_p2
% --> 007-cmrr_mbep2d_p3_mb2_1.4iso_VolCheck
% --> ...
% --> LM
% --> Physio
% --> ses-morning
% --> MRI
% --> 005-b1map_sag_p2
% --> 007-cmrr_mbep2d_p3_mb2_1.4iso_VolCheck
% --> ...
% --> LM
% --> Physio
% --> ses-struct
% --> MRI
% --> 002-B0_Shimming_4.0mm
% --> 003-B0_map
% --> ...
%
function Struct_Preprocessing_Synstrip(overwrite)
% Set environmental variables
% Paths
% load file 'your_paths.txt' this is a document containing:
% 1. path to the data
% 2. path to SPM
% 3. path to template
% 4. path to bash scripts
% 5. path to ANTs
% 6. path to edf toolbox
paths_txt = fileread('/home/islay/Documents/HILIGHT_CODES/MRI_processing/Structural_Prepoc/Prepared/your_paths.txt') ;
% split data by new lines
% the data is structured as follows:
% path to your data (example):
% corresponding path, for example: /mnt/data/Prepared
% ...
paths_txt = char(strsplit(paths_txt, '\n')) ;
% create a struct to store the paths we are interested in
paths = [];
paths.data = (strtrim(paths_txt(2,:)));
paths.spm = (strtrim(paths_txt(4,:)));
paths.templateANTs = (strtrim(paths_txt(6,:)));
paths.bash = (strtrim(paths_txt(8,:)));
paths.ANTs = (strtrim(paths_txt(10,:)));
paths.synthstrip = (strtrim(paths_txt(14,:)));
% add SPM path to matlab's search path
addpath(paths.spm)
% go to the path where the data is stored
cd(paths.data);
[subdir] = spm_select(Inf,'dir','Select subject directory');% select directories of subjects to analyses / directory has to be 'sub-0ID'
% get the i of subject storing the last 2 digits
IDs = str2num(subdir(:,end-2:end));
% add toolboxes
fprintf('\n path initialization done \n')
%% Start preprocessing
spm fmri % open progress window
for id = 1:length(IDs)
% change the format of the ID to 2 digits: For example: 30------>030
subject_ID = sprintf('%03d', IDs(id)) ;
% go to the data and extract information about:
% 1. name of the file which contains the data
% 2. location
% 3. Last modification date (timestamp)
% 4. Size of the file
% 5. Logical variable, 1 if it's folder, 0 if it's file
T = dir(fullfile(paths.data, 'HILIGHT_participants.xlsx')) ;
% store the data of the file
T = readtable(char([paths.data '/' T.name])) ;
% define the subject ID
entry_table = ['HILIGHT_', subject_ID ] ;
% create a logical array to identify the row location of the subject
index = contains(T.subject_Ids, entry_table) ;
% extract the information of the subjects we are interested in
patient_data = T(index,:) ;
% extract information of MP2RAGE acquisition and store it in a matrix
% 1st column: UNI session
% 2nd column: INV 1 session
% 3rd column: INV 2 session
mp2rage = [str2double(patient_data(1,:).MP2RAGE{:}) str2double(patient_data(1,:).Var19{:}) str2double(patient_data(1,:).Var20{:})] ;
% use isempty command to see if the MP2RAGE data for the subjects we
% are interested in exists
structural_data = isempty(mp2rage) ;
% MP2RAGE data doesn't exist, or the path to the subjects data we are
% interested in doesn't exist, it will display a message saying the
% data we are interested in wasn't found, we will jump into the next
% subject
if structural_data==1 || ~exist(fullfile(char([paths.data '/sub-' subject_ID])))
message = ['(1) No data found for sub-', subject_ID, 'or (2) No corresponding UNI, INV1 and/or INV2 sequences found in the Excel table'];
disp(message)
continue
end
%this part of the code search for the UNI, INV1 and INV2 files of the
%current subjects being processed
% Change the format of the subject enumeration style. For example: 09
% -----> 009
subject_ID = sprintf('%03d', IDs(id)) ;
message = ['processing structural images of subject: sub-', subject_ID] ;
disp(message)
fprintf('Looking for NIFTI structural files')
% list all the files in /mnt/data/Preprocessed and store it in a
% variable
structural_data = dir(char([paths.data])) ;
% we store the path of the subject's data we are interested in
File = fullfile(char([paths.data]),['sub-', subject_ID]);
% list all the files of the folder containing the sessions of the subject and store it in a
% variable
structural_data = dir(char([File '/ses-struct' '/MRI'])) ;
% convert the struct into cell
structural_data = struct2cell(structural_data).';
% Define paths of INV1, INV2 and UNI images
INV1_string = 'INV1' ;
INV2_string = 'INV2' ;
UNI_string = 'UNI_Images' ;
INV1 = fullfile(char([File '/ses-struct' '/MRI']), char(structural_data(endsWith(structural_data(:,1),INV1_string))));
INV2 = fullfile(char([File '/ses-struct' '/MRI']), char(structural_data(endsWith(structural_data(:,1),INV2_string))));
UNI = fullfile(char([File '/ses-struct' '/MRI']), char(structural_data(endsWith(structural_data(:,1),UNI_string)))) ;
% Define paths of UNI nifti image
UNI_dir = char(UNI) ;
UNI = dir(UNI) ;
% convert the struct into cell
UNI = struct2cell(UNI).';
nifti = '.nii' ;
UNI = fullfile(UNI_dir, char(UNI(endsWith(UNI(:,1),nifti))));
fprintf('\n UNI file found \n ')
% Define paths of INV1 nifti image
INV1_dir = char(INV1) ;
INV1 = dir(INV1) ;
% convert the struct into cell
INV1 = struct2cell(INV1).';
INV1 = fullfile(INV1_dir, char(INV1(endsWith(INV1(:,1),nifti))));
fprintf('\n INV1 file found \n')
% Define paths of INV2 nifti image
INV2_dir = char(INV2) ;
INV2 = dir(INV2) ;
% convert the struct into cell
INV2 = struct2cell(INV2).';
INV2 = fullfile(INV2_dir, char(INV2(endsWith(INV2(:,1),nifti))));
fprintf('\n INV2 file found \n')
%If not created yet, a new folder "ouput_preprocessing" will be created
%in the /.../prepared/sub-00ID/ses-/MRI/ path and contains a folder
%"mp2rage_UNI" were the preprocessed structural UNI images will be
%saved
data_dir = char([File '/ses-struct' '/MRI' '/output_preprocessing']) ;%we save that path to create an output folder
output_preprocessed_UNI = char([data_dir '/mp2rage_UNI']) ;
%verifying if an output already exist. If no, preprocessing will be
%conducted. If yes : no further processing will be conducted unless
%"overwrite" is set to true.
if ~isempty(dir(output_preprocessed_UNI)) && overwrite == false
message = ['Output already exists and "overwrite" set to false. No preprocessing conducted for sub-', subject_ID] ;
disp(message)
continue
end
% Create folders to store output files in case they don't exist
if ~exist(data_dir)
mkdir(data_dir)
mkdir(output_preprocessed_UNI)
elseif ~exist(output_preprocessed_UNI)
mkdir(output_preprocessed_UNI)
end
%set up batch
clear matlabbatch
spm_jobman('initcfg') ;
matlabbatch{1}.spm.tools.mp2rage.rmbg.INV1 = {INV1};
matlabbatch{1}.spm.tools.mp2rage.rmbg.INV2 = {INV2};
matlabbatch{1}.spm.tools.mp2rage.rmbg.UNI = {UNI};
matlabbatch{1}.spm.tools.mp2rage.rmbg.regularization = 5;
matlabbatch{1}.spm.tools.mp2rage.rmbg.output.dirfile.dirpath = {output_preprocessed_UNI};
matlabbatch{1}.spm.tools.mp2rage.rmbg.output.dirfile.filename = 'clean_UNI';
matlabbatch{1}.spm.tools.mp2rage.rmbg.show = 'yes';
matlabbatch{2}.spm.spatial.autoreor.image(1) = cfg_dep('Remove background: Background free image', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','files'));
matlabbatch{2}.spm.spatial.autoreor.imgtype = 't1canonical';
matlabbatch{2}.spm.spatial.autoreor.other = {''};
matlabbatch{3}.spm.spatial.preproc.channel.vols(1) = cfg_dep('Auto-Reorient: Auto-reoriented Image(s)', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','afiles'));
matlabbatch{3}.spm.spatial.preproc.channel.biasreg = 0.001;
matlabbatch{3}.spm.spatial.preproc.channel.biasfwhm = 30;
matlabbatch{3}.spm.spatial.preproc.channel.write = [0 1];
matlabbatch{3}.spm.spatial.preproc.tissue(1).tpm = cellstr([paths.spm '/tpm/eTPM.nii,1']);
matlabbatch{3}.spm.spatial.preproc.tissue(1).ngaus = 2;
matlabbatch{3}.spm.spatial.preproc.tissue(1).native = [1 0];
matlabbatch{3}.spm.spatial.preproc.tissue(1).warped = [1 1];
matlabbatch{3}.spm.spatial.preproc.tissue(2).tpm = cellstr([paths.spm '/tpm/eTPM.nii,2']);
matlabbatch{3}.spm.spatial.preproc.tissue(2).ngaus = 2;
matlabbatch{3}.spm.spatial.preproc.tissue(2).native = [1 0];
matlabbatch{3}.spm.spatial.preproc.tissue(2).warped = [1 1];
matlabbatch{3}.spm.spatial.preproc.tissue(3).tpm = cellstr([paths.spm '/tpm/eTPM.nii,3']);
matlabbatch{3}.spm.spatial.preproc.tissue(3).ngaus = 2;
matlabbatch{3}.spm.spatial.preproc.tissue(3).native = [1 0];
matlabbatch{3}.spm.spatial.preproc.tissue(3).warped = [1 1];
matlabbatch{3}.spm.spatial.preproc.tissue(4).tpm = cellstr([paths.spm '/tpm/eTPM.nii,4']);
matlabbatch{3}.spm.spatial.preproc.tissue(4).ngaus = 3;
matlabbatch{3}.spm.spatial.preproc.tissue(4).native = [0 0];
matlabbatch{3}.spm.spatial.preproc.tissue(4).warped = [0 0];
matlabbatch{3}.spm.spatial.preproc.tissue(5).tpm = cellstr([paths.spm '/tpm/eTPM.nii,5']);
matlabbatch{3}.spm.spatial.preproc.tissue(5).ngaus = 4;
matlabbatch{3}.spm.spatial.preproc.tissue(5).native = [0 0];
matlabbatch{3}.spm.spatial.preproc.tissue(5).warped = [0 0];
matlabbatch{3}.spm.spatial.preproc.tissue(6).tpm = cellstr([paths.spm '/tpm/eTPM.nii,6']);
matlabbatch{3}.spm.spatial.preproc.tissue(6).ngaus = 2;
matlabbatch{3}.spm.spatial.preproc.tissue(6).native = [0 0];
matlabbatch{3}.spm.spatial.preproc.tissue(6).warped = [0 0];
matlabbatch{3}.spm.spatial.preproc.warp.mrf = 1;
matlabbatch{3}.spm.spatial.preproc.warp.cleanup = 0;
matlabbatch{3}.spm.spatial.preproc.warp.reg = [0 0.001 0.5 0.05 0.2];
matlabbatch{3}.spm.spatial.preproc.warp.affreg = 'mni';
matlabbatch{3}.spm.spatial.preproc.warp.fwhm = 0;
matlabbatch{3}.spm.spatial.preproc.warp.samp = 2.5;
matlabbatch{3}.spm.spatial.preproc.warp.write = [1 1];
matlabbatch{3}.spm.spatial.preproc.warp.vox = NaN;
matlabbatch{3}.spm.spatial.preproc.warp.bb = [NaN NaN NaN
NaN NaN NaN];
matlabbatch{4}.spm.spatial.normalise.write.subj.def(1) = cfg_dep('Segment: Forward Deformations', substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','fordef', '()',{':'}));
matlabbatch{4}.spm.spatial.normalise.write.subj.resample(1) = cfg_dep('Segment: Bias Corrected (1)', substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','channel', '()',{1}, '.','biascorr', '()',{':'}));
matlabbatch{4}.spm.spatial.normalise.write.woptions.bb = [-78 -112 -70
78 76 85];
matlabbatch{4}.spm.spatial.normalise.write.woptions.vox = [1 1 1];
matlabbatch{4}.spm.spatial.normalise.write.woptions.interp = 4;
matlabbatch{4}.spm.spatial.normalise.write.woptions.prefix = 'w';
spm_jobman('run', matlabbatch) % run batch
fprintf('________________________________________________________')
message = ['taking mclean for brain extraction for subject: sub-', subject_ID] ;
disp(message)
% Brain extraction using Synthstrip
input1= [output_preprocessed_UNI '/mclean_UNI.nii']; % input image
input2= [output_preprocessed_UNI '/BET_STR.nii']; % output image name
cmdStr=[paths.synthstrip '/synthstrip-docker -i ' input1 ' -o ' input2];
system(cmdStr)
end
fprintf('\n done : all subjects have been processed \n')
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment