使用matlab实现指定三维区域随机生成N个不相交球的函数

自从写了这篇博客之后，更多的人来问我这个问题了。

简单写一个随机在一个三维区域生成球的matlab函数_风一样的航哥的博客-CSDN博客_comsol随机生成球体

上次简单写的功能，可以实现一些比较稀疏的分布，当分布比较密集的时候，会存在计算时间过长的问题，针对这个情况，本文做了更新。感谢粉丝的讨论和交流，让实现能够不断优化。

实现思路：1.随机生成，要用到随机数生成函数，matlab有好几种，这里实现使用了rand和randi函数，不同的是一个实现的是double类型，一个生成的是int类型的值。

2.位置随机，当然就采用了球心坐标是随机数，随机数在给定区间之内。

3.不相交，既然要不相交，那么任意两个球的距离都要大于球的直径，所以每次生成一个点之后，都要判断所有的距离是否满足条件，如果满足就加入到球心坐标数组中，不满足就下一次迭代。

4.如何画球体呢，本人使用了surf函数，surf(X,Y,Z) 创建一个三维曲面图，它是一个具有实色边和实色面的三维曲面。该函数将矩阵 Z 中的值绘制为由 X 和 Y 定义的 x-y 平面中的网格上方的高度。曲面的颜色根据 Z 指定的高度而变化。

下面是具体代码部分：

function  PlotNsphere_in_cubeRect(X1,Y1,Z1,X2,Y2,Z2,N,R)
%plotNcircle_in_rect(X1,X2,Z1,Y1,Y2,Z2,N,R) 函数说明：生成随机N个互不相交的球体
%  参数说明：三维区域为（X1，Y1，Z1）->（X2，Y2，Z2）
%  N：表示N个球体，R表示球体半径
%  输出是一张plot绘制的图，
close all
MAXIRTER = 10000; % 最大迭代次数
i=1;
s=[];%矩阵s，第一行存储x，第二行存储y，第三行存储z
c=1; %控制循环次数
rng('shuffle')          %%每次随机种子不一样，防止每次结果都是一样的
while c<MAXIRTER
    %% 如果采用以下三行代码，那么生成的随机数为整数
    tempx=randi([int32(X1),int32(X2)],1);
    tempy=randi([int32(Y1),int32(Y2)],1);
    tempz=randi([int32(Z1),int32(Z2)],1);
    %% 如果采用以下三行代码，那么生成的随机数为double类型，也即是默认的取值
%     tempx=X1+(X2-X1)*rand();
%     tempy=Y1+(Y2-Y1)*rand();
%     tempz=Z1+(Z2-Z1)*rand()
    %如果矩阵s为空矩阵，则存储第一个随机数值
    if(isempty(s)==1)
        s(1,i)=tempx;
        s(2,i)=tempy;
        s(3,i)=tempz;
        i=i+1;
        c=c+1;
    end
    %如果矩阵s不为空矩阵，求出各个坐标之间的距离矩阵disl，将与所有点之间距离都大于10的坐标存储
    if(isempty(s)==0)
        for j=1:i-1
            disl(1,j)=sqrt(((tempx-s(1,j))^2)+((tempy-s(2,j))^2)+((tempz-s(3,j))^2));
            j=j+1;
        end
        if(isempty(find(disl<2*R,1))==1)  % 或size(find(l<n),1)==0 表示矩阵l中小于n的个数等于0
            s(1,i)=tempx;
            s(2,i)=tempy;
            s(3,i)=tempz;
            i=i+1; 
        end
        c=c+1;
    end
    if(i>=N+1)
        break;
    end
end
X=s(1,:);
Y=s(2,:);
Z=s(3,:);

%% 绘制球体
    figure
for i=1:1:N
    myplotsphere(X(i),Y(i),Z(i),R);
    hold on 
    axis equal
    grid on
end
end

%% 绘制球体的子函数
function myplotsphere(X,Y,Z,R)
% myplotsphere使用三点和半径进行绘制一个球体
%   x,y,z;三维坐标；R：半径
[x,y,z]=sphere;
surf(R*x+X,R*y+Y,R*z+Z)
% mesh(R*x+X,R*y+Y,R*z+Z)
hold on
% axis equal
end

具体函数实现中，有注释，应该能看得懂。

实例：

当运行命令PlotNsphere_in_cubeRect(0,0,0,50+3.77,50+3.77,20+2,680,1.65)之后，等待一段时间会输出一图：

当运行PlotNsphere_in_cubeRect(0,0,0,50+3,50+3,20+2,1000,1.2)，也就是要生成1000个球体，会有结果：

 上述球体有点多，凑合看哈，最大测试就是1000啦，再大没有测试。

最后谢谢关注，如果问题，请私信或者留言。