Find the center of the cluster made from random spheres having different diameters

Hi,

I am making clusters using spheres having different diameters. The sphere's diameters and their center points are generated using random number generator. The diameters of the spheres lie between 0.1 and 0.5. I have to find the center point of the cluster which I generate. One idea I have is to take the average of all the spheres center points and it will be the center of the agglomerate but I am not sur if it is correct.

My code till now:

clear all; close all;

%% Data for Agglomerates
N = 50; % number of spheres
a = 0.1 ; % lowest diameter of sphere
b = 0.5 ; % highest diameter of sphere
Diam = a + (b-a).*rand(N,1);

% box dimensions that centers will fall in
xMax = 1;
yMax = 1;
zMax = 1;

% Generate random center points
x = rand(N,1)*xMax;
y = rand(N,1)*yMax;
z = rand(N,1)*zMax;

Axes = [x y z];

Data_agglo = [Diam Axes];
Data_Label ={'Diameter','X axis','Y axis','Z axis'};

R = Diam ./2;


%% Algorithm to find connected spheres and eliminate outliers-2
Data_agglo2 = Data_agglo;

% Find the minimum separation between sphere centers at which they will
% touch

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;
G = graph(Touch,'omitselfloops');
plot(G);

components = conncomp(G);
[gcnt,grp] = groupcounts(components(:));
[~,idx] = max(gcnt); % index of group with the most spheres
grpMax = grp(idx); % group number of the group with the most spheres
idxKeep = find(components == grpMax);
Data_agglo(~idxKeep,:) = [];
Data_agglo = Data_agglo(idxKeep,:);

R(~idxKeep,:) = [];
R = R(idxKeep,:);

%% To find the center of the agglomerates


%% 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(linspace(xmin,xmax,Nn), linspace(ymin,ymax,Nn), linspace(zmin,zmax,Nn));
[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];


%% 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
figure()
plot3(Xgrid(active),Ygrid(active),Zgrid(active),'o','MarkerFaceColor','red');
% hold on
% plot3(Xgrid(:),Ygrid(:), Zgrid(:),'.','MarkerFaceColor','blue');

axis equal
grid on

toc

Elapsed time is 0.368642 seconds.

active_dip = sum(sum(sum(active)));
volume = active_dip * dipS^3;

Before I generate the mesh, I want to find the center of the agglomerate whcih i generate

I also attach the sample file in .mat on how data for agglomerate looks like

Does anyone knows how it can be done?

[Walter Roberson]

Do you need the centroid, or the center of mass? (The two are the same if you have uniform density)

[Chris Dan]

All my spheres are of different diamters and they also overlap

[Chris Dan]

I think in this case center of mass and centroid will be differnet , my volume for each sphere is different but they are of the same material, in my case it is silica

[Walter Roberson]

If they overlap I think you need to subtract out the area of overlap, which could get a bit messy.

[Walter Roberson]

In some cases, it can be easier to "render" the spheres into an occupancy grid and work from that.

[Image Analyst]

Which is why, below, I suggested building a 3-D image and calling regionprops3. You might also call imfill depending on how you want to define centroid.

[Jon]

I think the approach also depends upon what is physically occurring when the spheres "overlap". Do they still maintain their spherical geometry and somehow have twice the density in the overlap region?

If so, it is simple you can just compute the center of mass using the same equations you would use if they were not overlapping.

If the do not maintain their spherical geometry with twice the density in the overlap, but the original spheres conserve their mass, they must assume some new, non spherical shape. In that case wouldn't you need to know what these shapes were in order to compute the overall center of mass?

[Walter Roberson]

I wonder if it can be treated as accretion around seeds? Morphological dilation.

I thought this was an aggregate situation though, like concrete? Because in aggregates the density is not typically uniform, instead having denser centers such as stone, with the interstices filled in with a different material.

Answer 1

[Walter Roberson]

Your proposal is not correct.

Consider the case of a circle of radius 1 at (0, 0) and a second circle of negligible radius at (10, 0). You propose to average the centers of the circles, which would be (0 0) (10 0) giving a centroid at (5 0)

Now replace the first circle with the approximation (-1,0), (0,-1), (1,0), (0,1) each of negligible radius. The top/bottom, left/right of the original. The centroid of the left stays the same, so you propose that the overall centroid should be the same. Sum of the x coordinates of the 5 circles is 10, sum of the y coordinates is 0, but the number of circles is 5 so the average would be (2,0) not (5,0)

You have to multiply the centroid by the area and divide by the total area. https://www.geogebra.org/m/nbexg5v9

Answer 2

[Image Analyst]

If you wanted to do it numerically instead of analytically by creating a 3-D volumetric image then you could easily get the centroid with a call to regionprops3.