fsl_mrs 13.6 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
13
# Quick imports
#import argparse
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
81
82
    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
83
    optional.add_argument('--report',action="store_true",
84
                          help='output html report')
Saad Jbabdi's avatar
Saad Jbabdi committed
85
    optional.add_argument('--verbose',action="store_true",
86
                          help='spit out verbose info')
Saad Jbabdi's avatar
Saad Jbabdi committed
87
88
    optional.add_argument('--phase_correct',action="store_true",
                          help='do phase correction')
Saad Jbabdi's avatar
Saad Jbabdi committed
89
    optional.add_argument('--overwrite',action="store_true",
90
                          help='overwrite existing output folder')
91
92
93
94
95
    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)
96
    optional.add_argument('--no_rescale',action="store_false",help='Forbid rescaling of FID/basis/H2O.')
97
    optional.add('--config', required=False, is_config_file=True, help='configuration file')
Saad Jbabdi's avatar
Saad Jbabdi committed
98
99


100
101
    #hide_args([h1,h2])
    
Saad Jbabdi's avatar
Saad Jbabdi committed
102
103
104
    # Parse command-line arguments
    args = p.parse_args()

105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
    # 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
123
124
    if not args.verbose:
        warnings.filterwarnings("ignore")
125
126

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

143
144
145
146
147
148
149
150

    # 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
151
152
    #######  Do the work #######

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


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

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

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

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

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

238

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


    # Parse metabolite groups
    metab_groups = args.metab_groups
247
248
249
250
251
252
253
254
255
256
257
258
259
    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
260
261
    elif metab_groups == 0:
        metab_groups = [0]*mrs.numBasis
262
263
264

    # Include Macromolecules? These should have their own metab groups
    if args.add_MM is not None:
Saad Jbabdi's avatar
Saad Jbabdi committed
265
266
        if not args.verbose:
            print('Adding macromolecules')
267
268
269
270
271
        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
272
    Fitargs = {'ppmlim':ppmlim,
273
274
275
               'method':args.algo,'baseline_order':args.baseline_order,
               'metab_groups':metab_groups}

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

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

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


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

Saad Jbabdi's avatar
Saad Jbabdi committed
325
326
327
328
329
330
    # 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
331

332
333
334
335
    # 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
336
337
338
339

    # Save quick summary figure
    report.fitting_summary_fig(mrs,res,
                               filename=os.path.join(args.output,'fit_summary.png'))
340
341
    
    
Saad Jbabdi's avatar
Saad Jbabdi committed
342
    # Create interactive HTML report
Saad Jbabdi's avatar
Saad Jbabdi committed
343
344
345
346
347
348
349
350
    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,
                             outdir=args.output,
                             date=datetime.datetime.now().strftime("%Y-%m-%d %H:%M"))
351
                         
352

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

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