GM-WM-CSF --> WM-GM-CSF

scripts/src/dwi2response/msmt_5tt.py

@@ -1,25 +1,25 @@
 def initParser(subparsers, base_parser):
   import argparse
-  parser = subparsers.add_parser('msmt_5tt', parents=[base_parser], help='Derive MSMT CSD responses based on a co-registered 5TT image')
+  parser = subparsers.add_parser('msmt_5tt', parents=[base_parser], help='Derive MSMT-CSD responses based on a co-registered 5TT image')
   arguments = parser.add_argument_group('Positional arguments specific to the \'msmt_5tt\' algorithm')
   arguments.add_argument('in_5tt', help='Input co-registered 5TT image')
-  arguments.add_argument('out_gm', help='Output GM response text file')
   arguments.add_argument('out_wm', help='Output WM response text file')
+  arguments.add_argument('out_gm', help='Output GM response text file')
   arguments.add_argument('out_csf', help='Output CSF response text file')
   options = parser.add_argument_group('Options specific to the \'msmt_5tt\' algorithm')
   options.add_argument('-dirs', help='Manually provide the fibre direction in each voxel (a tensor fit will be used otherwise)')
-  options.add_argument('-fa', type=float, default=0.2, help='Upper fractional anisotropy threshold for isotropic tissue (i.e. GM and CSF) voxel selection')
+  options.add_argument('-fa', type=float, default=0.2, help='Upper fractional anisotropy threshold for GM and CSF voxel selection')
   options.add_argument('-pvf', type=float, default=0.95, help='Partial volume fraction threshold for tissue voxel selection')
-  options.add_argument('-wm_algo', metavar='algorithm', default='tournier', help='dwi2response algorithm to use for white matter single-fibre voxel selection')
+  options.add_argument('-wm_algo', metavar='algorithm', default='tournier', help='dwi2response algorithm to use for WM single-fibre voxel selection')
   parser.set_defaults(algorithm='msmt_5tt')
   parser.set_defaults(single_shell=False)
 
 
 
 def checkOutputFiles():
   import lib.app
-  lib.app.checkOutputFile(lib.app.args.out_gm)
   lib.app.checkOutputFile(lib.app.args.out_wm)
+  lib.app.checkOutputFile(lib.app.args.out_gm)
   lib.app.checkOutputFile(lib.app.args.out_csf)
 
 
@@ -56,7 +56,7 @@ def execute():
   # Get shell information
   shells = [ int(round(float(x))) for x in getHeaderInfo('dwi.mif', 'shells').split() ]
   if len(shells) < 3:
-    warnMessage('Less than three b-value shells; response functions will not be applicable in MSMT CSD algorithm')
+    warnMessage('Less than three b-value shells; response functions will not be applicable in MSMT-CSD algorithm')
 
   # Get lmax information (if provided)
   wm_lmax = [ ]
@@ -75,9 +75,9 @@ def execute():
     shutil.copy('vector.mif', 'dirs.mif')
   runCommand('mrtransform 5tt.mif 5tt_regrid.mif -template fa.mif -interp linear')
 
-  # Tissue masks
-  runCommand('mrconvert 5tt_regrid.mif - -coord 3 0 -axes 0,1,2 | mrcalc - ' + str(lib.app.args.pvf) + ' -gt fa.mif ' + str(lib.app.args.fa) + ' -lt -mult mask.mif -mult gm_mask.mif')
+  # Basic tissue masks
   runCommand('mrconvert 5tt_regrid.mif - -coord 3 2 -axes 0,1,2 | mrcalc - ' + str(lib.app.args.pvf) + ' -gt mask.mif -mult wm_mask.mif')
+  runCommand('mrconvert 5tt_regrid.mif - -coord 3 0 -axes 0,1,2 | mrcalc - ' + str(lib.app.args.pvf) + ' -gt fa.mif ' + str(lib.app.args.fa) + ' -lt -mult mask.mif -mult gm_mask.mif')
   runCommand('mrconvert 5tt_regrid.mif - -coord 3 3 -axes 0,1,2 | mrcalc - ' + str(lib.app.args.pvf) + ' -gt fa.mif ' + str(lib.app.args.fa) + ' -lt -mult mask.mif -mult csf_mask.mif')
 
   # Revise WM mask to only include single-fibre voxels
@@ -88,14 +88,14 @@ def execute():
   runCommand('dwi2response -quiet -tempdir ' + lib.app.tempDir + recursive_cleanup_option + ' ' + lib.app.args.wm_algo + ' dwi.mif wm_ss_response.txt -mask wm_mask.mif -voxels wm_sf_mask.mif')
 
   # Check for empty masks
-  gm_voxels  = int(getImageStat('gm_mask.mif',    'count', 'gm_mask.mif'))
   wm_voxels  = int(getImageStat('wm_sf_mask.mif', 'count', 'wm_sf_mask.mif'))
+  gm_voxels  = int(getImageStat('gm_mask.mif',    'count', 'gm_mask.mif'))
   csf_voxels = int(getImageStat('csf_mask.mif',   'count', 'csf_mask.mif'))
   empty_masks = [ ]
-  if not gm_voxels:
-    empty_masks.append('GM')
   if not wm_voxels:
     empty_masks.append('WM')
+  if not gm_voxels:
+    empty_masks.append('GM')
   if not csf_voxels:
     empty_masks.append('CSF')
   if empty_masks:
@@ -111,44 +111,44 @@ def execute():
   # For each of the three tissues, generate a multi-shell response
   # Since here we're guaranteeing that GM and CSF will be isotropic in all shells, let's use mrstats rather than sh2response (seems a bit weird passing a directions file to sh2response with lmax=0...)
 
-  gm_responses  = [ ]
   wm_responses  = [ ]
+  gm_responses  = [ ]
   csf_responses = [ ]
   max_length = 0
 
   for index, b in enumerate(shells):
     dwi_path = 'dwi_b' + str(b) + '.mif'
-    mean_dwi_path = 'dwi_b' + str(b) + '_mean.mif'
-    # dwiextract will now yield a 4D image, even if there's only a single volume in a shell
+    # dwiextract will yield a 4D image, even if there's only a single volume in a shell
     runCommand('dwiextract dwi.mif -shell ' + str(b) + ' ' + dwi_path)
-    runCommand('mrmath ' + dwi_path + ' mean ' + mean_dwi_path + ' -axis 3')
-    gm_mean  = float(getImageStat(mean_dwi_path, 'mean', 'gm_mask.mif'))
-    csf_mean = float(getImageStat(mean_dwi_path, 'mean', 'csf_mask.mif'))
-    gm_responses .append( str(gm_mean  * math.sqrt(4.0 * math.pi)) )
-    csf_responses.append( str(csf_mean * math.sqrt(4.0 * math.pi)) )
     this_b_lmax_option = ''
     if wm_lmax:
       this_b_lmax_option = ' -lmax ' + str(wm_lmax[index])
     runCommand('amp2sh ' + dwi_path + ' - | sh2response - wm_sf_mask.mif dirs.mif wm_response_b' + str(b) + '.txt' + this_b_lmax_option)
     wm_response = open('wm_response_b' + str(b) + '.txt', 'r').read().split()
     wm_responses.append(wm_response)
     max_length = max(max_length, len(wm_response))
+    mean_dwi_path = 'dwi_b' + str(b) + '_mean.mif'
+    runCommand('mrmath ' + dwi_path + ' mean ' + mean_dwi_path + ' -axis 3')
+    gm_mean  = float(getImageStat(mean_dwi_path, 'mean', 'gm_mask.mif'))
+    csf_mean = float(getImageStat(mean_dwi_path, 'mean', 'csf_mask.mif'))
+    gm_responses .append( str(gm_mean  * math.sqrt(4.0 * math.pi)) )
+    csf_responses.append( str(csf_mean * math.sqrt(4.0 * math.pi)) )
 
-  with open('gm.txt', 'w') as f:
-    for line in gm_responses:
-      f.write(line + '\n')
   with open('wm.txt', 'w') as f:
     for line in wm_responses:
       line += ['0'] * (max_length - len(line))
       f.write(' '.join(line) + '\n')
+  with open('gm.txt', 'w') as f:
+    for line in gm_responses:
+      f.write(line + '\n')
   with open('csf.txt', 'w') as f:
     for line in csf_responses:
       f.write(line + '\n')
 
-  shutil.copyfile('gm.txt',  os.path.join(lib.app.workingDir, lib.app.args.out_gm))
   shutil.copyfile('wm.txt',  os.path.join(lib.app.workingDir, lib.app.args.out_wm))
+  shutil.copyfile('gm.txt',  os.path.join(lib.app.workingDir, lib.app.args.out_gm))
   shutil.copyfile('csf.txt', os.path.join(lib.app.workingDir, lib.app.args.out_csf))
 
-  # Generate output 4D binary image with voxel selections
-  runCommand('mrcat gm_mask.mif wm_sf_mask.mif csf_mask.mif voxels.mif -axis 3')
+  # Generate output 4D binary image with voxel selections; RGB as in MSMT-CSD paper
+  runCommand('mrcat csf_mask.mif gm_mask.mif wm_sf_mask.mif voxels.mif -axis 3')