Commit ac774d86 authored by William Clarke's avatar William Clarke
Browse files

Add padd option to the basis_tools add script.

parent 9ed8e49e
......@@ -3,8 +3,9 @@ This document contains the FSL-MRS release history in reverse chronological orde
1.1.3 (TBC)
------------------------------
- Added mrs_tools script. Replaces mrs_vis and mrs_info. Adds split/merge/reorder functionality.
- Added basis_tools script. Tools for manipulating (shifting, scaling, converting, differencing, and adding to) basis sets.
- Added basis_tools script. Tools for manipulating (shifting, scaling, converting, differencing, conjugating, and adding to) basis sets.
- Improved display of basis sets using mrs_tools or basis_tools.
- Added 'default' MEGA-PRESS MM option to fsl_mrs and mrs class.
- Under the hood refactor of basis, MRS, and MRSI classes.
- Updated density matrix simulator. Added some automatic testing.
......
......@@ -72,6 +72,8 @@ def main():
help='Rescale amplitude to mean of target basis.')
addparser.add_argument('--conj', action="store_true",
help='Conjugate added basis.')
addparser.add_argument('--pad', action="store_true",
help='Pad added FID if required.')
addparser.set_defaults(func=add)
# Shift tool
......@@ -107,7 +109,7 @@ def main():
# Differenceing (add/subtract)
diffparser = sp.add_parser(
'diff',
help='Form difference basis sepctra from two basis spectra.')
help='Form difference basis spectra from two basis spectra.')
diffparser.add_argument('file1', type=Path,
help='Basis file one')
diffparser.add_argument('file2', type=Path,
......@@ -121,6 +123,18 @@ def main():
help='Ignore missing basis sets in one basis.')
diffparser.set_defaults(func=diff)
# Conjugate
conjparser = sp.add_parser(
'conj',
help='Conjugate (reverse frequency axis) all or single basis.')
conjparser.add_argument('file', type=Path,
help='Basis file')
conjparser.add_argument('output', type=Path,
help='Output location, can overwrite.')
conjparser.add_argument('--metabolite', type=str, default=None,
help='Specify metabolite to conjugate')
conjparser.set_defaults(func=conj)
# Parse command-line arguments
args = p.parse_args()
......@@ -190,6 +204,7 @@ def add(args):
scale=args.scale,
width=width,
conj=args.conj,
pad=args.pad,
sim_info=args.info)
......@@ -238,5 +253,14 @@ def diff(args):
new.save(args.output)
def conj(args):
from fsl_mrs.utils import basis_tools
from fsl_mrs.utils.mrs_io import read_basis
basis_tools.conjugate_basis(
read_basis(args.file),
name=args.metabolite
).save(args.output, overwrite=True)
if __name__ == '__main__':
main()
......@@ -114,3 +114,15 @@ def test_diff(tmp_path):
assert out_loc.is_dir()
assert (out_loc / 'NAA.json').is_file()
def test_conj(tmp_path):
out_loc = tmp_path / 'test_basis'
subprocess.check_call(['basis_tools', 'conj',
'--metabolite', 'NAA',
str(fsl),
str(out_loc)])
assert out_loc.is_dir()
assert (out_loc / 'NAA.json').is_file()
......@@ -49,6 +49,11 @@ def test_add_basis(tmp_path):
assert exc_info.type is basis_tools.IncompatibleBasisError
assert exc_info.value.args[0] == "The new basis FID covers too little time, try padding."
basis_tools.add_basis(fid, 'mac1', cf, bw, out_loc, pad=True)
new_basis = mrs_io.read_basis(out_loc)
index = new_basis.names.index('mac1')
assert 'mac1' in new_basis.names
fid_pad = np.pad(fid, (0, fid.size))
basis_tools.add_basis(fid_pad, 'mac2', cf, bw, out_loc)
new_basis = mrs_io.read_basis(out_loc)
......@@ -128,3 +133,13 @@ def test_add_sub():
new = basis_tools.difference_basis_sets(basis_1, basis_2, missing_metabs='raise')
assert exc_info.type is basis_tools.IncompatibleBasisError
assert exc_info.value.args[0] == "NAA does not occur in basis_1."
def test_conj():
basis = mrs_io.read_basis(basis_off)
basis_conj = basis_tools.conjugate_basis(mrs_io.read_basis(basis_off))
assert np.allclose(basis_conj.original_basis_array, basis.original_basis_array.conj())
basis_conj = basis_tools.conjugate_basis(mrs_io.read_basis(basis_off), name='NAA')
index = basis_conj.names.index('NAA')
assert np.allclose(basis_conj.original_basis_array[:, index], basis.original_basis_array[:, index].conj())
......@@ -45,7 +45,7 @@ def convert_lcm_basis(path_to_basis, output_location=None):
basis.save(output_location, info_str=sim_info)
def add_basis(fid, name, cf, bw, target, scale=False, width=None, conj=False, sim_info='Manually added'):
def add_basis(fid, name, cf, bw, target, scale=False, width=None, conj=False, pad=False, sim_info='Manually added'):
"""Add an additional basis spectrum to an existing FSL formatted basis set.
Optionally rescale the norm of the new FID to the mean of the existing ones.
......@@ -67,6 +67,8 @@ def add_basis(fid, name, cf, bw, target, scale=False, width=None, conj=False, si
:type width: float, optional
:param conj: Conjugate added FID, defaults to False.
:type conj: bool, optional
:param pad: Pad input FID to target length if required, defaults to False.
:type pad: bool, optional
:param sim_info: String added to the meta.SimVersion field, defaults to 'Manually added'
:type sim_info: str, optional
"""
......@@ -85,7 +87,19 @@ def add_basis(fid, name, cf, bw, target, scale=False, width=None, conj=False, si
target_dt,
target_basis.original_points)
except InsufficentTimeCoverageError:
raise IncompatibleBasisError('The new basis FID covers too little time, try padding.')
if not pad:
raise IncompatibleBasisError('The new basis FID covers too little time, try padding.')
else:
# Pad fid to sufficent length
required_time = target_basis.original_points * target_dt
fid_dt = 1 / bw
required_points = int(np.ceil(required_time / fid_dt))
fid = np.pad(fid, (0, required_points - fid.size), constant_values=complex(0.0))
resampled_fid = ts_to_ts(fid,
1 / bw,
target_dt,
target_basis.original_points)
# 3. Scale if requested
if scale:
......@@ -164,6 +178,27 @@ def rescale_basis(basis, name, target_scale=None):
return basis
def conjugate_basis(basis: Basis, name=None):
"""Conjugate all FIDs or just a selected FID in a basis set.
This reverses the frequency axis.
:param basis: Basis object containing FID(s) to conjugate
:type basis: fsl_mrs.core.basis.Basis
:param name: Metabolite name to conjugate, defaults to None which will conjugate all.
:type name: str, optional
:return: Modified basis object
:rtype: fsl_mrs.core.basis.Basis
"""
if name is not None\
and name in basis.names:
b = basis.original_basis_array[:, basis.names.index(name)]
basis.update_fid(b.conj(), name)
else:
for b, name in zip(basis.original_basis_array.T, basis.names):
basis.update_fid(b.conj(), name)
return basis
def difference_basis_sets(basis_1, basis_2, add_or_subtract='add', missing_metabs='raise'):
"""Add or subtract basis sets to form a set of difference spectra
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment