Commit 42f0efe8 authored by William Clarke's avatar William Clarke
Browse files

Tests for new quantification and results as well as fitting classes.

parent 2f356e83
......@@ -357,24 +357,6 @@ class MRS(object):
if metabs is not None:
metabs = [m for m in self.names if m not in metabs]
self.ignore(metabs)
def combine(self,metabs):
"""
Create combined metabolite groups
Parameters
----------
metabs: list of lists
"""
if metabs is not None:
self.metab_groups = []
for mgroups in metabs:
group = []
for m in mgroups:
group.append(m)
self.metab_groups.append(group)
def add_peak(self,ppm,name,gamma=0,sigma=0):
......
......@@ -49,7 +49,7 @@ def main():
# FITTING ARGUMENTS
fitting_args.add_argument('--algo',default='Newton',type=str,
help='algorithm [Newton (fast) or MH (slow)]')
help='algorithm [Newton (fast, default) or MH (slow)]')
fitting_args.add_argument('--ignore',type=str,nargs='+',metavar='METAB',
help='ignore certain metabolites [repeatable]')
fitting_args.add_argument('--keep',type=str,nargs='+',metavar='METAB',
......@@ -61,7 +61,7 @@ def main():
fitting_args.add_argument('--h2o',default=None,type=str,metavar='H2O',
help='input .H2O file for quantification')
fitting_args.add_argument('--baseline_order',default=2,type=int,metavar=('ORDER'),
help='order of baseline polynomial (default=2)')
help='order of baseline polynomial (default=2,-1 disables)')
fitting_args.add_argument('--metab_groups',default=0,nargs='+',type=str_or_int_arg,
help="metabolite groups: list of groups or list of names for indept groups.")
fitting_args.add_argument('--add_MM',type=bool,
......@@ -71,6 +71,10 @@ def main():
# ADDITONAL OPTIONAL ARGUMENTS
optional.add_argument('--t1',type=str,default=None,metavar='IMAGE',
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('--tissue_fractions',type=float,nargs=3,default=[0.5,0.5,0.0],metavar='GM WM CSF',
help='Fractional tissue volumes for WM, GM, CSF. Defaults to 0.5, 0.5, 0.0.')
optional.add_argument('--central_frequency',default=None,type=float,
help='central frequency in Hz')
optional.add_argument('--dwell_time',default=None,type=float,
......@@ -217,8 +221,7 @@ def main():
# Keep/Ignore/Combine metabolites
mrs.keep(args.keep)
mrs.ignore(args.ignore)
mrs.combine(args.combine)
mrs.ignore(args.ignore)
# Do the fitting here
if args.verbose:
......@@ -226,10 +229,7 @@ def main():
print(' Algorithm = [{}]\n'.format(args.algo))
start = time.time()
# Do the fitting
if args.verbose:
print(' --- Run fitting --- ')
ppmlim=args.ppmlim
if ppmlim is not None:
ppmlim=(ppmlim[0],ppmlim[1])
......@@ -274,21 +274,42 @@ def main():
res = fitting.fit_FSLModel(mrs,**Fitargs)
# Quantification
# Echo time
if args.TE is not None:
echotime = args.TE*1E-3
elif 'TE' in basisheader['meta']:
echotime = basisheader['meta']['TE']
if echotime > 1.0: # Assume in ms.
echotime *= 1E-3
else:
echotime = None
# Internal and Water quantification if requested
if (mrs.H2O is None) or (echotime is None):
res.calculateConcScaling(mrs,referenceMetab=['Cr','PCr'])
else:
tfrac = {'WM':args.tissue_fractions[0],'GM':args.tissue_fractions[1],'CSF':args.tissue_fractions[2]}
res.calculateConcScaling(mrs,
referenceMetab=['Cr','PCr'],
waterRefFID=mrs.H2O,
tissueFractions=tfrac,
TE=echotime,
verbose=args.verbose)
# Combine metabolites.
if args.combine is not None:
res.combine(args.combine)
stop = time.time()
# Report on the fitting
if args.verbose:
print(' Fitting lasted : {:.3f} secs.\n'.format(stop-start))
#print(' Fitted concentrations : {}\n'.format(mrs.concentrations))
# Save output files
if args.verbose:
print('--->> Saving output files to {}\n'.format(args.output))
res.to_file(filename=os.path.join(args.output,'results_table.csv'),mrs=mrs,what='concentrations')
res.to_file(filename=os.path.join(args.output,'results_table.csv'),what='concentrations')
res.to_file(filename=os.path.join(args.output,'qc.csv'),what='qc')
res.to_file(filename=os.path.join(args.output,'all_parameters.csv'),what='parameters')
......
# Test the main svs fitting script
# Imports
import subprocess
import pytest
import os.path as op
# Files
testsPath = op.dirname(__file__)
data = {'metab':op.join(testsPath,'testdata/fsl_mrs/metab.nii'),
'water':op.join(testsPath,'testdata/fsl_mrs/water.nii'),
'basis':op.join(testsPath,'testdata/fsl_mrs/steam_basis')}
def test_fsl_mrs(tmp_path):
retcode = subprocess.check_call(['fsl_mrs',
'--data', data['metab'],
'--basis', data['basis'],
'--output', tmp_path,
'--h2o', data['water'],
'--TE', '11',
'--tissue_fractions', '0.45', '0.45', '0.1',
'--overwrite',
'--combine', 'Cr', 'PCr',
'--report'])
assert op.exists(op.join(tmp_path,'report.html'))
assert op.exists(op.join(tmp_path,'all_parameters.csv'))
assert op.exists(op.join(tmp_path,'qc.csv'))
assert op.exists(op.join(tmp_path,'results_table.csv'))
assert op.exists(op.join(tmp_path,'options.txt'))
\ No newline at end of file
from fsl_mrs.utils.synthetic import syntheticFID
from fsl_mrs.core import MRS
from fsl_mrs.utils.fitting import fit_FSLModel
from pytest import fixture
import numpy as np
@fixture
def data():
noiseCov = 0.01
amplitude = np.asarray([0.5,0.5,1.0])*10
chemshift = np.asarray([3.0,3.05,2.0])-4.65
lw = [10,10,10]
phases = [0,0,0]
g = [0,0,0]
basisNames = ['Cr','PCr','NAA']
basisFIDs = []
for idx,_ in enumerate(amplitude):
tmp,basisHdr = syntheticFID(noisecovariance=[[0.0]],
chemicalshift=[chemshift[idx]],
amplitude=[1.0],
linewidth=[lw[idx]/5],
phase=[phases[idx]],
g=[g[idx]])
basisFIDs.append(tmp[0])
basisFIDs = np.asarray(basisFIDs)
synFID,synHdr = syntheticFID(noisecovariance=[[noiseCov]],
chemicalshift=chemshift,
amplitude=amplitude,
linewidth=lw,
phase=phases,
g=g)
synMRS = MRS(FID =synFID[0],header=synHdr,basis =basisFIDs,basis_hdr=basisHdr,names=basisNames)
return synMRS,amplitude
def test_fit_FSLModel_Newton(data):
mrs = data[0]
amplitudes = data[1]
metab_groups = [0]*mrs.numBasis
Fitargs = {'ppmlim':[0.2,4.2],
'method':'Newton','baseline_order':-1,
'metab_groups':metab_groups}
res = fit_FSLModel(mrs,**Fitargs)
fittedconcs = res.getConc()
fittedRelconcs = res.getConc(scaling='internal')
assert np.allclose(fittedconcs,amplitudes,atol=1E-1)
assert np.allclose(fittedRelconcs,amplitudes/(amplitudes[0]+amplitudes[1]),atol=1E-1)
def test_fit_FSLModel_MH(data):
mrs = data[0]
amplitudes = data[1]
metab_groups = [0]*mrs.numBasis
Fitargs = {'ppmlim':[0.2,4.2],
'method':'MH','baseline_order':-1,
'metab_groups':metab_groups,
'MHSamples':100}
res = fit_FSLModel(mrs,**Fitargs)
fittedconcs = res.getConc()
fittedRelconcs = res.getConc(scaling='internal')
assert np.allclose(fittedconcs,amplitudes,atol=1E-1)
assert np.allclose(fittedRelconcs,amplitudes/(amplitudes[0]+amplitudes[1]),atol=1E-1)
\ No newline at end of file
......@@ -42,10 +42,10 @@ def test_calcQC():
'metab_groups':[0]}
res = fit_FSLModel(synMRS_basis,**Fitargs)
fwhm_test,snrSpec_test,snrPeaks_test = calcQC(synMRS_basis,res,ppmlim=[2.65,6.65])
print(f'Measured FWHM: {fwhm_test[0]:0.1f}')
print(f'Measured spec SNR: {snrSpec_test:0.1f}')
print(f'Measured peak SNR: {snrPeaks_test[0]:0.1f}')
assert np.isclose(fwhm_test,trueLW,atol=1E0)
assert np.isclose(snrSpec_test,SNR_noApod,atol=1E1)
assert np.isclose(snrPeaks_test,SNR,atol=2E1)
fwhm_test,SNRObj = calcQC(synMRS_basis,res,ppmlim=[2.65,6.65])
print(f'Measured FWHM: {fwhm_test.mean().to_numpy()[0]:0.1f}')
print(f'Measured spec SNR: {SNRObj.spectrum:0.1f}')
print(f'Measured peak SNR: {SNRObj.peaks.mean().to_numpy()[0]:0.1f}')
assert np.isclose(fwhm_test.mean().to_numpy(),trueLW,atol=1E0)
assert np.isclose(SNRObj.spectrum,SNR_noApod,atol=1E1)
assert np.isclose(SNRObj.peaks.mean().to_numpy(),SNR,atol=2E1)
# Tests for the quantify module.
# Utilise the independently constructed MRS fitting challenge data to test against
import os.path as op
from fsl_mrs.core import MRS
import fsl_mrs.utils.mrs_io as mrsio
from fsl_mrs.utils.fitting import fit_FSLModel
import numpy as np
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_quantifyWater():
basis,names,headerb = mrsio.read_basis(basisfile)
crIndex = names.index('Cr')
data,header = mrsio.read_FID(metabfile)
dataw,headerw = mrsio.read_FID(h2ofile)
mrs = MRS(FID=data,header=header,basis=basis[:,crIndex],names=['Cr'],basis_hdr=headerb[crIndex],H2O=dataw)
mrs.check_FID(repare=True)
mrs.check_Basis(repare=True)
Fitargs = {'ppmlim':[0.2,5.2],
'method':'MH','baseline_order':-1,
'metab_groups':[0]}
res = fit_FSLModel(mrs,**Fitargs)
tissueFractions = {'GM':0.6,'WM':0.4,'CSF':0.0}
TE = 0.03
T2dict = {'H2O_GM':0.110,
'H2O_WM':0.080,
'H2O_CSF':2.55,
'METAB':0.160}
res.calculateConcScaling(mrs,
referenceMetab=['Cr'],
waterRefFID=mrs.H2O,
tissueFractions=tissueFractions,
TE=TE,
T2= T2dict,
waterReferenceMetab='Cr',
wRefMetabProtons=5,
reflimits=(2,5),
verbose=False)
print(res.getConc(scaling='raw'))
print(res.getConc(scaling='internal'))
print(res.getConc(scaling='molality'))
print(res.getConc(scaling='molarity'))
assert np.allclose(res.getConc(scaling='internal'),1.0)
assert np.allclose(res.getConc(scaling='molarity'),10.72,atol=1E-1)
\ No newline at end of file
# Test features of the results class
# Imports
from fsl_mrs.utils.synthetic import syntheticFID
from fsl_mrs.core import MRS
from fsl_mrs.utils.fitting import fit_FSLModel
from pytest import fixture
import numpy as np
# Set up some synthetic data to use
@fixture(scope='module')
def data():
noiseCov = 0.01
amplitude = np.asarray([0.5,0.5,1.0])*10
chemshift = np.asarray([3.0,3.05,2.0])-4.65
lw = [10,10,10]
phases = [0,0,0]
g = [0,0,0]
basisNames = ['Cr','PCr','NAA']
begintime = 0.0001
basisFIDs = []
for idx,_ in enumerate(amplitude):
tmp,basisHdr = syntheticFID(noisecovariance=[[0.0]],
chemicalshift=[chemshift[idx]+0.1],
amplitude=[1.0],
linewidth=[lw[idx]/5],
phase=[phases[idx]],
g=[g[idx]],
begintime=begintime)
basisFIDs.append(tmp[0])
basisFIDs = np.asarray(basisFIDs)
synFID,synHdr = syntheticFID(noisecovariance=[[noiseCov]],
chemicalshift=chemshift,
amplitude=amplitude,
linewidth=lw,
phase=phases,
g=g)
synMRS = MRS(FID =synFID[0],header=synHdr,basis =basisFIDs,basis_hdr=basisHdr,names=basisNames)
metab_groups = [0]*synMRS.numBasis
Fitargs = {'ppmlim':[0.2,4.2],
'method':'MH','baseline_order':-1,
'metab_groups':metab_groups,
'MHSamples':100}
res = fit_FSLModel(synMRS,**Fitargs)
return res,amplitude
def test_peakcombination(data):
res = data[0]
amplitudes = data[1]
res.combine([['Cr','PCr']])
fittedconcs = res.getConc()
fittedRelconcs = res.getConc(scaling='internal')
amplitudes = np.append(amplitudes,amplitudes[0]+amplitudes[1])
assert 'Cr+PCr' in res.metabs
assert np.allclose(fittedconcs,amplitudes,atol=1E-1)
assert np.allclose(fittedRelconcs,amplitudes/(amplitudes[0]+amplitudes[1]),atol=1E-1)
def test_units(data):
res = data[0]
# Phase
p0,p1 = res.getPhaseParams(phi0='degrees',phi1='seconds')
assert np.isclose(p0,0.0,atol=1E-1)
assert np.isclose(p1,0.0001,atol=1E-5)
# Shift
shift = res.getShiftParams(units='ppm')
shift_hz = res.getShiftParams(units='Hz')
assert np.isclose(shift,0.1,atol=1E-2)
assert np.isclose(shift_hz,0.1*123.0,atol=1E-1)
# Linewidth
lw = res.getLineShapeParams(units='Hz')
lw_ppm = res.getLineShapeParams(units='ppm')
assert np.isclose(lw,8.0,atol=1E-1) #10-2
assert np.isclose(lw_ppm,8.0/123.0,atol=1E-1)
def test_qcOutput(data):
res = data[0]
SNR,FWHM = res.getQCParams()
assert np.allclose(FWHM,10.0,atol=1E0)
assert SNR.size == 3
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
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