Commit 9b400fcf authored by William Clarke's avatar William Clarke
Browse files

Modification to allow quantification without tissue fractions. Not happy with...

Modification to allow quantification without tissue fractions. Not happy with this quantification package yet but it will do.
parent a8a70b26
......@@ -74,7 +74,7 @@ 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('--tissue_fractions',type=float,nargs=3,default=[0.5,0.5,0.0],metavar='GM WM CSF',
optional.add_argument('--tissue_fractions',type=float,nargs=3,default=None,metavar='GM WM CSF',
help='Fractional tissue volumes for WM, GM, CSF. Defaults to 0.5, 0.5, 0.0.')
optional.add_argument('--internal_ref',type=str,default=['Cr','PCr'],nargs='+',
help='Metabolite(s) used as an internal reference. Defaults to tCr (Cr+PCr).')
......@@ -298,6 +298,13 @@ def main():
# Internal and Water quantification if requested
if (mrs.H2O is None) or (echotime is None):
res.calculateConcScaling(mrs,referenceMetab=args.internal_ref)
elif args.tissue_fractions is None:
res.calculateConcScaling(mrs,
referenceMetab=args.internal_ref,
waterRefFID=mrs.H2O,
tissueFractions=None,
TE=echotime,
verbose=args.verbose)
else:
tfrac = {'WM':args.tissue_fractions[0],'GM':args.tissue_fractions[1],'CSF':args.tissue_fractions[2]}
res.calculateConcScaling(mrs,
......
......@@ -19,8 +19,8 @@ def test_quantifyWater():
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)
mrs.check_FID(repair=True)
mrs.check_Basis(repair=True)
Fitargs = {'ppmlim':[0.2,5.2],
'method':'MH','baseline_order':-1,
......
......@@ -51,7 +51,7 @@ def quantifyWater(mrs,H2OFID,refFID,referenceName,concentrations,metaboliteNames
referenceName (str): Name of reference metabolite
concentrations (np.array): All metabolite raw concentrations
metaboliteNames (list:str): All metabolite names
refProtons (): Number of protons contributing to reference spectrum between reflimits
refProtons (int): Number of protons contributing to reference spectrum between reflimits
Q (QuantificationInfo object): Contains tissue information
reflimits (tuple:float): Limits of integration for reference metabolite
verbose (bool): Verbose output
......@@ -67,17 +67,20 @@ def quantifyWater(mrs,H2OFID,refFID,referenceName,concentrations,metaboliteNames
# Calculate concnetration scalings
# EQ 4 and 6 in https://doi.org/10.1002/nbm.4257
conc_molal = (SMObs *(Q.f_GM*Q.R_H2O_GM + Q.f_WM*Q.R_H2O_WM + Q.f_CSF*Q.R_H2O_CSF)/(SH2OObs*(1-Q.f_CSF)*Q.R_M)) \
* (H2O_PROTONS/refProtons)\
* H2O_MOLALITY
# conc_molal = (SMObs *(Q.f_GM*Q.R_H2O_GM + Q.f_WM*Q.R_H2O_WM + Q.f_CSF*Q.R_H2O_CSF)/(SH2OObs*(1-Q.f_CSF)*Q.R_M)) \
# * (H2O_PROTONS/refProtons)\
# * H2O_MOLALITY
conc_molar = (SMObs *(Q.f_GM*Q.d_GM*Q.R_H2O_GM + Q.f_WM*Q.d_WM*Q.R_H2O_WM + Q.f_CSF*Q.d_CSF*Q.R_H2O_CSF)/(SH2OObs*(1-Q.f_CSF)*Q.R_M))\
* (H2O_PROTONS/refProtons)\
* H2O_MOLALITY
# conc_molar = (SMObs *(Q.f_GM*Q.d_GM*Q.R_H2O_GM + Q.f_WM*Q.d_WM*Q.R_H2O_WM + Q.f_CSF*Q.d_CSF*Q.R_H2O_CSF)/(SH2OObs*(1-Q.f_CSF)*Q.R_M))\
# * (H2O_PROTONS/refProtons)\
# * H2O_MOLALITY
conc_molal = (SMObs/SH2OObs) * (H2O_PROTONS/refProtons) * Q.relaxationCorrWater_molal * Q.additionalCorr * Q.relaxationCorrMetab
conc_molar = (SMObs/SH2OObs) * (H2O_PROTONS/refProtons) * Q.relaxationCorrWater_molar * Q.additionalCorr * Q.relaxationCorrMetab
if verbose:
rcorwaterconc = (Q.f_GM*Q.d_GM*Q.R_H2O_GM + Q.f_WM*Q.d_WM*Q.R_H2O_WM + Q.f_CSF*Q.d_CSF*Q.R_H2O_CSF)*H2O_MOLALITY
metabRelaxCorr = 1/Q.R_M
rcorwaterconc = Q.relaxationCorrWater_molar
metabRelaxCorr = Q.relaxationCorrMetab
print(f'Metabolite area = {SMObs:0.2e}')
print(f'Water area = {SH2OObs:0.2e}')
print(f'Relaxation corrected water concentration = {rcorwaterconc:0.2e}')
......@@ -99,35 +102,54 @@ def quantifyWater(mrs,H2OFID,refFID,referenceName,concentrations,metaboliteNames
class QuantificationInfo(object):
""" Class encapsulating the information required to run internal water quantification scaling."""
def __init__(self,TE,T2,tissueFractions,tissueDensity=None):
def __init__(self,TE,T2,tissueFractions=None,tissueDensity=None):
"""Constructor for QuantificationInfo. Requires sequence and tissue information.
Args:
TE (float): Sequence echo time in seconds
T2 (dict:float): Tissue T2 values. Must contain 'H2O_GM', 'H2O_WM','H2O_CSF' and 'METAB' fields. In seconds.
T2 (dict:float): Tissue T2 values. Must contain ('H2O_GM', 'H2O_WM','H2O_CSF') or 'H2O' and 'METAB' fields. In seconds.
tissueFractions (dict:float): Tissue volume fractions, must contain 'WM', 'GM', 'CSF' fields.
tissueDensity (dict:float, optional): Tissue volume fractions, must contain 'WM', 'GM', 'CSF' fields. In units of g/ml.
"""
if tissueDensity is None:
self.d_GM = TISSUE_WATER_DENSITY['GM']
self.d_WM = TISSUE_WATER_DENSITY['WM']
self.d_CSF = TISSUE_WATER_DENSITY['CSF']
if tissueFractions is not None:
if tissueDensity is None:
self.d_GM = TISSUE_WATER_DENSITY['GM']
self.d_WM = TISSUE_WATER_DENSITY['WM']
self.d_CSF = TISSUE_WATER_DENSITY['CSF']
else:
self.d_GM = tissueDensity['GM']
self.d_WM = tissueDensity['WM']
self.d_CSF = tissueDensity['CSF']
# tissue fractions
self.f_GM = tissueFractions['GM']
self.f_WM = tissueFractions['WM']
self.f_CSF = tissueFractions['CSF']
# Calculate T2 relaxation terms
self.R_H2O_GM = np.exp(-TE/T2['H2O_GM'])
self.R_H2O_WM = np.exp(-TE/T2['H2O_WM'])
self.R_H2O_CSF = np.exp(-TE/T2['H2O_CSF'])
self.R_M = np.exp(-TE/T2['METAB'])
# Calculate parts of eqn
self.relaxationCorrWater_molal = (self.f_GM*self.R_H2O_GM + self.f_WM*self.R_H2O_WM + self.f_CSF*self.R_H2O_CSF)*H2O_MOLALITY
self.relaxationCorrWater_molar = (self.f_GM*self.d_GM*self.R_H2O_GM + self.f_WM*self.d_WM*self.R_H2O_WM + self.f_CSF*self.d_CSF*self.R_H2O_CSF)*H2O_MOLALITY
self.relaxationCorrMetab = 1/self.R_M
self.additionalCorr = 1/(1-self.f_CSF)
else:
self.d_GM = tissueDensity['GM']
self.d_WM = tissueDensity['WM']
self.d_CSF = tissueDensity['CSF']
# tissue fractions
self.f_GM = tissueFractions['GM']
self.f_WM = tissueFractions['WM']
self.f_CSF = tissueFractions['CSF']
# Calculate T2 relaxation terms
self.R_H2O_GM = np.exp(-TE/T2['H2O_GM'])
self.R_H2O_WM = np.exp(-TE/T2['H2O_WM'])
self.R_H2O_CSF = np.exp(-TE/T2['H2O_CSF'])
self.R_M = np.exp(-TE/T2['METAB'])
if 'H2O' in T2:
self.R_H2O = np.exp(-TE/T2['H2O'])
else:
self.R_H2O = 1.0
self.R_M = np.exp(-TE/T2['METAB'])
self.relaxationCorrWater_molal = self.R_H2O*H2O_MOLALITY
self.relaxationCorrWater_molar = self.R_H2O*H2O_MOLALITY
self.relaxationCorrMetab = 1/self.R_M
self.additionalCorr = 1.0
def selectT2Values(centralFreq):
""" Select T2 values based on field strength (7T or 3T)"""
......
......@@ -94,9 +94,7 @@ class FitRes(object):
internalRefScaling = quant.quantifyInternal(referenceMetab,self.getConc(),self.metabs)
if (waterRefFID is not None) and\
(tissueFractions is not None) and\
(TE is not None):
if (waterRefFID is not None) and (TE is not None):
refFID = self.predictedFID(mrs,mode=waterReferenceMetab,noBaseline=True)
if T2 == 'Default':
Q = quant.loadDefaultQuantificationInfo(TE,tissueFractions,mrs.centralFrequency/1E6)
......
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