91 lines
4.0 KiB
Python
91 lines
4.0 KiB
Python
# Synthetise new data and analyse it directly to return fit parameters
|
|
|
|
import numpy as np
|
|
import panda as pd
|
|
|
|
from fsl_mrs.core import MRS
|
|
from fsl_mrs.utils import plotting, mrs_io
|
|
from fsl_mrs.utils.synthetic import synthetic_from_basis as synth
|
|
from fsl_mrs.utils.misc import parse_metab_groups
|
|
from fsl_mrs.utils.fitting import fit_FSLModel
|
|
|
|
|
|
def synth_and_ana(noise_cov,
|
|
fit_parameters, fit_snr, fit_sdnoise, fit_varnoise,
|
|
df_parameter_synth, basis_path, output_path):
|
|
"""Synthetise spectra with given noise-covariance, analyse the data and
|
|
return table of fitting parameters with a fit-plot.
|
|
We ignore Gly. Independent scaling of Macro Molecules."""
|
|
|
|
# Load input
|
|
synth_parameter = df_parameter_synth.mean().to_dict()
|
|
basis, names, Bheader = mrs_io.read_basis(basis_path)
|
|
|
|
# Set the adjustment lists from our input
|
|
broadening = [(synth_parameter["gamma_0"], synth_parameter["sigma_0"]),
|
|
(synth_parameter["gamma_1"], synth_parameter["sigma_1"])]
|
|
shifting = [(synth_parameter["eps_0"]), (synth_parameter["eps_1"])]
|
|
baseline = [synth_parameter["B_real_0"], synth_parameter["B_imag_0"],
|
|
synth_parameter["B_real_1"], synth_parameter["B_imag_1"],
|
|
synth_parameter["B_real_2"], synth_parameter["B_imag_2"]]
|
|
coilphase = [synth_parameter["Phi0"]]
|
|
|
|
# Generate synthetic data
|
|
fidS, headerS, concentrationsS = synth.syntheticFromBasisFile(basis_path,
|
|
concentrations=synth_parameter,
|
|
ignore=['Gly'], ind_scaling=['mm'],
|
|
metab_groups='mm', broadening=broadening, shifting=shifting,
|
|
# correct for complex noise
|
|
noisecovariance=np.divide(
|
|
noise_cov, 2),
|
|
# CAVE: baseline chosen manually
|
|
baseline=baseline, baseline_ppm=(
|
|
.2, 4.2),
|
|
coilphase=coilphase,
|
|
bandwidth=6000)
|
|
# Create mrs object for further use
|
|
mrsA = MRS(FID=fidS, header=headerS, basis=basis,
|
|
basis_hdr=Bheader[0], names=names)
|
|
|
|
mrsA.ignore(['Gly'])
|
|
mrsA.processForFitting(ind_scaling=['mm'])
|
|
|
|
# Scale it to Input data
|
|
# rescaled_FID, __ = rescale_FID(mrsA.FID, 1/mrsA.scaling["FID"]*100)
|
|
mrsA.set_FID(fidS)
|
|
|
|
metab_groups = parse_metab_groups(mrsA, 'mm')
|
|
# Use Voigt line broadening, fit between .2 and 4.2,
|
|
# with a 2nd order polynomial baseline
|
|
FitArgs = {
|
|
'model': 'voigt',
|
|
'metab_groups': metab_groups,
|
|
'ppmlim': (.2, 4.2),
|
|
'baseline_order': 2}
|
|
res = fit_FSLModel(mrsA, **FitArgs)
|
|
|
|
# Combine highly correlated metabolites
|
|
combinationList = [['Glu', 'Gln'],
|
|
['GPC', 'PCho'],
|
|
['Cr', 'PCr'],
|
|
['Glc', 'Tau'],
|
|
["NAA", "NAAG"]]
|
|
res.combine(combinationList)
|
|
|
|
# Store parameters from fit
|
|
fit_parameters.append(res.fitResults)
|
|
fit_snr.append(res.SNR.spectrum)
|
|
# Old version of noise
|
|
# noise_sd=np.max(np.real(mrsA.get_spec(ppmlim=(.2,4.2))))/res.SNR.spectrum
|
|
noise_sd = np.std(mrsA.FID[1000:1600])
|
|
fit_sdnoise.append(noise_sd)
|
|
noise_var = np.var(mrsA.FID[1000:1600])
|
|
fit_varnoise.append(noise_var)
|
|
|
|
df_params = pd.concat(fit_parameters, ignore_index=True)
|
|
df_params['SNR'] = fit_snr
|
|
df_params['noise_sd'] = fit_sdnoise
|
|
df_params['noise_var'] = fit_varnoise
|
|
|
|
return df_params
|