From b89caaee3b361e38d393ba9f8cd8f3d7d35938f1 Mon Sep 17 00:00:00 2001 From: gagnonanthony <79757265+gagnonanthony@users.noreply.github.com> Date: Wed, 6 Dec 2023 02:35:59 +0000 Subject: [PATCH 1/2] add perp_diff to priors --- scripts/scil_compute_NODDI_priors.py | 53 ++++++++++++++++++++-------- 1 file changed, 38 insertions(+), 15 deletions(-) diff --git a/scripts/scil_compute_NODDI_priors.py b/scripts/scil_compute_NODDI_priors.py index 02f455376..1a119c175 100755 --- a/scripts/scil_compute_NODDI_priors.py +++ b/scripts/scil_compute_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. """ import argparse @@ -35,6 +36,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.') @@ -64,17 +67,21 @@ 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.') + '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.') g3.add_argument('--out_txt_ventricles', metavar='FILE', help='Output path for the text file containing the ' - 'ventricles average value of MD.\nIf not set, the file ' - 'will not be saved.') + 'ventricles average value of MD.\nIf not set, the ' + 'file will not be saved.') g3.add_argument('--out_mask_ventricles', metavar='FILE', help='Output path for the ventricule mask.\nIf not set, ' 'the mask will not be saved.') @@ -89,14 +96,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) @@ -107,6 +116,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) @@ -127,6 +137,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])] @@ -142,8 +155,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 @@ -169,14 +185,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)) From a5f67a7e7837554ea4216857039c6f1e6ca84020 Mon Sep 17 00:00:00 2001 From: gagnonanthony <79757265+gagnonanthony@users.noreply.github.com> Date: Tue, 19 Dec 2023 15:42:06 -0500 Subject: [PATCH 2/2] update test + pep8 --- scripts/scil_NODDI_priors.py | 6 +++--- scripts/tests/test_NODDI_priors.py | 15 +++++++++------ 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/scripts/scil_NODDI_priors.py b/scripts/scil_NODDI_priors.py index d21f28b34..297c730bb 100755 --- a/scripts/scil_NODDI_priors.py +++ b/scripts/scil_NODDI_priors.py @@ -2,7 +2,7 @@ # -*- coding: utf-8 -*- """ -Compute the axial (para_diff), radial (perp_diff), +Compute the axial (para_diff), radial (perp_diff), and mean (iso_diff) diffusivity priors for NODDI. """ @@ -157,7 +157,7 @@ def main(): 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]) @@ -187,7 +187,7 @@ def main(): 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') diff --git a/scripts/tests/test_NODDI_priors.py b/scripts/tests/test_NODDI_priors.py index 7ce35481f..577ee2f9e 100644 --- a/scripts/tests/test_NODDI_priors.py +++ b/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')