求助：椭球随机分布程序

[fusijie1]

各位大神，谁可以指导小弟这个程序该怎么写呢，谢谢

[fusijie1]

谁能看看这个程序啊，是生成随机球体的，可总是不对，希望得到帮助，谢谢

xyzr=[];
n=10000;
while n
x=100*rand;
y=100*rand;
z=100*rand;
r=3+2*rand;
if ~isempty(xyzr)
d1=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y).^2+(xyzr(:,3)-z).^2);
d2=xyzr(:,4)+r;
if all(d1>=d2)
if x>=r && x<=100-r && y>=r && y<=100-r && z>=r && z<=100-r
xyzr=[xyzr;x y z r];
end
if x<r && y>=r && y<=100-r%x<rµÄÇé¿ö
d1=sqrt((xyzr(:,1)-x-100).^2+(xyzr(:,2)-y).^2+(xyzr(:,3)-z).^2);
d2=xyzr(:,4)+r;
%% if y>=r && y<=100-r %x<rʱ,y>=r
if z>=r && z<=100-r %
if all(d1>=d2)
xyzr=[xyzr;x y z r;x+100 y z r];
end
else if z<r
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y).^2+(xyzr(:,2)-z-100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x+100 y z r;x y z r];
end
else
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y).^2+(xyzr(:,2)-z+100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x+100 y z r;x y z-100 r];
end
end
elseif x<r && y<r %x<r,y<r
if z>=r && z<=100-r %
if all(d1>=d2)
xyzr=[xyzr;x y z r;x+100 y+100 z r];
end
else if z<r
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y-100).^2+(xyzr(:,2)-z-100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x+100 y+100 z r;x y z r];
end
else
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y+100).^2+(xyzr(:,2)-z+100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x+100 y z r;x y-100 z-100 r];
end
end
end

if x>100-r %x>100-r
d1=sqrt((xyzr(:,1)-x+100).^2+(xyzr(:,2)-y).^2+(xyzr(:,2)-z).^2);
d2=xyzr(:,4)+r;
if y>=r && y<=100-r
if z>=r && z<=100-r
if all(d1>=d2)
xzyr=[xyzr;x y z r;x-100 y z r];
end
elseif z<r
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y).^2+(xyzr(:,2)-z-100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x-100 y z r;x y z+100 r];
end
else
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y).^2+(xyzr(:,2)-z+100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x-100 y z r;x y z-100 r];
end
end
elseif y<r
if z>=r && z<=100-r
if all(d1>=d2)
xzyr=[xyzr;x y z r;x-100 y-100 z r];
end
elseif z<r
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y-100).^2+(xyzr(:,2)-z-100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x-100 y z r;x y+100 z+100 r];
end
else
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y+100).^2+(xyzr(:,2)-z+100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4)
xyzr=[xyzr;x y z r;x-100 y z r;x y-100 z-100 r];
end
end
end

else
n=n-1;
end
else
xyzr=[x y z r];
if x<r
xyzr=[xyzr;x+100 y z r];
end
if x>100-r
xyzr=[xyzr;x-100 y z r];
end
if y<r
xyzr=[xyzr;x y+100 z r];
end
if y>100-r
xyzr=[xyzr;x y-100 z r];
end
if z<r
xyzr=[xyzr;x y z+100 r];
end
if z>100-r
xyzr=[xyzr;x y z-100 r];
end
end

%%t=[0:0.1:2*pi 0]';
delta=pi/12;
theta=0:delta:pi;
phi=0:delta:2*pi;

x=sin(theta)*cos(phi)*xyzr(:,[1 4])';
y=sin(theta)*sin(phi)*xyzr(:,[2 4])';
z=cos(theta)*xyzr(:,[3 4])';
h=plot(x,y,z);
hold on
plot(xyzr(:,1),xyzr(:,2),xyzr(:,3),'.')
set(h,'color','r')
axis equal
axis([0 100 0 100 0 100])
hold off

[fusijie1]

现在可以生成随机点，可是应该怎么生成球体呢？


%%%%%%%%%%%%%%%%%%edited by wen in 16th, July, 2013%%%%%%%%%%%%%%%%%%%
clear all;clc
xyzr=[];
n=0;
while (n<10) %(1) while
n=n+1;
x=100*rand;
y=100*rand;
z=100*rand;
r=3+2*rand;
if ~isempty(xyzr) %(2.1)if
d1=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y).^2+(xyzr(:,3)-z).^2);
d2=xyzr(:,4)+r;
else %(2.1)else
d1=sqrt((x)^2+(y)^2+(z)^2);
d2=d1+r;
end %(2.1)end
if all(d1>=d2) %(2.2)if
if (x>=r && x<=100-r && y>=r && y<=100-r && z>=r && z<=100-r) %(3.1)if
xyzr=[xyzr;x y z r];
elseif (x<r && y>=r && y<=100-r)%x<r的情况 %(3.1)elseif
d1=sqrt((xyzr(:,1)-x-100).^2+(xyzr(:,2)-y).^2+(xyzr(:,3)-z).^2);
d2=xyzr(:,4)+r;
if (z>=r && z<=100-r) % %(4.1)if
if all(d1>=d2) %(5.1)if
xyzr=[xyzr;x y z r;x+100 y z r];
end %(5.1)end
elseif (z<r) %(4.1) elseif
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y).^2+(xyzr(:,2)-z-100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4) %(5.2)if
xyzr=[xyzr;x y z r;x+100 y z r;x y z r];
end %(5.2)end
else %(4.1)else
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y).^2+(xyzr(:,2)-z+100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4) %(5.3)if
xyzr=[xyzr;x y z r;x+100 y z r;x y z-100 r];
end %(5.3)end
end %(4.1)end
elseif (x<r && y<r) %x<r,y<r %(3.1)elseif
if (z>=r && z<=100-r) % %(4.2)if
if all(d1>=d2) %(5.4)if
xyzr=[xyzr;x y z r;x+100 y+100 z r];
end %(5.4)end
elseif z<r %(4.2)elseif
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y-100).^2+(xyzr(:,2)-z-100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4) %(5.5)if
xyzr=[xyzr;x y z r;x+100 y+100 z r;x y z r];
end %(5.5)end
else %(4.2)else
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y+100).^2+(xyzr(:,2)-z+100).^2);
d4=xyzr(:,4)+r;
if (all(d1>=d2) && all(d3>=d4)) %(5.6)if
xyzr=[xyzr;x y z r;x+100 y z r;x y-100 z-100 r];
end %(5.6)end
end %(4.2)end
elseif (x>100-r) %x>100-r %(3.1)elseif
d1=sqrt((xyzr(:,1)-x+100).^2+(xyzr(:,2)-y).^2+(xyzr(:,2)-z).^2);
d2=xyzr(:,4)+r;
if ((y>=r)&& (y<=(100-r))) %(4.3)if
if ((z>=r)&&(z<=(100-r))) %(5.7)if
if all(d1>=d2) %(6.1)if
xzyr=[xyzr;x y z r;x-100 y z r];
end %(6.1)end
elseif (z<r) %(5.7)elseif
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y).^2+(xyzr(:,2)-z-100).^2);
d4=xyzr(:,4)+r;
if (all(d1>=d2) && all(d3>=d4)) %(6.2)if
xyzr=[xyzr;x y z r;x-100 y z r;x y z+100 r];
end %(6.2)end
else %(5.7)else
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y).^2+(xyzr(:,2)-z+100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4) %(6.3)if
xyzr=[xyzr;x y z r;x-100 y z r;x y z-100 r];
end %(6.3)end
end %(5.7)end
elseif (y<r) %(4.3)elseif
if ((z>=r)&&(z<=(100-r))) %(5.8)if
if all(d1>=d2) %(6.4)if
xzyr=[xyzr;x y z r;x-100 y-100 z r];
end %(6.5)end
elseif z<r %(5.8)elseif
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y-100).^2+(xyzr(:,2)-z-100).^2);
d4=xyzr(:,4)+r;
if (all(d1>=d2))&&(all(d3>=d4)) %(6.6)if
xyzr=[xyzr;x y z r;x-100 y z r;x y+100 z+100 r];
end %(6.6)end
else %(5.8)else
d3=sqrt((xyzr(:,1)-x).^2+(xyzr(:,2)-y+100).^2+(xyzr(:,2)-z+100).^2);
d4=xyzr(:,4)+r;
if all(d1>=d2) && all(d3>=d4) %(6.7)if
xyzr=[xyzr;x y z r;x-100 y z r;x y-100 z-100 r];
end %(6.7)end
end %(5.8)end
end %(4.3)end
else %(3.1)else
n=n-1;
end %(3.1)end


else %(2.1)else
if (x<r)
xyzr=[xyzr;x+100 y z r];
elseif (x>100-r)
xyzr=[xyzr;x-100 y z r];
elseif (y<r)
xyzr=[xyzr;x y+100 z r];
elseif (y>100-r)
xyzr=[xyzr;x y-100 z r];
elseif (z<r)
xyzr=[xyzr;x y z+100 r];
elseif (z>100-r)
xyzr=[xyzr;x y z-100 r];
else
xyzr=[x y z r];
end
end %(2.1)end

end
%(1)end
xyzr
figure(1)


theta=[0:0.1:pi 0]';
phi=[0:0.1:2*pi 0]';
%%x=[ones(size(theta)) sin(theta)]*[ones(size(phi)) cos(phi)]*xyzr(:,:,[1 4])';
%%y=[ones(size(theta)) sin(theta)]*[ones(size(phi)) sin(phi)]*xyzr(:,:,[2 4])';
%%z=[ones(size(theta)) cos(theta)]*xyzr(:,:,[3 4])';
%%h=plot3(x,y,z);
%%hold on

plot3(xyzr(:,1),xyzr(:,2),xyzr(:,3),'.');
% set(h,'color','r');
axis equal
axis([0 100 0 100 0 100])
% hold off