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

1 view (last 30 days)
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');
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
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);
% hold on
% plot3(Xgrid(:),Ygrid(:), Zgrid(:),'.','MarkerFaceColor','blue');
axis equal
grid on
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
Walter Roberson on 7 Jun 2022
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.

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 3 Jun 2022
Edited: Walter Roberson on 3 Jun 2022
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.

More Answers (1)

Image Analyst
Image Analyst on 3 Jun 2022
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.

Community Treasure Hunt

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

Start Hunting!