fsl_mrs 13.7 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
Saad Jbabdi's avatar
Saad Jbabdi committed
17

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



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

26
27
28
29
30
31
32
33
34

    # 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
35
    
36
    required     = p.add_argument_group('required arguments')
Saad Jbabdi's avatar
Saad Jbabdi committed
37
    fitting_args = p.add_argument_group('fitting options')
38
    optional     = p.add_argument_group('additional options')
Saad Jbabdi's avatar
Saad Jbabdi committed
39
40

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

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


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

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

145
146
147
148
149
150
151
152

    # 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
153
154
    #######  Do the work #######

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


162
163
    FID,dataheader            = mrs_io.read_FID(args.data)
    basis, names, basisheader = mrs_io.read_basis(args.basis)
164
        
165
166
167
168
    if args.h2o is not None:
        H2O,_                 = mrs_io.read_FID(args.h2o)
    else:
        H2O = None
169
      
170
171
172
173
174
    # 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
175
        if args.verbose:
176
            print(f'      Detected central frequency in header info cf = {cf:0.6f} MHz')
177
    else:
Saad Jbabdi's avatar
Saad Jbabdi committed
178
        raise(Exception('Cannot determine central frequency. Please either set it or include it in data header'))
179
180
181
182
183

    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
184
        if args.verbose:
185
            print(f'      Detected bandwidth in header info bw = {bw:0.1f} Hz')
186
    else:
Saad Jbabdi's avatar
Saad Jbabdi committed
187
        raise(Exception('Cannot determine bandwidth. Please either set it or include it in data header'))
Saad Jbabdi's avatar
Saad Jbabdi committed
188
189
190
191
192
193

    # 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]
194
        
Saad Jbabdi's avatar
Saad Jbabdi committed
195
    # Instantiate MRS object
Saad Jbabdi's avatar
Saad Jbabdi committed
196
    MRSargs = {'FID':FID,'basis':basis,'basis_hdr':basisheader,'names':names,'H2O':H2O,'cf':cf,'bw':bw}
197
198
    mrs     = MRS(**MRSargs)
    
199
200
201
202
203
    # 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
204
        conjugated = mrs.check_FID(repair=True)
205
206
207
208
209
210
211
212
        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
213
        conjugated = mrs.check_Basis(repair=True)
214
215
216
        if args.verbose:
            if conjugated == 1:
                warnings.warn('Basis has been checked and conjugated. Please check!',UserWarning)
217
218
219
220
    
    # Rescale FID, H2O and basis to have nice range
    if not args.no_rescale:
        mrs.rescaleForFitting()
Saad Jbabdi's avatar
Saad Jbabdi committed
221

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

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

240

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


    # Parse metabolite groups
    metab_groups = args.metab_groups
249
250
251
252
253
254
255
256
257
258
259
260
261
    if isinstance(metab_groups,list):
        if isinstance(metab_groups[0],str):
            tmp = [0]*mrs.numBasis
            grpcounter = 0
            for n in metab_groups:
                grpcounter += 1
                tmp[mrs.names.index(n)] = grpcounter
            metab_groups = tmp
        elif isinstance(metab_groups[0],int):
            if metab_groups == [0]:
                metab_groups = [0]*mrs.numBasis
            elif len(metab_groups) != mrs.numBasis:
                raise(Exception('Found {} metab_groups but there are {} basis functions'.format(len(metab_groups),mrs.numBasis)))
Saad Jbabdi's avatar
Saad Jbabdi committed
262
263
    elif metab_groups == 0:
        metab_groups = [0]*mrs.numBasis
264
265
266

    # Include Macromolecules? These should have their own metab groups
    if args.add_MM is not None:
Saad Jbabdi's avatar
Saad Jbabdi committed
267
268
        if not args.verbose:
            print('Adding macromolecules')
269
270
271
272
273
        nMM = mrs.add_MM_peaks()
        G   = [i+max(metab_groups)+1 for i in range(nMM)]
        metab_groups += G


Saad Jbabdi's avatar
Saad Jbabdi committed
274
    Fitargs = {'ppmlim':ppmlim,
275
276
277
               'method':args.algo,'baseline_order':args.baseline_order,
               'metab_groups':metab_groups}

Saad Jbabdi's avatar
Saad Jbabdi committed
278
279
280
281
282
283
    if args.verbose:
        print(mrs)
        print('Fitting args:')
        print(Fitargs)

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

285
286
287
288
    # Quantification
    # Echo time
    if args.TE is not None:
        echotime = args.TE*1E-3
289
290
291
292
293
294
295
    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
296
297
298
299
    else:
        echotime = None
    # Internal and Water quantification if requested
    if (mrs.H2O is None) or (echotime is None):
300
        res.calculateConcScaling(mrs,referenceMetab=args.internal_ref)
301
302
303
    else:        
        tfrac = {'WM':args.tissue_fractions[0],'GM':args.tissue_fractions[1],'CSF':args.tissue_fractions[2]}
        res.calculateConcScaling(mrs,
304
                                referenceMetab=args.internal_ref,
305
306
307
308
309
310
311
                                waterRefFID=mrs.H2O,
                                tissueFractions=tfrac,
                                TE=echotime,
                                verbose=args.verbose)
    # Combine metabolites.
    if args.combine is not None:
        res.combine(args.combine)
Saad Jbabdi's avatar
Saad Jbabdi committed
312
313
314
315
316
317
318
319
320
    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))
321
322


323
    res.to_file(filename=os.path.join(args.output,'results_table.csv'),what='concentrations')
Saad Jbabdi's avatar
Saad Jbabdi committed
324
325
    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')
Saad Jbabdi's avatar
Saad Jbabdi committed
326

Saad Jbabdi's avatar
Saad Jbabdi committed
327
328
329
330
331
332
    # Save image of MRS voxel
    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)
            fig.savefig(os.path.join(args.output,'voxel_location.png'))
Saad Jbabdi's avatar
Saad Jbabdi committed
333

334
335
336
337
    # Create short HTML report
    #fig = plotting.plotly_fit(mrs,res,ppmlim=ppmlim,proj='abs')
    #plotly.io.write_html(fig, file=os.path.join(args.output,'short_report.html'))
    
Saad Jbabdi's avatar
Saad Jbabdi committed
338
339
340
341

    # Save quick summary figure
    report.fitting_summary_fig(mrs,res,
                               filename=os.path.join(args.output,'fit_summary.png'))
342
343
    
    
Saad Jbabdi's avatar
Saad Jbabdi committed
344
    # Create interactive HTML report
Saad Jbabdi's avatar
Saad Jbabdi committed
345
346
347
348
349
350
351
    if args.report:
        report.create_report(mrs,res,
                             filename=os.path.join(args.output,'report.html'),
                             fidfile=args.data,
                             basisfile=args.basis,
                             h2ofile=args.h2o,
                             date=datetime.datetime.now().strftime("%Y-%m-%d %H:%M"))
352
                         
353

Saad Jbabdi's avatar
Saad Jbabdi committed
354
355
356
357
     
    if args.verbose:
        print('\n\n\nDone.')

358
359
360
361
362
def str_or_int_arg(x):    
    try:                           
        return int(x)
    except:
        return x
Saad Jbabdi's avatar
Saad Jbabdi committed
363
364
365
        
if __name__ == '__main__':
    main()