diff --git a/ci/integration_osf.sh b/ci/integration_osf.sh index f089f45..fd0f7ef 100755 --- a/ci/integration_osf.sh +++ b/ci/integration_osf.sh @@ -23,12 +23,16 @@ curl -L --fail -o "$DATA/mono_t2.zip" "https://osf.io/kujp3/download?version=3" curl -L --fail -o "$DATA/mtr.zip" "https://osf.io/erm2s/download?version=2" curl -L --fail -o "$DATA/mtsat.zip" "https://osf.io/c5wdb/download?version=4" curl -L --fail -o "$DATA/vfa_t1.zip" "https://osf.io/7wcvh/download?version=3" +curl -L --fail -o "$DATA/b1_dam.zip" "https://osf.io/mw3sq/download?version=3" +curl -L --fail -o "$DATA/b1_afi.zip" "https://osf.io/csjgx/download?version=9" unzip -o -q "$DATA/ir.zip" -d "$DATA/ir" unzip -o -q "$DATA/qmt.zip" -d "$DATA/qmt" unzip -o -q "$DATA/mono_t2.zip" -d "$DATA/mono_t2" unzip -o -q "$DATA/mtr.zip" -d "$DATA/mtr" unzip -o -q "$DATA/mtsat.zip" -d "$DATA/mtsat" unzip -o -q "$DATA/vfa_t1.zip" -d "$DATA/vfa_t1" +unzip -o -q "$DATA/b1_dam.zip" -d "$DATA/b1_dam" +unzip -o -q "$DATA/b1_afi.zip" -d "$DATA/b1_afi" # Locate the datasets by their key files (robust to archive folder layout). IR_MAT="$(find "$DATA/ir" -name 'IRData.mat' | head -1)" @@ -82,6 +86,25 @@ VFA_REF_M0="$(find "$DATA/vfa_t1" -path '*FitResults*' -name 'M0.nii.gz' | head [ -n "$VFA_REF_T1" ] || { echo "FitResults/T1.nii.gz not found in VFA archive"; exit 1; } [ -n "$VFA_REF_M0" ] || { echo "FitResults/M0.nii.gz not found in VFA archive"; exit 1; } +# b1_dam ships as two separate 3D NIfTIs — one per flip angle, alpha and +# 2*alpha — and no mask, with a qMRLab reference (FitResults/B1map.nii.gz). +B1DAM_A="$(find "$DATA/b1_dam" -name 'SFalpha.nii.gz' | head -1)" +B1DAM_2A="$(find "$DATA/b1_dam" -name 'SF2alpha.nii.gz' | head -1)" +B1DAM_REF="$(find "$DATA/b1_dam" -path '*FitResults*' -name 'B1map.nii.gz' | head -1)" +[ -n "$B1DAM_A" ] || { echo "SFalpha.nii.gz not found in b1_dam archive"; exit 1; } +[ -n "$B1DAM_2A" ] || { echo "SF2alpha.nii.gz not found in b1_dam archive"; exit 1; } +[ -n "$B1DAM_REF" ] || { echo "FitResults/B1map.nii.gz not found in b1_dam archive"; exit 1; } + +# b1_afi ships as two separate 3D NIfTIs — one per interleaved repetition time +# — and no mask, with a qMRLab reference (FitResults/B1map_raw.nii.gz; the +# `_filtered` map is FilterClass smoothing, which is not part of this model). +B1AFI_TR1="$(find "$DATA/b1_afi" -name 'AFIData1.nii.gz' | head -1)" +B1AFI_TR2="$(find "$DATA/b1_afi" -name 'AFIData2.nii.gz' | head -1)" +B1AFI_REF="$(find "$DATA/b1_afi" -path '*FitResults*' -name 'B1map_raw.nii.gz' | head -1)" +[ -n "$B1AFI_TR1" ] || { echo "AFIData1.nii.gz not found in b1_afi archive"; exit 1; } +[ -n "$B1AFI_TR2" ] || { echo "AFIData2.nii.gz not found in b1_afi archive"; exit 1; } +[ -n "$B1AFI_REF" ] || { echo "FitResults/B1map_raw.nii.gz not found in b1_afi archive"; exit 1; } + echo "Running IR fit..." "$BIN" fit --mat-data "$IR_MAT" --mask "$IR_MASK" \ --config recipes/non-bids/irt1_config.yaml --output-dir "$DATA/out_ir" @@ -134,6 +157,29 @@ echo "Running vfa_t1 bidsify + BIDS-path fit..." "$BIN" fit --bids-dir "$DATA/vfa_t1_bids" \ --config recipes/bids/vfa_t1_config.yaml --output-dir "$DATA/out_vfa_t1" +echo "Running b1_dam bidsify + BIDS-path fit..." +# The two flip-angle volumes ship as separate 3D files, so --nii-data is +# repeated once per volume in acquisition order; the non-BIDS recipe supplies +# the flip angles, written into the TB1DAM sidecars. The BIDS-path fit folds +# those sidecars back in. +"$BIN" bidsify --model b1_dam \ + --nii-data "$B1DAM_A" --nii-data "$B1DAM_2A" \ + --config recipes/non-bids/b1_dam_config.yaml --subject 01 --out "$DATA/b1_dam_bids" +"$BIN" fit --bids-dir "$DATA/b1_dam_bids" \ + --config recipes/bids/b1_dam_config.yaml --output-dir "$DATA/out_b1_dam" + +echo "Running b1_afi bidsify + BIDS-path fit..." +# The two interleaved repetition times ship as separate 3D files, so +# --nii-data is repeated once per volume in acquisition order; the non-BIDS +# recipe supplies the repetition times and the nominal flip angle, written +# into the TB1AFI sidecars. The BIDS-path fit folds those sidecars back in and +# reassembles the acq-tr1/acq-tr2 collection. +"$BIN" bidsify --model b1_afi \ + --nii-data "$B1AFI_TR1" --nii-data "$B1AFI_TR2" \ + --config recipes/non-bids/b1_afi_config.yaml --subject 01 --out "$DATA/b1_afi_bids" +"$BIN" fit --bids-dir "$DATA/b1_afi_bids" \ + --config recipes/bids/b1_afi_config.yaml --output-dir "$DATA/out_b1_afi" + echo "Asserting outputs..." for f in "$DATA/out_ir/T1.nii.gz" "$DATA/out_ramani/F.nii.gz" "$DATA/out_srp/F.nii.gz" \ "$DATA/out_mono_t2/qmrust/sub-01/anat/sub-01_T2map.nii.gz" \ @@ -141,7 +187,9 @@ for f in "$DATA/out_ir/T1.nii.gz" "$DATA/out_ramani/F.nii.gz" "$DATA/out_srp/F.n "$DATA/out_mtsat/qmrust/sub-01/anat/sub-01_MTsat.nii.gz" \ "$DATA/out_mtsat/qmrust/sub-01/anat/sub-01_T1map.nii.gz" \ "$DATA/out_vfa_t1/qmrust/sub-01/anat/sub-01_T1map.nii.gz" \ - "$DATA/out_vfa_t1/qmrust/sub-01/anat/sub-01_M0map.nii.gz"; do + "$DATA/out_vfa_t1/qmrust/sub-01/anat/sub-01_M0map.nii.gz" \ + "$DATA/out_b1_dam/qmrust/sub-01/fmap/sub-01_TB1map.nii.gz" \ + "$DATA/out_b1_afi/qmrust/sub-01/fmap/sub-01_TB1map.nii.gz"; do test -s "$f" || { echo "MISSING or empty: $f"; exit 1; } done @@ -189,6 +237,26 @@ python3 ci/compare_maps.py \ "$DATA/out_vfa_t1/qmrust/sub-01/anat/sub-01_M0map.nii.gz" "$VFA_REF_M0" \ --rel-tol 0.001 --min-frac 0.999 --min-corr 0.999 --label vfa_t1-M0 +# b1_dam is a closed-form arc-cosine of a signal ratio, dimensionless in both +# implementations, so every voxel agrees to double-precision round-off — this +# dataset ships no mask, so that holds over the whole image including the +# out-of-domain voxels where the ratio exceeds 1 and both take the magnitude of +# the complex principal value. +echo "Comparing b1_dam B1 map to qMRLab FitResults..." +python3 ci/compare_maps.py \ + "$DATA/out_b1_dam/qmrust/sub-01/fmap/sub-01_TB1map.nii.gz" "$B1DAM_REF" \ + --rel-tol 0.001 --min-frac 0.999 --min-corr 0.999 --label b1_dam-B1 + +# b1_afi is likewise closed form and dimensionless in both implementations. Its +# estimator pins an unphysical signal ratio (r > 1) to B1 = 0 while leaving a +# non-finite ratio NaN, and reproducing that arithmetic exactly is what keeps +# the zero and NaN footprints equal to qMRLab's rather than merely the +# well-behaved voxels agreeing. +echo "Comparing b1_afi B1 map to qMRLab FitResults..." +python3 ci/compare_maps.py \ + "$DATA/out_b1_afi/qmrust/sub-01/fmap/sub-01_TB1map.nii.gz" "$B1AFI_REF" \ + --rel-tol 0.001 --min-frac 0.999 --min-corr 0.999 --label b1_afi-B1 + # The BIDS path and the .mat path must produce identical maps from the same # input. CLAUDE.md names this as the definition of a behaviour-preserving # change, and the two tests that assert it are `#[ignore]`d because they need diff --git a/crates/qmrust-cli/src/bidsify.rs b/crates/qmrust-cli/src/bidsify.rs index 92a73a6..7a4b2bd 100644 --- a/crates/qmrust-cli/src/bidsify.rs +++ b/crates/qmrust-cli/src/bidsify.rs @@ -24,7 +24,8 @@ pub struct BidsifyArgs { pub model: String, pub mat_data: Option, pub mat_dir: Option, - pub nii_data: Option, + /// Either one 4D NIfTI, or one 3D NIfTI per volume in acquisition order. + pub nii_data: Vec, pub nii_dir: Option, pub nii_mask: Option, pub mask: Option, @@ -54,7 +55,7 @@ pub fn run_bidsify(args: BidsifyArgs) -> Result<()> { // mask flag rather than silently ignoring it. let (data, mask, mut aux, source_header) = if let Some(dir) = args.nii_dir.as_ref() { anyhow::ensure!( - args.nii_data.is_none() && args.mat_data.is_none() && args.mat_dir.is_none(), + args.nii_data.is_empty() && args.mat_data.is_none() && args.mat_dir.is_none(), "--nii-dir is mutually exclusive with --nii-data/--mat-data/--mat-dir" ); anyhow::ensure!( @@ -62,7 +63,7 @@ pub fn run_bidsify(args: BidsifyArgs) -> Result<()> { "--mask is for a .mat source; pass --nii-mask with --nii-dir" ); read_named_nifti_source(dir, args.nii_mask.as_deref(), model.as_ref())? - } else if let Some(nii) = args.nii_data.as_ref() { + } else if !args.nii_data.is_empty() { anyhow::ensure!( args.mat_data.is_none() && args.mat_dir.is_none(), "--nii-data is mutually exclusive with --mat-data/--mat-dir" @@ -71,7 +72,7 @@ pub fn run_bidsify(args: BidsifyArgs) -> Result<()> { args.mask.is_none(), "--mask is for a .mat source; pass --nii-mask with --nii-data" ); - read_nifti_source(nii, args.nii_mask.as_deref(), model.as_ref())? + read_nifti_source(&args.nii_data, args.nii_mask.as_deref(), model.as_ref())? } else { anyhow::ensure!( args.nii_mask.is_none(), @@ -204,11 +205,34 @@ fn read_aux_args(specs: &[String], model: &dyn Model) -> Result, model: &dyn Model) -> Result { - let (data, header) = io::nifti::read_4d_nifti(nii, model.n_volumes())?; +/// Read a NIfTI measurement (+ optional NIfTI mask), preserving its spatial +/// header. `paths` is either one 4D file whose 4th axis is the acquisition +/// axis, or one 3D file per volume in acquisition order — the layouts qMRLab +/// demo datasets ship in. Auxiliary maps come from `--aux` for a NIfTI source; +/// there is no directory convention to discover them from. +fn read_nifti_source( + paths: &[PathBuf], + nii_mask: Option<&Path>, + model: &dyn Model, +) -> Result { + let (data, header) = match paths { + // One file: its own shape is ambiguous for a single-slice series, which + // is why the reader is told how many volumes the model expects. + [one] => io::nifti::read_4d_nifti(one, model.n_volumes())?, + many => { + // Per-volume files carry no acquisition axis of their own, so their + // order *is* the protocol order the recipe declares; a count + // mismatch means the two disagree and the sidecars would be + // misassigned. + anyhow::ensure!( + many.len() == model.n_volumes(), + "{} --nii-data volumes given, but the model's protocol has {}", + many.len(), + model.n_volumes() + ); + io::nifti::stack_nii_volumes(many)? + } + }; let mask = match nii_mask { Some(p) => Some(io::nifti::read_mask_nifti(p)?), None => None, @@ -314,13 +338,18 @@ pub fn write_bids_tree( write_participants_row(out, subject)?; sync_bidsignore(out, spec.suffix)?; - let anat_dir = out.join(format!("sub-{subject}")).join("anat"); - std::fs::create_dir_all(&anat_dir)?; + // The acquisition's datatype directory follows its BIDS suffix, so a + // transmit-field series (TB1*) lands in `fmap/` and a weighted series in + // `anat/` without either the model or this writer choosing. + let acq_dir = out + .join(format!("sub-{subject}")) + .join(rust_bids::datatype_for_suffix(spec.suffix)); + std::fs::create_dir_all(&acq_dir)?; // A NIfTI source's spatial header is preserved; a `.mat` source carries // none, so emit qMRLab's make_nii-compatible minimal header (matches how // `qmrust fit --mat-data` treats .mat input). The volume writers force // datatype 64 regardless, so voxel data is always written as f64. - let header = match source_header { + let header = match source_header.filter(|h| has_spatial_transform(h)) { Some(h) => h.clone(), None => make_minimal_header(nx, ny, nz), }; @@ -329,9 +358,9 @@ pub fn write_bids_tree( let vol = data.index_axis(Axis(3), i).to_owned(); let bv = model.bids_volume(i); let base = filename_stem(subject, &bv.entities, spec.suffix); - write_inv_volume(&vol, &header, &anat_dir.join(format!("{base}.nii.gz")))?; + write_inv_volume(&vol, &header, &acq_dir.join(format!("{base}.nii.gz")))?; std::fs::write( - anat_dir.join(format!("{base}.json")), + acq_dir.join(format!("{base}.json")), serde_json::to_string_pretty(&bv.sidecar)?, )?; } @@ -364,16 +393,20 @@ fn filename_stem(subject: &str, entities: &[(&'static str, String)], suffix: &st stem } -/// The datatype directory (`fmap` or `anat`) an auxiliary BIDS suffix belongs -/// in, per BIDS field-map convention: `*B1map`/`*B0map` suffixes are -/// estimated field maps (`fmap`); everything else (R1 maps, brain masks, …) -/// is an anatomical derivative (`anat`). -fn aux_datatype(suffix: &str) -> &'static str { - if suffix.ends_with("B1map") || suffix.ends_with("B0map") { - "fmap" - } else { - "anat" - } +/// Whether a source header actually defines where its voxels sit in space. +/// +/// NIfTI encodes the mapping twice, in the qform and the sform, each gated by +/// its own code; `0` in both is the ANALYZE-compatible "unknown" case, which +/// leaves `pixdim` as the only hint and the origin and orientation undefined. +/// Such a header has nothing to preserve, and propagating it writes a BIDS +/// dataset whose images cannot be placed — in practice a viewer falls back to +/// `pixdim`, which in an unset header is as arbitrary as the transform it +/// accompanies (qMRLab's own `b1_dam` example declares a 2 micrometre slice). +/// A synthesized minimal header is the honest substitute: it says "voxel grid, +/// unit spacing" rather than implying a geometry nobody recorded. Voxel data is +/// untouched either way. +fn has_spatial_transform(header: &NiftiHeader) -> bool { + header.qform_code > 0 || header.sform_code > 0 } /// The physical unit a known field-map BIDS suffix is conventionally @@ -392,7 +425,8 @@ fn aux_units(suffix: &str) -> Option<&'static str> { /// Write one auxiliary map into the `preprocessed` derivatives pipeline, /// under the datatype directory its BIDS suffix implies (see -/// `aux_datatype`), with a `Units` sidecar for field-map-datatype suffixes. +/// [`rust_bids::datatype_for_suffix`]), with a `Units` sidecar for +/// field-map-datatype suffixes. fn write_aux_map( out: &Path, subject: &str, @@ -400,7 +434,7 @@ fn write_aux_map( vol: &Array3, header: &NiftiHeader, ) -> Result<()> { - let datatype = aux_datatype(suffix); + let datatype = rust_bids::datatype_for_suffix(suffix); let dir = preprocessed_datatype_dir(out, subject, datatype)?; let base = format!("sub-{subject}_{suffix}"); write_inv_volume(vol, header, &dir.join(format!("{base}.nii.gz")))?; @@ -716,6 +750,42 @@ mod tests { let _ = std::fs::remove_dir_all(&dir); } + /// One 3D file per volume: the count is the whole contract, because a + /// per-volume file carries no acquisition axis to check against. Given the + /// wrong number, every sidecar after the first mismatch would describe a + /// different volume than the one beside it. + #[test] + fn a_per_volume_source_must_supply_exactly_the_protocol_s_volumes() { + let dir = tmp_dir("nii-count"); + let model = ir_model(&[0.350, 0.650, 0.950], None); + let paths: Vec = (0..2) + .map(|i| { + let p = dir.join(format!("v{i}.nii")); + let data = ndarray::Array3::::zeros((2, 2, 1)); + nifti::writer::WriterOptions::new(&p) + .write_nifti(&data) + .unwrap(); + p + }) + .collect(); + + // Two files, three inversion times. + let err = read_nifti_source(&paths, None, model.as_ref()) + .unwrap_err() + .to_string(); + assert!( + err.contains("2 --nii-data volumes given") && err.contains("has 3"), + "unhelpful count error: {err}" + ); + + // The matching count reads, and the volume axis is the file axis. + let three = ir_model(&[0.350, 0.650], None); + let (data, _, _, _) = read_nifti_source(&paths, None, three.as_ref()).unwrap(); + assert_eq!(data.dim(), (2, 2, 1, 2)); + + let _ = std::fs::remove_dir_all(&dir); + } + /// A NIfTI source's spatial header must be carried into the written volumes /// (the orientation-preservation contract). With `Some(source_header)`, the /// distinctive affine survives to the output; the `None` path (a synthesized @@ -753,6 +823,54 @@ mod tests { let _ = std::fs::remove_dir_all(&dir); } + /// A source header with neither transform code set defines no geometry at + /// all, so there is nothing to preserve and its `pixdim` is as arbitrary as + /// the missing transform. Carrying it through writes images a viewer cannot + /// place: qMRLab's own `b1_dam` example declares a 2 micrometre slice, which + /// renders as an empty sliver. The synthesized minimal header replaces it, + /// while the voxel data stays byte-identical. + #[test] + fn bidsify_replaces_a_source_header_that_declares_no_transform() { + let dir = tmp_dir("no-transform"); + let ir_data = Array4::from_shape_fn((2, 2, 1, 3), |(i, j, _k, t)| { + (i * 10 + j) as f64 + t as f64 * 0.5 + }); + let model = ir_model(&[0.350, 0.650, 0.950], None); + + // Exactly the shape of the qMRLab b1_dam source: no qform, no sform, + // and a placeholder pixdim describing a degenerate slab. + let mut src = make_minimal_header(2, 2, 1); + src.qform_code = 0; + src.sform_code = 0; + src.srow_x = [0.0; 4]; + src.srow_y = [0.0; 4]; + src.srow_z = [0.0; 4]; + src.pixdim = [0.0, 0.1406, 0.1406, 0.002, 0.0, 0.0, 0.0, 0.0]; + + write_bids_tree(model.as_ref(), &ir_data, None, &[], "01", &dir, Some(&src)).unwrap(); + + let nii = dir + .join("sub-01") + .join("anat") + .join("sub-01_inv-1_IRT1.nii.gz"); + let (data, out) = io::nifti::read_map_nifti_with_header(&nii).unwrap(); + assert!( + out.qform_code > 0 || out.sform_code > 0, + "written volume still declares no spatial transform" + ); + assert_eq!(out.pixdim[1], 1.0, "degenerate pixdim was carried through"); + assert_eq!(out.pixdim[2], 1.0); + assert_eq!(out.pixdim[3], 1.0); + // The geometry was substituted; the voxels were not. + for i in 0..2 { + for j in 0..2 { + assert_eq!(data[[i, j, 0]], (i * 10 + j) as f64); + } + } + + let _ = std::fs::remove_dir_all(&dir); + } + /// Re-running bidsify for the same subject must not duplicate the /// participants.tsv row. #[test] diff --git a/crates/qmrust-cli/src/catalog.rs b/crates/qmrust-cli/src/catalog.rs index 83a7263..30a07ec 100644 --- a/crates/qmrust-cli/src/catalog.rs +++ b/crates/qmrust-cli/src/catalog.rs @@ -20,8 +20,19 @@ pub struct ModelCard { pub name: String, pub bids_suffix: String, pub title: String, + /// Documentation directory slug. Several categories may share one. pub category: String, + /// Top-level heading a reader browses by. + pub family: String, + /// Lucide icon name for the family. + pub family_icon: String, + /// Heading within the family, when the family is subdivided. + pub subgroup: Option, + /// The most specific heading: the subgroup when there is one, else family. pub category_title: String, + /// Position in the taxonomy's reading order; sort by this and group + /// consecutive runs to rebuild the family/subgroup tree. + pub category_order: u8, pub summary: String, pub equation: String, pub symbols: Vec, @@ -79,6 +90,10 @@ pub struct Output { pub name: String, pub bids_suffix: Option, pub unit: Option, + /// Datatype directory this map is written into, from the same rule the + /// writers use. `None` for a diagnostic, which is not written as a + /// derivative at all. + pub datatype: Option, pub diagnostic: bool, /// The declared display window `[min, max]` in this output's unit, when the /// model declares one; `None` leaves the scale to the data. @@ -168,6 +183,7 @@ fn card(entry: &ModelEntry, repo_root: &Path) -> Result { Some((_, suffix, unit)) => Output { display_range: window(&name), name, + datatype: Some(rust_bids::datatype_for_suffix(suffix).to_string()), bids_suffix: Some(suffix.to_string()), unit: Some(unit.to_string()), diagnostic: false, @@ -177,6 +193,7 @@ fn card(entry: &ModelEntry, repo_root: &Path) -> Result { name, bids_suffix: None, unit: None, + datatype: None, diagnostic: true, }, }, @@ -249,7 +266,11 @@ fn card(entry: &ModelEntry, repo_root: &Path) -> Result { bids_suffix: entry.bids_suffix.to_string(), title: doc.title.to_string(), category: doc.category.slug().to_string(), + family: doc.category.family().to_string(), + family_icon: doc.category.icon().to_string(), + subgroup: doc.category.subgroup().map(str::to_string), category_title: doc.category.title().to_string(), + category_order: doc.category.order(), summary: doc.summary.to_string(), equation: doc.equation.to_string(), symbols: doc diff --git a/crates/qmrust-cli/src/commands.rs b/crates/qmrust-cli/src/commands.rs index e2e4917..ea26817 100644 --- a/crates/qmrust-cli/src/commands.rs +++ b/crates/qmrust-cli/src/commands.rs @@ -507,7 +507,8 @@ fn collection_sources(c: &Collection, bids_dir: &Path) -> Vec { /// Write `model`'s declared BIDS output maps (`Model::bids_outputs()`) from a /// fitted `results` into a BIDS-derivatives tree rooted at `deriv_root`: -/// `deriv_root/qmrust/[/]/anat/[_]_.nii.gz` +/// `deriv_root/qmrust/[/]//[_]_.nii.gz`, +/// where the datatype is the one that suffix implies (`rust_bids::datatype_for_suffix`) /// plus a full provenance JSON sidecar (`prov.sidecar(units)`) next to each. /// Outputs `bids_outputs()` doesn't declare (diagnostics like /// `res`/`idx`/`kf`/`resnorm`) are never written — @@ -546,9 +547,6 @@ fn write_derivatives( Some(ses) => qmrust_root.join(subject).join(ses), None => qmrust_root.join(subject), }; - let anat_dir = subject_dir.join("anat"); - std::fs::create_dir_all(&anat_dir)?; - let entity_stem = match session { Some(ses) => format!("{subject}_{ses}"), None => subject.to_string(), @@ -558,10 +556,15 @@ fn write_derivatives( let Some(map) = results.get(output_name) else { continue; }; + // Each output map lands in the datatype directory its own BIDS suffix + // implies, so a fitted transmit map (TB1map) is written where the next + // model's B1 input resolver already looks for one. + let dir = subject_dir.join(rust_bids::datatype_for_suffix(suffix)); + std::fs::create_dir_all(&dir)?; let base = format!("{entity_stem}_{suffix}"); - let nii_path = anat_dir.join(format!("{base}.nii.gz")); + let nii_path = dir.join(format!("{base}.nii.gz")); io::nifti::write_map_nifti(map, header, &nii_path)?; - let json_path = anat_dir.join(format!("{base}.json")); + let json_path = dir.join(format!("{base}.json")); std::fs::write( &json_path, serde_json::to_string_pretty(&prov.sidecar(units))?, @@ -848,6 +851,22 @@ fn load_aux_and_mask( bids_dir: &Path, ) -> Result<(AuxMaps, Option>, Vec)> { let paths = rust_bids::resolve_input_paths(table, model, identity, mask_spec)?; + // An input the recipe asked for and the dataset does not hold: the fit + // still runs, so this is said out loud rather than raised. Named by + // collection, since a run over many subjects would otherwise repeat an + // unattributable line per subject. + let who = ["subject", "session", "run"] + .iter() + .filter_map(|k| identity.get(*k).map(|v| format!("{k}-{v}"))) + .collect::>() + .join(" "); + for w in &paths.warnings { + if who.is_empty() { + eprintln!("warning: {w}"); + } else { + eprintln!("warning ({who}): {w}"); + } + } let mut maps: Vec<(String, Option>)> = Vec::new(); for (name, path) in &paths.aux { let map = match path { @@ -1493,7 +1512,7 @@ mod tests { model: "inversion_recovery".to_string(), mat_data: Some(ir_mat), mat_dir: None, - nii_data: None, + nii_data: Vec::new(), nii_dir: None, nii_mask: None, mask: Some(ir_mask), @@ -1604,7 +1623,7 @@ mod tests { model: "qmt_spgr".to_string(), mat_data: Some(qmt_mat), mat_dir: None, - nii_data: None, + nii_data: Vec::new(), nii_dir: None, nii_mask: None, mask: None, diff --git a/crates/qmrust-cli/src/io/nifti.rs b/crates/qmrust-cli/src/io/nifti.rs index 96f7810..7926a37 100644 --- a/crates/qmrust-cli/src/io/nifti.rs +++ b/crates/qmrust-cli/src/io/nifti.rs @@ -4,7 +4,7 @@ use anyhow::{bail, Context, Result}; use ndarray::{Array3, Array4}; use nifti::writer::WriterOptions; use nifti::{IntoNdArray, NiftiHeader, NiftiObject, ReaderOptions}; -use std::path::Path; +use std::path::{Path, PathBuf}; /// Read a NIfTI file and return the raw dynamic-dimension array + header. fn read_nifti_raw(path: &Path) -> Result<(ndarray::ArrayD, NiftiHeader)> { @@ -151,19 +151,32 @@ pub fn read_map_nifti_with_header(path: &Path) -> Result<(Array3, NiftiHead /// is preserved for the output geometry. Every role file must exist and share /// the same spatial dims. pub fn read_named_nii_volumes(dir: &Path, roles: &[&str]) -> Result<(Array4, NiftiHeader)> { - let mut vols: Vec> = Vec::with_capacity(roles.len()); + let paths: Vec = roles + .iter() + .map(|r| dir.join(format!("{r}.nii.gz"))) + .collect(); + stack_nii_volumes(&paths) +} + +/// Read one 3D scalar NIfTI per path and stack them into a 4D array in the +/// given order — column `i` is `paths[i]`. Used wherever a measurement ships +/// as one file per volume rather than one 4D file: a `Named` model's role +/// files, or a `Series` model's per-volume files in acquisition order. The +/// first file's spatial header is preserved for the output geometry. Every +/// file must exist and share the same spatial dims. +pub fn stack_nii_volumes(paths: &[PathBuf]) -> Result<(Array4, NiftiHeader)> { + let mut vols: Vec> = Vec::with_capacity(paths.len()); let mut header: Option = None; let mut dims: Option<(usize, usize, usize)> = None; - for &role in roles { - let path = dir.join(format!("{role}.nii.gz")); - let (v, h) = read_map_nifti_with_header(&path) - .with_context(|| format!("reading named role '{role}' from {:?}", path))?; + for path in paths { + let (v, h) = read_map_nifti_with_header(path) + .with_context(|| format!("reading volume from {path:?}"))?; let d = v.dim(); match dims { None => dims = Some(d), Some(expected) if expected != d => bail!( - "role '{}' has spatial dims {:?}, expected {:?} (from the first role)", - role, + "{:?} has spatial dims {:?}, expected {:?} (from the first volume)", + path, d, expected ), @@ -174,8 +187,8 @@ pub fn read_named_nii_volumes(dir: &Path, roles: &[&str]) -> Result<(Array4 } vols.push(v); } - let (nx, ny, nz) = dims.with_context(|| "a named model must declare at least one role")?; - let mut out = Array4::::zeros((nx, ny, nz, roles.len())); + let (nx, ny, nz) = dims.with_context(|| "at least one volume is required")?; + let mut out = Array4::::zeros((nx, ny, nz, paths.len())); for (t, v) in vols.iter().enumerate() { out.index_axis_mut(ndarray::Axis(3), t).assign(v); } @@ -300,6 +313,50 @@ mod tests { path } + /// Per-volume files carry the acquisition axis *between* them, so the + /// stack's 4th axis has to follow `paths` exactly. A silent reorder here + /// would pair every volume with the wrong sidecar. + #[test] + fn stacking_follows_the_order_the_paths_were_given() { + let dir = TempDir::new("stack-order"); + let a = write_nifti(&dir.0, "a.nii", &[2, 3, 1]); + let b = write_nifti(&dir.0, "b.nii", &[2, 3, 1]); + let (fwd, _) = stack_nii_volumes(&[a.clone(), b.clone()]).unwrap(); + let (rev, _) = stack_nii_volumes(&[b, a]).unwrap(); + assert_eq!(fwd.dim(), (2, 3, 1, 2)); + // Same two files, opposite order: every voxel swaps between volumes. + for x in 0..2 { + for y in 0..3 { + assert_eq!(fwd[[x, y, 0, 0]], rev[[x, y, 0, 1]]); + assert_eq!(fwd[[x, y, 0, 1]], rev[[x, y, 0, 0]]); + } + } + } + + /// Volumes of different sizes cannot be one series, and stacking them would + /// otherwise write whichever geometry came first over all of them. + #[test] + fn stacking_rejects_volumes_of_different_dimensions() { + let dir = TempDir::new("stack-dims"); + let a = write_nifti(&dir.0, "a.nii", &[2, 3, 1]); + let b = write_nifti(&dir.0, "b.nii", &[2, 4, 1]); + let err = stack_nii_volumes(&[a, b]).unwrap_err().to_string(); + assert!(err.contains("expected"), "unhelpful dims error: {err}"); + } + + /// The output geometry is the first volume's, matching what the 4D reader + /// returns for a single file. + #[test] + fn stacking_keeps_the_first_volumes_header() { + let dir = TempDir::new("stack-header"); + let a = write_nifti(&dir.0, "a.nii", &[2, 3, 1]); + let b = write_nifti(&dir.0, "b.nii", &[2, 3, 1]); + let (_, first) = read_map_nifti_with_header(&a).unwrap(); + let (_, stacked) = stack_nii_volumes(&[a, b]).unwrap(); + assert_eq!(stacked.srow_x, first.srow_x); + assert_eq!(stacked.pixdim, first.pixdim); + } + /// NIfTI's `dim[0]` counts dimensions, and a value of 3 declares three /// *spatial* axes — there is no temporal axis in the file at all. Reading /// the third one as time turns a 143x218x220 brain into one slice imaged diff --git a/crates/qmrust-cli/src/main.rs b/crates/qmrust-cli/src/main.rs index f7e6378..885b317 100644 --- a/crates/qmrust-cli/src/main.rs +++ b/crates/qmrust-cli/src/main.rs @@ -126,11 +126,12 @@ enum Commands { #[arg(long)] mat_dir: Option, - /// Path to a 4D NIfTI measurement (echoes/TIs in the 4th axis), for - /// datasets that ship as NIfTI rather than qMLab .mat. Mutually - /// exclusive with --mat-data/--mat-dir. + /// NIfTI measurement for a Series model, for datasets that ship as + /// NIfTI rather than qMRLab .mat: either one 4D file (echoes/TIs/flip + /// angles in the 4th axis) or, repeated, one 3D file per volume in + /// acquisition order. Mutually exclusive with --mat-data/--mat-dir. #[arg(long)] - nii_data: Option, + nii_data: Vec, /// Directory of per-role NIfTIs (`.nii.gz`, e.g. MTS's /// MTw/PDw/T1w) for a Named model. Mutually exclusive with the other diff --git a/crates/qmrust-core/src/models/b1_afi/config.rs b/crates/qmrust-core/src/models/b1_afi/config.rs new file mode 100644 index 0000000..9944ed5 --- /dev/null +++ b/crates/qmrust-core/src/models/b1_afi/config.rs @@ -0,0 +1,149 @@ +//! b1_afi config, parsed from the top-level YAML keys. +//! +//! Actual Flip-Angle Imaging has no fit options: its estimate is closed form. +//! The acquisition is the nominal flip angle plus the two interleaved +//! excitation repetition times. + +use anyhow::{bail, Result}; +use serde::{Deserialize, Serialize}; + +#[derive(Debug, Clone, Default, Deserialize, Serialize)] +pub struct B1AfiConfig { + /// The two excitation repetition times in seconds (BIDS/SI), in + /// acquisition order. Their ratio is what the estimate depends on. + #[serde(default)] + pub repetition_times: Vec, + /// Nominal excitation flip angle in degrees, shared by both volumes. + #[serde(default)] + pub flip_angle: Option, +} + +impl B1AfiConfig { + /// Config-intrinsic validation. + /// + /// The model has no fit options, but the shape of the acquisition array is + /// intrinsic: `describe` runs only this, and a half-stated protocol would + /// reach the fitter, which reads both repetition times by index. Empty is + /// legitimate (nothing has been resolved yet); anything other than the pair + /// the method is defined on is not. + pub fn validate_options(&self) -> Result<()> { + if !matches!(self.repetition_times.len(), 0 | 2) { + bail!( + "b1_afi is defined on exactly 2 repetition times [TR1, TR2], got {}", + self.repetition_times.len() + ); + } + // Stated positively so a non-finite value is rejected too: `tr <= 0.0` + // is false for NaN, which would otherwise fit an all-NaN map, or be + // written into a sidecar as `null` by a `bidsify` that never reaches + // `validate_protocol`. + if !self + .repetition_times + .iter() + .all(|&tr| tr > 0.0 && tr.is_finite()) + { + bail!( + "b1_afi repetition times must be finite and > 0 seconds, got {:?}", + self.repetition_times + ); + } + if let Some(fa) = self.flip_angle { + if !(fa > 0.0 && fa.is_finite()) { + bail!("b1_afi nominal flip angle must be a finite value > 0 degrees, got {fa}"); + } + } + Ok(()) + } + + /// Protocol-completeness validation: run once the acquisition is final + /// (from `--config` for non-BIDS, or composed from sidecars for BIDS). + pub fn validate_protocol(&mut self) -> Result<()> { + if self.repetition_times.len() != 2 { + bail!( + "b1_afi needs exactly 2 repetition times [TR1, TR2], got {}", + self.repetition_times.len() + ); + } + // The estimator divides by `n - r` with n = TR2/TR1; equal repetition + // times collapse it to 0/0 for every voxel, and the acquisition would + // not be an AFI pair in the first place. + if self.repetition_times[0] == self.repetition_times[1] { + bail!( + "b1_afi needs two distinct repetition times, both are {}", + self.repetition_times[0] + ); + } + // The value itself is checked by `validate_options`; what a resolved + // protocol adds is that there has to be one. + if self.flip_angle.is_none() { + bail!("b1_afi requires a nominal 'flip_angle' in degrees"); + } + Ok(()) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn parses_top_level_keys() { + let v: serde_yaml::Value = + serde_yaml::from_str("model: b1_afi\nrepetition_times: [0.02, 0.1]\nflip_angle: 60\n") + .unwrap(); + let mut cfg: B1AfiConfig = serde_yaml::from_value(v).unwrap(); + cfg.validate_options().unwrap(); + cfg.validate_protocol().unwrap(); + assert_eq!(cfg.repetition_times, vec![0.02, 0.1]); + assert_eq!(cfg.flip_angle, Some(60.0)); + } + + #[test] + fn validate_options_passes_without_an_acquisition() { + B1AfiConfig::default().validate_options().unwrap(); + } + + #[test] + fn validate_options_rejects_values_no_protocol_could_make_sense_of() { + // These are properties of the values themselves, so they are checked by + // the gate `describe` runs: `bidsify` never reaches `validate_protocol`, + // and would otherwise write a NaN into a sidecar as `null`, or describe + // a one-volume collection whose fitter reads two by index. + let cfg = |trs: Vec, fa: Option| B1AfiConfig { + repetition_times: trs, + flip_angle: fa, + }; + assert!(cfg(vec![0.02], Some(60.0)).validate_options().is_err()); + assert!(cfg(vec![0.0, 0.1], Some(60.0)).validate_options().is_err()); + assert!(cfg(vec![f64::NAN, 0.1], Some(60.0)) + .validate_options() + .is_err()); + assert!(cfg(vec![0.02, f64::INFINITY], Some(60.0)) + .validate_options() + .is_err()); + assert!(cfg(vec![0.02, 0.1], Some(f64::NAN)) + .validate_options() + .is_err()); + assert!(cfg(vec![0.02, 0.1], Some(0.0)).validate_options().is_err()); + // Nothing resolved yet is not an error: that is what `describe` sees. + assert!(cfg(vec![], None).validate_options().is_ok()); + } + + #[test] + fn validate_protocol_requires_two_distinct_times_and_a_flip_angle() { + let cfg = |trs: Vec, fa: Option| B1AfiConfig { + repetition_times: trs, + flip_angle: fa, + }; + assert!(cfg(vec![0.02], Some(60.0)).validate_protocol().is_err()); + assert!(cfg(vec![0.05, 0.05], Some(60.0)) + .validate_protocol() + .is_err()); + assert!(cfg(vec![0.02, 0.1], None).validate_protocol().is_err()); + // Either acquisition order is accepted; the fitter reads TR1/TR2 by + // value, not by position. + cfg(vec![0.1, 0.02], Some(60.0)) + .validate_protocol() + .unwrap(); + } +} diff --git a/crates/qmrust-core/src/models/b1_afi/fit.rs b/crates/qmrust-core/src/models/b1_afi/fit.rs new file mode 100644 index 0000000..ce03e27 --- /dev/null +++ b/crates/qmrust-core/src/models/b1_afi/fit.rs @@ -0,0 +1,222 @@ +//! Actual Flip-Angle Imaging (AFI) B1+ mapping — Yarnykh (2007). +//! +//! One spoiled gradient-echo sequence interleaves two excitation repetition +//! times, TR1 and TR2 = n·TR1, at a single nominal flip angle. In the pulsed +//! steady state the ratio of the two signals depends on the achieved flip +//! angle `θ` alone, so proton density and receive sensitivity cancel: +//! +//! ```text +//! r = |S(TR2) / S(TR1)| cos θ = (r·n − 1) / (n − r) +//! B1 = θ / θ_nominal +//! ``` +//! +//! `B1` is dimensionless: the achieved flip angle as a fraction of the nominal +//! one, the same convention the `B1map` aux input of other models expects. +//! +//! The closed form above assumes TR ≪ T1, while [`B1AfiFitter::forward`] uses +//! the exact steady-state signal, which also depends on T1. Recovery is +//! therefore not exact — the residual is the method's own systematic bias, and +//! grows as T1 falls toward TR (see `forward_then_fit_recovers_b1_to_within_the_methods_own_bias`). + +use crate::models::b1_afi::config::B1AfiConfig; + +/// Physically plausible transmit-field range, dimensionless. The estimate is +/// closed form and unconstrained; these bounds are descriptive only. +pub const B1_BOUNDS: (f64, f64) = (0.0, 2.0); + +/// Longitudinal relaxation time range in seconds, for the T1 the forward +/// signal depends on. Never fitted (see `Model::fixed_mask`). +pub const T1_BOUNDS: (f64, f64) = (0.0, 10.0); + +/// Pre-computed fitter for AFI data. +pub struct B1AfiFitter { + /// Excitation repetition times in seconds, in acquisition order. + repetition_times: Vec, + /// Nominal excitation flip angle in degrees. + flip_angle: f64, +} + +impl B1AfiFitter { + pub fn new(cfg: &B1AfiConfig) -> Self { + Self { + repetition_times: cfg.repetition_times.clone(), + // `validate_protocol` guarantees a nominal flip angle on the fit + // path. `describe` skips it to interrogate a model's structure + // before any sidecar is resolved, and reads neither `forward` nor + // `fit_voxel`. + flip_angle: cfg.flip_angle.unwrap_or(f64::NAN), + } + } + + /// Ground-truth parameters in `forward` order. T1 is not recovered by the + /// fit; it exists because the exact steady-state signal depends on it. + pub fn param_names() -> [&'static str; 2] { + ["B1", "T1"] + } + + pub fn output_names() -> [&'static str; 1] { + ["B1"] + } + + /// Excitation repetition times in seconds, in the order + /// `forward`/`fit_voxel` expect them. + pub fn repetition_times(&self) -> &[f64] { + &self.repetition_times + } + + /// The nominal flip angle the estimate is expressed relative to, degrees. + pub fn flip_angle(&self) -> f64 { + self.flip_angle + } + + /// The shorter repetition time, in seconds. Identified by value rather + /// than by position, so either acquisition order works and the BIDS + /// `acq-tr1`/`acq-tr2` labels always name the right volume. + pub fn tr1(&self) -> f64 { + self.repetition_times[0].min(self.repetition_times[1]) + } + + /// The longer repetition time, in seconds. + pub fn tr2(&self) -> f64 { + self.repetition_times[0].max(self.repetition_times[1]) + } + + /// Exact steady-state AFI signal at unit amplitude, one value per + /// repetition time in `repetition_times()` order — qMRLab's + /// `afi_equation`. `b1` scales the nominal flip angle; `t1` is in seconds. + pub fn forward(&self, b1: f64, t1: f64) -> Vec { + let theta = (b1 * self.flip_angle).to_radians(); + let (cos, sin) = (theta.cos(), theta.sin()); + let (tr1, tr2) = (self.tr1(), self.tr2()); + let (e1, e2) = ((-tr1 / t1).exp(), (-tr2 / t1).exp()); + let den = 1.0 - e1 * e2 * cos * cos; + // The volume acquired at TR1 leads with the *other* interval's + // recovery, and vice versa — the asymmetry the estimate reads. + let mz1 = ((1.0 - e2) + (1.0 - e1) * e2 * cos) / den * sin; + let mz2 = ((1.0 - e1) + (1.0 - e2) * e1 * cos) / den * sin; + self.repetition_times + .iter() + .map(|&tr| if tr == tr1 { mz1 } else { mz2 }) + .collect() + } + + /// Fit a single voxel from the signals in `repetition_times()` order. + /// Returns values in `output_names()` order. + pub fn fit_voxel(&self, signal: &[f64]) -> Vec { + let (s1, s2) = if self.repetition_times[0] == self.tr1() { + (signal[0], signal[1]) + } else { + (signal[1], signal[0]) + }; + let n = self.tr2() / self.tr1(); + vec![afi_b1(s2 / s1, n, self.flip_angle)] + } +} + +/// The AFI estimate for one voxel: `acos((r·n − 1)/(n − r)) / θ_nominal`, +/// dimensionless, from the signal ratio `r` (signed; magnitude is taken here) +/// and the repetition-time ratio `n`. +/// +/// A ratio above 1 is unphysical — the longer-TR volume cannot out-recover the +/// shorter one — so qMRLab attributes it to noise and pins the estimate at +/// zero. It expresses that as the arithmetic `cos_arg·[r ≤ 1] + 1·[r > 1]`, +/// which is reproduced literally below because its IEEE behaviour is part of +/// the contract: `Inf·0` and `NaN·0` are both NaN, so a voxel whose ratio is +/// non-finite (an empty denominator, or `r` landing exactly on `n`) stays NaN +/// instead of being pinned to zero. Writing the pin as an `if`/`else` clamp +/// would silently turn those voxels into `B1 = 0` and change the NaN +/// footprint of the output map. +fn afi_b1(ratio: f64, n: f64, nominal_flip_angle_deg: f64) -> f64 { + let r = ratio.abs(); + let indicator = |b: bool| if b { 1.0 } else { 0.0 }; + let cos_arg = (r * n - 1.0) / (n - r) * indicator(r <= 1.0) + indicator(r > 1.0); + cos_arg.acos().to_degrees() / nominal_flip_angle_deg +} + +#[cfg(test)] +mod tests { + use super::*; + + fn fitter(trs: Vec) -> B1AfiFitter { + B1AfiFitter::new(&B1AfiConfig { + repetition_times: trs, + flip_angle: Some(60.0), + }) + } + + #[test] + fn forward_then_fit_recovers_b1_to_within_the_methods_own_bias() { + // `forward` is the exact steady-state signal; the estimator assumes + // TR << T1. The gap is AFI's real systematic error, not a porting + // defect: it is under 1% here and shrinks as T1 grows. + let f = fitter(vec![0.02, 0.1]); + for &t1 in &[0.5, 0.9, 2.0] { + for &truth in &[0.6, 0.8, 1.0, 1.2, 1.4] { + let got = f.fit_voxel(&f.forward(truth, t1))[0]; + let rel = (got - truth).abs() / truth; + assert!(rel < 0.01, "B1={truth} T1={t1}: got {got} (rel {rel})"); + // The bias is one-signed: the estimate always undershoots. + assert!( + got <= truth, + "B1={truth} T1={t1}: got {got}, expected undershoot" + ); + } + } + } + + #[test] + fn the_bias_vanishes_as_repetition_times_shrink_against_t1() { + // Same B1, ever shorter TRs: recovery tends to exact, which is the + // TR << T1 limit the closed form is derived in. The error falls in + // proportion to TR/T1, so each tenfold shortening buys a decade. + let mut previous = f64::INFINITY; + let mut first = 0.0; + for (i, scale) in [1.0, 0.1, 0.01].into_iter().enumerate() { + let f = fitter(vec![0.02 * scale, 0.1 * scale]); + let err = (f.fit_voxel(&f.forward(1.2, 0.9))[0] - 1.2).abs(); + assert!(err < previous, "error did not shrink at scale {scale}"); + if i == 0 { + first = err; + } + previous = err; + } + assert!( + previous < first / 50.0, + "shortening TR 100x only improved {first} to {previous}" + ); + } + + #[test] + fn a_common_amplitude_cancels() { + let f = fitter(vec![0.02, 0.1]); + let plain = f.forward(0.9, 0.9); + let scaled: Vec = plain.iter().map(|s| s * 4321.0).collect(); + assert!((f.fit_voxel(&plain)[0] - f.fit_voxel(&scaled)[0]).abs() < 1e-12); + } + + #[test] + fn either_acquisition_order_gives_the_same_estimate() { + // TR1/TR2 are identified by value, so a series listed longest-first + // must not invert the ratio. + let ascending = fitter(vec![0.02, 0.1]); + let descending = fitter(vec![0.1, 0.02]); + let sig = ascending.forward(1.1, 0.9); + let flipped: Vec = sig.iter().rev().copied().collect(); + let a = ascending.fit_voxel(&sig)[0]; + let b = descending.fit_voxel(&flipped)[0]; + assert!((a - b).abs() < 1e-12, "{a} vs {b}"); + } + + #[test] + fn an_unphysical_ratio_is_pinned_to_zero_but_a_non_finite_one_stays_nan() { + // r > 1 cannot happen physically, so qMRLab reads it as noise and + // reports B1 = 0 rather than an out-of-domain arc cosine... + assert_eq!(afi_b1(1.5, 5.0, 60.0), 0.0); + // ...but only where the ratio is finite. `r` landing exactly on `n` + // makes the quotient infinite, and an empty voxel makes it NaN; both + // stay NaN, which is what keeps the NaN footprint equal to qMRLab's. + assert!(afi_b1(5.0, 5.0, 60.0).is_nan()); + assert!(afi_b1(f64::NAN, 5.0, 60.0).is_nan()); + assert!(afi_b1(f64::INFINITY, 5.0, 60.0).is_nan()); + } +} diff --git a/crates/qmrust-core/src/models/b1_afi/mod.rs b/crates/qmrust-core/src/models/b1_afi/mod.rs new file mode 100644 index 0000000..336417f --- /dev/null +++ b/crates/qmrust-core/src/models/b1_afi/mod.rs @@ -0,0 +1,7 @@ +//! Actual Flip-Angle Imaging B1+ mapping. BIDS suffix: `TB1AFI`. + +pub mod config; +pub mod fit; +pub mod model; + +pub use model::{build, describe, dump, effective}; diff --git a/crates/qmrust-core/src/models/b1_afi/model.rs b/crates/qmrust-core/src/models/b1_afi/model.rs new file mode 100644 index 0000000..654b030 --- /dev/null +++ b/crates/qmrust-core/src/models/b1_afi/model.rs @@ -0,0 +1,386 @@ +//! b1_afi adapter onto the core `Model` trait. +//! +//! A `Series` of exactly two spoiled gradient-echo volumes acquired at +//! interleaved excitation repetition times, combined by the closed-form AFI +//! ratio — no iterative fit. BIDS indexes the pair by the `acq` entity, whose +//! values name the shorter and longer repetition time (`tr1`/`tr2`). + +use crate::core::model::{ + Aux, BidsSpec, BidsVolume, EntityRole, FitStrategy, InputSpec, Measurement, MeasurementKind, + Model, ProtoParam, Protocol, Sample, Scope, Source, +}; +use crate::models::b1_afi::config::B1AfiConfig; +use crate::models::b1_afi::fit::{B1AfiFitter, B1_BOUNDS, T1_BOUNDS}; +use anyhow::Result; +use serde_json::json; +use std::collections::BTreeMap; + +pub struct B1AfiModel { + fitter: B1AfiFitter, +} + +const B1AFI_ENTITIES: &[EntityRole] = &[EntityRole::Other("acq")]; + +/// One `{"RepetitionTimeExcitation": s}` identity row per volume, in +/// acquisition order. +fn b1_afi_rows(fitter: &B1AfiFitter) -> Vec> { + fitter + .repetition_times() + .iter() + .map(|&tr| BTreeMap::from([("RepetitionTimeExcitation".to_string(), tr)])) + .collect() +} + +impl B1AfiModel { + pub fn new(cfg: B1AfiConfig) -> Self { + Self { + fitter: B1AfiFitter::new(&cfg), + } + } +} + +impl Model for B1AfiModel { + fn param_names(&self) -> Vec<&'static str> { + B1AfiFitter::param_names().to_vec() + } + fn output_names(&self) -> Vec { + B1AfiFitter::output_names() + .iter() + .map(|s| s.to_string()) + .collect() + } + fn param_bounds(&self) -> Vec<(f64, f64)> { + vec![B1_BOUNDS, T1_BOUNDS] + } + fn fixed_mask(&self) -> Vec { + // T1 is a ground-truth parameter of the forward signal that the closed + // form never recovers — it assumes TR << T1 and divides it out. + vec![false, true] + } + fn required_inputs(&self) -> Vec { + // B1+ is what this model measures; it consumes no auxiliary map. + vec![] + } + fn measurement(&self) -> MeasurementKind { + MeasurementKind::Series { + rows: b1_afi_rows(&self.fitter), + } + } + fn strategy(&self) -> FitStrategy { + FitStrategy::Voxelwise + } + fn forward(&self, params: &[f64], _aux: &Aux) -> Measurement { + let values = self.fitter.forward(params[0], params[1]); + let samples = self + .fitter + .repetition_times() + .iter() + .zip(values) + .map(|(&tr, value)| Sample { + params: BTreeMap::from([("RepetitionTimeExcitation".to_string(), tr)]), + value, + }) + .collect(); + Measurement::Series(samples) + } + fn fit(&self, m: &Measurement, _aux: &Aux) -> Vec { + // Assemble in the fitter's own repetition-time order by matching each + // expected time to its sample by identity — never positionally. + // Swapping the two volumes would otherwise invert the ratio. + let samples = m.series(); + let signal: Vec = self + .fitter + .repetition_times() + .iter() + .map(|&tr| { + samples + .iter() + .find(|s| s.params.get("RepetitionTimeExcitation") == Some(&tr)) + .map(|s| s.value) + .unwrap_or_else(|| { + panic!("measurement has no sample with RepetitionTimeExcitation={tr}") + }) + }) + .collect(); + self.fitter.fit_voxel(&signal) + } + fn n_volumes(&self) -> usize { + self.fitter.repetition_times().len() + } + fn bids_volume(&self, index: usize) -> BidsVolume { + // The BIDS spec names the shorter repetition time `tr1` and the longer + // `tr2`, so the label follows the value, not the acquisition order. + let tr = self.fitter.repetition_times()[index]; + let label = if tr == self.fitter.tr1() { + "tr1" + } else { + "tr2" + }; + BidsVolume { + entities: vec![("acq", label.to_string())], + sidecar: BTreeMap::from([ + ("RepetitionTimeExcitation".to_string(), json!(tr)), + ("FlipAngle".to_string(), json!(self.fitter.flip_angle())), + ]), + } + } + fn bids(&self) -> Option { + Some(BidsSpec { + suffix: "TB1AFI", + entities: B1AFI_ENTITIES, + }) + } + fn protocol_schema(&self) -> Vec { + // The repetition time is the acquisition axis, so it alone identifies a + // volume: a `Series`' per-volume protocol rows *are* its volume + // identities (`engine::build_volume_ids`), and `forward` tags its + // samples the same way. The nominal flip angle is one value for the + // pair, so it is `Global` — as a per-volume param it would join the + // identity and no forward sample would match its volume. + vec![ + ProtoParam { + name: "RepetitionTimeExcitation", + source: Source::Field("RepetitionTimeExcitation"), + scope: Scope::PerVolume, + required: true, + }, + ProtoParam { + name: "FlipAngle", + source: Source::Field("FlipAngle"), + scope: Scope::Global, + required: true, + }, + ] + } + fn bids_outputs(&self) -> Vec<(&'static str, &'static str, &'static str)> { + // Dimensionless: the achieved flip angle as a fraction of the nominal + // one, which is the scaling every `B1map` aux consumer expects. + vec![("B1", "TB1map", "")] + } +} + +impl crate::core::model::ModelConfig for B1AfiConfig { + const NAME: &'static str = "b1_afi"; + const SUBKEY: Option<&'static str> = None; + const PROTOCOL_KEYS: &'static [&'static str] = &["repetition_times", "flip_angle"]; + + fn validate_options(&mut self) -> Result<()> { + B1AfiConfig::validate_options(self) + } + + fn ingest_protocol(&mut self, proto: &Protocol) -> Result<()> { + let times: Vec = proto + .volumes + .iter() + .filter_map(|m| m.get("RepetitionTimeExcitation").copied()) + .collect(); + if !times.is_empty() { + self.repetition_times = times; + } + // Global scope, so this is resolved once for the collection rather than + // per volume — read it independently of the per-volume rows. + if let Some(&fa) = proto.global.get("FlipAngle") { + self.flip_angle = Some(fa); + } + Ok(()) + } + + fn validate_protocol(&mut self) -> Result<()> { + B1AfiConfig::validate_protocol(self) + } + + fn into_model(self) -> Box { + Box::new(B1AfiModel::new(self)) + } +} + +/// Structural interrogation entry point (see [`describe_model`](crate::core::model::describe_model)). +pub fn describe(v: &serde_yaml::Value) -> Result> { + crate::core::model::describe_model::(v) +} + +/// Registry builder (see [`build_model`](crate::core::model::build_model)). +pub fn build(v: &serde_yaml::Value, proto: &Protocol) -> Result> { + crate::core::model::build_model::(v, proto) +} + +/// Registry dumper (see [`dump_model`](crate::core::model::dump_model)). +pub fn dump(v: &serde_yaml::Value) -> Result { + crate::core::model::dump_model::(v) +} + +/// Registry option-surface entry point (see +/// [`effective_model`](crate::core::model::effective_model)): every option this +/// model accepts, at its effective value, plus any validation complaint. +pub fn effective( + v: &serde_yaml::Value, + proto: &Protocol, +) -> Result { + crate::core::model::effective_model::(v, proto) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn b1_afi_value() -> serde_yaml::Value { + serde_yaml::from_str("model: b1_afi\nrepetition_times: [0.02, 0.1]\nflip_angle: 60\n") + .unwrap() + } + + /// A resolved protocol as the shell hands it over: one row per volume + /// carrying the identifying axis, plus the collection-wide flip angle. + fn resolved_proto(times: &[f64], flip_angle: f64) -> Protocol { + Protocol { + volumes: times + .iter() + .map(|&tr| BTreeMap::from([("RepetitionTimeExcitation".to_string(), tr)])) + .collect(), + global: BTreeMap::from([("FlipAngle".to_string(), flip_angle)]), + } + } + + #[test] + fn build_and_roundtrip_via_trait() { + let m = build(&b1_afi_value(), &Protocol::default()).unwrap(); + assert_eq!(m.param_names(), vec!["B1", "T1"]); + assert_eq!(m.output_names(), vec!["B1".to_string()]); + assert_eq!(m.n_volumes(), 2); + // The estimator assumes TR << T1, so recovery carries the method's own + // bias rather than being exact (see the fitter's own tests). + let sig = m.forward(&[0.85, 0.9], &Aux::new()); + let fitted = m.fit(&sig, &Aux::new()); + assert!((fitted[0] - 0.85).abs() / 0.85 < 0.01, "B1: {}", fitted[0]); + } + + #[test] + fn t1_is_a_forward_parameter_the_fit_never_recovers() { + let m = build(&b1_afi_value(), &Protocol::default()).unwrap(); + assert_eq!(m.fixed_mask(), vec![false, true]); + // It is real, not decorative: the forward signal moves with it. + let a = m.forward(&[1.0, 0.5], &Aux::new()); + let b = m.forward(&[1.0, 2.0], &Aux::new()); + assert!((a.series()[0].value - b.series()[0].value).abs() > 1e-6); + // ...and it is absent from the fitted outputs. + assert_eq!(m.fit(&a, &Aux::new()).len(), 1); + } + + #[test] + fn fit_assembles_by_identity_not_position() { + let m = build(&b1_afi_value(), &Protocol::default()).unwrap(); + let sig = m.forward(&[0.85, 0.9], &Aux::new()); + let mut reversed: Vec = sig + .series() + .iter() + .map(|s| Sample { + params: s.params.clone(), + value: s.value, + }) + .collect(); + reversed.reverse(); + assert_eq!( + m.fit(&sig, &Aux::new()), + m.fit(&Measurement::Series(reversed), &Aux::new()) + ); + } + + #[test] + #[should_panic(expected = "no sample with RepetitionTimeExcitation")] + fn fit_panics_on_unmatched_identity() { + let m = build(&b1_afi_value(), &Protocol::default()).unwrap(); + let bogus = Measurement::Series(vec![Sample { + params: BTreeMap::from([("RepetitionTimeExcitation".to_string(), 7.0)]), + value: 1.0, + }]); + let _ = m.fit(&bogus, &Aux::new()); + } + + #[test] + fn bids_folds_the_acquisition_from_protocol_and_labels_by_repetition_time() { + let proto = resolved_proto(&[0.03, 0.15], 55.0); + // Config carries no acquisition; the sidecars supply it. + let v: serde_yaml::Value = serde_yaml::from_str("model: b1_afi\n").unwrap(); + let m = build(&v, &proto).unwrap(); + assert_eq!(m.n_volumes(), 2); + let short = m.bids_volume(0); + assert_eq!(short.entities, vec![("acq", "tr1".to_string())]); + assert_eq!(short.sidecar["RepetitionTimeExcitation"], json!(0.03)); + assert_eq!(short.sidecar["FlipAngle"], json!(55.0)); + assert_eq!(m.bids_volume(1).entities, vec![("acq", "tr2".to_string())]); + } + + #[test] + fn the_tr_labels_follow_the_value_not_the_acquisition_order() { + // A series listed longest-first must still label the shorter time + // `tr1`, as the BIDS spec defines it. + let proto = resolved_proto(&[0.15, 0.03], 55.0); + let v: serde_yaml::Value = serde_yaml::from_str("model: b1_afi\n").unwrap(); + let m = build(&v, &proto).unwrap(); + assert_eq!(m.bids_volume(0).entities, vec![("acq", "tr2".to_string())]); + assert_eq!(m.bids_volume(1).entities, vec![("acq", "tr1".to_string())]); + } + + #[test] + fn forward_samples_carry_the_volume_identities_the_bids_path_builds() { + // A `Series` model's volume identities come from the resolved + // per-volume protocol (`engine::build_volume_ids`), while `forward` + // tags its samples from the model's own rows. Any protocol param the + // model does not also emit joins the identity on one side only, and + // every predicted sample then fails to match its volume — the fit + // still works (it queries one key), so the only symptom is a silently + // missing forward curve. Both sides must agree exactly. + let proto = resolved_proto(&[0.02, 0.1], 60.0); + let v: serde_yaml::Value = serde_yaml::from_str("model: b1_afi\n").unwrap(); + let m = build(&v, &proto).unwrap(); + + let ids = crate::engine::build_volume_ids(m.measurement(), &proto, m.n_volumes()).unwrap(); + let sig = m.forward(&[0.9, 0.9], &Aux::new()); + let samples = sig.series(); + assert_eq!(samples.len(), ids.len()); + for (id, sample) in ids.iter().zip(samples) { + let crate::core::model::VolumeId::Params(row) = id else { + panic!("b1_afi is a Series model; expected param-row identities") + }; + assert_eq!( + *row, sample.params, + "volume identity {row:?} has no matching forward sample identity {:?}", + sample.params + ); + } + } + + #[test] + fn declares_bids_tb1afi_and_no_aux() { + let m = build(&b1_afi_value(), &Protocol::default()).unwrap(); + assert_eq!(m.bids().unwrap().suffix, "TB1AFI"); + assert!(m.required_inputs().is_empty()); + } + + #[test] + fn bids_outputs_reference_real_output_names() { + let m = build(&b1_afi_value(), &Protocol::default()).unwrap(); + let names = m.output_names(); + for (out, _suffix, _unit) in m.bids_outputs() { + assert!(names.iter().any(|n| n == out), "{out} not in {names:?}"); + } + } + + #[test] + fn describe_succeeds_without_an_acquisition_and_exposes_schema() { + let v: serde_yaml::Value = serde_yaml::from_str("model: b1_afi\n").unwrap(); + let m = describe(&v).unwrap(); + let schema = m.protocol_schema(); + assert_eq!(schema[0].name, "RepetitionTimeExcitation"); + assert_eq!(schema[1].name, "FlipAngle"); + // The repetition time identifies a volume; the flip angle is one value + // for the collection. + assert!(matches!(schema[0].scope, Scope::PerVolume)); + assert!(matches!(schema[1].scope, Scope::Global)); + } + + #[test] + fn build_still_requires_an_acquisition_when_protocol_empty() { + let v: serde_yaml::Value = serde_yaml::from_str("model: b1_afi\n").unwrap(); + assert!(build(&v, &Protocol::default()).is_err()); + } +} diff --git a/crates/qmrust-core/src/models/b1_dam/config.rs b/crates/qmrust-core/src/models/b1_dam/config.rs new file mode 100644 index 0000000..06fd0e8 --- /dev/null +++ b/crates/qmrust-core/src/models/b1_dam/config.rs @@ -0,0 +1,131 @@ +//! b1_dam config, parsed from the top-level YAML keys. +//! +//! The double-angle method has no fit options: its estimate is closed form. +//! The only acquisition parameter is the pair of nominal flip angles, in +//! degrees, in acquisition order — `[alpha, 2*alpha]`. + +use anyhow::{bail, Result}; +use serde::{Deserialize, Serialize}; + +/// Relative tolerance on the `2*alpha` identity between the two flip angles. +/// Loose enough for a sidecar that rounds (e.g. 59.9998 / 120.0), tight enough +/// to reject a series that is not a double-angle pair at all. +const DOUBLE_ANGLE_RTOL: f64 = 1e-3; + +#[derive(Debug, Clone, Default, Deserialize, Serialize)] +pub struct B1DamConfig { + /// Nominal flip angles in degrees, in acquisition order: `[alpha, 2*alpha]`. + #[serde(default)] + pub flip_angles: Vec, +} + +impl B1DamConfig { + /// Config-intrinsic validation. + /// + /// The model has no fit options, but the shape of the acquisition array is + /// intrinsic: `describe` runs only this, and a half-stated protocol would + /// otherwise describe a one-volume collection this model can never fit. + /// Empty is legitimate (nothing has been resolved yet); anything other than + /// the pair the method is defined on is not. + pub fn validate_options(&self) -> Result<()> { + if !matches!(self.flip_angles.len(), 0 | 2) { + bail!( + "b1_dam is defined on exactly 2 flip angles [alpha, 2*alpha], got {}", + self.flip_angles.len() + ); + } + // Stated positively so a non-finite value is rejected too: `fa <= 0.0` + // is false for NaN, which would otherwise fit an all-NaN map, or be + // written into a sidecar as `null` by a `bidsify` that never reaches + // `validate_protocol`. + if !self + .flip_angles + .iter() + .all(|&fa| fa > 0.0 && fa.is_finite()) + { + bail!( + "b1_dam flip angles must be finite and > 0 degrees, got {:?}", + self.flip_angles + ); + } + Ok(()) + } + + /// Protocol-completeness validation: run once the flip angles are final + /// (from `--config` for non-BIDS, or composed from sidecars for BIDS). + /// The double-angle identity `S(2a)/(2*S(a)) = cos(a)` only holds when the + /// second angle really is twice the first, so a series that violates it + /// would yield a plausible-looking but wrong map. + pub fn validate_protocol(&mut self) -> Result<()> { + if self.flip_angles.len() != 2 { + bail!( + "b1_dam needs exactly 2 flip angles [alpha, 2*alpha], got {}", + self.flip_angles.len() + ); + } + // Both values are already known finite and positive (`validate_options`); + // what a resolved protocol adds is the relation between them. + let (alpha, alpha2) = (self.flip_angles[0], self.flip_angles[1]); + if (alpha2 - 2.0 * alpha).abs() > DOUBLE_ANGLE_RTOL * 2.0 * alpha { + bail!( + "b1_dam requires the second flip angle to be twice the first, \ + got alpha={alpha} and {alpha2} degrees" + ); + } + Ok(()) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn parses_top_level_keys() { + let v: serde_yaml::Value = + serde_yaml::from_str("model: b1_dam\nflip_angles: [60, 120]\n").unwrap(); + let mut cfg: B1DamConfig = serde_yaml::from_value(v).unwrap(); + cfg.validate_options().unwrap(); + cfg.validate_protocol().unwrap(); + assert_eq!(cfg.flip_angles, vec![60.0, 120.0]); + } + + #[test] + fn validate_options_passes_without_an_acquisition() { + B1DamConfig::default().validate_options().unwrap(); + } + + #[test] + fn validate_options_rejects_values_no_protocol_could_make_sense_of() { + // Properties of the values themselves, so they are checked by the gate + // `describe` runs: `bidsify` never reaches `validate_protocol`, and + // would otherwise write a NaN into a sidecar as `null`, or describe a + // one-volume collection this model can never fit. + let cfg = |fas: Vec| B1DamConfig { flip_angles: fas }; + assert!(cfg(vec![60.0]).validate_options().is_err()); + assert!(cfg(vec![0.0, 0.0]).validate_options().is_err()); + assert!(cfg(vec![f64::NAN, 120.0]).validate_options().is_err()); + assert!(cfg(vec![60.0, f64::INFINITY]).validate_options().is_err()); + // Nothing resolved yet is not an error: that is what `describe` sees. + assert!(cfg(vec![]).validate_options().is_ok()); + } + + #[test] + fn validate_protocol_requires_a_double_angle_pair() { + let mut wrong_count = B1DamConfig { + flip_angles: vec![60.0], + }; + assert!(wrong_count.validate_protocol().is_err()); + + let mut not_doubled = B1DamConfig { + flip_angles: vec![60.0, 90.0], + }; + assert!(not_doubled.validate_protocol().is_err()); + + // A sidecar that rounds still passes. + let mut rounded = B1DamConfig { + flip_angles: vec![59.9998, 120.0], + }; + rounded.validate_protocol().unwrap(); + } +} diff --git a/crates/qmrust-core/src/models/b1_dam/fit.rs b/crates/qmrust-core/src/models/b1_dam/fit.rs new file mode 100644 index 0000000..5104d0e --- /dev/null +++ b/crates/qmrust-core/src/models/b1_dam/fit.rs @@ -0,0 +1,179 @@ +//! Double-Angle Method (DAM) B1+ mapping — Insko & Bolinger (1993). +//! +//! Two spoiled gradient-echo volumes are acquired at nominal flip angles +//! `alpha` and `2*alpha`. With a common (unknown) amplitude the signals are +//! proportional to `sin(theta)` and `sin(2*theta)`, where `theta = B1 * alpha` +//! is the achieved flip angle, so the amplitude cancels in +//! +//! ```text +//! S(2a) / (2 * S(a)) = sin(2*theta) / (2 * sin(theta)) = cos(theta) +//! B1 = |acos(S(2a) / (2 * S(a)))| / alpha_radians +//! ``` +//! +//! `B1` is dimensionless: the achieved flip angle as a fraction of the nominal +//! one, so `1.0` means the transmit field is exactly on target. This is the +//! same convention the `B1map` aux input of other models expects. +//! +//! The amplitude the ratio divides out is then recovered from either signal as +//! `A = S(a) / sin(theta)`. It is receive-weighted and uncalibrated, so it is a +//! diagnostic rather than a written map — but two measurements determine two +//! unknowns exactly, and without it `forward` could only ever emit the shape of +//! the signal at unit scale, never the signal itself. + +use crate::models::b1_dam::config::B1DamConfig; + +/// Physically plausible transmit-field range, dimensionless. The estimate is +/// closed form and unconstrained; these bounds are descriptive only. +pub const B1_BOUNDS: (f64, f64) = (0.0, 2.0); + +/// Pre-computed fitter for double-angle B1+ data. +pub struct B1DamFitter { + /// Nominal flip angles in degrees, `[alpha, 2*alpha]`. + flip_angles: Vec, +} + +impl B1DamFitter { + pub fn new(cfg: &B1DamConfig) -> Self { + Self { + flip_angles: cfg.flip_angles.clone(), + } + } + + pub fn param_names() -> [&'static str; 2] { + ["B1", "A"] + } + + pub fn output_names() -> [&'static str; 2] { + ["B1", "A"] + } + + /// Nominal flip angles in degrees, in the order `forward`/`fit_voxel` + /// expect them: `[alpha, 2*alpha]`. + pub fn flip_angles(&self) -> &[f64] { + &self.flip_angles + } + + /// The nominal angle the estimate is expressed relative to, in degrees. + pub fn alpha(&self) -> f64 { + self.flip_angles[0] + } + + /// Noise-free DAM signal: `S(theta) = A * sin(B1 * theta)` for each nominal + /// angle, with `A` the receive-weighted equilibrium amplitude. + pub fn forward(&self, b1: f64, a: f64) -> Vec { + self.flip_angles + .iter() + .map(|&fa| a * (b1 * fa.to_radians()).sin()) + .collect() + } + + /// Fit a single voxel from `[S(alpha), S(2*alpha)]`. Returns values in + /// `output_names()` order. + /// + /// Two measurements, two unknowns: the ratio fixes the achieved angle with + /// the amplitude divided out, and the amplitude then follows from either + /// signal. Recovering it is what lets `forward` reproduce the measured + /// data rather than a unit-scaled shape of it. + pub fn fit_voxel(&self, signal: &[f64]) -> Vec { + let theta = acos_abs(signal[1] / (2.0 * signal[0])); + vec![theta / self.alpha().to_radians(), signal[0] / theta.sin()] + } +} + +/// `abs(acos(r))` evaluated over the complex plane, as MATLAB's `acos` does for +/// a real argument outside `[-1, 1]`. Noise and imperfect spoiling routinely +/// push the ratio out of domain, and qMRLab reports the magnitude of the +/// complex principal value there rather than a NaN, so voxels stay finite: +/// +/// ```text +/// r > 1: acos(r) = -i*acosh(r) -> |.| = acosh(r) +/// r < -1: acos(r) = pi - i*acosh(-r) -> |.| = hypot(pi, acosh(-r)) +/// ``` +/// +/// A non-finite ratio (a zero or NaN denominator) stays NaN, and the engine +/// records the voxel as unfitted. +fn acos_abs(r: f64) -> f64 { + if r.is_nan() { + f64::NAN + } else if r > 1.0 { + r.acosh() + } else if r < -1.0 { + std::f64::consts::PI.hypot((-r).acosh()) + } else { + r.acos().abs() + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn fitter(alpha: f64) -> B1DamFitter { + B1DamFitter::new(&B1DamConfig { + flip_angles: vec![alpha, 2.0 * alpha], + }) + } + + #[test] + fn forward_then_fit_recovers_b1_and_amplitude() { + let f = fitter(60.0); + for &truth in &[0.6, 0.85, 1.0, 1.15, 1.4] { + for & in &[1.0, 13_500.0] { + let sig = f.forward(truth, amp); + let out = f.fit_voxel(&sig); + assert!((out[0] - truth).abs() < 1e-12, "B1={truth}: got {}", out[0]); + assert!( + (out[1] - amp).abs() / amp < 1e-12, + "A={amp}: got {}", + out[1] + ); + } + } + } + + #[test] + fn the_fitted_pair_reproduces_the_measured_signal() { + // What the playground's forward curve plots against the measured + // points: with the amplitude recovered, the prediction lands on the + // data rather than on a unit-scaled shape of it near zero. + let f = fitter(60.0); + let measured = [13_500.0, 15_600.0]; + let out = f.fit_voxel(&measured); + let predicted = f.forward(out[0], out[1]); + for (m, p) in measured.iter().zip(&predicted) { + assert!((m - p).abs() / m < 1e-12, "measured {m} vs predicted {p}"); + } + } + + #[test] + fn a_common_amplitude_cancels() { + // Only the ratio of the two signals matters, so receive sensitivity and + // proton density drop out. + let f = fitter(60.0); + let plain = f.forward(0.9, 1.0); + let scaled: Vec = plain.iter().map(|s| s * 4321.0).collect(); + assert!((f.fit_voxel(&plain)[0] - f.fit_voxel(&scaled)[0]).abs() < 1e-12); + // The amplitude, by contrast, tracks the scale exactly. + assert!((f.fit_voxel(&scaled)[1] / f.fit_voxel(&plain)[1] - 4321.0).abs() < 1e-6); + } + + #[test] + fn out_of_domain_ratios_follow_matlabs_complex_acos_magnitude() { + // acos(2) = 1.3170i; acos(-2) = pi - 1.3170i. + assert!((acos_abs(2.0) - 2.0_f64.acosh()).abs() < 1e-12); + let expected = (std::f64::consts::PI.powi(2) + 2.0_f64.acosh().powi(2)).sqrt(); + assert!((acos_abs(-2.0) - expected).abs() < 1e-12); + // In domain it is the ordinary principal value. + assert!((acos_abs(0.5) - 0.5_f64.acos()).abs() < 1e-12); + assert!(acos_abs(f64::NAN).is_nan()); + } + + #[test] + fn a_zero_reference_signal_yields_nan() { + // 0/0 is NaN; a nonzero numerator over zero is an infinite ratio, whose + // acosh is also infinite — neither is a usable estimate. + let f = fitter(60.0); + assert!(f.fit_voxel(&[0.0, 0.0])[0].is_nan()); + assert!(!f.fit_voxel(&[0.0, 1.0])[0].is_finite()); + } +} diff --git a/crates/qmrust-core/src/models/b1_dam/mod.rs b/crates/qmrust-core/src/models/b1_dam/mod.rs new file mode 100644 index 0000000..ca99c58 --- /dev/null +++ b/crates/qmrust-core/src/models/b1_dam/mod.rs @@ -0,0 +1,7 @@ +//! Double-Angle Method B1+ mapping. BIDS suffix: `TB1DAM`. + +pub mod config; +pub mod fit; +pub mod model; + +pub use model::{build, describe, dump, effective}; diff --git a/crates/qmrust-core/src/models/b1_dam/model.rs b/crates/qmrust-core/src/models/b1_dam/model.rs new file mode 100644 index 0000000..056cd31 --- /dev/null +++ b/crates/qmrust-core/src/models/b1_dam/model.rs @@ -0,0 +1,337 @@ +//! b1_dam adapter onto the core `Model` trait. +//! +//! A `Series` of exactly two spoiled gradient-echo volumes, indexed by flip +//! angle, combined by the closed-form double-angle ratio — there is no +//! acquisition parameter beyond the angles themselves and no iterative fit. + +use crate::core::model::{ + Aux, BidsSpec, BidsVolume, EntityRole, FitStrategy, InputSpec, Measurement, MeasurementKind, + Model, ProtoParam, Protocol, Sample, Scope, Source, +}; +use crate::models::b1_dam::config::B1DamConfig; +use crate::models::b1_dam::fit::{B1DamFitter, B1_BOUNDS}; +use anyhow::Result; +use serde_json::json; +use std::collections::BTreeMap; + +pub struct B1DamModel { + fitter: B1DamFitter, +} + +const B1DAM_ENTITIES: &[EntityRole] = &[EntityRole::Flip]; + +/// One `{"FlipAngle": deg}` identity row per volume, in acquisition order. +fn b1_dam_rows(fitter: &B1DamFitter) -> Vec> { + fitter + .flip_angles() + .iter() + .map(|&fa| BTreeMap::from([("FlipAngle".to_string(), fa)])) + .collect() +} + +impl B1DamModel { + pub fn new(cfg: B1DamConfig) -> Self { + Self { + fitter: B1DamFitter::new(&cfg), + } + } +} + +impl Model for B1DamModel { + fn param_names(&self) -> Vec<&'static str> { + B1DamFitter::param_names().to_vec() + } + fn output_names(&self) -> Vec { + B1DamFitter::output_names() + .iter() + .map(|s| s.to_string()) + .collect() + } + fn param_bounds(&self) -> Vec<(f64, f64)> { + // The amplitude is an uncalibrated signal level, so it is unbounded. + vec![B1_BOUNDS, (f64::NEG_INFINITY, f64::INFINITY)] + } + fn fixed_mask(&self) -> Vec { + vec![false; 2] + } + fn required_inputs(&self) -> Vec { + // B1+ is what this model measures; it consumes no auxiliary map. + vec![] + } + fn measurement(&self) -> MeasurementKind { + MeasurementKind::Series { + rows: b1_dam_rows(&self.fitter), + } + } + fn strategy(&self) -> FitStrategy { + FitStrategy::Voxelwise + } + fn forward(&self, params: &[f64], _aux: &Aux) -> Measurement { + let values = self.fitter.forward(params[0], params[1]); + let samples = self + .fitter + .flip_angles() + .iter() + .zip(values) + .map(|(&fa, value)| Sample { + params: BTreeMap::from([("FlipAngle".to_string(), fa)]), + value, + }) + .collect(); + Measurement::Series(samples) + } + fn fit(&self, m: &Measurement, _aux: &Aux) -> Vec { + // Assemble in the fitter's own angle order by matching each expected + // angle to its sample by identity — never positionally. Swapping the + // two volumes would otherwise invert the ratio silently. + let samples = m.series(); + let signal: Vec = self + .fitter + .flip_angles() + .iter() + .map(|&fa| { + samples + .iter() + .find(|s| s.params.get("FlipAngle") == Some(&fa)) + .map(|s| s.value) + .unwrap_or_else(|| panic!("measurement has no sample with FlipAngle={fa}")) + }) + .collect(); + self.fitter.fit_voxel(&signal) + } + fn n_volumes(&self) -> usize { + self.fitter.flip_angles().len() + } + fn bids_volume(&self, index: usize) -> BidsVolume { + BidsVolume { + entities: vec![("flip", (index + 1).to_string())], + sidecar: BTreeMap::from([( + "FlipAngle".to_string(), + json!(self.fitter.flip_angles()[index]), + )]), + } + } + fn bids(&self) -> Option { + Some(BidsSpec { + suffix: "TB1DAM", + entities: B1DAM_ENTITIES, + }) + } + fn protocol_schema(&self) -> Vec { + // Flip angle is the whole acquisition axis, so it alone identifies a + // volume: a `Series`' per-volume protocol rows *are* its volume + // identities (`engine::build_volume_ids`), and `forward` tags its + // samples the same way. + vec![ProtoParam { + name: "FlipAngle", + source: Source::Field("FlipAngle"), + scope: Scope::PerVolume, + required: true, + }] + } + fn bids_outputs(&self) -> Vec<(&'static str, &'static str, &'static str)> { + // Dimensionless: the achieved flip angle as a fraction of the nominal + // one, which is the scaling every `B1map` aux consumer expects. `A` is + // deliberately absent: it is a receive-weighted signal level with no + // BIDS suffix of its own, so it stays a diagnostic output rather than + // becoming a written map. + vec![("B1", "TB1map", "")] + } +} + +impl crate::core::model::ModelConfig for B1DamConfig { + const NAME: &'static str = "b1_dam"; + const SUBKEY: Option<&'static str> = None; + const PROTOCOL_KEYS: &'static [&'static str] = &["flip_angles"]; + + fn validate_options(&mut self) -> Result<()> { + B1DamConfig::validate_options(self) + } + + fn ingest_protocol(&mut self, proto: &Protocol) -> Result<()> { + let angles: Vec = proto + .volumes + .iter() + .filter_map(|m| m.get("FlipAngle").copied()) + .collect(); + if !angles.is_empty() { + self.flip_angles = angles; + } + Ok(()) + } + + fn validate_protocol(&mut self) -> Result<()> { + B1DamConfig::validate_protocol(self) + } + + fn into_model(self) -> Box { + Box::new(B1DamModel::new(self)) + } +} + +/// Structural interrogation entry point (see [`describe_model`](crate::core::model::describe_model)). +pub fn describe(v: &serde_yaml::Value) -> Result> { + crate::core::model::describe_model::(v) +} + +/// Registry builder (see [`build_model`](crate::core::model::build_model)). +pub fn build(v: &serde_yaml::Value, proto: &Protocol) -> Result> { + crate::core::model::build_model::(v, proto) +} + +/// Registry dumper (see [`dump_model`](crate::core::model::dump_model)). +pub fn dump(v: &serde_yaml::Value) -> Result { + crate::core::model::dump_model::(v) +} + +/// Registry option-surface entry point (see +/// [`effective_model`](crate::core::model::effective_model)): every option this +/// model accepts, at its effective value, plus any validation complaint. +pub fn effective( + v: &serde_yaml::Value, + proto: &Protocol, +) -> Result { + crate::core::model::effective_model::(v, proto) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn b1_dam_value() -> serde_yaml::Value { + serde_yaml::from_str("model: b1_dam\nflip_angles: [60, 120]\n").unwrap() + } + + /// A resolved protocol as the shell hands it over: one row per volume + /// carrying the identifying axis. + fn resolved_proto(angles: &[f64]) -> Protocol { + Protocol { + volumes: angles + .iter() + .map(|&fa| BTreeMap::from([("FlipAngle".to_string(), fa)])) + .collect(), + global: BTreeMap::new(), + } + } + + #[test] + fn build_and_roundtrip_via_trait() { + let m = build(&b1_dam_value(), &Protocol::default()).unwrap(); + assert_eq!(m.param_names(), vec!["B1", "A"]); + assert_eq!(m.output_names(), vec!["B1".to_string(), "A".to_string()]); + assert_eq!(m.n_volumes(), 2); + let sig = m.forward(&[0.85, 13_500.0], &Aux::new()); + let fitted = m.fit(&sig, &Aux::new()); + assert!((fitted[0] - 0.85).abs() < 1e-12, "B1: {}", fitted[0]); + assert!((fitted[1] - 13_500.0).abs() < 1e-6, "A: {}", fitted[1]); + } + + #[test] + fn fit_assembles_by_identity_not_position() { + // Reversing the samples must not invert the ratio. + let m = build(&b1_dam_value(), &Protocol::default()).unwrap(); + let sig = m.forward(&[0.85, 1.0], &Aux::new()); + let mut reversed: Vec = sig + .series() + .iter() + .map(|s| Sample { + params: s.params.clone(), + value: s.value, + }) + .collect(); + reversed.reverse(); + assert_eq!( + m.fit(&sig, &Aux::new()), + m.fit(&Measurement::Series(reversed), &Aux::new()) + ); + } + + #[test] + #[should_panic(expected = "no sample with FlipAngle")] + fn fit_panics_on_unmatched_identity() { + let m = build(&b1_dam_value(), &Protocol::default()).unwrap(); + let bogus = Measurement::Series(vec![Sample { + params: BTreeMap::from([("FlipAngle".to_string(), 77.0)]), + value: 1.0, + }]); + let _ = m.fit(&bogus, &Aux::new()); + } + + #[test] + fn bids_folds_flip_angles_from_protocol() { + let proto = resolved_proto(&[45.0, 90.0]); + // Config carries no acquisition; the sidecars supply it. + let v: serde_yaml::Value = serde_yaml::from_str("model: b1_dam\n").unwrap(); + let m = build(&v, &proto).unwrap(); + assert_eq!(m.n_volumes(), 2); + let second = m.bids_volume(1); + assert_eq!(second.entities, vec![("flip", "2".to_string())]); + assert_eq!(second.sidecar["FlipAngle"], json!(90.0)); + } + + #[test] + fn forward_samples_carry_the_volume_identities_the_bids_path_builds() { + // A `Series` model's volume identities come from the resolved + // per-volume protocol (`engine::build_volume_ids`), while `forward` + // tags its samples from the model's own rows. Any drift between the two + // leaves the fit working (it queries one key by name) while the app's + // forward curve silently disappears, so assert they agree exactly. + let proto = resolved_proto(&[60.0, 120.0]); + let v: serde_yaml::Value = serde_yaml::from_str("model: b1_dam\n").unwrap(); + let m = build(&v, &proto).unwrap(); + + let ids = crate::engine::build_volume_ids(m.measurement(), &proto, m.n_volumes()).unwrap(); + let sig = m.forward(&[0.9, 1.0], &Aux::new()); + let samples = sig.series(); + assert_eq!(samples.len(), ids.len()); + for (id, sample) in ids.iter().zip(samples) { + let crate::core::model::VolumeId::Params(row) = id else { + panic!("b1_dam is a Series model; expected param-row identities") + }; + assert_eq!( + *row, sample.params, + "volume identity {row:?} has no matching forward sample identity {:?}", + sample.params + ); + } + } + + #[test] + fn declares_bids_tb1dam_and_no_aux() { + let m = build(&b1_dam_value(), &Protocol::default()).unwrap(); + assert_eq!(m.bids().unwrap().suffix, "TB1DAM"); + assert!(m.required_inputs().is_empty()); + } + + #[test] + fn bids_outputs_reference_real_output_names() { + let m = build(&b1_dam_value(), &Protocol::default()).unwrap(); + let names = m.output_names(); + for (out, _suffix, _unit) in m.bids_outputs() { + assert!(names.iter().any(|n| n == out), "{out} not in {names:?}"); + } + } + + #[test] + fn describe_succeeds_without_an_acquisition_and_exposes_schema() { + let v: serde_yaml::Value = serde_yaml::from_str("model: b1_dam\n").unwrap(); + let m = describe(&v).unwrap(); + let schema = m.protocol_schema(); + assert_eq!(schema.len(), 1); + assert_eq!(schema[0].name, "FlipAngle"); + assert!(matches!(schema[0].scope, Scope::PerVolume)); + } + + #[test] + fn build_rejects_a_series_that_is_not_a_double_angle_pair() { + let proto = resolved_proto(&[60.0, 90.0]); + let v: serde_yaml::Value = serde_yaml::from_str("model: b1_dam\n").unwrap(); + assert!(build(&v, &proto).is_err()); + } + + #[test] + fn build_still_requires_an_acquisition_when_protocol_empty() { + let v: serde_yaml::Value = serde_yaml::from_str("model: b1_dam\n").unwrap(); + assert!(build(&v, &Protocol::default()).is_err()); + } +} diff --git a/crates/qmrust-core/src/models/mod.rs b/crates/qmrust-core/src/models/mod.rs index e52b6ca..b1ce559 100644 --- a/crates/qmrust-core/src/models/mod.rs +++ b/crates/qmrust-core/src/models/mod.rs @@ -1,3 +1,5 @@ +pub mod b1_afi; +pub mod b1_dam; pub mod inversion_recovery; pub mod mono_t2; pub mod mt_ratio; diff --git a/crates/qmrust-core/src/registry.rs b/crates/qmrust-core/src/registry.rs index d7b46b3..23c2476 100644 --- a/crates/qmrust-core/src/registry.rs +++ b/crates/qmrust-core/src/registry.rs @@ -14,26 +14,89 @@ pub type Effective = fn(&serde_yaml::Value, &Protocol) -> Result &'static str { match self { - Category::Relaxometry => "relaxometry", - Category::MagnetizationTransfer => "magnetization-transfer", + Category::T1Relaxometry => "t1-relaxometry", + Category::T2Relaxometry => "t2-relaxometry", + Category::FieldMapping => "field-mapping", + Category::SemiQuantitativeMt | Category::QuantitativeMt => "magnetization-transfer", } } - /// Human-readable heading for the model gallery. + /// Top-level heading a reader browses by. + pub fn family(&self) -> &'static str { + match self { + Category::T1Relaxometry => "T1 Relaxometry", + Category::T2Relaxometry => "T2 Relaxometry", + Category::FieldMapping => "Field Mapping", + Category::SemiQuantitativeMt | Category::QuantitativeMt => "Magnetization Transfer", + } + } + + /// Heading within the family, for families worth subdividing. `None` when + /// the family's models sit directly under it. + pub fn subgroup(&self) -> Option<&'static str> { + match self { + Category::SemiQuantitativeMt => Some("Semi-quantitative MT"), + Category::QuantitativeMt => Some("Quantitative MT"), + _ => None, + } + } + + /// The most specific heading this category carries — its subgroup when it + /// has one, its family otherwise. pub fn title(&self) -> &'static str { + self.subgroup().unwrap_or_else(|| self.family()) + } + + /// Lucide icon name for this category's family, for surfaces that show one. + /// + /// Named here for the same reason the title and the order are: a family is + /// a registry fact, and a picker that mapped headings to glyphs itself + /// would need editing every time the taxonomy moved. Categories sharing a + /// family share its icon. + pub fn icon(&self) -> &'static str { match self { - Category::Relaxometry => "Relaxometry", - Category::MagnetizationTransfer => "Magnetization transfer", + Category::T1Relaxometry => "spline-pointer", + Category::T2Relaxometry => "spline-pointer-mirrored", + Category::FieldMapping => "circle-dot", + Category::SemiQuantitativeMt | Category::QuantitativeMt => "waves-arrow-up", + } + } + + /// Where this category sits in the order a reader is offered them. + /// + /// The registry lists models alphabetically, which is the wrong order to + /// read a method taxonomy in, and every surface that groups models — the + /// documentation gallery, the playground's picker — wants the same one. + /// Sorting by this and then grouping consecutive runs yields the whole + /// two-level tree, so the order lives here rather than being restated + /// wherever a list is drawn. + pub fn order(&self) -> u8 { + match self { + Category::T1Relaxometry => 0, + Category::T2Relaxometry => 1, + Category::FieldMapping => 2, + Category::SemiQuantitativeMt => 3, + Category::QuantitativeMt => 4, } } } @@ -102,6 +165,82 @@ pub struct ModelEntry { pub fn all() -> &'static [ModelEntry] { &[ + ModelEntry { + name: "b1_afi", + bids_suffix: "TB1AFI", + build: models::b1_afi::build, + describe: models::b1_afi::describe, + dump: models::b1_afi::dump, + effective: models::b1_afi::effective, + doc: ModelDoc { + title: "Actual Flip Angle B1+", + category: Category::FieldMapping, + summary: "Maps the transmit (B1+) field from a single spoiled \ + gradient-echo sequence that interleaves two excitation \ + repetition times at one nominal flip angle. In the pulsed \ + steady state the ratio of the two signals depends on the \ + achieved flip angle alone, so proton density and receive \ + sensitivity cancel. The result is dimensionless — the \ + achieved flip angle as a fraction of the nominal one — and \ + feeds the B1 correction of models such as MTsat and VFA. \ + The closed form assumes the repetition times are short \ + against T1, and underestimates B1 as that assumption \ + weakens.", + equation: r"\cos\theta = \frac{r n - 1}{n - r},\quad r = \left|\frac{S_{\mathrm{TR}_2}}{S_{\mathrm{TR}_1}}\right|,\quad n = \frac{\mathrm{TR}_2}{\mathrm{TR}_1},\quad B_1^{+} = \frac{\theta}{\theta_\mathrm{nom}}", + symbols: &[ + ("B1", "Relative transmit field", ""), + ("T1", "Longitudinal relaxation time", "s"), + ], + citations: &["yarnykh2007"], + source_dir: "crates/qmrust-core/src/models/b1_afi", + recipes: Recipes { + bids: "recipes/bids/b1_afi_config.yaml", + non_bids: "recipes/non-bids/b1_afi_config.yaml", + sim: None, + }, + enums: &[], + display_ranges: &[ + // A well-shimmed transmit coil sits within a few tens of + // percent of nominal; this window keeps 1.0 mid-scale. + ("B1", 0.5, 1.5), + ], + }, + }, + ModelEntry { + name: "b1_dam", + bids_suffix: "TB1DAM", + build: models::b1_dam::build, + describe: models::b1_dam::describe, + dump: models::b1_dam::dump, + effective: models::b1_dam::effective, + doc: ModelDoc { + title: "Double Angle B1+", + category: Category::FieldMapping, + summary: "Maps the transmit (B1+) field from two spoiled \ + gradient-echo volumes acquired at flip angles alpha and \ + twice alpha. The ratio of the two signals depends only on \ + the achieved flip angle, so proton density and receive \ + sensitivity cancel and the map needs no calibration. The \ + result is dimensionless — the achieved flip angle as a \ + fraction of the nominal one — and feeds the B1 correction \ + of models such as MTsat and VFA.", + equation: r"B_1^{+} = \frac{\left|\arccos\left(\frac{S_{2\alpha}}{2 S_{\alpha}}\right)\right|}{\alpha}", + symbols: &[("B1", "Relative transmit field", "")], + citations: &["insko1993"], + source_dir: "crates/qmrust-core/src/models/b1_dam", + recipes: Recipes { + bids: "recipes/bids/b1_dam_config.yaml", + non_bids: "recipes/non-bids/b1_dam_config.yaml", + sim: None, + }, + enums: &[], + display_ranges: &[ + // A well-shimmed transmit coil sits within a few tens of + // percent of nominal; this window keeps 1.0 mid-scale. + ("B1", 0.5, 1.5), + ], + }, + }, ModelEntry { name: "mt_ratio", bids_suffix: "MTR", @@ -110,8 +249,8 @@ pub fn all() -> &'static [ModelEntry] { dump: models::mt_ratio::dump, effective: models::mt_ratio::effective, doc: ModelDoc { - title: "Magnetization transfer ratio", - category: Category::MagnetizationTransfer, + title: "MT Ratio", + category: Category::SemiQuantitativeMt, summary: "Computes the magnetization transfer ratio from two \ images: one acquired with an off-resonance saturation pulse \ and one without. MTR is a semi-quantitative percentage that \ @@ -143,8 +282,8 @@ pub fn all() -> &'static [ModelEntry] { dump: models::mt_sat::dump, effective: models::mt_sat::effective, doc: ModelDoc { - title: "MT saturation", - category: Category::MagnetizationTransfer, + title: "MT Saturation", + category: Category::SemiQuantitativeMt, summary: "Derives the MT saturation parameter from three spoiled \ gradient-echo volumes — MT-weighted, PD-weighted and \ T1-weighted. Unlike MTR, MTsat removes the leading-order \ @@ -183,8 +322,8 @@ pub fn all() -> &'static [ModelEntry] { dump: models::mono_t2::dump, effective: models::mono_t2::effective, doc: ModelDoc { - title: "Mono-exponential T2", - category: Category::Relaxometry, + title: "Monoexp T2", + category: Category::T2Relaxometry, summary: "Fits the transverse relaxation time T2 from a \ multi-echo spin-echo series as a mono-exponential decay with \ a free amplitude. The fit runs either as a log-linear \ @@ -222,8 +361,8 @@ pub fn all() -> &'static [ModelEntry] { dump: models::inversion_recovery::dump, effective: models::inversion_recovery::effective, doc: ModelDoc { - title: "Inversion recovery T1", - category: Category::Relaxometry, + title: "Inversion Recovery", + category: Category::T1Relaxometry, summary: "Fits the longitudinal relaxation time T1 from a series \ of inversion-recovery images acquired at different inversion \ times. The magnitude signal is modelled as an exponential \ @@ -257,8 +396,8 @@ pub fn all() -> &'static [ModelEntry] { dump: models::vfa_t1::dump, effective: models::vfa_t1::effective, doc: ModelDoc { - title: "Variable flip angle T1", - category: Category::Relaxometry, + title: "Variable Flip Angle", + category: Category::T1Relaxometry, summary: "Fits the longitudinal relaxation time T1 from spoiled \ gradient-echo images acquired at two or more excitation flip \ angles and a single repetition time. Dividing the \ @@ -296,7 +435,7 @@ pub fn all() -> &'static [ModelEntry] { effective: models::qmt_spgr::effective, doc: ModelDoc { title: "qMT-SPGR", - category: Category::MagnetizationTransfer, + category: Category::QuantitativeMt, summary: "Two-pool quantitative magnetization transfer from a \ spoiled gradient-echo sequence with off-resonance saturation \ sampled across a grid of saturation flip angles and \ @@ -416,12 +555,53 @@ mod tests { for e in all() { let slug = e.doc.category.slug(); assert!( - slug.chars().all(|c| c.is_ascii_lowercase() || c == '-'), + slug.chars() + .all(|c| c.is_ascii_lowercase() || c.is_ascii_digit() || c == '-'), "{}: category slug '{}' is not URL-safe", e.name, slug ); assert!(!e.doc.category.title().is_empty()); + assert!(!e.doc.category.family().is_empty()); + } + } + + /// Several categories may share a documentation directory — the two + /// magnetization-transfer ones do — but only within a single family. Two + /// families writing into one directory would interleave their pages under + /// a heading naming just one of them, and the gallery links (built from the + /// slug) would point into the wrong section. + #[test] + fn categories_sharing_a_slug_belong_to_one_family() { + let mut family_of: std::collections::BTreeMap<&str, &str> = Default::default(); + for e in all() { + let (slug, family) = (e.doc.category.slug(), e.doc.category.family()); + if let Some(seen) = family_of.insert(slug, family) { + assert_eq!( + seen, family, + "slug '{slug}' is claimed by both '{seen}' and '{family}'" + ); + } + } + } + + /// Reading order must be a total order over the categories in play: two + /// categories at the same rank would group and sort unpredictably in both + /// the gallery and the picker, which sort by exactly this. + #[test] + fn category_order_is_unambiguous() { + let mut by_order: std::collections::BTreeMap = Default::default(); + for e in all() { + let c = e.doc.category; + if let Some(seen) = by_order.insert(c.order(), c.title()) { + assert_eq!( + seen, + c.title(), + "order {} is shared by '{seen}' and '{}'", + c.order(), + c.title() + ); + } } } } diff --git a/crates/qmrust-wasm/src/bids.rs b/crates/qmrust-wasm/src/bids.rs index f39a79d..f2d0773 100644 --- a/crates/qmrust-wasm/src/bids.rs +++ b/crates/qmrust-wasm/src/bids.rs @@ -183,7 +183,15 @@ pub fn resolve_bids( .collect(), volume_ids_json: emit_volume_ids(&volume_ids)?, protocol_json: emit_protocol(&proto)?, - warnings: c.warnings.iter().map(|w| w.message.clone()).collect(), + // Grouping warnings and input-resolution warnings reach the reader + // through the same channel: both say the dataset did not hold + // something the recipe asked for. + warnings: c + .warnings + .iter() + .map(|w| w.message.clone()) + .chain(paths.warnings.iter().cloned()) + .collect(), files, }); } diff --git a/crates/rust-bids/src/default_grouping.yaml b/crates/rust-bids/src/default_grouping.yaml index 8209100..fe1db1a 100644 --- a/crates/rust-bids/src/default_grouping.yaml +++ b/crates/rust-bids/src/default_grouping.yaml @@ -33,3 +33,9 @@ MTR: VFA: sequential_set: by: [flip] +TB1DAM: + sequential_set: + by: [flip] +TB1AFI: + sequential_set: + by: [acq] diff --git a/crates/rust-bids/src/inputs.rs b/crates/rust-bids/src/inputs.rs index 808dda5..f04e3e3 100644 --- a/crates/rust-bids/src/inputs.rs +++ b/crates/rust-bids/src/inputs.rs @@ -31,6 +31,14 @@ impl MaskSpec { /// Parse a recipe's `mask:` value, normalizing its entity keys (e.g. `desc` /// → `description`) to the form the file table stores them in. `None` when /// the recipe declares no `mask:` — no masking is then applied. + /// + /// An entity left blank is dropped rather than matched. A recipe is also a + /// template, and a blank field is the honest way to offer a knob for a + /// dataset whose mask nobody has named yet: writing a value there that the + /// example data does not have would state a fact about the dataset that is + /// not true. Blank therefore means "unconstrained" — with every entity + /// blank, `mask:` selects whatever single mask the dataset holds, and + /// nothing when it holds none. pub fn from_recipe(raw: &serde_yaml::Value, vocab: &Vocabulary) -> Result> { let Some(v) = raw.get("mask") else { return Ok(None); @@ -42,6 +50,7 @@ impl MaskSpec { entities: spec .entities .into_iter() + .filter(|(_, v)| !v.trim().is_empty()) .map(|(k, v)| (vocab.normalize_entity_key(&k), v)) .collect(), })) @@ -56,6 +65,10 @@ pub struct InputPaths { /// on its own default. pub aux: Vec<(String, Option)>, pub mask: Option, + /// What the caller should be told about inputs it asked for and did not + /// get. Absence is legitimate here — an optional aux and an unmatched mask + /// both leave the fit runnable — so these are reported, never fatal. + pub warnings: Vec, } impl InputPaths { @@ -143,16 +156,58 @@ pub fn resolve_input_paths( } } + let mut warnings = Vec::new(); let mask = match mask_spec { Some(spec) => { let mut extra: Vec<(&str, &str)> = vec![("suffix", spec.suffix.as_str())]; extra.extend(spec.entities.iter().map(|(k, v)| (k.as_str(), v.as_str()))); - find_row(table, identity, &extra)?.map(|row| row.path.clone()) + let found = find_row(table, identity, &extra)?.map(|row| row.path.clone()); + if found.is_none() { + // Declaring a mask and getting none is not an error — the fit + // runs over the whole image — but it must not pass in silence. + // The recipe is echoed verbatim into the provenance + // `Parameters`, so an unmatched `mask:` block reads there + // exactly like one that resolved, and only the `Sources` list + // (which omits it) reveals that nothing was applied. + warnings.push(match describe_mask(spec) { + Some(what) => format!( + "the recipe selects a mask ({what}) that this dataset has no match \ + for; fitting the whole image unmasked" + ), + None => format!( + "the recipe asks for a '{}' file and this dataset has none for \ + this collection; fitting the whole image unmasked", + spec.suffix + ), + }); + } + found } None => None, }; - Ok(InputPaths { aux, mask }) + Ok(InputPaths { + aux, + mask, + warnings, + }) +} + +/// A `mask:` block as a reader wrote it — `desc-brain mask` — for naming the +/// thing that was asked for in a message about not finding it. `None` when the +/// block names nothing in particular, which has no useful description and reads +/// as an empty parenthesis if forced into one. +fn describe_mask(spec: &MaskSpec) -> Option { + if spec.entities.is_empty() { + return None; + } + let mut parts: Vec = spec + .entities + .iter() + .map(|(k, v)| format!("{}-{}", crate::entities::short_key(k), v)) + .collect(); + parts.push(spec.suffix.clone()); + Some(parts.join(" ")) } #[cfg(test)] @@ -252,4 +307,119 @@ mod tests { // Suffix alone is too loose — both masks match, so it errors. assert!(find_row(&rows, &identity, &[("suffix", "mask")]).is_err()); } + + /// A recipe that selects a mask the dataset does not hold must say so. The + /// fit is still valid — it simply covers the whole image — but the recipe + /// is echoed verbatim into the provenance `Parameters`, where an unmatched + /// `mask:` block is indistinguishable from one that resolved. Only the + /// `Sources` list, by omission, records that nothing was applied, so + /// silence here reads as a masked fit to anyone who does not cross-check. + #[test] + fn a_mask_the_dataset_does_not_hold_is_reported_not_silently_skipped() { + use crate::vocab::Vocabulary; + use std::collections::BTreeMap; + + struct NoAuxModel; + impl Model for NoAuxModel { + fn param_names(&self) -> Vec<&'static str> { + vec!["x"] + } + fn output_names(&self) -> Vec { + vec!["x".into()] + } + fn param_bounds(&self) -> Vec<(f64, f64)> { + vec![(0.0, 1.0)] + } + fn fixed_mask(&self) -> Vec { + vec![false] + } + fn required_inputs(&self) -> Vec { + vec![] + } + fn measurement(&self) -> qmrust_core::core::model::MeasurementKind { + qmrust_core::core::model::MeasurementKind::Series { rows: vec![] } + } + fn forward( + &self, + _p: &[f64], + _a: &qmrust_core::core::model::Aux, + ) -> qmrust_core::core::model::Measurement { + qmrust_core::core::model::Measurement::Series(vec![]) + } + fn fit( + &self, + _m: &qmrust_core::core::model::Measurement, + _a: &qmrust_core::core::model::Aux, + ) -> Vec { + vec![0.0] + } + fn n_volumes(&self) -> usize { + 0 + } + fn bids_volume(&self, _i: usize) -> qmrust_core::core::model::BidsVolume { + unreachable!() + } + } + + let vocab = Vocabulary::bids(); + let raw: serde_yaml::Value = serde_yaml::from_str("mask:\n desc: brain\n").unwrap(); + let spec = MaskSpec::from_recipe(&raw, &vocab).unwrap().unwrap(); + let identity = BTreeMap::from([("subject".to_string(), "01".to_string())]); + + // A dataset holding no mask at all. + let paths = resolve_input_paths(&[], &NoAuxModel, &identity, Some(&spec)).unwrap(); + assert!(paths.mask.is_none()); + assert_eq!(paths.warnings.len(), 1, "{:?}", paths.warnings); + let msg = &paths.warnings[0]; + assert!( + msg.contains("desc-brain mask"), + "warning should name what was asked for: {msg}" + ); + + // Declaring no mask is not a warning: nothing was asked for. + let quiet = resolve_input_paths(&[], &NoAuxModel, &identity, None).unwrap(); + assert!(quiet.warnings.is_empty()); + } + + /// A blank entity is a knob nobody has set, not a demand for the empty + /// string. Recipes for datasets that ship no mask leave `desc` empty rather + /// than naming one the data does not have, so blank has to mean + /// "unconstrained" — otherwise those recipes would match nothing even on a + /// dataset that does have a mask. + #[test] + fn a_blank_mask_entity_is_unconstrained_not_a_match_on_the_empty_string() { + use crate::vocab::Vocabulary; + use std::collections::BTreeMap; + + let vocab = Vocabulary::bids(); + let raw: serde_yaml::Value = serde_yaml::from_str("mask:\n desc: \"\"\n").unwrap(); + let spec = MaskSpec::from_recipe(&raw, &vocab).unwrap().unwrap(); + assert_eq!(spec.suffix, "mask"); + assert!( + spec.entities.is_empty(), + "blank entity was kept: {:?}", + spec.entities + ); + + // With nothing constrained, the dataset's one mask is the one selected. + let row = BidsRow { + path: "derivatives/preprocessed/sub-01/anat/sub-01_desc-brain_mask.nii.gz".into(), + derivatives: Some("preprocessed".into()), + datatype: Some("anat".into()), + suffix: "mask".into(), + extension: ".nii.gz".into(), + entities: BTreeMap::from([ + ("subject".to_string(), "01".to_string()), + ("description".to_string(), "brain".to_string()), + ]), + sidecar_path: None, + }; + let identity = BTreeMap::from([("subject".to_string(), "01".to_string())]); + let rows = [row]; + let hit = find_row(&rows, &identity, &[("suffix", "mask")]).unwrap(); + assert!( + hit.is_some(), + "a blank spec should accept the only mask present" + ); + } } diff --git a/crates/rust-bids/src/lib.rs b/crates/rust-bids/src/lib.rs index 6e8bb83..3f05364 100644 --- a/crates/rust-bids/src/lib.rs +++ b/crates/rust-bids/src/lib.rs @@ -24,4 +24,4 @@ pub use resolve::{collections_for, resolve_set}; pub use scan::scan_dataset; pub use sidecar::{sidecar_for, sidecar_sources_for, Sidecar}; pub use table::{parse_to_table, row_column, table_filter, BidsRow}; -pub use vocab::Vocabulary; +pub use vocab::{datatype_for_suffix, Vocabulary}; diff --git a/crates/rust-bids/src/vocab.rs b/crates/rust-bids/src/vocab.rs index 60c2cd6..e6f9201 100644 --- a/crates/rust-bids/src/vocab.rs +++ b/crates/rust-bids/src/vocab.rs @@ -243,6 +243,33 @@ impl Vocabulary { } } +/// The datatype directory a BIDS suffix belongs in. +/// +/// BIDS files a suffix under one datatype wherever it appears, so this is a +/// property of the suffix alone — the same answer for a raw acquisition, a +/// preprocessed input, and a derivative output. Transmit- and receive-field +/// suffixes (`TB1*`, `RB1*`) and B0 field maps are `fmap`; everything else a +/// qMRI model reads or writes (weighted series, relaxation and MT maps, brain +/// masks) is `anat`. +/// +/// Scope: the suffixes qmrust's own models declare, via `Model::bids`, +/// `Model::bids_outputs` and `Model::required_inputs`. It is not a general BIDS +/// classifier, and the generic field-map suffixes (`epi`, `phasediff`, +/// `phase1`/`phase2`, `magnitude*`) deliberately have no branch here because no +/// registered model names one; the `anat` default would be wrong for them, so +/// a model that starts using one must add it rather than inherit the fallback. +pub fn datatype_for_suffix(suffix: &str) -> &'static str { + if suffix.starts_with("TB1") + || suffix.starts_with("RB1") + || suffix.ends_with("B0map") + || suffix == "fieldmap" + { + "fmap" + } else { + "anat" + } +} + #[cfg(test)] mod tests { use super::*; @@ -258,7 +285,7 @@ mod tests { // Registered model suffixes are built in (compile-time), so known // with no config — but a registered suffix that BIDS itself defines // stays canonical, so it needs no `.bidsignore` exemption. - for canonical in ["IRT1", "MESE", "MTR", "MTS", "VFA"] { + for canonical in ["IRT1", "MESE", "MTR", "MTS", "VFA", "TB1DAM"] { assert!(v.is_known_suffix(canonical)); assert!(!v.is_custom_suffix(canonical), "{canonical} is canonical"); } @@ -270,6 +297,22 @@ mod tests { assert!(!v.is_datatype("notadatatype")); } + #[test] + fn field_map_suffixes_are_filed_under_fmap_everything_else_under_anat() { + // The rule is keyed on the suffix alone, so a transmit-field + // acquisition, its fitted map, and a supplied B1 map all agree — a + // model that writes TB1map into `anat/` would be invisible to the B1 + // input resolver that reads `fmap/`. + for fieldmap in ["TB1DAM", "TB1map", "TB1EPI", "RB1COR", "B0map", "fieldmap"] { + assert_eq!(datatype_for_suffix(fieldmap), "fmap", "{fieldmap}"); + } + for anat in [ + "IRT1", "MESE", "MTR", "MTS", "VFA", "QMTSPGR", "T1map", "MTsat", + ] { + assert_eq!(datatype_for_suffix(anat), "anat", "{anat}"); + } + } + #[test] fn from_config_folds_registered_and_declared_customs() { let v = Vocabulary::from_config(&default_config()); diff --git a/docs/figures/b1_afi/inputs.webp b/docs/figures/b1_afi/inputs.webp new file mode 100644 index 0000000..0f12074 Binary files /dev/null and b/docs/figures/b1_afi/inputs.webp differ diff --git a/docs/figures/b1_afi/outputs.webp b/docs/figures/b1_afi/outputs.webp new file mode 100644 index 0000000..88fda8a Binary files /dev/null and b/docs/figures/b1_afi/outputs.webp differ diff --git a/docs/figures/b1_dam/curve.webp b/docs/figures/b1_dam/curve.webp new file mode 100644 index 0000000..9b79e55 Binary files /dev/null and b/docs/figures/b1_dam/curve.webp differ diff --git a/docs/figures/b1_dam/inputs.webp b/docs/figures/b1_dam/inputs.webp new file mode 100644 index 0000000..fb5a96a Binary files /dev/null and b/docs/figures/b1_dam/inputs.webp differ diff --git a/docs/figures/b1_dam/outputs.webp b/docs/figures/b1_dam/outputs.webp new file mode 100644 index 0000000..c522b4d Binary files /dev/null and b/docs/figures/b1_dam/outputs.webp differ diff --git a/docs/figures/inversion_recovery/inputs.webp b/docs/figures/inversion_recovery/inputs.webp index 67e7a96..07bc8dd 100644 Binary files a/docs/figures/inversion_recovery/inputs.webp and b/docs/figures/inversion_recovery/inputs.webp differ diff --git a/docs/figures/mono_t2/inputs.webp b/docs/figures/mono_t2/inputs.webp index 6c432c0..5222bf1 100644 Binary files a/docs/figures/mono_t2/inputs.webp and b/docs/figures/mono_t2/inputs.webp differ diff --git a/docs/figures/mt_ratio/inputs.webp b/docs/figures/mt_ratio/inputs.webp index 58f4cee..f42417a 100644 Binary files a/docs/figures/mt_ratio/inputs.webp and b/docs/figures/mt_ratio/inputs.webp differ diff --git a/docs/figures/mt_sat/inputs.webp b/docs/figures/mt_sat/inputs.webp index 467bead..faf888e 100644 Binary files a/docs/figures/mt_sat/inputs.webp and b/docs/figures/mt_sat/inputs.webp differ diff --git a/docs/figures/vfa_t1/inputs.webp b/docs/figures/vfa_t1/inputs.webp index 7f902b0..68cd4c5 100644 Binary files a/docs/figures/vfa_t1/inputs.webp and b/docs/figures/vfa_t1/inputs.webp differ diff --git a/docs/index.md b/docs/index.md index 4d3c609..a4e3512 100644 --- a/docs/index.md +++ b/docs/index.md @@ -17,61 +17,94 @@ data in this page, in your browser. ## Methods -### Magnetization transfer +### T1 Relaxometry ::::{grid} 1 1 2 2 :::{card} -:header: Magnetization transfer ratio -:link: models/magnetization-transfer/mt_ratio.md +:header: Inversion Recovery +:link: models/t1-relaxometry/inversion_recovery.md -Computes the magnetization transfer ratio from two images: one acquired with an off-resonance saturation pulse and one without. +Fits the longitudinal relaxation time T1 from a series of inversion-recovery images acquired at different inversion times. -`MTR` · `MTRmap` +`IRT1` · `T1map` ::: :::{card} -:header: MT saturation -:link: models/magnetization-transfer/mt_sat.md +:header: Variable Flip Angle +:link: models/t1-relaxometry/vfa_t1.md -Derives the MT saturation parameter from three spoiled gradient-echo volumes — MT-weighted, PD-weighted and T1-weighted. +Fits the longitudinal relaxation time T1 from spoiled gradient-echo images acquired at two or more excitation flip angles and a single repetition time. -`MTS` · `MTsat`, `T1map`, `MTRmap` +`VFA` · `T1map`, `M0map` ::: +:::: + +### T2 Relaxometry + +::::{grid} 1 1 2 2 :::{card} -:header: qMT-SPGR -:link: models/magnetization-transfer/qmt_spgr.md +:header: Monoexp T2 +:link: models/t2-relaxometry/mono_t2.md -Two-pool quantitative magnetization transfer from a spoiled gradient-echo sequence with off-resonance saturation sampled across a grid of saturation flip angles and frequency offsets. +Fits the transverse relaxation time T2 from a multi-echo spin-echo series as a mono-exponential decay with a free amplitude. -`QMTSPGR` · `Fmap`, `kRmap`, `R1Fmap`, `R1Rmap`, `T2Fmap`, `T2Rmap` +`MESE` · `T2map`, `M0map` ::: :::: -### Relaxometry +### Field Mapping ::::{grid} 1 1 2 2 :::{card} -:header: Mono-exponential T2 -:link: models/relaxometry/mono_t2.md +:header: Actual Flip Angle B1+ +:link: models/field-mapping/b1_afi.md -Fits the transverse relaxation time T2 from a multi-echo spin-echo series as a mono-exponential decay with a free amplitude. +Maps the transmit (B1+) field from a single spoiled gradient-echo sequence that interleaves two excitation repetition times at one nominal flip angle. -`MESE` · `T2map`, `M0map` +`TB1AFI` · `TB1map` ::: :::{card} -:header: Inversion recovery T1 -:link: models/relaxometry/inversion_recovery.md +:header: Double Angle B1+ +:link: models/field-mapping/b1_dam.md -Fits the longitudinal relaxation time T1 from a series of inversion-recovery images acquired at different inversion times. +Maps the transmit (B1+) field from two spoiled gradient-echo volumes acquired at flip angles alpha and twice alpha. -`IRT1` · `T1map` +`TB1DAM` · `TB1map` ::: +:::: + +### Magnetization Transfer + +#### Semi-quantitative MT + +::::{grid} 1 1 2 2 :::{card} -:header: Variable flip angle T1 -:link: models/relaxometry/vfa_t1.md +:header: MT Ratio +:link: models/magnetization-transfer/mt_ratio.md -Fits the longitudinal relaxation time T1 from spoiled gradient-echo images acquired at two or more excitation flip angles and a single repetition time. +Computes the magnetization transfer ratio from two images: one acquired with an off-resonance saturation pulse and one without. -`VFA` · `T1map`, `M0map` +`MTR` · `MTRmap` +::: +:::{card} +:header: MT Saturation +:link: models/magnetization-transfer/mt_sat.md + +Derives the MT saturation parameter from three spoiled gradient-echo volumes — MT-weighted, PD-weighted and T1-weighted. + +`MTS` · `MTsat`, `T1map`, `MTRmap` +::: +:::: + +#### Quantitative MT + +::::{grid} 1 1 2 2 +:::{card} +:header: qMT-SPGR +:link: models/magnetization-transfer/qmt_spgr.md + +Two-pool quantitative magnetization transfer from a spoiled gradient-echo sequence with off-resonance saturation sampled across a grid of saturation flip angles and frequency offsets. + +`QMTSPGR` · `Fmap`, `kRmap`, `R1Fmap`, `R1Rmap`, `T2Fmap`, `T2Rmap` ::: :::: diff --git a/docs/models/field-mapping/b1_afi.md b/docs/models/field-mapping/b1_afi.md new file mode 100644 index 0000000..77fcfe4 --- /dev/null +++ b/docs/models/field-mapping/b1_afi.md @@ -0,0 +1,196 @@ +--- +title: Actual Flip Angle B1+ +subtitle: Field Mapping · BIDS suffix TB1AFI +keywords: + - Field Mapping + - TB1AFI +thumbnail: ../../figures/b1_afi/outputs.webp +--- + + + +Config name +: `model: b1_afi` +BIDS suffix +: `TB1AFI` +Method family +: Field Mapping +Acquired volumes +: 2 +Fit granularity +: voxelwise + +Maps the transmit (B1+) field from a single spoiled gradient-echo sequence that interleaves two excitation repetition times at one nominal flip angle. In the pulsed steady state the ratio of the two signals depends on the achieved flip angle alone, so proton density and receive sensitivity cancel. The result is dimensionless — the achieved flip angle as a fraction of the nominal one — and feeds the B1 correction of models such as MTsat and VFA. The closed form assumes the repetition times are short against T1, and underestimates B1 as that assumption weakens. + +## Signal model + +```{math} +:label: eq-b1_afi +\cos\theta = \frac{r n - 1}{n - r},\quad r = \left|\frac{S_{\mathrm{TR}_2}}{S_{\mathrm{TR}_1}}\right|,\quad n = \frac{\mathrm{TR}_2}{\mathrm{TR}_1},\quad B_1^{+} = \frac{\theta}{\theta_\mathrm{nom}} +``` + +| Symbol | Quantity | Unit | +|---|---|---| +| `B1` | Relative transmit field | — | +| `T1` | Longitudinal relaxation time | s | + +## What the model reads + +A series of 2 volumes indexed by `RepetitionTimeExcitation`. Volumes are matched to these identities by value, so reordering the acquisition yields an identical fit. + +| # | RepetitionTimeExcitation | +|---|---| +| 1 | 0.02 | +| 2 | 0.1 | + +### Where the protocol comes from + +In a BIDS dataset these values are resolved for you: + +| Protocol parameter | Source | Key | Resolved | +|---|---|---|---| +| `RepetitionTimeExcitation` | sidecar field | `RepetitionTimeExcitation` | per volume | +| `FlipAngle` | sidecar field | `FlipAngle` | whole collection | + +## Inputs + +```{figure} ../../figures/b1_afi/inputs.webp +:label: fig-b1_afi-inputs +The 2 acquired volumes, each labeled with the identity the model reads it by. +``` + +## Outputs + +| Map | BIDS suffix | Unit | Fit bounds | Fitted | +|---|---|---|---|---| +| `B1` | `TB1map` | — | 0 … 2 | free | + +```{figure} ../../figures/b1_afi/outputs.webp +:label: fig-b1_afi-outputs +Fitted maps for one example subject. +``` + +## Usage + +::::{tab-set} +:::{tab-item} BIDS +:sync: bids + +The acquisition comes from the sidecars, so `recipes/bids/b1_afi_config.yaml` carries only algorithm options: + +```yaml +model: b1_afi + +# B1 = acos((r*n - 1) / (n - r)) / nomFA, with r the signal ratio of the two +# interleaved repetition times, n = TR2/TR1, and the map dimensionless +# (1.0 = the nominal flip angle was achieved). BIDS input: +# RepetitionTimeExcitation and FlipAngle are read from the JSON sidecars via +# the model's protocol_schema(), so neither is declared here. The TB1AFI +# collection is assembled by the `acq` entity (acq-tr1 / acq-tr2). Example: +# +# qmrust fit --bids-dir \ +# --config recipes/bids/b1_afi_config.yaml --output-dir out + +# BIDS mask selection. Left blank: this example dataset ships no mask, and +# naming one here would claim otherwise. Blank means "whatever mask is +# present", so a dataset of your own that has one is masked by it. +mask: + desc: "" +``` + +```bash +qmrust fit --bids-dir ds-mydata \ + --config recipes/bids/b1_afi_config.yaml --output-dir out +``` + +Fitted maps are written as BIDS derivatives: + +```text +out/derivatives/qmrust/sub-01/fmap/ + sub-01_TB1map.nii.gz + sub-01_TB1map.json +``` +::: + +:::{tab-item} Without BIDS +:sync: non-bids + +No sidecars, so `recipes/non-bids/b1_afi_config.yaml` declares the acquisition itself: + +```yaml +model: b1_afi + +# non-BIDS: protocol from this file; mask via --mask. +# Times in seconds (BIDS/SI), flip angle in degrees (BIDS-MRI convention). +# qMRLab's b1_afi default protocol is nomFA = 60 deg with TR1/TR2 of 20/100 ms; +# its Sequence matrix records those two in milliseconds, so they are converted +# once here at the shell boundary. Only their ratio enters the estimate. +repetition_times: +- 0.02 +- 0.10 +flip_angle: 60 + +# There are no fit options: the estimate is closed form. +# +# Two uses: +# - bidsify a qMRLab dataset (AFIData1.nii.gz / AFIData2.nii.gz) by naming +# one --nii-data per volume, in acquisition order: +# +# qmrust bidsify --model b1_afi \ +# --nii-data AFIData1.nii.gz --nii-data AFIData2.nii.gz \ +# --config recipes/non-bids/b1_afi_config.yaml --subject 01 --out ds +# +# - fit a stacked 4D NIfTI whose two volumes are [TR1, TR2]: +# +# qmrust fit --data afi.nii.gz --mask Mask.nii.gz \ +# --config recipes/non-bids/b1_afi_config.yaml --output-dir out_b1 +``` + +```bash +qmrust fit --data data.nii.gz --mask mask.nii.gz \ + --config recipes/non-bids/b1_afi_config.yaml \ + --output-dir out +``` + +See [Fitting without BIDS](../../guide/non-bids.md). +::: + +:::{tab-item} Browser +:sync: wasm + +The same core, compiled to WebAssembly: + +```js +import init, { fit_volume } from './pkg/qmrust_wasm.js'; +await init(); +const maps = fit_volume(configYaml, data, [nx, ny, nz, nt], + JSON.stringify(volumeIds), mask, '{}'); +``` + +Try it on real data in the [playground](../../playground.md), or read [Browser & wasm](../../guide/browser.md). +::: +:::: + +## Implementation + +- [`model.rs`](https://github.com/qmrlab/qmrust/blob/main/crates/qmrust-core/src/models/b1_afi/model.rs) — the `Model` implementation: measurement contract, protocol schema, BIDS outputs. +- [`fit.rs`](https://github.com/qmrlab/qmrust/blob/main/crates/qmrust-core/src/models/b1_afi/fit.rs) — the numerical fit. +- [`config.rs`](https://github.com/qmrlab/qmrust/blob/main/crates/qmrust-core/src/models/b1_afi/config.rs) — configuration and validation. + +:::{dropdown} Fully-resolved configuration +Every default materialized, as `qmrust dump-config` prints it. + +```yaml +model: b1_afi +repetition_times: +- 0.02 +- 0.1 +flip_angle: 60.0 +``` +::: + +## References + +@yarnykh2007 diff --git a/docs/models/field-mapping/b1_dam.md b/docs/models/field-mapping/b1_dam.md new file mode 100644 index 0000000..a34d594 --- /dev/null +++ b/docs/models/field-mapping/b1_dam.md @@ -0,0 +1,202 @@ +--- +title: Double Angle B1+ +subtitle: Field Mapping · BIDS suffix TB1DAM +keywords: + - Field Mapping + - TB1DAM +thumbnail: ../../figures/b1_dam/outputs.webp +--- + + + +Config name +: `model: b1_dam` +BIDS suffix +: `TB1DAM` +Method family +: Field Mapping +Acquired volumes +: 2 +Fit granularity +: voxelwise + +Maps the transmit (B1+) field from two spoiled gradient-echo volumes acquired at flip angles alpha and twice alpha. The ratio of the two signals depends only on the achieved flip angle, so proton density and receive sensitivity cancel and the map needs no calibration. The result is dimensionless — the achieved flip angle as a fraction of the nominal one — and feeds the B1 correction of models such as MTsat and VFA. + +## Signal model + +```{math} +:label: eq-b1_dam +B_1^{+} = \frac{\left|\arccos\left(\frac{S_{2\alpha}}{2 S_{\alpha}}\right)\right|}{\alpha} +``` + +| Symbol | Quantity | Unit | +|---|---|---| +| `B1` | Relative transmit field | — | + +## What the model reads + +A series of 2 volumes indexed by `FlipAngle`. Volumes are matched to these identities by value, so reordering the acquisition yields an identical fit. + +| # | FlipAngle | +|---|---| +| 1 | 60 | +| 2 | 120 | + +### Where the protocol comes from + +In a BIDS dataset these values are resolved for you: + +| Protocol parameter | Source | Key | Resolved | +|---|---|---|---| +| `FlipAngle` | sidecar field | `FlipAngle` | per volume | + +## Inputs + +```{figure} ../../figures/b1_dam/inputs.webp +:label: fig-b1_dam-inputs +The 2 acquired volumes, each labeled with the identity the model reads it by. +``` + +## Outputs + +| Map | BIDS suffix | Unit | Fit bounds | Fitted | +|---|---|---|---|---| +| `B1` | `TB1map` | — | 0 … 2 | free | + +:::{dropdown} Diagnostic outputs +Reported by the fit but not written as BIDS derivatives: they describe the fit rather than the tissue, and have no BIDS suffix to be written under. + +| Output | +|---| +| `A` | +::: + +```{figure} ../../figures/b1_dam/outputs.webp +:label: fig-b1_dam-outputs +Fitted maps for one example subject. +``` + +## Usage + +::::{tab-set} +:::{tab-item} BIDS +:sync: bids + +The acquisition comes from the sidecars, so `recipes/bids/b1_dam_config.yaml` carries only algorithm options: + +```yaml +model: b1_dam + +# B1 = |acos(S(2a) / (2 * S(a)))| / a, with the map dimensionless (1.0 = the +# nominal flip angle was achieved). BIDS input: the two flip angles are read +# from the JSON sidecars via the model's protocol_schema(), so they are not +# declared here. The TB1DAM collection is assembled by flip index. Example: +# +# qmrust fit --bids-dir \ +# --config recipes/bids/b1_dam_config.yaml --output-dir out + +# BIDS mask selection. Left blank: this example dataset ships no mask, and +# naming one here would claim otherwise. Blank means "whatever mask is +# present", so a dataset of your own that has one is masked by it. +mask: + desc: "" +``` + +```bash +qmrust fit --bids-dir ds-mydata \ + --config recipes/bids/b1_dam_config.yaml --output-dir out +``` + +Fitted maps are written as BIDS derivatives: + +```text +out/derivatives/qmrust/sub-01/fmap/ + sub-01_TB1map.nii.gz + sub-01_TB1map.json +``` +::: + +:::{tab-item} Without BIDS +:sync: non-bids + +No sidecars, so `recipes/non-bids/b1_dam_config.yaml` declares the acquisition itself: + +```yaml +model: b1_dam + +# non-BIDS: protocol from this file; mask via --mask. +# Nominal flip angles in degrees (BIDS-MRI convention), in acquisition order. +# The double-angle identity only holds when the second angle is twice the +# first, so the pair is validated before the fit runs. qMRLab's b1_dam default +# protocol is alpha = 60 degrees. +flip_angles: +- 60 +- 120 + +# There are no fit options: the estimate is closed form. +# +# Two uses: +# - bidsify a qMRLab dataset (SFalpha.mat / SF2alpha.mat plus Mask.mat) from +# --mat-dir; this config is the --config input: +# +# qmrust bidsify --model b1_dam --mat-dir \ +# --config recipes/non-bids/b1_dam_config.yaml --subject 01 --out ds +# +# - fit a stacked 4D NIfTI whose two volumes are [alpha, 2*alpha]: +# +# qmrust fit --data sfalpha_sf2alpha.nii.gz --mask Mask.mat \ +# --config recipes/non-bids/b1_dam_config.yaml --output-dir out_b1 +``` + +```bash +qmrust fit --data data.nii.gz --mask mask.nii.gz \ + --config recipes/non-bids/b1_dam_config.yaml \ + --output-dir out +``` + +See [Fitting without BIDS](../../guide/non-bids.md). +::: + +:::{tab-item} Browser +:sync: wasm + +The same core, compiled to WebAssembly: + +```js +import init, { fit_volume } from './pkg/qmrust_wasm.js'; +await init(); +const maps = fit_volume(configYaml, data, [nx, ny, nz, nt], + JSON.stringify(volumeIds), mask, '{}'); +``` + +Try it on real data in the [playground](../../playground.md), or read [Browser & wasm](../../guide/browser.md). +::: +:::: + +```{figure} ../../figures/b1_dam/curve.webp +:label: fig-b1_dam-curve +One voxel's measured data with the model's forward signal at the fitted parameters. The curve is computed by the same Rust code that does the fitting. +``` + +## Implementation + +- [`model.rs`](https://github.com/qmrlab/qmrust/blob/main/crates/qmrust-core/src/models/b1_dam/model.rs) — the `Model` implementation: measurement contract, protocol schema, BIDS outputs. +- [`fit.rs`](https://github.com/qmrlab/qmrust/blob/main/crates/qmrust-core/src/models/b1_dam/fit.rs) — the numerical fit. +- [`config.rs`](https://github.com/qmrlab/qmrust/blob/main/crates/qmrust-core/src/models/b1_dam/config.rs) — configuration and validation. + +:::{dropdown} Fully-resolved configuration +Every default materialized, as `qmrust dump-config` prints it. + +```yaml +model: b1_dam +flip_angles: +- 60.0 +- 120.0 +``` +::: + +## References + +@insko1993 diff --git a/docs/models/index.md b/docs/models/index.md index a00bd31..adb240b 100644 --- a/docs/models/index.md +++ b/docs/models/index.md @@ -6,60 +6,93 @@ Every model qmrust fits, grouped by method family. Each page documents what the model reads, what it produces, and how to run it with and without BIDS. -## Magnetization transfer +## T1 Relaxometry ::::{grid} 1 1 2 2 :::{card} -:header: Magnetization transfer ratio -:link: magnetization-transfer/mt_ratio.md +:header: Inversion Recovery +:link: t1-relaxometry/inversion_recovery.md -Computes the magnetization transfer ratio from two images: one acquired with an off-resonance saturation pulse and one without. +Fits the longitudinal relaxation time T1 from a series of inversion-recovery images acquired at different inversion times. -`MTR` · `MTRmap` +`IRT1` · `T1map` ::: :::{card} -:header: MT saturation -:link: magnetization-transfer/mt_sat.md +:header: Variable Flip Angle +:link: t1-relaxometry/vfa_t1.md -Derives the MT saturation parameter from three spoiled gradient-echo volumes — MT-weighted, PD-weighted and T1-weighted. +Fits the longitudinal relaxation time T1 from spoiled gradient-echo images acquired at two or more excitation flip angles and a single repetition time. -`MTS` · `MTsat`, `T1map`, `MTRmap` +`VFA` · `T1map`, `M0map` ::: +:::: + +## T2 Relaxometry + +::::{grid} 1 1 2 2 :::{card} -:header: qMT-SPGR -:link: magnetization-transfer/qmt_spgr.md +:header: Monoexp T2 +:link: t2-relaxometry/mono_t2.md -Two-pool quantitative magnetization transfer from a spoiled gradient-echo sequence with off-resonance saturation sampled across a grid of saturation flip angles and frequency offsets. +Fits the transverse relaxation time T2 from a multi-echo spin-echo series as a mono-exponential decay with a free amplitude. -`QMTSPGR` · `Fmap`, `kRmap`, `R1Fmap`, `R1Rmap`, `T2Fmap`, `T2Rmap` +`MESE` · `T2map`, `M0map` ::: :::: -## Relaxometry +## Field Mapping ::::{grid} 1 1 2 2 :::{card} -:header: Mono-exponential T2 -:link: relaxometry/mono_t2.md +:header: Actual Flip Angle B1+ +:link: field-mapping/b1_afi.md -Fits the transverse relaxation time T2 from a multi-echo spin-echo series as a mono-exponential decay with a free amplitude. +Maps the transmit (B1+) field from a single spoiled gradient-echo sequence that interleaves two excitation repetition times at one nominal flip angle. -`MESE` · `T2map`, `M0map` +`TB1AFI` · `TB1map` ::: :::{card} -:header: Inversion recovery T1 -:link: relaxometry/inversion_recovery.md +:header: Double Angle B1+ +:link: field-mapping/b1_dam.md -Fits the longitudinal relaxation time T1 from a series of inversion-recovery images acquired at different inversion times. +Maps the transmit (B1+) field from two spoiled gradient-echo volumes acquired at flip angles alpha and twice alpha. -`IRT1` · `T1map` +`TB1DAM` · `TB1map` ::: +:::: + +## Magnetization Transfer + +### Semi-quantitative MT + +::::{grid} 1 1 2 2 :::{card} -:header: Variable flip angle T1 -:link: relaxometry/vfa_t1.md +:header: MT Ratio +:link: magnetization-transfer/mt_ratio.md -Fits the longitudinal relaxation time T1 from spoiled gradient-echo images acquired at two or more excitation flip angles and a single repetition time. +Computes the magnetization transfer ratio from two images: one acquired with an off-resonance saturation pulse and one without. -`VFA` · `T1map`, `M0map` +`MTR` · `MTRmap` +::: +:::{card} +:header: MT Saturation +:link: magnetization-transfer/mt_sat.md + +Derives the MT saturation parameter from three spoiled gradient-echo volumes — MT-weighted, PD-weighted and T1-weighted. + +`MTS` · `MTsat`, `T1map`, `MTRmap` +::: +:::: + +### Quantitative MT + +::::{grid} 1 1 2 2 +:::{card} +:header: qMT-SPGR +:link: magnetization-transfer/qmt_spgr.md + +Two-pool quantitative magnetization transfer from a spoiled gradient-echo sequence with off-resonance saturation sampled across a grid of saturation flip angles and frequency offsets. + +`QMTSPGR` · `Fmap`, `kRmap`, `R1Fmap`, `R1Rmap`, `T2Fmap`, `T2Rmap` ::: :::: diff --git a/docs/models/magnetization-transfer/mt_ratio.md b/docs/models/magnetization-transfer/mt_ratio.md index b186535..073cb13 100644 --- a/docs/models/magnetization-transfer/mt_ratio.md +++ b/docs/models/magnetization-transfer/mt_ratio.md @@ -1,8 +1,8 @@ --- -title: Magnetization transfer ratio -subtitle: Magnetization transfer · BIDS suffix MTR +title: MT Ratio +subtitle: Semi-quantitative MT · BIDS suffix MTR keywords: - - Magnetization transfer + - Semi-quantitative MT - MTR thumbnail: ../../figures/mt_ratio/outputs.webp --- @@ -16,7 +16,7 @@ Config name BIDS suffix : `MTR` Method family -: Magnetization transfer +: Semi-quantitative MT Acquired volumes : 2 Fit granularity diff --git a/docs/models/magnetization-transfer/mt_sat.md b/docs/models/magnetization-transfer/mt_sat.md index 617e2ec..ef90ca9 100644 --- a/docs/models/magnetization-transfer/mt_sat.md +++ b/docs/models/magnetization-transfer/mt_sat.md @@ -1,8 +1,8 @@ --- -title: MT saturation -subtitle: Magnetization transfer · BIDS suffix MTS +title: MT Saturation +subtitle: Semi-quantitative MT · BIDS suffix MTS keywords: - - Magnetization transfer + - Semi-quantitative MT - MTS thumbnail: ../../figures/mt_sat/outputs.webp --- @@ -16,7 +16,7 @@ Config name BIDS suffix : `MTS` Method family -: Magnetization transfer +: Semi-quantitative MT Acquired volumes : 3 Fit granularity @@ -102,10 +102,11 @@ model: mt_sat b1_correction_factor: 0.4 export_mtr: true -# BIDS mask selection — disambiguates which mask to apply. If the dataset has -# no matching mask, no masking is applied (as in this example). +# BIDS mask selection. Left blank: this example dataset ships no mask, and +# naming one here would claim otherwise. Blank means "whatever mask is +# present", so a dataset of your own that has one is masked by it. mask: - desc: brain + desc: "" ``` ```bash diff --git a/docs/models/magnetization-transfer/qmt_spgr.md b/docs/models/magnetization-transfer/qmt_spgr.md index 29ad7e7..9df7731 100644 --- a/docs/models/magnetization-transfer/qmt_spgr.md +++ b/docs/models/magnetization-transfer/qmt_spgr.md @@ -1,8 +1,8 @@ --- title: qMT-SPGR -subtitle: Magnetization transfer · BIDS suffix QMTSPGR +subtitle: Quantitative MT · BIDS suffix QMTSPGR keywords: - - Magnetization transfer + - Quantitative MT - QMTSPGR thumbnail: ../../figures/qmt_spgr/outputs.webp --- @@ -16,7 +16,7 @@ Config name BIDS suffix : `QMTSPGR` Method family -: Magnetization transfer +: Quantitative MT Acquired volumes : 10 Fit granularity @@ -102,7 +102,7 @@ Auxiliary maps and the mask applied to this fit. | `T2r` | `T2Rmap` | s | 3e-06 … 5e-05 | free | :::{dropdown} Diagnostic outputs -Written alongside the quantitative maps, but not quantitative measurements — they describe the fit, not the tissue. +Reported by the fit but not written as BIDS derivatives: they describe the fit rather than the tissue, and have no BIDS suffix to be written under. | Output | |---| diff --git a/docs/models/relaxometry/inversion_recovery.md b/docs/models/t1-relaxometry/inversion_recovery.md similarity index 95% rename from docs/models/relaxometry/inversion_recovery.md rename to docs/models/t1-relaxometry/inversion_recovery.md index 805c5c7..14256b5 100644 --- a/docs/models/relaxometry/inversion_recovery.md +++ b/docs/models/t1-relaxometry/inversion_recovery.md @@ -1,8 +1,8 @@ --- -title: Inversion recovery T1 -subtitle: Relaxometry · BIDS suffix IRT1 +title: Inversion Recovery +subtitle: T1 Relaxometry · BIDS suffix IRT1 keywords: - - Relaxometry + - T1 Relaxometry - IRT1 thumbnail: ../../figures/inversion_recovery/outputs.webp --- @@ -16,7 +16,7 @@ Config name BIDS suffix : `IRT1` Method family -: Relaxometry +: T1 Relaxometry Acquired volumes : 9 Fit granularity @@ -81,7 +81,7 @@ Auxiliary maps and the mask applied to this fit. | `T1` | `T1map` | s | unbounded … unbounded | free | :::{dropdown} Diagnostic outputs -Written alongside the quantitative maps, but not quantitative measurements — they describe the fit, not the tissue. +Reported by the fit but not written as BIDS derivatives: they describe the fit rather than the tissue, and have no BIDS suffix to be written under. | Output | |---| diff --git a/docs/models/relaxometry/vfa_t1.md b/docs/models/t1-relaxometry/vfa_t1.md similarity index 98% rename from docs/models/relaxometry/vfa_t1.md rename to docs/models/t1-relaxometry/vfa_t1.md index 38cd100..0ec8a3a 100644 --- a/docs/models/relaxometry/vfa_t1.md +++ b/docs/models/t1-relaxometry/vfa_t1.md @@ -1,8 +1,8 @@ --- -title: Variable flip angle T1 -subtitle: Relaxometry · BIDS suffix VFA +title: Variable Flip Angle +subtitle: T1 Relaxometry · BIDS suffix VFA keywords: - - Relaxometry + - T1 Relaxometry - VFA thumbnail: ../../figures/vfa_t1/outputs.webp --- @@ -16,7 +16,7 @@ Config name BIDS suffix : `VFA` Method family -: Relaxometry +: T1 Relaxometry Acquired volumes : 2 Fit granularity diff --git a/docs/models/relaxometry/mono_t2.md b/docs/models/t2-relaxometry/mono_t2.md similarity index 98% rename from docs/models/relaxometry/mono_t2.md rename to docs/models/t2-relaxometry/mono_t2.md index ef8f703..f083448 100644 --- a/docs/models/relaxometry/mono_t2.md +++ b/docs/models/t2-relaxometry/mono_t2.md @@ -1,8 +1,8 @@ --- -title: Mono-exponential T2 -subtitle: Relaxometry · BIDS suffix MESE +title: Monoexp T2 +subtitle: T2 Relaxometry · BIDS suffix MESE keywords: - - Relaxometry + - T2 Relaxometry - MESE thumbnail: ../../figures/mono_t2/outputs.webp --- @@ -16,7 +16,7 @@ Config name BIDS suffix : `MESE` Method family -: Relaxometry +: T2 Relaxometry Acquired volumes : 30 Fit granularity diff --git a/docs/playground/app.css b/docs/playground/app.css index 5b27181..7e67c9c 100644 --- a/docs/playground/app.css +++ b/docs/playground/app.css @@ -291,6 +291,21 @@ button.primary:hover:not(:disabled) { filter: brightness(1.08); } .menu-item:focus-visible { outline: 2px solid var(--accent); outline-offset: -2px; } .menu-item .icon { color: var(--accent); flex: none; } +/* Shrink/grow sits in the same menu as the segmentation methods but is not a + command: it holds two buttons and a running description, so it opts out of + the row hover and keeps its label from pushing the steppers apart. */ +.menu-item.mask-steps { cursor: default; gap: 6px; } +.menu-item.mask-steps:hover { background: none; color: var(--ink-2); } +.menu-item.mask-steps span:not(.icon) { flex: 1; } +.step-btn { + flex: none; width: 20px; height: 20px; padding: 0; + border: var(--border-w) solid var(--line); border-radius: 5px; + background: var(--field); color: var(--ink); + font: 600 13px/1 var(--font-ui); cursor: pointer; +} +.step-btn:hover { border-color: var(--accent); color: var(--accent); } +.step-btn:focus-visible { outline: 2px solid var(--accent); outline-offset: 1px; } + /* Fitting progress: a permanent fixture above the Fit button (idle = empty track), not a topbar element that appears/disappears and shifts the rest of the page. The stripe is a white sheen rather than a token: it reads as a @@ -467,6 +482,102 @@ button.primary:hover:not(:disabled) { filter: brightness(1.08); } .model-line { margin: 12px 0; color: var(--ink-2); } .model-line select { margin-left: 4px; } +/* Model picker: a disclosure tree over the registry's taxonomy. Anchored + relative so the tree hangs under the button without a fixed-position + measurement, since this one always has room below it inside the panel. */ +.model-label { + display: flex; align-items: center; gap: 6px; margin-bottom: 6px; + color: var(--muted); font: 700 10px var(--font-ui); + text-transform: uppercase; letter-spacing: .08em; +} +.model-label .icon { color: var(--accent); } + +.model-picker { position: relative; } +.model-button { + display: flex; align-items: center; gap: 8px; width: 100%; + padding: 8px 10px; border: var(--border-w) solid var(--line); border-radius: 8px; + background: var(--field); color: var(--ink); + font: 600 13px var(--font-ui); text-align: left; cursor: pointer; +} +.model-button:hover { border-color: var(--accent); } +.model-button:focus-visible { outline: 2px solid var(--accent); outline-offset: 1px; } +.model-button .icon { color: var(--muted); flex: none; transition: transform .15s ease; } +.model-button[aria-expanded="true"] .icon { transform: rotate(90deg); color: var(--accent); } +.model-button-title { flex: 1; } +.model-suffix { + flex: none; font: 600 10px var(--font-mono); letter-spacing: .04em; + padding: 2px 6px; border-radius: 999px; + background: var(--panel-2); border: 1px solid var(--line); color: var(--muted); +} + +.model-tree { + position: absolute; z-index: 60; left: 0; right: 0; top: calc(100% + 4px); + max-height: 60vh; overflow-y: auto; + display: flex; flex-direction: column; gap: 1px; padding: 4px; + /* Opaque, unlike the shared panel token. This hangs directly over the recipe + form, and every theme's `--panel` is translucent (verdigris is .80), so the + fields underneath read straight through the list. The tint is kept rather + than dropped: the panel colour is layered over an opaque base instead of + being the base. */ + background-color: var(--bg); + background-image: linear-gradient(var(--panel), var(--panel)), var(--grad-soft); + border: var(--border-w) solid var(--line); border-radius: 10px; box-shadow: var(--shadow); +} +.model-tree[hidden] { display: none; } + +.model-family { + display: flex; align-items: center; gap: 6px; + padding: 6px 8px; border: 0; border-radius: 6px; + background: none; color: var(--ink); cursor: pointer; + font: 700 11px var(--font-ui); text-transform: uppercase; letter-spacing: .05em; + text-align: left; +} +.model-family:hover { background: var(--field); } +.model-family:focus-visible { outline: 2px solid var(--accent); outline-offset: -2px; } +.model-family span:not(.model-count) { flex: 1; } +/* The chevron is the fold's whole state: right is closed, down is open. */ +/* Two glyphs: the fold chevron turns, the family's own icon does not. */ +.model-family .icon { flex: none; } +.model-family .icon:first-child { color: var(--muted); transition: transform .15s ease; } +.model-family.open .icon:first-child { transform: rotate(90deg); color: var(--accent); } +.model-family .icon:nth-child(2) { color: var(--accent); } +.model-count { + flex: none; min-width: 16px; text-align: center; + font: 600 10px var(--font-mono); color: var(--muted); +} + +/* Nesting is shown the way the file tree shows it: one indent step per level + with a guide line down its left edge, so the two trees on this page read the + same way rather than each inventing its own depth cue. */ +.model-family-body { + display: flex; flex-direction: column; gap: 1px; + margin-left: 13px; padding-left: 10px; border-left: 1px solid var(--line); +} +.model-subgroup-body { + display: flex; flex-direction: column; gap: 1px; + margin-left: 6px; padding-left: 10px; border-left: 1px solid var(--line); +} +.model-subgroup { + padding: 5px 8px 2px; color: var(--muted); + font: 600 10px var(--font-ui); text-transform: uppercase; letter-spacing: .05em; +} +/* The family holding the loaded model, while everything is folded shut. */ +.model-family.holds-active { color: var(--accent); } +.model-family.holds-active .model-count { color: var(--accent); } +.model-option { + display: flex; align-items: center; gap: 8px; + padding: 6px 10px; border: 0; border-left: 2px solid transparent; border-radius: 6px; + background: none; color: var(--ink-2); cursor: pointer; + font: 500 12px var(--font-ui); text-align: left; +} +.model-option span:first-child { flex: 1; } +.model-option:hover { background: var(--field); color: var(--ink); } +.model-option:focus-visible { outline: 2px solid var(--accent); outline-offset: -2px; } +.model-option.active { + color: var(--ink); background: var(--field); border-left-color: var(--accent); +} +.model-option.active .model-suffix { border-color: var(--accent); color: var(--accent); } + .yaml-wrap { position: relative; } .yaml-pill { position: absolute; top: 8px; right: 8px; font: 700 10px var(--font-mono); diff --git a/docs/playground/app.js b/docs/playground/app.js index ee2b96d..825e4f0 100644 --- a/docs/playground/app.js +++ b/docs/playground/app.js @@ -48,6 +48,7 @@ import { THEMES, initTheme, onThemeChange, setFamily, setMode } from "./themes.j import { createLevelControl, repaintLevels } from "./level.js"; import { fitSlice } from "./fit.js"; import { loadModel } from "./model.js"; +import { buildModelTree, compareByTaxonomy } from "./picker.js"; // --------------------------------------------------------------------------- // Tooltips for the panel info buttons. One element, moved and refilled — a @@ -373,6 +374,7 @@ async function populateModels() { } app.modelNames = index.models; const select = $("model"); + const metas = {}; for (const name of index.models) { let meta; try { @@ -382,12 +384,26 @@ async function populateModels() { status(`Could not load ${name}: ${e.message}`, "error"); continue; } + metas[name] = meta; + } + // Taxonomy order, not the index's alphabetical one: the first option is what + // loads on arrival, and that should be the first model a reader is offered + // rather than whichever name sorts first. + const ordered = Object.keys(metas) + .map((name) => ({ name, meta: metas[name] })) + .sort(compareByTaxonomy) + .map((e) => e.name); + for (const name of ordered) { const opt = document.createElement("option"); opt.value = name; - opt.textContent = meta.title; + opt.textContent = metas[name].title; select.append(opt); } - return select.options.length > 0; + if (select.options.length === 0) return false; + // The tree is a view over the select, so it is built once the options are + // there and reads the taxonomy out of the payloads just loaded. + buildModelTree(metas, (name) => loadModel(name)); + return true; } async function main() { diff --git a/docs/playground/data/b1_afi.json b/docs/playground/data/b1_afi.json new file mode 100644 index 0000000..f7d3619 --- /dev/null +++ b/docs/playground/data/b1_afi.json @@ -0,0 +1,109 @@ +{ + "model": "b1_afi", + "title": "Actual Flip Angle B1+", + "bids_suffix": "TB1AFI", + "family": "Field Mapping", + "family_icon": "circle-dot", + "subgroup": null, + "category_order": 2, + "dims": [ + 88, + 128, + 1, + 2 + ], + "factor": 1, + "volume_ids": [ + { + "RepetitionTimeExcitation": 0.02 + }, + { + "RepetitionTimeExcitation": 0.1 + } + ], + "labels": [ + "RepetitionTimeExcitation=0.02", + "RepetitionTimeExcitation=0.1" + ], + "params": [ + "B1", + "T1" + ], + "outputs": [ + { + "name": "B1", + "unit": "", + "display_range": [ + 0.5, + 1.5 + ] + } + ], + "enums": [], + "config": "model: b1_afi\n\n# non-BIDS: protocol from this file; mask via --mask.\n# Times in seconds (BIDS/SI), flip angle in degrees (BIDS-MRI convention).\n# qMRLab's b1_afi default protocol is nomFA = 60 deg with TR1/TR2 of 20/100 ms;\n# its Sequence matrix records those two in milliseconds, so they are converted\n# once here at the shell boundary. Only their ratio enters the estimate.\nrepetition_times:\n- 0.02\n- 0.10\nflip_angle: 60\n\n# There are no fit options: the estimate is closed form.\n#\n# Two uses:\n# - bidsify a qMRLab dataset (AFIData1.nii.gz / AFIData2.nii.gz) by naming\n# one --nii-data per volume, in acquisition order:\n#\n# qmrust bidsify --model b1_afi \\\n# --nii-data AFIData1.nii.gz --nii-data AFIData2.nii.gz \\\n# --config recipes/non-bids/b1_afi_config.yaml --subject 01 --out ds\n#\n# - fit a stacked 4D NIfTI whose two volumes are [TR1, TR2]:\n#\n# qmrust fit --data afi.nii.gz --mask Mask.nii.gz \\\n# --config recipes/non-bids/b1_afi_config.yaml --output-dir out_b1\n", + "config_bids": "model: b1_afi\n\n# B1 = acos((r*n - 1) / (n - r)) / nomFA, with r the signal ratio of the two\n# interleaved repetition times, n = TR2/TR1, and the map dimensionless\n# (1.0 = the nominal flip angle was achieved). BIDS input:\n# RepetitionTimeExcitation and FlipAngle are read from the JSON sidecars via\n# the model's protocol_schema(), so neither is declared here. The TB1AFI\n# collection is assembled by the `acq` entity (acq-tr1 / acq-tr2). Example:\n#\n# qmrust fit --bids-dir \\\n# --config recipes/bids/b1_afi_config.yaml --output-dir out\n\n# BIDS mask selection. Left blank: this example dataset ships no mask, and\n# naming one here would claim otherwise. Blank means \"whatever mask is\n# present\", so a dataset of your own that has one is masked by it.\nmask:\n desc: \"\"\n", + "archive": "ds-tb1afi.zip", + "files": { + "data": "b1_afi.nii.gz", + "mask": null, + "aux": {} + }, + "probes": [ + { + "x": 44, + "y": 64, + "expected": { + "B1": 1.1303190829700478 + } + }, + { + "x": 41, + "y": 61, + "expected": { + "B1": 1.1335957094117102 + } + }, + { + "x": 41, + "y": 67, + "expected": { + "B1": 1.1318493144155053 + } + }, + { + "x": 47, + "y": 61, + "expected": { + "B1": 1.0979337368652193 + } + }, + { + "x": 47, + "y": 67, + "expected": { + "B1": 1.1429699650570844 + } + }, + { + "x": 44, + "y": 70, + "expected": { + "B1": 1.0617585161170247 + } + }, + { + "x": 50, + "y": 64, + "expected": { + "B1": 1.0750715539956834 + } + }, + { + "x": 38, + "y": 58, + "expected": { + "B1": 0.9942873505087643 + } + } + ] +} \ No newline at end of file diff --git a/docs/playground/data/b1_afi.nii.gz b/docs/playground/data/b1_afi.nii.gz new file mode 100644 index 0000000..04567cc Binary files /dev/null and b/docs/playground/data/b1_afi.nii.gz differ diff --git a/docs/playground/data/b1_dam.json b/docs/playground/data/b1_dam.json new file mode 100644 index 0000000..60440be --- /dev/null +++ b/docs/playground/data/b1_dam.json @@ -0,0 +1,109 @@ +{ + "model": "b1_dam", + "title": "Double Angle B1+", + "bids_suffix": "TB1DAM", + "family": "Field Mapping", + "family_icon": "circle-dot", + "subgroup": null, + "category_order": 2, + "dims": [ + 64, + 64, + 1, + 2 + ], + "factor": 1, + "volume_ids": [ + { + "FlipAngle": 60.0 + }, + { + "FlipAngle": 120.0 + } + ], + "labels": [ + "FlipAngle=60", + "FlipAngle=120" + ], + "params": [ + "B1", + "A" + ], + "outputs": [ + { + "name": "B1", + "unit": "", + "display_range": [ + 0.5, + 1.5 + ] + } + ], + "enums": [], + "config": "model: b1_dam\n\n# non-BIDS: protocol from this file; mask via --mask.\n# Nominal flip angles in degrees (BIDS-MRI convention), in acquisition order.\n# The double-angle identity only holds when the second angle is twice the\n# first, so the pair is validated before the fit runs. qMRLab's b1_dam default\n# protocol is alpha = 60 degrees.\nflip_angles:\n- 60\n- 120\n\n# There are no fit options: the estimate is closed form.\n#\n# Two uses:\n# - bidsify a qMRLab dataset (SFalpha.mat / SF2alpha.mat plus Mask.mat) from\n# --mat-dir; this config is the --config input:\n#\n# qmrust bidsify --model b1_dam --mat-dir \\\n# --config recipes/non-bids/b1_dam_config.yaml --subject 01 --out ds\n#\n# - fit a stacked 4D NIfTI whose two volumes are [alpha, 2*alpha]:\n#\n# qmrust fit --data sfalpha_sf2alpha.nii.gz --mask Mask.mat \\\n# --config recipes/non-bids/b1_dam_config.yaml --output-dir out_b1\n", + "config_bids": "model: b1_dam\n\n# B1 = |acos(S(2a) / (2 * S(a)))| / a, with the map dimensionless (1.0 = the\n# nominal flip angle was achieved). BIDS input: the two flip angles are read\n# from the JSON sidecars via the model's protocol_schema(), so they are not\n# declared here. The TB1DAM collection is assembled by flip index. Example:\n#\n# qmrust fit --bids-dir \\\n# --config recipes/bids/b1_dam_config.yaml --output-dir out\n\n# BIDS mask selection. Left blank: this example dataset ships no mask, and\n# naming one here would claim otherwise. Blank means \"whatever mask is\n# present\", so a dataset of your own that has one is masked by it.\nmask:\n desc: \"\"\n", + "archive": "ds-tb1dam.zip", + "files": { + "data": "b1_dam.nii.gz", + "mask": null, + "aux": {} + }, + "probes": [ + { + "x": 32, + "y": 32, + "expected": { + "B1": 0.9143969256837542 + } + }, + { + "x": 29, + "y": 29, + "expected": { + "B1": 0.9100976261270164 + } + }, + { + "x": 29, + "y": 35, + "expected": { + "B1": 0.9081585700225492 + } + }, + { + "x": 35, + "y": 29, + "expected": { + "B1": 0.9401362108666856 + } + }, + { + "x": 35, + "y": 35, + "expected": { + "B1": 0.9576968485584598 + } + }, + { + "x": 32, + "y": 38, + "expected": { + "B1": 0.9165485297676551 + } + }, + { + "x": 38, + "y": 32, + "expected": { + "B1": 0.923119926577423 + } + }, + { + "x": 26, + "y": 26, + "expected": { + "B1": 0.9323702485199705 + } + } + ] +} \ No newline at end of file diff --git a/docs/playground/data/b1_dam.nii.gz b/docs/playground/data/b1_dam.nii.gz new file mode 100644 index 0000000..84472bd Binary files /dev/null and b/docs/playground/data/b1_dam.nii.gz differ diff --git a/docs/playground/data/index.json b/docs/playground/data/index.json index 3f472df..d3b3155 100644 --- a/docs/playground/data/index.json +++ b/docs/playground/data/index.json @@ -1,5 +1,7 @@ { "models": [ + "b1_afi", + "b1_dam", "inversion_recovery", "mono_t2", "mt_ratio", diff --git a/docs/playground/data/inversion_recovery.json b/docs/playground/data/inversion_recovery.json index 208bddd..5924908 100644 --- a/docs/playground/data/inversion_recovery.json +++ b/docs/playground/data/inversion_recovery.json @@ -1,6 +1,11 @@ { "model": "inversion_recovery", - "title": "Inversion recovery T1", + "title": "Inversion Recovery", + "bids_suffix": "IRT1", + "family": "T1 Relaxometry", + "family_icon": "spline-pointer", + "subgroup": null, + "category_order": 0, "dims": [ 128, 128, diff --git a/docs/playground/data/inversion_recovery.nii.gz b/docs/playground/data/inversion_recovery.nii.gz index 0d05d63..28097d8 100644 Binary files a/docs/playground/data/inversion_recovery.nii.gz and b/docs/playground/data/inversion_recovery.nii.gz differ diff --git a/docs/playground/data/inversion_recovery_mask.nii.gz b/docs/playground/data/inversion_recovery_mask.nii.gz index 13e1e9a..c3a50b5 100644 Binary files a/docs/playground/data/inversion_recovery_mask.nii.gz and b/docs/playground/data/inversion_recovery_mask.nii.gz differ diff --git a/docs/playground/data/mono_t2.json b/docs/playground/data/mono_t2.json index 3f08837..dcbc6c2 100644 --- a/docs/playground/data/mono_t2.json +++ b/docs/playground/data/mono_t2.json @@ -1,6 +1,11 @@ { "model": "mono_t2", - "title": "Mono-exponential T2", + "title": "Monoexp T2", + "bids_suffix": "MESE", + "family": "T2 Relaxometry", + "family_icon": "spline-pointer-mirrored", + "subgroup": null, + "category_order": 1, "dims": [ 87, 107, diff --git a/docs/playground/data/mono_t2.nii.gz b/docs/playground/data/mono_t2.nii.gz index df0d237..345951b 100644 Binary files a/docs/playground/data/mono_t2.nii.gz and b/docs/playground/data/mono_t2.nii.gz differ diff --git a/docs/playground/data/mono_t2_mask.nii.gz b/docs/playground/data/mono_t2_mask.nii.gz index 23e7894..a2d26b7 100644 Binary files a/docs/playground/data/mono_t2_mask.nii.gz and b/docs/playground/data/mono_t2_mask.nii.gz differ diff --git a/docs/playground/data/mt_ratio.json b/docs/playground/data/mt_ratio.json index 021a8e0..e98c3db 100644 --- a/docs/playground/data/mt_ratio.json +++ b/docs/playground/data/mt_ratio.json @@ -1,6 +1,11 @@ { "model": "mt_ratio", - "title": "Magnetization transfer ratio", + "title": "MT Ratio", + "bids_suffix": "MTR", + "family": "Magnetization Transfer", + "family_icon": "waves-arrow-up", + "subgroup": "Semi-quantitative MT", + "category_order": 3, "dims": [ 128, 135, diff --git a/docs/playground/data/mt_ratio.nii.gz b/docs/playground/data/mt_ratio.nii.gz index 010e5f8..173291f 100644 Binary files a/docs/playground/data/mt_ratio.nii.gz and b/docs/playground/data/mt_ratio.nii.gz differ diff --git a/docs/playground/data/mt_ratio_mask.nii.gz b/docs/playground/data/mt_ratio_mask.nii.gz index db5232e..7a20ecc 100644 Binary files a/docs/playground/data/mt_ratio_mask.nii.gz and b/docs/playground/data/mt_ratio_mask.nii.gz differ diff --git a/docs/playground/data/mt_sat.json b/docs/playground/data/mt_sat.json index 31a7d8f..8af5b6d 100644 --- a/docs/playground/data/mt_sat.json +++ b/docs/playground/data/mt_sat.json @@ -1,6 +1,11 @@ { "model": "mt_sat", - "title": "MT saturation", + "title": "MT Saturation", + "bids_suffix": "MTS", + "family": "Magnetization Transfer", + "family_icon": "waves-arrow-up", + "subgroup": "Semi-quantitative MT", + "category_order": 3, "dims": [ 128, 128, @@ -51,7 +56,7 @@ ], "enums": [], "config": "model: mt_sat\n\n# Non-BIDS / bidsify protocol source: per-weighting nominal flip angle (deg)\n# and repetition time (s). bidsify writes these into the MTS sidecars\n# (FlipAngle / RepetitionTimeExcitation); a non-BIDS fit reads them directly.\nmtw: { flip_angle: 6, repetition_time: 0.028 }\npdw: { flip_angle: 6, repetition_time: 0.028 }\nt1w: { flip_angle: 20, repetition_time: 0.018 }\n\n# Options (always used; not BIDS metadata).\n# B1 correction is applied only when a B1 map (TB1map) is supplied.\nb1_correction_factor: 0.4\nexport_mtr: true\n\n# Non-BIDS input is a named set: one .nii.gz per role under --nii-dir\n# (MTw.nii.gz, PDw.nii.gz, T1w.nii.gz); mask via --mask/--nii-mask. Example:\n#\n# qmrust bidsify --model mt_sat --nii-dir \\\n# --config recipes/non-bids/mt_sat_config.yaml --subject 01 --out ds\n", - "config_bids": "model: mt_sat\n\n# BIDS input: each MTS role's FlipAngle and RepetitionTimeExcitation come from\n# its JSON sidecar (via the model's protocol_schema), so the per-role\n# acquisition is NOT declared here.\n\n# Options (always used; not BIDS metadata).\n# B1 correction is applied only when a B1 map (TB1map) is present in the dataset.\nb1_correction_factor: 0.4\nexport_mtr: true\n\n# BIDS mask selection \u2014 disambiguates which mask to apply. If the dataset has\n# no matching mask, no masking is applied (as in this example).\nmask:\n desc: brain\n", + "config_bids": "model: mt_sat\n\n# BIDS input: each MTS role's FlipAngle and RepetitionTimeExcitation come from\n# its JSON sidecar (via the model's protocol_schema), so the per-role\n# acquisition is NOT declared here.\n\n# Options (always used; not BIDS metadata).\n# B1 correction is applied only when a B1 map (TB1map) is present in the dataset.\nb1_correction_factor: 0.4\nexport_mtr: true\n\n# BIDS mask selection. Left blank: this example dataset ships no mask, and\n# naming one here would claim otherwise. Blank means \"whatever mask is\n# present\", so a dataset of your own that has one is masked by it.\nmask:\n desc: \"\"\n", "archive": "ds-mts.zip", "files": { "data": "mt_sat.nii.gz", diff --git a/docs/playground/data/mt_sat.nii.gz b/docs/playground/data/mt_sat.nii.gz index 1dbe8ca..12f0a78 100644 Binary files a/docs/playground/data/mt_sat.nii.gz and b/docs/playground/data/mt_sat.nii.gz differ diff --git a/docs/playground/data/qmt_spgr.json b/docs/playground/data/qmt_spgr.json index d1d3572..a481083 100644 --- a/docs/playground/data/qmt_spgr.json +++ b/docs/playground/data/qmt_spgr.json @@ -1,6 +1,11 @@ { "model": "qmt_spgr", "title": "qMT-SPGR", + "bids_suffix": "QMTSPGR", + "family": "Magnetization Transfer", + "family_icon": "waves-arrow-up", + "subgroup": "Quantitative MT", + "category_order": 4, "dims": [ 88, 128, diff --git a/docs/playground/data/qmt_spgr.nii.gz b/docs/playground/data/qmt_spgr.nii.gz index 4330448..fe621f1 100644 Binary files a/docs/playground/data/qmt_spgr.nii.gz and b/docs/playground/data/qmt_spgr.nii.gz differ diff --git a/docs/playground/data/qmt_spgr_B0map.nii.gz b/docs/playground/data/qmt_spgr_B0map.nii.gz index 652523f..a7dc4e8 100644 Binary files a/docs/playground/data/qmt_spgr_B0map.nii.gz and b/docs/playground/data/qmt_spgr_B0map.nii.gz differ diff --git a/docs/playground/data/qmt_spgr_B1map.nii.gz b/docs/playground/data/qmt_spgr_B1map.nii.gz index f518fd1..434a9bc 100644 Binary files a/docs/playground/data/qmt_spgr_B1map.nii.gz and b/docs/playground/data/qmt_spgr_B1map.nii.gz differ diff --git a/docs/playground/data/qmt_spgr_R1map.nii.gz b/docs/playground/data/qmt_spgr_R1map.nii.gz index a452f8a..9b54b1f 100644 Binary files a/docs/playground/data/qmt_spgr_R1map.nii.gz and b/docs/playground/data/qmt_spgr_R1map.nii.gz differ diff --git a/docs/playground/data/qmt_spgr_mask.nii.gz b/docs/playground/data/qmt_spgr_mask.nii.gz index b9d9d9c..89a10bd 100644 Binary files a/docs/playground/data/qmt_spgr_mask.nii.gz and b/docs/playground/data/qmt_spgr_mask.nii.gz differ diff --git a/docs/playground/data/sources.json b/docs/playground/data/sources.json index f027ce7..eeb186e 100644 --- a/docs/playground/data/sources.json +++ b/docs/playground/data/sources.json @@ -3,8 +3,8 @@ "_url_form": "Must be Zenodo's /api/ content URL. The human-facing /records//files/?download=1 form sends no Access-Control-Allow-Origin and is blocked by the browser; the /api/ form sends 'ACAO: *'.", "_versioning": "'record' must be a concrete version id, not the concept id (10.5281/zenodo.21611551 / recid 21611551) — the concept id 404s on the /api/ files route. Adding an archive to Zenodo mints a NEW version id, so publishing a new dataset is the one case that does require editing 'record', 'doi' and 'base' here.", "host": "zenodo", - "record": "21696048", - "doi": "10.5281/zenodo.21696048", - "base": "https://zenodo.org/api/records/21696048/files", + "record": "21753826", + "doi": "10.5281/zenodo.21753826", + "base": "https://zenodo.org/api/records/21753826/files", "suffix": "/content" } diff --git a/docs/playground/data/vfa_t1.json b/docs/playground/data/vfa_t1.json index 82a8fa9..3976fb7 100644 --- a/docs/playground/data/vfa_t1.json +++ b/docs/playground/data/vfa_t1.json @@ -1,6 +1,11 @@ { "model": "vfa_t1", - "title": "Variable flip angle T1", + "title": "Variable Flip Angle", + "bids_suffix": "VFA", + "family": "T1 Relaxometry", + "family_icon": "spline-pointer", + "subgroup": null, + "category_order": 0, "dims": [ 128, 128, diff --git a/docs/playground/data/vfa_t1.nii.gz b/docs/playground/data/vfa_t1.nii.gz index be99aef..6475f55 100644 Binary files a/docs/playground/data/vfa_t1.nii.gz and b/docs/playground/data/vfa_t1.nii.gz differ diff --git a/docs/playground/data/vfa_t1_B1map.nii.gz b/docs/playground/data/vfa_t1_B1map.nii.gz index b66608e..878be0c 100644 Binary files a/docs/playground/data/vfa_t1_B1map.nii.gz and b/docs/playground/data/vfa_t1_B1map.nii.gz differ diff --git a/docs/playground/data/vfa_t1_mask.nii.gz b/docs/playground/data/vfa_t1_mask.nii.gz index 871955e..08cb278 100644 Binary files a/docs/playground/data/vfa_t1_mask.nii.gz and b/docs/playground/data/vfa_t1_mask.nii.gz differ diff --git a/docs/playground/index.html b/docs/playground/index.html index 627efcd..b10b794 100644 --- a/docs/playground/index.html +++ b/docs/playground/index.html @@ -62,7 +62,18 @@
Recipe
-
model:
+ +
+
Models
+ + +