Added baseline to synthetic from basis.

fsl_mrs/tests/test_utils_baseline.py

+from fsl_mrs.utils import baseline
+from fsl_mrs.core import MRS
+import numpy as np
+
+
+def test_regressor_creation():
+
+    # Create dummy mrs
+    mrs = MRS(FID=np.zeros((1000,)), cf=100, bw=2000, nucleus='1H')
+    baseline_order = 3
+    ppmlim = (-1, 1)
+    B = baseline.prepare_baseline_regressor(mrs, baseline_order, ppmlim)
+
+    # Need some better tests...
+    assert B.shape == (1000, 8)
+    assert np.sum(B[:, 0]) == 100
+    assert np.sum(B[:, 1]) == 100j

fsl_mrs/tests/test_utils_synthetic.py

 from fsl_mrs.utils import synthetic as syn
 from fsl_mrs.utils.misc import FIDToSpec
+from fsl_mrs.core import MRS
 import numpy as np
+from pathlib import Path
+
+testsPath = Path(__file__).parent
+basis_path = testsPath / 'testdata/fsl_mrs/steam_basis'
+
 
 def test_noisecov():
     # Create positive semi-definite noise covariance
@@ -36,4 +42,30 @@ def test_syntheticFID():
     spec /= np.max(np.abs(spec))
     testSpec /= np.max(np.abs(testSpec))
 
-    assert np.isclose(spec,FIDToSpec(testFID[0]),atol = 1E-2,rtol = 1E0).all()
\ No newline at end of file
+    assert np.isclose(spec,FIDToSpec(testFID[0]),atol = 1E-2,rtol = 1E0).all()
+
+
+def test_syntheticFromBasis():
+    # TO DO
+    pass
+
+
+def test_syntheticFromBasis_baseline():
+
+    fid, header, _ = syn.syntheticFromBasisFile(str(basis_path),
+                                                baseline=((0.0, 0.0),),
+                                                concentrations={'Mac': 2.0},
+                                                noisecovariance=[[0.0]])
+
+    mrs = MRS(FID=fid, header=header)
+    mrs.conj_FID()
+
+    fid, header, _ = syn.syntheticFromBasisFile(str(basis_path),
+                                                baseline=((1.0, 1.0),),
+                                                concentrations={'Mac': 2.0},
+                                                noisecovariance=[[0.0]])
+
+    mrs2 = MRS(FID=fid, header=header)
+    mrs2.conj_FID()
+
+    assert np.allclose(mrs2.get_spec(), mrs.get_spec() + np.complex(1.0, -1.0))

fsl_mrs/utils/baseline.py

+#!/usr/bin/env python
+
+# baseline.py - Functions associated with the baseline description
+#
+# Author: Saad Jbabdi <saad@fmrib.ox.ac.uk>
+#         William Clarke <william.clarke@ndcn.ox.ac.uk>
+#
+# Copyright (C) 2019 University of Oxford
+# SHBASECOPYRIGHT
+
+import numpy as np
+from fsl_mrs.utils.misc import regress_out
+
+
+def prepare_baseline_regressor(mrs, baseline_order, ppmlim):
+    """
+       Complex baseline is polynomial
+
+    Parameters:
+    -----------
+    mrs            : MRS object
+    baseline_order : degree of polynomial (>=1)
+    ppmlim         : interval over which baseline is non-zero
+
+    Returns:
+    --------
+
+    2D numpy array
+    """
+
+    first, last = mrs.ppmlim_to_range(ppmlim)
+
+    B = []
+    x = np.zeros(mrs.numPoints, np.complex)
+    x[first:last] = np.linspace(-1, 1, last - first)
+
+    for i in range(baseline_order + 1):
+        regressor  = x**i
+        if i > 0:
+            regressor  = regress_out(regressor, B, keep_mean=False)
+
+        B.append(regressor.flatten())
+        B.append(1j * regressor.flatten())
+    B = np.asarray(B).T
+    tmp = B.copy()
+    B   = 0 * B
+    B[first:last, :] = tmp[first:last, :].copy()
+
+    return B

fsl_mrs/utils/fitting.py

@@ -13,6 +13,7 @@ import numpy as np
 from fsl_mrs.utils import models, misc
 from fsl_mrs.utils.stats import mh, vb, dist
 from fsl_mrs.utils.results import FitRes
+from fsl_mrs.utils.baseline import prepare_baseline_regressor
 
 from scipy.optimize import minimize
 
@@ -151,44 +152,6 @@ def init_FSLModel_Voigt(mrs,metab_groups,baseline,ppmlim):
 # ####################################################################################
 
 
-def prepare_baseline_regressor(mrs,baseline_order,ppmlim):
-    """
-       Complex baseline is polynomial
-
-    Parameters:
-    -----------
-    mrs            : MRS object
-    baseline_order : degree of polynomial (>=1)
-    ppmlim         : interval over which baseline is non-zero
-
-    Returns:
-    --------
-    
-    2D numpy array
-    """
-
-    first,last = mrs.ppmlim_to_range(ppmlim)
-    
-    B = []
-    x = np.zeros(mrs.numPoints,np.complex) 
-    x[first:last] = np.linspace(-1,1,last-first)
-    
-    for i in range(baseline_order+1):
-        regressor  = x**i
-        if i>0:
-            #regressor  = regressor - np.mean(regressor)
-            regressor  = misc.regress_out(regressor,B,keep_mean=False)
-            
-        B.append(regressor.flatten())
-        B.append(1j*regressor.flatten())
-    B = np.asarray(B).T
-    tmp = B.copy()
-    B   = 0*B
-    B[first:last,:] = tmp[first:last,:].copy()
-    
-    return B
-
-
 def get_bounds(num_basis,num_metab_groups,baseline_order,model,method,disableBaseline=False):
     if method == 'Newton':
         # conc

fsl_mrs/utils/synthetic/__init__.py

-from fsl_mrs.utils.synthetic.synthetic import syntheticFID
\ No newline at end of file
+from fsl_mrs.utils.synthetic.synthetic import syntheticFID
+from fsl_mrs.utils.synthetic.synthetic_from_basis import syntheticFromBasisFile

fsl_mrs/utils/synthetic/synthetic_from_basis.py

@@ -3,87 +3,95 @@
 # Author: Will Clarke <william.clarke@ndcn.ox.ac.uk>
 #         Saad Jbabdi <saad@fmrib.ox.ac.uk>
 #
-# Copyright (C) 2020 University of Oxford 
+# Copyright (C) 2020 University of Oxford
 # SHBASECOPYRIGHT
 
 import numpy as np
-from fsl_mrs.utils.misc import ts_to_ts,FIDToSpec,SpecToFID,rescale_FID
+from fsl_mrs.utils.misc import ts_to_ts, FIDToSpec, SpecToFID
 from fsl_mrs.utils import mrs_io
+from fsl_mrs.utils.baseline import prepare_baseline_regressor
+from fsl_mrs.core import MRS
+
+
 def standardConcentrations(basisNames):
     """Return standard concentrations for 1H MRS brain metabolites for those which match basis set names."""
     # These defaults are from the average of the MRS fitting challenge
-    standardconcs = {'Ala':0.60,
-                    'Asc':1.20,
-                    'Asp':2.40,
-                    'Cr':4.87,
-                    'GABA':1.20,
-                    'Glc':1.20,
-                    'Gln':3.37,
-                    'Glu':12.41,
-                    'GPC':0.74,
-                    'GSH':1.20,
-                    'Gly':1.20,
-                    'Ins':7.72,
-                    'Lac':0.60,
-                    'NAA':13.80,
-                    'NAAG':1.20,
-                    'PCho':0.85,
-                    'PCh':0.85,
-                    'PCr':4.87,
-                    'PE':1.80,
-                    'sIns':0.30,
-                    'Scyllo':0.30,
-                    'Tau':1.80}
+    standardconcs = {'Ala': 0.60,
+                     'Asc': 1.20,
+                     'Asp': 2.40,
+                     'Cr': 4.87,
+                     'GABA': 1.20,
+                     'Glc': 1.20,
+                     'Gln': 3.37,
+                     'Glu': 12.41,
+                     'GPC': 0.74,
+                     'GSH': 1.20,
+                     'Gly': 1.20,
+                     'Ins': 7.72,
+                     'Lac': 0.60,
+                     'NAA': 13.80,
+                     'NAAG': 1.20,
+                     'PCho': 0.85,
+                     'PCh': 0.85,
+                     'PCr': 4.87,
+                     'PE': 1.80,
+                     'sIns': 0.30,
+                     'Scyllo': 0.30,
+                     'Tau': 1.80}
     concs = []
     for name in basisNames:
         if name in standardconcs:
             concs.append(standardconcs[name])
         else:
             print(f'{name} not in standard concentrations. Setting to random between 1 and 5.')
-            concs.append(np.random.random()*(5.0-1.0)+1.0)
+            concs.append(np.random.random() * (5.0 - 1.0) + 1.0)
 
     return concs
 
-def syntheticFromBasisFile(basisFile,                        
-                            concentrations=None,
-                            broadening = (9.0,0.0),
-                            shifting= 0.0,
-                            baseline = [0,0],
-                            coilamps = [1.0],
-                            coilphase = [0.0],
-                            noisecovariance =[[0.1]],
-                            bandwidth = 4000.0,
-                            points = 2048):
+
+def syntheticFromBasisFile(basisFile,
+                           concentrations=None,
+                           broadening=(9.0, 0.0),
+                           shifting=0.0,
+                           baseline=None,
+                           coilamps=[1.0],
+                           coilphase=[0.0],
+                           noisecovariance=[[0.1]],
+                           bandwidth=4000.0,
+                           points=2048):
     """ Create synthetic data from a set of FSL-MRS basis files.
 
     Args:
             basisFile (str): path to directory containg basis spectra json files
             concentrations (list or dict or None, optional ): If None, standard concentrations will be used.
-                                                        If list of same length as basis spectra, then these will be used.
-                                                        Pass dict to overide standard values for specific metabolites. Key should be metabolite name.
-            broadening (list of tuples or tuple:floats, optional): Tuple containg a gamma and sigma or a list of tuples for 
-                                                                    each basis.
+                                                        If list of same length as basis spectra, then these will be
+                                                        used.
+                                                        Pass dict to overide standard values for specific metabolites.
+                                                        Key should be metabolite name.
+            broadening (list of tuples or tuple:floats, optional): Tuple containg a gamma and sigma or a list of tuples
+            for each basis.
             shifting (list of floats or float, optional): Eps shift value or a list for each basis.
             baseline (list of floats, optional): Baseline parameters. Not yet implemented
             coilamps (list of floats, optional): If multiple coils, specify magnitude scaling.
             coilphase (list of floats, optional): If multiple coils, specify phase.
             noisecovariance (list of floats, optional): N coils x N coils array of noise variance/covariance.
             bandwidth (float,optional): Bandwidth of output spectrum in Hz
-            points (int,optional): Number of points in output spectrum. 
+            points (int,optional): Number of points in output spectrum.
 
     Returns:
         FIDs: Numpy array of synthetic FIDs
         outHeader: Header suitable for loading FIDs into MRS object.
         concentrations: Final concentration scalings
     """
-    basis,names,header = mrs_io.read_basis(basisFile)
+    basis, names, header = mrs_io.read_basis(basisFile)
 
     if concentrations is None:
         concentrations = standardConcentrations(names)
-    elif isinstance(concentrations,(list,np.ndarray)):
+    elif isinstance(concentrations, (list, np.ndarray)):
         if len(concentrations) != len(names):
-            raise ValueError(f'Concentrations must have the same number of elements as basis spectra. {len(concentrations)} concentrations, {len(names)} basis spectra.')
-    elif isinstance(concentrations,dict):
+            raise ValueError(f'Concentrations must have the same number of elements as basis spectra.'
+                             f'{len(concentrations)} concentrations, {len(names)} basis spectra.')
+    elif isinstance(concentrations, dict):
         newconcs = []
         for name in names:
             if name in concentrations:
@@ -92,97 +100,129 @@ def syntheticFromBasisFile(basisFile,
                 newconcs.extend(standardConcentrations([name]))
         concentrations = newconcs
     else:
-        raise ValueError('Concentrations must be None, a list or a dict containing overides for particular metabolites.')
+        raise ValueError('Concentrations must be None, a list,'
+                         'or a dict containing overides for particular metabolites.')
 
-    
     FIDs = syntheticFromBasis(basis,
-                        header[0]['bandwidth'],
-                        concentrations,
-                        broadening = broadening,
-                        shifting= shifting,
-                        baseline = baseline,
-                        coilamps = coilamps,
-                        coilphase = coilphase,
-                        noisecovariance =noisecovariance,
-                        bandwidth = bandwidth,
-                        points = points)
-
-    outHeader = {'bandwidth':bandwidth,'centralFrequency':header[0]['centralFrequency']}
-    return FIDs,outHeader,concentrations
+                              header[0]['bandwidth'],
+                              header[0]['centralFrequency'],
+                              concentrations,
+                              broadening=broadening,
+                              shifting=shifting,
+                              baseline=baseline,
+                              coilamps=coilamps,
+                              coilphase=coilphase,
+                              noisecovariance=noisecovariance,
+                              bandwidth=bandwidth,
+                              points=points)
+
+    return FIDs, \
+        {'bandwidth': bandwidth, 'centralFrequency': header[0]['centralFrequency']}, \
+        concentrations
+
 
 def syntheticFromBasis(basis,
-                        basis_bandwidth,
-                        concentrations,
-                        broadening = (9.0,0.0),
-                        shifting= 0.0,
-                        baseline = [0,0],
-                        coilamps = [1.0],
-                        coilphase = [0.0],
-                        noisecovariance =[[0.1]],
-                        bandwidth = 4000.0,
-                        points = 2048):
+                       basis_bandwidth,
+                       basis_cf,
+                       concentrations,
+                       broadening=(9.0, 0.0),
+                       shifting=0.0,
+                       baseline=None,
+                       coilamps=[1.0],
+                       coilphase=[0.0],
+                       noisecovariance=[[0.1]],
+                       bandwidth=4000.0,
+                       points=2048):
     """ Create synthetic spectra from basis FIDs. Use syntheticFromBasisFile interface."""
     # sort out inputs
     numMetabs = basis.shape[1]
-    if len(concentrations)!=numMetabs:
+    if len(concentrations) != numMetabs:
         raise ValueError('Provide a concentration for each basis spectrum.')
 
-    if isinstance(broadening,list):
-        if len(broadening)!=numMetabs:
-            raise ValueError('Broadening values must be either a single tuple or list of tuples with the same number of elements as basis sets.')
+    if isinstance(broadening, list):
+        if len(broadening) != numMetabs:
+            raise ValueError('Broadening values must be either a single tuple,'
+                             'or list of tuples with the same number of elements as basis sets.')
         gammas = [b[0] for b in broadening]
         sigmas = [b[1] for b in broadening]
-    elif isinstance(broadening,tuple):
-        gammas = [broadening[0]]*numMetabs
-        sigmas = [broadening[1]]*numMetabs
+    elif isinstance(broadening, tuple):
+        gammas = [broadening[0]] * numMetabs
+        sigmas = [broadening[1]] * numMetabs
     else:
-        raise ValueError('Broadening values must be either a single tuple or list of tuples with the same number of elements as basis sets.')
+        raise ValueError('Broadening values must be either a single tuple,'
+                         ' or list of tuples with the same number of elements as basis sets.')
 
-    if isinstance(shifting,list):
-        if len(shifting)!=numMetabs:
-            raise ValueError('shifting values must be either a float or list with the same number of elements as basis sets.')
+    if isinstance(shifting, list):
+        if len(shifting) != numMetabs:
+            raise ValueError('shifting values must be either a float.'
+                             ' or list with the same number of elements as basis sets.')
         eps = shifting
-    elif isinstance(shifting,float):
-        eps = [shifting]*numMetabs
+    elif isinstance(shifting, float):
+        eps = [shifting] * numMetabs
     else:
-        raise ValueError('shifting values must be either a float or list with the same number of elements as basis sets.')
+        raise ValueError('shifting values must be either a float,'
+                         ' or list with the same number of elements as basis sets.')
 
     # Form noise vectors
     ncoils = len(coilamps)
     noisecovariance = np.asarray(noisecovariance)
     if len(coilphase) != ncoils:
         raise ValueError('Length of coilamps and coilphase must match.')
-    if noisecovariance.shape != (ncoils,ncoils):
+    if noisecovariance.shape != (ncoils, ncoils):
         raise ValueError('noisecovariance must be ncoils x ncoils.')
 
-    noise = np.random.multivariate_normal(np.zeros((ncoils)), noisecovariance, points) + 1j*np.random.multivariate_normal(np.zeros((ncoils)), noisecovariance, points)
+    noise = np.random.multivariate_normal(np.zeros((ncoils)), noisecovariance, points) \
+        + 1j * np.random.multivariate_normal(np.zeros((ncoils)), noisecovariance, points)
 
     # Interpolate basis
-    dwelltime = 1/bandwidth
-    basis_dwelltime= 1/basis_bandwidth
-    basis = ts_to_ts(basis,basis_dwelltime,dwelltime,points)
+    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)
-    dwellTime   = 1/bandwidth
-    timeAxis    = np.linspace(dwellTime,
-                              dwellTime*points,
-                              points)
-    for b,c,e,g,s in zip(basis.T,concentrations,eps,gammas,sigmas):
-        tmp   = b*np.exp(-(1j*e+g+timeAxis*s**2)*timeAxis)
+    baseFID    = np.zeros((points), np.complex128)
+    dwellTime  = 1 / bandwidth
+    timeAxis   = np.linspace(dwellTime,
+                             dwellTime * points,
+                             points)
+    for b, c, e, g, s in zip(basis.T, concentrations, eps, gammas, sigmas):
+        tmp   = b * np.exp(-(1j * e + g + timeAxis * s**2) * timeAxis)
         M     = FIDToSpec(tmp)
-        baseFID += SpecToFID(M*c)
-        
+        baseFID += SpecToFID(M * c)
 
     # Add baseline
-    # TO DO
+    if baseline is not None:
+        baseFID += calculate_baseline(baseline, bandwidth, basis_cf, points)
 
-    # Add noise and write tot output list
+    # Add noise and write to output list
     FIDs = []
-    for cDx,(camp,cphs) in enumerate(zip(coilamps,coilphase)):
-        FIDs.append((camp*np.exp(1j*cphs)*baseFID)+noise[:,cDx])
+    for cDx, (camp, cphs) in enumerate(zip(coilamps, coilphase)):
+        FIDs.append((camp * np.exp(1j * cphs) * baseFID) + noise[:, cDx])
     FIDs = np.asarray(FIDs).T
     FIDs = np.squeeze(FIDs)
 
-    return FIDs
\ No newline at end of file
+    return FIDs
+
+
+def calculate_baseline(bline_param, bandwidth, cf, points):
+    """Create a synthetic polynomial baseline specified by pairs of constants"""
+
+    order = len(bline_param) - 1
+    for opair in bline_param:
+        if len(opair) != 2:
+            raise ValueError('The baseline parameter must be a list of 2-tuples'
+                             ' containing real and imaginary coeeficients.')
+
+    dummy_mrs = MRS(FID=np.zeros((points,)),
+                    cf=cf,
+                    bw=bandwidth,
+                    nucleus='1H')
+
+    ppmlim = None
+    B = prepare_baseline_regressor(dummy_mrs, order, ppmlim)
+
+    bline_param = np.asarray(bline_param).flatten()
+
+    spec = np.sum(B * bline_param, axis=1)
+    return SpecToFID(spec)