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

Hi,

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);
    Axes(:,1)=aaa+(bbb-aaa)*rand(N,1);
    for k=2:M
    aaa=randperm(N);
    Axes(:,k)=Axes(aaa,1);
    end
    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 );
    else
    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 );
    end
    end
    end
    
    end
    
    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');
    end
    
    %% 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
    tic
    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);
    
    
    end
    
    plot3(Xgrid(active),Ygrid(active),Zgrid(active),'o','MarkerFaceColor','red');
    axis equal 
    grid on
    toc

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]

https://en.m.wikipedia.org/wiki/Equivalence_class build Equivalence_class relationships between the spheres according to whether they intersect. Afterwards, look for the largest equivalence class.

[Chris Dan]

Hii Walter Roberson

Do you know any sample code which i can look at ?

[Walter Roberson]

Untested.

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)) = {};
end

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 !)

[Chris Dan]

Thank you so much :-)

Answer 1

[Jon]

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, 
plot(G)

% 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]

Nice approach!

[Chris Dan]

Thank you so much Jon

Youd answer is really nice and solves the problem :-)

[Jon]

Glad to hear this worked. It is an interesting problem

[Chris Dan]

Hi Jon,

Do you know if I have to find out the cente of the agglomerate, how it can be done using graph theory?

[Jon]

How are your defining "center of the agglomerate"? Do you mean the center of mass (volume) of the collection?

[Jon]

For the center of mass, noting that the 4/3*pi*rho, where rho is the density, terms drop out, you could use:

rCubed = r(idxKeep).^3;
rCubedTotal = sum(rCubed);
xc = sum(rCubed .* x(idxKeep))/rCubedTotal;
yc = sum(rCubed .* y(idxKeep))/rCubedTotal;
zc = sum(rCubed .* z(idxKeep))/rCubedTotal;

[Chris Dan]

Hi Jon

Iam doing like this

Nx =size(active,1);
Ny =size(active,2);
Nz = size(active,3);
[Xi,Yi,Zi] = ndgrid(0:Nx-1,0:Ny-1,0:Nz-1);
ind = find(active(:));
Xfile = Xi(ind);
Yfile = Yi(ind);
Zfile = Zi(ind);
center_agglo_1 = [mean(Xfile) mean(Yfile) mean(Zfile)]; % using dipole positions of the whole object

[Chris Dan]

Hi Jon

is this correct approach?

[Jon]

I don't know enough about your application to know what the appropriate definition of "center of agglomerate" should be.

Also, are you looking for the "Center" of all of the spheres, or just the ones in the largest clump.

Assuming your index Xfile, are the indices of the spheres you want included in the calculation , then your calculation above, would give the number mean center, rather than my calculation which gives the mass mean center. Another words, your calculation weights the mean only according to how many spheres are in various postion without regard for their size. My calculation weights the center according to the mass of each sphere.

So for example suppose you had only two spheres, one at x = 0 the other at x = 1, and the one at x = 1 was much larger than the one at x = 0. Your calculation would give the center of the collection as x = 1/2, whereas my approach (center of mass) would put it to a value close to 1.

[Walter Roberson]

It is a continuation of https://www.mathworks.com/matlabcentral/answers/1733060-find-the-center-of-the-cluster-made-from-random-spheres-having-different-diameters 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]

@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 https://www.mathworks.com/matlabcentral/answers/1733060-find-the-center-of-the-cluster-made-from-random-spheres-having-different-diameters I would suggest not continuing here and just moving to your new question