grnodes=xyz;
radiusmin=0.25;
radiusmax=1.5;        %initial value
radiusstep=0.5;
NodeinMesh=0;
distances=zeros(length(grnodes));
closedistances=zeros(length(grnodes));
conn=zeros(length(grnodes));
percentinmesh=0;

while NodeinMesh<0.95
tic        
for k=1:length(grnodes)
grnodescelli{k,1} = grnodes(k,:);
grnodescellj{1,k} = grnodes(k,:);
end
for i=1:length(grnodes)
    for j=1:length(grnodes)
distances(i,j)= pdist([grnodescelli{i,1};grnodescellj{1,j}]);
    end
end

close=find(distances>=radiusmin & distances<=radiusmax);
far=find(distances<radiusmin & distances>radiusmax);
closedistances(close(:,:))=distances(close(:,:));
closedistances(far(:,:))=0;
nearby=closedistances;

for i=1:length(grnodes)
    f1=find(nearby(i,:)~=0);
    if ~isempty(f1)
    g1=find(nearby(i,:)==min(nearby(i,f1)));  
    conn(i,g1)=1;
    nearby(g1,i)=0;
    nearby(i,g1)=0;
    else
    end
    f2=find(nearby(i,:)~=0);
    if ~isempty(f2)
    g2=find(nearby(i,:)==min(nearby(i,f2)));  
    conn(i,g2)=1;
    nearby(g2,i)=0;
    nearby(i,g2)=0;
    else
    end
end

[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', grnodes(:,1),'YData', grnodes(:,2),'ZData', grnodes(:,3))

NodeDeg=degree(G);
% figure
% histogram(NodeDeg);

NodeTotal=length(grnodes);
NodeConnDeg=NodeDeg(NodeDeg>1);
NodeinMesh=length(NodeConnDeg)/NodeTotal;
percentinmesh=[percentinmesh;NodeinMesh];
radiusmax=radiusmax+radiusstep;
toc
end
Rmax=radiusmax;

% NodeMeanConnDeg=mean(NodeConnDeg);
% NodeSDConnDeg=std(NodeConnDeg);
% NodeDist=nonzeros(closedistances);
% NodeMeanDist=mean(NodeDist);
% NodeSDDist=std(NodeDist);

% figure
% P=double(grnodes);
% plot3(P(:,1),P(:,2),P(:,3),'.');
% grid on
% hold on
%shp = alphaShape(P(:,1),P(:,2),P(:,3),0.5);
% plot(shp);
%Surf=surfaceArea(shp);

%SummaryTable=table(NodeTotal,NodeinMesh,NodeMeanDist,NodeSDDist,NodeMeanConnDeg,NodeSDConnDeg,Surf,'VariableNames',{'NodeTotal','NodeinMesh','NodeMeanDist','NodeSDDist','NodeMeanConnDeg','NodeSDConnDeg','Surf'});
 
   
% mvol=mean(V(find(e1==name(i))));
% sdvol=std(V(find(e1==name(i))));
% results{i,1}=name(i);
% results{i,2}=mvol;
% results{i,3}=sdvol;