Putting two modalities in one coordinate frame

A CODEX panel and the H&E section beside it are two pictures of the same tissue in two coordinate systems. Every downstream claim — this marker is in that morphological compartment, this cell is inside that annotated region — is a claim about the map between them. Getting that map wrong does not produce an error. It produces a table.

That is what makes registration different from the steps around it. A segmentation that fails looks like nothing; a registration that fails looks like a slightly different answer. So the work splits in two, and the second half is the harder one: compute a transform, then prove it is the right transform without a ground truth to compare against.

ACCREDIT (Zhou, Zhao, Ren et al., OHSU) is the framework this page works through. It is MIT-licensed, runs six modality-specific pipelines behind one entry point — plus a senior variant for hard partial-overlap samples — scores every result with a reference-free composite, and escalates through a fixed ladder of recovery strategies when the score is low. The scoring is the interesting part and it is also where the traps are; the sections below are mostly about reading those numbers correctly.

This page is about the transform and the score. Getting the pixels and marker names out of a multiplex file is the multiplex-imaging-io skill; turning aligned pixels into a cell-by-marker table is codex-phenocycler. Both of those assume the frames already agree, which is what this page delivers.

What you have to obtain first

Nothing on the fluorescence-to-H&E path needs an account, a token or a GPU:

pip install numpy opencv-python tifffile imagecodecs SimpleITK scikit-image
git clone --depth 1 https://github.com/LeeZhou-bearway/ACCREDIT.git

imagecodecs is not in the project's own install line and you will need it. The H&E scans in the benchmark dataset used below are JPEG-compressed OME-TIFFs, and without it the pipeline's very first read fails:

ValueError: <COMPRESSION.JPEG: 7> requires the 'imagecodecs' package

openslide-python is needed only for the brightfield whole-slide path, which reads SVS through OpenSlide. The DAPI, CODEX and general-fluorescence pipelines import numpy, OpenCV, tifffile, skimage.color and SimpleITK and nothing else.

The H&E-to-IHC path additionally needs torch and LightGlue, and one of its two pretrained models is non-commercial. The project states this itself, in pipelines/he_ihc/METHODS.md, alongside the substitution:

| SuperPoint weights | MagicLeap research-only | No — use LightGlue's alternative | | LightGlue weights | Apache-2.0 | Yes |

For commercial deployment, replace SuperPoint with LightGlue's DISK backend (Apache-2.0) or use ALIKED (MIT).

So the restriction is on the feature extractor in one of six pipelines, not on the framework. Read it as three routes. Fluorescence to H&E — DAPI, CODEX, Xenium boundaries, spatial transcriptomics — never loads SuperPoint and is unrestricted. H&E to IHC as shipped loads SuperPoint and is academic-research-only. H&E to IHC with DISK or ALIKED substituted for the extractor is unrestricted and is the one to take for commercial work. The same file marks its reimplementation of the DeeperHistReg non-rigid core as "Check DHR repo (likely GPL-style) — Requires checking", which is upstream telling you the HE-IHC path has one more licence question open than the fluorescence path does. That is a second reason to lead with fluorescence where the science allows.

The rescue agent needs a metered API key that is billed separately from any subscription to a chat product. The project's README says an existing subscription covers it. It does not — API access is its own billed product, and a seat on a chat plan carries no API credit. Budget for per-token cost before you turn the rescue on. Set ANTHROPIC_API_KEY in the environment, or pass --api-key, which writes the same variable. You also need the anthropic package, which the top-level install line omits and which only pipelines/he_ihc/README.md mentions. Nothing about the deterministic pipelines requires any of this — the rescue is an optional branch and the sections below say exactly when it fires.

One representation both modalities can be compared in

Registration algorithms compare intensities. H&E is dark tissue on a light field; fluorescence is bright signal on a dark field. Nothing aligns until both are the same way up.

Two conversions, and they are not interchangeable:

import cv2, numpy as np
from skimage.color import rgb2hed

def he_inverted(rgb):
    """Whole-tissue proxy — everything stained shows up."""
    return (255 - cv2.cvtColor(rgb, cv2.COLOR_RGB2GRAY)).astype(np.uint8)

def he_hematoxylin(rgb):
    """Nuclear proxy — comparable to a DAPI channel. Normalise on percentiles."""
    h = rgb2hed(rgb.astype(np.float64) / 255.0)[:, :, 0]
    lo, hi = np.percentile(h, [1, 99])
    return (np.clip((h - lo) / (hi - lo + 1e-8), 0, 1) * 255).astype(np.uint8)

Inverted grayscale is a proxy for total optical density and responds to both stains; the hematoxylin channel isolates the nuclear component, so it is the side to compare a nuclear stain against. Measured on the pair in ## Try it, at the transform that block recovers, the effect is real and small:

moving channel	NCC vs inverted grayscale	NCC vs hematoxylin
fluorescence red	0.2546	0.2059
fluorescence green	0.2984	0.3060
fluorescence blue	0.3074	0.3150

The nuclear-looking channels prefer hematoxylin and the non-nuclear one prefers inverted grayscale, in the direction the mechanism predicts, by a few hundredths. Treat it as a tie-breaker once a registration is working, not as the thing that makes one work.

Normalise the hematoxylin channel on percentiles, not on min and max. ACCREDIT contains both — registration/preprocess.py uses the 1st and 99th percentiles, and run_fluo2he.he_to_h_channel uses h.min() and h.max(). On the brightfield thumbnail in ## Try it the two produce images with means of 42.1 and 3.4 out of 255; the min-max version puts 0.0% of pixels above intensity 32 against 35.4% for the percentile version, because one extreme pixel sets the whole range. It survives inside register_orb only because an Otsu threshold immediately follows it and recovers most of the mask (16.4% of the frame against 18.7%). Anywhere the H channel is used as an intensity rather than as input to a threshold, the min-max form is a bug waiting for a dust particle.

The fluorescence side needs the same treatment for the same reason — stretch between percentiles of the non-zero pixels, because a fluorescence frame is mostly true zeros and including them puts the 1st percentile at 0:

def stretch(plane):
    f = plane.astype(np.float32)
    pos = f[f > 0]
    lo, hi = np.percentile(pos, [1, 99])
    return np.clip((f - lo) / (hi - lo + 1e-10) * 255, 0, 255).astype(np.uint8)

Tissue masks come from Otsu on those representations, closed and dilated so that a villus lumen or a necrotic core does not read as background:

def tissue_mask(gray):
    _, m = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    k = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (15, 15))
    m = cv2.morphologyEx(m, cv2.MORPH_CLOSE, k, iterations=2)
    return cv2.morphologyEx(m, cv2.MORPH_DILATE, k, iterations=1)

The structuring element is 15 pixels at whatever resolution you hand it. ACCREDIT builds the fixed and moving masks with this same kernel at two different working scales, so the dilation is not the same physical distance on both. Resize both images to a common working size before masking, or the masks disagree by construction.

Never let it pick the modality for you

run_unified_agent.py detects the modality from file content and then, in single-pair mode, blocks on input() waiting for confirmation. An unattended agent hangs there forever. Pass -y, or use batch mode, which is non-interactive by construction.

More importantly, pass --mode explicitly, because the detector misroutes both of the commonest real inputs. Its fluorescence branch decides between the single-channel DAPI pipeline and the CODEX pipeline on len(tifffile.TiffFile(path).pages) > 5 — a page count, not a channel count. That distinction is the subject of the multiplex-imaging-io skill; here is what it costs. Four files, one detector:

file	pages	real channels	detected
one DAPI plane, plain 2-D grayscale TIFF	1	1	AxisError — crashes
3 fluorescence channels written as a 6-page pyramid or z-stack	6	3	codex
a 40-marker CODEX panel written as one 3-D page	1	40	dapi
brightfield H&E RGB named block_2.tif	1	3	he-ihc

Three of those four are wrong and the first does not survive at all. The crash is AxisError: axis 2 is out of bounds for array of dimension 2, raised inside the thumbnail loader on any 2-D image — which is the ordinary shape of a single-plane fluorescence export.

The last row is a substring match. ck is in the detector's IHC filename keyword list, for cytokeratin, and it is tested with in, so it matches inside block, stack, check and back. Rename that file section_2.tif and the same pixels are detected as he-he. Filename-driven pipeline selection is fine when the filenames were designed for it and catastrophic when they were not.

--mode takes boundary, dapi, codex, he-ihc, he-he, st-he or senior. Choose it from what you know about the acquisition, and let detection be a cross-check you read rather than a decision you delegate.

The transform ladder

Cross-modal pairs from adjacent sections or different instruments differ in orientation, scale, translation and — after that — in local deformation. Solve them in that order, because a refinement step cannot recover an orientation error.

Orientation and scale first, by search. There is no gradient to descend from a 180° error. ACCREDIT's DAPI pipeline enumerates 4 flips × 4 rotations × 6 scale ratios and keeps the best; the CODEX pipeline searches channels instead, scoring each one by cv2.matchTemplate with TM_CCOEFF_NORMED against the inverted H&E and registering the winner. Either way, the outer loop is exhaustive and the inner one is not.

Check where in the grid the winner landed. Those six ratios are literally [0.5, 0.7, 0.8, 1.0, 1.2, 1.5], applied after both images are scaled to their own working targets. The fine similarity fit that follows can move off the winning ratio, so the grid is not a hard bound — but it is the starting point for a local optimisation, and a winner sitting on the edge of the range means the optimum may be outside it. On the pair in ## Try it the pipeline's grid bottoms out at 0.5, the fine step settles at a full-resolution scale of 0.526, and an exhaustive full-resolution search finds 0.470 with 2.4× the structural correlation — 0.3151 against 0.1289. Nothing in the output says the winner was on a boundary. Compute the expected ratio from the two pixel sizes before you start, and pre-scale one image if it lands near or outside the ends of the range.

Translation by phase correlation, per candidate. cv2.phaseCorrelate on the two masks gives a translation in one FFT, so the search grid only has to cover scale and rotation.

Then a similarity or affine fit. ACCREDIT branches on whether the two images are within 20% of each other in size — ORB features on Otsu-binarised masks if so, SimpleITK Similarity2DTransform with Mattes mutual information if not.

That branch is worth knowing about, because the quantity it tests is not the quantity it is named after. res_ratio is the mean of the two per-axis pixel-dimension ratios, and the variable it sets is called same_resolution. On the pair in ## Try it the width ratio is 1.061 and the height ratio is 0.809 — a badly mismatched pair whose mean is 0.935, inside the ORB window. Two images at the same micrometres-per-pixel but different crops fail the test; two images at wildly different micrometres-per-pixel but similar pixel counts pass it. Neither image's actual pixel size is consulted. Read it off the file yourself.

Refine last, and check what the refinement actually applied. ecc_refine runs cv2.findTransformECC with MOTION_EUCLIDEAN, which estimates a rotation and a translation — and then adds only the two translation entries back into the transform. Feed it a pair differing by a 4° rotation and the linear part comes back untouched:

before: [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]
after : [[1.0, 0.0, -15.8562], [0.0, 1.0, 30.2891]]

That is a translation standing in for a rotation. It pulls the tissue centroid together and pushes the periphery further apart, and the tissue-overlap metrics reward it. Only accept a refinement that improved a metric sensitive to internal structure — see below.

The twelve numbers, and what each one is actually measuring

compute_quality_metrics_v2 returns twelve keys. Five of them, weighted, form the composite quality score (QCS); two more become multipliers on it; four are computed and never read; the twelfth is the composite itself.

key	what the name suggests	what it computes	in the composite
tissue_dice	the tissue overlaps	Dice of the two Otsu masks	weight 0.30
boundary_dice	the outlines agree	Dice restricted to a ±8 px band around either outline	weight 0.30
hausdorff_mean	Hausdorff distance	mean symmetric nearest-neighbour distance between the two largest contours	weight 0.15, as 1 − hd/HD_MAX
hausdorff_95	95th-percentile Hausdorff	95th percentile of the same distance set	unused
inlier_ratio	the images match	fraction of ORB matches kept by a freshly estimated RANSAC affine	weight 0.10, saturating at 0.15
n_inliers, n_matches	counts	counts	unused
ncc	the pixels agree	Pearson correlation inside the overlap	unused in the default composite
nmi	mutual information	(H(a)+H(b))/H(a,b) over pixels where both exceed zero	weight 0.15, saturating at 1.15
coverage	how much tissue is covered	warped mask area ÷ total frame pixels	multiplier min(1, cov/0.20)
bg_coverage	signal landing on background	warped mask ∩ H&E background ÷ H&E background	multiplier 1 − bg
qcs	the composite	the weighted sum × both multipliers	—

eval_px is a free parameter that moves the score, and it does not resize anything. It scales HD_MAX — 500 px at the default 4000 — and nothing else, so it is a promise about what resolution you warped at rather than an instruction. Pass the wrong one and the Hausdorff term silently changes. The same transform on the same pixels from ## Try it, whose images are 848 px on the long side:

eval_px	HD_MAX	Hausdorff sub-score	QCS
512	64	0.0000	0.5250
848, the true long side	106	0.3201	0.5729
4000, the default	500	0.8559	0.6531

A 0.13 swing on identical pixels, which is more than twice the gap between the pipeline's "accept this" threshold and its "call the rescue agent" threshold.

The DAPI pipeline does not pass its own working size. It builds its masks at WORK_SIZE = 2000 on the long side and calls the scorer with eval_px=4000, so HD_MAX is twice what it should be and the Hausdorff term is systematically over-generous. On the run in ## Try it that is not academic — the pipeline's own metrics rescored at 2000 give:

eval_px	QCS
4000, as the pipeline passes	0.6177
2000, its actual working size	0.5893

QCS_DONE is 0.60, so its own result moves from accepted to below threshold. Set eval_px to the long side you actually evaluated at, and do not compare a published HE-IHC score against a published DAPI one.

Four places the number is not what its name says

Each of these is reproduced by the Tier 0 block in ## Try it, which needs no download.

A perfect boundary match is scored as the worst possible one. The composite reads the Hausdorff term as metrics.get('hausdorff_mean', HD_MAX) or HD_MAX. Python's or treats 0.0 as absent, so the one value that means the contours coincide exactly is replaced by the value that means they could not be further apart. Holding everything else fixed:

hausdorff_mean	QCS
0.0	0.75
0.005	0.90
1.0	0.8977
500.0	0.75

The better of the two — the one whose contours coincide exactly — scores 0.15 lower than one that misses by a hundredth of a pixel. hausdorff_mean is rounded to two decimals before this, so anything under 0.005 px trips it, which in practice means an identity transform. That is exactly the self-registration sanity check people run first, and it reports 0.75 where the arithmetic says 0.90.

Only the largest contour of each mask is measured. Both Hausdorff values come from max(contours, key=len) on each mask, so a specimen with two tissue fragments is scored on one of them. Move the smaller fragment 340 px and leave the larger one identical: tissue_dice correctly drops to 0.9120 while hausdorff_mean and hausdorff_95 both report 0.0 — a perfect boundary match, which then hits the falsy-zero above and is scored as the worst one. Two defects in series, and the visible symptom is a good registration losing 0.15.

The inlier ratio does not measure your transform. orb_inlier_ratio detects ORB features in both images and runs cv2.estimateAffinePartial2D with RANSAC — estimating a new transform from scratch and reporting what fraction of matches it keeps. RANSAC absorbs whatever translation is present, so the number barely moves as you displace the result. Measured across the twelve shifted and rotated variants the ## Try it block scores — every shift from 0 to 160 px and every rotation up to 180° — it ranges from 0.326 to 0.565, and because it saturates at 0.15 its sub-score is 1.000 for all twelve. A tenth of the composite is a constant. On smooth synthetic images ORB finds no matches at all and it is a constant 0.000 instead. Neither constant is a measurement.

ncc is computed and then discarded. The module docstring calls normalized cross-correlation "most discriminating metric" and the default composite does not read it. Set it to −1.0, 0.0 or 1.0 with everything else fixed and the QCS is 0.6467 all three times. It is used only by composite_quality_score_cycif, the full-field-of-view variant, at weight 0.20 and normalised by a ceiling of 0.40 — so even there, the NCC of 0.3151 that the correct registration below achieves contributes 0.158 of 1.00, and an NCC of 0.11 contributes 0.055. If you want NCC to gate anything, read it directly.

The finding that matters — a silhouette is not an alignment

Four of the five weighted terms are computed on binary tissue masks. Masks describe the outline. Two outlines can agree while the structure inside them is displaced by tens of micrometres, which is the whole failure mode a cell-resolution claim depends on excluding.

Measured on the brightfield-versus-fluorescence pair in ## Try it — same section, two modalities, the fluorescence covering part of the brightfield field at a different scale. Three transforms, all scored at eval_px = 848, the images' true long side:

transform	scale	tissue_dice	ncc	QCS	what it looks like
the grid point with the best tissue_dice	0.860	0.8283	−0.0035	0.5620	structure visibly offset from structure
the grid point with the best ncc	0.470	0.5925	0.3151	0.5729	every crypt on its own crypt
ACCREDIT's dapi pipeline, 597 s	0.526	0.6760	0.1289	0.5951	correct, slightly loose

Two things to take from that table. The metric carrying the single largest weight ranks the visibly wrong transform first, by 40%. And the composite puts all three inside a band of 0.033, ordering them in a way uncorrelated with structural agreement — it correctly puts the worst last, and it inverts the top two, preferring the transform whose NCC is less than half the other's. NCC spans −0.0035 to 0.3151 over the same three.

Which grid point comes out on top depends on the search bounds, and so does the composite's ordering of the first two rows; narrow the grid and the wrong transform scores higher instead. The 2-to-3% margin is the finding. The ranking is not stable enough to be one.

The mechanism is not subtle once stated. The moving image covers only part of the fixed one, so Dice compares a small warped mask against the whole fixed tissue mask and every uncovered pixel counts against it. Scaling the moving image up covers more fixed tissue and raises Dice, while destroying the internal correspondence. Dice rewards covering more tissue. Correct registration does not always cover more tissue. ACCREDIT's coverage multiplier is aimed at the opposite failure — a warp that collapses — and does nothing about this one.

That mechanism also says where the problem stops, and the sweep at the end of ## Try it confirms it: crop the fixed image to the region the moving one covers, so both describe the same tissue, and tissue_dice on the same transform goes from 0.5925 to 0.9725 and the composite becomes decisive. The composite is valid on full-field pairs and unreliable on partial-field ones, and nothing in its output tells you which you have — which is what the null control below is for.

What the composite does do is separate gross failure from success. Displace the correct transform and it falls monotonically:

shift	0 px	5 px	10 px	20 px	40 px	80 px	160 px
QCS	0.5729	0.5317	0.4998	0.4608	0.4073	0.3236	0.1909
ncc	0.3151	0.1213	0.0391	−0.0494	−0.0280	0.0168	−0.0280

So it works as a catastrophe detector and not as a ranking. Notice where the pipeline's own thresholds sit against that table: QCS_DONE is 0.60, the rescue gate is 0.65 and the give-up threshold is 0.25. A 10-pixel displacement costs 0.073 — larger than the entire gap between "accept this" and "call the rescue agent". The operating thresholds are inside the noise band.

The null control — the cheapest honest test on this page

Score the registration you obtained. Then score the same registration with the moving image rotated 180°, and compare.

M180 = (np.vstack([M, [0, 0, 1]])
        @ np.vstack([cv2.getRotationMatrix2D((w / 2, h / 2), 180, 1.0), [0, 0, 1]]))[:2]

A 180° flip is the largest orientation error possible and the commonest gross failure in this domain. Whatever metric you are gating on must separate the answer from it by a wide margin. On the ## Try it pair:

	correct	rotated 180°	ratio
QCS	0.5729	0.4031	1.4×
ncc	0.3151	0.0255	12×
tissue_dice	0.5925	0.4754	1.2×

A compact specimen overlaps itself under rotation, so silhouette metrics cannot see orientation at all. Run this control on your pair before you trust any threshold on it. If your score ranks the answer within 1.5× of a 180° flip, it is not measuring alignment on your data and it must not be the gate. On a full-field version of the same pair the composite clears that bar by 90×, so this is a question about your framing rather than about the tool, and one warp answers it.

What to gate on instead

Read ncc on the overlap directly and require it to clear a floor you set from the null control. Report the composite alongside as context, not as the decision. The overlap mask is one line and the call is already in the library:

from quality_metrics import ncc, compute_quality_metrics_v2

overlap = ((fixed_mask > 0) & (warped_mask > 0)).astype(np.uint8) * 255
structure = ncc(fixed_inverted, warped_moving, mask=overlap)
metrics = compute_quality_metrics_v2(fixed_mask, warped_mask, fixed_inverted, warped_moving,
                                     eval_px=max(fixed_mask.shape))

Constrain the overlap when you search on NCC, or the optimiser will find a tiny high-correlation patch and stop. Requiring the overlap to be at least half the warped tissue's own area is enough, and it is what the ## Try it search does.

Two honest limits on this recommendation. NCC assumes a monotone intensity relationship between the modalities, which cross-modal pairs only approximately have — a realistic good value here is 0.32, not 0.9, and the floor has to be calibrated per pair rather than carried over. And nmi is the more principled cross-modal choice, except that the composite normalises it as (nmi − 1.0)/0.15 and saturates at 1.15, so anything genuinely informative is clipped to the same number; if you use NMI, use the raw value, not the sub-score.

The CODEX pipeline mis-scales the transform it scores

Verified against commit 9293473 on 2026-08-29, and worth re-checking before you read any CODEX QCS.

pipelines/codex/run_codex_pipeline.py holds two functions that convert a full-resolution transform down to the working resolution. _save_overlay multiplies the linear part by he_scale / fl_scale. _compute_qcs, in the same file, multiplies it by fl_scale / he_scale. They are reciprocals, so they cannot both be right, and the derivation says the overlay is: a point at he_scale · (A · p + t) in the working frame is (he_scale / fl_scale) · A · p_work + he_scale · t.

The consequence only appears when the H&E and the CODEX image differ in pixel dimensions, which is the normal case. Build a pair that is exactly a 2× rescaling of itself — a 6000 × 4500 H&E and a 3000 × 2250 CODEX — so the true transform is known to be [[2,0,0],[0,2,0]], and warp the nuclear mask both ways:

linear factor applied	mask Dice against the H&E
he_scale / fl_scale = 0.667, as _save_overlay uses	0.9609
fl_scale / he_scale = 1.500, as _compute_qcs uses	0.0822

Handed the exact correct transform, _compute_qcs returns QCS 0.1972. The same call with the ratio cancelled returns 0.5163. All six QCS evaluations in the CODEX decision cascade go through that function, and every threshold in it — 0.60 to proceed, 0.55, 0.35, 0.25 to give up — is read from that number.

On real data it changes the verdict, and the pipeline's own log shows it. Tier 2 of ## Try it runs the shipped CODEX pipeline on a 43-channel GBM panel against its name-matched H&E. At step 4 the region-of-interest module registers the pair, scores its own warp in its own frame, and prints:

  QCS=0.842  td=0.874  bd=0.874  ncc=0.555  hd=17.4
  Final dice: 0.8739  QCS: 0.8416
  ROI: QCS=0.2209

Two lines apart, one registration, 0.8416 and 0.2209. The first is the module scoring the warp it actually performed; the second is the parent re-scoring the same transform through _compute_qcs. The pipeline keeps the second, and four steps later prints STEP 6: All steps failed (QCS=0.2209 < 0.25) and writes a user_guidance.txt telling the reader to check orientation and scale by hand.

Putting the same transform back through the same function with the ratio cancelled gives QCS 0.6908, tissue_dice 0.9013, ncc 0.4475 — above the pipeline's own 0.60 acceptance threshold — and the overlay it wrote shows the CODEX signal sitting on the H&E tissue. The ROI module's own 0.8416 is itself inflated, by the eval_px mismatch described above, so 0.6908 is the number to believe. Both are more than three times the verdict the pipeline recorded.

Three practical consequences. Check whether it bites you: the applied factor is (fl_scale / he_scale)² times the correct one, where each scale is min(1, 4000 / long side). That is 1.0 in exactly two situations — the two long sides are equal, or both are at or below 4000 px so both scales clamp to 1.0. Whole-slide pairs are neither. In the synthetic pair above the factor is 2.25; where both images exceed 4000 px it is (L_fixed / L_moving)². Do not compare a CODEX QCS to a DAPI or HE-IHC one — they are not the same quantity. And score the CODEX result yourself with compute_quality_metrics_v2, warping the moving image into the fixed frame with your own scaling, which is what the ## Try it Tier 2 block does. The overlay image the pipeline writes is scaled correctly, so a good-looking overlay next to a poor QCS is the signature of this and not of a bad registration.

When the score is low

The CODEX pipeline escalates through a fixed ladder: register the best protein channel, then a nuclear-channel mutual-information rescue, then a rotation-and-phase-correlation rescue on the protein channel, then a region-of-interest registration with a channel and scale search, then the exhaustive DAPI orientation search as a fallback, and finally a user_guidance.txt file if the score is still under 0.25. Each step is entered on a threshold and the best result across all of them is kept, so a later step never makes things worse.

Above that sits an optional LLM rescue agent, gated in run_unified_agent.py on qcs < 0.65 and qcs > 0, for the dapi, codex, boundary and st-he modalities. Four things to know before you rely on it.

It cannot fire on a total failure. The qcs > 0 half of the condition means a registration scoring exactly zero — the worst outcome the pipeline can produce — skips the rescue entirely.

Its failures are swallowed. The call sits inside a bare except Exception that logs one line and continues. With no key, no package, or a stale model identifier, the run completes normally and returns the unrescued result. The only evidence is that line, which on a host without the SDK reads exactly:

Rescue agent failed: No module named 'anthropic'

Nothing about this reaches result.json or pipeline_summary.json. If you are consuming those files programmatically, grep the log for Rescue agent before concluding the rescue ran.

The model identifier is hard-coded, in two places in agent/agent.py. Read it out of that file rather than assuming it — model names are retired on the provider's schedule, and a stale one fails at the first API call and lands in the swallowed-exception path above.

It is a search, not an oracle. The agent selects among the same enhancement, rotation, channel-selection and re-registration tools the deterministic ladder already exposes, and it decides using the composite score. Everything in the sections above about that score applies to its decisions too. A rescue that improves the composite without improving the overlap NCC has improved the score, not the registration.

When a fetch returns 403 or 429, vary the request before concluding anything

This belongs in a registration page because assembling a cross-modal pair means fetching from several hosts, and a wrong conclusion here silently narrows what you will attempt.

assets.hubmapconsortium.org answers 403 to any User-Agent beginning curl/ or Wget/ and reaches the application for every other client string. Confirmed on 2026-08-29:

for ua in "curl/8.7.1" "Wget/1.21.4" "Mozilla/5.0"; do
  printf "%-16s " "$ua"
  curl -s -o /dev/null -w "%{http_code}\n" -I -A "$ua" "https://assets.hubmapconsortium.org/"
done
# curl/8.7.1       403
# Wget/1.21.4      403
# Mozilla/5.0      404

The 403 is a client-string block. The 404 is the application answering. Reported as-is, the first reads as "this data is not publicly accessible", which is false and would have been recorded as fact. A 429 from a data portal is usually a bot check on the landing page rather than a limit on the file host, and the same applies.

So: on 403 or 429, change the User-Agent, try the direct file URL rather than the landing page, and try one other client before writing anything down. Set -A on any curl aimed at a HuBMAP host. This is not about being clever with headers — it is that "I could not reach it" and "it is not available" are different findings, and only one of them is worth acting on.

What to write into the project

The transform is the artefact. Everything downstream is a re-derivation from it, so it has to be recoverable without re-running anything:

import json, numpy as np

np.save("moving_to_fixed_affine.npy", M)        # 2x3, maps moving pixels -> fixed pixels
json.dump({
    "fixed": {"path": fixed_path, "width": W, "height": H},
    "moving": {"path": moving_path, "width": w, "height": h},
    "transform_maps": "moving pixel coordinates -> fixed pixel coordinates",
    "search": {"scale": float(scale), "angle_deg": float(angle)},
    "metrics": {k: float(v) for k, v in metrics.items()},
    "ncc_overlap": float(structure),
    "ncc_overlap_rotated_180": float(structure_180),
}, open("registration.json", "w"), indent=2)

Alongside it, an overlay PNG. Every number on this page can be produced by a transform that is visibly wrong, and the two minutes it takes to look at the overlay is the only check that never misreports. Save it at a resolution where individual structures are distinguishable — a 2000 px long side is enough to see whether crypts land on crypts, and a 500 px thumbnail is not.

Write down which of the two H&E representations you used, which channel drove the registration, the search bounds, the working resolution, and the NCC floor you set from the null control. None of those is recoverable from the transform matrix, and each of them changes the answer.

Where this stops

Resolving a channel to a marker name, reading a region instead of a whole slide, and recovering the pixel size from a file's own tags are the multiplex-imaging-io skill. Choosing segmentation channels, running a segmenter over the aligned stack, measuring per cell and assembling the AnnData are codex-phenocycler. The spatial statistics that follow are spatial-phenotyping.

This page hands over one 2×3 matrix and an honest statement of how well it holds. It does not do non-rigid deformation beyond what the pipelines apply internally, it does not propagate a transform to segmentation masks or point sets, and it does not register three or more modalities into a common frame — for the last, register each moving image to the same fixed H&E rather than chaining, because chained transforms compound their errors and the composite score cannot see the compounding.

Try it

Three tiers. Tier 0 needs no data at all and reproduces the four metric defects. Tier 1 adds 0.2 MB and does a real cross-modal registration with the null control. Tier 2 is 849 MB and runs ACCREDIT's own CODEX pipeline on a published benchmark pair.

Data. Tier 1 is DAB-thumbnail.jpg and Fluo-thumbnail.jpg from Test Dataset for Whole Slide Image Registration (EPFL BIOP, Zenodo 10.5281/zenodo.5675686), CC BY 4.0, no account. They are one mouse duodenum section imaged twice — 848 × 484 brightfield with a DAB chromogen, and 799 × 598 fluorescence covering part of the same field at a different magnification. They are thumbnails, not the dataset: the record's 1,075 MB warpy-demo-project.zip holds the real slides. So Tier 1 exercises the representation, search, scoring and null control, and does not exercise OME-TIFF pyramid loading, channel selection or memory behaviour.

Tier 2 is ZH811_INF_v6.ome.tif (718 MB CODEX) and its 1:1 name-matched ZH811_INF_v6_HE.ome.tif (6.2 MB, inside the 131 MB HE_scans.zip) from Glioma spatialomics dataset (Zenodo 10.5281/zenodo.12624860), CC BY 4.0, no account. That is the CODEX-plus-H&E resource behind Greenwald, Galili-Darnell, Hoefflin et al., Cell 187(10) 2485–2501, 2024. Zenodo does not honour range requests on these files, so a dropped connection restarts the download from zero — use a retry loop rather than curl -C -. All four URLs confirmed reachable 2026-08-29.

Run. Cold, in an empty directory.

python3 -m venv .venv
. .venv/bin/activate
pip install numpy opencv-python tifffile imagecodecs SimpleITK scikit-image
git clone --depth 1 https://github.com/LeeZhou-bearway/ACCREDIT.git

Tier 0 — the twelve numbers and the four defects, no download

import sys, numpy as np, cv2
sys.path.insert(0, "ACCREDIT/quality")
from quality_metrics import (compute_quality_metrics_v2, composite_quality_score,
                             hausdorff_distance, tissue_dice)

H = W = 512
fixed = np.zeros((H, W), np.uint8)
cv2.ellipse(fixed, (256, 250), (170, 115), 20, 0, 360, 210, -1)
cv2.circle(fixed, (430, 440), 42, 210, -1)          # a second tissue fragment
fixed = cv2.GaussianBlur(fixed, (0, 0), 3)
fixed_mask = ((fixed > 0) * 255).astype(np.uint8)

def shift(img, d):
    return cv2.warpAffine(img, np.float32([[1, 0, d], [0, 1, d]]), (W, H))

print("displacement -> the twelve numbers, and the composite")
print(f"{'px':>4} {'dice':>7} {'bnd':>7} {'ncc':>8} {'hd_mean':>8} {'IR':>6} {'nmi':>7} {'cov':>7} {'bg':>7} {'QCS':>7}")
for d in (0, 5, 20, 60, 150):
    mv, mm = shift(fixed, d), shift(fixed_mask, d)
    m = compute_quality_metrics_v2(fixed_mask, mm, fixed, mv, eval_px=512)
    print(f"{d:>4} {m['tissue_dice']:>7.4f} {m['boundary_dice']:>7.4f} {m['ncc']:>8.4f} "
          f"{m['hausdorff_mean']:>8.2f} {m['inlier_ratio']:>6.3f} {m['nmi']:>7.4f} "
          f"{m['coverage']:>7.4f} {m['bg_coverage']:>7.4f} {m['qcs']:>7.4f}")
print("keys returned:", len(compute_quality_metrics_v2(fixed_mask, fixed_mask, fixed, fixed, eval_px=512)))

base = dict(tissue_dice=1.0, boundary_dice=1.0, inlier_ratio=0.0, nmi=2.0,
            coverage=0.5, bg_coverage=0.0)
print("\nA. a perfect boundary match is scored as the worst possible one")
for hd in (0.0, 0.005, 0.01, 1.0, 500.0):
    print(f"   hausdorff_mean={hd:<7} QCS={composite_quality_score({**base, 'hausdorff_mean': hd}, eval_px=512)}")

print("\nB. ncc is computed and then discarded by the default composite")
b = dict(tissue_dice=0.8, boundary_dice=0.6, hausdorff_mean=10.0, inlier_ratio=0.05,
         nmi=1.08, coverage=0.5, bg_coverage=0.02)
for v in (-1.0, 0.0, 1.0):
    print(f"   ncc={v:<5} QCS={composite_quality_score({**b, 'ncc': v}, eval_px=512)}")

print("\nC. only the largest contour is measured")
moved = fixed_mask.copy()
moved[380:500, 380:500] = 0                                    # delete the small fragment
cv2.circle(moved, (90, 450), 42, 255, -1)                      # and put it 340 px away
print(f"   tissue_dice          {tissue_dice(fixed_mask, moved):.4f}")
print(f"   hausdorff_mean, _95  {hausdorff_distance(fixed_mask, moved)}")

assert composite_quality_score({**base, 'hausdorff_mean': 0.0}, eval_px=512) \
     < composite_quality_score({**base, 'hausdorff_mean': 1.0}, eval_px=512)
assert (composite_quality_score({**b, 'ncc': -1.0}, eval_px=512)
        == composite_quality_score({**b, 'ncc': 1.0}, eval_px=512))
assert hausdorff_distance(fixed_mask, moved) == (0.0, 0.0) and tissue_dice(fixed_mask, moved) < 0.95
print("\nOK")

Expect — Tier 0. Invariants; a failure here means this page is wrong:

• Twelve keys come back from compute_quality_metrics_v2.
• hausdorff_mean = 0.0 scores strictly lower than hausdorff_mean = 1.0, everything else held fixed.
• The composite is identical for ncc = −1.0 and ncc = +1.0.
• The two-fragment mask reports (0.0, 0.0) for both Hausdorff values while tissue_dice is below 0.95 — the displaced fragment is invisible to one metric and visible to the other.
• inlier_ratio is 0.000 on every row. ORB finds no features on smooth synthetic shapes, so the 0.10 inlier term contributes nothing here — the mirror image of the real-data case below, where it is pinned at its ceiling instead.

Observed on 2026-08-29 with OpenCV 5.0.0 and numpy 2.5.2; these move if OpenCV changes its contour or ORB implementations:

  px    dice     bnd      ncc  hd_mean     IR     nmi     cov      bg     QCS
   0  1.0000  1.0000   1.0000     0.00  0.000  2.0000  0.2939  0.0000  0.7500
   5  0.9659  0.7845   0.7897     3.97  0.000  1.1986  0.2939  0.0142  0.8042
  20  0.8645  0.5289   0.2090    16.05  0.000  1.0335  0.2939  0.0564  0.5321
  60  0.6598  0.4834   0.0845    47.63  0.000  1.0200  0.2786  0.1271  0.3503
 150  0.3586  0.4615  -0.1121   103.21  0.000  1.0168  0.2191  0.1800  0.2155

The first two rows are the falsy-zero in the wild: a 5-pixel misalignment scores 0.8042 and perfect alignment scores 0.7500.

Tier 1 — a real cross-modal registration and the null control

curl -sL -o DAB-thumbnail.jpg  "https://zenodo.org/records/5675686/files/DAB-thumbnail.jpg?download=1"
curl -sL -o Fluo-thumbnail.jpg "https://zenodo.org/records/5675686/files/Fluo-thumbnail.jpg?download=1"

import sys, numpy as np, cv2
sys.path.insert(0, "ACCREDIT/quality")
from quality_metrics import compute_quality_metrics_v2, ncc, tissue_dice

def tissue(gray):
    _, m = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    k = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (15, 15))
    m = cv2.morphologyEx(m, cv2.MORPH_CLOSE, k, iterations=2)
    return cv2.morphologyEx(m, cv2.MORPH_DILATE, k, iterations=1)

def norm8(a):
    f = a.astype(np.float32); pos = f[f > 0]
    lo, hi = np.percentile(pos, [1, 99])
    return np.clip((f - lo) / (hi - lo + 1e-10) * 255, 0, 255).astype(np.uint8)

dab = cv2.cvtColor(cv2.imread("DAB-thumbnail.jpg"), cv2.COLOR_BGR2RGB)
flu = cv2.cvtColor(cv2.imread("Fluo-thumbnail.jpg"), cv2.COLOR_BGR2RGB)
H, W = dab.shape[:2]; h, w = flu.shape[:2]
print(f"fixed  brightfield DAB  {W}x{H}")
print(f"moving fluorescence     {w}x{h}")
print(f"per-axis size ratios: width {W/w:.3f} height {H/h:.3f} -> mean {(W/w + H/h)/2:.3f}")

he_inv = (255 - cv2.cvtColor(dab, cv2.COLOR_RGB2GRAY)).astype(np.uint8)
mv = norm8(cv2.cvtColor(flu, cv2.COLOR_RGB2GRAY))
he_m, mv_m = tissue(he_inv), tissue(mv)
cy, cx = h / 2, w / 2
mv_area = float((mv_m > 0).sum())

def build(s, th, dx, dy):
    M = cv2.getRotationMatrix2D((cx, cy), th, s); M[0, 2] += dx; M[1, 2] += dy; return M
def warp(M):
    return cv2.warpAffine(mv_m, M, (W, H)), cv2.warpAffine(mv, M, (W, H))
def score(M):
    wm, wf = warp(M)
    return compute_quality_metrics_v2(he_m, wm, he_inv, wf, eval_px=max(H, W))
def ncc_of(M):
    wm, wf = warp(M)
    ov = ((he_m > 0) & (wm > 0)).astype(np.uint8) * 255
    s2 = (M[0, 0] ** 2 + M[1, 0] ** 2)
    return ncc(he_inv, wf, mask=ov) if (ov > 0).sum() > 0.5 * mv_area * s2 else -1.0

cands = []
for s in np.arange(0.42, 0.90, 0.01):
    for th in np.arange(-6, 10.1, 2.0):
        R = cv2.getRotationMatrix2D((cx, cy), th, s)
        (dx, dy), _ = cv2.phaseCorrelate(np.float32(cv2.warpAffine(mv_m, R, (W, H))) / 255.,
                                         np.float32(he_m) / 255.)
        cands.append((s, th, dx, dy))

by_dice = max(cands, key=lambda p: tissue_dice(he_m, warp(build(*p))[0]))
by_ncc  = max(cands, key=lambda p: ncc_of(build(*p)))

print()
print(f"{'chosen by':>11} {'scale':>6} {'angle':>6} {'dice':>7} {'ncc':>8} {'IR':>6} {'nmi':>7} {'cov':>6} {'QCS':>7}")
for label, p in (("tissue dice", by_dice), ("ncc", by_ncc)):
    m = score(build(*p))
    print(f"{label:>11} {p[0]:>6.3f} {p[1]:>6.1f} {m['tissue_dice']:>7.4f} {m['ncc']:>8.4f} "
          f"{m['inlier_ratio']:>6.3f} {m['nmi']:>7.4f} {m['coverage']:>6.4f} {m['qcs']:>7.4f}")

M = build(*by_ncc); M3 = np.vstack([M, [0, 0, 1]])
print("\ndegrading the ncc-chosen transform by a pure shift")
print(f"{'px':>5} {'dice':>7} {'ncc':>8} {'QCS':>7}")
for d in (0, 5, 10, 20, 40, 80, 160):
    Md = M.copy(); Md[0, 2] += d; Md[1, 2] += d
    m = score(Md); print(f"{d:>5} {m['tissue_dice']:>7.4f} {m['ncc']:>8.4f} {m['qcs']:>7.4f}")

print("\nthe null control: rotate the moving image and re-score")
print(f"{'deg':>5} {'dice':>7} {'ncc':>8} {'QCS':>7}")
rot = {}
for th in (0, 10, 45, 90, 180):
    Mr = (M3 @ np.vstack([cv2.getRotationMatrix2D((cx, cy), th, 1.0), [0, 0, 1]]))[:2]
    rot[th] = score(Mr)
    print(f"{th:>5} {rot[th]['tissue_dice']:>7.4f} {rot[th]['ncc']:>8.4f} {rot[th]['qcs']:>7.4f}")

m0, m180 = rot[0], rot[180]
print(f"\nseparation between the answer and a 180-degree flip")
print(f"   QCS {m0['qcs']:.4f} / {m180['qcs']:.4f} = {m0['qcs']/m180['qcs']:.2f}x")
print(f"   NCC {m0['ncc']:.4f} / {abs(m180['ncc']):.4f} = {abs(m0['ncc']/m180['ncc']):.0f}x")

wf = warp(M)[1]; ovl = dab.copy(); k = wf > 30
ovl[k] = (ovl[k] * 0.55 + np.array([0, 255, 0], np.uint8) * 0.45).astype(np.uint8)
cv2.imwrite("overlay_ncc.png", cv2.cvtColor(ovl, cv2.COLOR_RGB2BGR))
np.save("M_moving_to_fixed.npy", M)
print("\nwrote overlay_ncc.png and M_moving_to_fixed.npy")

d_dice, d_ncc = score(build(*by_dice)), score(build(*by_ncc))
assert d_dice['tissue_dice'] > 1.3 * d_ncc['tissue_dice']          # dice prefers the wrong transform
assert abs(d_dice['qcs'] - d_ncc['qcs']) / d_ncc['qcs'] < 0.10     # the composite cannot separate them
assert d_ncc['ncc'] > 10 * abs(d_dice['ncc'])                      # ncc separates them by an order of magnitude
assert m0['ncc'] > 10 * abs(m180['ncc'])                           # ncc sees a 180-degree flip
assert m0['qcs'] / m180['qcs'] < 1.5                               # the composite barely does
assert m0['inlier_ratio'] > 0.15                                   # the inlier term sits on its ceiling
print("OK")

Expect — Tier 1. About ten seconds. Invariants:

• tissue_dice prefers the transform with essentially zero structural correlation, by more than 30%.
• The composite separates those two transforms by less than 10% — it cannot rank them.
• ncc separates them by more than an order of magnitude, and separates the answer from a 180° flip by more than 10×, while the composite stays under 1.5×.
• inlier_ratio exceeds 0.15 on the correct transform, so its sub-score is saturated at 1.000 and the 0.10 inlier term is constant.
• Open overlay_ncc.png. Every crypt cross-section in the green channel sits on its own brightfield counterpart. That is what an NCC of 0.32 looks like on a cross-modal pair, and it is the calibration the number has no meaning without.

Observed on 2026-08-29:

fixed  brightfield DAB  848x484
moving fluorescence     799x598
per-axis size ratios: width 1.061 height 0.809 -> mean 0.935

  chosen by  scale  angle    dice      ncc     IR     nmi    cov     QCS
tissue dice  0.860    8.0  0.8283  -0.0035  0.281  1.0222 0.4850  0.5620
        ncc  0.470    0.0  0.5925   0.3151  0.559  1.0584 0.2070  0.5729

degrading the ncc-chosen transform by a pure shift
   px    dice      ncc     QCS
    0  0.5925   0.3151  0.5729
    5  0.5783   0.1213  0.5317
   10  0.5609   0.0391  0.4998
   20  0.5304  -0.0494  0.4608
   40  0.4801  -0.0280  0.4073
   80  0.3947   0.0168  0.3236
  160  0.1888  -0.0280  0.1909

the null control: rotate the moving image and re-score
  deg    dice      ncc     QCS
    0  0.5925   0.3151  0.5729
   10  0.5698   0.0511  0.5089
   45  0.4997   0.0010  0.4089
   90  0.4792  -0.0205  0.3908
  180  0.4754   0.0255  0.4031

separation between the answer and a 180-degree flip
   QCS 0.5729 / 0.4031 = 1.42x
   NCC 0.3151 / 0.0255 = 12x

Here the composite happens to put the right transform first, by 0.0109. On a narrower search grid — 0.55 to 0.75 in steps of 0.005 — it puts the wrong one first instead, by 0.011. The margin is what holds; the ordering is not stable enough to be a finding, and a gate whose sign flips with a search bound is not a gate.

Across other inputs. The same reader was run against the ACCREDIT dapi pipeline end-to-end on this pair, and against a synthetic 6000 × 4500 H&E paired with a 3000 × 2250 CODEX stack:

input	what it exercised	outcome
this pair, through the shipped entry point in --mode dapi -y	the orientation search, ECC and BSpline steps	597 s on a 0.4-megapixel pair; scale 0.526, correct overlay; QCS reported as 0.6177 at eval_px=4000, 0.5893 at its true 2000 px working size, 0.5951 when its transform is rescored in the frame above; rescue gate fired at qcs < 0.65 and failed on the missing SDK
synthetic 2× rescaled pair, through _compute_qcs	the CODEX full-to-work transform scaling	true transform scored 0.1972; warping with the overlay function's ratio instead gives mask Dice 0.9609 against 0.0822
four TIFFs differing only in page layout and filename	detect_modality	one crash, two misroutes, one filename-driven misroute
a JPEG-compressed OME-TIFF H&E without imagecodecs	the first read in every pipeline	ValueError: <COMPRESSION.JPEG: 7> requires the 'imagecodecs' package

The three code defects generalised — the falsy-zero, the largest-contour restriction and the discarded NCC are properties of the source and appeared identically on synthetic and real data. The silhouette problem did not, and how it varies is the most useful thing in this section. Running the same transform and the same null control across four variants of the pair:

variant	tissue_dice	ncc	QCS	QCS vs its own 180° flip	NCC vs 180°
as published, 848 × 484	0.5925	0.3151	0.5729	1.42×	12×
both downsampled 2×, 424 × 242	0.5870	0.4374	0.5494	1.34×	57×
both downsampled 4×, 212 × 121	0.5556	0.4856	0.5281	1.17×	23×
moving cast to uint16	0.5925	0.3151	0.5729	1.42×	12×
fixed cropped to the moving's footprint	0.9725	0.3117	0.8146	90×	14×

Three conclusions, and the last one is the operative one:

Bit depth is a non-event. Percentile stretching over the non-zero pixels makes 8-bit and 16-bit inputs produce byte-identical scores. Nothing to guard against.

The composite discriminates worse at lower resolution. The 15-pixel morphology kernel, the 8-pixel boundary band and HD_MAX are all fixed pixel counts, so downsampling shrinks the real errors relative to them while leaving the tolerances put. NCC, having no length scale, sharpens instead. Do not validate a registration on a downsampled proxy.

The silhouette problem is a partial-field problem. Crop the fixed image to the region the moving one actually covers and both images describe the same tissue — tissue_dice goes to 0.9725 and the composite separates the answer from a 180° flip by 90× instead of 1.4×. So the composite is not broken; it is valid on full-field pairs and unreliable on partial-field ones, and nothing in its output says which you have. ## Try it's null control is how you find out, and it takes one warp.

One trap in the other direction. composite_quality_score_cycif exists for full-field pairs and is the right variant there — but on a partial-field pair its stricter coverage floor and area-ratio penalty both reward over-scaling. Feeding this pair to it, the score rises monotonically from 0.0055 at the correct transform to 0.1033 at 40% too large — nineteen times better for a registration whose NCC has collapsed from 0.3151 to 0.0614. Choose the variant from the framing, not from the modality name.

Tier 2 — ACCREDIT's own CODEX pipeline on a published benchmark pair

849 MB and roughly an hour of download at the rate Zenodo serves these files.

curl -sL --retry 5 --retry-delay 15 -o HE_scans.zip \
  "https://zenodo.org/records/12624860/files/HE_scans.zip?download=1"
unzip -q HE_scans.zip
until [ "$(stat -f%z ZH811_INF_v6.ome.tif 2>/dev/null || stat -c%s ZH811_INF_v6.ome.tif 2>/dev/null || echo 0)" -ge 718135498 ]; do
  curl -sL --retry 5 --retry-delay 15 -o ZH811_INF_v6.ome.tif \
    "https://zenodo.org/records/12624860/files/ZH811_INF_v6.ome.tif?download=1"
done

. .venv/bin/activate
cat > tier2.py <<'EOF'
import json, os, sys, numpy as np, cv2, tifffile
for sd in ("", "quality", "registration", "pipelines/codex",
           "pipelines/fluo_general", "pipelines/dapi"):
    sys.path.insert(0, os.path.join("ACCREDIT", sd))
from run_codex_pipeline import run_codex_pipeline, _compute_qcs
from quality_metrics import compute_quality_metrics_v2, ncc

HE, MV = "HE/ZH811_INF_v6_HE.ome.tif", "ZH811_INF_v6.ome.tif"
with tifffile.TiffFile(HE) as t:
    he_shape, he_levels = t.series[0].shape, len(t.series[0].levels)
with tifffile.TiffFile(MV) as t:
    mv_shape, mv_pages = t.series[0].shape, len(t.pages)
print("fixed  H&E  ", he_shape, "pyramid levels", he_levels)
print("moving CODEX", mv_shape, "pages", mv_pages)

he_long, mv_long = max(he_shape[:2]), max(mv_shape[-2:])
he_scale, mv_scale = min(1.0, 4000 / he_long), min(1.0, 4000 / mv_long)
print(f"long sides {he_long} and {mv_long}")
print(f"_compute_qcs applies {mv_scale/he_scale:.4f}; the correct factor is {he_scale/mv_scale:.4f}")
print(f"the scaling defect bites here: {abs(mv_scale/he_scale - he_scale/mv_scale) > 1e-9}")

res = run_codex_pipeline(HE, MV, "out_codex")
summary = json.load(open("out_codex/pipeline_summary.json"))
print("\nreported stage", summary["final_stage"], "reported QCS", summary["final_qcs"])

# Re-score the same transform with the correct full-to-work ratio.
M = np.array(res["M_full"] if "M_full" in res else json.load(open("out_codex/result.json"))["M_full"])
M_fixed = M.copy(); M_fixed[:, :2] *= (he_scale / mv_scale) / (mv_scale / he_scale)
qcs_corrected, m_corrected = _compute_qcs(HE, MV, M_fixed, "out_codex")
print("QCS with the ratio cancelled  ", qcs_corrected)
print("ncc on the overlap            ", m_corrected["ncc"])
print("tissue_dice                   ", m_corrected["tissue_dice"])
EOF
python tier2.py

Expect — Tier 2. About 18 minutes and a peak of roughly 9 GB of resident memory on ten cores; the 43-channel × 7-ratio scan in step 4 alone takes 694 s. The pipeline reads the whole 43-channel stack into memory, so budget for it. Nothing is printed until the process exits if you redirect stdout, because Python buffers it.

Invariants — a failure here means this page is wrong:

• The pipeline writes result.json, pipeline_summary.json, pipeline.log, user_guidance.txt and overlay_final.png into out_codex, and pipeline_summary.json carries final_stage, final_qcs and the six thresholds.
• The two images do not share a long-side pixel count, so the scaling defect applies: the block prints True for that line, and the two ratios it prints are reciprocals.
• The QCS printed with the ratio cancelled is more than three times the QCS the pipeline reports, and the ncc that comes with it is above 0.4 — a registration the pipeline declared a failure.
• overlay_final.png is scaled by the correct ratio and shows the CODEX signal filling the H&E tissue outline. Where the overlay looks right and the reported QCS looks bad, believe the overlay.

Observed on 2026-08-29, first and last lines of the run:

fixed  H&E   (10275, 8334, 3) pyramid levels 3
moving CODEX (43, 8885, 7194) pages 43
long sides 10275 and 8885
_compute_qcs applies 1.1564; the correct factor is 0.8647
the scaling defect bites here: True
...
STEP 6: All steps failed  (QCS=0.2209 < 0.25)
Total time: 1047s

reported stage failed reported QCS 0.2209
QCS with the ratio cancelled   0.6908
ncc on the overlap             0.4475
tissue_dice                    0.9013

Two more observed values worth knowing, both from the middle of the log. Channel selection scores all 43 channels by template match and picks ch39 at 0.5482; the region-of-interest step in the ladder independently picks ch31 by its own QCS scan. Different channel, same panel, two selection criteria — so "the best channel" is a property of the scoring rule, not of the panel, and it belongs in whatever you write down about the run.