% Created by Ari Segel, for any questions, email me at ari@wustl.edu
% This function was designed to facilitate a cleaner top level processing script
% This function currently makes GLM events and no events maps for the following paradigms
    % HW, OV, CV, GV, MOT, DFS, AC, OBS, IM, ISS, SLP, SS, and MS
% This function displays maps in three views for all of the above paradigms, 
    % typically view setting is interactable with NeuroDOT, but this function
    % manually sets the views for ease of use purposes


function [map, b, e, DM, EDM] = make_GLM_maps_260313_AS(dataIn, anatomyIn, hrf, affine, FOV, params, run)

    %% Parameters and Initialization
    if ~exist('params', 'var'), params=struct; end
    if ~isfield(params, 'GVTD_censor'), params.GVTD_censor = 1; end % do you want to threshold motion? 1 = yes, 0 = no
    if ~isfield(params, 'GVTD_th'), params.GVTD_th = 0.007; end % GVTD threshold for motion censoring
    if ~isfield(params, 'DoFilter'), params.DoFilter = 1; end % Filtering used in GLM
    if ~isfield(params, 'ctx'), params.ctx = 'std'; end % Cortex inflation
    if ~isfield(params, 'scale_factor'), params.scale_factor = 0.8; end % Multiply this number by std of data to get colorscale (happens below)
    if ~isfield(params, 'Use_hardcode_scale'), params.Use_hardcode_scale = 0; end % Flag for if user hardcodes colorscale outside of this function 
    if ~isfield(params, 'OL'), params.OL = 0; end % Overlap, advanced parameter, see NeuroDOT Data Visualization Tutorial for more on how to use this
    if ~isfield(params, 'Th')  ||  isempty(params.Th) % Color Scale thresholds
        params.Th.P = 0;
        params.Th.N = -params.Th.P;
    end
    if ~isfield(params, 'save_map'), params.save_map = 1; end % Save the figure? 1=yes 0=no
    if ~isfield(params, 'kill'), params.kill = 1; end % Close the figure after saving? 1=yes 0=no
    
    
    %% Run GLM
    [b,e,DM,EDM]=GLM_181206(squeeze(dataIn.cortex_Hb(:,:,1)),hrf,dataIn.info,params); % GLM code
    
    % File saving
    datafile=[params.flags.subj,'_',params.flags.dateA,'_',run]; % String for file name
    if params.GVTD_censor == 1
        if params.save_map == 1
            SaveGCF_to_PNG([datafile,'_GLM_DM_',params.flags.tag],params.kill);
        end
        GVTD_thresh_string = num2str(params.GVTD_th);
    end
    if params.GVTD_censor == 0 && params.kill == 1 %close GLM DM when not thresholding data - added option with kill to prevent closing
        close
        GVTD_thresh_string = 'none';
    elseif params.GVTD_censor == 0
        GVTD_thresh_string = 'none';
    end

    
    
    %% Get run_type
    run_name = run(isstrprop(run,'alpha')); % Get string for run Name
    
    if strcmp(run_name, 'HW') || strcmp(run_name, 'OV') || strcmp(run_name, 'CV') || strcmp(run_name, 'GV')...
            || strcmp(run_name, 'MS') || strcmp(run_name, 'SS')
        run_type = 'HW_GV_OV_CV';
    elseif strcmp(run_name, 'MOT') || strcmp(run_name, 'DFS') || strcmp(run_name, 'AC')...
            || strcmp(run_name, 'MOTshort') || strcmp(run_name, 'DFSshort')
        run_type = 'MOT_DFS_AC';
    elseif strcmp(run_name, 'OBS') || strcmp(run_name, 'IM') || strcmp(run_name, 'ISS')
        run_type = 'OBS_IM_ISS';
    elseif strcmp(run_name, 'SLP') || strcmp(run_name, 'MV')
        run_type = 'SLP';
%     elseif strcmp(run_name, 'MS') || strcmp(run_name, 'SS')
%         run_type = 'MS_SS';
    else
        disp('Unsupported run type used as an input. No maps will be made.')
    end
        
    
    
    %% Switch Case for Map Making
    switch run_type
        case 'HW_GV_OV_CV'
            %% Hearing Words, Verb Generation
            
            % Sort betas by event type
            if strcmp(params.type, 'no_events')
                b0 = b(:,2:end);
                b_rest = b0(:, dataIn.info.paradigm.Pulse_1);
                b_words = b0(:, dataIn.info.paradigm.Pulse_2);
            end
            
            switch params.type
                case 'events'
                    % Get Beta Map
                    bL = (squeeze(b(:,3))-squeeze(b(:,2))); % Task - Rest
                    bLvol = Good_Vox2vol(bL,dataIn.info.tissue.dim);
                    map = affine3d_img(bLvol,dataIn.info.tissue.dim,anatomyIn.info,affine); % Affine transform data to MNI space
                    if contains(params.flags.cap, 'LUMO')
                        map = map.*FOV; % Multiply by FOV, only for LUMO bc A was made on scaled mesh (may change this section in the future)
                    end

                    % Visualization Parameters
                    if length(params.scale_factor) > 1 % For these paradigms, scale factor should never be more than 1 value since we're only making one map
                        disp(['Scale factor set improperly, please only set one scale factor for ' run ' runs'])
                        params.scale_factor = 3.5; % Set to 3.5 bc we're using std now
                    end
                    params.Scale = params.scale_factor*std(bL(bL(:)~=0));
                    if params.Use_hardcode_scale == 1 %use hardcoded scale set by user if this flag is turned on
                            params.Scale = params.Scale_hardcoded(j);
                    end

                case 'no_events'
                    % Do T-test and get cohen's D
                    [~,~,~,d ]= TtestWelch(b_words, b_rest);
                    dvol=Good_Vox2vol(d,dataIn.info.tissue.dim);
                    map = affine3d_img(dvol,dataIn.info.tissue.dim,anatomyIn.info,affine); % Affine transform data to MNI space
                    if contains(params.flags.cap, 'LUMO')
                        map = map.*FOV; % Multiply by FOV, only for LUMO bc A was made on scaled mesh (may change this section in the future)
                    end

                    % Visualization Parameters
                    if length(params.scale_factor) > 1 % For these paradigms, scale factor should never be more than 1 value since we're only making one map
                        disp(['Scale factor set improperly, please only set one scale factor for ' run ' runs'])
                        params.scale_factor = 3.5; % Set to 3.5 bc we're using std now
                    end
                    params.Scale=params.scale_factor(j)*std(dvol(dvol(:)~=0)); %color scale parameters, you might want to adjust these  
                    if params.Use_hardcode_scale == 1 %use hardcoded scale set by user if this flag is turned on
                            params.Scale = params.Scale_hardcoded(j);
                    end
            end
            
            % Visualize
            % This will make a single figure with all 3 views of the brain
            figure('Color','k','Position',[100,100,1200,800]);
            %lateral view 
            subplot(2,2,[1,2]);
            params.fig_handle=gca;
            params.view='lat'; 
            PlotInterpSurfMesh(map,anatomyIn.MNIl,anatomyIn.MNIr,anatomyIn.info,params);
            title([{[run, ' GLM with ', params.type]};{['GVTD.Th=' GVTD_thresh_string ' Scale=' num2str(params.Scale)]};...
                     {[' # time pts kept=' num2str(size(DM,1)) ]}],'Color','w', 'FontSize', 8);
            %dorsal view 
            subplot(2,2,[3]);
            params.fig_handle=gca;
            params.view='dorsal';
            PlotInterpSurfMesh(map,anatomyIn.MNIl,anatomyIn.MNIr,anatomyIn.info,params);
            title([{[run, ' GLM with ', params.type]};{['GVTD.Th=' GVTD_thresh_string ' Scale=' num2str(params.Scale)]};...
                    {[' # time pts kept=' num2str(size(DM,1)) ]}],'Color','w', 'FontSize', 8);
            %posterior view 
            subplot(2,2,[4]);
            params.fig_handle=gca;
            params.view='post';
            PlotInterpSurfMesh(map,anatomyIn.MNIl,anatomyIn.MNIr,anatomyIn.info,params);
            title([{[run, ' GLM with ', params.type]};{['GVTD.Th=' GVTD_thresh_string ' Scale=' num2str(params.Scale)]};...
                    {[' # time pts kept=' num2str(size(DM,1)) ]}],'Color','w', 'FontSize', 8);
        
        case 'MOT_DFS_AC'
            %% Motor and Alternating Checkerboards
         
            
            % Setup
            figure('Color','k','Position',[100,0,1200,800]);
            type = {'left', 'right', 'left-right'};
            if length(params.scale_factor) < 3 % Make sure scale factor is not long bc this uses the same scale for all 3 views
                disp(['Only ' num2str(length(params.scale_factor)) ' scale factor given. Using the same scale factor for each map'])
                params.scale_factor = repmat(params.scale_factor, [1,length(type)]); % set to 3.5 because for these paradigms we use std of data instead of std of data to set colorscale
            end

            % Sort betas by event type
            if strcmp(params.type, 'no_events')
                b0 = b(:,2:end);
                b_rest = b0(:, dataIn.info.paradigm.Pulse_1);
                b_left = b0(:, dataIn.info.paradigm.Pulse_2);
                b_right = b0(:, dataIn.info.paradigm.Pulse_3);
            end

            % Map making loop
            for j = 1:length(type)

                stim = type{j};

                switch params.type
                    case 'events'
                        switch stim
                            case 'left'
                                % Make maps of left stimulus
                                bL=squeeze(b(:,3)-b(:,2));
                                bvol=Good_Vox2vol(bL,dataIn.info.tissue.dim);
                                params.Scale=params.scale_factor(j)*std(bvol(bvol(:)~=0)); %color scale parameters, you might want to adjust these
                            case 'right'
                                % Make maps of right stimulus
                                bR=squeeze(b(:,4)-b(:,2));
                                bvol=Good_Vox2vol(bR,dataIn.info.tissue.dim);
                                params.Scale= params.scale_factor(j)*std(bvol(bvol(:)~=0)); %color scale parameters, you might want to adjust these
                            case 'left-right'
                                % Make maps of contrast (left - right)
                                bC=squeeze(b(:,3)-b(:,4));
                                bvol=Good_Vox2vol(bC,dataIn.info.tissue.dim);
                                params.Scale=params.scale_factor(j)*std(bvol(bvol(:)~=0)); %color scale parameters, you might want to adjust these
                        end

                        % Affine transform data to MNI space
                        map=affine3d_img(bvol,dataIn.info.tissue.dim,anatomyIn.info,affine);

                    case 'no_events'
                            switch stim
                                case 'left'
                                    % Make maps of left stimulus
                                    [~,~,~,d ]= TtestWelch(b_left, b_rest);
                                    dvol=Good_Vox2vol(d,dataIn.info.tissue.dim);
                                    params.Scale=params.scale_factor(j)*std(dvol(dvol(:)~=0)); %color scale parameters, you might want to adjust these  
                                case 'right'
                                    % Make maps of right stimulus
                                    [~,~,~,d ]= TtestWelch(b_right, b_rest);
                                    dvol=Good_Vox2vol(d,dataIn.info.tissue.dim);
                                    params.Scale=params.scale_factor(j)*std(dvol(dvol(:)~=0)); %color scale parameters, you might want to adjust these  

                                case 'left-right'
                                    % Make maps of contrast (left - right)
                                    [~,~,~,d ]= TtestWelch(b_left, b_right);
                                    dvol=Good_Vox2vol(d,dataIn.info.tissue.dim);
                                    params.Scale=params.scale_factor(j)*std(dvol(dvol(:)~=0)); %color scale parameters, you might want to adjust these  
                            end

                            % Affine transform data to MNI space
                            map=affine3d_img(dvol,dataIn.info.tissue.dim,anatomyIn.info,affine);
                end

                if params.Use_hardcode_scale == 1 %use hardcoded scale set by user if this flag is turned on
                    params.Scale = params.Scale_hardcoded(j);
                end

                % Multiply by FOV, only for LUMO bc A was made on scaled mesh (may change this section in the future)
                if contains(params.flags.cap, 'LUMO')
                    map = map.*FOV; 
                end

                % Plot each view
                views = {'lat', 'post', 'dorsal'};
                num_views = length(views);
                subplot(3,num_views+1,j*4-3:j*4-2);
                params.fig_handle=gca;
                params.view = 'lat';
                PlotInterpSurfMesh(map,anatomyIn.MNIl,anatomyIn.MNIr,anatomyIn.info,params);
                title([{[run,' ', stim, ' GLM with ', params.type]};{['GVTD.Th=' GVTD_thresh_string ' Scale=' num2str(params.Scale)]};...
                     {[' # time pts kept=' num2str(size(DM,1)) ]}],'Color','w', 'FontSize', 8, 'interpreter', 'none');
                subplot(3,num_views+1,j.*4-1);
                params.fig_handle=gca;
                params.view = 'post';
                PlotInterpSurfMesh(map,anatomyIn.MNIl,anatomyIn.MNIr,anatomyIn.info,params);
                title([{[run,' ', stim, ' GLM with ', params.type]};{['GVTD.Th=' GVTD_thresh_string ' Scale=' num2str(params.Scale)]};...
                     {[' # time pts kept=' num2str(size(DM,1)) ]}],'Color','w', 'FontSize', 8, 'interpreter', 'none');
                subplot(3,num_views+1,j.*4);
                params.fig_handle=gca;
                params.view = 'dorsal';
                PlotInterpSurfMesh(map,anatomyIn.MNIl,anatomyIn.MNIr,anatomyIn.info,params);
                title([{[run,' ', stim, ' GLM with ', params.type]};{['GVTD.Th=' GVTD_thresh_string ' Scale=' num2str(params.Scale)]};...
                     {[' # time pts kept=' num2str(size(DM,1)) ]}],'Color','w', 'FontSize', 8, 'interpreter', 'none');
            end
            
            
        case 'OBS_IM_ISS'
            %% Observation / Immitation
            
            % Sort betas by event type
            if strcmp(params.type, 'no_events')
                b0 = b(:,2:end);
                b_rest = b0(:, dataIn.info.paradigm.Pulse_1);
                b_immitation = b0(:, dataIn.info.paradigm.Pulse_2);
            end
            
            switch params.type
                case 'events'
                    % Get map
                    bL=(squeeze(b(:,3))-squeeze(b(:,2)));
                    bLsmooth=smooth_Amat(bL,dataIn.info.tissue.dim,params.flags.gsigma); %smoothing, must run this while map is still 1D, before Good_vox2vol
                    bLvol = Good_Vox2vol(bLsmooth,dataIn.info.tissue.dim);
                    map = affine3d_img(bLvol,dataIn.info.tissue.dim,anatomyIn.info,affine); %affine transform data to MNI space
                    if contains(params.flags.cap, 'LUMO') % Multiply by FOV, only for LUMO bc A was made on scaled mesh (may change this section in the future)
                        map = map.*FOV; 
                    end

                    % Visualization parameters
                    if length(params.scale_factor) > 1 % Make sure scale factor is not long bc this uses the same scale for all 3 views
                        disp(['Scale factor set improperly, please only set one scale factor for ' run ' runs'])
                        params.scale_factor = 3.5; % set to 3.5 because for these paradigms we use std of data instead of std of data to set colorscale
                    end
                    params.Scale=params.scale_factor*std(bL(bL(:)~=0)); %scale parameters, you might need to adjust these to see activations, normally set to 5*std(bL(bL(:)~=0)) for immitation
                    if params.Use_hardcode_scale == 1 %use hardcoded scale set by user if this flag is turned on
                            params.Scale = params.Scale_hardcoded(j);
                    end
                case 'no_events'
                    % Do T-test and get cohen's D
                    [~,~,~,d ]= TtestWelch(b_words, b_rest);
                    dvol=Good_Vox2vol(d,dataIn.info.tissue.dim);
                    map = affine3d_img(dvol,dataIn.info.tissue.dim,anatomyIn.info,affine); % Affine transform data to MNI space
                    if contains(params.flags.cap, 'LUMO')
                        map = map.*FOV; % Multiply by FOV, only for LUMO bc A was made on scaled mesh (may change this section in the future)
                    end

                    % Visualization Parameters
                    if length(params.scale_factor) > 1 % For these paradigms, scale factor should never be more than 1 value since we're only making one map
                        disp(['Scale factor set improperly, please only set one scale factor for ' run ' runs'])
                        params.scale_factor = 3.5; % Set to 3.5 bc we're using std now
                    end
                    params.Scale=params.scale_factor(j)*std(dvol(dvol(:)~=0)); %color scale parameters, you might want to adjust these  
                    if params.Use_hardcode_scale == 1 %use hardcoded scale set by user if this flag is turned on
                            params.Scale = params.Scale_hardcoded(j);
                    end
            end
                    
            
            % Visualize
            % This will make a single figure with all 3 views of the brain
            figure('Color','k','Position',[100,100,1200,800]);
            %lateral view 
            subplot(2,2,[1,2]);
            params.fig_handle=gca;
            params.view='lat'; 
            PlotInterpSurfMesh(map,anatomyIn.MNIl,anatomyIn.MNIr,anatomyIn.info,params);
            title([run ' GLM with ', params.type, ' GVTD.Th=' GVTD_thresh_string ' Scale=' num2str(params.Scale) ' Thresh(p,n)='...
                    num2str(params.Th.P) ',' num2str(params.Th.N) ' # time pts kept=' num2str(size(DM,1)) ],'Color','w');
            %dorsal view 
            subplot(2,2,[3]);
            params.fig_handle=gca;
            params.view='dorsal';
            PlotInterpSurfMesh(map,anatomyIn.MNIl,anatomyIn.MNIr,anatomyIn.info,params);

            %posterior view 
            subplot(2,2,[4]);
            params.fig_handle=gca;
            params.view='post';
            PlotInterpSurfMesh(map,anatomyIn.MNIl,anatomyIn.MNIr,anatomyIn.info,params);
            
            
        case 'SLP'
            %% Semantic Language Processing
            
            if strcmp(params.type,'no_events')
                % Sort betas by event type
                b0 = b(:,2:end);
                b_rest = b0(:, dataIn.info.paradigm.Pulse_1);
                b_norm = b0(:, dataIn.info.paradigm.Pulse_2);
                b_noisy = b0(:, dataIn.info.paradigm.Pulse_3);
            end
            
            % Make figure window, this will get filled in using the below loop
            figure('Color','k','Position',[100,100,1200,800]);
            type = {'Norm - Rest', 'Noisy - Rest', 'Norm - Noisy'};
            for j = 1:length(type)

                stim = type{j};
                
                if length(params.scale_factor) < 3 % Make sure scale factor is not long bc this uses the same scale for all 3 views
                    disp(['Only ' num2str(length(params.scale_factor)) ' scale factor given. Using the same scale factor for each map'])
                    params.scale_factor = repmat(params.scale_factor, [1,length(type)]); % set to 3.5 because for these paradigms we use std of data instead of std of data to set colorscale
                end
                
                switch params.type
                    case 'events'
                
                        % Make volume and set scale
                        switch stim
                            case 'Norm - Rest'
                                % Make maps of normal audio stimulus
                                bL=squeeze(b(:,3)-b(:,2));
                                bvol=Good_Vox2vol(bL,dataIn.info.tissue.dim);
                                params.Scale=params.scale_factor(j)*std(bvol(bvol(:)~=0)); % color scale parameter
                            case 'Noisy - Rest'
                                % Make maps of noisy audio stimulus
                                bR=squeeze(b(:,4)-b(:,2));
                                bvol=Good_Vox2vol(bR,dataIn.info.tissue.dim);
                                params.Scale=params.scale_factor(j)*std(bvol(bvol(:)~=0)); %color scale parameter  
                            case 'Norm - Noisy'
                                % Make maps of contrast (normal - noisy)
                                bC=squeeze(b(:,3)-b(:,4));
                                bvol=Good_Vox2vol(bC,dataIn.info.tissue.dim);
                                params.Scale=params.scale_factor(j)*std(bvol(bvol(:)~=0)); %color scale parameter 
                        end

                        % Affine transform data to MNI space
                        map=affine3d_img(bvol,dataIn.info.tissue.dim,anatomyIn.info,affine);
                    
                    case 'no_events'
                        % Make volume and set scale
                        switch stim
                            case 'Norm - Rest'
                                % Make maps of left stimulus
                                [~,~,~,d ]= TtestWelch(b_norm, b_rest);
                                dvol=Good_Vox2vol(d,dataIn.info.tissue.dim);
                                params.Scale=params.scale_factor(j)*std(dvol(dvol(:)~=0)); 
                            case 'Noisy - Rest'
                                % Make maps of right stimulus
                                [~,~,~,d] = TtestWelch(b_noisy, b_rest);
                                dvol=Good_Vox2vol(d,dataIn.info.tissue.dim);
                                params.Scale=params.scale_factor(j)*std(dvol(dvol(:)~=0));  
                            case 'Norm - Noisy'
                                % Make maps of contrast (left - right)
                                [~,~,~,d] = TtestWelch(b_norm, b_noisy);
                                dvol=Good_Vox2vol(d,dataIn.info.tissue.dim);
                                params.Scale=params.scale_factor(j)*std(dvol(dvol(:)~=0)); 
                        end
                        
                        % Affine transform data to MNI space
                        map=affine3d_img(dvol,dataIn.info.tissue.dim,anatomyIn.info,affine);
                end
                
                if params.Use_hardcode_scale == 1 %use hardcoded scale set by user if this flag is turned on
                    params.Scale = params.Scale_hardcoded(j);
                end
                
                %Multiply by FOV
                map = map.*FOV;

                % Plot each view
                views = {'lat', 'post', 'dorsal'};
                num_views = length(views);
                subplot(3,num_views+1,j*4-3:j*4-2);
                params.fig_handle=gca;
                params.view = 'lat';
                PlotInterpSurfMesh(map,anatomyIn.MNIl,anatomyIn.MNIr,anatomyIn.info,params);
                title([{[run,' ', stim, ' GLM with ', params.type]};{['GVTD.Th=' GVTD_thresh_string ' Scale=' num2str(params.Scale)]};...
                     {[' # time pts kept=' num2str(size(DM,1)) ]}],'Color','w', 'FontSize', 8);
                subplot(3,num_views+1,j.*4-1);
                params.fig_handle=gca;
                params.view = 'post';
                PlotInterpSurfMesh(map,anatomyIn.MNIl,anatomyIn.MNIr,anatomyIn.info,params);
                title([{[run,' ', stim, ' GLM with ', params.type]};{['GVTD.Th=' GVTD_thresh_string ' Scale=' num2str(params.Scale)]};...
                     {[' # time pts kept=' num2str(size(DM,1)) ]}],'Color','w', 'FontSize', 8);
                subplot(3,num_views+1,j.*4);
                params.fig_handle=gca;
                params.view = 'dorsal';
                PlotInterpSurfMesh(map,anatomyIn.MNIl,anatomyIn.MNIr,anatomyIn.info,params);
                title([{[run,' ', stim, ' GLM with ', params.type]};{['GVTD.Th=' GVTD_thresh_string ' Scale=' num2str(params.Scale)]};...
                     {[' # time pts kept=' num2str(size(DM,1)) ]}],'Color','w', 'FontSize', 8, 'interpreter', 'none');
            end
            
    end
    
    %% Save Map
    if params.save_map == 1
        if params.GVTD_censor == 1
            SaveGCF_to_PNG([datafile,'_thresholded_map_',params.flags.tag],params.kill);
        end
        if params.GVTD_censor == 0
            SaveGCF_to_PNG([datafile,'_no_threshold_map_',params.flags.tag],params.kill);
        end
    else
        %do nothing
    end
    
    
end
