我需要为我的模拟创建一个具有或多或少等距点的球体。
我查看了一些代码，但它们都只给出了一个球体的表面，而我需要用从中心开始的等距点填充球。
之前，我创建了一个立方体，然后切掉了半径，但这给了我平坦的边缘，而我需要它更圆，所以它实际上类似于一个球。
如何在球中创建等距网格？
【问题讨论】：
我使用以下方法在球壳上创建一个近似等距点的网格。它基于选择适当的delta_theta，即（共）纬度 bin 宽度。然后将纵向宽度phi 选择为delta_theta./sind(theta_c)，其中theta_c 是箱的中心纬度。换句话说：当远离（朝向）赤道时，箱子的经度会变宽（变小）。
注意：我使用纬度和度数，即theta = [0 180) 和phi = [0 360)。它的最终输出是两个变量：theta 和 phi。它们的长度相等，(theta(i),phi(i)) 是 bin 的中心点。
delta_theta = 1; % latitudinal bin width
theta_c = 0:delta_theta:89; % Set up bin centres
theta_c = [-fliplr(theta_c) theta_c].'+90; % Centre the centres on the equator
delta_phi = delta_theta./sind(theta_c); % Grab longitudinal bin width per latitudinal bin
bins_per_theta = round(360./delta_phi); % Ensure an integer number of bins
theta = zeros(sum(bins_per_theta),1);
phi = theta;
kk=1;
for ii=1:numel(theta_c) % We need a loop here to deal with the changing number of bins
phi_c = linspace(0,360,bins_per_theta(ii));
kk_add = numel(phi_c);
theta(kk:kk+kk_add-1) = theta_c(ii);
phi(kk:kk+kk_add-1) = phi_c;
kk = kk+kk_add;
end
在此代码上进行一个简单的循环以根据半径适当减小delta_theta 会将其填充到一个完整的球。我会说delta_theta = 1/r，具有适当的常数，会随着半径的增加而减小纬度箱宽度。也许您还需要在 theta_c 和 phi_c 中进行偏移，以防止出现纯径向线。
【讨论】：
plot3 这样具有一些固定半径的例程，然后给出纬度和经度（可能必须将它们转换为弧度）。对于一个球，应该在柱子末端提到的延伸；将整个代码包装在半径r 上的循环中，并根据一些适当的标准更改delta_theta。再说一遍，您将拥有三个数组：r, theta, phi，它们是球坐标，可以这样绘制。
对于体积采样，最直接的解决方案始终是规则网格。随机选择点会导致采样效率降低（您需要更多点才能获得相同的表示，这是 Stereology 的结果）。
但是立方网格并不是我们可以用来对体积进行采样的唯一规则网格！晶体中的原子可以排列成各种不同的网格，最适合我们的是立方网格、面心立方网格 (FCC) 和体心立方网格 (BCC)——其他网格的各向同性较小。已经发表了许多关于使用这些网格对体积图像进行采样的研究论文，事实证明，在对体积进行采样时，FCC 和 BCC 网格的效率明显高于立方网格。这意味着，需要更少的样本来表示具有相同误差的体积。请注意，例如，水果供应商将使用 FCC 或 BCC 网格堆叠橙子和苹果，因为它们会导致尽可能密集的球体堆积。
是使用 FCC 还是 BCC 网格采样，还是像 Adriaan's 这样的方法取决于最重要的是：对体积本身进行良好采样，还是对体积边缘进行良好采样。
简单的三次方格由以下矩阵表示：
1 0 0
0 1 0
0 0 1
FCC 网格由以下矩阵表示：
1 1 0
1 0 1
0 1 1
BCC 网格由以下矩阵表示：
1 1 -1
1 -1 1
-1 1 1
这些矩阵的行代表网格的单位向量。任何网格点的坐标都可以通过矩阵乘以索引向量得到：
M * [i;j;k]
要找到球中的点集，我们可以这样做：
d = 0.3; % grid spacing
r = 2.3; % ball size
M = [1 1 0
1 0 1
0 1 1]; % FCC grid
p = -2*r/d : 2*r/d; % take a larger area around the ball, so we're sure to sample the whole thing
[i, j, k] = ndgrid(p, p, p);
p = [i(:).'; j(:).'; k(:).'];
p = (d * M) * p;
I = sqrt(sum(p.^2, 1)) < r; % can use `vecnorm` in newer MATLAB releases
p = p(:, I);
scatter3(p(1,:), p(2,:), p(3,:))
axis equal
【讨论】：