Merge branch 'bf/basis_shift_conj' into 'master'

BF: basis_tools all_shift conjugation

See merge request !129

CHANGELOG.rst

 This document contains the FSL-MRS release history in reverse chronological order.
 
-2.4.2 (WIP)
--------------------------------
+2.4.2 (Monday 17th March 2025)
+------------------------------
 - Added a `freeshift_lorentzian` model with inter-metabolite shifts (freeshift) and lorentzian lineshape.
 - Added `freeshift` option to mrsi `fsl_mrsi` and dynamic `fsl_dynmrs` fitting scripts.
 - Added option to show phase corrected spectum in the `fsl_mrs` HTML fitting report.
 - Added option to suppress alignment step in `fsl_mrs_preproc{_edit}`.
 - Changed syntax for water removal in `fsl_mrs_preproc{_edit}` from `--hlsvd` to `--remove-water`
 - Changed syntax for truncation in `fsl_mrs_preproc{_edit}` from `--leftshift {N}` to `--truncate-fid {N}`
+- Basis sets are automatically conjugated when using `basis_tools shift_all`.
 
 2.4.1 (Saturday 1st March 2025)
 -------------------------------

fsl_mrs/core/mrs.py

@@ -587,12 +587,10 @@ class MRS(object):
         if ppmlim is None:
             ppmlim = self.default_ppm_range
 
-        first, last = self.ppmlim_to_range(ppmlim)
-        Spec1 = np.real(misc.FIDToSpec(self.FID))[first:last]
-        Spec2 = np.real(misc.FIDToSpec(np.conj(self.FID)))[first:last]
-
-        if np.linalg.norm(misc.detrend(Spec1, deg=4)) < \
-           np.linalg.norm(misc.detrend(Spec2, deg=4)):
+        if misc.detect_conjugation(
+                self.FID,
+                self.ppmAxisShift,
+                ppmlim):
             if repair is False:
                 warnings.warn('YOU MAY NEED TO CONJUGATE YOUR FID!!!')
                 return -1
@@ -620,27 +618,16 @@ class MRS(object):
         if ppmlim is None:
             ppmlim = self.default_ppm_range
 
-        first, last = self.ppmlim_to_range(ppmlim)
-
-        conjOrNot = []
-        basis = self.basis
-        for b in basis.T:
-            Spec1 = np.real(misc.FIDToSpec(b))[first:last]
-            Spec2 = np.real(misc.FIDToSpec(np.conj(b)))[first:last]
-            if np.linalg.norm(misc.detrend(Spec1, deg=4)) < \
-               np.linalg.norm(misc.detrend(Spec2, deg=4)):
-                conjOrNot.append(1.0)
-            else:
-                conjOrNot.append(0.0)
-
-        if (sum(conjOrNot) / len(conjOrNot)) > 0.5:
+        if misc.detect_conjugation(
+                self.basis.T,
+                self.ppmAxisShift,
+                ppmlim):
             if repair is False:
                 warnings.warn('YOU MAY NEED TO CONJUGATE YOUR BASIS!!!')
                 return -1
             else:
                 self.conj_Basis = True
                 return 1
-
         return 0
 
     # Plotting functions

fsl_mrs/scripts/basis_tools.py

@@ -345,14 +345,42 @@ def shift(args):
 def all_shift(args):
     from fsl_mrs.utils import basis_tools
     from fsl_mrs.utils.mrs_io import read_basis
+    from fsl_mrs.utils.misc import detect_conjugation
     import pandas as pd
     import numpy as np
+    import json
+    from fsl_mrs.utils.constants import PPM_RANGE
 
     basis = read_basis(args.file)
+
+    # Extract settings from the saved fitting arguments
+    with open(args.results_dir / 'options.txt', "r") as f:
+        fit_args = json.loads(f.readline())
+    ppmlims = fit_args['ppmlim'] if fit_args['ppmlim'] else PPM_RANGE['1H']
+
+    # Load fitting results
     all_results = pd.read_csv(args.results_dir / 'all_parameters.csv', index_col=0)
     shift_res = all_results.filter(regex='eps', axis=0)['mean'] / (2 * np.pi * basis.cf)
-    for name, shift in zip(basis.names, shift_res.to_list()):
-        basis = basis_tools.shift_basis(basis, name, shift)
+
+    # Apply shifts to each basis
+    # Currently this assumes that all basis spectra are used.
+    assert len(basis.names) == len(shift_res.to_list())
+
+    def apply_shifts(x):
+        for name, shift in zip(x.names, shift_res.to_list()):
+            x = basis_tools.shift_basis(x, name, shift)
+        return x
+
+    if detect_conjugation(
+            basis.original_basis_array.T,
+            basis.original_ppm_shift_axis,
+            ppmlims):
+        from fsl_mrs.utils.basis_tools import conjugate_basis
+        basis = conjugate_basis(basis)
+        basis = apply_shifts(basis)
+        basis = conjugate_basis(basis)
+    else:
+        basis = apply_shifts(basis)
 
     basis.save(args.output, overwrite=True)
 

fsl_mrs/tests/test_utils_misc.py

@@ -222,3 +222,31 @@ def test_create_peak():
         gamma=10)
 
     assert np.allclose(one_ppm + zero_ppm, both_peaks)
+
+
+def test_detect_conjugation():
+    FIDs, headers = synth.syntheticFID(
+        chemicalshift=[-2, -3])
+    FID = FIDs[0]
+    FIDs = np.stack((FIDs[0], FIDs[0], FIDs[0].conj()))
+
+    assert not misc.detect_conjugation(
+        FID,
+        headers['ppmaxis'],
+        (0, -4))
+    assert misc.detect_conjugation(
+        FID.conj(),
+        headers['ppmaxis'],
+        (0, -4))
+
+    assert misc.detect_conjugation(
+        FIDs.conj(),
+        headers['ppmaxis'],
+        (0, -4))
+
+    # Test wrong orientation of data raises error
+    with pytest.raises(ValueError):
+        misc.detect_conjugation(
+            FIDs.T,
+            headers['ppmaxis'],
+            (0, -4))

fsl_mrs/utils/misc.py

@@ -792,6 +792,44 @@ def parse_metab_groups(mrs, metab_groups):
     return out
 
 
+def detect_conjugation(
+        data: np.ndarray,
+        ppmaxis: np.ndarray,
+        ppmlim: tuple) -> bool:
+    """Detect whether data should be conjugated based on
+    the amount of information content in ppm range.
+
+    :param data: FID or stack of FIDS (last dimension is time).
+    :type data: np.ndarray
+    :param ppmaxis: ppmaxis to match FFT of FID
+    :type ppmaxis: np.ndarray
+    :param ppmlim: Limits that define region of expected signal
+    :type ppmlim: tuple
+    :return: True if FIDs should be conjugated to maximise signal in limits.
+    :rtype: bool
+    """
+    if data.shape[-1] != ppmaxis.size:
+        raise ValueError("data's last dimension must matcht he size of ppmaxis.")
+
+    first, last = limit_to_range(ppmaxis, ppmlim)
+
+    def conj_or_not(x):
+        Spec1 = np.real(FIDToSpec(x))[first:last]
+        Spec2 = np.real(FIDToSpec(np.conj(x)))[first:last]
+        if np.linalg.norm(detrend(Spec1, deg=4)) < \
+                np.linalg.norm(detrend(Spec2, deg=4)):
+            return 1
+        else:
+            return 0
+
+    if data.ndim == 2:
+        return sum([conj_or_not(fid) for fid in data]) >= 0.5
+    elif data.ndim == 1:
+        return bool(conj_or_not(data))
+    else:
+        raise ValueError(f'data must have one or two dimensions, not {data.ndim}.')
+
+
 # ----- MRSI stuff ---- #
 def volume_to_list(data, mask):
     """