How to Generate 100 Random points inside the figure?

[x,y,z]=sphere;
mesh(x,y,z+9);

hold on;

% cylinder
r = 1.0;
h = 9.0;
m = h/r;
[R,A] = meshgrid(linspace(0,r,11),linspace(0,2*pi,41));
X = R .* cos(A);
Y = R .* sin(A);
Z = m*R;
% Cone around the z-axis, point at the origin
mesh(X,Y,Z)

hold on
axis equal

phi = -pi/3;
X1 = X*cos(phi) - Z*sin(phi);
Y1 = Y;
Z1 = X*sin(phi) + Z*cos(phi);

[x,y,z]=sphere;
mesh(x*0.6,y*0.6,z*0.6);

[Walter Roberson]

Which random distribution should be used? For example the slice at z = 6 has a larger circle than at z = 3, so if the distribution is to be according to slice cross-section or according to volume, more of the points would end up near the top.

[Nathaniel Datinguinoo]

We don't really mind if most of the points are on the top of the figure for as long as the 100 generated random points would be inside the figure.

[Walter Roberson]

rejection method might be easiest -- generate points within the bounding box, test if they are inside the object, accept if so.

[Nathaniel Datinguinoo]

How would it be done in the codes? I'm not that knowledgeable about MATLAB.

[Steven Lord]

Maybe build an alphaShape A from the sets of points that define your object then use inShape to determine if the randomly generated points are inside A?

Answer 1

[Jan]

Create the 3 shapes at first:

warning('off', 'MATLAB:alphaShape:DupPointsBasicWarnId');
Shape = cell(1, 3);

% Top sphere:
[x, y, z] = sphere;
Alpha     = 2;
Shape{1}  = alphaShape(x(:), y(:), z(:) + 9, Alpha);

% Cone:
r = 1.0;
h = 9.0;
m = h/r;
[R,A] = meshgrid(linspace(0,r,11), linspace(0,2*pi,41));
X2 = R .* cos(A);
Y2 = R .* sin(A);
Z2 = m * R;
Shape{2} = alphaShape(X2(:), Y2(:), Z2(:), Alpha);

% Bottom sphere:
Shape{3} = alphaShape(x(:) * 0.6, y(:) * 0.6, z(:) * 0.6, Alpha);
warning('on', 'MATLAB:alphaShape:DupPointsBasicWarnId');

% Surrounding cuboid:
coorMin   = [-1, -1, -0.6];  % Minimal x,y,z value
coorRange = [2, 2, 10.4];    % Max - Min of x,y,z value

Now choose random points inside the surrounding cuboid and collect the points, which are in any of the shapes:

% Get points inside the shapes:
maxIter = 1e6;  % Security limit
iter    = 0;
n       = 100;  % Number of output points
P       = zeros(n, 3);
iP      = 0;
while iP < n && iter <= maxIter
   Q = rand(1, 3) .* coorRange + coorMin;
   if inShape(Shape{1}, Q) || inShape(Shape{2}, Q) || inShape(Shape{3}, Q)
      iP       = iP + 1;
      P(iP, :) = Q;
   end
   iter = iter + 1;
end
if iter == maxIter  % Increase maxIter on demand:
   error('Cannot find enough points.');
end

This approach is simpler than creating one shape only, because you do not have to consider the overlap or know anything about the geometry or relative position of the objects. Even the sharp kink between the lower sphere and the cone is not problem here.

The smaller the volume of the shape containing the accepted points compared to the surrounding cuboid, the higher is the number of rejected points. Set maxIter accordingly.

Display the shape and the points:

% Show the shape and the points:
figure;
subplot(1,2,1, 'NextPlot', 'add');
for k = 1:3
   plot(Shape{k});
end
axis([-1,1, -1,1, -1,10], 'equal')
view(3)

subplot(1,2,2, 'NextPlot', 'add', 'XLim', [-1,1], 'YLim', [-1,1]);
plot3(P(:,1), P(:,2), P(:,3), '.');
axis([-1,1, -1,1, -1,10], 'equal')
view(3)

Here with n=1000 points:

[John BG]

Jan Simon

this is John BG <mailto:[e-mail] [e-mail]>

the lower half of the upper ball volume overlaps the inverted cone volume.

the tip of the cone, inside the lower ball, also overlaps volume with the lower ball volume.

If one really wants to scatter points uniformly, or with any other probability function, like it's mentioned if possible to scatter more points across higher levels, then, shouldn't all points of the volume have same weight?

The way you scatter the points you are giving double weight to that lower half of the upper ball, and the tip of the cone inside the lower ball.

Also, the way you scatter points, it may be possible that 2 or more random points are the physically the same point, because the points already selected are not removed from the pool.

I am sure you remember the discussion about similar problem, 2D scattering of large number of points in enclosed rectangle: once a point is selected as part of the random selection, it has to be removed from the pool for random selection to be consistent.

[Jan]

@John BG: Your arguments are not applicable here:

• The shown procedure accepts a point, if it is inside any of the 3 shapes. Then it does not matter if the shapes overlap: The points in the overlapping volume are not included multiple times and therefore there is no weighting effect.
• The original question asks for "random points". The author did not specify, if "random" includes or excludes multiple equal points. Remember that e.g. rand(1,2) produces equal values also with a chance of 2^-52 (this is 2.22e-16), and this is intended.
• Note that the chance to get the same random position in 3D is in the magnitude of 2^-104 / number of points. The code can run for trillions of years, before this happens.
• In opposite to this question, the mentioned 2D scattering problem required a minimal distance between the selected points. The solution posted there considers this restriction in the selection of points. If Nathaniel needs a unique set of points, he will mention this and it is very easy to consider:

if (inShape(Shape{1}, Q) || inShape(Shape{2}, Q) || inShape(Shape{3}, Q) && ...
   ~ismember(Q, P(1:iP, :), 'rows')

Answer 2

[John BG]

Hi Nathaniel

this is John BG <mailto:[e-mail] [e-mail]>

The following script, attached test_125.m, scatters points uniformly inside the volume enclosed by the figure you show in the the question

Although Stephen is right regarding that alphaShape and inShape being a way towards the solution, there's a problem when attempting to capture the points of this particular volume with a single alphaShape that I have solved using 2 alphaShapes.

1.

close all;clear all;clc

Na=41 % arc angle resolution
da=2*pi/Na;a=[3*pi/2:-da:pi/2]';
r01=.6;
ro1=r01*cos(a);
z1=r01*sin(a);
p1=[ro1 z1];
% figure(1);plot(p1(:,1),p1(:,2))   % check
% axis equal

r02=1;h1=9;
Nb=30;     % stick resolution
% dx1=.1;dz1=.2;p2=[[0:-dx1:-r2]' [0:dz1:h1]'];
p2=[linspace(0,-r02,Nb);linspace(0,h1,Nb)]';
ro2=p2(:,1);z2=p2(:,2);
% hold all
% figure(1);plot(p2(:,1),p2(:,2))   % check

% 1st correction: get intersection and trim
[ro_int1 z_int1]=intersections(p1(:,1),p1(:,2),p2(:,1),p2(:,2))
p1=[ro1(~and(ro1>=ro_int1,z1>=z_int1)) z1(~and(ro1>=ro_int1,z1>=z_int1))]
p2=[ro2(and(ro2<=ro_int1,z2>=z_int1)) z2(and(ro2<=ro_int1,z2>=z_int1))]

a2=[pi:-da:pi/2]';
ro3=r02*cos(a2);
z3=r02*sin(a2)+h1;

p3=[ro3 z3];
% hold all
% figure(1);plot(p3(:,1),p3(:,2))   % check

ro=[p1(:,1)' ro_int1 p2(:,1)' p3(:,1)'];
z=[p1(:,2)' z_int1 p2(:,2)' p3(:,2)'];
n_ro=numel(ro);

plot(ro,z)
axis equal

.

2.

generating the closed surface, not with 3 different shapes, but with cylindrical revolution of previous profile

% surface points
a3=[0:da:2*pi-da];  % phi revolving
z=repmat(z,Na,1);z=z(:);
ro=repmat(ro,Na,1);ro=ro(:);
ph=repmat(a3,n_ro,1)';ph=ph(:);
[xs,ys,zs]=pol2cart(ph,ro,z);  % surface points [x y z] format

% same problem chopped bat head when attempting to use alphaShape
shp1=alphaShape(xs,ys,zs);
figure(2);plot(shp1)

.

.

the problem is that despite alpha radius small enough, there's a horizontal plane splitting the bat cap

.

.

dr=.1 % 1D defining differential for Volume resolution
x_box_range=[-1.2:dr:1.2];
y_box_range=[-1.2:dr:1.2];
z_box_range=[-1:dr:11];

%%%%%%%%%%%%%%problem

[Xb,Yb,Zb]=meshgrid(x_box_range,y_box_range,z_box_range); % defining box points
in_shp1=inShape(shp1,Xb,Yb,Zb);  % inShape returning 1 across box volume points [Xb Yb Zb] within baseball bat surface

% figure(4);plot3(Xb,Yb,Zb,'go')

x1=Xb(in_shp1);y1=Yb(in_shp1);z1=Zb(in_shp1); % volume points, expected all, but top half ball points, most missing

% hold all
figure(3);plot3(x1,y1,z1,'b*')   % check
axis equal

%%%%%%%%%%%%%%

3.

way around this too few points in upper cap, with a 2nd alphaShape

.

% 2nd correction: getting the top half ball volume points
z0=9;  % plane halving top ball horizontaly
% radius 1st ball is 1
B1=71   
[X_mesh1,Y_mesh1,Z_mesh1]=sphere(B1);   
hm1=mesh(X_mesh1,Y_mesh1,Z_mesh1+z0);
X1=hm1.XData;Y1=hm1.YData;Z1=hm1.ZData;
figure(4);plot3(X1,Y1,Z1,'r*');axis equal;hold on
P1_3layer=zeros([size(X1),3]);     
P1_3layer(:,:,1)=X1;P1_3layer(:,:,2)=Y1;P1_3layer(:,:,3)=Z1;
% don't use reshape(P1_3layer,size(X1,1)*size(X1,2),3); 
P1=[];
for s=1:1:size(X1,1)
     for p=1:1:size(X1,2) 
         L0=squeeze(P1_3layer(s,p,:))';
         P1=[P1;L0];
     end
end
P1_below_z0=[];    % removing points 1st ball surface below z0=9
for k=1:1:length(P1)
    L1=P1(k,:);
    if L1(3)<z0-.1  % allow a bit of overlap otherwise alphaShape leaves empty entire 360  belt around upper ball
        P1_below_z0=[P1_below_z0 k];
    end
end
P1(P1_below_z0,:)=[];
x_1=P1(:,1);y_1=P1(:,2);z_1=P1(:,3);  

figure(4);plot3(x_1,y_1,z_1,'g*')  % check
% [th1,rh1,z1]=cart2pol(x_1,y_1,z_1);     % [x y z] 1st ball points to cylindrical format [theta rho z]

shp2=alphaShape(P1); 
figure(5);plot(shp2);
shp2.Alpha  % jagged because this radius=1
shp2.Alpha=2
figure(5);plot(shp2);
% volume

dr=.1 % 1D defining differential for Volume resolution
x_box_range=[-1.2:dr:1.2];
y_box_range=[-1.2:dr:1.2];
z_box_range=[8:dr:11];

[Xbox,Ybox,Zbox]=meshgrid(x_box_range,y_box_range,z_box_range);

P2=zeros(3,numel(x_box_range)*numel(y_box_range)*numel(z_box_range));  

while k<=numel(x_box_range)*numel(y_box_range)*numel(z_box_range)
    for s=1:1:numel(x_box_range)
        for p=1:1:numel(y_box_range)
            for t=1:1:numel(z_box_range)
                P2(:,k)=[x_box_range(p);y_box_range(s);z_box_range(t)];
                k=k+1;
            end
        end
    end
end

xbox=P2(1,:);ybox=P2(2,:);zbox=P2(3,:);

in_shp2=inShape(shp2,xbox,ybox,zbox);

x_top=xbox(in_shp2);y_top=ybox(in_shp2);z_top=zbox(in_shp2);
figure(7);plot3(x_top,y_top,z_top,'b*') 
axis equal

% all points in same matrix
x_vol=[x1' x_top];y_vol=[y1' y_top];z_vol=[z1' z_top];
Pvol=[x_vol;y_vol;z_vol];

figure;plot3(Pvol(1,:),Pvol(2,:),Pvol(3,:),'c*')
axis equal

Pv=Pvol;  % copy to gradually remove selected points

% scattering points inside volume
Pscatter=[];
Nscatter=100   
s=1;
while s<=Nscatter
    nx=randi(length(Pv),1,1);
    L4=[Pv(1,nx);Pv(2,nx);Pv(3,nx)];
    Pscatter=[Pscatter L4];
    Pv(:,nx)=[];   % remove selected point
    s=s+1;
end

xscatter=Pscatter(1,:);yscatter=Pscatter(2,:);zscatter=Pscatter(3,:);
figure;plot3(xscatter,yscatter,zscatter,'r*') % check
axis  equal

.

Nathaniel, thanks for your attention

if you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?

To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link

thanks in advance

John BG

<mailto:[e-mail] [e-mail]>