On lines

import itertools

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, Patch
from matplotlib.colors import LinearSegmentedColormap, ListedColormap

from skimage.draw import line, line_nd

# Comparators.
import cv2
from PIL import Image, ImageDraw

# Slots 1 and 2 of the reference categorical palette, validated as a pair:
# CVD dE 24.7, normal-vision dE 33.6, both well clear of the floors.
C_LINE = "#2a78d6"   # skimage.draw.line
C_ND = "#eb6834"     # skimage.draw.line_nd
C_BOTH = "#c9c8c1"   # pixels that both algorithms choose
C_OFF = "#f2f1ec"    # pixels neither chooses
INK = "#0b0b0b"
MUTED = "#52514e"
GRID = "#dedcd5"

plt.rcParams.update(
    {
        "figure.dpi": 110,
        "font.size": 9,
        "axes.titlesize": 9,
        "axes.titlecolor": MUTED,
        "figure.facecolor": "white",
    }
)

def pixel_axes(ax, shape, title=None):
    """Draw an empty pixel grid, axis 0 downwards and axis 1 rightwards."""
    n_i, n_j = shape
    ax.set_xlim(-0.5, n_j - 0.5)
    ax.set_ylim(n_i - 0.5, -0.5)
    ax.set_xticks(range(n_j))
    ax.set_yticks(range(n_i))
    ax.set_xticks(np.arange(n_j + 1) - 0.5, minor=True)
    ax.set_yticks(np.arange(n_i + 1) - 0.5, minor=True)
    ax.grid(which="minor", color=GRID, linewidth=0.8)
    ax.tick_params(which="both", length=0, labelsize=7, colors=MUTED)
    ax.set_aspect("equal")
    for spine in ax.spines.values():
        spine.set_visible(False)
    if title is not None:
        ax.set_title(title)
    return ax


def fill(ax, pixels, color, alpha=1.0):
    """Fill each ``(axis 0, axis 1)`` pixel in `pixels`."""
    for i, j in sorted(pixels):
        ax.add_patch(
            Rectangle(
                (j - 0.5, i - 0.5),
                1,
                1,
                facecolor=color,
                edgecolor="white",
                linewidth=1.0,
                alpha=alpha,
                zorder=1,
            )
        )
    return ax


def exact(ax, p, q, color=INK):
    """Overlay the ideal segment from `p` to `q`, with its endpoints."""
    ax.plot([p[1], q[1]], [p[0], q[0]], color=color, linewidth=1.4, zorder=3)
    ax.plot([p[1], q[1]], [p[0], q[0]], "o", color=color, markersize=4, zorder=4)
    return ax


def as_set(coords):
    """A tuple of index arrays, as a set of ``(axis 0, axis 1)`` pairs."""
    return set(zip(*(np.asarray(a).tolist() for a in coords)))


def sk_line(p, q):
    return as_set(line(p[0], p[1], q[0], q[1]))


def sk_nd(p, q):
    return as_set(line_nd(p, q, endpoint=True))

p, q = (0, 0), (1, 4)
shape = (3, 6)
# Pixel squares centred on integer indices: the segment meets these six.
touches = {(0, 0), (0, 1), (0, 2), (1, 2), (1, 3), (1, 4)}

fig, ax = plt.subplots(figsize=(3.6, 2.0))
pixel_axes(ax, shape)
fill(ax, touches, C_OFF)
exact(ax, p, q)
ax.set_title("the segment (0, 0) to (1, 4) and the pixels it touches")
fig.tight_layout()

def bresenham(start, stop):
    """Bresenham's line, in the same steps as `skimage.draw.line`."""
    start, stop = np.array(start), np.array(stop)
    delta = np.abs(stop - start)
    step = np.sign(stop - start)

    major = int(np.argmax(delta))    # the axis with further to travel
    minor = 1 - major

    # `error` is the midpoint test of the previous section, cleared of
    # fractions by the factor 2 * delta[major]:
    #
    #     error(k) = 2 * delta[major] * [ y(k + 1) - (m(k) + 1/2) ]
    #
    # The 2 clears the half, delta[major] clears the division, and both are
    # positive so the sign is untouched. At k = 0 the minor axis has not
    # moved, which leaves the expression below.
    error = 2 * delta[minor] - delta[major]

    at = start.copy()
    pixels = []
    for _ in range(delta[major]):
        pixels.append(tuple(at))

        # error(k + 1) - error(k) = 2 * delta[minor] - 2 * delta[major] * s,
        # where s is 1 when the minor axis steps and 0 when it does not. So
        # the subtraction is conditional and the addition is not.
        if error >= 0:
            at[minor] += step[minor]
            error -= 2 * delta[major]
        at[major] += step[major]
        error += 2 * delta[minor]

    pixels.append(tuple(stop))       # the endpoint is written, never computed
    return pixels

R = range(-7, 8)
all_pairs = [((a, b), (c, d)) for a, b, c, d in itertools.product(R, repeat=4)]


def sk_sequence(p, q):
    ii, jj = line(p[0], p[1], q[0], q[1])
    return list(zip(ii.tolist(), jj.tolist()))


matches = sum(bresenham(p, q) == sk_sequence(p, q) for p, q in all_pairs)
print(f"identical to skimage.draw.line on {matches}/{len(all_pairs)} pairs"
      f"  ({matches / len(all_pairs):.1%})")

def discrepancies(start, stop):
    """Per step: k, the exact position, the drawn integer, the gap, the error."""
    start, stop = np.array(start), np.array(stop)
    delta = np.abs(stop - start)
    major = int(np.argmax(delta))
    minor = 1 - major
    d_major, d_minor = int(delta[major]), int(delta[minor])
    if d_major == 0:
        return []

    error, drawn, out = 2 * d_minor - d_major, 0, []
    for k in range(d_major + 1):
        y_exact = k * d_minor / d_major
        out.append((k, y_exact, drawn, y_exact - drawn, error))
        if k == d_major:
            break
        if error >= 0:
            drawn += 1
            error -= 2 * d_major
        error += 2 * d_minor
    return out


gaps = [row[3] for p, q in all_pairs for row in discrepancies(p, q)]
print(f"discrepancy over {len(all_pairs)} lines: [{min(gaps)}, {max(gaps)})")
print(f"   never below -1/2 : {min(gaps) >= -0.5}")
print(f"   always below 1/2 : {max(gaps) < 0.5}")

print(f"\nthe tie case, (0, 0) to (1, 4):")
print(f"{'k':>4}{'y(k) exact':>13}{'m(k) drawn':>13}{'gap':>8}{'error':>8}")
for k, y_exact, drawn, gap, err in discrepancies((0, 0), (1, 4)):
    print(f"{k:>4}{y_exact:>13.2f}{drawn:>13}{gap:>8.2f}{err:>8}")

def error_matches_closed_form(start, stop):
    """Is `error` equal to 2 * D * [y(k + 1) - (m(k) + 1/2)] at every step?"""
    start, stop = np.array(start), np.array(stop)
    delta = np.abs(stop - start)
    major = int(np.argmax(delta))
    minor = 1 - major
    d_major, d_minor = int(delta[major]), int(delta[minor])

    error, drawn = 2 * d_minor - d_major, 0
    for k in range(d_major):
        closed_form = 2 * d_major * ((k + 1) * d_minor / d_major - drawn - 0.5)
        if abs(error - closed_form) > 1e-9:
            return False
        if error >= 0:
            error -= 2 * d_major
            drawn += 1
        error += 2 * d_minor
    return True


agree = sum(error_matches_closed_form(p, q) for p, q in all_pairs)
print(f"error matches the closed form on {agree}/{len(all_pairs)} pairs")

def decisions(start, stop):
    """The error sequence alone, with no positions."""
    start, stop = np.array(start), np.array(stop)
    delta = np.abs(stop - start)
    d_major, d_minor = int(delta.max()), int(delta.min())
    error, out = 2 * d_minor - d_major, []
    for _ in range(d_major):
        out.append(error)
        if error >= 0:
            error -= 2 * d_major
        error += 2 * d_minor
    return out


base = decisions((0, 0), (2, 7))
elsewhere = {
    "up and left, (-2, -7)": ((0, 0), (-2, -7)),
    "down and left, (2, -7)": ((0, 0), (2, -7)),
    "transposed, (7, 2)": ((0, 0), (7, 2)),
    "translated by (5, 5)": ((5, 5), (7, 12)),
}
for name, (start, stop) in elsewhere.items():
    print(f"{name:<24} same error sequence as (2, 7): "
          f"{decisions(start, stop) == base}")

def never_steps_twice(start, stop):
    """Would a second pass of the `while` body ever be taken?"""
    start, stop = np.array(start), np.array(stop)
    delta = np.abs(stop - start)
    d_major, d_minor = int(delta.max()), int(delta.min())
    error = 2 * d_minor - d_major
    for _ in range(d_major):
        if error >= 0:
            error -= 2 * d_major
            if error >= 0:
                return False
        error += 2 * d_minor
    return True


single = sum(never_steps_twice(p, q) for p, q in all_pairs)
print(f"one minor step per major step suffices on {single}/{len(all_pairs)} pairs")

fig, ax = plt.subplots(figsize=(3.6, 2.0))
pixel_axes(ax, shape)
fill(ax, sk_line(p, q), C_LINE)
exact(ax, p, q, color="white")
ax.set_title("line(0, 0, 1, 4)")
fig.tight_layout()

print(f"{'sample':>8}{'rounds to':>12}{'direction':>12}{'error':>9}"
      f"{'because':>28}")
for value in (0.5, 1.5, 2.5, 3.5):
    got = np.round(value)
    below, above = int(np.floor(value)), int(np.ceil(value))
    even = below if below % 2 == 0 else above
    print(f"{value:>8}{got:>12.0f}{'down' if got < value else 'up':>12}"
          f"{got - value:>+9.1f}{f'{even} is the even neighbour':>28}")

gap_start, gap_stop = (0.5, 0.0), (3.5, 3.0)
gap_samples = np.linspace(gap_start, gap_stop, 4)
gap_rows, gap_cols = gap_samples[:, 0], gap_samples[:, 1]

to_even = [(int(np.round(r)), int(c)) for r, c in gap_samples]
to_up = [(int(np.floor(r + 0.5)), int(c)) for r, c in gap_samples]
SHAPE_ROWS, SHAPE_COLS = 5, 4
SHAPE = (SHAPE_ROWS, SHAPE_COLS)

drawn_rows = {row for row, _ in to_even}
skipped_rows = [r for r in range(SHAPE_ROWS) if r not in drawn_rows]

fig, axes = plt.subplots(1, 3, figsize=(9.6, 2.9))

pixel_axes(axes[0], SHAPE, "every sample lands on a row boundary")
for boundary in (0.5, 1.5, 2.5, 3.5):
    axes[0].axhline(boundary, color=C_LINE, lw=1.4, ls="--", zorder=2)
exact(axes[0], gap_start, gap_stop, color=INK)
axes[0].plot(gap_cols, gap_rows, "o", color=INK, markersize=5, zorder=6)
for row, col in gap_samples:
    axes[0].annotate(f"row {row:g}", (col, row), textcoords="offset points",
                     xytext=(8, 9), fontsize=7, color=C_LINE, zorder=7,
                     bbox=dict(boxstyle="round,pad=0.15", fc="white", ec="none"))

pixel_axes(axes[1], SHAPE, "half to even: rows 0, 2, 2, 4")
fill(axes[1], to_even, C_ND)
for row in skipped_rows:
    axes[1].add_patch(Rectangle((-0.5, row - 0.5), SHAPE_COLS, 1, facecolor="none",
                                edgecolor=C_LINE, lw=1.6, ls=":", zorder=5))
exact(axes[1], gap_start, gap_stop, color=INK)

pixel_axes(axes[2], SHAPE, "half up: rows 1, 2, 3, 4")
fill(axes[2], to_up, C_ND)
exact(axes[2], gap_start, gap_stop, color=INK)

fig.suptitle("dotted band: a row the line crosses and no pixel is drawn in",
             y=1.03)
fig.tight_layout()

from skimage.measure import label

for name, pixels in (("half to even", to_even), ("half up", to_up)):
    canvas = np.zeros(SHAPE, int)
    for row, col in pixels:
        canvas[row, col] = 1
    print(f"   {name:<14}{label(canvas, connectivity=2).max()} connected "
          f"component(s) for {len(pixels)} pixels")

print(f"{'samples on one axis':<44}{'rounded':<22}{'worst step':>11}")
for label, seq in (
    ("spacing 1, starting on a half", np.arange(0.5, 4.5, 1.0)),
    ("spacing 1/2, starting on a half", np.arange(0.5, 3.0, 0.5)),
    ("spacing 1, starting on an integer", np.arange(0.0, 4.0, 1.0)),
    ("spacing 3/4, starting on a half", 0.5 + 0.75 * np.arange(5)),
):
    rounded = np.round(seq).astype(int)
    print(f"   {label:<41}{str(rounded):<22}{np.abs(np.diff(rounded)).max():>11}")

def guard_fires(first, spacing=1.0, n=5):
    """The `_round_safe` test, transcribed: a half fraction and unit spacing."""
    coords = first + np.arange(n) * spacing
    return bool(coords[0] % 1 == 0.5 and coords[1] - coords[0] == 1)


print("spacing 1, varying the first coordinate")
for first in (0.0, 0.25, 0.5, 1.5, 3.5, 7.5, -0.5, -2.5):
    print(f"   first = {first:>5}  -> {'floor' if guard_fires(first) else 'round'}")

print("\nfirst coordinate 3.5, varying the spacing")
for spacing in (1.0, 0.999999, 0.5, 2.0):
    print(f"   spacing = {spacing:<9} -> "
          f"{'floor' if guard_fires(3.5, spacing) else 'round'}")

def plain_round(start, stop):
    """What `line_nd` would give if it always used `np.round`."""
    npoints = int(np.ceil(np.max(np.abs(np.subtract(stop, start))))) + 1
    samples = np.linspace(start, stop, npoints, endpoint=True).T
    return [tuple(int(x) for x in c) for c in np.round(samples).astype(int).T]


def nd_pixels(start, stop):
    return [tuple(int(x) for x in c)
            for c in zip(*line_nd(start, stop, endpoint=True))]


distinct = [(a, b) for a, b in all_pairs if a != b]
same = sum(nd_pixels(a, b) == plain_round(a, b) for a, b in distinct)
print(f"integer endpoints: line_nd equals plain np.round on "
      f"{same}/{len(distinct)} pairs  ({same / len(distinct):.1%})")

half_start, half_stop = (0.5, 0.0), (4.5, 2.0)
print(f"from {half_start} to {half_stop}")
print(f"   line_nd    {nd_pixels(half_start, half_stop)}")
print(f"   np.round   {plain_round(half_start, half_stop)}")

ends = ((3.5, 0.0), (0.5, 3.0))

for name, (a, b) in (("named down", ends), ("named up", ends[::-1])):
    pixels = nd_pixels(a, b)
    steps = [int(np.max(np.abs(np.subtract(t, s))))
             for s, t in zip(pixels, pixels[1:])]
    print(f"{name}:  {a} -> {b}")
    print(f"   samples on axis 0   {np.linspace(a, b, 4)[:, 0]}")
    print(f"   line_nd             {pixels}")
    print(f"   steps between them  {steps}")

fig, axes = plt.subplots(1, 2, figsize=(6.0, 2.6))
for ax, (a, b) in zip(axes, (ends, ends[::-1])):
    pixel_axes(ax, (5, 4), f"{a} to {b}")
    fill(ax, sk_nd(a, b), C_ND)
    exact(ax, a, b, color=INK)
fig.suptitle("line_nd on one segment, named each way round", y=1.02)
fig.tight_layout()

rng = np.random.default_rng(0)
counts = {"every axis counts up": [0, 0], "some axis counts down": [0, 0]}

for _ in range(10000):
    a = rng.integers(0, 20, 2) + 0.5
    b = rng.integers(0, 20, 2) + 0.5
    if np.array_equal(a, b):
        continue
    key = "every axis counts up" if np.all(b >= a) else "some axis counts down"
    pixels = nd_pixels(a, b)
    counts[key][0] += any(np.max(np.abs(np.subtract(t, s))) > 1
                          for s, t in zip(pixels, pixels[1:]))
    counts[key][1] += 1

print("random half-integer endpoints in a 20x20 box")
for key, (bad, total) in counts.items():
    print(f"   {key:<24}{bad:>6} gapped /{total:>6}")

samples = np.linspace(p, q, 5, endpoint=True)

fig, ax = plt.subplots(figsize=(3.6, 2.0))
pixel_axes(ax, shape)
fill(ax, sk_nd(p, q), C_ND)
exact(ax, p, q, color="white")
ax.plot(samples[:, 1], samples[:, 0], "o", color=INK, markersize=5, zorder=5)
for i, j in samples:
    ax.annotate(
        f"{i:g}",
        (j, i),
        textcoords="offset points",
        xytext=(0, 9),
        ha="center",
        fontsize=7,
        color=INK,
    )
ax.set_title("line_nd((0, 0), (1, 4)) with its sample points, labelled by exact row")
fig.tight_layout()

print(f"{'column':>7}{'exact row':>11}{'line':>7}{'line_nd':>9}")
lr = dict(zip(*[a.tolist() for a in line(0, 0, 1, 4)][::-1]))
nr = dict(zip(*[a.tolist() for a in line_nd(p, q, endpoint=True)][::-1]))
for j in range(5):
    print(f"{j:>7}{0 + j * 0.25:>11.2f}{lr[j]:>7}{nr[j]:>9}")

fig, axes = plt.subplots(1, 2, figsize=(7.2, 2.1))
for ax, pix, color, name in (
    (axes[0], sk_line(p, q), C_LINE, "line"),
    (axes[1], sk_nd(p, q), C_ND, "line_nd"),
):
    pixel_axes(ax, shape, name)
    fill(ax, pix, color)
    exact(ax, p, q, color="white")
fig.suptitle("at the tie, line steps early and line_nd steps late", y=1.02)
fig.tight_layout()

n = 13
agree = np.zeros((n, n), dtype=int)
for i in range(n):
    for j in range(n):
        agree[i, j] = sk_line((0, 0), (i, j)) == sk_nd((0, 0), (i, j))

fig, ax = plt.subplots(figsize=(3.4, 3.4))
ax.imshow(agree, cmap=ListedColormap([C_ND, C_OFF]), vmin=0, vmax=1)
ax.set_xticks(range(0, n, 2))
ax.set_yticks(range(0, n, 2))
ax.tick_params(length=0, labelsize=7, colors=MUTED)
for spine in ax.spines.values():
    spine.set_visible(False)
ax.set_xlabel("end, axis 1", fontsize=8, color=MUTED)
ax.set_ylabel("end, axis 0", fontsize=8, color=MUTED)
ax.set_title("start fixed at (0, 0)")
fig.legend(
    handles=[
        Patch(facecolor=C_OFF, label="same pixels"),
        Patch(facecolor=C_ND, label="different pixels"),
    ],
    loc="lower center",
    ncols=2,
    frameon=False,
    fontsize=8,
)
fig.tight_layout(rect=(0, 0.08, 1, 1))

# Every ordered pair of endpoints in a 9x9 box, kept clear of the canvas edge
# so that Pillow and OpenCV have room to draw in section 7.
LO, HI = 4, 13
pts = list(itertools.product(range(LO, HI), repeat=2))
box = [(a, b) for a in pts for b in pts]

diff = sum(sk_line(a, b) != sk_nd(a, b) for a, b in box)
print(f"{len(box)} endpoint pairs, {diff} differ  ({diff / len(box):.1%})")

# The rate depends on the box: longer segments have more chances to hit a tie.
wide = range(-6, 7)
wide_pairs = [((a, b), (c, d)) for a, b, c, d in itertools.product(wide, repeat=4)]
wdiff = sum(sk_line(a, b) != sk_nd(a, b) for a, b in wide_pairs)
print(f"{len(wide_pairs)} pairs in a 13x13 box, {wdiff} differ  ({wdiff / len(wide_pairs):.1%})")

fig, axes = plt.subplots(2, 4, figsize=(9.6, 4.0))
for col, t in enumerate(range(4)):
    a, b = (t, 0), (1 + t, 4)
    for row, (f, color, name) in enumerate(
        ((sk_line, C_LINE, "line"), (sk_nd, C_ND, "line_nd"))
    ):
        ax = axes[row, col]
        pixel_axes(ax, (6, 6), f"{name}, shifted by {t}")
        fill(ax, f(a, b), color)
        exact(ax, a, b, color="white")
fig.suptitle(
    "shift the same segment down one row at a time: line keeps its shape, "
    "line_nd does not",
    y=1.01,
)
fig.tight_layout()

fig, axes = plt.subplots(1, 4, figsize=(9.6, 2.1))
panels = [
    (sk_line, p, q, C_LINE, "line, start to stop"),
    (sk_line, q, p, C_LINE, "line, stop to start"),
    (sk_nd, p, q, C_ND, "line_nd, start to stop"),
    (sk_nd, q, p, C_ND, "line_nd, stop to start"),
]
for ax, (f, a, b, color, name) in zip(axes, panels):
    pixel_axes(ax, shape, name)
    fill(ax, f(a, b), color)
    exact(ax, a, b, color="white")
fig.suptitle("line changes when you swap the ends; line_nd does not", y=1.04)
fig.tight_layout()

def symmetry(f, pairs, kind, shift=3):
    ok = 0
    for a, b in pairs:
        if kind == "reversal":
            ok += f(a, b) == f(b, a)
        else:
            moved = {(i + shift, j + shift) for i, j in f(a, b)}
            at = (a[0] + shift, a[1] + shift), (b[0] + shift, b[1] + shift)
            ok += moved == f(*at)
    return ok / len(pairs)


print(f"{'':16}{'reversal':>12}{'translation':>14}")
for name, f in (("line", sk_line), ("line_nd", sk_nd)):
    print(
        f"{name:<16}{symmetry(f, box, 'reversal'):>11.1%}"
        f"{symmetry(f, box, 'translation'):>13.1%}"
    )

CANVAS = 24


def pil_line(a, b):
    im = Image.new("L", (CANVAS, CANVAS), 0)
    ImageDraw.Draw(im).line([(a[1], a[0]), (b[1], b[0])], fill=255, width=1)
    return set(map(tuple, np.argwhere(np.array(im))))


def cv_line(a, b, connectivity=8):
    arr = np.zeros((CANVAS, CANVAS), np.uint8)
    cv2.line(arr, (a[1], a[0]), (b[1], b[0]), 255, 1, lineType=connectivity)
    return {(int(i), int(j)) for i, j in np.argwhere(arr)}

off = (4, 4)  # keep the segment inside the canvas
pa, qa = (off[0] + p[0], off[1] + p[1]), (off[0] + q[0], off[1] + q[1])

panels = [("skimage line", sk_line(pa, qa), C_LINE), ("pillow", pil_line(pa, qa), C_BOTH)]
panels.append(("opencv, 8-connected", cv_line(pa, qa), C_BOTH))
panels.append(("skimage line_nd", sk_nd(pa, qa), C_ND))

fig, axes = plt.subplots(1, len(panels), figsize=(2.4 * len(panels), 2.2))
for ax, (name, pix, color) in zip(axes, panels):
    sub = {(i - off[0] + 0, j - off[1] + 0) for i, j in pix}
    pixel_axes(ax, shape, name)
    fill(ax, sub, color)
    exact(ax, p, q, color="white")
fig.suptitle("the same tie case in four rasterisers", y=1.04)
fig.tight_layout()

fns = {
    "sk.line": sk_line,
    "sk.line_nd": sk_nd,
    "pillow": pil_line,
    "opencv8": cv_line,
    "opencv4": lambda a, b: cv_line(a, b, 4),
}

drawn = {k: [f(a, b) for a, b in box] for k, f in fns.items()}
names = list(fns)
print(f"agreement over {len(box)} endpoint pairs")
print(f"{'':12}" + "".join(f"{n:>11}" for n in names))
for a in names:
    row = "".join(
        f"{sum(x == y for x, y in zip(drawn[a], drawn[b])) / len(box):>10.1%} "
        for b in names
    )
    print(f"{a:<12}{row}")

# Every eleventh pair: Pillow and OpenCV each build an image per call, so the
# full corpus is slow here. The skimage figures above use all 6561 pairs.
sample = box[::11]
print(f"{'':14}{'reversal':>12}{'translation':>14}")
for name, f in fns.items():
    print(
        f"{name:<14}{symmetry(f, sample, 'reversal'):>11.1%}"
        f"{symmetry(f, sample, 'translation'):>13.1%}"
    )

def column_decreasing(a, b):
    """Order the endpoints so the column index decreases."""
    return (a, b) if b[1] <= a[1] else (b, a)


pillow_same = sum(pil_line(a, b) == sk_line(a, b) for a, b in box)
opencv_same = sum(cv_line(a, b) == sk_line(*column_decreasing(a, b)) for a, b in box)
print(f"over {len(box)} endpoint pairs")
print(f"   pillow  == sk.line, same endpoint order      : {pillow_same / len(box):.1%}")
print(f"   opencv8 == sk.line, column decreasing        : {opencv_same / len(box):.1%}")

from collections import Counter

pattern = Counter()
for a, b in box:
    if sk_line(a, b) == sk_line(b, a):
        continue                      # direction makes no difference here
    drawn_by_cv = cv_line(a, b)
    match = "forward" if drawn_by_cv == sk_line(a, b) else "reversed"
    d_col = b[1] - a[1]
    pattern[(match, "column increases" if d_col > 0 else "column decreases")] += 1

for (match, direction), n in sorted(pattern.items()):
    print(f"   opencv matched our {match:<8} run when the {direction}: {n:>5}")

from skimage.draw import line_aa

AA = 16
aa_start, aa_stop = (3, 2), (9, 13)

pillow_image = Image.new("L", (AA, AA), 0)
ImageDraw.Draw(pillow_image).line(
    [(aa_start[1], aa_start[0]), (aa_stop[1], aa_stop[0])], fill=255, width=1
)
pillow_aa = np.array(pillow_image) / 255

opencv_aa = np.zeros((AA, AA), np.uint8)
cv2.line(opencv_aa, (aa_start[1], aa_start[0]), (aa_stop[1], aa_stop[0]),
         255, 1, lineType=cv2.LINE_AA)
opencv_aa = opencv_aa / 255

skimage_aa = np.zeros((AA, AA))
rr, cc, value = line_aa(aa_start[0], aa_start[1], aa_stop[0], aa_stop[1])
skimage_aa[rr, cc] = value

for name, img in (("skimage line_aa", skimage_aa), ("opencv LINE_AA", opencv_aa),
                  ("pillow", pillow_aa)):
    print(f"{name:<18}{len(np.unique(img)):>3} distinct levels,"
          f"{int((img > 0).sum()):>4} pixels touched")

fig, axes = plt.subplots(1, 3, figsize=(9.0, 2.8))
for ax, (name, img) in zip(axes, (("skimage line_aa", skimage_aa),
                                  ("opencv LINE_AA", opencv_aa),
                                  ("pillow, no anti-aliasing", pillow_aa))):
    ax.imshow(img, cmap=LinearSegmentedColormap.from_list("c", [C_OFF, C_LINE]),
              vmin=0, vmax=1, interpolation="nearest")
    ax.plot([aa_start[1], aa_stop[1]], [aa_start[0], aa_stop[0]],
            color=INK, linewidth=1.0)
    ax.set_xticks([])
    ax.set_yticks([])
    for spine in ax.spines.values():
        spine.set_visible(False)
    ax.set_title(name)
fig.suptitle("coverage, not just which pixels", y=1.04)
fig.tight_layout()

import scipy.ndimage as ndi

# Slot 3 of the reference categorical palette; validates all-pairs with the
# blue and orange already in use.
C_FOUR = "#1baf7a"

S4 = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
S8 = np.ones((3, 3), int)

c8 = c4 = 0
for a, b in box:
    di, dj = abs(a[0] - b[0]), abs(a[1] - b[1])
    c8 += len(cv_line(a, b, 8)) == max(di, dj) + 1
    c4 += len(cv_line(a, b, 4)) == di + dj + 1
print(f"over {len(box)} endpoint pairs")
print(f"   8-connected count == max(|di|, |dj|) + 1 : {c8 / len(box):.1%}")
print(f"   4-connected count == |di| + |dj| + 1     : {c4 / len(box):.1%}")

a, b = (4, 4), (7, 13)
fig, axes = plt.subplots(1, 2, figsize=(7.2, 2.4))
for ax, conn, color, name in (
    (axes[0], 8, C_BOTH, "LINE_8, 8-connected"),
    (axes[1], 4, C_FOUR, "LINE_4, 4-connected"),
):
    pixel_axes(ax, (5, 12), f"{name}  ({len(cv_line(a, b, conn))} pixels)")
    fill(ax, {(i - 3, j - 3) for i, j in cv_line(a, b, conn)}, color)
    exact(ax, (a[0] - 3, a[1] - 3), (b[0] - 3, b[1] - 3), color="white")
fig.suptitle("the same segment, drawn twice", y=1.04)
fig.tight_layout()

def path_order(pixels, start, stop):
    """The pixels in the order the line visits them, major axis first."""
    delta = np.abs(np.subtract(stop, start))
    major = int(np.argmax(delta))
    minor = 1 - major
    sign = np.sign(np.subtract(stop, start))
    return sorted(pixels,
                  key=lambda ij: (sign[major] * ij[major], sign[minor] * ij[minor]))


def steps_taken(pixels, start, stop):
    """The distinct moves between consecutive pixels along the path."""
    walk = path_order(pixels, start, stop)
    return sorted({(abs(u[0] - v[0]), abs(u[1] - v[1]))
                   for u, v in zip(walk, walk[1:])})


for conn in (8, 4):
    pix = cv_line(a, b, conn)
    print(f"LINE_{conn}: {len(pix):>3} pixels, steps {steps_taken(pix, a, b)}")

diagonal_in_4 = sum((1, 1) in steps_taken(cv_line(u, v, 4), u, v) for u, v in box)
print(f"LINE_4 lines containing a (1, 1) move, out of {len(box)}: {diagonal_in_4}")

def as_mask(pixels, shape=(CANVAS, CANVAS)):
    """The pixel set as a boolean image."""
    m = np.zeros(shape, bool)
    for i, j in pixels:
        m[i, j] = True
    return m


print(f"{'':10}{'one component under S4':>26}{'under S8':>12}")
for conn in (8, 4):
    n4 = sum(ndi.label(as_mask(cv_line(a, b, conn)), structure=S4)[1] == 1
             for a, b in box)
    n8 = sum(ndi.label(as_mask(cv_line(a, b, conn)), structure=S8)[1] == 1
             for a, b in box)
    print(f"LINE_{conn:<5}{n4 / len(box):>25.1%}{n8 / len(box):>12.1%}")

axis_aligned = sum(a[0] == b[0] or a[1] == b[1] for a, b in box)
print(f"\npairs with no diagonal step at all: {axis_aligned / len(box):>10.1%}")

nested = sum(cv_line(a, b, 8) <= cv_line(a, b, 4) for a, b in box)
print(f"8-connected set contained in the 4-connected set: {nested / len(box):.1%}")

example = next((a, b) for a, b in box if not cv_line(a, b, 8) <= cv_line(a, b, 4))
s8, s4 = cv_line(*example, 8), cv_line(*example, 4)
print(f"\nfirst counter-example: {example[0]} to {example[1]}")
print(f"   only in LINE_8: {sorted(s8 - s4)}")
print(f"   only in LINE_4: {sorted(s4 - s8)}")

a, b = example
lo = (min(a[0], b[0]) - 1, min(a[1], b[1]) - 1)
span = (abs(a[0] - b[0]) + 3, abs(a[1] - b[1]) + 3)
fig, axes = plt.subplots(1, 2, figsize=(7.2, 2.2))
for ax, pix, color, name in (
    (axes[0], s8, C_BOTH, "LINE_8"), (axes[1], s4, C_FOUR, "LINE_4")
):
    pixel_axes(ax, span, name)
    fill(ax, {(i - lo[0], j - lo[1]) for i, j in pix}, color)
    exact(ax, (a[0] - lo[0], a[1] - lo[1]), (b[0] - lo[0], b[1] - lo[1]),
          color="white")
fig.suptitle("neither set contains the other", y=1.06)
fig.tight_layout()

def barrier(conn, shape=(28, 40), a=(2, 0), b=(25, 39)):
    """Draw a line across the array, then label what it leaves free."""
    arr = np.zeros(shape, np.uint8)
    cv2.line(arr, (a[1], a[0]), (b[1], b[0]), 255, 1, lineType=conn)
    labels, n = ndi.label(~(arr > 0), structure=S8)
    return arr > 0, labels, n


for conn in (8, 4):
    _, _, n = barrier(conn)
    verdict = "sealed" if n > 1 else "the fill leaks through"
    print(f"LINE_{conn}: the free area is {n} component(s)  -> {verdict}")

fig, axes = plt.subplots(2, 1, figsize=(5.4, 4.6))
for ax, conn in zip(axes, (8, 4)):
    m, labels, n = barrier(conn)
    ax.imshow(np.where(m, 0, labels),
              cmap=ListedColormap([C_FOUR, C_OFF, C_LINE]), vmin=0, vmax=2)
    ax.set_xticks([])
    ax.set_yticks([])
    for spine in ax.spines.values():
        spine.set_visible(False)
    ax.set_title(f"LINE_{conn}: {n} free component(s)")
fig.legend(
    handles=[
        Patch(facecolor=C_FOUR, label="the line"),
        Patch(facecolor=C_OFF, label="one side"),
        Patch(facecolor=C_LINE, label="the other side"),
    ],
    loc="lower center", ncols=3, frameon=False, fontsize=8,
)
fig.suptitle("only the 4-connected line divides the array in two", y=1.0)
fig.tight_layout(rect=(0, 0.07, 1, 1))

def line_4(start, stop):
    """Bresenham, with a corner pixel inserted at each diagonal step."""
    ii, jj = line(start[0], start[1], stop[0], stop[1])
    origin = np.asarray(start, float)
    direction = np.asarray(stop, float) - origin

    def offset(pixel):
        v = np.asarray(pixel, float) - origin
        return abs(v[0] * direction[1] - v[1] * direction[0])

    out = [(int(ii[0]), int(jj[0]))]
    for i, j in zip(ii[1:].tolist(), jj[1:].tolist()):
        pi, pj = out[-1]
        if i != pi and j != pj:
            out.append(min(((pi, j), (i, pj)), key=offset))
        out.append((i, j))
    return out

conn_ok = count_ok = 0
for a, b in box:
    pix = line_4(a, b)
    conn_ok += ndi.label(as_mask(pix), structure=S4)[1] == 1
    count_ok += len(set(pix)) == abs(a[0] - b[0]) + abs(a[1] - b[1]) + 1
print(f"over {len(box)} endpoint pairs")
print(f"   4-connected           {conn_ok / len(box):.1%}")
print(f"   Manhattan pixel count {count_ok / len(box):.1%}")

arr = np.zeros((28, 40), np.uint8)
for i, j in line_4((2, 0), (25, 39)):
    arr[i, j] = 255
print(f"   seals the array       {ndi.label(~(arr > 0), structure=S8)[1] > 1}")
same = sum(set(line_4(a, b)) == cv_line(a, b, 4) for a, b in box)
print(f"   same pixels as LINE_4 {same / len(box):.1%}")

a, b = (4, 4), (7, 13)
fig, axes = plt.subplots(1, 2, figsize=(7.2, 2.4))
for ax, pix, color, name in (
    (axes[0], set(line_4(a, b)), C_ND, "candidate line_4"),
    (axes[1], cv_line(a, b, 4), C_FOUR, "opencv LINE_4"),
):
    pixel_axes(ax, (5, 12), name)
    fill(ax, {(i - 3, j - 3) for i, j in pix}, color)
    exact(ax, (a[0] - 3, a[1] - 3), (b[0] - 3, b[1] - 3), color="white")
fig.suptitle("same guarantees, different corners", y=1.04)
fig.tight_layout()

def bresenham_nd(start, stop):
    """Bresenham in any number of dimensions: one error term per minor axis."""
    start, stop = np.array(start), np.array(stop)
    delta = np.abs(stop - start)
    step = np.sign(stop - start)

    major = int(np.argmax(delta))
    n_steps = int(delta[major])
    error = 2 * delta - n_steps          # one accumulator per axis

    at = start.copy()
    out = []
    for _ in range(n_steps):
        out.append(tuple(int(x) for x in at))
        for axis in range(len(delta)):
            if axis == major:
                continue
            if error[axis] >= 0:
                at[axis] += step[axis]
                error[axis] -= 2 * n_steps
            error[axis] += 2 * delta[axis]
        at[major] += step[major]

    out.append(tuple(int(x) for x in stop))
    return out

matches = sum(bresenham_nd(a, b) == sk_sequence(a, b) for a, b in all_pairs)
print(f"identical to skimage.draw.line in 2-D on {matches}/{len(all_pairs)} pairs"
      f"  ({matches / len(all_pairs):.1%})")

def nd_sequence(start, stop):
    """`line_nd` as a list of integer tuples, endpoint included."""
    return [tuple(int(x) for x in c)
            for c in zip(*line_nd(start, stop, endpoint=True))]


rng = np.random.default_rng(0)
grid = list(itertools.product(range(-4, 5), repeat=3))
choice = rng.choice(len(grid), size=(1500, 2))
triples = [(grid[a], grid[b]) for a, b in choice if grid[a] != grid[b]]

chebyshev = connected = 0
for start, stop in triples:
    pix = bresenham_nd(start, stop)
    chebyshev += len(pix) == max(abs(np.subtract(stop, start))) + 1
    connected += all(max(abs(np.subtract(v, u))) <= 1
                     for u, v in zip(pix, pix[1:]))

agree = sum(bresenham_nd(a, b) == nd_sequence(a, b) for a, b in triples)
print(f"over {len(triples)} random 3-D segments")
print(f"   point count is max(abs(delta)) + 1   : {chebyshev / len(triples):.1%}")
print(f"   every step moves at most 1 per axis  : {connected / len(triples):.1%}")
print(f"   identical to line_nd                 : {agree / len(triples):.1%}")

def holds(fn, kind, cases):
    ok = 0
    for start, stop in cases:
        if kind == "reversal":
            ok += set(fn(start, stop)) == set(fn(stop, start))
        else:
            shift = np.array([3, 1, 2])
            moved = {tuple(np.add(c, shift)) for c in fn(start, stop)}
            ok += moved == set(fn(tuple(np.add(start, shift)),
                                  tuple(np.add(stop, shift))))
    return ok / len(cases)


sample = triples[:600]
print(f"{'':16}{'reversal':>11}{'translation':>14}")
for name, fn in (("bresenham_nd", bresenham_nd), ("line_nd", nd_sequence)):
    print(f"{name:<16}{holds(fn, 'reversal', sample):>10.1%}"
          f"{holds(fn, 'translation', sample):>14.1%}")

# All three axes vary, and the two algorithms disagree on part of the run.
start3, stop3 = (0, 0, 0), (2, 4, 8)
walks = {"bresenham_nd": bresenham_nd(start3, stop3),
         "line_nd": nd_sequence(start3, stop3)}

fig, axes = plt.subplots(2, 2, figsize=(8.4, 4.0))
for row, (name, walk) in enumerate(walks.items()):
    colour = C_LINE if name == "bresenham_nd" else C_ND
    for col, axis in enumerate((0, 1)):
        ax = axes[row, col]
        shared = {(v[axis], v[2]) for v in walks["bresenham_nd"]} & \
                 {(v[axis], v[2]) for v in walks["line_nd"]}
        mine = {(v[axis], v[2]) for v in walk}
        pixel_axes(ax, (5 if axis else 3, 9),
                   f"{name}, axis {axis} against axis 2")
        fill(ax, shared, C_BOTH)
        fill(ax, mine - shared, colour)
        exact(ax, (start3[axis], start3[2]), (stop3[axis], stop3[2]),
              color="white")
fig.suptitle("grey where the projections coincide, coloured where they do not", y=1.02)
fig.tight_layout()

differ = [(u, v) for u, v in zip(walks["bresenham_nd"], walks["line_nd"]) if u != v]
print(f"{len(differ)} of {len(walks['bresenham_nd'])} positions differ:")
for u, v in differ:
    print(f"   bresenham_nd {u}   line_nd {v}")

print("line_nd on the true float segment:")
print("  ", nd_sequence((0.4, 0.2), (3.6, 9.8)))
print("bresenham on the rounded endpoints:")
print("  ", bresenham_nd((0, 0), (4, 10)))

spaced = np.array([0.5, 1.25, 2.0, 2.75, 3.5])   # the example from the PR
for name, rule in (("np.round", np.round),
                   ("np.floor", np.floor),
                   ("floor(x) + 0.5, as written in the PR",
                    lambda a: np.floor(a) + 0.5),
                   ("floor(x + 0.5), half up", lambda a: np.floor(a + 0.5))):
    got = rule(spaced).astype(int)
    print(f"   {name:<38}{str(got):<18}steps {np.diff(got)}")

def line_nd_halfup_pixels(a, b):
    """Ordered pixels from `line_nd` with half-up rounding."""
    a, b = np.asarray(a, float), np.asarray(b, float)
    npoints = int(np.ceil(np.max(np.abs(b - a)))) + 1
    coords = np.floor(np.linspace(a, b, npoints, endpoint=True).T + 0.5).astype(int)
    return [tuple(int(v) for v in c) for c in coords.T]


print("the PR's spacing, drawn as an actual line")
for a, b in (((0.5, 0.0), (3.5, 4.0)), ((0, 0), (3, 4))):
    now = nd_pixels(a, b)
    up = line_nd_halfup_pixels(a, b)
    hop = lambda pix: [int(np.max(np.abs(np.subtract(y, x))))
                       for x, y in zip(pix, pix[1:])]
    print(f"   {a} -> {b}")
    print(f"      line_nd now   {now}  steps {hop(now)}")
    print(f"      half up       {up}  steps {hop(up)}")

rules = {
    "half-up": lambda x: np.floor(x + 0.5),
    "half-to-even": np.round,
    "half-away": lambda x: np.sign(x) * np.floor(np.abs(x) + 0.5),
}

fine = np.arange(-4, 4, 1 / 512)
print(f"{'rule':<14}{'error range':>22}{'worst jump, step 1':>22}")
for name, rule in rules.items():
    errors = rule(fine) - fine
    steps = np.arange(-4, 4, 0.25)
    jump = max(abs(float(rule(np.float64(a + 1.0)) - rule(np.float64(a))))
               for a in steps)
    print(f"{name:<14}[{errors.min():+.3f}, {errors.max():+.3f}]{jump:>21.0f}")

def line_nd_halfup(a, b):
    """`line_nd` with half-up rounding instead of half-to-even."""
    a, b = np.asarray(a, float), np.asarray(b, float)
    npoints = int(np.ceil(np.max(np.abs(b - a)))) + 1
    coords = np.floor(np.linspace(a, b, npoints, endpoint=True).T + 0.5).astype(int)
    return set(zip(*(c.tolist() for c in coords)))


fig, axes = plt.subplots(2, 4, figsize=(9.6, 4.0))
for col, t in enumerate(range(4)):
    a, b = (t, 0), (1 + t, 4)
    for row, (f, name) in enumerate(
        ((sk_nd, "line_nd, now"), (line_nd_halfup, "line_nd, half up"))
    ):
        ax = axes[row, col]
        pixel_axes(ax, (6, 6), f"{name} (+{t})")
        fill(ax, f(a, b), C_ND)
        exact(ax, a, b, color="white")
fig.suptitle("half-up rounding makes line_nd keep its shape when it moves", y=1.01)
fig.tight_layout()

def line_sorted(a, b):
    """`line` with the endpoint order normalised."""
    a, b = sorted([tuple(a), tuple(b)])
    return sk_line(a, b)


fig, axes = plt.subplots(1, 4, figsize=(9.6, 2.1))
panels = [
    (sk_line, p, q, "line, forward"),
    (sk_line, q, p, "line, backward"),
    (line_sorted, p, q, "sorted, forward"),
    (line_sorted, q, p, "sorted, backward"),
]
for ax, (f, a, b, name) in zip(axes, panels):
    pixel_axes(ax, shape, name)
    fill(ax, f(a, b), C_LINE)
    exact(ax, a, b, color="white")
fig.suptitle("normalising the endpoints makes line direction-independent", y=1.04)
fig.tight_layout()

changed_nd = sum(sk_nd(a, b) != line_nd_halfup(a, b) for a, b in box)
changed_ln = sum(sk_line(a, b) != line_sorted(a, b) for a, b in box)
agree_now = sum(sk_line(a, b) == sk_nd(a, b) for a, b in box)
agree_fix = sum(line_sorted(a, b) == line_nd_halfup(a, b) for a, b in box)
n_box = len(box)

print(f"line_nd output changes under half-up rounding : {changed_nd / n_box:6.1%}")
print(f"line output changes under sorted endpoints    : {changed_ln / n_box:6.1%}")
print(f"line and line_nd agree, now                   : {agree_now / n_box:6.1%}")
print(f"line and line_nd agree, both fixed            : {agree_fix / n_box:6.1%}")

for name, f in (("line, sorted", line_sorted), ("line_nd, half up", line_nd_halfup)):
    print(
        f"\n{name}: reversal {symmetry(f, box, 'reversal'):.1%},"
        f" translation {symmetry(f, box, 'translation'):.1%}"
    )

resid = [(a, b) for a, b in box if line_sorted(a, b) != line_nd_halfup(a, b)]
a, b = resid[0]
lo = (min(a[0], b[0]) - 1, min(a[1], b[1]) - 1)

fig, axes = plt.subplots(1, 2, figsize=(7.2, 2.4))
span = (abs(a[0] - b[0]) + 3, abs(a[1] - b[1]) + 3)
for ax, (f, color, name) in zip(
    axes,
    ((line_sorted, C_LINE, "line, sorted"), (line_nd_halfup, C_ND, "line_nd, half up")),
):
    pixel_axes(ax, span, name)
    fill(ax, {(i - lo[0], j - lo[1]) for i, j in f(a, b)}, color)
    exact(ax, (a[0] - lo[0], a[1] - lo[1]), (b[0] - lo[0], b[1] - lo[1]), color="white")
fig.suptitle(f"both fixed, still different: {a} to {b}", y=1.04)
fig.tight_layout()

import PIL
import skimage

print(f"scikit-image {skimage.__version__}")
print(f"Pillow       {PIL.__version__}")
print(f"OpenCV       {cv2.__version__}")