Port b1_dam and b1_afi, and group models by method family

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,14 +157,39 @@ 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" \
          "$DATA/out_mtr/qmrust/sub-01/anat/sub-01_MTRmap.nii.gz" \
          "$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

crates/qmrust-cli/src/bidsify.rs

@@ -24,7 +24,8 @@ pub struct BidsifyArgs {
     pub model: String,
     pub mat_data: Option<PathBuf>,
     pub mat_dir: Option<PathBuf>,
-    pub nii_data: Option<PathBuf>,
+    /// Either one 4D NIfTI, or one 3D NIfTI per volume in acquisition order.
+    pub nii_data: Vec<PathBuf>,
     pub nii_dir: Option<PathBuf>,
     pub nii_mask: Option<PathBuf>,
     pub mask: Option<PathBuf>,
@@ -54,15 +55,15 @@ 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!(
             args.mask.is_none(),
             "--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<Vec<(String, Arr
         .collect()
 }
 
-/// Read a 4D NIfTI measurement (+ optional NIfTI mask), preserving its spatial
-/// header. Auxiliary maps come from `--aux` for a NIfTI source; there is no
-/// directory convention to discover them from.
-fn read_nifti_source(nii: &Path, nii_mask: Option<&Path>, model: &dyn Model) -> Result<Source> {
-    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<Source> {
+    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,15 +425,16 @@ 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,
     suffix: &str,
     vol: &Array3<f64>,
     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<std::path::PathBuf> = (0..2)
+            .map(|i| {
+                let p = dir.join(format!("v{i}.nii"));
+                let data = ndarray::Array3::<f64>::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]