在 matlab 中有效地矢量化逐元素操作

Question

我有一个代表 n 个球体的 nx4 矩阵 A 和一个代表 m 个点的 mx3 矩阵 B。我需要测试这些 m 点是否在任何球体内。我可以使用for 循环来做到这一点，但是对于较大的n 和m，这种方法效率非常低。如何矢量化此操作？我目前的方法是

A = [0.8622    1.1594    0.7457    0.6925;
     1.4325    0.2559    0.0520    0.4687;
     1.8465    0.3979    0.2850    0.4259;
     1.4387    0.8713    1.6585    0.4616;
     0.2383    1.5208    0.5415    0.9417;
     1.6812    0.2045    0.1290    0.1972];

B = [0.5689    0.9696    0.8196;
     0.5211    0.4462    0.6254;
     0.9000    0.4894    0.2202;
     0.4192    0.9229    0.4639];

for i=1:size(B,1)

    mask = vecnorm(A(:, 1:3) - B(i,:), 2, 2) < A(:, 4);

    if sum(mask) > 0
        C(i) = true;
    else
        C(i) = false;
    end %if

end %for

我测试了@LuisMendo 建议的方法，它似乎只加快了相当小的m 和n 的计算速度，但是对于大的m 和n，比如说，我的问题大约是10000 ，改善非常有限。但是@NickyMattsson 给了我一些提示。因为matlab中的逻辑运算比vecnorm要快，所以我先用粗检，找到点附近的球体，然后细检：

A = [0.8622    1.1594    0.7457    0.6925;
     1.4325    0.2559    0.0520    0.4687;
     1.8465    0.3979    0.2850    0.4259;
     1.4387    0.8713    1.6585    0.4616;
     0.2383    1.5208    0.5415    0.9417;
     1.6812    0.2045    0.1290    0.1972];

B = [0.5689    0.9696    0.8196;
     0.5211    0.4462    0.6254;
     0.9000    0.4894    0.2202;
     0.4192    0.9229    0.4639];

ids = 1:size(A, 1);

for i=1:size(B,1)

    % first a rough check
    xbound = abs(A(:, 1) - B(i, 1)) < A(:, 4);
    ybound = abs(A(:, 2) - B(i, 2)) < A(:, 4);
    zbound = abs(A(:, 3) - B(i, 3)) < A(:, 4);
    nears = ids(xbound & ybound & zbound);
    if isempty(nears)
        C(i) = false;
    else 
        % then a fine check
        mask = vecnorm(A(nears, 1:3) - B(i,:), 2, 2) < A(nears, 4);

        if sum(mask) > 0
            C(i) = true;
        else
            C(i) = false;
        end 
    end

end 

这可能会将时间减少到 1/2 或 1/3，这是可以接受的，如果我将 m 和 n 分成批次，它可能会更快，而不会造成太大的内存负担。 @CrisLuengo 提到了 R*-tree 的方法，但是实现起来好像挺复杂的XD

Answer 1

这使用implicit expansion 计算点和球心之间的所有距离，然后将它们与球半径进行比较：

C = any(vecnorm(permute(B, [1 3 2]) - permute(A(:,1:3), [3 1 2]), 2, 3) < A(:,4).', 2);

这可能比循环方法更快，但也更占用内存，因为会计算一个中间 m×n×3 数组。