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

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


当前回答

GLSL版:

// line (a -> b ) point p[enter image description here][1]
float distanceToLine(vec2 a, vec2 b, vec2 p) {
    float aside = dot((p - a),(b - a));
    if(aside< 0.0) return length(p-a);
    float bside = dot((p - b),(a - b));
    if(bside< 0.0) return length(p-b);
    vec2 pointOnLine = (bside*a + aside*b)/pow(length(a-b),2.0);
    return length(p - pointOnLine);
}

其他回答

这里是与c++答案相同的东西,但移植到pascal。点参数的顺序已经改变,以适应我的代码,但还是一样的东西。

function Dot(const p1, p2: PointF): double;
begin
  Result := p1.x * p2.x + p1.y * p2.y;
end;
function SubPoint(const p1, p2: PointF): PointF;
begin
  result.x := p1.x - p2.x;
  result.y := p1.y - p2.y;
end;

function ShortestDistance2(const p,v,w : PointF) : double;
var
  l2,t : double;
  projection,tt: PointF;
begin
  // Return minimum distance between line segment vw and point p
  //l2 := length_squared(v, w);  // i.e. |w-v|^2 -  avoid a sqrt
  l2 := Distance(v,w);
  l2 := MPower(l2,2);
  if (l2 = 0.0) then begin
    result:= Distance(p, v);   // v == w case
    exit;
  end;
  // 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 := Dot(SubPoint(p,v),SubPoint(w,v)) / l2;
  if (t < 0.0) then begin
    result := Distance(p, v);       // Beyond the 'v' end of the segment
    exit;
  end
  else if (t > 1.0) then begin
    result := Distance(p, w);  // Beyond the 'w' end of the segment
    exit;
  end;
  //projection := v + t * (w - v);  // Projection falls on the segment
  tt.x := v.x + t * (w.x - v.x);
  tt.y := v.y + t * (w.y - v.y);
  result := Distance(p, tt);
end;

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

这是一个为有限线段而做的实现,而不是像这里的大多数其他函数那样的无限线(这就是为什么我做这个)。

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

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;

    }

本想在GLSL中这样做,但如果可能的话,最好避免所有这些条件。使用clamp()可以避免两种端点情况:

// find closest point to P on line segment AB:
vec3 closest_point_on_line_segment(in vec3 P, in vec3 A, in vec3 B) {
    vec3 AP = P - A, AB = B - A;
    float l = dot(AB, AB);
    if (l <= 0.0000001) return A;    // A and B are practically the same
    return AP - AB*clamp(dot(AP, AB)/l, 0.0, 1.0);  // do the projection
}

如果您可以确定A和B彼此不会非常接近,则可以简化为删除If()。事实上,即使A和B是相同的,我的GPU仍然给出了这个无条件版本的正确结果(但这是使用pre-OpenGL 4.1,其中GLSL除零是未定义的):

// find closest point to P on line segment AB:
vec3 closest_point_on_line_segment(in vec3 P, in vec3 A, in vec3 B) {
    vec3 AP = P - A, AB = B - A;
    return AP - AB*clamp(dot(AP, AB)/dot(AB, AB), 0.0, 1.0);
}

计算距离是很简单的——GLSL提供了一个distance()函数,你可以在这个最近的点和P。

灵感来自Iñigo Quilez的胶囊距离函数代码

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