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

Merge branch 'multi_dim_ecc'

Handling of multi-dimensional data for dynamic preprocessing in ecc and align tool.
parents 8ff66716 db43d37c
This document contains the FSL-MRS release history in reverse chronological order.
1.1.5 (WIP)
-------------------------------
- fsl_mrs_proc align can now align across all higher dimension FIDs. Pass 'all' as dimension tag.
1.1.4 (Tuesday 3rd August 2021)
-------------------------------
- Fixed bug in calculation of molality concentration. Tissue mole fractions had been swapped for tissue volume fractions. Molar concentrations unaffected.
......
......@@ -428,6 +428,21 @@ class NIFTI_MRS(Image):
else:
raise TypeError('dim should be int or a string matching one of the dim tags.')
def iterate_over_spatial(self):
"""Iterate over spatial voxels yeilding a data array the shape of the FID and any higher dimensions + index.
:yield: Complex FID data with any higher dimensions. Index to data.
:rtype: tuple
"""
data = self.data
def calc_slice_idx(idx):
slice_obj = list(idx[:3]) + [slice(None), ] * (self.data.ndim - 3)
return tuple(slice_obj)
for idx in np.ndindex(data.shape[:3]):
yield self.data[idx], calc_slice_idx(idx)
def generate_mrs(self, dim=None, basis_file=None, basis=None, ref_data=None):
"""Generator for MRS or MRSI objects from the data, optionally returning a whole dimension as a list.
......
......@@ -85,6 +85,7 @@ def main():
help='List of files to align')
align_group.add_argument('--dim', type=str, default='DIM_DYN',
help='NIFTI-MRS dimension tag to align across.'
'Or "all" to align over all spectra in higer dimensions.'
'Default = DIM_DYN')
align_group.add_argument('--ppm', type=float, nargs=2,
metavar='<lower-limit upper-limit>',
......@@ -473,7 +474,9 @@ def average(dataobj, args):
def align(dataobj, args):
if args['dim'] not in dataobj.data.dim_tags:
if args['dim'].lower() == 'all':
pass
elif args['dim'] not in dataobj.data.dim_tags:
raise InappropriateDataError(f'Data ({dataobj.datafilename}) has no {args["dim"]} dimension.'
f' Dimensions are is {dataobj.data.dim_tags}.')
......
......@@ -91,6 +91,18 @@ def test_nifti_mrs_generator():
break
def test_nifti_mrs_spatial_generator():
obj = NIFTI_MRS(data['unprocessed'])
for gen_data, slice_idx in obj.iterate_over_spatial():
assert gen_data.shape == (4096, 32, 64)
assert slice_idx == (0, 0, 0,
slice(None, None, None),
slice(None, None, None),
slice(None, None, None))
break
def test_nifti_mrs_gen_mrs():
obj = NIFTI_MRS(data['unprocessed'])
......
......@@ -195,11 +195,42 @@ def mrsi_data_diff(tmp_path):
return testfile, nmrs
@pytest.fixture
def svs_data_uncomb_reps(tmp_path):
coils = 2
noiseconv = 0.1 * np.eye(coils)
coilamps = np.random.randn(coils)
coilphs = np.random.random(coils) * 2 * np.pi
FID, hdr = syntheticFID(noisecovariance=noiseconv,
coilamps=coilamps,
coilphase=coilphs,
points=512)
FID2, hdr = syntheticFID(noisecovariance=noiseconv,
coilamps=coilamps,
coilphase=coilphs,
points=512)
FID = np.stack((np.asarray(FID).T, np.asarray(FID2).T), axis=2)
FID = FID.reshape((1, 1, 1) + FID.shape)
nmrs = gen_new_nifti_mrs(FID,
hdr['dwelltime'],
hdr['centralFrequency'],
dim_tags=['DIM_COIL', 'DIM_DYN', None])
testname = 'svsdata_uncomb_reps.nii'
testfile = op.join(tmp_path, testname)
nmrs.save(testfile)
return testfile, nmrs
def splitdata(svs, mrsi):
return svs[0], mrsi[0], svs[1], mrsi[1]
def test_filecreation(svs_data, mrsi_data, svs_data_uncomb, mrsi_data_uncomb):
def test_filecreation(svs_data, mrsi_data, svs_data_uncomb, mrsi_data_uncomb, svs_data_uncomb_reps):
svsfile, mrsifile, svsdata, mrsidata = splitdata(svs_data, mrsi_data)
data = read_FID(svsfile)
......@@ -221,6 +252,11 @@ def test_filecreation(svs_data, mrsi_data, svs_data_uncomb, mrsi_data_uncomb):
assert data.shape == (2, 2, 2, 512, 4)
assert np.allclose(data.data, mrsidata.data)
# Uncombined and reps
data = read_FID(svs_data_uncomb_reps[0])
assert data.shape == (1, 1, 1, 512, 2, 2)
assert np.allclose(data.data, svs_data_uncomb_reps[1].data)
def test_coilcombine(svs_data_uncomb, mrsi_data_uncomb, tmp_path):
svsfile, mrsifile, svsdata, mrsidata = splitdata(svs_data_uncomb,
......@@ -313,7 +349,6 @@ def test_align(svs_data, mrsi_data, tmp_path):
assert np.allclose(data.data, directRun.data)
# Run coil combination on both sets of data using the command line
subprocess.check_call(['fsl_mrs_proc',
'align',
'--dim', 'DIM_DYN',
......@@ -331,6 +366,27 @@ def test_align(svs_data, mrsi_data, tmp_path):
assert np.allclose(data.data, directRun.data, atol=1E-1, rtol=1E-1)
def test_align_all(svs_data_uncomb_reps, tmp_path):
svsfile, svsdata = svs_data_uncomb_reps[0], svs_data_uncomb_reps[1]
# Run align on both sets of data using the command line
subprocess.check_call(['fsl_mrs_proc',
'align',
'--dim', 'all',
'--file', svsfile,
'--ppm', '-10', '10',
'--output', tmp_path,
'--filename', 'tmp'])
# Load result for comparison
data = read_FID(op.join(tmp_path, 'tmp.nii.gz'))
# Run directly
directRun = preproc.align(svsdata, 'all', ppmlim=(-10, 10))
assert np.allclose(data.data, directRun.data)
def test_ecc(svs_data, mrsi_data, tmp_path):
svsfile, mrsifile, svsdata, mrsidata = splitdata(svs_data, mrsi_data)
......
......@@ -75,6 +75,14 @@ def test_align():
== 'fsl_mrs.utils.preproc.nifti_mrs_proc.align, dim=DIM_DYN, '\
'target=None, ppmlim=(1.0, 4.0), niter=1, apodize=5.'
# Align across all spectra
aligned3 = nproc.align(with_coils, 'all', ppmlim=(1.0, 4.0), niter=1, apodize=5)
assert aligned3.hdr_ext['ProcessingApplied'][0]['Method'] == 'Frequency and phase correction'
assert aligned3.hdr_ext['ProcessingApplied'][0]['Details']\
== 'fsl_mrs.utils.preproc.nifti_mrs_proc.align, dim=all, '\
'target=None, ppmlim=(1.0, 4.0), niter=1, apodize=5.'
def test_aligndiff():
# For want of data this is a bizzare way of using this function.
......
......@@ -161,11 +161,11 @@ def average(data, dim, figure=False, report=None, report_all=False):
def align(data, dim, target=None, ppmlim=None, niter=2, apodize=10, figure=False, report=None, report_all=False):
'''Align frequencies of spectra across a dimension
specified by a tag.
'''Align frequency and phase of spectra. Can be run across a dimension (specified by a tag), or all spectra
stored in higher dimensions.
:param NIFTI_MRS data: Data to align
:param str dim: NIFTI-MRS dimension tag
:param str dim: NIFTI-MRS dimension tag, or 'all'
:param target: Optional target FID
:param ppmlim: ppm search limits.
:param int niter: Number of total iterations
......@@ -178,9 +178,19 @@ def align(data, dim, target=None, ppmlim=None, niter=2, apodize=10, figure=False
'''
aligned_obj = data.copy()
for dd, idx in data.iterate_over_dims(dim=dim,
iterate_over_space=True,
reduce_dim_index=False):
if dim.lower() == 'all':
generator = data.iterate_over_spatial()
else:
generator = data.iterate_over_dims(dim=dim,
iterate_over_space=True,
reduce_dim_index=False)
for dd, idx in generator:
if dim == 'all':
# flatten
original_shape = dd.shape
dd = dd.reshape(original_shape[0], -1)
out = preproc.phase_freq_align(
dd.T,
......@@ -193,12 +203,15 @@ def align(data, dim, target=None, ppmlim=None, niter=2, apodize=10, figure=False
verbose=False,
target=target)
aligned_obj[idx], phi, eps = out[0].T, out[1], out[2]
if dim == 'all':
aligned_obj[idx], phi, eps = out[0].T.reshape(original_shape), out[1], out[2]
else:
aligned_obj[idx], phi, eps = out[0].T, out[1], out[2]
if (figure or report) and (report_all or first_index(idx)):
from fsl_mrs.utils.preproc.align import phase_freq_align_report
fig = phase_freq_align_report(dd.T,
aligned_obj[idx].T,
out[0],
phi,
eps,
data.bandwidth,
......@@ -333,8 +346,10 @@ def ecc(data, reference, figure=False, report=None, report_all=False):
for dd, idx in data.iterate_over_dims(iterate_over_space=True):
if data.shape == reference.shape:
# Reference is the same shape as data, voxel-wise and spectrum-wise iteration
ref = reference[idx]
else:
# Only one reference FID, only iterate over spatial voxels.
ref = reference[idx[0], idx[1], idx[2], :]
corrected_obj[idx] = preproc.eddy_correct(dd, ref)
......
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