%This program calculates spike triggered landmark plot
%in v5 we take the median spike trig LM and divide by the median position
%triggered LM
clear all
close all
exp_date=input('please type in the experiment date : ','s');
load(strcat(exp_date,'_Stage1_output.mat'));
load(strcat(exp_date,'_Stage2_output.mat'));
ppc=7.5;%pixels per cm
for tn=2:length(hd_crct) %for all trials
    tn
    %% load the landmark locations data
    %% load the landmark locations data
    if filename{tn}(11)=='1' || filename{tn}(11)=='3' || filename{tn}(11)=='4' || filename{tn}(11)=='5'
        load(strcat('landmark_coordinates_',filename{tn}(1:8),'_t3.mat'));
        tankfilt=filter;
        clear filter
    else
        load(strcat('landmark_coordinates_',filename{tn}));
        tankfilt=filter;
        clear filter
    end
    y=0:10:1600;
    x=0:10:1200;
    resy=1600/10;
    resx=1200/10;
    LMoccup=zeros(resx,resy);
    filtgrid=polygrid(1200-tankfilt(:,2),tankfilt(:,1),0.01);
    hcgrid=polygrid(1200-hc(:,2),hc(:,1),0.01);
    stpgrid=polygrid(1200-stp(:,2),stp(:,1),0.01);
    ltpgrid=polygrid(1200-ltp(:,2),ltp(:,1),0.01);
    sqpgrid=polygrid(1200-sqp(:,2),sqp(:,1),0.01);
%%
    %for each spike find the center of the fish and set the grid around
    %that
    ndata=length(nuraw_count{tn});        
     for i=1:(resx)
        for j=1:(resy)
            %set the value of the grid points that are occupied by the tank
            %walls to 1
            temp1=((find((1200-tank_circle(2,:))>=x(i)&(1200-tank_circle(2,:))<x(i+1)&(tank_circle(1,:))>=y(j)&(tank_circle(1,:)<y(j+1)))));
            if ~isempty(temp1)
                LMoccup(i,j)=1;
            end
            %generate points that are inside the polygons belonging to
            %landmarks
            temp1=((find((filtgrid(:,1))>=x(i)&(filtgrid(:,1))<x(i+1)&(filtgrid(:,2))>=y(j)&(filtgrid(:,2)<y(j+1)))));
            if ~isempty(temp1)
                LMoccup(i,j)=1;
            end
            
            %check for other landmarks
            if filename{tn}(11)=='1'
                temp1=((find((hcgrid(:,1))>=x(i)&(hcgrid(:,1))<x(i+1)&(hcgrid(:,2))>=y(j)&(hcgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
            elseif filename{tn}(11)=='3'
                temp1=((find((stpgrid(:,1))>=x(i)&(stpgrid(:,1))<x(i+1)&(stpgrid(:,2))>=y(j)&(stpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((ltpgrid(:,1))>=x(i)&(ltpgrid(:,1))<x(i+1)&(ltpgrid(:,2))>=y(j)&(ltpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((sqpgrid(:,1))>=x(i)&(sqpgrid(:,1))<x(i+1)&(sqpgrid(:,2))>=y(j)&(sqpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((1200-sc(2,:))>=x(i)&(1200-sc(2,:))<x(i+1)&(sc(1,:))>=y(j)&(sc(1,:)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((1200-lc(2,:))>=x(i)&(1200-lc(2,:))<x(i+1)&(lc(1,:))>=y(j)&(lc(1,:)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
            elseif filename{tn}(11)=='4'
                temp1=((find((stpgrid(:,1))>=x(i)&(stpgrid(:,1))<x(i+1)&(stpgrid(:,2))>=y(j)&(stpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((sqpgrid(:,1))>=x(i)&(sqpgrid(:,1))<x(i+1)&(sqpgrid(:,2))>=y(j)&(sqpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((1200-sc(2,:))>=x(i)&(1200-sc(2,:))<x(i+1)&(sc(1,:))>=y(j)&(sc(1,:)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((1200-lc(2,:))>=x(i)&(1200-lc(2,:))<x(i+1)&(lc(1,:))>=y(j)&(lc(1,:)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
            elseif filename{tn}(11)== '5'
                temp1=((find((stpgrid(:,1))>=x(i)&(stpgrid(:,1))<x(i+1)&(stpgrid(:,2))>=y(j)&(stpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((ltpgrid(:,1))>=x(i)&(ltpgrid(:,1))<x(i+1)&(ltpgrid(:,2))>=y(j)&(ltpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((sqpgrid(:,1))>=x(i)&(sqpgrid(:,1))<x(i+1)&(sqpgrid(:,2))>=y(j)&(sqpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((1200-lc(2,:))>=x(i)&(1200-lc(2,:))<x(i+1)&(lc(1,:))>=y(j)&(lc(1,:)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
            elseif filename{tn}(11)== '6' || filename{tn}(11)== '7'
                temp1=((find((stpgrid(:,1))>=x(i)&(stpgrid(:,1))<x(i+1)&(stpgrid(:,2))>=y(j)&(stpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((ltpgrid(:,1))>=x(i)&(ltpgrid(:,1))<x(i+1)&(ltpgrid(:,2))>=y(j)&(ltpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((sqpgrid(:,1))>=x(i)&(sqpgrid(:,1))<x(i+1)&(sqpgrid(:,2))>=y(j)&(sqpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((1200-lc(2,:))>=x(i)&(1200-lc(2,:))<x(i+1)&(lc(1,:))>=y(j)&(lc(1,:)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
            else
                temp1=((find((stpgrid(:,1))>=x(i)&(stpgrid(:,1))<x(i+1)&(stpgrid(:,2))>=y(j)&(stpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
                temp1=((find((ltpgrid(:,1))>=x(i)&(ltpgrid(:,1))<x(i+1)&(ltpgrid(:,2))>=y(j)&(ltpgrid(:,2)<y(j+1)))));
                if ~isempty(temp1)
                    LMoccup(i,j)=1;
                end
            end
        end
    end
    %%calculate the position triggered landmark average. This will be
    %%calculated for each trial and used to normalize spike triggered LM to
    %%account for time spent at a given location'
    head=downsample(hd_crct{tn},20);
    tail=downsample(COM_crct{tn},20);%the mid point between head and tail is set at the center.
    %but since the tail is not really tail,
    clear dist
    clear index
    for i=1:length(head)
            dist(i)=pdist([head(i,1) head(i,2);tail(i,1) tail(i,2)],'euclidean');
            mid(i,:)=[(head(i,1)+tail(i,1))/2,(head(i,2)+tail(i,2))/2];
            h2tvec=[(head(i,1)-tail(i,1)),(head(i,2)-tail(i,2)),0];
            %calculate the translation vector that brings the midpoint to
            %construct a grid around the center point
            %calculate the heading angle of the fish relative to the north
            %pointing arrow
            n0=[800,400];
            n1=[800,600];
            N=[n1(1)-n0(1),n1(2)-n0(2),0];
            crossvec=cross(h2tvec,N);
            ThetaInDegreesL(i) =atan2d(norm(crossvec),dot(N,h2tvec))*sign(crossvec(3));
            x_center=800;
            y_center=600;
            center = repmat([x_center; y_center], 1, 1);
            R = [cosd(ThetaInDegreesL(i)) -sind(ThetaInDegreesL(i)); sind(ThetaInDegreesL(i)) cosd(ThetaInDegreesL(i))];
            midshift=mid(i,:)'-center;
            s0=R*midshift;
            midshift=s0+center;
            transvect=[(800-midshift(1))/10,(600-midshift(2))/10,];
            %construct the landmark configuration around the fish
            if ~isnan(ThetaInDegreesL(i))
                RotateImage=imrotate(LMoccup, -1*ThetaInDegreesL(i),'bilinear','crop');
                TranslateImageL{i}=imtranslate(RotateImage, transvect);
            else
                TranslateImageL{i}=NaN(size(LMoccup));
            end
    end
    %remove points where the distance between head and tail is too much,
    %this is taking out the points where tracking failed.
    range1=nanmean(dist)+3*nanstd(dist);
    range2=nanmean(dist)-3*nanstd(dist);
    index=find(dist(dist<range1 & dist>range2));
    M=zeros(length(index),resx,resy);
    for i=1:length(index)
            M(i,:,:)=cell2mat(TranslateImageL(index(i)));
    end
    
    PosTrigLM{tn}=squeeze(nansum(M));
    %only take pixels where the fish spent more than 1/10th of the tital
    %frames
    

    
    figure
    imagesc(PosTrigLM{tn})
    %colorbar
    set(gca,'YDir','normal')
    hold on
    %saveas(gcf,strcat(filename{tn},'_Postriglm_median.pdf'))
 
%    for i=1:120
%         for j=1:160
%             if PosTrigLM{tn}(i,j)<(length(index)/10)%~6 min or 1/10th of the trial length
%                 PosTrigLM{tn}(i,j)=NaN;
%             end
%         end
%     end   
        max_val=max(reshape(PosTrigLM{tn},[1,120*160]));
    for i=1:120
        for j=1:160
            %if PosTrigLM{tn}(i,j)<(length(index)*.05) || i<40 || i>70 || j<70 || j>90 %has to have spent more than 5% of the trial length
             if  i<40 || i>70 || j<70 || j>90
                PosTrigLM{tn}(i,j)=NaN;
            end
        end
    end   
        
    
        figure
    imagesc(PosTrigLM{tn})
    colorbar
    set(gca,'YDir','normal')
    hold on
    for w=1:ndata
        clear dist
        for k=1:nuraw_count{tn}(w)
            %figure
            %for each spike find the mid point between head and tail
            dist(k)=pdist([Place_int_unit{tn,w,1}(k) Place_int_unit{tn,w,2}(k);Place_int_unit_tail{tn,w,1}(k) Place_int_unit_tail{tn,w,2}(k)]);
            mid=[(Place_int_unit{tn,w,1}(k)+Place_int_unit_tail{tn,w,1}(k))/2,(Place_int_unit{tn,w,2}(k)+(Place_int_unit_tail{tn,w,2}(k)))/2];
            h2tvec=[Place_int_unit{tn,w,1}(k)-Place_int_unit_tail{tn,w,1}(k),Place_int_unit{tn,w,2}(k)-(Place_int_unit_tail{tn,w,2}(k)),0];
            %calculate the translation vector that brings the midpoint to
            %construct a grid around the center point
            %calculate the heading angle of the fish relative to the north
            %pointing arrow
            n0=[800,400];
            n1=[800,600];
            N=[n1(1)-n0(1),n1(2)-n0(2),0];
            crossvec=cross(h2tvec,N);
            ThetaInDegrees(k) =atan2d(norm(crossvec),dot(N,h2tvec))*sign(crossvec(3));
            x_center=800;
            y_center=600;
            center = repmat([x_center; y_center], 1, 1);
            R = [cosd(ThetaInDegrees(k)) -sind(ThetaInDegrees(k)); sind(ThetaInDegrees(k)) cosd(ThetaInDegrees(k))];
            midshift=mid'-center;
            s0=R*midshift;
            midshift=s0+center;
            transvect=[(800-midshift(1))/10,(600-midshift(2))/10,];
            %construct the landmark configuration around the fish
            RotateImage=imrotate(LMoccup, -1*ThetaInDegrees(k),'bilinear','crop');
            TranslateImage{k,w}=imtranslate(RotateImage, transvect);
        end
        range1=nanmean(dist)+3*nanstd(dist);
        range2=nanmean(dist)-3*nanstd(dist);
        index=find(dist(dist<range1 & dist>range2));

         M=zeros(length(index),resx,resy);
        for i=1:length(index)
            M(i,:,:)=cell2mat(TranslateImage(index(i),w));
        end
        SpkTrigLM{tn,w}=squeeze(nansum(M));
        figure
        imagesc(SpkTrigLM{tn,w})
        colorbar
        set(gca,'YDir','normal')
        hold on
        line([80 80],[47.5 62.5])
%         axis equal
        title(strcat('trial number',num2str(tn),'cell',uname{w,tn}))
        %saveas(gcf,strcat(filename{tn},'_',uname{w,tn},'_spktriglm_median.pdf'))
        %plot normalized spike triggered LM
        
        figure
        NormSpkTrigLM{tn,w}=(SpkTrigLM{tn,w}./(PosTrigLM{tn}));
        %plot a up to 5 cm distance fom the fish
        rad=3*ppc;
        imagesc(NormSpkTrigLM{tn,w})
        colorbar
        set(gca,'YDir','normal')
        hold on
        title(strcat('trial number',num2str(tn),'cell',uname{w,tn},'Normalized'))
        %saveas(gcf,strcat(filename{tn},'_',uname{w,tn},'_Normspktriglm_median.pdf'))

        
        %calculate the location of the peak in the lmtriggered average
        [maxval,idx]=(max(SpkTrigLM{tn,w}(:)));
        [row,col]=ind2sub(size(SpkTrigLM{tn,w}), idx);
        %calculate distance to the head of the fish in pixels
        %a point between fish's head and its center of mass is rowf=80,colf=60

        dist_cm{tn,w}=sqrt((row-60)^2+(col-80)^2)/ppc
        
        %calculate the location of the peak in the normalized lmtriggered average
        [maxval,idx]=(max(NormSpkTrigLM{tn,w}(:)));
        [row,col]=ind2sub(size(NormSpkTrigLM{tn,w}), idx);
        %calculate distance to the head of the fish in pixels
        %a point between fish's head and its center of mass is rowf=80,colf=60

        dist_cm_norm{tn,w}=sqrt((row-60)^2+(col-80)^2)/ppc

    

close all

    end

end
close all
%unit 1
for i=1:7
    if i~=1
        SpkTrigAllU1(i-1,:,:)=cell2mat(SpkTrigLM(i,1));
        PosTrigAllU1(i-1,:,:)=cell2mat(PosTrigLM(i));

        
    end
    
end
sumSpkTrigLMU1=(squeeze(nansum(SpkTrigAllU1)));
sumPosTrigLMU1=(squeeze(nansum(PosTrigAllU1)));

PosMean=prctile(reshape(sumPosTrigLMU1(40:70,70:90),1,31*21),10);
for i=1:120
    for j=1:160
        %if sumPosTrigLMU1(i,j)<(length(LtotalU1)*0.01) || i<40 || i>70 || j<70 || j>90
        if sumPosTrigLMU1(i,j)<PosMean || i<40 || i>70 || j<70 || j>90            
            sumPosTrigLMU1(i,j)=NaN;
        end
    end
end


NormSpkTrigLMU1=(sumSpkTrigLMU1./sumPosTrigLMU1);


temp=imgaussfilt(NormSpkTrigLMU1);
[maxval,idx]=max(temp(:))
[row,col]=ind2sub(size(temp), idx);
MaxLoc(1,:)=[row-57,col-80];%subtract the fish com
stlmUname{1}='U1';

reg_max(1,1)=max(reshape(temp(58:70,81:90),[1,numel(temp(58:70,81:90))]))
reg_max(1,2)=max(reshape(temp(58:70,70:80),[1,numel(temp(58:70,70:80))]));
reg_max(1,3)=max(reshape(temp(40:57,70:80),[1,numel(temp(40:57,70:80))]));
reg_max(1,4)=max(reshape(temp(40:57,81:90),[1,numel(temp(40:57,81:90))]));
%calculate an index head vs tail bias
index_side_pref(1,1)=(max(reg_max(1,1:2))-max(reg_max(1,3:4)))/(max(reg_max(1,1:2))+max(reg_max(1,3:4)))
index_side_pref(1,2)=(max(reg_max(1,[1,4]))-max(reg_max(1,[2,3])))/nansum([max(reg_max(1,[1,4])),max(reg_max(1,[2,3]))]);

temp=temp/maxval;
figure
imagesc(temp)
hold on
title('U1')
colorbar
set(gca,'YDir','normal')
hold on
line([80 80],[47.5 62.5])
xlim([71,89])
ylim([41,69])
axis equal

%saveas(gcf,'NormU1spktriglm.pdf')



%unit 4
for i=3:7
   SpkTrigAllU4(i-2,:,:)=cell2mat(SpkTrigLM(i,2));
   PosTrigAllU4(i-2,:,:)=cell2mat(PosTrigLM(i));

end
sumSpkTrigLMU4=(squeeze(nansum(SpkTrigAllU4)));
sumPosTrigLMU4=(squeeze(nansum(PosTrigAllU4)));


PosMean=prctile(reshape(sumPosTrigLMU4(40:70,70:90),1,31*21),10);
for i=1:120
    for j=1:160
        %if sumPosTrigLMU1(i,j)<(length(LtotalU1)*0.01) || i<40 || i>70 || j<70 || j>90
        if sumPosTrigLMU4(i,j)<PosMean || i<40 || i>70 || j<70 || j>90            
            sumPosTrigLMU4(i,j)=NaN;
        end
    end
end


NormSpkTrigLMU4=(sumSpkTrigLMU4./sumPosTrigLMU4);


temp=imgaussfilt(NormSpkTrigLMU4);

[maxval,idx]=max(temp(:))
[row,col]=ind2sub(size(temp), idx);
MaxLoc(2,:)=[row-57,col-80];%subtract the fish com
stlmUname{2}='U4';
temp=temp/maxval;
reg_max(2,1)=max(reshape(temp(58:70,81:90),[1,numel(temp(58:70,81:90))]))
reg_max(2,2)=max(reshape(temp(58:70,70:80),[1,numel(temp(58:70,70:80))]));
reg_max(2,3)=max(reshape(temp(40:57,70:80),[1,numel(temp(40:57,70:80))]));
reg_max(2,4)=max(reshape(temp(40:57,81:90),[1,numel(temp(40:57,81:90))]));
%calculate an index head vs tail bias
index_side_pref(2,1)=(max(reg_max(2,1:2))-max(reg_max(2,3:4)))/(max(reg_max(2,1:2))+max(reg_max(2,3:4)))
index_side_pref(2,2)=(max(reg_max(2,[1,4]))-max(reg_max(2,[2,3])))/nansum([max(reg_max(2,[1,4])),max(reg_max(2,[2,3]))]);

figure
imagesc(temp)
hold on
title('U4')
colorbar
set(gca,'YDir','normal')
hold on
line([80 80],[47.5 62.5])
xlim([71,89])
ylim([41,69])
axis equal

%saveas(gcf,'NormU4spktriglm.pdf')
%unit 5


for i=3:7
   SpkTrigAllU5(i-2,:,:)=cell2mat(SpkTrigLM(i,3));
   PosTrigAllU5(i-2,:,:)=cell2mat(PosTrigLM(i));

end
sumSpkTrigLMU5=(squeeze(nansum(SpkTrigAllU5)));
sumPosTrigLMU5=(squeeze(nansum(PosTrigAllU5)));


PosMean=prctile(reshape(sumPosTrigLMU5(40:70,70:90),1,31*21),10);
for i=1:120
    for j=1:160
        %if sumPosTrigLMU1(i,j)<(length(LtotalU1)*0.01) || i<40 || i>70 || j<70 || j>90
        if sumPosTrigLMU5(i,j)<PosMean || i<40 || i>70 || j<70 || j>90            
            sumPosTrigLMU5(i,j)=NaN;
        end
    end
end


NormSpkTrigLMU5=(sumSpkTrigLMU5./sumPosTrigLMU5);


temp=imgaussfilt(NormSpkTrigLMU5);
[maxval,idx]=max(temp(:))
[row,col]=ind2sub(size(temp), idx);
MaxLoc(3,:)=[row-57,col-80];%subtract the fish com
stlmUname{3}='U5';
temp=temp/maxval;
reg_max(3,1)=max(reshape(temp(58:70,81:90),[1,numel(temp(58:70,81:90))]))
reg_max(3,2)=max(reshape(temp(58:70,70:80),[1,numel(temp(58:70,70:80))]));
reg_max(3,3)=max(reshape(temp(40:57,70:80),[1,numel(temp(40:57,70:80))]));
reg_max(3,4)=max(reshape(temp(40:57,81:90),[1,numel(temp(40:57,81:90))]));
%calculate an index head vs tail bias
index_side_pref(3,1)=(max(reg_max(3,1:2))-max(reg_max(3,3:4)))/(max(reg_max(3,1:2))+max(reg_max(3,3:4)))
index_side_pref(3,2)=(max(reg_max(3,[1,4]))-max(reg_max(3,[2,3])))/nansum([max(reg_max(3,[1,4])),max(reg_max(3,[2,3]))]);

figure
imagesc(temp)
hold on
title('U5')
colorbar
set(gca,'YDir','normal')
hold on
line([80 80],[47.5 62.5])
xlim([71,89])
ylim([41,69])
axis equal

%saveas(gcf,'NormU5spktriglm.pdf')

%unit A
SpkTrigAllUA(1,:,:)=cell2mat(SpkTrigLM(7,4))
sumSpkTrigLMUA=(squeeze(SpkTrigAllUA(1,:,:)));
sumPosTrigLMUA=cell2mat(PosTrigLM(7));
PosMean=prctile(reshape(sumPosTrigLMUA(40:70,70:90),1,31*21),10);
for i=1:120
    for j=1:160
        %if sumPosTrigLMU1(i,j)<(length(LtotalU1)*0.01) || i<40 || i>70 || j<70 || j>90
        if sumPosTrigLMUA(i,j)<PosMean || i<40 || i>70 || j<70 || j>90            
            sumPosTrigLMUA(i,j)=NaN;
        end
    end
end


NormSpkTrigLMUA=(sumSpkTrigLMUA./sumPosTrigLMUA);


temp=imgaussfilt(NormSpkTrigLMUA);


[maxval,idx]=max(temp(:))
[row,col]=ind2sub(size(temp), idx);
MaxLoc(4,:)=[row-57,col-80];%subtract the fish com
stlmUname{4}='UA';
temp=temp/maxval
reg_max(4,1)=max(reshape(temp(58:70,81:90),[1,numel(temp(58:70,81:90))]))
reg_max(4,2)=max(reshape(temp(58:70,70:80),[1,numel(temp(58:70,70:80))]));
reg_max(4,3)=max(reshape(temp(40:57,70:80),[1,numel(temp(40:57,70:80))]));
reg_max(4,4)=max(reshape(temp(40:57,81:90),[1,numel(temp(40:57,81:90))]));
%calculate an index head vs tail bias
index_side_pref(4,1)=(max(reg_max(4,1:2))-max(reg_max(4,3:4)))/(max(reg_max(4,1:2))+max(reg_max(4,3:4)))
index_side_pref(4,2)=(max(reg_max(4,[1,4]))-max(reg_max(4,[2,3])))/nansum([max(reg_max(4,[1,4])),max(reg_max(4,[2,3]))]);

figure
imagesc(temp)
hold on
title('UA')
colorbar
set(gca,'YDir','normal')
hold on
line([80 80],[47.5 62.5])
xlim([71,89])
ylim([41,69])
axis equal

%saveas(gcf,'NormUAspktriglm.pdf')

%unit B
SpkTrigAllUB(1,:,:)=cell2mat(SpkTrigLM(6,4))
sumSpkTrigLMUB=(squeeze(SpkTrigAllUB(1,:,:)));
sumPosTrigLMUB=(cell2mat(PosTrigLM(6)));


PosMean=prctile(reshape(sumPosTrigLMUB(40:70,70:90),1,31*21),5);
for i=1:120
    for j=1:160
        %if sumPosTrigLMU1(i,j)<(length(LtotalU1)*0.01) || i<40 || i>70 || j<70 || j>90
        if sumPosTrigLMUB(i,j)<PosMean || i<40 || i>70 || j<70 || j>90            
            sumPosTrigLMUB(i,j)=NaN;
        end
    end
end


NormSpkTrigLMUB=(sumSpkTrigLMUB./sumPosTrigLMUB);


temp=imgaussfilt(NormSpkTrigLMUB);
[maxval,idx]=max(temp(:))
[row,col]=ind2sub(size(temp), idx);
MaxLoc(5,:)=[row-57,col-80];%subtract the fish com
stlmUname{5}='UB';
temp=temp/maxval;
reg_max(5,1)=max(reshape(temp(58:70,81:90),[1,numel(temp(58:70,81:90))]))
reg_max(5,2)=max(reshape(temp(58:70,70:80),[1,numel(temp(58:70,70:80))]));
reg_max(5,3)=max(reshape(temp(40:57,70:80),[1,numel(temp(40:57,70:80))]));
reg_max(5,4)=max(reshape(temp(40:57,81:90),[1,numel(temp(40:57,81:90))]));
%calculate an index head vs tail bias
index_side_pref(5,1)=(max(reg_max(5,1:2))-max(reg_max(5,3:4)))/(max(reg_max(5,1:2))+max(reg_max(5,3:4)))
index_side_pref(5,2)=(max(reg_max(5,[1,4]))-max(reg_max(5,[2,3])))/nansum([max(reg_max(5,[1,4])),max(reg_max(5,[2,3]))]);

    
figure
imagesc(temp)
hold on
title('UB')
colorbar
set(gca,'YDir','normal')
hold on
line([80 80],[47.5 62.5])
xlim([71,89])
ylim([41,69])
axis equal

%saveas(gcf,'NormUBpktriglm.pdf')


figure
histogram(cell2mat(dist_cm(:)),10, 'Normalization','Probability')
xlabel('Distance to LM peak (cm)')
ylabel('Probability')
figure
histogram(cell2mat(dist_cm_norm(:)),10, 'Normalization','Probability')
xlabel('Distance to LM peak of Normalized plot(cm)')
ylabel('Probability')
[t,r]=cart2pol(MaxLoc(:,2),MaxLoc(:,1));
ThetaLM=t*180/pi
RadLM=r;

%save(strcat(exp_date,'_Stage9_output_v6.mat'),'index_side_pref','reg_max','SpkTrigLM','dist_cm','dist_cm_norm','MaxLoc','ThetaLM','stlmUname','RadLM');