Commit 0b37f084 authored by Julius Welzel's avatar Julius Welzel
Browse files

"Almost" everything improved

parent 005acd0b

05_scripts/ana07_TOPOdiff.m

+245 −117
Original line number Diff line number Diff line
@@ -37,7 +37,8 @@ young = (SUBJ<80); 
% define trails of interest

ME_trials = [cfg.blck_idx{[1 4],1}];

MI_trials = [cfg.blck_idx{[2 3],1}];



cfg.AREA.channel = cfg.n_chan(~ismember(cfg.n_chan,[19,20]));

cfg.AREA.chanlocs = cfg.chanlocs(~ismember(cfg.n_chan,[19,20]));



FB = {'mue','beta','broad'};
@@ -50,7 +51,7 @@ load([PATHIN_ERD 'ERD_' FB{fb} '.mat']); 
load([PATHIN_EMG 'emg_all.mat']);



%% Zscore data over all participants only for MI trials

z_top = zscore(nanmean(m_rERD_trial(:,:,MI_trials),3),[],2);

z_top = zscore(nanmean(m_rERD_trial(:,cfg.AREA.channel,MI_trials),3),[],2);

top_old = nanmean(z_top(old,:));

top_young = nanmean(z_top(young,:));
@@ -60,10 +61,10 @@ top_young = nanmean(z_top(young,:)); 
%

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

%%

p_interp = 100; % points to interpolate

p_interp = 200; % points to interpolate



% get cartesian coordinates by MIKE X COHEN :D and grid data

[elocsX,elocsY] = pol2cart(pi/180*[cfg.chanlocs.theta],[cfg.chanlocs.radius]);

[elocsX,elocsY] = pol2cart(pi/180*[cfg.AREA.chanlocs.theta],[cfg.AREA.chanlocs.radius]);

elocsX = elocsX*-1; %swap left/ right

xlin = linspace(min(elocsX),max(elocsX),p_interp);

ylin = linspace(min(elocsY),max(elocsY),p_interp);
@@ -75,43 +76,21 @@ d_min = min(d_mat(d_mat > 0)); 
d_max = max(d_mat(d_mat > 0));



% calculate spatial frequency

spat_freq = d_max/(2*d_min);

n_level = 6; %number of rings from equiDist layout EASYCAP

n_intersect = ceil(n_level*spat_freq+(n_level*spat_freq)*0.2); % round up plus extra 20% for extra resolution





%% calculate surface area  of topoplot for each trial

%{

surf_old = {};

surf_young = {};



for s = find(young)

    disp(['SUBJ: ' num2str(SUBJ(s))])

   for t = MI_trials

       

       topo = zscore((m_rERD_trial(s,:,t)));

       surf_st = griddata(elocsX,elocsY,topo,Xgrid,Ygrid,'cubic');

       top_thresh = min(min(topo))*0.5;

       surf_st(surf_st>top_thresh) = NaN;

       surf_young{end+1} = SurfArea(surf_st,Xgrid,Ygrid);

       

   end

end



for s = find(old)

    disp(['SUBJ: ' num2str(SUBJ(s))])

   for t = MI_trials

spat_freq = (d_max/d_min)*2;

n_level = 4; %number of rings from equiDist layout EASYCAP

n_intersect = ceil(spat_freq); % round up 



       topo = zscore((m_rERD_trial(s,:,t)));

       surf_st = griddata(elocsX,elocsY,topo,Xgrid,Ygrid,'cubic');

       top_thresh = min(min(topo))*0.5;

       surf_st(surf_st>top_thresh) = NaN;

       surf_old{end+1} = SurfArea(surf_st,Xgrid,Ygrid);

% calculate intersect area

maps_st = zscore(nanmean(m_rERD_trial(:,cfg.AREA.channel,MI_trials),3),[],2);



   end

end

%}

load([PATHOUT_TOODIFF 'inter_area_noeye_' FB{fb} '.mat'])



% for s = 1:length(SUBJ)

%     if SUBJ(s)<50;c_map_ind = c_y_map;else;c_map_ind = c_o_map;end

%     area_z_t_avg(s,:) = inter_area(cfg.AREA.chanlocs,maps_st(s,:),n_intersect,SUBJ(s),PATHOUT_IND_PLOT,c_map_ind,FB{fb});

% end

% 

% save([PATHOUT_TOODIFF 'inter_area_noeye_' FB{fb} '.mat'],'area_z_t_avg')



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

%
@@ -121,139 +100,193 @@ end 


% only 3D topoplots over both groups with topoplot

c_lines = 200;

cuts = n_intersect/4;



ax1 = figure;

% ax1 = subplot(1,2,1);

min_z = min(min(top_young));

max_z = max(max(top_young));

lvl_inter = linspace(max_z,min_z,cuts);



figure;

ax1 = subplot(1,2,1);

Z_young = griddata(elocsX,elocsY,top_young,Xgrid,Ygrid,'cubic');

mesh(Xgrid,Ygrid,Z_young) % interpolate

colormap(ax1,c_map)

axis tight;hold on;    

plot3(elocsX,elocsY,top_young,'.k','MarkerSize',15) %plot electrodes

title '3D Topoplot Young'

% title '3D Topoplot Young'

zlabel 'z-Score [\sigma]'

set(gca,'YTickLabel',[]);

set(gca,'XTickLabel',[]);

hold on

%intersect lines

% surf(Xgrid,Ygrid,0*ones(size(Z_young)),'FaceColor','k','FaceAlpha',0.08,'EdgeColor','none')

for ii = 1:cuts

    surf(Xgrid,Ygrid,lvl_inter(ii)*ones(size(Z_young)),'FaceColor','k','FaceAlpha',0.08,'EdgeColor','none')

    hold on

end



%topoplot

[~,h_o] = contourf(Xgrid,Ygrid,Z_young,c_lines,'LineStyle','none');       %# get handle to contourgroup object

colormap (ax1,c_map)

hold on

[~,h_o_c] = contour(Xgrid,Ygrid,Z_young,cuts,'k');       %# get handle to contourgroup object



%# change the ZData property of the inner patches

h_o.ContourZLevel = 1.1*min(min(Z_young));

h_o.ContourZLevel = 1.3*min(min(Z_young));

h_o_c.ContourZLevel = 1.299*min(min(Z_young));



pic_r = 1/1;

save_fig(gcf,PATHOUT_TOODIFF,['3D_top_young_' FB{fb}],'figtype','.png','figsize',[0 0 22.7 pic_r*22.7],'fontsize',20)



% OLD 3D PLOT

ax2 = subplot(1,2,2);

ax2 = figure

% ax2 = subplot(1,2,2);

min_z = min(min(top_old));

max_z = max(max(top_old));

lvl_inter = linspace(max_z,min_z,cuts);



Z_old = griddata(elocsX,elocsY,top_old,Xgrid,Ygrid,'cubic');

mesh(Xgrid,Ygrid,Z_old) % interpolate

colormap(ax2,c_map)

axis tight;hold on;

plot3(elocsX,elocsY,top_old,'.k','MarkerSize',15)

title '3D Topolot Old'

% title '3D Topolot Old'

zlabel 'z-Score [\sigma]'

set(gca,'YTickLabel',[]);

set(gca,'XTickLabel',[]);

hold on

%intetsect planes

% surf(Xgrid,Ygrid,0*ones(size(Z_old)),'FaceColor','k','FaceAlpha',0.08,'EdgeColor','none')

% hold on

for ii = 1:cuts

    surf(Xgrid,Ygrid,lvl_inter(ii)*ones(size(Z_old)),'FaceColor','k','FaceAlpha',0.08,'EdgeColor','none')

    hold on

end



%topoplot

[~,h_o] = contourf(Xgrid,Ygrid,Z_old,c_lines,'LineStyle','none');       %# get handle to contourgroup object

colormap (ax1,c_map)

colormap (ax2,c_map)

hold on

[~,h_o_c] = contour(Xgrid,Ygrid,Z_old,cuts,'k');       %# get handle to contourgroup object



%# change the ZData property of the inner patches

h_o.ContourZLevel = 1.3*min(min(Z_old));

h_o_c.ContourZLevel = 1.299*min(min(Z_old));



save_fig(gcf,PATHOUT_TOODIFF,['3D_tops_' FB{fb}],'figtype','.png','figsize',[0 0 55 25],'fontsize',20)

pic_r = 1/1;

save_fig(gcf,PATHOUT_TOODIFF,['3D_top_old_' FB{fb}],'figtype','.png','figsize',[0 0 22.7 pic_r*22.7],'fontsize',20)



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

%

%                  Calculate area of intersection

%

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





gd(1).Z= griddata(elocsX,elocsY,top_young,Xgrid,Ygrid,'cubic'); %young

gd(2).Z = griddata(elocsX,elocsY,top_old,Xgrid,Ygrid,'cubic'); %old

group = {'Young','Old'};



plane_sect = linspace(0,-1.2,6);



for g = 1:length(gd)

    

    figure

    for i = 1:length(plane_sect)

        subplot(2,3,i)

        mesh(Xgrid,Ygrid,gd(g).Z) % interpolate

        hold on

        surf(Xgrid,Ygrid,plane_sect(i)*ones(size(gd(g).Z)),'FaceColor','k','FaceAlpha',0.08,'EdgeColor','none')

        hold on

% save_fig(gcf,PATHOUT_TOODIFF,['3D_tops_' FB{fb}],'figtype','.png','figsize',[0 0 55 25],'fontsize',20)



        % Take the difference between the two surface heights and find the contour

        % where that surface is zero.

        zdiff = gd(g).Z - (plane_sect(i)*ones(size(gd(g).Z)));

        C = contour(Xgrid, Ygrid, zdiff,[0,0]);



        [a, ca] = cont_area(C,Xgrid,Ygrid,gd(g).Z);

        % Visualize the line.

        cl = line(ca.xL, ca.yL, ca.zL, 'Color', 'k', 'LineWidth', 3);

        title(['Area : ' num2str(a) ' at z = ' num2str(plane_sect(i))])

    end

    sgtitle ([group{g} ' participants'])

    save_fig(gcf,PATHOUT_TOODIFF,['SurfArea_' FB{fb} '_' group{g}],'figtype','.fig')



end

%% Plot topographie interpolated with surface area fun

c_idx = ceil(linspace(56,9,n_intersect));



maps_st = zscore(nanmean(m_rERD_trial(:,:,MI_trials),3),[],2);

[v i] = max(mean(area_z_t_avg(old,:))-mean(area_z_t_avg(young,:))); % max diff in area

% find area where old is bigger



for s = 1:length(SUBJ)

    if SUBJ(s)<50;c_map_ind = c_y_map;else;c_map_ind = c_o_map;end

    area_z_t_avg(s,:) = inter_area(cfg.chanlocs,maps_st(s,:),n_intersect,SUBJ(s),PATHOUT_IND_PLOT,c_map_ind,FB{fb});

end

[m_a_old, se_a_old] = mean_SEM(area_z_t_avg(old,:)');

[m_a_young, se_a_young] = mean_SEM(area_z_t_avg(young,:)');

min_z_y = min(min(top_young));

max_z_y = max(max(top_young));

z_space_y = linspace(max_z_y,min_z_y,n_intersect);

min_z_o = min(min(top_old));

max_z_o = max(max(top_old));

z_space_o = linspace(max_z_o,min_z_o,n_intersect);

lev_diff_dist = [0 0.0015 0.001]

% i = ([m_a_old-se_a_old]-[m_a_young+se_a_young])>lev_diff_dist(fb);

full_w_y = fwhm(1:n_intersect,m_a_young);

full_w_o = fwhm(1:n_intersect,m_a_old);

i_y = [round(find(m_a_young==max(m_a_young))-0.5*full_w_y):round(find(m_a_young==max(m_a_young))+0.5*full_w_y)];

i_o = [round(find(m_a_old==max(m_a_old))-0.5*full_w_o):round(find(m_a_old==max(m_a_old))+0.5*full_w_o)];



[m_a_old se_a_old] = mean_SEM(area_z_t_avg(old,:)');

[m_a_young se_a_young] = mean_SEM(area_z_t_avg(young,:)');



save([PATHOUT_TOODIFF 'inter_area_' FB{fb} '.mat'],'area_z_t_avg')



%% Plot topographie interpolated with surface area fun

[v i] = max(mean(area_z_t_avg(old,:))-mean(area_z_t_avg(young,:))); % max diff in area



figure

ax1 = subplot(2,7,[1 2]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Young Topo%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ax1 = subplot(2,9,[6 7]);

Z_young = griddata(elocsX,elocsY,top_young,Xgrid,Ygrid,'cubic');

sy = mesh(Xgrid,Ygrid,Z_young); % interpolate

view([-90,1])

py = surf2patch(sy);

colormap(ax1,c_y_map)

axis tight;hold on;

plot3(elocsX,elocsY,top_young,'.k','MarkerSize',15)

title 'Interpolated Topoplot Young'

% title 'Interpolated Topoplot Young'

hold on

surf(Xgrid,Ygrid,v*ones(size(Z_young)),'FaceAlpha',0.2,'EdgeColor','none')

% up_i = find(i_y,'first');

surf(Xgrid,Ygrid,z_space_y(i_y(1))*ones(size(Z_young)),'FaceAlpha',0.2,'EdgeColor','none')

hold on

% down_i = find(i_y,1,'last');

surf(Xgrid,Ygrid,z_space_y(i_y(end))*ones(size(Z_young)),'FaceAlpha',0.2,'EdgeColor','none')



zlabel 'z-Score [\sigma]'

set(gca,'YTickLabel',[]);

set(gca,'XTickLabel',[]);



%side

ax1 = subplot(2,9,[8 9]);

Z_young = griddata(elocsX,elocsY,top_young,Xgrid,Ygrid,'cubic');

sy = mesh(Xgrid,Ygrid,Z_young); % interpolate

view([0,1])

py = surf2patch(sy);

colormap(ax1,c_y_map)

axis tight;hold on;

plot3(elocsX,elocsY,top_young,'.k','MarkerSize',15)

% title 'Interpolated Topoplot Young'

hold on

% up_i = find(i_y,1,'first');

surf(Xgrid,Ygrid,z_space_y(i_y(1))*ones(size(Z_young)),'FaceAlpha',0.2,'EdgeColor','none')

hold on

% down_i = find(i_y,1,'last');

surf(Xgrid,Ygrid,z_space_y(i_y(end))*ones(size(Z_young)),'FaceAlpha',0.2,'EdgeColor','none')



% zlabel 'z-Score [\sigma]'

set(gca,'YTickLabel',[]);

set(gca,'XTickLabel',[]);



ax2 = subplot(2,7,[7 8]+1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Old topo%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ax2 = subplot(2,9,[15 16]);

Z_old = griddata(elocsX,elocsY,top_old,Xgrid,Ygrid,'cubic');

mesh(Xgrid,Ygrid,Z_old) % interpolate

view([-90,1])

colormap(ax2,c_o_map)

axis tight;hold on;

plot3(elocsX,elocsY,top_old,'.k','MarkerSize',15)

title 'Interpolated Topoplot Old'

% title 'Interpolated Topoplot Old'

hold on

% up_i = find(i_o,1,'first');

surf(Xgrid,Ygrid,z_space_y(i_o(1))*ones(size(Z_old)),'FaceAlpha',0.2,'EdgeColor','none')

hold on

surf(Xgrid,Ygrid,v*ones(size(Z_old)),'FaceAlpha',0.2,'EdgeColor','none')

% down_i = find(i_o,1,'last');

surf(Xgrid,Ygrid,z_space_y(i_o(end))*ones(size(Z_young)),'FaceAlpha',0.2,'EdgeColor','none')



zlabel 'z-Score [\sigma]'

set(gca,'YTickLabel',[]);

set(gca,'XTickLabel',[]);



%side

ax2 = subplot(2,9,[17 18]);

Z_old = griddata(elocsX,elocsY,top_old,Xgrid,Ygrid,'cubic');

mesh(Xgrid,Ygrid,Z_old) % interpolate

view([0,1])

colormap(ax2,c_o_map)

axis tight;hold on;

plot3(elocsX,elocsY,top_old,'.k','MarkerSize',15)

% title 'Interpolated Topoplot Old'

hold on

% up_i = find(i_o,1,'first');

surf(Xgrid,Ygrid,z_space_y(i_o(1))*ones(size(Z_old)),'FaceAlpha',0.2,'EdgeColor','none')

hold on

% down_i = find(i_o,1,'last');

surf(Xgrid,Ygrid,z_space_y(i_o(end))*ones(size(Z_young)),'FaceAlpha',0.2,'EdgeColor','none')



% zlabel 'z-Score [\sigma]'

set(gca,'YTickLabel',[]);

set(gca,'XTickLabel',[]);



ax3 = subplot(2,7,[4:7,11:14]); % area

[m_a_old, se_a_old] = mean_SEM(area_z_t_avg(old,:)');

[m_a_young, se_a_young] = mean_SEM(area_z_t_avg(young,:)');



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%DIST%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ax3 = subplot(2,9,[1:4,10:13]); % area

% plot(1:n_intersect,area_z_t_avg(old,:),'o','MarkerFaceColor',c_old,'Color',c_old);hold on;

% plot(1:n_intersect,area_z_t_avg(young,:),'o','MarkerFaceColor',c_young,'Color',c_young) 

% scatter(1:n_intersect,m_a_old,[],c_o_map(c_idx,:),'fill')

bo = boundedline(flip(m_a_old),1:n_intersect,flip(se_a_old),'orientation', 'horiz','alpha','cmap',c_old);

bo.Color = c_old;
@@ -262,15 +295,25 @@ by = boundedline(flip(m_a_young),1:n_intersect,flip(se_a_young),'orientation', ' 
by.Color = c_young;



ylabel 'z-Score of intersect [\sigma]'

xlabel 'Area at intersect [u]'

xlabel 'Area at intersect [u^2]'

yticks([1  size(area_z_t_avg,2)])

yticklabels({'min Z','max Z'})

hline(n_intersect-i+1,'--k')

hl_o = hline([n_intersect-i_o(1)+1,n_intersect-i_o(end)+1],'--');

% line([m_a_old(n_intersect-i_o(1)) m_a_old(n_intersect-i_o(end))],[n_intersect-i_o(1)+1 ,n_intersect-i_o(end)+1],'Color',c_old)

set(findobj(hl_o,'Type', 'line'),'Color',c_old)

hl_y = hline([n_intersect-i_y(1)+1,n_intersect-i_y(end)+1],'--');

set(findobj(hl_y,'Type', 'line'),'Color',c_young);

% line([m_a_young(n_intersect-i_y(1)) m_a_young(n_intersect-i_y(end))],[n_intersect-i_y(1)+1 ,n_intersect-i_y(end)+1],'Color',c_young)







legend ([bo,by],{'Old','Young'})



legend ([bo,by,hl_o(1),hl_y(1)],{'Old','Young','FMHW','FMHW'},'NumColumns',2,'Location','southeast')

legend boxoff 

axis tight



save_fig(gcf,PATHOUT_TOODIFF,['surf_area_overview_' FB{fb}],'fontsize',20,'figsize',[0 0 40 20]);

pic_r = 1/2;

save_fig(gcf,PATHOUT_TOODIFF,['surf_area_overview_' FB{fb}],'fontsize',10,'figsize',[0 0 22.7 pic_r*22.7]);



%% 3D topo only

figure
@@ -281,7 +324,7 @@ py = surf2patch(sy); 
colormap(ax1,c_y_map)

axis tight;hold on;

plot3(elocsX,elocsY,top_young,'.k','MarkerSize',15)

title 'Interpolated Topoplot Young'

% title 'Interpolated Topoplot Young'

zlabel 'z-Score [\sigma]'

set(gca,'YTickLabel',[]);

set(gca,'XTickLabel',[]);
@@ -297,7 +340,7 @@ mesh(Xgrid,Ygrid,Z_old) % interpolate 
colormap(ax2,c_o_map)

axis tight;hold on;

plot3(elocsX,elocsY,top_old,'.k','MarkerSize',15)

title 'Interpolated Topoplot Old'

% title 'Interpolated Topoplot Old'

zlabel 'z-Score [\sigma]'

set(gca,'YTickLabel',[]);

set(gca,'XTickLabel',[]);
@@ -307,7 +350,8 @@ stlwrite([PATHOUT_TOODIFF 'mod_old3D.stl'],Xgrid,Ygrid,Z_old*0.2); 
fv = stlread([PATHOUT_TOODIFF 'mod_old3D.stl']);

stlwrite([PATHOUT_TOODIFF 'mod_old3D_' FB{fb} '.stl'],fv,'FaceColor',c_old); 



save_fig(gcf,PATHOUT_TOODIFF,['3D_topo_plain_' FB{fb}],'fontsize',18,'figsize',[0 0 12 20]);

pic_r = 2/1;

save_fig(gcf,PATHOUT_TOODIFF,['3D_topo_plain_' FB{fb}],'fontsize',20,'figsize',[0 0 22.7 pic_r*22.7]);



%% Disttribution params

if exist ('area_z_t_avg') ~= 1; load([PATHOUT_TOODIFF 'fun_surf_area.mat']);end
@@ -316,7 +360,14 @@ p = numSubplots(size(area_z_t_avg,2)); % for subplot config 
x_val = linspace(min(min(z_top)),max(max(z_top)),length(area_z_t_avg));

kurt = round(kurtosis(area_z_t_avg,[],1),1);

skew = round(skewness(area_z_t_avg,[],1),1);

c_idx = 56:-1:9;



%age diff

kurt_o = round(kurtosis(area_z_t_avg(old,i_o),[],1),1);

skew_o = round(skewness(area_z_t_avg(old,i_o),[],1),1);

kurt_y = round(kurtosis(area_z_t_avg(young,i_y),[],1),1);

skew_y = round(skewness(area_z_t_avg(young,i_y),[],1),1);



c_idx = ceil(linspace(56,9,n_intersect));



figure

for i = 1:size(area_z_t_avg,2)
@@ -338,12 +389,89 @@ for i = 1:size(area_z_t_avg,2) 
    hold on

    

    %add title & legend info

    legend([h_o(2) h_y(2)],{'Old','Young'},'Location','northeast');

    legend([h_o(2) h_y(2)],{'Old','Young'},'Location','best');

    legend boxoff

    title (['Kurt: ' num2str(kurt(i)) ' // Skew: ' num2str(skew(i))]);

    title (['K: ' num2str(kurt(i)) ' / S: ' num2str(skew(i))]);



end



save_fig(gcf,PATHOUT_TOODIFF,['dist_surf_area_' FB{fb}],'fontsize',8,'figsize',[0 0 50 40]);

pic_r = 1/1.3;

save_fig(gcf,PATHOUT_TOODIFF,['dist_surf_area_' FB{fb}],'fontsize',8,'figsize',[0 0 22.7 pic_r*22.7]);

save([PATHOUT_TOODIFF 'dist_param_' FB{fb} '.mat'],'kurt_o','skew_o','kurt_y','skew_y')



end



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

%

%                               Archive

%

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



%% calculate surface area  of topoplot for each trial

%{

surf_old = {};

surf_young = {};



for s = find(young)

    disp(['SUBJ: ' num2str(SUBJ(s))])

   for t = MI_trials

       

       topo = zscore((m_rERD_trial(s,:,t)));

       surf_st = griddata(elocsX,elocsY,topo,Xgrid,Ygrid,'cubic');

       top_thresh = min(min(topo))*0.5;

       surf_st(surf_st>top_thresh) = NaN;

       surf_young{end+1} = SurfArea(surf_st,Xgrid,Ygrid);

       

   end

end



for s = find(old)

    disp(['SUBJ: ' num2str(SUBJ(s))])

   for t = MI_trials

       

       topo = zscore((m_rERD_trial(s,:,t)));

       surf_st = griddata(elocsX,elocsY,topo,Xgrid,Ygrid,'cubic');

       top_thresh = min(min(topo))*0.5;

       surf_st(surf_st>top_thresh) = NaN;

       surf_old{end+1} = SurfArea(surf_st,Xgrid,Ygrid);

       

   end

end

%}



%% Calculate area of intersection



%{



gd(1).Z= griddata(elocsX,elocsY,top_young,Xgrid,Ygrid,'cubic'); %young

gd(2).Z = griddata(elocsX,elocsY,top_old,Xgrid,Ygrid,'cubic'); %old

group = {'Young','Old'};



plane_sect = linspace(0,-1.2,round(spat_freq));



for g = 1:length(gd)

    

    figure

    for i = 1:length(plane_sect)

        subplot(2,3,i)

        mesh(Xgrid,Ygrid,gd(g).Z) % interpolate

        hold on

        surf(Xgrid,Ygrid,plane_sect(i)*ones(size(gd(g).Z)),'FaceColor','k','FaceAlpha',0.08,'EdgeColor','none')

        hold on



        % Take the difference between the two surface heights and find the contour

        % where that surface is zero.

        zdiff = gd(g).Z - (plane_sect(i)*ones(size(gd(g).Z)));

        C = contour(Xgrid, Ygrid, zdiff,[0,0]);



        [a, ca] = cont_area(C,Xgrid,Ygrid,gd(g).Z);

        % Visualize the line.

        cl = line(ca.xL, ca.yL, ca.zL, 'Color', 'k', 'LineWidth', 3);

        title(['Area : ' num2str(a) ' at z = ' num2str(plane_sect(i))])

    end

    sgtitle ([group{g} ' participants'])

    save_fig(gcf,PATHOUT_TOODIFF,['SurfArea_' FB{fb} '_' group{g}],'figtype','.fig')



end



%}