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.6pi_std3d1.5_std2d1.5A.mat")
    load('field_elong_stats_TL.mat')
    neuron_end = 47;
    neurons = [ 1 3 9 15 16 17 20 23 24 25 27 29 34 35 36 40 41 42 43 44 45 47 48 49 50];
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*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

%%