
[ind,name]=grp2idx(OriginalImageName);

for n=1:length(name) %loop over single cells
tic
radiusmin=0.25;
radiusmax=1.5;        %initial value
radiusstep=0.1;
%preassign matrices for sinle cell
cellind=find(OriginalImageName==name(n));
nodexyz=xyz(cellind,:);
dist=zeros(length(cellind));
closedist=zeros(length(cellind));
conn=zeros(length(cellind));
%compute all paiwise distances between the mesh verteces
for k=1:length(cellind)
nodexyzi{k,1} = nodexyz(k,:);
nodexyzj{1,k} = nodexyz(k,:);
end
for i=1:length(cellind)
    for j=1:length(cellind)
dist(i,j)= pdist([nodexyzi{i,1};nodexyzj{1,j}]);
    end
end

NodeinMesh=0;          %need to pre-define for conditioning
while NodeinMesh<0.95  %run analysis until 95% of vertices are in the mesh

close=find(dist>=radiusmin & dist<=radiusmax);
far=find(dist<radiusmin | dist>radiusmax);
closedist(close(:,:))=dist(close(:,:));
closedist(far(:,:))=0;
nearby=closedist;

%create bond with 2 closest neibour nodes within the radius, if they exist
for i=1:length(cellind)
    neighbour1=find(nearby(i,:)~=0);
    if ~isempty(neighbour1)
    bond1=find(nearby(i,:)==min(nearby(i,neighbour1)));  
    conn(i,bond1)=1;
    nearby(bond1,i)=0;
    nearby(i,bond1)=0;
    else
    end
    neighbour2=find(nearby(i,:)~=0);
    if ~isempty(neighbour2)
    bond2=find(nearby(i,:)==min(nearby(i,neighbour2)));  
    conn(i,bond2)=1;
    nearby(bond2,i)=0;
    nearby(i,bond2)=0;
    else
    end
end
%create undirected graph
[row,col]=find(conn==1);
EdgeTable=table([row,col],'VariableNames',{'EndNodes'});    
G = graph(EdgeTable);
% figure
% plot(G,'EdgeColor','r','EdgeAlpha',0.3,'NodeColor','none','XData', nodexyz(:,1),'YData', nodexyz(:,2),'ZData', nodexyz(:,3))
NodeDeg=degree(G);
% figure
% histogram(NodeDeg);
NodeTotal=length(cellind);
NodeConnDeg=NodeDeg(NodeDeg>1);
NodeinMesh=length(NodeConnDeg)/NodeTotal;
% NodeMeanConnDeg=mean(NodeConnDeg);
% NodeSDConnDeg=std(NodeConnDeg);
% NodeDist=nonzeros(closedist);
% NodeMeanDist=mean(NodeDist);
% NodeSDDist=std(NodeDist);
radiusmax=radiusmax+radiusstep; %change max radius with a predefined step
end %end for the "while"

Rmax=radiusmax;
% figure
% plot3(nodexyz(:,1),nodexyz(:,2),nodexyz(:,3),'.');
% grid on
% hold on
% shp = alphaShape(nodexyz(:,1),nodexyz(:,2),nodexyz(:,3),radiusmax);
% plot(shp);
% Surf=surfaceArea(shp);

% mVol=mean(V(find(OriginalImageName==name(n))));
% sdVol=std(V(find(OriginalImageName==name(n))));
% 
results{n,1}=name(n);
results{n,2}=Rmax;
% results{n,2}=mVol;
% results{n,3}=sdVol;
% results{n,4}=NodeTotal;
% results{n,5}=NodeinMesh;
% results{n,6}=NodeMeanDist;
% results{n,7}=NodeSDDist;
% results{n,8}=NodeMeanConnDeg;
% results{n,9}=NodeSDConnDeg;
% results{n,10}=Surf;
clearvars -except results n OriginalImageName V xyz name
toc
end

%SummaryTable=array2table(results, 'VariableNames',{'name','mVol','sdVol','NodeTotal','NodeinMesh','NodeMeanDist','NodeSDDist','NodeMeanConnDeg','NodeSDConnDeg','Surf'});
SummaryTable=array2table(results, 'VariableNames',{'name','Rmax'});
