Module hc_chisep 

HC-ChiSep: hollow-cylinder χ-separation with signal-derived fiber orientation.

HC-ChiSep separates paramagnetic (χ+, iron) and diamagnetic (χ−, myelin) susceptibility from a conventional QSM (χ_total), an R2’ map, and multi-echo GRE magnitude — deriving white-matter fibre orientation from the magnitude’s multi-compartment interference pattern (Wharton & Bowtell hollow-cylinder model) rather than requiring DTI.

The pipeline (headline mode) is:

1. Dr+ self-calibration — the paramagnetic relaxivity is estimated from the 5th percentile of R2'/χ_total over confidently-paramagnetic voxels, falling back to the field-scaled empirical 137·B0/3 Hz/ppm.
2. Closed-form two-source solve — χ+ = R2'/Dr+, χ− = χ+ − χ_total (Ridani convention: no diamagnetic relaxivity outside myelinated WM), clamped to keep χ+ ≥ 0, |χ−| ≥ 0. This is the non-WM branch.
3. WM-likeness (beat) weighting — a soft weight w from a hollow-cylinder vs mono-exponential model-selection test on the magnitude decay, gated by χ_total sign. Voxels whose magnitude shows the multi-compartment “beat” and are diamagnetic-leaning get w → 1.
4. θ + MWF grid fit — for supported voxels, fibre angle θ and myelin-water fraction (MWF) are found by an R2’-anchored grid search over the hollow- cylinder magnitude library, then spatially regularised.
5. Separation — in WM, |χ−| = K_χ·MWF (myelin-content ↔ MWF anchor, 0.038 ppm ↔ MWF 0.12), χ+ = χ_total + |χ−| shrunk toward a self-calibrated median; the two branches are blended by w.

All the reference’s env-gated optional regularisers (TV/guided joint inversion, R2’ denoise, hard constraint) are off in headline mode and are omitted here.

Inputs are in the project’s standard units — χ_total in ppm, R2’ in Hz, TE in seconds, B0 in Tesla — and outputs follow the χ-separation sign convention (χ+ ≥ 0, χ− ≤ 0, χ_total = χ+ + χ−).

Reference: Wharton, S. & Bowtell, R. (2012). “Fiber orientation-dependent white matter contrast in gradient echo MRI.” PNAS 109(45):18559-18564. HC-ChiSep is the QSM-CI submission building on this biophysical model.

Structs

AnchoredFitter 🔒

Precomputed hollow-cylinder magnitude library over the (θ, MWF) grid.

HcChisepParams

Parameters for hc_chisep.

Constants

CHI_A 🔒

CHI_GATE_C 🔒

CHI_GATE_S 🔒

CHI_I 🔒

CHI_NEG_REF 🔒

E_EXCH 🔒

F_AXON 🔒

GAMMA_BAR 🔒

Reduced gyromagnetic ratio (Hz/T) used by the hollow-cylinder model.

G_RATIO 🔒

K_CHI 🔒

K_χ: ppm of |χ−| per unit MWF (= CHI_NEG_REF / MWF_REF).

LAM 🔒

MWF_MAX 🔒

MWF_MIN 🔒

MWF_REF 🔒

NCONV_SIGMA 🔒

NM 🔒

NO_BEAT_FRAC 🔒

NT 🔒

SMOOTH_CLS 🔒

SMOOTH_FIT 🔒

T2_A 🔒

T2_E 🔒

T2_M 🔒

W_CENTER 🔒

W_SCALE 🔒

W_SUPPORT 🔒

Functions

finish 🔒

Emit the closed-form branch as the final result (χ+ ≥ 0, χ− ≤ 0 signed).

gather_normalised 🔒

Gather (len, E) first-echo-normalised signals for the given voxel indices from an echo-major volume stack stack[e*n + i].

gaussian_filter_3d 🔒

3D separable Gaussian filter, scipy gaussian_filter semantics: reflect boundary, truncate = 4.0.

hc_chisep

HC-ChiSep source separation from a QSM, an R2’ map and multi-echo magnitude.

hc_compartment_freqs 🔒

Compartment frequencies (Hz) for fibre-to-B0 angle theta (radians). Returns (Δf_myelin, Δf_axon, Δf_extra); extra-axonal is exactly 0.

hc_wm_r2prime 🔒

Mono-exp-equivalent reversible rate the beat itself contributes (Hz): weighted-LS slope of log(|GRE|/SE) vs centred TE, negated and clamped ≥ 0.

hc_wm_se_signal 🔒

Spin-echo WM magnitude factor (offsets refocused → real pool-T2 mixture).

hc_wm_signal_mag 🔒

|complex hollow-cylinder GRE magnitude| at echo time te (s), angle theta (rad), myelin-water fraction mwf, and mesoscopic reversible rate r2p_meso.

median 🔒

median_filter_3 🔒

3D median filter with a 3×3×3 window (scipy median_filter size 3, reflect).

mono_sse 🔒

Mono-exponential fit SSE for one voxel’s first-echo-normalised signal.

percentile 🔒

numpy-style linear-interpolation percentile (q in [0, 100]).

reflect_index 🔒

Reflect an index into [0, len) using scipy’s reflect mode (edge repeated).

smooth_stack 🔒

Smooth each echo of an echo-major stack stack[e*n + i] by a 3D Gaussian.

to_echo_major 🔒

Convert a voxel-major multi-echo stack vm[i*ne + e] to echo-major (contiguous per-echo volumes) em[e*n + i].

volmap 🔒

Scatter per-voxel values back into a full volume, filling the rest with fill.