如何在 matlab/octave 中矢量化这个丑陋的嵌套循环？

Question

不幸的是，我无法加快以下函数的速度，并且阅读各种“如何在 Matlab/Octave 中进行矢量化”并没有帮助我解决这个特定主题。

这是我想要实现的目标： 我有一组二维样本点（以 x-y 坐标对给出）和一组线段（也是 x-y 坐标对）。这些点大致靠近线，我想得到每个样本点到最近线段的距离，但前提是我可以将样本垂直投影到线段上。

因此，虽然我已经设法将算法放入嵌套的 for 循环中，并且它为附加示例提供了正确的结果，但对于实际数据集（大约 4000 个样本点和 6000 个线段）来说，它的速度非常慢（正如预期的那样） ) 而且，看起来真的很丑。

谁能帮我为这段代码创建一个更复杂的版本？

编辑：此算法中使用的数学可以在这里查找： http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html

clc;
clear all;
close all;

ptsx = [0.5 3 5];
ptsy= [1 -1.5 0.5];
points = [ptsx; ptsy];

linesx = [0 2 4 6];
linesy = [0 0 0 0];
lines = [linesx; linesy];

% for each point in the sample dataset
for k=1:1:length(points(1,:))
  clear 'distvec'
  % calculate the distance to each line segment in the model dataset
  for l=1:1:length(lines(1,:))-1

    % vector of the line segment
    a = [lines(1,l+1)-lines(1,l), lines(2,l+1)-lines(2,l)];
    % vector from the start of the line segment to the sample point
    b = [points(1,k) - lines(1,l), points(2,k) - lines(2,l)];

    % check if the sample point can be projected onto the line segment
    if norm(a) ~= 0
      lba = dot(a,b)/norm(a)^2;
    else
      lba = -1;
    end

    if (lba >= 0) && (lba <= 1)
      % calculate distance from sample point to single line segment
      x1 = [lines(1,l) lines(2,l)];
      x2 = [lines(1,l+1) lines(2,l+1)];
      x0 = [points(1,k) points(2,k)];
      dist = abs(det([x2-x1; x1-x0]))/norm(x2-x1);
      distvec(end+1) = dist;
    end
  end
  dist(end+1) = min(distvec);

end

figure;
hold on;
plot(ptsx, ptsy, 'bo');
plot(linesx, linesy, 'r-o');
xlim([-1 7]);
ylim([-2 2]);

Answer 1

@Edit Variant 没有循环，但我可以管理最真实的矢量化。目前我无法测试它（它需要 r2016b 和更新版本），但根据 documentation 它应该可以工作

ptsx = [3 0.5 5];
ptsy= [-1.5 1 0.5];
%order row is [x,y]
points = [ptsx; ptsy]';

linesx = [0 2 4 6];
linesy = [0 0 0 0];
lines = kron([linesx;linesy],[1 1]);
lines = lines(:,2:size(lines,2)-1);
%each row is [x1 y1 x2 y2]
lines=reshape(lines,4,[])';

lenx = lines(:,3) - lines(:,1);
leny = lines(:,4) - lines(:,2);
%remove degenerated lines
sq_norm = [lenx.*lenx + leny.*leny]';
rem_idx = sq_norm < eps;
sq_norm(rem_idx) = 1;

%Starting from r2016b no need of dirty hacks with bsxfun. Finally!
diff_startx = points(:,1) - lines(:,1)';
diff_starty = points(:,2) - lines(:,2)';

pos = (diff_startx .* lenx + diff_starty .* leny) ./ sq_norm;
pos(pos < 0) = 0;
pos(pos > 1) = 1;
dist = hypot(pos .* lenx - diff_startx, pos .* leny - diff_starty);

---

我没有循环的变体

pts = [ 3   0.5 5;... %reordered for test purposes
       -1.5 1   0.5];
cpts = num2cell(pts,1);

lines = [0 2 4 6;...
         0 0 0 0];
%convert and split intercepts into two x-y pairs
lines = kron(lines,[1 1]);
lines = lines(:,2:size(lines,2)-1); %suggesting three intercepts
clines = mat2cell(lines, 2, repmat(2,1,size(lines,2)/2));

% [lines(1,2) - lines(1,1), lines(2,2) - lines(2,1)]
diff_lines = diff(lines,1,2); 
diff_lines = diff_lines(:,1:2:size(diff_lines,2));
diff_lines = num2cell(diff_lines,1);
% norm for each line
norms = cellfun(@(x) norm(x), diff_lines,'un',0);
idx = cellfun(@(x) x~=0, norms);

% [lines(1,2) - lines(1,1), lines(2,2) - lines(2,1)]
diff_start = cellfun(@(x,y) cellfun(@(z) {x(:,1)-z(:,1)},y),...
             cpts,num2cell(repmat(clines, numel(cpts),1),2)','un',0);
%choose corresponding to nonzero norm
lba = cellfun(@(x,y,n) cellfun(@(b) dot(x,b)/norm(n)^2, y, 'un', 0) ...
              ,diff_lines(idx),diff_start(idx),norms(idx),'un',0);
%find first matching lba for each point
idx_line = cellfun(@(x) find(cellfun(@(y) (y >= 0) & (y <= 1), x),1,'first'), lba);
%reorder intermediate results
clines = clines(idx_line);
norms = norms(idx_line);
diff_lines = diff_lines(idx_line);
% [lines(1,2) - pounts(1,1), lines(2,2) - pounts(2,1)]
diff_start = cellfun(@(x,y) y(:,1)-x(:,1),...
             cpts,clines,'un',0);
dists = cellfun(@(x,y,n) abs(det([x'; y']))/n, diff_lines, diff_start, norms)

---

编辑 TS 代码以在 MATLAB 中工作

lines = [linesx; linesy];
% Edit1: distance storage preallocation
distmat = nan(length(ptsx),length(linesx)-1);
for k=1:1:length(points(1,:))
...
dist = abs(det([x2-x1; x1-x0]))/norm(x2-x1);
% Edit 2: filling the storage
      distmat(k,l) = dist;
    end
...
% Edit 3: getting the distance
dists_vals = cellfun(@(x) min(x(~isnan(x))), num2cell(distmat,2))

Answer 2

GNU Octave 的解决方案，也应该在 MATLAB 中工作。如果不是，请报告。

ptsx = [0.5 3 5];
ptsy= [1 -1.5 0.5];
linesx = [0 2 4 6];
linesy = [0 0 0 0];

p = [ptsx;ptsy];
% start of lines
l = [linesx(1:end-1);linesy(1:end-1)];
% vector perpendicular on line
v = [diff(linesy);-diff(linesx)];
% make unit vector
v = v ./ hypot (v(1,:),v(2,:));
v = repmat (v, 1, 1, size (p, 2));
% vector from points (in third dimension) to start of lines (second dimension)
r = bsxfun (@minus, permute (p, [1 3 2]), l);
d = abs (dot (v, r));
dist = squeeze (min (d, [], 2))

给予

dist =

   1.00000
   1.50000
   0.50000