我如何确定两条直线是否相交,如果相交,在x,y点处?
当前回答
基于@Gareth Rees的回答,Python版本:
import numpy as np
def np_perp( a ) :
b = np.empty_like(a)
b[0] = a[1]
b[1] = -a[0]
return b
def np_cross_product(a, b):
return np.dot(a, np_perp(b))
def np_seg_intersect(a, b, considerCollinearOverlapAsIntersect = False):
# https://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect/565282#565282
# http://www.codeproject.com/Tips/862988/Find-the-intersection-point-of-two-line-segments
r = a[1] - a[0]
s = b[1] - b[0]
v = b[0] - a[0]
num = np_cross_product(v, r)
denom = np_cross_product(r, s)
# If r x s = 0 and (q - p) x r = 0, then the two lines are collinear.
if np.isclose(denom, 0) and np.isclose(num, 0):
# 1. If either 0 <= (q - p) * r <= r * r or 0 <= (p - q) * s <= * s
# then the two lines are overlapping,
if(considerCollinearOverlapAsIntersect):
vDotR = np.dot(v, r)
aDotS = np.dot(-v, s)
if (0 <= vDotR and vDotR <= np.dot(r,r)) or (0 <= aDotS and aDotS <= np.dot(s,s)):
return True
# 2. If neither 0 <= (q - p) * r = r * r nor 0 <= (p - q) * s <= s * s
# then the two lines are collinear but disjoint.
# No need to implement this expression, as it follows from the expression above.
return None
if np.isclose(denom, 0) and not np.isclose(num, 0):
# Parallel and non intersecting
return None
u = num / denom
t = np_cross_product(v, s) / denom
if u >= 0 and u <= 1 and t >= 0 and t <= 1:
res = b[0] + (s*u)
return res
# Otherwise, the two line segments are not parallel but do not intersect.
return None
其他回答
以下是对加文回答的改进。马普的解决方案也类似,但都没有推迟分割。
这实际上也是Gareth Rees的答案的一个实际应用,因为向量积在2D中的等价是补点积,这段代码用了其中的三个。切换到3D并使用叉积,在最后插入s和t,结果是3D中直线之间的两个最近点。 不管怎样,2D解:
int get_line_intersection(float p0_x, float p0_y, float p1_x, float p1_y,
float p2_x, float p2_y, float p3_x, float p3_y, float *i_x, float *i_y)
{
float s02_x, s02_y, s10_x, s10_y, s32_x, s32_y, s_numer, t_numer, denom, t;
s10_x = p1_x - p0_x;
s10_y = p1_y - p0_y;
s32_x = p3_x - p2_x;
s32_y = p3_y - p2_y;
denom = s10_x * s32_y - s32_x * s10_y;
if (denom == 0)
return 0; // Collinear
bool denomPositive = denom > 0;
s02_x = p0_x - p2_x;
s02_y = p0_y - p2_y;
s_numer = s10_x * s02_y - s10_y * s02_x;
if ((s_numer < 0) == denomPositive)
return 0; // No collision
t_numer = s32_x * s02_y - s32_y * s02_x;
if ((t_numer < 0) == denomPositive)
return 0; // No collision
if (((s_numer > denom) == denomPositive) || ((t_numer > denom) == denomPositive))
return 0; // No collision
// Collision detected
t = t_numer / denom;
if (i_x != NULL)
*i_x = p0_x + (t * s10_x);
if (i_y != NULL)
*i_y = p0_y + (t * s10_y);
return 1;
}
基本上,它将除法延迟到最后一刻,并将大多数测试移动到某些计算完成之前,从而增加了早期退出。最后,它还避免了直线平行时的除零情况。
您可能还想考虑使用ε检验,而不是与零比较。非常接近平行的线会产生稍微偏离的结果。这不是一个bug,这是浮点数学的一个限制。
问题可以简化成这样一个问题:从A到B和从C到D的两条直线相交吗?然后你可以问它四次(在直线和矩形的四条边之间)。
这是做这个的矢量数学。假设A到B的直线就是问题中的直线C到D的直线是其中一条矩形直线。我的表示法是Ax是A的x坐标Cy是c的y坐标“*”表示点积,例如A*B = Ax*Bx + Ay*By。
E = B-A = ( Bx-Ax, By-Ay )
F = D-C = ( Dx-Cx, Dy-Cy )
P = ( -Ey, Ex )
h = ( (A-C) * P ) / ( F * P )
h是键。如果h在0和1之间,两条线相交,否则不相交。如果F*P为零,当然不能进行计算,但在这种情况下,直线是平行的,因此只有在明显的情况下才相交。
交点是C + F*h。
更多的乐趣:
如果h恰好等于0或1,两条直线的端点相交。你可以认为这是一个“交集”,也可以认为不是。
具体来说,h是直线长度乘以多少才能恰好与另一条直线相交。
因此,如果h<0,这意味着矩形线在给定直线的“后面”(“方向”是“从A到B”),如果h>1,矩形线在给定直线的“前面”。
推导:
A和C是指向直线起点的向量;E和F是由A和C端点组成的直线。
对于平面上任意两条不平行线,必须恰好有一对标量g和h,使得这个方程成立:
A + E*g = C + F*h
为什么?因为两条不平行线必须相交,这意味着你可以将这两条线按一定比例缩放并相互接触。
(起初,这看起来像一个有两个未知数的方程!但当你考虑到这是一个二维矢量方程时,它就不是,这意味着这是一对x和y的方程)
我们必须消去其中一个变量。一个简单的方法是使E项为零。要做到这一点,用一个向量对方程两边做点积这个向量与E点乘到0,我把上面的向量称为P,我做了E的明显变换。
你现在有:
A*P = C*P + F*P*h
(A-C)*P = (F*P)*h
( (A-C)*P ) / (F*P) = h
C和Objective-C
基于Gareth Rees的回答
const AGKLine AGKLineZero = (AGKLine){(CGPoint){0.0, 0.0}, (CGPoint){0.0, 0.0}};
AGKLine AGKLineMake(CGPoint start, CGPoint end)
{
return (AGKLine){start, end};
}
double AGKLineLength(AGKLine l)
{
return CGPointLengthBetween_AGK(l.start, l.end);
}
BOOL AGKLineIntersection(AGKLine l1, AGKLine l2, CGPoint *out_pointOfIntersection)
{
// http://stackoverflow.com/a/565282/202451
CGPoint p = l1.start;
CGPoint q = l2.start;
CGPoint r = CGPointSubtract_AGK(l1.end, l1.start);
CGPoint s = CGPointSubtract_AGK(l2.end, l2.start);
double s_r_crossProduct = CGPointCrossProductZComponent_AGK(r, s);
double t = CGPointCrossProductZComponent_AGK(CGPointSubtract_AGK(q, p), s) / s_r_crossProduct;
double u = CGPointCrossProductZComponent_AGK(CGPointSubtract_AGK(q, p), r) / s_r_crossProduct;
if(t < 0 || t > 1.0 || u < 0 || u > 1.0)
{
if(out_pointOfIntersection != NULL)
{
*out_pointOfIntersection = CGPointZero;
}
return NO;
}
else
{
if(out_pointOfIntersection != NULL)
{
CGPoint i = CGPointAdd_AGK(p, CGPointMultiply_AGK(r, t));
*out_pointOfIntersection = i;
}
return YES;
}
}
CGFloat CGPointCrossProductZComponent_AGK(CGPoint v1, CGPoint v2)
{
return v1.x * v2.y - v1.y * v2.x;
}
CGPoint CGPointSubtract_AGK(CGPoint p1, CGPoint p2)
{
return (CGPoint){p1.x - p2.x, p1.y - p2.y};
}
CGPoint CGPointAdd_AGK(CGPoint p1, CGPoint p2)
{
return (CGPoint){p1.x + p2.x, p1.y + p2.y};
}
CGFloat CGPointCrossProductZComponent_AGK(CGPoint v1, CGPoint v2)
{
return v1.x * v2.y - v1.y * v2.x;
}
CGPoint CGPointMultiply_AGK(CGPoint p1, CGFloat factor)
{
return (CGPoint){p1.x * factor, p1.y * factor};
}
许多函数和结构都是私有的,但是你应该很容易就能知道发生了什么。 这是公开的在这个回购https://github.com/hfossli/AGGeometryKit/
我试过其中一些答案,但它们对我不起作用(对不起伙计们);在网上搜索之后,我找到了这个。
对他的代码做了一点修改,我现在有了这个函数,它将返回交点,如果没有找到交点,它将返回- 1,1。
Public Function intercetion(ByVal ax As Integer, ByVal ay As Integer, ByVal bx As Integer, ByVal by As Integer, ByVal cx As Integer, ByVal cy As Integer, ByVal dx As Integer, ByVal dy As Integer) As Point
'// Determines the intersection point of the line segment defined by points A and B
'// with the line segment defined by points C and D.
'//
'// Returns YES if the intersection point was found, and stores that point in X,Y.
'// Returns NO if there is no determinable intersection point, in which case X,Y will
'// be unmodified.
Dim distAB, theCos, theSin, newX, ABpos As Double
'// Fail if either line segment is zero-length.
If ax = bx And ay = by Or cx = dx And cy = dy Then Return New Point(-1, -1)
'// Fail if the segments share an end-point.
If ax = cx And ay = cy Or bx = cx And by = cy Or ax = dx And ay = dy Or bx = dx And by = dy Then Return New Point(-1, -1)
'// (1) Translate the system so that point A is on the origin.
bx -= ax
by -= ay
cx -= ax
cy -= ay
dx -= ax
dy -= ay
'// Discover the length of segment A-B.
distAB = Math.Sqrt(bx * bx + by * by)
'// (2) Rotate the system so that point B is on the positive X axis.
theCos = bx / distAB
theSin = by / distAB
newX = cx * theCos + cy * theSin
cy = cy * theCos - cx * theSin
cx = newX
newX = dx * theCos + dy * theSin
dy = dy * theCos - dx * theSin
dx = newX
'// Fail if segment C-D doesn't cross line A-B.
If cy < 0 And dy < 0 Or cy >= 0 And dy >= 0 Then Return New Point(-1, -1)
'// (3) Discover the position of the intersection point along line A-B.
ABpos = dx + (cx - dx) * dy / (dy - cy)
'// Fail if segment C-D crosses line A-B outside of segment A-B.
If ABpos < 0 Or ABpos > distAB Then Return New Point(-1, -1)
'// (4) Apply the discovered position to line A-B in the original coordinate system.
'*X=Ax+ABpos*theCos
'*Y=Ay+ABpos*theSin
'// Success.
Return New Point(ax + ABpos * theCos, ay + ABpos * theSin)
End Function
FWIW,下面的函数(在C中)既检测线的交点,又确定交点。这是基于Andre LeMothe的“Tricks of the Windows Game Programming Gurus”中的一个算法。这与其他答案(例如Gareth的答案)中的一些算法并没有什么不同。然后LeMothe使用克莱默法则(不要问我)来解这些方程。
我可以证明它在我的小行星克隆中起作用,并且似乎正确地处理了Elemental, Dan和Wodzu在其他答案中描述的边缘情况。它也可能比KingNestor发布的代码快,因为它都是乘法和除法,没有平方根!
我想这里有一些除以0的可能性,尽管在我的例子中这不是问题。很容易修改以避免崩溃。
// Returns 1 if the lines intersect, otherwise 0. In addition, if the lines
// intersect the intersection point may be stored in the floats i_x and i_y.
char get_line_intersection(float p0_x, float p0_y, float p1_x, float p1_y,
float p2_x, float p2_y, float p3_x, float p3_y, float *i_x, float *i_y)
{
float s1_x, s1_y, s2_x, s2_y;
s1_x = p1_x - p0_x; s1_y = p1_y - p0_y;
s2_x = p3_x - p2_x; s2_y = p3_y - p2_y;
float s, t;
s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / (-s2_x * s1_y + s1_x * s2_y);
t = ( s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / (-s2_x * s1_y + s1_x * s2_y);
if (s >= 0 && s <= 1 && t >= 0 && t <= 1)
{
// Collision detected
if (i_x != NULL)
*i_x = p0_x + (t * s1_x);
if (i_y != NULL)
*i_y = p0_y + (t * s1_y);
return 1;
}
return 0; // No collision
}
顺便说一句,我必须说,在LeMothe的书中,虽然他显然得到了正确的算法,但他展示的具体示例插入了错误的数字,并且计算错误。例如:
(4 * (4-1) + 12 * (7-1))/(17 * 4 + 12 * 10) = 844/0.88 = 0.44
这让我困惑了好几个小时。:(