from math import hypot, acos, cos, sin, pi, atan2
from collections import defaultdict
import heapq as hq
import time
EPS = 1e-7
# returns True if the point corresponding to the angle x lies between the points for angles a and b
def is_between_angles(x, a, b):
a %= 2 * pi
b %= 2 * pi
x %= 2 * pi
delta = (b - a) % (2 * pi)
if delta <= pi:
return (x - a) % (2 * pi) <= delta
else:
return (a - x) % (2 * pi) <= (2 * pi - delta)
def is_between_angles_biggest(x, a, b):
a %= 2 * pi
b %= 2 * pi
x %= 2 * pi
delta = (b - a) % (2 * pi)
if delta <= pi:
# smallest arc is from a to b (counterclockwise)
return (x - a) % (2 * pi) > delta
else:
# smallest arc is from b to a (counterclockwise)
return (a - x) % (2 * pi) > (2 * pi - delta)
# Returns the angle between the x-axis and vector(center, pt)
def angle(center, pt):
return atan2(pt.y - center.y, pt.x - center.x)
class _Point:
def __init__(self, x, y):
self.x = x
self.y = y
def __str__(self):
return f"Point({{{self.x}, {self.y}}})"
return f"({self.x}, {self.y})"
def __add__(self, other):
return _Point(self.x + other.x, self.y + other.y)
def __sub__(self, other):
return _Point(self.x - other.x, self.y - other.y)
def __neg__(self):
return _Point(-self.x, -self.y)
def __rmul__(self, k):
return _Point(k * self.x, k * self.y)
def __eq__(self, other):
return self.x == other.x and self.y == other.y
def __lt__(self, other):
return (self.x, self.y) <= (other.x, other.y)
class Segment:
def __init__(self, p1, p2):
self.p1 = p1
self.p2 = p2
def length(self):
return hypot(self.p1.x - self.p2.x, self.p1.y - self.p2.y)
def __repr__(self):
return f"Segment {self.p1} --- {self.p2}"
# returns True if the segment intersect a Circle
def intersect(self, circle):
u = (
(self.p2.x - self.p1.x) * (circle.center.x - self.p1.x)
+ (self.p2.y - self.p1.y) * (circle.center.y - self.p1.y)
) / self.length() ** 2
# pt: the shortest point of segment [p1, p2] to the center of the circle
pt = self.p1 + max(0, min(u, 1)) * (self.p2 - self.p1)
seg = Segment(circle.center, pt)
# Some threshold to prevent approximation errors
r, l = circle.radius, seg.length()
return l + EPS < r
class _Circle:
def __init__(self, center, radius):
self.center = center
self.radius = radius
def __repr__(self):
return f"Circle(Point({{{self.center.x}, {self.center.y}}}), {self.radius})"
return f"_Circle of center {self.center} and radius {self.radius}"
# finds the two points (if they exist) of the circle that have tangents passing through point p
def tangents_to_point(self, p):
c, r = self.center, self.radius
d = Segment(c, p).length()
if r > d:
return (False, [])
u = max(-1, min(r / d, 1))
t = acos(u)
q = c + r / d * (p - c)
p1 = _Point(
c.x + cos(t) * (q.x - c.x) - sin(t) * (q.y - c.y),
c.y + sin(t) * (q.x - c.x) + cos(t) * (q.y - c.y),
)
p2 = _Point(
c.x + cos(t) * (q.x - c.x) + sin(t) * (q.y - c.y),
c.y - sin(t) * (q.x - c.x) + cos(t) * (q.y - c.y),
)
return (True, [p1, p2])
# Given two circles c1 and c2, finds (if they exist) the points p1, p2 in c1, q1 and q2 in c2 such that (p1, q2) and (p2, q1) are internal bi tangents
def internal_bitangents(self, other):
r1, r2 = self.radius, other.radius
x1, y1 = self.center.x, self.center.y
x2, y2 = other.center.x, other.center.y
dis = hypot(x2 - x1, y2 - y1)
if r1 + r2 >= dis + EPS:
return (False, [])
u = max(-1, min((r1 + r2) / dis, 1))
t = acos(u)
p = self.center + r1 / dis * (other.center - self.center)
p1 = _Point(
x1 + cos(t) * (p.x - x1) - sin(t) * (p.y - y1),
y1 + sin(t) * (p.x - x1) + cos(t) * (p.y - y1),
)
p2 = _Point(
x1 + cos(t) * (p.x - x1) + sin(t) * (p.y - y1),
y1 - sin(t) * (p.x - x1) + cos(t) * (p.y - y1),
)
#########
q = other.center + r2 / dis * (self.center - other.center)
q1 = _Point(
x2 + cos(t) * (q.x - x2) - sin(t) * (q.y - y2),
y2 + sin(t) * (q.x - x2) + cos(t) * (q.y - y2),
)
q2 = _Point(
x2 + cos(t) * (q.x - x2) + sin(t) * (q.y - y2),
y2 - sin(t) * (q.x - x2) + cos(t) * (q.y - y2),
)
return (True, [p1, p2, q1, q2])
# Given two circles c1 and c2, finds (if they exist) the points p1, p2 in c1, q1 and q2 in c2 such that (p1, q2) and (p2, q1) are external bi tangents
def external_bitangents(self, other):
r1, r2 = self.radius, other.radius
x1, y1 = self.center.x, self.center.y
x2, y2 = other.center.x, other.center.y
dis = hypot(x2 - x1, y2 - y1)
if abs(r1 - r2) > dis:
return (False, [])
u = max(-1, min((r1 - r2) / dis, 1))
t = acos(u)
p = self.center + r1 / dis * (other.center - self.center)
p1 = _Point(
x1 + cos(t) * (p.x - x1) - sin(t) * (p.y - y1),
y1 + sin(t) * (p.x - x1) + cos(t) * (p.y - y1),
)
p2 = _Point(
x1 + cos(t) * (p.x - x1) + sin(t) * (p.y - y1),
y1 - sin(t) * (p.x - x1) + cos(t) * (p.y - y1),
)
#########
q = self.center + (dis + r2) / dis * (other.center - self.center)
q1 = _Point(
x2 + cos(t) * (q.x - x2) - sin(t) * (q.y - y2),
y2 + sin(t) * (q.x - x2) + cos(t) * (q.y - y2),
)
q2 = _Point(
x2 + cos(t) * (q.x - x2) + sin(t) * (q.y - y2),
y2 - sin(t) * (q.x - x2) + cos(t) * (q.y - y2),
)
return (True, [p1, p2, q1, q2])
# gets the intersection with another circle
def intersect(self, other):
assert type(other) == _Circle, "Element is not a circle"
ca, cb = self.center, other.center
ra, rb = self.radius, other.radius
d = Segment(ca, cb).length()
# No intersection
if ra + rb < d or abs(ra - rb) > d:
return (False, [])
a = (ra**2 + d**2 - rb**2) / (2 * d)
# 1 intersection point
if abs(a - ra) <= EPS:
c = ca + a / d * (cb - ca)
return (False, [c])
# 2 intersection points
d = ca + ra / d * (cb - ca)
u = max(-1, min(a / ra, 1))
t = acos(u)
p = _Point(
ca.x + cos(t) * (d.x - ca.x) - sin(t) * (d.y - ca.y),
ca.y + sin(t) * (d.x - ca.x) + cos(t) * (d.y - ca.y),
)
q = _Point(
ca.x + cos(t) * (d.x - ca.x) + sin(t) * (d.y - ca.y),
ca.y - sin(t) * (d.x - ca.x) + cos(t) * (d.y - ca.y),
)
return (True, [p, q])
class PathFinder:
def __init__(self, pa, pb, circles):
self.start = pa
self.end = pb
self.circles = circles
self.elements = [self.start, self.end]
self.adj = defaultdict(list)
t1 = time.time()
n = len(circles)
# Direct path
seg = Segment(self.start, self.end)
if all(not seg.intersect(circle) for circle in self.circles):
l = seg.length()
self.adj[0].append((1, l))
self.adj[1].append((0, l))
num = 2
cs = defaultdict(list)
for i in range(n):
ci = self.circles[i]
for j in range(i + 1, n):
cj = self.circles[j]
b, ps = ci.internal_bitangents(cj)
self.elements.extend(ps if b else [None] * 4)
if b:
s1 = Segment(ps[0], ps[2])
s2 = Segment(ps[1], ps[3])
if not self.segment_intersect_circles(s1, i, j):
l = s1.length()
self.adj[num].append((num + 2, l))
self.adj[num + 2].append((num, l))
cs[i].append(num)
cs[j].append(num + 2)
if not self.segment_intersect_circles(s2, i, j):
l = s2.length()
self.adj[num + 1].append((num + 3, l))
self.adj[num + 3].append((num + 1, l))
cs[i].append(num + 1)
cs[j].append(num + 3)
num += 4
########################
c, qs = ci.external_bitangents(cj)
self.elements.extend(qs if c else [None] * 4)
if c:
s1 = Segment(qs[0], qs[2])
s2 = Segment(qs[1], qs[3])
_c1 = all(not s1.intersect(circle) for circle in self.circles)
# print(i, j, _c1, [s1.intersect(circle) for k, circle in enumerate(self.circles) if not k in (i, j)], s1)
if _c1:
l = s1.length()
self.adj[num].append((num + 2, l))
self.adj[num + 2].append((num, l))
cs[i].append(num)
cs[j].append(num + 2)
_c2 = all(not s2.intersect(circle) for circle in self.circles)
# print(i, j, _c2, [s2.intersect(circle) for k, circle in enumerate(self.circles) if not k in (i, j)], s2)
if _c2:
l = s2.length()
self.adj[num + 1].append((num + 3, l))
self.adj[num + 3].append((num + 1, l))
cs[i].append(num + 1)
cs[j].append(num + 3)
num += 4
t2 = time.time()
print("Phase 1", t2 - t1)
for i, c in enumerate(circles):
b, ps = c.tangents_to_point(pa)
if b:
for q in ps:
s = Segment(pa, q)
if all(not s.intersect(circle) for k, circle in enumerate(circles) if k != i):
l = s.length()
self.adj[0].append((num, l))
self.adj[num].append((0, l))
cs[i].append(num)
self.elements.append(q)
num += 1
b, ps = c.tangents_to_point(pb)
if b:
for q in ps:
s = Segment(pb, q)
if all(not s.intersect(circle) for k, circle in enumerate(circles) if k != i):
l = s.length()
self.adj[1].append((num, l))
self.adj[num].append((1, l))
cs[i].append(num)
self.elements.append(q)
num += 1
t3 = time.time()
print("Phase 2", t3 - t2)
for i, _s in cs.items():
k = len(_s)
for a in range(k):
for b in range(a + 1, k):
p, q = self.elements[_s[a]], self.elements[_s[b]]
m, M = PathFinder.arc_distance(p, q, circles[i])
if self.arc_is_good(p, q, i):
self.adj[_s[a]].append((_s[b], m))
self.adj[_s[b]].append((_s[a], m))
if self.arc_is_good_big(p, q, i):
self.adj[_s[a]].append((_s[b], M))
self.adj[_s[b]].append((_s[a], M))
t4 = time.time()
print("Phase 3 :", t4 - t3)
def __repr__(self):
sep = "\n" + "#" * 20 + "\n"
arr = []
for i, (p, t) in enumerate(self.elements):
arr.append(f"_Point{i}:{p} of type {t}")
n = f"Number of elements : {len(self.elements)}"
m = "\n".join(
"{1} ({2}): {3}".format(key, self.elements[key], value)
for key, value in self.adj
)
return (
f"Tangent Visibility Graph\n{n}\n{sep.join(arr)}\nAdjacency matrix :\n{m}"
)
# checks if the smallest arc p, q lying on circle of index circle_idx intersect any of the other circles
def arc_is_good(self, p, q, circle_idx):
circle = self.circles[circle_idx]
for k, c in enumerate(self.circles):
if k == circle_idx:
continue
b, pts = circle.intersect(c)
if not b:
continue
center = circle.center
t1 = angle(center, p)
t2 = angle(center, q)
# if any of the intersection lies between t1 and t2, the segment is not valid
if any(is_between_angles(angle(center, p), t1, t2) for p in pts):
return False
return True
# checks if the biggest arc p, q lying on circle of index circle_idx intersect any of the other circles
def arc_is_good_big(self, p, q, circle_idx):
circle = self.circles[circle_idx]
for k, c in enumerate(self.circles):
if k == circle_idx:
continue
b, pts = circle.intersect(c)
if not b:
continue
center = circle.center
t1 = angle(center, p)
t2 = angle(center, q)
# if any of the intersection lies between t1 and t2, the segment is not valid
if any(is_between_angles_biggest(angle(center, p), t1, t2) for p in pts):
return False
return True
def segment_intersect_circles(self, s, i, j):
for k in range(len(self.circles)):
if k in (i, j):
continue
if s.intersect(self.circles[k]):
return True
return False
def tangent_intersect_circles(self, s, i):
for k in range(len(self.circles)):
if k == i:
continue
if s.intersect(self.circles[k]):
return True
return False
# finds the shortest distance between 2 points on a circle
@staticmethod
def arc_distance(p, q, c):
center, radius = (
c.center,
c.radius,
)
u = (
(p.x - center.x) * (q.x - center.x) + (p.y - center.y) * (q.y - center.y)
) / radius**2
# clamp to 0..1 to prevent math domain errors
u = max(-1, min(u, 1))
t = acos(u)
l1, l2 = radius * t, radius * (2 * pi - t)
return min(l1, l2), max(l1, l2)
# heuristic for the a star implementation
@staticmethod
def heuristic(p, q):
return hypot(p.x - q.x, p.y - q.y)
# classic a star implementation to find the shrtest path from start to end (elements 0 and 1 of self.elements)
def a_star_search(self):
t0 = time.time()
N = len(self.elements)
G = {
i: float("inf") for i in range(N)
} # Actual movement cost to each position from the start position
F = {
i: float("inf") for i in range(N)
} # Estimated movement cost of start to end going via this position
start, end = 0, 1
# Initialize starting values
G[start] = 0
p_start = self.elements[start]
p_end = self.elements[end]
d_init = PathFinder.heuristic(p_start, p_end)
F[start] = d_init
open_vertices = [(d_init, start)]
seen = set([start])
hq.heapify(open_vertices)
came_from = {}
while len(open_vertices) > 0:
_, current = hq.heappop(open_vertices)
# Check if we have reached the goal
if current == end:
# Retrace our route backward
# print("path:", came_from)
path = [current]
while current in came_from:
current = came_from[current]
path.append(current)
path.reverse()
t1 = time.time()
print("Phase 4 (A star) :", t1 - t0)
return path, F[end]
# Update scores for vertices near the current position
for neighbor, dis in self.adj[current]:
p_end = self.elements[end]
tentative_gscore = G[current] + dis
if tentative_gscore < G[neighbor]:
came_from[neighbor] = current
G[neighbor] = tentative_gscore
p_neighbor = self.elements[neighbor]
F[neighbor] = tentative_gscore + PathFinder.heuristic(
p_neighbor, p_end
)
if not neighbor in seen:
seen.add(neighbor)
hq.heappush(open_vertices, (F[neighbor], neighbor))
return [], -1
# Random Tests 1 - seeded (0 of 5 Assertions)
# Random Tests 1 - not seeded (0 of 15 Assertions)
# Random Tests 2 - seeded (0 of 3 Assertions)
# Random Tests 2 - not seeded (0 of 7 Assertions)
def shortest_path_length(a, b, cs):
print(a, b, cs)
pa = _Point(a.x, a.y)
pb = _Point(b.x, b.y)
circles = [_Circle(_Point(c.ctr.x, c.ctr.y), c.r) for c in cs]
t0 = time.time()
if len(circles) < 20:
print(pa, pb, '\n'.join(str(c) for c in circles))
pf = PathFinder(pa, pb, circles)
# print('\n'.join(str(p) for p in pf.elements if p is not None))
try:
path, dis = pf.a_star_search()
t1 = time.time()
print("A star (end) :", t1 - t0)
return dis
except:
return -1
return 0
To embed this program on your website, copy the following code and paste it into your website's HTML: