% load your file

% variable name in the .mat file for WT = ["tf_wt_global", "tf_wt_delta",

% "tf_wt_theta", "tf_wt_alpha", "tf_wt_beta", "tf_wt_gamma"," tf_wt_HFO"]

% variable name in the .mat file for KO = ["tf_ko_global", "tf_ko_delta",

% "tf_ko_theta", "tf_ko_alpha", "tf_ko_beta", "tf_ko_gamma"," tf_ko_HFO"]

% example

load("f04_D_001_001_Spatial-power-spectrum-global-WT.mat");

signal = tf_wt_global;

% obtain montage and boudaries for the plot

[boundaries_func, montage_func] = get_boundary();

%run topoplot function

topoplot_visualisation(boundaries_func, montage_func, signal)

%% function to obtain boundaries and montage to pass as argument in the topoplot function

function [boundaries, montage] = get_boundary()

boundaries = [-4.400, 0.030; -4.180, 0.609; -3.960, 1.148; -3.740, 1.646; -3.520, 2.105; -3.300, 2.525; -3.080, 2.908; -2.860, 3.255;...

    -2.640, 3.566; -2.420, 3.843; -2.200, 4.086; -1.980, 4.298; -1.760, 4.4799; -1.540, 4.6321; -1.320, 4.7567; -1.100, 4.8553;...

    -0.880, 4.9298; -0.660, 4.9822; -0.440, 5.0150; -0.220, 5.0312; 0, 5.035; 0.220, 5.0312; 0.440, 5.0150; 0.660, 4.9822;...

    0.880, 4.9298; 1.100, 4.8553; 1.320, 4.7567; 1.540, 4.6321; 1.760, 4.4799; 1.980, 4.2986; 2.200, 4.0867; 2.420, 3.8430;...

    2.640, 3.5662; 2.860, 3.2551; 3.080, 2.9087; 3.300, 2.5258; 3.520, 2.1054; 3.740, 1.6466; 3.960, 1.1484; 4.180, 0.6099;...

    4.400, 0.0302; 4.400, 0.0302; 4.467, -0.1597; 4.5268, -0.3497; 4.5799, -0.5397; 4.6266, -0.7297; 4.6673, -0.9197; 4.7025, -1.1097;...

    4.7326, -1.2997; 4.7579, -1.4897; 4.7789, -1.6797; 4.7960, -1.8697; 4.8095, -2.0597; 4.8199, -2.2497; 4.8277, -2.4397; 4.8331, -2.6297;...

    4.8366, -2.8197; 4.8387, -3.0097; 4.8396, -3.1997; 4.8399, -3.3897; 4.8384, -3.5797; 4.8177, -3.7697; 4.7776, -3.9597; 4.7237, -4.1497;...

    4.6620, -4.3397; 4.5958, -4.5297; 4.5021, -4.7197; 4.400, -4.8937; 4.1800, -5.1191; 3.9600, -5.3285; 3.7400, -5.5223; 3.5200, -5.7007;...

    3.3000, -5.8642; 3.0800, -6.0131; 2.8600, -6.1478; 2.6400, -6.2688; 2.4200, -6.3764; 2.2000, -6.4712; 1.9800, -6.5536; 1.7600, -6.6241;...

    1.5400, -6.6833; 1.3200, -6.7317; 1.1000, -6.7701; 0.8800, -6.7991; 0.6600, -6.8194; 0.4400, -6.8322; 0.2200, -6.8385; 0, -6.840;...

    -0.220, -6.8385; -0.440, -6.8322; -0.660, -6.8194; -0.880, -6.7991; -1.100, -6.7701; -1.320, -6.7317; -1.540, -6.6833; -1.760, -6.6241;...

    -1.980, -6.5536; -2.200, -6.4712; -2.420, -6.3764; -2.640, -6.2688; -2.860, -6.1478; -3.080, -6.0131; -3.300, -5.8642; -3.520, -5.7007;...

    -3.740, -5.5223; -3.960, -5.3285; -4.180, -5.1191; -4.400, -4.89370; -4.5021, -4.7197; -4.5958, -4.5297; -4.6620, -4.3397; -4.7237, -4.1497;...

    -4.7776, -3.9597; -4.8177, -3.7697; -4.8384, -3.5797; -4.8399, -3.3897; -4.8397, -3.1997; -4.8387, -3.0097; -4.8367, -2.8197; -4.8331, -2.6297;...

    -4.8277, -2.4397; -4.8200, -2.2497; -4.8095, -2.0597; -4.7960, -1.8697; -4.7789, -1.6797; -4.7579, -1.4897; -4.7326, -1.2997; -4.7025, -1.1097;...

    -4.6673, -0.9197; -4.6266, -0.7297; -4.5799, -0.5397; -4.5268, -0.3497; -4.4670, -0.1597; -4.4000, 0.03025];

xmax = max(boundaries(:,1)); xmin = min(boundaries(:,1));

ymax = max(boundaries(:,2)); ymin = min(boundaries(:,2));

montage = [4.05, -4.14; 2.24, -4.14; 1.0, -4.14; 4.13, -3.04; 2.88, -3.04; 1.13, -3.04; 4.05, -1.96; 2.88, -1.96; 1.12, -1.96; 3.5, -0.48;...

    2.12, -0.48; 1.93, 1.04; 0.55, 1.04; 1.5, 2.3; 0.5, 2.3;...

    -0.5, 2.3; -1.5, 2.3; -0.55, 1.04; -1.93, 1.04; -2.12, -0.48; -3.5, -0.48; -1.12, -1.96; -2.88, -1.96; -4.05, -1.96; -1.13, -3.04;...

    -2.88, -3.04; -4.13, -3.04; -1.0, -4.10; -2.24, -4.14; -4.05, -4.14];

x = montage(:,1); y = montage(:,2);

end

%% Function to make the plot

function topoplot_visualisation(boundaries, montage, data)

%% preparing the important variables

GRID_SCALE = 60;

SHADING = 'interp';

CONTOURNUM = 30;

signal = data;

xmax = max(boundaries(:,1)); xmin = min(boundaries(:,1));

ymax = max(boundaries(:,2)); ymin = min(boundaries(:,2));

x = montage(:,1); y = montage(:,2);

xi   = linspace(xmin,xmax,GRID_SCALE);   % x-axis description (row vector)

yi   = linspace(ymin,ymax,GRID_SCALE);   % y-axis description (row vector)

[Xi,Yi,Zi] = griddata(x,y,signal,xi,yi','v4');

%% plot

cla  % clear current axis

hold on

h = gca; % uses current axes

AXHEADFAC = 1.05;     % do not leave room for external ears if head cartoon

set(gca,'Xlim',[xmin xmax]*AXHEADFAC,'Ylim',[ymin ymax]*AXHEADFAC);

unsh = (GRID_SCALE+1)/GRID_SCALE; % un-shrink the effects of 'interp' SHADING

handle = surface(Xi*unsh,Yi*unsh,zeros(size(Zi)),Zi,'EdgeColor','none','FaceColor',SHADING);

set(gca, "Clim", [-10 20])

colormap("jet")

[C,h] = contour(Xi,Yi,Zi,CONTOURNUM,'k','hittest','off', 'ShowText','off',"LabelFormat","%0.1f");

pgonX = [xmin*unsh xmax*unsh xmax*unsh xmin*unsh nan boundaries(:,1)'];

pgonY = [ymin*unsh ymin*unsh ymax*unsh ymax*unsh nan boundaries(:,2)'];

pgon = polyshape(pgonX,pgonY);

patch_coordinate = pgon.Vertices;

[F,V] = poly2fv(patch_coordinate(:,1), patch_coordinate(:,2));

patchH = patch('Faces', F, 'Vertices', V);

patchH.EdgeColor = "w";

patchH.FaceColor = "w";

plot(boundaries(:,1), boundaries(:,2), "k-")

hold on

plot(montage(:,1), montage(:,2), "ro", MarkerFaceColor="k", MarkerSize=10)

ref = [2.20, -6.0];

plot(ref(1), ref(2), "ks", MarkerSize = 20,MarkerFaceColor="k" )

bregma = [0,0];

plot(bregma(1), bregma(2), "k+", MarkerSize=17,LineWidth=4 )

xlabel( "ML coordinate (mm)");

ylabel( "AP coordinate(mm)");

axis equal

end