How to Generate 100 Random points inside the figure?

10 次查看（过去 30 天）
Nathaniel Datinguinoo
编辑： Jan 2017-8-28
[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);
  5 个评论
显示 3更早的评论
Nathaniel Datinguinoo
Nathaniel Datinguinoo 2017-8-18
How would it be done in the codes? I'm not that knowledgeable about MATLAB.
Steven Lord
Steven Lord 2017-8-18
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?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Jan
Jan 2017-8-20
编辑：Jan 2017-8-20
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:
  2 个评论
John BG
John BG 2017-8-20
Jan Simon
this is John BG <mailto:jgb2012@sky.com jgb2012@sky.com>
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
Jan 2017-8-21
编辑：Jan 2017-8-28
@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')

请先登录，再进行评论。

更多回答（1 个）

John BG
John BG 2017-8-19
编辑：John BG 2017-8-20
Hi Nathaniel
this is John BG <mailto:jgb2012@sky.com jgb2012@sky.com>
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:jgb2012@sky.com jgb2012@sky.com>

类别

在 Help Center 和 File Exchange 中查找有关 Bounding Regions 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by