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

QC code. Not yet integrated, but tested. Minor speed improvements by using...

QC code. Not yet integrated, but tested. Minor speed improvements by using scipy in misc.FIDToSpec. Some cleanups in fitting and models code.
parent 75dadedf
......@@ -63,11 +63,14 @@ class MRS(object):
# Set Basis info
if basis is not None:
self.basis = basis.copy()
# Handle single basis spectra
if self.basis.ndim==1:
self.basis = self.basis[:,np.newaxis]
# Assume that there will always be more timepoints than basis spectra.
if self.basis.shape[0] < self.basis.shape[1]:
self.basis = self.basis.T
self.numBasis = basis.shape[1]
self.numBasisPoints = basis.shape[0]
self.numBasis = self.basis.shape[1]
self.numBasisPoints = self.basis.shape[0]
if (names is not None) and (basis_hdr is not None):
self.names = names.copy()
......
from fsl_mrs.utils.synthetic import syntheticFID
from fsl_mrs.utils.qc import specApodise,calcQC
from fsl_mrs.utils.fitting import fit_FSLModel
from fsl_mrs.core import MRS
import numpy as np
def test_calcQC():
# Syntetic data
synFID,synHdr = syntheticFID(noisecovariance=[[0.1]],points=2*2048,chemicalshift=[0],amplitude=[6.0],linewidth=[10])
synFIDNoise,synHdrNoise = syntheticFID(noisecovariance=[[0.1]],points=2*2048,chemicalshift=[0],amplitude=[0],linewidth=[10])
basisFID,basisHdr = syntheticFID(noisecovariance=[[0.0]],points=2*2048,chemicalshift=[0],amplitude=[0.1],linewidth=[2])
synMRS = MRS(FID =synFID[0],header=synHdr)
synMRSNoise = MRS(FID =synFIDNoise[0],header=synHdrNoise)
synMRSNoNoise = MRS(FID =synHdr['noiseless'],header=synHdr)
synMRS_basis = MRS(FID =synFID[0],header=synHdr,basis =basisFID[0] ,basis_hdr=basisHdr,names=['Peak1'])
truenoiseSD = np.sqrt(synHdrNoise['cov'][0,0])
pureNoiseMeasured = np.std(synMRSNoise.getSpectrum())
realnoise = np.std(np.real(synMRSNoise.getSpectrum()))
imagNoise = np.std(np.imag(synMRSNoise.getSpectrum()))
print(f'True cmplx noise = {truenoiseSD:0.3f}, pure noise measured = {pureNoiseMeasured:0.3f} (real/imag = {realnoise:0.3f}/{imagNoise:0.3f})')
# Calc SNR without apodisation from the no noise and pure noise spectra
truePeakHeight = np.max(np.real(synMRSNoNoise.getSpectrum()))
SNR_noApod = truePeakHeight/pureNoiseMeasured
print(f'SNR no apod: {SNR_noApod:0.1f} ({truePeakHeight:0.2e}/{pureNoiseMeasured:0.2e})')
# Calc SNR with apodisation from the no noise and pure noise spectra
trueLW = synHdr['inputopts']['linewidth'][0]
trueApodSpec_Noise = specApodise(synMRSNoise,trueLW)
apodNoise = np.std(trueApodSpec_Noise)
trueApodSpec_noNoise = specApodise(synMRSNoNoise,trueLW)
peakHeigtApod = np.max(np.real(trueApodSpec_noNoise))
SNR = peakHeigtApod/apodNoise
print(f'SNR w. apod: {SNR:0.1f} ({peakHeigtApod:0.2e}/{apodNoise:0.2e})')
metab_groups = [0]
Fitargs = {'ppmlim':[2.65,6.65],
'method':'Newton','baseline_order':-1,
'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)
......@@ -102,15 +102,15 @@ class FitRes(object):
import pandas as pd
df = pd.DataFrame()
if what is 'concentrations':
if what == 'concentrations':
df['Metab'] = mrs.names
df['mMol/kg'] = self.conc_h2o
df['%CRLB'] = self.perc_SD[:mrs.numBasis]
df['/Cr'] = self.conc_cr
elif what is 'qc':
elif what == 'qc':
df['Measure'] = ['SNR']
df['Value'] = [self.snr]
elif what is 'parameters':
elif what == 'parameters':
df['Parameter'] = self.params_names
df['Value'] = self.params
......@@ -360,7 +360,7 @@ def prepare_baseline_regressor(mrs,baseline_order,ppmlim):
return B
def get_bounds(num_basis,num_metab_groups,baseline_order,model,method):
def get_bounds(num_basis,num_metab_groups,baseline_order,model,method,disableBaseline=False):
if method == 'Newton':
# conc
bnds = [(0,None)]*num_basis
......@@ -373,7 +373,10 @@ def get_bounds(num_basis,num_metab_groups,baseline_order,model,method):
bnds.extend([(None,None)]*2)
# baseline
n = (1+baseline_order)*2
bnds.extend([(None,None)]*n)
if disableBaseline:
bnds.extend([(0.0,0.0)]*n)
else:
bnds.extend([(None,None)]*n)
return bnds
elif method == 'MH':
......@@ -395,8 +398,12 @@ def get_bounds(num_basis,num_metab_groups,baseline_order,model,method):
UB.extend([MAX]*2)
# baseline
n = (1+baseline_order)*2
LB.extend([MIN]*n)
UB.extend([MAX]*n)
if disableBaseline:
LB.extend([0.0]*n)
UB.extend([0.0]*n)
else:
LB.extend([MIN]*n)
UB.extend([MAX]*n)
return LB,UB
......@@ -440,7 +447,7 @@ def fit_FSLModel(mrs,
"""
A simplified version of LCModel
"""
err_func,grad_func,forward,x2p,_ = model.getModelFunctions(model)
err_func,grad_func,forward,x2p,_ = models.getModelFunctions(model)
if model == 'lorentzian':
init_func = init_FSLModel # initilisation of params
elif model == 'voigt':
......@@ -455,6 +462,12 @@ def fit_FSLModel(mrs,
# shorter names for some of the useful stuff
freq,time,basis=mrs.frequencyAxis,mrs.timeAxis,mrs.basis
# Handle completely disabling basline
if baseline_order < 0:
baseline_order = 0 # Generate one order of baseline parameters
disableBaseline = True # But diable by setting bounds to 0
else:
disableBaseline = False
# Results object
results = FitRes(model)
......@@ -476,7 +489,7 @@ def fit_FSLModel(mrs,
# Fitting
if method == 'Newton':
# Bounds
bounds = get_bounds(mrs.numBasis,g,baseline_order,model,method)
bounds = get_bounds(mrs.numBasis,g,baseline_order,model,method,disableBaseline=disableBaseline)
res = minimize(err_func, x0, args=constants,
method='TNC',jac=grad_func,bounds=bounds)
# collect results
......@@ -526,7 +539,7 @@ def fit_FSLModel(mrs,
results.pred = misc.SpecToFID(results.pred_spec) # predict FID not Spec
# baseline
baseline = model.getFittedModel(model,results.params,results.base_poly,metab_groups,mrs,baselineOnly= True)
baseline = models.getFittedModel(model,results.params,results.base_poly,metab_groups,mrs,baselineOnly= True)
baseline = misc.SpecToFID(baseline)
results.baseline = baseline
p = x2p(results.params,mrs.numBasis,g)
......@@ -552,8 +565,8 @@ def fit_FSLModel(mrs,
results.mcmc_cor = np.ma.corrcoef(results.mcmc_samples.T)
results.mcmc_var = np.var(results.mcmc_samples,axis=0)
results.perc_SD = np.sqrt(results.crlb) / results.params*100
with np.errstate(divide='ignore', invalid='ignore'):
results.perc_SD = np.sqrt(results.crlb) / results.params*100
results.perc_SD[results.perc_SD>999] = 999 # Like LCModel :)
results.perc_SD[np.isnan(results.perc_SD)] = 999
......@@ -575,9 +588,10 @@ def fit_FSLModel(mrs,
# nuisance parameters in sensible units
if model == 'lorentzian':
con,gamma,eps,phi0,phi1,b = models.FSLModel_x2param(results.params,mrs.numBasis,g)
results.inv_gamma_sec = 1/gamma
con,gamma,eps,phi0,phi1,b = models.FSLModel_x2param(results.params,mrs.numBasis,g)
results.gamma = gamma
results.gamma_hz = gamma/np.pi
results.inv_gamma_sec = 1/results.gamma_hz
elif model == 'voigt':
con,gamma,sigma,eps,phi0,phi1,b = models.FSLModel_x2param_Voigt(results.params,mrs.numBasis,g)
results.inv_gamma_sec = 1/gamma
......@@ -589,7 +603,8 @@ def fit_FSLModel(mrs,
results.phi1_deg_per_ppm = phi1*np.pi/180.0 * mrs.centralFrequency / 1E6
results.eps = eps
results.eps_ppm = eps / mrs.centralFrequency * 1E6
results.b_norm = b/b[0]
with np.errstate(divide='ignore', invalid='ignore'):
results.b_norm = b/b[0]
# QC parameters (like LCModel)
......
......@@ -8,6 +8,7 @@
# SHBASECOPYRIGHT
import numpy as np
import scipy.fft
from scipy.signal import butter, lfilter
from scipy.interpolate import interp1d
import itertools as it
......@@ -45,17 +46,11 @@ def FIDToSpec(FID,axis=0):
Returns:
x (np.array) : array of spectra
"""
# By convention the first point of the fid is special cased
def scaleFID(x):
x[0,...] *= 0.5
return x
# move indicated axis into first dimension
# copy so we don't modify the first fid point
FID = np.moveaxis(FID,axis,0).copy()
x = np.fft.fftshift(np.fft.fft(scaleFID(FID),axis=0),axes=0)/FID.shape[0]
x = np.moveaxis(x,0,axis)
return x
# By convention the first point of the fid is special cased
FID[0] *=0.5
out = scipy.fft.fftshift(scipy.fft.fft(FID,axis=axis,norm='ortho'),axes=axis)
FID[0] *=2
return out
def SpecToFID(spec,axis=0):
""" Convert spectrum to FID
......@@ -67,14 +62,10 @@ def SpecToFID(spec,axis=0):
Returns:
x (np.array) : array of FIDs
"""
def scaleFID(x):
x[0] *= 2
return x
spec = np.moveaxis(spec,axis,0).copy()
x = scaleFID(np.fft.ifft(np.fft.ifftshift(spec,axes=0),axis=0)*spec.shape[0])
x = np.moveaxis(x,0,axis)
return x
"""
fid = scipy.fft.ifft(scipy.fft.ifftshift(spec,axes=axis),axis=axis,norm='ortho')
fid[0] *= 2
return fid
def calculateAxes(bandwidth,centralFrequency,points):
dwellTime = 1/bandwidth
......
......@@ -214,7 +214,8 @@ def FSLModel_forward(x,nu,t,m,B,G,g):
E = np.zeros((m.shape[0],g),dtype=np.complex)
for gg in range(g):
E[:,gg] = np.exp(-(1j*eps[gg]+gamma[gg])*t).flatten()
# E = np.exp(-(1j*eps+gamma)*t) # THis is actually slower! But maybe more optimisable longterm with numexpr or numba
tmp = np.zeros(m.shape,dtype=np.complex)
for i,gg in enumerate(G):
tmp[:,i] = m[:,i]*E[:,gg]
......
......@@ -89,8 +89,9 @@ def plot_fit(mrs,pred=None,ppmlim=(0.40,4.2),out=None,baseline=None,proj='real')
axis = mrs.ppmAxisShift
first = np.argmin(np.abs(axis[0:int(mrs.numPoints/2)]-ppmlim[0]))
last = np.argmin(np.abs(axis[0:int(mrs.numPoints/2)]-ppmlim[1]))
first,last = mrs.ppmlim_to_range(ppmlim=ppmlim,shift=True)
# first = np.argmin(np.abs(axis[0:int(mrs.numPoints/2)]-ppmlim[0]))
# last = np.argmin(np.abs(axis[0:int(mrs.numPoints/2)]-ppmlim[1]))
# turn to real numbers
......
......@@ -33,6 +33,9 @@ def calcQC(mrs,res,ppmlim=(0.2,4.2)):
return fwhm,snrSpec,snrPeaks
class NoiseNotFoundError(ValueError):
pass
def calcQCOnResults(mrs,res,resparams,ppmlim):
""" Calculate QC metrics on single instance of fitting results
......@@ -45,7 +48,16 @@ def calcQCOnResults(mrs,res,resparams,ppmlim):
for basemrs in basisMRS:
combinedSpectrum += np.real(basemrs.getSpectrum())
normCombSpec = combinedSpectrum/np.max(combinedSpectrum)
noiseRegion = np.abs(normCombSpec)<0.015
noiseThreshold = 0.001
noiseRegion = np.abs(normCombSpec)<noiseThreshold
# print(np.sum(noiseRegion))
while np.sum(noiseRegion)<100:
if noiseThreshold>0.1:
raise NoiseNotFoundError(f'Unable to identify suitable noise area. Only {np.sum(noiseRegion)} points of {normCombSpec.size} found. Minimum of 100 needed.')
noiseThreshold += 0.001
noiseRegion = np.abs(normCombSpec)<noiseThreshold
# print(np.sum(noiseRegion))
# Noise region OS masks
noiseOSMask = detectOS(mrs,noiseRegion)
combinedMask = noiseRegion&noiseOSMask
......@@ -60,7 +72,7 @@ def calcQCOnResults(mrs,res,resparams,ppmlim):
basisSNR = []
for basemrs in basisMRS:
#FWHM
fwhm_curr,_,_ = idPeaksCalcFWHM(basemrs,estimatedFWHM=res.eps[0],ppmlim=ppmlim)
fwhm_curr,_,_ = idPeaksCalcFWHM(basemrs,estimatedFWHM=res.gamma_hz[0],ppmlim=ppmlim)
fwhm.append(fwhm_curr)
#Basis SNR
......@@ -118,6 +130,10 @@ def matchedFilterSNR(mrs,basismrs,lw,noisemask,ppmlim):
f,l = mrs.ppmlim_to_range(ppmlim=ppmlim)
peakHeight = np.max(np.real(apodbasis[f:l]))
currSNR= peakHeight/currNoise
# import matplotlib.pyplot as plt
# plt.plot(np.real(apodSpec))
# plt.plot(noisemask*np.max(np.real(apodSpec)))
# plt.show()
# print(f'SNR: {currSNR:0.1f} ({peakHeight:0.2e}/{currNoise:0.2e}), LW: {lw:0.1f}')
return peakHeight/currNoise
......
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