first commit

fsl_mrs/scripts/fsl_mrs_preproc

+#!/usr/bin/env python
+
+# fsl_mrs_preproc - wrapper script for MRS preprocessing
+#
+# Author: Saad Jbabdi <saad@fmrib.ox.ac.uk>
+#
+# Copyright (C) 2019 University of Oxford 
+# SHBASECOPYRIGHT
+
+# Quick imports
+#import argparse
+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)
+
+
+def main():
+    # Parse command-line arguments
+    p = configargparse.ArgParser(add_config_file_help=False,description="FSL Magnetic Resonance Spectroscopy - Preprocessing")
+
+    p.add_argument('-v','--version', action='version', version=__version__)
+
+    required     = p.add_argument_group('required arguments')
+    optional     = p.add_argument_group('additional options')
+
+
+    # REQUIRED ARGUMENTS
+    required.add_argument('--data',
+                          required=True,type=str,nargs='+',metavar='<str>',
+                          help='list of input files')    
+    required.add_argument('--output',
+                          required=True,type=str,metavar='<str>',
+                          help='output folder')
+
+    # ADDITONAL OPTIONAL ARGUMENTS
+
+    optional.add_argument('--t1',type=str,default=None,metavar='IMAGE',
+                          help='structural image (for report)')
+    optional.add_argument('--verbose',action="store_true",
+                          help='spit out verbose info')
+    optional.add_argument('--conjugate',action="store_true",
+                          help='apply conjugate to FID')
+    optional.add_argument('--overwrite',action="store_true",
+                          help='overwrite existing output folder')
+    optional.add('--config', required=False, is_config_file=True, help='configuration file')
+
+
+    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')
+
+    # Parse command-line arguments
+    args = p.parse_args()
+    
+    # Output kickass splash screen
+    if args.verbose:
+        splash(logo='mrs')
+    
+
+    # ######################################################
+    # 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 preproc
+    from fsl_mrs.utils import plotting
+    from fsl_mrs.utils import misc
+    import datetime
+    # ######################################################
+
+    
+    
+    # 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 #######
+    if args.verbose:
+        print('Load the data....')
+    # Read all data
+    FIDlist= []
+    shape    = None
+    datatype = None
+    for filename in args.data:
+        dtype      = mrs_io.check_datatype(filename)
+        fid,header = mrs_io.read_FID(filename)
+        shape_pre  = fid.shape
+        fid        = np.squeeze(fid)
+        if args.conjugate:
+            fid = np.conj(fid)
+        if shape is not None:
+            if fid.shape != shape:
+                raise(Exception('FIDs are incompatible in shapes'))
+            if dtype != datatype:
+                raise(Exception('FID files must all be of the same data type (e.g. RAW or NIFTI)'))
+        else:
+            shape    = fid.shape
+            datatype = dtype
+        
+        FIDlist.append(fid)
+
+    if args.verbose:
+        print('.... Found {} repeats of the FIDs with shape {} each.\n\n'.format(len(FIDlist),shape_pre))
+        
+    # Determine if coils have been combined already
+    if args.verbose:
+        print('.... Determine if coil combination is needed')
+    if len(shape) == 1:
+        do_coil_combine = False
+        if args.verbose:
+            print('   ----> NO.\n')
+    else:
+        do_coil_combine = True
+        if args.verbose:
+            print('  ----> YES.\n')
+
+    # Do preproc
+    if args.verbose:
+            print('Begin proprocessing.... ')
+            
+    if do_coil_combine:
+        FIDcomblist = []
+        if args.verbose:
+            print('... Coil Combination ...')
+        for idx,fid in enumerate(FIDlist):
+            if args.verbose:
+                print('{}/{}'.format(idx,len(FIDlist)),end='\r')
+            FIDcomb = preproc.combine_FIDs(list(fid.T),'svd',do_dephase=True,do_prewhiten=True)
+            FIDcomblist.append(FIDcomb)
+        FIDcomblist = preproc.dephase(FIDcomblist)
+    else:
+        FIDcomblist = FIDlist
+        
+    # Average the data (if asked to do so - this might not be the case for FMRS)
+    if args.verbose:
+        print('... Average FIDs ...')
+        FID_final = preproc.combine_FIDs(FIDcomblist,'mean',do_dephase=True,do_prewhiten=False)
+
+    
+    # Save the data
+    if datatype == 'RAW':
+        filename = os.path.join(args.output,'metab.RAW')
+        mrs_io.saveRAW(filename,FID_final,header)
+    elif datatype == 'NIFTI':
+        filename = os.path.join(args.output,'metab.nii')
+        data = np.reshape(FID_final,shape_pre[:4])
+        mrs_io.saveNIFTI(filename,data,header.affine)
+
+
+    # Produce html report
+    if args.verbose:
+        print('Create report')
+    if args.t1 is not None:
+        if datatype == 'NIFTI':
+            fig = plotting.plot_world_orient(args.t1,args.data[0])
+            fig.savefig(os.path.join(args.output,'voxel_location.png'))
+
+    centralFrequency = 123E6
+    if args.central_frequency is not None:
+        centralFrequency = args.central_frequency
+
+
+    if args.dwell_time is not None:
+        dwellTime = args.dwell_time
+    else:
+        import nibabel as nib
+        img = nib.load(args.data[0])
+        dwellTime = img.header.get_zooms()[3]
+    bandwidth = 1/dwellTime
+
+    
+    fig = plotting.plot_spectra(FIDlist=FIDcomblist,
+                                bandwidth=bandwidth,
+                                centralFrequency=centralFrequency,
+                                ppmlim=(.2,4.2),
+                                single_FID=FID_final,plot_avg=False)
+    fig.savefig(os.path.join(args.output,'coil_combined.png'))
+
+    # Produce QC measures?
+
+
+
+    
+
+        
+if __name__ == '__main__':
+    main()

fsl_mrs/scripts/mrs_vis

+#!/usr/bin/env python
+
+# mrs_vis - quick MRS visualisation
+#
+# Author: Saad Jbabdi <saad@fmrib.ox.ac.uk>
+#
+# Copyright (C) 2019 University of Oxford 
+# SHBASECOPYRIGHT
+
+import sys
+
+def main():
+    nargs = len(sys.argv) - 1
+    if nargs < 3:
+        print('\nUsage : mrs_vis <FIDfile> <centralFrequency (Hz)> <dwellTime (sec)> [True/False]\n')
+        exit()
+
+    FIDfile = sys.argv[1]
+    cf      = float(sys.argv[2])
+    dt      = float(sys.argv[3])
+    conj    = False
+    if nargs==4:
+        if sys.argv[4] == 'True':
+            conj = True
+
+    from fsl_mrs.utils.plotting import plot_spectrum
+    from fsl_mrs.utils.mrs_io import read_FID
+    import matplotlib.pyplot as plt
+
+    FID,_     = read_FID(FIDfile)
+    if len(FID.shape) > 1:
+        FID = FID.flatten()
+    if conj:
+        from numpy import conjugate as conj
+        FID = conj(FID)
+
+    plt.figure(figsize=(8,8))
+    fig     = plot_spectrum(FID,1/dt,cf,ppmlim=(.2,4.2))
+    plt.show()
+    
+
+if __name__ == '__main__':
+    main()

fsl_mrs/utils/preproc.py

+#!/usr/bin/env python
+
+# preproc.py - Preprocessing utilities
+#
+# Author: Saad Jbabdi <saad@fmrib.ox.ac.uk>
+#
+# Copyright (C) 2019 University of Oxford 
+# SHBASECOPYRIGHT
+
+import numpy as np
+from scipy.signal import tukey
+
+# tmp
+from scipy.sparse.linalg import svds
+
+from fsl_mrs.utils import misc
+
+# Functions
+#    apodize --> check!
+#    phase_freq_align --> check!
+#    channel_combine --> check!
+#    eddy_correct : subtract water FID phase from metab FID
+
+# TODO
+# residual water removal
+# outlier removal (of individual repeats)
+
+# QC measures?
+#    Peak signal / noise variance (between coils / channels)
+# 
+
+
+def truncate(FID,k,first_or_last='last'):
+    """
+    Truncate parts of a FID
+    
+    Parameters:
+    -----------
+    FID           : array-like
+    k             : int (number of timepoints to remove)
+    first_or_last : either 'first' or 'last' (which bit to truncate)
+
+    Returns:
+    --------
+    array-like
+    """
+    FID_trunc = FID.copy()
+    
+    if first_or_last == 'first':
+        return FID_trunc[k:]
+    elif first_or_last == 'last':
+        return FID_trunc[:-k]
+    else:
+        raise(Exception("Last parameter must either be 'first' or 'last'"))
+
+
+def apodize(FID):
+    """
+       FID Apodization
+
+    TODO: optional choice of filter
+
+    Parameters:
+    -----------
+     FID : array-like
+    Returns:
+    --------
+    Array-like
+    """
+    window = tukey(FID.size * 2)[FID.size:]
+    return window*FID
+
+
+# Phase-Freq alignment functions
+def align_FID(mrs,src_FID,tgt_FID,ppmlim=None):
+    """
+       Phase and frequency alignment
+
+    Parameters
+    ----------
+    mrs : MRS Object
+    src_FID : array-like
+    tgt_FID : array-like
+    ppmlim : tuple
+
+    Returns
+    -------
+    array-like
+    """
+
+    # Internal functions so they can see globals
+    def shift_phase_freq(FID,phi,eps,extract=True):        
+        sFID = np.exp(-1j*phi)*misc.shift_FID(mrs,FID,eps)
+        if extract:
+            sFID = misc.extract_spectrum(mrs,sFID,ppmlim=ppmlim)
+        return sFID
+    def cf(p):
+        phi    = p[0] #phase shift
+        eps    = p[1] #freq shift
+        FID    = shift_phase_freq(src_FID,phi,eps)
+        target = misc.extract_spectrum(mrs,tgt_FID,ppmlim=ppmlim)
+        xx     = np.linalg.norm(FID-target)
+        return xx
+    x0  = np.array([0,0])
+    res = minimize(cf, x0, method='Powell')
+    phi = res.x[0]
+    eps = res.x[1]
+    
+    return shift_phase_freq(src_FID,phi,eps,extract=False)
+
+def get_target_FID(FIDlist,target='mean'):
+    """
+    target can be 'mean' or 'first' or 'nearest_to_mean'
+    """
+    if target == 'mean':
+        return sum(FIDlist) / len(FIDlist)
+    elif target == 'first':
+        return FIDlist[0].copy()
+    elif target == 'nearest_to_mean':
+        avg = sum(FIDlist) / len(FIDlist)
+        d   = [np.linalg.norm(fid-avg) for fid in FIDlist]
+        return FIDlist[np.argmin(d)].copy()
+    else:
+        raise(Exception('Unknown target type {}'.format(target)))
+
+
+def phase_freq_align(FIDlist,bandwidth,centralFrequency,ppmlim=None,niter=2,verbose=False):
+    """
+    Algorithm:
+       Average spectra
+       Loop over all spectra and find best phase/frequency shifts
+       Iterate
+
+    TODO: 
+       test the code
+       add dedrifting?
+
+    Parameters:
+    -----------
+    FIDlist          : list
+    bandwidth        : float (unit=Hz)
+    centralFrequency : float (unit=Hz)
+    ppmlim           : tuple
+    niter            : int
+    verbose          : bool
+    
+    Returns:
+    --------
+    list of FID aligned to each other
+    """
+    all_FIDs = FIDlist.copy()
+    for iter in range(niter):
+        if verbose:
+            print(' ---- iteration {} ----\n'.format(iter))
+        target = get_target_FID(FIDlist,target='nearest_to_mean')
+        MRSargrgs = {'FID':target,'bw':bandwidth,'cf':centralFrequency}
+        mrs = MRS(**MRSargs)
+        for idx,FID in enumerate(all_FIDs):
+            if verbose:
+                print('... aligning FID number {}'.format(idx),end='\r')
+            all_FIDs[idx] = align_FID(mrs,FID,target,ppmlim=ppmlim)
+        if verbose:
+            print('\n')
+    return all_FIDs
+
+
+# Methods for FID combination
+
+def dephase(FIDlist):
+    """
+      Uses first data point of each FID to dephase each FID in list
+      Returns a list of dephased FIDs
+    """
+    return [fid*np.exp(-1j*np.angle(fid[0])) for fid in FIDlist]
+
+def prewhiten(FIDlist,prop=.1):
+    """
+       Uses noise covariance to prewhiten data
+
+    Parameters:
+    -----------
+    FIDlist : list of FIDs
+    prop    : proportion of data used to estimate noise covariance
+
+    Returns : 
+    list of FIDs
+    pre-whitening matrix
+    """
+    FIDs = np.asarray(FIDlist,dtype=np.complex)
+    # Estimate noise covariance
+    start = int((1-prop)*FIDs.shape[1])
+    # Raise warning if not enough samples
+    if (FIDs.shape[1]-start)<1.5*FIDs.shape[0]:
+        raise(Warning('You may not have enough samples to robustly estimate the noise covariance'))
+    C     = np.cov(FIDs[:,start:])
+    D,V   = np.linalg.eigh(C)
+    # Pre-whitening matrix
+    W     = np.diag(1/np.sqrt(D))@np.conj(V.T)
+    # Pre-whitened data
+    FIDs = W@FIDs
+    return list(FIDs),W
+
+
+def svd_reduce(FIDlist,W=None,return_alpha=False):
+    """
+    Combine different channels by SVD method
+    Based on C.T. Rodgers and M.D. Robson, Magn Reson Med 63:881–891, 2010
+
+    Parameters:
+    -----------
+    FIDlist      : list of FIDs
+    W            : pre-whitening matrix (only used to calculate sensitivities)
+    return_alpha : return sensitivities?
+
+    Returns:
+    --------
+    array-like (FID)
+    array-like (sensitivities) - optional
+    """
+    FIDs  = np.asarray(FIDlist)
+    #U,S,V = np.linalg.svd(FIDs)
+    U,S,V = svds(FIDs,k=1)  #this is much faster
+
+    # get arbitrary amplitude
+    iW = np.eye(FIDs.shape[0])
+    if W is not None:
+        iW = np.linalg.inv(W)
+    amp  = np.linalg.norm(np.mean(iW@FIDs,axis=0))
+    
+    # combined channels
+    FID = amp*V[0,:].T
+
+    if return_alpha:
+        # sensitivities per channel
+        if W is not None:
+            alpha = iW@U[:,0]
+        else:
+            alpha = U[:,0]
+
+        return FID,alpha
+    else:
+        return FID
+
+def combine_FIDs(FIDlist,method,do_prewhiten=False,do_dephase=False,do_phase_correct=False):
+    """
+       Combine FIDs (either from multiple coils or multiple averages
+    
+    Parameters:
+    -----------
+    FIDlist   : list of FIDs
+    method    : one of 'mean', 'svd'
+    prewhiten : bool
+    dephase   : bool
+
+    Returns:
+    --------
+    array-like
+
+    """
+
+    # Pre-whitening
+    W = None
+    if do_prewhiten:
+        FIDlist,W = prewhiten(FIDlist)
+
+    # Dephasing
+    if do_dephase:
+        FIDlist   = dephase(FIDlist)
+
+    # Combining channels
+    if method == 'mean':
+        return sum(FIDlist) / len(FIDlist)
+    elif method == 'svd':
+        return svd_reduce(FIDlist,W)
+    else:
+        raise(Exception("Unknown method '{}'. Should be either 'mean' or 'svd'".format(method)))
+
+
+def eddy_correct(FIDmet,FIDw):
+    """
+    Subtract water phase from metab phase
+    Typically do this after coil combination
+
+    Parameters:
+    -----------
+    FIDmet : array-like (FID for the metabolites)
+    FIDw   : array-like (water FID)
+
+    Returns:
+    --------
+    array-like (corrected metabolite FID)
+    """
+
+    return FIDmet / np.exp(1j*np.angle(FIDw))