Merge branch 'enh/lcm_basis_conv' into 'master'

ENH: Update to peak removal in basis conversion.

See merge request fsl/fsl_mrs!47

CHANGELOG.rst

@@ -3,7 +3,8 @@ This document contains the FSL-MRS release history in reverse chronological orde
 1.1.11 (Monday 4th April 2022)
 ------------------------------
 - Now able to choose the number of workers in fsl_mrs_sim.
-- HSLVD routines now use dense (non-sparse) routine for better numerical stability
+- Basis conversion now can remove reference peaks in a single step.
+- Peak removal in basis set now defaults to zeroing rather than HLSVD for greater numerical stability. Mimics LCModel approach.
 - Updates and corrections to documentation, references to new FSL Course MRS section added.
 - Fixed bugs in LCModel basis set handling.
 - Removed divide by zero warnings in quantification of voxels where fitting has failed.

fsl_mrs/scripts/basis_tools

@@ -54,6 +54,10 @@ def main():
                                help='Input basis file or folder')
     convertparser.add_argument('output', type=Path,
                                help='Output fsl formatted folder, will be created if needed.')
+    convertparser.add_argument('--remove_reference', action="store_true",
+                               help='Remove LCModel reference peak.')
+    convertparser.add_argument('--hlsvd', action="store_true",
+                               help='Use HLSVD peak removal, rather than zeroing.')
     convertparser.set_defaults(func=convert)
 
     # Add tool - add a json formatted fid to a basis set
@@ -148,6 +152,8 @@ def main():
                         help='Specify metabolite to edit')
     mutual.add_argument('--all', action="store_true",
                         help='Edit all basis spectra.')
+    rempeakparser.add_argument('--hlsvd', action="store_true",
+                               help='Use HLSVD peak removal, rather than zeroing.')
     rempeakparser.add_argument('--ppmlim', default=(-.2, .2), type=float,
                                nargs=2, metavar=('LOW', 'HIGH'),
                                help='Peak removal ppm range (default=(-.2, .2))')
@@ -194,8 +200,20 @@ def convert(args):
     :type args: Namespace
     """
     from fsl_mrs.utils import basis_tools
+    from fsl_mrs.utils.mrs_io import read_basis
     basis_tools.convert_lcm_basis(args.input, args.output)
 
+    if args.remove_reference:
+        # TODO sort this conjugation mess out.
+        basis = read_basis(args.output)
+        basis = basis_tools.conjugate_basis(basis)
+        basis = basis_tools.remove_peak(
+            basis,
+            (-.2, .2),
+            all=True,
+            use_hlsvd=args.hlsvd)
+        basis_tools.conjugate_basis(basis).save(args.output, overwrite=True)
+
 
 def add(args):
     from fsl_mrs.utils.mrs_io import read_basis
@@ -298,7 +316,8 @@ def rem_peak(args):
         read_basis(args.file),
         args.ppmlim,
         name=args.metabolite,
-        all=args.all
+        all=args.all,
+        use_hlsvd=args.hlsvd
     ).save(args.output, overwrite=True)
 
 

fsl_mrs/tests/test_scripts_basis_tools.py

@@ -43,6 +43,16 @@ def test_convert(tmp_path):
     assert (tmp_path / 'new' / 'NAA.json').is_file()
 
 
+def test_convert_with_remove(tmp_path):
+    subprocess.check_call(['basis_tools', 'convert',
+                           '--remove_reference',
+                           str(lcm),
+                           str(tmp_path / 'new')])
+
+    assert (tmp_path / 'new').is_dir()
+    assert (tmp_path / 'new' / 'NAA.json').is_file()
+
+
 def test_add(tmp_path):
     out_loc = tmp_path / 'test_basis'
     copytree(fsl, out_loc)
@@ -126,3 +136,30 @@ def test_conj(tmp_path):
 
     assert out_loc.is_dir()
     assert (out_loc / 'NAA.json').is_file()
+
+
+def test_remove(tmp_path):
+    out_loc = tmp_path / 'test_basis'
+
+    subprocess.check_call(['basis_tools', 'remove_peak',
+                           '--all',
+                           '--ppmlim', '1.8', '2.0',
+                           str(fsl),
+                           str(out_loc)])
+
+    assert out_loc.is_dir()
+    assert (out_loc / 'NAA.json').is_file()
+
+
+def test_remove_hlsvd(tmp_path):
+    out_loc = tmp_path / 'test_basis'
+
+    subprocess.check_call(['basis_tools', 'remove_peak',
+                           '--all',
+                           '--ppmlim', '1.8', '2.0',
+                           '--hlsvd',
+                           str(fsl),
+                           str(out_loc)])
+
+    assert out_loc.is_dir()
+    assert (out_loc / 'NAA.json').is_file()

fsl_mrs/tests/test_utils_preproc.py

@@ -160,6 +160,36 @@ def test_apodize():
     assert apodFID[0] == 1.0
 
 
+# Test zero_spectrum by removing one peak from a two peak fid
+def test_remove():
+    # low noise
+    testFIDs, testHdrs = syn.syntheticFID(
+        noisecovariance=[[1E-6]],
+        amplitude=[1.0, 1.0],
+        chemicalshift=[-6, 0],
+        damping=[3, 3])
+    limits = [-7, -5]
+
+    testFIDs_single, _ = syn.syntheticFID(
+        noisecovariance=[[0]],
+        amplitude=[0.0, 1.0],
+        chemicalshift=[-6, 0],
+        damping=[3, 3])
+    # (FID,dwelltime,centralFrequency,limits,limitUnits = 'ppm',numSingularValues=50)
+
+    removedFID = preproc.zero_spectrum(
+        testFIDs[0],
+        testHdrs['dwelltime'],
+        testHdrs['centralFrequency'],
+        limits,
+        limitUnits='ppm')
+
+    assert np.allclose(
+        np.real(FIDToSpec(removedFID)),
+        np.real(FIDToSpec(testFIDs_single[0])),
+        atol=1E-2, rtol=1E-2)
+
+
 # Test hlsvd by removing one peak from a two peak fid
 def test_hlsvd():
     # low noise

fsl_mrs/utils/basis_tools.py

@@ -14,7 +14,7 @@ from fsl_mrs.core import MRS
 from fsl_mrs.core.basis import Basis
 from fsl_mrs.utils.qc import idPeaksCalcFWHM
 from fsl_mrs.utils.misc import ts_to_ts, InsufficentTimeCoverageError
-from fsl_mrs.utils.preproc import hlsvd
+from fsl_mrs.utils.preproc import hlsvd, zero_spectrum
 
 
 class IncompatibleBasisError(Exception):
@@ -246,26 +246,31 @@ def difference_basis_sets(basis_1, basis_2, add_or_subtract='add', missing_metab
     return Basis(np.asarray(difference_spec), names, headers)
 
 
-def remove_peak(basis, limits, name=None, all=False):
+def remove_peak(basis, limits, name=None, all=False, use_hlsvd=False):
+
+    def removal_func(fid):
+        if use_hlsvd:
+            return hlsvd(
+                fid,
+                basis.original_dwell,
+                basis.cf * 1E6,
+                limits,
+                limitUnits='ppm+shift',
+                numSingularValues=5)
+        else:
+            return zero_spectrum(
+                fid,
+                basis.original_dwell,
+                basis.cf * 1E6,
+                limits,
+                limitUnits='ppmshift')
 
     if all:
         for b, n in zip(basis.original_basis_array.T, basis.names):
-            corrected_obj = hlsvd(b,
-                                  basis.original_dwell,
-                                  basis.cf * 1E6,
-                                  limits,
-                                  limitUnits='ppm+shift',
-                                  numSingularValues=5)
-            basis.update_fid(corrected_obj, n)
+            basis.update_fid(removal_func(b), n)
     else:
         index = basis.names.index(name)
         indexed_fid = basis.original_basis_array[:, index]
-        corrected_obj = hlsvd(indexed_fid,
-                              basis.original_dwell,
-                              basis.cf * 1E6,
-                              limits,
-                              limitUnits='ppm+shift',
-                              numSingularValues=5)
-        basis.update_fid(corrected_obj, name)
+        basis.update_fid(removal_func(indexed_fid), name)
 
     return basis

fsl_mrs/utils/preproc/__init__.py

@@ -5,6 +5,6 @@ from fsl_mrs.utils.preproc.phasing import phaseCorrect, applyPhase
 from fsl_mrs.utils.preproc.eddycorrect import eddy_correct
 from fsl_mrs.utils.preproc.shifting import truncate, pad, timeshift, freqshift, shiftToRef
 from fsl_mrs.utils.preproc.filtering import apodize
-from fsl_mrs.utils.preproc.remove import hlsvd, model_fid_hlsvd
+from fsl_mrs.utils.preproc.remove import hlsvd, model_fid_hlsvd, zero_spectrum
 from fsl_mrs.utils.preproc.general import add, subtract
 from fsl_mrs.utils.preproc.unlike import identifyUnlikeFIDs

fsl_mrs/utils/preproc/remove.py

-# remove.py - HLSVD-based correction
+# remove.py - Removal of peaks, zeroing or HLSVD-based correction
 #
 # Author: Saad Jbabdi <saad@fmrib.ox.ac.uk>
 #         William Clarke <william.clarke@ndcn.ox.ac.uk>
@@ -8,10 +8,38 @@
 
 import numpy as np
 import hlsvdpropy
-from fsl_mrs.utils.misc import checkCFUnits
+from fsl_mrs.utils.misc import checkCFUnits, limit_to_range, calculateAxes, FIDToSpec, SpecToFID
 from fsl_mrs.utils.constants import H2O_PPM_TO_TMS
 
 
+def zero_spectrum(FID, dwelltime, centralFrequency, limits, limitUnits='ppmshift'):
+    """Zero part of a spectrum between limits
+
+    :param FID: FID to modify
+    :type FID: nump.array
+    :param dwelltime: Dwelltime in seconds
+    :type dwelltime: float
+    :param centralFrequency: Central imaging frequency in Hz or MHz
+    :type centralFrequency: float
+    :param limits: ppm limits
+    :type limits: tuple
+    :param limitUnits: Whether limits include shift ('ppm' or 'ppmshift), defaults to 'ppmshift'
+    :type limitUnits: str, optional
+    :return: Modified FID
+    :rtype: np.array
+    """
+    ppm_axes = calculateAxes(
+        1 / dwelltime,
+        checkCFUnits(centralFrequency),
+        len(FID),
+        H2O_PPM_TO_TMS)[limitUnits]
+    first, last = limit_to_range(ppm_axes, limits)
+    spec = FIDToSpec(FID)
+    spec[first:last] = 0.0 + 1j * 0.0
+    mod_fid = SpecToFID(spec)
+    return mod_fid
+
+
 def model_fid_hlsvd(FID, dwelltime, centralFrequency, limits=None,
                     limitUnits='ppm', numSingularValues=20):
     """Model a section of the FID using HLSVD. Optionally retain components
@@ -92,7 +120,7 @@ def _hlsvd(FID, dwelltime, centralFrequency, limits,
     :return: HLSVD modeled FID
     """
     m = FID.size // 2
-    r = hlsvdpropy.hlsvdpro(FID, numSingularValues, m=m, sparse=False)
+    r = hlsvdpropy.hlsvdpro(FID, numSingularValues, m=m, sparse=True)
     r = hlsvdpropy.convert_hlsvd_result(r, dwelltime)
     nsv_found, singular_values, frequencies, damping_factors, amplitudes, \
         phases = r[0:6]