general updates

README.md

 # FSL-MRS
 
-### Installation instructions
 
-```
-	git clone https://git.fmrib.ox.ac.uk/saad/fsl_mrs.git
-	cd fsl_mrs
-	pip install .
-```
+### Description
 
+FSL-MRS is a collection of python modules and wrapper scripts for pre-processing and model fitting of Magnetic Resonance Spectroscopy (MRS) data.
+
+---
+
+### Installation 
+
+    git clone https://git.fmrib.ox.ac.uk/saad/fsl_mrs.git
+    cd fsl_mrs
+    pip install .
+
+---
 
 ### Content
 
+#### Scripts:
+
+- **fsl\_mrs**
+: fit a single spectrum 
+- **fsl\_mrsi**
+: fit a 4D volume of spectra 
+- **fsl\_mrs\_preproc**
+: pre-processing (coil combination, averaging, eddy-current correction) 
+- **fsl\_mrs\_sim**
+: simulate basis
+- **mrs_vis**
+: quick visualisation of the spectrum
+
+---
+
+#### Usage
+
+For each of the wrapper scripts above, simply type `<name_of_script> --help` to get the usage.
+
+
+#### File types
+
+FSL-MRS accepts FID data in NIFTI format. It can also read .RAW format (like LCModel). 
+
+#### Working in python
+
+If you don't want to use the wrapper scripts, you can use the python modules directly in your own python scripts/programs. Here are some examples below:
+
+- Pre-processing
+- Model fitting - single voxel
+- Model fitting - MRSImaging
+
+
 
 
 

fsl_mrs/core.py

@@ -36,20 +36,23 @@ class MRS(object):
         # (now copying the data - looks ugly but better than referencing.
         # now I can run multiple times with different setups)
         self.FID               = FID.copy()
-        self.basis             = basis.copy()
+        self.basis             = basis
+        if basis is not None:
+            self.basis         = basis.copy()
         if H2O is not None:
             self.H2O           = H2O.copy()
         else:
             self.H2O           = None
         self.centralFrequency  = cf     # Hz
         self.bandwidth         = bw     # Hz
-        self.names             = names.copy()
+        self.names             = names
+        if names is not None:
+            self.names         = names.copy()
 
 
         # Set remaining class attributes
         self.Spec              = None        
         self.numPoints         = None
-        self.Spec              = None
         self.numBasis          = None
         
         # Constants 
@@ -169,13 +172,13 @@ class MRS(object):
         return int(first),int(last)
 
 
-    def resample_basis(self):
+    def resample_basis(self,dwelltime):
         """
            Sometimes the basis is simulated using different timings (dwelltime)
            This interpolates the basis to match the FID
         """
         # RESAMPLE BASIS FUNCTION
-        bdt    = self.basis_dwellTime
+        bdt    = dwelltime
         bbw    = 1/bdt
         bn     = self.basis.shape[0]
         
@@ -205,10 +208,11 @@ class MRS(object):
         0 if check successful and -1 if not (also issues warning)
 
         """
-        Spec1 = np.fft.fft(self.FID)
-        Spec2 = np.fft.fft(np.conj(self.FID))
         first,last = self.ppmlim_to_range(ppmlim)
-        if np.linalg.norm(Spec1[first:last]) < np.linalg.norm(Spec2[first:last]):
+        Spec1 = np.real(np.fft.fft(self.FID))[first:last]
+        Spec2 = np.real(np.fft.fft(np.conj(self.FID)))[first:last]
+        
+        if np.linalg.norm(misc.detrend(Spec1,deg=4)) < np.linalg.norm(misc.detrend(Spec2,deg=4)):
             if repare is False:
                 warnings.warn('YOU MAY NEED TO CONJUGATE YOU FID!!!')
                 return -1
@@ -219,6 +223,9 @@ class MRS(object):
         return 0
 
     def conj_FID(self):
+        """
+        Conjugate FID and recalculate spectrum
+        """
         self.FID  = np.conj(self.FID)
         self.Spec = np.fft.fft(self.FID)
 
@@ -254,7 +261,6 @@ class MRS(object):
         """
         if metabs is not None:
             metabs = [m for m in self.names if m not in metabs]
-            #metabs = list(set(self.names)-set(metabs))
             self.ignore(metabs)
 
 
@@ -306,7 +312,7 @@ class MRS(object):
 
         return len(ppmlist)
                 
-    # I/O functions
+    # I/O functions  [NOW OBSOLETE?]
     @staticmethod 
     def read(filename,TYPE='RAW'):
         """

fsl_mrs/scripts/fsl_mrs

@@ -58,7 +58,7 @@ def main():
                    help='combine certain metabolites [repeatable]')
     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',
+    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)')
@@ -68,8 +68,9 @@ def main():
                               help="include default macromolecule peaks")
 
     
-    
     # ADDITONAL OPTIONAL ARGUMENTS
+    optional.add_argument('--t1',type=str,default=None,metavar='IMAGE',
+                          help='structural image (for report)')
     optional.add_argument('--central_frequency',default=None,type=float,
                            help='central frequency in Hz')
     optional.add_argument('--dwell_time',default=None,type=float,
@@ -78,6 +79,8 @@ def main():
                           help='output html report')
     optional.add_argument('--verbose',action="store_true",
                           help='spit out verbose info')
+    optional.add_argument('--phase_correct',action="store_true",
+                          help='do phase correction')
     optional.add_argument('--overwrite',action="store_true",
                           help='overwrite existing output folder')
     optional.add('--config', required=False, is_config_file=True, help='configuration file')
@@ -106,6 +109,8 @@ def main():
     import datetime
     import plotly
     # ######################################################
+    if not args.verbose:
+        warnings.filterwarnings("ignore")
 
     
     # Check if output folder exists
@@ -120,9 +125,9 @@ def main():
             exit()
         else:
             shutil.rmtree(args.output)
-            os.mkdir(args.output)
+            os.makedirs(args.output,exist_ok=True)
     else:
-        os.mkdir(args.output)
+        os.makedirs(args.output,exist_ok=True)
 
 
     # Save chosen arguments
@@ -134,7 +139,13 @@ def main():
         
     #######  Do the work #######
 
-    # Read data/h2o/basis    
+    # Read data/h2o/basis
+    if args.verbose:        
+        print('--->> Read input data and basis\n')
+        print('  {}'.format(args.data))
+        print('  {}\n'.format(args.basis))
+
+
     FID,dataheader            = mrs_io.read_FID(args.data)
     basis, names, basisheader = mrs_io.read_basis(args.basis)
     if args.h2o is not None:
@@ -142,27 +153,37 @@ def main():
     else:
         H2O = None
 
+    # Squeeze the data/h2o (e.g. if from NIFTI single voxel)
+    FID = np.squeeze(FID)
+    if H2O is not None:
+        H2O = np.squeeze(H2O)
+        
     # Collect useful info
     if args.central_frequency is not None:
         cf = args.central_frequency
     elif dataheader['centralFrequency'] is not None:
         cf = dataheader['centralFrequency']
+        if args.verbose:
+            print('      Detected central frequency in header info cf = {} MHz'.format(cf*1E-6))
     else:
-        raise(Exception('Cannot determine central frequency. Please either set it on include it in data header'))
+        raise(Exception('Cannot determine central frequency. Please either set it or 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']
+        if args.verbose:
+            print('      Detected bandwidth in header info bw = {} Hz'.format(bw))
     else:
-        raise(Exception('Cannot determine central frequency. Please either set it on include it in data header'))
+        raise(Exception('Cannot determine bandwidth. Please either set it or 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)
+    if basisheader is not None:
+        if bw != basisheader['bandwidth']:
+            dwell     = 1/basisheader['bandwidth']
+            new_dwell = 1/bw
+            basis     = misc.resample_ts(basis,dwell,new_dwell)
 
 
     # Instantiate MRS object
@@ -170,25 +191,24 @@ def main():
     mrs     = MRS(**MRSargs)
     
     # Check the FID
-    mrs.check_FID(repare=True)
+    conjugated = mrs.check_FID(repare=True)
+    if args.verbose:
+        if conjugated == 1:
+            raise(Warning('Warning :: FID has been checked and conjugated. Please check!'))
 
     # 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))
-
-            
+    if args.phase_correct:
+        if args.verbose:
+            print('--->>  Phase correction\n')
+        mrs.FID  = misc.phase_correct(mrs,mrs.FID)
+        mrs.Spec = np.fft.fft(mrs.FID)
+              
         
     # Keep/Ignore/Combine metabolites
     mrs.keep(args.keep)
     mrs.ignore(args.ignore)
     mrs.combine(args.combine)
 
-    
     # Do the fitting here
     if args.verbose:
         print('--->> Start fitting\n\n')
@@ -197,12 +217,12 @@ def main():
 
         
     # Do the fitting
-    print('     --- Run fitting --- ')
+    if args.verbose:
+        print('     --- Run fitting --- ')
     ppmlim=args.ppmlim
     if ppmlim is not None:
         ppmlim=(ppmlim[0],ppmlim[1])
  
-    #mrs.fit(model = args.model, method=args.algo,ppmlim=args.ppmlim)
 
 
     # Parse metabolite groups
@@ -214,16 +234,23 @@ def main():
 
     # Include Macromolecules? These should have their own metab groups
     if args.add_MM is not None:
+        if not args.verbose:
+            print('Adding macromolecules')
         nMM = mrs.add_MM_peaks()
         G   = [i+max(metab_groups)+1 for i in range(nMM)]
         metab_groups += G
 
 
-    Fitargs = {'mrs':mrs,'ppmlim':ppmlim,
+    Fitargs = {'ppmlim':ppmlim,
                'method':args.algo,'baseline_order':args.baseline_order,
                'metab_groups':metab_groups}
 
-    res = fitting.fit_FSLModel(**Fitargs)
+    if args.verbose:
+        print(mrs)
+        print('Fitting args:')
+        print(Fitargs)
+
+    res = fitting.fit_FSLModel(mrs,**Fitargs)
 
 
     
@@ -239,9 +266,16 @@ def main():
         print('--->> Saving output files to {}\n'.format(args.output))
 
 
-    res.to_file(mrs,filename=os.path.join(args.output,'results_table.csv'))
-    
+    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,'qc.csv'),what='qc')
+    res.to_file(filename=os.path.join(args.output,'all_parameters.csv'),what='parameters')
 
+    # Save image of MRS voxel
+    if args.t1 is not None:
+        datatype = mrs_io.check_datatype(args.data)
+        if datatype == 'NIFTI':            
+            fig = plotting.plot_world_orient(args.t1,args.data)
+            fig.savefig(os.path.join(args.output,'voxel_location.png'))
 
     # Create short HTML report
     #fig = plotting.plotly_fit(mrs,res,ppmlim=ppmlim,proj='abs')
@@ -250,13 +284,14 @@ def main():
     
     
     # Creat HTML report
-    report.create_report(mrs,res,
-                         filename=os.path.join(args.output,'report.html'),
-                         fidfile=args.data,
-                         basisfile=args.basis,
-                         h2ofile=args.h2o,
-                         outdir=args.output,
-                         date=datetime.datetime.now().strftime("%Y-%m-%d %H:%M"))
+    if args.report:
+        report.create_report(mrs,res,
+                             filename=os.path.join(args.output,'report.html'),
+                             fidfile=args.data,
+                             basisfile=args.basis,
+                             h2ofile=args.h2o,
+                             outdir=args.output,
+                             date=datetime.datetime.now().strftime("%Y-%m-%d %H:%M"))
                          
 
     report.fitting_summary_fig(mrs,res,

fsl_mrs/scripts/fsl_mrsi

@@ -60,6 +60,8 @@ fitting_args.add_argument('--add_MM',action="store_true",
 
 
 # ADDITONAL OPTIONAL ARGUMENTS
+optional.add_argument('--conjugate',action="store_true",
+                      help='conjugate the data')
 optional.add_argument('--single_proc',action="store_true",
                       help='do not run in parallel')
 optional.add_argument('--report',action="store_true",
@@ -143,6 +145,8 @@ numBasis                  = basis.shape[1]
 
 # Get array data
 data = np.asanyarray(data_hdr.dataobj)
+if args.conjugate:
+    data = np.conj(data)
 if args.h2o is not None:
     h2o  = np.asanyarray(h2o_hdr.dataobj) 
 mask = np.asanyarray(mask_hdr.dataobj) 
@@ -191,7 +195,7 @@ 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)
+    #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)]
@@ -222,9 +226,9 @@ if args.single_proc:
     results = []
     for idx,FID in enumerate(fid_list):
         if args.verbose:
-            print('{}/{} voxels fitted'.format(idx,len(fid_list)),end='\r')
+            print('{}/{} voxels fitted'.format(idx,len(fid_list))) #,end='\r')
         mrs = MRS(FID=FID,H2O=h2o_list[idx],**MRSargs)
-        mrs.check_FID(repare=True)
+        #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
@@ -306,16 +310,18 @@ if results[0].conc_h2o is not None:
 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')
+if args.conjugate:
+    pred_vol = np.conjugate(pred_vol)
+#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)

fsl_mrs/utils/fitting.py

@@ -14,7 +14,7 @@ from fsl_mrs.utils.constants import *
 from fsl_mrs.utils import mh
 from fsl_mrs.core import MRS
 
-from scipy.optimize import minimize
+from scipy.optimize import minimize,nnls
 
 
 class FitRes(object):
@@ -83,20 +83,44 @@ class FitRes(object):
         for i in range(g):
             self.params_names.extend(["eps_{}".format(i)])
 
+        self.params_names.extend(['Phi0','Phi1'])
+        
+        for i in range(baseline_order+1):
+            self.params_names.extend(["B_real_{}".format(i)])
+
         for i in range(baseline_order+1):
-            self.params_names.extend(["B_{}".format(i)])
+            self.params_names.extend(["B_imag_{}".format(i)])
 
-    def to_file(self,mrs,filename):
+    def to_file(self,filename,mrs=None,what='concentrations'):
+        """
+        Save results to a csv file
+
+        Parameters:
+        -----------
+        filename : str
+        mrs      : MRS obj (only used if what = 'concentrations')
+        what     : one of 'concentrations, 'qc', 'parameters'
+        """
         import pandas as pd
-        df = pd.DataFrame()
-        df['Metab']          = mrs.names
-        df['mMol/kg']        = self.conc_h2o
-        df['%CRLB']          = self.perc_SD[:mrs.numBasis]
-        df['/Cr']            = self.conc_cr
-        df.to_csv(filename)
+        df                   = pd.DataFrame()
 
-        
+        if what is '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':
+            df['Measure'] = ['SNR']
+            df['Value']   = [self.snr]
+        elif what is 'parameters':
+            df['Parameter'] = self.params_names
+            df['Value']     = self.params
 
+        df.to_csv(filename,index=False)
+
+         
+
+        
 def print_params(x,mrs,metab_groups,ref_metab='Cr',scale_factor=1):
     """
        Print parameters 
@@ -152,7 +176,55 @@ def quantify(mrs,concentrations,ref='Cr',to_h2o=False,scale=1):
 
     return res
 
+
+
 def init_gamma_eps(mrs):
+    """
+       Initialise gamma/epsilon parameters
+       
+    """
+
+    ppmlim = (.2,4.2)
+    first,last = mrs.ppmlim_to_range(ppmlim)
+    
+    y = np.fft.fft(mrs.FID)[first:last]
+    #y = np.real(y).flatten() #
+    y = np.concatenate((np.real(y),np.imag(y)),axis=0).flatten()
+
+    def modify_basis(mrs,gamma,eps):
+        bs = mrs.basis * np.exp(-(gamma+1j*eps)*mrs.timeAxis)        
+        bs = np.fft.fft(bs,axis=0)
+        bs = bs[first:last,:]
+        #return np.real(bs) #
+        return np.concatenate((np.real(bs),np.imag(bs)),axis=0)
+            
+    def cf(p):
+        gamma,eps = np.exp(p[0]),p[1]
+        bs        = modify_basis(mrs,gamma,eps)
+        #print(gamma)
+        #print(eps)
+        #beta,_    = nnls(bs,y)
+        beta      = np.linalg.pinv(bs)@y
+        pred      = bs@beta                
+        loss      = np.mean((pred-y)**2)
+        return loss
+
+    x0     = np.array([-10,0])
+    #bounds = [(0,1e3),(-1e3,1e3)]
+    #res    = minimize(cf, x0, method='TNC',bounds=bounds)
+
+    
+    res    = minimize(cf, x0, method='Powell')
+    res    = minimize(cf, x0=res.x, method='Nelder-Mead')
+    res    = minimize(cf, x0=res.x, method='TNC')
+    
+
+    g,e    = np.exp(res.x[0]),res.x[1]
+    
+    #print([g,e])
+    return g,e
+
+def init_gamma_eps_old(mrs):
     """
        Initialise gamma/epsilon parameters
        This is done by summing all the basis FIDs and
@@ -180,6 +252,7 @@ def init_gamma_eps(mrs):
     
     return g,e
 
+
 def init_FSLModel(mrs,metab_groups,baseline_order):
     """
        Initialise params of FSLModel
@@ -189,13 +262,16 @@ def init_FSLModel(mrs,metab_groups,baseline_order):
     # 3. How about phi0 and phi1?
     
     # 1. gamma/eps
+
+
     gamma,eps = init_gamma_eps(mrs)
     new_basis = mrs.basis*np.exp(-(gamma+1j*eps)*mrs.timeAxis)    
     
     data   = np.append(np.real(mrs.FID),np.imag(mrs.FID),axis=0)
     desmat = np.append(np.real(new_basis),np.imag(new_basis),axis=0)
-    con    = np.real(np.linalg.pinv(desmat)@data)   
-                
+    #con    = np.real(np.linalg.pinv(desmat)@data)   
+    con,_  = nnls(desmat,data.flatten())
+    
     # Append 
     x0  = con   # concentrations
         
@@ -233,8 +309,7 @@ def init_gamma_sigma_eps(mrs):
     g   = res.x[0]
     s   = res.x[1]
     e   = res.x[2]
-    
-    
+        
     return g,s,e
 
 def init_FSLModel_Voigt(mrs,metab_groups,baseline_order):
@@ -287,7 +362,9 @@ def prepare_baseline_regressor(mrs,baseline_order,first,last):
     
     return B
 
-def fit_FSLModel(mrs,method='Newton',ppmlim=None,baseline_order=5,metab_groups=None,model='lorentzian'):
+
+
+def fit_FSLModel(mrs,method='Newton',ppmlim=None,baseline_order=5,metab_groups=None,model='lorentzian',x0=None):
     """
         A simplified version of LCModel
     """
@@ -317,18 +394,27 @@ def fit_FSLModel(mrs,method='Newton',ppmlim=None,baseline_order=5,metab_groups=N
     # Results object
     results = FitRes(model)
     results.fill_names(mrs,baseline_order,metab_groups)
-    
-    # Initialise all params
-    x0= init_func(mrs,metab_groups,baseline_order)
 
     # Prepare baseline
     B                 = prepare_baseline_regressor(mrs,baseline_order,first,last)
     results.base_poly = B
-
+    
     # Constants
     g         = results.g
     constants = (freq,time,basis,B,metab_groups,g,data,first,last)    
+
+    if x0 is None:
+        # Initialise all params
+        x0= init_func(mrs,metab_groups,baseline_order)