Clear Filters
Clear Filters

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

1 view (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?
Walter Roberson
Walter Roberson on 4 May 2022
equivalence_class_numbers = (1:number_of_spheres).';
equivalence_class_members = num2cell(equilvance_class_numbers);
for idx = 2 : number_of_spheres
sp = Index_of_conflicting_spheres(idx, sphere_coordinates, sphere_radii);
if isempty(sp); continue; end
sp = [sp; idx];
new_class = min(sp);
new_members = unique(cell2mat(equivalence_class_members(sp)));
equivalance_class_members(new_class) = new_members;
equivalance_class_numbers(sp) = new_class;
equivalance_class_members(setdiff(sp, new_class)) = {};
At the end of this, equivalance_class_numbers will be a vector of length number_of_spheres that gives the cluster number each sphere belongs to. The numbers will likely have gaps -- for example if sphere #3 is connected to sphere #1 then the class number (1) and at (3) would both be 1, and 3 would not appear in any class number. For some purposes you would want to renumber without gaps (findgroups() would help for that), but for the purpose of your particular need, you do not care about the cluster number in absolute terms, you just want to know what is joined to what.
At the end of the code, equivalance_class_members is a cell, and for any cluster number that does not exist (like the 3 in the above example, because sphere 3 is part of cluster 1), the entry in equivalance_class_members will be empty. For each cluster number that is populated, equivalance_class_members will be a list of all the sphere numbers that belong to the same cluster.
If you see a equivalance_class_members entry that is a scalar (but not empty), that means that the sphere is isolated, and in such a case equivalance_class_numbers(K) == K .
You can cellfun(@length, equivalance_class_members) and use max() to find the largest cluster (at least in terms of numbers of members -- not in terms of volume !)

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
Walter Roberson
Walter Roberson on 7 Jun 2022
It is a continuation of in which it turns out that the spheres are overlapping but they want the centroid. If not for the overlap one could weight the center of the sphere by its volume, but with the overlap it gets nasty to solve analytically.
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 2-D and 3-D Plots 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!