python
def cross(a, b):
return a[0] * b[1] - a[1] * b[0]
def proper_intersection(seg1, seg2):
"""
返回两条线段的严格内部交点。
为突出扫描线算法机制，此函数只处理
不平行、没有端点接触的普通相交情况。
seg1, seg2: ((x1, y1), (x2, y2))
"""
a, b = seg1
c, d = seg2
r = (b[0] - a[0], b[1] - a[1])
s = (d[0] - c[0], d[1] - c[1])
q = (c[0] - a[0], c[1] - a[1])
denominator = cross(r, s)
if abs(denominator) < 1e-12:
return None
t = cross(q, s) / denominator
u = cross(q, r) / denominator
eps = 1e-10
if not (eps < t < 1-eps):
return None
if not (eps < u < 1-eps):
return None
return (
a[0] + t * r[0],
a[1] + t * r[1]
)
python
from heapq import heappush, heappop
from itertools import count
def cross(a, b):
return a[0] * b[1] - a[1] * b[0]
def proper_intersection(seg1, seg2):
a, b = seg1
c, d = seg2
r = (b[0]-a[0], b[1]-a[1])
s = (d[0]-c[0], d[1]-c[1])
q = (c[0]-a[0], c[1]-a[1])
den = cross(r, s)
if abs(den) < 1e-12:
return None
t = cross(q, s) / den
u = cross(q, r) / den
if not (1e-10 < t < 1-1e-10):
return None
if not (1e-10 < u < 1-1e-10):
return None
return (
a[0] + t*r[0],
a[1] + t*r[1]
)
def sweep_intersections(segments):
"""
教学版事件驱动扫描线算法。
输入：
{id: ((x1,y1), (x2,y2)), ...}
要求：
非竖直线段；
没有端点接触或共线重叠；
没有三线共点；
关键事件具有不同的x坐标。
使用普通列表维护活动状态，
不具备生产级Bentley-Ottmann时间界。
"""
S = {}
for name, (a, b) in segments.items():
if a[0] == b[0]:
raise ValueError("教学版不支持竖直线段")
S[name] = (
(a, b) if a[0] < b[0] else (b, a)
)
events = []
serial = count()
active = []
scheduled = set()
reported = []
def add_event(x, kind, a, b=None):
heappush(
events,
(x, kind, next(serial), a, b)
)
def y_at(name, x):
(x1, y1), (x2, y2) = S[name]
return (
y1 + (x-x1)*(y2-y1)/(x2-x1)
)
def arrange(a, b, current_x):
"""检查一对新相邻线段的未来交点。"""
if a is None or b is None:
return
key = tuple(sorted((a, b)))
if key in scheduled:
return
point = proper_intersection(
S[a], S[b]
)
if point is None:
return
if point[0] <= current_x + 1e-8:
return
scheduled.add(key)
add_event(point[0], 1, a, b)
def reorder(x):
active.sort(
key=lambda name: y_at(name, x+1e-6)
)
for name, (left, right) in S.items():
add_event(left[0], 0, name)
add_event(right[0], 2, name)
while events:
x, kind, _, a, b = heappop(events)
if kind == 0:
active.append(a)
reorder(x)
i = active.index(a)
if i > 0:
arrange(active[i-1], a, x)
if i+1 < len(active):
arrange(a, active[i+1], x)
elif kind == 2:
i = active.index(a)
lower = (
active[i-1] if i > 0 else None
)
upper = (
active[i+1]
if i+1 < len(active)
else None
)
active.pop(i)
arrange(lower, upper, x)
else:
i = active.index(a)
j = active.index(b)
if abs(i-j) != 1:
continue
point = proper_intersection(
S[a], S[b]
)
reported.append((a, b, point))
reorder(x)
for name in (a, b):
i = active.index(name)
if i > 0:
arrange(
active[i-1], name, x
)
if i+1 < len(active):
arrange(
name, active[i+1], x
)
return reported
segments = {
"A": ((0.8, 1.5), (9.0, 8.4)),
"B": ((1.2, 8.2), (9.5, 1.2)),
"C": ((2.0, 4.2), (8.4, 4.4)),
"D": ((3.2, 0.8), (7.8, 6.6)),
}
results = sweep_intersections(segments)
for a, b, point in results:
print(
a, b,
tuple(round(v, 3) for v in point)
)
print("总交点数量:", len(results))
python
from itertools import combinations
def brute_force(segments):
result = []
for (a, sa), (b, sb) in combinations(
segments.items(), 2
):
point = proper_intersection(sa, sb)
if point is not None:
result.append((a, b, point))
return result
brute = brute_force(segments)
actual_pairs = {
tuple(sorted((a, b)))
for a, b, _ in results
}
expected_pairs = {
tuple(sorted((a, b)))
for a, b, _ in brute
}
assert actual_pairs == expected_pairs
print("扫描线结果与穷举结果一致")
python
from shapely.geometry import LineString
from shapely.strtree import STRtree
segments = [
LineString([(0.8, 1.5), (9.0, 8.4)]),
LineString([(1.2, 8.2), (9.5, 1.2)]),
LineString([(2.0, 4.2), (8.4, 4.4)]),
LineString([(3.2, 0.8), (7.8, 6.6)]),
]
tree = STRtree(segments)
seen = set()
for i, segment in enumerate(segments):
candidates = tree.query(
segment,
predicate="intersects"
)
for j in candidates:
j = int(j)
if i >= j:
continue
pair = (i, j)
if pair in seen:
continue
seen.add(pair)
geometry = segment.intersection(segments[j])
print(
"线段", i, "与", j,
"相交，交集：",
geometry
)