提高MATLAB逻辑索引的性能

Question

我正在编写一个函数，该函数从大型坐标集中选择所有在一条线周围的边界框内的点（以加速RANSAC拟合）。昨晚我运行了包含该函数的脚本，在总共55k秒的时间内，有30k秒（在55M次调用中）花费在以下代码行上（X是3×N笛卡尔坐标，xmin、xmax等是与X长度相同且对应的边界向量）：

inlierIdx = ( X(1,:) >= xmin & X(1,:) <= xmax & X(2,:) >= ymin & X(2,:) <= ymax );

你能帮我加速这个吗？短路计算可以很好地提高性能，但似乎无法与索引一起使用。

如果有完全不同的更好方法，这是函数代码的其余部分（p1，p2是定义线的点，a是边界框的偏移量；额外步骤需要另外12k秒钟处理5500万次调用）：

function [inlierX, xmin, xmax, ymin, ymax] = selectbox(X, L, a) % xmin ... ymax only for plotting purposes
% X: 3xN coords to check; L: 3x2 vector defining line; a: offset of bounding box
p1 = L(:,1);
p2 = L(:,2);

constfactor = (X(3,:)-p1(3))/(p2(3)-p1(3)); % precompute for following steps
xmin = p1(1) + (p2(1)-p1(1))*constfactor - a; % line parallel to x-component of defined line, with offset
xmax = xmin + 2*a;
ymin = p1(2) + (p2(2)-p1(2))*constfactor - a;
ymax = ymin + 2*a; 

inlierIdx = ( X(1,:) >= xmin & X(1,:) <= xmax & X(2,:) >= ymin & X(2,:) <= ymax );

if p1(3) == p2(3) % singularity if line is horizontal, discard then
    inlierIdx = [];
end

inlierX = X(:,inlierIdx); % save & return inlier coordinates

我突然想到应该将那个奇异点if语句移动，这样如果条件成立就可以跳过计算，但这会稍微降低性能。

该函数针对每个RANSAC样本进行调用，以选择仅距离采样线合理距离的点来计算它们的距离以进行阈值处理。

MATLAB版本为R2016a。

并行计算工具箱可用，我尝试将步骤移入arrayfuns，并使用gpuArrays调用整个函数，但速度要慢得多。

整个上下文如下：

prepare coordinates

while trialcount < expectedNumTrialsNeeded

draw random sample (generated array of randomly sampled coordinates and index into it)

check if sample is not degenerate and has allowed angle

compute number of inliers:

    this calls my selectbox function, to select a smaller set of points, which are near enough to check - takes 30k sec of 55k

    take those points and compute their distance to the line, all points within threshold are inliers - takes 12k secs of 55k

if number of inliers > best number of inliers yet, this is new best set

update expected number of trials needed to find sample of only inliers
increment trialcount

end

return best set

我处理包含1万到10万个点的数组，对于每条线我想要拟合，我运行1E+5...1E+6次试验来找到最佳组合。我有大约50条线需要进行拟合，通过删除找到的线上的异常值并在其余数据上运行算法来完成。

- LimaKilo

---

你能发一些更多的代码吗？有时候分析器会有点误导人，特别是当循环没有JIT编译时...另外，你使用的MATLAB版本是什么？ - Rody Oldenhuis

@Rody 已完成，感谢您的关注。 - LimaKilo

更多的代码？比如说，这是一个函数还是脚本？是什么在调用这个函数？只有 inlierX 重要，还是 inlierIdx 也很重要？请把你所拥有的一切都发一个最小化可运行示例！ :) - Rody Oldenhuis

@Rody 完成了，请再次检查 :) - LimaKilo

不确定这会产生多大影响，你可以检查 nargout，如果你只需要 inlierX，那么不要实例化 xmin 等（我假设你并不总是绘图）。 - Steve Heim

@Steve：看起来没有什么效果，不过还是谢谢。 - LimaKilo

Answer 1

你的代码中是否考虑使用列索引？我的意思是将X作为Nx3的笛卡尔坐标系，而不是3xN。据我所知，由于MATLAB数组以FORTRAN方式存储，因此使用列索引可能会更快。

- Alex Chudinov

---

我曾考虑过按列索引，但直觉告诉我，3xN 应该使用按列索引来寻址单独的列。所以实际上是相反的，按列索引意味着在一列内进行索引？ - LimaKilo

虽然这种方法的加速通常不是惊人的 - 大约只有1.5倍左右的因素。 - Rody Oldenhuis

@RobyOldenhuis 同意。关于线性索引，您有什么看法？在这种情况下，另一种我会使用的方式是 C++ dll 或 mexfunction。 - Alex Chudinov

Answer 2

我怀疑那一行并不是真正的问题所在。你说你要调用这个函数5500万次，我们能看到做这件事的代码吗？因为我强烈怀疑问题实际上出在这段代码中。

我怀疑你在循环中调用了这个函数，这将使得MATLAB无法有效地使用其JIT编译器加速调用。如果是这种情况，那么确实，显示的最大负担就是这一行，但它之所以如此缓慢，是因为你的代码结构迫使MATLAB继续调用解释器，而不是执行机器代码。如果确实是这种情况，那么将这个小函数直接复制粘贴到循环体中将极大地加速所有操作（前提是循环体的其余部分也可以被加速...）

总之，这就是我能够理解的内容。它减少了比较的次数，但增加了索引的数量；我也怀疑这样做是否更快...好吧，只需进行1000次调用并进行比较即可。

function [inlierX, xmin, xmax, ymin, ymax] = selectbox(X, L, a)

    % X: 3xN coords to check; L: 3x2 vector defining line; a: offset of bounding box
    p1 = L(:,1);
    p2 = L(:,2);

    constfactor = (X(3,:)-p1(3)) / (p2(3)-p1(3));

    % line parallel to x-component of defined line, with offset
    xmin = p1(1) + (p2(1)-p1(1))*constfactor - a; 
    xmax = xmin + 2*a;
    ymin = p1(2) + (p2(2)-p1(2))*constfactor - a;
    ymax = ymin + 2*a;

    % singularity if line is horizontal, discard then
    if p1(3) ~= p2(3)        
        inlierIdx            = X(1,:)         >= xmin;        
        inlierIdx(inlierIdx) = X(1,inlierIdx) <= xmax(inlierIdx);
        inlierIdx(inlierIdx) = X(2,inlierIdx) >= ymin(inlierIdx);
        inlierIdx(inlierIdx) = X(2,inlierIdx) <= ymax(inlierIdx);                
    else
        inlierIdx = [];
    end

    inlierX = X(:,inlierIdx); % save & return inlier coordinates

end

- Rody Oldenhuis

---

你的代码明显比我的原始代码慢得多，210秒对比60秒。整个数组上的第一个布尔运算花费了5秒，其他的花费了约45秒，因此索引似乎是问题所在。 - LimaKilo

@LukasK. 你试过将函数体直接复制粘贴到循环体中吗？ - Rody Oldenhuis

您IP地址为143.198.54.68，由于运营成本限制，当前对于免费用户的使用频率限制为每个IP每72小时10次对话，如需解除限制，请点击左下角设置图标按钮（手机用户先点击左上角菜单按钮）。 - LimaKilo

@LukasK. 嗯，太好了！你无意中解决了一个你不知道存在的问题 :) 实际上，你确实改变了一些东西 - 你记得引入奇异条件了吗？ - Rody Oldenhuis

@LukasK。好的，很高兴能帮忙。我仍然觉得我可以在这里有用，所以当他回来并同意分享一些玩具数据集和代码时，请让我知道 :) - Rody Oldenhuis