Preproc_finalver1.m

%% Preproc final
% Make sure the functional files do not start with m
clear;

addpath('/mnt/d/toolboxes/spm12/spm12/')
path = uigetdir('/mnt/d/codes/preprocessing related/');
cd (path);
Subj_list=dir('*0*');
nt = 404;
other_name_func = 1;
other_name_anat = 1;

tic
%%  Realignment
sub_not_processed = zeros(1,length(Subj_list));
for subji =1:length(Subj_list)                                              % number of subjects
    cd(fullfile(Subj_list(subji).folder,Subj_list(subji).name))             % Change directory to the subject folder
    Subj_dir = dir();                                                       % Subject folder
    Subj_dir([1,2]) = [];                                                   % Delete first two as they are not files
    if exist('func','dir')                                                  % if the func dir exits then make that as the functional folder
        func_fl_name = 'func';
        path_func = fullfile(Subj_dir(1).folder,'func');
    end

    if exist('anat','dir')                                                  % if the anat dir exits then make that as the anatomical folder
        anat_fl_name = 'anat';
        path_anat = fullfile(Subj_dir(1).folder,'anat');
    end

    if other_name_func && ~exist('func_fl_name','var')                         % if the path_func variable does not exits that means the name of functional folder is different
        % Specify the functional folder
        path_func = string(uigetdir(fullfile(Subj_list(subji).folder,Subj_list(subji).name),'Specify the Functional scan folder')); % Convert Char to str so we can use erase
        temp = erase(path_func,string(fullfile(Subj_list(subji).folder,Subj_list(subji).name)));             %Get the path of the func folder from the subject directory  %strsplit(path_func,'/');
        method=questdlg('Do the other subjects also have the same name', 'Same Name?','Yes','No','Yes');
        switch method
            case 'Yes'
                func_fl_name = string(temp);
                other_name_func = 0;
            case 'No'
                other_name_func = 1;
                clear func_fl_name
        end
    end    
    if other_name_anat && ~exist('anat_fl_name','var')                         % if the path_anat variable does not exits that means the name of anatomical folder is different
        % Specify the anatomical folder
        path_anat = string(uigetdir(fullfile(Subj_list(subji).folder,Subj_list(subji).name),'Specify the Anatomical scan folder'));
        temp = erase(path_anat,string(fullfile(Subj_list(subji).folder,Subj_list(subji).name))); %strsplit(path_anat,'/');
        method=questdlg('Do the other subjects also have the same name', 'Same Name?','Yes','No','Yes');
        switch method
            case 'Yes'
                anat_fl_name = string(temp);
                other_name_anat = 0;
            case 'No'
                other_name_anat = 1;
                clear anat_fl_name
        end
    end
    func_dir = dir(fullfile(Subj_list(subji).folder,Subj_list(subji).name,func_fl_name));                                              % Get the functional image folder
    func_dir([1,2]) = [];
    if exist('rest_name','var') 
        if strcmp(rest_name(end-2:end),'.gz')                          % if .gz file unzip it first
            rest_image = dir(fullfile(func_dir(1).folder,rest_name));
            gunzip(fullfile(func_dir(1).folder,rest_image.name)) 
            rest_image = dir(fullfile(func_dir(1).folder,rest_name(1:end-3)));
            temp = niftiread(fullfile(func_dir(1).folder,rest_image.name));
            nt = size(temp,4);
        elseif strcmp(rest_name(end-3:end),'.nii')
            rest_image = dir(fullfile(func_dir(1).folder,rest_name));
            temp = niftiread(fullfile(func_dir(1).folder,rest_image.name));
            nt = size(temp,4);
        end
        if isempty(rest_image)                                       % if the mentioned file does not exist
            sub_not_processed(subji) = 1;
            continue
        end
    else
        rest_image = dir(fullfile(func_dir(1).folder,'*.nii'));            
        if isempty(rest_image)
            rest_image = dir(fullfile(func_dir(1).folder,'*.gz'));
            if (length(rest_image) > 1)                                                                                             % If there are more than 1 gz files then
                rest_image_name = uigetfile({'*.gz'},'More than one gz files in the Functional folder Please choose 1',path_func);  % ask for which gz file to choose
                gunzip(fullfile(func_dir(1).folder,rest_image_name))                                                                   % Unzip the gz file
                rest_image = dir(fullfile(func_dir(1).folder,rest_image_name(1:end-3)));                                               % Get the .nii file
                temp = niftiread(fullfile(func_dir(1).folder,rest_image.name));
                nt = size(temp,4);
                method=questdlg('Do the other subjects also have the similar name', 'Similar Name?','Yes','No','Yes');
                switch method
                case 'Yes'
                    rest_name = char(inputdlg(['Please type the letters similar in other subjects. You selected ',rest_image_name]));
                    other_name_rest = 0;
                case 'No'
                    other_name_rest = 1;
                    clear rest_name
                end
            end
        elseif (length(rest_image) > 1)                                                                                           % If there are more than 1 nifti files then
        clear rest_image
        rest_image_name = uigetfile({'*.nii'},'More than one nifti files in the Functional folder Please choose 1',path_func);
        rest_image = dir(fullfile(func_dir(1).folder,rest_image_name));                                               % Get the .nii file
                method=questdlg('Do the other subjects also have the similar name', 'Similar Name?','Yes','No','Yes');
                switch method
                case 'Yes'
                    rest_name = inputdlg(['Please type the letters similar in other subjects. You selected ',rest_image_name]);
                    other_name_rest = 0;
                case 'No'
                    other_name_rest = 1;
                    clear rest_name
                end
        end
        
    end

    for ii = 1:nt                                                           % number of timepoints
        pathlist{ii,1} = [char(fullfile(Subj_list(subji).folder,Subj_list(subji).name,func_fl_name,rest_image.name)),',',num2str(ii)];
    end

    matlabbatch{1}.spm.spatial.realign.estwrite.data = {pathlist}';

    matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.quality = 0.9;           % Quality 1--> max quality
    matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.sep = 4;                 % separation (in mm) between the points sampled in the reference image. Smaller speration gives more accurate results, but will be slower. 
    matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.fwhm = 5;                % FWHM Gaussian smothing kernel parameter in mm
    matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.rtm = 1;                 % Register to mean
    matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.interp = 2;              % The method by which the images are sampled when estimating the optimum transformation. 2nd degree bspline
    matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.wrap = [0 0 0];          % Directions in the volumes the values should wrap around in.
    matlabbatch{1}.spm.spatial.realign.estwrite.eoptions.weight = '';             % Optional weighting image to weight each voxel of the reference image differently when estimating the realignment parameters.
    matlabbatch{1}.spm.spatial.realign.estwrite.roptions.which = [2 1];           % 
    matlabbatch{1}.spm.spatial.realign.estwrite.roptions.interp = 4;              % 4th degree bspline
    matlabbatch{1}.spm.spatial.realign.estwrite.roptions.wrap = [0 0 0];          % no wraping
    matlabbatch{1}.spm.spatial.realign.estwrite.roptions.mask = 1;                % Because  of  subject  motion,  different images are likely to have different patterns of zeros from where it was not possible to sample data. With masking enabled, the program searches through the whole time series looking for voxels which need to be sampled from outside the original images. Where this occurs, that voxel is set to zero for the whole set of images (unless the image format can represent NaN, in which case NaNs are used where possible).
    matlabbatch{1}.spm.spatial.realign.estwrite.roptions.prefix = 'r';            % prefix

    spm('defaults', 'FMRI');
    spm_jobman('run', matlabbatch);
%     cfg_util('run',matlabbatch);
    clear matlabbatch;
    clearvars -except path Subj_list nt other_name_func other_name_anat func_fl_name anat_fl_name path_func path_anat sc_st rest_name sub_not_processed
     
end
	
%% Framewise Displacement 
% check the framewize displacement and exclude any if the criteria is not
% matched

excluded = zeros(1,length(Subj_list));
p = 1;
for subji = 1:length(Subj_list)                                % For loop on number od subjects
    % Get rest folder
    rest_folder = fullfile(Subj_list(subji).folder,Subj_list(subji).name,func_fl_name);            % folder path
    text_file = dir(fullfile(rest_folder,'*.txt'));
    rp_rest = load(fullfile(text_file.folder,text_file.name));              % input the text file generated at the Realignment stage
   
    rp_diff_trans = diff(rp_rest(:,1:3));                      % The first 3 parameters tell the displacement in x,y, and z direction in mm. Here the vector difference operator is used to get the derivative of vector. (framewise difference)
    rp_diff_rotat = diff(rp_rest(:,4:6)*180/pi);               % The last  3 parameters tell the rotation values pitch, yaw and roll in radians (here we convert them to degrees)
   
    multp = rp_diff_trans*rp_diff_trans';                      % diagonals will represent dx^2 + dy^2 + dz^2 
    fd_trans = sqrt(diag(multp));                              % sqrt (dx^2 + dy^2 + dz^2)
   
    multp = rp_diff_rotat*rp_diff_rotat';                      % diagonals will represent dpitch^2 + dyaw^2 + droll^2
    fd_rotat = sqrt(diag(multp));                              % sqrt (dpitch^2 + dyaw^2 + droll^2)
   
    fd_max_trans(subji,1) = max(fd_trans);                     % Get the maximum traslational framewize displacement
    fd_max_rotat(subji,1) = max(fd_rotat);                     % Get the maximum rotational framewize displacement
   
    fd_mean_trans(subji,1) = mean(fd_trans);                   % Get the mean traslational framewize displacement
    fd_mean_rotat(subji,1) = mean(fd_rotat);                   % Get the mean rotational framewize displacement
    
    if fd_max_trans(subji)>1.5 || fd_max_rotat(subji)>1.5 || fd_mean_trans(subji)>0.2 || fd_mean_rotat(subji)>0.2
        excluded(subji) = 1;
        excluded_subjects(p) = {Subj_list(subji).name};
        p = p+1;
    end
end

excluded = find(fd_max_trans>1.5|fd_max_rotat>1.5|fd_mean_trans>0.2|fd_mean_rotat>0.2);
% Try less stringent threshold for children data:
% excluded = find(fd_max_trans>2|fd_max_rotat>2|fd_mean_trans>0.3|fd_mean_rotat>0.3);
excludedID = Subj_list(excluded);


if nnz(excluded)
    ask = questdlg(['A total of ',num2str(nnz(excluded)),' Subjects were discarded out of ',num2str(length(Subj_list)),' subjects. Do you want to Proceed'],'Proceed?','Yes','No','No');

    switch ask
        case 'No'
            return
    end   
end
%%
skull_strip=questdlg('Do you want to Skull Strip','Skull Strip', ...
                                                    'Yes AFNI','Yes FSL','No','No');

switch skull_strip 
    case 'Yes AFNI'
        sc_st = 2;         % 2: AFNI
    case 'Yes FSL'
        sc_st = 1;         % 1: FSL
    case 'No'
        sc_st = 0;         % 0: No Skull Stripping
end

%% Preprocess the non-excluded subjects
subji_ex = 1;
T = zeros(length(Subj_list) - nnz(excluded),1);
for subji = 1:length(Subj_list)                                                       % For loop on number of subjects

    if excluded(subji) == 0
    
        anat_folder = fullfile(Subj_list(subji).folder,...                                % Get the anatomical image folder
                        Subj_list(subji).name,anat_fl_name);
%% Scull stripping
% You may perform skull stripping at this stage for better coregisteration
% output. For skull stripping you may either use AFNI or FSL
        if sc_st
            cd(anat_folder)                                                                   % Change Directory to the anat image folder
            anat_image = dir('*.nii');                                                        % get the anatomical image
            if ~(length(anat_image) == 1)                                                     % If there are more than 1 nifti files
                clear anat_image
                anat_image_name = uigetfile({'*.nii'},...                                     % Ask the user to choose 1
                    'More than one nifti files in the anatomical folder Please choose 1',...
                    anat_folder);
            end
            if sc_st == 1
                system(['bet2 ',anat_image.name, [' b',anat_image.name]])                         % FSL command to skull strip the image
            elseif sc_st == 2
                system(['3dSkullStrip -input ',anat_image.name, ' -prefix' [' b',anat_image.name]]) % AFNI command to skull strip the image
            end
        end


%% COREGISTRATION REST
% Coregister anatomical image to the functional image

        func_folder = fullfile(Subj_list(subji).folder,...                                % Get the functional image folder
        Subj_list(subji).name,func_fl_name);
        cd(anat_folder)                                                                   % Change Directory to the anat image folder
        if sc_st == 1
            anat_image = dir('b*.nii.gz');     
            gunzip(fullfile(anat_image.folder,anat_image.name))
            anat_image = dir('b*.nii');     
            matlabbatch{1}.spm.spatial.coreg.estimate.ref = {[fullfile(anat_image.folder,anat_image.name),',1']};           % Reference image (this does not change) (Here its the anatomical image)
        elseif sc_st == 2
            anat_image = dir('b*.nii'); 
            matlabbatch{1}.spm.spatial.coreg.estimate.ref = {[fullfile(anat_image.folder,anat_image.name),',1']};       % Reference image (this does not change) (Here its the anatomical image)

        else
            anat_image = dir('*.nii'); 
            if length(anat_image) == 1
                matlabbatch{1}.spm.spatial.coreg.estimate.ref = {[fullfile(anat_image.folder,anat_image.name),',1']};       % Reference image (this does not change) (Here its the anatomical image)
            else 
                anat_image_name = uigetfile({'*.nii'},...                                     % Ask the user to choose 1
                    'More than one nifti files in the anatomical folder Please choose 1',...
                    anat_folder);
                anat_image = dir(anat_image_name);
                matlabbatch{1}.spm.spatial.coreg.estimate.ref = {[fullfile(anat_image.folder,anat_image.name),',1']};
            end
        end
        cd(func_folder)
        mean_func = dir('m*.nii');
        matlabbatch{1}.spm.spatial.coreg.estimate.source = {[fullfile(mean_func.folder,mean_func.name),',1']}; % Source image (This will change) (here its the functional image)

        func = dir('r*.nii');
        for ii= 1:nt    % Number of timepoints
            pathlist{ii,1} = strcat([fullfile(func.folder,func.name),',',num2str(ii)]);
        end

        matlabbatch{1}.spm.spatial.coreg.estimate.other = (pathlist);         % These are any images that need to remain in alignment with the moved image

        matlabbatch{1}.spm.spatial.coreg.estimate.eoptions.cost_fun = 'nmi';  % Cost function is normalized mutual information
        matlabbatch{1}.spm.spatial.coreg.estimate.eoptions.sep = [4 2];       % The average distance between sampled points (in mm).
        matlabbatch{1}.spm.spatial.coreg.estimate.eoptions.tol = [0.02 0.02 0.02 0.001 0.001 0.001 0.01 0.01 0.01 0.001 0.001 0.001]; % Iterations  stop  when differences between successive estimates are less than the required tolerance.
        matlabbatch{1}.spm.spatial.coreg.estimate.eoptions.fwhm = [5 5];      % Histogram smoothing by Gaussian smoothing to apply to the 256x256 joint histogram. 
        
%         cfg_util('run',matlabbatch);
        spm('defaults', 'FMRI');
        spm_jobman('run', matlabbatch);
        clear matlabbatch;
        
%% Segmentation 
        matlabbatch{1}.spm.spatial.preproc.channel.vols = {fullfile(anat_image.folder,anat_image.name)};  % Select Volumes for processing
        matlabbatch{1}.spm.spatial.preproc.channel.biasreg = 0.001;
        matlabbatch{1}.spm.spatial.preproc.channel.biasfwhm = 60;
        matlabbatch{1}.spm.spatial.preproc.channel.write = [0 1];
        matlabbatch{1}.spm.spatial.preproc.tissue(1).tpm = {'/mnt/d/iit_mandi/Research/spm12/spm12/tpm/TPM.nii,1'};
        matlabbatch{1}.spm.spatial.preproc.tissue(1).ngaus = 1;
        matlabbatch{1}.spm.spatial.preproc.tissue(1).native = [1 0];
        matlabbatch{1}.spm.spatial.preproc.tissue(1).warped = [0 0];
        matlabbatch{1}.spm.spatial.preproc.tissue(2).tpm = {'/mnt/d/iit_mandi/Research/spm12/spm12/tpm/TPM.nii,2'};
        matlabbatch{1}.spm.spatial.preproc.tissue(2).ngaus = 1;
        matlabbatch{1}.spm.spatial.preproc.tissue(2).native = [1 0];
        matlabbatch{1}.spm.spatial.preproc.tissue(2).warped = [0 0];
        matlabbatch{1}.spm.spatial.preproc.tissue(3).tpm = {'/mnt/d/iit_mandi/Research/spm12/spm12/tpm/TPM.nii,3'};
        matlabbatch{1}.spm.spatial.preproc.tissue(3).ngaus = 2;
        matlabbatch{1}.spm.spatial.preproc.tissue(3).native = [1 0];
        matlabbatch{1}.spm.spatial.preproc.tissue(3).warped = [0 0];
        matlabbatch{1}.spm.spatial.preproc.tissue(4).tpm = {'/mnt/d/iit_mandi/Research/spm12/spm12/tpm/TPM.nii,4'};
        matlabbatch{1}.spm.spatial.preproc.tissue(4).ngaus = 3;
        matlabbatch{1}.spm.spatial.preproc.tissue(4).native = [0 0];
        matlabbatch{1}.spm.spatial.preproc.tissue(4).warped = [0 0];
        matlabbatch{1}.spm.spatial.preproc.tissue(5).tpm = {'/mnt/d/iit_mandi/Research/spm12/spm12/tpm/TPM.nii,5'};
        matlabbatch{1}.spm.spatial.preproc.tissue(5).ngaus = 4;
        matlabbatch{1}.spm.spatial.preproc.tissue(5).native = [0 0];
        matlabbatch{1}.spm.spatial.preproc.tissue(5).warped = [0 0];
        matlabbatch{1}.spm.spatial.preproc.tissue(6).tpm = {'/mnt/d/iit_mandi/Research/spm12/spm12/tpm/TPM.nii,6'};
        matlabbatch{1}.spm.spatial.preproc.tissue(6).ngaus = 2;
        matlabbatch{1}.spm.spatial.preproc.tissue(6).native = [0 0];
        matlabbatch{1}.spm.spatial.preproc.tissue(6).warped = [0 0];
        matlabbatch{1}.spm.spatial.preproc.warp.mrf = 1;
        matlabbatch{1}.spm.spatial.preproc.warp.cleanup = 1;
        matlabbatch{1}.spm.spatial.preproc.warp.reg = [0 0.001 0.5 0.05 0.2];
        matlabbatch{1}.spm.spatial.preproc.warp.affreg = 'mni';
        matlabbatch{1}.spm.spatial.preproc.warp.fwhm = 0;
        matlabbatch{1}.spm.spatial.preproc.warp.samp = 3;
        matlabbatch{1}.spm.spatial.preproc.warp.write = [0 1];
%         cfg_util('run',matlabbatch);
        spm('defaults', 'FMRI');
        spm_jobman('run', matlabbatch);
        clear matlabbatch;

%% Making CSF_MASK for REST

        func = dir('r*.nii');
        cd(anat_folder) 
        c3 = dir('c3*.nii');
        matlabbatch{1}.spm.util.imcalc.input = {
            [fullfile(func.folder,func.name),',1']%[path Subj_list(subji).name '/wrREST1.nii,1']
             fullfile(c3.folder,c3.name) %[path Subj_list(subji).name '/wc3MPRAGE.nii']
            };
        matlabbatch{1}.spm.util.imcalc.output = 'CSF_MASK99.nii';
        matlabbatch{1}.spm.util.imcalc.outdir = {anat_folder};%{[path Subj_list(subji).name '/']};
        matlabbatch{1}.spm.util.imcalc.expression = 'i2>0.99';
        matlabbatch{1}.spm.util.imcalc.var = struct('name', {}, 'value', {});
        matlabbatch{1}.spm.util.imcalc.options.dmtx = 0;
        matlabbatch{1}.spm.util.imcalc.options.mask = 0;
        matlabbatch{1}.spm.util.imcalc.options.interp = 1;
        matlabbatch{1}.spm.util.imcalc.options.dtype = 4;

%         cfg_util('run',matlabbatch); 
        spm('defaults', 'FMRI');
        spm_jobman('run', matlabbatch);
        clear matlabbatch;


%% Making WM_MASK for REST

        func_folder = fullfile(Subj_list(subji).folder,...                                % Get the functional image folder
                              Subj_list(subji).name,func_fl_name);
        cd(func_folder)
        func = dir('r*.nii');

        anat_folder = fullfile(Subj_list(subji).folder,...                                % Get the anatomical image folder
                              Subj_list(subji).name,anat_fl_name);
        cd(anat_folder) 
        c2 = dir('c2*.nii');
        matlabbatch{1}.spm.util.imcalc.input = {
        [fullfile(func.folder,func.name),',1']%[path Subj_list(subji).name '/wrREST1.nii,1']
         fullfile(c2.folder,c2.name) %[path Subj_list(subji).name '/wc3MPRAGE.nii']
        };
        matlabbatch{1}.spm.util.imcalc.output = 'WM_MASK99.nii';
        matlabbatch{1}.spm.util.imcalc.outdir = {anat_folder};%{[path Subj_list(subji).name '/']};
        matlabbatch{1}.spm.util.imcalc.expression = 'i2>0.99';
        matlabbatch{1}.spm.util.imcalc.var = struct('name', {}, 'value', {});
        matlabbatch{1}.spm.util.imcalc.options.dmtx = 0;
        matlabbatch{1}.spm.util.imcalc.options.mask = 0;
        matlabbatch{1}.spm.util.imcalc.options.interp = 1;
        matlabbatch{1}.spm.util.imcalc.options.dtype = 4;

%         cfg_util('run',matlabbatch);  
        spm('defaults', 'FMRI');
        spm_jobman('run', matlabbatch);
        clear matlabbatch;

%% Extracting MEAN CSF and WM time-series

        cd(anat_folder) 
        % loading CSF MASKS
        CSF_MASK = spm_vol('CSF_MASK99.nii');%spm_vol([path Subj_list(subji).name '/CSF_MASK98.nii']);
        CSF_files  = spm_read_vols(CSF_MASK);
        image_dim = CSF_MASK.dim;
        CSF_Files_RS  = reshape(CSF_files,image_dim(1)*image_dim(2)*image_dim(3),1);
        clear image_dim

        % loading WM MASKS
        WM_MASK = spm_vol('WM_MASK99.nii');%spm_vol([path Subj_list(subji).name '/WM_MASK98.nii']);
        WM_files  = spm_read_vols(WM_MASK);
        image_dim = WM_MASK.dim;
        WM_Files_RS  = reshape(WM_files,image_dim(1)*image_dim(2)*image_dim(3),1);

        clear WM_REST CSF_REST REST_list REST_files image_dim2 REST_RS
        %%

    %     REST_folder = 'D:\NJIT\Research\Preprocessing\Practice\with_normalization\sub-10225\func';
        for imagei = 1:nt

            REST_list = spm_vol([fullfile(func.folder,func.name),',',num2str(imagei)]); %spm_vol([path Subj_list(subji).name '/wrREST1.nii,',num2str(imagei)]);
            REST_files = spm_read_vols(REST_list);
            image_dim2 = REST_list.dim;
            REST_RS = reshape(REST_files,image_dim2(1)*image_dim2(2)*image_dim2(3),1);

            WM_REST(:,imagei)  = REST_RS(find(WM_Files_RS > 0.5));
            CSF_REST(:,imagei) = REST_RS(find(CSF_Files_RS > 0.5));

        end
  %%  
        clear MEAN_WM_REST MEAN_CSF_REST VARname1 z*

        n1=isnan(CSF_REST);
        z1=find(n1>0);
        CSF_REST(z1)=0;

        n2=isnan(WM_REST);
        z2=find(n2>0);
        WM_REST(z2)=0;

        MEAN_WM  = mean(WM_REST);
        MEAN_CSF = mean(CSF_REST);


        MEAN_WM_REST  = MEAN_WM';
        MEAN_CSF_REST = MEAN_CSF';

        VARname1 = fullfile(func.folder,'covariance_wm_csf_REST.mat') ;%[path Subj_list(subji).name '/covariance_wm_csf_REST.mat'];
        save (VARname1, 'MEAN_WM_REST','MEAN_CSF_REST');
    

%% TEMPORAL REGRESSION: removing nusisance signals and motion correction

        cd(func_folder)
%         func = dir('r*.nii');
        % Loading One Image from MASKED_wrREST
        clear v_pair_* y_pair_* image_dim_* v_image_* y_image_* rp_* pca_CSF_* pca_WM_* mean_image_*
        v_pair_REST = spm_vol([fullfile(func.folder,func.name),',15']);%spm_vol([REST_folder 'wrREST1.nii,15']);
        y_pair_REST = spm_read_vols(v_pair_REST);
        image_dim_REST = v_pair_REST.dim;
        
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		%%% *** Please Note: The REST_mask.nii is a mask to perform regression only within the brain ...
		... This script assumes that such a mask has already been created. One can use AFNI's 3dAutomask or ...
		... a modified version of FSL's bet module to do so on the functional data or any other suitable method that works ...
		... If the mask is not present in each subject's folder this part of the script will fail of course ...
		... Modify the script accordingly if you do not want to use any mask, which will obviously take much longer ....
		
        system(['3dAutomask -prefix rest_mask.nii ',  func.name])
        % Loading REST_MASK
        MASK_list=spm_vol(fullfile(func.folder,'rest_mask.nii,1'));
        REST_MASK=spm_read_vols(MASK_list);
		
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
        % Loading Time Series of MASKED_wrREST
        clear imagei
        for imagei = 1:nt
            v_image_REST = spm_vol(fullfile(func.folder,[func.name,',',num2str(imagei)]));%spm_vol([REST_folder 'wrREST1.nii,' num2str(imagei)]);
            y_image_REST(:,:,:,imagei) = spm_read_vols(v_image_REST);
        end
        
        fprintf('loading covariates for REST... \n')
        
        fprintf('First loading the motion parameters for REST... \n')
        txt_file = dir(fullfile(func.folder,'*.txt'));
        rp_REST1=load(fullfile(txt_file.folder,txt_file.name));
        rp_REST_temp=rp_REST1(1:nt,:);
        rp_REST = zscore(rp_REST_temp);
        rp_previous_REST = [0 0 0 0 0 0; rp_REST(1:end-1,:)];
        rp_auto_REST = [rp_REST rp_REST.^2 rp_previous_REST rp_previous_REST.^2];
        
        fprintf('Now loading csf and wm pca components for REST... \n')
        
        fprintf('Loading covariates from: \n')
        fprintf(fullfile(func.folder,'covariance_wm_csf_REST.mat'))
        cov_file=load(fullfile(func.folder,'covariance_wm_csf_REST.mat'));
        MEAN_WM_REST=cov_file.MEAN_WM_REST;
        MEAN_CSF_REST=cov_file.MEAN_CSF_REST;
        
        fprintf('Loading of covariates complete! \n')
        
        fprintf('regressing out covariates for REST... \n')
        mean_image_REST = mean(y_image_REST,4);   % calculate mean across time for all voxels | can be added back after regression to improve ICA performance
        
        %%
        clear vi vj vk
        for vi = 1:image_dim_REST(1)
            fprintf('%d',vi)
            for vj = 1:image_dim_REST(2)
                for vk = 1:image_dim_REST(3)
                    if REST_MASK(vi,vj,vk)==0
                        
                        y_image_REST_regressed(vi,vj,vk,:) = zeros(1,1,1,nt);
                        
                    else
                        
                        % CSF, WM and MOTION Regressors and Regressors from OTHER SESSIONs for the same subject
                        
                        b = zscore([MEAN_WM_REST MEAN_CSF_REST rp_auto_REST]);
                        
                        [b,bint,y_image_REST_regressed(vi,vj,vk,:)] = regress(shiftdim(y_image_REST(vi,vj,vk,:),3),[b ones(nt,1)]);
                    end
                end
                
            end
        end
        
        fprintf('REST done! \n')
        fprintf('Release the Oompa Loompas!!\n')
        
        fprintf('saving REST image... \n')
        
        clear V
        V.dim   = v_pair_REST.dim;
        V.dt    = v_pair_REST.dt;
        V.mat   = v_pair_REST.mat;
        V.pinfo = v_pair_REST.pinfo;
        
        clear imagei
        for imagei = 1:nt
            V.n=[imagei 1];
            V.fname = fullfile(func.folder,['REG_',func.name]);
            V = spm_write_vol(V,y_image_REST_regressed(:,:,:,imagei)+mean_image_REST);
        end
        
%% Nomralization

        matlabbatch{1}.spm.spatial.normalise.write.subj.def = {fullfile(anat_folder,['y_',anat_image.name])};%{'/mnt/d/NJIT/Research/Preprocessing/Practice/new/sub-10471/anat/y_sub-10471_T1w.nii'};
        clear pathlist
        for ii= 1:nt    % Number of timepoints
            pathlist{ii,1} = strcat([fullfile(func.folder,['REG_',func.name]),',',num2str(ii)]);
        end
        matlabbatch{1}.spm.spatial.normalise.write.subj.resample = (pathlist);
        %%
        matlabbatch{1}.spm.spatial.normalise.write.woptions.bb = [-90 -126 -72;  %[-78 -112 -70;78 76 85];
                                                                     90 90 108];
        matlabbatch{1}.spm.spatial.normalise.write.woptions.vox = [2 2 2];
        matlabbatch{1}.spm.spatial.normalise.write.woptions.interp = 4;
        matlabbatch{1}.spm.spatial.normalise.write.woptions.prefix = 'w';
        
        cfg_util('run',matlabbatch); 
%         spm('defaults', 'FMRI');
%         spm_jobman('run', matlabbatch);
        clear matlabbatch
%% Smoothing

        % Load Regressed Images
        for imagei = 1:190
            reg_images{imagei, 1} = (fullfile(func.folder,['wREG_',func.name,',',num2str(imagei)]));
        end

        matlabbatch{1}.spm.spatial.smooth.data = reg_images;
        matlabbatch{1}.spm.spatial.smooth.fwhm = [8 8 8];
        matlabbatch{1}.spm.spatial.smooth.dtype = 0;
        matlabbatch{1}.spm.spatial.smooth.im = 0;
        matlabbatch{1}.spm.spatial.smooth.prefix = 's';

%         cfg_util('run',matlabbatch);  
        spm('defaults', 'FMRI');
        spm_jobman('run', matlabbatch);
        clear matlabbatch

    end
    subji_ex = subji_ex + 1;
    disp('################################################################################################################')
    disp(['Completed Preprocessing Subject ',Subj_list(subji).name])
    T(subji_ex) = toc;

    disp(['Time Taken for subject ',Subj_list(subji).name, ' is ',num2str(T(subji_ex) - T(subji_ex - 1)),' seconds'])
    disp('################################################################################################################')
end

disp('Preprocessing is Finally over, I can now rest and you can get to work :-) ')

%%
% To format the participants.tsv file

% input: file path---> the path of xlsx file


function T_formated = format_participants(filepath)

     T = readtable(filepath);
     T(:,6) = [];
     T_formated = unique(T,'rows');
end