Add remove peak functionality, tests needed.

fsl_mrs/scripts/basis_tools

@@ -135,6 +135,24 @@ def main():
                             help='Specify metabolite to conjugate')
     conjparser.set_defaults(func=conj)
 
+    # Remove peak
+    rempeakparser = sp.add_parser(
+        'remove_peak',
+        help='Edit a basis set by removing a peak (e.g. TMS reference).')
+    rempeakparser.add_argument('file', type=Path,
+                               help='Basis file')
+    rempeakparser.add_argument('output', type=Path,
+                               help='Output location, can overwrite.')
+    mutual = rempeakparser.add_mutually_exclusive_group(required=True)
+    mutual.add_argument('--metabolite', type=str, default=None,
+                        help='Specify metabolite to edit')
+    mutual.add_argument('--all', action="store_true",
+                        help='Edit all basis spectra.')
+    rempeakparser.add_argument('--ppmlim', default=(-.2, .2), type=float,
+                               nargs=2, metavar=('LOW', 'HIGH'),
+                               help='Peak removal ppm range (default=(-.2, .2))')
+    rempeakparser.set_defaults(func=rem_peak)
+
     # Parse command-line arguments
     args = p.parse_args()
 
@@ -262,5 +280,16 @@ def conj(args):
     ).save(args.output, overwrite=True)
 
 
+def rem_peak(args):
+    from fsl_mrs.utils import basis_tools
+    from fsl_mrs.utils.mrs_io import read_basis
+    basis_tools.remove_peak(
+        read_basis(args.file),
+        args.ppmlim,
+        name=args.metabolite,
+        all=args.all
+    ).save(args.output, overwrite=True)
+
+
 if __name__ == '__main__':
     main()

fsl_mrs/utils/basis_tools.py

@@ -14,6 +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
 
 
 class IncompatibleBasisError(Exception):
@@ -244,3 +245,26 @@ def difference_basis_sets(basis_1, basis_2, add_or_subtract='add', missing_metab
                             'fwhm': basis_2.basis_fwhm[index]})
 
     return Basis(np.asarray(difference_spec), names, headers)
+
+
+def remove_peak(basis, limits, name=None, all=False):
+
+    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')
+            basis.update_fid(corrected_obj, 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')
+        basis.update_fid(corrected_obj, name)
+
+    return basis