From ef4f075dc377027c5c31919920b354a876e489da Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Wed, 6 Mar 2024 10:51:05 -0500 Subject: [PATCH 01/12] adjust_b1_header: Ok! --- scilpy/reconst/mti.py | 6 ++++-- scripts/scil_mti_adjust_B1_header.py | 6 +++--- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/scilpy/reconst/mti.py b/scilpy/reconst/mti.py index 066457e2e..976c53ad7 100644 --- a/scilpy/reconst/mti.py +++ b/scilpy/reconst/mti.py @@ -363,8 +363,10 @@ def adjust_B1_map_header(B1_img, slope): Parameters ---------- - B1_img: B1 nifti image object. - slope: Slope value, obtained from the image header. + B1_img: nifti image object + The B1 map. + slope: float + The slope value, obtained from the image header. Returns ---------- diff --git a/scripts/scil_mti_adjust_B1_header.py b/scripts/scil_mti_adjust_B1_header.py index 85812a40c..2553946ae 100644 --- a/scripts/scil_mti_adjust_B1_header.py +++ b/scripts/scil_mti_adjust_B1_header.py @@ -1,8 +1,8 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ -Correct B1 map header problem. - +Correct B1 map header problem, by applying the scaling (slope) and setting +the slope to 1. """ import argparse @@ -40,7 +40,7 @@ def main(): logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_outputs_exist(parser, args, args.out_B1_map) - assert_inputs_exist(parser, (args.in_B1_map, args.in_B1_json)) + assert_inputs_exist(parser, [args.in_B1_map, args.in_B1_json]) with open(args.in_B1_json) as curr_json: b1_json = json.load(curr_json) From f495ab7490213d237c93f036f2c9c9b2aea2d52c Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Wed, 6 Mar 2024 10:59:32 -0500 Subject: [PATCH 02/12] qball_metrics: ok! --- scripts/scil_qball_metrics.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/scripts/scil_qball_metrics.py b/scripts/scil_qball_metrics.py index 5c9b8776e..eba9aab20 100755 --- a/scripts/scil_qball_metrics.py +++ b/scripts/scil_qball_metrics.py @@ -8,8 +8,8 @@ By default, will output all possible files, using default names. Specific names can be specified using the file flags specified in the "File flags" section. -If --not_all is set, only the files specified explicitly by the flags -will be output. +If --not_all is set, only the files specified explicitly by the flags will be +output. See [Descoteaux et al MRM 2007, Aganj et al MRM 2009] for details and [Cote et al MEDIA 2013] for quantitative comparisons. @@ -69,21 +69,21 @@ def _build_arg_parser(): 'following flags.') g = p.add_argument_group(title='File flags') - g.add_argument('--gfa', default='', + g.add_argument('--gfa', help='Output filename for the generalized fractional ' 'anisotropy [gfa.nii.gz].') - g.add_argument('--peaks', default='', + g.add_argument('--peaks', help='Output filename for the extracted peaks ' '[peaks.nii.gz].') - g.add_argument('--peak_indices', default='', + g.add_argument('--peak_indices', help='Output filename for the generated peaks ' 'indices on the sphere [peaks_indices.nii.gz].') - g.add_argument('--sh', default='', + g.add_argument('--sh', help='Output filename for the spherical harmonics ' 'coefficients [sh.nii.gz].') - g.add_argument('--nufo', default='', + g.add_argument('--nufo', help='Output filename for the NUFO map [nufo.nii.gz].') - g.add_argument('--a_power', default='', + g.add_argument('--a_power', help='Output filename for the anisotropic power map' '[anisotropic_power.nii.gz].') @@ -112,14 +112,13 @@ def main(): arglist = [args.gfa, args.peaks, args.peak_indices, args.sh, args.nufo, args.a_power] if args.not_all and not any(arglist): - parser.error('When using --not_all, you need to specify at least ' + - 'one file to output.') + parser.error('When using --not_all, you need to specify at least one ' + 'file to output.') assert_inputs_exist(parser, [args.in_dwi, args.in_bval, args.in_bvec]) - assert_outputs_exist(parser, args, arglist) - validate_nbr_processes(parser, args) + assert_outputs_exist(parser, args, [], optional=arglist) - nbr_processes = args.nbr_processes + nbr_processes = validate_nbr_processes(parser, args) parallel = nbr_processes > 1 # Load data @@ -129,7 +128,8 @@ def main(): bvals, bvecs = read_bvals_bvecs(args.in_bval, args.in_bvec) if not is_normalized_bvecs(bvecs): - logging.warning('Your b-vectors do not seem normalized...') + logging.warning('Your b-vectors do not seem normalized... Normalizing ' + 'now.') bvecs = normalize_bvecs(bvecs) # Usage of gtab.b0s_mask in dipy's models is not very well documented, but From 0f43091765a8f585cd8584b0eaa5d1393f98805c Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Wed, 6 Mar 2024 11:16:56 -0500 Subject: [PATCH 03/12] sh_convert: ok! (Clarifying sh_basis stuff) --- scilpy/io/utils.py | 84 ++++++++++++++++++++------------------ scilpy/reconst/sh.py | 3 +- scripts/scil_sh_convert.py | 2 +- 3 files changed, 47 insertions(+), 42 deletions(-) diff --git a/scilpy/io/utils.py b/scilpy/io/utils.py index c2c8ee5e3..ce75e03a7 100644 --- a/scilpy/io/utils.py +++ b/scilpy/io/utils.py @@ -180,12 +180,12 @@ def add_processes_arg(parser): parser.add_argument('--processes', dest='nbr_processes', metavar='NBR', type=int, default=1, help='Number of sub-processes to start. \n' - 'Default: [%(default)s]') + 'Default: [%(default)s]') def add_reference_arg(parser, arg_name=None): if arg_name: - parser.add_argument('--'+arg_name+'_ref', + parser.add_argument('--' + arg_name + '_ref', help='Reference anatomy for {} (if tck/vtk/fib/dpy' ') file\n' 'support (.nii or .nii.gz).'.format(arg_name)) @@ -302,8 +302,8 @@ def add_sh_basis_args(parser, mandatory=False, input_output=False): if input_output: nargs = 2 def_val = ['descoteaux07_legacy', 'tournier07'] - input_output_msg = '\nBoth the input and output bases are ' +\ - 'required, in that order.' + input_output_msg = ('\nBoth the input and output bases are ' + 'required, in that order.') else: nargs = 1 def_val = ['descoteaux07_legacy'] @@ -311,25 +311,25 @@ def add_sh_basis_args(parser, mandatory=False, input_output=False): choices = ['descoteaux07', 'tournier07', 'descoteaux07_legacy', 'tournier07_legacy'] - help_msg = 'Spherical harmonics basis used for the SH coefficients. ' +\ - input_output_msg +\ - '\nMust be either \'descoteaux07\', \'tournier07\', \n' +\ - '\'descoteaux07_legacy\' or \'tournier07_legacy\'' +\ - ' [%(default)s]:\n' +\ - ' \'descoteaux07\' : SH basis from the Descoteaux ' +\ - 'et al.\n' +\ - ' MRM 2007 paper\n' +\ - ' \'tournier07\' : SH basis from the new ' +\ - 'Tournier et al.\n' +\ - ' NeuroImage 2019 paper, as in ' +\ - 'MRtrix 3.\n' +\ - ' \'descoteaux07_legacy\': SH basis from the legacy Dipy ' +\ - 'implementation\n' +\ - ' of the ' +\ - 'Descoteaux et al. MRM 2007 paper\n' +\ - ' \'tournier07_legacy\' : SH basis from the legacy ' +\ - 'Tournier et al.\n' +\ - ' NeuroImage 2007 paper.' + help_msg = ("Spherical harmonics basis used for the SH coefficients. " + "{}\n" + "Must be either descoteaux07', 'tournier07', \n" + "'descoteaux07_legacy' or 'tournier07_legacy' [%(default)s]:\n" + " 'descoteaux07' : SH basis from the Descoteaux " + "et al.\n" + " MRM 2007 paper\n" + " 'tournier07' : SH basis from the new " + "Tournier et al.\n" + " NeuroImage 2019 paper, as in " + "MRtrix 3.\n" + " 'descoteaux07_legacy': SH basis from the legacy Dipy " + "implementation\n" + " of the Descoteaux et al. MRM 2007 " + "paper\n" + " 'tournier07_legacy' : SH basis from the legacy " + "Tournier et al.\n" + " NeuroImage 2007 paper." + .format(input_output_msg)) if mandatory: arg_name = 'sh_basis' @@ -353,10 +353,14 @@ def parse_sh_basis_arg(args): Returns ------- - sh_basis : string - Spherical harmonic basis name. - is_legacy : bool - Whether or not the SH basis is in its legacy form. + if args.sh_basis is a list of one string: + sh_basis : string + Spherical harmonic basis name. + is_legacy : bool + Whether the SH basis is in its legacy form. + else: (args:sh_basis is a list of two strings) + Returns a Tuple of 4 values: + (sh_basis_in, is_legacy_in, sh_basis_out, is_legacy_out) """ sh_basis_name = args.sh_basis[0] sh_basis = 'descoteaux07' if 'descoteaux07' in sh_basis_name \ @@ -373,7 +377,7 @@ def parse_sh_basis_arg(args): def add_nifti_screenshot_default_args( - parser, slice_ids_mandatory=True, transparency_mask_mandatory=True + parser, slice_ids_mandatory=True, transparency_mask_mandatory=True ): _mask_prefix = "" if transparency_mask_mandatory else "--" @@ -385,8 +389,8 @@ def add_nifti_screenshot_default_args( "the transparency mask are selected." _output_help = "Name of the output image(s). If multiple slices are " \ "provided (or none), their index will be append to " \ - "the name (e.g. volume.jpg, volume.png becomes " \ - "volume_slice_0.jpg, volume_slice_0.png)." + "the name (e.g. volume.jpg, volume.png becomes " \ + "volume_slice_0.jpg, volume_slice_0.png)." # Positional arguments parser.add_argument( @@ -422,12 +426,12 @@ def add_nifti_screenshot_default_args( def add_nifti_screenshot_overlays_args( - parser, labelmap_overlay=True, mask_overlay=True, - transparency_is_overlay=False + parser, labelmap_overlay=True, mask_overlay=True, + transparency_is_overlay=False ): if labelmap_overlay: parser.add_argument( - "--in_labelmap", help="Labelmap 3D Nifti image (.nii/.nii.gz).") + "--in_labelmap", help="Labelmap 3D Nifti image (.nii/.nii.gz).") parser.add_argument( "--labelmap_cmap_name", default="viridis", help="Colormap name for the labelmap image data. [%(default)s]") @@ -540,6 +544,7 @@ def assert_inputs_exist(parser, required, optional=None): optional: string or list of paths Optional paths to be checked. """ + def check(path): if not os.path.isfile(path): parser.error('Input file {} does not exist'.format(path)) @@ -576,6 +581,7 @@ def assert_outputs_exist(parser, args, required, optional=None, check_dir_exists: bool Test if output directory exists. """ + def check(path): if os.path.isfile(path) and not args.overwrite: parser.error('Output file {} exists. Use -f to force ' @@ -620,6 +626,7 @@ def assert_output_dirs_exist_and_empty(parser, args, required, create_dir: bool If true, create the directory if it does not exist. """ + def check(path): if not os.path.isdir(path): if not create_dir: @@ -751,18 +758,18 @@ def read_info_from_mb_bdo(filename): geometry = root.attrib['type'] center_tag = root.find('origin') flip = [-1, -1, 1] - center = [flip[0]*float(center_tag.attrib['x'].replace(',', '.')), - flip[1]*float(center_tag.attrib['y'].replace(',', '.')), - flip[2]*float(center_tag.attrib['z'].replace(',', '.'))] + center = [flip[0] * float(center_tag.attrib['x'].replace(',', '.')), + flip[1] * float(center_tag.attrib['y'].replace(',', '.')), + flip[2] * float(center_tag.attrib['z'].replace(',', '.'))] row_list = tree.iter('Row') radius = [None, None, None] for i, row in enumerate(row_list): for j in range(0, 3): if j == i: - key = 'col' + str(j+1) + key = 'col' + str(j + 1) radius[i] = float(row.attrib[key].replace(',', '.')) else: - key = 'col' + str(j+1) + key = 'col' + str(j + 1) value = float(row.attrib[key].replace(',', '.')) if abs(value) > 0.01: raise ValueError('Does not support rotation, for now \n' @@ -912,7 +919,6 @@ def range_checker(arg: str): def get_default_screenshotting_data(args): - volume_img = nib.load(args.in_volume) transparency_mask_img = None diff --git a/scilpy/reconst/sh.py b/scilpy/reconst/sh.py index b0cb5b72f..b86d99119 100644 --- a/scilpy/reconst/sh.py +++ b/scilpy/reconst/sh.py @@ -577,8 +577,7 @@ def convert_sh_basis(shm_coeff, sphere, mask=None, mask = np.sum(shm_coeff, axis=3).astype(bool) nbr_processes = multiprocessing.cpu_count() \ - if nbr_processes is None or nbr_processes < 0 \ - else nbr_processes + if nbr_processes is None or nbr_processes < 0 else nbr_processes # Ravel the first 3 dimensions while keeping the 4th intact, like a list of # 1D time series voxels. Then separate it in chunks of len(nbr_processes). diff --git a/scripts/scil_sh_convert.py b/scripts/scil_sh_convert.py index 5cd192521..1251221c1 100755 --- a/scripts/scil_sh_convert.py +++ b/scripts/scil_sh_convert.py @@ -2,7 +2,7 @@ # -*- coding: utf-8 -*- """ -Convert a SH file between the two of the following bases choices: +Convert a SH file between the two of the following base choices: 'descoteaux07', 'descoteaux07_legacy', 'tournier07' or 'tournier07_legacy'. Using the sh_basis argument, both the input and the output SH bases must be given, in the order. For more information about the bases, see From c7a3039594f3673de5877ccae69c66f81cdac566 Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Wed, 6 Mar 2024 11:19:24 -0500 Subject: [PATCH 04/12] sh_fusion: ok! --- scripts/scil_sh_fusion.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/scripts/scil_sh_fusion.py b/scripts/scil_sh_fusion.py index 748835754..ced8acdd8 100755 --- a/scripts/scil_sh_fusion.py +++ b/scripts/scil_sh_fusion.py @@ -3,8 +3,8 @@ """ Merge a list of Spherical Harmonics files. -This merges the coefficients of multiple Spherical Harmonics files -by taking, for each coefficient, the one with the largest magnitude. +This merges the coefficients of multiple Spherical Harmonics files by taking, +for each coefficient, the one with the largest magnitude. Can be used to merge fODFs computed from different shells into 1, while conserving the most relevant information. @@ -67,8 +67,7 @@ def main(): out_coeffs = first_im.get_fdata(dtype=np.float32) for sh_file in args.in_shs[1:]: - im = nib.load(sh_file) - im_dat = im.get_fdata(dtype=np.float32) + im_dat = nib.load(sh_file).get_fdata(dtype=np.float32) out_coeffs = np.where(np.abs(im_dat) > np.abs(out_coeffs), im_dat, out_coeffs) From bfe5ee1d66e4582a8d5b367c379f2a4b4705bc2e Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Wed, 6 Mar 2024 11:26:15 -0500 Subject: [PATCH 05/12] Sh_to_rish: ok! --- scripts/scil_sh_to_rish.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/scripts/scil_sh_to_rish.py b/scripts/scil_sh_to_rish.py index efccc079a..59bb20bf4 100755 --- a/scripts/scil_sh_to_rish.py +++ b/scripts/scil_sh_to_rish.py @@ -2,12 +2,11 @@ # -*- coding: utf-8 -*- """ -Compute the RISH (Rotationally Invariant Spherical Harmonics) features -of an SH signal [1]. +Compute the RISH (Rotationally Invariant Spherical Harmonics) features of an SH +signal [1]. -Each RISH feature map is the total energy of its -associated order. Mathematically, it is the sum of the squared SH -coefficients of the SH order. +Each RISH feature map is the total energy of its associated order. +Mathematically, it is the sum of the squared SH coefficients of the SH order. This script supports both symmetrical and asymmetrical SH images as input, of any SH order. @@ -37,9 +36,12 @@ def _build_arg_parser(): p = argparse.ArgumentParser( description=__doc__, formatter_class=argparse.RawTextHelpFormatter) p.add_argument('in_sh', - help='Path of the sh image. Must be a symmetric SH file.') + help='Path of the sh image. They can be formatted in any ' + 'sh basis, but we \nexpect it to be a symmetrical ' + 'one. Else, provide --full_basis.') p.add_argument('out_prefix', - help='Prefix of the output RISH files to save.') + help='Prefix of the output RISH files to save. Suffixes ' + 'will be \nbased on the sh orders.') p.add_argument('--full_basis', action="store_true", help="Input SH image uses a full SH basis (asymmetrical).") p.add_argument('--mask', @@ -81,6 +83,7 @@ def main(): # Save each RISH feature as a separate file for i, fname in enumerate(output_fnames): + logging.info("Saving {}".format(fname)) nib.save(nib.Nifti1Image(rish[..., i], sh_img.affine), fname) From 213ef028391c40884c2f35808ebc0a29825a501e Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Wed, 6 Mar 2024 13:36:25 -0500 Subject: [PATCH 06/12] sh_to_sf: still a little big but ok! Added option in test. --- scilpy/reconst/sh.py | 2 +- scripts/scil_sh_to_sf.py | 2 +- scripts/tests/test_sh_to_sf.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/scilpy/reconst/sh.py b/scilpy/reconst/sh.py index b86d99119..f6ff22bcc 100644 --- a/scilpy/reconst/sh.py +++ b/scilpy/reconst/sh.py @@ -646,7 +646,7 @@ def convert_sh_to_sf(shm_coeff, sphere, mask=None, dtype="float32", If True, use a full SH basis (even and odd orders) for the input SH coefficients. is_input_legacy : bool, optional - Whether or not the input basis is in its legacy form. + Whether the input basis is in its legacy form. nbr_processes: int, optional The number of subprocesses to use. Default: multiprocessing.cpu_count() diff --git a/scripts/scil_sh_to_sf.py b/scripts/scil_sh_to_sf.py index 02f11b979..875cba13c 100755 --- a/scripts/scil_sh_to_sf.py +++ b/scripts/scil_sh_to_sf.py @@ -151,7 +151,7 @@ def main(): # Sample SF from SH if args.sphere: sphere = get_sphere(args.sphere) - elif args.in_bvec: + else: # args.in_bvec is set. gtab = gradient_table(bvals, bvecs, b0_threshold=args.b0_threshold) # Remove bvecs corresponding to b0 images bvecs = bvecs[np.logical_not(gtab.b0s_mask)] diff --git a/scripts/tests/test_sh_to_sf.py b/scripts/tests/test_sh_to_sf.py index 6473bc85a..879d9ec91 100755 --- a/scripts/tests/test_sh_to_sf.py +++ b/scripts/tests/test_sh_to_sf.py @@ -61,7 +61,7 @@ def test_execution_no_bval(script_runner): # --sphere but no --bval ret = script_runner.run('scil_sh_to_sf.py', in_sh, 'sf_724.nii.gz', '--in_b0', in_b0, - '--out_bvec', 'sf_724.bvec', + '--out_bvec', 'sf_724.bvec', '--b0_scaling', '--sphere', 'symmetric724', '--dtype', 'float32', '-f') assert ret.success From c54be779f2862e23c0d17cd7d7480c4c1f45f740 Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Wed, 6 Mar 2024 14:31:01 -0500 Subject: [PATCH 07/12] NODDI_priors: clarifying --- scilpy/reconst/noddi.py | 0 scripts/scil_NODDI_priors.py | 120 +++++++++++++++++------------------ 2 files changed, 58 insertions(+), 62 deletions(-) create mode 100644 scilpy/reconst/noddi.py diff --git a/scilpy/reconst/noddi.py b/scilpy/reconst/noddi.py new file mode 100644 index 000000000..e69de29bb diff --git a/scripts/scil_NODDI_priors.py b/scripts/scil_NODDI_priors.py index 9d74b3393..9847adec0 100755 --- a/scripts/scil_NODDI_priors.py +++ b/scripts/scil_NODDI_priors.py @@ -20,20 +20,18 @@ add_overwrite_arg, add_verbose_arg) - EPILOG = """ Reference: [1] Zhang H, Schneider T, Wheeler-Kingshott CA, Alexander DC. - NODDI: practical in vivo neurite orientation dispersion - and density imaging of the human brain. - NeuroImage. 2012 Jul 16;61:1000-16. + NODDI: practical in vivo neurite orientation dispersion and density + imaging of the human brain. NeuroImage. 2012 Jul 16;61:1000-16. """ def _build_arg_parser(): p = argparse.ArgumentParser( description=__doc__, epilog=EPILOG, - formatter_class=argparse.RawDescriptionHelpFormatter) + formatter_class=argparse.RawTextHelpFormatter) p.add_argument('in_FA', help='Path to the FA volume.') p.add_argument('in_AD', @@ -45,28 +43,29 @@ def _build_arg_parser(): g1 = p.add_argument_group('Metrics options') g1.add_argument( - '--fa_min', type=float, default='0.7', + '--fa_single_fiber', type=float, default=0.7, help='Minimal threshold of FA (voxels above that threshold are ' - 'considered in the single fiber mask). [%(default)s]') + 'considered in \nthe single fiber mask). [%(default)s]') g1.add_argument( - '--fa_max', type=float, default='0.1', + '--fa_ventricles', type=float, default=0.1, help='Maximal threshold of FA (voxels under that threshold are ' - 'considered in the ventricles). [%(default)s]') + 'considered in \nthe ventricles). [%(default)s]') g1.add_argument( - '--md_min', dest='md_min', type=float, default='0.003', + '--md_ventricles', type=float, default=0.003, help='Minimal threshold of MD in mm2/s (voxels above that threshold ' - 'are considered for in the ventricles). [%(default)s]') + 'are considered \nfor in the ventricles). [%(default)s]') g2 = p.add_argument_group('Regions options') g2.add_argument( '--roi_radius', type=int, default=20, help='Radius of the region used to estimate the priors. The roi will ' - 'be a cube spanning from ROI_CENTER in each direction. ' + 'be a cube spanning \nfrom ROI_CENTER in each direction. ' '[%(default)s]') g2.add_argument( - '--roi_center', metavar='tuple(3)', nargs="+", type=int, + '--roi_center', metavar='pos', nargs=3, type=int, help='Center of the roi of size roi_radius used to estimate the ' - 'priors. [center of the 3D volume]') + 'priors; a 3-value coordinate. \nIf not set, uses the center of ' + 'the 3D volume.') g3 = p.add_argument_group('Outputs') g3.add_argument('--out_txt_1fiber_para', metavar='FILE', @@ -99,6 +98,7 @@ def main(): args = parser.parse_args() logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + # Verifications assert_inputs_exist(parser, [args.in_AD, args.in_FA, args.in_MD, args.in_RD]) assert_outputs_exist(parser, args, [], @@ -110,6 +110,7 @@ def main(): assert_same_resolution([args.in_AD, args.in_FA, args.in_MD, args.in_RD]) + # Loading fa_img = nib.load(args.in_FA) fa_data = fa_img.get_fdata(dtype=np.float32) affine = fa_img.affine @@ -118,10 +119,7 @@ def main(): 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) - - # center + # Finding ROI center if args.roi_center is None: ci, cj, ck = np.array(fa_data.shape[:3]) // 2 else: @@ -132,76 +130,74 @@ def main(): else: ci, cj, ck = args.roi_center + # Get values in the ROI w = args.roi_radius - fa_shape = fa_data.shape - 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])] - roi_fa = fa_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])] - - logging.info('fa_min, fa_max, md_min: {}, {}, {}'.format( - args.fa_min, args.fa_max, args.md_min)) - - indices = np.where((roi_fa > args.fa_min) & (roi_fa < 0.95)) - N = roi_ad[indices].shape[0] - - logging.info('Number of voxels found in single fiber area: {}'.format(N)) - - 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]) - + roi_posx = slice(max(int(ci - w), 0), min(int(ci + w), fa_data.shape[0])) + roi_posy = slice(max(int(cj - w), 0), min(int(cj + w), fa_data.shape[1])) + roi_posz = slice(max(int(ck - w), 0), min(int(ck + w), fa_data.shape[2])) + roi_ad = ad_data[roi_posx, roi_posy, roi_posz] + roi_rd = rd_data[roi_posx, roi_posy, roi_posz] + roi_md = md_data[roi_posx, roi_posy, roi_posz] + roi_fa = fa_data[roi_posx, roi_posy, roi_posz] + + # Get information in single fiber voxels + # Taking voxels with FA < 0.95 just to avoid using weird broken voxels. + indices = np.where((roi_fa > args.fa_single_fiber) & (roi_fa < 0.95)) + nb_voxels = roi_fa[indices].shape[0] + logging.info('Number of voxels found in single fiber area (FA in range ' + '{}-{}]: {}'.format(args.fa_single_fiber, 0.95, nb_voxels)) + single_fiber_ad_mean = np.mean(roi_ad[indices]) + single_fiber_ad_std = np.std(roi_ad[indices]) + single_fiber_rd_mean = np.mean(roi_rd[indices]) + single_fiber_rd_std = np.std(roi_rd[indices]) + + # Create mask of single fiber in ROI indices[0][:] += ci - w indices[1][:] += cj - w indices[2][:] += ck - w - mask_cc[indices] = 1 - - indices = np.where((roi_md > args.md_min) & (roi_fa < args.fa_max)) - N = roi_md[indices].shape[0] + mask_single_fiber = np.zeros(fa_data.shape, dtype=np.uint8) + mask_single_fiber[indices] = 1 - logging.info('Number of voxels found in ventricles: {}'.format(N)) + # Get information in ventricles + indices = np.where((roi_md > args.md_ventricles) & + (roi_fa < args.fa_ventricles)) + nb_voxels = roi_md[indices].shape[0] + logging.info('Number of voxels found in ventricles (FA < {} and MD > {}): ' + '{}'.format(args.fa_ventricles, args.md_ventricles, nb_voxels)) vent_avg = np.mean(roi_md[indices]) vent_std = np.std(roi_md[indices]) + # Create mask of ventricle in ROI indices[0][:] += ci - w indices[1][:] += cj - w indices[2][:] += ck - w + mask_vent = np.zeros(fa_data.shape, dtype=np.uint8) mask_vent[indices] = 1 + # Saving if args.out_mask_1fiber: - nib.save(nib.Nifti1Image(mask_cc, affine), args.out_mask_1fiber) + nib.save(nib.Nifti1Image(mask_single_fiber, affine), + args.out_mask_1fiber) if args.out_mask_ventricles: nib.save(nib.Nifti1Image(mask_vent, affine), args.out_mask_ventricles) if args.out_txt_1fiber_para: - np.savetxt(args.out_txt_1fiber_para, [cc_avg_para], fmt='%f') + np.savetxt(args.out_txt_1fiber_para, [single_fiber_ad_mean], fmt='%f') if args.out_txt_1fiber_perp: - np.savetxt(args.out_txt_1fiber_perp, [cc_avg_perp], fmt='%f') + np.savetxt(args.out_txt_1fiber_perp, [single_fiber_rd_mean], 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_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)) + logging.info("Average AD in single fiber areas: {} +- {}" + .format(single_fiber_ad_mean, single_fiber_ad_std)) + logging.info("Average RD in single fiber areas: {} +- {}" + .format(single_fiber_rd_mean, single_fiber_rd_std)) + logging.info("Average MD in ventricles: {} +- {}" + .format(vent_avg, vent_std)) if __name__ == "__main__": From 62343d72ab213979321fe98d655e50c02aa6f081 Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Thu, 7 Mar 2024 09:12:24 -0500 Subject: [PATCH 08/12] Noddi maps: ok. Use official tol arg. --- scilpy/reconst/noddi.py | 0 scripts/scil_NODDI_maps.py | 54 ++++++++++++++++++-------------- scripts/tests/test_NODDI_maps.py | 2 +- 3 files changed, 31 insertions(+), 25 deletions(-) delete mode 100644 scilpy/reconst/noddi.py diff --git a/scilpy/reconst/noddi.py b/scilpy/reconst/noddi.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/scripts/scil_NODDI_maps.py b/scripts/scil_NODDI_maps.py index b5eb9853d..abe7cd231 100755 --- a/scripts/scil_NODDI_maps.py +++ b/scripts/scil_NODDI_maps.py @@ -26,8 +26,11 @@ add_verbose_arg, assert_inputs_exist, assert_output_dirs_exist_and_empty, - redirect_stdout_c) -from scilpy.gradients.bvec_bval_tools import identify_shells + redirect_stdout_c, add_tolerance_arg, + add_skip_b0_check_arg) +from scilpy.gradients.bvec_bval_tools import (check_b0_threshold, + identify_shells) + EPILOG = """ Reference: @@ -41,7 +44,7 @@ def _build_arg_parser(): p = argparse.ArgumentParser( description=__doc__, epilog=EPILOG, - formatter_class=argparse.RawDescriptionHelpFormatter) + formatter_class=argparse.RawTextHelpFormatter) p.add_argument('in_dwi', help='DWI file acquired with a NODDI compatible protocol ' @@ -56,10 +59,9 @@ def _build_arg_parser(): p.add_argument('--out_dir', default="results", help='Output directory for the NODDI results. ' '[%(default)s]') - p.add_argument('--b_thr', type=int, default=40, - help='Limit value to consider that a b-value is on an ' - 'existing shell. Above this limit, the b-value is ' - 'placed on a new shell. This includes b0s values.') + add_tolerance_arg(p) + add_skip_b0_check_arg(p, will_overwrite_with_min=False, + b0_tol_name='--tolerance') g1 = p.add_argument_group(title='Model options') g1.add_argument('--para_diff', type=float, default=1.7e-3, @@ -101,37 +103,40 @@ def main(): logging.getLogger().setLevel(logging.getLevelName(args.verbose)) redirected_stdout = redirect_stdout(sys.stdout) + # Verifications if args.compute_only and not args.save_kernels: parser.error('--compute_only must be used with --save_kernels.') assert_inputs_exist(parser, [args.in_dwi, args.in_bval, args.in_bvec], args.mask) - assert_output_dirs_exist_and_empty(parser, args, - args.out_dir, + assert_output_dirs_exist_and_empty(parser, args, args.out_dir, optional=args.save_kernels) - # Generage a scheme file from the bvals and bvecs files + # Generate a scheme file from the bvals and bvecs files + bvals, _ = read_bvals_bvecs(args.in_bval, args.in_bvec) + _ = check_b0_threshold(bvals.min(), b0_thr=args.tolerance, + skip_b0_check=args.skip_b0_check) + shells_centroids, indices_shells = identify_shells(bvals, args.tolerance, + round_centroids=True) + + # Save the resulting bvals to a temporary file tmp_dir = tempfile.TemporaryDirectory() tmp_scheme_filename = os.path.join(tmp_dir.name, 'gradients.b') tmp_bval_filename = os.path.join(tmp_dir.name, 'bval') - bvals, _ = read_bvals_bvecs(args.in_bval, args.in_bvec) - shells_centroids, indices_shells = identify_shells(bvals, - args.b_thr, - round_centroids=True) np.savetxt(tmp_bval_filename, shells_centroids[indices_shells], newline=' ', fmt='%i') fsl2mrtrix(tmp_bval_filename, args.in_bvec, tmp_scheme_filename) - logging.info('Compute NODDI with AMICO on {} shells at found ' - 'at {}.'.format(len(shells_centroids), shells_centroids)) + + logging.info('Will compute NODDI with AMICO on {} shells at found at {}.' + .format(len(shells_centroids), np.sort(shells_centroids))) with redirected_stdout: # Load the data amico.core.setup() ae = amico.Evaluation('.', '.') - ae.load_data(args.in_dwi, - tmp_scheme_filename, - mask_filename=args.mask) + ae.load_data(args.in_dwi, tmp_scheme_filename, mask_filename=args.mask) + # Compute the response functions ae.set_model("NODDI") @@ -139,18 +144,18 @@ def main(): intra_orient_distr = np.hstack((np.array([0.03, 0.06]), np.linspace(0.09, 0.99, 10))) - ae.model.set(args.para_diff, args.iso_diff, - intra_vol_frac, intra_orient_distr, - False) + ae.model.set(dPar=args.para_diff, dIso=args.iso_diff, + IC_VFs=intra_vol_frac, IC_ODs=intra_orient_distr, + isExvivo=False) ae.set_solver(lambda1=args.lambda1, lambda2=args.lambda2) # The kernels are, by default, set to be in the current directory # Depending on the choice, manually change the saving location if args.save_kernels: - kernels_dir = os.path.join(args.save_kernels) + kernels_dir = args.save_kernels regenerate_kernels = True elif args.load_kernels: - kernels_dir = os.path.join(args.load_kernels) + kernels_dir = args.load_kernels regenerate_kernels = False else: kernels_dir = os.path.join(tmp_dir.name, 'kernels', ae.model.id) @@ -166,6 +171,7 @@ def main(): # Model fit ae.fit() + # Save the results ae.save_results() diff --git a/scripts/tests/test_NODDI_maps.py b/scripts/tests/test_NODDI_maps.py index 946656adc..dfbce2763 100644 --- a/scripts/tests/test_NODDI_maps.py +++ b/scripts/tests/test_NODDI_maps.py @@ -28,7 +28,7 @@ def test_execution_commit_amico(script_runner): 'mask.nii.gz') ret = script_runner.run('scil_NODDI_maps.py', in_dwi, in_bval, in_bvec, '--mask', mask, - '--out_dir', 'noddi', '--b_thr', '30', + '--out_dir', 'noddi', '--tol', '30', '--para_diff', '0.0017', '--iso_diff', '0.003', '--lambda1', '0.5', '--lambda2', '0.001', '--processes', '1') From a8d6ec806588348828e638357d47b26adcf8a9c7 Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Thu, 7 Mar 2024 09:15:35 -0500 Subject: [PATCH 09/12] Answering Max comment: NODDI fails on single shell --- scripts/scil_NODDI_maps.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/scripts/scil_NODDI_maps.py b/scripts/scil_NODDI_maps.py index abe7cd231..9c29ea487 100755 --- a/scripts/scil_NODDI_maps.py +++ b/scripts/scil_NODDI_maps.py @@ -120,6 +120,16 @@ def main(): shells_centroids, indices_shells = identify_shells(bvals, args.tolerance, round_centroids=True) + nb_shells = len(shells_centroids) + if nb_shells <= 1: + raise ValueError("Amico's NODDI works with data with more than one " + "shell, but you seem to have single-shell data (we " + "found shells {}). Change tolerance if necessary." + .format(np.sort(shells_centroids))) + + logging.info('Will compute NODDI with AMICO on {} shells at found at {}.' + .format(len(shells_centroids), np.sort(shells_centroids))) + # Save the resulting bvals to a temporary file tmp_dir = tempfile.TemporaryDirectory() tmp_scheme_filename = os.path.join(tmp_dir.name, 'gradients.b') @@ -128,9 +138,6 @@ def main(): newline=' ', fmt='%i') fsl2mrtrix(tmp_bval_filename, args.in_bvec, tmp_scheme_filename) - logging.info('Will compute NODDI with AMICO on {} shells at found at {}.' - .format(len(shells_centroids), np.sort(shells_centroids))) - with redirected_stdout: # Load the data amico.core.setup() From 1d04bbcb37b767474f4d94ec972fdfe77494df81 Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Thu, 7 Mar 2024 09:21:20 -0500 Subject: [PATCH 10/12] Fix pep8 --- scripts/scil_NODDI_maps.py | 2 +- scripts/scil_NODDI_priors.py | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/scripts/scil_NODDI_maps.py b/scripts/scil_NODDI_maps.py index 9c29ea487..bb94bb58f 100755 --- a/scripts/scil_NODDI_maps.py +++ b/scripts/scil_NODDI_maps.py @@ -126,7 +126,7 @@ def main(): "shell, but you seem to have single-shell data (we " "found shells {}). Change tolerance if necessary." .format(np.sort(shells_centroids))) - + logging.info('Will compute NODDI with AMICO on {} shells at found at {}.' .format(len(shells_centroids), np.sort(shells_centroids))) diff --git a/scripts/scil_NODDI_priors.py b/scripts/scil_NODDI_priors.py index 9847adec0..14591cd6e 100755 --- a/scripts/scil_NODDI_priors.py +++ b/scripts/scil_NODDI_priors.py @@ -23,7 +23,7 @@ EPILOG = """ Reference: [1] Zhang H, Schneider T, Wheeler-Kingshott CA, Alexander DC. - NODDI: practical in vivo neurite orientation dispersion and density + NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain. NeuroImage. 2012 Jul 16;61:1000-16. """ @@ -163,7 +163,8 @@ def main(): (roi_fa < args.fa_ventricles)) nb_voxels = roi_md[indices].shape[0] logging.info('Number of voxels found in ventricles (FA < {} and MD > {}): ' - '{}'.format(args.fa_ventricles, args.md_ventricles, nb_voxels)) + '{}' + .format(args.fa_ventricles, args.md_ventricles, nb_voxels)) vent_avg = np.mean(roi_md[indices]) vent_std = np.std(roi_md[indices]) From 8fb1182f3ef76c621176168924ac41a668baf01d Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Thu, 7 Mar 2024 13:21:27 -0500 Subject: [PATCH 11/12] Answer Phil comments --- scripts/scil_NODDI_priors.py | 17 +++++++++-------- scripts/scil_sh_convert.py | 2 +- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/scripts/scil_NODDI_priors.py b/scripts/scil_NODDI_priors.py index 14591cd6e..7158d643c 100755 --- a/scripts/scil_NODDI_priors.py +++ b/scripts/scil_NODDI_priors.py @@ -43,15 +43,15 @@ def _build_arg_parser(): g1 = p.add_argument_group('Metrics options') g1.add_argument( - '--fa_single_fiber', type=float, default=0.7, + '--fa_min_single_fiber', type=float, default=0.7, help='Minimal threshold of FA (voxels above that threshold are ' 'considered in \nthe single fiber mask). [%(default)s]') g1.add_argument( - '--fa_ventricles', type=float, default=0.1, + '--fa_max_ventricles', type=float, default=0.1, help='Maximal threshold of FA (voxels under that threshold are ' 'considered in \nthe ventricles). [%(default)s]') g1.add_argument( - '--md_ventricles', type=float, default=0.003, + '--md_min_ventricles', type=float, default=0.003, help='Minimal threshold of MD in mm2/s (voxels above that threshold ' 'are considered \nfor in the ventricles). [%(default)s]') @@ -121,6 +121,7 @@ def main(): # Finding ROI center if args.roi_center is None: + # Using the center of the image: single-fiber region should be the CC ci, cj, ck = np.array(fa_data.shape[:3]) // 2 else: if len(args.roi_center) != 3: @@ -142,10 +143,10 @@ def main(): # Get information in single fiber voxels # Taking voxels with FA < 0.95 just to avoid using weird broken voxels. - indices = np.where((roi_fa > args.fa_single_fiber) & (roi_fa < 0.95)) + indices = np.where((roi_fa > args.fa_min_single_fiber) & (roi_fa < 0.95)) nb_voxels = roi_fa[indices].shape[0] logging.info('Number of voxels found in single fiber area (FA in range ' - '{}-{}]: {}'.format(args.fa_single_fiber, 0.95, nb_voxels)) + '{}-{}]: {}'.format(args.fa_min_single_fiber, 0.95, nb_voxels)) single_fiber_ad_mean = np.mean(roi_ad[indices]) single_fiber_ad_std = np.std(roi_ad[indices]) single_fiber_rd_mean = np.mean(roi_rd[indices]) @@ -159,12 +160,12 @@ def main(): mask_single_fiber[indices] = 1 # Get information in ventricles - indices = np.where((roi_md > args.md_ventricles) & - (roi_fa < args.fa_ventricles)) + indices = np.where((roi_md > args.md_min_ventricles) & + (roi_fa < args.fa_max_ventricles)) nb_voxels = roi_md[indices].shape[0] logging.info('Number of voxels found in ventricles (FA < {} and MD > {}): ' '{}' - .format(args.fa_ventricles, args.md_ventricles, nb_voxels)) + .format(args.fa_max_ventricles, args.md_min_ventricles, nb_voxels)) vent_avg = np.mean(roi_md[indices]) vent_std = np.std(roi_md[indices]) diff --git a/scripts/scil_sh_convert.py b/scripts/scil_sh_convert.py index 1251221c1..51d383d40 100755 --- a/scripts/scil_sh_convert.py +++ b/scripts/scil_sh_convert.py @@ -2,7 +2,7 @@ # -*- coding: utf-8 -*- """ -Convert a SH file between the two of the following base choices: +Convert a SH file between the two of the following basis choices: 'descoteaux07', 'descoteaux07_legacy', 'tournier07' or 'tournier07_legacy'. Using the sh_basis argument, both the input and the output SH bases must be given, in the order. For more information about the bases, see From 91e58561319a34a7b54bc1d8c0bb2831688fe735 Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Thu, 7 Mar 2024 13:22:23 -0500 Subject: [PATCH 12/12] Pep8 again --- scripts/scil_NODDI_priors.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/scripts/scil_NODDI_priors.py b/scripts/scil_NODDI_priors.py index 7158d643c..417f845f4 100755 --- a/scripts/scil_NODDI_priors.py +++ b/scripts/scil_NODDI_priors.py @@ -146,7 +146,8 @@ def main(): indices = np.where((roi_fa > args.fa_min_single_fiber) & (roi_fa < 0.95)) nb_voxels = roi_fa[indices].shape[0] logging.info('Number of voxels found in single fiber area (FA in range ' - '{}-{}]: {}'.format(args.fa_min_single_fiber, 0.95, nb_voxels)) + '{}-{}]: {}' + .format(args.fa_min_single_fiber, 0.95, nb_voxels)) single_fiber_ad_mean = np.mean(roi_ad[indices]) single_fiber_ad_std = np.std(roi_ad[indices]) single_fiber_rd_mean = np.mean(roi_rd[indices]) @@ -164,8 +165,8 @@ def main(): (roi_fa < args.fa_max_ventricles)) nb_voxels = roi_md[indices].shape[0] logging.info('Number of voxels found in ventricles (FA < {} and MD > {}): ' - '{}' - .format(args.fa_max_ventricles, args.md_min_ventricles, nb_voxels)) + '{}'.format(args.fa_max_ventricles, args.md_min_ventricles, + nb_voxels)) vent_avg = np.mean(roi_md[indices]) vent_std = np.std(roi_md[indices])