将点矢量化为线/段/光线距离函数（matlab）

Question

我正在使用here中的函数来计算多个光线到空间中多个点的距离。我现在的方式是通过形成两个for循环，你可以在下面的代码中注意到，因为distancePoint2Line()函数没有矢量化。然而，当我必须通过大约2000-3000光线的光线并在空间中找到6000-7000点的距离时，这需要花费很多时间。因此，我想尽可能地矢量化distancePoint2Line()以优化处理时间。

现有distancePoint2Line()函数的输入是起点，线的终点，线段或光线，空间中的独立点和线的类型，简要给出如下：

function [d, C, t0] = distancePoint2Line(A, B, P, varargin)
% - lineType definition
if (nargin < 4)
    lineType = 'line';
lineType = 'line';
else
    lineType = varargin{1};
end
%% Algorithm

% Direction vector 
% M = B - A;
M = B;

% Running parameter t0 defines the intersection point of line through A and B
% and the perpendicular through P
t0  = dot(M, P - A) / dot(M, M);

% Intersection point of the perpendicular and line through A and B
intersectPnt = A + t0 * M;

switch lower(lineType)
    case 'line'
        % Line: intersection point is always closest.
        C   = intersectPnt;

    case 'segment'
        % Line segment
        if t0 < 0
            % Start point is closest.
            C   = A;
        elseif t0 > 1
            % End point is closest.
            C   = B;
        else
            % Intersection point is closest.
            C   = intersectPnt;
        end     

    case 'ray'
        % Ray
        if t0 < 0
            % Start point is closest.
            C   = A;
        else
            % Intersection point is closest.
            C   = intersectPnt;
        end

end

% Distance between independent point and closest point
d   = norm(P-C);

end

所以如果我认为我有下图：

绿点是我的起源，蓝色光线是从原点投射到空间的光线，红点是空间中的各个点，我想测量它们与每条光线的距离，就像我到目前为止所做的那样。下列：

origin = [10 10 0];
rays = [ 0.0420891179448099 -0.246431544552697  0.968245836551854;
         0.310536647204174   0.682690992135673  0.661437827766148;
        -0.746496065726383  -0.0724128707829756 0.661437827766148;
         0.435959418522208  -0.610278121352698 0.661437827766148];
mypts = [4579 4246 1962;
         3961 4670 3277;
         4796 3393  856;
         3277 3787 3529;
          179 3713  160];

for j=1:size(rays,1)
    for i=1:size(mypts,1)
        [D(i,j), C(i,j,:), t0(i,j)] = distancePoint2Line(origin, rays(j,:), mypts(j,:), 'ray');
    end
end

但是，正如我所说，当光线数量和个别点数增加时，上述实施方案太慢了。

因此，我尝试修改distancePoint2Line()函数，以便我对处理进行矢量化并将原点，所有结束点和所有单个点一次作为参数传递。

[d, c, tt0] = distancePoint2LineVec(origin, rays, mypts, 'ray');

以下是我的修改，直到获得交叉点：

function [d, C, t0] = distancePoint2LineVec(A, B, P, varargin)
% - lineType definition
if (nargin < 4)
    lineType = 'line';
lineType = 'line';
else
    lineType = varargin{1};
end
%% Algorithm

% Direction vector 
% M = B - A;
% M = bsxfun(@minus, B, A); 
M = B;

% Running parameter t0 defines the intersection point of line through A and B
% and the perpendicular through P
% t0  = dot(M, P - A) / dot(M, M);
K = bsxfun(@minus, P, A);
V = M*K.';
U = M*M.';
t0 = bsxfun(@rdivide, V, diag(U));
t0 = t0(:);

% Intersection point of the perpendicular and line through A and B
% intersectPnt = A + t0 * M;
mul = M.*repmat(t0, [1 size(M,1)]); <---- Getting matrix dimension error here
intersectPnt = bsxfun(@plus,A,mul);

switch lower(lineType)
    case 'line'
        % Line: intersection point is always closest.
        C   = intersectPnt;

    case 'segment'
        % Line segment
        if t0 < 0
            % Start point is closest.
            C   = A;
        elseif t0 > 1
            % End point is closest.
            C   = B;
        else
            % Intersection point is closest.
            C   = intersectPnt;
        end     

    case 'ray'
        % Ray
        if t0 < 0
            % Start point is closest.
            C   = A;
        else
            % Intersection point is closest.
            C   = intersectPnt;
        end

end

% Distance between independent point and closest point
d   = norm(P-C);

end

但重点是我在那里得到一个错误，因为尺寸与乘法不匹配。我搜索了一下，似乎没有直接的方法来做到这一点。任何想法如何有效地解决它？

谢谢。

Answer 1

好的，这是我结束的第一个解决方案，看起来非常快：

function [d, C, t0] = distancePoint2LineVec(A, B, P, varargin)
% - lineType definition
if (nargin < 4)
    lineType = 'line';
else
    lineType = varargin{1};
end
%% Algorithm

% find if the origin point is withing the points and initialize it to Inf
[~,indx]=ismember(A,P,'rows');

if indx ~= 0
    P(indx,:) = Inf;
end

% Direction vector 
% M = B - A;
% M = bsxfun(@minus, B, A); 
M = B;

% Running parameter t0 defines the intersection point of line through A and B
% and the perpendicular through P
% t0  = dot(M, P - A) / dot(M, M);
K = bsxfun(@minus, P, A);
V = M*K.';
U = M*M.';
t0 = bsxfun(@rdivide, V, diag(U));
t0 = t0(:);

% Intersection point of the perpendicular and line through A and B
% intersectPnt = A + t0 * M;
mul = bsxfun(@times, repmat(M, [size(P,1) 1]), t0);
intersectPnt = bsxfun(@plus, repmat(A, [size(t0,1) 1]), mul);

C = intersectPnt;

idx = find(t0<0);
idx2 = find(t0>1);
% C(idx,:) = repmat(A,[size(idx,1) 1]);

switch lower(lineType)
    case 'line'
        % Line: intersection point is always closest.
        C   = intersectPnt;

    case 'segment'
        % Line segment
%         if t0 < 0
            % Start point is closest.
%             C   = A;
            C(idx,:) = repmat(A,[size(idx,1) 1]);
%         elseif t0 > 1
            % End point is closest.
%             C   = B;
            C(idx2,:) = repmat(B,[size(idx2,1) 1]);
%         else
            % Intersection point is closest.
%             C   = intersectPnt;
%         end     

    case 'ray'
        % Ray
%         if t0 < 0
            % Start point is closest.
%             C   = A;
%             C   = Inf;
%         else
            % Intersection point is closest.
%             C   = intersectPnt;
%         end
        C(idx,:) = repmat([Inf Inf Inf],[size(idx,1) 1]);

end

% Distance between independent point and closest point
% d   = norm(P-C);
d = pdist2(P,C);
d = d(sub2ind(size(d), repmat((1:size(d,1)).', 1, size(B,1)), (1:size(B,1))+((1:size(B,1):size(d,2))-1).'));


C = permute(reshape(C, size(B,1), size(P,1), []), [2 1 3]);
t0 = reshape(t0, size(B,1), size(P,1))';

end

如果对更多改进有任何意见或建议，我将很高兴听到。再次感谢。