Find closest sphere in a cluster with random spheres having different diameters-2

3 views (last 30 days)
I posted this question before and i got good help in how to find closest sphere in a cluster with random spheres having different diameters. I am making clusters of spheres using random number generator. Some of the spheres are connected and some of them are not connected. I want to exclude the not connected ones, so I only get the clusters/connected spheres .I have random spheres having random diameters between 0.1 and 0.5.
here is my code till now:
%% Data for Agglomerates
N = 15; % number of spheres
a = 0.1 ; % lowest diameter of sphere
b = 0.5 ; % highest diameter of sphere
Diam = a + (b-a).*rand(N,1);
aaa= 0; % minimum center x and y coordinate limit for spheres
bbb= 1.5 ; % maximum center x and y coordinat limit for sphere
M=3 ;
Axes= zeros(N,M);
for k=2:M
Axes_Label ={'X axis','Y axis','Z axis'};
Data_agglo = [Diam Axes];
Data_Label ={'Diameter','X axis','Y axis','Z axis'};
R = Diam ./2;
%% Algorithm to find connected spheres and eliminate outliers
distances = zeros(size(Data_agglo,1));
active_spheres = zeros(size(Data_agglo,1));
delta = 1e-4;
for p =1:1:size(Data_agglo,1)
for q = 1:1:size(Data_agglo,1)
if size(Data_agglo,1)==1
distances(p,q) = sqrt((Data_agglo(q,2)-Data_agglo(p,2)).^2 + (Data_agglo(q,3)-Data_agglo(p,3)).^2 + (Data_agglo(q,4)-Data_agglo(p,4)).^2);
active_spheres(p,q) = (distances(p,q) <= (Data_agglo(p,1)/2)+(Data_agglo(q,1)/2) - delta );
if p ~= q
distances(p,q) = sqrt((Data_agglo(q,2)-Data_agglo(p,2)).^2 + (Data_agglo(q,3)-Data_agglo(p,3)).^2 + (Data_agglo(q,4)-Data_agglo(p,4)).^2);
active_spheres(p,q) = (distances(p,q) <= (Data_agglo(p,1)/2)+(Data_agglo(q,1)/2) - delta );
rem_particle = sum(active_spheres);
rem_particle = find(rem_particle == 0);
Data_agglo2 = Data_agglo;
Data_agglo(rem_particle,:) = [];
R(rem_particle) = [];
if isempty(R)
error('No spheres are connected!\n');
%% generate mesh
dipS = 0.01;
xmin = min(Data_agglo(:,2)-R);
xmax = max(Data_agglo(:,2)+R);
ymin = min(Data_agglo(:,3)-R);
ymax = max(Data_agglo(:,3)+R);
zmin = min(Data_agglo(:,4)-R);
zmax = max(Data_agglo(:,4)+R);
[Xgrid,Ygrid,Zgrid]= ndgrid((xmin:dipS:xmax)-(xmin+xmax)/2,(ymin:dipS:ymax)-(ymin+ymax)/2,(zmin:dipS:zmax)-(zmin+zmax)/2);
Data_agglo(:,2:4) = Data_agglo(:,2:4) - [(xmin+xmax)/2,(ymin+ymax)/2,(zmin+zmax)/2];
% plot3(Xgrid(:),Ygrid(:), Zgrid(:),'.','MarkerFaceColor','blue');
% hold all
%% get active dipoles
active = false(size(Xgrid));
for i =1:1:size(Data_agglo,1)
active = active | (((Xgrid - Data_agglo(i,2)).^2 + (Ygrid - Data_agglo(i,3)).^2 + (Zgrid - Data_agglo(i,4)).^2) <= R(i).^2);
axis equal
grid on
The problem I am getting is :
I eliminate the spheres, but I am left with two sets of connected spheres. I want to choose one of them and delete the other one. there is no pirority which one to use.
I am attaching the figure file as well. Other problem I get is that when i give the number of spheres =1, my program doesnt work. Also if it is possible to find the largest set of connected spheres as shown in the picture and keep them?
Does anyone knows how it can be done?

Sign in to comment.

Accepted Answer

Jon on 4 May 2022
Edited: Jon on 4 May 2022
Here is another approach, in which we consider the spheres to be nodes in an undirected graph with edges connecting the nodes whose corresponding spheres touch. We can then use some graph theory algorithms to find the connected components, which I think are just what you are looking for. Here is some example code, untested except to see that it runs (I didn't follow your notation, or scaling of the locations etc, just wanted to get the main idea across. I think you can adapt what I did here to the specifics of your problem)
% parameters
N = 15; % number of spheres
% box dimensions that centers will fall in
xMax = 1
yMax = 2
zMax = 1
% minimum and maximum diameters
a = 0.1;
b = 0.5;
% Generate random center points
x = rand(N,1)*xMax;
y = rand(N,1)*yMax;
z = rand(N,1)*zMax;
% Generate random diameters,
d = a + (b-a)*rand(N,1);
% Find the minimum separation between sphere centers at which they will
% touch
r = d/2; % radius
Stouch = r + r'; % use scalar expansion between column and row to get every pair
% Find the distance between every pair of points, making use of scalar
% expansion of difference between a column and row to get every pair
S = sqrt((x-x').^2 + (y-y').^2 + (z -z').^2)
% Generate logical matrix whose elements are set to 1 (true) if the ith and
% jth sphere touch
Touch = S <= Stouch;
% Consider as a graph problem with a node for each sphere and edges
% connecting the vertices if the corresponding spheres touch. So our Touch
% matrix is exactly the adjacency matrix of this graph
% Obviously each sphere touches itself
% so don't include self loops (edges that connect the node to itself)
G = graph(Touch,'omitselfloops');
% Plot the graph just to visualize it,
% Find the connected components, there is a path between any two nodes in a
% connected component, this means the connected components are exactly the
% clusters of connected spheres you are looking for
components = conncomp(G); % components(i) gives the cluster number the ith sphere belongs to
% Find the connected component with the most spheres in it using groupcounts
% gc is the number of spheres in each group
[gcnt,grp] = groupcounts(components(:)); % use (:) to make sure it is a column
[~,idx] = max(gcnt); % index of group with the most spheres
grpMax = grp(idx); % group number of the group with the most spheres
% Get a indices of the spheres that are in the largest cluster that we wish
% to keep
idxKeep = find(components == grpMax);
>> idxKeep
idxKeep =
2 3 5 7 11 12
Jon on 7 Jun 2022
@Walter Roberson, thanks for catching that. It has been long enough since looking at this, that I forgot that the spheres could overlap. Also I was unaware that a new question had been opened on this topic.
@hamzah khan looks like you have a lot of suggestions to work with already in I would suggest not continuing here and just moving to your new question

Sign in to comment.

More Answers (0)


Find more on Statistics and Machine Learning Toolbox in Help Center and File Exchange




Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!