ocr.py

#!/usr/bin/env python3
"""
Cardinal Chains OCR
Extracts puzzle structure AND solution paths from screenshot PNGs.

Image encoding:
  - Colored cell       = active grid cell; color = which chain it belongs to
  - White cell         = wall (not part of the puzzle)
  - X marker           = chain starting position
  - Two-dash symbol    = number 1 (pixel-art)
  - Digit 2/3/4...     = that number in pixel-art font
  - No border between same-chain adjacent cells that are CONSECUTIVE in the path
  - Border present between same-chain adjacent cells NOT consecutive in the path

Output JSON per level:
  {
    rows, cols,
    grid: [[{value, chain, start?} | null, ...], ...],
    chains: {
      "0": {start: [r,c], path: [[r,c], [r,c], ...]},
      ...
    },
    palette: [[R,G,B], ...]
  }
"""

import sys
import json
import numpy as np
from pathlib import Path
from PIL import Image
from scipy import ndimage


# ── Color helpers ──────────────────────────────────────────────────────────────

def color_distance(a, b):
    return float(np.sqrt(np.sum((np.array(a, float) - np.array(b, float)) ** 2)))

def dominant_bg_color(region_rgb):
    """Median non-dark, non-white color of a region."""
    pixels = region_rgb.reshape(-1, 3)
    valid = pixels[
        ~np.all(pixels < 80, axis=1) &
        ~np.all(pixels > 220, axis=1)
    ]
    if not len(valid):
        return None
    return tuple(int(v) for v in np.median(valid, axis=0))


# ── Step 1: Find the grid region ───────────────────────────────────────────────

def find_nonwhite_bounds(img):
    mask = ~np.all(img > 240, axis=2)
    rows = np.where(mask.any(axis=1))[0]
    cols = np.where(mask.any(axis=0))[0]
    if not len(rows) or not len(cols):
        raise ValueError("No non-white pixels found")
    return int(rows[0]), int(rows[-1]), int(cols[0]), int(cols[-1])

def find_border_width(dark_mask):
    """Outer border thickness on each side (rows/cols that are ≥60% dark).
    Tolerates up to 5px of non-dark gap before the actual border line."""
    def bw(arr):
        n = len(arr)
        # Skip up to 5 non-dark pixels (e.g. 1px rendering gap before grid border)
        start = 0
        while start < min(5, n) and arr[start] < 0.60:
            start += 1
        count = 0
        for i in range(start, n):
            if arr[i] >= 0.60:
                count += 1
            else:
                break
        return start + count if count > 0 else 0
    t = bw(dark_mask.mean(axis=1))
    b = bw(dark_mask.mean(axis=1)[::-1])
    l = bw(dark_mask.mean(axis=0))
    r = bw(dark_mask.mean(axis=0)[::-1])
    return t, b, l, r


# ── Step 2: Separator detection ────────────────────────────────────────────────

def find_color_change_seps(scan, min_dark_run=3, color_thr=90, large_run=40,
                           check_edges=True):
    """
    Find separator positions in a 1-D color scan (N×3).

    A dark run qualifies when the cells on each side have clearly different
    visual categories (different chain colors, or active vs wall).

    color_thr=90: high enough to ignore same-chain rendering variation (~50-80)
    while catching real chain-to-chain transitions (always ≥100 apart).

    For runs ≥ large_run px (merged wall-border): both edges are emitted.
    """
    N = len(scan)
    dark = np.all(scan < 80, axis=1)
    runs = []
    in_run = False
    for i, d in enumerate(dark):
        if d and not in_run:
            run_start = i; in_run = True
        elif not d and in_run:
            runs.append((run_start, i)); in_run = False
    if in_run:
        runs.append((run_start, N))

    def sample_side(pixels):
        """
        Returns (color_or_None, is_wall).
        is_wall=True if all non-dark pixels are clearly white (> 200 all channels).
        color_or_None is the median colored background, or None if indeterminate.
        Skip the first 3px (far-end border shadow) and the last 2px (near-end
        transition pixels that are often anti-aliased shadow of the dark run itself).
        """
        if len(pixels) > 5:
            p = pixels[3:-2]
        elif len(pixels) > 3:
            p = pixels[3:]
        else:
            p = pixels
        if len(p) == 0:
            return None, False
        non_dark = p[~np.all(p < 80, axis=1)]
        if len(non_dark) == 0:
            return None, False
        # Is it a wall? All non-dark pixels are white-ish (> 200 on all channels)
        if np.all(non_dark > 200):
            return None, True
        # Colored background: exclude near-white pixels too
        colored = non_dark[~np.all(non_dark > 200, axis=1)]
        if len(colored) == 0:
            return None, False
        return tuple(np.median(colored, axis=0)), False

    positions = []
    for rs, re in runs:
        run_len = re - rs
        if run_len < min_dark_run:
            continue

        at_edge = check_edges and ((rs <= 10) or (re >= N - 10))

        if at_edge:
            is_sep = True
        else:
            left_win  = scan[max(0, rs - 25): rs]
            right_win = scan[re: min(N, re + 25)]
            lc, l_wall = sample_side(left_win)
            rc, r_wall = sample_side(right_win[::-1])  # read from run outward

            if l_wall and r_wall:
                # Both walls: true sep lines are thin (≤ 8px); symbol dashes are
                # much longer (≥ 12px), so gate on run length.
                is_sep = (run_len <= 8)
            elif l_wall or r_wall:
                # Exactly one side is a wall cell → definite separator.
                is_sep = True
            elif lc is None or rc is None:
                # Cannot determine color (symbol interior) → not a separator
                is_sep = False
            else:
                is_sep = color_distance(lc, rc) >= color_thr

        if not is_sep:
            continue

        if run_len >= large_run:
            positions.append(rs + 2)
            positions.append(re - 2)
        else:
            positions.append((rs + re) // 2)

    return positions


def detect_all_seps(region, scan_step=5, border=8):
    """
    Densely scan the inner region (outer borders already stripped) in both
    axes; vote on separator positions.  Working on the inner image prevents
    anti-aliased outer-border pixels from contaminating sample_side().
    Returns (row_sep_positions, col_sep_positions) in *region* coordinates.
    """
    H, W = region.shape[:2]
    dark_mask = np.all(region < 80, axis=2)
    bt, bb, bl, br = find_border_width(dark_mask)
    inner = region[bt: H - bb if bb else H, bl: W - br if br else W]
    iH, iW = inner.shape[:2]

    row_votes = {}   # position (region-coord) → vote count
    col_votes = {}   # position (region-coord) → vote count

    for c in range(1, iW - 1, scan_step):
        for p in find_color_change_seps(inner[:, c], check_edges=False):
            rp = p + bt   # convert inner → region coord
            row_votes[rp] = row_votes.get(rp, 0) + 1

    for r in range(1, iH - 1, scan_step):
        for p in find_color_change_seps(inner[r, :], check_edges=False):
            cp = p + bl
            col_votes[cp] = col_votes.get(cp, 0) + 1

    def cluster_votes(votes, tol=8, min_votes=2):
        if not votes:
            return []
        positions = sorted(votes)
        groups = [[positions[0]]]
        for p in positions[1:]:
            if p - groups[-1][-1] <= tol:
                groups[-1].append(p)
            else:
                groups.append([p])
        result = []
        for g in groups:
            total = sum(votes[p] for p in g)
            if total >= min_votes:
                result.append(int(np.median(g)))
        return sorted(result)

    return cluster_votes(row_votes), cluster_votes(col_votes)


def infer_complete_seps(detected, extent, min_cell=30, cell_size=None):
    """
    Given a list of detected separator positions (in [0, extent]) and the
    grid extent, fill in missing separators by finding the dominant cell
    spacing and repeating it.

    Always includes 0 and extent as outer boundaries.
    If cell_size is given, it is used directly instead of being inferred.
    """
    if not detected:
        return [0, extent]

    # Ensure outer boundaries are present
    detected = sorted(set([0] + list(detected) + [extent]))

    if cell_size is None:
        # Find spacings between consecutive separators
        spacings = [detected[i + 1] - detected[i] for i in range(len(detected) - 1)]
        spacings = [s for s in spacings if s > min_cell]
        if not spacings:
            return detected

        # The most common spacing (mode) is the cell size
        from collections import Counter
        rounded = [round(s / 5) * 5 for s in spacings]  # round to 5px
        cell_size = Counter(rounded).most_common(1)[0][0]

    if cell_size == 0:
        return detected

    # Build a complete set of separators starting from the first boundary
    seps = [detected[0]]
    pos = detected[0] + cell_size
    while pos < extent - cell_size * 0.3:
        # Snap to nearest detected separator if within tolerance
        near = [d for d in detected if abs(d - pos) < cell_size * 0.25]
        if near:
            pos = min(near, key=lambda d: abs(d - pos))
        seps.append(int(round(pos)))
        pos += cell_size
    seps.append(detected[-1])

    # Deduplicate and sort
    return sorted(set(seps))


def dark_profile_sep_positions(profile, extent, min_frac=0.15):
    """
    Find separator positions from a dark-fraction profile using a natural-break
    threshold.  Returns midpoints of runs above the threshold, in profile coords.
    """
    vals = sorted(v for v in profile if v >= min_frac)
    if not vals:
        return []
    if len(vals) < 2:
        return [int(np.argmax(profile))]
    diffs = [vals[i + 1] - vals[i] for i in range(len(vals) - 1)]
    max_gap = max(diffs)
    if max_gap < 0.05:
        threshold = max(vals) * 0.7
    else:
        gap_idx = max(range(len(diffs)), key=lambda i: diffs[i])
        threshold = (vals[gap_idx] + vals[gap_idx + 1]) / 2
    seps = []
    i = 0
    while i < extent:
        if profile[i] >= threshold:
            j = i
            while j < extent and profile[j] >= threshold:
                j += 1
            seps.append((i + j - 1) // 2)
            i = j
        else:
            i += 1
    return seps


# ── Step 3: Fall-back blob detection for same-color grids ─────────────────────

def find_symbol_centers(region, bt, bb, bl, br):
    """
    Find dark-blob centers in the inner region (excluding outer borders).
    Returns list of (cy, cx) in inner-region coordinates.
    """
    H, W = region.shape[:2]
    inner = region[bt: H - bb if bb else H, bl: W - br if br else W]
    iH, iW = inner.shape[:2]

    dark_mask = np.all(inner < 80, axis=2)
    labeled, n = ndimage.label(dark_mask)

    blobs = []
    for i in range(1, n + 1):
        ys, xs = np.where(labeled == i)
        size = len(ys)
        if size < 8:
            continue
        bh = int(ys.max() - ys.min() + 1)
        bw = int(xs.max() - xs.min() + 1)
        blobs.append({
            'cy': float(np.mean(ys)),
            'cx': float(np.mean(xs)),
            'size': size,
            'h': bh,
            'w': bw,
        })

    if not blobs:
        return []

    # Filter out very large structural blobs (borders, wall outlines)
    max_size = max(b['size'] for b in blobs)
    threshold = min(max_size * 0.25, 1500)
    symbol_blobs = [b for b in blobs if b['size'] <= threshold] or blobs

    # Filter out separator-line artifacts: very wide, very short blobs (aspect ratio > 8:1)
    symbol_blobs = [b for b in symbol_blobs if b['w'] <= 8 * b['h']] or symbol_blobs

    # Return one center per blob; infer_grid_from_blobs will cluster these
    return [(b['cy'], b['cx']) for b in symbol_blobs]


def natural_break_threshold(vals, fallback_frac=0.4):
    """
    Find the threshold that separates "within-symbol" distances from
    "between-cell" distances using the largest gap in the sorted diff list.
    """
    if len(vals) < 2:
        return fallback_frac * (max(vals) - min(vals)) if vals else 10
    s = sorted(vals)
    diffs = [s[i + 1] - s[i] for i in range(len(s) - 1)]
    if not diffs:
        return 10
    sdiffs = sorted(diffs)
    gaps = [sdiffs[i + 1] - sdiffs[i] for i in range(len(sdiffs) - 1)]
    if not gaps or max(gaps) < 3:
        return float(np.mean(diffs)) * 0.7
    split_idx = int(np.argmax(gaps))
    return (sdiffs[split_idx] + sdiffs[split_idx + 1]) / 2


def cluster_1d_natural(vals):
    """Cluster 1-D values using natural-break threshold."""
    if not vals:
        return []
    tol = natural_break_threshold(vals)
    tol = max(tol, 5)
    vals = sorted(vals)
    groups = [[vals[0]]]
    for v in vals[1:]:
        if v - groups[-1][-1] <= tol:
            groups[-1].append(v)
        else:
            groups.append([v])
    return [float(np.mean(g)) for g in groups]


def _grid_from_spacing(centers_1d, extent, cell_size):
    """
    Given blob positions along one axis, a cell_size, and the axis extent,
    return separator positions by detecting the grid phase and snapping to a
    regular lattice.
    """
    vals = np.array(centers_1d, float)
    n_cells = max(1, round(extent / cell_size))

    # Phase detection via circular mean of (val mod cell_size)
    angles = (vals % cell_size) / cell_size * 2 * np.pi
    sin_mean = float(np.mean(np.sin(angles)))
    cos_mean = float(np.mean(np.cos(angles)))
    offset = float(np.arctan2(sin_mean, cos_mean) / (2 * np.pi) * cell_size) % cell_size

    # Normalise offset into [0, cell_size): first centre is inside the image.
    offset = offset % cell_size

    # n_cells cell centres → n_cells−1 interior separators + 0 and extent
    centres = [offset + k * cell_size for k in range(n_cells)]
    seps = [0]
    for i in range(len(centres) - 1):
        mid = int(round((centres[i] + centres[i + 1]) / 2))
        if 0 < mid < extent:
            seps.append(mid)
    seps.append(extent)
    return sorted(set(seps))


def infer_grid_from_blobs(centers, iH, iW):
    """
    Fit a regular grid to symbol centers using row spacing as the cell-size
    anchor, then derive columns via circular-phase detection.
    Returns (row_seps, col_seps) in inner coordinates, or None.
    """
    if not centers:
        return None

    # Row clustering via natural break (works well — within-row y-spread is
    # much smaller than between-row gap)
    row_c = cluster_1d_natural([c[0] for c in centers])
    n_rows = len(row_c)
    if n_rows == 0:
        return None

    if n_rows > 1:
        row_spacings = [row_c[i + 1] - row_c[i] for i in range(n_rows - 1)]
        cell_size = float(np.median(row_spacings))
    else:
        cell_size = float(iH)

    # Row separators: midpoints between successive row centres + outer bounds
    row_seps = [0]
    for i in range(n_rows - 1):
        row_seps.append(int(round((row_c[i] + row_c[i + 1]) / 2)))
    row_seps.append(iH)

    # Column positions: use phase detection with estimated cell_size
    col_seps = _grid_from_spacing([c[1] for c in centers], iW, cell_size)

    return sorted(set(row_seps)), sorted(set(col_seps))


# ── Step 4: Build cell grid ────────────────────────────────────────────────────

def seps_to_cells(row_seps, col_seps):
    cells = []
    for i in range(len(row_seps) - 1):
        row = []
        for j in range(len(col_seps) - 1):
            row.append((row_seps[i], row_seps[i + 1],
                        col_seps[j], col_seps[j + 1]))
        cells.append(row)
    return cells


# ── Step 5: Cell classification ────────────────────────────────────────────────

def classify_cell(inner, t, b, l, r):
    mh = max(1, (b - t) // 6)
    mw = max(1, (r - l) // 6)
    cell = inner[t + mh: b - mh, l + mw: r - mw]
    if cell.size == 0:
        return {'type': 'wall', 'color': None, 'has_x': False, 'number': None}

    color = dominant_bg_color(cell)
    if color is None:
        return {'type': 'wall', 'color': None, 'has_x': False, 'number': None}

    dark_frac = float(np.mean(np.all(cell < 80, axis=2)))
    # X marks fill ~16-30% of the cell; digit symbols fill <14%
    has_x = dark_frac > 0.15

    return {'type': 'active', 'color': color, 'has_x': has_x, 'number': None}


# ── Step 6: Number reading ─────────────────────────────────────────────────────

def read_symbol_fingerprint(inner, t, b, l, r):
    mh = max(1, (b - t) // 8)
    mw = max(1, (r - l) // 8)
    cell = inner[t + mh: b - mh, l + mw: r - mw]
    if cell.size == 0:
        return None
    dark = (np.all(cell < 80, axis=2).astype(np.uint8)) * 255
    thumb = np.array(Image.fromarray(dark).resize((8, 8), Image.LANCZOS)) > 128
    return thumb


def assign_numbers(inner, cells_bounds, classifications):
    fingerprints = {}
    dark_fracs = {}

    for ri, row in enumerate(cells_bounds):
        for ci, (t, b, l, r) in enumerate(row):
            info = classifications[ri][ci]
            if info['type'] != 'active':
                continue
            if info['has_x']:
                info['number'] = 'x'
                continue
            fp = read_symbol_fingerprint(inner, t, b, l, r)
            if fp is None:
                info['number'] = 1
                continue
            mh = max(1, (b - t) // 6)
            mw = max(1, (r - l) // 6)
            cell = inner[t + mh: b - mh, l + mw: r - mw]
            dark_fracs[(ri, ci)] = float(np.mean(np.all(cell < 80, axis=2)))
            fingerprints[(ri, ci)] = fp

    # Cluster fingerprints
    templates = []
    pos_to_tmpl = {}
    for pos, fp in fingerprints.items():
        best_idx, best_dist = None, 64
        for idx, (tmpl, _) in enumerate(templates):
            d = int(np.sum(fp != tmpl))
            if d < best_dist:
                best_dist = d; best_idx = idx
        if best_dist <= 18:
            templates[best_idx][1].append(pos)
            pos_to_tmpl[pos] = best_idx
        else:
            pos_to_tmpl[pos] = len(templates)
            templates.append((fp, [pos]))

    # Sort templates by avg dark fraction → digit 1, 2, 3, …
    def avg_frac(t):
        return np.mean([dark_fracs.get(p, 0) for p in t[1]]) if t[1] else 0

    templates_sorted = sorted(range(len(templates)), key=lambda i: avg_frac(templates[i]))
    digit_map = {tidx: digit for digit, tidx in enumerate(templates_sorted, 1)}

    for pos, tidx in pos_to_tmpl.items():
        ri, ci = pos
        classifications[ri][ci]['number'] = digit_map.get(tidx, 1)

    return classifications


# ── Step 7: Color → chain ID ───────────────────────────────────────────────────

def assign_chain_ids(classifications):
    palette = []
    for row in classifications:
        for info in row:
            if info['type'] != 'active' or info['color'] is None:
                continue
            c = info['color']
            match = next((i for i, rep in enumerate(palette) if color_distance(c, rep) < 30), None)
            if match is None:
                match = len(palette)
                palette.append(c)
            info['chain_id'] = match
    return classifications, palette


# ── Step 8: Path extraction ────────────────────────────────────────────────────

def check_border_between(inner, r1, c1, r2, c2, cells_bounds, dark_threshold=0.05):
    """
    Return True if there is a visible dark separator between adjacent cells
    (r1,c1) and (r2,c2).  False means they are connected (consecutive in path).

    Adjacent cells share a boundary coordinate, so we sample a band of ±half
    pixels around that boundary.  For a 2-4 px separator line this gives a
    reliable dark-fraction signal without cutting into cell content.
    """
    t1, b1, l1, r_1 = cells_bounds[r1][c1]
    t2, b2, l2, r_2 = cells_bounds[r2][c2]
    iH, iW = inner.shape[:2]

    if r1 == r2:
        # Horizontally adjacent: look for a vertical dark line at x = r_1 ≈ l2
        sep_x = (r_1 + l2) // 2
        half = max(3, min(r_1 - l1, r_2 - l2) // 8)
        sep_l = max(0, sep_x - half)
        sep_r = min(iW, sep_x + half)
        if sep_r <= sep_l:
            return False
        cell_h = min(b1, b2) - max(t1, t2)
        q = max(1, cell_h // 5)
        top = max(t1, t2) + q
        bot = min(b1, b2) - q
        strip = inner[top:bot, sep_l:sep_r]
    else:
        # Vertically adjacent: look for a horizontal dark line at y = b1 ≈ t2
        sep_y = (b1 + t2) // 2
        half = max(3, min(b1 - t1, b2 - t2) // 8)
        sep_t = max(0, sep_y - half)
        sep_b = min(iH, sep_y + half)
        if sep_b <= sep_t:
            return False
        cell_w = min(r_1, r_2) - max(l1, l2)
        q = max(1, cell_w // 5)
        left = max(l1, l2) + q
        right = min(r_1, r_2) - q
        strip = inner[sep_t:sep_b, left:right]

    if strip.size == 0:
        return False

    dark_frac = float(np.mean(np.all(strip < 80, axis=2)))
    return dark_frac >= dark_threshold


def extract_chain_paths(inner, cells_bounds, classifications):
    """
    For each chain, determine the ordered path from X to end by finding
    which adjacent same-chain cell pairs have NO border between them
    (those are consecutive in the path).

    Returns dict: chain_id → {'start': [r,c], 'path': [[r,c], ...]}
    """
    nrows = len(cells_bounds)
    ncols = max(len(row) for row in cells_bounds)

    # Build adjacency: for each cell, list all same-chain neighbors + border status
    adj = {}  # (r,c) → list of (r2,c2, has_border)

    for ri in range(nrows):
        for ci in range(len(cells_bounds[ri])):
            info = classifications[ri][ci]
            if info['type'] != 'active':
                continue
            chain_id = info.get('chain_id', -1)
            neighbors = []
            for dr, dc in [(-1, 0), (1, 0), (0, -1), (0, 1)]:
                nr, nc = ri + dr, ci + dc
                if nr < 0 or nr >= nrows or nc < 0 or nc >= len(cells_bounds[nr]):
                    continue
                ninfo = classifications[nr][nc]
                if ninfo['type'] != 'active' or ninfo.get('chain_id') != chain_id:
                    continue
                has_border = check_border_between(inner, ri, ci, nr, nc, cells_bounds)
                neighbors.append((nr, nc, has_border))
            adj[(ri, ci)] = neighbors

    # For each chain, find the start cell and traverse the no-border path
    chains_info = {}
    start_cells = {}  # chain_id → (r, c)

    for ri, row in enumerate(cells_bounds):
        for ci in range(len(row)):
            info = classifications[ri][ci]
            if info['type'] == 'active' and info.get('has_x'):
                cid = info.get('chain_id', -1)
                start_cells[cid] = (ri, ci)

    for cid, start in start_cells.items():
        # BFS/DFS along no-border edges from start
        path = []
        visited = set()
        stack = [start]
        # The path graph should be a simple path (no branching)
        # Use iterative DFS following only no-border edges
        current = start
        visited.add(current)
        path.append(list(current))

        while True:
            r, c = current
            # Find unvisited no-border neighbors
            no_border_neighbors = [
                (nr, nc) for nr, nc, hb in adj.get((r, c), [])
                if not hb and (nr, nc) not in visited
            ]
            if not no_border_neighbors:
                break
            # Prefer neighbor with lower value first (non-decreasing constraint)
            def cell_val(pos):
                v = classifications[pos[0]][pos[1]].get('number', 1)
                if v == 'x':
                    return 0
                return v if isinstance(v, int) else 0

            next_cell = min(no_border_neighbors, key=cell_val)
            visited.add(next_cell)
            path.append(list(next_cell))
            current = next_cell

        chains_info[cid] = {'start': list(start), 'path': path}

    return chains_info


# ── Main pipeline ──────────────────────────────────────────────────────────────

def parse_image(path):
    img = np.array(Image.open(path).convert('RGB'))

    # 1. Find grid region
    top, bottom, left, right = find_nonwhite_bounds(img)
    region = img[top: bottom + 1, left: right + 1]
    H, W = region.shape[:2]

    # 2. Outer border
    dark_mask = np.all(region < 80, axis=2)
    bt, bb, bl, br = find_border_width(dark_mask)
    inner = region[bt: H - bb if bb else H, bl: W - br if br else W]
    iH, iW = inner.shape[:2]

    # 3. Detect separators (densely, with voting)
    raw_row_seps, raw_col_seps = detect_all_seps(region, scan_step=5, border=bt + 1)

    def to_inner(seps, offset, extent):
        return sorted(set(
            [0]
            + [s - offset for s in seps if offset < s < offset + extent]
            + [extent]
        ))

    row_seps = infer_complete_seps(to_inner(raw_row_seps, bt, iH), iH)
    col_seps = infer_complete_seps(to_inner(raw_col_seps, bl, iW), iW)

    def spacing_from_seps(seps):
        gaps = [seps[i+1] - seps[i] for i in range(len(seps)-1)]
        return float(np.median(gaps)) if gaps else None

    # ── Fallback when scan gives implausibly small cell spacing ────────────────
    # Happens when "1"-symbol dash edges flood col/row voting with false positives.
    dark_inner_mask = np.all(inner < 80, axis=2)

    col_spacing = spacing_from_seps(col_seps) or iW
    if col_spacing < 30:
        col_dp = dark_profile_sep_positions(dark_inner_mask.mean(axis=0), iW)
        if col_dp:
            col_seps = infer_complete_seps(col_dp, iW)
            col_spacing = spacing_from_seps(col_seps) or iW

    row_spacing = spacing_from_seps(row_seps) or iH
    if col_spacing > row_spacing * 1.3 and col_spacing >= 30:
        # Col detection is much more reliable; rebuild rows using col cell size.
        row_raw_inner = [s - bt for s in raw_row_seps if 0 < s - bt < iH]
        row_seps = infer_complete_seps(row_raw_inner, iH,
                                       cell_size=int(round(col_spacing)))
        row_spacing = spacing_from_seps(row_seps) or iH
    # ── End fallback ───────────────────────────────────────────────────────────

    n_row_int = sum(1 for s in row_seps if 0 < s < iH)
    n_col_int = sum(1 for s in col_seps if 0 < s < iW)

    def build_seps_from_spacing(spacing, extent):
        seps = [0]
        pos = spacing
        while pos < extent - spacing * 0.3:
            seps.append(int(round(pos)))
            pos += spacing
        seps.append(extent)
        return sorted(set(seps))

    if n_row_int >= 1 and n_col_int == 0:
        # Use row spacing to infer column grid
        col_seps = build_seps_from_spacing(spacing_from_seps(row_seps), iW)
        n_col_int = len(col_seps) - 2
    elif n_col_int >= 1 and n_row_int == 0:
        # Use col spacing to infer row grid
        row_seps = build_seps_from_spacing(spacing_from_seps(col_seps), iH)
        n_row_int = len(row_seps) - 2

    use_seps = n_row_int >= 1 and n_col_int >= 1

    # Validate: if cell size is implausibly large (>120px), the separator
    # detection caught only a chain boundary, not the full grid — fall back.
    if use_seps:
        row_spacing = spacing_from_seps(row_seps) or iH
        col_spacing = spacing_from_seps(col_seps) or iW
        if row_spacing > 120 or col_spacing > 120:
            use_seps = False

    if not use_seps:
        # Fall back to blob detection
        centers = find_symbol_centers(region, bt, bb, bl, br)
        inferred = infer_grid_from_blobs(centers, iH, iW)
        if inferred is None:
            raise ValueError(f"Cannot infer grid layout from {path}")
        row_seps, col_seps = inferred

    # 4. Build cells
    cells_bounds = seps_to_cells(row_seps, col_seps)

    # 5. Classify cells
    classifications = [
        [classify_cell(inner, t, b, l, r) for (t, b, l, r) in row]
        for row in cells_bounds
    ]

    # 6. Read numbers
    classifications = assign_numbers(inner, cells_bounds, classifications)

    # 7. Assign chain IDs
    classifications, palette = assign_chain_ids(classifications)

    # 8. Extract paths
    chains_info = extract_chain_paths(inner, cells_bounds, classifications)

    # 9. Build output grid
    grid = []
    for ri, row in enumerate(cells_bounds):
        grow = []
        for ci in range(len(row)):
            info = classifications[ri][ci]
            if info['type'] == 'wall':
                grow.append(None)
            else:
                entry = {
                    'value': info['number'],
                    'chain': info.get('chain_id', -1),
                }
                if info.get('has_x'):
                    entry['start'] = True
                grow.append(entry)
        grid.append(grow)

    return {
        'rows': len(cells_bounds),
        'cols': max(len(r) for r in cells_bounds) if cells_bounds else 0,
        'grid': grid,
        'chains': {str(k): v for k, v in chains_info.items()},
        'palette': [list(c) for c in palette],
    }


class NumpyEncoder(json.JSONEncoder):
    def default(self, obj):
        if isinstance(obj, np.integer):
            return int(obj)
        if isinstance(obj, np.floating):
            return float(obj)
        if isinstance(obj, np.ndarray):
            return obj.tolist()
        return super().default(obj)


if __name__ == '__main__':
    path = sys.argv[1] if len(sys.argv) > 1 else '.-001.png'
    result = parse_image(path)
    print(json.dumps(result, indent=2, cls=NumpyEncoder))