fsl_mrs 15 KB
Newer Older
Saad Jbabdi's avatar
Saad Jbabdi committed
1
2
3
4
5
#!/usr/bin/env python

# fsl_mrs - wrapper script for MRS fitting
#
# Author: Saad Jbabdi <saad@fmrib.ox.ac.uk>
6
#         William Carke <william.clarke@ndcn.ox.ac.uk>
Saad Jbabdi's avatar
Saad Jbabdi committed
7
8
9
10
#
# Copyright (C) 2019 University of Oxford 
# SHBASECOPYRIGHT

11
12
# Quick imports
#import argparse
13
from fsl_mrs.aux import configargparse
Saad Jbabdi's avatar
Saad Jbabdi committed
14

15
16
from fsl_mrs import __version__
from fsl_mrs.utils.splash import splash
17
import json
Saad Jbabdi's avatar
Saad Jbabdi committed
18

19
# NOTE!!!! THERE ARE MORE IMPORTS IN THE CODE BELOW (AFTER ARGPARSING)
Saad Jbabdi's avatar
Saad Jbabdi committed
20
21
22
23
24



def main():
    # Parse command-line arguments
25
    p = configargparse.ArgParser(add_config_file_help=False,description="FSL Magnetic Resonance Spectroscopy Wrapper Script")
Saad Jbabdi's avatar
Saad Jbabdi committed
26

27
28
29
30
31
32
33
34
35

    # utility for hiding certain arguments
    def hide_args(arglist):
        for action in arglist:
            action.help=p.SUPPRESS

        
    #p = argparse.ArgumentParser(description='FSL Magnetic Resonance Spectroscopy Tool')
    p.add_argument('-v','--version', action='version', version=__version__)
Saad Jbabdi's avatar
Saad Jbabdi committed
36
    
37
    required     = p.add_argument_group('required arguments')
Saad Jbabdi's avatar
Saad Jbabdi committed
38
    fitting_args = p.add_argument_group('fitting options')
39
    optional     = p.add_argument_group('additional options')
Saad Jbabdi's avatar
Saad Jbabdi committed
40
41

    # REQUIRED ARGUMENTS
42
    required.add_argument('--data',
43
                          required=True,type=str,metavar='<str>',
Saad Jbabdi's avatar
Saad Jbabdi committed
44
                          help='input FID file')
45
    required.add_argument('--basis',
Saad Jbabdi's avatar
Saad Jbabdi committed
46
                          required=True,type=str,metavar='<str>',
47
                          help='.BASIS file or folder containing basis spectra (will read all files within)')
48
    required.add_argument('--output',
Saad Jbabdi's avatar
Saad Jbabdi committed
49
50
51
52
53
                          required=True,type=str,metavar='<str>',
                          help='output folder')
    
    # FITTING ARGUMENTS
    fitting_args.add_argument('--algo',default='Newton',type=str,
54
                   help='algorithm [Newton (fast, default) or MH (slow)]')
Saad Jbabdi's avatar
Saad Jbabdi committed
55
56
    fitting_args.add_argument('--ignore',type=str,nargs='+',metavar='METAB',
                   help='ignore certain metabolites [repeatable]')
57
58
    fitting_args.add_argument('--keep',type=str,nargs='+',metavar='METAB',
                   help='only keep these metabolites')
Saad Jbabdi's avatar
Saad Jbabdi committed
59
    fitting_args.add_argument('--combine',type=str,nargs='+',action='append',metavar='METAB',
60
61
62
                   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))')
Saad Jbabdi's avatar
Saad Jbabdi committed
63
    fitting_args.add_argument('--h2o',default=None,type=str,metavar='H2O',
Saad Jbabdi's avatar
Saad Jbabdi committed
64
                   help='input .H2O file for quantification')
65
    fitting_args.add_argument('--baseline_order',default=2,type=int,metavar=('ORDER'),
66
                   help='order of baseline polynomial (default=2,-1 disables)')
67
68
    fitting_args.add_argument('--metab_groups',default=0,nargs='+',type=str_or_int_arg,
                              help="metabolite groups: list of groups or list of names for indept groups.")
69
    fitting_args.add_argument('--add_MM',action="store_true",
70
                              help="include default macromolecule peaks")
71
72
    fitting_args.add_argument('--lorentzian',action="store_true",
                              help="Enable purely lorentzian broadening (default is Voigt)")                          
73

74
    
Saad Jbabdi's avatar
Saad Jbabdi committed
75
    # ADDITONAL OPTIONAL ARGUMENTS
Saad Jbabdi's avatar
Saad Jbabdi committed
76
77
    optional.add_argument('--t1',type=str,default=None,metavar='IMAGE',
                          help='structural image (for report)')
78
79
    optional.add_argument('--TE',type=float,default=None,metavar='TE',
                          help='Echo time for relaxation correction (ms)')
80
    optional.add_argument('--tissue_frac',type=tissue_frac_arg,action=TissueFracAction,nargs='+',default=None,metavar='WM GM CSF OR json',
81
                          help='Fractional tissue volumes for WM, GM, CSF or json segmentation file. Defaults to pure water scaling.')
82
83
    optional.add_argument('--internal_ref',type=str,default=['Cr','PCr'],nargs='+',
                          help='Metabolite(s) used as an internal reference. Defaults to tCr (Cr+PCr).')
84
85
86
87
    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')
Saad Jbabdi's avatar
Saad Jbabdi committed
88
    optional.add_argument('--report',action="store_true",
89
                          help='output html report')
Saad Jbabdi's avatar
Saad Jbabdi committed
90
    optional.add_argument('--verbose',action="store_true",
91
                          help='spit out verbose info')
Saad Jbabdi's avatar
Saad Jbabdi committed
92
93
    optional.add_argument('--phase_correct',action="store_true",
                          help='do phase correction')
Saad Jbabdi's avatar
Saad Jbabdi committed
94
    optional.add_argument('--overwrite',action="store_true",
95
                          help='overwrite existing output folder')
96
97
98
99
100
    optional.add_argument('--conj_fid',dest='conjfid',action="store_true",help='Force conjugation of FID')
    optional.add_argument('--no_conj_fid',dest='conjfid',action="store_false",help='Forbid automatic conjugation of FID')
    optional.add_argument('--conj_basis',dest='conjbasis',action="store_true",help='Force conjugation of basis')
    optional.add_argument('--no_conj_basis',dest='conjbasis',action="store_false",help='Forbid automatic conjugation of basis')
    optional.set_defaults(conjfid=None,conjbasis=None)
101
    optional.add_argument('--no_rescale',action="store_true",help='Forbid rescaling of FID/basis/H2O.')
102
    optional.add('--config', required=False, is_config_file=True, help='configuration file')
Saad Jbabdi's avatar
Saad Jbabdi committed
103
104


105
106
    #hide_args([h1,h2])
    
Saad Jbabdi's avatar
Saad Jbabdi committed
107
108
    # Parse command-line arguments
    args = p.parse_args()
109
    
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
    # 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 fitting
    from fsl_mrs.utils import plotting
    from fsl_mrs.utils import misc
    import datetime
    import plotly
    # ######################################################
Saad Jbabdi's avatar
Saad Jbabdi committed
128
129
    if not args.verbose:
        warnings.filterwarnings("ignore")
130
131

    
Saad Jbabdi's avatar
Saad Jbabdi committed
132
133
134
135
136
137
138
139
140
141
142
143
    # 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)
Saad Jbabdi's avatar
Saad Jbabdi committed
144
            os.makedirs(args.output,exist_ok=True)
Saad Jbabdi's avatar
Saad Jbabdi committed
145
    else:
Saad Jbabdi's avatar
Saad Jbabdi committed
146
        os.makedirs(args.output,exist_ok=True)
Saad Jbabdi's avatar
Saad Jbabdi committed
147

148
149
150
151
152
153
154
155

    # 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())
        
        
Saad Jbabdi's avatar
Saad Jbabdi committed
156
157
    #######  Do the work #######

Saad Jbabdi's avatar
Saad Jbabdi committed
158
159
160
161
162
163
164
    # Read data/h2o/basis
    if args.verbose:        
        print('--->> Read input data and basis\n')
        print('  {}'.format(args.data))
        print('  {}\n'.format(args.basis))


165
166
    FID,dataheader            = mrs_io.read_FID(args.data)
    basis, names, basisheader = mrs_io.read_basis(args.basis)
167
        
168
169
170
171
    if args.h2o is not None:
        H2O,_                 = mrs_io.read_FID(args.h2o)
    else:
        H2O = None
172
      
173
174
175
176
177
    # Collect useful info
    if args.central_frequency is not None:
        cf = args.central_frequency
    elif dataheader['centralFrequency'] is not None:
        cf = dataheader['centralFrequency']
Saad Jbabdi's avatar
Saad Jbabdi committed
178
        if args.verbose:
179
            print(f'      Detected central frequency in header info cf = {cf:0.6f} MHz')
180
    else:
Saad Jbabdi's avatar
Saad Jbabdi committed
181
        raise(Exception('Cannot determine central frequency. Please either set it or include it in data header'))
182
183
184
185
186

    if args.dwell_time is not None:
        bw = 1/args.dwell_time
    elif dataheader['bandwidth'] is not None:
        bw = dataheader['bandwidth']
Saad Jbabdi's avatar
Saad Jbabdi committed
187
        if args.verbose:
188
            print(f'      Detected bandwidth in header info bw = {bw:0.1f} Hz')
189
    else:
Saad Jbabdi's avatar
Saad Jbabdi committed
190
        raise(Exception('Cannot determine bandwidth. Please either set it or include it in data header'))
Saad Jbabdi's avatar
Saad Jbabdi committed
191
192
193
194
195
196

    # Fix case where basis file contains no header info (e.g. .RAW)
    if basisheader is None:
        basisheader = {'bandwidth':bw,'dwelltime':1/bw,'centralFrequency':cf}
    else:
        basisheader = basisheader[0]
197
        
Saad Jbabdi's avatar
Saad Jbabdi committed
198
    # Instantiate MRS object
Saad Jbabdi's avatar
Saad Jbabdi committed
199
    MRSargs = {'FID':FID,'basis':basis,'basis_hdr':basisheader,'names':names,'H2O':H2O,'cf':cf,'bw':bw}
200
201
    mrs     = MRS(**MRSargs)
    
202
203
204
205
206
    # Check the FID and basis / conjugate
    if args.conjfid is not None:
        if args.conjfid:
            mrs.conj_FID()
    else:
William Clarke's avatar
William Clarke committed
207
        conjugated = mrs.check_FID(repair=True)
208
209
210
211
212
213
214
215
        if args.verbose:
            if conjugated == 1:
                warnings.warn('FID has been checked and conjugated. Please check!',UserWarning)

    if args.conjbasis is not None:
        if args.conjbasis:
            mrs.conj_Basis()
    else:
William Clarke's avatar
William Clarke committed
216
        conjugated = mrs.check_Basis(repair=True)
217
218
219
        if args.verbose:
            if conjugated == 1:
                warnings.warn('Basis has been checked and conjugated. Please check!',UserWarning)
220
221
222
223
    
    # Rescale FID, H2O and basis to have nice range
    if not args.no_rescale:
        mrs.rescaleForFitting()
Saad Jbabdi's avatar
Saad Jbabdi committed
224

225
    # Do phase correction
Saad Jbabdi's avatar
Saad Jbabdi committed
226
227
228
229
    if args.phase_correct:
        if args.verbose:
            print('--->>  Phase correction\n')
        mrs.FID  = misc.phase_correct(mrs,mrs.FID)
230
        mrs.Spec = misc.FIDToSpec(mrs.FID)
Saad Jbabdi's avatar
Saad Jbabdi committed
231
              
Saad Jbabdi's avatar
Saad Jbabdi committed
232
        
233
    # Keep/Ignore metabolites
234
    mrs.keep(args.keep)
235
    mrs.ignore(args.ignore)    
236

Saad Jbabdi's avatar
Saad Jbabdi committed
237
238
239
    # Do the fitting here
    if args.verbose:
        print('--->> Start fitting\n\n')
240
        print('    Algorithm = [{}]\n'.format(args.algo))
Saad Jbabdi's avatar
Saad Jbabdi committed
241
242
        start = time.time()

243

244
245
246
    ppmlim=args.ppmlim
    if ppmlim is not None:
        ppmlim=(ppmlim[0],ppmlim[1])
247
248
249
250
 


    # Parse metabolite groups
251
252
    metab_groups = misc.parse_metab_groups(mrs,args.metab_groups)
    
253
254

    # Include Macromolecules? These should have their own metab groups
255
    if args.add_MM:
Saad Jbabdi's avatar
Saad Jbabdi committed
256
257
        if not args.verbose:
            print('Adding macromolecules')
258
        nMM = mrs.add_MM_peaks(gamma=10,sigma=20)
259
260
261
        G   = [i+max(metab_groups)+1 for i in range(nMM)]
        metab_groups += G

262
263
264
265
266
267
268
    if args.lorentzian:
        Fitargs = {'ppmlim':ppmlim,
               'method':args.algo,'baseline_order':args.baseline_order,
               'metab_groups':metab_groups,
               'model':'lorentzian'}
    else:
        Fitargs = {'ppmlim':ppmlim,
269
               'method':args.algo,'baseline_order':args.baseline_order,
270
271
               'metab_groups':metab_groups,
               'model':'voigt'}
272

Saad Jbabdi's avatar
Saad Jbabdi committed
273
274
275
276
277
278
    if args.verbose:
        print(mrs)
        print('Fitting args:')
        print(Fitargs)

    res = fitting.fit_FSLModel(mrs,**Fitargs)
279

280
281
282
283
    # Quantification
    # Echo time
    if args.TE is not None:
        echotime = args.TE*1E-3
284
285
286
287
288
289
290
    elif 'meta' in basisheader:
        if 'TE' in basisheader['meta']:
            echotime = basisheader['meta']['TE']
            if echotime > 1.0: # Assume in ms.
                echotime *= 1E-3
            else:
                echotime = None
291
292
    elif 'TE' in dataheader:
        echotime = dataheader['TE']
293
294
295
296
    else:
        echotime = None
    # Internal and Water quantification if requested
    if (mrs.H2O is None) or (echotime is None):
297
298
        if echotime is None:
            warnings.warn('H2O file provided but could not determine TE: no absolute quantification will be performed.',UserWarning)
299
        res.calculateConcScaling(mrs,referenceMetab=args.internal_ref)
300
    elif args.tissue_frac is None:
301
302
303
304
305
306
        res.calculateConcScaling(mrs,
                                referenceMetab=args.internal_ref,
                                waterRefFID=mrs.H2O,
                                tissueFractions=None,
                                TE=echotime,
                                verbose=args.verbose)
307
    else:
308
        res.calculateConcScaling(mrs,
309
                                referenceMetab=args.internal_ref,
310
                                waterRefFID=mrs.H2O,
311
                                tissueFractions=args.tissue_frac,
312
313
314
315
316
                                TE=echotime,
                                verbose=args.verbose)
    # Combine metabolites.
    if args.combine is not None:
        res.combine(args.combine)
Saad Jbabdi's avatar
Saad Jbabdi committed
317
318
319
320
321
322
323
324
325
    stop = time.time()

    
    # Report on the fitting
    if args.verbose:
        print('    Fitting lasted          : {:.3f} secs.\n'.format(stop-start))
    # Save output files
    if args.verbose:
        print('--->> Saving output files to {}\n'.format(args.output))
326
327


328
329
    res.to_file(filename=os.path.join(args.output,'summary.csv'),what='summary')
    res.to_file(filename=os.path.join(args.output,'concentrations.csv'),what='concentrations')
Saad Jbabdi's avatar
Saad Jbabdi committed
330
331
    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')
332
333
334
    if args.algo == 'MH':
        res.to_file(filename=os.path.join(args.output,'concentration_samples.csv'),what='concentrations-mh')
        res.to_file(filename=os.path.join(args.output,'all_samples.csv'),what='parameters-mh')
Saad Jbabdi's avatar
Saad Jbabdi committed
335

Saad Jbabdi's avatar
Saad Jbabdi committed
336
    # Save image of MRS voxel
337
    location_fig = None
Saad Jbabdi's avatar
Saad Jbabdi committed
338
339
340
341
    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)
342
            fig.tight_layout()
343
344
            location_fig = os.path.join(args.output,'voxel_location.png')
            fig.savefig(location_fig,bbox_inches='tight',facecolor='k')
Saad Jbabdi's avatar
Saad Jbabdi committed
345
346
347
348

    # Save quick summary figure
    report.fitting_summary_fig(mrs,res,
                               filename=os.path.join(args.output,'fit_summary.png'))
349
350
    
    
Saad Jbabdi's avatar
Saad Jbabdi committed
351
    # Create interactive HTML report
Saad Jbabdi's avatar
Saad Jbabdi committed
352
353
354
355
356
    if args.report:
        report.create_report(mrs,res,
                             filename=os.path.join(args.output,'report.html'),
                             fidfile=args.data,
                             basisfile=args.basis,
357
358
359
                             h2ofile=args.h2o,                             
                             date=datetime.datetime.now().strftime("%Y-%m-%d %H:%M"),
                             location_fig = location_fig)
360
                         
361

Saad Jbabdi's avatar
Saad Jbabdi committed
362
363
364
365
     
    if args.verbose:
        print('\n\n\nDone.')

366
367
368
369
370
def str_or_int_arg(x):    
    try:                           
        return int(x)
    except:
        return x
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386

def tissue_frac_arg(x):
    try:
        with open(x) as jsonFile:
            jsonString = jsonFile.read()                   
        return json.loads(jsonString)
    except:
        return float(x)

class TissueFracAction(configargparse.Action):
    """Sort out tissue fraction types. Should return dict"""
    def __call__(self, parser, namespace, values, option_string=None):
        if isinstance(values[0],dict):
            setattr(namespace,self.dest, values[0])
        else:
            setattr(namespace,self.dest, {'WM':values[0],'GM':values[1],'CSF':values[2]})
Saad Jbabdi's avatar
Saad Jbabdi committed
387
388
389
        
if __name__ == '__main__':
    main()