%%
clear all; clc;
%%
traj_file='Trajectory_interpolated_tilted_lattice_no_diagAzra.csv';
traj_data=csvread(traj_file);
cent_file = 'Tilted_data_props_b0.8pi_std3d1.5_std2d1.5.mat';
load(cent_file);
spike_neu_file = 'encoded_TL_t20p3_b0.8pi_3d_testA.csv';
spike_dat = csvread(spike_neu_file);
% spike_file = 'G:\.shortcut-targets-by-id\1_HPNBtvegcfjF9AB6-jAJHtWwqwC4YCF\3d\3d_lattice_maze\Field_elongation\pf_spike';
% load(spike_file);
voxel_file = 'CC_voxel.mat';
load(voxel_file);

voxel_cord_f = 'voxel_coords.mat';
load(voxel_cord_f);

%%
vox_dat = CC_voxel(~cellfun('isempty',CC_voxel));
neurons = [5,6,7,10,14,18,25,26,33,36,39,41,46]; %With aligned phi case
spike_data = [];
for neu=neurons
    spike_data = [spike_data, spike_dat(:,neu)];
end



%%
% seperating the spikes for each place field
vl = 0.25;
hvl = vl/2;
pf_spike = {};
pf_traj = {};
for neuron = 1:numel(vox_dat)
% for neuron = 3:3
    neuron
    spike_neu = spike_data(:,neuron);
    %neuron
    for j = 1:numel(vox_dat{1,neuron})
        pf = vox_dat{1,neuron}{1,j};
        pf_sz = size(pf);
        if pf_sz(1) > 50
            spike_temp = []; traj_temp =[];
            for vox = pf'
                vx = voxel_coords(vox,1); vy = voxel_coords(vox,2); vz = voxel_coords(vox,3);
                for i = 1:(length(traj_data)-1)
                    px=traj_data(i,1); py=traj_data(i,2); pz=traj_data(i,3);
                    if abs(vx-px)<hvl & abs(vy-py)<hvl & abs(vz-pz)<hvl
                        spike_temp = [spike_temp, spike_neu(i)];
                        traj_temp = [traj_temp; traj_data(i)];
                    end
                end
            end
            pf_spike{1,numel(pf_spike)+1} = spike_temp;
            pf_traj{1,numel(pf_traj)+1} = traj_temp;
        end
    end
end

%%
% Centroids and Volumes of all the place fields
minv = 50;
Centroids = [];
Vols = [];
for neuron = 1:46
    neuron
    propstable=FIRR_3dprops(neuron).props;
    for i = 1:height(propstable)
        vol=propstable{i,"Volume"};
        if vol > minv
            convol=propstable{i,"ConvexVolume"};
            Vols = [Vols, convol];
            Centroids = [Centroids; propstable{i,"Centroid"}];
        end
    end
end

Centroids = Centroids*0.25 + 1;

%%
% create a sphere
shuffl_data{1,8100} = [];
n_shuffle = 100;
[x,y,z] = sphere;
shuffl_data = {};
for i = 1:length(Vols)
% for i = 1:1
    rad = 3*((Vols(i))/(pi))^(1/3);
    rad = rad*0.25;
    x_ = x*rad;
    y_ = y*rad;
    z_ = z*rad;
    c = Centroids(i,:);
    cx=c(1); cy=c(2); cz=c(3);
    inside = [];
    for i2 = 1:length(traj_data)
        i2;
        p = traj_data(i2,:);
        px=p(1); py=p(2); pz=p(3);
        d = ((cx-px)^2+(cy-py)^2+(cz-pz)^2)^0.5; 
        if d<rad
            inside = [inside; p];
        end
    end
    n_spike = length(pf_spike{1,i});
    for i3=1:n_shuffle
        rand_x = normrnd(cx, 1.5*rad, [1,n_spike]);
        rand_y = normrnd(cy, 1.5*rad, [1,n_spike]);
        rand_z = normrnd(cz, 1*rad, [1,n_spike]);

        rand_pts = [rand_x', rand_y', rand_z'];

        rand_spike_id = knnsearch(inside, rand_pts);
        rand_spike_traj = inside(rand_spike_id,:); 
        shuffl_data{1,((i-1)*n_shuffle) + i3} = rand_spike_traj;
    end
    
    i
end

%%
% data = inside;
% x=data(:,1);
% y=data(:,2);
% z=data(:,3);
% T=delaunay(x,y);
% trisurf(T,x,y,z);
% xy = [x,y];
% a = xy(T(:,2),:)-xy(T(:,1),:);
% b = xy(T(:,3),:)-xy(T(:,1),:);
% V = ((a(:,1).*b(:,2)-a(:,2).*b(:,1))' * sum(z(T),2))/6


%