clc
clear all
close all
%%
%inputs
%CHANGES TO MAKE-values of field length are in voxels ,have to convert
% firr3dprops_file="props_TL_t20p3_b1.8.mat";
%
% load(firr3dprops_file);
% rng('default');
vl=0.25;%voxel length
firr_thresh=0.2;%0.2
tilted=0;%check if tilted
%calculate stats from which neuron to which neuron?
neuron_start=1;
if tilted
load("Tilted_data_props.mat")
load("field_elongstats_tilted.mat")
neuron_end = 50;
neurons = [];
else
load("Aligned_data_props_b0.8pi_std3d1.5_std2d1.5A.mat")
load('field_elong_stats_TL.mat')
neuron_end = 47;
neurons = [ 1 5 8 9 10 11 14 18 20 25 26 27 33 36 39];
end
minv=50;%minimum voxels to be considered as place field
%CHOOSE PLOTS TO PLOT
identify_placef=1;%set this to 1 to identify place fields ,check min voxels
calc_orientation_elongation=1;
plot_orientation=1; %plot orientation on sphere
plot_elongation=1;
plot_centroids=1;
plot_binary=1;%binary morphology
fields_per_cell=0;%plot fields per cell
plot_pf=0;%plot place field firing rate map
calc_proportion=0; %calculate fields parallel to x,y,z
plot_flength=0;
plotfirr=0;
neuron_firr=22;%neuron for which to plot firr
rotm = [[ 0.70710678, -0.57922797, -0.40557979],
[ 0. , 0.57357644, -0.81915204],
[ 0.70710678, 0.57922797, 0.40557979]];
XYZ = [1 0 0; 0 1 0; 0 0 1];
ABC = rotm*XYZ;
axes_cords = [XYZ; ABC];
% set(gca,'FontSize',30);
%% %%%%%%% IDENTIFY PLACE CELLS %%%%%%%%%%%%
x1=1.125:0.25:5.875 ;
y1=1.125:0.25:5.875;
z1=1.125:0.25:5.875;
[X,Y,Z]=meshgrid(x1, y1 ,z1);
voxel_coords=[Y(:) X(:) Z(:)];%voxel coords increase along x then y then z
orientation=double.empty;
orientationplot=double.empty;
elongation=double.empty;
yfl=double.empty;
xfl=double.empty;
zfl=double.empty;
centroids=double.empty;
%%identify place fields from all neurons and store in pf structure
if identify_placef
pf=struct;
pf_counter=1;
for neuron=neurons
propstable=FIRR_3dprops(neuron).props; %get 3d properties table of the neuron
firrtable=FIRR_3dprops(neuron).firr;
for i=1:height(propstable)
vol=propstable{i,"Volume"};
if vol>minv %if volume of cluster is greater than minvoxels
pf(pf_counter).props=propstable(i,:);
pf(pf_counter).neuron=neuron;
pf(pf_counter).firr=firrtable;
pf_counter=pf_counter+1;
end
end
end
end
%% %%%%%%% PLOT CENTROIDS FIG 4 %%%%%%%%%%%%
if plot_centroids
%%%%%%%%%%%%% FIG 4 a1) Aligned or FIG 4 b1) Tilted %%%%%%%%%%%%%%
set(0, 'DefaultFigureRenderer', 'painters');
for i=1:size(pf,2)
centroid=pf(i).props.Centroid;
centroids=[centroids;centroid];
centroids=centroids;
end
% centroids = rotm * centroids';
centroids_scaled=(5/20)*(centroids)+1;
if tilted == 1
centroids_scaled1 = rotm * centroids_scaled';
centroids_scaled1 = centroids_scaled1';
end
%%hold off
figure
if tilted == 1
V1 = [1 1 1];
V2 = [6 1 1];
V3 = [6 1 6];
V4 = [1 1 6];
V5 = [1 6 6];
V6 = [1 6 1];
V7 = [1 6 6];
V8 = [6 6 6];
V9 = [6 6 1];
V10 = [1 6 1];
V11 = [1 1 1];
V12 = [1 1 6];
V13 = [6 1 6];
V14 = [6 6 6];
V15 = [6 6 1];
V16 = [6 1 1];
C = [V1;V2;V3;V4;V5;V6;V7;V8;V9;V10;V11;V12;V13;V14;V15;V16];
E = rotm*C';
E = E';
figure
plot3(E(:,1),E(:,2),E(:,3),'black')
hold on
end
if tilted
scatter3(centroids_scaled1(:,1) ,centroids_scaled1(:,2), centroids_scaled1(:,3),30, centroids_scaled1(:,3),"filled")
title("Centroids of place fields- Tilted lattice")
hold on
else
scatter3(centroids_scaled(:,1) ,centroids_scaled(:,2), centroids_scaled(:,3),30, centroids_scaled(:,3),"filled")
title("Centroids of place fields- Aligned lattice")
hold on
end
xlabel('x','Fontsize',15)
ylabel('y','Fontsize',15)
zlabel('z','Fontsize',15)
if tilted
medx=mean(centroids_scaled1(:,1));
medy=mean(centroids_scaled1(:,2));
medz=mean(centroids_scaled1(:,3));
else
medx=mean(centroids_scaled(:,1));
medy=mean(centroids_scaled(:,2));
medz=mean(centroids_scaled(:,3));
end
if tilted == 0
hold on
%
scatter3(medx,medy,medz,100,'filled','red')
hold on
end
if tilted == 1
hold on
scatter3(medx,medy,medz,100,'filled','red')% scatter3(medx,3.5,3.5,100,'filled','green')
end
x2=3.5*ones(length(x1),1);% x2=3.5*ones(length(x1),1);
y2=3.5*ones(length(y1),1);
z2=3.5*ones(length(z1),1);
if tilted ==0
plot3(x1,y2,z2,'Color','black');
hold on
plot3(x2,y1,z2,'Color','black');
hold on
plot3(x2,y2,z1,'Color','black');
hold on
end
hold off
%%%%%%%%%%% Fig 4a2) Aligned or Fig 4b2) Tilted %%%%%%%%%%%
if tilted
cc = inv(rotm)*centroids_scaled';
cc = cc';
end
cc = centroids_scaled;
center = [3.5, 3.5, 3.5];
cc = cc - center;
bw = 0.7;
figure;
xi = -5:0.1:5;
[~,~,~,f1] = histDENSITY(cc(:,1),xi,bw);
[~,~,~,f2] = histDENSITY(cc(:,2),xi,bw);
[~,~,~,f3] = histDENSITY(cc(:,3),xi,bw);
leng_line = max([max(f1) max(f2) max(f3)]);
p1 = plot(xi(:),f1,'LineWidth',1.5, 'color', 'red'); hold on;
p2 = plot(xi(:),f2,'LineWidth',1.5, 'Color', 'green');
p3 = plot(xi(:),f3,'LineWidth',1.5, 'Color', 'blue');
line([nanmedian(cc(:,1)),nanmedian(cc(:,1))], [0 leng_line], 'color', 'red');
line([nanmedian(cc(:,2)),nanmedian(cc(:,2))], [0 leng_line], 'color', 'green');
line([nanmedian(cc(:,3)),nanmedian(cc(:,3))], [0 leng_line], 'color', 'blue');
legend('X', 'Y', 'Z');
xlabel('Distance from center'); ylabel('Proportion of fields');
%%%%%%% Fig 4a3) Aligned or Fig 4b3) Tilted %%%%%%%%
xi = -0.5:.001:0.5;
med_rand_x =[]; med_rand_y =[]; med_rand_z=[];
for i = 1:1000
med_rand_x = [med_rand_x nanmedian((5*rand(1,length(cc)))-2.5)];
med_rand_y = [med_rand_y nanmedian((5*rand(1,length(cc)))-2.5)];
med_rand_z = [med_rand_z nanmedian((5*rand(1,length(cc)))-2.5)];
end
med_rand = [med_rand_x; med_rand_y; med_rand_z]';
colors = ['r'; 'g'; 'b'];
fig = figure;
for j =1:3
subplot(3,1,j);
temp = med_rand(:,j);
[~,~,~,f] = histDENSITY(temp,xi,0.01);
p1 = plot(xi(:),f,'LineWidth',1.5, 'color', colors(j)); %hold on;
line([nanmedian(cc(:,j)),nanmedian(cc(:,j))], [0 max(f)], 'color', 'black', 'LineWidth', 1.5);
line([prctile(temp, 97.5), prctile(temp, 97.5)], [0 max(f)], 'color', 'red', 'LineWidth', 1.5);
line([prctile(temp, 2.5), prctile(temp, 2.5)], [0 max(f)], 'color', 'red','LineWidth', 1.5);
xlabel("Median field position");
ylabel("Proportion of fields");
end
end
%% %%%%%%% CALCULATE ORIENTATION AND ELONGATION %%%%%%%%%%%%
%orientation ,elongation
if calc_orientation_elongation
orientation=double.empty;
for i=1:length(pf)%calculate stats
% i
%calculate orientation
b=pf(i).props.ConvexHull{1,1};
a=pca(b);
if tilted
a = rotm * a;
orientation=[orientation a(:,1)];
orientationplot=[orientationplot a(:,1) -a(:,1)];
else
orientation=[orientation a(:,1)];
orientationplot=[orientationplot a(:,1) -a(:,1)];
end
%calculate elongation
t=pf(i).props;
arr=t.PrincipalAxisLength;
elong=2*arr(1)/(arr(2)+arr(3));
elongation=[elongation ;elong];
end
end
extent_field = [];
for neuron=neurons
propstable=FIRR_3dprops(neuron).props;
for i=1:height(propstable)
vol=propstable{i,"Volume"};
if vol>minv
extent_field = [extent_field; propstable.BoundingBox(i,4:6)];
end
end
end
xi = 0:1:25;
[~,~,~,f1] = histDENSITY(extent_field(:,1),xi,1);
[~,~,~,f2] = histDENSITY(extent_field(:,2),xi,1);
[~,~,~,f3] = histDENSITY(extent_field(:,3),xi,1);
figure;
p1 = plot(xi(:), f1, 'color', 'r');hold on;
p2 = plot(xi(:), f2, 'color', 'g');
p3 = plot(xi(:), f3, 'color', 'b');
%% %%%%%%% BINARY MORPHOLOGY FIG 7 %%%%%%%%%%%%
if plot_binary
total=0;
%create table like array that has SE length and correponding
%proportions
Px=double.empty;
Py=double.empty;
Pz=double.empty;
% Px(1,:)=[3:2:19];
% Py(1,:)=[3:2:19];
% Pz(1,:)=[3:2:19];
Px(1,:)=1:1:19;
Py(1,:)=1:1:19;
Pz(1,:)=1:1:19;
Px(2,:)=zeros(1,19);
Py(2,:)=zeros(1,19);
Pz(2,:)=zeros(1,19);
ind=1 ; %index simply represents index of array l which contains structuring element length ;
for l=1:1:19 %
ex = strel('cuboid',[l 1 1]); %elementx
ey = strel('cuboid',[1 l 1]);
ez = strel('cuboid',[1 1 l]);
total=0;
for neuron=neurons %binarize rate map of each neuron using 10% of maximum as thresh and make binary volumetric image
%neuron
bw=zeros(20,20,20);
firr=FIRR_3dprops(neuron).firr;
nanIndices = isnan(firr);
infIndices = isinf(firr);
firr(nanIndices | infIndices) = 0;
counter=1;
threshold=0.1*max(firr);
for k=1:20 % this loop is to create binary volumetric image for a neuron (with threshold)
j=1;
for j=1:20
i=1;
for i=1:20
if firr(counter)>threshold
bw(i,j,k)=1;
end
counter=counter+1;
end
end
end
%pause;
Bx=imerode(bw,ex);
By=imerode(bw,ey);
Bz=imerode(bw,ez);
total1=sum(Bx,"all")+sum(By,"all")+sum(Bz,"all");%total number of voxels with a particular connectivity for a particular neuron's firing rate
total=total+total1; %total ,after complete iteration of all neurons,is the sum of all voxels remaining for that particular SE length
if neuron==neuron_start
Px(2,ind)=sum(Bx,"all");
Py(2,ind)=sum(By,"all");
Pz(2,ind)=sum(Bz,"all");
else
Px(2,ind)=Px(2,ind)+sum(Bx,"all");%sum of remaining voxels of every iterated neuron along x axis
Py(2,ind)=Py(2,ind)+sum(By,"all");
Pz(2,ind)=Pz(2,ind)+sum(Bz,"all");
end
end%neuron is inner loop
Px(2,ind)=Px(2,ind)/total;
Py(2,ind)=Py(2,ind)/total;
Pz(2,ind)=Pz(2,ind)/total;
ind=ind+1;
end %strel is outer loop
figure
plot(Px(1,:),Px(2,:), 'Color', 'r', 'LineWidth', 1.5)
hold on
plot(Py(1,:),Py(2,:), 'Color', 'g', 'LineWidth', 1.5)
hold on
plot(Pz(1,:),Pz(2,:), 'Color', 'b', 'LineWidth', 1.5)
ylabel('Proportion of voxels','FontSize',15)
xlabel('Structuring element length','FontSize',15)
if tilted
legend('A','B','C')
else
legend('X','Y','Z')
end
xlim([0 20]); ylim([0 1]);
title('Binary morphology','Fontsize',16)
end
%% %%%%%%% PLOT ELONGATION FIG 5 %%%%%%%%%%%%
if plot_elongation
%%%%%%%%%%% Fig 5a1) Aligned or Fig 5b1) Tilted %%%%%%%%%%%
hold off
figure
scatter(ones(size(elongation)),elongation,'filled','blue','jitter','on','jitterAmount',0.05);
hold on
yline(1,'Color','black')
hold on
box1=boxplot(elongation','Whisker',5); hold on;
set(box1,'Color','black');
ylim([0,5]); yticks([1, 2, 3, 4, 5]); xticks([1]);
sph_fields = num2str(round((sum(field_elong_stats(:,5))*100)/length(field_elong_stats),2)) + " %";
t = text(0.85,0.5,sph_fields); t.FontSize = 15;
t=title("Elongation of place fields-Aligned",'FontSize',16);
set(t);
set(gca,'linewidth',2);
ylabel("Field elongation",'FontSize',15);
hold off
%%%%%%%%%%% Fig 5a2) Aligned or Fig 5b2) Tilted %%%%%%%%%%%
end
%% %%%%%%% PLOT ORIENTATION FIG 6 %%%%%%%%%%%%
if plot_orientation
if tilted==0
figure
x1=-1.2:0.1:1.2;
y1=-1.2:0.1:1.2;
z1=-1.2:0.1:1.2;
x2=zeros(25,1);
y2=zeros(25,1);
z2=zeros(25,1);
plot3(x1,y2,z2,'Color','black'); hold on
plot3(x2,y1,z2,'Color','black'); hold on
plot3(x2,y2,z1,'Color','black'); hold on
V=double.empty;
V =orientationplot.';
[X,Y,Z,F2] = projectEIGS(V,256,10);
mx = nanmax(F2(:));
F2m = (F2 ./ mx);
s = surf(X,Y,Z,F2m);
s.EdgeColor = 'k';
s.EdgeAlpha = 0.0;
colormap(gca,'jet');
axis on vis3d;
rotate3d on;
title("Place field orientation",'FontSize',16)
xlabel('Y'); ylabel('X'); zlabel('Z'); hold off
elseif tilted==1
scale=-1.5:1:1.5; %scale of reference axes (tilted)
%axis A
i=1;
j=0;
k=0;
i=i*scale;
j=j*scale;
k=k*scale;
x_axis = [i; j; k];
x_axis = rotm * x_axis;
i = x_axis(1,:);
j = x_axis(2,:);
k = x_axis(3,:);
figure()
plot3(i,j,k);
hold on
%axis B
i = 0;
j = 1;
k = 0;
i=i*scale;
j=j*scale;
k=k*scale;
y_axis = [i; j; k];
y_axis = rotm * y_axis;
i = y_axis(1,:);
j = y_axis(2,:);
k = y_axis(3,:);
plot3(i,j,k);
hold on
%axis C
i = 0;
j = 0;
k = 1;
i=i*scale;
j=j*scale;
k=k*scale;
z_axis = [i; j; k];
z_axis = rotm * z_axis;
i = z_axis(1,:);
j = z_axis(2,:);
k = z_axis(3,:);
plot3(i,j,k);
V=double.empty;
V =orientationplot';
[X,Y,Z,F2] = projectEIGS(V,256,10);
mx = nanmax(F2(:));
F2m = (F2 ./ mx);
s = surf(X,Y,Z,F2m);
s.EdgeColor = 'k';
s.EdgeAlpha = 0.0;
colormap(gca,'jet');
axis on vis3d;
rotate3d on;
title("Place field orientation")
xlabel('Y')
ylabel('X')
zlabel('Z')
end
%%Statistical analysis
theta_bool = [];
theta_bool1 = [];
%collect all primary eigenvectors
ev = [];
for neuron=neurons
Vol = FIRR_3dprops(neuron).props.Volume;
for i = 1:length(Vol)
if (Vol(i)>50)
temp2 = FIRR_3dprops(neuron).props.EigenVectors(i);
ev = [ev temp2{1,1}(:,1)];
end
end
end
%check if eigenvectors are in vicinity of axes
if tilted
ev_tilted = rotm*ev;
ev_tilted = ev_tilted';
for i = 1:length(ev_tilted)
ev_temp = ev_tilted(i,:);
theta = acosd(abs(ev_temp*XYZ'));
temp = theta<30;
theta_bool = [theta_bool; temp];
end
ev= ev';
for i = 1:length(ev)
ev_temp2 = ev(i,:);
theta = acosd(abs(ev_temp2*XYZ'));
temp = theta<30;
theta_bool1 = [theta_bool1; temp];
end
else
ev= ev';
for i = 1:length(ev)
ev_temp = ev(i,:);
theta = acosd(abs(ev_temp*axes_cords'));
temp = theta<30;
theta_bool1 = [theta_bool1; temp];
end
end
theta_bool = [theta_bool theta_bool1];
ratio_ev = sum(theta_bool,1)/length(theta_bool);
xyz_abc = sum(ratio_ev(1:3))/sum(ratio_ev(4:6));
%%error bars and confidence intervals
% randomly choose pfs from uniform distribution and create error bars
shuffled_ratio = [];
if tilted
for i=1:1000
rindx = randi(length(ev), 1,length(ev));
shuffl_pf_t = ev_tilted(rindx,:);
theta_dat = []; theta_dat1 = [];
for j = 1:length(shuffl_pf_t)
temp = shuffl_pf_t(j,:);
theta = acosd(abs(temp*XYZ'));
temp = theta<30;
theta_dat = [theta_dat; temp];
end
shuffl_pf = ev(rindx,:);
for j = 1:length(shuffl_pf)
ev_temp2 = shuffl_pf(j,:);
theta = acosd(abs(ev_temp2*XYZ'));
temp = theta<30;
theta_dat1 = [theta_dat1; temp];
end
theta_dat2 = [theta_dat theta_dat1];
shuffled_ratio = [shuffled_ratio; sum(theta_dat2,1)/length(theta_bool)];
end
else
for i=1:1000
rindx = randi(length(ev), 1,length(ev));
shuffl_pf = ev(rindx,:);
theta_dat2 = [];
for j = 1:length(shuffl_pf)
temp = shuffl_pf(j,:);
theta = acosd(abs(temp*axes_cords'));
temp = theta<30;
theta_dat2 = [theta_dat2; temp];
end
shuffled_ratio = [shuffled_ratio; sum(theta_dat2,1)/length(theta_bool)];
end
end
%%%%%%%%%%% Fig 6a1) Aligned or Fig 6b1) Tilted %%%%%%%%%%%
%%with random points
% randomly generate points to know the by chance line and also the XYZ/ABC
% distribution
figure;
total = [];
for i = 1:1000
TH = 2*pi*rand(1,1000);
PH = asin(-1+2*rand(1,1000));
[X,Y,Z] = sph2cart(TH,PH,1);
%for j = 1:1000
ev_rand = [X; Y; Z];
ev_rand = ev_rand';
theta = acosd(abs(ev_rand*axes_cords'));
temp = theta<30;
total = [total; sum(temp, 1)/length(temp)];
end
rand_shuffl_ratio = sum(total(:,1:3), 2)./sum(total(:,4:6), 2);
% plotting distribution
xi = 0:0.01:2;
dist = sum(total(:,1:3), 2)./sum(total(:,4:6), 2);
llev = prctile(rand_shuffl_ratio, 1);
hlev = prctile(rand_shuffl_ratio, 99);
[~, ~, ~, fn] = histDENSITY(dist, xi, 0.09);
% p = plot(xi(:),fn,'LineWidth',1.5, 'color', 'black');
area(xi,fn,'FaceColor', [0.9 0.9 0.9],'LineWidth',1.5);
line([xyz_abc xyz_abc], [min(fn), max(fn)], 'color', 'b', 'LineWidth', 1.5);
line([llev llev], [min(fn), max(fn)], 'color', 'r', 'LineWidth', 1.5);
line([hlev hlev], [min(fn), max(fn)], 'color', 'r', 'LineWidth', 1.5);
xlabel('XYZ/ABC'); ylabel('Probability')
%%%%%%%%%%% Fig 6a2) Aligned or Fig 6b1) Tilted %%%%%%%%%%%
%%plotting errorbar
figure;
prcr1 = prctile(total, 2.5); % percentiles for chance line
prcr2 = prctile(total, 50);
prcr3 = prctile(total, 97.5);
prcs1 = prctile(shuffled_ratio, 2.5); % percentile for error bar
prcs2 = prctile(shuffled_ratio, 97.5);
for i = 1:6
errorbar(i,ratio_ev(i),ratio_ev(i)-prcs1(i),ratio_ev(i)-prcs2(i), 'LineWidth', 1.0); hold on;
line([0.5 7], [prcr1(i) prcr1(i)], 'color', 'red', 'LineWidth', 0.5);
line([0.5 7], [prcr3(i) prcr3(i)], 'color', 'red', 'LineWidth', 0.5);
line([0.5 7], [prcr2(i) prcr2(i)], 'color', 'red', 'LineWidth', 1.5);
plot(i,ratio_ev(i), 'Marker','o', 'LineWidth', 1.5, 'MarkerSize', 15);
xlim([0.5,7]);
ylim([0,0.8]);
end
names = {'X'; 'Y'; 'Z'; 'A'; 'B'; 'C'};
set(gca,'xtick',1:6,'xticklabel',names)
xlabel('Maze Axis'); ylabel('Proportion of fields')
end
%% BINARY MORPHOLOGY ANOVA APPARENTLY??
%binary morphology ANOVA, creating a new loop coz there's a difference in
%the way ill be calculating proportion, here I'll take proportion for each
%neuron seperately and outer loop will be neurons
% if binary_anova
%
% binary_stat_table=[];
% Px_statb4=zeros(1,12);
% Py_statb4=zeros(1,12);
% Pz_statb4=zeros(1,12);
%
%
%
%
%
%
%
%
%
%
%
%
% for neuron=setdiff(neuron_start:neuron_end,avoid_neurons) %binarize rate map of each neuron using 10% of maximum as thresh and make binary volumetric image
% %neuron
%
%
%
% for l=8:1:19 %
% bw=zeros(20,20,20);
% firr=FIRR_3dprops(neuron).firr;
% counter=1;
% threshold=0.1*max(firr);
% for k=1:20 % this loop is to create binary volumetric image for a neuron (with threshold)
% j=1;
% for j=1:20
% i=1;
% for i=1:20
% if firr(counter)>threshold
% bw(i,j,k)=1;
% end
% counter=counter+1;
% end
% end
% end
%
% ex = strel('cuboid',[l 1 1]); %elementx
% ey = strel('cuboid',[1 l 1]);
% ez = strel('cuboid',[1 1 l]);
% total1=0;
%
%
% Bx=imerode(bw,ex);
% By=imerode(bw,ey);
% Bz=imerode(bw,ez);
%
% total1=sum(Bx,"all")+sum(By,"all")+sum(Bz,"all");%total number of voxels with a particular connectivity for a particular neuron's firing rate
%
%
%
% ind=l;
% Px_statb4(ind-7)=sum(Bx,"all");
% Py_statb4(ind-7)=sum(By,"all");
% Pz_statb4(ind-7)=sum(Bz,"all");
%
%
%
%
% Px_stat(ind-7)=Px_statb4(ind-7)/total1;
% Py_stat(ind-7)=Py_statb4(ind-7)/total1;
% Pz_stat(ind-7)=Pz_statb4(ind-7)/total1;
%
%
%
%
%
%
%
%
%
% end %l is inner loop for anova calculation
%
% binary_stat_singleneuron=[neuron Px_stat Py_stat Pz_stat];
% binary_stat_table=[binary_stat_table;binary_stat_singleneuron];
% binary_stat_table(isnan(binary_stat_table))=0;%NaN turned into 0s
% %
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
% end %neuron is outer loop for anova calculation
% % for i=1:1:58 %plot the means to check if they're good
% % binary_stat_table(33,i)=mean(binary_stat_table(:,i));
% %
% % end
% %
%
%
%
% end
%% %%%%%%% PLOT FIRING FIELDS %%%%%%%%%%%%
if plotfirr
firr_arr=FIRR_3dprops(neuron_firr).firr;
firr_arr(isnan(firr_arr))=0;
col=vals2colormap(firr_arr,'jet',[0.2 1]);
%set nans to zero
figure
plotcube_original([5 5 5],[1 1 1],0,[1 1 1])
hold on
for i=1:length(voxel_coords)%plot 3d firing rate voxel
P = [voxel_coords(i,1),voxel_coords(i,2),voxel_coords(i,3)] ; % your center point
L = [vl,vl,vl] ; % your cube dimensions
O = P-L/2 ;
if firr_arr(i)>firr_thresh*max(firr_arr)
plotcube(L,O,0.3,col(i,:)) %firr_col(i,:) % use function plotcube
hold on
end
end
colormap(jet)
title(["Firing rate of neuron ",num2str(neuron_firr)],'Fontsize',16)
xlabel("X",'Fontsize',16)
ylabel("Y",'Fontsize',16)
zlabel("Z",'Fontsize',16)
end
%% %%%%%%% I HAVE NO CLUE WHAT THIS SECTION DOES %%%%%%%%%%%%
if calc_proportion
orientation=[orientation(1,:);orientation(2,:);orientation(3,:)];
count_x=0
count_y=0
count_z=0
count_a=0
count_b=0
count_c=0
for n=1:len(pf)
i=[1 0 0]
j=[0 1 0]
k=[0 0 1]
a=[sin(deg2rad(35.26))*cos(deg2rad(60)) sin(deg2rad(35.26))*sin(deg2rad(60)) cos(deg2rad(35.26))];
b=[sin(deg2rad(35.26))*cos(deg2rad(-60)) sin(deg2rad(35.26))*sin(deg2rad(-60)) cos(deg2rad(35.26))];
c=[sin(deg2rad(35.26))*cos(deg2rad(180)) sin(deg2rad(35.26))*sin(deg2rad(180)) cos(deg2rad(35.26))];
if 0.95<dot(orientation(:,n),i)<1
count_x=count_x+1;
end
if 0.95<dot(orientation(:,n),j)<1
count_y=count_y+1;
end
if 0.95<dot(orientation(:,n),k)<1
count_z=count_z+1;
end
if 0.95<dot(orientation(:,n),a)<1
count_a=count_a+1;
end
if 0.95<dot(orientation(:,n),b)<1
count_b=count_b+1;
end
if 0.95<dot(orientation(:,n),c)<1
count_c=count_c+1;
end
end
end
%% %%%%%%% PLOT FIELD LENGTH %%%%%%%%%%%%
%field length
if plot_flength
for i=1:size(pf,2)
flz=pf(i).props.BoundingBox(6);
zfl=[zfl;flz];
flx=pf(i).props.BoundingBox(5);
xfl=[xfl;flx];
fly=pf(i).props.BoundingBox(4);
yfl=[yfl;fly];
end
figure
% histogram(zfl)
% hold on
% [H, P, h] = bootmode (zfl, 2, 1000,'gaussian')
histogram(xfl)
hold on
histogram(yfl)
title("Field lengths")
legend("X","Y")
end
%fields per cell
if fields_per_cell
count1=0
count2=0
count3=0
count4=0
count5=0
count6=0
count7=0
for neuron=1:400
neuron
counter=0
for i=1:size(pf,2)
if neuron==pf(i).neuron
counter=counter+1
end
end
if counter==1
count1=count1+1
end
if counter==2
count2=count2+1
end
if counter==3
count3=count3+1
end
if counter>3
count4=count4+1
end
end
y=[0 count1 count2 count3 count4 0]/400;
hold off
figure
plot([0 1 2 3 4 5],y);
title("Fields per cell-algined")
xlabel('No. of fields-4 implies >3')
ylabel('Proportion of fields')
hold off
end
%% %%%%%%% AGAIN NO CLUE %%%%%%%%%%%%
%place field firing rate
if plot_pf %plot place field
for m=1:length(pf)
firr_arr=pf(m).firr;
voxelid=pf(m).props.VoxelIdxList{1,1};
for i=1:length(voxel_coords) %looping over every voxel
if any(voxelid==i)
continue
else
firr_arr(i)=0;
end
end
firr_arr(isnan(firr_arr))=0;
col=vals2colormap(firr_arr,'jet',[0.2 1])
% firr_col=zeros(length(firr_arr),3)
% firr_arr(isnan(firr_arr))=0
% firr_col(:,1)=(1-firr_arr(:))/max(firr_arr)
% firr_col(:,2)=firr_col(:,1)
% firr_col(:,3)=firr_col(:,1)
hold off
figure(m);
title(['Placefield ',num2str(m)])
plotcube([5 5 5],[1 1 1],0,[1 1 1])
xlabel('x')
ylabel('y')
hold on
for i=1:size(voxel_coords,1)%plot 3d firing rate voxel
P = [voxel_coords(i,1),voxel_coords(i,2),voxel_coords(i,3)] ; % your center point
L = [vl,vl,vl] ; % your cube dimensions
O = P-L/2 ;
if firr_arr(i)>firr_thresh*max(firr_arr)
plotcube(L,O,0.4,col(i,:)) % use function plotcube
hold on
end
end
end
end
%%