我需要一个基本函数来求点到线段的最短距离。你可以随意用任何你想要的语言来编写解决方案;我可以把它翻译成什么我正在使用(Javascript)。
编辑:线段由两个端点定义。线段AB由两点A (x1,y1)和B (x2,y2)定义。我要求的是这条线段到点C (x3,y3)的距离。我的几何技能生疏了,所以我看到的例子让我很困惑,我很遗憾地承认。
我需要一个基本函数来求点到线段的最短距离。你可以随意用任何你想要的语言来编写解决方案;我可以把它翻译成什么我正在使用(Javascript)。
编辑:线段由两个端点定义。线段AB由两点A (x1,y1)和B (x2,y2)定义。我要求的是这条线段到点C (x3,y3)的距离。我的几何技能生疏了,所以我看到的例子让我很困惑,我很遗憾地承认。
当前回答
JavaScript中一个基于这个公式的更简洁的解决方案:
distToSegment: function (point, linePointA, linePointB){
var x0 = point.X;
var y0 = point.Y;
var x1 = linePointA.X;
var y1 = linePointA.Y;
var x2 = linePointB.X;
var y2 = linePointB.Y;
var Dx = (x2 - x1);
var Dy = (y2 - y1);
var numerator = Math.abs(Dy*x0 - Dx*y0 - x1*y2 + x2*y1);
var denominator = Math.sqrt(Dx*Dx + Dy*Dy);
if (denominator == 0) {
return this.dist2(point, linePointA);
}
return numerator/denominator;
}
其他回答
下面是devnullicus转换为c#的c++版本。对于我的实现,我需要知道交叉点,并找到他的解决方案。
public static bool PointSegmentDistanceSquared(PointF point, PointF lineStart, PointF lineEnd, out double distance, out PointF intersectPoint)
{
const double kMinSegmentLenSquared = 0.00000001; // adjust to suit. If you use float, you'll probably want something like 0.000001f
const double kEpsilon = 1.0E-14; // adjust to suit. If you use floats, you'll probably want something like 1E-7f
double dX = lineEnd.X - lineStart.X;
double dY = lineEnd.Y - lineStart.Y;
double dp1X = point.X - lineStart.X;
double dp1Y = point.Y - lineStart.Y;
double segLenSquared = (dX * dX) + (dY * dY);
double t = 0.0;
if (segLenSquared >= -kMinSegmentLenSquared && segLenSquared <= kMinSegmentLenSquared)
{
// segment is a point.
intersectPoint = lineStart;
t = 0.0;
distance = ((dp1X * dp1X) + (dp1Y * dp1Y));
}
else
{
// Project a line from p to the segment [p1,p2]. By considering the line
// extending the segment, parameterized as p1 + (t * (p2 - p1)),
// we find projection of point p onto the line.
// It falls where t = [(p - p1) . (p2 - p1)] / |p2 - p1|^2
t = ((dp1X * dX) + (dp1Y * dY)) / segLenSquared;
if (t < kEpsilon)
{
// intersects at or to the "left" of first segment vertex (lineStart.X, lineStart.Y). If t is approximately 0.0, then
// intersection is at p1. If t is less than that, then there is no intersection (i.e. p is not within
// the 'bounds' of the segment)
if (t > -kEpsilon)
{
// intersects at 1st segment vertex
t = 0.0;
}
// set our 'intersection' point to p1.
intersectPoint = lineStart;
// Note: If you wanted the ACTUAL intersection point of where the projected lines would intersect if
// we were doing PointLineDistanceSquared, then intersectPoint.X would be (lineStart.X + (t * dx)) and intersectPoint.Y would be (lineStart.Y + (t * dy)).
}
else if (t > (1.0 - kEpsilon))
{
// intersects at or to the "right" of second segment vertex (lineEnd.X, lineEnd.Y). If t is approximately 1.0, then
// intersection is at p2. If t is greater than that, then there is no intersection (i.e. p is not within
// the 'bounds' of the segment)
if (t < (1.0 + kEpsilon))
{
// intersects at 2nd segment vertex
t = 1.0;
}
// set our 'intersection' point to p2.
intersectPoint = lineEnd;
// Note: If you wanted the ACTUAL intersection point of where the projected lines would intersect if
// we were doing PointLineDistanceSquared, then intersectPoint.X would be (lineStart.X + (t * dx)) and intersectPoint.Y would be (lineStart.Y + (t * dy)).
}
else
{
// The projection of the point to the point on the segment that is perpendicular succeeded and the point
// is 'within' the bounds of the segment. Set the intersection point as that projected point.
intersectPoint = new PointF((float)(lineStart.X + (t * dX)), (float)(lineStart.Y + (t * dY)));
}
// return the squared distance from p to the intersection point. Note that we return the squared distance
// as an optimization because many times you just need to compare relative distances and the squared values
// works fine for that. If you want the ACTUAL distance, just take the square root of this value.
double dpqX = point.X - intersectPoint.X;
double dpqY = point.Y - intersectPoint.Y;
distance = ((dpqX * dpqX) + (dpqY * dpqY));
}
return true;
}
忍不住用python来编码:)
from math import sqrt, fabs
def pdis(a, b, c):
t = b[0]-a[0], b[1]-a[1] # Vector ab
dd = sqrt(t[0]**2+t[1]**2) # Length of ab
t = t[0]/dd, t[1]/dd # unit vector of ab
n = -t[1], t[0] # normal unit vector to ab
ac = c[0]-a[0], c[1]-a[1] # vector ac
return fabs(ac[0]*n[0]+ac[1]*n[1]) # Projection of ac to n (the minimum distance)
print pdis((1,1), (2,2), (2,0)) # Example (answer is 1.414)
fortran也是一样:)
real function pdis(a, b, c)
real, dimension(0:1), intent(in) :: a, b, c
real, dimension(0:1) :: t, n, ac
real :: dd
t = b - a ! Vector ab
dd = sqrt(t(0)**2+t(1)**2) ! Length of ab
t = t/dd ! unit vector of ab
n = (/-t(1), t(0)/) ! normal unit vector to ab
ac = c - a ! vector ac
pdis = abs(ac(0)*n(0)+ac(1)*n(1)) ! Projection of ac to n (the minimum distance)
end function pdis
program test
print *, pdis((/1.0,1.0/), (/2.0,2.0/), (/2.0,0.0/)) ! Example (answer is 1.414)
end program test
%Matlab solution by Tim from Cody
function ans=distP2S(x0,y0,x1,y1,x2,y2)
% Point is x0,y0
z=complex(x0-x1,y0-y1);
complex(x2-x1,y2-y1);
abs(z-ans*min(1,max(0,real(z/ans))));
快速实现http://paulbourke.net/geometry/pointlineplane/source.c
static func magnitude(p1: CGPoint, p2: CGPoint) -> CGFloat {
let vector = CGPoint(x: p2.x - p1.x, y: p2.y - p1.y)
return sqrt(pow(vector.x, 2) + pow(vector.y, 2))
}
/// http://paulbourke.net/geometry/pointlineplane/
/// http://paulbourke.net/geometry/pointlineplane/source.c
static func pointDistanceToLine(point: CGPoint, lineStart: CGPoint, lineEnd: CGPoint) -> CGFloat? {
let lineMag = magnitude(p1: lineEnd, p2: lineStart)
let u = (((point.x - lineStart.x) * (lineEnd.x - lineStart.x)) +
((point.y - lineStart.y) * (lineEnd.y - lineStart.y))) /
(lineMag * lineMag)
if u < 0 || u > 1 {
// closest point does not fall within the line segment
return nil
}
let intersectionX = lineStart.x + u * (lineEnd.x - lineStart.x)
let intersectionY = lineStart.y + u * (lineEnd.y - lineStart.y)
return magnitude(p1: point, p2: CGPoint(x: intersectionX, y: intersectionY))
}
c#版本
public static FP DistanceToLineSegment(FPVector3 a, FPVector3 b, FPVector3 point)
{
var d = b - a;
var s = d.SqrMagnitude;
var ds = d / s;
var lambda = FPVector3.Dot(point - a, ds);
var p = FPMath.Clamp01(lambda) * d;
return (a + p - point).Magnitude;
}