检查坐标是否在另一个矩阵中任何点的给定距离内

Question

在 Matlab 中，我有两个包含坐标的相当大的矩阵（A 和 B）。两条线代表 x 和 y，每列代表一个笛卡尔坐标 (x;y)。

现在我想将矩阵 B 中距离矩阵 A 中任意点小于 1（米）的所有点存储在一个新矩阵中。

我可以遍历所有数据，但是非常耗时（矩阵为 2x800000）。

有什么方法可以提高性能？

---

这是我当前代码的结构：

new_vec = [0;0];
for i=1:length(A)
    cur_x = A(1, i);
    cur_y = A(2, i);

    for j=1:length(B)
        if B(2, j) <= cur_y + 1 && B(2, j) >= cur_y - 1 && ...
               B(1, j) <= cur_x + 1 && B(1, j) >= cur_x - 1
            new_vec = [new_vec, [B(1, j); B(2, j)]];
        end
    end
end

Answer 1

另一种选择是使用pdist2，如下所示：

new_vec = B(:, any(pdist2(A', B', 'Chebychev') < 1, 1));

请注意，pdist2 总是比你的方法快，但可能比 Wolfie 的建议慢，因为 pdist2 总是计算 A 和 B 的所有点之间的所有距离。

比较

我会比较：

• 原文：您在回答中提供的代码
• 优化：Wolfie提供的代码
• pdist2：我使用pdist2 的解决方案
• bsxfun：rahnema1 的回答
• bsxfun (>=2016b)：rahnema1 使用 2016b 新功能的答案

使用以下示例数据

A = rand(2, N)*N*relativeAmplitude;
B = rand(2, N)*N*relativeAmplitude;

N 和relativeAmplitude=1 函数的执行时间：

relativeAmplitude 和N=10000 函数的执行时间：

结论

所有解决方案（Wolfie 的、rahnema1 的和我的）都比原始算法更快。

优化（Wolfie）与pdist2（我的）：如果B 的索引很可能会在A 中找到，那么 Wolfie 的答案可能会快 50 倍，但是如果不太可能，pdist 可能会快 50%。请注意，我的解决方案的执行时间与relativeAmplitude 无关，而 Wolfie 的则不是，但在某些情况下，Wolfie 的回答可能要快得多。

bsxfun (rahnema1) vs pdist2 (mine)：如果没有新的 R2016b 功能，bsxfun 总是比 pdist2 慢约 50%，否则这两种方法总是 (几乎）同样快。

Answer 2

性能改进基于您当前的实现：

% Appending is bad practise for memory management, you should initialise the 
% entire output array at first.
new_vec = NaN(size(B));
% You should not use i as a loop variable, since you are overwriting the default i=sqrt(-1)
% Also length(A)=max(size(A)), clearer to use size(A,2)

% Loops have been swapped as we want to exit the *A* looping when satisfied
for jj=1:size(B,2)
    % No need to re-assign current variables each loop, waste of time/memory

    % Same as before, j also is sqrt(-1) by default!
    % We could remove this loop entirely using vectorization, but it's likely quicker to
    % loop *until the condition is satisfied* then exit the loop early, avoiding many ops.
    for ii=1:size(A,2)
        % We can *half* the number of logical operations by using the abs distance!
        if abs(B(2,jj)-A(2,ii)) <= 1 && abs(B(1,jj) - A(1,ii)) <= 1
            % We pre-allocated, so no need to append - use direct indexing
            new_vec(:,jj) = B(:,jj);
            % Now the condition is satisfied for B(:,jj), exit the jj loop!
            break;
        end
    end
end
% We initialised an array of NaNs, remove the columns which are still NaN
new_vec(:, isnan(new_vec(1,:))) = [];

亮点：

• 最佳实践：不要使用i或j作为循环变量，它们的默认值为sqrt(-1)
• Pre-allocate memory, don't append results during looping。
• 减少逻辑检查次数，使用abs进行绝对距离检查。
• 不要在不必要的情况下为每个循环分配临时变量。
• 因为当给定的B 坐标在任意A 坐标的距离内时您会很高兴，所以尽早退出循环以避免进一步检查。

---

初始化NaNs 的二维数组的一种（可能对内存更友好的）替代方法是初始化一个布尔值false 数组，然后在每个满足的jj 索引处使其为真。最后我们会做new_vec = B(:,booleanVector);

Answer 3

这是一个矢量化的解决方案：

cur_x = A(1,:);
cur_y = A(2,:);
B1= reshape(B(1,:),[],1);
B2= reshape(B(2,:),[],1);
condition = abs(bsxfun(@minus,B2,cur_y))<=1 & ...
            abs(bsxfun(@minus,B1,cur_x))<=1;

[x ,~]=find(condition);
new_vec = [[0;0] B(:,x)];

从 MATLAB r2016b 开始，您可以将 condition 写为：

condition = abs(B2-cur_y)<=1 & ...
            abs(B1-cur_x)<=1;

或

condition = B2 <= cur_y + 1 & B2 >= cur_y - 1 & ...
               B1 <= cur_x + 1 & B1 >= cur_x - 1;

*abs(B(2,jj)-A(2,ii)) <= 1 的想法从@Wolfie 答案中窃取。