Commit 076399d2 authored by William Clarke's avatar William Clarke
Browse files

TNew quantification structure with fitted water and T1 relaxation.

parent 0443d171
......@@ -67,7 +67,7 @@ def main():
fitting_args.add_argument('--baseline_order', default=2, type=int,
metavar=('ORDER'),
help='order of baseline polynomial'
' (default=2,-1 disables)')
' (default=2, -1 disables)')
fitting_args.add_argument('--metab_groups', default=0, nargs='+',
type=str_or_int_arg,
help='metabolite groups: list of groups'
......@@ -89,6 +89,8 @@ def main():
help='structural image (for report)')
optional.add_argument('--TE', type=float, default=None, metavar='TE',
help='Echo time for relaxation correction (ms)')
optional.add_argument('--TR', type=float, default=None, metavar='TR',
help='Repetition time for relaxation correction (s)')
optional.add_argument('--tissue_frac', type=tissue_frac_arg,
action=TissueFracAction, nargs='+',
default=None, metavar='WM GM CSF OR json',
......@@ -105,8 +107,6 @@ def main():
help='output html report')
optional.add_argument('--verbose', action="store_true",
help='spit out verbose info')
optional.add_argument('--phase_correct', action="store_true",
help='do phase correction')
optional.add_argument('--overwrite', action="store_true",
help='overwrite existing output folder')
optional.add_argument('--conj_fid', dest='conjfid', action="store_true",
......@@ -147,6 +147,7 @@ def main():
from fsl_mrs.utils import fitting
from fsl_mrs.utils import plotting
from fsl_mrs.utils import misc
from fsl_mrs.utils import quantify
import datetime
# ######################################################
if not args.verbose:
......@@ -226,12 +227,6 @@ def main():
if not args.no_rescale:
mrs.rescaleForFitting(ind_scaling=args.ind_scale)
# Do phase correction
if args.phase_correct:
if args.verbose:
print('--->> Phase correction\n')
mrs.FID = misc.phase_correct(mrs, mrs.FID)
# Keep/Ignore metabolites
mrs.keep(args.keep)
mrs.ignore(args.ignore)
......@@ -284,40 +279,45 @@ def main():
# Echo time
if args.TE is not None:
echotime = args.TE * 1E-3
elif 'meta' in basisheader and 'TE' in basisheader['meta']:
echotime = basisheader['meta']['TE']
if echotime > 1.0: # Assume in ms.
echotime *= 1E-3
else:
echotime = None
elif 'EchoTime' in FID.hdr_ext:
echotime = FID.hdr_ext['EchoTime']
else:
echotime = None
# Repetition time
if args.TR is not None:
repetition_time = args.TR
elif 'RepetitionTime' in FID.hdr_ext:
repetition_time = FID.hdr_ext['RepetitionTime']
else:
repetition_time = None
# Internal and Water quantification if requested
if (mrs.H2O is None) or (echotime is None):
if (mrs.H2O is None) or (echotime is None) or (repetition_time is None):
if echotime is None:
warnings.warn('H2O file provided but could not determine TE:'
' no absolute quantification will be performed.',
UserWarning)
res.calculateConcScaling(mrs, referenceMetab=args.internal_ref)
elif args.tissue_frac is not None:
res.calculateConcScaling(mrs,
referenceMetab=args.internal_ref,
waterRefFID=mrs.H2O,
tissueFractions=args.tissue_frac,
TE=echotime,
verbose=args.verbose,
add_scale=args.h2o_scale)
if repetition_time is None:
warnings.warn('H2O file provided but could not determine TR:'
' no absolute quantification will be performed.',
UserWarning)
res.calculateConcScaling(mrs, internal_reference=args.internal_ref, verbose=args.verbose)
else:
# Form quantification information
q_info = quantify.QuantificationInfo(echotime,
repetition_time,
mrs.names,
mrs.centralFrequency / 1E6)
if args.tissue_frac:
q_info.set_fractions(args.tissue_frac)
if args.h2o_scale:
q_info.add_corr = args.h2o_scale
res.calculateConcScaling(mrs,
referenceMetab=args.internal_ref,
waterRefFID=mrs.H2O,
tissueFractions=None,
TE=echotime,
verbose=args.verbose,
add_scale=args.h2o_scale)
quant_info=q_info,
internal_reference=args.internal_ref,
verbose=args.verbose)
# Combine metabolites.
if args.combine is not None:
......
......@@ -7,12 +7,14 @@
#
# Copyright (C) 2019 University of Oxford
from fsl_mrs.auxiliary import configargparse
import time
import warnings
from fsl_mrs.auxiliary import configargparse
from fsl_mrs import __version__
from fsl_mrs.utils.splash import splash
from fsl_mrs.utils import fitting, misc, mrs_io
import time
from fsl_mrs.utils import fitting, misc, mrs_io, quantify
# NOTE!!!! THERE ARE MORE IMPORTS IN THE CODE BELOW (AFTER ARGPARSING)
......@@ -83,6 +85,8 @@ def main():
# ADDITONAL OPTIONAL ARGUMENTS
optional.add_argument('--TE', type=float, default=None, metavar='TE',
help='Echo time for relaxation correction (ms)')
optional.add_argument('--TR', type=float, default=None, metavar='TR',
help='Repetition time for relaxation correction (s)')
optional.add_argument('--tissue_frac', type=str, nargs=3, default=None,
help='Tissue fraction nifti files registered to'
' MRSI volume. Supplied in order: WM, GM, CSF.')
......@@ -122,7 +126,6 @@ def main():
# DO THE IMPORTS AFTER PARSING TO SPEED UP HELP DISPLAY
import os
import shutil
import warnings
import numpy as np
from fsl_mrs.utils import report
from fsl_mrs.core import NIFTI_MRS
......@@ -212,16 +215,17 @@ def main():
# Echo time
if args.TE is not None:
echotime = args.TE * 1E-3
elif 'meta' in mrsi.basis_hdr and 'TE' in mrsi.basis_hdr['meta']:
echotime = mrsi.basis_hdr['meta']['TE']
if echotime > 1.0: # Assume in ms.
echotime *= 1E-3
else:
echotime = None
elif 'TE' in mrsi.header:
echotime = mrsi.header['TE']
else:
echotime = None
# Repetition time
if args.TR is not None:
repetition_time = args.TR
elif 'RepetitionTime' in mrsi_data.hdr_ext:
repetition_time = mrsi_data.hdr_ext['RepetitionTime']
else:
repetition_time = None
# Fitting
if args.verbose:
......@@ -234,7 +238,7 @@ def main():
mrs = mrsi.mrs_from_average()
Fitargs_init = Fitargs.copy()
Fitargs_init['method'] = 'Newton'
res_init, _ = runvoxel([mrs, 0, None], args, Fitargs_init, echotime)
res_init, _ = runvoxel([mrs, 0, None], args, Fitargs_init, echotime, repetition_time)
Fitargs['x0'] = res_init.params
# quick summary figure
......@@ -260,7 +264,7 @@ def main():
print(' Running sequentially (are you sure about that?) ')
results = []
for idx, mrs in enumerate(mrsi):
res = runvoxel(mrs, args, Fitargs, echotime)
res = runvoxel(mrs, args, Fitargs, echotime, repetition_time)
results.append(res)
if args.verbose:
print(f'{idx+1}/{mrsi.num_masked_voxels} voxels completed')
......@@ -268,7 +272,7 @@ def main():
if args.verbose:
print(f' Parallelising over {mp.cpu_count()} workers ')
func = partial(runvoxel, args=args, Fitargs=Fitargs, echotime=echotime)
func = partial(runvoxel, args=args, Fitargs=Fitargs, echotime=echotime, repetition_time=repetition_time)
with mp.Pool(processes=mp.cpu_count()) as p:
results = p.map_async(func, mrsi)
if args.verbose:
......@@ -406,7 +410,7 @@ def main():
print('\n\n\nDone.')
def runvoxel(mrs_in, args, Fitargs, echotime):
def runvoxel(mrs_in, args, Fitargs, echotime, repetition_time):
mrs, index, tissue_seg = mrs_in
# Parse metabolite groups
......@@ -420,22 +424,33 @@ def runvoxel(mrs_in, args, Fitargs, echotime):
new_metab_groups = metab_groups.copy()
res = fitting.fit_FSLModel(mrs, **Fitargs, metab_groups=new_metab_groups)
# Quantification
# Internal and Water quantification if requested
if (mrs.H2O is None) or (echotime is None):
res.calculateConcScaling(mrs, referenceMetab=args.internal_ref)
elif tissue_seg is None:
res.calculateConcScaling(mrs,
referenceMetab=args.internal_ref,
waterRefFID=mrs.H2O,
tissueFractions=None,
TE=echotime)
if (mrs.H2O is None) or (echotime is None) or (repetition_time is None):
if echotime is None:
warnings.warn('H2O file provided but could not determine TE:'
' no absolute quantification will be performed.',
UserWarning)
if repetition_time is None:
warnings.warn('H2O file provided but could not determine TR:'
' no absolute quantification will be performed.',
UserWarning)
res.calculateConcScaling(mrs, internal_reference=args.internal_ref, verbose=args.verbose)
else:
# Form quantification information
q_info = quantify.QuantificationInfo(echotime,
repetition_time,
mrs.names,
mrs.centralFrequency / 1E6)
if tissue_seg:
q_info.set_fractions(tissue_seg)
if args.h2o_scale:
q_info.add_corr = args.h2o_scale
res.calculateConcScaling(mrs,
referenceMetab=args.internal_ref,
waterRefFID=mrs.H2O,
tissueFractions=tissue_seg,
TE=echotime)
quant_info=q_info,
internal_reference=args.internal_ref,
verbose=args.verbose)
# Combine metabolites.
if args.combine is not None:
res.combine(args.combine)
......
......@@ -7,15 +7,59 @@ Copyright Will Clarke, University of Oxford, 2021'''
import os.path as op
import numpy as np
import fsl_mrs.utils.mrs_io as mrsio
from fsl_mrs.utils.fitting import fit_FSLModel
import numpy as np
import fsl_mrs.utils.quantify as quant
from fsl_mrs.utils.constants import STANDARD_T1, STANDARD_T2
metabfile = op.join(op.dirname(__file__), 'testdata/quantify/Cr_10mM_test_water_scaling_WS.txt')
h2ofile = op.join(op.dirname(__file__), 'testdata/quantify/Cr_10mM_test_water_scaling_nWS.txt')
basisfile = op.join(op.dirname(__file__), 'testdata/quantify/basisset_JMRUI')
def test_QuantificationInfo():
qci = quant.QuantificationInfo(0.000, 40, ['Cr', 'NAA'], 298)
assert qci.relax_corr_water_molal > 55500
assert qci.relax_corr_water_molar > 55500
qci = quant.QuantificationInfo(0.000, 40, ['Cr', 'NAA'], 127)
assert qci.relax_corr_water_molal > 55500
assert qci.relax_corr_water_molar > 55500
qci = quant.QuantificationInfo(0.010, 3, ['Cr', 'NAA'], 127)
t2s = STANDARD_T2['3T']
t1s = STANDARD_T1['3T']
assert np.isclose(qci.R_H2O_WM, np.exp(-0.010 / t2s['H2O_WM']) * (1 - np.exp(-3 / t1s['H2O_WM'])))
assert np.isclose(qci.R_H2O_GM, np.exp(-0.010 / t2s['H2O_GM']) * (1 - np.exp(-3 / t1s['H2O_GM'])))
assert np.isclose(qci.R_H2O_CSF, np.exp(-0.010 / t2s['H2O_CSF']) * (1 - np.exp(-3 / t1s['H2O_CSF'])))
qci = quant.QuantificationInfo(0.010, 3, ['Cr', 'NAA'], 298)
t2s = STANDARD_T2['7T']
t1s = STANDARD_T1['7T']
assert np.isclose(qci.R_H2O_WM, np.exp(-0.010 / t2s['H2O_WM']) * (1 - np.exp(-3 / t1s['H2O_WM'])))
assert np.isclose(qci.R_H2O_GM, np.exp(-0.010 / t2s['H2O_GM']) * (1 - np.exp(-3 / t1s['H2O_GM'])))
assert np.isclose(qci.R_H2O_CSF, np.exp(-0.010 / t2s['H2O_CSF']) * (1 - np.exp(-3 / t1s['H2O_CSF'])))
assert qci.ref_metab == 'Cr'
assert qci.ref_protons == 5
assert qci.ref_limits == (2, 5)
qci = quant.QuantificationInfo(0.010, 3, ['NAA'], 298)
assert qci.ref_metab == 'NAA'
assert qci.ref_protons == 3
assert qci.ref_limits == (1.8, 2.2)
qci.set_fractions({'GM': 0.45, 'WM': 0.45, 'CSF': 0.1})
assert qci._fractions is not None
assert np.isclose(qci.csf_corr, 1 / 0.9)
qci.add_corr = 5.0
assert qci.add_corr == 5.0
def test_quantifyWater():
basis, names, headerb = mrsio.read_basis(basisfile)
crIndex = names.index('Cr')
......@@ -31,28 +75,32 @@ def test_quantifyWater():
Fitargs = {'ppmlim': [0.2, 5.2],
'method': 'MH',
'baseline_order': -1,
'baseline_order': 0,
'metab_groups': [0]}
res = fit_FSLModel(mrs, **Fitargs)
tissueFractions = {'GM': 0.6, 'WM': 0.4, 'CSF': 0.0}
TE = 0.03
TR = 20
T2dict = {'H2O_GM': 0.110,
'H2O_WM': 0.080,
'H2O_CSF': 2.55,
'METAB': 0.160}
q_info = quant.QuantificationInfo(
TE,
TR,
mrs.names,
mrs.centralFrequency / 1E6,
t2=T2dict)
q_info.set_fractions(tissueFractions)
res.calculateConcScaling(mrs,
referenceMetab=['Cr'],
waterRefFID=mrs.H2O,
tissueFractions=tissueFractions,
TE=TE,
T2=T2dict,
waterReferenceMetab='Cr',
wRefMetabProtons=5,
reflimits=(2, 5),
verbose=False)
q_info,
internal_reference=['Cr'],
verbose=True)
print(res.getConc(scaling='raw'))
print(res.getConc(scaling='internal'))
......@@ -60,4 +108,4 @@ def test_quantifyWater():
print(res.getConc(scaling='molarity'))
assert np.allclose(res.getConc(scaling='internal'), 1.0)
assert np.allclose(res.getConc(scaling='molarity'), 10.59, atol=1E-1)
assert np.allclose(res.getConc(scaling='molarity'), 10.78, atol=1E-1)
......@@ -37,6 +37,13 @@ H2O_MOLECULAR_MASS = 18.01528 # g/mol
H2O_MOLALITY = 55.51E3 # mmol/kg
H2O_PROTONS = 2
# Water referencing metabolites
# Define a list of sensible metabolites to use and
# the number of protons between the default limits of 2 and 5
WATER_SCALING_METAB = ['Cr', 'PCr', 'NAA']
WATER_SCALING_METAB_PROTONS = [5, 5, 3]
WATER_SCALING_DEFAULT_LIMITS = [(2, 5), (2, 5), (1.8, 2.2)]
# T1 parameters
# Derived from a survey of the literature. References listed below.
# Metabolite values derived from an average of NAA, Cr and Cho peaks.
......
......@@ -581,7 +581,8 @@ def getModelJac(model):
def getFittedModel(model, resParams, base_poly, metab_groups, mrs,
basisSelect=None, baselineOnly=False, noBaseline=False):
basisSelect=None, baselineOnly=False, noBaseline=False,
no_phase=False):
""" Return the predicted model given some fitting parameters
model (str) : Model string
......@@ -601,8 +602,15 @@ def getFittedModel(model, resParams, base_poly, metab_groups, mrs,
else:
bp = base_poly
if no_phase:
p = x2p(resParams, numBasis, numGroups)
p = p[0:-3] + (0.0, 0.0) + p[-1:]
params = p2x(*p)
else:
params = resParams
if basisSelect is None and not baselineOnly:
return forward(resParams,
return forward(params,
mrs.frequencyAxis,
mrs.timeAxis,
mrs.basis,
......@@ -610,7 +618,7 @@ def getFittedModel(model, resParams, base_poly, metab_groups, mrs,
metab_groups,
numGroups)
elif baselineOnly:
p = x2p(resParams, numBasis, numGroups)
p = x2p(params, numBasis, numGroups)
p = (np.zeros(numBasis),) + p[1:]
xx = p2x(*p)
return forward(xx,
......@@ -621,7 +629,7 @@ def getFittedModel(model, resParams, base_poly, metab_groups, mrs,
metab_groups,
numGroups)
elif basisSelect is not None:
p = x2p(resParams, numBasis, numGroups)
p = x2p(params, numBasis, numGroups)
tmp = np.zeros(numBasis)
basisIdx = mrs.names.index(basisSelect)
tmp[basisIdx] = p[0][basisIdx]
......
This diff is collapsed.
......@@ -77,12 +77,13 @@ def create_plotly_div(mrs, res):
divs['indiv'] = to_div(fig) # plotly.offline.plot(fig, output_type='div',include_plotlyjs='cdn')
# Quantification table
if (res.concScalings['molality'] is not None) and hasattr(res.concScalings['info']['q_info'], 'f_GM'):
Q = res.concScalings['info']['q_info']
if (res.concScalings['molality'] is not None) and hasattr(res.concScalings['quant_info'], 'f_GM'):
Q = res.concScalings['quant_info']
quant_df = pd.DataFrame()
quant_df['Tissue-water densities (g/cm^3)'] = [Q.d_GM, Q.d_WM, Q.d_CSF]
quant_df['Tissue volume fractions'] = [Q.f_GM, Q.f_WM, Q.f_CSF]
quant_df['Water T2 (ms)'] = [Q.T2['H2O_GM'] * 1000, Q.T2['H2O_WM'] * 1000, Q.T2['H2O_CSF'] * 1000]
quant_df['Water T2 (ms)'] = [Q.t2['H2O_GM'] * 1000, Q.t2['H2O_WM'] * 1000, Q.t2['H2O_CSF'] * 1000]
quant_df['Water T1 (s)'] = [Q.t1['H2O_GM'], Q.t1['H2O_WM'], Q.t1['H2O_CSF']]
quant_df.index = ['GM', 'WM', 'CSF']
quant_df.index.name = 'Tissue'
quant_df.reset_index(inplace=True)
......@@ -268,28 +269,38 @@ def create_report(mrs, res, filename, fidfile, basisfile, h2ofile, date, locatio
# Table of CSF,GM,WM
# Fractions
# T2 info (water)
# T1 info (water)
# T2 info (metab)
# T1 info (metab)
# Tissue water densities
# Relaxation corrected water
# Relaxation correction for metab
# Final scalings for molality & molarity
if res.concScalings['molality'] is not None:
Q = res.concScalings['info']['q_info']
relax_water_conc = Q.relaxationCorrWater_molar
metabRelaxCorr = Q.relaxationCorrMetab
Q = res.concScalings['quant_info']
relax_water_conc = Q.relax_corr_water_molar
metabRelaxCorr = Q.relax_corr_metab
if hasattr(res.concScalings['info']['q_info'], 'f_GM'):
if res.concScalings['quant_info'].f_GM is not None:
section = f"""
<h1><a name="quantification">Quantification information</a></h1>
<div style="width:70%">{divs['quant_table']}</div>
<table>
<tr>
<td class="titles">Metabolite T<sub>2</sub>:</td>
<td>{1000*Q.T2['METAB']} ms</td>
<td>{1000*Q.t2['METAB']} ms</td>
</tr>
<tr>
<td class="titles">Metabolite T<sub>1</sub>:</td>
<td>{Q.t1['METAB']} s</td>
</tr>
<tr>
<td class="titles">Sequence echo time (T<sub>E</sub>):</td>
<td>{1000*Q.TE} ms</td>
<td>{1000*Q.te} ms</td>
</tr>
<tr>
<td class="titles">Sequence repetition time (T<sub>R</sub>):</td>
<td>{Q.tr} s</td>
</tr>
<tr>
<td class="titles">T<sub>2</sub> relaxation corrected water concentration:</td>
......@@ -311,29 +322,32 @@ def create_report(mrs, res, filename, fidfile, basisfile, h2ofile, date, locatio
<hr>
"""
else:
if 'H2O' in Q.T2:
waterstr = f"""
<td class="titles">Water T<sub>2</sub></td>
<td> : {1000*Q.T2['H2O']} ms</td>
"""
else:
waterstr = """
<td class="titles">Water T<sub>2</sub></td>
<td> : Inf (no relaxation)</td>
"""
section = f"""
<h1><a name="quantification">Quantification information</a></h1>
<table>
<tr>
{waterstr}
<td class="titles">Water T<sub>2</sub></td>
<td> : {(1000*Q.t2['H2O_GM'] + 1000*Q.t2['H2O_WM']) / 2} ms</td>
</tr>
<tr>
<td class="titles">Metabolite T<sub>2</sub></td>
<td> : {1000*Q.T2['METAB']} ms</td>
<td class="titles">Water T<sub>1</sub></td>
<td> : {(1000*Q.t1['H2O_GM'] + 1000*Q.t1['H2O_WM']) / 2} ms</td>
</tr>
<tr>
<td class="titles">Metabolite T<sub>2</sub>:</td>
<td>{1000*Q.t2['METAB']} ms</td>
</tr>
<tr>
<td class="titles">Metabolite T<sub>1</sub>:</td>
<td>{Q.t1['METAB']} s</td>
</tr>
<tr>
<td class="titles">Sequence echo time (T<sub>E</sub>):</td>
<td>{1000*Q.te} ms</td>
</tr>
<tr>
<td class="titles">Sequence echo time (T<sub>E</sub>)</td>
<td> : {1000*Q.TE} ms</td>
<td class="titles">Sequence repetition time (T<sub>R</sub>):</td>
<td>{Q.tr} s</td>
</tr>
<tr>
<td class="titles">T<sub>2</sub> relaxation corrected water concentration</td>
......
......@@ -11,6 +11,7 @@ import fsl_mrs.utils.models as models
import fsl_mrs.utils.quantify as quant
import fsl_mrs.utils.qc as qc
from fsl_mrs.utils.misc import FIDToSpec, SpecToFID, calculate_lap_cov
import numpy as np
from copy import deepcopy
......@@ -96,61 +97,49 @@ class FitRes(object):
self.combine([['Cr', 'PCr']])
self.calculateConcScaling(mrs)
def calculateConcScaling(self, mrs,
referenceMetab=['Cr', 'PCr'],
waterRefFID=None,
tissueFractions=None,
TE=None,
T2='Default',
waterReferenceMetab='Cr',
wRefMetabProtons=5,
reflimits=(2, 5),
verbose=False,
Q=None,
add_scale=None):
self.intrefstr = '+'.join(referenceMetab)
self.referenceMetab = referenceMetab
self.waterReferenceMetab = waterReferenceMetab
self.waterReferenceMetabProtons = wRefMetabProtons
internalRefScaling = quant.quantifyInternal(referenceMetab, self.getConc(), self.metabs)
if (waterRefFID is not None) and ((TE is not None) or (Q is not None)):
refFID = self.predictedFID(mrs, mode=waterReferenceMetab, noBaseline=True)
if Q is None:
if T2 == 'Default':
Q = quant.loadDefaultQuantificationInfo(TE, tissueFractions, mrs.centralFrequency / 1E6)
else:
Q = quant.QuantificationInfo(TE, T2, tissueFractions)
if add_scale is not None:
Q.additionalCorr *= add_scale
molalityScaling, molarityScaling, quant_info = quant.quantifyWater(mrs,
waterRefFID,
refFID,
waterReferenceMetab,
self.getConc(),
self.metabs,
wRefMetabProtons,
Q,
reflimits=reflimits,
verbose=verbose)
def calculateConcScaling(self,
mrs,
quant_info=None,
internal_reference=['Cr', 'PCr'],
verbose=False):
"""Run calculation of internal and (if possible) water concentration scaling.
:param mrs: MRS object
:type mrs: fsl_mrs.core.mrs.MRS
:param quant_info: Quantification information for water scaling, defaults to None
:type quant_info: fsl_mrs.utils.quantify.QuantificationInfo, optional
:param internal_reference: Internal referencing emtabolite, defaults to ['Cr', 'PCr'] i.e. tCr
:type internal_reference: list, optional
:param verbose: Enable for verbose output, defaults to False
:type verbose: bool, optional