"""Line drawings -> one continuous theta-rho path, for logos and pictures.

The ball can't lift, so every stroke of a drawing has to be joined into one line. To get from
the end of one stroke to the start of the next, the path runs back along lines it has already
drawn (that leaves no new mark) and only crosses fresh sand where it must, as short as possible.

Input is an SVG file (paths, lines, polylines, polygons, rectangles, circles and ellipses; fills
are ignored, every shape is drawn as its outline) or a list of strokes. Draw the picture as you
want it in the sand: centre lines for thin parts, outlines for solid ones.
"""
import heapq
import math
import re
import xml.etree.ElementTree as ET

from .patterns import densify_xy, xy_to_polar

# ---------------------------------------------------------------- SVG reading
NUM = r'[-+]?(?:\d+\.?\d*|\.\d+)(?:[eE][-+]?\d+)?'
TOKEN = re.compile(r'[MmLlHhVvCcSsQqTtAaZz]|' + NUM)
SKIP = {'defs', 'clipPath', 'mask', 'symbol', 'marker', 'pattern', 'title', 'desc', 'metadata', 'style', 'text'}


def _mul(m, n):
    a, b, c, d, e, f = m
    a2, b2, c2, d2, e2, f2 = n
    return (a * a2 + c * b2, b * a2 + d * b2, a * c2 + c * d2, b * c2 + d * d2, a * e2 + c * f2 + e, b * e2 + d * f2 + f)


def _apply(m, p):
    a, b, c, d, e, f = m
    return (a * p[0] + c * p[1] + e, b * p[0] + d * p[1] + f)


def parse_transform(text):
    m = (1, 0, 0, 1, 0, 0)
    for name, args in re.findall(r'(\w+)\s*\(([^)]*)\)', text or ''):
        v = [float(x) for x in re.findall(NUM, args)]
        if name == 'translate':
            t = (1, 0, 0, 1, v[0], v[1] if len(v) > 1 else 0)
        elif name == 'scale':
            t = (v[0], 0, 0, v[1] if len(v) > 1 else v[0], 0, 0)
        elif name == 'rotate':
            r = math.radians(v[0])
            t = (math.cos(r), math.sin(r), -math.sin(r), math.cos(r), 0, 0)
            if len(v) == 3:
                t = _mul(_mul((1, 0, 0, 1, v[1], v[2]), t), (1, 0, 0, 1, -v[1], -v[2]))
        elif name == 'matrix' and len(v) == 6:
            t = tuple(v)
        elif name == 'skewX':
            t = (1, 0, math.tan(math.radians(v[0])), 1, 0, 0)
        elif name == 'skewY':
            t = (1, math.tan(math.radians(v[0])), 0, 1, 0, 0)
        else:
            continue
        m = _mul(m, t)
    return m


def _arc(p0, rx, ry, phi, large, sweep, p1, tol):
    """SVG elliptical arc as points (endpoint to centre form, SVG spec F.6)."""
    if rx == 0 or ry == 0 or p0 == p1:
        return [p1]
    rx, ry = abs(rx), abs(ry)
    cp, sp = math.cos(phi), math.sin(phi)
    dx, dy = (p0[0] - p1[0]) / 2, (p0[1] - p1[1]) / 2
    x1, y1 = cp * dx + sp * dy, -sp * dx + cp * dy
    lam = x1 * x1 / (rx * rx) + y1 * y1 / (ry * ry)
    if lam > 1:
        rx, ry = rx * math.sqrt(lam), ry * math.sqrt(lam)
    num = rx * rx * ry * ry - rx * rx * y1 * y1 - ry * ry * x1 * x1
    den = rx * rx * y1 * y1 + ry * ry * x1 * x1
    k = math.sqrt(max(0.0, num / den)) if den else 0.0
    if large == sweep:
        k = -k
    cx1, cy1 = k * rx * y1 / ry, -k * ry * x1 / rx
    cx, cy = cp * cx1 - sp * cy1 + (p0[0] + p1[0]) / 2, sp * cx1 + cp * cy1 + (p0[1] + p1[1]) / 2

    def ang(ux, uy, vx, vy):
        return math.atan2(ux * vy - uy * vx, ux * vx + uy * vy)
    t1 = ang(1, 0, (x1 - cx1) / rx, (y1 - cy1) / ry)
    dt = ang((x1 - cx1) / rx, (y1 - cy1) / ry, (-x1 - cx1) / rx, (-y1 - cy1) / ry)
    if not sweep and dt > 0:
        dt -= 2 * math.pi
    elif sweep and dt < 0:
        dt += 2 * math.pi
    n = max(4, int(math.ceil(abs(dt) * max(rx, ry) / tol)))
    out = []
    for i in range(1, n + 1):
        t = t1 + dt * i / n
        ex, ey = rx * math.cos(t), ry * math.sin(t)
        out.append((cp * ex - sp * ey + cx, sp * ex + cp * ey + cy))
    out[-1] = p1
    return out


def _bezier(pts, tol):
    L = sum(math.dist(a, b) for a, b in zip(pts, pts[1:]))
    n = max(2, min(400, int(math.ceil(L / tol))))
    out = []
    for i in range(1, n + 1):
        t = i / n
        q = list(pts)
        while len(q) > 1:
            q = [((1 - t) * a[0] + t * b[0], (1 - t) * a[1] + t * b[1]) for a, b in zip(q, q[1:])]
        out.append(q[0])
    return out


def path_polylines(d, tol=0.5):
    """An SVG path's d attribute -> list of polylines (one per subpath)."""
    toks = TOKEN.findall(d)
    i, cmd = 0, None
    cur = start = (0.0, 0.0)
    last_ctrl = None
    lines, poly = [], []

    def num():
        nonlocal i
        v = float(toks[i])
        i += 1
        return v

    def flag():
        # arc flags may be written packed together ("a5 5 0 105 5")
        nonlocal i
        t = toks[i]
        if t in ('0', '1'):
            i += 1
            return int(t)
        toks[i] = t[1:]
        return int(t[0])

    while i < len(toks):
        t = toks[i]
        if re.match(r'[A-Za-z]', t):
            cmd = t
            i += 1
            if cmd in 'Zz':
                if poly:
                    poly.append(start)
                    lines.append(poly)
                poly, cur, last_ctrl = [], start, None
                continue
        elif cmd is None:
            i += 1
            continue
        rel = cmd.islower()
        C = cmd.upper()
        ox, oy = cur if rel else (0.0, 0.0)
        if C == 'M':
            if len(poly) > 1:
                lines.append(poly)
            cur = (ox + num(), oy + num())
            start, poly, last_ctrl = cur, [cur], None
            cmd = 'l' if rel else 'L'          # further pairs are line-tos
            continue
        if not poly:
            poly = [cur]
        if C == 'L':
            cur = (ox + num(), oy + num())
            poly.append(cur)
            last_ctrl = None
        elif C == 'H':
            cur = ((cur[0] if rel else 0.0) + num(), cur[1])
            poly.append(cur)
            last_ctrl = None
        elif C == 'V':
            cur = (cur[0], (cur[1] if rel else 0.0) + num())
            poly.append(cur)
            last_ctrl = None
        elif C == 'C':
            c1 = (ox + num(), oy + num())
            c2 = (ox + num(), oy + num())
            p = (ox + num(), oy + num())
            poly += _bezier([cur, c1, c2, p], tol)
            cur, last_ctrl = p, ('C', c2)
        elif C == 'S':
            c1 = (2 * cur[0] - last_ctrl[1][0], 2 * cur[1] - last_ctrl[1][1]) if last_ctrl and last_ctrl[0] == 'C' else cur
            c2 = (ox + num(), oy + num())
            p = (ox + num(), oy + num())
            poly += _bezier([cur, c1, c2, p], tol)
            cur, last_ctrl = p, ('C', c2)
        elif C == 'Q':
            c1 = (ox + num(), oy + num())
            p = (ox + num(), oy + num())
            poly += _bezier([cur, c1, p], tol)
            cur, last_ctrl = p, ('Q', c1)
        elif C == 'T':
            c1 = (2 * cur[0] - last_ctrl[1][0], 2 * cur[1] - last_ctrl[1][1]) if last_ctrl and last_ctrl[0] == 'Q' else cur
            p = (ox + num(), oy + num())
            poly += _bezier([cur, c1, p], tol)
            cur, last_ctrl = p, ('Q', c1)
        elif C == 'A':
            rx, ry, rot = num(), num(), num()
            large, sweep = flag(), flag()
            p = (ox + num(), oy + num())
            poly += _arc(cur, rx, ry, math.radians(rot), large, sweep, p, tol)
            cur, last_ctrl = p, None
        else:
            i += 1
    if len(poly) > 1:
        lines.append(poly)
    return lines


def _ellipse(cx, cy, rx, ry, tol):
    n = max(12, int(math.ceil(2 * math.pi * max(rx, ry) / tol)))
    return [(cx + rx * math.cos(2 * math.pi * i / n), cy + ry * math.sin(2 * math.pi * i / n)) for i in range(n + 1)]


def svg_strokes(text):
    """Every shape in an SVG document as polylines, in the document's coordinates (y down)."""
    root = ET.fromstring(text)
    vb = [float(v) for v in re.findall(NUM, root.get('viewBox', ''))]
    size = max(vb[2], vb[3]) if len(vb) == 4 else max(
        float((re.findall(NUM, root.get('width', '')) or ['200'])[0]),
        float((re.findall(NUM, root.get('height', '')) or ['200'])[0]))
    tol = size / 500.0
    out = []

    def f(el, k, default=0.0):
        v = re.findall(NUM, el.get(k, ''))
        return float(v[0]) if v else default

    def walk(el, m):
        tag = el.tag.split('}')[-1]
        if tag in SKIP:
            return
        m = _mul(m, parse_transform(el.get('transform')))
        polys = []
        if tag == 'path':
            polys = path_polylines(el.get('d', ''), tol)
        elif tag == 'line':
            polys = [[(f(el, 'x1'), f(el, 'y1')), (f(el, 'x2'), f(el, 'y2'))]]
        elif tag in ('polyline', 'polygon'):
            v = [float(x) for x in re.findall(NUM, el.get('points', ''))]
            p = list(zip(v[0::2], v[1::2]))
            if tag == 'polygon' and p:
                p.append(p[0])
            polys = [p]
        elif tag == 'rect':
            x, y, w, h = f(el, 'x'), f(el, 'y'), f(el, 'width'), f(el, 'height')
            polys = [[(x, y), (x + w, y), (x + w, y + h), (x, y + h), (x, y)]]
        elif tag == 'circle':
            r = f(el, 'r')
            polys = [_ellipse(f(el, 'cx'), f(el, 'cy'), r, r, tol)]
        elif tag == 'ellipse':
            polys = [_ellipse(f(el, 'cx'), f(el, 'cy'), f(el, 'rx'), f(el, 'ry'), tol)]
        for p in polys:
            if len(p) > 1:
                out.append([_apply(m, q) for q in p])
        for ch in el:
            walk(ch, m)

    walk(root, (1, 0, 0, 1, 0, 0))
    return out


# ---------------------------------------------------------------- joining strokes into one line
def normalise(strokes, fit=0.9, flip_y=True):
    """Centre the drawing and scale it to fit inside a circle of radius `fit` (the table is 1)."""
    pts = [p for s in strokes for p in s]
    x0, x1 = min(p[0] for p in pts), max(p[0] for p in pts)
    y0, y1 = min(p[1] for p in pts), max(p[1] for p in pts)
    cx, cy = (x0 + x1) / 2, (y0 + y1) / 2
    R = max(math.hypot(p[0] - cx, p[1] - cy) for p in pts) or 1.0
    k = fit / R
    sy = -k if flip_y else k
    return [[((x - cx) * k, (y - cy) * sy) for x, y in s] for s in strokes]


class _Graph:
    """The lines drawn so far: points about `snap` apart, joined where the drawing joins or crosses."""

    def __init__(self, snap):
        self.snap = snap
        self.pos, self.adj, self.cells = [], [], {}

    def node(self, p):
        cx, cy = int(math.floor(p[0] / self.snap)), int(math.floor(p[1] / self.snap))
        near = []
        for gx in (cx - 1, cx, cx + 1):
            for gy in (cy - 1, cy, cy + 1):
                for j in self.cells.get((gx, gy), ()):
                    dd = math.dist(self.pos[j], p)
                    if dd < self.snap * 0.5:
                        return j                      # the same point
                    if dd < self.snap * 1.5:
                        near.append((j, dd))
        i = len(self.pos)
        self.pos.append(p)
        self.adj.append([])
        self.cells.setdefault((cx, cy), []).append(i)
        for j, dd in near:                            # lines that meet or cross here
            self.adj[i].append((j, dd))
            self.adj[j].append((i, dd))
        return i

    def add_line(self, pts):
        prev = None
        for p in pts:
            i = self.node(p)
            if prev is not None and prev != i:
                w = math.dist(self.pos[prev], self.pos[i])
                self.adj[prev].append((i, w))
                self.adj[i].append((prev, w))
            prev = i
        return prev

    def dijkstra(self, src):
        dist, prev = {src: 0.0}, {}
        heap = [(0.0, src)]
        while heap:
            d, u = heapq.heappop(heap)
            if d > dist[u]:
                continue
            for v, w in self.adj[u]:
                nd = d + w
                if nd < dist.get(v, 1e18):
                    dist[v], prev[v] = nd, u
                    heapq.heappush(heap, (nd, v))
        return dist, prev


def join_strokes(strokes, jump_cost=4.0, step=0.005, snap=0.006, stats=None):
    """Order the strokes and join them into one line (unit-disc coordinates). It starts and ends
    at the centre or the edge, whichever is closer to the drawing (a pattern must, so the erase
    spiral can hand over to it). Going back over drawn lines costs its length; crossing fresh
    sand costs `jump_cost` times its length."""
    strokes = [densify_xy(s, step) for s in strokes if len(s) > 1 and sum(math.dist(a, b) for a, b in zip(s, s[1:])) > step]
    g = _Graph(snap)
    cur_node = None
    drawn = retrace = jump = 0.0
    todo = list(range(len(strokes)))

    def entries(k):
        s = strokes[k]
        if math.dist(s[0], s[-1]) < 2 * snap:         # closed: start anywhere on it
            every = max(1, int(0.03 / step))
            return [(k, j, 'loop') for j in range(0, len(s) - 1, every)]
        return [(k, 0, 'fwd'), (k, len(s) - 1, 'rev')]

    # start: the centre, or the edge next to the stroke end that is nearest to it
    ends = [strokes[k][j] for k in todo for (_, j, _) in entries(k)]
    p_in = min(ends, key=lambda p: math.hypot(*p))
    p_out = max(ends, key=lambda p: math.hypot(*p))
    if 1 - math.hypot(*p_out) < 2 * math.hypot(*p_in):      # a short line at the edge shows less than one in the middle
        r = math.hypot(*p_out)
        start = (p_out[0] / r, p_out[1] / r)
    else:
        start = (0.0, 0.0)
    path = [start]

    def best_route(targets):
        """Cheapest way from where the ball is to one of the target points."""
        if cur_node is None:
            dist, prevs = {}, {}
        else:
            dist, prevs = g.dijkstra(cur_node)
        here = path[-1]
        cell = 0.05
        grid = {}
        for v in dist:
            p = g.pos[v]
            grid.setdefault((int(math.floor(p[0] / cell)), int(math.floor(p[1] / cell))), []).append(v)
        best = None
        for t, p in targets:
            c0 = jump_cost * math.dist(here, p)
            cand = (c0, None, t, p)
            r = c0 / jump_cost
            if dist:
                gx0, gx1 = int(math.floor((p[0] - r) / cell)), int(math.floor((p[0] + r) / cell))
                gy0, gy1 = int(math.floor((p[1] - r) / cell)), int(math.floor((p[1] + r) / cell))
                if (gx1 - gx0 + 1) * (gy1 - gy0 + 1) > len(grid):
                    cells = list(grid.values())
                else:
                    cells = [grid[(gx, gy)] for gx in range(gx0, gx1 + 1) for gy in range(gy0, gy1 + 1) if (gx, gy) in grid]
                for vs in cells:
                    for v in vs:
                        c = dist[v] + jump_cost * math.dist(g.pos[v], p)
                        if c < cand[0]:
                            cand = (c, v, t, p)
            if best is None or cand[0] < best[0]:
                best = cand
        _, v, t, p = best
        back = []
        if v is not None:
            while v != cur_node:
                back.append(g.pos[v])
                v = prevs[v]
            back.reverse()
        return t, back, p

    while todo:
        targets = []
        for k in todo:
            for e in entries(k):
                targets.append((e, strokes[e[0]][e[1]]))
        (k, j, how), back, p = best_route(targets)
        if back:
            retrace += sum(math.dist(a, b) for a, b in zip([path[-1]] + back, back))
            path += back
        jump += math.dist(path[-1], p)
        path += densify_xy([path[-1], p], step)[1:]
        s = strokes[k]
        if how == 'loop':
            seq = s[j:] + s[1:j + 1]
        elif how == 'rev':
            seq = s[::-1]
        else:
            seq = s
        drawn += sum(math.dist(a, b) for a, b in zip(seq, seq[1:]))
        path += seq[1:]
        cur_node = g.add_line(seq)
        todo.remove(k)
    # finish at the centre or the edge, whichever is cheaper from somewhere already drawn
    dist, prevs = g.dijkstra(cur_node)
    best = (jump_cost * min(2 * math.hypot(*path[-1]), 1 - math.hypot(*path[-1])), cur_node)
    for v, d in dist.items():
        c = d + jump_cost * min(2 * math.hypot(*g.pos[v]), 1 - math.hypot(*g.pos[v]))
        if c < best[0]:
            best = (c, v)
    v = best[1]
    back = []
    while v != cur_node:
        back.append(g.pos[v])
        v = prevs[v]
    back.reverse()
    q = back[-1] if back else path[-1]
    rq = math.hypot(*q)
    end = (0.0, 0.0) if 2 * rq < 1 - rq else (q[0] / rq, q[1] / rq)
    if back:
        retrace += sum(math.dist(a, b) for a, b in zip([path[-1]] + back, back))
        path += back
    jump += math.dist(path[-1], end)
    path += densify_xy([path[-1], end], step)[1:]
    if stats is not None:
        stats.update(start=start, end=end)
        stats.update(drawn=drawn, retrace=retrace, jump=jump)
    return path


def drawing_to_polar(strokes, fit=0.9, border=False, stats=None, **kw):
    """Strokes (any units, y down as in SVG) -> a theta-rho pattern that starts and ends at the centre or edge.
    border: also draw a circle round the picture, just inside the edge."""
    s = normalise(strokes, fit)
    if border:
        n = 720
        s.append([(0.97 * math.cos(2 * math.pi * i / n), 0.97 * math.sin(2 * math.pi * i / n)) for i in range(n + 1)])
    xy = join_strokes(s, stats=stats, **kw)
    return xy_to_polar(densify_xy(xy, 0.01))


def svg_to_polar(text, **kw):
    return drawing_to_polar(svg_strokes(text), **kw)