%% sFNIRS presentation %%
% Section 1: Generating a Light Model
% Section 2: Full Data Processing and Data Visualization


%% Section 1: Generating a Light Model %%
%% GENERATING A LIGHT MODEL OVERVIEW
% 1. start with Mask and Pad
% 2. make a Low Density mesh of the mask
% 3. use AlignMe to place the Pad on the LD mesh
% 4. remove unneeded components of the mask
% 5. make a HD mesh of the scooped mask
% 6. run AlignMe on the LD-aligned Pad and the HD mesh
% 7. update the info file and name it specifically for this light model
% 8. prepare mesh for nirfast
% 9. set flags and run makeAnirfast
% 10. package and save A with info etc
% 11. view single measurement, Gsd, and flat field reconstruction
% 12. view LD mesh, HD mesh, array, and cortex together


%% Pathing and Parameters
%Set Output Directory
%Change the path on the following line to whatever output directory you want to use
outputdir = '/Users/Ari/Desktop/Matlab/sfnirs_test1'; % put path to directory in between quotes
cd(outputdir)

%Parameter Initialization
p.Cmap='jet'; p.Scale=5; p.Th.P=0; p.Th.N=-p.Th.P; p.PD=1; p.BG=[0,0,0];
pM.orientation='coord'; pM.Cmap.P='gray'; pM.EdgeColor='none';
pS = pM; pS = rmfield(pS, 'orientation'); % make copy of params, to be used with PlotSlices


%% Load a Segmented Volume and generate head mesh
[mask,infoT1]=LoadVolumetricData(['Segmented_MNI152nl_on_MNI111_nifti'],[],'nii');
PlotSlices(mask,infoT1,p)       % Visualize the segmented mask: 
        % Note, PlotSlices is an interactive plot. 
        % Hit q or the middle mouse button to quit


%% Generate low density head mesh
% Parameters for generating your mesh
meshname=['LD_Mesh'];      % Provide a name for your mesh name here
param.facet_distance=5;    % Node position error tolerance at boundary, changed to 5 from .75 12/2/19
param.facet_size=3;        % boundary element size parameter
param.cell_size=5;         % Volume element size parameter
param.info=infoT1;
param.Offset=[0,0,0];
param.r0=5;                % nodes outside of mask must be set to scalp==5;
param.CheckMeshQuality=0;
param.Mode=1;              % make simple mesh with no region labels
tic;meshLD=NirfastMesh_Region(mask,meshname,param);toc

PlotMeshSurface(meshLD,pM) %visualize in index (non-negative) space
view([135,30]) %view the mesh where the nose of the head is pointing down and to the right
%this will give you a good view of where the center of the mesh is

%%%% IF you get an error, this is due to NIRFAST using an old version of
%%%% mode. Go into NIRFAST/toolbox, change mode.m to modeOLD.m. Then the
%%%% correct version of mode will be used and this visualization will work.

% Put coordinates back in true space
meshLD.nodes=change_space_coords(meshLD.nodes,infoT1,'coord');
PlotMeshSurface(meshLD,pM) %visualize in coordinate space
view([135,30])%view the mesh where the nose of the head is pointing down and to the right
%this will give you a good view of where the center of the mesh is
%now the mesh should be centered at (0,0,0)

% Look at the inside
m1.nodes=meshLD.nodes;m1.elements=meshLD.elements;
m2=CutMesh(m1,find(m1.nodes(:,2)>0)); %cut LD mesh in half
[Ia,Ib]=ismember(m2.nodes,meshLD.nodes,'rows');Ib(Ib==0)=[];
PlotMeshSurface(m2,pM); %visualize back half of LD mesh
view([135,30])%view the mesh where the nose of the head is pointing down and to the right
%this will give you a good view of where the center of the mesh is
%now the mesh should be centered at (0,0,0)


%% Load initial Pad file depending on data structure:
padname='AdultV24x28';
load(['Pad_',padname,'.mat']); % Pad file
tpos=cat(1,info.optodes.spos3,info.optodes.dpos3);
Ns=size(info.optodes.spos3,1);
Nd=size(info.optodes.dpos3,1);
rad=info;

paramsFoci.color=cat(1,repmat([1,0,0],Ns,1),repmat([0,0,1],Nd,1));
paramsFoci.color(1,:) = [1 0.4 0.6]; % pink for s1
paramsFoci.color(Ns+1,:) = [0.3010, 0.7450, 0.9330]; %light blue for d1
paramsFoci.radius = 2; %set radius of spheres to 2

%Visualize Pad
Draw_Foci_191203(tpos, paramsFoci);view([-40,30]) %3D plot of pad
PlotCap(info) %2D plot of pad

% Visualize LD mesh with optode positions 
% Optodes have not been relaxed onto mesh yet
PlotMeshSurface(meshLD,pM);Draw_Foci_191203(tpos, paramsFoci);
view([0,90])%view the mesh where the nose of the head is pointing down and to the right
%this will give you a good view of where the center of the mesh is
%you should be able to see a few optodes peeking out of the mesh
%the array has yet to be lined up correctly, which you'll do in the next section with hummingbird


%% AlignMe Section:  Move grid from arbitrary location to approximate target on mesh, Relax grid on head and view
% Location of atlas fiducials: Nasion; Left preauricular point; Right preauricular point; Cz; Inion
atlasFiducials = [- 0.65,  -84.1, -31.88; ... % If using mesh2EEG- EEGPts(1,:)
                   80.78,  15.95, -41.89; ... % EEGPts(155,:)
                  -80.78,  15.95, -41.89; ... % EEGPts(175,:)
                   0.233,   9.63, 97.296; ... % EEGPts(165,:)
                   -0.65, 117.69, -11.78];    % EEGPts(329,:)
               
% Create an instance of our custom DataStorage HANDLE class to store variables
ds = DataStorage(); 

% Input structure
ds.dI.tpos = tpos;       
ds.dI.mesh = meshLD;
ds.dI.pad = info;
ds.dI.pM = pM;
ds.dI.paramsFoci = paramsFoci;
ds.dI.Ns = Ns;
ds.dI.photoPath = ''; %no participant photos for this tutorial, so leave as empty string
ds.dI.imageType = '.jpeg'; %type of image file used, default is .jpeg
ds.dI.atlasFiducials = atlasFiducials;

% Create an instance of your App Designer application, passing variable 'ds' into the app.  
% The code is stuck at this line until your app closes, which destroys 'myapp'
% But the data is assigned to ds variable in workspace
% Run AlignMe
myapp=AlignMe_2020b(ds);
while isvalid(myapp); pause(0.1); end % Wait for app to close before continuing script

% get relaxed optode positions from AlignMe
tposNew = ds.dO.tpos2_relaxed;

% visualize mesh with relaxed optodes
pM.reg=0;
PlotMeshSurface(meshLD,pM);Draw_Foci_191203(tposNew,paramsFoci)
view(0,0) %posterior view

% Get affine matrix that can be used to transform FROM participant space TO MNI space
if isfield(ds.dO, 'affineTform') %if mesh scaled, affineTform field will exists, save it to workspace
    affine_Subj2MNI = [ds.dO.affineTform, zeros(3,1)];
    save('affine_matrix_Subject_to_MNI.mat', 'affine_Subj2MNI')
else %otherwise, set affine_Subj2MNI to eye(4)
    affine_Subj2MNI = eye(4);
end


%% Generate smaller mask based on optode locations
maskCrop=IceCreamScoopCap4Mask(infoT1,tposNew,mask,40); %40 = mm of mesh to scoop out
PlotSlices(maskCrop,infoT1,pS) %visualize smaller mask


%% Generate High Density Head Mesh (Can take up to 15 min)
%If you get an error when running NirfastMesh_Region
%try chagning the mesh name and clearing your output directory of all files
meshname=['_HD_Mesh1'];      % Provide a name for your mesh name here, if making multiple meshes, provide a different name for each mesh
param.facet_distance=2.0;   % Node position error tolerance at boundary
param.facet_size=0.8;       % boundary element size parameter
param.cell_size=1.0;        % Volume element size parameter
param.Mode=0;               % Set Mode=0 to make nirfast compliant mesh
param.info=infoT1;          % Make sure info is infoT1 like with LD mesh
param.Offset=[0,0,0];
tic;meshHD=NirfastMesh_Region(maskCrop,meshname,param);toc

%make copy of mesh that only contains nodes and elements for visualization purposes
visMeshHD.nodes = meshHD.nodes;
visMeshHD.elements = meshHD.elements;
PlotMeshSurface(visMeshHD,pM);view([70,60]) %Visualize in index space
%view should be as if you're looking at the mesh from above and on the right side
%mesh will also be angled so that the top left corner is in the top center of the plot

% Put coordinates back in coordinate space
meshHD.nodes=change_space_coords(meshHD.nodes,infoT1,'coord'); %for actual mesh
visMeshHD.nodes=change_space_coords(visMeshHD.nodes,infoT1,'coord'); %for mesh that's visualized
PlotMeshSurface(visMeshHD,pM);view([70,60]) %Visualize in coordinate space
%view should be as if you're looking at the mesh from above and on the right side
%mesh will also be angled so that the top left corner is in the top center of the plot


%% Relax optodes onto HD Mesh
% Create an instance of  DataStorage class
ds_HD = DataStorage();

% Input structure
ds_HD.dI.tpos = tposNew;       
ds_HD.dI.mesh = meshHD;
ds_HD.dI.pad = info;
ds_HD.dI.pM = pM;
ds_HD.dI.paramsFoci = paramsFoci;
ds_HD.dI.Ns = Ns;
ds_HD.dI.photoPath = ''; %no participant photos for this tutorial, so leave as empty string
ds_HD.dI.imageType = '.jpeg'; %type of image file used, default is .jpeg
ds_HD.dI.atlasFiducials = atlasFiducials;

% Create an instance of your App Designer application,
% passing variable 'ds' into the app.  
% The code is stuck at this line until your app closes, which destroys 'myapp'
% But the data is assigned to ds variable in workspace
% Run AlignMe
myapp = AlignMe_2020b(ds_HD);
while isvalid(myapp); pause(0.1); end % Wait for app to close before continuing script

% get relaxed optode positions from AlignMe
tposNew_HD = ds_HD.dO.tpos2_relaxed;

% visualize mesh with relaxed optodes
pM.reg=0;
PlotMeshSurface(visMeshHD,pM);Draw_Foci_191203(tposNew_HD,paramsFoci)
view(0,0) %posterior view


%% Update Pad info Structure 
info.optodes.spos3=tposNew_HD(1:Ns,:);
info.optodes.dpos3=tposNew_HD((Ns+1):end,:);
m=0;
for d=1:Nd
    for s=1:Ns
        m=m+1;
        info.pairs.r3d(m)=norm(tposNew_HD(s,:)-tposNew_HD(d+Ns,:));
        info.pairs.r3d(m+(Ns*Nd))=info.pairs.r3d(m);
    end
end

% histogram of number of measurements per SD vs radius between SD pairs
figure;histogram(info.pairs.r3d,1000);xlabel('R_S_D');ylabel('N_m_e_a_s');

% save pad file
gridname=padname;
save(['Pad',meshname,'_',gridname,'_', '.mat'],'info')


%% PREPARE! --> mesh with grid array in same file set for NIRFAST
mesh=meshHD;
mesh=PrepareMeshForNIRFAST(mesh,[meshname,'_',gridname],tposNew_HD);
PlotMeshSurface(mesh,pM);Draw_Foci_191203(tposNew_HD, paramsFoci);
view([0,0]) %posterior view


% One last visualization check...
% A portion of the mesh will be removed for this specific visualization
% This is so you can evaluate whether the optodes are placed too deeply within the mesh
m3=CutMesh(mesh,intersect(find(mesh.nodes(:,3)>0),find(mesh.nodes(:,1)>0)));
[Ia,Ib]=ismember(m3.nodes,mesh.nodes,'rows');Ib(Ib==0)=[];
m3.region=mesh.region(Ib);
PlotMeshSurface(m3,pM);Draw_Foci_191203(tposNew_HD, paramsFoci);
view([135,40]) %view mesh from behind with an angled top-down view



%% Calculate Sensitivity Profile
% Set flags
flags.tag=[padname,'_on',meshname,'_test'];
flags.gridname=gridname;
flags.meshname=meshname;
flags.head='info';
flags.info=infoT1;                  % Your T1 info file
flags.gthresh=1e-5;                 % Voxelation threshold in G
flags.voxmm=2;                      % Voxelation resolution (mm)
flags.labels.r1='csf';              % Regions for optical properties
flags.labels.r2='white';
flags.labels.r3='gray';
flags.labels.r4='bone';
flags.labels.r5='skin';
flags.op.lambda=[750,850];          % Wavelengths (nm)
flags.op.mua_skin=[0.0170,0.0190];  % Baseline absorption
flags.op.mua_bone=[0.0116,0.0139];
flags.op.mua_csf=[0.0040,0.0040];
flags.op.mua_gray=[0.0180,0.0192];
flags.op.mua_white=[0.0167,0.0208];
flags.op.musp_skin=[0.74,0.64];     % Baseline reduced scattering coeff
flags.op.musp_bone=[0.94,0.84];
flags.op.musp_csf=[0.3,0.3];
flags.op.musp_gray=[0.8359,0.6726];
flags.op.musp_white=[1.1908,1.0107];
flags.op.n_skin=[1.4,1.4];          % Index of refraction
flags.op.n_bone=[1.4,1.4];
flags.op.n_csf=[1.4,1.4];
flags.op.n_gray=[1.4,1.4];
flags.op.n_white=[1.4,1.4];
flags.srcnum=Ns;                    % Number of sources
flags.t4=affine_Subj2MNI;           % Affine matrix for going from subject-specific space to MNI space
flags.t4_target='MNI'; % string
flags.makeA=1; % don't make A, just make G
flags.Hz=0;
if flags.Hz, flags.tag = [flags.tag,'FD']; end

% Run makeAnirfast to get sensitivity matrix
Ti=tic;[A,dim,Gsd]=makeAnirfaster(mesh,flags); % size(A)= [Nwl, Nmeas, Nvox]
disp(['<makeAnirfast took ',num2str(toc(Ti))])

%Gsd vs Rsd provides a simulated light fall-off
figure;semilogy(info.pairs.r3d(1:(Ns*Nd)),Gsd,'*');
legend({num2str(flags.op.lambda')})
xlabel('R_S_D [mm]');ylabel('G_S_D');xlim([0,100])


%% Package data and save A
[Nwl,Nmeas,Nvox]=size(A);
A=reshape(permute(A,[2,1,3]),Nwl*Nmeas,Nvox);

% Place spatial information about light model in info.tissue structure
info.tissue.dim=dim;
info.tissue.affine=flags.t4;
info.tissue.infoT1=infoT1;
info.tissue.affine_target='MNI';
info.tissue.flags=flags;

save(['A_',flags.tag,'.mat'],'A','info','-v7.3') %save A


%% Visualize aspects of sensitivity profile
t1=affine3d_img(mask,infoT1,dim,eye(4)); % put anatomical volume in dim space

% Visualize Single SD pair (S1 D6)
keep=info.pairs.WL==2 & info.pairs.Src==1 & info.pairs.Det==6; % SD pair set here
foo=squeeze(A(keep,:));              % Single meas pair
fooV=Good_Vox2vol(foo',dim);
fooV=fooV./max(fooV(:));
fooV=log10(1e2.*fooV);                  % top 2 o.o.m.
pA.PD=1;pA.Scale=2;pA.Th.P=0;pA.Th.N=-pA.Th.P;
PlotSlices(t1,dim,pA,fooV)

% Visualize Single SD pair (S2 D2)
keep=info.pairs.WL==2 & info.pairs.Src==2 & info.pairs.Det==2; % SD pair set here
foo=squeeze(A(keep,:));              % Single meas pair
fooV=Good_Vox2vol(foo',dim);
fooV=fooV./max(fooV(:));
fooV=log10(1e2.*fooV);                  % top 2 o.o.m.
pA.PD=1;pA.Scale=2;pA.Th.P=0;pA.Th.N=-pA.Th.P;
PlotSlices(t1,dim,pA,fooV)

% FFR
keep=(info.pairs.WL==2 & info.pairs.r2d<=40);
a=squeeze(A(keep,:));
iA=Tikhonov_invert_Amat(a,0.01,0.1);
iA=smooth_Amat(iA,dim,5); %5 = smoothing parameter
ffr=makeFlatFieldRecon(a,iA);

fooV=Good_Vox2vol(ffr,dim);
fooV=fooV./max(fooV(:));
pA.PD=1;pA.Scale=1;pA.Th.P=1e-2;pA.Th.N=-pA.Th.P;
PlotSlices(t1,dim,pA,fooV)


%% Visualize alignment of LD mesh, HD mesh, array, and cortices
%  View alignment
Cortical_mesh = load(['MNI164k_big.mat']);
Anat.CtxL = Cortical_mesh.MNIl;
Anat.CtxR = Cortical_mesh.MNIr;
pA0l=struct;
figure('Color','k','Position',[500,100,1050,1000])
pM1.fig_handle=gca;pM1.FaceColor=[1,1,1].*0.5;pM1.EdgeColor = 'none';
pM1.SpecularStrength=0.25;pM1.DiffuseStrength=0.5;
PlotMeshSurface(meshHD,pM1)                        % HD mesh
pA0l.fig_handle=gca;
pA0l.FaceColor=[1,1,1].*0.5;pA0l.EdgeColor='none';
pA0l.AmbientStrength=0.25;pA0l.DiffuseStrength=0.25;
pA0l.SpecularStrength=0.025;
PlotMeshSurface(Anat.CtxL,pA0l)                 % Cortical Surfaces
PlotMeshSurface(Anat.CtxR,pA0l)
pA2.fig_handle=gca;
pA2.FaceAlpha=0;
pA2.EdgeAlpha=0.5;
pA2.EdgeColor=[1,1,1].*0.25;
PlotMeshSurface(meshLD,pA2)                     % LD mesh
set(gca,'Color','k')
view([90,0]) % display in Lateral view
Draw_Foci_191203(cat(1,info.optodes.spos3,info.optodes.dpos3),paramsFoci)  % Array
axis off

% For posterior view, uncomment the following line
% view([0,0])



%% Clear and Close everything before next section
clear all; close all


%% Section 2: Full Data Processing and Visualization %%

%% FULL DATA PROCESSING AND VISUALIZATION OVERVIEW
% Load data from NeuroDOT_Data_Sample_CCW2.mat. 
% Generate raw data quality figures
% Pre-process data
% Visualize effects of filtering and pre-processing
% Block-Average data and visualize
% Reconstruction
% Spectroscopy
% Visualize data in anatomical space with PlotSlices
% Visualize data in mesh space with PlotMeshSurface


%% Load Measurement data
dataset='CCW2'; % CCW1, CCW2, CW1, IN1, OUT1, GV1, HW1, HW2, HW3_Noisy,  RW1
load(['NeuroDOT_Data_Sample_',dataset,'.mat']); % data, info, flags

% Set parameters for A and block length for quick processing examples
A_fn='A_AdultV24x28_on_HD_Mesh1_test.mat';    % Sensitivity Matrix
dt=36;                      % Block length
tp=16;                      % Example (block averaged) time point



%% General data QC with synchpts if present
Plot_RawData_Time_Traces_Overview(data,info);           % Time traces
Plot_RawData_Cap_DQC(data,info);                        % Cap-relevant views
Plot_RawData_Metrics_II_DQC(data,info)                  % Raw data quality figs


%% Pre-process data
% Pre-processing pipeline, but only the part before filtering data
lmdata_b4_filt = logmean(data);                                           % Logmean Light Levels
info_b4_filt  = FindGoodMeas(lmdata_b4_filt , info, 0.075);               % 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 = logmean(data);                                                   % Logmean Light Levels
info = FindGoodMeas(lmdata, info, 0.075);                                 % 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);                                               % Superficial Signal Regression
[lmdata, ~] = regcorr(lmdata, info, hem);
lmdata = lowpass(lmdata, 0.5, info.system.framerate);                     % Low Pass Filter 2 (0.5 Hz)
[lmdata, info] = resample_tts(lmdata, info, 1, 1e-5);                     % 1 Hz Resampling (1 Hz)
info.GVTD = CalcGVTD(lmdata(info.MEAS.GI & info.pairs.r2d<20,:));         % Calculate GVTD


%% Visualize Effects of Pre-Proecssing
% Visualize Pre-filtering
% Measurements to include
keep_b4_filt = info_b4_filt.pairs.WL==2 & info_b4_filt.pairs.r2d < 40 & 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 < 40 & 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(mean(lmdata(keep,:),1)),info.system.framerate); % Generate average spectrum
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
Plot_TimeTrace_With_PowerSpectrum(lmdata,info); % Data, Grayplots organized by measurement #, Powerspectrum
GrayPlots_Rsd_by_Wavelength(lmdata,info); %Grayplots organized by SD separation


%% Block Averaging the measurement data and view
badata = BlockAverage(lmdata, info.paradigm.synchpts(info.paradigm.Pulse_2), dt);

badata=bsxfun(@minus,badata,mean(badata,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 #')


%% RECONSTRUCTION PIPELINE
if ~exist('A', 'var')       % In case running by hand or re-running script
    A=load([A_fn],'info','A');
    if length(size(A.A))>2  % A data structure [wl X meas X vox]-->[meas X vox]
        [Nwl,Nmeas,Nvox]=size(A.A);
        A.A=reshape(permute(A.A,[2,1,3]),Nwl*Nmeas,Nvox);
    end        
end
Nvox=size(A.A,2);
Nt=size(lmdata,2);
cortex_mu_a=zeros(Nvox,Nt,2);

for j = 1:2
    keep = (info.pairs.WL == j) & (info.pairs.r2d <= 40) & info.MEAS.GI;
    disp('> Inverting A')                
    iA = Tikhonov_invert_Amat(A.A(keep, :), 0.01, 0.1); % Invert A-Matrix (Ari sets lambda1 to 0.1 for better vis in pstt slides)
    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);
cortex_HbR = cortex_Hb(:, :, 2);
cortex_HbT = cortex_HbO + cortex_HbR;


%% Select Volumetric visualizations of block averaged data
if ~exist('MNI', 'var')
[MNI,infoB]=LoadVolumetricData('Segmented_MNI152nl_on_MNI111_nifti',[],'nii'); % Load MRI (same data set as in A matrix dim)
end
MNI_dim = affine3d_img(MNI,infoB,A.info.tissue.dim,eye(4),'nearest'); % Transform to DOT volume space

% Block Average Data
badata_HbO = BlockAverage(cortex_HbO, info.paradigm.synchpts(info.paradigm.Pulse_2), dt);
badata_HbO=bsxfun(@minus,badata_HbO,badata_HbO(:,1));
badata_HbOvol = Good_Vox2vol(badata_HbO,A.info.tissue.dim);
tp_Eg=squeeze(badata_HbOvol(:,:,:,tp));

% Explore PlotSlices - The basics
PlotSlices(MNI_dim)                             % Anatomy only
PlotSlices(MNI_dim,A.info.tissue.dim)           % Anatomy + volumetric data
PlotSlices(MNI_dim,A.info.tissue.dim,[],tp_Eg); % Anatomy + volumetric data + functional data

% Match your view of PlotSlices to the images in the tutorial PPT (slide 68)
PlotSlices(MNI_dim,A.info.tissue.dim,[],tp_Eg);
NoBA_HbO = Good_Vox2vol(cortex_HbO, A.info.tissue.dim); %get non block averaged data in time
PlotSlicesTimeTrace(MNI_dim, A.info.tissue.dim, [], NoBA_HbO, info)

% Explore PlotSlices - in depth (Slides 69-72)
% Defaults only, no params structure as input
PlotSlices(MNI_dim,A.info.tissue.dim,[],tp_Eg);
% Params 1
Params1.PD = 1; % Turn Positive Definite On
PlotSlices(MNI_dim,A.info.tissue.dim,Params1,tp_Eg);
% Params 2
Params2.Th.P = 0.001;
Params2.Th.N = - Params2.Th.P;
PlotSlices(MNI_dim,A.info.tissue.dim,Params2,tp_Eg);
% Params 3
Params3.Th.P = 0.001;
Params3.Th.N = -0.004;
PlotSlices(MNI_dim,A.info.tissue.dim,Params3,tp_Eg);
% Params 4
Params4.Scale = 2*max(abs(tp_Eg(:)));
PlotSlices(MNI_dim,A.info.tissue.dim,Params4,tp_Eg);
% Params 5
Params5.Cmap = 'cool';
PlotSlices(MNI_dim,A.info.tissue.dim,Params5,tp_Eg);
% Params 6
Params6.Cmap.P = 'winter';
Params6.Cmap.N = 'summer';
PlotSlices(MNI_dim,A.info.tissue.dim,Params6,tp_Eg);
% Params 7
Params7.cbmode = 1;
Params7.cbticks = [0,0.25,0.5,0.75,1];
PlotSlices(MNI_dim,A.info.tissue.dim,Params7,tp_Eg);
% Params 8
Params8.cbmode = 1;
Params8.cblabels = {'Hello', 'Goodbye'};
PlotSlices(MNI_dim,A.info.tissue.dim,Params8,tp_Eg);
% Params 9 
Params9.cbmode = 1;
Params9.cbticks = [0,0.25,0.5,0.75,1];
Params9.cblabels = {'1', '2', '3', '4', '5'};
Params9.Cmap.P = 'winter';
Params9.Cmap.N = 'summer';
Params9.Scale = 0.8*max(abs(tp_Eg(:)));
Params9.Th.P = 0.001;
Params9.Th.N = -0.004;
PlotSlices(MNI_dim,A.info.tissue.dim,Params9,tp_Eg);


% Visualize the data (slide 73)
PlotSlices(tp_Eg,A.info.tissue.dim);            % Data by itself
PlotSlices(MNI_dim,A.info.tissue.dim,[],tp_Eg); % Data with anatomical underlay
% Set parameters to visualize more specific aspects of data
Params.Scale=0.8*max(abs(tp_Eg(:)));     % Scale wrt/max of data
Params.Th.P=0.5*Params.Scale;            % Threshold to see strong activations
Params.Th.N=-Params.Th.P;                % Thresholds go both ways
Params.Cmap='jet';
PlotSlices(MNI_dim,A.info.tissue.dim,Params,tp_Eg);

% Explore the block-averaged data a bit more interactively
Params.Scale=0.8*max(abs(badata_HbOvol(:)));
Params.Th.P=0;
Params.Th.N=-Params.Th.P;
PlotSlicesTimeTrace(MNI_dim,A.info.tissue.dim,Params,badata_HbOvol,info)

% Explore the not-block-averaged data a bit more interactively
HbOvol = Good_Vox2vol(cortex_HbO,A.info.tissue.dim);
Params.Scale=4e-3;
Params.Th.P=1e-3;
Params.Th.N=-Params.Th.P;
PlotSlicesTimeTrace(MNI_dim,A.info.tissue.dim,Params,HbOvol,info)


%% Select Surface visualizations
if ~exist('MNIl', 'var'),load(['MNI164k_big.mat']);end
HbO_atlas = affine3d_img(badata_HbOvol,A.info.tissue.dim,infoB,eye(4));
tp_Eg_atlas=squeeze(HbO_atlas(:,:,:,tp));
pS=Params;
pS.view='post'; % Posterior view
pS.ctx='std'; % Standard pial cortical view
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, infoB, pS); %(Slide 34)


% Play around with visualization parameters
% View 
% Dorsal view
pS.view = 'dorsal';
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, infoB, pS);
% Lateral view
pS.view = 'lat';
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, infoB, pS);

% Inflation 
pS.view='post'; %reset to posterior view before visualizing
% Standard pial cortical view
pS.ctx='std'; 
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, infoB, pS);
% Inflated pial cortical view
pS.ctx='inf'; 
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, infoB, pS);
% Very Inflated pial cortical view
pS.ctx='vinf';
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, infoB, pS);

% Thresholds 
pS.ctx='std'; %reset to standard inflation before visualizing
% Set to 0
pS.Th.P = 0;
pS.Th.N = -pS.Th.P;
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, infoB, pS);
% Set a bit higher
pS.Th.P = 1e-4;
pS.Th.N = -pS.Th.P;
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, infoB, pS);
% Set even higher
pS.Th.P = 1e-3;
pS.Th.N = -pS.Th.P;
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, infoB, pS);
% Set very high
pS.Th.P = 1e-2;
pS.Th.N = -pS.Th.P;
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, infoB, pS);

% Scale
pS.Th.P = 0;
pS.Th.N = -pS.Th.P; %reset thresholds first
% Set to 4e-3
pS.Scale = 4e-3;
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, infoB, pS);
% Set lower
pS.Scale = 1e-3;
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, infoB, pS);
% Set higher
pS.Scale = 9e-3;
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, infoB, pS);

% Colormaps
% reset scale and thresholds first
pS.Scale = 4e-3; 
pS.Th.P = 1e-3;
pS.Th.N = -pS.Th.P;
% Same colormap for P/N
pS.Cmap = 'cool';
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, infoB, pS);
% Differnt colormaps for P/N
pS = rmfield(pS, 'Cmap'); %remove Cmap field so we can replace it with a struct containing both positive and negative colormaps
pS.Cmap.P = 'winter';
pS.Cmap.N = 'summer';
PlotInterpSurfMesh(tp_Eg_atlas, MNIl,MNIr, infoB, pS);


%% True Color Parameter
params.Cmap = 'jet';
% Get data at another timepoint
tp_Eg_atlas_2=squeeze(HbO_atlas(:,:,:,tp+9)); %adding 9 to tp gets us data when the checkerboard is in a different quadrant

% Make mask of timepoint = 16 data
mask_temp_1 = tp_Eg_atlas;
mask1 = zeros(size(mask_temp_1));
mask1(mask_temp_1 > pS.Th.P) = 1; %set tp1 positive activations to 1
PlotInterpSurfMesh(mask1, MNIl,MNIr, infoB, pS);
mask_temp_2 = tp_Eg_atlas_2;
mask2 = zeros(size(mask_temp_2));
mask2(mask_temp_2 > pS.Th.P) = 1; %set tp2 positive activations to 1
PlotInterpSurfMesh(mask2, MNIl,MNIr, infoB, pS);

% Combine masks to plot
mask3 = zeros(size(mask_temp_1));
mask3(mask1 == 1 & mask2 == 1) = 3;
mask3(mask1 == 1 & mask3 ~= 3) = 1;
mask3(mask2 == 1 & mask3 ~= 3) = 2;

% Set parameters
pS_TC = pS;
pS_TC.Cmap = [1,0,0; 0,0,1; 1,0,1];
pS_TC.TC = 1; %True Color On
pS_TC.OL = 1; %Overlap On
PlotInterpSurfMesh(mask3, MNIl,MNIr, infoB, pS_TC); 


%% Gif Making
%re-reconstruct + do spectroscopy
for j = 1:2
    keep1 = (info.pairs.WL == j) & (info.pairs.r2d <= 40) & info.MEAS.GI;
    disp('> Inverting A')                
    iA1 = Tikhonov_invert_Amat(A.A(keep1, :), 0.01, 0.1); % Invert A-Matrix
    disp('> Smoothing iA')
    iA1 = smooth_Amat(iA1, A.info.tissue.dim, 5);         % Smooth Inverted A-Matrix      
    cortex_mu_a1(:, :, j) = reconstruct_img(lmdata(keep1, :), iA1);% Reconstruct Image Volume
end
cortex_Hb1 = spectroscopy_img(cortex_mu_a1, E); %Spectroscopy

%Params and BA
pS_gif = pS;
pS_gif.GVTD_censor = 1; % do you want to remove motion? 1 = remove motion and 0 = don't remove motion
pS_gif.GVTD_th = 0.007; % motion above whatever number is here is removed
[BA_out,~,BT_out] = BlockAverage(squeeze(cortex_Hb1(:,:,1)),...
    info.paradigm.synchpts(info.paradigm.Pulse_2), dt); 
BAvol=Good_Vox2vol(BA_out,A.info.tissue.dim);
fooMNI=affine3d_img(squeeze(mean(BAvol(:,:,:,12:16),4)),... %data to the tissue, run this to set scale for gifs
    A.info.tissue.dim,infoB,info.tissue.affine);

%Actually make the gifs
% Please make sure you are in your own folder when running these.
filename = append('CCW2', '_','_test5', '.gif'); %sets file name
% Code that makes gif
for i = 1:1:dt %1:1:36 for retinotopy, 1:1:30 for auditory
    % Get map for each timepoint and visualize it
    fooMNI=affine3d_img(squeeze(mean(BAvol(:,:,:,i),4)),A.info.tissue.dim,infoB,A.info.tissue.affine);
    pS_gif.Scale=max(fooMNI(fooMNI(:)~=0)); %option to set scale for each indivdual image here
    pS_gif.Th.P=0.25*max(fooMNI(fooMNI(:)~=0));
    pS_gif.PD = 1;
    PlotInterpSurfMesh(fooMNI,MNIl,MNIr,infoB,pS_gif); %for Gates/BigCap data
    % change the title for each run/view
    title(['CCW2 Frame ' , num2str(i)], 'color','white')
    frame = getframe(gcf);
    clear axis
    im{i} = frame2im(frame);

end
% Change the filename for each run
figure;
for i = 1:dt %1:1:36 for retinotopy, 1:1:30 for auditory
    subplot(6,6,i) %(6,6,i) for retinotopy, (6,5,i) for auditory
    imshow(im{i});
end
for i = 1:dt %1:1:36 for retinotopy, 1:1:30 for auditory
    [Z,map] = rgb2ind(im{i},256);
    if i == 1
        imwrite(Z,map,filename,'gif','LoopCount',Inf,'DelayTime',0.2);
    else
        imwrite(Z,map,filename,'gif','WriteMode','append','DelayTime',0.2);
    end
end
close all





