我需要一个基本函数来求点到线段的最短距离。你可以随意用任何你想要的语言来编写解决方案;我可以把它翻译成什么我正在使用(Javascript)。
编辑:线段由两个端点定义。线段AB由两点A (x1,y1)和B (x2,y2)定义。我要求的是这条线段到点C (x3,y3)的距离。我的几何技能生疏了,所以我看到的例子让我很困惑,我很遗憾地承认。
我需要一个基本函数来求点到线段的最短距离。你可以随意用任何你想要的语言来编写解决方案;我可以把它翻译成什么我正在使用(Javascript)。
编辑:线段由两个端点定义。线段AB由两点A (x1,y1)和B (x2,y2)定义。我要求的是这条线段到点C (x3,y3)的距离。我的几何技能生疏了,所以我看到的例子让我很困惑,我很遗憾地承认。
当前回答
这里是与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;
其他回答
使用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))
下面是HSQLDB的SQL实现:
CREATE FUNCTION dist_to_segment(px double, py double, vx double, vy double, wx double, wy double)
RETURNS double
BEGIN atomic
declare l2 double;
declare t double;
declare nx double;
declare ny double;
set l2 =(vx - wx)*(vx - wx) + (vy - wy)*(vy - wy);
IF l2 = 0 THEN
RETURN sqrt((vx - px)*(vx - px) + (vy - py)*(vy - py));
ELSE
set t = ((px - vx) * (wx - vx) + (py - vy) * (wy - vy)) / l2;
set t = GREATEST(0, LEAST(1, t));
set nx=vx + t * (wx - vx);
set ny=vy + t * (wy - vy);
RETURN sqrt((nx - px)*(nx - px) + (ny - py)*(ny - py));
END IF;
END;
Postgres的实现:
CREATE FUNCTION dist_to_segment(px numeric, py numeric, vx numeric, vy numeric, wx numeric, wy numeric)
RETURNS numeric
AS $$
declare l2 numeric;
declare t numeric;
declare nx numeric;
declare ny numeric;
BEGIN
l2 := (vx - wx)*(vx - wx) + (vy - wy)*(vy - wy);
IF l2 = 0 THEN
RETURN sqrt((vx - px)*(vx - px) + (vy - py)*(vy - py));
ELSE
t := ((px - vx) * (wx - vx) + (py - vy) * (wy - vy)) / l2;
t := GREATEST(0, LEAST(1, t));
nx := vx + t * (wx - vx);
ny := vy + t * (wy - vy);
RETURN sqrt((nx - px)*(nx - px) + (ny - py)*(ny - py));
END IF;
END;
$$ LANGUAGE plpgsql;
忍不住用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
特征c++版本的3D线段和点
// Return minimum distance between line segment: head--->tail and point
double MinimumDistance(Eigen::Vector3d head, Eigen::Vector3d tail,Eigen::Vector3d point)
{
double l2 = std::pow((head - tail).norm(),2);
if(l2 ==0.0) return (head - point).norm();// head == tail case
// Consider the line extending the segment, parameterized as head + t (tail - point).
// We find projection of point onto the line.
// It falls where t = [(point-head) . (tail-head)] / |tail-head|^2
// We clamp t from [0,1] to handle points outside the segment head--->tail.
double t = max(0,min(1,(point-head).dot(tail-head)/l2));
Eigen::Vector3d projection = head + t*(tail-head);
return (point - projection).norm();
}
在javascript中使用几何:
var a = { x:20, y:20};//start segment
var b = { x:40, y:30};//end segment
var c = { x:37, y:14};//point
// magnitude from a to c
var ac = Math.sqrt( Math.pow( ( a.x - c.x ), 2 ) + Math.pow( ( a.y - c.y ), 2) );
// magnitude from b to c
var bc = Math.sqrt( Math.pow( ( b.x - c.x ), 2 ) + Math.pow( ( b.y - c.y ), 2 ) );
// magnitude from a to b (base)
var ab = Math.sqrt( Math.pow( ( a.x - b.x ), 2 ) + Math.pow( ( a.y - b.y ), 2 ) );
// perimeter of triangle
var p = ac + bc + ab;
// area of the triangle
var area = Math.sqrt( p/2 * ( p/2 - ac) * ( p/2 - bc ) * ( p/2 - ab ) );
// height of the triangle = distance
var h = ( area * 2 ) / ab;
console.log ("height: " + h);