一点与线段之间的最短距离

444

我需要一个基本函数来查找点与线段之间最短的距离。您可以使用任何编程语言编写解决方案,我可以将其转换为我正在使用的语言(Javascript)。

编辑:我的线段由两个端点定义。因此,我的线段AB由两个点A(x1,y1)B(x2,y2)定义。我正在尝试查找该线段与点C(x3,y3)之间的距离。我的几何技能有些生疏,所以我看到的示例很令人困惑,对此我感到抱歉。


4
@ArthurKalliokoski:那个链接失效了,但我找到了一份副本:http://paulbourke.net/geometry/pointline/ - Gunther Struyf
我不得不自己寻找并从谷歌上偶然发现了这个 -- 如果有人正在寻找一种实现方式,并且可以使用Python,那么Shapely是一个不错的选择。你需要寻找路径的LineString类。 - TC1
12
@GuntherStruyf说的链接也失效了,但这个类似的链接有效:http://paulbourke.net/geometry/pointlineplane/。 - Michael
1
如果有人想要查找点与直线之间的距离,而不是点与线段之间的距离,请查看此链接:https://gist.github.com/rhyolight/2846020 - Nick Budden
1
上面的链接已经失效。这里提供伪代码和C++示例,详细解释和推导,就像一本教科书一样,http://geomalgorithms.com/a02-_lines.html - Eric
显示剩余3条评论
56个回答

519

Eli,您选择的代码是错误的。在线段所在的直线附近但远离线段一端的点将被错误地判断为靠近线段。 更新:不再接受上述错误答案。

以下是正确的C++代码。它假定一个二维向量类class vec2 {float x,y;},具有添加、减去、缩放等运算符,以及距离和点积函数(即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实现,所以这里有一个没有依赖项的实现(没有注释,但是它是上面内容的直接移植)。点被表示为具有xy属性的对象。

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(假设已存在)可以找到两个点之间距离的平方。

4
谢谢 @Grumdrig,你的JavaScript解决方案非常棒,节省了很多时间。我将你的解决方案移植到Objective-C,并在下面添加了它。 - awolf
14
p在一条直线上的投影是距离p最近的点。(并且,在投影处与直线垂直的直线将通过点p。)数字t表示从向量vw的线段上,投影落在哪个位置。如果t为0,则投影正好在点v上;如果t为1,则投影在点w上;例如,如果t为0.5,则投影在中间位置。如果t小于0或大于1,则投影落在线段的一端或另一端。在这种情况下,到线段的距离将是到较近端点的距离。 - Grumdrig
1
我怎样才能在3D中实现这个? - clankill3r
2
糟糕 - 没有注意到有人提供了一个3D版本。@RogiSolorzano,您需要先将纬度和经度坐标转换为三维空间中的x、y、z坐标。 - Grumdrig
1
它能工作但效率低下。代码将无论如何进行分支,以求得最小值/最大值。当最小值或最大值夹紧t时,没有必要计算公式的其余部分,你已经知道需要计算pvw之间的距离。 - Soonts
显示剩余23条评论

222

这里是JavaScript中最简单的完整代码。

x、y是您的目标点,x1、y1到x2、y2是您的线段。

更新:根据评论修正了长度为0的线段问题。

function pDistance(x, y, x1, y1, x2, y2) {

  var A = x - x1;
  var B = y - y1;
  var C = x2 - x1;
  var D = y2 - y1;

  var dot = A * C + B * D;
  var len_sq = C * C + D * D;
  var param = -1;
  if (len_sq != 0) //in case of 0 length line
      param = dot / len_sq;

  var xx, yy;

  if (param < 0) {
    xx = x1;
    yy = y1;
  }
  else if (param > 1) {
    xx = x2;
    yy = y2;
  }
  else {
    xx = x1 + param * C;
    yy = y1 + param * D;
  }

  var dx = x - xx;
  var dy = y - yy;
  return Math.sqrt(dx * dx + dy * dy);
}

用于帮助可视化解决方案的图像

更新:Kotlin版本

fun getDistance(x: Double, y: Double, x1: Double, y1: Double, x2: Double, y2: Double): Double {
    val a = x - x1
    val b = y - y1
    val c = x2 - x1
    val d = y2 - y1

    val lenSq = c * c + d * d
    val param = if (lenSq != .0) { //in case of 0 length line
        val dot = a * c + b * d
        dot / lenSq
    } else {
        -1.0
    }

    val (xx, yy) = when {
        param < 0 -> x1 to y1
        param > 1 -> x2 to y2
        else -> x1 + param * c to y1 + param * d
    }

    val dx = x - xx
    val dy = y - yy
    return hypot(dx, dy)
}

12
在我看过的所有解决这个问题的代码中,我最喜欢这个。它非常清晰易读。然而,它背后的数学有点神秘。例如,点积除以长度平方实际上代表什么? - user1815201
2
点积除以长度的平方可以给出从(x1,y1)到投影距离。这是点(x,y)最接近的线段的一部分。请注意最终的else子句,其中计算了(xx,yy)- 这是点(x,y)在线段(x1,y1)-(x2,y2)上的投影。 - Logan Pickup
5
代码中对长度为0的线段的检查太靠下了。在到达安全检查之前,'len_sq' 就已经为零,代码将会除以0。请将该检查提前到代码的位置。 - HostedMetrics.com
36
米。返回的结果以米为单位。 - Joshua
2
@ViniBiso 如果你没有从其他评论中明显地看出来,那些开玩笑回答“米”的人,这个距离会以你输入的单位返回。所以如果你输入米,它就会返回米。如果你输入英里,它就会返回英里。如果你输入shklobmorphs,它就会返回shklobmorphs。几乎所有的几何数学都是如此。 - Clonkex
显示剩余15条评论

83

这是一个为有限线段而非无限线实现的函数,与其他大多数函数不同(这也是我创建它的原因)。

由Paul Bourke实现的理论

Python:

def dist(x1, y1, x2, y2, x3, y3): # x3,y3 is the point
    px = x2-x1
    py = y2-y1

    norm = px*px + py*py

    u =  ((x3 - x1) * px + (y3 - y1) * py) / float(norm)

    if u > 1:
        u = 1
    elif u < 0:
        u = 0

    x = x1 + u * px
    y = y1 + u * py

    dx = x - x3
    dy = y - y3

    # Note: If the actual distance does not matter,
    # if you only want to compare what this function
    # returns to other results of this function, you
    # can just return the squared distance instead
    # (i.e. remove the sqrt) to gain a little performance

    dist = (dx*dx + dy*dy)**.5

    return dist

AS3:

public static function segmentDistToPoint(segA:Point, segB:Point, p:Point):Number
{
    var p2:Point = new Point(segB.x - segA.x, segB.y - segA.y);
    var something:Number = p2.x*p2.x + p2.y*p2.y;
    var u:Number = ((p.x - segA.x) * p2.x + (p.y - segA.y) * p2.y) / something;

    if (u > 1)
        u = 1;
    else if (u < 0)
        u = 0;

    var x:Number = segA.x + u * p2.x;
    var y:Number = segA.y + u * p2.y;

    var dx:Number = x - p.x;
    var dy:Number = y - p.y;

    var dist:Number = Math.sqrt(dx*dx + dy*dy);

    return dist;
}

Java


(翻译:Java)
private double shortestDistance(float x1,float y1,float x2,float y2,float x3,float y3)
    {
        float px=x2-x1;
        float py=y2-y1;
        float temp=(px*px)+(py*py);
        float u=((x3 - x1) * px + (y3 - y1) * py) / (temp);
        if(u>1){
            u=1;
        }
        else if(u<0){
            u=0;
        }
        float x = x1 + u * px;
        float y = y1 + u * py;

        float dx = x - x3;
        float dy = y - y3;
        double dist = Math.sqrt(dx*dx + dy*dy);
        return dist;

    }

2
抱歉,但我尝试了这个方法,它仍然给出了无限延伸的结果。不过,我发现Grumdig的答案可以解决问题。 - Frederik
1
在那种情况下,你可能使用不当或者非无穷的含义与其他意思不同。这里有一个代码示例:http://boomie.se/upload/Drawdebug.swf - quano
奇怪。好吧,在 SWF 示例中我使用的正是这段代码(复制粘贴过来的),如果你想要的是那个 SWF 的行为,那么这就是代码。 - quano
35
变量名称的选择远非良好(p2,something,u...)。 - miguelSantirso
3
我已经尝试了该函数的Python版本,并发现如果参数是整数,则会显示不正确的结果。distAnother(0, 0, 4, 0, 2, 2) 显示为 2.8284271247461903(不正确)。 distAnother(0., 0., 4., 0., 2., 2.)显示为 2.0(正确)。请注意此问题,我认为代码可以改进,在某个位置进行浮点数转换。 - Vladimir Obrizan
显示剩余8条评论

23

在我的问题讨论串 如何在C、C# / .NET 2.0或Java中计算点和线段之间的最短2D距离?中,有人要求我在找到C#的答案后将其放在这里:因此,这是修改自http://www.topcoder.com/tc?d1=tutorials&d2=geometry1&module=Static的答案:

//Compute the dot product AB . BC
private double DotProduct(double[] pointA, double[] pointB, double[] pointC)
{
    double[] AB = new double[2];
    double[] BC = new double[2];
    AB[0] = pointB[0] - pointA[0];
    AB[1] = pointB[1] - pointA[1];
    BC[0] = pointC[0] - pointB[0];
    BC[1] = pointC[1] - pointB[1];
    double dot = AB[0] * BC[0] + AB[1] * BC[1];

    return dot;
}

//Compute the cross product AB x AC
private double CrossProduct(double[] pointA, double[] pointB, double[] pointC)
{
    double[] AB = new double[2];
    double[] AC = new double[2];
    AB[0] = pointB[0] - pointA[0];
    AB[1] = pointB[1] - pointA[1];
    AC[0] = pointC[0] - pointA[0];
    AC[1] = pointC[1] - pointA[1];
    double cross = AB[0] * AC[1] - AB[1] * AC[0];

    return cross;
}

//Compute the distance from A to B
double Distance(double[] pointA, double[] pointB)
{
    double d1 = pointA[0] - pointB[0];
    double d2 = pointA[1] - pointB[1];

    return Math.Sqrt(d1 * d1 + d2 * d2);
}

//Compute the distance from AB to C
//if isSegment is true, AB is a segment, not a line.
double LineToPointDistance2D(double[] pointA, double[] pointB, double[] pointC, 
    bool isSegment)
{
    double dist = CrossProduct(pointA, pointB, pointC) / Distance(pointA, pointB);
    if (isSegment)
    {
        double dot1 = DotProduct(pointA, pointB, pointC);
        if (dot1 > 0) 
            return Distance(pointB, pointC);

        double dot2 = DotProduct(pointB, pointA, pointC);
        if (dot2 > 0) 
            return Distance(pointA, pointC);
    }
    return Math.Abs(dist);
} 

我在@SO上不是来回答问题而是提问,所以我希望我不会因为批判而得到无数的踩,但是我只想(并且受到鼓励)分享别人的想法,因为这个线程中的解决方案要么使用了一些奇怪的语言(Fortran、Mathematica),要么被某些人标记为有错误。唯一对我有用的一个解决方案(来自Grumdrig)是用C++编写的,没有人标记它有错误。但是它缺少被调用的方法(例如点等)。


1
谢谢您发布这篇文章。但是最后一种方法中似乎存在一个明显的优化:在确定需要计算距离之前,不要计算距离。 - RenniePet
2
DotProduct的注释说它计算AB.AC,但实际上它计算的是AB.BC。 - J23
根据定义,向量叉积返回一个向量,但这里返回一个标量。 - Steak Overflow

21

在Mathematica中

它使用线段的参数化描述,并将点投影到由该线段定义的直线上。当参数在线段中从0到1变化时,如果投影超出了这个范围,我们计算到相应端点的距离,而不是直角线。

Clear["Global`*"];
 distance[{start_, end_}, pt_] := 
   Module[{param},
   param = ((pt - start).(end - start))/Norm[end - start]^2; (*parameter. the "."
                                                       here means vector product*)

   Which[
    param < 0, EuclideanDistance[start, pt],                 (*If outside bounds*)
    param > 1, EuclideanDistance[end, pt],
    True, EuclideanDistance[pt, start + param (end - start)] (*Normal distance*)
    ]
   ];  

绘图结果:

Plot3D[distance[{{0, 0}, {1, 0}}, {xp, yp}], {xp, -1, 2}, {yp, -1, 2}]

alt text

绘制距离小于截止距离的点:

alt text

等高线图:

enter image description here


21

在 F# 中,从点 c 到线段 ab 之间的距离为:

let pointToLineSegmentDistance (a: Vector, b: Vector) (c: Vector) =
  let d = b - a
  let s = d.Length
  let lambda = (c - a) * d / s
  let p = (lambda |> max 0.0 |> min s) * d / s
  (a + p - c).Length

向量d在线段上从点a指向b。将d/sc-a的点积给出无限直线和点c之间最接近点的参数。使用minmax函数将该参数夹紧到0..s范围内,以使该点位于ab之间。最后,a+p-c长度为c到线段上最近点的距离。

示例用法:

pointToLineSegmentDistance (Vector(0.0, 0.0), Vector(1.0, 0.0)) (Vector(-1.0, 1.0))

1
我认为最后一行是不正确的,应该写成:(a + p - c).Length - Blair Holloway
这仍然不能完全解决问题。修正该函数的一种方法是重新定义 lambdap,如下所示: let lambda = (c - a) * d / (s * s)let p = a + (lambda |> max 0.0 |> min 1.0) * d。 之后,该函数将返回正确的距离,例如当 a = (0,1)b = (1,0)c = (1,1) 的情况。 - mikkoma
回复代码块很糟糕 :/ - Xelnath

21

如果有人感兴趣,这里是将Joshua的JavaScript代码简单转换为Objective-C的方法:

- (double)distanceToPoint:(CGPoint)p fromLineSegmentBetween:(CGPoint)l1 and:(CGPoint)l2
{
    double A = p.x - l1.x;
    double B = p.y - l1.y;
    double C = l2.x - l1.x;
    double D = l2.y - l1.y;

    double dot = A * C + B * D;
    double len_sq = C * C + D * D;
    double param = dot / len_sq;

    double xx, yy;

    if (param < 0 || (l1.x == l2.x && l1.y == l2.y)) {
        xx = l1.x;
        yy = l1.y;
    }
    else if (param > 1) {
        xx = l2.x;
        yy = l2.y;
    }
    else {
        xx = l1.x + param * C;
        yy = l1.y + param * D;
    }

    double dx = p.x - xx;
    double dy = p.y - yy;

    return sqrtf(dx * dx + dy * dy);
}

我需要这个解决方案与 MKMapPoint 一起使用,所以我将分享它,以防其他人需要。只需进行一些小修改,这将返回以米为单位的距离:

- (double)distanceToPoint:(MKMapPoint)p fromLineSegmentBetween:(MKMapPoint)l1 and:(MKMapPoint)l2
{
    double A = p.x - l1.x;
    double B = p.y - l1.y;
    double C = l2.x - l1.x;
    double D = l2.y - l1.y;

    double dot = A * C + B * D;
    double len_sq = C * C + D * D;
    double param = dot / len_sq;

    double xx, yy;

    if (param < 0 || (l1.x == l2.x && l1.y == l2.y)) {
        xx = l1.x;
        yy = l1.y;
    }
    else if (param > 1) {
        xx = l2.x;
        yy = l2.y;
    }
    else {
        xx = l1.x + param * C;
        yy = l1.y + param * D;
    }

    return MKMetersBetweenMapPoints(p, MKMapPointMake(xx, yy));
}

这对我来说似乎很有效。感谢您的转换。 - Gregir
值得注意的是,(xx, yy) 是最近点的位置。我稍微修改了你的代码,使其返回点和距离,重构了名称以描述它们的含义,并在以下链接提供了示例:http://stackoverflow.com/a/28028023/849616。 - Nat

11

嘿,我昨天刚写的。它是用Actionscript 3.0编写的,它基本上与Javascript相同,尽管你可能没有相同的Point类。

//st = start of line segment
//b = the line segment (as in: st + b = end of line segment)
//pt = point to test
//Returns distance from point to line segment.  
//Note: nearest point on the segment to the test point is right there if we ever need it
public static function linePointDist( st:Point, b:Point, pt:Point ):Number
{
    var nearestPt:Point; //closest point on seqment to pt

    var keyDot:Number = dot( b, pt.subtract( st ) ); //key dot product
    var bLenSq:Number = dot( b, b ); //Segment length squared

    if( keyDot <= 0 )  //pt is "behind" st, use st
    {
        nearestPt = st  
    }
    else if( keyDot >= bLenSq ) //pt is "past" end of segment, use end (notice we are saving twin sqrts here cuz)
    {
        nearestPt = st.add(b);
    }
    else //pt is inside segment, reuse keyDot and bLenSq to get percent of seqment to move in to find closest point
    {
        var keyDotToPctOfB:Number = keyDot/bLenSq; //REM dot product comes squared
        var partOfB:Point = new Point( b.x * keyDotToPctOfB, b.y * keyDotToPctOfB );
        nearestPt = st.add(partOfB);
    }

    var dist:Number = (pt.subtract(nearestPt)).length;

    return dist;
}

此外,在这里有一个相当详尽且易于理解的讨论:notejot.com


谢谢 - 这正是我在寻找的代码类型。我已经发布了自己的答案,因为我设法编写出适用于当前浏览器Javascript的代码,但我已将您的答案标记为已接受,因为它简单、写得好、易于理解,并且非常感激。 - Eli Courtwright
这不是缺少了点方法吗?无论如何,计算很容易:vec1.x * vec2.x + vec1.y * vec2.y - quano

11

以下是我对 @Grumdrig 上面解决方案的Objective-C端口:

CGFloat sqr(CGFloat x) { return x*x; }
CGFloat dist2(CGPoint v, CGPoint w) { return sqr(v.x - w.x) + sqr(v.y - w.y); }
CGFloat distanceToSegmentSquared(CGPoint p, CGPoint v, CGPoint w)
{
    CGFloat l2 = dist2(v, w);
    if (l2 == 0.0f) return dist2(p, v);

    CGFloat t = ((p.x - v.x) * (w.x - v.x) + (p.y - v.y) * (w.y - v.y)) / l2;
    if (t < 0.0f) return dist2(p, v);
    if (t > 1.0f) return dist2(p, w);
    return dist2(p, CGPointMake(v.x + t * (w.x - v.x), v.y + t * (w.y - v.y)));
}
CGFloat distanceToSegment(CGPoint point, CGPoint segmentPointV, CGPoint segmentPointW)
{
    return sqrtf(distanceToSegmentSquared(point, segmentPointV, segmentPointW));
}

我从这行代码中得到了“nan”的返回值。你有什么想法吗?(顺便说一句,感谢你用Obj-C编写这段代码!)return dist2(p, CGPointMake(v.x + t * (w.x - v.x), v.y + t * (w.y - v.y))) - Gregir
sqrtf()是将x平方,而不是获取它的平方根。 - Senseful
@Senseful 不太确定你的意思。sqrtf 是平方根。https://developer.apple.com/library/mac/documentation/Darwin/Reference/Manpages/man3/sqrtf.3.html - awolf
@awolf:看一下上面的第一行代码。它定义了方法 sqrtf(x) = x*x - Senseful
@Senseful 谢谢,它的名称错误而不是执行了错误的操作。 - awolf

11

使用反正切的一行解决方案:

思路是将A移动到(0,0),并顺时针旋转三角形以使C位于X轴上, 这样,By将是距离。

  1. a角度 = Atan(Cy - Ay, Cx - Ax);
  2. b角度 = Atan(By - Ay, Bx - Ax);
  3. AB长度 = Sqrt( (Bx - Ax)^2 + (By - Ay)^2 )
  4. 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))

只有在旋转后b在线段上方(b.y不为负且b.x在a.x和c.x之间)时才有效。如果b在旋转后在线段下方,则会返回负值,例如,如果线段从0,0到10,0,并且点b在x,y处,则无论y或x是什么,距离都将为y,这是到通过这两个点的(无限长的)直线的距离,而不是到线段的距离。 - MaxSilvester

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接