Эффективно векторизовать поэлементную операцию в Matlab

Question

У меня есть матрица nx4 A, представляющая n сфер, и матрица mx3 B, представляющая m точек. Мне нужно проверить, находятся ли эти 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 * -дерева, но кажется, что реализация довольно сложная XD

Answer 1

При этом используется неявное расширение , чтобы вычислить все расстояния между точками и центрами сфер, а затем сравнить их с радиусами сфер:

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

Это, вероятно, быстрее, чем циклический подход, но также требует больше памяти, поскольку вычисляется промежуточный массив m × n × 3.