在我自己的问题线程如何计算在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);
} 

我不是要回答问题，而是要问问题，所以我希望我不会因为某些原因而得到数百万张反对票，而是批评。我只是想(并被鼓励)分享其他人的想法，因为这个帖子中的解决方案要么是用一些奇异的语言(Fortran, Mathematica)，要么被某人标记为错误。对我来说唯一有用的(由Grumdrig编写)是用c++编写的，没有人标记它有错误。但是它缺少被调用的方法(dot等)。

下面是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;
}

Lua解决方案

-- distance from point (px, py) to line segment (x1, y1, x2, y2)
function distPointToLine(px,py,x1,y1,x2,y2) -- point, start and end of the segment
    local dx,dy = x2-x1,y2-y1
    local length = math.sqrt(dx*dx+dy*dy)
    dx,dy = dx/length,dy/length -- normalization
    local p = dx*(px-x1)+dy*(py-y1)
    if p < 0 then
        dx,dy = px-x1,py-y1
        return math.sqrt(dx*dx+dy*dy), x1, y1 -- distance, nearest point
    elseif p > length then
        dx,dy = px-x2,py-y2
        return math.sqrt(dx*dx+dy*dy), x2, y2 -- distance, nearest point
    end
    return math.abs(dy*(px-x1)-dx*(py-y1)), x1+dx*p, y1+dy*p -- distance, nearest point
end

对于折线(有两条以上线段的线):

-- if the (poly-)line has several segments, just iterate through all of them:
function nearest_sector_in_line (x, y, line)
    local x1, y1, x2, y2, min_dist
    local ax,ay = line[1], line[2]
    for j = 3, #line-1, 2 do
        local bx,by = line[j], line[j+1]
        local dist = distPointToLine(x,y,ax,ay,bx,by)
        if not min_dist or dist < min_dist then
            min_dist = dist
            x1, y1, x2, y2 = ax,ay,bx,by
        end
        ax, ay = bx, by
    end
    return x1, y1, x2, y2
end

例子:

-- call it:
local x1, y1, x2, y2 = nearest_sector_in_line (7, 4, {0,0, 10,0, 10,10, 0,10})

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)
}



 

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

以下是Grumdrig解决方案的一个更完整的说明。这个版本还返回最近的点本身。

#include "stdio.h"
#include "math.h"

class Vec2
{
public:
    float _x;
    float _y;

    Vec2()
    {
        _x = 0;
        _y = 0;
    }

    Vec2( const float x, const float y )
    {
        _x = x;
        _y = y;
    }

    Vec2 operator+( const Vec2 &v ) const
    {
        return Vec2( this->_x + v._x, this->_y + v._y );
    }

    Vec2 operator-( const Vec2 &v ) const
    {
        return Vec2( this->_x - v._x, this->_y - v._y );
    }

    Vec2 operator*( const float f ) const
    {
        return Vec2( this->_x * f, this->_y * f );
    }

    float DistanceToSquared( const Vec2 p ) const
    {
        const float dX = p._x - this->_x;
        const float dY = p._y - this->_y;

        return dX * dX + dY * dY;
    }

    float DistanceTo( const Vec2 p ) const
    {
        return sqrt( this->DistanceToSquared( p ) );
    }

    float DotProduct( const Vec2 p ) const
    {
        return this->_x * p._x + this->_y * p._y;
    }
};

// return minimum distance between line segment vw and point p, and the closest point on the line segment, q
float DistanceFromLineSegmentToPoint( const Vec2 v, const Vec2 w, const Vec2 p, Vec2 * const q )
{
    const float distSq = v.DistanceToSquared( w ); // i.e. |w-v|^2 ... avoid a sqrt
    if ( distSq == 0.0 )
    {
        // v == w case
        (*q) = v;

        return v.DistanceTo( p );
    }

    // 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

    const float t = ( p - v ).DotProduct( w - v ) / distSq;
    if ( t < 0.0 )
    {
        // beyond the v end of the segment
        (*q) = v;

        return v.DistanceTo( p );
    }
    else if ( t > 1.0 )
    {
        // beyond the w end of the segment
        (*q) = w;

        return w.DistanceTo( p );
    }

    // projection falls on the segment
    const Vec2 projection = v + ( ( w - v ) * t );

    (*q) = projection;

    return p.DistanceTo( projection );
}

float DistanceFromLineSegmentToPoint( float segmentX1, float segmentY1, float segmentX2, float segmentY2, float pX, float pY, float *qX, float *qY )
{
    Vec2 q;

    float distance = DistanceFromLineSegmentToPoint( Vec2( segmentX1, segmentY1 ), Vec2( segmentX2, segmentY2 ), Vec2( pX, pY ), &q );

    (*qX) = q._x;
    (*qY) = q._y;

    return distance;
}

void TestDistanceFromLineSegmentToPoint( float segmentX1, float segmentY1, float segmentX2, float segmentY2, float pX, float pY )
{
    float qX;
    float qY;
    float d = DistanceFromLineSegmentToPoint( segmentX1, segmentY1, segmentX2, segmentY2, pX, pY, &qX, &qY );
    printf( "line segment = ( ( %f, %f ), ( %f, %f ) ), p = ( %f, %f ), distance = %f, q = ( %f, %f )\n",
            segmentX1, segmentY1, segmentX2, segmentY2, pX, pY, d, qX, qY );
}

void TestDistanceFromLineSegmentToPoint()
{
    TestDistanceFromLineSegmentToPoint( 0, 0, 1, 1, 1, 0 );
    TestDistanceFromLineSegmentToPoint( 0, 0, 20, 10, 5, 4 );
    TestDistanceFromLineSegmentToPoint( 0, 0, 20, 10, 30, 15 );
    TestDistanceFromLineSegmentToPoint( 0, 0, 20, 10, -30, 15 );
    TestDistanceFromLineSegmentToPoint( 0, 0, 10, 0, 5, 1 );
    TestDistanceFromLineSegmentToPoint( 0, 0, 0, 10, 1, 5 );
}