New basis tools and tests.

fsl_mrs/scripts/basis_tools

+#!/usr/bin/env python
+
+""" basis_tools - top level script for calling general mrs basis set handling tools
+
+Author: William Clarke <william.clarke@ndcn.ox.ac.uk>
+
+Copyright (C) 2021 University of Oxford
+SHBASECOPYRIGHT
+"""
+
+import argparse
+from pathlib import Path
+
+from fsl_mrs import __version__
+
+
+def main():
+    # Parse command-line arguments
+    p = argparse.ArgumentParser(description="FSL Magnetic Resonance Spectroscopy - Basis Set Tools")
+
+    p.add_argument('-v', '--version', action='version', version=__version__)
+
+    sp = p.add_subparsers(title='subcommands',
+                          description='Availible tools',
+                          required=True,
+                          dest='subcommand')
+
+    # Info tool
+    infoparser = sp.add_parser(
+        'info',
+        help='Information about the basis set.')
+    infoparser.add_argument(
+        'file',
+        type=Path,
+        help='Basis file or folder')
+    infoparser.set_defaults(func=info)
+
+    # Vis tool
+    visparser = sp.add_parser(
+        'vis',
+        help='Quick visualisation of a basis set.')
+    visparser.add_argument('file', type=Path, metavar='FILE or DIR',
+                           help='NIfTI file or directory of basis sets')
+    visparser.add_argument('--ppmlim', default=(.2, 4.2), type=float,
+                           nargs=2, metavar=('LOW', 'HIGH'),
+                           help='limit the fit to a freq range (default=(.2,4.2))')
+    visparser.set_defaults(func=vis)
+
+    # Convert tool - Convert lcm or jmrui format to fsl format
+    convertparser = sp.add_parser(
+        'convert',
+        help='Convert LCModel or jMRUI formated basis to FSL format.')
+    convertparser.add_argument('input', type=Path,
+                               help='Input basis file or folder')
+    convertparser.add_argument('output', type=Path,
+                               help='Output fsl formatted folder, will be created if needed.')
+    convertparser.set_defaults(func=convert)
+
+    # Add tool - add a json formatted fid to a basis set
+    addparser = sp.add_parser(
+        'add',
+        help='Add a json formatted basis spectrum to a basis set.')
+    addparser.add_argument('new', type=Path,
+                           help='Basis to add')
+    addparser.add_argument('target', type=str,
+                           help='Basis to add to.')
+    addparser.add_argument('--info', type=str, default=None,
+                           help='Optional info string added to basis file.')
+    addparser.add_argument('--name', type=str, default=None,
+                           help='Optionally provide a different name for the basis.')
+    addparser.add_argument('--scale', action="store_true",
+                           help='Rescale amplitude to mean of target basis.')
+    addparser.set_defaults(func=add)
+
+    # Shift tool
+    shiftparser = sp.add_parser(
+        'shift',
+        help='Frequency shift a single basis spectrum.')
+    shiftparser.add_argument('file', type=Path,
+                             help='Basis file to modify')
+    shiftparser.add_argument('metabolite', type=str,
+                             help='Metabolite to shift')
+    shiftparser.add_argument('ppm_shift', type=float,
+                             help='Shift to apply (in ppm).')
+    shiftparser.add_argument('output', type=Path,
+                             help='Output location, can overwrite.')
+    shiftparser.set_defaults(func=shift)
+
+    # Rescale tool
+    scaleparser = sp.add_parser(
+        'scale',
+        help='Rescale a single basis spectrum to the mean of the other spectra.')
+    scaleparser.add_argument('file', type=Path,
+                             help='Basis file to modify')
+    scaleparser.add_argument('metabolite', type=str,
+                             help='Metabolite to shift')
+    scaleparser.add_argument('--target_scale', type=float,
+                             default=None,
+                             help='Use to specifiy an explicit scaling.')
+    scaleparser.add_argument('output', type=Path,
+                             help='Output location, can overwrite.')
+    scaleparser.set_defaults(func=rescale)
+
+    # Parse command-line arguments
+    args = p.parse_args()
+
+    # Call function
+    args.func(args)
+
+
+def info(args):
+    """Prints basic information about Basis sets
+    :param args: Argparse interpreted arguments
+    :type args: Namespace
+    """
+    from fsl_mrs.utils.mrs_io import read_basis
+    print(read_basis(args.file))
+
+
+def vis(args):
+    """Visualiser for Basis set
+    :param args: Argparse interpreted arguments
+    :type args: Namespace
+    """
+    from fsl_mrs.utils.mrs_io import read_basis
+    import matplotlib.pyplot as plt
+
+    # Some heuristics
+    if args.file.is_dir():
+        conj = True
+    else:
+        conj = False
+
+    basis = read_basis(args.file)
+    _ = basis.plot(ppmlim=args.ppmlim, conjugate=conj)
+    plt.show()
+
+
+def convert(args):
+    """Converter for lcm/jmrui basis sets
+    :param args: Argparse interpreted arguments
+    :type args: Namespace
+    """
+    from fsl_mrs.utils import basis_tools
+    basis_tools.convert_lcm_basis(args.input, args.output)
+
+
+def add(args):
+    from fsl_mrs.utils import basis_tools
+    import json
+    import numpy as np
+
+    with open(args.new, 'r') as fp:
+        json_dict = json.load(fp)
+
+    if 'basis' in json_dict:
+        json_dict = json_dict['basis']
+
+    fid = np.asarray(json_dict['basis_re']) + 1j * np.asarray(json_dict['basis_im'])
+    if args.name:
+        name = args.name
+    else:
+        name = json_dict['basis_name']
+    cf = json_dict['basis_centre']
+    bw = 1 / json_dict['basis_dwell']
+    width = json_dict['basis_width']
+
+    basis_tools.add_basis(fid, name, cf, bw, args.target, scale=args.scale, width=width, sim_info=args.info)
+
+
+def shift(args):
+    from fsl_mrs.utils import basis_tools
+    from fsl_mrs.utils.mrs_io import read_basis
+
+    basis = read_basis(args.file)
+    basis = basis_tools.shift_basis(basis, args.metabolite, args.ppm_shift)
+
+    basis.save(args.output, overwrite=True)
+
+
+def rescale(args):
+    from fsl_mrs.utils import basis_tools
+    from fsl_mrs.utils.mrs_io import read_basis
+
+    basis = read_basis(args.file)
+    basis = basis_tools.rescale_basis(basis, args.metabolite, args.target_scale)
+
+    basis.save(args.output, overwrite=True)
+
+
+if __name__ == '__main__':
+    main()

fsl_mrs/tests/test_scripts_basis_tools.py

+'''FSL-MRS test script
+
+Test the basis tools script
+
+Copyright Will Clarke, University of Oxford, 2021'''
+
+# Imports
+import subprocess
+from pathlib import Path
+from shutil import copytree
+# from unittest.mock import patch
+
+# Files
+testsPath = Path(__file__).parent
+
+jmrui = testsPath / 'testdata/mrs_io/basisset_JMRUI'
+lcm = testsPath / 'testdata/mrs_io/basisset_LCModel.BASIS'
+fsl = testsPath / 'testdata/mrs_io/basisset_FSL'
+mac = testsPath / 'testdata/basis_tools/macSeed.json'
+
+
+def test_info():
+    subprocess.check_call(['basis_tools', 'info', str(jmrui)])
+    subprocess.check_call(['basis_tools', 'info', str(lcm)])
+    subprocess.check_call(['basis_tools', 'info', str(fsl)])
+
+
+# To figure out one day
+# @patch("matplotlib.pyplot.show")
+# def test_vis(tmp_path):
+#     subprocess.check_call(['basis_tools', 'vis',
+#                            '--ppmlim', '0.2', '5.2',
+#                            str(jmrui)])
+
+def test_convert(tmp_path):
+    subprocess.check_call(['basis_tools', 'convert',
+                           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)
+
+    subprocess.check_call(['basis_tools', 'add',
+                           '--info', 'some info',
+                           '--scale',
+                           '--name', 'new_basis',
+                           str(mac),
+                           str(out_loc)])
+
+    assert (out_loc / 'new_basis.json').is_file()
+
+
+def test_shift(tmp_path):
+    out_loc = tmp_path / 'test_basis'
+
+    subprocess.check_call(['basis_tools', 'shift',
+                           str(fsl),
+                           'NAA',
+                           '1.0',
+                           str(out_loc)])
+
+    assert out_loc.is_dir()
+    assert (out_loc / 'NAA.json').is_file()
+
+
+def test_rescale(tmp_path):
+    out_loc = tmp_path / 'test_basis1'
+
+    subprocess.check_call(['basis_tools', 'scale',
+                           str(fsl),
+                           'NAA',
+                           str(out_loc)])
+    assert out_loc.is_dir()
+    assert (out_loc / 'NAA.json').is_file()
+
+    out_loc = tmp_path / 'test_basis2'
+    subprocess.check_call(['basis_tools', 'scale',
+                           str(fsl),
+                           'NAA',
+                           '--target_scale', '1.0',
+                           str(out_loc)])
+
+    assert out_loc.is_dir()
+    assert (out_loc / 'NAA.json').is_file()

fsl_mrs/tests/test_utils_basis_tools.py

@@ -73,3 +73,21 @@ def test_shift():
     t = basis.original_time_axis
     shifted_fid = basis.original_basis_array[:, index] * np.exp(-1j * 2 * np.pi * t * amount_in_hz)
     assert np.allclose(shifted.original_basis_array[:, index], shifted_fid)
+
+
+def test_rescale():
+    basis = mrs_io.read_basis(fsl_basis_path)
+
+    index = basis.names.index('Mac')
+    indexed_fid = basis.original_basis_array[:, index]
+    original_scale = np.linalg.norm(indexed_fid)
+
+    basis = basis_tools.rescale_basis(basis, 'Mac')
+    indexed_fid = basis.original_basis_array[:, index]
+    new_scale = np.linalg.norm(indexed_fid)
+    assert new_scale != original_scale
+
+    basis = basis_tools.rescale_basis(basis, 'Mac', target_scale=1.0)
+    indexed_fid = basis.original_basis_array[:, index]
+    new_scale = np.linalg.norm(indexed_fid)
+    assert np.isclose(new_scale, 1.0)

fsl_mrs/utils/basis_tools.py

@@ -62,6 +62,10 @@ def add_basis(fid, name, cf, bw, target, scale=False, width=None, sim_info='Manu
     :param scale: Rescale the fid so its norm is the mean of the norms of the
      other basis spectra, defaults to False
     :type scale: bool, optional
+    :param width: Width (fwhm in Hz) of added basis, defaults to None. If None it will be estimated.
+    :type width: float, optional
+    :param sim_info: String added to the meta.SimVersion field, defaults to 'Manually added'
+    :type sim_info: str, optional
     """
     # 1. Check that target exists
     target = Path(target)
@@ -97,7 +101,7 @@ def add_basis(fid, name, cf, bw, target, scale=False, width=None, sim_info='Manu
     target_basis.add_fid_to_basis(resampled_fid, name, width=width)
 
     # 6. Write to json without overwriting existing files
-    target_basis.save(target)
+    target_basis.save(target, info_str=sim_info)
 
 
 def shift_basis(basis, name, amount):
@@ -126,4 +130,26 @@ def shift_basis(basis, name, amount):
 
 
 def rescale_basis(basis, name, target_scale=None):
-    pass
\ No newline at end of file
+    """Rescale a single basis FID in place to either the mean of the set or a particular target
+
+    :param basis: Unmodified basis set
+    :type basis: fsl_mrs.core.basis.Basis
+    :param name: Metabolite name to SCALE
+    :type name: str
+    :param target_scale: Norm of slected basis will be scaled to this target, defaults to None = mean of other FIDs
+    :type target_scale: float, optional
+    :return: Basis object with new scaling
+    :rtype: fsl_mrs.core.basis.Basis
+    """
+    index = basis.names.index(name)
+    indexed_fid = basis.original_basis_array[:, index]
+    all_except_index = np.delete(basis.original_basis_array, index, axis=1)
+
+    if target_scale:
+        indexed_fid *= target_scale / np.linalg.norm(indexed_fid)
+    else:
+        indexed_fid *= np.linalg.norm(np.mean(all_except_index, axis=1), axis=0) / np.linalg.norm(indexed_fid)
+
+    basis.update_fid(indexed_fid, name)
+
+    return basis

setup.py

@@ -42,6 +42,7 @@ setup(name='fsl_mrs',
                'fsl_mrs/scripts/fsl_mrs_proc',
                'fsl_mrs/scripts/fsl_mrs_sim',
                'fsl_mrs/scripts/mrs_tools',
+               'fsl_mrs/scripts/basis_tools',
                'fsl_mrs/scripts/merge_mrs_reports',
                'fsl_mrs/scripts/svs_segment',
                'fsl_mrs/scripts/mrsi_segment',