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

New hlsvd modelling command added for dwmrs ecc at high b values.

parent 29634598
......@@ -3,6 +3,7 @@ This document contains the FSL-MRS release history in reverse chronological orde
1.1.5 (WIP)
-------------------------------
- fsl_mrs_proc align can now align across all higher dimension FIDs. Pass 'all' as dimension tag.
- New command "fsl_mrs_proc model". HSLVD modelling of peaks in defined region. Number of components settable.
1.1.4 (Tuesday 3rd August 2021)
-------------------------------
......
......@@ -144,6 +144,22 @@ def main():
hlsvdparser.set_defaults(func=remove)
add_common_args(hlsvdparser)
# model subcommand - model peaks using HLSVD
modelparser = sp.add_parser('model', add_help=False,
help='Model peaks with HLSVD.')
model_group = modelparser.add_argument_group('HLSVD modelling arguments')
model_group.add_argument('--file', type=str, required=True,
help='Data file(s)')
model_group.add_argument('--ppm', type=float, nargs=2,
metavar='<lower-limit upper-limit>',
default=[4.5, 4.8],
help='ppm limits of removal window')
model_group.add_argument('--components', type=int,
default=5,
help='Number of components to model.')
modelparser.set_defaults(func=model)
add_common_args(modelparser)
# tshift subcommand - shift/resample in timedomain
tshiftparser = sp.add_parser('tshift', add_help=False,
help='shift/resample in timedomain.')
......@@ -524,6 +540,16 @@ def remove(dataobj, args):
return datacontainer(corrected, dataobj.datafilename)
def model(dataobj, args):
modelled = preproc.hlsvd_model_peaks(dataobj.data,
limits=args['ppm'],
components=args['components'],
report=args['generateReports'],
report_all=args['allreports'])
return datacontainer(modelled, dataobj.datafilename)
def tshift(dataobj, args):
shifted = preproc.tshift(dataobj.data,
tshiftStart=args['tshiftStart'],
......
......@@ -459,6 +459,44 @@ def test_remove(svs_data, mrsi_data, tmp_path):
assert np.allclose(data.data, directRun.data)
def test_model(svs_data, mrsi_data, tmp_path):
svsfile, mrsifile, svsdata, mrsidata = splitdata(svs_data, mrsi_data)
# Run model on both sets of data using the command line
subprocess.check_call(['fsl_mrs_proc',
'model',
'--file', svsfile,
'--ppm', '-10', '10',
'--components', '5',
'--output', tmp_path,
'--filename', 'tmp'])
# Load result for comparison
data = read_FID(op.join(tmp_path, 'tmp.nii.gz'))
# Run directly
directRun = preproc.hlsvd_model_peaks(svsdata, (-10, 10), components=5)
assert np.allclose(data.data, directRun.data)
# Run model on both sets of data using the command line
subprocess.check_call(['fsl_mrs_proc',
'model',
'--file', mrsifile,
'--ppm', '-10', '10',
'--components', '5',
'--output', tmp_path,
'--filename', 'tmp'])
# Load result for comparison
data = read_FID(op.join(tmp_path, 'tmp.nii.gz'))
# Run directly
directRun = preproc.hlsvd_model_peaks(mrsidata, (-10, 10), components=5)
assert np.allclose(data.data, directRun.data)
def test_align_diff(svs_data_diff, mrsi_data_diff, tmp_path):
svsfile, svsdata = svs_data_diff[0], svs_data_diff[1]
......
......@@ -118,6 +118,18 @@ def test_remove():
== 'fsl_mrs.utils.preproc.nifti_mrs_proc.remove_peaks, limits=(4, 5.3), limit_units=ppm+shift.'
def test_hlsvd_model():
nmrs_obj = read_FID(wrefc)
nmrs_obj = nproc.average(nmrs_obj, 'DIM_DYN')
modeled = nproc.hlsvd_model_peaks(nmrs_obj, (4, 5.30), components=3)
assert modeled.hdr_ext['ProcessingApplied'][1]['Method'] == 'HLSVD modeling'
assert modeled.hdr_ext['ProcessingApplied'][1]['Details']\
== 'fsl_mrs.utils.preproc.nifti_mrs_proc.hlsvd_model_peaks,'\
' limits=(4, 5.3), limit_units=ppm+shift, components=3.'
def test_tshift():
nmrs_obj = read_FID(wrefc)
nmrs_obj = nproc.average(nmrs_obj, 'DIM_DYN')
......
......@@ -179,6 +179,24 @@ def test_hlsvd():
assert np.allclose(np.real(removedFID), np.real(onResFID), atol=1E-2, rtol=1E-2)
# Test hlsvd modelling by 'denoising'
def test_model_fid_hlsvd():
# low noise
testFIDs, testHdrs = syn.syntheticFID(noisecovariance=[[1E-4]], amplitude=[1.0], chemicalshift=[0])
noislessFIDs, _ = syn.syntheticFID(noisecovariance=[[0]], amplitude=[1.0], chemicalshift=[0])
limits = [-1.5, 1.5]
modelledFID = preproc.model_fid_hlsvd(testFIDs[0],
testHdrs['dwelltime'],
testHdrs['centralFrequency'],
limits,
limitUnits='ppm',
numSingularValues=5)
assert np.allclose(np.real(modelledFID), np.real(noislessFIDs), atol=1E-2, rtol=1E-2)
def test_pad_truncate():
mockFID = np.full((1024), 1.0)
testFID = preproc.pad(mockFID, 1, first_or_last='last')
......
......@@ -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
from fsl_mrs.utils.preproc.remove import hlsvd, model_fid_hlsvd
from fsl_mrs.utils.preproc.general import add, subtract
from fsl_mrs.utils.preproc.unlike import identifyUnlikeFIDs
......@@ -419,6 +419,54 @@ def remove_peaks(data, limits, limit_units='ppm+shift', figure=False, report=Non
return corrected_obj
def hlsvd_model_peaks(data, limits,
limit_units='ppm+shift', components=5, figure=False, report=None, report_all=False):
'''Apply HLSVD to model spectum
:param NIFTI_MRS data: Data to model
:param limits: ppm limits between which spectrum will be modeled
:param str limit_units: Can be 'Hz', 'ppm' or 'ppm+shift'.
:param int components: Number of lorentzian components to model
:param figure: True to show figure.
:param report: Provide output location as path to generate report
:param report_all: True to output all indicies
:return: Corrected data in NIFTI_MRS format.
'''
corrected_obj = data.copy()
for dd, idx in data.iterate_over_dims(iterate_over_space=True):
corrected_obj[idx] = preproc.model_fid_hlsvd(
dd,
data.dwelltime,
data.spectrometer_frequency[0],
limits,
limitUnits=limit_units,
numSingularValues=components)
if (figure or report) and (report_all or first_index(idx)):
from fsl_mrs.utils.preproc.remove import hlsvd_report
fig = hlsvd_report(dd,
corrected_obj[idx],
limits,
data.bandwidth,
data.spectrometer_frequency[0],
nucleus=data.nucleus[0],
limitUnits=limit_units,
html=report)
if figure:
fig.show()
# Update processing prov
processing_info = f'{__name__}.hlsvd_model_peaks, '
processing_info += f'limits={limits}, '
processing_info += f'limit_units={limit_units}, '
processing_info += f'components={components}.'
update_processing_prov(corrected_obj, 'HLSVD modeling', processing_info)
return corrected_obj
def tshift(data, tshiftStart=0.0, tshiftEnd=0.0, samples=None, figure=False, report=None, report_all=False):
'''Apply time shift or resampling to each FID
:param NIFTI_MRS data: Data to shift
......
#!/usr/bin/env python
# remove.py - HLSVD-based correction
#
# Author: Saad Jbabdi <saad@fmrib.ox.ac.uk>
......@@ -14,6 +12,36 @@ from fsl_mrs.utils.misc import checkCFUnits
from fsl_mrs.utils.constants import H2O_PPM_TO_TMS
def model_fid_hlsvd(FID, dwelltime, centralFrequency, limits=None,
limitUnits='ppm', numSingularValues=20):
"""Model a section of the FID using HLSVD. Optionally retain components
only within the frequenccy/ppm limits.
:param FID: Time domain data
:type FID: numpy.array
:param dwelltime: dwell time in seconds
:type dwelltime: float
:param centralFrequency: Central frequency in Hz
:type centralFrequency: float
:param limits: Frequency/ppm limits, defaults to None
:type limits: tuple, optional
:param limitUnits: Axis that limits are given in. By Default
in ppm, relative to receiver frequency (no shift). Can be 'Hz', 'ppm'
or 'ppm+shift'. Defaults to 'ppm'
:type limitUnits: str, optional
:param numSingularValues: Max number of singular values, defaults to 20
:type numSingularValues: int, optional
"""
return _hlsvd(
FID,
dwelltime,
centralFrequency,
limits,
limitUnits=limitUnits,
numSingularValues=numSingularValues)
def hlsvd(FID, dwelltime, centralFrequency, limits,
limitUnits='ppm', numSingularValues=20):
""" Run HLSVDPRO on FID
......@@ -31,9 +59,38 @@ def hlsvd(FID, dwelltime, centralFrequency, limits,
Returns:
FID (ndarray): Modified FID
"""
# nsv_found, singular_values, frequencies, damping_factors, amplitudes, \
# phases = hlsvdpropy.hlsvd(FID, numSingularValues, dwelltime)
sumFID = _hlsvd(
FID,
dwelltime,
centralFrequency,
limits,
limitUnits=limitUnits,
numSingularValues=numSingularValues)
return FID - sumFID
def _hlsvd(FID, dwelltime, centralFrequency, limits,
limitUnits='ppm', numSingularValues=20):
"""Run hlsvdpro on FID and return spectrum modeled by HLSVD.
:param FID: Time domain data
:type FID: numpy.array
:param dwelltime: dwell time in seconds
:type dwelltime: float
:param centralFrequency: Central frequency in Hz
:type centralFrequency: float
:param limits: Frequency/ppm limits, defaults to None
:type limits: tuple, optional
:param limitUnits: Axis that limits are given in. By Default
in ppm, relative to receiver frequency (no shift). Can be 'Hz', 'ppm'
or 'ppm+shift'. Defaults to 'ppm'
:type limitUnits: str, optional
:param numSingularValues: Max number of singular values, defaults to 20
:type numSingularValues: int, optional
:return: HLSVD modeled FID
"""
m = FID.size // 2
r = hlsvdpropy.hlsvdpro(FID, numSingularValues, m=m, sparse=True)
r = hlsvdpropy.convert_hlsvd_result(r, dwelltime)
......@@ -72,8 +129,7 @@ def hlsvd(FID, dwelltime, centralFrequency, limits,
sumFID += a * np.exp((timeAxis / d)
+ 1j * 2 * np.pi
* (f * timeAxis + p / 360.0))
return FID - sumFID
return sumFID
def hlsvd_report(inFID,
......
Markdown is supported
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