Fill area with random circles having different diameters

89 Ansichten (letzte 30 Tage)
mikymike89
mikymike89 am 12 Jun. 2018
Kommentiert: RAJESH am 19 Jan. 2023
I should fill the area of a 500x500 square with random circles having random diameters between 10 and 50 (without overlap). Then, I need the output file of the generated coordinates. Can someone help me, please?
  2 Kommentare
KSSV
KSSV am 12 Jun. 2018
How many number of circles? Any further conditions?
Image Analyst
Image Analyst am 12 Jun. 2018
Is this homework? Do you want the circles in a digital image (if so, what resolution), or do you want them in an overlay (graphical plane) like you can do with plot(), scatter(), or viscircles()? Be sure to look at Walter's answer below.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Anton Semechko
Anton Semechko am 14 Jun. 2018
Hopefully this iteration of code (see below) is the last one.
% Unconstrained circle packing example. Only centroids constrained to the rectangle.
ab=[500 500]; % rectangle dimensions; [width height]
R_min=5; % minimum circle radius
R_max=25; % maximum circle radius
cnst=false;
[C,R]=random_circle_packing_rectangle(ab,R_min,R_max,cnst);
% Constrained circle packing example
ab=[500 500]; % rectangle dimensions
R_min=5; % minimum circle radius
R_max=25; % maximum circle radius
cnst=true;
[C,R]=random_circle_packing_rectangle(ab,R_min,R_max,cnst);
% Example of how to export data
FileName=sprintf('%s%sRandom circle packing data.csv',cd,filesep);
csvwrite(FileName,[C R]);
fprintf('Data saved to %s\n',FileName)
============================================================================================================
function [C,R]=random_circle_packing_rectangle(ab,R_min,R_max,cnst,vis)
% Random circle packing inside a rectangle.
%
% INPUT:
% - ab : 1-by-2 vector specify rectangle dimensions; ab=[width height].
% ab=[500 500] is the default setting. Coordinates of the
% lower-left and upper-right corners of the rectangle are
% (0,0) and (ab(1),ab(2)), respectively.
% - R_min : minimum circle radius. R_min=min(ab)/100 is the default setting.
% - R_max : maximum circle radius. R_max=min(ab)/20 is the default setting.
% - cnst : set cnst=true to ensure all circles fit into the rectangle.
% cnst=false is the default setting, meaning that only circle
% centroids will be constrained to the boundary and interior of
% the rectangle.
% - vis : set vis=false to suppress visualization. vis=true is the
% default setting.
%
% OUTPUT:
% - C : Q-by-2 array of sphere centroids
% - r : Q-by-1 array of sphere radii
%
% Default settings
if nargin<1 || isempty(ab)
ab=[500 500];
elseif numel(ab)~=2 || ~isnumeric(ab) || ~ismatrix(ab) || sum(ab(:)<1E-6 | isinf(ab(:)))>0
error('Invalid entry for 1st input argument (ab)')
else
ab=ab(:)';
end
if nargin<2 || isempty(R_min)
R_min=min(ab)/100;
elseif numel(R_min)~=1 || ~isnumeric(R_min) || R_min>(min(ab)/4-2E-12) || R_min<min(ab)/1E3
error('Invalid entry for 2nd input argument (D_min)')
end
if nargin<3 || isempty(R_max)
R_max=min(ab)/20;
elseif numel(R_max)~=1 || ~isnumeric(R_max) || R_max<R_min || R_max>(min(ab)/4-2E-12)
error('Invalid entry for 3rd input argument (D_max)')
end
if nargin<4 || isempty(cnst)
cnst=false;
elseif numel(cnst)~=1 || ~islogical(cnst)
error('Invalid entry for 4th input argument (cnst)')
end
if nargin<5 || isempty(vis)
vis=true;
elseif numel(vis)~=1 || ~islogical(vis)
error('Invalid entry for 5th input argument (vis)')
end
% Grid used to keep track of unoccupied space inside the rectangle
dx=max(min(ab)/2E3,R_min/50);
x=0:dx:ab(1);
y=0:dx:ab(2);
[x,y]=meshgrid(x,y);
G=[x(:) y(:)];
clear x y
% Start by placing circles along the edges if cnst=false
dR=R_max-R_min;
Cc=bsxfun(@times,ab,[0 0; 1 0; 1 1; 0 1]); % corner vertices
if ~cnst
[Xa,Xb]=deal(Cc,circshift(Cc,[-1 0]));
Rc=dR*rand(4,1)+R_min;
[Rc_a,Rc_b]=deal(Rc,circshift(Rc,[-1 0]));
[C,R]=deal(cell(4,1));
for i=1:4
[Ci,Ri]=SampleLineSegment(Xa(i,:),Xb(i,:),Rc_a(i),Rc_b(i),R_min,R_max);
Ci(end,:)=[];
C{i}=Ci;
Ri(end)=[];
R{i}=Ri;
end
C=cell2mat(C);
R=cell2mat(R);
% Update grid
for i=1:size(C,1), G=update_grid(C(i,:),R(i),G,R_min); end
else
% Remove all grid points less than R_min units from the boundary
G_max=G+R_min+1E-12;
G_min=G-R_min-1E-12;
chk_in=bsxfun(@le,G_max,ab) & bsxfun(@ge,G_min,[0 0]);
chk_in=sum(chk_in,2)==2;
G(~chk_in,:)=[];
clear G_max G_min chk_in
C=[]; R=[];
end
% Begin visualization
if vis
hf=figure('color','w');
axis equal
set(gca,'XLim',[0 ab(1)]+ab(1)*[-1/20 1/20],'YLim',[0 ab(2)]+ab(2)*[-1/20 1/20],'box','on')
hold on
drawnow
Hg=plot(G(:,1),G(:,2),'.k','MarkerSize',2);
t=linspace(0,2*pi,1E2)';
P=[cos(t) sin(t)]; % unit circle
if ~isempty(C)
for i=1:size(C,1)
Pm=bsxfun(@plus,R(i)*P,C(i,:));
h=fill(Pm(:,1),Pm(:,2),'r');
set(h,'FaceAlpha',0.25)
end
drawnow
end
end
f=5:5:100;
fprintf('Progress : ')
fprintf(' %-3u ',f)
fprintf(' (%%complete)\n')
fprintf(' ')
Ng=size(G,1);
f=size(G,1)-round((f/100)*Ng);
% Use rejection sampling to populate interior of the rectangle
flag=true;
n=0; cnt=0; m=0;
while ~isempty(G) && cnt<1E4 && (~vis || (vis && ishandle(hf)))
n=n+1;
% New circle
if flag && (cnt>500 || size(G,1)<0.95*Ng)
flag=false;
Rg=R_max*ones(size(G,1),1);
end
i=[];
if cnt<=500 && flag
X_new=ab.*rand(1,2); % centroid
else
i=randi(size(G,1));
X_new=G(i,:)+(dx/2)*(2*rand(1,2)-1);
X_new=min(max(X_new,[0 0]),ab);
if cnt>1E3
Rg(:)=max(0.95*Rg,R_min);
end
end
if isempty(i)
R_new=dR*rand(1)+R_min; % radius
else
R_new=(Rg(i)-R_min)*rand(1)+R_min;
end
% Check if the circle fits inside the rectangle when cnst=true
if cnst
X_new_max=X_new+R_new+1E-12;
X_new_min=X_new-R_new-1E-12;
chk_in=X_new_max<=ab & X_new_min>=[0 0];
if sum(chk_in)<2
cnt=cnt+1;
continue
end
end
% Does the new circle intersect with any other circles?
if ~isempty(C)
d_in=sqrt(sum(bsxfun(@minus,C,X_new).^2,2));
id=d_in<(R+R_new);
if sum(id)>0
cnt=cnt+1;
if ~isempty(i)
Rg(i)=min(0.95*Rg(i),min(0.99*(R_new+dx/2),R_max));
Rg(i)=max(Rg(i),R_min);
end
continue
end
end
% Accept new circle
cnt=0;
m=m+1;
C=cat(1,C,X_new);
R=cat(1,R,R_new);
[G,id]=update_grid(X_new,R_new,G,R_min);
if ~flag, Rg(id)=[]; end
% Visualization
if vis && ishandle(hf)
Pm=bsxfun(@plus,R_new*P,X_new);
h=fill(Pm(:,1),Pm(:,2),'r');
set(h,'FaceAlpha',0.25)
if mod(m,50)==0
delete(Hg)
Hg=plot(G(:,1),G(:,2),'.k','MarkerSize',2);
drawnow
end
end
% Progress update
if size(G,1)<=f(1)
f(1)=[];
fprintf('* ')
end
end
fprintf('\n')
% Show rectangle
if vis && ishandle(hf)
Cc=[Cc;Cc(1,:)];
plot(Cc(:,1),Cc(:,2),'--k','LineWidth',2)
delete(Hg)
end
if nargout<1, clear C R; end
function [G,id]=update_grid(X_new,R_new,G,R_min)
% Remove grid points within R_new+R_min units of new circle
D=sum(bsxfun(@minus,G,X_new).^2,2);
id=D<(R_new+R_min+1E-12)^2;
G(id,:)=[];
function [C,R]=SampleLineSegment(Xa,Xb,Ra,Rb,R_min,R_max)
% Place circles along line segment between points Xa and Xb
r=Xb-Xa;
L=norm(r);
r=r/norm(L);
dR=R_max-R_min;
C=Xa; R=Ra;
while true
R_new=dR*rand(1)+R_min;
C_new=C(end,:)+r*(R(end)+R_new+R_max*rand(1));
D=L - norm(C_new + r*(R_new+Rb) - Xa); % will there be enough space left for the end point with radius Rb?
if D<2*(R_min+1E-12)
C=cat(1,C,Xb);
R=cat(1,R,Rb);
break
else
C=cat(1,C,C_new);
R=cat(1,R,R_new);
end
end
  10 Kommentare
Philip Steinwascher
Philip Steinwascher am 12 Jan. 2023
This code is great. I was looking for something almost like that. How would you need to edit the code for fitting a variable number of different circles for which one sets the exact radii plus a fixed distance between them?
Walter Roberson
Walter Roberson am 12 Jan. 2023
Requiring a fixed distance between N-dimensional spheres in an N-dimensional space always has a maximum of (N+1) objects placed. After (N+1), you can place additional objects but they will not be the given distance from all other points.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (5)

Walter Roberson
Walter Roberson am 12 Jun. 2018
How do you know when to give up on the random placement?
N = 500;
minr = 10; maxr = 50;
maxtries = 500000;
intcoords = false;
sq = zeros(N, N, 'uint8');
[X, Y] = ndgrid(1:N, 1:N);
nc = 0;
iteration = 0;
trynum = 0;
maxgoodtry = 0;
fmt = 'iteration #%d, nc = %d, try #%d';
fig = gcf;
set(fig, 'Units', 'pixels', 'Position', [0 0 N+30, N+30]);
image(sq);
colormap(jet(2))
axis([0 N+1 0 N+1]);
title(sprintf(fmt, iteration, nc, trynum));
drawnow();
rx = []; ry = []; rr = [];
while trynum <= maxtries
iteration = iteration + 1;
if intcoords
r = randi([minr, maxr]);
cxy = randi([r+1, N-r], 1, 2);
else
r = minr + rand(1,1) * (maxr-minr);
cxy = r + 1 + rand(1,2) * (N - 2*r - 1);
end
mask = (X-cxy(1)).^2 + (Y-cxy(2)).^2 <= r^2;
if nnz(sq & mask) > 0
trynum = trynum + 1;
else
sq = sq | mask;
image(sq);
nc = nc + 1;
rx(nc) = cxy(1); ry(nc) = cxy(2); rr(nc) = r;
title(sprintf(fmt, iteration, nc, trynum));
drawnow();
maxgoodtry = max(maxgoodtry, trynum);
trynum = 1;
end
end
fprintf('finished at iteration %d, hardest success took %d tries\n', iteration, maxgoodtry);
fprintf('Number of circles: %d\n', length(rx));
  9 Kommentare
jaehee jeong
jaehee jeong am 11 Jun. 2020
@Walter
I really want to make it up to you.
Thank you
davide cademartori
davide cademartori am 19 Jan. 2023
Hello, thanks for all your work.
I have a couple of further questions related to the code:
1)Is it possible to implement the partial overlapping of the circles up to r/2 of the new circle introduced at every iteration?
2)Do you have any suggestion on how to implement a criterion which stops the surface filling when there is a continuous path of connected circles from the bottom to the top of rectangle/square?
Thanks again for all the work you have already done.

Melden Sie sich an, um zu kommentieren.


Anton Semechko
Anton Semechko am 12 Jun. 2018
Bearbeitet: Anton Semechko am 12 Jun. 2018
Here is a demo that uses rejection sampling and Delaunay triangulation
function [C,r]=random_circle_packing_demo
%
% - C : Q-by-2 array of sphere centroids
% - r : Q-by-1 array of sphere radii
%
% Problem settings
a=500; % size of the square
D_min=10; % min diameter of inscribed circle
D_max=50; % max diameter of inscribed circle
% Start by placing randomply spaced points along the boundary
d_min=D_min/a;
d_max=D_max/a;
d=d_max-d_min;
Xb=cell(4,1);
for i=1:4
t=[];
t_net=0;
while t_net<1
t=cat(1,t,d*rand(1)+d_min);
t_net=t_net+t(end);
end
t=t/t_net;
[t_min,id_min]=min(t);
while t_min<d_min
t_new=d*rand(1)+d_min;
if t_new>t_min
dt=t_new-t_min;
if t_new>t_min*(1+dt)
t(id_min)=t_new;
t=t/(1+dt);
[t_min,id_min]=min(t);
end
end
end
t=cat(1,0,t);
t(end)=[];
Xi=ones(numel(t),2);
if i==1
Xi(:,1)=a*cumsum(t);
Xi(:,2)=0;
elseif i==2
Xi(:,1)=a;
Xi(:,2)=a*cumsum(t);
elseif i==3
Xi(:,1)=a*(1-cumsum(t));
Xi(:,2)=a;
else
Xi(:,1)=0;
Xi(:,2)=a*(1-cumsum(t));
end
Xb{i}=Xi;
if i==1
figure('color','w')
axis equal
set(gca,'XLim',[0 a]+a*[-1/20 1/20],'YLim',[0 a]+a*[-1/20 1/20],'box','on')
hold on
drawnow
end
plot(Xi(:,1),Xi(:,2),'.k','MarkerSize',10)
end
Xb=cell2mat(Xb);
% Use rejection sampling to populate interior of the square
X=Xb;
D=D_max-D_min;
n=0; cnt=0;
while cnt<100 % terminate when no additional points can be inserted after 100 iterations
n=n+1;
%disp([n size(X,1) cnt])
% New point
X_new=a*rand(1,2);
% Randomly chosen distance threshold, D_thr, between D_min and D_max. The point
% process being simulated depends on how this parameter is chosen. If
% you want the spheres to be packed closer together, bias D_thr towards
% D_min.
D_thr=D*rand(1)+D_min;
% Reject new point if it is closer to X than D_thr
d_out=sum(bsxfun(@minus,X,X_new).^2,2);
if min(d_out)<D_thr^2
cnt=cnt+1;
continue
end
% Accept new point
cnt=0;
X=cat(1,X,X_new);
plot(X_new(:,1),X_new(:,2),'.k','MarkerSize',10)
if mod(n,10)==0, drawnow; end
end
% Delaunay triangulation
DT=delaunayTriangulation(X);
triplot(DT)
[C,r]=incenter(DT);
id=r<D_min/2 | r>D_max/2;
C(id,:)=[];
r(id)=[];
plot(C(:,1),C(:,2),'.r','MarkerSize',10)
N=size(C,1);
% Plot the spheres
figure('color','w')
axis equal
set(gca,'XLim',[0 a]+a*[-1/20 1/20],'YLim',[0 a]+a*[-1/20 1/20],'box','on')
hold on
t=linspace(0,2*pi,1E2)';
P=[cos(t) sin(t)];
for i=1:N
Pi=bsxfun(@plus,r(i)*P,C(i,:));
plot(Pi(:,1),Pi(:,2),'-k')
if mod(i,10)==0, drawnow; end
end
  4 Kommentare
Anton Semechko
Anton Semechko am 12 Jun. 2018
To better fill in the area, set 'D_thr' variable inside the while loop to D_min:
D_thr=D_min;
Simplest way to "export" results of the simulation is to call
[C,r]=random_circle_packing_demo;
from command window. Once the simulation is complete, right click on the 'C' variable (centroids) in the workspace, select "Copy", then paste "C" into spreadsheet. Repeat the same for 'r' variable (radii).
Anton Semechko
Anton Semechko am 12 Jun. 2018
Actually, setting D_thr to D_min will not produce tighter packing of the circles. After sampling is completed I think that the Voronoi cells corresponding to circles with radii below D_min will have to be modified locally using Lloyd's relaxation algorithm to get tighter packing. I don't really have time to experiment with this. Maybe someone else on here can help you.

Melden Sie sich an, um zu kommentieren.


Anton Semechko
Anton Semechko am 12 Jun. 2018
Here is another implementation of random circle packing. This version produces much tighter packing than my previous demo thanks to Walter's idea of using a grid to keep track of unoccupied space inside the square.
function [C,R]=random_circle_packing_demo2
%
% - C : Q-by-2 array of sphere centroids
% - r : Q-by-1 array of sphere radii
%
% Problem settings
a=500; % size of the square
D_min=10; % min circle diameter
D_max=50; % max circle diameter
% Unit circle
t=linspace(0,2*pi,1E2)';
P=[cos(t) sin(t)];
% Sampling grid
x=linspace(0,a,2E3);
dx=x(2)-x(1);
[x,y]=meshgrid(x);
G=[x(:) y(:)];
clear x y
% Begin visualization
hf=figure('color','w');
axis equal
set(gca,'XLim',[0 a]+a*[-1/20 1/20],'YLim',[0 a]+a*[-1/20 1/20],'box','on')
hold on
drawnow
f=5:5:100;
fprintf('Progress : ')
fprintf(' %-3u ',f)
fprintf(' (%%complete)\n')
fprintf(' ')
f=size(G,1)-round((f/100)*size(G,1));
% Use rejection sampling to populate interior of the square
C=[]; R=[];
D=D_max-D_min;
flag=true;
n=0; cnt=0; m=0;
while ~isempty(G) && ishandle(hf)
n=n+1;
% New circle
if cnt>500, flag=false; end
if cnt<500 && flag
X_new=a*rand(1,2); % centroid
else
i=randi(size(G,1));
X_new=G(i,:)+(dx/2)*(2*rand(1,2)-1);
end
R_new=(D*rand(1)+D_min)/2; % radius
% Does the new circle intersect with any other circles?
if n>1
d_in=sqrt(sum(bsxfun(@minus,C,X_new).^2,2));
id=d_in<(R+R_new);
if sum(id)>0
cnt=cnt+1;
continue
end
end
% Accept new circle
cnt=0;
m=m+1;
C=cat(1,C,X_new);
R=cat(1,R,R_new);
G=update_grid(X_new,R_new,G,D_min/2);
% Visualization
if ishandle(hf)
Pm=bsxfun(@plus,R_new*P,X_new);
h=fill(Pm(:,1),Pm(:,2),'r');
set(h,'FaceAlpha',0.25)
if mod(m,10)==0, drawnow; end
end
% Progress update
if size(G,1)<=f(1)
f(1)=[];
fprintf('* ')
end
end
fprintf('\n')
if nargout<1
clear C R
end
function G=update_grid(X_new,R_new,G,R_min)
% Remove grid points within R_new+R_min units of new circle
D=sum(bsxfun(@minus,G,X_new).^2,2);
id=D<(R_new+R_min+1E-12)^2;
G(id,:)=[];
  18 Kommentare
RAJESH
RAJESH am 19 Jan. 2023
@Anton Semechko, Could help me with Random Position of short Fiber(RVE) and Generate Fiber.

Melden Sie sich an, um zu kommentieren.


jahanzaib ahmad
jahanzaib ahmad am 16 Feb. 2019
Bearbeitet: jahanzaib ahmad am 16 Feb. 2019
can some one do such packing with polygons ? without using delauny triagulation
  7 Kommentare
Walter Roberson
Walter Roberson am 25 Feb. 2019
There is a notable difference between what the paper talks about and what you are talking about. In the paper, they are permitted to rotate the object at need, but your requirements appear to be that the polygons are to be at random orientations. The algorithms in the paper are permitted to use heuristics and computations to figure out where to best put a piece, whereas your requirements appear to be that pieces are added at random locations if they fit.
The paper deals with packing algorithms. Your question appears to deal with a situation similar to if you were to pour a bunch of objects into a box, remove the ones that overlap something beneath it, pour those removed ones on top, remove the ones that overlap, and so on, until finally no more pieces can possibly fit.
jahanzaib ahmad
jahanzaib ahmad am 25 Feb. 2019
yes.. Walter Roberson i cant see any material like polygon packing or random placement on mathworks .so i m doing that . i have made code just trying to improve it and i will post it ..
i m trying is randomly placing it and then checking overlap and placeing it to new position. but such a method is slow .can u please suggest any good way e.g i m trying now it that if a polygon has occupied some place in square i dont want to produce random number in that place so that new polygons has few chances to goto that place again so nimber of overlaps/intersections will reduce .? what do u say ? or any other method i can search for free space in a box and place new ploygono there in random way ?

Melden Sie sich an, um zu kommentieren.


Mohammad Faris
Mohammad Faris am 20 Feb. 2019
Hey everybody.
i have some issue in 3d dimension
i wanna fill a cylinder with some spherical glasses. can we rework this method for 3d dimension?
  7 Kommentare
Yi Liu
Yi Liu am 31 Dez. 2020
The following function does fill a cylinder space with radom sized ellipsoids without overlapping. The resulting ellipsoids are colour coded by their sizes (small-blue to large-orange). The following plots are generated by calling:
[C,R]=random_ellipsoid_packing_demo(500,600,0.3);
The function is:
function [C,R]=random_ellipsoid_packing_demo(d,h,ffactor)
% Randomly pack different size ellipsoids in a cylinder value with dim d,
% height h, and a filling factor ffactor (0.0 - 1.0).
%
% - C : Q-by-3 array of ellipsiod centroids
% - R : Q-by-3 array of ellipsiod radii
%
% Problem settings
%a=500; % size of the square
tvol = pi*(d/2)^2*h;
D_min=5; % min ellipsiod diameter
D_max=25; % max ellipsiod diameter
% Unit circle
t=linspace(0,2*pi,1E2)';
P=[cos(t) sin(t)];
% Unit ellipsiod (sphere)
[xs,ys,zs]=ellipsoid(0,0,0,1,1,1);
% Sampling grid
x=linspace(0,d,1E3);
y=linspace(0,d,1E3);
z=linspace(0,h,1E3);
dx=x(2)-x(1);
[x,y,z]=ndgrid(x,y,z);
G=[x(:) y(:) z(:)];
clear x y z
cmap=jet;
Ridx = linspace(D_min,D_max,256);
% Begin visualization
hf=figure('color','w');
plot3(d/2*P(:,1)+d/2,d/2*P(:,2)+d/2,h*ones(size(t)),':r')
hold on
plot3(d/2*P(:,1)+d/2,d/2*P(:,2)+d/2,0*ones(size(t)),':r')
axis equal
set(gca,'XLim',[0 d]+d*[-1/20 1/20],'YLim',[0 d]+d*[-1/20 1/20],'ZLim',[0 h]+h*[-1/20 1/20],'box','on')
drawnow
C=zeros(1,3); R=zeros(1,3);
D=D_max-D_min;
flag=true;
n=0; cnt=0; m=0;
fillfact=0;
fillper=0;
while ~isempty(G) && ishandle(hf) && fillfact<ffactor*tvol
n=n+1;
% New ellipsiod
if cnt>500, flag=false; end
if cnt<500 && flag
X_new=[d,d,h].*rand(1,3); % centroid
else
i=randi(size(G,1));
X_new=G(i,:)+[dx/2,dx/2,h/2].*(2*rand(1,3)-1);
end
R_new = D*rand(1,3)+D_min; % radius
% Does the ellipsiod intersect with boundary (cylinder)?
if ((X_new(1)-d/2-R_new(1))^2+(X_new(2)-d/2-R_new(2))^2) > (d/2)^2 ...
|| ((X_new(1)-d/2+R_new(1))^2+(X_new(2)-d/2+R_new(2))^2) > (d/2)^2 ... % partly outside boundary
|| (X_new(3)-R_new(3))<0 || (X_new(3)+R_new(3))>h % top or bottom
continue
end
% Does the new ellipsiod intersect with any other ellipsiods?
if n>1
d_in=sqrt(sum(bsxfun(@minus,C,X_new).^2,2));
id=d_in<(R+R_new);
if (sum(id)>0)
cnt=cnt+1;
continue
end
end
% Accept new ellipsiod
cnt=0;
m=m+1;
C=cat(1,C,X_new);
R=cat(1,R,R_new);
fillfact=fillfact+4/3*pi*R_new(1)*R_new(2)*R_new(3);
idx = find(mean(R_new)<Ridx,1);
ccmap=zs*0;
ccmap(:,:,1)=cmap(idx,1);
ccmap(:,:,2)=cmap(idx,2);
ccmap(:,:,3)=cmap(idx,3);
% Visualization
if ishandle(hf)
hh=surf(xs*R_new(1)+X_new(1),ys*R_new(2)+X_new(2),zs*R_new(3)+X_new(3),ccmap);
set(hh,'EdgeColor','None','FaceAlpha',0.3,'FaceColor',cmap(idx,:))
if mod(m,20)==0, drawnow; end
end
% Progress update
if mod(fillfact/tvol*100,5)<0.05 && ceil(fillfact/tvol*100)>fillper
fillper=ceil(fillfact/tvol*100);
fprintf('Progress (%% filled): %s\n',num2str(fillfact/tvol*100,'%2.1f'))
end
end
title(['Fill Factor = ',num2str(fillfact/tvol*100,'%.1f'),'%'])
figure
[np,bin]=hist(R(:),ceil(D));
yyaxis left
bar(bin,np)
ylabel('Number of Ellipsoids')
xlabel('Radius')
yyaxis right
plot(bin,np.*bin.^3*pi*4/3,'-*')
ylabel('Volume of Ellipsoids at Radius')
title('Statistics')
if nargout<1
clear C R
end
The resulting plots are:
Stats of the filling are:
RAJESH
RAJESH am 13 Jan. 2023
the above code is not running

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Graphics Performance finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by