创建一个大致等距点的球

问题描述

我需要创建一个球体，该球体在模拟中具有或多或少的等距点。
我看过一些代码，但是它们都只给出一个球体的表面，而我需要用从中心开始的等距点填充球。

之前，我创建了一个立方体，然后切出了半径，但是这给了我平坦的边缘，同时我需要使其更圆一些，因此实际上就像一个球。

如何在球中创建等距网格？

解决方法

我使用以下内容在球形壳体上创建近似等距点的网格。它基于适当的delta_theta的选择，即（共）横向箱宽度。然后将纵向phi宽度选择为delta_theta./sind(theta_c)，其中theta_c是分箱的中心纬度。换句话说：当离开（朝向）赤道时，垃圾箱的经度变宽（变小）。

注意：我使用的是姓氏和程度，即theta = [0 180)和phi = [0 360)。它的最终输出是两个变量：theta和phi。它们的长度相等，(theta(i),phi(i))是垃圾箱的中心点。

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中设置偏移量，以防止出现纯径向线。

,

对于采样体积，最直接的解决方案始终是规则的网格。随机选择点会导致采样效率降低（您需要更多点才能获得相同的表示，这是“立体学”的结果。）

但是立方网格并不是我们可以用来采样体积的唯一常规网格！晶体中的原子可以排列在各种不同的栅格中，最适合我们的是立方栅格，面心立方栅格（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