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

Embed on website

To embed this program on your website, copy the following code and paste it into your website's HTML: