我需要一个基本函数来求点到线段的最短距离。你可以随意用任何你想要的语言来编写解决方案;我可以把它翻译成什么我正在使用(Javascript)。

编辑:线段由两个端点定义。线段AB由两点A (x1,y1)和B (x2,y2)定义。我要求的是这条线段到点C (x3,y3)的距离。我的几何技能生疏了,所以我看到的例子让我很困惑,我很遗憾地承认。


当前回答

Consider this modification to Grumdrig's answer above. Many times you'll find that floating point imprecision can cause problems. I'm using doubles in the version below, but you can easily change to floats. The important part is that it uses an epsilon to handle the "slop". In addition, you'll many times want to know WHERE the intersection happened, or if it happened at all. If the returned t is < 0.0 or > 1.0, no collision occurred. However, even if no collision occurred, many times you'll want to know where the closest point on the segment to P is, and thus I use qx and qy to return this location.

double PointSegmentDistanceSquared( double px, double py,
                                    double p1x, double p1y,
                                    double p2x, double p2y,
                                    double& t,
                                    double& qx, double& qy)
{
    static const double kMinSegmentLenSquared = 0.00000001;  // adjust to suit.  If you use float, you'll probably want something like 0.000001f
    static const double kEpsilon = 1.0E-14;  // adjust to suit.  If you use floats, you'll probably want something like 1E-7f
    double dx = p2x - p1x;
    double dy = p2y - p1y;
    double dp1x = px - p1x;
    double dp1y = py - p1y;
    const double segLenSquared = (dx * dx) + (dy * dy);
    if (segLenSquared >= -kMinSegmentLenSquared && segLenSquared <= kMinSegmentLenSquared)
    {
        // segment is a point.
        qx = p1x;
        qy = p1y;
        t = 0.0;
        return ((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 (p1x, p1y).  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.
            qx = p1x;
            qy = p1y;
            // Note: If you wanted the ACTUAL intersection point of where the projected lines would intersect if
            // we were doing PointLineDistanceSquared, then qx would be (p1x + (t * dx)) and qy would be (p1y + (t * dy)).
        }
        else if (t > (1.0 - kEpsilon))
        {
            // intersects at or to the "right" of second segment vertex (p2x, p2y).  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.
            qx = p2x;
            qy = p2y;
            // Note: If you wanted the ACTUAL intersection point of where the projected lines would intersect if
            // we were doing PointLineDistanceSquared, then qx would be (p1x + (t * dx)) and qy would be (p1y + (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.
            qx = p1x + (t * dx);
            qy = p1y + (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 = px - qx;
        double dpqy = py - qy;
        return ((dpqx * dpqx) + (dpqy * dpqy));
    }
}

其他回答

嘿,我昨天才写的。它在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)

这是非常快的,因为我试图避免任何数学。战俘的功能。

如你所见,在函数的最后,我得到了直线的距离。

代码来自lib http://www.draw2d.org/graphiti/jsdoc/#!/例子

/**
 * Static util function to determine is a point(px,py) on the line(x1,y1,x2,y2)
 * A simple hit test.
 * 
 * @return {boolean}
 * @static
 * @private
 * @param {Number} coronaWidth the accepted corona for the hit test
 * @param {Number} X1 x coordinate of the start point of the line
 * @param {Number} Y1 y coordinate of the start point of the line
 * @param {Number} X2 x coordinate of the end point of the line
 * @param {Number} Y2 y coordinate of the end point of the line
 * @param {Number} px x coordinate of the point to test
 * @param {Number} py y coordinate of the point to test
 **/
graphiti.shape.basic.Line.hit= function( coronaWidth, X1, Y1,  X2,  Y2, px, py)
{
  // Adjust vectors relative to X1,Y1
  // X2,Y2 becomes relative vector from X1,Y1 to end of segment
  X2 -= X1;
  Y2 -= Y1;
  // px,py becomes relative vector from X1,Y1 to test point
  px -= X1;
  py -= Y1;
  var dotprod = px * X2 + py * Y2;
  var projlenSq;
  if (dotprod <= 0.0) {
      // px,py is on the side of X1,Y1 away from X2,Y2
      // distance to segment is length of px,py vector
      // "length of its (clipped) projection" is now 0.0
      projlenSq = 0.0;
  } else {
      // switch to backwards vectors relative to X2,Y2
      // X2,Y2 are already the negative of X1,Y1=>X2,Y2
      // to get px,py to be the negative of px,py=>X2,Y2
      // the dot product of two negated vectors is the same
      // as the dot product of the two normal vectors
      px = X2 - px;
      py = Y2 - py;
      dotprod = px * X2 + py * Y2;
      if (dotprod <= 0.0) {
          // px,py is on the side of X2,Y2 away from X1,Y1
          // distance to segment is length of (backwards) px,py vector
          // "length of its (clipped) projection" is now 0.0
          projlenSq = 0.0;
      } else {
          // px,py is between X1,Y1 and X2,Y2
          // dotprod is the length of the px,py vector
          // projected on the X2,Y2=>X1,Y1 vector times the
          // length of the X2,Y2=>X1,Y1 vector
          projlenSq = dotprod * dotprod / (X2 * X2 + Y2 * Y2);
      }
  }
    // Distance to line is now the length of the relative point
    // vector minus the length of its projection onto the line
    // (which is zero if the projection falls outside the range
    //  of the line segment).
    var lenSq = px * px + py * py - projlenSq;
    if (lenSq < 0) {
        lenSq = 0;
    }
    return Math.sqrt(lenSq)<coronaWidth;
};

一个2D和3D的解决方案

考虑基底的变化,使得线段变成(0,0,0)-(d, 0,0)和点(u, v, 0)。在这个平面上,最短的距离由

    u ≤ 0 -> d(A, C)
0 ≤ u ≤ d -> |v|
d ≤ u     -> d(B, C)

(到其中一个端点或到支撑线的距离,取决于到该线的投影。等距轨迹由两个半圆和两条线段组成。)

式中,d为AB线段的长度,u、v分别为AB/d (AB方向的单位矢量)与AC的标量积和外积的模量。

AB.AC ≤ 0             -> |AC|
    0 ≤ AB.AC ≤ AB²   -> |ABxAC|/|AB|
          AB² ≤ AB.AC -> |BC|

忍不住用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

使用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))