我需要一个基本函数来求点到线段的最短距离。你可以随意用任何你想要的语言来编写解决方案;我可以把它翻译成什么我正在使用(Javascript)。
编辑:线段由两个端点定义。线段AB由两点A (x1,y1)和B (x2,y2)定义。我要求的是这条线段到点C (x3,y3)的距离。我的几何技能生疏了,所以我看到的例子让我很困惑,我很遗憾地承认。
我需要一个基本函数来求点到线段的最短距离。你可以随意用任何你想要的语言来编写解决方案;我可以把它翻译成什么我正在使用(Javascript)。
编辑:线段由两个端点定义。线段AB由两点A (x1,y1)和B (x2,y2)定义。我要求的是这条线段到点C (x3,y3)的距离。我的几何技能生疏了,所以我看到的例子让我很困惑,我很遗憾地承认。
当前回答
伊莱,你选定的代码是错误的。在线段所在直线附近但远离线段一端的点将被错误地判断为接近线段。更新:上面提到的错误答案已不再被接受。
下面是一些正确的c++代码。它假设一个2d向量类vec2 {float x,y;},本质上,带有加法、subract、缩放等运算符,以及一个距离和点积函数(即x1 x2 + y1 y2)。
float minimum_distance(vec2 v, vec2 w, vec2 p) {
// Return minimum distance between line segment vw and point p
const float l2 = length_squared(v, w); // i.e. |w-v|^2 - avoid a sqrt
if (l2 == 0.0) return distance(p, v); // v == w case
// Consider the line extending the segment, parameterized as v + t (w - v).
// We find projection of point p onto the line.
// It falls where t = [(p-v) . (w-v)] / |w-v|^2
// We clamp t from [0,1] to handle points outside the segment vw.
const float t = max(0, min(1, dot(p - v, w - v) / l2));
const vec2 projection = v + t * (w - v); // Projection falls on the segment
return distance(p, projection);
}
编辑:我需要一个Javascript实现,所以在这里,没有依赖关系(或注释,但它是一个直接的端口以上)。点被表示为具有x和y属性的对象。
function sqr(x) { return x * x }
function dist2(v, w) { return sqr(v.x - w.x) + sqr(v.y - w.y) }
function distToSegmentSquared(p, v, w) {
var l2 = dist2(v, w);
if (l2 == 0) return dist2(p, v);
var t = ((p.x - v.x) * (w.x - v.x) + (p.y - v.y) * (w.y - v.y)) / l2;
t = Math.max(0, Math.min(1, t));
return dist2(p, { x: v.x + t * (w.x - v.x),
y: v.y + t * (w.y - v.y) });
}
function distToSegment(p, v, w) { return Math.sqrt(distToSegmentSquared(p, v, w)); }
编辑2:我需要一个Java版本,但更重要的是,我需要3d版本而不是2d版本。
float dist_to_segment_squared(float px, float py, float pz, float lx1, float ly1, float lz1, float lx2, float ly2, float lz2) {
float line_dist = dist_sq(lx1, ly1, lz1, lx2, ly2, lz2);
if (line_dist == 0) return dist_sq(px, py, pz, lx1, ly1, lz1);
float t = ((px - lx1) * (lx2 - lx1) + (py - ly1) * (ly2 - ly1) + (pz - lz1) * (lz2 - lz1)) / line_dist;
t = constrain(t, 0, 1);
return dist_sq(px, py, pz, lx1 + t * (lx2 - lx1), ly1 + t * (ly2 - ly1), lz1 + t * (lz2 - lz1));
}
这里,在函数参数中,<px,py,pz>是问题点,线段有端点<lx1,ly1,lz1>和<lx2,ly2,lz2>。函数dist_sq(假定存在)求两点之间距离的平方。
其他回答
in R
#distance beetween segment ab and point c in 2D space
getDistance_ort_2 <- function(a, b, c){
#go to complex numbers
A<-c(a[1]+1i*a[2],b[1]+1i*b[2])
q=c[1]+1i*c[2]
#function to get coefficients of line (ab)
getAlphaBeta <- function(A)
{ a<-Re(A[2])-Re(A[1])
b<-Im(A[2])-Im(A[1])
ab<-as.numeric()
ab[1] <- -Re(A[1])*b/a+Im(A[1])
ab[2] <-b/a
if(Im(A[1])==Im(A[2])) ab<- c(Im(A[1]),0)
if(Re(A[1])==Re(A[2])) ab <- NA
return(ab)
}
#function to get coefficients of line ortogonal to line (ab) which goes through point q
getAlphaBeta_ort<-function(A,q)
{ ab <- getAlphaBeta(A)
coef<-c(Re(q)/ab[2]+Im(q),-1/ab[2])
if(Re(A[1])==Re(A[2])) coef<-c(Im(q),0)
return(coef)
}
#function to get coordinates of interception point
#between line (ab) and its ortogonal which goes through point q
getIntersection_ort <- function(A, q){
A.ab <- getAlphaBeta(A)
q.ab <- getAlphaBeta_ort(A,q)
if (!is.na(A.ab[1])&A.ab[2]==0) {
x<-Re(q)
y<-Im(A[1])}
if (is.na(A.ab[1])) {
x<-Re(A[1])
y<-Im(q)
}
if (!is.na(A.ab[1])&A.ab[2]!=0) {
x <- (q.ab[1] - A.ab[1])/(A.ab[2] - q.ab[2])
y <- q.ab[1] + q.ab[2]*x}
xy <- x + 1i*y
return(xy)
}
intersect<-getIntersection_ort(A,q)
if ((Mod(A[1]-intersect)+Mod(A[2]-intersect))>Mod(A[1]-A[2])) {dist<-min(Mod(A[1]-q),Mod(A[2]-q))
} else dist<-Mod(q-intersect)
return(dist)
}
看起来几乎每个人都在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%。
下面是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;
}
特征c++版本的3D线段和点
// Return minimum distance between line segment: head--->tail and point
double MinimumDistance(Eigen::Vector3d head, Eigen::Vector3d tail,Eigen::Vector3d point)
{
double l2 = std::pow((head - tail).norm(),2);
if(l2 ==0.0) return (head - point).norm();// head == tail case
// Consider the line extending the segment, parameterized as head + t (tail - point).
// We find projection of point onto the line.
// It falls where t = [(point-head) . (tail-head)] / |tail-head|^2
// We clamp t from [0,1] to handle points outside the segment head--->tail.
double t = max(0,min(1,(point-head).dot(tail-head)/l2));
Eigen::Vector3d projection = head + t*(tail-head);
return (point - projection).norm();
}
I'm assuming you want to find the shortest distance between the point and a line segment; to do this, you need to find the line (lineA) which is perpendicular to your line segment (lineB) which goes through your point, determine the intersection between that line (lineA) and your line which goes through your line segment (lineB); if that point is between the two points of your line segment, then the distance is the distance between your point and the point you just found which is the intersection of lineA and lineB; if the point is not between the two points of your line segment, you need to get the distance between your point and the closer of two ends of the line segment; this can be done easily by taking the square distance (to avoid a square root) between the point and the two points of the line segment; whichever is closer, take the square root of that one.