Idea

How much area does a retinal neuron actually cover, and how does that compare to naive approximations like convex hulls?

Traditional geometry-based approaches (convex or concave hulls) wildly overestimate coverage because neurons are sparse: long, thin branches reach far from their neighbors. The real coverage region is where a neuron can actually make contact — its skeleton thickened by its inherent radius (measured from EM reconstruction) plus a small buffer for reach.

This project computes that region efficiently in two ways:

1. exact — high-precision geometric union of tapered-capsule polygons
2. raster — fast approximate grid-based method, within ~1% of exact

Both support removing holes (gaps within the coverage) and downsampling to a fixed vertex budget for compact storage.

---

Implementation

from neuron_hulls import simulate_skeleton, skeleton_coverage

skeleton = simulate_skeleton(seed=0)                  # nodes (n, 2), edges (m, 2), radii (n,)
result = skeleton_coverage(
    *skeleton,
    reach=3.0,          # extra buffer beyond each node's own radius
    fill_holes=True,    # or max_hole_area=... to drop only the small ones
    max_vertices=256,   # downsample the stored hull
)
result.area, result.polygon, result.n_vertices, result.n_holes

The region is

union over edges of { x : dist(x, edge) <= radius_scale * radius(x along edge) + reach }

i.e. every neurite thickened by its own (linearly interpolated) radius plus reach. radius_scale=0 drops the radii entirely, which is the right call for a cell far larger than its neurites are thick.

What the arbor actually covers, against a convex hull

Simulated cells, reach = 3:

cell	nodes	convex hull	coverage	holes filled	% of convex
dense	812	47,387	14,997	22,777	32%
sparse	288	34,598	6,314	6,263	18%
long_range	293	96,686	8,471	8,525	9%
large	3,788	554,703	98,581	180,692	18%

The convex hull overestimates by 3x on a densely branched cell and 11x on three long neurites — which was the complaint that started this.

Backends

Two, computing the same region, both reporting which one ran in .method:

• exact — one tapered-capsule polygon per edge (the convex hull of the two end disks, built as a single vectorised vertex array), dissolved with a cascaded union. Accurate to the polygon approximation of the round caps: 0.04% on an elongated capsule at the default quad_segs=8, verified against closed-form areas.
• raster — one Euclidean distance transform on a grid, then each pixel's nearest segment is recovered and the tapered distance test evaluated in closed form. Within ~1% of exact, and its cost follows the grid rather than the skeleton.

auto (the default) takes raster from 2,000 edges upwards, exact below.

Speed

Same cell, edges subdivided so only the node count changes:

nodes	exact	raster	raster, area only	speed-up	area vs exact
3,788	256 ms	101 ms	72 ms	2.5x	+0.96%
15,149	1,297 ms	101 ms	71 ms	12.8x	+0.59%
60,593	4,866 ms	131 ms	95 ms	37x	+0.19%
121,185	9,434 ms	168 ms	139 ms	56x	+0.12%

raster is essentially flat in node count — the grid, not the skeleton, sets the cost. grid (default 512) is the speed/detail dial; as_polygon=False skips tracing the outline when only the number is wanted.

Three things were needed to make the raster accurate enough to be worth using, each of which was worth about a percent or more:

1. Nearest segment, not nearest seed pixel. Storing one representative point per seed pixel caps the effective sampling along a neurite at the pixel size, and the resulting scalloping cost 3.5% on a dense skeleton. Storing the owning segment and evaluating point-to-segment distance in closed form removes it.
2. Let the projection overhang the segment ends by half a pixel. Each pixel is answered by one segment, but the truly nearest point often lies just past its end on the next one along; without the overhang, coverage came out short wherever segments are shorter than a pixel.
3. Maximise radius(t) - distance(t), not -distance(t). On a tapering segment the deepest coverage is not at the closest point — a fat proximal end reaches ground the thin end it points at cannot. Ignoring this got a disk nested inside a fatter one 44% wrong.

Downsampling the stored hull

A traced outline runs to ~1,200 points. max_vertices searches for the smallest Douglas-Peucker tolerance that meets the budget, and max_area_error (default 2%) caps how much area that may cost. When the two conflict the area wins and the polygon comes back over budget, so a too-tight budget degrades into a larger polygon rather than a wrong one:

budget	allow 2%	allow 20%
none	1,185 pts	1,185 pts
512	396 pts, −1.8%	396 pts, −1.8%
256	396 pts, −1.8%	191 pts, +1.9%
128	396 pts, −1.8%	90 pts, +13.8%

Pair it with fill_holes=True: every hole costs at least four vertices, so a budget fights the holes before it touches the detail. At 20% tolerance the shape starts shattering into spikes, which is what the area guard exists to stop by default.

Layout

neuron_hulls/simulate.py   growth model for 2-D skeletons; EXAMPLES presets
neuron_hulls/coverage.py   skeleton_coverage() and the two backends
neuron_hulls/geometry.py   tapered capsules, hole removal, vertex budgeting
neuron_hulls/plot.py       plot_skeleton(), plot_coverage()
scripts/demo.py            simulate -> measure -> plot, writes figures/
scripts/benchmark.py       the tables above

One thing that did not work

Quantising the buffer distances into a few bins and making one shapely.buffer call per bin on a MultiLineString of that bin's edges. It sounds like it should beat unioning thousands of separate polygons, but GEOS dissolves the whole MultiLineString in a single overlay and that degrades sharply with segment count: 11 s at 15k edges against 0.9 s for exact, and four minutes at 60k. Chunking the calls bought only a constant factor. Many small convex polygons plus a cascaded union is the shape of problem GEOS is good at, so the backend was removed rather than shipped.