Opt Art — Ariadne

AriadneTypes

The data contract shared by every Ariadne stage: DensityField (what the image became), PointSet (where the stipples are), TourSet (what order to draw them in).

Coordinates are image pixel space: origin top-left, y increasing downward, Float64. Density is ink, not brightness — 1.0 is maximum ink. Both conventions are fixed here and assumed everywhere downstream.

This module holds no algorithms. It exists so that the stages agree on shapes without depending on each other.

DensityField(ρ, γ, mass)

The ink field. ρ[i, j] ∈ [0, 1] with row 1 at the top of the image and 1.0 meaning maximum ink (black). γ records the exponent used to produce it, and mass is Σρ — the total ink, supplied by the caller.

ρ is Float32 deliberately: it is scanned repeatedly and the precision is irrelevant for a value in [0,1] that will be sampled stochastically. Anything accumulating over it should use Float64.

A 2-D point. SVector rather than a tuple or two parallel vectors, because the relaxation's inner loop must not allocate, and Vector{SVector{2,Float64}} is stored inline.

PointSet(coords, weights, bbox)

The stipple points. coords is in image pixel space. weights is per-point Voronoi mass, empty until LloydPolish computes it. bbox is (xmin, ymin, xmax, ymax) — the image extent, not the extent of the points.

Mutable, and mutated in place by relax! and polish!: at 100k points, copying per epoch is real cost for no benefit.

TourSet(tours, closed)

Index sequences into a PointSet's coords. closed means the last point of each tour connects back to the first, which is left implicit rather than duplicated in the index list.

Int32 rather than Int: the 2-opt inner loop in Part 2 is memory-bound, and halving the index width halves its traffic.

tour_indices_valid(ts, npoints) -> Bool

Whether every index in 1:npoints appears exactly once across all tours.

This is the parent design's checked invariant. A tour set that drops or duplicates points still renders — it just renders the wrong picture, with some stipples missing and others drawn twice — so it must be asserted rather than eyeballed.

Tone

Stage S0: an image file becomes a density field.

The one subtlety worth stating loudly: sRGB is linearized before luminance is computed. Skipping that step is the most common bug in stippling code — it makes every midtone systematically wrong, because sRGB 0.5 carries about 21% of the light of sRGB 1.0, not 50%.

Output density is ink: 1.0 is black. This inverts the brightness convention used by the Truchet modules, and the inversion happens here, once.

density_field(img; γ=2.0, unsharp_σ=0.0, unsharp_amount=0.0, floor=0.0) -> DensityField

Image to ink field: linearize, re-encode, optionally sharpen, invert, raise to γ, then floor.

Luminance is computed in linear light (Rec. 709 weights are only meaningful there) and then mapped back through linear_to_srgb before inversion, so that ink tracks perceived darkness. Inverting linear light directly compresses the whole picture into heavy ink and the subject never appears; see the note on linear_to_srgb.

γ = 2.0 is the default because for a line drawing the ink is the tour, whose length grows like √n (Beardwood–Halton–Hammersley), so darkness scales with the square root of point density and the density fed to the stippler should be roughly squared to compensate. That argument is asymptotic and unverified here — treat γ as a free parameter, not a settled constant.

floor, like γ, is an artistic control, not a settled constant. Any gamma-raised density strictly below floor is snapped to exactly 0.0f0 rather than left as sparse ink. "Exactly" matters: AliasSampler refuses to draw zero-weight indices, so a hard zero — not a small nonzero value — is what turns a light region into genuinely blank paper (e.g. Bosch's Father and Son background) instead of sparse stipple. mass is computed from ρ after flooring, so it always reports the true post-floor ink total and downstream mass-conservation invariants (LloydPolish) keep holding. The default 0.0 floors nothing, so it is a strict no-op against Part 1's renders.

Unsharp masking is off by default. It helps TSP art noticeably, because points must keep a minimum spacing and so fine detail is lost unless pre-emphasized, but it is a taste decision and belongs to the caller.

linear_to_srgb(c) -> Float64

The inverse sRGB transfer function: linear light in [0,1] back to a gamma-encoded value. Exact inverse of srgb_to_linear, same knee.

Why this round trip exists — do not "optimize" it away. luminance linearizes each channel before applying the Rec. 709 weights, because those weights describe how much light each primary contributes and are only meaningful in linear light. But density_field then needs a perceptual quantity: ink coverage should track how dark a region looks, not how little light it emits. Inverting linear light instead lifts every midtone toward heavy ink — on the cameraman test image it left the 5th-percentile density at 0.27, so even blank sky received a quarter of maximum stipple density and the subject never emerged from the mesh. So we weight in linear light, then re-encode before inverting. The round trip is the point.

The c >= 1 branch is not cosmetic. 1.055 * 1.0^(1/2.4) - 0.055 evaluates to 0.9999999999999999, not 1.0, because 0.055 is not representable in binary floating point. Without the guard, a pure white pixel comes back a half-ulp short of 1 and picks up 1.1e-16 of ink — nonzero, so blank paper is no longer exactly blank. Snapping the endpoint keeps "white means the pen never goes there" an exact statement rather than an approximate one.

load_density(path; kwargs...) -> DensityField

Read an image from disk and convert it. Keyword arguments pass through to density_field.

luminance(img) -> Matrix{Float64}

Rec. 709 relative luminance in linear light, in [0,1]. Row 1 is the top of the image, as in the source raster.

srgb_to_linear(c) -> Float64

The sRGB electro-optical transfer function: gamma-encoded channel value in [0,1] to linear light. Piecewise, with a linear segment below the 0.04045 knee to avoid the infinite slope of the pure power law at zero.

AliasSampler

Walker's alias method: O(1) draws from a fixed discrete distribution, after an O(n) build.

Both the initial stippling and the relaxation sample from the density field — on the order of 10⁷ to 10⁹ draws for a full run — so the per-draw cost is worth caring about. A binary search over a cumulative distribution is O(log n) and would dominate.

AliasTable(weights)

Build the alias table. weights must be finite, non-negative, and not all zero; they are normalized internally, so any positive scale works.

Zero-weight entries are never drawn.

draw(rng, t) -> Int

One sample, as an index into the original weight vector. O(1).

draw_pixel(rng, t, nrows) -> (row, col)

One sample from a table built over a flattened matrix, converted back to (row, column). Julia stores matrices column-major, so linear index k is row mod1(k, nrows), column cld(k, nrows).

Stipple

Stage S1: the density field becomes an initial set of points.

This is density-proportional importance sampling — draw pixels in proportion to their ink, jitter within the pixel. The result is clumpy, with visible voids: a white-noise spectrum rather than the blue-noise one the drawing wants. That is expected. MacQueen fixes it; a good initialization only saves epochs.

Pixel (i, j) occupies x ∈ [j-1, j), y ∈ [i-1, i), so the image extent is (0, 0, width, height) and coordinates are directly comparable to the density field's indices.

density_table(d) -> AliasTable

An alias table over the flattened density field. Built once and reused: the relaxation samples from the same distribution millions of times.

initial_points(d, n; rng, jitter=true) -> PointSet

Draw n points with probability proportional to ink.

With jitter=true (the default) each point is placed uniformly within its pixel; with jitter=false it sits at the pixel's top-left corner, which is useful for pinning the coordinate convention in tests but produces a visibly gridded point set.

BucketGrid

A uniform-cell spatial index over a point set that mutates on every iteration.

That mutation is why this exists rather than a k-d tree: the relaxation moves one point per sample and runs on the order of 10⁸ samples, so the index must support O(1) incremental update. A static tree would have to be rebuilt.

The structure works here specifically because of what the relaxation is doing: its whole purpose is to equidistribute the points, so the grid stays balanced. It would be a poor choice for a permanently clustered set — and dark image regions are genuinely denser — so the ring scan is written to stay correct, if slower, when balance fails.

Grid(coords, bbox, cellsize)

Buckets of point indices on a uniform lattice. cells[r, c] holds the indices of the points currently inside that cell.

build_grid(ps) -> Grid

Size the cells at roughly the mean nearest-neighbour spacing, √(area/n), which puts about one point per cell.

cell_of(g, p) -> (row, col)

The cell containing p, clamped to the grid. Clamping (rather than erroring) means a point that drifts outside the image extent still has a home, which matters because Lloyd centroids can land on the boundary.

move_point!(g, idx, from, to)

Update membership after point idx moved. A move within one cell touches nothing; a move across a boundary is one deletion and one push, both O(1) expected because cells hold ~1 point.

from must be the point's last known position — the one that produced the cell it is currently registered in. There is no runtime check for this (no @assert, no @boundscheck): the function is called on the order of 10⁸ times and the hot path stays clean. Passing a stale from silently corrupts the index — the deletion lookup finds nothing in the (wrong) old cell and is skipped, but the push! into the new cell still happens, so idx ends up registered in two cells permanently. Callers that track positions themselves (e.g. MacQueen's relaxation) must pass exactly the value they last handed to this function as to.

nearest(g, coords, q) -> Int

Index of the point nearest q, or 0 if the grid is empty.

Scans the home cell, then expanding square rings. Termination is distance- based: the scan stops when the closest possible point in the next ring is farther than the best found so far. Stopping on ring index instead gives subtly wrong neighbours near cell boundaries — an error that never crashes and shows up only as structured artifacts in the finished drawing.

MacQueen

Stage S2: MacQueen's sequential k-means, which Bosch calls the "tractor beam".

Repeatedly draw a point from the image density; whichever stipple is nearest gets tugged toward it by a shrinking amount. The fixed point is a centroidal Voronoi tessellation weighted by the density — the same target as Secord's weighted Voronoi stippling, reached by stochastic approximation instead of by constructing Voronoi cells.

This is the heart of the pipeline. The initial sampling is clumpy white noise; what makes the finished drawing read as an image rather than as static is that this stage turns it into something close to blue noise.

Mean nearest-neighbour distance and its coefficient of variation, via the grid. CV falling from ≈0.5 (Poisson) toward ≈0.2 indicates blue-noise-like spacing.

learning_rate(schedule, c; α=0.7, η₀=0.01) -> Float64

Step size for a stipple that has now won c samples.

• :harmonic — 1/c. Makes the stipple the exact running mean of its samples. Provably convergent, but early samples (drawn when the partition was still wrong) carry the same weight as late ones, and the step decays so fast that a point cannot escape a bad basin.
• :power — c^(−α). Robbins–Monro conditions hold for α ∈ (0.5, 1]. The default α = 0.7 keeps points moving long enough to correct early error.
• :constant — η₀. Never converges; tracks a stationary distribution. Useful as an annealing prefix before a harmonic tail.

relax!(ps, d; epochs=150, schedule=:power, α=0.7, η₀=0.01, rng, reseed_after=5)

Relax ps in place toward a density-weighted CVT. One epoch is length(ps.coords) samples. Returns (; mean_displacement, max_displacement, mean_nn, cv_nn, reseeded) for the final epoch.

reseed_after is a warm-up: from that epoch onward, any point that has still won no samples at all is treated as dead and re-seeded by drawing a fresh sample from the density. This keeps n exactly constant and prevents stranded stipples in blown-out highlights.

The deadness test is cumulative ("has never won"), not per-epoch ("missed the last few epochs"). One epoch is only n samples for n points, so by balls-in-bins occupancy a perfectly healthy point wins nothing in a given epoch with probability (1-1/n)^n → 1/e ≈ 0.37. A rule that reseeds after k consecutive idle epochs therefore fires on healthy points at a rate of about n·e^{-k} per epoch — roughly 3 points per epoch at n=400, k=5 — and each such "reseed" teleports a well-placed point away and re-inflates the learning rate of the point it split. Measured, that drives the nearest-neighbour CV up with epochs (0.32 → 0.55 over 400 epochs) instead of down. A point that has never won, by contrast, is genuinely in a region with no ink, and the count of such points saturates instead of growing.

LloydPolish

Stage S2's finishing pass: a few iterations of density-weighted Lloyd relaxation.

MacQueen leaves a slightly clumpier point set than full Lloyd, because stochastic approximation carries residual noise. Two or three explicit centroid steps clean that up.

The centroids are computed by rasterizing: assign every pixel to its nearest point, accumulate Σρ·x / Σρ per cell. That is O(pixels), exact at pixel resolution, and is what Secord's original stippling did. It also means this module needs no Voronoi library — which is why the project's dependency list stays as short as it does.

polish!(ps, d; iterations=3) -> ps

Move each point to its density-weighted Voronoi centroid, iterations times, then leave ps.weights holding the per-point Voronoi masses for the FINAL point positions – i.e. ps.weights is exactly what a fresh call to voronoi_accumulate(ps, d) would return after polish! has returned.

Those weights are the signal for variable stroke width at render time: a point owning more ink sits in a darker region and can carry a heavier line. Because they describe the final positions, they stay valid for that purpose – weights computed against an earlier iterate would misattribute mass to the wrong location.

iterations=0 computes weights for the points' current positions without moving them at all.

voronoi_accumulate(ps, d) -> (centroids, weights)

For each point, the density-weighted centroid of its Voronoi cell and the total density mass in that cell.

A pixel (i, j) contributes at its centre, x = (j−1)+0.5, y = (i−1)+0.5. Points owning no mass get their current position back as the centroid and a weight of zero — never NaN, which would propagate silently through every later stage.

Because every pixel is assigned to exactly one point, sum(weights) equals the field's total mass. That is the mass-conservation invariant.

Segment

Generates the region labels that Partition's :labels mode consumes, so subtours can follow subject boundaries rather than tone or geography.

Two modes, both returning a label matrix with labels exactly 1..k:

• segment_seeded — you supply k seed points, each grows into one region. The direct analogue of hand-painting a mask: a person decides what the subjects are, and the algorithm finds their boundaries.
• segment_auto — felzenszwalb produces however many regions its parameters yield, then they are merged down to at most k.

These are not two conveniences for the same result. Measured on the cameraman: seeded gives 1 connected component per group; automatic merge-down gives 28–46. Seeded is what reproduces Bosch's Figure 6.23's clean subject division; automatic is closer in character to tone banding, and should not be described as a substitute for a mask.

Segmentation is a parallel input to the pipeline, not a stage inside it: it consumes the source image and produces labels, while the density path consumes the same image independently. The two meet at Partition.

relabel_by_size(labels) -> Matrix{Int}

Renumber labels to exactly 1..k in descending pixel-count order: the largest region becomes group 1.

Not optional. Segmentation output is not contiguous — four of six merge trials produced label sets with gaps — and Partition would index groups that do not exist. Sorting by size (rather than by first appearance) also makes the numbering deterministic across runs, which is what makes "assign ink colour 2 to region 2" stable.

Ties are broken by the original label value, so the result never depends on scan or hash order.

segment_auto(img, k; scale=100, min_size=100) -> Matrix{Int}

Segment automatically with felzenszwalb, then merge down to at most k regions and renumber 1..k by descending size.

scale and min_size control how many regions felzenszwalb produces before merging — measured on the cameraman: (300, 500) gives 7, (100, 100) gives 18, (50, 50) gives 32. They are artistic parameters, like γ; their right values depend on the image.

The merged regions are not spatially contiguous. prune_segments merges only into adjacent regions, so every step is local, but felzenszwalb output is already fragmented and merging chains those fragments together — adjacency does not preserve global connectedness. Measured at k=2 on the cameraman: 42 and 28 connected components. Seeded segmentation gives 1 per group, which is why it is the primary path for subject-shaped subtours.

segment_seeded(img, seeds) -> Matrix{Int}

Grow one region from each seed. seeds is a vector of ((row, col), label) pairs; the number of distinct labels determines k.

Coordinates are the caller's business — chosen by eye, carried from a previous run, or hardcoded. There is deliberately no interactive picking, so the API stays scriptable and reproducible.

Returns labels renumbered 1..k by descending region size.

Partition

Stage S3: split the stipple points into k groups, each of which becomes its own closed subtour.

Bosch's Hands, Version Three (Figure 6.23) is two subtours in two colours, divided by subject — Adam's hand against God's hand and the background. No algorithm infers that division, which is why :mask exists: a person paints a second image, and its regions become the groups.

Four modes:

• :mask — a supplied image assigns each point by its coordinate. The only mode that can honour an arbitrary, subject-shaped boundary.
• :labels — a supplied integer label matrix assigns each point by its coordinate, the same way :mask does but without deriving the labels from colours. Independent of any segmentation package.
• :tone — split by density into equal-count tonal bands. Screenprint-like layering; cannot separate two objects of similar tone.
• :spatial — capacity-constrained k-means over positions. Fully automatic; boundaries fall where the clustering puts them, and adjacent tours show seams. Group sizes are capped from above at ceil(n/k) but not bounded below, so this is weaker than true balanced k-means — see _partition_spatial's docstring.

Every index appears in exactly one group — the invariant tour_indices_valid later checks on the tours themselves.

Assign each point by looking up its coordinate in an integer label matrix.

Takes the labels directly rather than deriving them from colours, because ImageSegmentation.labels_map already returns a Matrix{Int} — routing that through an RGB image would be a lossy round trip. Deliberately independent of any segmentation package: any source of an integer matrix works, including a hand-authored one.

Labels must be exactly 1..k. A gap would mean Partition returns a group index nothing maps to, which surfaces far from the cause.

Capacity-constrained k-means over point positions.

Standard Lloyd iteration, but assignment is capacity-constrained: each group may hold AT MOST ceil(n/k) points. Points are assigned in order of increasing regret — the gap between the nearest and second-nearest centroid — so the points with the least to lose from being displaced are the ones that get displaced.

The actual guarantee this gives is one-sided: cap = cld(n, k) bounds group sizes from above only. Nothing bounds a group's size from below, so a group can end up materially smaller than n/k while every group still respects the cap. This is a stable fixed point of the Lloyd iteration, not a transient — running more iterations does not fix it. Concretely, for n=101, k=11 (target size 9.18), measured group sizes come out with a spread of 3-5 between the smallest and largest group, e.g. [7,7,7,10,10,10,10,10,10,10,10]; this is typical, not a rare outlier, and gets worse when n/k is not close to an integer.

This is weaker than true balanced k-means, which would bound sizes both above and below. Do not rely on group sizes being close to n/k; only rely on no group exceeding ceil(n/k).

Equal-count tonal bands, darkest first. Band 1 holds the darkest points, band k the lightest.

Band sizes are within one of each other (floor(n/k) or ceil(n/k)), and band 1 always gets one of the larger sizes when there is a remainder, so band 1 is genuinely the darkest and largest band. But when n is not a multiple of k, the extra points do NOT land on a contiguous prefix of bands 1..rem — the assignment band = min(k, floor((rank-1)*k/n) + 1) can and does skip bands, so an earlier (darker) band can end up smaller than a later one. For example n=10, k=6 gives sizes [2,2,1,2,2,1] in band order: band 3 is smaller than band 4. Do not assume "darker bands are bigger or equal" beyond band 1.

Equal-count rather than equal-ink-mass: every layer gets (approximately) the same number of points, so each subtour takes comparable time to draw. Equal-mass banding is a reasonable alternative and is not implemented.

mask_groups(mask) -> (lookup, k)

Map a mask image to a matrix of group numbers and the number of distinct groups.

Distinct colours are sorted by (red, green, blue) and numbered 1..k in that order. Sorting rather than first-appearance order means the same mask always produces the same numbering, so "region 2 gets the second ink colour" is stable across runs and across scan orders.

The mask must have flat regions of exactly equal colour: save it as PNG rather than JPEG and disable antialiasing when you paint it, because colours are compared by exact equality. An antialiased or JPEG-compressed mask produces a distinct colour per blended edge pixel and fails loudly — "mask has 790 distinct colours but k=2 was requested" — rather than silently mis-grouping.

partition(ps, k; mode, kwargs...) -> Vector{Vector{Int32}}

Split ps into k groups of indices. Returns a vector of k index vectors; every index in 1:length(ps.coords) appears exactly once across them. Groups may be empty.

Modes and their keyword arguments:

• mode=:mask, mask::AbstractMatrix{<:Colorant} — assign by mask region.
• mode=:labels, labels::AbstractMatrix{<:Integer} — assign by an integer label matrix directly; labels must be exactly 1..k.
• mode=:tone, density::DensityField — equal-count tonal bands, darkest first.
• mode=:spatial, rng, iterations=25 — capacity-constrained k-means over positions; see _partition_spatial for the (weaker than true balance) guarantee this actually gives.

HilbertTour

A provisional stage S4: order the points along a Hilbert space-filling curve and call that the tour.

This is about 25% above optimal, which sounds bad and is in fact perfectly usable — the curve is locality-preserving, so the ordering is non-crossing and visually plausible. It produces a real picture immediately, which is worth more during development than a better tour that does not exist yet.

The near-optimal solver (greedy construction plus 2-opt/Or-opt on a Delaunay ∪ kNN candidate set) replaces this later. This module stays as the fast path.

hilbert_index(order, x, y) -> Int

Distance along the order-order Hilbert curve to lattice cell (x, y), where x, y ∈ [0, 2^order).

The standard iterative construction: walk from the most significant bit down, accumulating the quadrant contribution and rotating the coordinate frame so that the next level is expressed in the current quadrant's orientation.

hilbert_tour(ps, groups; order=16) -> TourSet

One closed subtour per group, each ordered along the Hilbert curve.

The curve is computed over the shared image bbox, not each group's own extent. Per-group extents would give each subtour a differently-scaled curve, so adjacent subtours would wander at different granularities and read as inconsistent. Indices stay global — they index ps.coords — so tour_indices_valid holds across the whole set.

Throws ArgumentError if the groups do not partition 1:length(ps.coords).

hilbert_tour(ps; order=16) -> TourSet

One closed tour visiting every point, ordered along the Hilbert curve.

order = 16 gives a 65536 × 65536 lattice, which is far finer than any point set we produce, so distinct points essentially never share a cell. When they do, the tie is broken by point index — arbitrary but deterministic, which is what reproducibility requires.

TourCandidates

The candidate edge set the TSP improvement loop is restricted to: for each point, its k nearest neighbours.

Every good 2-opt implementation limits moves to a candidate list — scanning all n-1 partners per city is O(n^2) per pass and pointless, since an improving move almost always connects near neighbours.

The parent design specifies Delaunay edges ∪ kNN. This builds kNN only, and the reason is specific to this pipeline: our points come out of MacQueen relaxation and Lloyd polish, whose whole purpose is blue-noise equidistribution. That is exactly the regime where a uniform-grid kNN is both fast and nearly complete. candidates is written as an interface so a Delaunay-backed method can be added later without touching the solver.

CandidateSet(lists)

lists[i] holds the indices of point i's candidate neighbours, sorted by increasing distance and never containing i itself.

Int32 to match TourSet.tours: these lists are scanned constantly in the improvement loop, which is memory-bound.

_knn(g, coords, i, k, heap) -> Vector{Int32}

The k points nearest coords[i], excluding i.

Expanding square rings around the home cell, with distance-based termination: the scan stops when the closest possible point in the next ring is farther than the current k-th best. BucketGrid.nearest documents why stopping on ring index instead is subtly wrong — it yields wrong neighbours near cell boundaries, an error that never crashes and shows up only as structured artifacts in the finished drawing. The same argument applies here with the bound taken against the k-th best rather than the best.

_offer!(heap, d2, idx, k)

Insert (d2, idx) into the sorted-ascending heap, keeping at most k entries. Linear insertion rather than a real binary heap: k is 8–10, so the constant factor of an array insert beats the pointer chasing, and keeping the vector sorted means heap[end] is the termination bound for free.

candidates(ps; k=8) -> CandidateSet

The k nearest neighbours of every point, via the bucket grid.

k = 8 is the parent design's recommendation (§5.2.2). Larger k finds more improving moves per pass at linear cost in scan time; the gain flattens quickly above 10 for equidistributed points.

neighbors(cs, i) -> Vector{Int32}

Point i's candidate neighbours, nearest first.

TourArray

The reference tour representation: a position array plus its inverse.

This is deliberately the obvious implementation. It exists to be believed without argument, so that TourTwoLevel — which is not obvious, and whose failure mode is a valid permutation with a wrong length — can be differential-tested against something trustworthy.

It is not a fallback and not dead code: the solver runs on it in the test suite, and optimize_tour(...; rep=:array) selects it.

Reversal is O(n) worst case, mitigated by always reversing the shorter of the two arcs. The parent design (§5.2.4) rates that acceptable to ~100k points.

ArrayTour(cities)

tour[p] is the city at position p; pos[c] is the position of city c. The two are maintained as exact inverses — check_invariants asserts it.

Cities must be exactly 1:n in some order.

_reverse!(t, from, to, len)

Reverse the len positions from from to to inclusive, walking forward with wrapping, repairing pos as it goes.

between(t, a, b, c) -> Bool

Walking forward from a, is b reached strictly before c?

false when b == a or b == c. Three integer comparisons on positions, after rotating so a sits at the origin — the rotation is what makes the wrapping case fall out without a special branch.

check_invariants(t) -> Bool

Whether pos is the exact inverse of tour. Called after every move in the test suite; this is the cheap check that catches a botched reversal immediately rather than as a wrong tour length much later.

move!(t, a, b, c, d)

The 2-opt move. Requires next(t, a) == b and next(t, c) == d. Removes edges (a,b) and (c,d), adds (a,c) and (b,d), by reversing the path b…c.

Reverses whichever of b…c or d…a is shorter — they are equivalent up to tour orientation, and picking the shorter one is what keeps this affordable.

next(t, a) -> Int32

The city following a, wrapping at the end.

prev(t, a) -> Int32

The city preceding a, wrapping at the start.

splice!(t, s, e, dest)

Or-opt: remove the forward path s…e and reinsert it immediately after dest, preserving the segment's internal orientation.

dest must not lie within s…e. Rebuilds the affected span rather than doing pointer surgery — the segment is 1–3 cities and the array is contiguous, so a rebuild is both simpler and faster than being clever.

tour_order(t) -> Vector{Int32}

The cities in visiting order. A copy, so callers cannot corrupt the state.

TourTwoLevel

The production tour representation: a two-level doubly-linked list.

The problem it solves: 2-opt reverses a tour segment, which is O(n) in a plain array. At 200k points with millions of accepted moves, that dominates everything. Splitting the tour into ~√n segments, each carrying a reversed bit, turns a reversal into flipping the order of a range of segments and toggling their bits — O(√n), touching no interior cities at all.

The cost is that every operation must now consult the reversed bit to know which way "forward" runs inside a segment. That is exactly the kind of detail that yields a valid permutation with a wrong length: the tour still renders, the picture still looks plausible, and cheap invariants still pass. This is why TourArray exists and why the differential test in test_ariadne_tourrep.jl is not optional.

Implementation note: this keeps city order inside each segment as a Vector rather than true linked-list pointers. At √n ≈ 450 for n = 200k, the vector is cache-resident and splitting it is a copy of a few hundred Int32s, which beats pointer chasing. The "linked list" is at the segment level.

Segment

A run of consecutive cities. seq orders segments around the tour; rev means the cities in items are visited back-to-front.

TwoLevelTour(cities)

Cities partitioned into ~√n segments.

segof[c] is the index into segs of the segment holding city c; idxof[c] is its position within that segment's items, always in storage order, not visiting order. Reading a visiting order therefore means consulting rev.

_position(t, c) -> Int

The city's 1-based position in visiting order. O(number of segments), so this is for between and correctness checks, not for the hot path.

_rebuild!(t, order)

Repartition order into ~√n segments and renumber everything. O(n).

Called at construction and whenever the segments have fragmented enough that their sizes no longer approximate √n.

Segment with the next seq, wrapping to the lowest.

Segment with the previous seq, wrapping to the highest.

between(t, a, b, c) -> Bool

Walking forward from a, is b reached strictly before c?

Same contract as TourArray.between, including b == a and b == c being false. Computed from visiting positions rotated so a is the origin.

check_invariants(t) -> Bool

Whether the structure is self-consistent: segof/idxof point where the segments say they do, segment sequence numbers are distinct, and the visiting order is a permutation of 1:n.

move!(t, a, b, c, d)

The 2-opt move: remove (a,b) and (c,d), add (a,c) and (b,d), reversing the path b…c. Requires next(t, a) == b and next(t, c) == d.

Reverses the shorter of b…c and d…a, matching TourArray — the two are equivalent up to orientation, and the differential test only passes if both representations produce the same cycle, which they do either way.

The reversal itself materializes the affected span, reverses it, and rebuilds. That is O(n) rather than the O(√n) the structure permits, and it is a deliberate first cut: correctness first, verified by the differential test, with the segment-range optimization to follow once the benchmark shows the reversal cost actually dominates. The structure — segments, seq, rev — is already in place to support it.

next(t, a) -> Int32

The city following a. Steps within the segment when possible — respecting the reversed bit — and otherwise moves to the head of the next segment by seq.

prev(t, a) -> Int32

The city preceding a. The mirror of next.

splice!(t, s, e, dest)

Or-opt: remove the forward path s…e (1–3 cities) and reinsert it immediately after dest, preserving its internal orientation.

Same contract as TourArray.splice!, including the error when dest lies inside the segment.

tour_order(t) -> Vector{Int32}

The cities in visiting order: segments by seq, and each segment's items forward or backward according to its rev bit.

TourRep

The tour-representation interface: five operations plus a distance helper, dispatched over ArrayTour and TwoLevelTour.

This exists because two modules need the same operations. TourSolve drives them for candidate-limited 2-opt/Or-opt; TourCrossings drives them for crossing repair, which deliberately ignores the candidate list. Naming the shared contract here keeps TourCrossings from reaching into TourSolve's privates — a coupling that would break quietly the first time either is refactored.

The dispatch is written out explicitly rather than via a common supertype because both representation modules export identically-named functions (next, move!, …) that collide with each other and with Base, so every call must be fully qualified.

Walking forward from a, is b reached strictly before c?

rep_dist(coords, i, j) -> Float64

Euclidean distance between two points, by index.

rep_move!(r, a, b, c, d)

The 2-opt move. Requires rep_next(r, a) == b and rep_next(r, c) == d. Removes edges (a,b) and (c,d), adds (a,c) and (b,d).

The city following a in the tour.

The cities in visiting order.

The city preceding a in the tour.

rep_splice!(r, s, e, dest)

Or-opt relocation: move the forward path s…e to just after dest, preserving its internal orientation.

TourSolve

Stage S4: a PointSet becomes a near-optimal TourSet.

Greedy-edge construction on the candidate set, then 2-opt and Or-opt driven by don't-look bits. The parent design (§5.2.3–4) puts greedy at ~15–20% above optimal and the improvement pass at ~4–6%, against the Hilbert curve's ~25%.

The solver is written against a tour representation interface — next, prev, between, move!, splice! — and never learns which representation it is driving. TourArray is simple and trusted; TourTwoLevel is fast. That separation is what lets the differential test verify the fast one against the simple one.

LockedEdges(n)

Edges the improvement pass must not break, as an unordered-pair set.

Empty by default and unused by P1's own callers: this exists because P4 (symmetric tours) must pin the seam edges where a fundamental domain's tour meets its mirrored copy, and retrofitting a lock check into a converged optimizer is far harder than having the move filter consult one from the start.

The unordered storage matters. A locked edge's endpoints swap adjacency direction whenever a 2-opt reverses the segment containing it, so a direction-sensitive check would silently stop protecting the edge after the first reversal.

_default_max_sweeps(m) -> Int

Default sweep cap for remove_crossings!, as a function of the SUBTOUR size m (not the whole point set — grouped mode solves each group independently, so a small group must not inherit a cap sized for the total n).

max(100, 2m). Justification: measured sweeps-to-convergence on shuffled (worst-case, adversarial) input grows roughly linearly in m —

2m comfortably covers the observed slope (roughly m/13 at worst) with better than 10x headroom at every measured size, and max(100, ...) keeps small subtours from getting an unnecessarily tight cap when 2m < 100. Actual solver output (not shuffled) converges in single digits of sweeps regardless of m — 5 sweeps at m=1000, 8 at m=3000 — so this cap is essentially never approached on the realistic path; it only matters for adversarial inputs, and each sweep on an already-clean tour costs one cheap detection pass (the loop returns immediately once a sweep finds no crossings), so a generous cap costs nothing when it is not needed.

_find(parent, x) -> Int32

Union-find with path halving. Used to reject an edge that would close a cycle before every city is on the path.

_solve_subset(ps, idxs, k, rep, locked, max_passes, jordan=true, max_sweeps=nothing) -> Vector{Int32}

Solve the TSP over the subset idxs and return global indices.

Builds a compact sub-PointSet so the candidate grid is sized to the subset's own density, then maps the resulting local tour back through idxs. When jordan is true, runs the crossing repair (via JORDAN_REPAIR) on this subset before mapping back — which is what makes grouped mode repair each subtour independently, for free.

max_sweeps, when nothing, defaults via _default_max_sweeps on the SUBSET's own size m = length(idxs) — deliberately not the whole point set's n, since grouped mode solves each group as an independent subtour.

Mark a city dirty and queue it for rescanning.

_try_or_opt!(rep, coords, cs, a, locked, dirty, queue) -> Bool

Look for an improving Or-opt move: relocate the segment of 1–3 cities starting at a to sit beside one of a's candidate neighbours.

Or-opt requires no segment reversal, so it is cheap in both representations, and it reaches improvements 2-opt structurally cannot — notably moving a single badly-placed city between two distant parts of the tour.

_try_two_opt!(rep, coords, cs, a, locked, dirty, queue) -> Bool

Look for an improving 2-opt move on an edge incident to a.

For each tour direction, take the edge (a, b) and try replacing it with (a, c) for each candidate c. The paired removal is (c, d) where d follows c in the same direction. The move is improving when d(a,b) + d(c,d) > d(a,c) + d(b,d).

Candidates are scanned nearest-first and the scan stops as soon as d(a,c) >= d(a,b): beyond that no candidate can produce a positive gain on this edge, since the remaining terms cannot compensate. This is the standard neighbour-list pruning and it is most of the speed.

_walk(adj, n) -> Vector{Int32}

Read the adjacency structure off as a visiting order, starting at city 1.

greedy_tour(ps, cs) -> Vector{Int32}

Greedy edge matching: consider candidate edges shortest-first, accepting one when both endpoints still have degree < 2 and it does not close a premature cycle.

The candidate edges alone will not complete a tour — greedy always strands some fragments, and short candidate lists guarantee it. The leftover fragment endpoints are then joined nearest-first, which is O(f²) in the fragment count f. That is acceptable because f is small (a few percent of n) and it avoids a second spatial index.

improve!(ps, rep, cs; locked=nothing, max_passes=50) -> Int

2-opt and Or-opt on the candidate set, driven by don't-look bits. Mutates rep. Returns the number of accepted moves.

Don't-look bits. After a move, only the cities whose incident edges changed can newly admit an improving move, so a queue of "dirty" cities replaces rescanning everything. A city is popped, scanned, and — if nothing improves — marked clean and left alone until a neighbouring move dirties it again. This is what turns 2-opt from O(n²) per pass into something near-linear in practice.

Why moves are restricted to candidate neighbours. An improving 2-opt move almost always shortens at least one edge to something in the candidate list; scanning all partners costs O(n) per city for a vanishing extra yield.

The consequence, spelled out because it looks like a bug otherwise: this does not eliminate all crossings. A crossing whose partner is in neither endpoint's candidate list is never considered. Guaranteed non-crossing is P2's job, and P1 deliberately does not assert it.

is_locked(le, a, b) -> Bool

Whether the undirected edge (a, b) is locked. Order-independent by construction.

lock_edge!(le, a, b)

Mark the undirected edge (a, b) as unbreakable.

optimize_tour(ps, groups; kwargs...) -> TourSet

One near-optimal closed subtour per group.

Each group is solved as an independent TSP — smaller subproblems land closer to optimal, per the parent design §6.2(a). No moves cross a group boundary, so the seams between adjacent subtours are inherited from the partitioning stage, not introduced here.

Indices stay global: they index ps.coords, so tour_indices_valid holds across the whole set. Throws ArgumentError if the groups do not partition 1:length(ps.coords), matching hilbert_tour.

max_sweeps (see the single-tour method's docstring) is forwarded to every group's repair unchanged; when nothing, each group gets a default sized to its own subtour length, not the whole point set — a small group should not pay for a cap computed from n.

optimize_tour(ps; k=8, rep=:twolevel, locked=nothing, max_passes=50, jordan=true,
             max_sweeps=nothing) -> TourSet

One near-optimal closed tour through every point.

Greedy-edge construction on a k-nearest-neighbour candidate set, then 2-opt and Or-opt to convergence. This is stage S4's default; hilbert_tour remains as the instant, ~25%-above-optimal preview path.

rep selects the tour representation — :twolevel (production) or :array (the reference implementation). Both produce the same tour; the option exists so the test suite can compare them.

jordan (default true) runs the crossing-repair pass from TourCrossings after improve! converges, certifying the returned tour has zero self-intersections. Pass jordan=false to skip it and reproduce P1 behaviour exactly. See JORDAN_REPAIR for how the repair is wired in without a module cycle.

max_sweeps, when jordan is true, caps remove_crossings!'s re-sweep loop. nothing (the default) lets _solve_subset pick a size-appropriate default — see its docstring for the exact formula and the measurements behind it. Pass an explicit integer to override, e.g. to budget more sweeps for an unusually adversarial input, or a small one to force hit_cap in a test.

tour_length(ps, order) -> Float64

Total Euclidean length of the closed cycle visiting order.

TourCrossings

Stage S4b: make the tour a simple closed curve — no self-intersections.

This is the project's success criterion #2, and P1's solver cannot deliver it. improve! restricts 2-opt moves to each city's k-nearest candidate list, and a crossing between two geometrically distant tour segments need not have either partner in the other's candidate list, so no improving move is ever considered. Repair here deliberately bypasses the candidate list.

Every crossing admits an improving 2-opt move — if (a,b) and (c,d) cross then |ac| + |bd| < |ab| + |cd| — so removing crossings also shortens the tour. The Jordan guarantee and tour quality point the same way.

Detection buckets each tour edge into every grid cell its bounding box overlaps, then tests only edge pairs sharing a cell. Measured crossing counts on real solver output are sparse (56 in 10,000 edges), which is why this suffices and the Bentley–Ottmann sweep line is deferred.

The correctness hazard is the bucketing. An edge omitted from any cell its bounding box overlaps means a crossing is missed silently — no crash, nothing visibly wrong, just a tour that is not a Jordan curve while claiming to be. The brute-force oracle in the test suite is what makes this trustworthy; it is not optional.

CrossingReport(repaired, remaining, blocked, swept, hit_cap)

The outcome of remove_crossings!.

A struct rather than a bare count because three distinct outcomes must be distinguishable: a clean fixpoint, a pass stopped by the sweep cap, and crossings left in place because a locked edge forbade the repair. Collapsing them into one integer would make an unmet guarantee look like a met one.

blocked is a running total accumulated across every sweep, not just the last one — it answers "how many repair attempts were skipped for a locked edge over the whole run?" regardless of whether the loop stopped at a blocked-only fixpoint or at max_sweeps. The two are independent: a run can be capped (hit_cap == true) while also having blocked crossings on every sweep (blocked > 0), and the field must say so rather than silently reporting only the final sweep's local count.

The certified-zero guarantee is exactly remaining == 0. It implies !hit_cap, by construction: remaining == 0 is only ever returned from the branch that finds a sweep with no crossings at all, and that branch always reports hit_cap == false.

It does not formally imply blocked == 0. Because blocked accumulates across sweeps, a run could in principle skip a locked crossing on one sweep and then have an unrelated repair elsewhere dissolve that same crossing later, ending clean with a nonzero count. A search over 4798 runs (seeds 1–400, n ∈ {12,20,30,40}, each of the first three crossings locked in turn) found no such case, so it is at most rare — but read remaining for the guarantee, not blocked.

_bucket_edges(ps, tour, cellsize) -> Dict{Tuple{Int,Int},Vector{Int}}

Map each grid cell to the indices of every edge whose bounding box overlaps it.

Conservative by construction: an edge is registered in every cell spanned by its bounding box, from floor of the minimum corner to floor of the maximum. An edge is a straight line so it cannot leave its own bounding box, and a crossing requires overlapping bounding boxes — therefore two crossing edges always share at least one cell. Narrowing this (registering only endpoint cells, say) would silently miss crossings on long edges.

_cellsize(ps, tour) -> Float64

Cell edge length, sized at roughly a typical tour edge length so a typical edge touches only a handful of cells. Falls back to the bbox diagonal for degenerate inputs, which puts everything in one cell — slow but correct.

Uses the median, not the mean, of edge lengths, and clamps the result from below relative to the bbox diagonal. Both guard the same failure mode: a tour with a skewed edge-length distribution (many tiny edges, a few enormous ones — e.g. a dense cluster plus one distant outlier) drags the mean toward zero even though the long edges' bounding boxes still span the full extent. Since _bucket_edges registers every edge in every cell its bbox overlaps, a near-zero cellsize turns each long edge's registration into O((length/cellsize)^2) cells — quadratic memory blowup on a tour with zero crossings. The median is robust to that skew because it tracks the typical edge rather than being dragged by outliers, and the diag/2048 floor is a second, independent backstop: even if the median itself is small (e.g. every edge is short except one outlier), cellsize is never allowed to fall far enough below the bbox scale to blow up a single long edge's cell count. This cannot affect correctness — _bucket_edges iterates the full bbox-derived cell range whatever cellsize is, so any positive finite cellsize still buckets every edge conservatively; only performance is at stake here.

_edge_endpoints(tour, i) -> (city_from, city_to)

Edge i of the closed tour: from tour[i] to tour[i+1], wrapping.

_poly_cellsize(pts) -> Float64

Cell edge length for a closed polyline, by the same rule _cellsize uses for a tour: the median edge length, floored at diag/2048.

The two share a rule rather than a code path because a polyline has no PointSet and no bbox — the extent is derived from the points themselves. See _cellsize for why the median and the floor are both load-bearing: they guard against a skewed edge-length distribution turning one long edge's registration into quadratic memory.

_poly_edge(pts, i) -> (P2, P2)

Edge i of the closed polyline: from pts[i] to pts[i+1], wrapping.

_side(a, b, p) -> Int

Which side of the directed line a→b does p lie on? +1 left, -1 right, 0 exactly collinear. The sign of the 2-D cross product.

count_crossings(ps, tour) -> Int

How many pairs of tour edges properly cross. Zero means the drawing is a simple closed curve.

count_crossings(pts::Vector{P2}) -> Int

How many pairs of edges of the closed polyline properly cross.

find_crossings(ps, tour) -> Vector{Tuple{Int,Int}}

Every properly-crossing pair of tour edges, as edge index pairs (positions in tour, not city ids), with i < j.

Pairs sharing an endpoint are skipped: adjacent tour edges always share one, and a shared endpoint can never be a proper crossing anyway.

find_crossings(pts::Vector{P2}) -> Vector{Tuple{Int,Int}}

Every properly-crossing pair of edges of the closed polyline through pts, as edge-index pairs with i < j.

This is the geometric core. find_crossings(ps, tour) is the same computation with the points reached through a tour's indices; a sampled Bézier curve is the same computation with the points reached directly. Naming one implementation means the brute-force oracle that protects the tour detector protects the curve detector too.

Index-adjacent edges are skipped: edges i and i+1 always share a point, and edge n shares one with edge 1. Coordinate equality between non-adjacent edges is NOT skipped — segments_properly_cross already rejects it, because a shared endpoint zeroes an orientation test.

remove_crossings!(ps, rep; locked=nothing, max_sweeps=100) -> CrossingReport

Repair tour self-intersections until a full sweep finds none, mutating rep.

Each repair is an ordinary 2-opt move — the crossing pair (a,b) × (c,d) becomes (a,c) and (b,d) — applied through the representation interface, so this works on either representation. The difference from improve! is that the pair comes from geometric detection rather than a candidate list, which is what lets it fix crossings improve! structurally cannot see.

Why it re-sweeps. Repairing one crossing can create a new one elsewhere, so this is a fixpoint loop, not a single pass. Each accepted repair strictly shortens the tour, so the same crossing cannot recur forever; but there is no clean a-priori bound on the number of sweeps, hence max_sweeps and the hit_cap flag.

How many sweeps to expect depends heavily on the regime. On real solver output (greedy construction + improve!, the actual P1→P2 pipeline this module exists to finish) convergence is fast — a handful of sweeps regardless of size (5 sweeps measured at n=1000, 8 at n=3000), because improve! already removes most crossings and leaves only the few improve!'s candidate-list restriction cannot see. On adversarial input — a shuffled tour, far denser in crossings than anything the solver produces — sweeps-to-convergence instead grow roughly linearly in tour size: 13 at n=50, 21 at n=100, 30 at n=200, 51 at n=400, 81 at n=800, 130 at n=1600. TourSolve.optimize_tour sizes its default max_sweeps off this adversarial slope (see TourSolve._default_max_sweeps) precisely because a flat default would run out of headroom on large, dense inputs even though it never does on realistic solver output.

Locked edges. A repair that would break a locked edge is skipped and counted in blocked rather than applied. Those crossings remain, so remaining > 0 and the Jordan guarantee legitimately does not hold — the caller is told rather than misled. This matters for P4, where seam edges are pinned.

blocked is a running total, accumulated across every sweep the loop performs — not just the last one. A caller asking "how many crossings could not be repaired because of locks?" wants the whole run's answer, whether the loop stopped because it reached a fixpoint of blocked-only crossings or because it hit max_sweeps while locks and open crossings were both still present on every sweep. Both return paths report the same accumulated quantity, so blocked means one consistent thing regardless of why the loop stopped.

segments_properly_cross(p1, p2, p3, p4) -> Bool

Whether segments p1p2 and p3p4 cross transversally.

Deliberately strict: a 0 from any orientation test means some endpoint is collinear with the other segment, which covers touching at a point, T-junctions, and collinear overlap. None of those is a crossing for our purpose — the tour is a closed polyline whose adjacent edges legitimately meet at shared points, and treating a touch as a crossing would send the repair loop chasing non-problems forever.

TourSmooth

Stage S5a: the tour becomes a smooth curve instead of a spiky polyline.

The raw tour looks angular at high zoom and is hard on a plotter's acceleration. Catmull-Rom is the right smoothing family here because it is interpolating — the curve passes through every stipple point, so it does not move ink and therefore does not change tone. A B-spline is smoother but approximating: it pulls the curve off the points and lightens the image.

Smoothing threatens the Jordan property P2 established. Catmull-Rom overshoots outside the corner on tight turns, and two overshoots can intersect even when the underlying polyline is simple — so a tour certified crossing-free can render as a curve that crosses itself. TourCrossings.count_crossings cannot see this, because it inspects the tour rather than the curve. Repair therefore happens here, by reducing tension where overshoot causes a crossing.

That repair has a proven floor: at tension 0 a segment IS the straight polyline edge, which P2 already guarantees is crossing-free.

BezierSegment(p0, c1, c2, p3)

One cubic Bézier: endpoints p0, p3 and control points c1, c2.

Cubic rather than a sampled polyline because SVG emits these directly as C commands — one compact instruction per segment instead of thousands of tiny line segments.

SmoothPath(segments, tension, closed)

A smoothed tour: one BezierSegment per input point (closed), plus the tension actually used for each.

tension records what was used, not what was requested, because the repair loop lowers it per segment. That turns "which corners needed relaxing, and how much" into data a test can assert on rather than invisible internal state — the same reason CrossingReport is a struct and not a bare count.

SmoothReport(relaxed, remaining, min_tension, iterations, hit_cap, degenerate)

The outcome of smooth_polyline.

A struct rather than a bare count because the outcomes are distinct and a caller needs to tell them apart:

degenerate is the artistically interesting one: the Jordan guarantee held, but the picture has a flat spot where a corner used to curve. That is a deliberate trade the caller may want to know about, so it is reported rather than silently accepted.

A third stopping condition, represented without a new field. The repair loop can also stop because a pass lowered no tension: every segment involved in a remaining crossing is already pinned at the 0.0 floor, so another pass would resample and detect the same crossings for free with no chance of changing anything. That is neither "clean" nor "gave up because the budget ran out" — it is "gave up because there is nothing left to try." It is still fully readable from the existing fields: remaining > 0 && !hit_cap && iterations < max_relax means the loop stopped early with crossings left, but not because max_relax was reached — i.e. every implicated segment had already saturated at tension 0. No new field is added for this because the combination is already unambiguous and every existing caller that checks remaining == 0 for the guarantee, or hit_cap for "should I raise the cap," continues to get the right answer.

_cr_segment(p0, p1, p2, p3, a) -> BezierSegment

The uniform Catmull-Rom segment from p1 to p2, as a cubic Bézier.

c1 = p1 + a * (p2 - p0) / 6
c2 = p2 - a * (p3 - p1) / 6

The /6 (rather than /3) makes a == 1 the standard Catmull-Rom tangent (p2 - p0)/2, since a cubic Bézier's start tangent is 3*(c1 - p0).

At a == 0 this returns exactly the straight segment c1 == p1, c2 == p2. The repair loop depends on that identity holding exactly.

_rebuild(pts, tens) -> SmoothPath

Rebuild every segment at its own current tension.

_segment_samples(seg) -> Int

How many samples this segment needs, from how far it bows off its chord.

Flat segments need almost none; tight turns need many. Overshoot — the thing detection is looking for — happens exactly where curvature is high, so concentrating samples there puts them where they matter instead of spreading them uniformly.

_self_crossing(seg, k) -> Bool

Whether the segment's own arc, sampled at k points, self-intersects.

A cubic Bézier can loop on itself at most once, so this is cheap and bounded — unlike the pairwise cross-segment check, it never needs a spatial index.

The sample sequence here is an open arc from t=0 to t=1, not the implicitly-closed polyline find_crossings(::Vector{P2}) assumes: there is no edge from the last sample back to the first, because nothing wraps a single segment around on itself. Reusing find_crossings directly would silently introduce a spurious closing edge, so this loops directly over segments_properly_cross and skips only genuine index-adjacent pairs (which always share an endpoint) rather than the closed-curve wrap skip.

bezier_point(seg, t) -> P2

The point at parameter t ∈ [0,1] along the cubic, by the Bernstein form.

catmullrom_path(pts, tension) -> SmoothPath

Smooth the closed polyline through pts. One segment per point; indices wrap, so the segment from pts[i] to pts[i+1] uses pts[i-1] and pts[i+2].

Fewer than three points cannot form a closed curve, so the result is empty.

path_crossings(path; scale=1) -> Vector{Tuple{Int,Int}}

Self-intersections of the smoothed curve, as segment index pairs.

Sampled crossings between different segments are mapped back through sample_path's owner vector and deduplicated, so a single overshoot spanning several samples reports once per segment pair.

A single segment can also loop on itself. At full tension, a cubic Catmull-Rom-derived Bézier self-intersects in roughly 0.3% of random 4-point trials — not the rare edge case the pairwise check alone would suggest — and _segment_samples already allocates enough resolution to see it. Silently dropping same-segment sample pairs (as an earlier version of this function did) makes that loop invisible to Task 4's repair loop regardless of scale, which is a silent hole in the Jordan guarantee: a looped corner would never get its tension lowered. So a same-segment loop is detected separately, with _self_crossing, and reported as (si, si) — a diagonal pair the repair loop can recognise as "relax this one segment" rather than "reconnect two edges."

sample_path(path; scale=1) -> (Vector{P2}, Vector{Int})

Sample the curve into a closed polyline, and report which segment each sample came from.

The second return value is what lets repair act: detection finds crossings between samples, but tension is a property of segments. Without the map, the repair loop cannot tell which Bézier to relax.

scale multiplies every segment's sample count. It exists so the sufficiency of the sampling can be measured — doubling it must not reveal new crossings — rather than assumed.

The closing point is not repeated: the polyline is implicitly closed, matching find_crossings(::Vector{P2}).

smooth_polyline(pts; tension=0.5, preserve_jordan=true, max_relax=20)
    -> (SmoothPath, SmoothReport)

Smooth the closed polyline through pts, optionally repairing any self-intersection the smoothing introduces.

The repair. Where the curve crosses itself, tension is halved on the four segments that produced the crossing — the two the samples belong to and their two predecessors — and those segments are rebuilt. Repeat until clean or max_relax passes.

Why all four and not the sharpest. A crossing between stretches far apart in the tour has no single guilty corner; the bulge is produced from both sides. Relaxing only one leaves the crossing in place and the loop can bounce between candidates without converging.

Why it terminates. Tension only decreases, halving each pass, and at 0 a segment is exactly the straight polyline edge — which the tour-level Jordan guarantee already certifies as crossing-free. The floor is known safe, so max_relax is a budget, not a safety net.

preserve_jordan=false smooths without repairing: maximum smoothness for a caller who does not care about self-intersection. The default preserves the guarantee.

Valid range for tension. Must be finite and non-negative: 0.0 is the straight polyline (the repair floor), 1.0 is standard Catmull-Rom, and values above 1.0 are accepted (an unusually pronounced overshoot, not a meaningless one) but are increasingly likely to need repair and to hit max_relax. Negative tension inverts the tangent direction and produces self-looping curves even on gentle input, which is never useful here, so it throws ArgumentError rather than silently producing a degenerate picture.

smooth_tours(ps, ts; tension=0.5, preserve_jordan=true, max_relax=20)
    -> (Vector{SmoothPath}, Vector{SmoothReport})

Smooth every subtour independently.

Per subtour, exactly as optimize_tour's jordan works: each subtour is a separate closed pen stroke, so each gets its own curve and its own guarantee. Crossings between distinct subtours are out of scope by design — they are separate strokes, and a plotter lifts the pen between them.

AriadneRender

Stage S5: a TourSet becomes an SVG, PDF, or PNG.

Each tour is drawn as one lines! call over a single point vector. This is the same lesson TruchetRender records for poly!: Makie treats every call as a separate plot object, so per-segment drawing is catastrophically slower (~67× on the Truchet grid) for pixel-identical output. Measured here: a 100,000-point tour saves in about 4 s at 2.35 MB of SVG, 0.61 MB of PDF.

Points arrive in image space, y downward. The flip into Makie's upward axis happens here and only here.

Expand a scalar to one value per tour, or check that a supplied vector has exactly one entry per tour.

Colour and width are artistic controls: a scalar means "the same for every subtour", a vector means "one per subtour, in tour order".

render_tours(ps, ts; figure_size=(1000,1000), background=:white, ink=:black,
             linewidth=0.6) -> Figure

Draw every tour. The axis is fixed to the point set's bbox with DataAspect(), so the drawing keeps the source image's proportions regardless of figure size.

ink and linewidth each accept either a scalar, applied to every tour, or a vector with one entry per tour (in tour order). Colour and width are artistic controls chosen by the person making the picture, not machinery for encoding tone: Bosch's Figure 6.23, Hands, Version Three, draws Adam's hand and God's hand as two subtours in two colours at the same line width, which is the motivating case for per-tour colour here.

preserve_jordan and max_relax are smooth_tours' escape hatch, exposed here so a caller can reach it without dropping to the smoothing API directly: preserve_jordan=true (the default) repairs any curve self-intersection smoothing introduces, up to max_relax passes per subtour; false skips repair for maximum smoothness. If the guarantee is not met – or is met only by driving some segment to a flat, degenerate straight line – this function warns rather than rendering the discrepancy silently. See the module's TourSolve._solve_subset for the same pattern applied to the tour-level Jordan guarantee.

save_tours(path, ps, ts; kwargs...) -> String

Render and write. The format follows the file extension: .svg, .pdf, or .png. Returns the path.

smooth_points(ps, path) -> Vector{Point2f}

The drawable polyline for a smoothed tour: the curve sampled densely, flipped into plot space, closed by repeating the first point.

Makie draws polylines, not Béziers, so the curve is sampled for the screen. The SmoothPath remains the real object — it is what SVG C commands would be emitted from, and what crossing detection runs on.

tour_points(ps, tour, closed) -> Vector{Point2f}

The drawable polyline for one tour: coordinates flipped into plot space, with the first point repeated at the end when closed.

The flip is y_plot = (ymin + ymax) − y_image, which maps the top of the image to the top of the axis.