Commit 57bea33f authored by William Clarke's avatar William Clarke
Browse files

Merge branch 'enh/dyn_based_align' into 'master'

Enh: dyn based align

See merge request !38
parents 9a8d16f7 0b35a1aa
Pipeline #12082 passed with stages
in 16 minutes and 54 seconds
This document contains the FSL-MRS release history in reverse chronological order.
1.1.9 (WIP)
-------------------------------
1.1.9 (Tuesday 30th November 2021)
----------------------------------
- Fixed typos in fsl_mrs_proc help.
- Fixed simulator bug for edited sequence coherence filters.
- Modified API of syntheticFromBasis function.
......@@ -11,6 +11,8 @@ This document contains the FSL-MRS release history in reverse chronological orde
- Added MH sample option to fsl_mrs, matches interactive python interface.
- Changes to the dynamic fitting results API.
- Allow tissue fractions with larger errors, but normalise. Error thrown if sum < 0.9.
- Specialist phase & frequency alignment via dynamic fitting added.
- Added fsl_mrs_preproc_edit as a script for preprocessing edited data.
1.1.8 (Tuesday 5th October 2021)
--------------------------------
......
#!/usr/bin/env python
# fsl_mrs_preproc_edit - wrapper script for (MEGA) edited single voxel MRS preprocessing
#
# Author: William Clarke <william.clarke@ndcn.ox.ac.uk>
#
# Copyright (C) 2021 University of Oxford
# SHBASECOPYRIGHT
# Quick imports
from fsl_mrs.auxiliary import configargparse
from fsl_mrs import __version__
from fsl_mrs.utils.splash import splash
import os.path as op
from os import mkdir
# NOTE!!!! THERE ARE MORE IMPORTS IN THE CODE BELOW (AFTER ARGPARSING)
def main():
# Parse command-line arguments
p = configargparse.ArgParser(
add_config_file_help=False,
description="FSL Magnetic Resonance Spectroscopy"
" - Complete edited SVS Preprocessing")
p.add_argument('-v', '--version', action='version', version=__version__)
required = p.add_argument_group('required arguments')
optional = p.add_argument_group('additional options')
# REQUIRED ARGUMENTS
required.add_argument('--data',
required=True, type=str, metavar='<str>',
help='input files')
required.add_argument('--reference',
required=True, type=str, metavar='<str>',
help='Reference non-water suppressed file')
required.add_argument('--output',
required=True, type=str, metavar='<str>',
help='output folder')
# ADDITONAL OPTIONAL ARGUMENTS
optional.add_argument('--quant', type=str,
default=None, metavar='<str>',
help='Water reference data for'
' quantification (Optional).')
optional.add_argument('--ecc', type=str,
default=None, metavar='<str>',
help='Water reference data for eddy'
' current correction (Optional).')
# option to not average, this will toggle the average property to False
optional.add_argument('--noaverage', action="store_false", dest='average',
help='Do not average repetitions.')
optional.add_argument('--hlsvd', action="store_true",
help='Apply HLSVD for residual water removal.')
optional.add_argument('--dynamic_align', action="store_true",
help='Align spectra based on dynamic fitting. Must specify dynamic basis sets.')
optional.add_argument('--dynamic_basis', type=str, nargs=2,
help='Paths to the two editing condition basis sets. '
'Only needed if --dynamic_align is specified.')
optional.add_argument('--leftshift', type=int, metavar='POINTS',
help='Remove points at the start of the fid.')
optional.add_argument('--dynamic_align_ppm', type=float, nargs=2, default=(1.9, 4.2),
help='PPM limit for dynamics dimension alignment. Default=(0.0, 4.2).')
optional.add_argument('--edit_align_ppm', type=float, nargs=2, default=(2.5, 4.2),
help='PPM limit for edit dimension alignment. Default=(2.5, 4.2).')
optional.add_argument('--t1', type=str, default=None, metavar='IMAGE',
help='structural image (for report)')
optional.add_argument('--verbose', action="store_true",
help='spit out verbose info')
optional.add_argument('--conjugate', action="store_true",
help='apply conjugate to FID')
optional.add_argument('--overwrite', action="store_true",
help='overwrite existing output folder')
optional.add_argument('--report', action="store_true",
help='Generate report in output folder')
optional.add('--config', required=False, is_config_file=True,
help='configuration file')
# Parse command-line arguments
args = p.parse_args()
# Output kickass splash screen
if args.verbose:
splash(logo='mrs')
# ######################################################
# DO THE IMPORTS AFTER PARSING TO SPEED UP HELP DISPLAY
import shutil
import numpy as np
from fsl_mrs.utils.preproc import nifti_mrs_proc
from fsl_mrs.utils.preproc import dyn_based_proc as dproc
from fsl_mrs.utils import nifti_mrs_tools as ntools
from fsl_mrs.utils import plotting
from fsl_mrs.utils import mrs_io
# ######################################################
# Check if output folder exists
overwrite = args.overwrite
if op.exists(args.output):
if not overwrite:
print(f"Folder '{args.output}' exists."
"Are you sure you want to delete it? [Y,N]")
response = input()
overwrite = response.upper() == "Y"
if not overwrite:
print('Early stopping...')
exit()
else:
shutil.rmtree(args.output)
mkdir(args.output)
else:
mkdir(args.output)
# Save chosen arguments
with open(op.join(args.output, "options.txt"), "w") as f:
f.write(str(args))
f.write("\n--------\n")
f.write(p.format_values())
if args.report:
report_dir = args.output
else:
# Will suppress report generation
report_dir = None
# ###### Do the work #######
if args.verbose:
print('Load the data....')
# Read the data
# Suppressed data
supp_data = mrs_io.read_FID(args.data)
if args.verbose:
print(f'.... Found data with shape {supp_data.shape}.\n\n')
# Reference data
ref_data = mrs_io.read_FID(args.reference)
if args.verbose:
print(f'.... Found reference with shape {ref_data.shape}.\n\n')
if args.quant is not None:
# Separate quant data
quant_data = mrs_io.read_FID(args.quant)
if args.verbose:
print(f'.... Found quant with shape {quant_data.shape}.\n\n')
else:
quant_data = None
if args.ecc is not None:
# Separate ecc data
ecc_data = mrs_io.read_FID(args.ecc)
if args.verbose:
print(f'.... Found ecc with shape {ecc_data.shape}.\n\n')
else:
ecc_data = None
# Data conjugation
if args.conjugate:
if args.verbose:
print('Conjugation explicitly set,'
'applying conjugation.')
supp_data = nifti_mrs_proc.conjugate(supp_data)
ref_data = nifti_mrs_proc.conjugate(ref_data)
if args.quant is not None:
quant_data = nifti_mrs_proc.conjugate(quant_data)
if args.ecc is not None:
ecc_data = nifti_mrs_proc.conjugate(ecc_data)
# Determine if coils have been combined already
if args.verbose:
print('.... Determine if coil combination is needed')
if 'DIM_COIL' in supp_data.dim_tags:
do_coil_combine = True
if args.verbose:
print(' ----> YES.\n')
else:
do_coil_combine = False
if args.verbose:
print(' ----> NO.\n')
# Do preproc
if args.verbose:
print('Begin proprocessing.... ')
if do_coil_combine:
if args.verbose:
print('... Coil Combination ...')
if 'DIM_DYN' in ref_data.dim_tags:
avg_ref_data = nifti_mrs_proc.average(ref_data, 'DIM_DYN')
else:
avg_ref_data = ref_data
if 'DIM_EDIT' in avg_ref_data.dim_tags:
avg_ref_data = nifti_mrs_proc.average(avg_ref_data, 'DIM_EDIT')
else:
avg_ref_data = avg_ref_data
supp_data = nifti_mrs_proc.coilcombine(supp_data, reference=avg_ref_data, report=report_dir)
ref_data = nifti_mrs_proc.coilcombine(ref_data, reference=avg_ref_data)
if args.quant is not None:
quant_data = nifti_mrs_proc.coilcombine(quant_data, reference=avg_ref_data)
if args.ecc is not None:
ecc_data = nifti_mrs_proc.coilcombine(ecc_data, reference=avg_ref_data)
# Frequency and phase align the FIDs in the ON/OFF condition
if args.verbose:
print('... Align Dynamics ...')
supp_data = nifti_mrs_proc.align(
supp_data,
'DIM_DYN',
ppmlim=args.dynamic_align_ppm,
niter=4,
apodize=0.0,
report=report_dir,
report_all=True)
if args.dynamic_align:
# Run dynamic fitting based alignment
edit_0, edit_1 = ntools.split(supp_data, 'DIM_EDIT', 0)
basis_0 = mrs_io.read_basis(args.dynamic_basis[0])
basis_1 = mrs_io.read_basis(args.dynamic_basis[1])
fitargs_0 = {'baseline_order': 1,
'model': 'lorentzian',
'ppmlim': (0.2, 4.2)}
fitargs_1 = {'baseline_order': 1,
'model': 'lorentzian',
'ppmlim': (2.0, 4.2)}
edit_0_aligned, eps0, phi0 = dproc.align_by_dynamic_fit(edit_0, basis_0, fitargs_0)
edit_1_aligned, eps1, phi1 = dproc.align_by_dynamic_fit(edit_1, basis_1, fitargs_1)
if report_dir is not None:
dproc.align_by_dynamic_fit_report(edit_0, edit_0_aligned, eps0, phi0, html=report_dir)
dproc.align_by_dynamic_fit_report(edit_1, edit_1_aligned, eps1, phi1, html=report_dir)
supp_data = ntools.merge([edit_0_aligned, edit_1_aligned], 'DIM_EDIT')
if 'DIM_DYN' in ref_data.dim_tags:
ref_data = nifti_mrs_proc.align(ref_data, 'DIM_DYN', ppmlim=(0, 8))
if args.quant is not None and ('DIM_DYN' in quant_data.dim_tags):
quant_data = nifti_mrs_proc.align(quant_data, 'DIM_DYN', ppmlim=(0, 8))
if args.ecc is not None and ('DIM_DYN' in ecc_data.dim_tags):
ecc_data = nifti_mrs_proc.align(ecc_data, 'DIM_DYN', ppmlim=(0, 8))
# Average the data (if asked to do so)
if args.average:
if args.verbose:
print('... Average FIDs ...')
supp_data = nifti_mrs_proc.average(supp_data, 'DIM_DYN', report=report_dir)
if 'DIM_DYN' in ref_data.dim_tags:
ref_data = nifti_mrs_proc.average(ref_data, 'DIM_DYN')
if args.quant is not None and ('DIM_DYN' in quant_data.dim_tags):
quant_data = nifti_mrs_proc.average(quant_data, 'DIM_DYN')
# Always average ecc if it exists as a separate scan
if args.ecc is not None and ('DIM_DYN' in ecc_data.dim_tags):
ecc_data = nifti_mrs_proc.average(ecc_data, 'DIM_DYN')
# Eddy current correction
if args.verbose:
print('... ECC correction ...')
if args.ecc is not None:
eccRef = ecc_data
else:
if 'DIM_DYN' in ref_data.dim_tags:
eccRef = nifti_mrs_proc.average(ref_data, 'DIM_DYN')
else:
eccRef = ref_data
supp_data = nifti_mrs_proc.ecc(supp_data, eccRef, report=report_dir)
ref_data = nifti_mrs_proc.ecc(ref_data, eccRef)
# Typically if a separate "quantification" water reference
# has been collected it will have most gradients removed
# (OVS and water suppression), therefore use it as it's own reference.
if args.quant is not None:
quant_data = nifti_mrs_proc.ecc(quant_data, quant_data)
# HLSVD
if args.hlsvd:
if args.verbose:
print('... Residual water removal ...')
hlsvdlimits = [-0.25, 0.25]
supp_data = nifti_mrs_proc.remove_peaks(supp_data, hlsvdlimits, limit_units='ppm', report=report_dir)
if args.leftshift:
if args.verbose:
print('... Truncation ...')
supp_data = nifti_mrs_proc.truncate_or_pad(supp_data, -args.leftshift, 'first', report=report_dir)
ref_data = nifti_mrs_proc.truncate_or_pad(ref_data, -args.leftshift, 'first')
if args.quant is not None:
quant_data = nifti_mrs_proc.truncate_or_pad(quant_data, -args.leftshift, 'first')
# Apply shift to reference
if args.verbose:
print('... Shift Cr to 3.027 ...')
supp_data = nifti_mrs_proc.shift_to_reference(supp_data, 3.027, (2.9, 3.1), report=report_dir)
# Apply phasing based on a single peak (Cr)
if args.verbose:
print('... Phasing on tCr ...')
supp_data = nifti_mrs_proc.phase_correct(supp_data, (2.9, 3.1), report=report_dir)
if args.quant is not None:
if 'DIM_EDIT' in quant_data.dim_tags:
quant_data = nifti_mrs_proc.align(quant_data, 'DIM_EDIT', ppmlim=(0, 8), niter=1, apodize=0.0)
quant_data = nifti_mrs_proc.average(quant_data, 'DIM_EDIT')
final_wref = nifti_mrs_proc.phase_correct(quant_data, (4.55, 4.7), hlsvd=False, report=report_dir)
else:
if 'DIM_EDIT' in ref_data.dim_tags:
ref_data = nifti_mrs_proc.align(ref_data, 'DIM_EDIT', ppmlim=(0, 8), niter=1, apodize=0.0)
ref_data = nifti_mrs_proc.average(ref_data, 'DIM_EDIT')
final_wref = nifti_mrs_proc.phase_correct(ref_data, (4.55, 4.7), hlsvd=False, report=report_dir)
# Align between edit dimension
if args.dynamic_align and args.average:
metab_edit_align, eps, phi = dproc.align_by_dynamic_fit(supp_data, [basis_0, basis_1], fitargs_1)
if report_dir is not None:
dproc.align_by_dynamic_fit_report(supp_data, metab_edit_align, eps, phi, html=report_dir)
elif args.dynamic_align:
# Do not average and dynamic alignment
# Work out alignement on averaged data
supp_data_avg = nifti_mrs_proc.average(supp_data, 'DIM_DYN')
metab_edit_align, eps, phi = dproc.align_by_dynamic_fit(supp_data_avg, [basis_0, basis_1], fitargs_1)
# Split, apply correction, recombine...
unavg_0 = ntools.split(supp_data, 'DIM_EDIT', 0)
unavg_1 = ntools.split(supp_data, 'DIM_EDIT', 0)
unavg_0_aligned = nifti_mrs_proc.fshift(unavg_0, eps[0] * 1 / (2 * np.pi))
unavg_0_aligned = nifti_mrs_proc.apply_fixed_phase(unavg_0_aligned, phi[0] * 180 / np.pi)
unavg_1_aligned = nifti_mrs_proc.fshift(unavg_1, eps[1] * 1 / (2 * np.pi))
unavg_1_aligned = nifti_mrs_proc.apply_fixed_phase(unavg_1_aligned, phi[1] * 180 / np.pi)
metab_edit_align = ntools.merge([unavg_0_aligned, unavg_1_aligned], 'DIM_EDIT')
else:
metab_edit_align = nifti_mrs_proc.align(
supp_data,
'DIM_EDIT',
ppmlim=args.edit_align_ppm,
niter=4,
apodize=0.0,
report=report_dir)
# Differencing
metab_diff = nifti_mrs_proc.subtract(metab_edit_align, dim='DIM_EDIT', figure=False)
metab_diff = nifti_mrs_proc.apply_fixed_phase(metab_diff, 180.0, report=report_dir)
cond_0, cond_1 = ntools.split(metab_edit_align, 'DIM_EDIT', 0)
# Save the data
if args.verbose:
print('... Saving data ...')
metab_diff.save(op.join(args.output, 'diff'))
cond_0.save(op.join(args.output, 'edit_0'))
cond_1.save(op.join(args.output, 'edit_1'))
final_wref.save(op.join(args.output, 'wref'))
# Produce full html report
if args.report:
import subprocess
import glob
if args.verbose:
print('Create report')
htmlfiles = glob.glob(op.join(args.output, '*.html'))
subprocess.call(['merge_mrs_reports', '-d',
op.join(args.output, 'metab'),
'-o', args.output,
'--delete'] + htmlfiles)
if args.t1 is not None:
fig = plotting.plot_world_orient(args.t1, args.data)
fig.savefig(op.join(args.output, 'voxel_location.png'))
if __name__ == '__main__':
main()
'''FSL-MRS test script
Test the edited svs preprocessing script
Copyright Will Clarke, University of Oxford, 2021'''
import subprocess
from pathlib import Path
from fsl_mrs.utils import mrs_io
testsPath = Path(__file__).parent
data = testsPath / 'testdata/fsl_mrs_preproc_edit'
t1 = str(testsPath / 'testdata/svs_segment/T1.anat/T1_biascorr.nii.gz')
def test_preproc(tmp_path):
metab = str(data / 'metab_raw.nii.gz')
wrefc = str(data / 'wref_internal.nii.gz')
wrefq = str(data / 'wref_quant.nii.gz')
ecc = str(data / 'wref_internal.nii.gz')
retcode = subprocess.check_call(
['fsl_mrs_preproc_edit',
'--output', str(tmp_path),
'--data', metab,
'--reference', wrefc,
'--quant', wrefq,
'--ecc', ecc,
'--t1', t1,
'--dynamic_align_ppm', '1.8', '4.2',
'--edit_align_ppm', '2.5', '4.2',
'--hlsvd',
'--leftshift', '2',
'--overwrite',
'--report',
'--verbose'])
assert retcode == 0
assert (tmp_path / 'diff.nii.gz').exists()
assert (tmp_path / 'edit_0.nii.gz').exists()
assert (tmp_path / 'edit_1.nii.gz').exists()
assert (tmp_path / 'wref.nii.gz').exists()
assert (tmp_path / 'options.txt').exists()
assert (tmp_path / 'mergedReports.html').exists()
assert (tmp_path / 'voxel_location.png').exists()
def test_preproc_noavg(tmp_path):
metab = str(data / 'metab_raw.nii.gz')
wrefc = str(data / 'wref_internal.nii.gz')
wrefq = str(data / 'wref_quant.nii.gz')
ecc = str(data / 'wref_internal.nii.gz')
retcode = subprocess.check_call(
['fsl_mrs_preproc_edit',
'--output', str(tmp_path),
'--data', metab,
'--reference', wrefc,
'--quant', wrefq,
'--ecc', ecc,
'--t1', t1,
'--dynamic_align_ppm', '1.8', '4.2',
'--edit_align_ppm', '2.5', '4.2',
'--leftshift', '2',
'--overwrite',
'--report',
'--verbose',
'--noaverage'])
assert retcode == 0
assert (tmp_path / 'diff.nii.gz').exists()
assert (tmp_path / 'edit_0.nii.gz').exists()
assert (tmp_path / 'edit_1.nii.gz').exists()
assert (tmp_path / 'wref.nii.gz').exists()
assert (tmp_path / 'options.txt').exists()
assert (tmp_path / 'mergedReports.html').exists()
assert (tmp_path / 'voxel_location.png').exists()
diff = mrs_io.read_FID(tmp_path / 'diff.nii.gz')
on = mrs_io.read_FID(tmp_path / 'edit_0.nii.gz')
off = mrs_io.read_FID(tmp_path / 'edit_1.nii.gz')
assert diff.shape[-1] == 16
assert on.shape[-1] == 16
assert off.shape[-1] == 16
......@@ -150,7 +150,8 @@ def test_dynMRS_fit(fixed_ratio_mrs):
baseline_order=0,
metab_groups=[0, 0],
rescale=False)
res, _ = dyn_obj.fit()
init = dyn_obj.initialise(indiv_init=None)
res, _ = dyn_obj.fit(init=init)
concs = res.data_frame.filter(like='conc').to_numpy()
assert np.allclose(concs, [1, 1, 1, 1], atol=0.1)
......
'''FSL-MRS test script
Test dynamic fitting based preprocessing
Copyright Will Clarke, University of Oxford, 2021
'''
from pathlib import Path
from fsl_mrs.utils.preproc import nifti_mrs_proc as nproc
from fsl_mrs.utils.preproc import dyn_based_proc as dproc
from fsl_mrs.utils.mrs_io import read_FID, read_basis
from fsl_mrs.utils.nifti_mrs_tools import split
from fsl_mrs.utils import basis_tools as btools
testsPath = Path(__file__).parent
data = testsPath / 'testdata' / 'fsl_mrs_preproc'
metab = data / 'metab_raw.nii.gz'
wrefc = data / 'wref_raw.nii.gz'
basis_path = testsPath / 'testdata' / 'fsl_mrs' / 'steam_basis'
def test_dyn_align(tmp_path):
nmrs_obj = read_FID(metab)
nmrs_ref_obj = read_FID(wrefc)
nmrs_ref_obj = nproc.average(nmrs_ref_obj, 'DIM_DYN')
combined = nproc.coilcombine(nmrs_obj, reference=nmrs_ref_obj)
reduced_data, _ = split(combined, 'DIM_DYN', 2)
aligned_1 = nproc.align(reduced_data, 'DIM_DYN', ppmlim=(0.2, 4.2), apodize=0.0)
basis = btools.conjugate_basis(read_basis(basis_path))
fitargs = {'ppmlim': (0.2, 4.2), 'baseline_order': 1}
aligned_2 = dproc.align_by_dynamic_fit(aligned_1, basis, fitargs=fitargs)
dproc.align_by_dynamic_fit_report(
aligned_1,
aligned_2[0],
aligned_2[1],
aligned_2[2],
ppmlim=fitargs['ppmlim'],
html=str(tmp_path / 'align_report.html'))
assert aligned_2[0].hdr_ext['ProcessingApplied'][2]['Method'] == 'Frequency and phase correction'
assert aligned_2[0].hdr_ext['ProcessingApplied'][2]['Details']\
== "fsl_mrs.utils.preproc.dyn_based_proc.align_by_dynamic_fit, "\
"fitargs={'ppmlim': (0.2, 4.2), 'baseline_order': 1}."
assert (tmp_path / 'align_report.html').is_file()
"""Model for dynamic fitting based alignment
Author: William Clarke <william.clarke@ndcn.ox.ac.uk>
Copyright (C) 2021 University of Oxford
"""
Parameters = {
'Phi_0': 'variable',
'Phi_1': 'fixed',
'conc': 'fixed',
'eps': 'variable',
'gamma': 'fixed',
'sigma': 'fixed',
'baseline': 'fixed'
}
Bounds = {
'gamma': (0, None),
'sigma': (0, None)
}
# dyn_based_proc.py - Module containing processing functions based on dynamic fitting
#
# Author: William Clarke <william.clarke@ndcn.ox.ac.uk>
#
# Copyright (C) 2021 University of Oxford
# SHBASECOPYRIGHT
import os.path as op