计算特殊情况下的点到线段的距离 [英] Calculate distance point to line segment with special cases
问题描述
我需要实现以下逻辑: 给定一组2d采样点(作为x-y坐标对)和一组线段(也为x-y坐标对). 如何计算(矢量化)点 pi 与线 Li 的距离?
I need to implement the following logic: Given a set of 2d sample points (given as x-y-coordinate pairs) and a set of line segments (also x-y-coordinate pairs). Edit 1: How to calculate (vectorized) the distance of the points pi to the lines Li?
这些点大约在直线附近,我想获得从每个样本点到最近的线段的距离.可能是点,有些关"(请参见第一张图片中的p6),这些点可以通过以下算法找到:
The points are approximately near the lines and I want to get the distance from each of the sample points to the nearest line segment. The may be points, which are a bit "off" (see p6 in the first picture) and these could be found by the following algorithm:
- 情况:到该线段的采样点投影为左外".在这种情况下,我需要从p1到x1的欧几里德距离.
- 情况:到该线段的采样点投影在该线段的内部".在这种情况下,我需要从p2到x1到x2的线的距离.
- 情况:到该线段的采样点投影为右外".在这种情况下,我需要从p3到x2的欧几里德距离.
有一个矢量化的解决方案(感谢用户 Andy )使用全局"投影,始终为每个点线段对假设情况2.但是,这将返回p1 ... p3的距离[1 1 1]
,其中所需距离将为[1.4142 1 1.4142]
.可以修改这些代码以满足这些需求吗?
There is a vectorized solution (thanks to User Andy) which uses the projection "globally", always assuming case 2 for each point-linesegment-pair. However, this returns the distance [1 1 1]
for p1 ... p3 , where the desired distance would be [1.4142 1 1.4142]
. Can this code be modified to these needs?
ptsx = [1 3 5];
ptsy= [1 1 1];
linesx = [2 4];
linesy = [0 0];
points = [ptsx;ptsy];
lines = [linesx;linesy];
% 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 (points, 2));
% vector from points (in third dimension) to start of lines (second dimension)
r = bsxfun (@minus, permute (points, [1 3 2]), l);
d = abs (dot (v, r));
dist = squeeze (min (d, [], 2))
从数学上讲,可以通过查看vec(pi-x1)
在vec(x2-x1)
上的投影长度来区分大小写.如果该长度因数<. 0,则该点位于线段的左侧",如果在0到1之间,则可以进行垂直投影;如果大于1,则该点位于线段的右侧".
Mathematically speaking, the cases can be separated by looking at the length of the projection of vec(pi-x1)
onto vec(x2-x1)
. If this length factor is < 0, then the point is "left" of the line segment, if it is between 0 and 1, perpendicular projection is possible, if it is > 1, then the point is "right" of the line segment.
我将添加一个伪代码来说明如何使用双for循环解决此问题,但是由于我有大约6000个样本和10000行,因此循环解决方案对我来说不是一个选择.
Edit 1: I will add a pseudocode to explain how this could be solved with a double for-loop, but since I have around 6000 samples and 10000 lines, the loop solution is not an option for me.
for each sample point pi
for each linesegment Li
a = vector from start of Li to end of Li
b = vector from pi to start of Li
relLength = dot(a,b)/norm(a)^2
if relLength < 0: distance = euclidean distance from start of Li to pi
if relLength > 1: distance = euclidean distance from end of Li to pi
else: distance = perpendicular distance from pi to Li
endfor
endfor
编辑2/2017-09-07:我设法向量化了该算法的第一部分. relLength 现在包含每个pi-startOfLi
在每个线段上的投影的相对长度.
Edit 2 / 2017-09-07: I managed to vectorize the first part of this algorithm. relLength now contains the relative length of the projection of each pi-startOfLi
onto each line segment.
ptsx = [0.5 2 3 5.5 8 11];
ptsy= [1 2 -1.5 0.5 4 5];
linesx = [0 2 4 6 10 10 0 0];
linesy = [0 0 0 0 0 4 4 0];
points = [ptsx;ptsy];
lines = [linesx;linesy];
% Start of all lines
L1 = [linesx(1:end-1); linesy(1:end-1)];
% Vector of each line segment
a = [diff(linesx); diff(linesy)];
b = bsxfun(@minus, permute(points, [1 3 2]), L1);
aRep = repmat(a, 1, 1, length(b(1,1,:)));
relLength = dot(aRep,b)./norm(a, 'cols').^2
推荐答案
在GNU Octave中:
In GNU Octave:
points = [1 4.3 3.7 2.9;0.8 0.8 2.1 -0.5];
lines = [0 2 4 3.6;0 -1 1 1.75];
% plot them
hold off
plot (points(1,:), points(2,:), 'or')
hold on
plot (lines(1,:), lines(2,:), '-xb')
text (points(1,:), points(2,:),...
arrayfun (@(x) sprintf(' p%i',x),...
1:columns(points),'UniformOutput', false))
axis ('equal')
grid on
zoom (0.9);
% some intermediate vars
s_lines = lines (:,1:end-1); % start of lines
d_lines = diff(lines, 1, 2); % vectors between line points
l_lines = hypot (d_lines(1,:),
d_lines(2,:)); % length of lines
现在做真实的"工作:
% vectors perpendicular on lines
v = [0 1;-1 0] * d_lines;
vn = v ./ norm (v, 2, 'cols'); %make unit vector
% extend to number of points
vn = repmat (vn, 1, 1, columns (points));
points_3 = permute (points, [1 3 2]);
% vectors from points (in third dimension) to start of lines (second dimension)
d = dot (vn, points_3 - s_lines);
% check if the projection is on line
tmp = dot (repmat (d_lines, 1, 1, columns (points)),...
points_3 - s_lines)./l_lines.^2;
point_hits_line = tmp > 0 & tmp < 1;
% set othogonal distance to Inf if there is no hit
d(~ point_hits_line) = Inf;
dist = squeeze (min (abs(d), [], 2));
% calculate the euclidean distance from points to line start/stop
% if the projection doesn't hit the line
nh = isinf (dist);
tmp = points_3(:,:,nh) - lines;
tmp = hypot(tmp(1,:,:),tmp(2,:,:));
tmp = min (tmp, [], 2);
% store the result back
dist (nh) = tmp
将结果绘制成点周围的黄色圆圈
plot the results as yellow circles around the points
% use for loops because this hasn't to be fast
t = linspace (0, 2*pi, 40);
for k=1:numel(dist)
plot (points (1, k) + cos (t) * dist(k),
points (2, k) + sin (t) * dist(k),
'-y')
end
这篇关于计算特殊情况下的点到线段的距离的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!