
%% Section 1: Load Data

% Choose dataset
% dataset = 'ND_LUMO_Data_Sample_MOT001.mat';                               % Motor finger tapping
% dataset = 'ND_LUMO_Data_Sample_AC001.mat';                                % Alternating checkerboards
% dataset = 'ND_LUMO_Data_Sample_HW001.mat';                                % Hearing words
% dataset = 'ND_LUMO_Data_Sample_GV001.mat';                                % Generating verbs
dataset = 'ND_LUMO_Data_Sample_CCW001.mat';                               % Retinotopy - counter clockwise checkerboard

% Set flags for preprocessing and image reconstruction
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
flags.bthresh = 0.075;                                                    % Good measurement threshold
flags.omega_hp1 = 0.02;                                                   % High pass filter
flags.omega_lp1 = 1;                                                      % First low pass filter
flags.omega_lp2 = 0.2;                                                    % Second low pass filter 
flags.lambda1 = 0.1;                                                      % Lambda 1 controls tradoff between image resolution and noise
flags.lambda2 = 0.1;                                                      % Lambda 2 controls spatially variant regularization
flags.gsigma = 3;                                                         % Spatial Smoothing
flags.rsd_max = 35;                                                       % Maximum S-D separationt to include in reconstruction


%% Section 2: Preprocessing and Image Reconstruction

% Load Measurement data
load(dataset)    

% Crop data to match A Matrix' measurement list
temp = load(['A_LUMO_51_tile_on_particpant_mesh_r3dlt50mm.mat'],'info');  % Load measurement list from A Matrix
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);


% 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


% Full Pre-processing Pipeline
[lmdata, info.MEAS.Phi_0] = logmean(data);                                % Logmean Light Levels
info = FindGoodMeas(lmdata, info, flags.bthresh, flags);                  % Detect Noisy Channels
lmdata = detrend_tts(lmdata);                                             % Detrend Data
lmdata = highpass(lmdata, flags.omega_hp1, info.system.framerate);        % High Pass Filter (0.02 Hz)
lmdata = lowpass(lmdata, flags.omega_lp1, 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, flags.omega_lp2, 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


% Reconstruction
if ~exist('A', 'var')                                                     % In case running by hand or re-running script
    A=load(['A_LUMO_51_tile_on_particpant_mesh_r3dlt50mm'],'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 <= flags.rsd_max) & info.MEAS.GI;
    disp('> Inverting A')                
    iA = Tikhonov_invert_Amat(A.A(keep, :), flags.lambda1, flags.lambda2);% Invert A-Matrix
    disp('> Smoothing iA')
    iA = smooth_Amat(iA, A.info.tissue.dim, flags.gsigma);                % 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



%% Section 3: Load files used in visualization
 
[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
LUMO_to_MNI = A.info.tissue.affine;                                                  % Affine matrix from A-matrix volume space to MNI volume space
load('FOV_LUMO_51_tile_on_particpant_mesh_r3dlt50mm')                                % FOV for LUMO 51 tile cap
load('GoodVox3mm.mat')                                                               % GoodVox from A matrix in 3mm space (needed for CCW surface visualization)


%% Section 4: Visualize Block Averaged Data - MOT and AC

% Parameters for visualization
dt = 34;                                                                  % Block length
tp = 14;                                                                  % Select time point for surface vis
params.GVTD_th = 0.001;                                                   % GVTD Threshold: motion above whatever number is here is removed
tkeep = info.GVTD < params.GVTD_th;                                       % Keep time points with low motion
scale_factor_vol = 0.8;                                                   % Color scale for volumetric visualizations - multiplied by max(data)
scale_factor_surf = 0.3;                                                  % Color scale for surface visualizations - multiplied by max(data)
threshold_factor = 0;                                                     % Color thresholds for visualizations - try setting this to 0.25 - multiplied by color scale
Params.view='dorsal';                                                     % Cortex view - 'dorsal' for MOT, 'post' for AC
Params.ctx='std';                                                         % Cortex inflation


% Volumetric Visualizations (LEFT SIDE)
[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=scale_factor_vol*max(BA_MNI_L(BA_MNI_L(:)~=0));              % Set color scale based off of data
params.Th.P = threshold_factor*params.Scale;                              % Color Scale Threshold
params.Th.N = -params.Th.P;
PlotSlicesTimeTrace(MNI_333,info333,params,BA_MNI_L,info)                 % Visualize


% 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=scale_factor_vol*max(BA_MNI_L(BA_MNI_L(:)~=0));              % Set color scale based off of data
params.Th.P = threshold_factor*params.Scale;                              % Color Scale Threshold
params.Th.N = -params.Th.P;
PlotSlicesTimeTrace(MNI_333,info333,params,BA_MNI_R,info)                 % Visualize


% Surface Map - Left Hand Side
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=scale_factor_surf*max(abs(BAvol_L(:)));                      % Color Scale
Params.Th.P=0;                                                            % Color Scale Threhsolds
Params.Th.N=-Params.Th.P;
PlotInterpSurfMesh(tp_Eg_atlas_L, MNIl,MNIr, info333, Params); 
title(['Left, Block Average, Time Point ', num2str(tp)], 'Color', 'w')


% Surface Map - Right Hand Side                                                                
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=scale_factor_surf*max(abs(BAvol_L(:)));                      % Color Scale
Params.Th.P=0;                                                            % Color Scale Threhsolds
Params.Th.N=-Params.Th.P;
PlotInterpSurfMesh(tp_Eg_atlas_L, MNIl,MNIr, info333, Params); 
title(['Right, Block Average, Time Point ', num2str(tp)], 'Color', 'w')



%% Section 5: Visualize Block Averaged Data - HW and GV

% Parameters for visualization
dt = 30;                                                                  % Block length
tp = 12;                                                                  % Select time point for surface vis
params.GVTD_th = 0.001;                                                   % GVTD Threshold: motion above whatever number is here is removed
tkeep = info.GVTD < params.GVTD_th;                                       % Keep time points with low motion
scale_factor_vol = 0.5;                                                   % Color scale for volumetric visualizations - multiplied by max(data)
scale_factor_surf = 0.3;                                                  % Color scale for surface visualizations - multiplied by max(data)
threshold_factor = 0;                                                     % Color thresholds for visualizations - try setting this to 0.25 - multiplied by color scale
Params.view='lat';                                                        % Cortex view
Params.ctx='std';                                                         % Cortex inflation


% Volumetric Visualizations 
[BA_out,~,BT_out] = BlockAverage(squeeze(cortex_Hb(:,:,1)),...
    info.paradigm.synchpts(info.paradigm.Pulse_2), dt, tkeep); 
BA_out=bsxfun(@minus,BA_out,mean(BA_out(:,1:4),2));
BAvol=Good_Vox2vol(BA_out,A.info.tissue.dim);                             % Turn the block average matrix into a volume
BA_MNI=affine3d_img(squeeze(BAvol(:,:,:,:)),...                           % Affine transform data to the tissue
    A.info.tissue.dim,info333,LUMO_to_MNI);
params.Scale=scale_factor_vol*max(BA_MNI(BA_MNI(:)~=0));                  % Set color scale based off of data
params.Th.P = threshold_factor*params.Scale;                              % Color Scale Threshold
params.Th.N = -params.Th.P;
PlotSlicesTimeTrace(MNI_333,info333,params,BA_MNI,info)                   % Visualize


% Surface Map
HbO_atlas = affine3d_img(BAvol,A.info.tissue.dim,info333,LUMO_to_MNI);    % Affine transform data to the tissue
tp_Eg_atlas=squeeze(HbO_atlas(:,:,:,tp));
Params.Scale=scale_factor_surf*max(abs(BAvol(:)));                        % Color Scale
Params.Th.P=0;                                                            % Color Scale Threhsolds
Params.Th.N=-Params.Th.P;
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, info333, Params); 
title(['Block Average, Time Point ', num2str(tp)], 'Color', 'w')



%% Section 6: Visualize Block Averaged Data - CCW

% Parameters for visualization
tp = 13;                                                                  % Time point to visualize on surface
dt = 36;                                                                  % Block length
params.GVTD_th = 0.001;                                                   % GVTD Threshold: motion above whatever number is here is removed
tkeep = info.GVTD < params.GVTD_th;                                       % Keep time points with low motion
scale_factor_vol = 0.8;                                                   % Color scale for volumetric visualizations - multiplied by max(data)
scale_factor_surf = 0.3;                                                  % Color scale for surface visualizations - multiplied by max(data)
threshold_factor = 0;                                                     % Color thresholds for visualizations - try setting this to 0.25 - multiplied by color scale
Params.view='post';                                                       % Cortex view
Params.ctx='std';                                                         % Cortex inflation  


% Volumetric Visualizations 
[BA_out,~,BT_out] = BlockAverage(squeeze(cortex_Hb(:,:,1)),...
    info.paradigm.synchpts(info.paradigm.Pulse_2), dt, tkeep); 
BA_out=bsxfun(@minus,BA_out,mean(BA_out(:,1:4),2));
BAvol=Good_Vox2vol(BA_out,A.info.tissue.dim);                             % Turn the block average matrix into a volume
BA_MNI=affine3d_img(squeeze(BAvol(:,:,:,:)),...                           % Affine transform data to the tissue
    A.info.tissue.dim,info333,LUMO_to_MNI);
params.Scale=scale_factor_vol*max(BA_MNI(BA_MNI(:)~=0));                  % Set color scale based off of data
params.Th.P = threshold_factor*params.Scale;                              % Color Scale Threshold
params.Th.N = -params.Th.P;
PlotSlicesTimeTrace(MNI_333,info333,params,BA_MNI,info)                   % Visualize


% Surface Visualization in quadrants (this is some expert level visualization)
dt = round(mean(diff(info.paradigm.synchpts(info.paradigm.Pulse_2))));
dtq = round(dt/4);
ti = 2;
tf = 4;
mask = MNI_333 < 4 & MNI_333 > 0;

% Visualization parameters
pSba.Scale = 1; pSba.Th.P = 0; pSba.Th.N = -pSba.Th.P; pSba.PD = 0;
pSba.Cmap = 'jet'; pSba.view = 'post'; pSba.ctx = 'std';

% Make the figures
figure('Color','k','Units','Normalized','Position',[0.05,0.5,0.6,0.6])
for j=1:4
    subplot(2,2,j)
    BA_mean_tp = mean(BA_out(:,(((j-1).*dtq)+[ti:tf])),2);
    BAvol_mean_tp=Good_Vox2vol(BA_mean_tp,A.info.tissue.dim);                             % Turn the block average matrix into a volume
    foo=affine3d_img(squeeze(BAvol_mean_tp(:,:,:,:)),...                           % Affine transform data to the tissue
        A.info.tissue.dim,info333,LUMO_to_MNI);
    foo = foo./max(foo(:).*mask(:));
    pSba.fig_handle = gca;
    PlotInterpSurfMesh(foo, MNIl,MNIr, info333, pSba);
    title(['\mu HbO t=',num2str(((j-1).*dtq)+ti),'-',...
        num2str(((j-1).*dtq)+tf)],'Color','w')
end


%% Section 7: Visualize GLM Data (all data except for CCW)

% What Paradigm are you making maps for
run_type = 'GV001';                                                      % Run name goes here, options: 'MOT001', 'AC001', 'HW001', 'GV001', 'CCW001' 


% Flags for visualization function
flags.subj = 'NeuroDOT_Example_Subject';                                  % Subject ID            
flags.dateA = '260315';                                                   % Scan Date
flags.cap = 'LUMO';                                                       % Flag for LUMO


% Visualization parameters
params.GVTD_censor = 1;                                                   % Do you want to remove motion? 1 = yes, 0 = no
params.GVTD_th = 0.001;                                                   % GVTD threshold, start at 0.007, then lower the threshold using the histogram to crop out noisy data
params.type = 'events';                                                   % Type of GLM you want to run - choose either 'events' or 'no_events'
params.scale_factor = 4.5; % *std(data)                                   % Color scale parameter - play around with this to make activation stand out more/less
params.Th.P = 0;                                                          % Color scale thresholds
params.Th.N = -params.Th.P;
params.DoFilter=  1;                                                      % Filter data during GLM? 1 = yes, 0 = no
params.ctx = 'std';                                                       % Cortex inflation parameter - options: 'std', 'inf', 'vinf', 'flat'
params.save_map = 0;                                                      % Do you want to save the figure? 1 = yes, 0 = no
params.kill = 0;                                                          % Do you want to close the figure after its saved? 1 = yes 0 = no;  
params.flags = flags;                                                     % Put flags in params for visualization function to read


% Get cortexHB_input and brain mesh input
dataIn.cortex_Hb = cortex_Hb;
dataIn.info = info;
dataIn.info.tissue = A.info.tissue;
anatomyIn.MNIl = MNIl;
anatomyIn.MNIr = MNIr;
anatomyIn.info = info333;


% Make the map
[map, b, e, DM, EDM] = make_GLM_maps_260313_AS(dataIn, anatomyIn, hrf, LUMO_to_MNI, FOV, params, run_type); 



%% Make Phase Encoded Retinotopy Maps (CCW only)

% PER maps for CCW
info.tissue = A.info.tissue;
MNI_dim = affine3d_img(MNI_333,info333,A.info.tissue.dim,inv(A.info.tissue.affine),'nearest'); % Transform to DOT volume space
View_Data_QA_Quick_PER_maps_AS(cortex_HbO, info, MNI_dim, info333, MNIl,MNIr,'ND_Workshop_Glasgow_CCW_data_sample',0);