Reordering EPI based on reco when supplied

brukerapi/dataset.py

@@ -4,6 +4,7 @@
 import os.path
 import re
 import warnings
+from contextlib import suppress
 from copy import deepcopy
 from pathlib import Path
 
@@ -445,7 +446,10 @@ def add_parameter_file(self, file):
             dataset["PVM_DwDir"].value
 
         """
-        path = traverse(self.path.parent, RELATIVE_PATHS[self.type][file])
+        if file in RELATIVE_PATHS[self.type]:
+            path = traverse(self.path.parent, RELATIVE_PATHS[self.type][file])
+        else:
+            path = as_path(file)
 
         if not hasattr(self, "_parameters") or self._parameters is None:
             self._parameters = {path.name: JCAMPDX(path)}
@@ -474,6 +478,16 @@ def _read_parameters(self):
                 # if jcampdx file is not found, but not required, pass
                 pass
 
+        # The vendor's first reconstruction is the conventional default. A
+        # fid can have several reconstructions, however, so callers can select
+        # the one whose input-order declaration applies with ``reco_path=``.
+        reco_path = self._state.get("reco_path")
+        if reco_path is not None:
+            self.add_parameter_file(reco_path)
+        elif self.type in {"fid", "fid_proc"}:
+            with suppress(FileNotFoundError):
+                self.add_parameter_file("reco")
+
     def _write_parameters(self, parent):
         for type_, jcampdx in self._parameters.items():
             jcampdx.write(parent / type_)

brukerapi/schemas.py

@@ -162,6 +162,45 @@ def acquisition_factor(self):
         """Number of stored scalar samples per logical sample."""
         return 1 if str(self._dataset._parameter_value("AQ_mod", "")).lower() == "qf" else 2
 
+    def _declared_or_inferred(self, parameter, declared_value, inferred, warning_flag):
+        """Prefer a loaded format declaration, retaining scheme inference."""
+        if declared_value is None:
+            return inferred
+        if declared_value != inferred and not getattr(self._dataset, warning_flag, False):
+            warnings.warn(
+                f"{parameter}={self._dataset._parameter_value(parameter)!r} disagrees with "
+                f"scheme_id={self._dataset.scheme_id!r} for {self._dataset.path}; "
+                f"using the declared value",
+                RuntimeWarning,
+                stacklevel=3,
+            )
+            setattr(self._dataset, warning_flag, True)
+        return declared_value
+
+    @property
+    def continuous_train(self):
+        """Whether §5.3 declares a phase-factor continuous acquisition train."""
+        value = self._dataset._parameter_value("ACQ_scan_size")
+        declared = None if value is None else str(value).strip("<>").casefold() == "acq_phase_factor_scans"
+        return self._declared_or_inferred(
+            "ACQ_scan_size",
+            declared,
+            "EPI" in self._dataset.scheme_id,
+            "_warned_scan_size_disagreement",
+        )
+
+    @property
+    def mirror_odd_lines(self):
+        """Whether a selected reconstruction requests §6.2 odd-line reversal."""
+        value = self._dataset._parameter_value("RECO_inp_order")
+        declared = None if value is None else str(value).strip("<>").upper() == "REV_ALT_ROWS"
+        return self._declared_or_inferred(
+            "RECO_inp_order",
+            declared,
+            "EPI" in self._dataset.scheme_id,
+            "_warned_reco_input_order_disagreement",
+        )
+
     @property
     def layouts(self):
         """Dictionary of possible logical layouts of data
@@ -193,7 +232,7 @@ def layouts(self):
 
         layouts["encoding_permuted"] = tuple(np.array(layouts["encoding_space"])[np.array(layouts["permute"])])
 
-        if "EPI" in self._dataset.scheme_id:
+        if self.continuous_train:
             discarded = self._dataset.block_size - self._dataset.acq_length
             block_format = self._dataset._parameter_value("GO_block_size")
             scan_shift = self._dataset._parameter_value("ACQ_scan_shift", 0)
@@ -221,7 +260,7 @@ def deserialize(self, data, layouts):
         # Typically for RARE, or EPI
         data = self._reorder_fid_lines(data, dir="FW")
 
-        if "EPI" in self._dataset.scheme_id:
+        if self.mirror_odd_lines:
             data = self._mirror_odd_lines(data)
 
         data = self._reorder_objects(data, dir="FW")
@@ -385,7 +424,7 @@ def _mirror_odd_lines(self, data):
     def serialize(self, data, layouts):
         data = self._reorder_objects(data, dir="BW")
 
-        if "EPI" in self._dataset.scheme_id:
+        if self.mirror_odd_lines:
             data = self._mirror_odd_lines(data)
 
         data = self._reorder_fid_lines(data, dir="BW")

test/test_dataset.py

@@ -691,7 +691,8 @@ def test_monotonic_phase_encode_steps_skip_reordering():
 
 @pytest.mark.skipif(not PV51_STUDY_PATH.is_dir(), reason="PV51 test data is not available")
 def test_epi_standard_kblock_trims_trailing_padding():
-    dataset = Dataset(PV51_STUDY_PATH / "13" / "fid")
+    with pytest.warns(RuntimeWarning, match="RECO_inp_order=.*disagrees with scheme_id"):
+        dataset = Dataset(PV51_STUDY_PATH / "13" / "fid")
     dataset["GO_block_size"].val_str = "Standard_KBlock_Format"
     dataset.block_size = 256
     dataset.acq_length = 200
@@ -706,7 +707,8 @@ def test_epi_standard_kblock_trims_trailing_padding():
 
 @pytest.mark.skipif(not PV51_STUDY_PATH.is_dir(), reason="PV51 test data is not available")
 def test_epi_standard_kblock_warns_on_nonzero_discarded_samples():
-    dataset = Dataset(PV51_STUDY_PATH / "13" / "fid")
+    with pytest.warns(RuntimeWarning, match="RECO_inp_order=.*disagrees with scheme_id"):
+        dataset = Dataset(PV51_STUDY_PATH / "13" / "fid")
     dataset["GO_block_size"].val_str = "Standard_KBlock_Format"
     dataset.block_size = 8
     dataset.acq_length = 6

test/test_raw_layouts.py

@@ -5,9 +5,10 @@
 """
 
 import numpy as np
+import pytest
 
 from brukerapi.dataset import LOAD_STAGES, Dataset
-from test.synthetic import Verbatim, write_binary, write_fid
+from test.synthetic import Verbatim, write_binary, write_fid, write_jcampdx
 
 # pv6 EPSI, shrunk: 6 read points, 4 phase steps, 2 segments, 4 spectral points
 # per read line -- the real scan is 96 x 64 x 4 x 64.
@@ -84,6 +85,48 @@ def test_epsi_keeps_every_stored_complex_sample(tmp_path):
     assert np.array_equal(np.sort_complex(dataset.data.reshape(-1)), np.sort_complex(expected))
 
 
+def test_fid_uses_pdata_one_reco_unless_a_different_one_is_selected(tmp_path):
+    """`pdata/1` is conventional; `reco_path` selects another reconstruction."""
+    path = write_fid(
+        tmp_path / "30",
+        acqp={
+            "GO_block_size": "continuous",
+            "GO_raw_data_format": "GO_32BIT_SGN_INT",
+            "BYTORDA": "little",
+            "ACQ_dim": 2,
+            "ACQ_dim_desc": Verbatim("( 2 )\nSpatial Spatial"),
+            "ACQ_size": np.array([8, 2]),
+            "ACQ_phase_factor": 1,
+            "ACQ_scan_size": "ACQ_phase_factor_scans",
+            "PULPROG": "<FLASH.ppg>",
+            "NI": 1,
+            "NR": 1,
+        },
+        method={"Method": "<Bruker:FLASH>", "PVM_DigNp": 4, "PVM_EncMatrix": np.array([4, 2]), "PVM_EncNReceivers": 1},
+        blocks=2,
+    )
+    write_jcampdx(path.parent / "pdata" / "1" / "reco", {"RECO_inp_order": "NO_REORDERING"})
+    reco = write_jcampdx(path.parent / "pdata" / "2" / "reco", {"RECO_inp_order": "REV_ALT_ROWS"})
+
+    default = Dataset(path, load=LOAD_STAGES["properties"])
+    default.load_schema()
+    assert "reco" in default.parameters
+    assert default["RECO_inp_order"].value == "NO_REORDERING"
+    assert not default._schema.mirror_odd_lines
+
+    dataset = Dataset(path, reco_path=reco, load=LOAD_STAGES["properties"])
+    dataset.load_schema()
+    assert "reco" in dataset.parameters
+    with pytest.warns(RuntimeWarning, match="ACQ_scan_size=.*disagrees with scheme_id"):
+        assert dataset._schema.continuous_train
+    with pytest.warns(RuntimeWarning, match="RECO_inp_order=.*disagrees with scheme_id"):
+        assert dataset._schema.mirror_odd_lines
+
+    lines = np.arange(8).reshape((4, 2), order="F")
+    expected = np.column_stack((lines[:, 0], lines[::-1, 1]))
+    assert np.array_equal(dataset._schema._mirror_odd_lines(lines), expected)
+
+
 def kblock_fid(tmp_path, acqp, method, *, blocks, samples_per_block):
     """A Standard_KBlock_Format fid: real samples first, then zero padding."""
     block_size = int(np.ceil(int(np.atleast_1d(acqp["ACQ_size"])[0]) * int(method["PVM_EncNReceivers"]) * 4 / 1024.0) * 1024 // 4)