Наименьшее расстояние между точкой и отрезком - PullRequest
323 голосов
/ 11 мая 2009

Мне нужна основная функция, чтобы найти кратчайшее расстояние между точкой и отрезком. Не стесняйтесь написать решение на любом языке, который вы хотите; Я могу перевести это в то, что я использую (Javascript).

РЕДАКТИРОВАТЬ: мой сегмент линии определяется двумя конечными точками. Таким образом, мой отрезок AB определяется двумя точками A (x1,y1) и B (x2,y2). Я пытаюсь найти расстояние между этим отрезком и точкой C (x3,y3). Мои навыки геометрии ржавые, поэтому примеры, которые я видел, сбивают с толку, извините, признаюсь.

Ответы [ 49 ]

8 голосов
/ 09 июня 2012

Вот более полное изложение решения Грумдрига. Эта версия также возвращает самую близкую точку.

#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 );
}
7 голосов
/ 19 ноября 2013

Однолинейное решение с использованием арктангенсов:

Идея состоит в том, чтобы переместить A в (0, 0) и повернуть треугольник по часовой стрелке, чтобы C лежал на оси X, когда это произойдет, По будет расстояние.

  1. угол = Атан (Cy - Ay, Cx - Ax);
  2. b угол = Атан (By - Ay, Bx - Axe);
  3. длина AB = Sqrt ((Bx - Axe) ^ 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))
7 голосов
/ 11 июля 2012

Рассмотрите эту модификацию ответа Грамдрига выше. Много раз вы обнаружите, что неточность с плавающей запятой может вызвать проблемы. Я использую двойники в версии ниже, но вы можете легко перейти на плавающие. Важной частью является то, что он использует эпсилон для обработки "помои". Кроме того, вы много раз захотите узнать, ГДЕ произошло пересечение или произошло ли оно вообще. Если возвращенное значение t <0,0 или> 1,0, столкновения не произошло. Тем не менее, даже если столкновения не произошло, много раз вы захотите узнать, где находится ближайшая точка на сегменте к P, и поэтому я использую qx и qy, чтобы вернуть это местоположение.

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));
    }
}
6 голосов
/ 11 мая 2009

Я предполагаю, что вы хотите найти кратчайшее расстояние между точкой и отрезком; Для этого вам нужно найти линию (lineA), которая перпендикулярна вашему отрезку (lineB), который проходит через вашу точку, определить пересечение между этой линией (lineA) и вашей линией, проходящей через ваш отрезок line (lineB) ; если эта точка находится между двумя точками вашего отрезка, то расстояние - это расстояние между вашей точкой и точкой, которую вы только что нашли, которая является пересечением линии A и линии B; если точка не находится между двумя точками вашего отрезка, вам нужно получить расстояние между вашей точкой и ближе к двум концам отрезка; это можно легко сделать, взяв квадратное расстояние (чтобы избежать квадратного корня) между точкой и двумя точками отрезка; в зависимости от того, что ближе, возьмите квадратный корень этого.

5 голосов
/ 11 августа 2014

Реализация Grumdrig на C ++ / JavaScript была для меня очень полезной, поэтому я предоставил прямой порт Python, который я использую. Полный код здесь .

class Point(object):
  def __init__(self, x, y):
    self.x = float(x)
    self.y = float(y)

def square(x):
  return x * x

def distance_squared(v, w):
  return square(v.x - w.x) + square(v.y - w.y)

def distance_point_segment_squared(p, v, w):
  # Segment length squared, |w-v|^2
  d2 = distance_squared(v, w) 
  if d2 == 0: 
    # v == w, return distance to v
    return distance_squared(p, v)
  # 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
  t = ((p.x - v.x) * (w.x - v.x) + (p.y - v.y) * (w.y - v.y)) / d2;
  if t < 0:
    # Beyond v end of the segment
    return distance_squared(p, v)
  elif t > 1.0:
    # Beyond w end of the segment
    return distance_squared(p, w)
  else:
    # Projection falls on the segment.
    proj = Point(v.x + t * (w.x - v.x), v.y + t * (w.y - v.y))
    # print proj.x, proj.y
    return distance_squared(p, proj)
5 голосов
/ 28 февраля 2012

Код Matlab со встроенным «самотестированием», если они вызывают функцию без аргументов:

function r = distPointToLineSegment( xy0, xy1, xyP )
% r = distPointToLineSegment( xy0, xy1, xyP )

if( nargin < 3 )
    selfTest();
    r=0;
else
    vx = xy0(1)-xyP(1);
    vy = xy0(2)-xyP(2);
    ux = xy1(1)-xy0(1);
    uy = xy1(2)-xy0(2);
    lenSqr= (ux*ux+uy*uy);
    detP= -vx*ux + -vy*uy;

    if( detP < 0 )
        r = norm(xy0-xyP,2);
    elseif( detP > lenSqr )
        r = norm(xy1-xyP,2);
    else
        r = abs(ux*vy-uy*vx)/sqrt(lenSqr);
    end
end


    function selfTest()
        %#ok<*NASGU>
        disp(['invalid args, distPointToLineSegment running (recursive)  self-test...']);

        ptA = [1;1]; ptB = [-1;-1];
        ptC = [1/2;1/2];  % on the line
        ptD = [-2;-1.5];  % too far from line segment
        ptE = [1/2;0];    % should be same as perpendicular distance to line
        ptF = [1.5;1.5];      % along the A-B but outside of the segment

        distCtoAB = distPointToLineSegment(ptA,ptB,ptC)
        distDtoAB = distPointToLineSegment(ptA,ptB,ptD)
        distEtoAB = distPointToLineSegment(ptA,ptB,ptE)
        distFtoAB = distPointToLineSegment(ptA,ptB,ptF)
        figure(1); clf;
        circle = @(x, y, r, c) rectangle('Position', [x-r, y-r, 2*r, 2*r], ...
            'Curvature', [1 1], 'EdgeColor', c);
        plot([ptA(1) ptB(1)],[ptA(2) ptB(2)],'r-x'); hold on;
        plot(ptC(1),ptC(2),'b+'); circle(ptC(1),ptC(2), 0.5e-1, 'b');
        plot(ptD(1),ptD(2),'g+'); circle(ptD(1),ptD(2), distDtoAB, 'g');
        plot(ptE(1),ptE(2),'k+'); circle(ptE(1),ptE(2), distEtoAB, 'k');
        plot(ptF(1),ptF(2),'m+'); circle(ptF(1),ptF(2), distFtoAB, 'm');
        hold off;
        axis([-3 3 -3 3]); axis equal;
    end

end
4 голосов
/ 16 мая 2012

А теперь и мое решение ...... (Javascript)

Это очень быстро, потому что я стараюсь избегать любых функций Math.pow.

Как видите, в конце функции у меня есть расстояние до линии.

код из библиотеки http://www.draw2d.org/graphiti/jsdoc/#!/example

/**
 * 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;
};
4 голосов
/ 17 октября 2012

закодировано в t-sql

точка равна (@px, @py), а отрезок проходит от (@ax, @ay) до (@bx, @by)

create function fn_sqr (@NumberToSquare decimal(18,10)) 
returns decimal(18,10)
as 
begin
    declare @Result decimal(18,10)
    set @Result = @NumberToSquare * @NumberToSquare
    return @Result
end
go

create function fn_Distance(@ax decimal (18,10) , @ay decimal (18,10), @bx decimal(18,10),  @by decimal(18,10)) 
returns decimal(18,10)
as
begin
    declare @Result decimal(18,10)
    set @Result = (select dbo.fn_sqr(@ax - @bx) + dbo.fn_sqr(@ay - @by) )
    return @Result
end
go

create function fn_DistanceToSegmentSquared(@px decimal(18,10), @py decimal(18,10), @ax decimal(18,10), @ay decimal(18,10), @bx decimal(18,10), @by decimal(18,10)) 
returns decimal(18,10)
as 
begin
    declare @l2 decimal(18,10)
    set @l2 = (select dbo.fn_Distance(@ax, @ay, @bx, @by))
    if @l2 = 0
        return dbo.fn_Distance(@px, @py, @ax, @ay)
    declare @t decimal(18,10)
    set @t = ((@px - @ax) * (@bx - @ax) + (@py - @ay) * (@by - @ay)) / @l2
    if (@t < 0) 
        return dbo.fn_Distance(@px, @py, @ax, @ay);
    if (@t > 1) 
        return dbo.fn_Distance(@px, @py, @bx, @by);
    return dbo.fn_Distance(@px, @py,  @ax + @t * (@bx - @ax),  @ay + @t * (@by - @ay))
end
go

create function fn_DistanceToSegment(@px decimal(18,10), @py decimal(18,10), @ax decimal(18,10), @ay decimal(18,10), @bx decimal(18,10), @by decimal(18,10)) 
returns decimal(18,10)
as 
begin
    return sqrt(dbo.fn_DistanceToSegmentSquared(@px, @py , @ax , @ay , @bx , @by ))
end
go

--example execution for distance from a point at (6,1) to line segment that runs from (4,2) to (2,1)
select dbo.fn_DistanceToSegment(6, 1, 4, 2, 2, 1) 
--result = 2.2360679775

--example execution for distance from a point at (-3,-2) to line segment that runs from (0,-2) to (-2,1)
select dbo.fn_DistanceToSegment(-3, -2, 0, -2, -2, 1) 
--result = 2.4961508830

--example execution for distance from a point at (0,-2) to line segment that runs from (0,-2) to (-2,1)
select dbo.fn_DistanceToSegment(0,-2, 0, -2, -2, 1) 
--result = 0.0000000000
4 голосов
/ 01 апреля 2013

Похоже, что почти все остальные в StackOverflow предоставили ответ (пока 23 ответа), так что вот мой вклад для C #. В основном это основано на ответе М. Каца, который, в свою очередь, основан на ответе Грумдрига.

   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%, но я даже не уверен в этом - различия в различных тестовых прогонах были одного порядка. Я также попытался измерить версию, опубликованную Матти (плюс очевидная оптимизация), и эта версия, кажется, примерно на 4% медленнее, чем версия, основанная на коде Каца / Грумдрига.

Редактировать: Между прочим, я также пытался измерить метод, который находит расстояние до бесконечной линии (не отрезка), используя перекрестное произведение (и Sqrt ()), и это примерно на 32% быстрее. *

3 голосов
/ 24 марта 2014

Вот версия devnullicus's 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;
}
...