% For in vivo spine data analyzed using MapManager V_20171028 
% Author: Richard Roth (rhroth@stanford.edu)
% Date: October 2019
% Based on MapManager-Matlab-1.0.1-alpha by Robert Cudmore, 20171016

%     Used for NearesNeighbor.m
%     Requires the following functions:
%     segmentCall_withDist.m
%     removeNanRows.m


function retPool = getmyMapListValuesCondDist_Neighbour(myMapList, mySegmentTable, Condition, SessionConditions, NeighborConditions, StatName, SecStatName, Channel, PersCond, MinSpines)

sessincl = ismember(SessionConditions,NeighborConditions);
nnball3D={};
nbsall3D={};
nnball=[];
numnbs = 8;

% only include segments that are not excluded
rows = mySegmentTable.Include == 1;
myGoodSegmentTable = mySegmentTable(rows,:);
clearvars rows

% only include segments that are part of condition to plot
rows=strcmp(myGoodSegmentTable.MouseCondition,Condition);
myCondSegmentTable = myGoodSegmentTable(rows,:);
clearvars SegmentNumber i rows values valuesgeo meanvalues stdev val valneurons valneuronsSegmentsmeans valmice valmiceSegmentsmeans valmiceNeuronsmeans valmiceNeuronsSegmentsmeans valneuronsSegmentsgeomeans
clearvars valmiceSegmentsgeomeans valmiceNeuronsSegmentsgeomeans valmiceNeuronsgeomeans se retPool NeuronNames NeuronName geomeanvalues currentvalues currentmousename count MapName MouseName
clearvars nbs nnb dist ind idx

for i = 1:size(myCondSegmentTable,1)
   
MapName = myCondSegmentTable.MapName{i};
MouseName = myCondSegmentTable.MouseName{i};
SegmentNumber = myCondSegmentTable.SegmentNumber(i);

[currentvalues,currentvaluesPDist] = segmentCall_withDist(myMapList, MapName, SegmentNumber, SessionConditions, StatName, SecStatName, Channel, PersCond, MinSpines);

% Invert pDist values if segmentorientation is inverted 
% Only applies if mySegmentTable has a column titled segmentoritentaion
% with -1 values for inverted orientations
if ismember('segmentoritentaion', myCondSegmentTable.Properties.VariableNames) == 1
SegmentOrientation = myCondSegmentTable.segmentoritentaion(i);
if strcmp(SecStatName, 'pDist') == 1
    if SegmentOrientation == -1
    currentvaluesPDist.values = -currentvaluesPDist.values;
    currentvaluesPDist.mean = -currentvaluesPDist.mean;
    currentvaluesPDist.geomean = -currentvaluesPDist.geomean;  
    end
end
end

[currentVals, currentValInds] = removeNanRows(currentvalues.values);

%pDist value of spine (average pDist of all time points)
nearestNeighbour = nanmean(currentvaluesPDist.values(currentValInds,:),2);
%Stat of interest value of spine (averaege of specified timepoints)
nearestNeighbour(:,2) = nanmean(currentVals(:,sessincl),2);

% sort by pDist
nearestNeighbourSorted=sortrows(nearestNeighbour,1);

%exclude segment if it only contains 3 or less spines
if size(nearestNeighbourSorted,1) <= 3
    nnb=NaN(length(nearestNeighbourSorted),3);
    nbs = NaN(size(nearestNeighbourSorted,1),numnbs);
else
    %%
    nnb=NaN(length(nearestNeighbourSorted),3);
row=2;
    nnb(row,1)=nearestNeighbourSorted(row,2);
    if abs(nearestNeighbourSorted(row,1)-nearestNeighbourSorted(row-1,1)) < abs(nearestNeighbourSorted(row+1,1)-nearestNeighbourSorted(row,1)) %Check which side is nearest neighbour
    nnb(row,2)=nearestNeighbourSorted(row-1,2); %Populate value for nearest neighbor 
    nnb(row,3)=abs(nearestNeighbourSorted(row,1)-nearestNeighbourSorted(row-1,1)); %Populate distance to nearest neighbor 
    else
    nnb(row,2)=nearestNeighbourSorted(row+1,2);
    nnb(row,3)=abs(nearestNeighbourSorted(row+1,1)-nearestNeighbourSorted(row,1));
    end
  
   
    %Remove spine pairs with distance over 3.5 um
    if nnb(row,3) >= 3.5
    nnb(row,1)=NaN;
    nnb(row,2)=NaN;
    nnb(row,3)=NaN;
    end
    if nnb(row,3) <= 1
    nnb(row,1)=NaN;
    nnb(row,2)=NaN;
    nnb(row,3)=NaN;
    end
for row=3:length(nearestNeighbourSorted)-1
    nnb(row,1)=nearestNeighbourSorted(row,2);
    if abs(nearestNeighbourSorted(row,1)-nearestNeighbourSorted(row-1,1)) < abs(nearestNeighbourSorted(row+1,1)-nearestNeighbourSorted(row,1))
    if nnb(row-1,2) == nearestNeighbourSorted(row,2)
    nnb(row,1)=NaN;
    nnb(row,2)=NaN;
    else        
    nnb(row,1)=nearestNeighbourSorted(row,2);    
    nnb(row,2)=nearestNeighbourSorted(row-1,2);
    nnb(row,3)=abs(nearestNeighbourSorted(row,1)-nearestNeighbourSorted(row-1,1));
    end
    
    else
    nnb(row,2)=nearestNeighbourSorted(row+1,2);
    nnb(row,3)=abs(nearestNeighbourSorted(row+1,1)-nearestNeighbourSorted(row,1));
    end
    
    %Remove spine pairs with distance over 3.5 um
    if nnb(row,3) >= 3.5
    nnb(row,1)=NaN;
    nnb(row,2)=NaN;
    nnb(row,3)=NaN;
    end
    if nnb(row,3) <= 1
    nnb(row,1)=NaN;
    nnb(row,2)=NaN;
    nnb(row,3)=NaN;
    end
end


%%

nbs = NaN(size(nearestNeighbourSorted,1),numnbs);
startrow = round(numnbs/2);
endrow = size(nearestNeighbourSorted,1)-round(numnbs/2);

for row = startrow:endrow
    [dist, ind] = sort(abs(nearestNeighbourSorted(:,1)-nearestNeighbourSorted(row,1)));
    
    
    for j = 1:numnbs
    nbs(row,j) = nearestNeighbourSorted(ind(j),2); 
    end
    if max(nbs(row,:)) == nbs(row,1) % only fill row if main spine was strongest in surrounding n
    else
    nbs(row,:) = NaN;
    end
    
end    

end

nnball=[nnball;nnb];

t=size(nnball3D,1);
nnball3D{t+1,1}=nnb;
nnball3D{t+1,2}=MapName;

nbsall3D{t+1,1}=nbs;
nbsall3D{t+1,2}=MapName;
nbsall3D{t+1,3}=nanmean(nearestNeighbourSorted(:,2));

retPool.nnball=nnball;
retPool.nnball3D=nnball3D;
retPool.nbsall3D=nbsall3D;

clearvars nbs nnb dist ind idx

end
retPool.nnball=nnball;
retPool.nnball3D=nnball3D;
retPool.nbsall3D=nbsall3D;
end