2023_Blaschke_Stroke_tDCS_rsfMRI/scripts/addendum/register_data.m

parent_path = "E:\Projects\CBIrep_Imaging\subjects";
cd (parent_path)
addpath("..\scripts\toolbox\DVARS-master\")
addpath("..\scripts\toolbox\DVARS-master\Nifti_Util\")
addpath("..\scripts")
addpath("..\scripts\toolbox\spm12")
addpath("..\scripts\toolbox\fMRI-Quality-Checker-master")
addpath("..\scripts\toolbox\spm12\matlabbatch")
addpath("..\scripts\toolbox\nifti_utils-master\")

%% Preps
% px = dir('sub*');
BIDS_path = "E:\Projects\CBIrep_Imaging\subjects";
% if ~exist(BIDS_path,'dir')
%     mkdir(BIDS_path)
% end
% physio_files = "F:\Projects\StrokePT_tDCS_rsfmri\physio";
% cd (physio_files)
% phx = dir("SB*");
% timeX = {'tpA','tpB','tpC','tpD'};

% master = zeros300,5};
% noPhysio = zeros(300);
% master_tag = {'subject','group','timepoint','dataset','physio','volumes','bad_volumes','meanDeltaDvar'};
tick = 1;
tack = 1;

% Check for registration quality (user input)?
quali = "off";


%% level1: group folders
cd (parent_path)
gx = dir("sub*"); 

for g = 4
    group = gx(g).name;
    cd (fullfile(gx(g).folder,gx(g).name))

%% level2: timepoints
    tx = dir('tp*');
    for t = 1:length(tx)
        time = tx(t).name;
        cd (fullfile(tx(t).folder,tx(t).name))
        session = strcat('ses-',num2str(find(strcmp(time,timeX))));

        tix = dir("SB*")

%% level3: individual datasets
        for m = 1:length(tix)
            cd (fullfile(parent_path,group,time))
            dataset = tix(m).name;

%% find corresponding Physio dataset
            physioID = find(~cellfun('isempty',regexp({phx.name},dataset,'once')));
            if ~isempty(physioID)
                physiofolder = fullfile(physio_files,phx(physioID).name);
                master{tick,5} = "y";
            else
                noPhysio{tack} = dataset;
                tack = tack +1; 
                master{tick,5} = "n";
            end
%%        
            sjx = strfind(dataset,'_');
            subject = dataset(1:sjx(1)-1);
            
            master {tick,1} = subject; 
            master {tick,2} = group;
            master {tick,3} = num2str(t); 
            master {tick,4} = dataset;
          
            disp(strcat("running dataset:",dataset));
            cd(dataset)

%% copy data to BIDS
            folder_from = pwd;
            folder_to = fullfile(BIDS_path,strcat("sub-",subject(3:end)),session);
% % % %               
% % % %                 if ~exist(folder_to,"dir")
% % % %                     mkdir(folder_to)
% % % %                 end
% % % % %                 copy MRI data
% % % %                 copyfile(folder_from,folder_to)
% % % %                 
% % % % %                 copy physio data
% % % %             if exist("physiofolder","var")
% % % %                 physio_to = fullfile(folder_to,"physio",phx(physioID).name);
% % % %                 if ~exist(physio_to,"dir")
% % % %                     mkdir(physio_to)
% % % %                 end
% % % %                 copyfile(physiofolder,physio_to);
% % % %             end

%%
cd (folder_to)
% movefile("fMRI","func")
% movefile("T2w","anat")
 % a) get structural 
        if exist(fullfile(folder_to,"anat/"))==7 && exist (fullfile(folder_to,"func/"))==7
            cd (fullfile(folder_to,"anat/"))
            struct_orig = dir("*.1.nii.gz");
            struct_name = struct_orig(length(struct_orig)).name;  
            struct_fullfile = fullfile(struct_orig(length(struct_orig)).folder,struct_orig(length(struct_orig)).name); 
            
            struct_reg = dir("*BiasBet_AnnorsfMRI.nii.gz");
            if ~isempty(struct_reg)
            structReg_fullfile = fullfile(struct_reg.folder,struct_reg.name);
            
            struct_BET = dir("*BiasBet.nii.gz");
            structBET_fullfile = fullfile(struct_BET.folder,struct_BET.name);
            end
%             copyfile(struct_fullfile,fullfile(folder_to,"anat/",strcat(struct_name(1:end-7),'_T2w.nii.gz')))
            struct_orig = dir("*_T2w.nii.gz");
            struct_name = struct_orig(length(struct_orig)).name;  
            struct_fullfile = fullfile(struct_orig(length(struct_orig)).folder,struct_orig(length(struct_orig)).name); 
            

        % b) get fMRI
            cd (fullfile(folder_to,"func/"))
            fmri_orig = dir("*.1.nii.gz");
            clear ('fx_rep')
            for jx = 1:length(fmri_orig)
               fx_scan = niftiinfo(fmri_orig(jx).name);
               fx_rep(jx) = fx_scan.ImageSize(4);
            end
            fx_id = find(fx_rep>100);
            fx = length(fx_id);
% % % % 
% % % %             if fx>1
% % % %                 fmri_conc = niftiread(fmri_orig(fx_id(1)).name);
% % % %                 info = niftiinfo(fmri_orig(fx_id(1)).name);
% % % %                 info.ImageSize = [96,96,16,(fx*105-((fx-1)*5))];
% % % %                 for imgx = 2:fx
% % % %                   fmri_add = niftiread(fmri_orig(fx_id(1)).name);
% % % %                   fmri_conc = cat(4,fmri_conc,fmri_add(:,:,:,6:end)); 
% % % %                 end
% % % %                 conc_fmri = strcat(dataset,'.conc.nii');
% % % %                 niftiwrite(fmri_conc,conc_fmri,info);
% % % %                 gzip(conc_fmri)
% % % %                 fMRI_fullfile = strcat(fullfile(folder_to,'func',conc_fmri),".gz");
% % % %                 fMRI_file = dir('*conc.nii.gz');
% % % %             elseif fx==1
% % % %                 fMRI_fullfile = fullfile(fmri_orig(fx_id).folder,fmri_orig(fx_id).name);
% % % %                 fMRI_file = fmri_orig(fx_id);
% % % %             end
              fMRI_fullfile = fullfile(fmri_orig(fx_id).folder,fmri_orig(fx_id).name);
                fMRI_file = fmri_orig(fx_id);
            fMRI_name = fMRI_file.name;
%             copyfile(fMRI_fullfile,fullfile(fMRI_file.folder,strcat(fMRI_name(1:end-7),'_bold.nii.gz')))


            info = niftiinfo(fMRI_fullfile);
            master{tick,6} = info.ImageSize(4);

            fmri_smooth = dir("*SmoothBet.nii.gz");
            if ~isempty(fmri_smooth)
            fmriSmooth_fullfile = fullfile(fmri_smooth.folder,fmri_smooth.name);
            end
%             if ~exist("analyses","dir")
%                 mkdir("analyses")
%             end
%             cd ("analyses")

% % % % %% get DVARS
% % % %         V1 = load_untouch_nii(fMRI_file.name);
% % % %         V2 = V1.img;
% % % %         X0 = size(V2,1); Y0 = size(V2,2); Z0 = size(V2,3); T0 = size(V2,4);
% % % %         I0 = prod([X0,Y0,Z0]);
% % % %         Y  = reshape(V2,[I0,T0]); clear V2 V1;
% % % %         [DVARS,DVARS_Stat]=DVARSCalc(Y,'scale',1/100,'TransPower',1/3,'RDVARS','verbose',1);
% % % %         [V,DSE_Stat]=DSEvars(Y,'scale',1/100);
% % % % %             fMRIDiag_plot(V,DVARS_Stat)
% % % %     %         saveas(gcf,strcat(fmri_name(1:end-7),'_',"DVARS_plot.png"))
% % % %     %         pause(2)
% % % %     %         close all
% % % %         
% % % %         bad_points = find(DVARS_Stat.pvals<0.001);
% % % %         bad_points (bad_points<5)=[];
% % % %         perc_bad = length(bad_points)/(info.ImageSize(4)-5)*100;
% % % %             disp("amount of bad points: ")
% % % %             disp(num2str(perc_bad))
% % % %             disp("%")
% % % % 
% % % % %         master{tick,6} = fmri_name;
% % % % %         bad_points_list{tick,2} = bad_points;
% % % %         master{tick,8} = mean(DSE_Stat.DeltapDvar(6:end),'omitnan');
% % % % %         bad_points_list{tick,5} = mean(DSE_Stat.DeltapSvar(6:end));
% % % %         master{tick,7} = perc_bad;

        % c) get Regression File
            cd (fullfile(folder_to,"func/"))
            ABAreg_orig = dir("*SmoothBet_AnnoSplit_rsfMRI.nii.gz");
            if ~isempty(ABAreg_orig)
            ABAreg_fullfile = fullfile(ABAreg_orig.folder,ABAreg_orig.name);
            end
            if exist("regr","dir")
                cd regr\
                SFRGR_orig = dir("*_SFRGR.nii.gz");
                if ~isempty(SFRGR_orig)
                SFRGR_fullfile = fullfile(SFRGR_orig.folder,SFRGR_orig.name);
                else
                    SFRGR_fullfile = [];
                end
            end
%% plot raw images to check registration
     if ~exist("quali",'var')
         quali = "off"
     end
     if quali == "on"                 
            fMRI_imag = niftiread(fmriSmooth_fullfile);
            fMRI_imag = imrotate(fMRI_imag,-90);
            
            if ~isempty(ABAreg_orig)
            fMRI_Reg = niftiread(ABAreg_fullfile);
            fMRI_Reg = imrotate(fMRI_Reg,-90);
            else
                fMRI_Reg = zeros(96,96,16);
            end

            T2raw = niftiread(structBET_fullfile);
            T2raw = imrotate(T2raw,-90);
            if ~isempty(struct_reg)
            T2Reg = niftiread(structReg_fullfile);
            T2Reg = imrotate(T2Reg,-90);
            else
                T2Reg = zeros(256,256,28);
            end
            
            Loc = [8 13];

            figure('units','normalized','outerposition',[0 0 1 1])
                        for ix = 1:2
                            % Raw T2w
                            subplot(2,4,1+(4*(ix-1)))
                                imagesc(T2raw(:,:,floor(Loc(ix)/16*48)))
                                colormap("hot")
                            % T2 ARA Reg
                            subplot(2,4,2+(4*(ix-1)))
                                imagesc(T2Reg(:,:,floor(Loc(ix)/16*48)))
%                                 colormap("hot") 
                            % Raw fMRI
                            subplot(2,4,3+(4*(ix-1)))
                                imagesc(fMRI_imag(:,:,Loc(ix)))
%                                 colormap("gray")
                            % fMRI ARA Reg
                            subplot(2,4,4+(4*(ix-1)))
                                imagesc(fMRI_Reg(:,:,Loc(ix)))
%                                 colormap("hot")
                        end
                        title(strcat(subject,'-',time))
             pause
             close
  %% user input to evaluate quality
             RegQuali{tick,1}=subject;
             RegQuali{tick,2}=time;
             
             % get input of image quality
             checkT2    = input("is T2 ok? Y/N [Y]");
             if isempty(checkT2)
                 checkT2 = "Y";
             else 
                 checkT2 = "N";
             end
             RegQuali{tick,3} = checkT2;
             
             checkfMRI  = input("is fMRI ok? Y/N [Y]");
             if isempty(checkfMRI)
                 checkfMRI = "Y";
             else 
                 checkfMRI = "N";
             end
             disp('')
             RegQuali{tick,4} = checkfMRI;
     end
%%

             tick=tick +1;




        end
               

    end

end
end



% %         % a) get structural 
% %         if exist(fullfile(parent_path,subject,"MRI",time,"T2w/"))==7 && exist (fullfile(parent_path,subject,"MRI",time,"fMRI/"))==7
% %             cd (fullfile(parent_path,subject,"MRI",time,"T2w/"))
% %             struct_orig = dir("*.1.nii.gz");
% %             struct_name = struct_orig(length(struct_orig)).name;  
% %             struct_fullfile = fullfile(struct_orig(length(struct_orig)).folder,struct_orig(length(struct_orig)).name); 
% %             
% %             struct_reg = dir("*BiasBet_AnnorsfMRI.nii.gz");
% %             if ~isempty(struct_reg)
% %             structReg_fullfile = fullfile(struct_reg.folder,struct_reg.name);
% %             end
% %             struct_BET = dir("*BiasBet.nii.gz");
% %             structBET_fullfile = fullfile(struct_BET.folder,struct_BET.name);
% % 
% % 
% %         % b) get fMRI
% %             cd (fullfile(parent_path,subject,"MRI",time,"fMRI/"))
% %             fmri_orig = dir("*.1.nii.gz");
% %             fmri_name = fmri_orig(length(fmri_orig)).name;
% %             fMRI_fullfile = fullfile(fmri_orig(length(fmri_orig)).folder,fmri_orig(length(fmri_orig)).name); 
% %             fmri_smooth = dir("*SmoothBet.nii.gz");
% %             fmriSmooth_fullfile = fullfile(fmri_smooth.folder,fmri_smooth.name);
% % 
% %             if ~exist("analyses","dir")
% %                 mkdir("analyses")
% %             end
% %             cd ("analyses")
% % 
% %         % c) get Regression File
% %             cd (fullfile(parent_path,subject,"MRI",time,"fMRI/"))
% %             ABAreg_orig = dir("*SmoothBet_AnnoSplit_rsfMRI.nii.gz");
% %             if ~isempty(ABAreg_orig)
% %             ABAreg_fullfile = fullfile(ABAreg_orig.folder,ABAreg_orig.name);
% %             end
% %             cd regr\
% %             SFRGR_orig = dir("*_SFRGR.nii.gz");
% %             if ~isempty(SFRGR_orig)
% %             SFRGR_fullfile = fullfile(SFRGR_orig.folder,SFRGR_orig.name);
% %             else
% %                 SFRGR_fullfile = [];
% %             end
% %         % plot raw images to check registration
% %             
% %              
% %             fMRI_imag = niftiread(fmriSmooth_fullfile);
% %             fMRI_imag = imrotate(fMRI_imag,-90);
% %             
% %             if ~isempty(ABAreg_orig)
% %             fMRI_Reg = niftiread(ABAreg_fullfile);
% %             fMRI_Reg = imrotate(fMRI_Reg,-90);
% %             else
% %                 fMRI_Reg = zeros(96,96,16);
% %             end
% % 
% %             T2raw = niftiread(structBET_fullfile);
% %             T2raw = imrotate(T2raw,-90);
% %             if ~isempty(struct_reg)
% %             T2Reg = niftiread(structReg_fullfile);
% %             T2Reg = imrotate(T2Reg,-90);
% %             else
% %                 T2Reg = zeros(256,256,28);
% %             end
% %             
% %             Loc = [8 13];
% % %%
% %                         figure('units','normalized','outerposition',[0 0 1 1])
% %                         for ix = 1:2
% %                             % Raw T2w
% %                             subplot(2,4,1+(4*(ix-1)))
% %                                 imagesc(T2raw(:,:,floor(Loc(ix)/16*28)))
% %                                 colormap("hot")
% %                             % T2 ARA Reg
% %                             subplot(2,4,2+(4*(ix-1)))
% %                                 imagesc(T2Reg(:,:,floor(Loc(ix)/16*28)))
% % %                                 colormap("hot") 
% %                             % Raw fMRI
% %                             subplot(2,4,3+(4*(ix-1)))
% %                                 imagesc(fMRI_imag(:,:,Loc(ix)))
% % %                                 colormap("gray")
% %                             % fMRI ARA Reg
% %                             subplot(2,4,4+(4*(ix-1)))
% %                                 imagesc(fMRI_Reg(:,:,Loc(ix)))
% % %                                 colormap("hot")
% %                         end
% %                         title(strcat(subject,'-',time))
% %              pause
% %              close
% %   %%
% %              RegQuali{tick,1}=subject;
% %              RegQuali{tick,2}=time;
% %              
% %              % get input of image quality
% %              checkT2    = input("is T2 ok? Y/N [Y]");
% %              if isempty(checkT2)
% %                  checkT2 = "Y";
% %              else 
% %                  checkT2 = "N";
% %              end
% %              RegQuali{tick,3} = checkT2;
% %              
% %              checkfMRI  = input("is fMRI ok? Y/N [Y]");
% %              if isempty(checkfMRI)
% %                  checkfMRI = "Y";
% %              else 
% %                  checkfMRI = "N";
% %              end
% %              disp('')
% %              RegQuali{tick,4} = checkfMRI;
% % 
% %              tick = tick + 1;