Merge branch '1.0.4' into 'master'

V1.0.4

See merge request fsl/fsl_mrs!10

.gitignore

@@ -18,3 +18,4 @@ example_usage/fsl_mrs_preproc/
 example_usage/fsl_mrs_preproc_simple/
 example_usage/fsl_mrs_proc/
 example_usage/report.html
+*.egg-info

CHANGELOG.rst

 This document contains the FSL-MRS release history in reverse chronological order.
 
+1.0.4 (Friday 14th August 2020)
+-------------------------------
+- Fixed bug in automatic conjugation facility of fsl_mrs_preproc
+- jmrui text file reader now handles files with both FID and spectra
+
 1.0.3 (Friday 10th July 2020)
 -----------------------------
 - Changed to pure python version of HLSVDPRO (hlsvdpropy). Slight speed penalty

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.3'
+version = '1.0.4'
 release = version
 
 # From PM's fsleyes doc

fsl_mrs/core/MRS.py

@@ -463,6 +463,13 @@ class MRS(object):
         self.FID         = FID.copy()
         self.numPoints   = self.FID.size
         self.Spec        = misc.FIDToSpec(self.FID)
+
+
+    def plot(self,ppmlim=(0.2,4.2)):
+        from fsl_mrs.utils.plotting import plot_spectrum
+
+        plot_spectrum(self,ppmlim=ppmlim)
+
         
 
         

fsl_mrs/scripts/fsl_mrs_preproc

@@ -65,6 +65,8 @@ def main():
                           help='spit out verbose info')
     optional.add_argument('--conjugate', action="store_true",
                           help='apply conjugate to FID')
+    optional.add_argument('--no_conjugate', action="store_true",
+                          help='Forbid automatic conjugation')                      
     optional.add_argument('--overwrite', action="store_true",
                           help='overwrite existing output folder')
     optional.add_argument('--report', action="store_true",
@@ -116,31 +118,53 @@ def main():
 
     # Read all data
     # Suppressed data
-    FIDlist, hdr, shape, datatype = readData(args.data, args.conjugate)
+    FIDlist, hdr, shape, datatype = readData(args.data)
     if args.verbose:
         print(f'.... Found {len(FIDlist)} repeats of the data'
               f' FIDs with shape {shape} each.\n\n')
     # Reference data
-    Reflist, refhdr, refshape, _ = readData(args.reference, args.conjugate)
+    Reflist, refhdr, refshape, _ = readData(args.reference)
     if args.verbose:
         print(f'.... Found {len(Reflist)} repeats of the reference'
               f' FIDs with shape {refshape} each.\n\n')
 
     if args.quant is not None:
         # Separate quant data
-        quantlist, quanthdr, quantshape, _ = readData(args.quant,
-                                                      args.conjugate)
+        quantlist, quanthdr, quantshape, _ = readData(args.quant)
         if args.verbose:
             print(f'.... Found {len(quantlist)} repeats of the quant'
                   f' FIDs with shape {quantshape} each.\n\n')
 
     if args.ecc is not None:
         # Separate ecc data
-        ecclist, ecchdr, eccshape, _ = readData(args.ecc, args.conjugate)
+        ecclist, ecchdr, eccshape, _ = readData(args.ecc)
         if args.verbose:
             print(f'.... Found {len(ecclist)} repeats of the ecc FIDs'
                   f' with shape {eccshape} each.\n\n')
 
+    # Data conjugation
+    if args.no_conjugate:
+        if args.verbose:
+            print('No conjugation explicitly set,'
+                  'skipping automatic determination.')
+        pass
+    elif args.conjugate or auto_conj(FIDlist,hdr):
+        if args.verbose and args.conjugate:
+            print('Conjugation explicitly set,'
+                  'applying conjugation.')
+        elif args.verbose:
+            print('Automatic conjugation applied')
+
+        def iterFunction(flist):
+            return [np.conj(fid) for fid in flist]
+
+        FIDlist = iterFunction(FIDlist)
+        Reflist = iterFunction(Reflist)
+        if args.quant is not None:
+            quantlist = iterFunction(quantlist)
+        if args.ecc is not None:
+            ecclist = iterFunction(ecclist)
+
     # Determine if coils have been combined already
     if args.verbose:
         print('.... Determine if coil combination is needed')
@@ -401,24 +425,15 @@ def main():
                 fig.savefig(op.join(args.output, 'voxel_location.png'))
 
 
-def readData(files, conjugate):
-    from fsl_mrs.core import MRS
+def readData(files):
     FIDlist = []
     shape = None
     datatype = None
-    auto_conj = 0
     for filename in files:
         dtype = mrs_io.check_datatype(filename)
         fid, header = mrs_io.read_FID(filename)
         fid = np.squeeze(fid)
 
-        if fid.ndim > 1:
-            mrs = MRS(FID=np.mean(fid, axis=1), header=header)
-            auto_conj += mrs.check_FID(repair=True)
-        else:
-            mrs = MRS(FID=fid, header=header)
-            auto_conj += mrs.check_FID(repair=True)            
-
         if shape is not None:
             if fid.shape != shape:
                 raise(Exception('FIDs are incompatible in shapes'))
@@ -431,10 +446,6 @@ def readData(files, conjugate):
 
         FIDlist.append(fid)
 
-    auto_conj /= len(FIDlist)
-    if conjugate or auto_conj > 0.5:
-        FIDlist = [np.conj(f) for f in FIDlist]
-
     return FIDlist, header, shape, datatype
 
 
@@ -454,5 +465,23 @@ def saveData(FIDlist, hdr, outputDir, basename, datatype='NIFTI'):
             mrs_io.fsl_io.saveNIFTI(filename, f, hdr)
 
 
+def auto_conj(FIDlist, hdr):
+    from fsl_mrs.core import MRS
+    auto_conj = 0
+    for fid in FIDlist:
+        if fid.ndim > 1:
+            mrs = MRS(FID=np.mean(fid, axis=1), header=hdr)
+            auto_conj += mrs.check_FID(repair=True)
+        else:
+            mrs = MRS(FID=fid, header=hdr)
+            auto_conj += mrs.check_FID(repair=True)
+
+    auto_conj /= len(FIDlist)
+    if auto_conj > 0.5:
+        return True
+    else:
+        return False
+
+
 if __name__ == '__main__':
     main()

fsl_mrs/utils/mrs_io/jmrui_io.py

@@ -44,7 +44,10 @@ def readjMRUItxt(filename,unpack_header=True):
                 continue
                 
             if recordData:
-                data.append(list(map(float,line.split())))
+                curr_data = line.split()
+                if len(curr_data) > 2:
+                    curr_data = curr_data[:2]
+                data.append(list(map(float, curr_data)))
 
     # Reshape data
     data = np.concatenate([np.array(i) for i in data])
@@ -53,7 +56,7 @@ def readjMRUItxt(filename,unpack_header=True):
     # Clean up header
     header = translateHeader(header)
 
-    return data,header
+    return data, header
 
 # Translate jMRUI header to mandatory fields
 def translateHeader(header):

fsl_mrs/utils/plotting.py

@@ -1156,4 +1156,80 @@ def plot_world_orient(t1file,voxfile):
 
     return plt.gcf()
 
+def plot_world_orient_ax(t1file,voxfile,ax,orientation='axial'):
+    """
+    Plot sagittal/coronal/axial T1 centered on voxel in axes passed
+    Overlay voxel in red
+    Plots in world coordinates with the 'MNI' convention
+    """
+    t1      = nib.load(t1file)
+    vox     = nib.load(voxfile)
+    t1_data = t1.get_fdata()
+
+    # transform t1 data into world coordinate system,
+    # resampling to 1mm^3.
+    #
+    # This is a touch fiddly because the affine_transform
+    # function uses the destination coordinate system as
+    # data indices. So we need to remove any translation
+    # in the original affine, to ensure that the origin
+    # of the destination coordinate system is forced to
+    # be (0, 0, 0).
+    #
+    # We figure all of this out by transforming the image
+    # bounding box coordinates into world coordinates,
+    # and then figuring out a suitable offset and resampling
+    # data shape from them.
+    extents       = zip((0, 0, 0), t1_data.shape)
+    extents       = np.asarray(list(it.product(*extents)))
+    extents       = ijk2xyz(extents, t1.affine)
+    offset        = extents.min(axis=0)
+    offaff        = np.eye(4)
+    offaff[:3, 3] = -offset
+    shape         = (extents.max(axis=0) - offset).astype(np.int)[:3]
+
+    t1_data = ndimage.affine_transform(t1_data,
+                                       np.dot(offaff, t1.affine),
+                                       output_shape=shape,
+                                       order=3,
+                                       mode='constant',
+                                       cval=0)
+
+    # centre of MRS voxel in (transformed) T1 voxel space
+    ijk = xyz2ijk(ijk2xyz([0,0,0],vox.affine),np.linalg.inv(offaff))
+    i,j,k = ijk
+
+    
+    si,sj,sk = np.array(vox.header.get_zooms())[:3]/2
+    # Do the plotting
+    plt.sca(ax)
+    if orientation == 'sagital':
+        slice = np.squeeze(t1_data[int(i),:,:]).T
+        rect = np.asarray([[j-sj,k-sk],
+                        [j+sj,k-sk],
+                        [j+sj,k+sk],
+                        [j-sj,k+sk],
+                        [j-sj,k-sk]])
+        do_plot_slice(slice,rect)
+    elif orientation == 'coronal':
+    
+        slice = np.squeeze(t1_data[:,int(j),:]).T
+        rect = np.asarray([[i-si,k-sk],
+                        [i+si,k-sk],
+                        [i+si,k+sk],
+                        [i-si,k+sk],
+                        [i-si,k-sk]])
+        do_plot_slice(slice,rect)
+        ax.invert_xaxis()
+    elif orientation == 'axial':
+        slice = np.squeeze(t1_data[:,:,int(k)]).T
+        rect = np.asarray([[i-si,j-sj],
+                        [i+si,j-sj],
+                        [i+si,j+sj],
+                        [i-si,j+sj],
+                        [i-si,j-sj]])
+        do_plot_slice(slice,rect)
+        ax.invert_xaxis()
+
+    return plt.gca()
 

fsl_mrs/utils/synthetic/synthetic_from_basis.py

@@ -7,7 +7,7 @@
 # SHBASECOPYRIGHT
 
 import numpy as np
-from fsl_mrs.utils.misc import ts_to_ts
+from fsl_mrs.utils.misc import ts_to_ts,FIDToSpec,SpecToFID,rescale_FID
 from fsl_mrs.utils import mrs_io
 def standardConcentrations(basisNames):
     """Return standard concentrations for 1H MRS brain metabolites for those which match basis set names."""
@@ -161,6 +161,7 @@ def syntheticFromBasis(basis,
     dwelltime = 1/bandwidth
     basis_dwelltime= 1/basis_bandwidth
     basis = ts_to_ts(basis,basis_dwelltime,dwelltime,points)
+    # basis = rescale_FID(basis,scale=100)
 
     # Create the spectrum
     baseFID    = np.zeros((points),np.complex128)
@@ -169,7 +170,10 @@ def syntheticFromBasis(basis,
                               dwellTime*points,
                               points)
     for b,c,e,g,s in zip(basis.T,concentrations,eps,gammas,sigmas):
-        baseFID += c*b*np.exp(-(1j*e+g+timeAxis*s**2)*timeAxis)
+        tmp   = b*np.exp(-(1j*e+g+timeAxis*s**2)*timeAxis)
+        M     = FIDToSpec(tmp)
+        baseFID += SpecToFID(M*c)
+        
 
     # Add baseline
     # TO DO