adding MRSI, complex baseline, and Voigt

fsl_mrs/scripts/fsl_mrs

@@ -14,10 +14,8 @@ import configargparse
 from fsl_mrs import __version__
 from fsl_mrs.utils.splash import splash
 
+# NOTE!!!! THERE ARE MORE IMPORTS IN THE CODE BELOW (AFTER ARGPARSING)
 
-# TODO:
-# Look here: http://blog.vwelch.com/2011/04/combining-configparser-and-argparse.html
-# For a good way to allow config files
 
 
 def main():
@@ -34,9 +32,9 @@ def main():
     #p = argparse.ArgumentParser(description='FSL Magnetic Resonance Spectroscopy Tool')
     p.add_argument('-v','--version', action='version', version=__version__)
     
-    required = p.add_argument_group('required arguments')
+    required     = p.add_argument_group('required arguments')
     fitting_args = p.add_argument_group('fitting options')
-    optional = p.add_argument_group('additional options')
+    optional     = p.add_argument_group('additional options')
 
     # REQUIRED ARGUMENTS
     required.add_argument('--data',
@@ -62,23 +60,26 @@ def main():
                    help='limit the fit to a freq range (default=(.2,4.2))')
     fitting_args.add_argument('--h2o',type=str,metavar='H2O',
                    help='input .H2O file for quantification')
-    fitting_args.add_argument('--T2s',type=float,nargs=4,metavar=('T2metab','T2GM','T2WM','T2CSF'),
-                   help='T2 values (in milliseconds) for metabolites, and water GM, WM, and CSF')
-    fitting_args.add_argument('--volfrac',type=float,nargs=3,metavar=('GM','WM','CSF'),
-                   help='volume fractions of GM, WM, and CSF (in most cases these should add up to one)')
     fitting_args.add_argument('--baseline_order',default=2,type=int,metavar=('ORDER'),
                    help='order of baseline polynomial (default=2)')
     fitting_args.add_argument('--metab_groups',default=0,nargs='+',type=int,
-                              help="metabolite groups.")
+                              help="metabolite groups")
+    fitting_args.add_argument('--add_MM',type=bool,
+                              help="include default macromolecule peaks")
+
     
     
     # ADDITONAL OPTIONAL ARGUMENTS
+    optional.add_argument('--central_frequency',default=None,type=float,
+                           help='central frequency in Hz')
+    optional.add_argument('--dwell_time',default=None,type=float,
+                           help='dwell time in seconds')
     optional.add_argument('--report',action="store_true",
-                   help='output html report')
+                          help='output html report')
     optional.add_argument('--verbose',action="store_true",
-                   help='spit out verbose info')
+                          help='spit out verbose info')
     optional.add_argument('--overwrite',action="store_true",
-                   help='overwrite existing output folder')
+                          help='overwrite existing output folder')
     optional.add('--config', required=False, is_config_file=True, help='configuration file')
 
 
@@ -133,43 +134,60 @@ def main():
         
     #######  Do the work #######
 
+    # Read data/h2o/basis    
+    FID,dataheader            = mrs_io.read_FID(args.data)
+    basis, names, basisheader = mrs_io.read_basis(args.basis)
+    if args.h2o is not None:
+        H2O,_                 = mrs_io.read_FID(args.h2o)
+    else:
+        H2O = None
+
+    # Collect useful info
+    if args.central_frequency is not None:
+        cf = args.central_frequency
+    elif dataheader['centralFrequency'] is not None:
+        cf = dataheader['centralFrequency']
+    else:
+        raise(Exception('Cannot determine central frequency. Please either set it on include it in data header'))
+
+    if args.dwell_time is not None:
+        bw = 1/args.dwell_time
+    elif dataheader['bandwidth'] is not None:
+        bw = dataheader['bandwidth']
+    else:
+        raise(Exception('Cannot determine central frequency. Please either set it on include it in data header'))
+
+        
+    # Resample basis?
+    if bw != basisheader['bandwidth']:
+        dwell     = 1/basisheader['bandwidth']
+        new_dwell = 1/bw
+        basis     = misc.resample_ts(basis,dwell,new_dwell)
+
+
     # Instantiate MRS object
-    mrs = MRS()
+    MRSargs = {'FID':FID,'basis':basis,'names':names,'H2O':H2O,'cf':cf,'bw':bw}
+    mrs     = MRS(**MRSargs)
+    
+    # Check the FID
+    mrs.check_FID(repare=True)
 
+    # Do phase correction
+    mrs.FID  = misc.phase_correct(mrs.FID)
+    mrs.Spec = np.fft.fft(mrs.FID)
+  
     if args.verbose:        
         print('--->> Read input data and basis\n')
         print('  {}\n'.format(args.data))
         print('  {}\n'.format(args.basis))
 
-    # Read data and Basis
-    mrs.read_data(args.data)
-    if os.path.isdir(args.basis):
-        mrs.read_basis_from_folder(args.basis)
-    else:
-        mrs.read_basis_from_file(args.basis)
-
-    # Read H2O file if provided
-    if args.h2o is not None:
-        if args.verbose:
-            print('---> Read h2o file {}\n'.format(args.h2o))
-        mrs.read_h2o(args.h2o)
-
-    # Check that H2O quantif can include T2 factor
-    if mrs.echotime is None:
-        warnings.warn('\n Cannot determine echo time. Will assume it is =30ms!!!!\n')
-        mrs.echotime = 30
-
-    if args.T2s is None:
-        T2s = np.asarray([160,110,80,2000])
-        warnings.warn('\n T2s not provided. Assuming 3T approx values\n')
-    else:
-        T2s = args.T2s
-        
+            
         
     # Keep/Ignore/Combine metabolites
     mrs.keep(args.keep)
     mrs.ignore(args.ignore)
     mrs.combine(args.combine)
+
     
     # Do the fitting here
     if args.verbose:
@@ -186,9 +204,6 @@ def main():
  
     #mrs.fit(model = args.model, method=args.algo,ppmlim=args.ppmlim)
 
-    # Do phase correction
-    mrs.FID  = misc.phase_correct(mrs.FID)
-    mrs.Spec = np.fft.fft(mrs.FID)
 
     # Parse metabolite groups
     metab_groups = args.metab_groups
@@ -196,20 +211,24 @@ def main():
         metab_groups = [0]*mrs.numBasis
     if len(metab_groups) != mrs.numBasis:
         raise(Exception('Found {} metab_groups but there are {} basis functions'.format(len(metab_groups),mrs.numBasis)))
-    
-    res = fitting.fit_FSLModel(mrs,method=args.algo,
-                               ppmlim=ppmlim,
-                               baseline_order=args.baseline_order,
-                               metab_groups=metab_groups)
+
+    # Include Macromolecules? These should have their own metab groups
+    if args.add_MM is not None:
+        nMM = mrs.add_MM_peaks()
+        G   = [i+max(metab_groups)+1 for i in range(nMM)]
+        metab_groups += G
+
+
+    Fitargs = {'mrs':mrs,'ppmlim':ppmlim,
+               'method':args.algo,'baseline_order':args.baseline_order,
+               'metab_groups':metab_groups}
+
+    res = fitting.fit_FSLModel(**Fitargs)
 
 
     
     stop = time.time()
 
-    # Do some extra bits for quantification
-    #mrs.post_process(metab=['Cr','PCr'],scale=8.0,T2s=T2s,volfrac=args.volfrac)
-    #mrs.reference_matebolite = 'Cr+PCr'
-    #print(mrs.all_con_names_h2o)
     
     # Report on the fitting
     if args.verbose:
@@ -239,7 +258,9 @@ def main():
                          outdir=args.output,
                          date=datetime.datetime.now().strftime("%Y-%m-%d %H:%M"))
                          
-    
+
+    report.fitting_summary_fig(mrs,res,
+                               filename=os.path.join(args.output,'fit_summary.png'))
      
     if args.verbose:
         print('\n\n\nDone.')

fsl_mrs/scripts/fsl_mrsi

+#!/usr/bin/env python
+
+# fsl_mrs - wrapper script for MRSI fitting
+#
+# Author: Saad Jbabdi <saad@fmrib.ox.ac.uk>
+#
+# Copyright (C) 2019 University of Oxford 
+
+import configargparse
+
+from fsl_mrs import __version__
+from fsl_mrs.utils.splash import splash
+
+
+
+
+# Parse command-line arguments
+p = configargparse.ArgParser(add_config_file_help=False,
+                             description="FSL Magnetic Resonance Spectroscopy Imaging Wrapper Script") 
+
+p.add_argument('-v','--version', action='version', version=__version__)
+    
+required     = p.add_argument_group('required arguments')
+fitting_args = p.add_argument_group('fitting options')
+optional     = p.add_argument_group('additional options')
+
+
+# REQUIRED ARGUMENTS
+required.add_argument('--data',
+                      required=True,type=str,metavar='<str>.NII',
+                      help='input NIFTI file')
+required.add_argument('--basis',
+                      required=True,type=str,metavar='<str>',
+                      help='.BASIS file or folder containing basis spectra (will read all .RAW files within)')
+required.add_argument('--mask',
+                      required=True,type=str,metavar='<str>',
+                      help='mask volume')			 
+required.add_argument('--output',
+                      required=True,type=str,metavar='<str>',
+                      help='output folder')
+required.add_argument('--central_frequency',required=True,type=float,
+                      help='central frequency in Hz')
+required.add_argument('--dwell_time',required=True,type=float,
+                      help='dwell time in seconds')
+
+# FITTING ARGUMENTS
+fitting_args.add_argument('--algo',default='Newton',type=str,
+                          help='algorithm [Newton (fast) or MH (slow)]')
+fitting_args.add_argument('--ppmlim',default=(.2,4.2),type=float,nargs=2,metavar=('LOW','HIGH'),
+                          help='limit the fit to a freq range (default=(.2,4.2))')
+fitting_args.add_argument('--h2o',type=str,metavar='H2O',
+                          help='input H2O reference NIFTI for quantification')
+fitting_args.add_argument('--baseline_order',default=2,type=int,metavar=('ORDER'),
+                          help='order of baseline polynomial (default=2)')
+fitting_args.add_argument('--metab_groups',default=0,nargs='+',type=int,
+                          help="metabolite groups")
+fitting_args.add_argument('--add_MM',action="store_true",
+                          help="include default macromolecule peaks")
+
+
+
+# ADDITONAL OPTIONAL ARGUMENTS
+optional.add_argument('--single_proc',action="store_true",
+                      help='do not run in parallel')
+optional.add_argument('--report',action="store_true",
+                      help='output html report')
+optional.add_argument('--verbose',action="store_true",
+                      help='spit out verbose info')
+optional.add_argument('--overwrite',action="store_true",
+                      help='overwrite existing output folder')
+optional.add('--config', required=False, is_config_file=True, help='configuration file')
+
+    
+# Parse command-line arguments
+args = p.parse_args()
+
+# Output kickass splash screen
+if args.verbose:
+    splash(logo='mrsi')
+
+
+# ######################################################
+# DO THE IMPORTS AFTER PARSING TO SPEED UP HELP DISPLAY
+import time, os, sys, shutil, warnings
+import numpy as np
+from fsl_mrs.core import MRS
+from fsl_mrs.utils import mrs_io
+from fsl_mrs.utils import report
+from fsl_mrs.utils import fitting
+from fsl_mrs.utils import plotting
+from fsl_mrs.utils import misc
+import datetime
+#import plotly
+import nibabel as nib
+from functools import partial
+# ######################################################
+
+    
+# Check if output folder exists
+overwrite = args.overwrite
+if os.path.exists(args.output):
+    if not overwrite:
+        print("Folder '{}' exists. Are you sure you want to delete it? [Y,N]".format(args.output))
+        response = input()
+        overwrite = response.upper() == "Y"
+    if not overwrite:
+        print('Early stopping...')
+        exit()
+    else:
+        shutil.rmtree(args.output)
+        os.mkdir(args.output)
+else:
+    os.mkdir(args.output)
+
+
+# Save chosen arguments
+with open(os.path.join(args.output,"options.txt"),"w") as f:
+    f.write(str(args))
+    f.write("\n--------\n")
+    f.write(p.format_values())
+        
+        
+#######  Do the work #######
+
+# Read data/h2o/basis
+if args.verbose:
+    print('--->> Load data\n')
+    print('      FID   : {}'.format(args.data))
+data_hdr = nib.load(args.data)
+if args.h2o is not None:
+    if args.verbose:
+        print('      H2O   : {}'.format(args.h2o))
+    h2o_hdr  = nib.load(args.h2o)
+if args.verbose:
+    print('      mask  : {}'.format(args.mask))
+mask_hdr = nib.load(args.mask)
+
+numPoints                 = data_hdr.shape[-1]
+if args.verbose:
+    print('      basis : {}'.format(args.basis))
+basis, names, basisheader = mrs_io.readLCModelBasis(args.basis,numPoints)
+numBasis                  = basis.shape[1]
+
+# Get array data
+data = np.asanyarray(data_hdr.dataobj)
+if args.h2o is not None:
+    h2o  = np.asanyarray(h2o_hdr.dataobj) 
+mask = np.asanyarray(mask_hdr.dataobj) 
+if data.shape[:3] != mask.shape:
+    mask = mask[:,:,None]
+
+# data as list (for parallel proc)
+fid_list = misc.volume_to_list(data,mask)
+if args.h2o is not None:
+    h2o_list = misc.volume_to_list(h2o,mask)
+
+
+
+# # Resample basis?
+cf = args.central_frequency
+bw = 1/args.dwell_time
+# if bw != basisheader['bandwidth']:
+#     dwell     = 1/basisheader['bandwidth']
+#     new_dwell = 1/bw
+#     basis     = misc.resample_ts(basis,dwell,new_dwell)
+
+# ppmlim for fitting
+ppmlim=args.ppmlim
+if ppmlim is not None:
+    ppmlim=(ppmlim[0],ppmlim[1])
+
+# Parse metabolite groups
+metab_groups = args.metab_groups
+if metab_groups == 0:
+    metab_groups = [0]*numBasis
+if len(metab_groups) != numBasis:
+    raise(Exception('Found {} metab_groups but there are {} basis functions'.format(len(metab_groups),numBasis)))
+
+
+# Store info in disctionaries to be passed to MRS and fitting
+MRSargs = {'basis':basis,'names':names,'cf':cf,'bw':bw}
+Fitargs = {'ppmlim':ppmlim,
+           'method':args.algo,'baseline_order':args.baseline_order,
+           'model':'lorentzian'}
+
+import multiprocessing as mp
+from functools import partial
+global_counter = mp.Value('L')
+
+
+# Define some ugly local functions for parallel processing
+def runvoxel(FIDH2O,MRSargs,Fitargs):
+    mrs = MRS(FID=FIDH2O[0],H2O=FIDH2O[1],**MRSargs)        
+    mrs.check_FID(repare=True)
+    if args.add_MM:
+        n = mrs.add_MM_peaks()
+        new_metab_groups = [i+max(metab_groups)+1 for i in range(n)]
+        new_metab_groups = metab_groups + new_metab_groups
+    else:
+        new_metab_groups = metab_groups.copy()
+    res = fitting.fit_FSLModel(mrs, **Fitargs,metab_groups=new_metab_groups)   
+    with global_counter.get_lock():
+        global_counter.value += 1
+    return res
+
+def parallel_runs(data_list):
+    pool    = mp.Pool(processes=mp.cpu_count())
+    func    = partial(runvoxel, MRSargs=MRSargs, Fitargs=Fitargs) 
+    results = pool.map_async(func, data_list)
+    return results
+
+    
+# Fitting
+if args.verbose:
+    print('\n--->> Start fitting\n\n')
+    print('    Algorithm = [{}]\n'.format(args.algo))
+
+warnings.filterwarnings("ignore")
+if args.single_proc:
+    if args.verbose:
+        print('    Running sequentially (are you sure about that?) ')
+    results = []
+    for idx,FID in enumerate(fid_list):
+        if args.verbose:
+            print('{}/{} voxels fitted'.format(idx,len(fid_list)),end='\r')
+        mrs = MRS(FID=FID,H2O=h2o_list[idx],**MRSargs)
+        mrs.check_FID(repare=True)
+        n = mrs.add_MM_peaks()
+        new_metab_groups = [i+max(metab_groups)+1 for i in range(n)]
+        new_metab_groups = metab_groups + new_metab_groups
+        res = fitting.fit_FSLModel(mrs, **Fitargs,metab_groups=new_metab_groups)   
+        results.append(res)
+else:
+    if args.verbose:
+        print('    Parallelising over {} workers '.format(mp.cpu_count()))
+    t0  = time.time()
+    results = parallel_runs(zip(fid_list,h2o_list))
+
+    while not results.ready():
+        if args.verbose:
+            print('{}/{} voxels completed'.format(global_counter.value,len(fid_list)), end='\r')
+        time.sleep(1)
+    if args.verbose:
+        print('{}/{} voxels completed'.format(global_counter.value,len(fid_list)), end='\r')
+        print('\n\nFitting done in {:0f} secs.'.format(time.time()-t0))
+
+
+    if not results.successful():
+        raise(Exception("Fitting unsuccessful :-(((((("))
+    else:
+        results = results.get()
+            
+# Save output files
+if args.verbose:
+    print('--->> Saving output files to {}\n'.format(args.output))
+    
+# ResFit --> Images
+# Store concentrations, uncertainties, residuals, predictions
+
+# All params
+folder = os.path.join(args.output,'all_params')
+os.mkdir(folder)
+params = np.asarray([r.params for r in results])
+names  = results[0].params_names
+report.save_params(params,names,data_hdr,mask,folder)
+
+# All params
+folder = os.path.join(args.output,'all_perc_crlb')
+os.mkdir(folder)
+params = np.asarray([r.perc_SD for r in results])
+names  = results[0].params_names
+report.save_params(params,names,data_hdr,mask,folder)
+
+# Concentrations
+folder = os.path.join(args.output,'conc')
+os.mkdir(folder)
+params = np.asarray([r.conc for r in results])
+names  = results[0].names
+report.save_params(params,names,data_hdr,mask,folder)
+
+# Concentrations/Cr
+if results[0].conc_cr is not None:
+    folder = os.path.join(args.output,'conc_cr')
+    os.mkdir(folder)
+    params = np.asarray([r.conc_cr for r in results])
+    names  = results[0].names
+    report.save_params(params,names,data_hdr,mask,folder)
+
+# Concentrations/Cr+Pcr
+if results[0].conc_cr_pcr is not None:
+    folder = os.path.join(args.output,'conc_cr_pcr')
+    os.mkdir(folder)
+    params = np.asarray([r.conc_cr_pcr for r in results])
+    names  = results[0].names
+    report.save_params(params,names,data_hdr,mask,folder)
+
+# Concentrations/H2O
+if results[0].conc_h2o is not None:
+    folder = os.path.join(args.output,'conc_h2o')
+    os.mkdir(folder)
+    params = np.asarray([r.conc_h2o for r in results])
+    names  = results[0].names
+    report.save_params(params,names,data_hdr,mask,folder)
+
+# Store predictions
+pred_list = [r.pred for r in results]
+pred_vol  = misc.list_to_volume(pred_list,mask,dtype=np.complex)
+# check if data has been conjugated and if so conjugate the predictions
+FID   = fid_list[0] 
+Spec1 = np.fft.fft(FID)
+Spec2 = np.fft.fft(np.conj(FID))
+mrs = MRS(FID=fid_list[0],**MRSargs)
+first,last = mrs.ppmlim_to_range(ppmlim)
+if np.linalg.norm(Spec1[first:last]) < np.linalg.norm(Spec2[first:last]):
+    print('Data has been conjugated')
+    pred_vol = np.conj(pred_vol)
+else:
+    print('Data has NOT been conjugated')
+img       = nib.Nifti1Image(pred_vol,affine=data_hdr.affine,header=data_hdr.header)
+filename  = os.path.join(os.path.join(args.output,'pred.nii.gz'))
+nib.save(img, filename)
+
+# MSE
+mse_list = [r.mse for r in results]
+mse_vol  = misc.list_to_volume(mse_list,mask)
+img       = nib.Nifti1Image(mse_vol,affine=data_hdr.affine)
+filename  = os.path.join(os.path.join(args.output,'mse.nii.gz'))
+nib.save(img, filename)
+
+# SNR
+snr_list = [r.snr for r in results]
+snr_vol  = misc.list_to_volume(snr_list,mask)
+img       = nib.Nifti1Image(snr_vol,affine=data_hdr.affine)
+filename  = os.path.join(os.path.join(args.output,'snr.nii.gz'))
+nib.save(img, filename)
+
+
+if args.verbose:
+    print('\n\n\nDone.')

fsl_mrs/utils/fitting.py

@@ -12,6 +12,7 @@ from fsl_mrs.utils.models import *
 from fsl_mrs.utils.misc import *
 from fsl_mrs.utils.constants import *
 from fsl_mrs.utils import mh
+from fsl_mrs.core import MRS
 
 from scipy.optimize import minimize
 
@@ -20,7 +21,8 @@ class FitRes(object):
     """
        Collects fitting results
     """
-    def __init__(self):
+    def __init__(self,model):
+        self.model        = model
         self.params       = None
         self.crlb         = None
         self.mcmc_samples = None
@@ -30,6 +32,7 @@ class FitRes(object):
         self.perc_SD      = None
         self.conc_h2o     = None
         self.conc_cr      = None
+        self.conc_cr_pcr  = None
         self.cov          = None
         self.params_names = None
         self.residuals    = None
@@ -49,7 +52,7 @@ class FitRes(object):
         #out += " cov            = {}\n".format(self.cov)
         out += " phi0 (deg)     = {}\n".format(self.phi0_deg)
         out += " phi1 (deg/ppm) = {}\n".format(self.phi1_deg_per_ppm)
-        out += "inv_gamma_ms    = {}\n".format(self.inv_gamma_sec)
+        out += "inv_gamma_sec   = {}\n".format(self.inv_gamma_sec)
         out += "eps_ppm         = {}\n".format(self.eps_ppm)
         out += "b_norm          = {}\n".format(self.b_norm)
 
@@ -58,7 +61,7 @@ class FitRes(object):
         
     def fill_names(self,mrs,baseline_order=0,metab_groups=None):
         """
-        mrs : MRS Object
+        mrs            : MRS Object
         baseline_order : int
         metab_groups   : list (by default assumes single metab group)
         """
@@ -127,15 +130,20 @@ def quantify(mrs,concentrations,ref='Cr',to_h2o=False,scale=1):
     """
         Quantification
     """
-    ref_con = concentrations[mrs.names.index(ref)]
+    if isinstance(ref,list):
+        ref_con = 0
+        for met in ref:
+            ref_con += concentrations[mrs.names.index(met)]
+    else:
+        ref_con = concentrations[mrs.names.index(ref)]
 
-    if to_h2o is not True:
+    if to_h2o is not True or mrs.H2O is None:
         QuantifFactor = scale/ref_con
     else:
         ref_fid       = mrs.basis[:,mrs.names.index(ref)]
         ref_area      = calculate_area(mrs,ref_con*ref_fid)
         H2O_area      = calculate_area(mrs,mrs.H2O) 
-        refH2O_ratio   = ref_area/H2O_area
+        refH2O_ratio  = ref_area/H2O_area
         H2O_protons   = 2
         ref_protons   = num_protons[ref]