Draw compound medial tires as separated cycles

papers/medial_tire_decompositions_of_plane_triangulations/lib/draw_random_medial_tire_decompositions.py

@@ -30,6 +30,8 @@ import matplotlib
 matplotlib.use("Agg")
 import matplotlib.pyplot as plt
 from matplotlib.lines import Line2D
+from matplotlib.patches import PathPatch
+from matplotlib.path import Path as MplPath
 import networkx as nx
 
 if str(PAPER_DIR) not in sys.path:
@@ -368,13 +370,13 @@ def draw_tire_tree(ax, nodes: list[TreadNode], tree_edges):
     ax.axis("off")
 
 
-def vertex_xy(k: int, n: int, radius: float) -> tuple[float, float]:
+def vertex_xy(k: int, n: int, radius: float, rotation: float = 0.0) -> tuple[float, float]:
-    angle = math.pi / 2 - 2 * math.pi * k / n
+    angle = math.pi / 2 + rotation - 2 * math.pi * k / n
     return radius * math.cos(angle), radius * math.sin(angle)
 
 
-def edge_midpoint_angle(i: int, n: int) -> float:
+def edge_midpoint_angle(i: int, n: int, rotation: float = 0.0) -> float:
-    return math.pi / 2 - 2 * math.pi * (i + 0.5) / n
+    return math.pi / 2 + rotation - 2 * math.pi * (i + 0.5) / n
 
 
 def annular_cycle_edges(node: TreadNode) -> set[tuple]:
@@ -386,99 +388,99 @@ def annular_cycle_edges(node: TreadNode) -> set[tuple]:
     return edges
 
 
-def draw_compound_tread_model(ax, node: TreadNode):
+def shared_up_apex_occurrences(node: TreadNode) -> dict[tuple, list[tuple[int, int]]]:
-    """Draw a compound tread using a planar layout of its actual medial graph."""
+    occurrences: dict[tuple, list[tuple[int, int]]] = defaultdict(list)
-    try:
-        pos = nx.planar_layout(node.medial)
-    except nx.NetworkXException:
-        pos = nx.spring_layout(node.medial, seed=node.idx)
-
-    cycle_edges = annular_cycle_edges(node)
-    non_cycle_edges = [
-        edge for edge in node.medial.edges()
-        if tuple(sorted(edge)) not in cycle_edges
-    ]
-
-    nx.draw_networkx_edges(
-        node.medial,
-        pos,
-        edgelist=non_cycle_edges,
-        ax=ax,
-        edge_color="#cbd5e1",
-        width=0.7,
-    )
-    nx.draw_networkx_edges(
-        node.medial,
-        pos,
-        edgelist=list(cycle_edges),
-        ax=ax,
-        edge_color="black",
-        width=1.4,
-    )
-
     annular = set(node.annular)
-    singleton_down = set(node.down) - set(node.bites)
+    for cycle_idx, order in enumerate(node.annular_cycles):
-    categories = [
+        n = len(order)
-        (annular, "black", 13, "none"),
+        for i, a in enumerate(order):
-        (set(node.up) - annular, "#2563eb", 18, "none"),
+            b = order[(i + 1) % n]
-        (singleton_down - annular, "#dc2626", 18, "none"),
+            apexes = [
-        (set(node.bites) - annular, "#7f1d1d", 28, "black"),
+                w for w in set(node.medial.neighbors(a)) & set(node.medial.neighbors(b))
-    ]
+                if w not in annular
-    for vertices, color, size, edgecolor in categories:
+            ]
-        drawn = [v for v in vertices if v in pos]
+            for apex in apexes:
-        if not drawn:
+                if apex in node.up:
-            continue
+                    occurrences[apex].append((cycle_idx, i))
-        ax.scatter(
+    return {apex: where for apex, where in occurrences.items() if len(where) > 1}
-            [pos[v][0] for v in drawn],
+
-            [pos[v][1] for v in drawn],
+
-            s=size,
+def compound_cycle_rotations(node: TreadNode) -> list[float]:
-            color=color,
+    rotations = [0.0] * len(node.annular_cycles)
-            edgecolors=edgecolor,
+    shared = shared_up_apex_occurrences(node)
-            linewidths=0.4,
+    by_cycle: dict[int, list[int]] = defaultdict(list)
-            zorder=3,
+    for occurrences in shared.values():
+        for cycle_idx, edge_idx in occurrences:
+            by_cycle[cycle_idx].append(edge_idx)
+
+    for cycle_idx, edge_indices in by_cycle.items():
+        n = len(node.annular_cycles[cycle_idx])
+        sx = sy = 0.0
+        for edge_idx in edge_indices:
+            angle = edge_midpoint_angle(edge_idx, n)
+            sx += math.cos(angle)
+            sy += math.sin(angle)
+        if sx or sy:
+            rotations[cycle_idx] = math.pi / 2 - math.atan2(sy, sx)
+    return rotations
+
+
+def compound_cycle_order(node: TreadNode) -> list[int]:
+    shared = shared_up_apex_occurrences(node)
+    adj: dict[int, set[int]] = defaultdict(set)
+    for occurrences in shared.values():
+        cycles = sorted({cycle_idx for cycle_idx, _edge_idx in occurrences})
+        for a, b in zip(cycles, cycles[1:]):
+            adj[a].add(b)
+            adj[b].add(a)
+
+    remaining = set(range(len(node.annular_cycles)))
+    order = []
+    while remaining:
+        candidates = sorted(
+            remaining,
+            key=lambda idx: (len(adj[idx]) if adj[idx] else 999, idx),
         )
-
+        start = candidates[0]
-    xs = [p[0] for p in pos.values()]
+        stack = [(start, None)]
-    ys = [p[1] for p in pos.values()]
+        while stack:
-    xpad = max(0.05, (max(xs) - min(xs)) * 0.12)
+            current, parent = stack.pop()
-    ypad = max(0.05, (max(ys) - min(ys)) * 0.12)
+            if current not in remaining:
-    ax.set_xlim(min(xs) - xpad, max(xs) + xpad)
+                continue
-    ax.set_ylim(min(ys) - ypad, max(ys) + ypad)
+            remaining.remove(current)
-    ax.set_aspect("equal")
+            order.append(current)
-    ax.axis("off")
+            children = sorted(
+                (neighbor for neighbor in adj[current] if neighbor != parent),
+                key=lambda idx: (len(adj[idx]), idx),
+                reverse=True,
+            )
+            for child in children:
+                stack.append((child, current))
+    return order
 
 
-def draw_tread_model(ax, node: TreadNode):
+def cycle_layout(
-    if len(node.annular_cycles) > 1:
+    node: TreadNode,
-        draw_compound_tread_model(ax, node)
+    rotations: list[float],
-        singleton_down = set(node.down) - set(node.bites)
+    display_order: list[int] | None = None,
-        ax.set_title(
+    cycle_spacing: float = 3.25,
-            f"T{node.idx} d={node.depth}: {len(node.annular_cycles)} annular cycle(s)\n"
+):
-            f"ann={len(node.annular)} up={len(node.up)} down={len(singleton_down)} "
-            f"bite={len(node.bites)}",
-            fontsize=6.4,
-            pad=1.5,
-        )
-        return
-
     cycle_count = len(node.annular_cycles)
-    offsets = [3.25 * (i - (cycle_count - 1) / 2) for i in range(cycle_count)]
+    if display_order is None:
-    apex_positions: dict[tuple, tuple[float, float]] = {}
+        display_order = list(range(cycle_count))
-    apex_corners: dict[tuple, list[tuple[float, float]]] = defaultdict(list)
+    rank = {cycle_idx: i for i, cycle_idx in enumerate(display_order)}
-    ann_positions: dict[tuple, tuple[float, float]] = {}
+    offsets = [cycle_spacing * (rank[i] - (cycle_count - 1) / 2) for i in range(cycle_count)]
+    ann_positions: dict[tuple[int, tuple], tuple[float, float]] = {}
+    apex_positions: dict[tuple[int, tuple], tuple[float, float]] = {}
+    apex_corners: dict[tuple[int, tuple], list[tuple[float, float]]] = defaultdict(list)
 
     for cycle_idx, order in enumerate(node.annular_cycles):
         n = len(order)
         dx = offsets[cycle_idx]
-        ann = {
+        rotation = rotations[cycle_idx]
-            vertex: (dx + x, y)
+        for k, vertex in enumerate(order):
-            for vertex, (x, y) in zip(order, [vertex_xy(k, n, 1.0) for k in range(n)])
+            x, y = vertex_xy(k, n, 1.0, rotation)
-        }
+            ann_positions[(cycle_idx, vertex)] = (dx + x, y)
-        ann_positions.update(ann)
-
-        cyc_x = [ann[v][0] for v in order] + [ann[order[0]][0]]
-        cyc_y = [ann[v][1] for v in order] + [ann[order[0]][1]]
-        ax.plot(cyc_x, cyc_y, color="black", lw=1.3, zorder=2)
 
         for i, a in enumerate(order):
             b = order[(i + 1) % n]
@@ -487,33 +489,157 @@ def draw_tread_model(ax, node: TreadNode):
                 if w not in node.annular
             ]
             for apex in apexes:
-                apex_corners[apex].extend([ann[a], ann[b]])
+                key = (cycle_idx, apex)
-                if apex in apex_positions:
+                apex_corners[key].extend([
+                    ann_positions[(cycle_idx, a)],
+                    ann_positions[(cycle_idx, b)],
+                ])
+                if key in apex_positions:
                     continue
-                angle = edge_midpoint_angle(i, n)
+                angle = edge_midpoint_angle(i, n, rotation)
-                if apex in node.up:
+                radius = 1.42 if apex in node.up else 0.58
-                    radius = 1.42
+                apex_positions[key] = (
-                else:
-                    radius = 0.58
-                apex_positions[apex] = (
                     dx + radius * math.cos(angle),
                     radius * math.sin(angle),
                 )
 
-    for apex, corners in apex_corners.items():
+    for key, corners in apex_corners.items():
+        _cycle_idx, apex = key
         if apex in node.bites and corners:
             cx = sum(p[0] for p in corners) / len(corners)
             cy = sum(p[1] for p in corners) / len(corners)
-            center_x = sum(offsets) / len(offsets) if offsets else 0.0
+            apex_positions[key] = (0.82 * cx, 0.82 * cy)
-            apex_positions[apex] = (
+
-                center_x + 0.82 * (cx - center_x),
+    return offsets, ann_positions, apex_positions, apex_corners
-                0.82 * cy,
+
-            )
+
-        pos = apex_positions[apex]
+def orientation(a, b, c) -> float:
+    return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])
+
+
+def segments_cross(a, b, c, d) -> bool:
+    eps = 1e-9
+    if max(a[0], b[0]) + eps < min(c[0], d[0]):
+        return False
+    if max(c[0], d[0]) + eps < min(a[0], b[0]):
+        return False
+    if max(a[1], b[1]) + eps < min(c[1], d[1]):
+        return False
+    if max(c[1], d[1]) + eps < min(a[1], b[1]):
+        return False
+    o1 = orientation(a, b, c)
+    o2 = orientation(a, b, d)
+    o3 = orientation(c, d, a)
+    o4 = orientation(c, d, b)
+    return (o1 * o2 < -eps) and (o3 * o4 < -eps)
+
+
+def polylines_cross(first, second) -> bool:
+    for i in range(len(first) - 1):
+        for j in range(len(second) - 1):
+            if segments_cross(first[i], first[i + 1], second[j], second[j + 1]):
+                return True
+    return False
+
+
+def quadratic_points(start, control, end, samples: int = 40):
+    points = []
+    for k in range(samples + 1):
+        t = k / samples
+        x = (1 - t) ** 2 * start[0] + 2 * (1 - t) * t * control[0] + t ** 2 * end[0]
+        y = (1 - t) ** 2 * start[1] + 2 * (1 - t) * t * control[1] + t ** 2 * end[1]
+        points.append((x, y))
+    return points
+
+
+def crossing_free_shared_arcs(
+    node: TreadNode,
+    apex_positions,
+    top_y: float,
+    bottom_y: float,
+):
+    shared = shared_up_apex_occurrences(node)
+    arc_specs = []
+    routed = []
+    pairs = []
+    for apex, occurrences in shared.items():
+        ordered = sorted(occurrences, key=lambda item: apex_positions[(item[0], apex)][0])
+        for left, right in zip(ordered, ordered[1:]):
+            pairs.append((apex, left, right))
+
+    pairs.sort(key=lambda item: (
+        abs(apex_positions[(item[1][0], item[0])][0] - apex_positions[(item[2][0], item[0])][0]),
+        min(item[1][0], item[2][0]),
+        item[0],
+    ))
+
+    for arc_idx, (apex, left, right) in enumerate(pairs):
+        start = apex_positions[(left[0], apex)]
+        end = apex_positions[(right[0], apex)]
+        if start[0] > end[0]:
+            start, end = end, start
+        for lane in range(64):
+            candidates = [
+                ((start[0] + end[0]) / 2, top_y + 0.42 + 0.18 * lane),
+                ((start[0] + end[0]) / 2, bottom_y - 0.42 - 0.18 * lane),
+            ]
+            for control in candidates:
+                points = quadratic_points(start, control, end)
+                if not any(polylines_cross(points, existing) for existing in routed):
+                    routed.append(points)
+                    arc_specs.append((start, control, end))
+                    break
+            else:
+                continue
+            break
+        else:
+            control = ((start[0] + end[0]) / 2, top_y + 0.42 + 0.18 * (64 + arc_idx))
+            routed.append(quadratic_points(start, control, end))
+            arc_specs.append((start, control, end))
+    return arc_specs
+
+
+def draw_tread_cycles(ax, node: TreadNode, connect_shared: bool):
+    rotations = compound_cycle_rotations(node) if connect_shared else [0.0] * len(node.annular_cycles)
+    display_order = compound_cycle_order(node) if connect_shared else None
+    offsets, ann_positions, apex_positions, apex_corners = cycle_layout(
+        node, rotations, display_order=display_order
+    )
+
+    for cycle_idx, order in enumerate(node.annular_cycles):
+        cyc_x = [ann_positions[(cycle_idx, v)][0] for v in order] + [
+            ann_positions[(cycle_idx, order[0])][0]
+        ]
+        cyc_y = [ann_positions[(cycle_idx, v)][1] for v in order] + [
+            ann_positions[(cycle_idx, order[0])][1]
+        ]
+        ax.plot(cyc_x, cyc_y, color="black", lw=1.3, zorder=2)
+
+    for key, corners in apex_corners.items():
+        pos = apex_positions[key]
         for corner in corners:
             ax.plot([pos[0], corner[0]], [pos[1], corner[1]], color="#9ca3af", lw=0.5)
 
-    for apex, pos in apex_positions.items():
+    if connect_shared and apex_positions:
+        all_positions = list(ann_positions.values()) + list(apex_positions.values())
+        top_y = max(p[1] for p in all_positions)
+        bottom_y = min(p[1] for p in all_positions)
+        for start, control, end in crossing_free_shared_arcs(
+            node, apex_positions, top_y, bottom_y
+        ):
+            path = MplPath([start, control, end], [MplPath.MOVETO, MplPath.CURVE3, MplPath.CURVE3])
+            ax.add_patch(
+                PathPatch(
+                    path,
+                    facecolor="none",
+                    edgecolor="#475569",
+                    lw=0.8,
+                    linestyle=(0, (1.2, 2.0)),
+                    zorder=1,
+                )
+            )
+
+    for (_cycle_idx, apex), pos in apex_positions.items():
         if apex in node.up:
             color, size, edgecolor = "#2563eb", 13, "none"
         elif apex in node.bites:
@@ -539,6 +665,34 @@ def draw_tread_model(ax, node: TreadNode):
             zorder=4,
         )
 
+    xs = [p[0] for p in list(ann_positions.values()) + list(apex_positions.values())]
+    ys = [p[1] for p in list(ann_positions.values()) + list(apex_positions.values())]
+    if connect_shared:
+        ys.append(max(ys) + 1.2)
+        ys.append(min(ys) - 1.2)
+    xpad = 1.1 if connect_shared else 1.7
+    ypad = 0.25 if connect_shared else 0.0
+    ax.set_xlim(min(xs, default=min(offsets, default=0.0)) - xpad, max(xs, default=max(offsets, default=0.0)) + xpad)
+    ax.set_ylim(min(ys, default=-1.65) - 0.25, max(ys, default=1.65) + ypad)
+    ax.set_aspect("equal")
+    ax.axis("off")
+
+
+def draw_tread_model(ax, node: TreadNode):
+    if len(node.annular_cycles) > 1:
+        draw_tread_cycles(ax, node, connect_shared=True)
+        singleton_down = set(node.down) - set(node.bites)
+        ax.set_title(
+            f"T{node.idx} d={node.depth}: {len(node.annular_cycles)} annular cycle(s)\n"
+            f"ann={len(node.annular)} up={len(node.up)} down={len(singleton_down)} "
+            f"bite={len(node.bites)}",
+            fontsize=6.4,
+            pad=1.5,
+        )
+        return
+
+    draw_tread_cycles(ax, node, connect_shared=False)
+
     singleton_down = set(node.down) - set(node.bites)
     ax.set_title(
         f"T{node.idx} d={node.depth}: {len(node.annular_cycles)} annular cycle(s)\n"
@@ -547,11 +701,6 @@ def draw_tread_model(ax, node: TreadNode):
         fontsize=6.4,
         pad=1.5,
     )
-    pad = 1.7
-    ax.set_xlim(min(offsets, default=0.0) - pad, max(offsets, default=0.0) + pad)
-    ax.set_ylim(-1.65, 1.65)
-    ax.set_aspect("equal")
-    ax.axis("off")
 
 
 def draw_medial_tire_grid(fig, outer_spec, nodes):
@@ -629,6 +778,8 @@ def draw_case(out_dir: Path, graph_idx: int, g: nx.Graph, source: int, augment:
         Line2D([0], [0], marker="o", color="w", label="inserted vertex",
                markerfacecolor="#fde68a", markeredgecolor="#7c3aed", markersize=8),
         Line2D([0], [0], color="black", lw=1.3, label="annular cycle A(T)"),
+        Line2D([0], [0], color="#475569", lw=0.8, linestyle=(0, (1.2, 2.0)),
+               label="shared up apex"),
         Line2D([0], [0], marker="o", color="w", label="up tooth",
                markerfacecolor="#2563eb", markersize=6),
         Line2D([0], [0], marker="o", color="w", label="down tooth",
@@ -636,7 +787,7 @@ def draw_case(out_dir: Path, graph_idx: int, g: nx.Graph, source: int, augment:
         Line2D([0], [0], marker="o", color="w", label="bite apex",
                markerfacecolor="#7f1d1d", markeredgecolor="black", markersize=6),
     ]
-    fig.legend(handles=legend, loc="lower center", ncol=5, fontsize=9)
+    fig.legend(handles=legend, loc="lower center", ncol=6, fontsize=9)
     fig.subplots_adjust(left=0.03, right=0.99, top=0.92, bottom=0.08, wspace=0.08, hspace=0.16)
 
     png = out_dir / f"random_c5_n30_medial_tire_decomposition_{graph_idx}.png"