%% NEURODOT VISUALIZATION TUTORIAL SCRIPT
% This script is designed to process and visualize a set of Motor finger tapping data collected on a LUMO system at WashU
% You will first process the data and visualize it at the measurement level
% Then you will perform image reconstruction and spectroscopy to visualize the data volumetrically and on a mesh surface
% The end of the script, lines 327-end provide a sandbox for playing around with NeuroDOT's many visualization parameters
% For any questions, please see Adam, Ari, or Emma! We will be walking around the room to help as we can :)


%% Load Measurement data
load('ND_LUMO_Data_Sample_MOT001.mat')

%% Crop data to match A Matrix' measurement list

% Load measurement list from A Matrix
temp = load(['A_LUMO_51tile_KKI_on_HD_Mesh0_r3dlt50mm.mat'],'info'); 

% Crop data
mA=cat(2,temp.info.pairs.Src,temp.info.pairs.Det,temp.info.pairs.WL);
mD=cat(2,info.pairs.Src,info.pairs.Det,info.pairs.WL);
keepInfo=ismember(mD,mA,'rows');
data=data(keepInfo,:);                                                    % Choose measurementps from data
info.pairs=temp.info.pairs;                                               % Grab relevant measurements from pairs with A
info.MEAS.Clipped = info.MEAS.Clipped(keepInfo);


%% Set flags for measurement level processing and visualization
flags.RsdMax = 35;                                                        % Maximum source-detector separation
flags.AbsDiff = 0.5;                                                      % For gradual drop-outs
flags.FloorTh = 1e-7;                                                     % For low light levels
flags.DtKeep = 1;                                                         % For sudden channel drop-outs
flags.LFO_GI = 1;                                                         % Only include good measurments in light falloff plot
    
    
%% General data QC with synchpts if present
Plot_RawData_Time_Traces_Overview(data, info, flags);                     % Time traces
info.DQ_metrics.WL = Plot_RawData_Metrics_II_DQC(data, info, flags);       % Raw data quality figs
info = Plot_RawData_Cap_DQC(data, info, flags);                           % Cap-relevant views


%% Pre-process data

% Pre-processing pipeline, but only the part before filtering data
[lmdata_b4_filt, info.MEAS.Phi_0] = logmean(data);                        % Logmean Light Levels
info_b4_filt  = FindGoodMeas(lmdata_b4_filt, info, 0.075, flags);        % Detect Noisy Channels
lmdata_b4_filt  = detrend_tts(lmdata_b4_filt );                           % Detrend Data
info_b4_filt.GVTD = CalcGVTD(lmdata_b4_filt(info_b4_filt.MEAS.GI & info_b4_filt.pairs.r2d<20,:));         % Calculate GVTD

% Full Pre-processing Pipeline
[lmdata, info.MEAS.Phi_0] = logmean(data);                                % Logmean Light Levels
info = FindGoodMeas(lmdata, info, 0.075, flags);                          % Detect Noisy Channels
lmdata = detrend_tts(lmdata);                                             % Detrend Data
lmdata = highpass(lmdata, .02, info.system.framerate);                    % High Pass Filter (0.02 Hz)
lmdata = lowpass(lmdata, 1, info.system.framerate);                       % Low Pass Filter 1 (1.0 Hz)
hem = gethem(lmdata, info, 'r2d', [1,15]);                                % Superficial Signal Regression
[lmdata, ~] = regcorr(lmdata, info, hem);
lmdata = lowpass(lmdata, 0.2, info.system.framerate);                     % Low Pass Filter 2 (0.2 Hz)
info.GVTD = CalcGVTD(lmdata(info.MEAS.GI & info.pairs.r2d<20,:));         % Calculate GVTD
[lmdata, info] = resample_tts(lmdata, info, 1, 1e-5);                     % 1 Hz Resampling (1 Hz)
info.GVTD_filt_rs = resample_tts(info.GVTD', info, 1, 1e-5,...
    info.system.init_framerate)';                                         % Resample GVTD to 1 Hz


%% Visualize Pre-filtered Data
% Measurements to include
keep_b4_filt = info_b4_filt.pairs.WL==2 & info_b4_filt.pairs.r2d < 35 & info_b4_filt.MEAS.GI; 

% Visualize
figure('Position',[100 100 550 780])
subplot(3,1,1); plot(lmdata_b4_filt(keep_b4_filt,:)'); 
set(gca,'XLimSpec','tight'), xlabel('Time (samples)'), 
ylabel('log(\phi/\phi_0)') 
m_b4_filt=max(max(abs(lmdata_b4_filt(keep_b4_filt,:))));
subplot(3,1,2); imagesc(lmdata_b4_filt(keep_b4_filt,:),[-1,1].*m_b4_filt); 
colorbar('Location','northoutside');
xlabel('Time (samples)');ylabel('Measurement #')
[ftmag_b4_filt,ftdomain_b4_filt] = fft_tts(squeeze(mean(lmdata_b4_filt(keep_b4_filt,:),1)),info_b4_filt.system.framerate); % Generate average spectrum
subplot(3,1,3); semilogx(ftdomain_b4_filt,ftmag_b4_filt);
xlabel('Frequency (Hz)');ylabel('|X(f)|');xlim([1e-3 1])

nlrGrayPlots_180818(lmdata_b4_filt,info_b4_filt);                         % Gray Plot with synch points


%% Visualize Filtered Data
% Measurements to include
keep = info.pairs.WL==2 & info.pairs.r2d < 35 & info.MEAS.GI;             % Measurements to include

% Visualize
figure('Position',[100 100 550 780])
subplot(3,1,1); plot(lmdata(keep,:)'); 
set(gca,'XLimSpec','tight'), xlabel('Time (samples)'), 
ylabel('log(\phi/\phi_0)') 
m=max(max(abs(lmdata(keep,:))));
subplot(3,1,2); imagesc(lmdata(keep,:),[-1,1].*m); 
colorbar('Location','northoutside');
xlabel('Time (samples)');ylabel('Measurement #')
[ftmag,ftdomain] = fft_tts(squeeze(lmdata(keep,:)),info.system.framerate); % Generate average spectrum
ftmag = rms(ftmag,1);
subplot(3,1,3); semilogx(ftdomain,ftmag);
xlabel('Frequency (Hz)');ylabel('|X(f)|');xlim([1e-3 1])

nlrGrayPlots_180818(lmdata,info);                                         % Gray Plot with synch points

% Compare 750nm and 850nm wavelengths
params.rlimits = [1,35];                                                  % Only look at data from meas with SD separation between 1mm and 35mm
Plot_TimeTrace_With_PowerSpectrum(lmdata,info,params);                    % Data, Grayplots organized by measurement #, Powerspectrum
GrayPlots_Rsd_by_Wavelength(lmdata,info);                                 % Grayplots organized by SD separation


%% Block Averaging the measurement data and view

% Initialization
LeftRight = {'Pulse_2', 'Pulse_3'};                                       % Pulse 2 = left, Pulse 3 = right
dt = 34;                                                                  % Total time for task + rest in seconds

% Loop to block average measurement data and visualize
for j = 1:length(LeftRight)
    badata = BlockAverage(lmdata, info.paradigm.synchpts(info.paradigm.(LeftRight{j})), dt);

    badata=bsxfun(@minus,badata,mean(badata(:,1:4),2));

    figure('Position',[100 100 550 780])
    subplot(2,1,1); plot(badata(keep,:)'); 
    set(gca,'XLimSpec','tight'), xlabel('Time (samples)'), 
    ylabel('log(\phi/\phi_0)') 
    m=max(max(abs(badata(keep,:))));
    subplot(2,1,2); imagesc(badata(keep,:),[-1,1].*m); 
    colorbar('Location','northoutside');
    xlabel('Time (samples)');
    ylabel('Measurement #')
end


%% RECONSTRUCTION PIPELINE
if ~exist('A', 'var')        % In case running by hand or re-running script
    A=load(['A_LUMO_51tile_KKI_on_HD_Mesh0_r3dlt50mm.mat'],'info','A');     
end
Nvox=size(A.A,2);
Nt=size(lmdata,2);
Nwl = length(unique(info.pairs.WL));
cortex_mu_a=zeros(Nvox,Nt,2);
for j = 1:Nwl
    keep = (A.info.pairs.WL == j) & (A.info.pairs.r2d <= 35) & info.MEAS.GI;
    disp('> Inverting A')                
    iA = Tikhonov_invert_Amat(A.A(keep, :), 0.1, 0.1);                    % Invert A-Matrix
    disp('> Smoothing iA')
    iA = smooth_Amat(iA, A.info.tissue.dim, 3);                           % Smooth Inverted A-Matrix      
    cortex_mu_a(:, :, j) = reconstruct_img(lmdata(keep, :), iA);          % Reconstruct Image Volume
end


%% Spectroscopy
if ~exist('E', 'var'),load('E.mat'),end
cortex_Hb = spectroscopy_img(cortex_mu_a, E);
cortex_HbO = cortex_Hb(:, :, 1);                                          % Oxy
cortex_HbR = cortex_Hb(:, :, 2);                                          % Deoxy
cortex_HbT = cortex_HbO + cortex_HbR;                                     % Total


%% Load files for Visualizing Reconstructed Data

[MNI_111,info111]=LoadVolumetricData('Segmented_MNI152nl_on_MNI111_nifti',[],'nii'); % MRI with 1mm voxellation - MNI T1
[MNI_333,info333]=LoadVolumetricData('Segmented_MNI152nl_on_333_nifti',[],'nii');    % MRI with 3mm voxellation - MNI T1
load(['hrf_DOT3.mat']);                                                              % HRF for GLM processing
load(['MNI164k_big.mat'])                                                            % Brain Mesh
% load('LUMO_51tile_affine_matrix_MNI_to_subj_vol.mat');                               % Affine matrix from A-matrix volume space to MNI 1mm volume space
LUMO_to_MNI = A.info.tissue.affine;                                                  % Affine matrix from A-matrix volume space to MNI 1mm volume space
load('LUMO_51tile_KKI_r3dlt50_FOV_on_MNI.mat')                                                   % LUMO cap field of view


%% Visualize Block-Averaged Data Volumetrically

% Threshold: motion above whatever number is here is removed
params.GVTD_th = 0.001; 

% LEFT HAND: Block Average, then affine transform data to tissue space
tkeep = info.GVTD < params.GVTD_th; 
[BA_out_L,~,BT_out] = BlockAverage(squeeze(cortex_Hb(:,:,1)),...
    info.paradigm.synchpts(info.paradigm.Pulse_2), dt, tkeep); 
BA_out_L=bsxfun(@minus,BA_out_L,mean(BA_out_L(:,1:4),2));
BAvol_L=Good_Vox2vol(BA_out_L,A.info.tissue.dim);                         % Turn the block average matrix into a volume
BA_MNI_L=affine3d_img(squeeze(BAvol_L(:,:,:,:)),...                       % Affine transform data to the tissue
    A.info.tissue.dim,info333,LUMO_to_MNI);
params.Scale=0.8*max(BA_MNI_L(BA_MNI_L(:)~=0));                           % Set color scale based off of data
params.Th.P = 0;                                                          % Color Scale Threshold, try changing to 0.25*params.Scale;
params.Th.N = -params.Th.P;

%Block averaging was done in the above section, so all we need to do is run PlotSlicesTimeTrace
PlotSlicesTimeTrace(MNI_333,info333,params,BA_MNI_L,info)


% RIGHT HAND: Block Average, then affine transform data to tissue space
tkeep = info.GVTD < params.GVTD_th; 
[BA_out_R,~,BT_out] = BlockAverage(squeeze(cortex_Hb(:,:,1)),...
    info.paradigm.synchpts(info.paradigm.Pulse_3), dt, tkeep); 
BA_out_R=bsxfun(@minus,BA_out_R,mean(BA_out_R(:,1:4),2));
BAvol_R=Good_Vox2vol(BA_out_R,A.info.tissue.dim);                         % Turn the block average matrix into a volume
BA_MNI_R=affine3d_img(squeeze(BAvol_R(:,:,:,:)),...                       % Affine transform data to the tissue
    A.info.tissue.dim,info333,LUMO_to_MNI);
params.Scale=0.8*max(BA_MNI_R(BA_MNI_R(:)~=0));                           % Set color scale based off of data
params.Th.P = 0;                                                          % Color Scale Threshold 0.25*params.Scale;
params.Th.N = -params.Th.P;

%Block averaging was done in the above section, so all we need to do is run PlotSlicesTimeTrace
PlotSlicesTimeTrace(MNI_333,info333,params,BA_MNI_R,info)


%% Visualize Block-Averaged Data on the Surface

% Left Hand Side
tp = 14;                                                                  % Select time point
HbO_atlas_L = affine3d_img(BAvol_L,A.info.tissue.dim,info333,LUMO_to_MNI);% Affine transform data to the tissue
tp_Eg_atlas_L=squeeze(HbO_atlas_L(:,:,:,tp));
Params.Scale=0.3*max(abs(BAvol_L(:)));                                    % Color Scale
Params.Th.P=0;                                                            % Color Scale Threhsolds
Params.Th.N=-Params.Th.P;
Params.view='dorsal'; % Posterior view
Params.ctx='std'; % Standard pial cortical view
PlotInterpSurfMesh(tp_Eg_atlas_L, MNIl,MNIr, info333, Params); 
title(['Left Hand Finger Tapping, Block Average, Time Point ', num2str(tp)], 'Color', 'w')

% Right Hand Side
tp = 13;                                                                  % Select time point
HbO_atlas_R = affine3d_img(BAvol_R,A.info.tissue.dim,info333,LUMO_to_MNI);% Affine transform data to the tissue
tp_Eg_atlas_L=squeeze(HbO_atlas_R(:,:,:,tp));
Params.Scale=0.3*max(abs(BAvol_L(:)));                                    % Color Scale
Params.Th.P=0;                                                            % Color Scale Threhsolds
Params.Th.N=-Params.Th.P;
Params.view='dorsal';                                                     % Posterior view
Params.ctx='std';                                                         % Standard pial cortical view
PlotInterpSurfMesh(tp_Eg_atlas_L, MNIl,MNIr, info333, Params); 
title(['Right Hand Finger Tapping, Block Average, Time Point ', num2str(tp)], 'Color', 'w')


%% Visualize Non Block-Averaged Data Volumetrically

HbOvol = Good_Vox2vol(cortex_HbO,A.info.tissue.dim);                      % Make Volume

% Spatial Normalization - Note: this section takes ~1-3 minutes to run
HbO_MNI111 = affine3d_img(squeeze(HbOvol(:,:,:,:)),...                    % Affine transform data to the tissue - takes a while
    A.info.tissue.dim,info111,LUMO_to_MNI);
HbO_MNI333 = affine3d_img(squeeze(HbO_MNI111(:,:,:,:)),...                % Affine transform data to the tissue - takes a while
    info111,info333,eye(4), 'nearest');

% Visualize on 1mm - NOTE: THIS WILL TAKE > 10min TO RUN, DO NOT UNCOMMENT RUN THIS SECTION UNLESS YOU WANT TO WAIT A WHILE
% Params.Scale = max(HbO_MNI111(HbO_MNI111(:)~=0));
% Params.Th.P = 0;
% Params.Th.N = -Params.Th.P;
% PlotSlicesTimeTrace(MNI_111,info111,Params,HbO_MNI111,info)               

% Visualize on 3mm - Note: This section doesn't take as much time to run
Params.Scale = 0.4*max(HbO_MNI333(HbO_MNI333(:)~=0));
Params.Th.P = 0;
Params.Th.N = -Params.Th.P;
PlotSlicesTimeTrace(MNI_333,info333,Params,HbO_MNI333,info)


%% View GLM Maps - run this section once, then change the GVTD threshold in line 356 and run it again

% GVTD params
% Most important thing here is the threshold, set in line 233
params_GLM = struct;
params_GLM.GVTD_censor = 1;                                        % Do you want to remove motion? 1 = yes, 0 = no
params_GLM.GVTD_th = 0.001;                                        % Motion above whatever number is here is removed (reference GVTD plot)
params_GLM.DoFilter = 0;
params_GLM.type = 'events';                                        % Try this again with 'no_events' as the type

% % Remove initial fixation - Optional, but useful!
info_no_IF = info;
info_OG = info;                                                    % Backup
info_no_IF.paradigm.synchpts = info_no_IF.paradigm.synchpts(2:end);
info_no_IF.paradigm.init_synchpts = info_no_IF.paradigm.init_synchpts(2:end);
info_no_IF.paradigm.synchtype = info_no_IF.paradigm.synchtype(2:end);
info_no_IF.paradigm.synchtimes = info_no_IF.paradigm.synchtimes(2:end);
info_no_IF.paradigm.Pulse_1 = info_no_IF.paradigm.Pulse_1(2:end);
info_no_IF.paradigm.Pulse_1 = info_no_IF.paradigm.Pulse_1 - 1;
info_no_IF.paradigm.Pulse_2 = info_no_IF.paradigm.Pulse_2 - 1;
info_no_IF.paradigm.Pulse_3 = info_no_IF.paradigm.Pulse_3 - 1;
info = info_no_IF;

% Run GLM
[b,e,DM,EDM]=GLM_181206(squeeze(cortex_Hb(:,:,1)),hrf,info,params_GLM); % GLM code - rerun after setting new threshold

% Make figure window, this will get filled in using the below loop
figure('Color','k','Position',[100,100,800,1000]);

% Map making loop
run = 'MOT001';
type = {'left', 'right', 'contrast'};
for j = 1:length(type)
    
    stim = type{j};
    
    switch stim
        case 'left'
            % Make maps of left stimulus
            bL=squeeze(b(:,3)-b(:,2));
            bvol=Good_Vox2vol(bL,A.info.tissue.dim);
            params_GLM.Scale=5.5*std(bvol(bvol(:)~=0));            % Color Scale
            
        case 'right'
            % Make maps of right stimulus
            bR=squeeze(b(:,4)-b(:,2));
            bvol=Good_Vox2vol(bR,A.info.tissue.dim);
            params_GLM.Scale=5.5*std(bvol(bvol(:)~=0));            % Color Scale
            
        case 'contrast'
            % Make maps of contrast (left - right)
            bC=squeeze(b(:,3)-b(:,4));
            bvol=Good_Vox2vol(bC,A.info.tissue.dim);
            params_GLM.Scale=5.5*std(bvol(bvol(:)~=0));            % Color Scale
    end

    % Affine transform data to MNI space
    fooMNI=affine3d_img(bvol,A.info.tissue.dim,info111,LUMO_to_MNI);      % Spatial Normalization

    %Multiply by FOV
    fooMNI = fooMNI.*FOV;

    % Visualization Parameters 
    params_GLM.Th.P = 0;                                           % Color scale thresholds
    params_GLM.Th.N = -params_GLM.Th.P;
    params_GLM.view = 'dorsal';                                    % View of brain
    params_GLM.ctx = 'std';                                        % Cortex Inflation
    params_GLM.OL = 0;                                             % Overlap

    %Visualize
    % This will make a single figure with all 3 pulses (stim)
    subplot(3,1,[j]);
    params_GLM.fig_handle=gca;
    PlotInterpSurfMesh(fooMNI,MNIl,MNIr,info111,params_GLM);
    title([run ' GLM GVTD.Th=' num2str(params_GLM.GVTD_th) ' Scale=' num2str(params_GLM.Scale)...
        ' Thresh(p,n)=' num2str(params_GLM.Th.P) ',' num2str(params_GLM.Th.N) ' # time pts kept=' num2str(size(DM,1)) ],'Color','w', 'FontSize', 8);
end


%% Play Around with Visualization Parameters - Volumetric Data
close all

% Setup
Left_tp = 14;
Left_data = BA_MNI_L(:,:,:,Left_tp);
Right_tp = 13;
Right_data = BA_MNI_R(:,:,:,Right_tp);

% Choose which data you want to visualize (uncomment the line you want to use)
data_V = Left_data;
% dataV = Right_data;


% Explore PlotSlices - in depth
%Anatomy Only
PlotSlices(MNI_111,info111)

% Params 1 - Defaults
Params1.Th.P = 0;
Params1.Th.N = - Params1.Th.P;
PlotSlices(MNI_111,info111,Params1,data_V)

% Params 2 - Positive Definite
Params2 = Params1;
Params2.PD = 1; 
PlotSlices(MNI_111,info111,Params2,data_V)

% Params 3 - Setting Different Thresholds
Params3.Th.P = 1e-9;
Params3.Th.N = -4e-8;
PlotSlices(MNI_111,info111,Params3,data_V)

% Params 4 - Color Scale
Params4 = Params1;
Params4.Scale = 2*max(abs(data_V(:)));
PlotSlices(MNI_111,info111,Params4,data_V)

% Params 5 - Colormap
Params5 = Params1;
Params5.Cmap = 'cool';
PlotSlices(MNI_111,info111,Params5,data_V)

% Params 6 - Different Colormaps for Positive and Negative
Params6 = Params1;
Params6.Cmap.P = 'cool';
Params6.Cmap.N = 'summer';
PlotSlices(MNI_111,info111,Params6,data_V)

% Params 7 - Colorbar Ticks
Params7 = Params1;
Params7.cbmode = 1;
Params7.cbticks = [0,0.25,0.5,0.75,1];
PlotSlices(MNI_111,info111,Params7,data_V)

% Params 8 - Colorbar Labels
Params8 = Params1;
Params8.cbmode = 1;
Params8.cblabels = {'Hello', 'Goodbye'};
PlotSlices(MNI_111,info111,Params8,data_V)

% Params 9 - Putting it All together
Params9.cbmode = 1;
Params9.cbticks = [0,0.25,0.5,0.75,1];
Params9.cblabels = {'1', '2', '3', '4', '5'};
Params9.Cmap.P = 'cool';
Params9.Cmap.N = 'summer';
Params9.Scale = 0.8*max(abs(data_V(:)));
Params9.Th.P = 1e-8;
Params9.Th.N = -4e-7;
PlotSlices(MNI_111,info111,Params9,data_V)

%% Play Around with Visualization Parameters - Reconstructed Data - GLM Contrast Map
close all

% Choose Left or Right Data to use here
% Left Data
tp = 14;                                                                  
HbO_atlas_L = affine3d_img(BAvol_L,A.info.tissue.dim,info111,LUMO_to_MNI);
data_surf=squeeze(HbO_atlas_L(:,:,:,tp));
% Right Data
% tp = 13;                                                                  
% HbO_atlas_R = affine3d_img(BAvol_R,A.info.tissue.dim,info111,LUMO_to_MNI);
% data_surf = squeeze(HbO_atlas_R(:,:,:,tp));


% Set Initial visualization params
Params.Scale=0.3*max(abs(data_surf(:)));                                  % Color Scale
Params.Th.P=0;                                                            % Color Scale Threhsolds
Params.Th.N=-Params.Th.P;

% Params as set in GLM section 
PlotInterpSurfMesh(data_surf,MNIl,MNIr,infoMNI,Params);

% Change the View
% Posterior view
Params.view = 'post';
PlotInterpSurfMesh(data_surf,MNIl,MNIr,infoMNI,Params);
% Lateral view
Params.view = 'lat';
PlotInterpSurfMesh(data_surf,MNIl,MNIr,infoMNI,Params);

% Change the Cortex Inflation
Params.view='dorsal';                                                 % Reset to dorsal view before visualizing
% Standard pial cortical view
Params.ctx='std'; 
PlotInterpSurfMesh(data_surf,MNIl,MNIr,infoMNI,Params);
% Inflated pial cortical view
Params.ctx='inf'; 
PlotInterpSurfMesh(data_surf,MNIl,MNIr,infoMNI,Params);
% Very Inflated pial cortical view
Params.ctx='vinf';
PlotInterpSurfMesh(data_surf,MNIl,MNIr,infoMNI,Params);

% Change the Color Scale Thresholds
Params.ctx='std';                                                     % Reset to standard inflation before visualizing
% Set to 0
Params.Th.P = 0;
Params.Th.N = -Params.Th.P;
PlotInterpSurfMesh(data_surf,MNIl,MNIr,infoMNI,Params);
% Set a bit higher
Params.Th.P = 1e-6;
Params.Th.N = -Params.Th.P;
PlotInterpSurfMesh(data_surf,MNIl,MNIr,infoMNI,Params);
% Set even higher
Params.Th.P = 4e-5;
Params.Th.N = -Params.Th.P;
PlotInterpSurfMesh(data_surf,MNIl,MNIr,infoMNI,Params);
% Set very high
Params.Th.P = 1e-3;
Params.Th.N = -Params.Th.P;
PlotInterpSurfMesh(data_surf,MNIl,MNIr,infoMNI,Params);

% Change the Color Scale
Params.Th.P = 0;
Params.Th.N = -Params.Th.P;                                       % Reset thresholds first
% Set to 5.5*std
Params.Scale = 3.5*std(data_surf(data_surf(:)~=0));
PlotInterpSurfMesh(data_surf,MNIl,MNIr,infoMNI,Params);
% Set lower
Params.Scale = 1.5*std(data_surf(data_surf(:)~=0));
PlotInterpSurfMesh(data_surf,MNIl,MNIr,infoMNI,Params);
% Set higher
Params.Scale = 6.5*std(data_surf(data_surf(:)~=0));
PlotInterpSurfMesh(data_surf,MNIl,MNIr,infoMNI,Params);

% Colormaps
% reset scale and thresholds first
Params.Scale = 3.5*std(data_surf(data_surf(:)~=0));
Params.Th.P = 0;
Params.Th.N = -Params.Th.P;
% Same colormap for P/N
Params.Cmap = 'cool';
PlotInterpSurfMesh(data_surf,MNIl,MNIr,infoMNI,Params);
% Differnt colormaps for P/N
Params = rmfield(Params, 'Cmap');                                 % Remove Cmap field so we can replace it with a struct containing both positive and negative colormaps
Params.Cmap.P = 'jet';
Params.Cmap.N = 'summer';
PlotInterpSurfMesh(data_surf,MNIl,MNIr,infoMNI,Params);
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%