我需要一个基本函数来求点到线段的最短距离。你可以随意用任何你想要的语言来编写解决方案;我可以把它翻译成什么我正在使用(Javascript)。
编辑:线段由两个端点定义。线段AB由两点A (x1,y1)和B (x2,y2)定义。我要求的是这条线段到点C (x3,y3)的距离。我的几何技能生疏了,所以我看到的例子让我很困惑,我很遗憾地承认。
我需要一个基本函数来求点到线段的最短距离。你可以随意用任何你想要的语言来编写解决方案;我可以把它翻译成什么我正在使用(Javascript)。
编辑:线段由两个端点定义。线段AB由两点A (x1,y1)和B (x2,y2)定义。我要求的是这条线段到点C (x3,y3)的距离。我的几何技能生疏了,所以我看到的例子让我很困惑,我很遗憾地承认。
当前回答
这是一个自成体系的Delphi / Pascal版本的函数,基于上面约书亚的答案。使用TPoint用于VCL屏幕图形,但应该易于根据需要进行调整。
function DistancePtToSegment( pt, pt1, pt2: TPoint): double;
var
a, b, c, d: double;
len_sq: double;
param: double;
xx, yy: double;
dx, dy: double;
begin
a := pt.x - pt1.x;
b := pt.y - pt1.y;
c := pt2.x - pt1.x;
d := pt2.y - pt1.y;
len_sq := (c * c) + (d * d);
param := -1;
if (len_sq <> 0) then
begin
param := ((a * c) + (b * d)) / len_sq;
end;
if param < 0 then
begin
xx := pt1.x;
yy := pt1.y;
end
else if param > 1 then
begin
xx := pt2.x;
yy := pt2.y;
end
else begin
xx := pt1.x + param * c;
yy := pt1.y + param * d;
end;
dx := pt.x - xx;
dy := pt.y - yy;
result := sqrt( (dx * dx) + (dy * dy))
end;
其他回答
这是我最后写的代码。这段代码假设一个点以{x:5, y:7}的形式定义。注意,这不是绝对最有效的方法,但它是我能想到的最简单、最容易理解的代码。
// a, b, and c in the code below are all points
function distance(a, b)
{
var dx = a.x - b.x;
var dy = a.y - b.y;
return Math.sqrt(dx*dx + dy*dy);
}
function Segment(a, b)
{
var ab = {
x: b.x - a.x,
y: b.y - a.y
};
var length = distance(a, b);
function cross(c) {
return ab.x * (c.y-a.y) - ab.y * (c.x-a.x);
};
this.distanceFrom = function(c) {
return Math.min(distance(a,c),
distance(b,c),
Math.abs(cross(c) / length));
};
}
看起来几乎每个人都在StackOverflow上贡献了一个答案(目前为止有23个答案),所以这里是我对c#的贡献。这主要是基于M. Katz的回答,而Katz的回答又基于Grumdrig的回答。
public struct MyVector
{
private readonly double _x, _y;
// Constructor
public MyVector(double x, double y)
{
_x = x;
_y = y;
}
// Distance from this point to another point, squared
private double DistanceSquared(MyVector otherPoint)
{
double dx = otherPoint._x - this._x;
double dy = otherPoint._y - this._y;
return dx * dx + dy * dy;
}
// Find the distance from this point to a line segment (which is not the same as from this
// point to anywhere on an infinite line). Also returns the closest point.
public double DistanceToLineSegment(MyVector lineSegmentPoint1, MyVector lineSegmentPoint2,
out MyVector closestPoint)
{
return Math.Sqrt(DistanceToLineSegmentSquared(lineSegmentPoint1, lineSegmentPoint2,
out closestPoint));
}
// Same as above, but avoid using Sqrt(), saves a new nanoseconds in cases where you only want
// to compare several distances to find the smallest or largest, but don't need the distance
public double DistanceToLineSegmentSquared(MyVector lineSegmentPoint1,
MyVector lineSegmentPoint2, out MyVector closestPoint)
{
// Compute length of line segment (squared) and handle special case of coincident points
double segmentLengthSquared = lineSegmentPoint1.DistanceSquared(lineSegmentPoint2);
if (segmentLengthSquared < 1E-7f) // Arbitrary "close enough for government work" value
{
closestPoint = lineSegmentPoint1;
return this.DistanceSquared(closestPoint);
}
// Use the magic formula to compute the "projection" of this point on the infinite line
MyVector lineSegment = lineSegmentPoint2 - lineSegmentPoint1;
double t = (this - lineSegmentPoint1).DotProduct(lineSegment) / segmentLengthSquared;
// Handle the two cases where the projection is not on the line segment, and the case where
// the projection is on the segment
if (t <= 0)
closestPoint = lineSegmentPoint1;
else if (t >= 1)
closestPoint = lineSegmentPoint2;
else
closestPoint = lineSegmentPoint1 + (lineSegment * t);
return this.DistanceSquared(closestPoint);
}
public double DotProduct(MyVector otherVector)
{
return this._x * otherVector._x + this._y * otherVector._y;
}
public static MyVector operator +(MyVector leftVector, MyVector rightVector)
{
return new MyVector(leftVector._x + rightVector._x, leftVector._y + rightVector._y);
}
public static MyVector operator -(MyVector leftVector, MyVector rightVector)
{
return new MyVector(leftVector._x - rightVector._x, leftVector._y - rightVector._y);
}
public static MyVector operator *(MyVector aVector, double aScalar)
{
return new MyVector(aVector._x * aScalar, aVector._y * aScalar);
}
// Added using ReSharper due to CodeAnalysis nagging
public bool Equals(MyVector other)
{
return _x.Equals(other._x) && _y.Equals(other._y);
}
public override bool Equals(object obj)
{
if (ReferenceEquals(null, obj)) return false;
return obj is MyVector && Equals((MyVector) obj);
}
public override int GetHashCode()
{
unchecked
{
return (_x.GetHashCode()*397) ^ _y.GetHashCode();
}
}
public static bool operator ==(MyVector left, MyVector right)
{
return left.Equals(right);
}
public static bool operator !=(MyVector left, MyVector right)
{
return !left.Equals(right);
}
}
这是一个小测试程序。
public static class JustTesting
{
public static void Main()
{
Stopwatch stopwatch = new Stopwatch();
stopwatch.Start();
for (int i = 0; i < 10000000; i++)
{
TestIt(1, 0, 0, 0, 1, 1, 0.70710678118654757);
TestIt(5, 4, 0, 0, 20, 10, 1.3416407864998738);
TestIt(30, 15, 0, 0, 20, 10, 11.180339887498949);
TestIt(-30, 15, 0, 0, 20, 10, 33.541019662496844);
TestIt(5, 1, 0, 0, 10, 0, 1.0);
TestIt(1, 5, 0, 0, 0, 10, 1.0);
}
stopwatch.Stop();
TimeSpan timeSpan = stopwatch.Elapsed;
}
private static void TestIt(float aPointX, float aPointY,
float lineSegmentPoint1X, float lineSegmentPoint1Y,
float lineSegmentPoint2X, float lineSegmentPoint2Y,
double expectedAnswer)
{
// Katz
double d1 = DistanceFromPointToLineSegment(new MyVector(aPointX, aPointY),
new MyVector(lineSegmentPoint1X, lineSegmentPoint1Y),
new MyVector(lineSegmentPoint2X, lineSegmentPoint2Y));
Debug.Assert(d1 == expectedAnswer);
/*
// Katz using squared distance
double d2 = DistanceFromPointToLineSegmentSquared(new MyVector(aPointX, aPointY),
new MyVector(lineSegmentPoint1X, lineSegmentPoint1Y),
new MyVector(lineSegmentPoint2X, lineSegmentPoint2Y));
Debug.Assert(Math.Abs(d2 - expectedAnswer * expectedAnswer) < 1E-7f);
*/
/*
// Matti (optimized)
double d3 = FloatVector.DistanceToLineSegment(new PointF(aPointX, aPointY),
new PointF(lineSegmentPoint1X, lineSegmentPoint1Y),
new PointF(lineSegmentPoint2X, lineSegmentPoint2Y));
Debug.Assert(Math.Abs(d3 - expectedAnswer) < 1E-7f);
*/
}
private static double DistanceFromPointToLineSegment(MyVector aPoint,
MyVector lineSegmentPoint1, MyVector lineSegmentPoint2)
{
MyVector closestPoint; // Not used
return aPoint.DistanceToLineSegment(lineSegmentPoint1, lineSegmentPoint2,
out closestPoint);
}
private static double DistanceFromPointToLineSegmentSquared(MyVector aPoint,
MyVector lineSegmentPoint1, MyVector lineSegmentPoint2)
{
MyVector closestPoint; // Not used
return aPoint.DistanceToLineSegmentSquared(lineSegmentPoint1, lineSegmentPoint2,
out closestPoint);
}
}
如您所见,我试图衡量使用避免Sqrt()方法的版本与使用普通版本之间的差异。我的测试表明你可能可以节省2.5%,但我甚至不确定——各种测试运行中的变化是相同的数量级。我还试着测量了Matti发布的版本(加上一个明显的优化),该版本似乎比基于Katz/Grumdrig代码的版本慢了大约4%。
编辑:顺便说一句,我还尝试过测量一种方法,该方法使用叉乘(和平方根())来查找到无限直线(不是线段)的距离,它大约快32%。
我制作了一个交互式Desmos图来演示如何实现这一点:
https://www.desmos.com/calculator/kswrm8ddum
红点是A点,绿点是B点,C点是蓝色点。 您可以拖动图形中的点来查看值的变化。 左边的值“s”是线段的参数(即s = 0表示点A, s = 1表示点B)。 值“d”是第三点到经过A和B的直线的距离。
编辑:
有趣的小见解:坐标(s, d)是坐标系中第三点C的坐标,AB是单位x轴,单位y轴垂直于AB。
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;
}
使用arctangents的一行解决方案:
思路是将A移动到(0,0),并顺时针旋转三角形,使C位于X轴上, 当这种情况发生时,By就是距离。
a角= Atan(Cy - Ay, Cx - Ax); b角= Atan(By - Ay, Bx - Ax); AB长度=平方根((Bx - Ax)²+ (By - Ay)²) By = Sin (bAngle - aAngle) * ABLength
C#
public double Distance(Point a, Point b, Point c)
{
// normalize points
Point cn = new Point(c.X - a.X, c.Y - a.Y);
Point bn = new Point(b.X - a.X, b.Y - a.Y);
double angle = Math.Atan2(bn.Y, bn.X) - Math.Atan2(cn.Y, cn.X);
double abLength = Math.Sqrt(bn.X*bn.X + bn.Y*bn.Y);
return Math.Sin(angle)*abLength;
}
一行c#(要转换为SQL)
double distance = Math.Sin(Math.Atan2(b.Y - a.Y, b.X - a.X) - Math.Atan2(c.Y - a.Y, c.X - a.X)) * Math.Sqrt((b.X - a.X) * (b.X - a.X) + (b.Y - a.Y) * (b.Y - a.Y))