 codex/radial-sampling-missed-branch.py

"""Original planar sampling example; write the PNG in the current directory.

Tested with Python 3.14.4, NumPy 2.3.5 and Matplotlib 3.10.7.
Synthetic line segments, not an image or a reconstruction of anatomy.
"""

from pathlib import Path

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np


branch_start = np.array([2.1, 0.0])
branch_end = np.array([2.4, 0.7])
direction = branch_end - branch_start
assert np.isclose(np.linalg.norm(branch_start), 2.1)
assert np.isclose(np.linalg.norm(branch_end), 2.5)
assert np.dot(branch_start, direction) > 0  # Radius increases along this segment.
probe = 2.25
aa = np.dot(direction, direction)
bb = 2 * np.dot(branch_start, direction)
cc = np.dot(branch_start, branch_start) - probe**2
parameter = (-bb + np.sqrt(bb**2 - 4 * aa * cc)) / (2 * aa)
assert 0 < parameter < 1
branch_hit = branch_start + parameter * direction
assert np.isclose(np.linalg.norm(branch_hit), probe)


def intersections(radius):
    return int(0 < radius <= 4) + int(2.1 < radius <= 2.5)


coarse = np.array([1, 2, 3, 4])
assert [intersections(r) for r in coarse] == [1, 1, 1, 1]
assert intersections(probe) == 2
assert intersections(2.1) == 1  # Shared endpoint is one distinct intersection.
assert intersections(2.5) == 2

plt.rcParams.update({"font.size": 10, "figure.facecolor": "white"})
fig, axes = plt.subplots(1, 2, figsize=(9.6, 4.2), dpi=100)
fig.subplots_adjust(left=0.065, right=0.985, top=0.78, bottom=0.20, wspace=0.24)
fig.suptitle("A short branch can fall between sampled radii", fontsize=15, y=0.96)
blue, purple, green = "#1766a0", "#982871", "#237b43"
ax = axes[0]
ax.plot([0, 4], [0, 0], color=blue, lw=2.5)
ax.plot([branch_start[0], branch_end[0]], [branch_start[1], branch_end[1]],
        color=blue, lw=2.5)
angle = np.linspace(0, np.pi / 2, 301)
ax.plot(probe * np.cos(angle), probe * np.sin(angle), color=green, ls="--", lw=1.5)
ax.scatter([0], [0], color="black", s=24, zorder=4)
ax.scatter([probe, branch_hit[0]], [0, branch_hit[1]], color=purple,
           marker="D", s=33, zorder=5)
ax.text(0, -0.20, "center", ha="center", fontsize=9)
ax.text(2.0, -0.20, "r = 2.1", ha="center", fontsize=9)
ax.text(2.55, 0.75, "r = 2.5", fontsize=9)
ax.text(0.04, 0.90, "Probe radius = 2.25", color=green, transform=ax.transAxes)
ax.set(xlim=(-0.3, 4.2), ylim=(-0.35, 1.5), xlabel="x (chosen units)",
       ylabel="y (chosen units)")
ax.set_aspect("equal", adjustable="box")
ax.set_title("Two intersections at an extra probe radius", fontsize=11)
ax.spines[["top", "right"]].set_visible(False)
ax.grid(alpha=0.12)

ax = axes[1]
for left, right, value in [(1, 2.1, 1), (2.1, 2.5, 2), (2.5, 4, 1)]:
    ax.plot([left, right], [value, value], color=blue, lw=2.5)
ax.scatter([2.1, 2.5], [2, 1], facecolors="white", edgecolors=blue, s=40, zorder=4)
ax.scatter([2.1, 2.5], [1, 2], color=blue, s=26, zorder=4)
ax.scatter(coarse, np.ones_like(coarse), color="black", s=32, zorder=5,
           label="Samples at radii 1, 2, 3, 4")
ax.scatter([probe], [2], color=purple, marker="D", s=35, zorder=6,
           label="Extra sample at radius 2.25")
ax.set(xlim=(0.8, 4.2), ylim=(0.65, 2.40), xlabel="Radius r (chosen units)",
       ylabel="Distinct intersections N(r)")
ax.set_xticks([1, 2, 3, 4])
ax.set_yticks([1, 2])
ax.set_title("The coarse samples all report one intersection", fontsize=11)
ax.legend(loc="upper right", frameon=False, fontsize=8)
ax.spines[["top", "right"]].set_visible(False)
ax.grid(alpha=0.12)
fig.text(0.5, 0.055, "Synthetic geometry; the shared branch endpoint is counted once.",
         ha="center", fontsize=10)
fig.savefig(Path("radial-sampling-missed-branch.png"), facecolor="white")
plt.close(fig)