我有一条从A到B的线段,以及一个位于C处、半径为R的圆。
如何使用良好的算法检查该线段是否与圆相交?并且在圆的边缘上发生交点的坐标是什么?
我有一条从A到B的线段,以及一个位于C处、半径为R的圆。
如何使用良好的算法检查该线段是否与圆相交?并且在圆的边缘上发生交点的坐标是什么?
我知道这个帖子已经有一段时间了。根据chmike的回答并由Aqib Mumtaz改进,他们给出了一个很好的答案,但只适用于无限线,正如Aqib所说。因此,我添加了一些比较来确定线段是否触碰到圆,我用Python编写了它。
def LineIntersectCircle(c, r, p1, p2):
#p1 is the first line point
#p2 is the second line point
#c is the circle's center
#r is the circle's radius
p3 = [p1[0]-c[0], p1[1]-c[1]]
p4 = [p2[0]-c[0], p2[1]-c[1]]
m = (p4[1] - p3[1]) / (p4[0] - p3[0])
b = p3[1] - m * p3[0]
underRadical = math.pow(r,2)*math.pow(m,2) + math.pow(r,2) - math.pow(b,2)
if (underRadical < 0):
print("NOT")
else:
t1 = (-2*m*b+2*math.sqrt(underRadical)) / (2 * math.pow(m,2) + 2)
t2 = (-2*m*b-2*math.sqrt(underRadical)) / (2 * math.pow(m,2) + 2)
i1 = [t1+c[0], m * t1 + b + c[1]]
i2 = [t2+c[0], m * t2 + b + c[1]]
if p1[0] > p2[0]: #Si el punto 1 es mayor al 2 en X
if (i1[0] < p1[0]) and (i1[0] > p2[0]): #Si el punto iX esta entre 2 y 1 en X
if p1[1] > p2[1]: #Si el punto 1 es mayor al 2 en Y
if (i1[1] < p1[1]) and (i1[1] > p2[1]): #Si el punto iy esta entre 2 y 1
print("Intersection")
if p1[1] < p2[1]: #Si el punto 2 es mayo al 2 en Y
if (i1[1] > p1[1]) and (i1[1] < p2[1]): #Si el punto iy esta entre 1 y 2
print("Intersection")
if p1[0] < p2[0]: #Si el punto 2 es mayor al 1 en X
if (i1[0] > p1[0]) and (i1[0] < p2[0]): #Si el punto iX esta entre 1 y 2 en X
if p1[1] > p2[1]: #Si el punto 1 es mayor al 2 en Y
if (i1[1] < p1[1]) and (i1[1] > p2[1]): #Si el punto iy esta entre 2 y 1
print("Intersection")
if p1[1] < p2[1]: #Si el punto 2 es mayo al 2 en Y
if (i1[1] > p1[1]) and (i1[1] < p2[1]): #Si el punto iy esta entre 1 y 2
print("Intersection")
if p1[0] > p2[0]: #Si el punto 1 es mayor al 2 en X
if (i2[0] < p1[0]) and (i2[0] > p2[0]): #Si el punto iX esta entre 2 y 1 en X
if p1[1] > p2[1]: #Si el punto 1 es mayor al 2 en Y
if (i2[1] < p1[1]) and (i2[1] > p2[1]): #Si el punto iy esta entre 2 y 1
print("Intersection")
if p1[1] < p2[1]: #Si el punto 2 es mayo al 2 en Y
if (i2[1] > p1[1]) and (i2[1] < p2[1]): #Si el punto iy esta entre 1 y 2
print("Intersection")
if p1[0] < p2[0]: #Si el punto 2 es mayor al 1 en X
if (i2[0] > p1[0]) and (i2[0] < p2[0]): #Si el punto iX esta entre 1 y 2 en X
if p1[1] > p2[1]: #Si el punto 1 es mayor al 2 en Y
if (i2[1] < p1[1]) and (i2[1] > p2[1]): #Si el punto iy esta entre 2 y 1
print("Intersection")
if p1[1] < p2[1]: #Si el punto 2 es mayo al 2 en Y
if (i2[1] > p1[1]) and (i2[1] < p2[1]): #Si el punto iy esta entre 1 y 2
print("Intersection")
我只需要这个,所以我想出了这个解决方案。这种语言是maxscript,但应该很容易翻译成任何其他语言。 sideA、sideB和CircleRadius都是标量,其余变量都是点,格式为[x,y,z]。我假设z=0以在平面XY上解决问题。
fn projectPoint p1 p2 p3 = --project p1 perpendicular to the line p2-p3
(
local v= normalize (p3-p2)
local p= (p1-p2)
p2+((dot v p)*v)
)
fn findIntersectionLineCircle CircleCenter CircleRadius LineP1 LineP2=
(
pp=projectPoint CircleCenter LineP1 LineP2
sideA=distance pp CircleCenter
--use pythagoras to solve the third side
sideB=sqrt(CircleRadius^2-sideA^2) -- this will return NaN if they don't intersect
IntersectV=normalize (pp-CircleCenter)
perpV=[IntersectV.y,-IntersectV.x,IntersectV.z]
--project the point to both sides to find the solutions
solution1=pp+(sideB*perpV)
solution2=pp-(sideB*perpV)
return #(solution1,solution2)
)
这个Java函数返回一个DVec2对象。它需要一个DVec2作为圆心,圆的半径和一条线。
public static DVec2 CircLine(DVec2 C, double r, Line line)
{
DVec2 A = line.p1;
DVec2 B = line.p2;
DVec2 P;
DVec2 AC = new DVec2( C );
AC.sub(A);
DVec2 AB = new DVec2( B );
AB.sub(A);
double ab2 = AB.dot(AB);
double acab = AC.dot(AB);
double t = acab / ab2;
if (t < 0.0)
t = 0.0;
else if (t > 1.0)
t = 1.0;
//P = A + t * AB;
P = new DVec2( AB );
P.mul( t );
P.add( A );
DVec2 H = new DVec2( P );
H.sub( C );
double h2 = H.dot(H);
double r2 = r * r;
if(h2 > r2)
return null;
else
return P;
}
这是我用TypeScript编写的解决方案,遵循了@Mizipzor提出的想法(使用投影):
/**
* Determines whether a line segment defined by a start and end point intersects with a sphere defined by a center point and a radius
* @param a the start point of the line segment
* @param b the end point of the line segment
* @param c the center point of the sphere
* @param r the radius of the sphere
*/
export function lineSphereIntersects(
a: IPoint,
b: IPoint,
c: IPoint,
r: number
): boolean {
// find the three sides of the triangle formed by the three points
const ab: number = distance(a, b);
const ac: number = distance(a, c);
const bc: number = distance(b, c);
// check to see if either ends of the line segment are inside of the sphere
if (ac < r || bc < r) {
return true;
}
// find the angle between the line segment and the center of the sphere
const numerator: number = Math.pow(ac, 2) + Math.pow(ab, 2) - Math.pow(bc, 2);
const denominator: number = 2 * ac * ab;
const cab: number = Math.acos(numerator / denominator);
// find the distance from the center of the sphere and the line segment
const cd: number = Math.sin(cab) * ac;
// if the radius is at least as long as the distance between the center and the line
if (r >= cd) {
// find the distance between the line start and the point on the line closest to
// the center of the sphere
const ad: number = Math.cos(cab) * ac;
// intersection occurs when the point on the line closest to the sphere center is
// no further away than the end of the line
return ad <= ab;
}
return false;
}
export function distance(a: IPoint, b: IPoint): number {
return Math.sqrt(
Math.pow(b.z - a.z, 2) + Math.pow(b.y - a.y, 2) + Math.pow(b.x - a.x, 2)
);
}
export interface IPoint {
x: number;
y: number;
z: number;
}
并且反过来x' = x * cos(theta) + y * sin(theta)
y' = - x * sin(theta) + y * cos(theta)
考虑线段x = x' * cos(theta) - y' * sin(theta)
y = x' * sin(theta) + y' * cos(theta)
P1-P2
的方向(在 -x 方向上为 45°),我们可以取 theta=45°
。将旋转的第二个方程代入 x-y 系的圆的方程:x^2+y^2=1
,经过简单的运算,我们得到了在 x'-y' 系中相同的圆的方程:x'^2+y'^2=1
。基于 @Joe Skeen 的 Python 解决方案
def check_line_segment_circle_intersection(line, point, radious):
""" Checks whether a point intersects with a line defined by two points.
A `point` is list with two values: [2, 3]
A `line` is list with two points: [point1, point2]
"""
line_distance = distance(line[0], line[1])
distance_start_to_point = distance(line[0], point)
distance_end_to_point = distance(line[1], point)
if (distance_start_to_point <= radious or distance_end_to_point <= radious):
return True
# angle between line and point with law of cosines
numerator = (math.pow(distance_start_to_point, 2)
+ math.pow(line_distance, 2)
- math.pow(distance_end_to_point, 2))
denominator = 2 * distance_start_to_point * line_distance
ratio = numerator / denominator
ratio = ratio if ratio <= 1 else 1 # To account for float errors
ratio = ratio if ratio >= -1 else -1 # To account for float errors
angle = math.acos(ratio)
# distance from the point to the line with sin projection
distance_line_to_point = math.sin(angle) * distance_start_to_point
if distance_line_to_point <= radious:
point_projection_in_line = math.cos(angle) * distance_start_to_point
# Intersection occurs whent the point projection in the line is less
# than the line distance and positive
return point_projection_in_line <= line_distance and point_projection_in_line >= 0
return False
def distance(point1, point2):
return math.sqrt(
math.pow(point1[1] - point2[1], 2) +
math.pow(point1[0] - point2[0], 2)
)
这里有一个用golang编写的解决方案。该方法与此处发布的其他答案类似,但并不完全相同。它易于实现,并已经过测试。以下是步骤:
这里推导出二次方程的A、B和C值,其中(n-et)和(m-dt)分别是线的x和y坐标的方程。r是圆的半径。
(n-et)(n-et) + (m-dt)(m-dt) = rr
nn - 2etn + etet + mm - 2mdt + dtdt = rr
(ee+dd)tt - 2(en + dm)t + nn + mm - rr = 0
因此A = ee + dd,B = -2(en + dm),C = nn + mm - rr。
以下是该函数的golang代码:
package geom
import (
"math"
)
// SegmentCircleIntersection return points of intersection between a circle and
// a line segment. The Boolean intersects returns true if one or
// more solutions exist. If only one solution exists,
// x1 == x2 and y1 == y2.
// s1x and s1y are coordinates for one end point of the segment, and
// s2x and s2y are coordinates for the other end of the segment.
// cx and cy are the coordinates of the center of the circle and
// r is the radius of the circle.
func SegmentCircleIntersection(s1x, s1y, s2x, s2y, cx, cy, r float64) (x1, y1, x2, y2 float64, intersects bool) {
// (n-et) and (m-dt) are expressions for the x and y coordinates
// of a parameterized line in coordinates whose origin is the
// center of the circle.
// When t = 0, (n-et) == s1x - cx and (m-dt) == s1y - cy
// When t = 1, (n-et) == s2x - cx and (m-dt) == s2y - cy.
n := s2x - cx
m := s2y - cy
e := s2x - s1x
d := s2y - s1y
// lineFunc checks if the t parameter is in the segment and if so
// calculates the line point in the unshifted coordinates (adds back
// cx and cy.
lineFunc := func(t float64) (x, y float64, inBounds bool) {
inBounds = t >= 0 && t <= 1 // Check bounds on closed segment
// To check bounds for an open segment use t > 0 && t < 1
if inBounds { // Calc coords for point in segment
x = n - e*t + cx
y = m - d*t + cy
}
return
}
// Since we want the points on the line distance r from the origin,
// (n-et)(n-et) + (m-dt)(m-dt) = rr.
// Expanding and collecting terms yeilds the following quadratic equation:
A, B, C := e*e+d*d, -2*(e*n+m*d), n*n+m*m-r*r
D := B*B - 4*A*C // discriminant of quadratic
if D < 0 {
return // No solution
}
D = math.Sqrt(D)
var p1In, p2In bool
x1, y1, p1In = lineFunc((-B + D) / (2 * A)) // First root
if D == 0.0 {
intersects = p1In
x2, y2 = x1, y1
return // Only possible solution, quadratic has one root.
}
x2, y2, p2In = lineFunc((-B - D) / (2 * A)) // Second root
intersects = p1In || p2In
if p1In == false { // Only x2, y2 may be valid solutions
x1, y1 = x2, y2
} else if p2In == false { // Only x1, y1 are valid solutions
x2, y2 = x1, y1
}
return
}
package geom_test
import (
"testing"
. "**put your package path here**"
)
func CheckEpsilon(t *testing.T, v, epsilon float64, message string) {
if v > epsilon || v < -epsilon {
t.Error(message, v, epsilon)
t.FailNow()
}
}
func TestSegmentCircleIntersection(t *testing.T) {
epsilon := 1e-10 // Something smallish
x1, y1 := 5.0, 2.0 // segment end point 1
x2, y2 := 50.0, 30.0 // segment end point 2
cx, cy := 100.0, 90.0 // center of circle
r := 80.0
segx, segy := x2-x1, y2-y1
testCntr, solutionCntr := 0, 0
for i := -100; i < 100; i++ {
for j := -100; j < 100; j++ {
testCntr++
s1x, s2x := x1+float64(i), x2+float64(i)
s1y, s2y := y1+float64(j), y2+float64(j)
sc1x, sc1y := s1x-cx, s1y-cy
seg1Inside := sc1x*sc1x+sc1y*sc1y < r*r
sc2x, sc2y := s2x-cx, s2y-cy
seg2Inside := sc2x*sc2x+sc2y*sc2y < r*r
p1x, p1y, p2x, p2y, intersects := SegmentCircleIntersection(s1x, s1y, s2x, s2y, cx, cy, r)
if intersects {
solutionCntr++
//Check if points are on circle
c1x, c1y := p1x-cx, p1y-cy
deltaLen1 := (c1x*c1x + c1y*c1y) - r*r
CheckEpsilon(t, deltaLen1, epsilon, "p1 not on circle")
c2x, c2y := p2x-cx, p2y-cy
deltaLen2 := (c2x*c2x + c2y*c2y) - r*r
CheckEpsilon(t, deltaLen2, epsilon, "p2 not on circle")
// Check if points are on the line through the line segment
// "cross product" of vector from a segment point to the point
// and the vector for the segment should be near zero
vp1x, vp1y := p1x-s1x, p1y-s1y
crossProd1 := vp1x*segy - vp1y*segx
CheckEpsilon(t, crossProd1, epsilon, "p1 not on line ")
vp2x, vp2y := p2x-s1x, p2y-s1y
crossProd2 := vp2x*segy - vp2y*segx
CheckEpsilon(t, crossProd2, epsilon, "p2 not on line ")
// Check if point is between points s1 and s2 on line
// This means the sign of the dot prod of the segment vector
// and point to segment end point vectors are opposite for
// either end.
wp1x, wp1y := p1x-s2x, p1y-s2y
dp1v := vp1x*segx + vp1y*segy
dp1w := wp1x*segx + wp1y*segy
if (dp1v < 0 && dp1w < 0) || (dp1v > 0 && dp1w > 0) {
t.Error("point not contained in segment ", dp1v, dp1w)
t.FailNow()
}
wp2x, wp2y := p2x-s2x, p2y-s2y
dp2v := vp2x*segx + vp2y*segy
dp2w := wp2x*segx + wp2y*segy
if (dp2v < 0 && dp2w < 0) || (dp2v > 0 && dp2w > 0) {
t.Error("point not contained in segment ", dp2v, dp2w)
t.FailNow()
}
if s1x == s2x && s2y == s1y { //Only one solution
// Test that one end of the segment is withing the radius of the circle
// and one is not
if seg1Inside && seg2Inside {
t.Error("Only one solution but both line segment ends inside")
t.FailNow()
}
if !seg1Inside && !seg2Inside {
t.Error("Only one solution but both line segment ends outside")
t.FailNow()
}
}
} else { // No intersection, check if both points outside or inside
if (seg1Inside && !seg2Inside) || (!seg1Inside && seg2Inside) {
t.Error("No solution but only one point in radius of circle")
t.FailNow()
}
}
}
}
t.Log("Tested ", testCntr, " examples and found ", solutionCntr, " solutions.")
}
以下是测试的输出结果:
=== RUN TestSegmentCircleIntersection
--- PASS: TestSegmentCircleIntersection (0.00s)
geom_test.go:105: Tested 40000 examples and found 7343 solutions.
虽然我认为使用线圆交点,然后检查交点是否在端点之间更好,而且可能更便宜,但我想添加这个更直观的解决方案。
我喜欢把问题看作是“香肠问题”,这应该适用于任何维度,而不改变算法。这种解决方案将无法找到交点。
以下是我的想法:
(我使用“小于”,但根据我们正在测试的内容,“小于或等于”也可以使用。)
图像:最粗的线段是所选线段,没有示例圆。有点粗糙,一些像素略微错位。
function boolean pointInSausage(sp1,sp2,r,c) {
if ( !(pointLineDist(c,sp1,sp2) < r) ) {
return false;
}
double a = dist(sp1,c);
double b = dist(sp2,c);
double l;
double s;
if (a>b) {
l = a;
s = b;
} else {
l = b;
s = a;
}
double segLength = dist(sp1,sp2);
if ( l < sqrt(segLength*segLength+r*r) ) {
return true;
}
return s < r;
}
如果发现任何问题,我会进行编辑或撤回。
#Line segment points (A0, Af) defined by xA0, yA0, xAf, yAf; circle center denoted by xB, yB; rB=radius of circle, rA = radius of point (set to zero for your application)
def staticCollision_f(xA0, yA0, xAf, yAf, rA, xB, yB, rB): #note potential speed up here by casting all variables to same type and/or using Cython
#Build equations of a line for linear agents (convert y = mx + b to ax + by + c = 0 means that a = -m, b = 1, c = -b
m_v = (yAf - yA0) / (xAf - xA0)
b_v = yAf - m_v * xAf
rEff = rA + rB #radii are added since we are considering the agent path as a thin line
#Check if points (circles) are within line segment (find center of line segment and check if circle is within radius of this point)
segmentMask = np.sqrt( (yB - (yA0+yAf)/2)**2 + (xB - (xA0+xAf)/2)**2 ) < np.sqrt( (yAf - yA0)**2 + (xAf - xA0)**2 ) / 2 + rEff
#Calculate perpendicular distance between line and a point
dist_v = np.abs(-m_v * xB + yB - b_v) / np.sqrt(m_v**2 + 1)
collisionMask = (dist_v < rEff) & segmentMask
#return True if collision is detected
return collisionMask, collisionMask.any()