Merge branch 'v1.0.1' into 'master'

V1.0.1

See merge request !5

CHANGELOG.rst

 This document contains the FSL-MRS release history in reverse chronological order.
 
+1.0.1 (Friday 19th June 2020)
+--------------------------------
+- Output folder in fsl_mrs_proc will now be created if it does not exist.
+- fsl_mrs_proc now handles data with a singleton coil dimension correctly.
+- --ind_scale and --disable_MH_priors options added to fsl_mrs and fsl_mrsi.
 
 1.0.0 (Wednesday 17th June 2020)
 --------------------------------

docs/user_docs/conf.py

@@ -25,7 +25,7 @@ copyright = f'{date.year}, Will Clarke & Saad Jbabdi, University of Oxford, Oxfo
 author = 'William Clarke'
 
 # The full version, including alpha/beta/rc tags
-version = '1.0.0'
+version = '1.0.1'
 release = version
 
 # From PM's fsleyes doc

docs/user_docs/fitting.rst

@@ -158,8 +158,10 @@ Below are detailed explanations of some of the optional arguments in the wrapper
     Add macromolecule peaks at the following frequencies: 0.9, 1.2, 1.4, 1.7 ppm and a doublet at 2.08 & 3.0 ppm
 :code:`--lorentzian`        
     By default the lineshape is a Voigt (lorentizian+gaussian). Use this flag to set to Lorentzian.
-
-
+:code:`--ind_scale`        
+    Allow independent scaling of specified basis spectra before fitting. For example this can be used to independently scale empirically measured macromolecules combined with simulated metabolite spectra.
+:code:`--disable_MH_priors`        
+    Disable the priors on the MH fitting. The priors are tuned for *in vivo* human brain spectroscopy. Use this option if your spectra has significantly different line widths, phases or large shifts. E.g. in liquid phase phantom or (potentially) pre-clinical systems. 
 
 The wrapper scripts can also take a configuration file as an input. For example, say we have a text file called :code:`config.txt` which contains the below:
 

docs/user_docs/troubleshooting.rst

@@ -9,4 +9,10 @@ MRS specific questions may be better answered on the `MRSHub forums <https://for
 
 1. Unable to find example data
     If you installed FSL-MRS through conda the example data can be downloaded directly from the GitLab repository `folder <https://git.fmrib.ox.ac.uk/fsl/fsl_mrs/-/tree/master/example_usage>`_.
- 
\ No newline at end of file
+ 
+2. Poor fits
+    Two problems are commonley diagnosed when poor fits are seen:
+
+    1)  Basis spectra are inconsistently scaled. For example empirically derived macromolecular basis spectra can be orders of magnitude larger than the other basis spectra. Before fitting, fsl_mrs(i) scales the magnitude of the data and basis spectra to a known range. Relative scales are preserved within the basis spectra. To permit fsl_mrs(i) to apply different scales to individual basis spectra use the :code:`--ind_scale` option with a list of basis names.
+
+    2)  The data might have parameters unlike a 7T or 3T human *in vivo* brain spectrum. I.e. the spectrum originates from a pre-clinical system or from phantom. In this case the MCMC priors which are suitable for *in vivo* human case can be disabled using the :code:`--disable_MH_priors` option.
\ No newline at end of file

fsl_mrs/core/MRS.py

@@ -252,13 +252,13 @@ class MRS(object):
 
         
     # Helper functions
-    def processForFitting(self,ppmlim=(.2,4.2),):
+    def processForFitting(self,ppmlim=(.2,4.2),ind_scaling=None):
         """ Apply rescaling and run the conjugation checks"""
         self.check_FID(ppmlim=ppmlim,repair=True)
         self.check_Basis(ppmlim=ppmlim,repair=True)
-        self.rescaleForFitting()
+        self.rescaleForFitting(ind_scaling=ind_scaling)
 
-    def rescaleForFitting(self,scale=100):
+    def rescaleForFitting(self,scale=100,ind_scaling=None):
         """ Apply rescaling across data, basis and H20"""
         
         scaledFID,scaling = misc.rescale_FID(self.FID,scale=scale)
@@ -267,7 +267,17 @@ class MRS(object):
             self.H2O *= scaling
 
         if self.basis is not None:
-            self.basis,scaling_basis = misc.rescale_FID(self.basis,scale=scale)
+            if ind_scaling is None:
+                self.basis,scaling_basis = misc.rescale_FID(self.basis,scale=scale)
+            else:
+                index = [self.names.index(n) for n in ind_scaling]
+                mask = np.zeros_like(self.names,dtype=bool)
+                mask[index] = True
+                self.basis[:,~mask],scaling_basis = misc.rescale_FID(self.basis[:,~mask],scale=scale)
+                scaling_basis = [scaling_basis]
+                for idx in index:
+                    self.basis[:,idx],tmp = misc.rescale_FID(self.basis[:,idx],scale=scale)
+                    scaling_basis.append(tmp)
         else:
             scaling_basis = None
 

fsl_mrs/core/MRSI.py

@@ -10,10 +10,9 @@
 
 import numpy as np
 from fsl_mrs.core import MRS
-from fsl_mrs.utils import mrs_io,plotting,fitting,misc
+from fsl_mrs.utils import mrs_io, misc
 import matplotlib.pyplot as plt
 import nibabel as nib
-import os.path as op
 from fsl_mrs.utils.mrs_io.fsl_io import saveNIFTI
 
 class MRSI(object):
@@ -67,6 +66,7 @@ class MRSI(object):
         self.rescale        = False
         self.keep           = None
         self.ignore         = None
+        self.ind_scaling    = None
 
         self._store_scalings = None
         
@@ -165,7 +165,7 @@ class MRSI(object):
             mrs.check_FID(repair=True)
 
         if self.rescale:
-            mrs.rescaleForFitting()      
+            mrs.rescaleForFitting(ind_scaling=self.ind_scaling)      
                 
     def plot(self,mask=True,ppmlim=(0.2,4.2)):
         if mask:

fsl_mrs/scripts/fsl_mrs

@@ -77,6 +77,12 @@ def main():
     fitting_args.add_argument('--lorentzian', action="store_true",
                               help='Enable purely lorentzian broadening'
                                    ' (default is Voigt)')
+    fitting_args.add_argument('--ind_scale', default=None, type=str,
+                              nargs='+',
+                              help='List of basis spectra to scale'
+                                   ' independently of other basis spectra.')
+    fitting_args.add_argument('--disable_MH_priors', action="store_true",
+                              help="Disable MH priors.")
 
     # ADDITIONAL OPTIONAL ARGUMENTS
     optional.add_argument('--t1', type=str, default=None, metavar='IMAGE',
@@ -244,7 +250,7 @@ def main():
 
     # Rescale FID, H2O and basis to have nice range
     if not args.no_rescale:
-        mrs.rescaleForFitting()
+        mrs.rescaleForFitting(ind_scaling=args.ind_scale)
 
     # Do phase correction
     if args.phase_correct:
@@ -278,18 +284,21 @@ def main():
         G = [i+max(metab_groups)+1 for i in range(nMM)]
         metab_groups += G
 
+    # Choose fitting lineshape model.
     if args.lorentzian:
         Fitargs = {'ppmlim': ppmlim,
                    'method': args.algo,
                    'baseline_order': args.baseline_order,
                    'metab_groups': metab_groups,
-                   'model': 'lorentzian'}
+                   'model': 'lorentzian',
+                   'disable_mh_priors': args.disable_MH_priors}
     else:
         Fitargs = {'ppmlim': ppmlim,
                    'method': args.algo,
                    'baseline_order': args.baseline_order,
                    'metab_groups': metab_groups,
-                   'model': 'voigt'}
+                   'model': 'voigt',
+                   'disable_mh_priors': args.disable_MH_priors}
 
     if args.verbose:
         print(mrs)

fsl_mrs/scripts/fsl_mrsi

@@ -73,6 +73,12 @@ def main():
     fitting_args.add_argument('--lorentzian', action="store_true",
                               help="Enable purely lorentzian broadening"
                                    " (default is Voigt)")
+    fitting_args.add_argument('--ind_scale', default=None, type=str,
+                              nargs='+',
+                              help='List of basis spectra to scale'
+                                   ' independently of other basis spectra.')
+    fitting_args.add_argument('--disable_MH_priors', action="store_true",
+                              help="Disable MH priors.")
 
     # ADDITONAL OPTIONAL ARGUMENTS
     optional.add_argument('--TE', type=float, default=None, metavar='TE',
@@ -195,6 +201,9 @@ def main():
     else:
         Fitargs['model'] = 'voigt'
 
+    if args.disable_MH_priors:
+        Fitargs['disable_mh_priors'] = True
+
     # Echo time
     if args.TE is not None:
         echotime = args.TE*1E-3

fsl_mrs/utils/preproc/align.py

@@ -286,7 +286,7 @@ def phase_freq_align_report(inFIDs,outFIDs,hdr,phi,eps,ppmlim=None,shift=True,ht
         elif op.isdir(op.dirname(html)) and op.splitext(html)[1]=='.html':
             htmlfile = html
         else:
-            raise ValueError('html must be file ')
+            raise ValueError('Report html path must be file or directory. ')
         
         opName = 'Align'
         timestr = datetime.now().strftime("%H:%M:%S")
@@ -412,7 +412,7 @@ def phase_freq_align_diff_report(inFIDs0,inFIDs1,outFIDs0,outFIDs1,hdr,eps,phi,p
         elif op.isdir(op.dirname(html)) and op.splitext(html)[1]=='.html':
             htmlfile = html
         else:
-            raise ValueError('html must be file ')
+            raise ValueError('Report html path must be file or directory. ')
         
         opName = 'AlignDiff'
         timestr = datetime.now().strftime("%H:%M:%S")
@@ -440,103 +440,3 @@ def phase_freq_align_diff_report(inFIDs0,inFIDs1,outFIDs0,outFIDs1,hdr,eps,phi,p
         return fig,fig2,fig3
     else:
         return fig,fig2,fig3
-
-# Matplotlib
-# def phase_freq_align_report(inFIDs,outFIDs,hdr,phi,eps,ppmlim=None):
-#     from matplotlib import pyplot as plt
-#     from fsl_mrs.core import MRS
-#     from fsl_mrs.utils.preproc.combine import combine_FIDs
-#     from fsl_mrs.utils.plotting import styleSpectrumAxes
-#     fig, axs = plt.subplots(2, 2,figsize=(10,10))
-#     axs[0,0].plot(np.array(phi)*(180.0/np.pi))
-#     axs[0,0].set_ylabel(r'$\phi$ (degrees)')
-#     axs[0,1].plot(eps)
-#     axs[0,1].set_ylabel('Shift (Hz)')
-
-#     meanIn = combine_FIDs(inFIDs,'mean')
-#     meanOut = combine_FIDs(outFIDs,'mean')
-#     toMRSobj = lambda fid : MRS(FID=fid,header=hdr)
-#     meanIn = toMRSobj(meanIn)
-#     meanOut = toMRSobj(meanOut)
-
-#     toPlotIn,toPlotOut = [],[]
-#     for fid in inFIDs:
-#         toPlotIn.append(toMRSobj(fid))
-#     for fid in outFIDs:
-#         toPlotOut.append(toMRSobj(fid))
-#     for fid in toPlotIn:
-#             axs[1,0].plot(fid.getAxes(ppmlim=ppmlim),np.real(fid.getSpectrum(ppmlim=ppmlim)))
-#     axs[1,0].plot(meanIn.getAxes(ppmlim=ppmlim),np.real(meanIn.getSpectrum(ppmlim=ppmlim)),'k')        
-#     styleSpectrumAxes(ax=axs[1,0])
-#     for fid in toPlotOut:
-#             axs[1,1].plot(fid.getAxes(ppmlim=ppmlim),np.real(fid.getSpectrum(ppmlim=ppmlim)))
-#     axs[1,1].plot(meanOut.getAxes(ppmlim=ppmlim),np.real(meanOut.getSpectrum(ppmlim=ppmlim)),'k')        
-#     styleSpectrumAxes(ax=axs[1,1])
-#     plt.tight_layout()
-#     plt.show()
-
-#     fig = plt.figure(figsize=(10,10))
-#     plt.plot(meanIn.getAxes(ppmlim=ppmlim),np.real(meanIn.getSpectrum(ppmlim=ppmlim)),'k',label='Unaligned', linewidth=2)
-#     plt.plot(meanOut.getAxes(ppmlim=ppmlim),np.real(meanOut.getSpectrum(ppmlim=ppmlim)),'r--',label='Aligned', linewidth=2)
-#     styleSpectrumAxes(ax=plt.gca())
-#     plt.legend()
-#     plt.rcParams.update({'font.size': 12})
-#     plt.show()
-
-# def phase_freq_align_diff_report(inFIDs0,inFIDs1,outFIDs0,outFIDs1,hdr,eps,phi,ppmlim=None,diffType='add',shift=True):
-#     from matplotlib import pyplot as plt
-#     from fsl_mrs.core import MRS
-#     from fsl_mrs.utils.preproc.combine import combine_FIDs
-#     from fsl_mrs.utils.plotting import styleSpectrumAxes
-#     fig, axs = plt.subplots(2, 2,figsize=(10,10))
-#     axs[0,0].plot(np.array(phi)*(180.0/np.pi))
-#     axs[0,0].set_ylabel(r'$\phi$ (degrees)')
-#     axs[0,1].plot(eps)
-#     axs[0,1].set_ylabel('Shift (Hz)')
-
-#     diffFIDListIn = []
-#     diffFIDListOut = []
-#     for fid0i,fid1i,fid0o,fid1o in zip(inFIDs0,inFIDs1,outFIDs0,outFIDs1):
-#         if diffType.lower() == 'add':
-#             diffFIDListIn.append(add(fid1i,fid0i))
-#             diffFIDListOut.append(add(fid1o,fid0o))
-#         elif diffType.lower() == 'sub':
-#             diffFIDListIn.append(subtract(fid1i,fid0i))
-#             diffFIDListOut.append(subtract(fid1o,fid0o))
-#         else:
-#             raise ValueError('diffType must be add or sub.') 
-
-#     meanIn = combine_FIDs(diffFIDListIn,'mean')
-#     meanOut = combine_FIDs(diffFIDListOut,'mean')
-#     toMRSobj = lambda fid : MRS(FID=fid,header=hdr)
-#     meanIn = toMRSobj(meanIn)
-#     meanOut = toMRSobj(meanOut)
-
-#     if shift:
-#         axis = 'ppmshift'
-#     else:
-#         axis = 'ppm'
-
-#     toPlotIn,toPlotOut = [],[]
-#     for fid in diffFIDListIn:
-#         toPlotIn.append(toMRSobj(fid))
-#     for fid in diffFIDListOut:
-#         toPlotOut.append(toMRSobj(fid))
-#     for fid in toPlotIn:
-#             axs[1,0].plot(fid.getAxes(ppmlim=ppmlim, axis=axis),np.real(fid.getSpectrum(ppmlim=ppmlim, shift=shift)))
-#     axs[1,0].plot(meanIn.getAxes(ppmlim=ppmlim, axis=axis),np.real(meanIn.getSpectrum(ppmlim=ppmlim, shift=shift)),'k')        
-#     styleSpectrumAxes(ax=axs[1,0])
-#     for fid in toPlotOut:
-#             axs[1,1].plot(fid.getAxes(ppmlim=ppmlim, axis=axis),np.real(fid.getSpectrum(ppmlim=ppmlim, shift=shift)))
-#     axs[1,1].plot(meanOut.getAxes(ppmlim=ppmlim, axis=axis),np.real(meanOut.getSpectrum(ppmlim=ppmlim, shift=shift)),'k')        
-#     styleSpectrumAxes(ax=axs[1,1])
-#     plt.tight_layout()
-#     plt.show()
-
-#     fig = plt.figure(figsize=(10,10))
-#     plt.plot(meanIn.getAxes(ppmlim=ppmlim, axis=axis),np.real(meanIn.getSpectrum(ppmlim=ppmlim, shift=shift)),'k',label='Unaligned', linewidth=2)
-#     plt.plot(meanOut.getAxes(ppmlim=ppmlim, axis=axis),np.real(meanOut.getSpectrum(ppmlim=ppmlim, shift=shift)),'r--',label='Aligned', linewidth=2)
-#     styleSpectrumAxes(ax=plt.gca())
-#     plt.legend()
-#     plt.rcParams.update({'font.size': 12})
-#     plt.show()

fsl_mrs/utils/preproc/combine.py

@@ -237,7 +237,7 @@ def combine_FIDs_report(inFIDs,outFID,hdr,ncha=2,ppmlim = (0.0,6.0),method='not
         elif op.isdir(op.dirname(html)) and op.splitext(html)[1]=='.html':
             htmlfile = html
         else:
-            raise ValueError('html must be file ')
+            raise ValueError('Report html path must be file or directory. ')
         
         opName = 'Combination'
         timestr = datetime.now().strftime("%H:%M:%S")

fsl_mrs/utils/preproc/eddycorrect.py

@@ -93,7 +93,7 @@ def eddy_correct_report(inFID,outFID,phsRef,hdr,ppmlim = (0.2,4.2),html=None):
         elif op.isdir(op.dirname(html)) and op.splitext(html)[1]=='.html':
             htmlfile = html
         else:
-            raise ValueError('html must be file ')
+            raise ValueError('Report html path must be file or directory. ')
         
         opName = 'ECC'
         timestr = datetime.now().strftime("%H:%M:%S")
@@ -116,21 +116,3 @@ def eddy_correct_report(inFID,outFID,phsRef,hdr,ppmlim = (0.2,4.2),html=None):
         return fig,fig2
     else:
         return fig,fig2
-
-# Matplotlib
-# def eddy_correct_report(inFID,outFID,hdr,ppmlim = (0.2,4.2)):
-    # from matplotlib import pyplot as plt
-    # from fsl_mrs.core import MRS
-    # from fsl_mrs.utils.plotting import styleSpectrumAxes
-
-    # toMRSobj = lambda fid : MRS(FID=fid,header=hdr)
-    # plotIn = toMRSobj(inFID)
-    # plotOut = toMRSobj(outFID)    
-    
-    # fig = plt.figure(figsize=(10,10))
-    # plt.plot(plotIn.getAxes(ppmlim=ppmlim),np.real(plotIn.getSpectrum(ppmlim=ppmlim)),'k',label='Uncorrected', linewidth=2)
-    # plt.plot(plotOut.getAxes(ppmlim=ppmlim),np.real(plotOut.getSpectrum(ppmlim=ppmlim)),'r--',label='Corrected', linewidth=2)
-    # styleSpectrumAxes(ax=plt.gca())
-    # plt.legend()
-    # plt.rcParams.update({'font.size': 12})
-    # plt.show()

fsl_mrs/utils/preproc/filtering.py

@@ -81,7 +81,7 @@ def apodize_report(inFID,outFID,hdr,plotlim = (0.2,6),html=None):
         elif op.isdir(op.dirname(html)) and op.splitext(html)[1]=='.html':
             htmlfile = html
         else:
-            raise ValueError('html must be file ')
+            raise ValueError('Report html path must be file or directory. ')
         
         opName = 'Apodization'
         timestr = datetime.now().strftime("%H:%M:%S")

fsl_mrs/utils/preproc/general.py

@@ -97,7 +97,7 @@ def add_subtract_report(inFID,inFID2,outFID,hdr,ppmlim=(0.2,4.2),function='Not s
         elif op.isdir(op.dirname(html)) and op.splitext(html)[1]=='.html':
             htmlfile = html
         else:
-            raise ValueError('html must be file ')
+            raise ValueError('Report html path must be file or directory. ')
         
         opName = function
         timestr = datetime.now().strftime("%H:%M:%S")

fsl_mrs/utils/preproc/phasing.py

@@ -116,7 +116,7 @@ def phaseCorrect_report(inFID,outFID,hdr,position,ppmlim=(2.8,3.2),html=None):
         elif op.isdir(op.dirname(html)) and op.splitext(html)[1]=='.html':
             htmlfile = html
         else:
-            raise ValueError('html must be file ')
+            raise ValueError('Report html path must be file or directory. ')
         
         opName = 'Phase correction'
         timestr = datetime.now().strftime("%H:%M:%S")

fsl_mrs/utils/preproc/remove.py

@@ -117,7 +117,7 @@ def hlsvd_report(inFID,outFID,hdr,limits,limitUnits = 'ppm',plotlim = (0.2,6),ht
         elif op.isdir(op.dirname(html)) and op.splitext(html)[1]=='.html':
             htmlfile = html
         else:
-            raise ValueError('html must be file ')
+            raise ValueError('Report html path must be file or directory. ')
         
         opName = 'HLSVD'
         timestr = datetime.now().strftime("%H:%M:%S")

fsl_mrs/utils/preproc/shifting.py

@@ -181,7 +181,7 @@ def shift_report(inFID,outFID,inHdr,outHdr,ppmlim = (0.2,4.2),html=None,function
         elif op.isdir(op.dirname(html)) and op.splitext(html)[1]=='.html':
             htmlfile = html
         else:
-            raise ValueError('html must be file ')
+            raise ValueError('Report html path must be file or directory. ')
         
         operation,function,description = reportStrings(function)
 

fsl_mrs/utils/preproc/unlike.py

@@ -139,7 +139,7 @@ def identifyUnlikeFIDs_report(goodFIDs,badFIDs,hdr,keepIndicies,rmIndicies,metri
         elif op.isdir(op.dirname(html)) and op.splitext(html)[1]=='.html':
             htmlfile = html
         else:
-            raise ValueError('html must be file ')
+            raise ValueError('Report html path must be file or directory. ')
         
         opName = 'BadAverageRemoval'
         timestr = datetime.now().strftime("%H:%M:%S")