PyMicroglia

Microglial imaging analysis: registration, filtering, segmentation, tracing and figures for bioluminescence time-lapses, over a keyed artefact store that makes re-analysis cheap.

Read the documentation

The store is the idea the package is built around. Registration is the expensive step, and its whole output is a 298 KB table of per-frame shifts against a 10.8 GB input — 0.003%. So the expensive steps save what they derived — transforms, masks, labels, traces — rather than a compressed copy of the pixels, and a changed parameter misses the cache on its own. Keep the recipe, not the cooked dish.

Install

pip install PyMicroglia

Optional extras: kit (run records and house style), figure, seg, video, rhythm, scn (the automatic SCN outline, which lives in PySCNSlice), test.

git clone https://github.com/Jay2owe/PyMicroglia.git
cd PyMicroglia
pip install -e ".[kit,test]"

Self-contained, by rule

PyMicroglia depends on nothing in Protocols/Analysis. Those scripts and macros are sources it was copied from — never imported, never called, never edited. They keep working exactly as they do today, and the package works with them absent. tests/test_self_contained.py enforces it, and the whole suite passes with the folder unavailable:

$env:PYMICROGLIA_PROTOCOLS = "/nowhere"; python -m pytest

analysis_kit is imported softly, through pymicroglia._optional.kit() only: if the audit layer is missing, the run record is skipped and the science carries on. circadian_workbench, by contrast, is a hard optional dependency — silently skipping a periodogram would be worse than failing.

What works today

Stages 01 to 10, bar three gates that need stage 11's pipeline: the package, the parameter reader, the action registry, the command line, the artefact store, the image series layer, registration, cosmic-ray removal, unmixing, the two display filters, segmentation, regions of interest, traces, controls, the rhythm handoff, the figures and the videos.

pymicroglia doctor                 # healthy? and which Python is answering?
pymicroglia discover               # 26 actions, what is implemented, what is pending
pymicroglia describe remove_cosmic_rays
pymicroglia validate remove_cosmic_rays seed_z=8

Twenty-six actions are declared and none are pending: every one binds to a function that exists. That count is how the plan reported its own progress while it was being built, and validate says which target is missing rather than failing obscurely. A bound action accepts exactly the parameters describe lists, checked in tests/test_registry.py, so a keyword an agent is shown is never one the call rejects six hours in.

The catalogue behind describe ships as a data file, generated once from the protocol scripts and then owned by this package:

python tools/regenerate_catalogue.py          # reseed from the sources
python tools/regenerate_catalogue.py --check  # fail if it has drifted

One parameter name means one thing, so the shared vocabulary holds each name's type, units and prose. Defaults are per action, because crf is honestly 18 in one video exporter and 20 in another; describe merges the two so an agent sees the default that actually applies.

The store

from pymicroglia import store

source = store.fingerprint("VID52_C1_timestack.tif")   # 24 MB read, not 10.8 GB
store.put("registration", source, params, kind="table", value=shifts,
          name="registration_shifts_and_qc", output_dir=run_folder,
          method_version=METHOD_VERSION)

hit = store.get("registration", source, params, method_version=METHOD_VERSION)

A miss says why, because "cache miss" tells somebody to wait six hours and this tells them whether they meant to:

Cache miss for stage 'registration' on VID52_C1_timestack.tif.
Nearest stored artefact differs in:
  margin_px   stored 2   requested 8
Everything else matches, including the source and the METHOD_VERSION.
  stored at ...\AI_Exports\registered_2026-08-18\registration_shifts_and_qc.csv

resolve finds the one upstream artefact a stage should build on, so nobody has to remember which dated AI_Exports folder held the right registration. It never guesses: an explicit path wins, and otherwise exactly one match resolves. Two matches is an error naming both and the parameters that separate them.

Three classes of stored thing, and they are not equivalent:

Tier B moved into the project folder on 2026-08-20. It used to be forced local and a synced path was refused outright, on the grounds that syncing rebuildable pixels is a waste — which it is, and so is re-deriving them for hours on a machine that has never seen them. Two consequences worth knowing:

• Pin the folder. Dropbox frees space by turning a file into an online-only placeholder, and memory-mapping a 21 GB placeholder faults it back down the network a page at a time. Mark PixelStore "Make available offline"; doctor reports the count when it finds one and says the same thing.
• The index did not move. manifest.json is small, rewritten after every run, and written by both machines — exactly what a sync client turns into a pair of conflicted copies. It stays in %LOCALAPPDATA% and rebuilds by scanning.

Writes go to a memory-mapped .partial.npy renamed into place, so an interrupted run leaves nothing the next one will trust.

pymicroglia scan "...\AI_Exports\registered_2026-07-23"  # index one folder
pymicroglia store --evict                                  # back under the cap
pymicroglia verify raw.tif --full                          # whole-file hash, once

scan is not recursive, by rule: a recursive search over a synced path hydrates online-only files by the gigabyte.

Opening a time-lapse

from pymicroglia import open_series

stack = open_series("VID52_C1_phase-green-red_timestack.tif")
stack.shape                  # (1018, 3, 1024, 1024) — no pixel decoded
stack.meta.times_h[-1]       # 8.48, from the OME per-plane timestamps
stack.frame(0, 2)            # one image
stack.registered(0, 2)       # the same frame, with the stored transform
stack.crop(164, 212, 608, 360).frame(0, 2)      # a view; nothing moves

registered() is where the store and the pixels meet, and it is what makes re-analysis cheap: given a 298 KB shifts table, any frame of a 10.8 GB stack can be produced registered on demand. The 21 GB registered copy is an optimisation, not a prerequisite. If two registrations match, it raises rather than choosing.

Two questions the file cannot answer go through the decision class, so a person is asked at most once:

from pymicroglia import assign_channels, usable_window

assign_channels(stack)                                # inferred
assign_channels(stack, override="dluc=1,bf=0,struct=2")   # and remembered
usable_window(stack.meta.times_h)                     # longest unbroken block

io carries the three things that make file handling on this machine different: paths past 260 characters, handles held briefly by Dropbox and the antivirus, and writes that must never leave a truncated file at the target.

Registering

from pymicroglia import registration

registration.estimate_and_apply("VID52_C1_timestack.tif", registration_channel=1)
registration.estimate_and_apply_three_channel("VID52_B6_timestack.tif")
registration.export_registered_stack("VID52_B6_timestack.tif")   # finds the shifts

Two methods, kept apart because the engines differ for a reason. Reference builds a temporal reference from coarsely aligned frames and re-estimates every frame against it. Red sequential registers each frame against the one before, on the red neuronal channel, falling back to a segmented centroid when the correlation is weak — neurons hold still and microglia do not, and the microglia are what is being measured.

Both write registration_shifts_and_qc.csv and registration_summary.csv, under those names, in the engines' exact column order and fixed-precision formatting. Those filenames appear in a methods paragraph that ends up in a manuscript, and the numbers are compared to the engine's stored output cell by cell rather than with a tolerance.

estimate_only=True writes the two tables and no TIFF — the tables are the analysis; the registered stack is a 21 GB convenience that rebuilds from them. A second run with the same settings reads the stored shifts and does no estimation, which on a real recording is eleven minutes saved.

Filtering, and the line it must not cross

Two modules, because AGENTS.md draws a line here and a line inside one file is one an import can cross by accident.

raw --> registration --> cosmic-ray removal --> unsmoothed measurement
                              |
                              +--> display.py   (looking at only, never measured)

cosmic/ and filtering.py are the measurement branch. Both change pixels a number will later be computed from, so both write a record of what they changed and key it on the registration it came from:

from pymicroglia import cosmic, filtering

cleaned = cosmic.remove_cosmic_rays("VID52_registered.tif", signal_channel=2)
cleaned.mask.sum()                       # which pixel-frames were replaced
cleaned.summary["pixel_frames_replaced"]
filtering.unmix("VID52_registered.tif", unmixing_coefficient=0.025)

An outlier is a temporal event, never an absolute count. dluc_pipeline.py records why: a count threshold tuned on dim data deleted the brightest cell in a brighter dataset. Every cut is a multiple of the recording's own noise, so the same stack at half the gain has the same pixels repaired. seed_z carries that CAUTION unchanged — it is the setting that decides what counts as data.

One rule is asked three times: of a pixel, of a run of pixels along a line, and of a pixel at the top of the camera's range. That last one has not measured anything, so it is censored rather than replaced — reduced by a fitted bleed along its own row. mirror_placebo=True fits the same bleed on the side nothing bleeds toward and writes numbers only; it must come back near zero.

The dLuc pipeline reaches the same method through cosmic.clean_stack_in_place, because at that point what it has is a registered array and not a file. Same rule, same settings, same METHOD_VERSION — tests/test_cosmic_in_place.py is the assertion that the two agree pixel for pixel.

The matched-line method this replaced on 2026-08-20 is in superseded/matched_line.py. Nothing calls it; it stays importable as filtering.remove_cosmic_rays so that a run record written before that date still reproduces its run.

Unmixing is a step a caller opts into by name, never something a series does to itself on load: microglia_red_only_video_export.py unmixes and microglia_raw_registered_stack_export.py deliberately does not, and that difference between two exports of one recording has to stay visible.

display.py is the display branch, and its two methods are not interchangeable:

from pymicroglia import display

display.remove_static_background("cleaned.tif", frame_interval_h=0.5)
display.bioluminescence_display("cleaned.tif")   # writes *_DISPLAY_ONLY.tif

Every output of that branch carries _DISPLAY_ONLY in its name and display_only: true in its sidecar, because either signal alone can be lost — the flag survives a rename, the name survives a deleted sidecar. Hand one to a measurement and it stops:

DisplayOnlyInput: ..._display_DISPLAY_ONLY.tif is display-only and cannot feed
a measurement. See AGENTS.md: VSI conversion -> registration -> cosmic-ray
removal -> unsmoothed measurement, with display-only smoothing on a separate
branch. Display smoothing must never feed masks, traces, amplitudes or
statistics.

There is no force. A keyword that switched the guard off would be used once, in a hurry, and the resulting number would be indistinguishable from a real one.

Both display methods reproduce their engines' stacks bit for bit on a real recording — 241 frames of 504x504, every pixel, both display points to the last bit of a float, and every quality-control number identical. Cosmic-ray removal does the same against the engine's stored output, from a 77 KB fixture that ships with the package.

Finding cells

Segmentation's output is a label image and a few numbers per cell against a gigabyte-scale input, so it is keyed and stored like registration.

from pymicroglia import segmentation, roi

found = segmentation.segment("cleaned.tif", channels="dluc=0,struct=2")
len(found.cells), len(found.candidates)
found.reference.sigma            # background, read off tissue
roi.export_roi("cleaned.tif")    # RoiSet.zip Fiji can open

The accepted whole-SCN outline is also a named action, and it is the one whose algorithm is not in this package. Outlining a suprachiasmatic nucleus has nothing to do with microglia, so the roughly 6,200 lines that do it moved to PySCNSlice on 2026-08-23; pymicroglia.scn_outline is now a delegate, and the action needs pip install "PyMicroglia[scn]" to run. The interface below did not change.

Its input can be a registered red-channel time mean or a registered ImageJ/OME hyperstack. For a multi-channel hyperstack, scn_channel selects the outline channel using one-based ImageJ numbering. If there are multiple depth planes, scn_z selects the outline depth the same way. scn_time="mean" averages over time and remains the default; scn_time="max" uses a per-pixel maximum projection, while an integer such as scn_time=320 uses one-based source frame 320. The selected outline image determines one orientation and crop, which are then streamed across every frame, channel and depth plane. A sibling validfield_<key>.tif is found automatically for a named mean image; otherwise pass valid_mask.

from pymicroglia import scn_outline

result = scn_outline.automatic_scn_outline(
    "meanred_recording_01.tif",  # orientation + standard square crop are defaults
)
result["output"]          # full label: 0 background, 1 output-left, 2 output-right
result["oriented_source"] # full source image transformed into the same geometry
result["cropped_output"]  # standard square label crop around the SCN centre
result["cropped_source"]  # matching square source crop
result["report"]          # settings, transforms, crop geometry, QC and hashes

The equivalent multi-channel call is:

result = scn_outline.automatic_scn_outline(
    "registered_hyperstack.ome.tif",
    scn_channel=3,  # one-based: the third ImageJ channel
    scn_time=320,   # or "mean" (default) or "max"
    valid_mask="registered_hyperstack_validfield.tif",
)
result["outline_input"]       # exact mean, maximum projection or frame used
result["cropped_source"]      # all frames/channels, transformed and square-cropped

crop_mode accepts "tight", "standard", "wide", "custom" or "none". The three presets are relative to the smallest outline-centred square that contains the complete outline. For an exact output resolution, pass crop_size_px=512; the crop remains centred on the complete SCN outline and a size that would cut the outline is refused.

Pass orient_scn=False only when the original source pose is required. With orientation disabled, the port reproduces all ten original accepted labels byte for byte and matches the A007 A4/A5 crop-stability and carved-channel corrections. With orientation enabled, it reproduces all 16 approved A006 core masks and the approved source-image transform. A new acquisition still returns an open_questions check until its outline has been compared with hand-drawn evidence; passing through the same code is not evidence that the anatomy or image scale is the same.

Two of the five refusals in dluc_pipeline.py's header live here, and the code is shaped so they cannot be undone by accident:

• No minimum cell area. A 40 px floor once discarded a real 9.3-sigma cell with 30 px above threshold — it failed on size, not brightness. So there is no min_area parameter on segment_cells or anything it calls, and a test asserts the absence. RELAX_BELOW is also 40 and is the opposite setting: a component below it is regrown, not dropped.
• Background comes from the structural channel, never the image corners. The corners sit on the instrumental gradient, which inflated the noise estimate by 28 % and took the cell count from 5 to 4. Every mask is drawn at a multiple of that sigma, so a corner estimate shrinks all of them at once.

The tissue mask is checked rather than trusted: each candidate channel is scored by whether the bioluminescence is actually inside its mask, because the statistics that name the structural channel can pick the wrong one and nothing downstream notices.

A hand-drawn ROI is a decision, keyed on the source alone — it survives a parameter change, a METHOD_VERSION bump and a full cache eviction. scn_roi automates when nobody has answered, records that it did, and defers to the person the moment one exists. roi.py carries its own ImageJ .roi reader and writer, and reads the RoiSet.zip the existing pipeline wrote in August.

Traces, and what has to be true before a period is reported

from pymicroglia import tracing, controls, rhythm

traces = tracing.extract_traces("cleaned.tif")
controls.run_controls("cleaned.tif")          # decoys + the instrumental control
result = rhythm.test_rhythm("cleaned.tif")
result.periods[0]["period_hours"], result.control_passes

dF/F divides by each trace's window mean, not by its instantaneous rolling baseline. One reference cell's rolling baseline falls to 5 % of its mean, and the textbook definition turned that into a +634 % spike that flattened its panel. The rolling baseline is still subtracted — that is the detrend — but the divisor is one number per trace. tests/test_tracing_parity.py rebuilds that cell and asserts both halves: that the textbook form still blows up on it, and that this one does not.

Decoys go on tissue and the verdict is read in absolute counts. Off tissue their mask-minus-ring baseline is about zero, so dF/F divides by nothing and a real cell has to beat a distribution of near-singular ratios. Asking for off-tissue decoys raises; there is no flag.

No control, no rhythm result. In the reference dataset the structural channel showed a 22.8 h sinusoid at Lomb-Scargle power 0.966 — and the same rhythm was present off tissue where there is no sample, in a second channel, and in image sharpness. A daily focus cycle, not biology. So test_rhythm refuses a source with no instrumental-control artefact:

ControlMissing: No instrumental control for cleaned.tif. In the reference
dataset the structural channel showed a 22.8 h sinusoid at Lomb-Scargle power
0.966 — and the same rhythm was present off tissue where there is no sample, in
a second channel, and in image sharpness. It was a daily focus cycle, not
biology. Run controls.run_controls(source) first — its result is stored, so
this costs once per source.

The verdict then travels attached to the result, so no figure or record can show a period without also showing whether the control passed.

rhythm.py implements no period statistic. Every one comes from circadian_workbench, and a test asserts the absence by searching the module. Unlike analysis_kit, the workbench is a hard optional dependency — a missing audit layer costs a run record, a missing periodogram would cost a result — so its absence raises a named ImportError and pymicroglia doctor says up front whether rhythm analysis is available.

Figures, and the four rules that keep them small

PyFLASH/plotting.py is 31,497 lines — half that codebase — holding 620 functions, 557 of them private, to serve about forty public plots. Roughly 790 lines per plot, because each carried its own layout, saving, labelling and statistics helpers. Nobody decided that; it happened one function at a time. The same had started here, with figures at 1,042 of dluc_pipeline.py's 3,643 lines before this package existed.

So visualisation/ is built around four structural facts, each checked in tests/test_visualisation_limits.py:

The second is load-bearing. A figure that cannot compute cannot grow a private helper stack, because there is nothing for the helpers to do — and it means every figure's exact plotted data already exists as a table when it is saved. When the trace panel outgrew the cap, its computing half moved to pymicroglia/trace_tables.py. That is the cap working.

from pymicroglia import trace_tables

trace_tables.trace_panel(
    "traces_24h.csv",
    output_formats=("svg", "pdf", "png", "jpg", "tif", "webp", "avif", "heif"),
    dpi=300,
)  # twelve panels, one ReproFig identity across every format

Every figure arrives as a plot-that bundle without being asked:

traces_24h_panels.png
traces_24h_panels_plotted.csv        every value actually drawn
traces_24h_panels_provenance.json    sources with SHA256, settings, artefacts drawn
traces_24h_panels_bundle/
  README.md
  data/sources.csv  data/sources.md  data/src/   data/der/figure_data.csv
  fig/traces_24h_panels.svg          fig/preview.png

A source over 32 MB is hashed and listed but not copied — a bundle beside every figure drawn from a 10 GB acquisition documents nothing and fills a disk, and the row and the README both say so.

The trace panel is a port of trace_panel_figure.py, which has real users: a PowerShell wrapper, a saved JSON spec, and another folder's workflow contract that names it as a stage's review artefact. It reproduces that engine's output exactly — every plotted value bit-identical across all twelve panels, all twelve standard-deviation annotations, and 0 of 247,572 pixels different. The style port changed where a colour is written down, not what it is:

DEFAULT_COLOUR = "#a340d1"   ->  colour("dluc")
COLOUR_CYCLE   = [six hex]   ->  cycle("semantic")
UNDERLAY_GREY  = "0.72"      ->  GREY_LEVEL["raw"]

All three resolve to the values the engine spelled out. "0.72" stays a string because a Matplotlib grey level and a float mean different things.

theme defaults to "engine" on the trace panel, which leaves Matplotlib's defaults alone so a saved run does not move under somebody. Pass "pyflash" for the house look. Every other figure uses the house theme and passes the kit's own conformance check inside this project's test suite — which is where it belongs, because four projects each drew their own version of the house style and the colours came apart without anyone seeing it.

The quality-control figures read what a stage stored, never pixels they recompute:

A figure that needs a stored artefact and cannot find one says which action to run. cell_overlay will not segment anything for you: it draws what a stage found.

Videos, and the line they must not cross

A movie is something to look at, never something to measure. Every export records itself as a display-only artefact, so the same refusal that stops a display-filtered stack feeding a trace also stops a movie doing it — there is no flag that turns it off.

Five engines render movies and carry 91 of the project's parameters between them. Two of them exported a registered TIFF and a movie from one call, which is the conflation this stage undoes:

before                                after
  red_only_video_export(...)            filtering.unmix(...)      -> artefact
    unmixes                             registration.apply(...)   -> artefact
    applies shifts                      video.red_only(...)       -> MP4
    writes a registered TIFF stack
    renders a red-on-black MP4        the caller composes the three

unmixing_coefficient is still accepted and still recorded — a movie has to be able to say what it was drawn from — and is never applied. A test asserts that nothing under video/ writes a TIFF.

Playback speed is stated in experimental hours per second:

video.encode.frame_rate(frame_interval_h=0.5, hours_per_second=12)   # 24.0

At 12, one biological day takes two seconds of screen time whatever the acquisition interval was. The frame rate is derived; no frame is dropped or duplicated. A rate outside 1–60 fps is refused as a wrong speed rather than rendered.

Nothing shells out. Three of the four engines called ffmpeg through subprocess; this encodes through imageio-ffmpeg, which ships its own binary, so pymicroglia doctor reports video_export_available before a six-hour run rather than after it.

Parity has a method, because a movie has no bit-identical comparison. ffmpeg output is not reproducible across versions, so the test separates the two sources of difference:

They agree to three decimals, so the render is exact and the residual is the encoder. Every number the engines' manifests quote is reproduced exactly: the tissue mask at 93,447 pixels, the endpoint gain at 1.945895522388, the display ranges to nine figures.

Two findings worth knowing:

• There are two dluc_purple maps. tiff_stack_to_mp4.py has five stops, cry1_dluc_photon_pipeline.py has three. Merging them would repaint one set of movies, so they are dluc_purple and dluc_purple_photon.
• There are two Otsu thresholds, and both are right. Registration clips the histogram to the 1st–99.8th percentile so a cosmic ray cannot flatten it; the organotypic exporter uses the full range because its footage has no such outliers. Using one for both moved a display range by 0.04 counts.

Reading a parameter block

from pymicroglia import harvest, harvest_many

block = harvest("Protocols/Analysis/microglia_cosmic_ray_removal.py")
block.method_version                       # '2026-08-21-selectable-replacement'
len(block)                                 # 22
doc = block.by_constant("DEFAULT_SEED_Z")
doc.name, doc.type, doc.default, doc.units # ('seed_z', 'float', 12.0, '-')
doc.description                            # the full CAUTION note, joined

harvest reads a script's PROTOCOL PARAMETERS block into ParamDoc entries without importing or executing it — the real engines pull numpy, tifffile, scipy, scikit-image and matplotlib at module scope, and reading a comment should not cost that. Python, ImageJ macro and PowerShell blocks are all handled; seven of the thirteen protocols take parameters from a macro and two from PowerShell, so a Python-only reader would miss a quarter of the schema.

harvest_many merges the files that make up one protocol, first occurrence winning, which is how protocol-that's register.py builds INDEX.md — so the counts agree with it.

Tests

python -m pytest -q

Source-verification tests compare the reader against the real protocol scripts and skip when those are unavailable. Everything else runs from synthetic fixtures — text ones in tests/fixtures/, generated OME-TIFFs from tests/fixtures.py.

tests/test_real_stack.py is the exception. Three of stage 04's exit gates are claims about real data — how fast a ten-gigabyte stack opens, what a real VSI conversion's OME block contains, and whether a registered frame matches the stack the existing engine wrote — and none can honestly be closed against a 32-by-24 stand-in. Those tests skip until pointed at real files, because every acquisition here is a Dropbox online-only placeholder and reading one downloads 10.1 GB. Make one available offline and they run on their own:

python -m pytest tests/test_real_stack.py -v -rs

tests/test_registration_parity.py is the same idea for stage 05, and the reason that stage exists: it re-estimates registration on a real recording and compares the result to registration_shifts_and_qc.csv cell by cell, every frame, as text — not with a tolerance, because the engine writes fixed-precision numbers and so does this package. No engine is executed; the comparison is against output it produced earlier. It is off by default because it takes about ten minutes, and it caches the estimate afterwards so a re-check is seconds:

$env:PYMICROGLIA_PARITY = "1"; python -m pytest tests/test_registration_parity.py -v -s

Stage 06's two parity suites are cheaper and both run on their own. tests/test_filtering_parity.py holds the engine's own five unit cases, translated but otherwise unchanged, plus a comparison against a cleaned stack and replacement mask the engine wrote — all from tests/fixtures/cosmic_ray_reference/, 77 KB, so it needs nothing else present. tests/test_display_parity.py compares both display methods against the engine outputs in Tmem_Cry_leaktest_magenta_display_2026-08-17; it reads 344 MB and takes about 45 seconds, so it runs when those files are on the machine and skips with an explanation when they are not.

No test runs an engine. Every reference is output an engine produced earlier, which is the only kind of parity claim that still holds when Protocols/Analysis is absent.

License

PyMicroglia is released under the BSD 3-Clause License; the license text is included in every source and wheel distribution.

Citation

Until an archived release is available, cite the repository URL and the PyMicroglia version reported by pymicroglia doctor.