2023_Blaschke_Stroke_tDCS_rsfMRI/scripts/addendum/checkMRI.m

% function MRI_quality(BIDS_path)
% Stefan J. Blaschke, Neurology Department, University Hospital Cologne
% VERSION 1.0., 31.12.2022
%
% matlab script to examine preprocessed fMRI dataset
% display of raw imaging
% overlay with registered atlas
% SNR and tSNR
% display of tSNR with mask outlines
% DVARS and FD
% -  Afyouni S. & Nichols T.E., Insights and inference for DVARS, 2017
%   http://www.biorxiv.org/content/early/2017/04/06/125021.1
% necessity of GSR (Chen et al., 2022, Magn Reson Med)


% header - set up path to BIDS folder
head_path = 'G:\Projects\StrokePT_tDCS_rsfmri'; % 
rate_data = 1;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% run script 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd (head_path)
    addpath(".\scripts\toolbox\DVARS-master\")
    addpath(".\scripts\toolbox\DVARS-master\Nifti_Util\")
    addpath(".\scripts\toolbox\BCT\")
    addpath(".\scripts\toolbox\spm12")
    addpath(".\scripts\toolbox\fMRI-Quality-Checker-master")
    addpath(".\scripts\toolbox\spm12\matlabbatch")
    addpath(".\scripts\toolbox\nifti_utils-master\")
    addpath(".\scripts\toolbox\nifti_utils-master\external\NIFTI")
    addpath(".\scripts\toolbox\FSLNets\")
    addpath(".\scripts\")
%%
% create data table to store features

varnames = {'animal','ses','ID','MRI date','file path','SNR','tSNR','DVARS','bad points','IOR','MNS','physio','qualy'};
vartype = repmat({'cellstr'},length(varnames),1);
Tmaster = table('Size',[100 length(varnames)],'VariableTypes',vartype,'VariableNames',varnames);

parent_path = fullfile(head_path,"subjects/"); % parent_path = BIDS MRI path
cd (parent_path)
px = dir('sub*');
%%
% loop through Data
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%
idx = 1;
for sub_idx = 1:length(px) % start loop subjects

    animal_path = fullfile(parent_path,px(sub_idx).name);
    cd (animal_path)
    tx = dir('ses*');
 
    for ses_idx = 1:length(tx) % start loop sessions
%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%

       session_path = fullfile(animal_path,tx(ses_idx).name);
        cd (session_path)

[path,session,~] = fileparts(session_path);
[~,animal,~] = fileparts(path);

Tmaster{idx,"animal"} = {animal};
Tmaster{idx,"ses"} = {session};
Mat(idx).animal = animal;
Mat(idx).ses = session;
fprintf('Run# %u: running dataset from %s, %s...',idx,animal,session)
 

% grab images
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if exist(fullfile(session_path,"anat/"))==7
% a) get structural 
cd (fullfile(session_path,"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_parts = strsplit(struct_orig(length(struct_orig)).name,'_');
    Tmaster(idx,'ID') = struct_parts(2);
    Mat(idx).ID = struct_parts(2);
    Tmaster(idx,'MRI date') = struct_parts(1);
    Mat(idx).date = struct_parts(1);
 
% b.1) get fMRI
if exist(fullfile(session_path,"func/"))==7
    cd (fullfile(session_path,"func/"))
        fmri_orig = dir("*.1.nii.gz");
        if length(fmri_orig)>1
%             for lx = 1:length(fmri_orig)
%                 sx = niftiinfo(fmri_orig(lx).name);
%                 volumes(lx) = sx.ImageSize(4);
%             end
            image_size = [fmri_orig.bytes];
%             image_size(volumes<100) = 0;
            fmri_orig = fmri_orig(find(image_size==max([image_size])));
        end
         
        ANNO_orig = dir("*AnnoSplit_rsfMRI.nii.gz");
        if ~isempty(fmri_orig) && ~isempty(ANNO_orig)
        fmri_name = fmri_orig(length(fmri_orig)).name;
        fMRI_fullfile = fullfile(fmri_orig(length(fmri_orig)).folder,fmri_orig(length(fmri_orig)).name); 
        ANNO_fullfile = fullfile(ANNO_orig(length(ANNO_orig)).folder,ANNO_orig(length(ANNO_orig)).name);
        atlas = load_untouch_nii(ANNO_fullfile);
        atlas_img = imrotate(atlas.img,-90);
        M = atlas_img;
        M(M>0) = 1;

% check physio regression
if exist("rawMonData","dir")==7
    cd ("rawMonData\")
    if ~isempty(dir('*I32'))
        Tmaster{idx,'physio'} = {'y'};
    else
        Tmaster{idx,'physio'} = {'n'};
    end
else
     Tmaster{idx,'physio'} = {'n'};
end

% b.2) get Regression File
    cd (fullfile(session_path,"func/"))
    save_path = fullfile(session_path,"func/","results/")
if ~exist(save_path)
    mkdir(save_path)
end

if exist('regr','dir')==7
        cd regr\
        SFRGR_orig = dir("*_SFRGR.nii.gz");
        TCfile = dir("MasksTCsSplit*.txt.mat");
        fmri_mask = dir("*thres_mask.nii.gz");
        
if ~isempty(TCfile) && ~isempty(SFRGR_orig)
        SFRGR_file = char(fullfile(SFRGR_orig.folder,SFRGR_orig.name));
        TC_fullfile = fullfile(TCfile.folder,TCfile.name);
        Tmaster(idx,'file path') = {TC_fullfile};

% load and plot SFRGR file        
SFRGR = load_untouch_nii(SFRGR_file);
    SFRGR_img = SFRGR.img;
    SFRGR_mean = imrotate(squeeze(mean(SFRGR_img,4)),-90);
    
    
% load and add boundaries of mask
mask = load_untouch_nii(fullfile(fmri_mask.folder,fmri_mask.name));
    mask_img = imrotate(mask.img,-90);
    mask_img(mask_img>0) = 1;

% % % Global signal regression
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RestData = gunzip(fullfile(SFRGR_orig.folder,SFRGR_orig.name));
%     WB_mask = gunzip(fullfile(fmri_mask.folder,fmri_mask.name));
%     GNI = Rest2GNI(char(RestData),char(WB_mask),100);
%     if GNI >3
%         need_GSR = "GSR";
%     else
%         need_GSR = "no GSR";
%     end
% 
% flag{1} = need_GSR;

% % % SNR (spatial, temporal)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   info = niftiinfo(fMRI_fullfile);
            V1 = load_untouch_nii(fMRI_fullfile);
            V2 = V1.img;
            fmri_data = flip(flip(imrotate(V2,-90),3),2);
            [Ni, Nj, Nk, Nt] = size(fmri_data);
            func_2Dimg = reshape(fmri_data, Ni*Nj*Nk, Nt);
            [Nv, Nt] = size(func_2Dimg);

            func_mean = mean(func_2Dimg, 2,'omitnan');
            func_mean3D = reshape(func_mean, Ni, Nj, Nk);
            func_stddev = std(double(func_2Dimg), 0, 2,'omitnan');
            func_stddev3D = reshape(func_stddev, Ni, Nj, Nk);
            
            tSNR_2Dimg = mean(func_2Dimg, 2,'omitnan')./std(double(func_2Dimg), 0, 2,'omitnan');
            tSNR_3Dimg = reshape(tSNR_2Dimg, Ni, Nj, Nk);
       
            tSNR = squeeze(mean2(tSNR_3Dimg(mask_img==1)));
            SNR = squeeze(mean(func_mean3D(mask_img==1)))/std(func_stddev3D(mask_img==0));
Tmaster{idx,'tSNR'} = {tSNR};
Tmaster{idx,'SNR'} = {SNR};
Mat(idx).tSNR = tSNR;
Mat(idx).SNR = SNR; 

%             cd ..\
%             niftiwrite(tSNR_image,strcat(fmri_name(1:end-10),'tSNR.nii'))

SFRGR_img(mask_img==0)=nan;
GSR = squeeze(mean(mean(mean(SFRGR_img,3,'omitnan'),2,'omitnan'),1,'omitnan'));
% % 
% % % get DVARS
        fmri_BET = fmri_data;
        fmri_BET(mask_img==0) = nan;
        [X0,Y0,Z0,T0] = size(fmri_BET); 
        I0 = prod([X0,Y0,Z0]);
        Y  = reshape(fmri_BET,[I0,T0]);
        [DVARS,DVARS_Stat]=DVARSCalc(Y,'scale',1/100,'TransPower',1/3,'RDVARS','verbose',0);
        [V,DSE_Stat]=DSEvars(Y,'scale',1/100);
        
        Tmaster{idx,'DVARS'} = {mean(DVARS)};
        Dx = find(DVARS_Stat.pvals<0.05./(T0-1) & DVARS_Stat.DeltapDvar>5); %print corrupted DVARS data-points
        
        bad_points = Dx;%(Dx > 3);

        if isempty(bad_points)
            badTP = 0;
        else
            badTP = length(bad_points);
        end

        scratch = zeros(1,4*length(bad_points));
        for ptx = 1:length(bad_points)
            scratch(((ptx-1)*4+1):((ptx-1)*4+4)) = bad_points(ptx)-2:bad_points(ptx)+1;
        end
        scratch = unique(scratch(scratch>0));
        scratch(scratch>info.ImageSize(4)-5) = [];

        Tmaster{idx,"bad points"} = {badTP*100/(info.ImageSize(4)-5)};
        Mat(idx).badvol = {badTP*100/(info.ImageSize(4)-5)};
        Mat(idx).DVARS = DVARS;
  
 % getFD
        cd (fullfile(session_path,"func/"))
        if exist("txtRegrPython","dir")
        cd txtRegrPython\
        txtRegr_files = dir("*mcf*");
            if size(txtRegr_files,1)>1
            txtRegr = txtRegr_files(floor(length(txtRegr_files)/2)).name;
            txtRegr_table = readtable(txtRegr,'VariableNamingRule','preserve');
                if ~isempty(txtRegr_table)
                    Regression_parameters = txtRegr_table(:,5:10);
                    MP = cell2mat(table2cell(Regression_parameters));
                    FD_measures = calculateFD(MP, .5, 0.5);
                    FD = FD_measures.FD;
                    Tmaster{idx,'FD'}=squeeze(mean(FD));
                    Mat(idx).FD = FD;
                   
                else
                    FD = [];
                    Mat(idx).FD = [];
                    disp('Dataset has no available regression file')
                end
            else 
                FD = [];
                Mat(idx).FD = [];
                disp('Dataset has no available regression file')
            end
        else
                FD = [];
                Mat(idx).FD = [];
                disp('Dataset has no available regression file')
        end
        
[mNet,mNet_scrub,mNet_scrub_GSR,IOR] = grabMRI(TC_fullfile,GSR,scratch,rate_data);

Tmaster{idx,'IOR'} = {IOR};
Mat(idx).IOR = IOR;
MNS = squeeze(mean(strengths_und(mNet)));
Tmaster{idx,'MNS'} = {MNS/(size(mNet,1)-1)};

Mat(idx).mNet = mNet;
Mat(idx).mNet_scrub = mNet_scrub;
Mat(idx).mNet_scrub_GSR = mNet_scrub_GSR;

T_sub = Tmaster(idx,["animal","ses","tSNR","SNR","DVARS","bad points","IOR","MNS"]);

%% plotting and rating data
%%%%%%%%%%%%%%%%%%%%%%%%%%
if rate_data == 1
pause()
saveas(gcf,fullfile(save_path,strcat(animal,'_',session,'_plots.png')))
close
    structural_fn = char(gunzip(struct_fullfile));
    vc= load_untouch_nii(structural_fn);
    vc = imrotate(vc.img,-90);
    figure(2)
    montage(vc)
    clim([0 max(max(max(vc)))])
pause()

functional4D_fn = char(gunzip(fMRI_fullfile));
    vc= load_untouch_nii(functional4D_fn);
    vc = imrotate(squeeze(mean(vc.img,4)),-90);
    figure(3)
    montage(vc)
    clim([0 max(max(max(vc)))])
pause()

plotMRI(tSNR_3Dimg,mask_img,'hot')
pause()

    plotMRI(SFRGR_mean,M)
pause()

figure('units','normalized','outerposition',[0 0 1 1])
   subplot(10,3,4:3:30)
   x = find(mNet~=0);
   [c,stat] = cdfplot(abs(mNet(x)));
   c.LineWidth = 2;
   c.Color = [0 0 1];
   hold on
   x = find(mNet_scrub_GSR~=0);
   [c,stat] = cdfplot(abs(mNet_scrub_GSR(x)));
   xlim([0 1])
   c.LineWidth = 2;
   c.Color = [1 0 0];
   legend({'raw Mat','GSR Mat'},'Location','southeast')
   title('')
   ylabel('cumulative distribution')
   xlabel('node strength')
   hold off
  
   subplot(10,3,5:3:30)
        mNetD = threshold_proportional(mNet,0.2);
%         nNetT = zeros(size(mNetD));
%         mNetT(mNetD~=0) = mNet_full(mNetD~=0)

        plotmat_values_rsfMRI(mNetD,[],[-(stat.mean*3) stat.mean*3],0,[]);
        hold on
        axis square;
            plot(get(gca,'xlim'),[mean(get(gca,'ylim')), mean(get(gca,'ylim'))],'k-','linewidth',2);
            plot([mean(get(gca,'xlim')), mean(get(gca,'xlim'))], get(gca,'ylim'),'k-','linewidth',2);
       
   subplot(10,3,6:3:30)
   GSRx=GSR;
   GSRx(scratch)=nan;
    plot(GSRx,'r-','LineWidth',2)
    xlabel('timepoint')
    ylabel('value')
    title('globel signal')


%Get the table in string form.
TString = evalc('disp(T_sub)');
% Use TeX Markup for bold formatting and underscores.
TString = strrep(TString,'<strong>','\bf');
TString = strrep(TString,'</strong>','\rm');
TString = strrep(TString,'_','\_');
% Get a fixed-width font.
FixedWidth = get(0,'FixedWidthFontName');
% Output the table using the annotation command.
annotation(gcf,'Textbox','String',TString,'Interpreter','Tex',...
    'FontName',FixedWidth,'Units','Normalized','Position',[0 0 1 1]);

    pause()
saveas(gcf,fullfile(save_path,strcat(animal,'_',session,'_parameters.png')))


    % dialog to rate preprocessing
    qm_set = 0;
    while qm_set == 0
    qm = input("was registration successful? y/n:",'s');
        if isempty(qm) | length(qm)>1
             disp('please print y or n' )
        elseif qm ~= 'n' && qm ~= 'y' && qm ~= 'd' &&  qm ~= 'g' && qm ~= 'x'
            disp('please print y or n' )
        else
           qm_set = 1;
        end
    end
    close all
end
                    
                    else
                        qm = 0;
                    end % if loop no Anno?
                else
                    qm = 0;
                end % if loop no anat?
            else
                qm = 0;
            end % if loop no func?  
        else 
            qm = 0;
        end % if loop no TCfile?
else
    qm = 0;
end % if loop no regression?
if ~exist('qm','var')
    qm = '';
end
Tmaster{idx,'qualy'} = {qm}; 
Mat(idx).qualy = qm;


%%
clearvars -except head_path idx Mat Tmaster parent_path px tx rate_data ses_idx sub_idx animal_path
idx = idx +1;
    end % end loop sessions
end % end loop subjects

%%
cd (head_path)
writetable(Tmaster,'PT_tDCS_quali_stat_rate20230201.xlsx')
save("PT_tDCS_Mat_rate_20230201","Mat")
%%
fprintf('idx = %u\n',idx)
fprintf('subject: %u\n',sub_idx)
fprintf('session: %u\n',ses_idx)
% 3.1.23 11:27: idx = 206, sub 46, ses 1