你可以做很多事情来改进你的程序——算法和代码。
在代码方面，真正让你慢下来的一件事是，你不仅使用for循环（这很慢），而且在行中
P = [P;P1];
您将元素附加到数组。每次发生这种情况时，Matlab 都需要找到一个新的地方来放置数据，复制过程中的所有点。这很快变得非常缓慢。预分配数组
P = zeros(1000000, 3);
跟踪到目前为止您找到的点数 N，并将距离计算更改为
D = pdist2(P1, P(1:N, :), 'euclidean');
至少会解决这个问题......
另一个问题是您检查所有先前找到的点的新点 - 因此，当您有 100 个点时，您检查大约 100x100，对于 1000，它是 1000x1000。然后你可以看到这个算法至少是 O(N^3) ......不计算随着密度的增加你会得到更多“未命中”的事实。N=10E6的AO(N^3)算法至少需要10E18个周期；如果您有一台 4 GHz 机器，每次比较有一个时钟周期，您将需要 2.5E9 秒 = 大约 80 年。您可以尝试并行处理，但这只是蛮力 - 谁想要呢？
我建议您考虑将问题分解成更小的部分（从字面上看）：例如，如果您将球体分成大约为最大距离大小的小盒子，并且对于每个盒子，您跟踪哪些点在它，那么您只需要检查“这个”框及其直接邻居中的点 - 总共 27 个框。如果您的盒子的宽度为 2.5 毫米，那么您将拥有 100x100x100 = 1M 的盒子。这看起来很多，但现在您的计算时间将大大减少，因为（在算法结束时）每个盒子平均只有 1 个点......当然，使用您使用的距离标准，您将在中心附近有更多点，但这是一个细节。
您需要的数据结构将是一个 100x100x100 的单元格数组，每个单元格都包含迄今为止“在该单元格中”找到的好点的索引。元胞数组的问题在于它不适合向量化。相反，如果您有内存，则可以将其分配为 10x100x100x100 的 4D 数组，假设每个单元格不超过 10 个点（如果这样做，则必须单独处理；在这里与我合作...） . 使用-1尚未找到的点的索引
那么您的支票将是这样的：
% initializing:
bigList = zeros(10,102,102,102)-1; % avoid hitting the edge...
NPlist = zeros(102, 102, 102); % track # valid points in each box
bottomcorner = [-25.5, -25.5, -25.5]; % boxes span from -25.5 to +25.5
cellSize = 0.5;
.
% in your loop:
P1= [x, y, z];
cellCoords = ceil(P1/cellSize);
goodFlag = true;
pointsSoFar = bigList(:, cellCoords(1)+(-1:1), cellCoords(2)+(-1:1), cellCoords(3)+(-1:1));
pointsToCheck = find(pointsSoFar>0); % this is where the big gains come...
r=sum(P1.^2);
D = pdist2(P1,P(pointsToCheck, :),'euclidean'); % euclidean distance
if D>0.146*r^(2/3)
P(k,:) = P1;
% check the syntax of this line...
cci = ind2sub([102 102 102], cellCoords(1), cellCoords(2), cellCoords(3));
NP(cci)=NP(cci)+1; % increasing number of points in this box
% you want to handle the case where this > 10!!!
bigList(NP(cci), cci) = k;
k=k+1;
end
....
我不知道你能不能从这里拿走；如果你不能，请在笔记中说明，我可能在这个周末有时间更详细地编写代码。有一些方法可以通过一些矢量化来加快速度，但很快就会变得难以管理。
我认为在空间中随机放置大量点，然后使用上面的进行巨大的矢量化剔除，可能是要走的路。但我建议先采取一些小步骤......如果你能让上述工作正常，你可以进一步优化（数组大小等）。
