diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 18a8bb6b79984407701c08b18b2c9c2a83dfbb96..6ee756a17c0c7a0d9f1288ba79937d92ee462bc0 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -1,13 +1,14 @@ 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) ------------------------------- diff --git a/fsl_mrs/core/mrs.py b/fsl_mrs/core/mrs.py index 0904c529e913907b4a7dd48f702e98bb6ff1f176..4f692b2f0c2d686a15cf98389d60edcc75d8bbe5 100644 --- a/fsl_mrs/core/mrs.py +++ b/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 diff --git a/fsl_mrs/scripts/basis_tools.py b/fsl_mrs/scripts/basis_tools.py index 76e436dc4f80ab790929a994550f2f08700f0d3b..5e73c101581ebc9507e6dc97c4754c36a5546d96 100644 --- a/fsl_mrs/scripts/basis_tools.py +++ b/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) diff --git a/fsl_mrs/tests/test_utils_misc.py b/fsl_mrs/tests/test_utils_misc.py index 374c7cb4e05ceb0d1e2c88a4aeb70dd14c3447d3..61395388ccfe60160eb1fab5a6dcd12c117a14f4 100644 --- a/fsl_mrs/tests/test_utils_misc.py +++ b/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)) diff --git a/fsl_mrs/utils/misc.py b/fsl_mrs/utils/misc.py index 2e1f03bf151ddc20568ce9f8392555a31a0fb5e1..448d876037f01663c7dcca5aac6ad606d6106df3 100644 --- a/fsl_mrs/utils/misc.py +++ b/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): """