Merge pull request #872 from gagnonanthony/perpdiff

Add perpendicular diffusivity in NODDI priors

scripts/scil_NODDI_priors.py

@@ -2,7 +2,8 @@
 # -*- coding: utf-8 -*-
 
 """
-Compute the axial (para_diff) and mean (iso_diff) diffusivity priors for NODDI.
+Compute the axial (para_diff), radial (perp_diff),
+and mean (iso_diff) diffusivity priors for NODDI.
 
 Formerly: scil_compute_NODDI_priors.py
 """
@@ -37,6 +38,8 @@ def _build_arg_parser():
                    help='Path to the FA volume.')
     p.add_argument('in_AD',
                    help='Path to the axial diffusivity (AD) volume.')
+    p.add_argument('in_RD',
+                   help='Path to the radial diffusivity (RD) volume.')
     p.add_argument('in_MD',
                    help='Path to the mean diffusivity (MD) volume.')
 
@@ -66,10 +69,14 @@ def _build_arg_parser():
              'priors. [center of the 3D volume]')
 
     g3 = p.add_argument_group('Outputs')
-    g3.add_argument('--out_txt_1fiber', metavar='FILE',
+    g3.add_argument('--out_txt_1fiber_para', metavar='FILE',
                     help='Output path for the text file containing the single '
                          'fiber average value of AD.\nIf not set, the file '
                          'will not be saved.')
+    g3.add_argument('--out_txt_1fiber_perp', metavar='FILE',
+                    help='Output path for the text file containing the single '
+                         'fiber average value of RD.\nIf not set, the file '
+                         'will not be saved.')
     g3.add_argument('--out_mask_1fiber', metavar='FILE',
                     help='Output path for single fiber mask. If not set, the '
                          'mask will not be saved.')
@@ -91,14 +98,16 @@ def main():
     parser = _build_arg_parser()
     args = parser.parse_args()
 
-    assert_inputs_exist(parser, [args.in_AD, args.in_FA, args.in_MD])
+    assert_inputs_exist(parser, [args.in_AD, args.in_FA, args.in_MD,
+                                 args.in_RD])
     assert_outputs_exist(parser, args, [],
                          [args.out_mask_1fiber,
                           args.out_mask_ventricles,
                           args.out_txt_ventricles,
-                          args.out_txt_1fiber])
+                          args.out_txt_1fiber_para,
+                          args.out_txt_1fiber_perp])
 
-    assert_same_resolution([args.in_AD, args.in_FA, args.in_MD])
+    assert_same_resolution([args.in_AD, args.in_FA, args.in_MD, args.in_RD])
 
     log_level = logging.DEBUG if args.verbose else logging.INFO
     logging.getLogger().setLevel(log_level)
@@ -109,6 +118,7 @@ def main():
 
     md_data = nib.load(args.in_MD).get_fdata(dtype=np.float32)
     ad_data = nib.load(args.in_AD).get_fdata(dtype=np.float32)
+    rd_data = nib.load(args.in_RD).get_fdata(dtype=np.float32)
 
     mask_cc = np.zeros(fa_data.shape, dtype=np.uint8)
     mask_vent = np.zeros(fa_data.shape, dtype=np.uint8)
@@ -129,6 +139,9 @@ def main():
     roi_ad = ad_data[max(int(ci - w), 0): min(int(ci + w), fa_shape[0]),
                      max(int(cj - w), 0): min(int(cj + w), fa_shape[1]),
                      max(int(ck - w), 0): min(int(ck + w), fa_shape[2])]
+    roi_rd = rd_data[max(int(ci - w), 0): min(int(ci + w), fa_shape[0]),
+                     max(int(cj - w), 0): min(int(cj + w), fa_shape[1]),
+                     max(int(ck - w), 0): min(int(ck + w), fa_shape[2])]
     roi_md = md_data[max(int(ci - w), 0): min(int(ci + w), fa_shape[0]),
                      max(int(cj - w), 0): min(int(cj + w), fa_shape[1]),
                      max(int(ck - w), 0): min(int(ck + w), fa_shape[2])]
@@ -144,8 +157,11 @@ def main():
 
     logging.debug('Number of voxels found in single fiber area: {}'.format(N))
 
-    cc_avg = np.mean(roi_ad[indices])
-    cc_std = np.std(roi_ad[indices])
+    cc_avg_para = np.mean(roi_ad[indices])
+    cc_std_para = np.std(roi_ad[indices])
+
+    cc_avg_perp = np.mean(roi_rd[indices])
+    cc_std_perp = np.std(roi_rd[indices])
 
     indices[0][:] += ci - w
     indices[1][:] += cj - w
@@ -171,14 +187,21 @@ def main():
     if args.out_mask_ventricles:
         nib.save(nib.Nifti1Image(mask_vent, affine), args.out_mask_ventricles)
 
-    if args.out_txt_1fiber:
-        np.savetxt(args.out_txt_1fiber, [cc_avg], fmt='%f')
+    if args.out_txt_1fiber_para:
+        np.savetxt(args.out_txt_1fiber_para, [cc_avg_para], fmt='%f')
+
+    if args.out_txt_1fiber_perp:
+        np.savetxt(args.out_txt_1fiber_perp, [cc_avg_perp], fmt='%f')
 
     if args.out_txt_ventricles:
         np.savetxt(args.out_txt_ventricles, [vent_avg], fmt='%f')
 
-    logging.info("Average AD in single fiber areas: {} +- {}".format(cc_avg,
-                                                                     cc_std))
+    logging.info("Average AD in single fiber areas: {} +- {}".format(
+                                                                cc_avg_para,
+                                                                cc_std_para))
+    logging.info("Average RD in single fiber areas: {} +- {}".format(
+                                                                cc_avg_perp,
+                                                                cc_std_perp))
     logging.info("Average MD in ventricles: {} +- {}".format(vent_avg,
                                                              vent_std))
 

scripts/tests/test_NODDI_priors.py

@@ -7,7 +7,7 @@
 from scilpy.io.fetcher import fetch_data, get_home, get_testing_files_dict
 
 # If they already exist, this only takes 5 seconds (check md5sum)
-fetch_data(get_testing_files_dict(), keys=['commit_amico.zip'])
+fetch_data(get_testing_files_dict(), keys=['processing.zip'])
 tmp_dir = tempfile.TemporaryDirectory()
 
 
@@ -18,14 +18,17 @@ def test_help_option(script_runner):
 
 def test_execution_commit_amico(script_runner):
     os.chdir(os.path.expanduser(tmp_dir.name))
-    in_fa = os.path.join(get_home(), 'commit_amico',
+    in_fa = os.path.join(get_home(), 'processing',
                          'fa.nii.gz')
-    in_ad = os.path.join(get_home(), 'commit_amico',
+    in_ad = os.path.join(get_home(), 'processing',
                          'ad.nii.gz')
-    in_md = os.path.join(get_home(), 'commit_amico',
+    in_md = os.path.join(get_home(), 'processing',
                          'md.nii.gz')
-    ret = script_runner.run('scil_NODDI_priors.py', in_fa, in_ad, in_md,
-                            '--out_txt_1fiber', '1fiber.txt',
+    in_rd = os.path.join(get_home(), 'processing',
+                         'rd.nii.gz')
+    ret = script_runner.run('scil_NODDI_priors.py', in_fa, in_ad, in_rd, in_md,
+                            '--out_txt_1fiber_para', '1fiber_para.txt',
+                            '--out_txt_1fiber_perp', '1fiber_perp.txt',
                             '--out_mask_1fiber', '1fiber.nii.gz',
                             '--out_txt_ventricles', 'ventricules.txt',
                             '--out_mask_ventricles', 'ventricules.nii.gz')