这段代码可以进一步矢量化以消除循环吗？

Question

我正在研究MATLAB中的光线跟踪几何问题,并且在我的程序中遇到了瓶颈.

该函数接收光线的起点和终点(lStart和lEnd),一组平面点和法线(pPoint和规范).然后,该函数计算沿着与每个平面相交的光线的距离.

以下是对原始等式的引用:https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection#Algebraic_form

我到目前为止的代码如下:

dists = (diag(bsxfun(@minus, pPoint, lStart) * norms')) ./ ((lEnd - lStart) * norms')';

这是在循环中调用的,例如:

nRays   = size(lStart, 1);
nPlanes = size(pPoint, 1);
dists   = zeros(nPlanes, nRays);

for rayCtr = 1:nRays

    dists(:, rayCtr) = (diag(bsxfun(@minus, pPoint, lStart(rayCtr, :)) * norms')) ./...
         ((lEnd(rayCtr, :) - lStart(rayCtr, :)) * norms')';

end

这对于单个光线非常有效.

给定一个射线[1 x 3]和300个平面[300 x 3],我得到一个[300 x 1]矩阵,每个平面交点的距离.

我正在努力的是,将其矢量化以在光线列表上工作.

典型数据集中的大小为:

lStart, lEnd  = [14e6, 3];
pPoint, norms = [300,  3];

光线处理通常分为数万个 - 以适应内存.对于每个批次,我想消除rayCtr循环.使用此方法,整个程序只需8个小时(32位,Windows,2GB RAM).

以下是六个平面(形成一个长方体)和两个光线作为MWE的一些坐标:

pPoint = [-0.5000   -0.5000   -0.5000;
          -0.5000   -0.5000    0.5000;
          -0.5000   -0.5000   -0.5000;
          -0.5000    0.5000   -0.5000;
          -0.5000   -0.5000   -0.5000;
           0.5000   -0.5000   -0.5000]

norms = [0  0   1;
         0  0   1;
         0  1   0;
         0  1   0;
         1  0   0;
         1  0   0]

lStart = [-1 0 0;
          -1 0.25 0]

lEnd   = [1 0 0;
          1 0.25 0]

该示例的预期输出是:

dists = [-Inf -Inf;
          Inf  Inf; 
         -Inf -Inf; 
          Inf  Inf; 
          0.25 0.25; 
          0.75 0.75]

非常感谢.

更新:通过在接受的答案中提出的解决方案,运行时间缩短到大约30分钟 - 现在受到读写操作和体素查找的限制.

Answer 1

我想你需要的是

dists=sum(bsxfun(@times,bsxfun(@minus,...
                               permute(pPoint,[1 3 2]),permute(lStart,[3 1 2])),...
                 permute(norms,[1 3 2])),3)...
        ./(sum(bsxfun(@times,...
                      permute(lEnd-lStart,[3 1 2]),permute(norms,[1 3 2])),3))

这假定pPoint和norms是大小[nlay 3],lStart而且lEnd是大小[nray 3].结果是大小[nlay nray],每个对应于(层,光线)对.

这为您的示例提供了正确的结果:

dists =

      -Inf      -Inf
       Inf       Inf
      -Inf      -Inf
       Inf       Inf
    0.2500    0.2500
    0.7500    0.7500

这是另一种介绍fast matrix-multiplication分母部分计算的方法 -

p1 = sum(bsxfun(@times,bsxfun(@minus,pPoint,permute(lStart,[3 2 1])),norms),2)
p2 = norms*(lEnd - lStart).'
dists = squeeze(p1)./p2

由于lStart被列为重型数据集,因此最好保持原样并置换周围的事物.因此,获得的另一种方法squeeze(p1)是 -

squeeze(sum(bsxfun(@times,bsxfun(@minus,permute(pPoint,[3 2 1]),lStart),permute(norms,[3 2 1])),2)).'