From 2ff000221886e1fe30f661dac0e08615cbd835ed Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 16 Nov 2023 15:16:37 -0500 Subject: [PATCH 01/97] project maps to streamline points and do math --- scripts/scil_project_map_to_streamlines.py | 158 ++++++++++++ scripts/scil_streamlines_point_math.py | 286 +++++++++++++++++++++ 2 files changed, 444 insertions(+) create mode 100755 scripts/scil_project_map_to_streamlines.py create mode 100755 scripts/scil_streamlines_point_math.py diff --git a/scripts/scil_project_map_to_streamlines.py b/scripts/scil_project_map_to_streamlines.py new file mode 100755 index 000000000..be45a93b1 --- /dev/null +++ b/scripts/scil_project_map_to_streamlines.py @@ -0,0 +1,158 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +Projects metrics extracted from a map onto the endpoints of streamlines. + +The default options will take data from a nifti image (3D or ND) and +project it onto the endpoints of streamlines. +""" + +import argparse +import logging +from copy import deepcopy + +from numpy import asarray, reshape, repeat + +from scilpy.io.image import assert_same_resolution, load_img +from scilpy.io.streamlines import load_tractogram_with_reference +from scilpy.io.utils import (add_overwrite_arg, + assert_inputs_exist, + assert_outputs_exist, + add_reference_arg, + add_verbose_arg) + +from scilpy.image.volume_space_management import DataVolume + +from dipy.io.streamline import save_tractogram, StatefulTractogram + +def project_metric_to_streamlines(sft, metric, endpoints_only=False): + """ + Projects a metric onto the points of streamlines. + + Parameters + ---------- + sft: StatefulTractogram + Input tractogram. + metric: DataVolume + Input metric. + + Optional: + --------- + endpoints_only: bool + If True, will only project the metric onto the endpoints of the + streamlines (all values along streamlines set to zero). If False, + will project the metric onto all points of the streamlines. + + Returns + ------- + streamline_data: + metric projected to each point of the streamlines. + """ + if len(metric.data.shape) == 4: + dimension = metric.data.shape[3] + else: + dimension = 1 + + streamline_data = [] + if endpoints_only: + for s in sft.streamlines: + p1_data = metric.get_value_at_coordinate(s[0][0], s[0][1], s[0][2], space=sft.space, origin=sft.origin) + p2_data = metric.get_value_at_coordinate(s[-1][0], s[-1][1], s[-1][2], space=sft.space, origin=sft.origin) + thisstreamline_data = [] + for p in s: + thisstreamline_data.append(asarray(repeat(0, p1_data.shape[0]))) + + thisstreamline_data[0] = p1_data + thisstreamline_data[-1] = p2_data + thisstreamline_data = asarray(thisstreamline_data) + + streamline_data.append(reshape(thisstreamline_data, (len(thisstreamline_data), dimension))) + else: + for s in sft.streamlines: + thisstreamline_data = [] + for p in s: + thisstreamline_data.append(metric.get_value_at_coordinate(p[0], p[1], p[2], space=sft.space, origin=sft.origin)) + + streamline_data.append(reshape(thisstreamline_data, (len(thisstreamline_data), dimension))) + + return streamline_data + +def _build_arg_parser(): + p = argparse.ArgumentParser( + description=__doc__, + formatter_class=argparse.RawTextHelpFormatter) + + # Mandatory arguments input and output tractogram must be in trk format + p.add_argument('in_tractogram', + help='Fiber bundle file.') + p.add_argument('in_metric', + help='Nifti metric to project onto streamlines.') + p.add_argument('out_tractogram', + help='Output file.') + + # Optional arguments + p.add_argument('--trilinear', action='store_true', + help='If set, will use trilinear interpolation \n' + 'else will use nearest neighbor interpolation \n' + 'by default.') + + p.add_argument('--endpoints_only', action='store_true', + help='If set, will only project the metric onto the endpoints \n' + 'of the streamlines (all other values along streamlines set to zero). \n' + 'If not set, will project the metric onto all points of the streamlines.') + + p.add_argument('--dpp_name', default='metric', + help='Name of the data_per_point to be saved in the output tractogram. \n' + '(Default: %(default)s)') + add_reference_arg(p) + add_overwrite_arg(p) + add_verbose_arg(p) + return p + + +def main(): + parser = _build_arg_parser() + args = parser.parse_args() + + assert_inputs_exist(parser, [args.in_tractogram, args.in_metric]) + + assert_outputs_exist(parser, args, [args.out_tractogram]) + + logging.debug("Loading the tractogram...") + sft = load_tractogram_with_reference(parser, args, args.in_tractogram) + sft.to_vox() + sft.to_corner() + + if len(sft.streamlines) == 0: + logging.warning('Empty bundle file {}. Skipping'.format(args.bundle)) + return + + logging.debug("Loading the metric...") + metric_img,metric_dtype = load_img(args.in_metric) + assert_same_resolution((args.in_tractogram, args.in_metric)) + metric_data = metric_img.get_fdata(caching='unchanged', dtype=float) + metric_res = metric_img.header.get_zooms()[:3] + + if args.trilinear: + interp = "trilinear" + else: + interp = "nearest" + + metric = DataVolume(metric_data, metric_res, interp) + + logging.debug("Projecting metric onto streamlines") + streamline_data = project_metric_to_streamlines(sft, metric, + endpoints_only=args.endpoints_only) + + logging.debug("Saving the tractogram...") + data_per_point = {} + data_per_point[args.dpp_name] = streamline_data + out_sft = StatefulTractogram(sft.streamlines, metric_img, + sft.space, sft.origin, data_per_point=data_per_point) + save_tractogram(out_sft, args.out_tractogram) + print(out_sft.data_per_point[args.dpp_name][0][0]) + print(out_sft.data_per_point[args.dpp_name][0][-1]) + +if __name__ == '__main__': + main() diff --git a/scripts/scil_streamlines_point_math.py b/scripts/scil_streamlines_point_math.py new file mode 100755 index 000000000..b70cbe47f --- /dev/null +++ b/scripts/scil_streamlines_point_math.py @@ -0,0 +1,286 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" +Performs an operation on data per point from streamlines +resulting in data per streamline. The supported +operations are: + +If the input is singular per point then: +mean: mean across points per streamline +sum: sum across points per streamline +min: min value across points per streamline +max: max value across points per streamline + +endpoints_only is set, min/max/mean/sum will only +be calculated using the streamline endpoints + +If the input is multivalued per point then: +mean: mean calculated per point for each streamline +sum: sum calculated per point for each streamline +min: min value calculated per point for each streamline +max: max value calculated per point for each streamline + +endpoints_only is set, min/max/mean/sum will only +be calculated at the streamline endpoints + +Endpoint only operations: +correlation: correlation calculated between arrays extracted from +streamline endpoints (data must be multivalued per point) +""" + +import argparse +import logging + +from copy import deepcopy + +from dipy.io.streamline import save_tractogram, StatefulTractogram +import numpy as np +from numpy import asarray, reshape + +from scilpy.io.streamlines import load_tractogram_with_reference +from scilpy.io.utils import (add_bbox_arg, + add_overwrite_arg, + add_reference_arg, + add_verbose_arg, + assert_inputs_exist, + assert_outputs_exist) + +def perform_operation_per_point(op_name, sft, dpp_name='metric', endpoints_only=False): + """Peforms an operation per point for all streamlines. + + Parameters + ---------- + op_name: str + A callable that takes a list of streamline data per point (4D) and + returns a list of streamline data per point. + sft: StatefulTractogram + The streamlines used in the operation. + dpp_name: str + The name of the data per point to be used in the operation. + endpoints_only: bool + If True, will only perform operation on endpoints + + Returns + ------- + new_sft: StatefulTractogram + sft with data per streamline resulting from the operation. + """ + # Performing operation + call_op = OPERATIONS[op_name] + if endpoints_only: + new_data_per_point = [] + for s in sft.data_per_point[dpp_name]: + this_data_per_point = [] + for p in s: + this_data_per_point.append(asarray(0.)) + this_data_per_point[0] = call_op(s[0]) + this_data_per_point[-1] = call_op(s[-1]) + new_data_per_point.append(reshape(this_data_per_point, (len(this_data_per_point), 1))) + else: + new_data_per_point = [] + for s in sft.data_per_point[dpp_name]: + this_data_per_point = [] + for p in s: + this_data_per_point.append(call_op(p)) + new_data_per_point.append(reshape(this_data_per_point, (len(this_data_per_point), 1))) + + # Extracting streamlines + return new_data_per_point + +def perform_operation_per_streamline(op_name, sft, dpp_name='metric', endpoints_only=False): + """Peforms an operation across points for each streamline. + + Parameters + ---------- + op_name: str + A callable that takes a list of streamline data per streamline and + returns a list of data per streamline. + sft: StatefulTractogram + The streamlines used in the operation. + dpp_name: str + The name of the data per point to be used in the operation. + endpoints_only: bool + If True, will only perform operation on endpoints + + Returns + ------- + new_sft: StatefulTractogram + sft with data per streamline resulting from the operation. + """ + # Performing operation + call_op = OPERATIONS[op_name] + if endpoints_only: + new_data_per_streamline = [] + for s in sft.data_per_point[dpp_name]: + new_data_per_streamline.append(call_op([s[0], s[-1]])) + else: + new_data_per_streamline = [] + for s in sft.data_per_point[dpp_name]: + new_data_per_streamline.append(call_op(s)) + + return new_data_per_streamline + +def perform_operation_across_endpoints(op_name, sft, dpp_name='metric'): + """Peforms an operation across endpoints for each streamline. + + Parameters + ---------- + op_name: str + A callable that takes a list of streamline data per streamline and + returns a list of data per streamline. + sft: StatefulTractogram + The streamlines used in the operation. + dpp_name: str + The name of the data per point to be used in the operation. + + Returns + ------- + new_sft: StatefulTractogram + sft with data per streamline resulting from the operation. + """ + # Performing operation + call_op = OPERATIONS[op_name] + new_data_per_streamline = [] + for s in sft.data_per_point[dpp_name]: + new_data_per_streamline.append(call_op(s[0], s[-1])[0,1]) + + return new_data_per_streamline + +def mean(array): + return np.mean(array) + +def sum(array): + return np.sum(array) + +def min(array): + return np.min(array) + +def max(array): + return np.max(array) + +def correlation(array1, array2): + return np.corrcoef(array1, array2) + +OPERATIONS = { + 'mean': mean, + 'sum': sum, + 'min': min, + 'max': max, + 'correlation': correlation, +} + +def _build_arg_parser(): + p = argparse.ArgumentParser( + formatter_class=argparse.RawTextHelpFormatter, + description=__doc__) + + p.add_argument('operation', metavar='OPERATION', + choices=['mean', 'sum', 'min', + 'max', 'correlation'], + help='The type of operation to be performed on the ' + 'streamlines. Must\nbe one of the following: ' + '%(choices)s.') + p.add_argument('in_tractogram', metavar='INPUT_FILE', + help='Input tractogram containing streamlines and ' + 'metadata.') + p.add_argument('out_tractogram', metavar='OUTPUT_FILE', + help='The file where the remaining streamlines ' + 'are saved.') + + p.add_argument('--endpoints_only', action='store_true', default=False, + help='If set, will only perform operation on endpoints \n' + 'If not set, will perform operation on all streamline points.') + p.add_argument('--dpp_name', default='metric', + help='Name of the data_per_point for operation to be ' + 'performed on. (Default: %(default)s)') + + p.add_argument('--output_dpp_name', default='metric_math', + help='Name of the resulting data_per_point to be saved in the output ' + 'tractogram. (Default: %(default)s)') + + add_reference_arg(p) + add_verbose_arg(p) + add_overwrite_arg(p) + add_bbox_arg(p) + + return p + +def main(): + parser = _build_arg_parser() + args = parser.parse_args() + + if args.verbose: + logging.getLogger().setLevel(logging.INFO) + + assert_inputs_exist(parser, args.in_tractogram) + assert_outputs_exist(parser, args, args.out_tractogram) + + # Load the input files. + logging.info("Loading file {}".format(args.in_tractogram)) + sft = load_tractogram_with_reference(parser, args, args.in_tractogram) + + if len(sft.streamlines) == 0: + logging.info("Input tractogram contains no streamlines. Exiting.") + return + + # Check to see if the data per point exists. + if args.dpp_name not in sft.data_per_point: + logging.info('Data per point {} not found in input tractogram.' + .format(args.dpp_name)) + return + + # check size of first streamline data_per_point + if sft.data_per_point[args.dpp_name][0].shape[1] == 1: + is_singular = True + else: + is_singular = False + + if args.operation == 'correlation' and is_singular: + logging.info('Correlation operation requires multivalued data per ' + 'point. Exiting.') + return + + # Perform the requested operation. + if not is_singular: + if args.operation == 'correlation': + logging.info('Performing {}} across endpoint data') + new_data_per_streamline = perform_operation_across_endpoints(args.operation,sft, args.dpp_name) + + # Adding data per streamline to new_sft + new_sft = StatefulTractogram(sft.streamlines, sft.space_attributes, + sft.space, sft.origin, + data_per_streamline={args.output_dpp_name: new_data_per_streamline}) + else: + # Results in new data per point + logging.info('Performing {} on data from each streamine point.'.format(args.operation)) + new_data_per_point = perform_operation_per_point( + args.operation, sft, args.dpp_name, args.endpoints_only) + + # Adding data per point to new_sft + new_sft = StatefulTractogram(sft.streamlines, sft.space_attributes, + sft.space, sft.origin, + data_per_point={args.output_dpp_name: new_data_per_point}) + else: + # Results in new data per streamline + logging.info('Performing {} across each streamline.'.format(args.operation)) + new_data_per_streamline = perform_operation_per_streamline( + args.operation, sft, args.output_dpp_name, args.endpoints_only) + + # Adding data per streamline to new_sft + new_sft = StatefulTractogram(sft.streamlines, sft.space_attributes, + sft.space, sft.origin, + data_per_streamline={args.dpp_name: new_data_per_streamline}) + + if len(new_sft) == 0 and not args.save_empty: + logging.info("Empty resulting tractogram. Not saving results.") + return + + # Save the new streamlines (and metadata) + logging.info('Saving {} streamlines to {}.'.format(len(new_sft), + args.out_tractogram)) + save_tractogram(new_sft, args.out_tractogram, + bbox_valid_check=args.bbox_check) + +if __name__ == "__main__": + main() From 7eef429d0f9da20e226c1a3354bffc765c81ea14 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 23 Nov 2023 13:34:49 -0500 Subject: [PATCH 02/97] fixed pep8 --- scripts/scil_project_map_to_streamlines.py | 55 +++++++++++++--------- scripts/scil_streamlines_point_math.py | 55 ++++++++++++++-------- 2 files changed, 68 insertions(+), 42 deletions(-) diff --git a/scripts/scil_project_map_to_streamlines.py b/scripts/scil_project_map_to_streamlines.py index be45a93b1..c9f0497d5 100755 --- a/scripts/scil_project_map_to_streamlines.py +++ b/scripts/scil_project_map_to_streamlines.py @@ -26,6 +26,7 @@ from dipy.io.streamline import save_tractogram, StatefulTractogram + def project_metric_to_streamlines(sft, metric, endpoints_only=False): """ Projects a metric onto the points of streamlines. @@ -43,7 +44,7 @@ def project_metric_to_streamlines(sft, metric, endpoints_only=False): If True, will only project the metric onto the endpoints of the streamlines (all values along streamlines set to zero). If False, will project the metric onto all points of the streamlines. - + Returns ------- streamline_data: @@ -57,38 +58,45 @@ def project_metric_to_streamlines(sft, metric, endpoints_only=False): streamline_data = [] if endpoints_only: for s in sft.streamlines: - p1_data = metric.get_value_at_coordinate(s[0][0], s[0][1], s[0][2], space=sft.space, origin=sft.origin) - p2_data = metric.get_value_at_coordinate(s[-1][0], s[-1][1], s[-1][2], space=sft.space, origin=sft.origin) + p1_data = metric.get_value_at_coordinate( + s[0][0], s[0][1], s[0][2], space=sft.space, origin=sft.origin) + p2_data = metric.get_value_at_coordinate( + s[-1][0], s[-1][1], s[-1][2], space=sft.space, origin=sft.origin) thisstreamline_data = [] for p in s: - thisstreamline_data.append(asarray(repeat(0, p1_data.shape[0]))) + thisstreamline_data.append( + asarray(repeat(0, p1_data.shape[0]))) thisstreamline_data[0] = p1_data thisstreamline_data[-1] = p2_data thisstreamline_data = asarray(thisstreamline_data) - streamline_data.append(reshape(thisstreamline_data, (len(thisstreamline_data), dimension))) + streamline_data.append( + reshape(thisstreamline_data, (len(thisstreamline_data), dimension))) else: for s in sft.streamlines: thisstreamline_data = [] for p in s: - thisstreamline_data.append(metric.get_value_at_coordinate(p[0], p[1], p[2], space=sft.space, origin=sft.origin)) + thisstreamline_data.append(metric.get_value_at_coordinate( + p[0], p[1], p[2], space=sft.space, origin=sft.origin)) - streamline_data.append(reshape(thisstreamline_data, (len(thisstreamline_data), dimension))) + streamline_data.append( + reshape(thisstreamline_data, (len(thisstreamline_data), dimension))) return streamline_data - + + def _build_arg_parser(): p = argparse.ArgumentParser( description=__doc__, formatter_class=argparse.RawTextHelpFormatter) # Mandatory arguments input and output tractogram must be in trk format - p.add_argument('in_tractogram', + p.add_argument('in_tractogram', help='Fiber bundle file.') p.add_argument('in_metric', help='Nifti metric to project onto streamlines.') - p.add_argument('out_tractogram', + p.add_argument('out_tractogram', help='Output file.') # Optional arguments @@ -96,11 +104,11 @@ def _build_arg_parser(): help='If set, will use trilinear interpolation \n' 'else will use nearest neighbor interpolation \n' 'by default.') - + p.add_argument('--endpoints_only', action='store_true', - help='If set, will only project the metric onto the endpoints \n' - 'of the streamlines (all other values along streamlines set to zero). \n' - 'If not set, will project the metric onto all points of the streamlines.') + help='If set, will only project the metric onto the endpoints \n' + 'of the streamlines (all other values along streamlines set to zero). \n' + 'If not set, will project the metric onto all points of the streamlines.') p.add_argument('--dpp_name', default='metric', help='Name of the data_per_point to be saved in the output tractogram. \n' @@ -116,7 +124,7 @@ def main(): args = parser.parse_args() assert_inputs_exist(parser, [args.in_tractogram, args.in_metric]) - + assert_outputs_exist(parser, args, [args.out_tractogram]) logging.debug("Loading the tractogram...") @@ -127,23 +135,23 @@ def main(): if len(sft.streamlines) == 0: logging.warning('Empty bundle file {}. Skipping'.format(args.bundle)) return - + logging.debug("Loading the metric...") - metric_img,metric_dtype = load_img(args.in_metric) + metric_img, metric_dtype = load_img(args.in_metric) assert_same_resolution((args.in_tractogram, args.in_metric)) metric_data = metric_img.get_fdata(caching='unchanged', dtype=float) metric_res = metric_img.header.get_zooms()[:3] - + if args.trilinear: - interp = "trilinear" + interp = "trilinear" else: interp = "nearest" - + metric = DataVolume(metric_data, metric_res, interp) - + logging.debug("Projecting metric onto streamlines") - streamline_data = project_metric_to_streamlines(sft, metric, - endpoints_only=args.endpoints_only) + streamline_data = project_metric_to_streamlines(sft, metric, + endpoints_only=args.endpoints_only) logging.debug("Saving the tractogram...") data_per_point = {} @@ -154,5 +162,6 @@ def main(): print(out_sft.data_per_point[args.dpp_name][0][0]) print(out_sft.data_per_point[args.dpp_name][0][-1]) + if __name__ == '__main__': main() diff --git a/scripts/scil_streamlines_point_math.py b/scripts/scil_streamlines_point_math.py index b70cbe47f..d3e30446d 100755 --- a/scripts/scil_streamlines_point_math.py +++ b/scripts/scil_streamlines_point_math.py @@ -46,6 +46,7 @@ assert_inputs_exist, assert_outputs_exist) + def perform_operation_per_point(op_name, sft, dpp_name='metric', endpoints_only=False): """Peforms an operation per point for all streamlines. @@ -60,7 +61,7 @@ def perform_operation_per_point(op_name, sft, dpp_name='metric', endpoints_only= The name of the data per point to be used in the operation. endpoints_only: bool If True, will only perform operation on endpoints - + Returns ------- new_sft: StatefulTractogram @@ -76,18 +77,21 @@ def perform_operation_per_point(op_name, sft, dpp_name='metric', endpoints_only= this_data_per_point.append(asarray(0.)) this_data_per_point[0] = call_op(s[0]) this_data_per_point[-1] = call_op(s[-1]) - new_data_per_point.append(reshape(this_data_per_point, (len(this_data_per_point), 1))) + new_data_per_point.append( + reshape(this_data_per_point, (len(this_data_per_point), 1))) else: new_data_per_point = [] for s in sft.data_per_point[dpp_name]: this_data_per_point = [] for p in s: this_data_per_point.append(call_op(p)) - new_data_per_point.append(reshape(this_data_per_point, (len(this_data_per_point), 1))) + new_data_per_point.append( + reshape(this_data_per_point, (len(this_data_per_point), 1))) # Extracting streamlines return new_data_per_point + def perform_operation_per_streamline(op_name, sft, dpp_name='metric', endpoints_only=False): """Peforms an operation across points for each streamline. @@ -102,7 +106,7 @@ def perform_operation_per_streamline(op_name, sft, dpp_name='metric', endpoints_ The name of the data per point to be used in the operation. endpoints_only: bool If True, will only perform operation on endpoints - + Returns ------- new_sft: StatefulTractogram @@ -121,6 +125,7 @@ def perform_operation_per_streamline(op_name, sft, dpp_name='metric', endpoints_ return new_data_per_streamline + def perform_operation_across_endpoints(op_name, sft, dpp_name='metric'): """Peforms an operation across endpoints for each streamline. @@ -133,7 +138,7 @@ def perform_operation_across_endpoints(op_name, sft, dpp_name='metric'): The streamlines used in the operation. dpp_name: str The name of the data per point to be used in the operation. - + Returns ------- new_sft: StatefulTractogram @@ -143,25 +148,31 @@ def perform_operation_across_endpoints(op_name, sft, dpp_name='metric'): call_op = OPERATIONS[op_name] new_data_per_streamline = [] for s in sft.data_per_point[dpp_name]: - new_data_per_streamline.append(call_op(s[0], s[-1])[0,1]) + new_data_per_streamline.append(call_op(s[0], s[-1])[0, 1]) return new_data_per_streamline + def mean(array): return np.mean(array) + def sum(array): return np.sum(array) + def min(array): return np.min(array) + def max(array): return np.max(array) + def correlation(array1, array2): return np.corrcoef(array1, array2) + OPERATIONS = { 'mean': mean, 'sum': sum, @@ -170,6 +181,7 @@ def correlation(array1, array2): 'correlation': correlation, } + def _build_arg_parser(): p = argparse.ArgumentParser( formatter_class=argparse.RawTextHelpFormatter, @@ -187,14 +199,14 @@ def _build_arg_parser(): p.add_argument('out_tractogram', metavar='OUTPUT_FILE', help='The file where the remaining streamlines ' 'are saved.') - + p.add_argument('--endpoints_only', action='store_true', default=False, - help='If set, will only perform operation on endpoints \n' - 'If not set, will perform operation on all streamline points.') + help='If set, will only perform operation on endpoints \n' + 'If not set, will perform operation on all streamline points.') p.add_argument('--dpp_name', default='metric', help='Name of the data_per_point for operation to be ' 'performed on. (Default: %(default)s)') - + p.add_argument('--output_dpp_name', default='metric_math', help='Name of the resulting data_per_point to be saved in the output ' 'tractogram. (Default: %(default)s)') @@ -206,6 +218,7 @@ def _build_arg_parser(): return p + def main(): parser = _build_arg_parser() args = parser.parse_args() @@ -227,9 +240,9 @@ def main(): # Check to see if the data per point exists. if args.dpp_name not in sft.data_per_point: logging.info('Data per point {} not found in input tractogram.' - .format(args.dpp_name)) + .format(args.dpp_name)) return - + # check size of first streamline data_per_point if sft.data_per_point[args.dpp_name][0].shape[1] == 1: is_singular = True @@ -238,22 +251,24 @@ def main(): if args.operation == 'correlation' and is_singular: logging.info('Correlation operation requires multivalued data per ' - 'point. Exiting.') + 'point. Exiting.') return - + # Perform the requested operation. if not is_singular: if args.operation == 'correlation': logging.info('Performing {}} across endpoint data') - new_data_per_streamline = perform_operation_across_endpoints(args.operation,sft, args.dpp_name) - + new_data_per_streamline = perform_operation_across_endpoints( + args.operation, sft, args.dpp_name) + # Adding data per streamline to new_sft new_sft = StatefulTractogram(sft.streamlines, sft.space_attributes, sft.space, sft.origin, data_per_streamline={args.output_dpp_name: new_data_per_streamline}) else: # Results in new data per point - logging.info('Performing {} on data from each streamine point.'.format(args.operation)) + logging.info( + 'Performing {} on data from each streamine point.'.format(args.operation)) new_data_per_point = perform_operation_per_point( args.operation, sft, args.dpp_name, args.endpoints_only) @@ -263,10 +278,11 @@ def main(): data_per_point={args.output_dpp_name: new_data_per_point}) else: # Results in new data per streamline - logging.info('Performing {} across each streamline.'.format(args.operation)) + logging.info( + 'Performing {} across each streamline.'.format(args.operation)) new_data_per_streamline = perform_operation_per_streamline( args.operation, sft, args.output_dpp_name, args.endpoints_only) - + # Adding data per streamline to new_sft new_sft = StatefulTractogram(sft.streamlines, sft.space_attributes, sft.space, sft.origin, @@ -282,5 +298,6 @@ def main(): save_tractogram(new_sft, args.out_tractogram, bbox_valid_check=args.bbox_check) + if __name__ == "__main__": main() From 4d64eb6f74b8f8ddc7a2730e14879b701d900f90 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 23 Nov 2023 13:36:18 -0500 Subject: [PATCH 03/97] fixed logging line --- scripts/scil_streamlines_point_math.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/scil_streamlines_point_math.py b/scripts/scil_streamlines_point_math.py index d3e30446d..52bab6c8a 100755 --- a/scripts/scil_streamlines_point_math.py +++ b/scripts/scil_streamlines_point_math.py @@ -257,7 +257,7 @@ def main(): # Perform the requested operation. if not is_singular: if args.operation == 'correlation': - logging.info('Performing {}} across endpoint data') + logging.info('Performing {} across endpoint data.'.format(args.operation)) new_data_per_streamline = perform_operation_across_endpoints( args.operation, sft, args.dpp_name) From 2fb63cd2f16c93c3ef0564ad27844b74b47d9e86 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 23 Nov 2023 13:49:59 -0500 Subject: [PATCH 04/97] fixed trailing whitespace --- scripts/scil_project_map_to_streamlines.py | 10 +++++----- scripts/scil_streamlines_point_math.py | 6 +++--- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/scripts/scil_project_map_to_streamlines.py b/scripts/scil_project_map_to_streamlines.py index c9f0497d5..59a9588bd 100755 --- a/scripts/scil_project_map_to_streamlines.py +++ b/scripts/scil_project_map_to_streamlines.py @@ -2,9 +2,9 @@ # -*- coding: utf-8 -*- """ -Projects metrics extracted from a map onto the endpoints of streamlines. +Projects metrics extracted from a map onto the endpoints of streamlines. -The default options will take data from a nifti image (3D or ND) and +The default options will take data from a nifti image (3D or ND) and project it onto the endpoints of streamlines. """ @@ -41,13 +41,13 @@ def project_metric_to_streamlines(sft, metric, endpoints_only=False): Optional: --------- endpoints_only: bool - If True, will only project the metric onto the endpoints of the - streamlines (all values along streamlines set to zero). If False, + If True, will only project the metric onto the endpoints of the + streamlines (all values along streamlines set to zero). If False, will project the metric onto all points of the streamlines. Returns ------- - streamline_data: + streamline_data: metric projected to each point of the streamlines. """ if len(metric.data.shape) == 4: diff --git a/scripts/scil_streamlines_point_math.py b/scripts/scil_streamlines_point_math.py index 52bab6c8a..a1487d733 100755 --- a/scripts/scil_streamlines_point_math.py +++ b/scripts/scil_streamlines_point_math.py @@ -2,7 +2,7 @@ # -*- coding: utf-8 -*- """ -Performs an operation on data per point from streamlines +Performs an operation on data per point from streamlines resulting in data per streamline. The supported operations are: @@ -12,10 +12,10 @@ min: min value across points per streamline max: max value across points per streamline -endpoints_only is set, min/max/mean/sum will only +endpoints_only is set, min/max/mean/sum will only be calculated using the streamline endpoints -If the input is multivalued per point then: +If the input is multivalued per point then: mean: mean calculated per point for each streamline sum: sum calculated per point for each streamline min: min value calculated per point for each streamline From ba5a516387654475b2df9edfd154677acc9131a8 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 23 Nov 2023 14:21:07 -0500 Subject: [PATCH 05/97] removed print statements --- scripts/scil_project_map_to_streamlines.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/scripts/scil_project_map_to_streamlines.py b/scripts/scil_project_map_to_streamlines.py index 59a9588bd..1e05cdad3 100755 --- a/scripts/scil_project_map_to_streamlines.py +++ b/scripts/scil_project_map_to_streamlines.py @@ -159,8 +159,6 @@ def main(): out_sft = StatefulTractogram(sft.streamlines, metric_img, sft.space, sft.origin, data_per_point=data_per_point) save_tractogram(out_sft, args.out_tractogram) - print(out_sft.data_per_point[args.dpp_name][0][0]) - print(out_sft.data_per_point[args.dpp_name][0][-1]) if __name__ == '__main__': From 797bd95451954e0ca96201c5d3ff3f2d9e2ee0f2 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 23 Nov 2023 14:37:24 -0500 Subject: [PATCH 06/97] removed useless assert resolution --- scripts/scil_project_map_to_streamlines.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/scripts/scil_project_map_to_streamlines.py b/scripts/scil_project_map_to_streamlines.py index 1e05cdad3..ddc22358f 100755 --- a/scripts/scil_project_map_to_streamlines.py +++ b/scripts/scil_project_map_to_streamlines.py @@ -14,7 +14,7 @@ from numpy import asarray, reshape, repeat -from scilpy.io.image import assert_same_resolution, load_img +from scilpy.io.image import load_img from scilpy.io.streamlines import load_tractogram_with_reference from scilpy.io.utils import (add_overwrite_arg, assert_inputs_exist, @@ -138,7 +138,6 @@ def main(): logging.debug("Loading the metric...") metric_img, metric_dtype = load_img(args.in_metric) - assert_same_resolution((args.in_tractogram, args.in_metric)) metric_data = metric_img.get_fdata(caching='unchanged', dtype=float) metric_res = metric_img.header.get_zooms()[:3] From cf12dc7bcc8f4a5d5d06a2181d000686adc3cca1 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Fri, 24 Nov 2023 09:55:30 -0500 Subject: [PATCH 07/97] changed to voxel world coords for interpolation --- scripts/scil_project_map_to_streamlines.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/scil_project_map_to_streamlines.py b/scripts/scil_project_map_to_streamlines.py index ddc22358f..2ef109064 100755 --- a/scripts/scil_project_map_to_streamlines.py +++ b/scripts/scil_project_map_to_streamlines.py @@ -129,7 +129,7 @@ def main(): logging.debug("Loading the tractogram...") sft = load_tractogram_with_reference(parser, args, args.in_tractogram) - sft.to_vox() + sft.to_voxmm() sft.to_corner() if len(sft.streamlines) == 0: From 5d014a82b6b9cccb204e2d86328a7e95cdc0ee58 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Fri, 24 Nov 2023 10:57:14 -0500 Subject: [PATCH 08/97] fixed dimension error on endpoints only single val --- scripts/scil_project_map_to_streamlines.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/scripts/scil_project_map_to_streamlines.py b/scripts/scil_project_map_to_streamlines.py index 2ef109064..bbb64acf1 100755 --- a/scripts/scil_project_map_to_streamlines.py +++ b/scripts/scil_project_map_to_streamlines.py @@ -64,8 +64,11 @@ def project_metric_to_streamlines(sft, metric, endpoints_only=False): s[-1][0], s[-1][1], s[-1][2], space=sft.space, origin=sft.origin) thisstreamline_data = [] for p in s: - thisstreamline_data.append( - asarray(repeat(0, p1_data.shape[0]))) + if dimension == 1: + thisstreamline_data.append(0) + else: + thisstreamline_data.append( + asarray(repeat(0, p1_data.shape[0]))) thisstreamline_data[0] = p1_data thisstreamline_data[-1] = p2_data From 6de542cf5726faa8efd2c177ac4d935c1e03c6b4 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Fri, 24 Nov 2023 10:58:07 -0500 Subject: [PATCH 09/97] remove deepcopy --- scripts/scil_project_map_to_streamlines.py | 1 - 1 file changed, 1 deletion(-) diff --git a/scripts/scil_project_map_to_streamlines.py b/scripts/scil_project_map_to_streamlines.py index bbb64acf1..e1764599d 100755 --- a/scripts/scil_project_map_to_streamlines.py +++ b/scripts/scil_project_map_to_streamlines.py @@ -10,7 +10,6 @@ import argparse import logging -from copy import deepcopy from numpy import asarray, reshape, repeat From ba839fb5cdacb6a1714fa58387d9d3f01ec17f3e Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Fri, 24 Nov 2023 11:09:40 -0500 Subject: [PATCH 10/97] fixed too long lines --- scripts/scil_project_map_to_streamlines.py | 28 +++++++++------- scripts/scil_streamlines_point_math.py | 37 +++++++++++++++------- 2 files changed, 42 insertions(+), 23 deletions(-) diff --git a/scripts/scil_project_map_to_streamlines.py b/scripts/scil_project_map_to_streamlines.py index e1764599d..fe4a05c3a 100755 --- a/scripts/scil_project_map_to_streamlines.py +++ b/scripts/scil_project_map_to_streamlines.py @@ -58,9 +58,11 @@ def project_metric_to_streamlines(sft, metric, endpoints_only=False): if endpoints_only: for s in sft.streamlines: p1_data = metric.get_value_at_coordinate( - s[0][0], s[0][1], s[0][2], space=sft.space, origin=sft.origin) + s[0][0], s[0][1], s[0][2], + space=sft.space, origin=sft.origin) p2_data = metric.get_value_at_coordinate( - s[-1][0], s[-1][1], s[-1][2], space=sft.space, origin=sft.origin) + s[-1][0], s[-1][1], s[-1][2], + space=sft.space, origin=sft.origin) thisstreamline_data = [] for p in s: if dimension == 1: @@ -74,7 +76,8 @@ def project_metric_to_streamlines(sft, metric, endpoints_only=False): thisstreamline_data = asarray(thisstreamline_data) streamline_data.append( - reshape(thisstreamline_data, (len(thisstreamline_data), dimension))) + reshape(thisstreamline_data, + (len(thisstreamline_data), dimension))) else: for s in sft.streamlines: thisstreamline_data = [] @@ -83,7 +86,8 @@ def project_metric_to_streamlines(sft, metric, endpoints_only=False): p[0], p[1], p[2], space=sft.space, origin=sft.origin)) streamline_data.append( - reshape(thisstreamline_data, (len(thisstreamline_data), dimension))) + reshape(thisstreamline_data, + (len(thisstreamline_data), dimension))) return streamline_data @@ -108,13 +112,14 @@ def _build_arg_parser(): 'by default.') p.add_argument('--endpoints_only', action='store_true', - help='If set, will only project the metric onto the endpoints \n' - 'of the streamlines (all other values along streamlines set to zero). \n' - 'If not set, will project the metric onto all points of the streamlines.') + help='If set, will only project the metric onto the \n' + 'endpoints of the streamlines (all other values along \n' + ' streamlines set to zero). If not set, will project \n' + ' the metric onto all points of the streamlines.') p.add_argument('--dpp_name', default='metric', - help='Name of the data_per_point to be saved in the output tractogram. \n' - '(Default: %(default)s)') + help='Name of the data_per_point to be saved in the \n' + 'output tractogram. (Default: %(default)s)') add_reference_arg(p) add_overwrite_arg(p) add_verbose_arg(p) @@ -152,13 +157,14 @@ def main(): logging.debug("Projecting metric onto streamlines") streamline_data = project_metric_to_streamlines(sft, metric, - endpoints_only=args.endpoints_only) + endpoints_only=args.endpoints_only) logging.debug("Saving the tractogram...") data_per_point = {} data_per_point[args.dpp_name] = streamline_data out_sft = StatefulTractogram(sft.streamlines, metric_img, - sft.space, sft.origin, data_per_point=data_per_point) + sft.space, sft.origin, + data_per_point=data_per_point) save_tractogram(out_sft, args.out_tractogram) diff --git a/scripts/scil_streamlines_point_math.py b/scripts/scil_streamlines_point_math.py index a1487d733..2b42ea0ea 100755 --- a/scripts/scil_streamlines_point_math.py +++ b/scripts/scil_streamlines_point_math.py @@ -47,7 +47,8 @@ assert_outputs_exist) -def perform_operation_per_point(op_name, sft, dpp_name='metric', endpoints_only=False): +def perform_operation_per_point(op_name, sft, dpp_name='metric', + endpoints_only=False): """Peforms an operation per point for all streamlines. Parameters @@ -92,7 +93,8 @@ def perform_operation_per_point(op_name, sft, dpp_name='metric', endpoints_only= return new_data_per_point -def perform_operation_per_streamline(op_name, sft, dpp_name='metric', endpoints_only=False): +def perform_operation_per_streamline(op_name, sft, dpp_name='metric', + endpoints_only=False): """Peforms an operation across points for each streamline. Parameters @@ -202,14 +204,15 @@ def _build_arg_parser(): p.add_argument('--endpoints_only', action='store_true', default=False, help='If set, will only perform operation on endpoints \n' - 'If not set, will perform operation on all streamline points.') + 'If not set, will perform operation on all streamline \n' + 'points.') p.add_argument('--dpp_name', default='metric', help='Name of the data_per_point for operation to be ' 'performed on. (Default: %(default)s)') p.add_argument('--output_dpp_name', default='metric_math', - help='Name of the resulting data_per_point to be saved in the output ' - 'tractogram. (Default: %(default)s)') + help='Name of the resulting data_per_point to be saved \n' + 'in the output tractogram. (Default: %(default)s)') add_reference_arg(p) add_verbose_arg(p) @@ -257,25 +260,33 @@ def main(): # Perform the requested operation. if not is_singular: if args.operation == 'correlation': - logging.info('Performing {} across endpoint data.'.format(args.operation)) + logging.info('Performing {} across endpoint data.'.format( + args.operation)) new_data_per_streamline = perform_operation_across_endpoints( args.operation, sft, args.dpp_name) # Adding data per streamline to new_sft - new_sft = StatefulTractogram(sft.streamlines, sft.space_attributes, + new_sft = StatefulTractogram(sft.streamlines, + sft.space_attributes, sft.space, sft.origin, - data_per_streamline={args.output_dpp_name: new_data_per_streamline}) + data_per_streamline={ + args.output_dpp_name: + new_data_per_streamline}) else: # Results in new data per point logging.info( - 'Performing {} on data from each streamine point.'.format(args.operation)) + 'Performing {} on data from each streamine point.'.format( + args.operation)) new_data_per_point = perform_operation_per_point( args.operation, sft, args.dpp_name, args.endpoints_only) # Adding data per point to new_sft - new_sft = StatefulTractogram(sft.streamlines, sft.space_attributes, + new_sft = StatefulTractogram(sft.streamlines, + sft.space_attributes, sft.space, sft.origin, - data_per_point={args.output_dpp_name: new_data_per_point}) + data_per_point={ + args.output_dpp_name: + new_data_per_point}) else: # Results in new data per streamline logging.info( @@ -286,7 +297,9 @@ def main(): # Adding data per streamline to new_sft new_sft = StatefulTractogram(sft.streamlines, sft.space_attributes, sft.space, sft.origin, - data_per_streamline={args.dpp_name: new_data_per_streamline}) + data_per_streamline={ + args.dpp_name: + new_data_per_streamline}) if len(new_sft) == 0 and not args.save_empty: logging.info("Empty resulting tractogram. Not saving results.") From a4af4acf6c2a741621fed66c8c0e48b9b080ca8a Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Fri, 24 Nov 2023 11:11:14 -0500 Subject: [PATCH 11/97] fixed indent --- scripts/scil_project_map_to_streamlines.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/scripts/scil_project_map_to_streamlines.py b/scripts/scil_project_map_to_streamlines.py index fe4a05c3a..04e08fc43 100755 --- a/scripts/scil_project_map_to_streamlines.py +++ b/scripts/scil_project_map_to_streamlines.py @@ -157,7 +157,8 @@ def main(): logging.debug("Projecting metric onto streamlines") streamline_data = project_metric_to_streamlines(sft, metric, - endpoints_only=args.endpoints_only) + endpoints_only= + args.endpoints_only) logging.debug("Saving the tractogram...") data_per_point = {} From fb0d9afe206c35993f1ec20b9531bd0424f2c6bf Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Fri, 24 Nov 2023 11:16:35 -0500 Subject: [PATCH 12/97] fixexed indents --- scripts/scil_project_map_to_streamlines.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/scripts/scil_project_map_to_streamlines.py b/scripts/scil_project_map_to_streamlines.py index 04e08fc43..818700f84 100755 --- a/scripts/scil_project_map_to_streamlines.py +++ b/scripts/scil_project_map_to_streamlines.py @@ -156,9 +156,9 @@ def main(): metric = DataVolume(metric_data, metric_res, interp) logging.debug("Projecting metric onto streamlines") - streamline_data = project_metric_to_streamlines(sft, metric, - endpoints_only= - args.endpoints_only) + streamline_data = project_metric_to_streamlines( + sft, metric, + endpoints_only=args.endpoints_only) logging.debug("Saving the tractogram...") data_per_point = {} From 76a78c44fea7f89435e3522277e852d4ebd20a47 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Wed, 29 Nov 2023 14:13:55 -0500 Subject: [PATCH 13/97] update output dpp name --- scripts/scil_streamlines_point_math.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/scil_streamlines_point_math.py b/scripts/scil_streamlines_point_math.py index 2b42ea0ea..5bbad7deb 100755 --- a/scripts/scil_streamlines_point_math.py +++ b/scripts/scil_streamlines_point_math.py @@ -292,13 +292,13 @@ def main(): logging.info( 'Performing {} across each streamline.'.format(args.operation)) new_data_per_streamline = perform_operation_per_streamline( - args.operation, sft, args.output_dpp_name, args.endpoints_only) + args.operation, sft, args.dpp_name, args.endpoints_only) # Adding data per streamline to new_sft new_sft = StatefulTractogram(sft.streamlines, sft.space_attributes, sft.space, sft.origin, data_per_streamline={ - args.dpp_name: + args.output_dpp_name: new_data_per_streamline}) if len(new_sft) == 0 and not args.save_empty: 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 14/97] 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 0e6fb5a5d406122f537a3fe6c1881a11c39e3e16 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 14 Dec 2023 09:56:18 -0500 Subject: [PATCH 15/97] moved project_metrics func to module --- scilpy/tractograms/streamline_operations.py | 65 +++++++++++++++++++ scripts/scil_project_map_to_streamlines.py | 70 +-------------------- 2 files changed, 66 insertions(+), 69 deletions(-) diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index c97936f44..0fea389e3 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -50,6 +50,71 @@ def _get_point_on_line(first_point, second_point, vox_lower_corner): return first_point + ray * (t0 + t1) / 2. +def project_metric_to_streamlines(sft, metric, endpoints_only=False): + """ + Projects a metric onto the points of streamlines. + + Parameters + ---------- + sft: StatefulTractogram + Input tractogram. + metric: DataVolume + Input metric. + + Optional: + --------- + endpoints_only: bool + If True, will only project the metric onto the endpoints of the + streamlines (all values along streamlines set to zero). If False, + will project the metric onto all points of the streamlines. + + Returns + ------- + streamline_data: + metric projected to each point of the streamlines. + """ + if len(metric.data.shape) == 4: + dimension = metric.data.shape[3] + else: + dimension = 1 + + streamline_data = [] + if endpoints_only: + for s in sft.streamlines: + p1_data = metric.get_value_at_coordinate( + s[0][0], s[0][1], s[0][2], + space=sft.space, origin=sft.origin) + p2_data = metric.get_value_at_coordinate( + s[-1][0], s[-1][1], s[-1][2], + space=sft.space, origin=sft.origin) + thisstreamline_data = [] + for p in s: + if dimension == 1: + thisstreamline_data.append(0) + else: + thisstreamline_data.append( + np.asarray(np.repeat(0, p1_data.shape[0]))) + + thisstreamline_data[0] = p1_data + thisstreamline_data[-1] = p2_data + thisstreamline_data = np.asarray(thisstreamline_data) + + streamline_data.append( + np.reshape(thisstreamline_data, + (len(thisstreamline_data), dimension))) + else: + for s in sft.streamlines: + thisstreamline_data = [] + for p in s: + thisstreamline_data.append(metric.get_value_at_coordinate( + p[0], p[1], p[2], space=sft.space, origin=sft.origin)) + + streamline_data.append( + np.reshape(thisstreamline_data, + (len(thisstreamline_data), dimension))) + + return streamline_data + def filter_streamlines_by_length(sft, min_length=0., max_length=np.inf): """ Filter streamlines using minimum and max length. diff --git a/scripts/scil_project_map_to_streamlines.py b/scripts/scil_project_map_to_streamlines.py index 818700f84..724ad5dd1 100755 --- a/scripts/scil_project_map_to_streamlines.py +++ b/scripts/scil_project_map_to_streamlines.py @@ -11,10 +11,9 @@ import argparse import logging -from numpy import asarray, reshape, repeat - from scilpy.io.image import load_img from scilpy.io.streamlines import load_tractogram_with_reference +from scilpy.tractograms.streamline_operations import project_metric_to_streamlines from scilpy.io.utils import (add_overwrite_arg, assert_inputs_exist, assert_outputs_exist, @@ -25,73 +24,6 @@ from dipy.io.streamline import save_tractogram, StatefulTractogram - -def project_metric_to_streamlines(sft, metric, endpoints_only=False): - """ - Projects a metric onto the points of streamlines. - - Parameters - ---------- - sft: StatefulTractogram - Input tractogram. - metric: DataVolume - Input metric. - - Optional: - --------- - endpoints_only: bool - If True, will only project the metric onto the endpoints of the - streamlines (all values along streamlines set to zero). If False, - will project the metric onto all points of the streamlines. - - Returns - ------- - streamline_data: - metric projected to each point of the streamlines. - """ - if len(metric.data.shape) == 4: - dimension = metric.data.shape[3] - else: - dimension = 1 - - streamline_data = [] - if endpoints_only: - for s in sft.streamlines: - p1_data = metric.get_value_at_coordinate( - s[0][0], s[0][1], s[0][2], - space=sft.space, origin=sft.origin) - p2_data = metric.get_value_at_coordinate( - s[-1][0], s[-1][1], s[-1][2], - space=sft.space, origin=sft.origin) - thisstreamline_data = [] - for p in s: - if dimension == 1: - thisstreamline_data.append(0) - else: - thisstreamline_data.append( - asarray(repeat(0, p1_data.shape[0]))) - - thisstreamline_data[0] = p1_data - thisstreamline_data[-1] = p2_data - thisstreamline_data = asarray(thisstreamline_data) - - streamline_data.append( - reshape(thisstreamline_data, - (len(thisstreamline_data), dimension))) - else: - for s in sft.streamlines: - thisstreamline_data = [] - for p in s: - thisstreamline_data.append(metric.get_value_at_coordinate( - p[0], p[1], p[2], space=sft.space, origin=sft.origin)) - - streamline_data.append( - reshape(thisstreamline_data, - (len(thisstreamline_data), dimension))) - - return streamline_data - - def _build_arg_parser(): p = argparse.ArgumentParser( description=__doc__, From 53ec9cc1ea79ff8bf21c9b1781e4421c1c168c98 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 14 Dec 2023 10:17:12 -0500 Subject: [PATCH 16/97] moved streamline ops to module --- scilpy/tractograms/streamline_operations.py | 269 +++++++++++++++----- scripts/scil_streamlines_point_math.py | 148 +---------- 2 files changed, 208 insertions(+), 209 deletions(-) diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index 0fea389e3..7b7dcab36 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -49,72 +49,6 @@ def _get_point_on_line(first_point, second_point, vox_lower_corner): return first_point + ray * (t0 + t1) / 2. - -def project_metric_to_streamlines(sft, metric, endpoints_only=False): - """ - Projects a metric onto the points of streamlines. - - Parameters - ---------- - sft: StatefulTractogram - Input tractogram. - metric: DataVolume - Input metric. - - Optional: - --------- - endpoints_only: bool - If True, will only project the metric onto the endpoints of the - streamlines (all values along streamlines set to zero). If False, - will project the metric onto all points of the streamlines. - - Returns - ------- - streamline_data: - metric projected to each point of the streamlines. - """ - if len(metric.data.shape) == 4: - dimension = metric.data.shape[3] - else: - dimension = 1 - - streamline_data = [] - if endpoints_only: - for s in sft.streamlines: - p1_data = metric.get_value_at_coordinate( - s[0][0], s[0][1], s[0][2], - space=sft.space, origin=sft.origin) - p2_data = metric.get_value_at_coordinate( - s[-1][0], s[-1][1], s[-1][2], - space=sft.space, origin=sft.origin) - thisstreamline_data = [] - for p in s: - if dimension == 1: - thisstreamline_data.append(0) - else: - thisstreamline_data.append( - np.asarray(np.repeat(0, p1_data.shape[0]))) - - thisstreamline_data[0] = p1_data - thisstreamline_data[-1] = p2_data - thisstreamline_data = np.asarray(thisstreamline_data) - - streamline_data.append( - np.reshape(thisstreamline_data, - (len(thisstreamline_data), dimension))) - else: - for s in sft.streamlines: - thisstreamline_data = [] - for p in s: - thisstreamline_data.append(metric.get_value_at_coordinate( - p[0], p[1], p[2], space=sft.space, origin=sft.origin)) - - streamline_data.append( - np.reshape(thisstreamline_data, - (len(thisstreamline_data), dimension))) - - return streamline_data - def filter_streamlines_by_length(sft, min_length=0., max_length=np.inf): """ Filter streamlines using minimum and max length. @@ -445,3 +379,206 @@ def smooth_line_spline(streamline, sigma, nb_ctrl_points): smoothed_streamline[-1] = streamline[-1] return smoothed_streamline + + +def perform_streamline_operation_per_point(op_name, sft, dpp_name='metric', + endpoints_only=False): + """Peforms an operation per point for all streamlines. + + Parameters + ---------- + op_name: str + A callable that takes a list of streamline data per point (4D) and + returns a list of streamline data per point. + sft: StatefulTractogram + The streamlines used in the operation. + dpp_name: str + The name of the data per point to be used in the operation. + endpoints_only: bool + If True, will only perform operation on endpoints + + Returns + ------- + new_sft: StatefulTractogram + sft with data per streamline resulting from the operation. + """ + # Performing operation + call_op = OPERATIONS[op_name] + if endpoints_only: + new_data_per_point = [] + for s in sft.data_per_point[dpp_name]: + this_data_per_point = [] + for p in s: + this_data_per_point.append(np.asarray(0.)) + this_data_per_point[0] = call_op(s[0]) + this_data_per_point[-1] = call_op(s[-1]) + new_data_per_point.append( + np.reshape(this_data_per_point, (len(this_data_per_point), 1))) + else: + new_data_per_point = [] + for s in sft.data_per_point[dpp_name]: + this_data_per_point = [] + for p in s: + this_data_per_point.append(call_op(p)) + new_data_per_point.append( + np.reshape(this_data_per_point, (len(this_data_per_point), 1))) + + # Extracting streamlines + return new_data_per_point + + +def perform_operation_per_streamline(op_name, sft, dpp_name='metric', + endpoints_only=False): + """Peforms an operation across points for each streamline. + + Parameters + ---------- + op_name: str + A callable that takes a list of streamline data per streamline and + returns a list of data per streamline. + sft: StatefulTractogram + The streamlines used in the operation. + dpp_name: str + The name of the data per point to be used in the operation. + endpoints_only: bool + If True, will only perform operation on endpoints + + Returns + ------- + new_sft: StatefulTractogram + sft with data per streamline resulting from the operation. + """ + # Performing operation + call_op = OPERATIONS[op_name] + if endpoints_only: + new_data_per_streamline = [] + for s in sft.data_per_point[dpp_name]: + new_data_per_streamline.append(call_op([s[0], s[-1]])) + else: + new_data_per_streamline = [] + for s in sft.data_per_point[dpp_name]: + new_data_per_streamline.append(call_op(s)) + + return new_data_per_streamline + + +def perform_streamline_operation_on_endpoints(op_name, sft, dpp_name='metric'): + """Peforms an operation across endpoints for each streamline. + + Parameters + ---------- + op_name: str + A callable that takes a list of streamline data per streamline and + returns a list of data per streamline. + sft: StatefulTractogram + The streamlines used in the operation. + dpp_name: str + The name of the data per point to be used in the operation. + + Returns + ------- + new_sft: StatefulTractogram + sft with data per streamline resulting from the operation. + """ + # Performing operation + call_op = OPERATIONS[op_name] + new_data_per_streamline = [] + for s in sft.data_per_point[dpp_name]: + new_data_per_streamline.append(call_op(s[0], s[-1])[0, 1]) + + return new_data_per_streamline + + +def mean(array): + return np.mean(array) + + +def sum(array): + return np.sum(array) + + +def min(array): + return np.min(array) + + +def max(array): + return np.max(array) + + +def correlation(array1, array2): + return np.corrcoef(array1, array2) + + +OPERATIONS = { + 'mean': mean, + 'sum': sum, + 'min': min, + 'max': max, + 'correlation': correlation, +} + + +def project_metric_to_streamlines(sft, metric, endpoints_only=False): + """ + Projects a metric onto the points of streamlines. + + Parameters + ---------- + sft: StatefulTractogram + Input tractogram. + metric: DataVolume + Input metric. + + Optional: + --------- + endpoints_only: bool + If True, will only project the metric onto the endpoints of the + streamlines (all values along streamlines set to zero). If False, + will project the metric onto all points of the streamlines. + + Returns + ------- + streamline_data: + metric projected to each point of the streamlines. + """ + if len(metric.data.shape) == 4: + dimension = metric.data.shape[3] + else: + dimension = 1 + + streamline_data = [] + if endpoints_only: + for s in sft.streamlines: + p1_data = metric.get_value_at_coordinate( + s[0][0], s[0][1], s[0][2], + space=sft.space, origin=sft.origin) + p2_data = metric.get_value_at_coordinate( + s[-1][0], s[-1][1], s[-1][2], + space=sft.space, origin=sft.origin) + thisstreamline_data = [] + for p in s: + if dimension == 1: + thisstreamline_data.append(0) + else: + thisstreamline_data.append( + np.asarray(np.repeat(0, p1_data.shape[0]))) + + thisstreamline_data[0] = p1_data + thisstreamline_data[-1] = p2_data + thisstreamline_data = np.asarray(thisstreamline_data) + + streamline_data.append( + np.reshape(thisstreamline_data, + (len(thisstreamline_data), dimension))) + else: + for s in sft.streamlines: + thisstreamline_data = [] + for p in s: + thisstreamline_data.append(metric.get_value_at_coordinate( + p[0], p[1], p[2], space=sft.space, origin=sft.origin)) + + streamline_data.append( + np.reshape(thisstreamline_data, + (len(thisstreamline_data), dimension))) + + return streamline_data diff --git a/scripts/scil_streamlines_point_math.py b/scripts/scil_streamlines_point_math.py index 5bbad7deb..463de8504 100755 --- a/scripts/scil_streamlines_point_math.py +++ b/scripts/scil_streamlines_point_math.py @@ -32,11 +32,7 @@ import argparse import logging -from copy import deepcopy - from dipy.io.streamline import save_tractogram, StatefulTractogram -import numpy as np -from numpy import asarray, reshape from scilpy.io.streamlines import load_tractogram_with_reference from scilpy.io.utils import (add_bbox_arg, @@ -46,143 +42,9 @@ assert_inputs_exist, assert_outputs_exist) - -def perform_operation_per_point(op_name, sft, dpp_name='metric', - endpoints_only=False): - """Peforms an operation per point for all streamlines. - - Parameters - ---------- - op_name: str - A callable that takes a list of streamline data per point (4D) and - returns a list of streamline data per point. - sft: StatefulTractogram - The streamlines used in the operation. - dpp_name: str - The name of the data per point to be used in the operation. - endpoints_only: bool - If True, will only perform operation on endpoints - - Returns - ------- - new_sft: StatefulTractogram - sft with data per streamline resulting from the operation. - """ - # Performing operation - call_op = OPERATIONS[op_name] - if endpoints_only: - new_data_per_point = [] - for s in sft.data_per_point[dpp_name]: - this_data_per_point = [] - for p in s: - this_data_per_point.append(asarray(0.)) - this_data_per_point[0] = call_op(s[0]) - this_data_per_point[-1] = call_op(s[-1]) - new_data_per_point.append( - reshape(this_data_per_point, (len(this_data_per_point), 1))) - else: - new_data_per_point = [] - for s in sft.data_per_point[dpp_name]: - this_data_per_point = [] - for p in s: - this_data_per_point.append(call_op(p)) - new_data_per_point.append( - reshape(this_data_per_point, (len(this_data_per_point), 1))) - - # Extracting streamlines - return new_data_per_point - - -def perform_operation_per_streamline(op_name, sft, dpp_name='metric', - endpoints_only=False): - """Peforms an operation across points for each streamline. - - Parameters - ---------- - op_name: str - A callable that takes a list of streamline data per streamline and - returns a list of data per streamline. - sft: StatefulTractogram - The streamlines used in the operation. - dpp_name: str - The name of the data per point to be used in the operation. - endpoints_only: bool - If True, will only perform operation on endpoints - - Returns - ------- - new_sft: StatefulTractogram - sft with data per streamline resulting from the operation. - """ - # Performing operation - call_op = OPERATIONS[op_name] - if endpoints_only: - new_data_per_streamline = [] - for s in sft.data_per_point[dpp_name]: - new_data_per_streamline.append(call_op([s[0], s[-1]])) - else: - new_data_per_streamline = [] - for s in sft.data_per_point[dpp_name]: - new_data_per_streamline.append(call_op(s)) - - return new_data_per_streamline - - -def perform_operation_across_endpoints(op_name, sft, dpp_name='metric'): - """Peforms an operation across endpoints for each streamline. - - Parameters - ---------- - op_name: str - A callable that takes a list of streamline data per streamline and - returns a list of data per streamline. - sft: StatefulTractogram - The streamlines used in the operation. - dpp_name: str - The name of the data per point to be used in the operation. - - Returns - ------- - new_sft: StatefulTractogram - sft with data per streamline resulting from the operation. - """ - # Performing operation - call_op = OPERATIONS[op_name] - new_data_per_streamline = [] - for s in sft.data_per_point[dpp_name]: - new_data_per_streamline.append(call_op(s[0], s[-1])[0, 1]) - - return new_data_per_streamline - - -def mean(array): - return np.mean(array) - - -def sum(array): - return np.sum(array) - - -def min(array): - return np.min(array) - - -def max(array): - return np.max(array) - - -def correlation(array1, array2): - return np.corrcoef(array1, array2) - - -OPERATIONS = { - 'mean': mean, - 'sum': sum, - 'min': min, - 'max': max, - 'correlation': correlation, -} - +from scilpy.tractograms.streamline_operations import (perform_streamline_operation_on_endpoints, + perform_streamline_operation_per_point, + perform_operation_per_streamline) def _build_arg_parser(): p = argparse.ArgumentParser( @@ -262,7 +124,7 @@ def main(): if args.operation == 'correlation': logging.info('Performing {} across endpoint data.'.format( args.operation)) - new_data_per_streamline = perform_operation_across_endpoints( + new_data_per_streamline = perform_streamline_operation_on_endpoints( args.operation, sft, args.dpp_name) # Adding data per streamline to new_sft @@ -277,7 +139,7 @@ def main(): logging.info( 'Performing {} on data from each streamine point.'.format( args.operation)) - new_data_per_point = perform_operation_per_point( + new_data_per_point = perform_streamline_operation_per_point( args.operation, sft, args.dpp_name, args.endpoints_only) # Adding data per point to new_sft From 061f3998123fa16ea34e48c1fd953989aa3426ea Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 14 Dec 2023 10:28:17 -0500 Subject: [PATCH 17/97] added blank unit tests --- scilpy/tractograms/streamline_operations.py | 1 + .../tests/tests_streamline_operations.py | 21 +++++++++++++++++++ 2 files changed, 22 insertions(+) diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index 7b7dcab36..552609d49 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -402,6 +402,7 @@ def perform_streamline_operation_per_point(op_name, sft, dpp_name='metric', new_sft: StatefulTractogram sft with data per streamline resulting from the operation. """ + # Performing operation call_op = OPERATIONS[op_name] if endpoints_only: diff --git a/scilpy/tractograms/tests/tests_streamline_operations.py b/scilpy/tractograms/tests/tests_streamline_operations.py index e0714f420..4029cd4bf 100644 --- a/scilpy/tractograms/tests/tests_streamline_operations.py +++ b/scilpy/tractograms/tests/tests_streamline_operations.py @@ -58,3 +58,24 @@ def test_resample_streamlines_step_size(): def compute_streamline_segment(): # toDo pass + + +def test_project_metric_to_streamlines(): + # toDo + pass + + +def test_perform_streamline_operation_per_point(): + # toDo + pass + + +def test_perform_operation_per_streamline(): + # toDo + pass + + +def test_perform_streamline_operation_on_endpoints(): + # toDo + pass + From c8b38895a1d076f16e041f8139203ee836a354c4 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 14 Dec 2023 11:24:05 -0500 Subject: [PATCH 18/97] added script test --- .../tests/test_project_map_to_streamlines.py | 84 +++++++++++++++++++ 1 file changed, 84 insertions(+) create mode 100644 scripts/tests/test_project_map_to_streamlines.py diff --git a/scripts/tests/test_project_map_to_streamlines.py b/scripts/tests/test_project_map_to_streamlines.py new file mode 100644 index 000000000..9f2f417f3 --- /dev/null +++ b/scripts/tests/test_project_map_to_streamlines.py @@ -0,0 +1,84 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import os +import tempfile + +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=['others.zip']) +tmp_dir = tempfile.TemporaryDirectory() + + +##### Deprecated file but it should still be running. +##### For more exhaustive tests, see test_tractogram_math.py + +def test_help_option(script_runner): + ret = script_runner.run('scil_project_map_to_streamlines.py', '--help') + assert ret.success + + +def test_execution_3D_map(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + in_fa = os.path.join(get_home(), 'others', 'fa.nii.gz') + in_tracto_1 = os.path.join(get_home(), 'others', + 'IFGWM_sub.trk') + + ret = script_runner.run('scil_project_map_to_streamlines.py', + in_tracto_1, in_fa, 'fa_on_streamlines.trk') + assert ret.success + +def test_execution_4D_map(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + in_rgb = os.path.join(get_home(), 'others', 'rgb.nii.gz') + in_tracto_1 = os.path.join(get_home(), 'others', + 'IFGWM_sub.trk') + + ret = script_runner.run('scil_project_map_to_streamlines.py', + in_tracto_1, in_rgb, 'rgb_on_streamlines.trk') + assert ret.success + +def test_execution_3D_map_endpoints_only(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + in_fa = os.path.join(get_home(), 'others', 'fa.nii.gz') + in_tracto_1 = os.path.join(get_home(), 'others', + 'IFGWM_sub.trk') + + ret = script_runner.run('scil_project_map_to_streamlines.py', + in_tracto_1, in_fa, 'fa_on_streamlines.trk', + '--endpoints_only') + assert ret.success + +def test_execution_4D_map_endpoints_only(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + in_rgb = os.path.join(get_home(), 'others', 'rgb.nii.gz') + in_tracto_1 = os.path.join(get_home(), 'others', + 'IFGWM_sub.trk') + + ret = script_runner.run('scil_project_map_to_streamlines.py', + in_tracto_1, in_rgb, 'rgb_on_streamlines.trk', + '--endpoints_only') + assert ret.success + +def test_execution_3D_map_dpp_name(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + in_fa = os.path.join(get_home(), 'others', 'fa.nii.gz') + in_tracto_1 = os.path.join(get_home(), 'others', + 'IFGWM_sub.trk') + + ret = script_runner.run('scil_project_map_to_streamlines.py', + in_tracto_1, in_fa, 'fa_on_streamlines.trk', + '--dpp_name', 'fa') + assert ret.success + +def test_execution_3D_map_trilinear(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + in_fa = os.path.join(get_home(), 'others', 'fa.nii.gz') + in_tracto_1 = os.path.join(get_home(), 'others', + 'IFGWM_sub.trk') + + ret = script_runner.run('scil_project_map_to_streamlines.py', + in_tracto_1, in_fa, 'fa_on_streamlines.trk', + '--trilinear') + assert ret.success \ No newline at end of file From 1f0ac45cab3b56fdcf162d497f35f235cd8e10a2 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 14 Dec 2023 11:56:09 -0500 Subject: [PATCH 19/97] renamed scripts --- ...il_streamlines_point_math.py => scil_tractogram_point_math.py} | 0 ...reamlines.py => scil_tractogram_project_map_to_streamlines.py} | 0 ...reamlines.py => test_tractogram_project_map_to_streamlines.py} | 0 3 files changed, 0 insertions(+), 0 deletions(-) rename scripts/{scil_streamlines_point_math.py => scil_tractogram_point_math.py} (100%) rename scripts/{scil_project_map_to_streamlines.py => scil_tractogram_project_map_to_streamlines.py} (100%) rename scripts/tests/{test_project_map_to_streamlines.py => test_tractogram_project_map_to_streamlines.py} (100%) diff --git a/scripts/scil_streamlines_point_math.py b/scripts/scil_tractogram_point_math.py similarity index 100% rename from scripts/scil_streamlines_point_math.py rename to scripts/scil_tractogram_point_math.py diff --git a/scripts/scil_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py similarity index 100% rename from scripts/scil_project_map_to_streamlines.py rename to scripts/scil_tractogram_project_map_to_streamlines.py diff --git a/scripts/tests/test_project_map_to_streamlines.py b/scripts/tests/test_tractogram_project_map_to_streamlines.py similarity index 100% rename from scripts/tests/test_project_map_to_streamlines.py rename to scripts/tests/test_tractogram_project_map_to_streamlines.py From 44fe6976bcdddbaffa03728272ed2f7865d02da9 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 14 Dec 2023 11:58:33 -0500 Subject: [PATCH 20/97] update naming in tests --- .../test_tractogram_project_map_to_streamlines.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/scripts/tests/test_tractogram_project_map_to_streamlines.py b/scripts/tests/test_tractogram_project_map_to_streamlines.py index 9f2f417f3..d6c744191 100644 --- a/scripts/tests/test_tractogram_project_map_to_streamlines.py +++ b/scripts/tests/test_tractogram_project_map_to_streamlines.py @@ -15,7 +15,7 @@ ##### For more exhaustive tests, see test_tractogram_math.py def test_help_option(script_runner): - ret = script_runner.run('scil_project_map_to_streamlines.py', '--help') + ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', '--help') assert ret.success @@ -25,7 +25,7 @@ def test_execution_3D_map(script_runner): in_tracto_1 = os.path.join(get_home(), 'others', 'IFGWM_sub.trk') - ret = script_runner.run('scil_project_map_to_streamlines.py', + ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', in_tracto_1, in_fa, 'fa_on_streamlines.trk') assert ret.success @@ -35,7 +35,7 @@ def test_execution_4D_map(script_runner): in_tracto_1 = os.path.join(get_home(), 'others', 'IFGWM_sub.trk') - ret = script_runner.run('scil_project_map_to_streamlines.py', + ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', in_tracto_1, in_rgb, 'rgb_on_streamlines.trk') assert ret.success @@ -45,7 +45,7 @@ def test_execution_3D_map_endpoints_only(script_runner): in_tracto_1 = os.path.join(get_home(), 'others', 'IFGWM_sub.trk') - ret = script_runner.run('scil_project_map_to_streamlines.py', + ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', in_tracto_1, in_fa, 'fa_on_streamlines.trk', '--endpoints_only') assert ret.success @@ -56,7 +56,7 @@ def test_execution_4D_map_endpoints_only(script_runner): in_tracto_1 = os.path.join(get_home(), 'others', 'IFGWM_sub.trk') - ret = script_runner.run('scil_project_map_to_streamlines.py', + ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', in_tracto_1, in_rgb, 'rgb_on_streamlines.trk', '--endpoints_only') assert ret.success @@ -67,7 +67,7 @@ def test_execution_3D_map_dpp_name(script_runner): in_tracto_1 = os.path.join(get_home(), 'others', 'IFGWM_sub.trk') - ret = script_runner.run('scil_project_map_to_streamlines.py', + ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', in_tracto_1, in_fa, 'fa_on_streamlines.trk', '--dpp_name', 'fa') assert ret.success @@ -78,7 +78,7 @@ def test_execution_3D_map_trilinear(script_runner): in_tracto_1 = os.path.join(get_home(), 'others', 'IFGWM_sub.trk') - ret = script_runner.run('scil_project_map_to_streamlines.py', + ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', in_tracto_1, in_fa, 'fa_on_streamlines.trk', '--trilinear') assert ret.success \ No newline at end of file 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 21/97] 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') From 16e308464e85d178c45d5e6dccf12a1aa5bb9ee9 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Wed, 10 Jan 2024 13:20:38 -0500 Subject: [PATCH 22/97] add blank unit tests to streamline operations --- .../tests/test_streamline_operations.py | 20 +++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/scilpy/tractograms/tests/test_streamline_operations.py b/scilpy/tractograms/tests/test_streamline_operations.py index faeb4e57c..b37c4aec7 100644 --- a/scilpy/tractograms/tests/test_streamline_operations.py +++ b/scilpy/tractograms/tests/test_streamline_operations.py @@ -367,3 +367,23 @@ def test_smooth_line_spline(): dist_2 = np.linalg.norm(noisy_streamline - smoothed_streamline) assert dist_1 < dist_2 + + +def test_project_metric_to_streamlines(): + # toDo + pass + + +def test_perform_streamline_operation_per_point(): + # toDo + pass + + +def test_perform_operation_per_streamline(): + # toDo + pass + + +def test_perform_streamline_operation_on_endpoints(): + # toDo + pass From b5660fde0aea3a94438ddabf7777f7c7598d7df9 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 11 Jan 2024 13:27:59 -0500 Subject: [PATCH 23/97] pep errors --- scilpy/tractograms/streamline_operations.py | 6 +++--- scripts/scil_tractogram_point_math.py | 16 +++++++++------- ...scil_tractogram_project_map_to_streamlines.py | 3 ++- ...test_tractogram_project_map_to_streamlines.py | 14 ++++++++------ 4 files changed, 22 insertions(+), 17 deletions(-) diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index 6fb3197a4..212ff0ef4 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -391,7 +391,7 @@ def smooth_line_spline(streamline, smoothing_parameter, nb_ctrl_points): def perform_streamline_operation_per_point(op_name, sft, dpp_name='metric', - endpoints_only=False): + endpoints_only=False): """Peforms an operation per point for all streamlines. Parameters @@ -411,7 +411,7 @@ def perform_streamline_operation_per_point(op_name, sft, dpp_name='metric', new_sft: StatefulTractogram sft with data per streamline resulting from the operation. """ - + # Performing operation call_op = OPERATIONS[op_name] if endpoints_only: @@ -438,7 +438,7 @@ def perform_streamline_operation_per_point(op_name, sft, dpp_name='metric', def perform_operation_per_streamline(op_name, sft, dpp_name='metric', - endpoints_only=False): + endpoints_only=False): """Peforms an operation across points for each streamline. Parameters diff --git a/scripts/scil_tractogram_point_math.py b/scripts/scil_tractogram_point_math.py index 463de8504..c9594d97c 100755 --- a/scripts/scil_tractogram_point_math.py +++ b/scripts/scil_tractogram_point_math.py @@ -42,9 +42,11 @@ assert_inputs_exist, assert_outputs_exist) -from scilpy.tractograms.streamline_operations import (perform_streamline_operation_on_endpoints, - perform_streamline_operation_per_point, - perform_operation_per_streamline) +from scilpy.tractograms.streamline_operations import ( + perform_streamline_operation_on_endpoints, + perform_streamline_operation_per_point, + perform_operation_per_streamline) + def _build_arg_parser(): p = argparse.ArgumentParser( @@ -124,7 +126,7 @@ def main(): if args.operation == 'correlation': logging.info('Performing {} across endpoint data.'.format( args.operation)) - new_data_per_streamline = perform_streamline_operation_on_endpoints( + new_dps= perform_streamline_operation_on_endpoints( args.operation, sft, args.dpp_name) # Adding data per streamline to new_sft @@ -133,13 +135,13 @@ def main(): sft.space, sft.origin, data_per_streamline={ args.output_dpp_name: - new_data_per_streamline}) + new_dps}) else: # Results in new data per point logging.info( 'Performing {} on data from each streamine point.'.format( args.operation)) - new_data_per_point = perform_streamline_operation_per_point( + new_dpp = perform_streamline_operation_per_point( args.operation, sft, args.dpp_name, args.endpoints_only) # Adding data per point to new_sft @@ -148,7 +150,7 @@ def main(): sft.space, sft.origin, data_per_point={ args.output_dpp_name: - new_data_per_point}) + new_dpp}) else: # Results in new data per streamline logging.info( diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index 724ad5dd1..78f11dc3d 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -13,7 +13,8 @@ from scilpy.io.image import load_img from scilpy.io.streamlines import load_tractogram_with_reference -from scilpy.tractograms.streamline_operations import project_metric_to_streamlines +from scilpy.tractograms.streamline_operations import ( + project_metric_to_streamlines) from scilpy.io.utils import (add_overwrite_arg, assert_inputs_exist, assert_outputs_exist, diff --git a/scripts/tests/test_tractogram_project_map_to_streamlines.py b/scripts/tests/test_tractogram_project_map_to_streamlines.py index d6c744191..95d99000c 100644 --- a/scripts/tests/test_tractogram_project_map_to_streamlines.py +++ b/scripts/tests/test_tractogram_project_map_to_streamlines.py @@ -10,12 +10,9 @@ fetch_data(get_testing_files_dict(), keys=['others.zip']) tmp_dir = tempfile.TemporaryDirectory() - -##### Deprecated file but it should still be running. -##### For more exhaustive tests, see test_tractogram_math.py - def test_help_option(script_runner): - ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', '--help') + ret = script_runner.run( + 'scil_tractogram_project_map_to_streamlines.py', '--help') assert ret.success @@ -29,6 +26,7 @@ def test_execution_3D_map(script_runner): in_tracto_1, in_fa, 'fa_on_streamlines.trk') assert ret.success + def test_execution_4D_map(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) in_rgb = os.path.join(get_home(), 'others', 'rgb.nii.gz') @@ -39,6 +37,7 @@ def test_execution_4D_map(script_runner): in_tracto_1, in_rgb, 'rgb_on_streamlines.trk') assert ret.success + def test_execution_3D_map_endpoints_only(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) in_fa = os.path.join(get_home(), 'others', 'fa.nii.gz') @@ -50,6 +49,7 @@ def test_execution_3D_map_endpoints_only(script_runner): '--endpoints_only') assert ret.success + def test_execution_4D_map_endpoints_only(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) in_rgb = os.path.join(get_home(), 'others', 'rgb.nii.gz') @@ -61,6 +61,7 @@ def test_execution_4D_map_endpoints_only(script_runner): '--endpoints_only') assert ret.success + def test_execution_3D_map_dpp_name(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) in_fa = os.path.join(get_home(), 'others', 'fa.nii.gz') @@ -72,6 +73,7 @@ def test_execution_3D_map_dpp_name(script_runner): '--dpp_name', 'fa') assert ret.success + def test_execution_3D_map_trilinear(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) in_fa = os.path.join(get_home(), 'others', 'fa.nii.gz') @@ -81,4 +83,4 @@ def test_execution_3D_map_trilinear(script_runner): ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', in_tracto_1, in_fa, 'fa_on_streamlines.trk', '--trilinear') - assert ret.success \ No newline at end of file + assert ret.success From a0fc9929e8ac21215e1f52a916d9eea08e46615a Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 11 Jan 2024 13:38:45 -0500 Subject: [PATCH 24/97] pep errors 2 --- scilpy/tractograms/streamline_operations.py | 4 ++-- scripts/scil_tractogram_point_math.py | 2 +- scripts/scil_tractogram_project_map_to_streamlines.py | 1 + scripts/tests/test_tractogram_project_map_to_streamlines.py | 1 + 4 files changed, 5 insertions(+), 3 deletions(-) diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index 212ff0ef4..20d8be343 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -391,7 +391,7 @@ def smooth_line_spline(streamline, smoothing_parameter, nb_ctrl_points): def perform_streamline_operation_per_point(op_name, sft, dpp_name='metric', - endpoints_only=False): + endpoints_only=False): """Peforms an operation per point for all streamlines. Parameters @@ -438,7 +438,7 @@ def perform_streamline_operation_per_point(op_name, sft, dpp_name='metric', def perform_operation_per_streamline(op_name, sft, dpp_name='metric', - endpoints_only=False): + endpoints_only=False): """Peforms an operation across points for each streamline. Parameters diff --git a/scripts/scil_tractogram_point_math.py b/scripts/scil_tractogram_point_math.py index c9594d97c..7adf7d61a 100755 --- a/scripts/scil_tractogram_point_math.py +++ b/scripts/scil_tractogram_point_math.py @@ -126,7 +126,7 @@ def main(): if args.operation == 'correlation': logging.info('Performing {} across endpoint data.'.format( args.operation)) - new_dps= perform_streamline_operation_on_endpoints( + new_dps = perform_streamline_operation_on_endpoints( args.operation, sft, args.dpp_name) # Adding data per streamline to new_sft diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index 78f11dc3d..34ef2152c 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -25,6 +25,7 @@ from dipy.io.streamline import save_tractogram, StatefulTractogram + def _build_arg_parser(): p = argparse.ArgumentParser( description=__doc__, diff --git a/scripts/tests/test_tractogram_project_map_to_streamlines.py b/scripts/tests/test_tractogram_project_map_to_streamlines.py index 95d99000c..c87c31b93 100644 --- a/scripts/tests/test_tractogram_project_map_to_streamlines.py +++ b/scripts/tests/test_tractogram_project_map_to_streamlines.py @@ -10,6 +10,7 @@ fetch_data(get_testing_files_dict(), keys=['others.zip']) tmp_dir = tempfile.TemporaryDirectory() + def test_help_option(script_runner): ret = script_runner.run( 'scil_tractogram_project_map_to_streamlines.py', '--help') From e1795fb3c2782501a1c05d3b6ec0a56e8030ab5f Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 15 Jan 2024 11:53:23 -0500 Subject: [PATCH 25/97] tests passed --- ...t_tractogram_project_map_to_streamlines.py | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/scripts/tests/test_tractogram_project_map_to_streamlines.py b/scripts/tests/test_tractogram_project_map_to_streamlines.py index c87c31b93..13dcf4fcb 100644 --- a/scripts/tests/test_tractogram_project_map_to_streamlines.py +++ b/scripts/tests/test_tractogram_project_map_to_streamlines.py @@ -19,12 +19,12 @@ def test_help_option(script_runner): def test_execution_3D_map(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) - in_fa = os.path.join(get_home(), 'others', 'fa.nii.gz') + in_t1 = os.path.join(get_home(), 'tractometry', 'mni_masked.nii.gz') in_tracto_1 = os.path.join(get_home(), 'others', 'IFGWM_sub.trk') ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_tracto_1, in_fa, 'fa_on_streamlines.trk') + in_tracto_1, in_t1, 't1_on_streamlines.trk') assert ret.success @@ -41,12 +41,12 @@ def test_execution_4D_map(script_runner): def test_execution_3D_map_endpoints_only(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) - in_fa = os.path.join(get_home(), 'others', 'fa.nii.gz') + in_t1 = os.path.join(get_home(), 'tractometry', 'mni_masked.nii.gz') in_tracto_1 = os.path.join(get_home(), 'others', 'IFGWM_sub.trk') ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_tracto_1, in_fa, 'fa_on_streamlines.trk', + in_tracto_1, in_t1, 't1_on_streamlines_endpoints.trk', '--endpoints_only') assert ret.success @@ -58,30 +58,30 @@ def test_execution_4D_map_endpoints_only(script_runner): 'IFGWM_sub.trk') ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_tracto_1, in_rgb, 'rgb_on_streamlines.trk', + in_tracto_1, in_rgb, 'rgb_on_streamlines_endpoints.trk', '--endpoints_only') assert ret.success def test_execution_3D_map_dpp_name(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) - in_fa = os.path.join(get_home(), 'others', 'fa.nii.gz') + in_t1 = os.path.join(get_home(), 'tractometry', 'mni_masked.nii.gz') in_tracto_1 = os.path.join(get_home(), 'others', 'IFGWM_sub.trk') ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_tracto_1, in_fa, 'fa_on_streamlines.trk', - '--dpp_name', 'fa') + in_tracto_1, in_t1, 't1_on_streamlines_dpp_name.trk', + '--dpp_name', 't1') assert ret.success def test_execution_3D_map_trilinear(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) - in_fa = os.path.join(get_home(), 'others', 'fa.nii.gz') + in_t1 = os.path.join(get_home(), 'tractometry', 'mni_masked.nii.gz') in_tracto_1 = os.path.join(get_home(), 'others', 'IFGWM_sub.trk') ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_tracto_1, in_fa, 'fa_on_streamlines.trk', + in_tracto_1, in_t1, 't1_on_streamlines_trilinear.trk', '--trilinear') assert ret.success From dded388d1b0474d9d08673dcfbb1cf185189b84a Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 15 Jan 2024 12:05:29 -0500 Subject: [PATCH 26/97] added test for point math --- scripts/tests/test_tractogram_point_math.py | 35 +++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 scripts/tests/test_tractogram_point_math.py diff --git a/scripts/tests/test_tractogram_point_math.py b/scripts/tests/test_tractogram_point_math.py new file mode 100644 index 000000000..db240cceb --- /dev/null +++ b/scripts/tests/test_tractogram_point_math.py @@ -0,0 +1,35 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import os +import tempfile + +from scilpy.io.fetcher import get_testing_files_dict, fetch_data, get_home + + +# If they already exist, this only takes 5 seconds (check md5sum) +fetch_data(get_testing_files_dict(), keys=['tractometry.zip']) +tmp_dir = tempfile.TemporaryDirectory() + + +def test_help_option(script_runner): + ret = script_runner.run('scil_tractogram_point_math.py', '--help') + assert ret.success + + +def test_execution_tractogram_point_math_mean_3D_defaults(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + in_bundle = os.path.join(get_home(), 'tractometry', + 'IFGWM_uni.trk') + in_t1 = os.path.join(get_home(), 'tractometry', + 'mni_masked.nii.gz') + + t1_on_bundle = 't1_on_streamlines.trk' + + script_runner.run('scil_tractogram_project_map_to_streamlines.py', + in_bundle, in_t1, t1_on_bundle) + + ret = script_runner.run('scil_tractogram_point_math.py', 'mean', t1_on_bundle, + 't1_mean_on_streamlines.trk') + + assert ret.success From 7592d60fd7ceb7255d10e8713ac0ad510569a196 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 15 Jan 2024 12:39:18 -0500 Subject: [PATCH 27/97] changed arg name for point math output --- scripts/scil_tractogram_point_math.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/scripts/scil_tractogram_point_math.py b/scripts/scil_tractogram_point_math.py index 7adf7d61a..6a22e6b3b 100755 --- a/scripts/scil_tractogram_point_math.py +++ b/scripts/scil_tractogram_point_math.py @@ -74,9 +74,10 @@ def _build_arg_parser(): help='Name of the data_per_point for operation to be ' 'performed on. (Default: %(default)s)') - p.add_argument('--output_dpp_name', default='metric_math', - help='Name of the resulting data_per_point to be saved \n' - 'in the output tractogram. (Default: %(default)s)') + p.add_argument('--output_name', default='metric_math', + help='Name of the resulting data_per_point or ' + 'data_per_streamline to be saved in the output ' + 'tractogram. (Default: %(default)s)') add_reference_arg(p) add_verbose_arg(p) @@ -134,7 +135,7 @@ def main(): sft.space_attributes, sft.space, sft.origin, data_per_streamline={ - args.output_dpp_name: + args.output_name: new_dps}) else: # Results in new data per point @@ -149,7 +150,7 @@ def main(): sft.space_attributes, sft.space, sft.origin, data_per_point={ - args.output_dpp_name: + args.output_name: new_dpp}) else: # Results in new data per streamline @@ -162,7 +163,7 @@ def main(): new_sft = StatefulTractogram(sft.streamlines, sft.space_attributes, sft.space, sft.origin, data_per_streamline={ - args.output_dpp_name: + args.output_name: new_data_per_streamline}) if len(new_sft) == 0 and not args.save_empty: From 8a6973d2a3d1705020eabd860a5e18a14f65e88e Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 15 Jan 2024 13:57:08 -0500 Subject: [PATCH 28/97] pep 8 long lines --- scripts/tests/test_tractogram_point_math.py | 4 +++- .../test_tractogram_project_map_to_streamlines.py | 12 ++++++++---- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/scripts/tests/test_tractogram_point_math.py b/scripts/tests/test_tractogram_point_math.py index db240cceb..90f1dc135 100644 --- a/scripts/tests/test_tractogram_point_math.py +++ b/scripts/tests/test_tractogram_point_math.py @@ -29,7 +29,9 @@ def test_execution_tractogram_point_math_mean_3D_defaults(script_runner): script_runner.run('scil_tractogram_project_map_to_streamlines.py', in_bundle, in_t1, t1_on_bundle) - ret = script_runner.run('scil_tractogram_point_math.py', 'mean', t1_on_bundle, + ret = script_runner.run('scil_tractogram_point_math.py', + 'mean', + t1_on_bundle, 't1_mean_on_streamlines.trk') assert ret.success diff --git a/scripts/tests/test_tractogram_project_map_to_streamlines.py b/scripts/tests/test_tractogram_project_map_to_streamlines.py index 13dcf4fcb..83adef48f 100644 --- a/scripts/tests/test_tractogram_project_map_to_streamlines.py +++ b/scripts/tests/test_tractogram_project_map_to_streamlines.py @@ -46,7 +46,8 @@ def test_execution_3D_map_endpoints_only(script_runner): 'IFGWM_sub.trk') ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_tracto_1, in_t1, 't1_on_streamlines_endpoints.trk', + in_tracto_1, in_t1, + 't1_on_streamlines_endpoints.trk', '--endpoints_only') assert ret.success @@ -58,7 +59,8 @@ def test_execution_4D_map_endpoints_only(script_runner): 'IFGWM_sub.trk') ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_tracto_1, in_rgb, 'rgb_on_streamlines_endpoints.trk', + in_tracto_1, in_rgb, + 'rgb_on_streamlines_endpoints.trk', '--endpoints_only') assert ret.success @@ -70,7 +72,8 @@ def test_execution_3D_map_dpp_name(script_runner): 'IFGWM_sub.trk') ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_tracto_1, in_t1, 't1_on_streamlines_dpp_name.trk', + in_tracto_1, in_t1, + 't1_on_streamlines_dpp_name.trk', '--dpp_name', 't1') assert ret.success @@ -82,6 +85,7 @@ def test_execution_3D_map_trilinear(script_runner): 'IFGWM_sub.trk') ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_tracto_1, in_t1, 't1_on_streamlines_trilinear.trk', + in_tracto_1, in_t1, + 't1_on_streamlines_trilinear.trk', '--trilinear') assert ret.success From 62c103c48b85699777a3749609fdf8df44eec836 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 15 Jan 2024 15:16:45 -0500 Subject: [PATCH 29/97] updated function names in stream_ops --- scilpy/tractograms/streamline_operations.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index 20d8be343..fe727426d 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -499,32 +499,32 @@ def perform_streamline_operation_on_endpoints(op_name, sft, dpp_name='metric'): return new_data_per_streamline -def mean(array): +def stream_mean(array): return np.mean(array) -def sum(array): +def stream_sum(array): return np.sum(array) -def min(array): +def stream_min(array): return np.min(array) -def max(array): +def stream_max(array): return np.max(array) -def correlation(array1, array2): +def stream_correlation(array1, array2): return np.corrcoef(array1, array2) OPERATIONS = { - 'mean': mean, - 'sum': sum, - 'min': min, - 'max': max, - 'correlation': correlation, + 'mean': stream_mean, + 'sum': stream_sum, + 'min': stream_min, + 'max': stream_max, + 'correlation': stream_correlation, } From ffc2c88e07e7c4f4027a35a2a6df890f49fa15e2 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 25 Jan 2024 12:56:08 -0500 Subject: [PATCH 30/97] Fixed for loops and added multi input --- scilpy/tractograms/streamline_operations.py | 14 +- scripts/scil_tractogram_point_math.py | 128 +++++++++--------- ...l_tractogram_project_map_to_streamlines.py | 72 +++++----- scripts/tests/.python-version | 1 + 4 files changed, 111 insertions(+), 104 deletions(-) create mode 100644 scripts/tests/.python-version diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index fe727426d..becc00029 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -417,9 +417,7 @@ def perform_streamline_operation_per_point(op_name, sft, dpp_name='metric', if endpoints_only: new_data_per_point = [] for s in sft.data_per_point[dpp_name]: - this_data_per_point = [] - for p in s: - this_data_per_point.append(np.asarray(0.)) + this_data_per_point = np.nan * np.ones((len(s), 1)) this_data_per_point[0] = call_op(s[0]) this_data_per_point[-1] = call_op(s[-1]) new_data_per_point.append( @@ -566,12 +564,10 @@ def project_metric_to_streamlines(sft, metric, endpoints_only=False): s[-1][0], s[-1][1], s[-1][2], space=sft.space, origin=sft.origin) thisstreamline_data = [] - for p in s: - if dimension == 1: - thisstreamline_data.append(0) - else: - thisstreamline_data.append( - np.asarray(np.repeat(0, p1_data.shape[0]))) + if dimension == 1: + thisstreamline_data = np.ones((len(s), 1)) * np.nan + else: + thisstreamline_data = np.ones((len(s), p1_data.shape[0])) * np.nan thisstreamline_data[0] = p1_data thisstreamline_data[-1] = p2_data diff --git a/scripts/scil_tractogram_point_math.py b/scripts/scil_tractogram_point_math.py index 6a22e6b3b..a7f1b3abe 100755 --- a/scripts/scil_tractogram_point_math.py +++ b/scripts/scil_tractogram_point_math.py @@ -62,6 +62,15 @@ def _build_arg_parser(): p.add_argument('in_tractogram', metavar='INPUT_FILE', help='Input tractogram containing streamlines and ' 'metadata.') + p.add_argument('--in_dpp_name', nargs='+', action='append', + required=True, + help='Name of the data_per_point for operation to be ' + 'performed on.') + p.add_argument('--out_name', nargs='+', action='append', + required=True, + help='Name of the resulting data_per_point or ' + 'data_per_streamline to be saved in the output ' + 'tractogram. (Default: %(default)s)') p.add_argument('out_tractogram', metavar='OUTPUT_FILE', help='The file where the remaining streamlines ' 'are saved.') @@ -70,14 +79,7 @@ def _build_arg_parser(): help='If set, will only perform operation on endpoints \n' 'If not set, will perform operation on all streamline \n' 'points.') - p.add_argument('--dpp_name', default='metric', - help='Name of the data_per_point for operation to be ' - 'performed on. (Default: %(default)s)') - p.add_argument('--output_name', default='metric_math', - help='Name of the resulting data_per_point or ' - 'data_per_streamline to be saved in the output ' - 'tractogram. (Default: %(default)s)') add_reference_arg(p) add_verbose_arg(p) @@ -105,66 +107,66 @@ def main(): logging.info("Input tractogram contains no streamlines. Exiting.") return - # Check to see if the data per point exists. - if args.dpp_name not in sft.data_per_point: - logging.info('Data per point {} not found in input tractogram.' - .format(args.dpp_name)) - return - - # check size of first streamline data_per_point - if sft.data_per_point[args.dpp_name][0].shape[1] == 1: - is_singular = True - else: - is_singular = False - - if args.operation == 'correlation' and is_singular: - logging.info('Correlation operation requires multivalued data per ' - 'point. Exiting.') - return + if len(args.in_dpp_name[0]) != len(args.out_name[0]): + parser.error('The number of in_dpp_names and out_names must be ' + 'the same.') + + data_per_point = {} + data_per_streamline = {} + for in_dpp_name, out_name in zip(args.in_dpp_name[0], + args.out_name[0]): + # Check to see if the data per point exists. + if in_dpp_name not in sft.data_per_point: + logging.info('Data per point {} not found in input tractogram.' + .format(in_dpp_name)) + return + + # check size of first streamline data_per_point + if sft.data_per_point[in_dpp_name][0].shape[1] == 1: + is_singular = True + else: + is_singular = False - # Perform the requested operation. - if not is_singular: - if args.operation == 'correlation': - logging.info('Performing {} across endpoint data.'.format( - args.operation)) - new_dps = perform_streamline_operation_on_endpoints( - args.operation, sft, args.dpp_name) + if args.operation == 'correlation' and is_singular: + logging.info('Correlation operation requires multivalued data per ' + 'point. Exiting.') + return - # Adding data per streamline to new_sft - new_sft = StatefulTractogram(sft.streamlines, - sft.space_attributes, - sft.space, sft.origin, - data_per_streamline={ - args.output_name: - new_dps}) + # Perform the requested operation. + if not is_singular: + if args.operation == 'correlation': + logging.info('Performing {} across endpoint data.'.format( + args.operation)) + new_dps = perform_streamline_operation_on_endpoints( + args.operation, sft, in_dpp_name) + + # Adding data per streamline to new_sft + data_per_streamline[out_name] = new_dps + else: + # Results in new data per point + logging.info( + 'Performing {} on data from each streamine point.'.format( + args.operation)) + new_dpp = perform_streamline_operation_per_point( + args.operation, sft, in_dpp_name, args.endpoints_only) + + # Adding data per point to new_sft + data_per_point[out_name] = new_dpp else: - # Results in new data per point + # Results in new data per streamline logging.info( - 'Performing {} on data from each streamine point.'.format( - args.operation)) - new_dpp = perform_streamline_operation_per_point( - args.operation, sft, args.dpp_name, args.endpoints_only) - - # Adding data per point to new_sft - new_sft = StatefulTractogram(sft.streamlines, - sft.space_attributes, - sft.space, sft.origin, - data_per_point={ - args.output_name: - new_dpp}) - else: - # Results in new data per streamline - logging.info( - 'Performing {} across each streamline.'.format(args.operation)) - new_data_per_streamline = perform_operation_per_streamline( - args.operation, sft, args.dpp_name, args.endpoints_only) - - # Adding data per streamline to new_sft - new_sft = StatefulTractogram(sft.streamlines, sft.space_attributes, - sft.space, sft.origin, - data_per_streamline={ - args.output_name: - new_data_per_streamline}) + 'Performing {} across each streamline.'.format(args.operation)) + new_data_per_streamline = perform_operation_per_streamline( + args.operation, sft, in_dpp_name, args.endpoints_only) + + # Adding data per streamline to new_sft + data_per_streamline[out_name] = new_data_per_streamline + + new_sft = StatefulTractogram(sft.streamlines, + sft.space_attributes, + sft.space, sft.origin, + data_per_point=data_per_point, + data_per_streamline=data_per_streamline) if len(new_sft) == 0 and not args.save_empty: logging.info("Empty resulting tractogram. Not saving results.") diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index 34ef2152c..ef001c9c6 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -11,20 +11,18 @@ import argparse import logging +from dipy.io.streamline import save_tractogram, StatefulTractogram + from scilpy.io.image import load_img from scilpy.io.streamlines import load_tractogram_with_reference -from scilpy.tractograms.streamline_operations import ( - project_metric_to_streamlines) from scilpy.io.utils import (add_overwrite_arg, - assert_inputs_exist, - assert_outputs_exist, add_reference_arg, - add_verbose_arg) - + add_verbose_arg, + assert_inputs_exist, + assert_outputs_exist) from scilpy.image.volume_space_management import DataVolume - -from dipy.io.streamline import save_tractogram, StatefulTractogram - +from scilpy.tractograms.streamline_operations import ( + project_metric_to_streamlines) def _build_arg_parser(): p = argparse.ArgumentParser( @@ -34,10 +32,13 @@ def _build_arg_parser(): # Mandatory arguments input and output tractogram must be in trk format p.add_argument('in_tractogram', help='Fiber bundle file.') - p.add_argument('in_metric', - help='Nifti metric to project onto streamlines.') p.add_argument('out_tractogram', help='Output file.') + p.add_argument('--in_metric', nargs='+', action='append', required=True, + help='Nifti metric to project onto streamlines.') + p.add_argument('--out_dpp_name', nargs='+', action='append', required=True, + help='Name of the data_per_point to be saved in the \n' + 'output tractogram.') # Optional arguments p.add_argument('--trilinear', action='store_true', @@ -48,12 +49,9 @@ def _build_arg_parser(): p.add_argument('--endpoints_only', action='store_true', help='If set, will only project the metric onto the \n' 'endpoints of the streamlines (all other values along \n' - ' streamlines set to zero). If not set, will project \n' + ' streamlines will be NaN). If not set, will project \n' ' the metric onto all points of the streamlines.') - p.add_argument('--dpp_name', default='metric', - help='Name of the data_per_point to be saved in the \n' - 'output tractogram. (Default: %(default)s)') add_reference_arg(p) add_overwrite_arg(p) add_verbose_arg(p) @@ -64,7 +62,10 @@ def main(): parser = _build_arg_parser() args = parser.parse_args() - assert_inputs_exist(parser, [args.in_tractogram, args.in_metric]) + input_file_list = [args.in_tractogram] + for fmetric in args.in_metric: + input_file_list += fmetric + assert_inputs_exist(parser, input_file_list) assert_outputs_exist(parser, args, [args.out_tractogram]) @@ -77,26 +78,33 @@ def main(): logging.warning('Empty bundle file {}. Skipping'.format(args.bundle)) return - logging.debug("Loading the metric...") - metric_img, metric_dtype = load_img(args.in_metric) - metric_data = metric_img.get_fdata(caching='unchanged', dtype=float) - metric_res = metric_img.header.get_zooms()[:3] + # Check to see if the number of metrics and dpp_names are the same + if len(args.in_metric) != len(args.out_dpp_name): + parser.error('The number of metrics and dpp_names must be the same.') - if args.trilinear: - interp = "trilinear" - else: - interp = "nearest" + data_per_point = {} + for fmetric, dpp_name in zip(args.in_metric[0], args.out_dpp_name[0]): + logging.debug("Loading the metric...") + metric_img, metric_dtype = load_img(fmetric) + metric_data = metric_img.get_fdata(caching='unchanged', dtype=float) + metric_res = metric_img.header.get_zooms()[:3] - metric = DataVolume(metric_data, metric_res, interp) + if args.trilinear: + interp = "trilinear" + else: + interp = "nearest" - logging.debug("Projecting metric onto streamlines") - streamline_data = project_metric_to_streamlines( - sft, metric, - endpoints_only=args.endpoints_only) + metric = DataVolume(metric_data, metric_res, interp) + + logging.debug("Projecting metric onto streamlines") + streamline_data = project_metric_to_streamlines( + sft, metric, + endpoints_only=args.endpoints_only) + + logging.debug("Saving the tractogram...") + + data_per_point[dpp_name] = streamline_data - logging.debug("Saving the tractogram...") - data_per_point = {} - data_per_point[args.dpp_name] = streamline_data out_sft = StatefulTractogram(sft.streamlines, metric_img, sft.space, sft.origin, data_per_point=data_per_point) diff --git a/scripts/tests/.python-version b/scripts/tests/.python-version new file mode 100644 index 000000000..c8cfe3959 --- /dev/null +++ b/scripts/tests/.python-version @@ -0,0 +1 @@ +3.10 From 89f40d2f5dbdb219effe42589d066fdc5090f0d3 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 25 Jan 2024 13:02:38 -0500 Subject: [PATCH 31/97] updated projection tests with new required args --- ...t_tractogram_project_map_to_streamlines.py | 34 ++++++++----------- 1 file changed, 15 insertions(+), 19 deletions(-) diff --git a/scripts/tests/test_tractogram_project_map_to_streamlines.py b/scripts/tests/test_tractogram_project_map_to_streamlines.py index 83adef48f..2ed5b9b56 100644 --- a/scripts/tests/test_tractogram_project_map_to_streamlines.py +++ b/scripts/tests/test_tractogram_project_map_to_streamlines.py @@ -24,7 +24,9 @@ def test_execution_3D_map(script_runner): 'IFGWM_sub.trk') ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_tracto_1, in_t1, 't1_on_streamlines.trk') + in_tracto_1, 't1_on_streamlines.trk', + '--in_metric', in_t1, + '--out_dpp_name', 't1') assert ret.success @@ -35,7 +37,9 @@ def test_execution_4D_map(script_runner): 'IFGWM_sub.trk') ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_tracto_1, in_rgb, 'rgb_on_streamlines.trk') + in_tracto_1, 'rgb_on_streamlines.trk', + '--in_metric', in_rgb, + '--out_dpp_name', 'rgb') assert ret.success @@ -46,8 +50,10 @@ def test_execution_3D_map_endpoints_only(script_runner): 'IFGWM_sub.trk') ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_tracto_1, in_t1, + in_tracto_1, 't1_on_streamlines_endpoints.trk', + '--in_metric', in_t1, + '--out_dpp_name', 't1', '--endpoints_only') assert ret.success @@ -59,25 +65,13 @@ def test_execution_4D_map_endpoints_only(script_runner): 'IFGWM_sub.trk') ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_tracto_1, in_rgb, + in_tracto_1, 'rgb_on_streamlines_endpoints.trk', + '--in_metric', in_rgb, + '--out_dpp_name', 'rgb', '--endpoints_only') assert ret.success - -def test_execution_3D_map_dpp_name(script_runner): - os.chdir(os.path.expanduser(tmp_dir.name)) - in_t1 = os.path.join(get_home(), 'tractometry', 'mni_masked.nii.gz') - in_tracto_1 = os.path.join(get_home(), 'others', - 'IFGWM_sub.trk') - - ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_tracto_1, in_t1, - 't1_on_streamlines_dpp_name.trk', - '--dpp_name', 't1') - assert ret.success - - def test_execution_3D_map_trilinear(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) in_t1 = os.path.join(get_home(), 'tractometry', 'mni_masked.nii.gz') @@ -85,7 +79,9 @@ def test_execution_3D_map_trilinear(script_runner): 'IFGWM_sub.trk') ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_tracto_1, in_t1, + in_tracto_1, 't1_on_streamlines_trilinear.trk', + '--in_metric', in_t1, + '--out_dpp_name', 't1', '--trilinear') assert ret.success From e8bcaddfec51c92a8566f4d0c338a5c1f3c9e7a7 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 25 Jan 2024 13:28:16 -0500 Subject: [PATCH 32/97] updated point math test with req args --- scripts/tests/test_tractogram_point_math.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/scripts/tests/test_tractogram_point_math.py b/scripts/tests/test_tractogram_point_math.py index 90f1dc135..fa1b871da 100644 --- a/scripts/tests/test_tractogram_point_math.py +++ b/scripts/tests/test_tractogram_point_math.py @@ -27,11 +27,15 @@ def test_execution_tractogram_point_math_mean_3D_defaults(script_runner): t1_on_bundle = 't1_on_streamlines.trk' script_runner.run('scil_tractogram_project_map_to_streamlines.py', - in_bundle, in_t1, t1_on_bundle) + in_bundle, t1_on_bundle, + '--in_metric', in_t1, + '--out_dpp_name', 't1') ret = script_runner.run('scil_tractogram_point_math.py', 'mean', t1_on_bundle, - 't1_mean_on_streamlines.trk') + 't1_mean_on_streamlines.trk', + '--in_dpp_name', 't1', + '--out_name', 't1_mean') assert ret.success From 212975bb36d2b02e85d08e6d41a77d1bb3db19c8 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 25 Jan 2024 15:25:00 -0500 Subject: [PATCH 33/97] added test for point math correlation --- scripts/tests/test_tractogram_point_math.py | 25 +++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/scripts/tests/test_tractogram_point_math.py b/scripts/tests/test_tractogram_point_math.py index fa1b871da..fdf049a51 100644 --- a/scripts/tests/test_tractogram_point_math.py +++ b/scripts/tests/test_tractogram_point_math.py @@ -9,6 +9,7 @@ # If they already exist, this only takes 5 seconds (check md5sum) fetch_data(get_testing_files_dict(), keys=['tractometry.zip']) +fetch_data(get_testing_files_dict(), keys=['tracking.zip']) tmp_dir = tempfile.TemporaryDirectory() @@ -39,3 +40,27 @@ def test_execution_tractogram_point_math_mean_3D_defaults(script_runner): '--out_name', 't1_mean') assert ret.success + + +def test_execution_tractogram_point_math_mean_4D_correlation(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + in_bundle = os.path.join(get_home(), 'tracking', + 'local_split_0.trk') + + in_fodf = os.path.join(get_home(), 'tracking', + 'fodf.nii.gz') + fodf_on_bundle = 'fodf_on_streamlines.trk' + + script_runner.run('scil_tractogram_project_map_to_streamlines.py', + in_bundle, fodf_on_bundle, + '--in_metric', in_fodf, in_fodf, + '--out_dpp_name', 'fodf', 'fodf2') + + ret = script_runner.run('scil_tractogram_point_math.py', + 'correlation', + fodf_on_bundle, + 'fodf_correlation_on_streamlines.trk', + '--in_dpp_name', 'fodf', + '--out_name', 'fodf_correlation') + + assert ret.success From bdcf75079ea16353c16217f3bc599577efb82035 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Fri, 26 Jan 2024 09:31:04 -0500 Subject: [PATCH 34/97] pep8 fix --- scilpy/tractograms/streamline_operations.py | 3 ++- scripts/scil_tractogram_point_math.py | 7 +++---- scripts/scil_tractogram_project_map_to_streamlines.py | 1 + scripts/tests/test_tractogram_point_math.py | 2 +- .../tests/test_tractogram_project_map_to_streamlines.py | 1 + 5 files changed, 8 insertions(+), 6 deletions(-) diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index becc00029..750531224 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -567,7 +567,8 @@ def project_metric_to_streamlines(sft, metric, endpoints_only=False): if dimension == 1: thisstreamline_data = np.ones((len(s), 1)) * np.nan else: - thisstreamline_data = np.ones((len(s), p1_data.shape[0])) * np.nan + thisstreamline_data = np.ones( + (len(s), p1_data.shape[0])) * np.nan thisstreamline_data[0] = p1_data thisstreamline_data[-1] = p2_data diff --git a/scripts/scil_tractogram_point_math.py b/scripts/scil_tractogram_point_math.py index a7f1b3abe..1862e3a4c 100755 --- a/scripts/scil_tractogram_point_math.py +++ b/scripts/scil_tractogram_point_math.py @@ -80,7 +80,6 @@ def _build_arg_parser(): 'If not set, will perform operation on all streamline \n' 'points.') - add_reference_arg(p) add_verbose_arg(p) add_overwrite_arg(p) @@ -114,11 +113,11 @@ def main(): data_per_point = {} data_per_streamline = {} for in_dpp_name, out_name in zip(args.in_dpp_name[0], - args.out_name[0]): + args.out_name[0]): # Check to see if the data per point exists. if in_dpp_name not in sft.data_per_point: logging.info('Data per point {} not found in input tractogram.' - .format(in_dpp_name)) + .format(in_dpp_name)) return # check size of first streamline data_per_point @@ -129,7 +128,7 @@ def main(): if args.operation == 'correlation' and is_singular: logging.info('Correlation operation requires multivalued data per ' - 'point. Exiting.') + 'point. Exiting.') return # Perform the requested operation. diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index ef001c9c6..89a113201 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -24,6 +24,7 @@ from scilpy.tractograms.streamline_operations import ( project_metric_to_streamlines) + def _build_arg_parser(): p = argparse.ArgumentParser( description=__doc__, diff --git a/scripts/tests/test_tractogram_point_math.py b/scripts/tests/test_tractogram_point_math.py index fdf049a51..99ab6c966 100644 --- a/scripts/tests/test_tractogram_point_math.py +++ b/scripts/tests/test_tractogram_point_math.py @@ -48,7 +48,7 @@ def test_execution_tractogram_point_math_mean_4D_correlation(script_runner): 'local_split_0.trk') in_fodf = os.path.join(get_home(), 'tracking', - 'fodf.nii.gz') + 'fodf.nii.gz') fodf_on_bundle = 'fodf_on_streamlines.trk' script_runner.run('scil_tractogram_project_map_to_streamlines.py', diff --git a/scripts/tests/test_tractogram_project_map_to_streamlines.py b/scripts/tests/test_tractogram_project_map_to_streamlines.py index 2ed5b9b56..c4b171490 100644 --- a/scripts/tests/test_tractogram_project_map_to_streamlines.py +++ b/scripts/tests/test_tractogram_project_map_to_streamlines.py @@ -72,6 +72,7 @@ def test_execution_4D_map_endpoints_only(script_runner): '--endpoints_only') assert ret.success + def test_execution_3D_map_trilinear(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) in_t1 = os.path.join(get_home(), 'tractometry', 'mni_masked.nii.gz') From 4c2d8448618499bb73e9d6ac8fd4b359838d0f51 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 29 Jan 2024 11:53:47 -0500 Subject: [PATCH 35/97] Updated doc and logging in point_math --- scripts/scil_tractogram_point_math.py | 46 +++++++++++++-------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/scripts/scil_tractogram_point_math.py b/scripts/scil_tractogram_point_math.py index 1862e3a4c..b5e6fa485 100755 --- a/scripts/scil_tractogram_point_math.py +++ b/scripts/scil_tractogram_point_math.py @@ -6,23 +6,24 @@ resulting in data per streamline. The supported operations are: -If the input is singular per point then: +If the chosen dpp key refers to numerical values (one single number per point): mean: mean across points per streamline sum: sum across points per streamline min: min value across points per streamline max: max value across points per streamline -endpoints_only is set, min/max/mean/sum will only +if endpoints_only is set, min/max/mean/sum will only be calculated using the streamline endpoints -If the input is multivalued per point then: +If the chosen dpp key refers to data stored as arrays (multiple numbers per point): mean: mean calculated per point for each streamline sum: sum calculated per point for each streamline min: min value calculated per point for each streamline max: max value calculated per point for each streamline -endpoints_only is set, min/max/mean/sum will only -be calculated at the streamline endpoints +if endpoints_only is set, min/max/mean/sum will only +be calculated at the streamline endpoints the rest of the +values along the streamline will be NaN Endpoint only operations: correlation: correlation calculated between arrays extracted from @@ -41,7 +42,6 @@ add_verbose_arg, assert_inputs_exist, assert_outputs_exist) - from scilpy.tractograms.streamline_operations import ( perform_streamline_operation_on_endpoints, perform_streamline_operation_per_point, @@ -64,13 +64,16 @@ def _build_arg_parser(): 'metadata.') p.add_argument('--in_dpp_name', nargs='+', action='append', required=True, - help='Name of the data_per_point for operation to be ' - 'performed on.') + help='Name or list of names of the data_per_point for ' + 'operation to be performed on. If more than one dpp ' + 'is selected, the same operation will be applied ' + 'separately to each one.') p.add_argument('--out_name', nargs='+', action='append', required=True, help='Name of the resulting data_per_point or ' 'data_per_streamline to be saved in the output ' - 'tractogram. (Default: %(default)s)') + 'tractogram. If more than one --in_dpp_name was used, ' + 'enter the same number of --out_name values.') p.add_argument('out_tractogram', metavar='OUTPUT_FILE', help='The file where the remaining streamlines ' 'are saved.') @@ -110,6 +113,10 @@ def main(): parser.error('The number of in_dpp_names and out_names must be ' 'the same.') + # Check to see if there are duplicates in the out_names + if len(args.out_name[0]) != len(set(args.out_name[0])): + parser.error('The output names (out_names) must be unique.') + data_per_point = {} data_per_streamline = {} for in_dpp_name, out_name in zip(args.in_dpp_name[0], @@ -134,31 +141,28 @@ def main(): # Perform the requested operation. if not is_singular: if args.operation == 'correlation': - logging.info('Performing {} across endpoint data.'.format( - args.operation)) + logging.info('Performing {} across endpoint data and saving as new dpp {}'.format( + args.operation, out_name)) new_dps = perform_streamline_operation_on_endpoints( args.operation, sft, in_dpp_name) - # Adding data per streamline to new_sft data_per_streamline[out_name] = new_dps else: # Results in new data per point logging.info( - 'Performing {} on data from each streamine point.'.format( - args.operation)) + 'Performing {} on data from each streamine point ' + 'and saving as new dpp {}'.format( + args.operation, out_name)) new_dpp = perform_streamline_operation_per_point( args.operation, sft, in_dpp_name, args.endpoints_only) - - # Adding data per point to new_sft data_per_point[out_name] = new_dpp else: # Results in new data per streamline logging.info( - 'Performing {} across each streamline.'.format(args.operation)) + 'Performing {} across each streamline and saving resulting ' + 'data per streamline {}'.format(args.operation, out_name)) new_data_per_streamline = perform_operation_per_streamline( args.operation, sft, in_dpp_name, args.endpoints_only) - - # Adding data per streamline to new_sft data_per_streamline[out_name] = new_data_per_streamline new_sft = StatefulTractogram(sft.streamlines, @@ -167,10 +171,6 @@ def main(): data_per_point=data_per_point, data_per_streamline=data_per_streamline) - if len(new_sft) == 0 and not args.save_empty: - logging.info("Empty resulting tractogram. Not saving results.") - return - # Save the new streamlines (and metadata) logging.info('Saving {} streamlines to {}.'.format(len(new_sft), args.out_tractogram)) From 67873b99dbb210ff2bdcd1971b23c0f196fe9923 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 29 Jan 2024 15:15:48 -0500 Subject: [PATCH 36/97] added dpp and dpp mode argument --- scilpy/tractograms/streamline_operations.py | 11 +- scripts/scil_tractogram_point_math.py | 109 +++++++++++--------- 2 files changed, 65 insertions(+), 55 deletions(-) diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index 750531224..8e053fe08 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -437,7 +437,7 @@ def perform_streamline_operation_per_point(op_name, sft, dpp_name='metric', def perform_operation_per_streamline(op_name, sft, dpp_name='metric', endpoints_only=False): - """Peforms an operation across points for each streamline. + """Performs an operation across all data points for each streamline. Parameters ---------- @@ -461,16 +461,19 @@ def perform_operation_per_streamline(op_name, sft, dpp_name='metric', if endpoints_only: new_data_per_streamline = [] for s in sft.data_per_point[dpp_name]: - new_data_per_streamline.append(call_op([s[0], s[-1]])) + start = s[0] + end = s[-1] + concat = np.concatenate((start[:], end[:]), axis=0) + new_data_per_streamline.append(call_op(concat)) else: new_data_per_streamline = [] for s in sft.data_per_point[dpp_name]: - new_data_per_streamline.append(call_op(s)) + new_data_per_streamline.append(call_op(s[:])) return new_data_per_streamline -def perform_streamline_operation_on_endpoints(op_name, sft, dpp_name='metric'): +def perform_pairwise_streamline_operation_on_endpoints(op_name, sft, dpp_name='metric'): """Peforms an operation across endpoints for each streamline. Parameters diff --git a/scripts/scil_tractogram_point_math.py b/scripts/scil_tractogram_point_math.py index b5e6fa485..4803e6747 100755 --- a/scripts/scil_tractogram_point_math.py +++ b/scripts/scil_tractogram_point_math.py @@ -2,32 +2,29 @@ # -*- coding: utf-8 -*- """ -Performs an operation on data per point from streamlines -resulting in data per streamline. The supported -operations are: - -If the chosen dpp key refers to numerical values (one single number per point): -mean: mean across points per streamline -sum: sum across points per streamline -min: min value across points per streamline -max: max value across points per streamline - -if endpoints_only is set, min/max/mean/sum will only -be calculated using the streamline endpoints - -If the chosen dpp key refers to data stored as arrays (multiple numbers per point): -mean: mean calculated per point for each streamline -sum: sum calculated per point for each streamline -min: min value calculated per point for each streamline -max: max value calculated per point for each streamline - -if endpoints_only is set, min/max/mean/sum will only +Performs an operation on data per point from input streamlines. + +Two modes of operation are supported: data per streamline (dps) and data per +point (dpp). + +In dps mode, the operation is performed across all dimensions of the data +resulting in a single value per streamline. + +In dpp mode the operation is performed on each point separately, resulting in +a single value per point. + +If endpoints_only and dpp mode is set operation will only be calculated at the streamline endpoints the rest of the values along the streamline will be NaN -Endpoint only operations: +If endpoints_only and dps mode is set operation will be calculated +across the data at the endpoints and stored as a +single value per streamline. + +Endpoint only operation: correlation: correlation calculated between arrays extracted from -streamline endpoints (data must be multivalued per point) +streamline endpoints (data must be multivalued per point) and dps +mode must be set. """ import argparse @@ -43,7 +40,7 @@ assert_inputs_exist, assert_outputs_exist) from scilpy.tractograms.streamline_operations import ( - perform_streamline_operation_on_endpoints, + perform_pairwise_streamline_operation_on_endpoints, perform_streamline_operation_per_point, perform_operation_per_streamline) @@ -59,6 +56,13 @@ def _build_arg_parser(): help='The type of operation to be performed on the ' 'streamlines. Must\nbe one of the following: ' '%(choices)s.') + p.add_argument('dpp_or_dps', metavar='DPP_OR_DPS', + choices=['dpp', 'dps'], + help='Set to dps if the operation is to be performed ' + 'across all dimensions resulting in a single value per ' + 'streamline. Set to dpp if the operation is to be ' + 'performed on each point separately resulting in a single ' + 'value per point.') p.add_argument('in_tractogram', metavar='INPUT_FILE', help='Input tractogram containing streamlines and ' 'metadata.') @@ -117,46 +121,49 @@ def main(): if len(args.out_name[0]) != len(set(args.out_name[0])): parser.error('The output names (out_names) must be unique.') - data_per_point = {} - data_per_streamline = {} - for in_dpp_name, out_name in zip(args.in_dpp_name[0], - args.out_name[0]): + # Input name checks + for in_dpp_name in args.in_dpp_name[0]: # Check to see if the data per point exists. if in_dpp_name not in sft.data_per_point: logging.info('Data per point {} not found in input tractogram.' .format(in_dpp_name)) return - # check size of first streamline data_per_point - if sft.data_per_point[in_dpp_name][0].shape[1] == 1: - is_singular = True - else: - is_singular = False - - if args.operation == 'correlation' and is_singular: + data_shape = sft.data_per_point[in_dpp_name][0].shape + if args.operation == 'correlation' and len(data_shape) == 1: logging.info('Correlation operation requires multivalued data per ' 'point. Exiting.') return + if args.operation == 'correlation' and args.dpp_or_dps == 'dpp': + logging.info('Correlation operation requires dps mode. Exiting.') + return + + data_per_point = {} + data_per_streamline = {} + for in_dpp_name, out_name in zip(args.in_dpp_name[0], + args.out_name[0]): + + # Perform the requested operation. - if not is_singular: - if args.operation == 'correlation': - logging.info('Performing {} across endpoint data and saving as new dpp {}'.format( + if args.operation == 'correlation': + logging.info('Performing {} across endpoint data and saving as ' + 'new dpp {}'.format( + args.operation, out_name)) + new_dps = perform_pairwise_streamline_operation_on_endpoints( + args.operation, sft, in_dpp_name) + + data_per_streamline[out_name] = new_dps + elif args.dpp_or_dps == 'dpp': + # Results in new data per point + logging.info( + 'Performing {} on data from each streamine point ' + 'and saving as new dpp {}'.format( args.operation, out_name)) - new_dps = perform_streamline_operation_on_endpoints( - args.operation, sft, in_dpp_name) - - data_per_streamline[out_name] = new_dps - else: - # Results in new data per point - logging.info( - 'Performing {} on data from each streamine point ' - 'and saving as new dpp {}'.format( - args.operation, out_name)) - new_dpp = perform_streamline_operation_per_point( - args.operation, sft, in_dpp_name, args.endpoints_only) - data_per_point[out_name] = new_dpp - else: + new_dpp = perform_streamline_operation_per_point( + args.operation, sft, in_dpp_name, args.endpoints_only) + data_per_point[out_name] = new_dpp + elif args.dpp_or_dps == 'dps': # Results in new data per streamline logging.info( 'Performing {} across each streamline and saving resulting ' From 154f8ed0ea41f26d68f25cb36dc1ab6a3e8cb2f8 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 29 Jan 2024 16:17:35 -0500 Subject: [PATCH 37/97] added overwrite functionality --- scripts/scil_tractogram_point_math.py | 53 +++++++++++++++++++++++---- 1 file changed, 46 insertions(+), 7 deletions(-) diff --git a/scripts/scil_tractogram_point_math.py b/scripts/scil_tractogram_point_math.py index 4803e6747..85d7fe4e1 100755 --- a/scripts/scil_tractogram_point_math.py +++ b/scripts/scil_tractogram_point_math.py @@ -86,7 +86,10 @@ def _build_arg_parser(): help='If set, will only perform operation on endpoints \n' 'If not set, will perform operation on all streamline \n' 'points.') - + p.add_argument('--overwrite_data', action='store_true', default=False, + help='If set, will overwrite the data_per_point or ' + 'data_per_streamline in the output tractogram, otherwise ' + 'previous data will be preserved in the output tractogram.') add_reference_arg(p) add_verbose_arg(p) add_overwrite_arg(p) @@ -139,6 +142,13 @@ def main(): logging.info('Correlation operation requires dps mode. Exiting.') return + if args.overwrite_data: + if in_dpp_name in args.out_name[0]: + logging.info('out_name {} already exists in input tractogram. ' + 'Unset overwrite_data or choose a different ' + 'out_name. Exiting.'.format(in_dpp_name)) + return + data_per_point = {} data_per_streamline = {} for in_dpp_name, out_name in zip(args.in_dpp_name[0], @@ -172,11 +182,41 @@ def main(): args.operation, sft, in_dpp_name, args.endpoints_only) data_per_streamline[out_name] = new_data_per_streamline - new_sft = StatefulTractogram(sft.streamlines, - sft.space_attributes, - sft.space, sft.origin, - data_per_point=data_per_point, - data_per_streamline=data_per_streamline) + + if args.overwrite_data: + new_sft = sft.from_sft(sft.streamlines, sft, + data_per_point=data_per_point, + data_per_streamline=data_per_streamline) + else: + old_data_per_streamline = {} + for key, value in sft.data_per_streamline.items(): + old_data_per_streamline[key] = value + + old_data_per_point = {} + for key, value in sft.data_per_point.items(): + old_data_per_point[key] = value + + if data_per_point is not None: + for key, value in data_per_point.items(): + old_data_per_point[key] = value + + if data_per_streamline is not None: + for key, value in data_per_streamline.items(): + old_data_per_streamline[key] = value + + new_sft = sft.from_sft(sft.streamlines, sft, + data_per_point=old_data_per_point, + data_per_streamline=old_data_per_streamline) + + # Print DPP names + logging.info('New data per point names:') + for key in new_sft.data_per_point.keys(): + logging.info(key) + + # Print DPS names + logging.info('New data per streamline names:') + for key in new_sft.data_per_streamline.keys(): + logging.info(key) # Save the new streamlines (and metadata) logging.info('Saving {} streamlines to {}.'.format(len(new_sft), @@ -184,6 +224,5 @@ def main(): save_tractogram(new_sft, args.out_tractogram, bbox_valid_check=args.bbox_check) - if __name__ == "__main__": main() From e7b5f0be81c177ed6d26107db984a277f18923de Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 29 Jan 2024 20:14:18 -0500 Subject: [PATCH 38/97] simplify overwrite --- scripts/scil_tractogram_point_math.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/scripts/scil_tractogram_point_math.py b/scripts/scil_tractogram_point_math.py index 85d7fe4e1..eaf429e94 100755 --- a/scripts/scil_tractogram_point_math.py +++ b/scripts/scil_tractogram_point_math.py @@ -188,13 +188,8 @@ def main(): data_per_point=data_per_point, data_per_streamline=data_per_streamline) else: - old_data_per_streamline = {} - for key, value in sft.data_per_streamline.items(): - old_data_per_streamline[key] = value - - old_data_per_point = {} - for key, value in sft.data_per_point.items(): - old_data_per_point[key] = value + old_data_per_streamline = sft.data_per_streamline + old_data_per_point = sft.data_per_point if data_per_point is not None: for key, value in data_per_point.items(): From 39f65ef17d7d26108d22aa111053d9ae263e89d4 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 29 Jan 2024 20:17:04 -0500 Subject: [PATCH 39/97] remove append option --- scripts/scil_tractogram_point_math.py | 4 ++-- scripts/scil_tractogram_project_map_to_streamlines.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/scripts/scil_tractogram_point_math.py b/scripts/scil_tractogram_point_math.py index eaf429e94..e6224fc1e 100755 --- a/scripts/scil_tractogram_point_math.py +++ b/scripts/scil_tractogram_point_math.py @@ -66,13 +66,13 @@ def _build_arg_parser(): p.add_argument('in_tractogram', metavar='INPUT_FILE', help='Input tractogram containing streamlines and ' 'metadata.') - p.add_argument('--in_dpp_name', nargs='+', action='append', + p.add_argument('--in_dpp_name', nargs='+', required=True, help='Name or list of names of the data_per_point for ' 'operation to be performed on. If more than one dpp ' 'is selected, the same operation will be applied ' 'separately to each one.') - p.add_argument('--out_name', nargs='+', action='append', + p.add_argument('--out_name', nargs='+', required=True, help='Name of the resulting data_per_point or ' 'data_per_streamline to be saved in the output ' diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index 89a113201..a7a76ac20 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -35,9 +35,9 @@ def _build_arg_parser(): help='Fiber bundle file.') p.add_argument('out_tractogram', help='Output file.') - p.add_argument('--in_metric', nargs='+', action='append', required=True, + p.add_argument('--in_metric', nargs='+', required=True, help='Nifti metric to project onto streamlines.') - p.add_argument('--out_dpp_name', nargs='+', action='append', required=True, + p.add_argument('--out_dpp_name', nargs='+', required=True, help='Name of the data_per_point to be saved in the \n' 'output tractogram.') From 6b33f7c325fd0c2340c3a30aef887340edbccb7d Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Tue, 30 Jan 2024 13:22:52 -0500 Subject: [PATCH 40/97] arg multiinput just narg+ --- scripts/scil_tractogram_point_math.py | 24 +++++----- ...l_tractogram_project_map_to_streamlines.py | 47 ++++++++++++++----- 2 files changed, 46 insertions(+), 25 deletions(-) diff --git a/scripts/scil_tractogram_point_math.py b/scripts/scil_tractogram_point_math.py index e6224fc1e..d058bb4ea 100755 --- a/scripts/scil_tractogram_point_math.py +++ b/scripts/scil_tractogram_point_math.py @@ -66,14 +66,12 @@ def _build_arg_parser(): p.add_argument('in_tractogram', metavar='INPUT_FILE', help='Input tractogram containing streamlines and ' 'metadata.') - p.add_argument('--in_dpp_name', nargs='+', - required=True, + p.add_argument('--in_dpp_name', nargs='+', required=True, help='Name or list of names of the data_per_point for ' 'operation to be performed on. If more than one dpp ' 'is selected, the same operation will be applied ' 'separately to each one.') - p.add_argument('--out_name', nargs='+', - required=True, + p.add_argument('--out_name', nargs='+', required=True, help='Name of the resulting data_per_point or ' 'data_per_streamline to be saved in the output ' 'tractogram. If more than one --in_dpp_name was used, ' @@ -90,6 +88,7 @@ def _build_arg_parser(): help='If set, will overwrite the data_per_point or ' 'data_per_streamline in the output tractogram, otherwise ' 'previous data will be preserved in the output tractogram.') + add_reference_arg(p) add_verbose_arg(p) add_overwrite_arg(p) @@ -116,22 +115,23 @@ def main(): logging.info("Input tractogram contains no streamlines. Exiting.") return - if len(args.in_dpp_name[0]) != len(args.out_name[0]): + if len(args.in_dpp_name) != len(args.out_name): parser.error('The number of in_dpp_names and out_names must be ' 'the same.') # Check to see if there are duplicates in the out_names - if len(args.out_name[0]) != len(set(args.out_name[0])): + if len(args.out_name) != len(set(args.out_name)): parser.error('The output names (out_names) must be unique.') # Input name checks - for in_dpp_name in args.in_dpp_name[0]: + for in_dpp_name in args.in_dpp_name: # Check to see if the data per point exists. if in_dpp_name not in sft.data_per_point: logging.info('Data per point {} not found in input tractogram.' .format(in_dpp_name)) return + # Check if first data_per_point is multivalued data_shape = sft.data_per_point[in_dpp_name][0].shape if args.operation == 'correlation' and len(data_shape) == 1: logging.info('Correlation operation requires multivalued data per ' @@ -142,17 +142,17 @@ def main(): logging.info('Correlation operation requires dps mode. Exiting.') return - if args.overwrite_data: - if in_dpp_name in args.out_name[0]: + if not args.overwrite_data: + if in_dpp_name in args.out_name: logging.info('out_name {} already exists in input tractogram. ' - 'Unset overwrite_data or choose a different ' + 'Set overwrite_data or choose a different ' 'out_name. Exiting.'.format(in_dpp_name)) return data_per_point = {} data_per_streamline = {} - for in_dpp_name, out_name in zip(args.in_dpp_name[0], - args.out_name[0]): + for in_dpp_name, out_name in zip(args.in_dpp_name, + args.out_name): # Perform the requested operation. diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index a7a76ac20..c90cc9456 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -5,7 +5,7 @@ Projects metrics extracted from a map onto the endpoints of streamlines. The default options will take data from a nifti image (3D or ND) and -project it onto the endpoints of streamlines. +project it onto the points of streamlines. """ import argparse @@ -35,7 +35,7 @@ def _build_arg_parser(): help='Fiber bundle file.') p.add_argument('out_tractogram', help='Output file.') - p.add_argument('--in_metric', nargs='+', required=True, + p.add_argument('--in_maps', nargs='+', required=True, help='Nifti metric to project onto streamlines.') p.add_argument('--out_dpp_name', nargs='+', required=True, help='Name of the data_per_point to be saved in the \n' @@ -46,12 +46,15 @@ def _build_arg_parser(): help='If set, will use trilinear interpolation \n' 'else will use nearest neighbor interpolation \n' 'by default.') - p.add_argument('--endpoints_only', action='store_true', help='If set, will only project the metric onto the \n' 'endpoints of the streamlines (all other values along \n' ' streamlines will be NaN). If not set, will project \n' ' the metric onto all points of the streamlines.') + p.add_argument('--overwrite_data', action='store_true', default=False, + help='If set, will overwrite data_per_point in the ' + 'output tractogram, otherwise previous data will be ' + 'preserved in the output tractogram.') add_reference_arg(p) add_overwrite_arg(p) @@ -63,10 +66,7 @@ def main(): parser = _build_arg_parser() args = parser.parse_args() - input_file_list = [args.in_tractogram] - for fmetric in args.in_metric: - input_file_list += fmetric - assert_inputs_exist(parser, input_file_list) + assert_inputs_exist(parser, [args.in_tractogram] + args.in_maps) assert_outputs_exist(parser, args, [args.out_tractogram]) @@ -76,15 +76,28 @@ def main(): sft.to_corner() if len(sft.streamlines) == 0: - logging.warning('Empty bundle file {}. Skipping'.format(args.bundle)) + logging.warning('Empty bundle file {}. Skipping'.format(args.in_tractogram)) return # Check to see if the number of metrics and dpp_names are the same - if len(args.in_metric) != len(args.out_dpp_name): + if len(args.in_maps) != len(args.out_dpp_name): parser.error('The number of metrics and dpp_names must be the same.') + # Check to see if there are duplicates in the out_dpp_names + if len(args.out_dpp_name) != len(set(args.out_dpp_name)): + parser.error('The output names (out_dpp_names) must be unique.') + + # Check to see if the output names already exist in the input tractogram + if not args.overwrite_data: + for out_dpp_name in args.out_dpp_name: + if out_dpp_name in sft.data_per_point: + logging.info('out_name {} already exists in input tractogram. ' + 'Set overwrite_data or choose a different ' + 'out_name. Exiting.'.format(out_dpp_name)) + return + data_per_point = {} - for fmetric, dpp_name in zip(args.in_metric[0], args.out_dpp_name[0]): + for fmetric, dpp_name in zip(args.in_maps, args.out_dpp_name): logging.debug("Loading the metric...") metric_img, metric_dtype = load_img(fmetric) metric_data = metric_img.get_fdata(caching='unchanged', dtype=float) @@ -106,9 +119,17 @@ def main(): data_per_point[dpp_name] = streamline_data - out_sft = StatefulTractogram(sft.streamlines, metric_img, - sft.space, sft.origin, - data_per_point=data_per_point) + if args.overwrite_data: + out_sft = sft.from_sft(sft.streamlines, sft, + data_per_point=data_per_point) + else: + old_data_per_point = sft.data_per_point + for dpp_name in data_per_point: + old_data_per_point[dpp_name] = data_per_point[dpp_name] + out_sft = sft.from_sft(sft.streamlines, sft, + data_per_point=old_data_per_point) + + print(out_sft.data_per_point) save_tractogram(out_sft, args.out_tractogram) From e679fa13272d484ae3dc1551620372e7102d4ec9 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Tue, 30 Jan 2024 14:48:12 -0500 Subject: [PATCH 41/97] fixed project map tests --- .../test_tractogram_project_map_to_streamlines.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/scripts/tests/test_tractogram_project_map_to_streamlines.py b/scripts/tests/test_tractogram_project_map_to_streamlines.py index c4b171490..ce3b125fe 100644 --- a/scripts/tests/test_tractogram_project_map_to_streamlines.py +++ b/scripts/tests/test_tractogram_project_map_to_streamlines.py @@ -25,7 +25,7 @@ def test_execution_3D_map(script_runner): ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', in_tracto_1, 't1_on_streamlines.trk', - '--in_metric', in_t1, + '--in_maps', in_t1, '--out_dpp_name', 't1') assert ret.success @@ -38,7 +38,7 @@ def test_execution_4D_map(script_runner): ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', in_tracto_1, 'rgb_on_streamlines.trk', - '--in_metric', in_rgb, + '--in_maps', in_rgb, '--out_dpp_name', 'rgb') assert ret.success @@ -52,7 +52,7 @@ def test_execution_3D_map_endpoints_only(script_runner): ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', in_tracto_1, 't1_on_streamlines_endpoints.trk', - '--in_metric', in_t1, + '--in_maps', in_t1, '--out_dpp_name', 't1', '--endpoints_only') assert ret.success @@ -67,7 +67,7 @@ def test_execution_4D_map_endpoints_only(script_runner): ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', in_tracto_1, 'rgb_on_streamlines_endpoints.trk', - '--in_metric', in_rgb, + '--in_maps', in_rgb, '--out_dpp_name', 'rgb', '--endpoints_only') assert ret.success @@ -82,7 +82,7 @@ def test_execution_3D_map_trilinear(script_runner): ret = script_runner.run('scil_tractogram_project_map_to_streamlines.py', in_tracto_1, 't1_on_streamlines_trilinear.trk', - '--in_metric', in_t1, + '--in_maps', in_t1, '--out_dpp_name', 't1', '--trilinear') assert ret.success From 4046cc11b5f7fcb5cd80fc5f8649aa46450258ae Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Tue, 30 Jan 2024 14:49:50 -0500 Subject: [PATCH 42/97] rm doc reference to endpoints --- scripts/scil_tractogram_project_map_to_streamlines.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index c90cc9456..f702dd2d0 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -2,7 +2,7 @@ # -*- coding: utf-8 -*- """ -Projects metrics extracted from a map onto the endpoints of streamlines. +Projects metrics extracted from a map onto the points of streamlines. The default options will take data from a nifti image (3D or ND) and project it onto the points of streamlines. From 0c5cdd3fa2319c7d906918ba61de72c33992e12a Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Tue, 30 Jan 2024 14:53:16 -0500 Subject: [PATCH 43/97] update test tractogram point math new args --- scripts/tests/test_tractogram_point_math.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/scripts/tests/test_tractogram_point_math.py b/scripts/tests/test_tractogram_point_math.py index 99ab6c966..f819f2523 100644 --- a/scripts/tests/test_tractogram_point_math.py +++ b/scripts/tests/test_tractogram_point_math.py @@ -29,11 +29,12 @@ def test_execution_tractogram_point_math_mean_3D_defaults(script_runner): script_runner.run('scil_tractogram_project_map_to_streamlines.py', in_bundle, t1_on_bundle, - '--in_metric', in_t1, + '--in_maps', in_t1, '--out_dpp_name', 't1') ret = script_runner.run('scil_tractogram_point_math.py', 'mean', + 'dps', t1_on_bundle, 't1_mean_on_streamlines.trk', '--in_dpp_name', 't1', @@ -53,11 +54,12 @@ def test_execution_tractogram_point_math_mean_4D_correlation(script_runner): script_runner.run('scil_tractogram_project_map_to_streamlines.py', in_bundle, fodf_on_bundle, - '--in_metric', in_fodf, in_fodf, + '--in_maps', in_fodf, in_fodf, '--out_dpp_name', 'fodf', 'fodf2') ret = script_runner.run('scil_tractogram_point_math.py', 'correlation', + 'dps', fodf_on_bundle, 'fodf_correlation_on_streamlines.trk', '--in_dpp_name', 'fodf', From 0963198c2f3eb2cc5adc60a4adc4d58681536bb6 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Tue, 30 Jan 2024 14:57:33 -0500 Subject: [PATCH 44/97] updated name of to dpp_math --- ...tractogram_point_math.py => scil_tractogram_dpp_math.py} | 0 ...tractogram_point_math.py => test_tractogram_dpp_math.py} | 6 +++--- 2 files changed, 3 insertions(+), 3 deletions(-) rename scripts/{scil_tractogram_point_math.py => scil_tractogram_dpp_math.py} (100%) rename scripts/tests/{test_tractogram_point_math.py => test_tractogram_dpp_math.py} (92%) diff --git a/scripts/scil_tractogram_point_math.py b/scripts/scil_tractogram_dpp_math.py similarity index 100% rename from scripts/scil_tractogram_point_math.py rename to scripts/scil_tractogram_dpp_math.py diff --git a/scripts/tests/test_tractogram_point_math.py b/scripts/tests/test_tractogram_dpp_math.py similarity index 92% rename from scripts/tests/test_tractogram_point_math.py rename to scripts/tests/test_tractogram_dpp_math.py index f819f2523..ef89d6eca 100644 --- a/scripts/tests/test_tractogram_point_math.py +++ b/scripts/tests/test_tractogram_dpp_math.py @@ -14,7 +14,7 @@ def test_help_option(script_runner): - ret = script_runner.run('scil_tractogram_point_math.py', '--help') + ret = script_runner.run('scil_tractogram_dpp_math.py', '--help') assert ret.success @@ -32,7 +32,7 @@ def test_execution_tractogram_point_math_mean_3D_defaults(script_runner): '--in_maps', in_t1, '--out_dpp_name', 't1') - ret = script_runner.run('scil_tractogram_point_math.py', + ret = script_runner.run('scil_tractogram_dpp_math.py', 'mean', 'dps', t1_on_bundle, @@ -57,7 +57,7 @@ def test_execution_tractogram_point_math_mean_4D_correlation(script_runner): '--in_maps', in_fodf, in_fodf, '--out_dpp_name', 'fodf', 'fodf2') - ret = script_runner.run('scil_tractogram_point_math.py', + ret = script_runner.run('scil_tractogram_dpp_math.py', 'correlation', 'dps', fodf_on_bundle, From ff8c3779e3614d378eacf4e4043ed452518e2c8c Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Tue, 30 Jan 2024 15:12:42 -0500 Subject: [PATCH 45/97] update doc for clearer on what happens in 4D case --- scripts/scil_tractogram_project_map_to_streamlines.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index f702dd2d0..1be51496d 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -4,8 +4,10 @@ """ Projects metrics extracted from a map onto the points of streamlines. -The default options will take data from a nifti image (3D or ND) and -project it onto the points of streamlines. +The default options will take data from a nifti image (3D or 4D) and +project it onto the points of streamlines. If the image is 4D, the data +is stored as a list of 1D arrays per streamline. If the image is 3D, the data is stored +as a list of values per streamline. """ import argparse From d3fd65ac0c6361fda0697f64cbc4d418e598251b Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Tue, 30 Jan 2024 15:44:28 -0500 Subject: [PATCH 46/97] remove print --- scripts/scil_tractogram_project_map_to_streamlines.py | 1 - 1 file changed, 1 deletion(-) diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index 1be51496d..eb37f2c02 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -131,7 +131,6 @@ def main(): out_sft = sft.from_sft(sft.streamlines, sft, data_per_point=old_data_per_point) - print(out_sft.data_per_point) save_tractogram(out_sft, args.out_tractogram) From 0778e9613a899addfcf1101b2f3bed1566f512c5 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 5 Feb 2024 11:00:55 -0500 Subject: [PATCH 47/97] moved project func to new file --- scilpy/tractograms/dps_and_dpp_management.py | 65 +++++++++++++++++++ scilpy/tractograms/streamline_operations.py | 65 ------------------- ...l_tractogram_project_map_to_streamlines.py | 2 +- 3 files changed, 66 insertions(+), 66 deletions(-) create mode 100644 scilpy/tractograms/dps_and_dpp_management.py diff --git a/scilpy/tractograms/dps_and_dpp_management.py b/scilpy/tractograms/dps_and_dpp_management.py new file mode 100644 index 000000000..a365b1a79 --- /dev/null +++ b/scilpy/tractograms/dps_and_dpp_management.py @@ -0,0 +1,65 @@ +import numpy as np + +def project_metric_to_streamlines(sft, metric, endpoints_only=False): + """ + Projects a metric onto the points of streamlines. + + Parameters + ---------- + sft: StatefulTractogram + Input tractogram. + metric: DataVolume + Input metric. + + Optional: + --------- + endpoints_only: bool + If True, will only project the metric onto the endpoints of the + streamlines (all values along streamlines set to zero). If False, + will project the metric onto all points of the streamlines. + + Returns + ------- + streamline_data: + metric projected to each point of the streamlines. + """ + if len(metric.data.shape) == 4: + dimension = metric.data.shape[3] + else: + dimension = 1 + + streamline_data = [] + if endpoints_only: + for s in sft.streamlines: + p1_data = metric.get_value_at_coordinate( + s[0][0], s[0][1], s[0][2], + space=sft.space, origin=sft.origin) + p2_data = metric.get_value_at_coordinate( + s[-1][0], s[-1][1], s[-1][2], + space=sft.space, origin=sft.origin) + thisstreamline_data = [] + if dimension == 1: + thisstreamline_data = np.ones((len(s), 1)) * np.nan + else: + thisstreamline_data = np.ones( + (len(s), p1_data.shape[0])) * np.nan + + thisstreamline_data[0] = p1_data + thisstreamline_data[-1] = p2_data + thisstreamline_data = np.asarray(thisstreamline_data) + + streamline_data.append( + np.reshape(thisstreamline_data, + (len(thisstreamline_data), dimension))) + else: + for s in sft.streamlines: + thisstreamline_data = [] + for p in s: + thisstreamline_data.append(metric.get_value_at_coordinate( + p[0], p[1], p[2], space=sft.space, origin=sft.origin)) + + streamline_data.append( + np.reshape(thisstreamline_data, + (len(thisstreamline_data), dimension))) + + return streamline_data diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index 8e053fe08..2d0bdf2c7 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -527,68 +527,3 @@ def stream_correlation(array1, array2): 'max': stream_max, 'correlation': stream_correlation, } - - -def project_metric_to_streamlines(sft, metric, endpoints_only=False): - """ - Projects a metric onto the points of streamlines. - - Parameters - ---------- - sft: StatefulTractogram - Input tractogram. - metric: DataVolume - Input metric. - - Optional: - --------- - endpoints_only: bool - If True, will only project the metric onto the endpoints of the - streamlines (all values along streamlines set to zero). If False, - will project the metric onto all points of the streamlines. - - Returns - ------- - streamline_data: - metric projected to each point of the streamlines. - """ - if len(metric.data.shape) == 4: - dimension = metric.data.shape[3] - else: - dimension = 1 - - streamline_data = [] - if endpoints_only: - for s in sft.streamlines: - p1_data = metric.get_value_at_coordinate( - s[0][0], s[0][1], s[0][2], - space=sft.space, origin=sft.origin) - p2_data = metric.get_value_at_coordinate( - s[-1][0], s[-1][1], s[-1][2], - space=sft.space, origin=sft.origin) - thisstreamline_data = [] - if dimension == 1: - thisstreamline_data = np.ones((len(s), 1)) * np.nan - else: - thisstreamline_data = np.ones( - (len(s), p1_data.shape[0])) * np.nan - - thisstreamline_data[0] = p1_data - thisstreamline_data[-1] = p2_data - thisstreamline_data = np.asarray(thisstreamline_data) - - streamline_data.append( - np.reshape(thisstreamline_data, - (len(thisstreamline_data), dimension))) - else: - for s in sft.streamlines: - thisstreamline_data = [] - for p in s: - thisstreamline_data.append(metric.get_value_at_coordinate( - p[0], p[1], p[2], space=sft.space, origin=sft.origin)) - - streamline_data.append( - np.reshape(thisstreamline_data, - (len(thisstreamline_data), dimension))) - - return streamline_data diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index eb37f2c02..c6578bbc9 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -23,7 +23,7 @@ assert_inputs_exist, assert_outputs_exist) from scilpy.image.volume_space_management import DataVolume -from scilpy.tractograms.streamline_operations import ( +from scilpy.tractograms.dps_and_dpp_management import ( project_metric_to_streamlines) From 918d4c279f74937cecdc7fd8359aa63f46c3cfc1 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 5 Feb 2024 11:03:00 -0500 Subject: [PATCH 48/97] changed metric to map --- scilpy/tractograms/dps_and_dpp_management.py | 24 ++++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/scilpy/tractograms/dps_and_dpp_management.py b/scilpy/tractograms/dps_and_dpp_management.py index a365b1a79..cf79189be 100644 --- a/scilpy/tractograms/dps_and_dpp_management.py +++ b/scilpy/tractograms/dps_and_dpp_management.py @@ -1,40 +1,40 @@ import numpy as np -def project_metric_to_streamlines(sft, metric, endpoints_only=False): +def project_map_to_streamlines(sft, map, endpoints_only=False): """ - Projects a metric onto the points of streamlines. + Projects a map onto the points of streamlines. Parameters ---------- sft: StatefulTractogram Input tractogram. - metric: DataVolume - Input metric. + map: DataVolume + Input map. Optional: --------- endpoints_only: bool - If True, will only project the metric onto the endpoints of the + If True, will only project the map onto the endpoints of the streamlines (all values along streamlines set to zero). If False, - will project the metric onto all points of the streamlines. + will project the map onto all points of the streamlines. Returns ------- streamline_data: - metric projected to each point of the streamlines. + map projected to each point of the streamlines. """ - if len(metric.data.shape) == 4: - dimension = metric.data.shape[3] + if len(map.data.shape) == 4: + dimension = map.data.shape[3] else: dimension = 1 streamline_data = [] if endpoints_only: for s in sft.streamlines: - p1_data = metric.get_value_at_coordinate( + p1_data = map.get_value_at_coordinate( s[0][0], s[0][1], s[0][2], space=sft.space, origin=sft.origin) - p2_data = metric.get_value_at_coordinate( + p2_data = map.get_value_at_coordinate( s[-1][0], s[-1][1], s[-1][2], space=sft.space, origin=sft.origin) thisstreamline_data = [] @@ -55,7 +55,7 @@ def project_metric_to_streamlines(sft, metric, endpoints_only=False): for s in sft.streamlines: thisstreamline_data = [] for p in s: - thisstreamline_data.append(metric.get_value_at_coordinate( + thisstreamline_data.append(map.get_value_at_coordinate( p[0], p[1], p[2], space=sft.space, origin=sft.origin)) streamline_data.append( From ce94ab07dd428c2f60ad158c8a3a32680d1fdaf6 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 5 Feb 2024 11:04:18 -0500 Subject: [PATCH 49/97] changed metric to map --- ...l_tractogram_project_map_to_streamlines.py | 32 +++++++++---------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index c6578bbc9..d31e31d35 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -2,7 +2,7 @@ # -*- coding: utf-8 -*- """ -Projects metrics extracted from a map onto the points of streamlines. +Projects maps extracted from a map onto the points of streamlines. The default options will take data from a nifti image (3D or 4D) and project it onto the points of streamlines. If the image is 4D, the data @@ -24,7 +24,7 @@ assert_outputs_exist) from scilpy.image.volume_space_management import DataVolume from scilpy.tractograms.dps_and_dpp_management import ( - project_metric_to_streamlines) + project_map_to_streamlines) def _build_arg_parser(): @@ -38,7 +38,7 @@ def _build_arg_parser(): p.add_argument('out_tractogram', help='Output file.') p.add_argument('--in_maps', nargs='+', required=True, - help='Nifti metric to project onto streamlines.') + help='Nifti map to project onto streamlines.') p.add_argument('--out_dpp_name', nargs='+', required=True, help='Name of the data_per_point to be saved in the \n' 'output tractogram.') @@ -49,10 +49,10 @@ def _build_arg_parser(): 'else will use nearest neighbor interpolation \n' 'by default.') p.add_argument('--endpoints_only', action='store_true', - help='If set, will only project the metric onto the \n' + help='If set, will only project the map onto the \n' 'endpoints of the streamlines (all other values along \n' ' streamlines will be NaN). If not set, will project \n' - ' the metric onto all points of the streamlines.') + ' the map onto all points of the streamlines.') p.add_argument('--overwrite_data', action='store_true', default=False, help='If set, will overwrite data_per_point in the ' 'output tractogram, otherwise previous data will be ' @@ -81,9 +81,9 @@ def main(): logging.warning('Empty bundle file {}. Skipping'.format(args.in_tractogram)) return - # Check to see if the number of metrics and dpp_names are the same + # Check to see if the number of maps and dpp_names are the same if len(args.in_maps) != len(args.out_dpp_name): - parser.error('The number of metrics and dpp_names must be the same.') + parser.error('The number of maps and dpp_names must be the same.') # Check to see if there are duplicates in the out_dpp_names if len(args.out_dpp_name) != len(set(args.out_dpp_name)): @@ -99,22 +99,22 @@ def main(): return data_per_point = {} - for fmetric, dpp_name in zip(args.in_maps, args.out_dpp_name): - logging.debug("Loading the metric...") - metric_img, metric_dtype = load_img(fmetric) - metric_data = metric_img.get_fdata(caching='unchanged', dtype=float) - metric_res = metric_img.header.get_zooms()[:3] + for fmap, dpp_name in zip(args.in_maps, args.out_dpp_name): + logging.debug("Loading the map...") + map_img, map_dtype = load_img(fmap) + map_data = map_img.get_fdata(caching='unchanged', dtype=float) + map_res = map_img.header.get_zooms()[:3] if args.trilinear: interp = "trilinear" else: interp = "nearest" - metric = DataVolume(metric_data, metric_res, interp) + map = DataVolume(map_data, map_res, interp) - logging.debug("Projecting metric onto streamlines") - streamline_data = project_metric_to_streamlines( - sft, metric, + logging.debug("Projecting map onto streamlines") + streamline_data = project_map_to_streamlines( + sft, map, endpoints_only=args.endpoints_only) logging.debug("Saving the tractogram...") From 75a54bff072bd88551ad345fd7de4cdb714f807a Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 5 Feb 2024 11:08:11 -0500 Subject: [PATCH 50/97] added dummy test for dpp management --- scilpy/tractograms/tests/test_dps_and_dpp_management.py | 3 +++ scilpy/tractograms/tests/test_streamline_operations.py | 5 ----- 2 files changed, 3 insertions(+), 5 deletions(-) create mode 100644 scilpy/tractograms/tests/test_dps_and_dpp_management.py diff --git a/scilpy/tractograms/tests/test_dps_and_dpp_management.py b/scilpy/tractograms/tests/test_dps_and_dpp_management.py new file mode 100644 index 000000000..fb93d736a --- /dev/null +++ b/scilpy/tractograms/tests/test_dps_and_dpp_management.py @@ -0,0 +1,3 @@ +def test_project_map_to_streamlines(): + # toDo + pass diff --git a/scilpy/tractograms/tests/test_streamline_operations.py b/scilpy/tractograms/tests/test_streamline_operations.py index b37c4aec7..f8e0d4e8a 100644 --- a/scilpy/tractograms/tests/test_streamline_operations.py +++ b/scilpy/tractograms/tests/test_streamline_operations.py @@ -369,11 +369,6 @@ def test_smooth_line_spline(): assert dist_1 < dist_2 -def test_project_metric_to_streamlines(): - # toDo - pass - - def test_perform_streamline_operation_per_point(): # toDo pass From 74d6a8714487119f545aa11631944fda71c87e98 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 5 Feb 2024 11:59:43 -0500 Subject: [PATCH 51/97] 4D dps math - now outputs array per streamline --- scilpy/tractograms/streamline_operations.py | 14 +++++++------- .../scil_tractogram_project_map_to_streamlines.py | 2 +- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index 2d0bdf2c7..8324d0e94 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -463,12 +463,13 @@ def perform_operation_per_streamline(op_name, sft, dpp_name='metric', for s in sft.data_per_point[dpp_name]: start = s[0] end = s[-1] - concat = np.concatenate((start[:], end[:]), axis=0) + concat = np.concatenate((start[:], end[:])) new_data_per_streamline.append(call_op(concat)) else: new_data_per_streamline = [] for s in sft.data_per_point[dpp_name]: - new_data_per_streamline.append(call_op(s[:])) + s_np = np.asarray(s) + new_data_per_streamline.append(call_op(s_np)) return new_data_per_streamline @@ -501,19 +502,18 @@ def perform_pairwise_streamline_operation_on_endpoints(op_name, sft, dpp_name='m def stream_mean(array): - return np.mean(array) + return np.squeeze(np.mean(array,axis=0)) def stream_sum(array): - return np.sum(array) + return np.squeeze(np.sum(array,axis=0)) def stream_min(array): - return np.min(array) - + return np.squeeze(np.min(array,axis=0)) def stream_max(array): - return np.max(array) + return np.squeeze(np.max(array,axis=0)) def stream_correlation(array1, array2): diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index d31e31d35..10e047f20 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -130,7 +130,7 @@ def main(): old_data_per_point[dpp_name] = data_per_point[dpp_name] out_sft = sft.from_sft(sft.streamlines, sft, data_per_point=old_data_per_point) - + save_tractogram(out_sft, args.out_tractogram) From 507df9ff9977247216803c4550d6d893625e397a Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 5 Feb 2024 12:15:25 -0500 Subject: [PATCH 52/97] pep8 issues --- scilpy/tractograms/dps_and_dpp_management.py | 1 + scilpy/tractograms/streamline_operations.py | 3 ++- scripts/scil_tractogram_dpp_math.py | 6 ++---- scripts/scil_tractogram_project_map_to_streamlines.py | 9 +++++---- 4 files changed, 10 insertions(+), 9 deletions(-) diff --git a/scilpy/tractograms/dps_and_dpp_management.py b/scilpy/tractograms/dps_and_dpp_management.py index cf79189be..373add8a5 100644 --- a/scilpy/tractograms/dps_and_dpp_management.py +++ b/scilpy/tractograms/dps_and_dpp_management.py @@ -1,5 +1,6 @@ import numpy as np + def project_map_to_streamlines(sft, map, endpoints_only=False): """ Projects a map onto the points of streamlines. diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index 8324d0e94..1b5655cac 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -474,7 +474,8 @@ def perform_operation_per_streamline(op_name, sft, dpp_name='metric', return new_data_per_streamline -def perform_pairwise_streamline_operation_on_endpoints(op_name, sft, dpp_name='metric'): +def perform_pairwise_streamline_operation_on_endpoints(op_name, sft, + dpp_name='metric'): """Peforms an operation across endpoints for each streamline. Parameters diff --git a/scripts/scil_tractogram_dpp_math.py b/scripts/scil_tractogram_dpp_math.py index d058bb4ea..92584a04d 100755 --- a/scripts/scil_tractogram_dpp_math.py +++ b/scripts/scil_tractogram_dpp_math.py @@ -153,13 +153,11 @@ def main(): data_per_streamline = {} for in_dpp_name, out_name in zip(args.in_dpp_name, args.out_name): - - # Perform the requested operation. if args.operation == 'correlation': logging.info('Performing {} across endpoint data and saving as ' 'new dpp {}'.format( - args.operation, out_name)) + args.operation, out_name)) new_dps = perform_pairwise_streamline_operation_on_endpoints( args.operation, sft, in_dpp_name) @@ -182,7 +180,6 @@ def main(): args.operation, sft, in_dpp_name, args.endpoints_only) data_per_streamline[out_name] = new_data_per_streamline - if args.overwrite_data: new_sft = sft.from_sft(sft.streamlines, sft, data_per_point=data_per_point, @@ -219,5 +216,6 @@ def main(): save_tractogram(new_sft, args.out_tractogram, bbox_valid_check=args.bbox_check) + if __name__ == "__main__": main() diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index 10e047f20..18d0a64ed 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -6,8 +6,8 @@ The default options will take data from a nifti image (3D or 4D) and project it onto the points of streamlines. If the image is 4D, the data -is stored as a list of 1D arrays per streamline. If the image is 3D, the data is stored -as a list of values per streamline. +is stored as a list of 1D arrays per streamline. If the image is 3D, +the data is stored as a list of values per streamline. """ import argparse @@ -78,7 +78,8 @@ def main(): sft.to_corner() if len(sft.streamlines) == 0: - logging.warning('Empty bundle file {}. Skipping'.format(args.in_tractogram)) + logging.warning('Empty bundle file {}. Skipping'.format( + args.in_tractogram)) return # Check to see if the number of maps and dpp_names are the same @@ -130,7 +131,7 @@ def main(): old_data_per_point[dpp_name] = data_per_point[dpp_name] out_sft = sft.from_sft(sft.streamlines, sft, data_per_point=old_data_per_point) - + save_tractogram(out_sft, args.out_tractogram) From 3a59f6a2f7662497f689a14808ffd41dc65e1326 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Wed, 7 Feb 2024 14:14:07 -0500 Subject: [PATCH 53/97] pep8 --- scilpy/tractograms/streamline_operations.py | 9 +++++---- scripts/scil_tractogram_dpp_math.py | 3 +-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index 1b5655cac..3a98316e0 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -503,18 +503,19 @@ def perform_pairwise_streamline_operation_on_endpoints(op_name, sft, def stream_mean(array): - return np.squeeze(np.mean(array,axis=0)) + return np.squeeze(np.mean(array, axis=0)) def stream_sum(array): - return np.squeeze(np.sum(array,axis=0)) + return np.squeeze(np.sum(array, axis=0)) def stream_min(array): - return np.squeeze(np.min(array,axis=0)) + return np.squeeze(np.min(array, axis=0)) + def stream_max(array): - return np.squeeze(np.max(array,axis=0)) + return np.squeeze(np.max(array, axis=0)) def stream_correlation(array1, array2): diff --git a/scripts/scil_tractogram_dpp_math.py b/scripts/scil_tractogram_dpp_math.py index 92584a04d..12034bf8f 100755 --- a/scripts/scil_tractogram_dpp_math.py +++ b/scripts/scil_tractogram_dpp_math.py @@ -156,8 +156,7 @@ def main(): # Perform the requested operation. if args.operation == 'correlation': logging.info('Performing {} across endpoint data and saving as ' - 'new dpp {}'.format( - args.operation, out_name)) + 'new dpp {}'.format(args.operation, out_name)) new_dps = perform_pairwise_streamline_operation_on_endpoints( args.operation, sft, in_dpp_name) From 1f69f940a152f95625e7950afbe3069e2763ff6f Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Wed, 7 Feb 2024 14:14:47 -0500 Subject: [PATCH 54/97] remove py version file --- .python-version | 1 - 1 file changed, 1 deletion(-) delete mode 100644 .python-version diff --git a/.python-version b/.python-version deleted file mode 100644 index c8cfe3959..000000000 --- a/.python-version +++ /dev/null @@ -1 +0,0 @@ -3.10 From 47d4a21fe4509382c0ffe298e37f9d024ee840c1 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Wed, 7 Feb 2024 14:15:19 -0500 Subject: [PATCH 55/97] added utf header --- scilpy/tractograms/dps_and_dpp_management.py | 1 + 1 file changed, 1 insertion(+) diff --git a/scilpy/tractograms/dps_and_dpp_management.py b/scilpy/tractograms/dps_and_dpp_management.py index 373add8a5..d919e4788 100644 --- a/scilpy/tractograms/dps_and_dpp_management.py +++ b/scilpy/tractograms/dps_and_dpp_management.py @@ -1,3 +1,4 @@ +# -*- coding: utf-8 -*- import numpy as np From 8e8bb69838430c9e05c0d893b6ff0f46e6540576 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 8 Feb 2024 14:24:01 -0500 Subject: [PATCH 56/97] removed py version file --- scripts/tests/.python-version | 1 - 1 file changed, 1 deletion(-) delete mode 100644 scripts/tests/.python-version diff --git a/scripts/tests/.python-version b/scripts/tests/.python-version deleted file mode 100644 index c8cfe3959..000000000 --- a/scripts/tests/.python-version +++ /dev/null @@ -1 +0,0 @@ -3.10 From 9ea5ec0bdaa099b161a0f6f4e910332e65e639c2 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Thu, 8 Feb 2024 14:29:37 -0500 Subject: [PATCH 57/97] added back py ver file in root --- .python-version | 1 + 1 file changed, 1 insertion(+) create mode 100644 .python-version diff --git a/.python-version b/.python-version new file mode 100644 index 000000000..7c7a975f4 --- /dev/null +++ b/.python-version @@ -0,0 +1 @@ +3.10 \ No newline at end of file From 2e801b7045795f24a7c4936fa7e712522da4f398 Mon Sep 17 00:00:00 2001 From: karp2601 Date: Tue, 13 Feb 2024 09:57:45 -0500 Subject: [PATCH 58/97] Adding msmt support --- scripts/scil_frf_set_diffusivities.py | 22 ++++++++++++++++---- scripts/tests/test_frf_set_diffusivities.py | 23 ++++++++++++++++++--- 2 files changed, 38 insertions(+), 7 deletions(-) diff --git a/scripts/scil_frf_set_diffusivities.py b/scripts/scil_frf_set_diffusivities.py index 6ae138701..dbec3dbec 100755 --- a/scripts/scil_frf_set_diffusivities.py +++ b/scripts/scil_frf_set_diffusivities.py @@ -6,7 +6,8 @@ Use this script when you want to use a fixed response function and keep the mean b0. -The FRF file is obtained from scil_frf_ssst.py +The FRF file is obtained from scil_frf_ssst.py or scil_frf_msmt.py in the case +of multi-shell data. Formerly: scil_set_response_function.py """ @@ -30,7 +31,9 @@ def _build_arg_parser(): p.add_argument('new_frf', metavar='tuple', help='Replace the response function with\n' 'this fiber response function x 10**-4 (e.g. ' - '15,4,4).') + '15,4,4). \nIf multi-shell, write the first shell' + ', then the second shell, \nand the third, etc ' + '(e.g. 15,4,4,13,5,5,12,5,5).') p.add_argument('output_frf_file', metavar='output', help='Path of the new FRF file.') p.add_argument('--no_factor', action='store_true', @@ -51,13 +54,24 @@ def main(): assert_inputs_exist(parser, args.frf_file) assert_outputs_exist(parser, args, args.output_frf_file) - frf_file = np.array(np.loadtxt(args.frf_file)) + frf_file = np.array(np.loadtxt(args.frf_file)).T new_frf = np.array(literal_eval(args.new_frf), dtype=np.float64) if not args.no_factor: new_frf *= 10 ** -4 b0_mean = frf_file[3] - response = np.concatenate([new_frf, [b0_mean]]) + if new_frf.shape[0] % 3 != 0: + raise ValueError('Inputed new frf is not valid. There should be ' + 'three values per shell, and thus the total number ' + 'of values should be a multiple of three.') + + nb_shells = int(new_frf.shape[0] / 3) + new_frf = new_frf.reshape((nb_shells, 3)) + + response = np.empty((nb_shells, 4)) + response[:, 0:3] = new_frf + response[:, 3] = b0_mean + np.savetxt(args.output_frf_file, response) diff --git a/scripts/tests/test_frf_set_diffusivities.py b/scripts/tests/test_frf_set_diffusivities.py index 631796cde..9ceaeb0fa 100644 --- a/scripts/tests/test_frf_set_diffusivities.py +++ b/scripts/tests/test_frf_set_diffusivities.py @@ -7,7 +7,8 @@ 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=['processing.zip']) +fetch_data(get_testing_files_dict(), + keys=['processing.zip', 'commit_amico.zip']) tmp_dir = tempfile.TemporaryDirectory() @@ -16,10 +17,26 @@ def test_help_option(script_runner): assert ret.success -def test_execution_processing(script_runner): +def test_execution_processing_ssst(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) in_frf = os.path.join(get_home(), 'processing', 'frf.txt') ret = script_runner.run('scil_frf_set_diffusivities.py', in_frf, - '15,4,4', 'nfrf.txt') + '15,4,4', 'new_frf.txt', '-f') assert ret.success + +def test_execution_processing_msmt(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + in_frf = os.path.join(get_home(), 'commit_amico', + 'wm_frf.txt') + ret = script_runner.run('scil_frf_set_diffusivities.py', in_frf, + '15,4,4,13,4,4,12,5,5', 'new_frf.txt', '-f') + assert ret.success + +def test_execution_processing__wrong_input(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + in_frf = os.path.join(get_home(), 'commit_amico', + 'wm_frf.txt') + ret = script_runner.run('scil_frf_set_diffusivities.py', in_frf, + '15,4,4,13,4,4', 'new_frf.txt', '-f') + assert not ret.success From 828b90d99aa44f09386b281f01843c8e62df8821 Mon Sep 17 00:00:00 2001 From: karp2601 Date: Tue, 13 Feb 2024 10:01:42 -0500 Subject: [PATCH 59/97] Pep 8 --- scripts/tests/test_frf_set_diffusivities.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/scripts/tests/test_frf_set_diffusivities.py b/scripts/tests/test_frf_set_diffusivities.py index 9ceaeb0fa..673df6e74 100644 --- a/scripts/tests/test_frf_set_diffusivities.py +++ b/scripts/tests/test_frf_set_diffusivities.py @@ -25,6 +25,7 @@ def test_execution_processing_ssst(script_runner): '15,4,4', 'new_frf.txt', '-f') assert ret.success + def test_execution_processing_msmt(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) in_frf = os.path.join(get_home(), 'commit_amico', @@ -33,6 +34,7 @@ def test_execution_processing_msmt(script_runner): '15,4,4,13,4,4,12,5,5', 'new_frf.txt', '-f') assert ret.success + def test_execution_processing__wrong_input(script_runner): os.chdir(os.path.expanduser(tmp_dir.name)) in_frf = os.path.join(get_home(), 'commit_amico', From be217c7345cce2a1db2968accbc1cddfa738c2e2 Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Thu, 14 Dec 2023 11:58:32 -0500 Subject: [PATCH 60/97] normalize_bvecs: Remove unused variable (verified in all of our scripts) --- scilpy/gradients/bvec_bval_tools.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/scilpy/gradients/bvec_bval_tools.py b/scilpy/gradients/bvec_bval_tools.py index 88c13f61e..8893fa1e1 100644 --- a/scilpy/gradients/bvec_bval_tools.py +++ b/scilpy/gradients/bvec_bval_tools.py @@ -37,7 +37,7 @@ def is_normalized_bvecs(bvecs): bvecs_norm == 0)) -def normalize_bvecs(bvecs, filename=None): +def normalize_bvecs(bvecs): """ Normalize b-vectors @@ -45,8 +45,6 @@ def normalize_bvecs(bvecs, filename=None): ---------- bvecs : (N, 3) array input b-vectors (N, 3) array - filename : string - output filename where to save the normalized bvecs Returns ------- @@ -58,10 +56,6 @@ def normalize_bvecs(bvecs, filename=None): idx = bvecs_norm != 0 bvecs[idx] /= bvecs_norm[idx, None] - if filename is not None: - logging.info('Saving new bvecs: {}'.format(filename)) - np.savetxt(filename, np.transpose(bvecs), "%.8f") - return bvecs From b2b93f81f65597c4c56a8545697a666733a314aa Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Thu, 14 Dec 2023 12:22:22 -0500 Subject: [PATCH 61/97] Prepare tests --- scilpy/gradients/tests/test_bvec_bval_tools.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/scilpy/gradients/tests/test_bvec_bval_tools.py b/scilpy/gradients/tests/test_bvec_bval_tools.py index 92d005ee7..5f0db7498 100644 --- a/scilpy/gradients/tests/test_bvec_bval_tools.py +++ b/scilpy/gradients/tests/test_bvec_bval_tools.py @@ -1,17 +1,22 @@ # -*- coding: utf-8 -*- import numpy as np -from scilpy.gradients.bvec_bval_tools import round_bvals_to_shell +from scilpy.gradients.bvec_bval_tools import (is_normalized_bvecs, + round_bvals_to_shell, normalize_bvecs) + +bvecs = np.asarray([[1.0, 1.0, 1.0], + [1.0, 0.0, 1.0], + [0.0, 1.0, 0.0]]) def test_is_normalized_bvecs(): - # toDO - pass + assert not is_normalized_bvecs(bvecs) + assert is_normalized_bvecs( + bvecs / np.linalg.norm(bvecs, axis=1, keepdims=True)) def test_normalize_bvecs(): - # toDo - pass + assert is_normalized_bvecs(normalize_bvecs(bvecs)) def test_check_b0_threshold(): From 2e9c3b007f6fbe60a9ef3f7aabec8bda136bb91d Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Tue, 23 Jan 2024 14:11:40 -0500 Subject: [PATCH 62/97] mrtrix2fsl and fsl2mrtrix: Moved to io.gradients. They load and save. --- scilpy/gradients/bvec_bval_tools.py | 72 ------------------- .../gradients/tests/test_bvec_bval_tools.py | 17 ++--- scilpy/io/gradients.py | 63 ++++++++++++++++ scripts/scil_NODDI_maps.py | 3 +- scripts/scil_freewater_maps.py | 3 +- scripts/scil_gradients_convert.py | 2 +- scripts/scil_tractogram_commit.py | 3 +- 7 files changed, 74 insertions(+), 89 deletions(-) diff --git a/scilpy/gradients/bvec_bval_tools.py b/scilpy/gradients/bvec_bval_tools.py index 8893fa1e1..c0fadebe4 100644 --- a/scilpy/gradients/bvec_bval_tools.py +++ b/scilpy/gradients/bvec_bval_tools.py @@ -6,9 +6,6 @@ from dipy.core.gradients import get_bval_indices import numpy as np -from scilpy.io.gradients import (save_gradient_sampling_fsl, - save_gradient_sampling_mrtrix) - DEFAULT_B0_THRESHOLD = 20 @@ -114,75 +111,6 @@ def check_b0_threshold( return b0_thr -def fsl2mrtrix(fsl_bval_filename, fsl_bvec_filename, mrtrix_filename): - """ - Convert a fsl dir_grad.bvec/.bval files to mrtrix encoding.b file. - - Parameters - ---------- - fsl_bval_filename: str - path to input fsl bval file. - fsl_bvec_filename: str - path to input fsl bvec file. - mrtrix_filename : str - path to output mrtrix encoding.b file. - - Returns - ------- - """ - - shells = np.loadtxt(fsl_bval_filename) - points = np.loadtxt(fsl_bvec_filename) - bvals = np.unique(shells).tolist() - - # Remove .bval and .bvec if present - mrtrix_filename = mrtrix_filename.replace('.b', '') - - if not points.shape[0] == 3: - points = points.transpose() - logging.warning('WARNING: Your bvecs seem transposed. ' + - 'Transposing them.') - - shell_idx = [int(np.where(bval == bvals)[0]) for bval in shells] - save_gradient_sampling_mrtrix(points, - shell_idx, - bvals, - mrtrix_filename + '.b') - - -def mrtrix2fsl(mrtrix_filename, fsl_filename): - """ - Convert a mrtrix encoding.b file to fsl dir_grad.bvec/.bval files. - - Parameters - ---------- - mrtrix_filename : str - path to mrtrix encoding.b file. - fsl_bval_filename: str - path to the output fsl files. Files will be named - fsl_bval_filename.bval and fsl_bval_filename.bvec. - """ - # Remove .bval and .bvec if present - fsl_filename = fsl_filename.replace('.bval', '') - fsl_filename = fsl_filename.replace('.bvec', '') - - mrtrix_b = np.loadtxt(mrtrix_filename) - if not len(mrtrix_b.shape) == 2 or not mrtrix_b.shape[1] == 4: - raise ValueError('mrtrix file must have 4 columns') - - points = np.array([mrtrix_b[:, 0], mrtrix_b[:, 1], mrtrix_b[:, 2]]) - shells = np.array(mrtrix_b[:, 3]) - - bvals = np.unique(shells).tolist() - shell_idx = [int(np.where(bval == bvals)[0]) for bval in shells] - - save_gradient_sampling_fsl(points, - shell_idx, - bvals, - filename_bval=fsl_filename + '.bval', - filename_bvec=fsl_filename + '.bvec') - - def identify_shells(bvals, threshold=40.0, roundCentroids=False, sort=False): """ Guessing the shells from the b-values. Returns the list of shells and, for diff --git a/scilpy/gradients/tests/test_bvec_bval_tools.py b/scilpy/gradients/tests/test_bvec_bval_tools.py index 5f0db7498..f23ce19e0 100644 --- a/scilpy/gradients/tests/test_bvec_bval_tools.py +++ b/scilpy/gradients/tests/test_bvec_bval_tools.py @@ -1,8 +1,9 @@ # -*- coding: utf-8 -*- import numpy as np -from scilpy.gradients.bvec_bval_tools import (is_normalized_bvecs, - round_bvals_to_shell, normalize_bvecs) +from scilpy.gradients.bvec_bval_tools import ( + check_b0_threshold, is_normalized_bvecs, normalize_bvecs, + round_bvals_to_shell) bvecs = np.asarray([[1.0, 1.0, 1.0], [1.0, 0.0, 1.0], @@ -20,17 +21,7 @@ def test_normalize_bvecs(): def test_check_b0_threshold(): - # toDo - pass - - -def test_fsl2mrtrix(): - # toDo - pass - - -def test_mrtrix2fsl(): - # toDo + # toDo To be modified (see PR#867). pass diff --git a/scilpy/io/gradients.py b/scilpy/io/gradients.py index 3c168c2e7..c72103e61 100644 --- a/scilpy/io/gradients.py +++ b/scilpy/io/gradients.py @@ -6,6 +6,69 @@ import numpy as np +def fsl2mrtrix(fsl_bval_filename, fsl_bvec_filename, mrtrix_filename): + """ + Convert a fsl dir_grad.bvec/.bval files to mrtrix encoding.b file. + Saves the result. + + Parameters + ---------- + fsl_bval_filename: str + path to input fsl bval file. + fsl_bvec_filename: str + path to input fsl bvec file. + mrtrix_filename : str + path to output mrtrix encoding.b file. + """ + shells = np.loadtxt(fsl_bval_filename) + points = np.loadtxt(fsl_bvec_filename) + bvals = np.unique(shells).tolist() + + # Remove .bval and .bvec if present + mrtrix_filename = mrtrix_filename.replace('.b', '') + + if not points.shape[0] == 3: + points = points.transpose() + logging.warning('WARNING: Your bvecs seem transposed. ' + + 'Transposing them.') + + shell_idx = [int(np.where(bval == bvals)[0]) for bval in shells] + save_gradient_sampling_mrtrix(points, shell_idx, bvals, + mrtrix_filename + '.b') + + +def mrtrix2fsl(mrtrix_filename, fsl_filename): + """ + Convert a mrtrix encoding.b file to fsl dir_grad.bvec/.bval files. + Saves the result. + + Parameters + ---------- + mrtrix_filename : str + path to mrtrix encoding.b file. + fsl_filename: str + path to the output fsl files. Files will be named + fsl_bval_filename.bval and fsl_bval_filename.bvec. + """ + # Remove .bval and .bvec if present + fsl_filename = fsl_filename.replace('.bval', '') + fsl_filename = fsl_filename.replace('.bvec', '') + + mrtrix_b = np.loadtxt(mrtrix_filename) + if not len(mrtrix_b.shape) == 2 or not mrtrix_b.shape[1] == 4: + raise ValueError('mrtrix file must have 4 columns') + + points = np.array([mrtrix_b[:, 0], mrtrix_b[:, 1], mrtrix_b[:, 2]]) + shells = np.array(mrtrix_b[:, 3]) + + bvals = np.unique(shells).tolist() + shell_idx = [int(np.where(bval == bvals)[0]) for bval in shells] + + save_gradient_sampling_fsl(points, shell_idx, bvals, + filename_bval=fsl_filename + '.bval', + filename_bvec=fsl_filename + '.bvec') + + def save_gradient_sampling_mrtrix(bvecs, shell_idx, bvals, filename): """ Save table gradient (MRtrix format) diff --git a/scripts/scil_NODDI_maps.py b/scripts/scil_NODDI_maps.py index bb720a85c..b2d890bf8 100755 --- a/scripts/scil_NODDI_maps.py +++ b/scripts/scil_NODDI_maps.py @@ -20,13 +20,14 @@ from dipy.io.gradients import read_bvals_bvecs import numpy as np +from scilpy.io.gradients import fsl2mrtrix from scilpy.io.utils import (add_overwrite_arg, add_processes_arg, add_verbose_arg, assert_inputs_exist, assert_output_dirs_exist_and_empty, redirect_stdout_c) -from scilpy.gradients.bvec_bval_tools import fsl2mrtrix, identify_shells +from scilpy.gradients.bvec_bval_tools import identify_shells EPILOG = """ Reference: diff --git a/scripts/scil_freewater_maps.py b/scripts/scil_freewater_maps.py index b3fb44ba6..643c75380 100755 --- a/scripts/scil_freewater_maps.py +++ b/scripts/scil_freewater_maps.py @@ -20,13 +20,14 @@ from dipy.io.gradients import read_bvals_bvecs import numpy as np +from scilpy.io.gradients import fsl2mrtrix from scilpy.io.utils import (add_overwrite_arg, add_processes_arg, add_verbose_arg, assert_inputs_exist, assert_output_dirs_exist_and_empty, redirect_stdout_c) -from scilpy.gradients.bvec_bval_tools import fsl2mrtrix, identify_shells +from scilpy.gradients.bvec_bval_tools import identify_shells EPILOG = """ diff --git a/scripts/scil_gradients_convert.py b/scripts/scil_gradients_convert.py index f3cfd0025..f2c2d2766 100755 --- a/scripts/scil_gradients_convert.py +++ b/scripts/scil_gradients_convert.py @@ -14,7 +14,7 @@ from scilpy.io.utils import (assert_gradients_filenames_valid, assert_inputs_exist, assert_outputs_exist, add_overwrite_arg, add_verbose_arg) -from scilpy.gradients.bvec_bval_tools import fsl2mrtrix, mrtrix2fsl +from scilpy.io.gradients import fsl2mrtrix, mrtrix2fsl def _build_arg_parser(): diff --git a/scripts/scil_tractogram_commit.py b/scripts/scil_tractogram_commit.py index 736a62229..b52c6b36b 100755 --- a/scripts/scil_tractogram_commit.py +++ b/scripts/scil_tractogram_commit.py @@ -84,6 +84,7 @@ import numpy as np import nibabel as nib +from scilpy.io.gradients import fsl2mrtrix from scilpy.io.streamlines import (reconstruct_streamlines, reconstruct_streamlines_from_hdf5) from scilpy.io.utils import (add_overwrite_arg, @@ -92,7 +93,7 @@ assert_inputs_exist, assert_output_dirs_exist_and_empty, redirect_stdout_c) -from scilpy.gradients.bvec_bval_tools import fsl2mrtrix, identify_shells +from scilpy.gradients.bvec_bval_tools import identify_shells EPILOG = """ From 21edbd31ff5746b8de79690fdb2294a5b5fea0c1 Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Tue, 23 Jan 2024 14:52:14 -0500 Subject: [PATCH 63/97] Test identify_shells. Added a warning. Renamed tol --- scilpy/gradients/bvec_bval_tools.py | 29 ++++++++++----- .../gradients/tests/test_bvec_bval_tools.py | 37 +++++++++++++++++-- scilpy/image/volume_operations.py | 4 +- scripts/scil_NODDI_maps.py | 2 +- scripts/scil_freewater_maps.py | 2 +- scripts/scil_tractogram_commit.py | 2 +- 6 files changed, 58 insertions(+), 18 deletions(-) diff --git a/scilpy/gradients/bvec_bval_tools.py b/scilpy/gradients/bvec_bval_tools.py index c0fadebe4..252fd9d74 100644 --- a/scilpy/gradients/bvec_bval_tools.py +++ b/scilpy/gradients/bvec_bval_tools.py @@ -111,7 +111,7 @@ def check_b0_threshold( return b0_thr -def identify_shells(bvals, threshold=40.0, roundCentroids=False, sort=False): +def identify_shells(bvals, tol=40.0, round_centroids=False, sort=False): """ Guessing the shells from the b-values. Returns the list of shells and, for each b-value, the associated shell. @@ -128,10 +128,11 @@ def identify_shells(bvals, threshold=40.0, roundCentroids=False, sort=False): ---------- bvals: array (N,) Array of bvals - threshold: float - 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. - roundCentroids: bool + tol: float + Limit difference to centroid to consider that a b-value is on an + existing shell. On or above this limit, the b-value is placed on a new + shell. + round_centroids: bool If true will round shell values to the nearest 10. sort: bool Sort centroids and shell_indices associated. @@ -151,7 +152,7 @@ def identify_shells(bvals, threshold=40.0, roundCentroids=False, sort=False): bval_centroids = [bvals[0]] for bval in bvals[1:]: diffs = np.abs(np.asarray(bval_centroids, dtype=float) - bval) - if not len(np.where(diffs < threshold)[0]): + if not len(np.where(diffs < tol)[0]): # Found no bval in bval centroids close enough to the current one. # Create new centroid (i.e. new shell) bval_centroids.append(bval) @@ -163,15 +164,23 @@ def identify_shells(bvals, threshold=40.0, roundCentroids=False, sort=False): shell_indices = np.argmin(np.abs(bvals_for_diffs - centroids), axis=1) - if roundCentroids: + if round_centroids: centroids = np.round(centroids, decimals=-1) + # Ex: with bvals [0, 5], threshold 5, we get centroids 0, 5. + # Rounded, we get centroids 0, 0. + if len(np.unique(centroids)) != len(centroids): + logging.warning("With option to round the centroids to the " + "nearest 10, with tolerance {}, we get unclear " + "division of the shells. Use this data carefully." + .format(tol)) + if sort: sort_index = np.argsort(centroids) - sorted_centroids = np.zeros(centroids.shape) - sorted_indices = np.zeros(shell_indices.shape) + sorted_centroids = centroids[sort_index] + + sorted_indices = np.zeros(shell_indices.shape, dtype=int) for i in range(len(centroids)): - sorted_centroids[i] = centroids[sort_index[i]] sorted_indices[shell_indices == i] = sort_index[i] return sorted_centroids, sorted_indices diff --git a/scilpy/gradients/tests/test_bvec_bval_tools.py b/scilpy/gradients/tests/test_bvec_bval_tools.py index f23ce19e0..2f900fd78 100644 --- a/scilpy/gradients/tests/test_bvec_bval_tools.py +++ b/scilpy/gradients/tests/test_bvec_bval_tools.py @@ -3,7 +3,7 @@ from scilpy.gradients.bvec_bval_tools import ( check_b0_threshold, is_normalized_bvecs, normalize_bvecs, - round_bvals_to_shell) + round_bvals_to_shell, identify_shells) bvecs = np.asarray([[1.0, 1.0, 1.0], [1.0, 0.0, 1.0], @@ -26,8 +26,39 @@ def test_check_b0_threshold(): def test_identify_shells(): - # toDo - pass + def _subtest_identify_shells(bvals, threshold, + expected_raw_centroids, expected_raw_shells, + expected_round_sorted_centroids, + expected_round_sorted_shells): + bvals = np.asarray(bvals) + + # 1) Not rounded, not sorted + c, s = identify_shells(bvals, threshold) + assert np.array_equal(c, expected_raw_centroids) + assert np.array_equal(s, expected_raw_shells) + + # 2) Rounded, sorted + c, s = identify_shells(bvals, threshold, round_centroids=True, + sort=True) + assert np.array_equal(c, expected_round_sorted_centroids) + assert np.array_equal(s, expected_round_sorted_shells) + + # Test 1. All easy. Over the limit for 0, 5, 15. Clear difference for + # 100, 2000. + _subtest_identify_shells(bvals=[0, 0, 5, 15, 2000, 100], threshold=50, + expected_raw_centroids=[0, 2000, 100], + expected_raw_shells=[0, 0, 0, 0, 1, 2], + expected_round_sorted_centroids=[0, 100, 2000], + expected_round_sorted_shells=[0, 0, 0, 0, 2, 1]) + + # Test 2. Threshold on the limit. + # Additional difficulty with option rounded: two shells with the same + # value, but a warning is printed. Should it raise an error? + _subtest_identify_shells(bvals=[0, 0, 5, 2000, 100], threshold=5, + expected_raw_centroids=[0, 5, 2000, 100], + expected_raw_shells=[0, 0, 1, 2, 3], + expected_round_sorted_centroids=[0, 0, 100, 2000], + expected_round_sorted_shells=[0, 0, 1, 3, 2]) def test_str_to_axis_index(): diff --git a/scilpy/image/volume_operations.py b/scilpy/image/volume_operations.py index 463c39c8d..16fd8b04f 100644 --- a/scilpy/image/volume_operations.py +++ b/scilpy/image/volume_operations.py @@ -279,8 +279,8 @@ def compute_snr(dwi, bval, bvec, b0_thr, mask, mask = get_data_as_mask(mask, dtype=bool) if split_shells: - centroids, shell_indices = identify_shells(bval, threshold=40.0, - roundCentroids=False, + centroids, shell_indices = identify_shells(bval, tol=40.0, + round_centroids=False, sort=False) bval = centroids[shell_indices] diff --git a/scripts/scil_NODDI_maps.py b/scripts/scil_NODDI_maps.py index b2d890bf8..82b12e40a 100755 --- a/scripts/scil_NODDI_maps.py +++ b/scripts/scil_NODDI_maps.py @@ -119,7 +119,7 @@ def main(): bvals, _ = read_bvals_bvecs(args.in_bval, args.in_bvec) shells_centroids, indices_shells = identify_shells(bvals, args.b_thr, - roundCentroids=True) + 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) diff --git a/scripts/scil_freewater_maps.py b/scripts/scil_freewater_maps.py index 643c75380..7e108148a 100755 --- a/scripts/scil_freewater_maps.py +++ b/scripts/scil_freewater_maps.py @@ -126,7 +126,7 @@ def main(): bvals, _ = read_bvals_bvecs(args.in_bval, args.in_bvec) shells_centroids, indices_shells = identify_shells(bvals, args.b_thr, - roundCentroids=True) + 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) diff --git a/scripts/scil_tractogram_commit.py b/scripts/scil_tractogram_commit.py index b52c6b36b..f03b9aba2 100755 --- a/scripts/scil_tractogram_commit.py +++ b/scripts/scil_tractogram_commit.py @@ -385,7 +385,7 @@ def main(): 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, - roundCentroids=True) + 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) From 514e277ada99063710789d6cecdda34572c9ee00 Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Tue, 23 Jan 2024 15:02:21 -0500 Subject: [PATCH 64/97] Add test flip and swap axis --- scilpy/gradients/bvec_bval_tools.py | 21 +++++++++------- .../gradients/tests/test_bvec_bval_tools.py | 25 +++++++++++++------ scripts/scil_gradients_modify_axes.py | 12 ++++++++- 3 files changed, 41 insertions(+), 17 deletions(-) diff --git a/scilpy/gradients/bvec_bval_tools.py b/scilpy/gradients/bvec_bval_tools.py index 252fd9d74..1f99d95a1 100644 --- a/scilpy/gradients/bvec_bval_tools.py +++ b/scilpy/gradients/bvec_bval_tools.py @@ -28,7 +28,6 @@ def is_normalized_bvecs(bvecs): ------- True/False """ - bvecs_norm = np.linalg.norm(bvecs, axis=1) return np.all(np.logical_or(np.abs(bvecs_norm - 1) < 1e-3, bvecs_norm == 0)) @@ -48,7 +47,7 @@ def normalize_bvecs(bvecs): bvecs : (N, 3) normalized b-vectors """ - + bvecs = bvecs.copy() # Avoid in-place modification. bvecs_norm = np.linalg.norm(bvecs, axis=1) idx = bvecs_norm != 0 bvecs[idx] /= bvecs_norm[idx, None] @@ -217,12 +216,13 @@ def flip_gradient_sampling(bvecs, axes, sampling_type): Parameters ---------- bvecs: np.ndarray - Loaded bvecs. In the case 'mrtrix' the bvecs actually also contain the - bvals. + bvecs loaded directly, not re-formatted. Careful! Must respect the + format (not verified here). axes: list of int - List of axes to flip (e.g. [0, 1]) + List of axes to flip (e.g. [0, 1]). See str_to_axis_index. sampling_type: str - Either 'mrtrix' or 'fsl'. + Either 'mrtrix': bvecs are of shape (N, 4) or + 'fsl': bvecs are of shape (3, N) Returns ------- @@ -230,6 +230,8 @@ def flip_gradient_sampling(bvecs, axes, sampling_type): The final bvecs. """ assert sampling_type in ['mrtrix', 'fsl'] + + bvecs = bvecs.copy() # Avoid in-place modification. if sampling_type == 'mrtrix': for axis in axes: bvecs[:, axis] *= -1 @@ -246,12 +248,13 @@ def swap_gradient_axis(bvecs, final_order, sampling_type): Parameters ---------- bvecs: np.array - Loaded bvecs. In the case 'mrtrix' the bvecs actually also contain the - bvals. + bvecs loaded directly, not re-formatted. Careful! Must respect the + format (not verified here). final_order: new order Final order (ex, 2 1 0) sampling_type: str - Either 'mrtrix' or 'fsl'. + Either 'mrtrix': bvecs are of shape (N, 4) or + 'fsl': bvecs are of shape (3, N) Returns ------- diff --git a/scilpy/gradients/tests/test_bvec_bval_tools.py b/scilpy/gradients/tests/test_bvec_bval_tools.py index 2f900fd78..64ef15f8f 100644 --- a/scilpy/gradients/tests/test_bvec_bval_tools.py +++ b/scilpy/gradients/tests/test_bvec_bval_tools.py @@ -3,11 +3,13 @@ from scilpy.gradients.bvec_bval_tools import ( check_b0_threshold, is_normalized_bvecs, normalize_bvecs, - round_bvals_to_shell, identify_shells) + round_bvals_to_shell, identify_shells, str_to_axis_index, flip_gradient_sampling, + swap_gradient_axis) bvecs = np.asarray([[1.0, 1.0, 1.0], [1.0, 0.0, 1.0], - [0.0, 1.0, 0.0]]) + [0.0, 1.0, 0.0], + [8.0, 1.0, 1.0]]) def test_is_normalized_bvecs(): @@ -63,17 +65,26 @@ def _subtest_identify_shells(bvals, threshold, def test_str_to_axis_index(): # Very simple, nothing to do - pass + assert str_to_axis_index('x') == 0 def test_flip_gradient_sampling(): - # toDo - pass + fsl_bvecs = bvecs.T + b = flip_gradient_sampling(fsl_bvecs, axes=[0], sampling_type='fsl') + assert np.array_equal(b, np.asarray([[-1.0, 1.0, 1.0], + [-1.0, 0.0, 1.0], + [-0.0, 1.0, 0.0], + [-8.0, 1.0, 1.0]]).T) def test_swap_gradient_axis(): - # toDo - pass + fsl_bvecs = bvecs.T + final_order = [1, 0, 2] + b = swap_gradient_axis(fsl_bvecs, final_order, sampling_type='fsl') + assert np.array_equal(b, np.asarray([[1.0, 1.0, 1.0], + [0.0, 1.0, 1.0], + [1.0, 0.0, 0.0], + [1.0, 8.0, 1.0]]).T) def test_round_bvals_to_shell(): diff --git a/scripts/scil_gradients_modify_axes.py b/scripts/scil_gradients_modify_axes.py index 99f90c73f..a642725c6 100755 --- a/scripts/scil_gradients_modify_axes.py +++ b/scripts/scil_gradients_modify_axes.py @@ -83,14 +83,24 @@ def main(): if len(np.unique(swapped_order)) != 3: parser.error("final_order should contain the three axis.") + # Could use io.gradients method, but we don't have a bvals file here, only + # treating with the bvecs. Loading directly. bvecs = np.loadtxt(args.in_gradient_sampling_file) if ext_in == '.bvec': - # Supposing FSL format + # Supposing FSL format: Per columns + if bvecs.shape[0] != 3: + parser.error("b-vectors format for a .b file should be FSL, " + "and contain 3 lines (x, y, z), but got {}" + .format(bvecs.shape[0])) bvecs = flip_gradient_sampling(bvecs, axes_to_flip, 'fsl') bvecs = swap_gradient_axis(bvecs, swapped_order, 'fsl') np.savetxt(args.out_gradient_sampling_file, bvecs, "%.8f") else: # ext == '.b': # Supposing mrtrix format + if bvecs.shape[1] != 4: + parser.error("b-vectors format for a .b file should be mrtrix, " + "and contain 4 columns (x, y, z, bval), but got {}" + .format(bvecs.shape[1])) bvecs = flip_gradient_sampling(bvecs, axes_to_flip, 'mrtrix') bvecs = swap_gradient_axis(bvecs, swapped_order, 'mrtrix') np.savetxt(args.out_gradient_sampling_file, bvecs, From 0d6e7a77b7ab8320919cb9fd589d23269b6bd3ab Mon Sep 17 00:00:00 2001 From: Emmanuelle Renauld Date: Wed, 24 Jan 2024 15:40:25 -0500 Subject: [PATCH 65/97] Answer Arnaud's offline comment --- scilpy/gradients/tests/test_bvec_bval_tools.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/scilpy/gradients/tests/test_bvec_bval_tools.py b/scilpy/gradients/tests/test_bvec_bval_tools.py index 64ef15f8f..3f413f3b1 100644 --- a/scilpy/gradients/tests/test_bvec_bval_tools.py +++ b/scilpy/gradients/tests/test_bvec_bval_tools.py @@ -66,6 +66,9 @@ def _subtest_identify_shells(bvals, threshold, def test_str_to_axis_index(): # Very simple, nothing to do assert str_to_axis_index('x') == 0 + assert str_to_axis_index('y') == 1 + assert str_to_axis_index('z') == 2 + def test_flip_gradient_sampling(): From f0bc75a698be98b56ac4e555e673b068c570e29d Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Wed, 14 Feb 2024 10:14:02 -0500 Subject: [PATCH 66/97] Add test axis 'v' and fix pep8 --- scilpy/gradients/tests/test_bvec_bval_tools.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/scilpy/gradients/tests/test_bvec_bval_tools.py b/scilpy/gradients/tests/test_bvec_bval_tools.py index 3f413f3b1..3d82cc188 100644 --- a/scilpy/gradients/tests/test_bvec_bval_tools.py +++ b/scilpy/gradients/tests/test_bvec_bval_tools.py @@ -2,8 +2,8 @@ import numpy as np from scilpy.gradients.bvec_bval_tools import ( - check_b0_threshold, is_normalized_bvecs, normalize_bvecs, - round_bvals_to_shell, identify_shells, str_to_axis_index, flip_gradient_sampling, + identify_shells, is_normalized_bvecs, flip_gradient_sampling, + normalize_bvecs, round_bvals_to_shell, str_to_axis_index, swap_gradient_axis) bvecs = np.asarray([[1.0, 1.0, 1.0], @@ -68,7 +68,7 @@ def test_str_to_axis_index(): assert str_to_axis_index('x') == 0 assert str_to_axis_index('y') == 1 assert str_to_axis_index('z') == 2 - + assert str_to_axis_index('v') is None def test_flip_gradient_sampling(): From dedb657c5474c6fd4eca27734a12dec0dcf835e0 Mon Sep 17 00:00:00 2001 From: karp2601 Date: Wed, 14 Feb 2024 11:49:19 -0500 Subject: [PATCH 67/97] Fixing pep8 --- scripts/scil_frf_set_diffusivities.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/scil_frf_set_diffusivities.py b/scripts/scil_frf_set_diffusivities.py index dbec3dbec..be564e3af 100755 --- a/scripts/scil_frf_set_diffusivities.py +++ b/scripts/scil_frf_set_diffusivities.py @@ -40,10 +40,10 @@ def _build_arg_parser(): help='If supplied, the fiber response function is\n' 'evaluated without the x 10**-4 factor. [%(default)s].' ) - + add_verbose_arg(p) add_overwrite_arg(p) - + return p From 1cf2a0f937134599f5f98b0a44b4659bf15f5383 Mon Sep 17 00:00:00 2001 From: AlexVCaron Date: Wed, 14 Feb 2024 14:57:19 -0500 Subject: [PATCH 68/97] disable test run on forks --- .github/workflows/test.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 082d54c0a..740cf4692 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -19,9 +19,10 @@ env: jobs: test: runs-on: scilus-runners + if: github.repository == 'scilus/scilpy' steps: - name: Checkout repository for merge - uses: actions/checkout@v4.1.1 + uses: actions/checkout@v4 - name: Fetch python version from repository id: python-selector From a7952308d2aed5060fbb108202981d3bad9e9734 Mon Sep 17 00:00:00 2001 From: AlexVCaron Date: Thu, 15 Feb 2024 12:59:31 -0500 Subject: [PATCH 69/97] write sources relative to repository --- .coveragerc | 1 + 1 file changed, 1 insertion(+) diff --git a/.coveragerc b/.coveragerc index 10530e8c5..e748ff743 100644 --- a/.coveragerc +++ b/.coveragerc @@ -2,6 +2,7 @@ branch = True concurrency = multiprocessing data_file = .test_reports/.coverage +relative_files = True source = scilpy/ scripts/ From 62c678841c08db98445352a6e0a71984b3ea6560 Mon Sep 17 00:00:00 2001 From: AlexVCaron Date: Thu, 15 Feb 2024 13:01:24 -0500 Subject: [PATCH 70/97] first try running only simple fast test to test upload and effect on codecov.io --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 740cf4692..9273974b6 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -53,7 +53,7 @@ jobs: - name: Run tests run: | export C_INCLUDE_PATH=$pythonLocation/include/python${{ steps.python-selector.outputs.python-version }}:$C_INCLUDE_PATH - pytest --cov-report term-missing:skip-covered + pytest --cov-report term-missing:skip-covered scripts/tests/test_bids_validate.py - name: Upload coverage reports to Codecov uses: codecov/codecov-action@v3 From 62e71e1500ff3860b3c1a9b18df2efac367a8c18 Mon Sep 17 00:00:00 2001 From: AlexVCaron Date: Thu, 15 Feb 2024 13:19:07 -0500 Subject: [PATCH 71/97] switch back to v4, to see if they fixed something --- .github/workflows/test.yml | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 9273974b6..5fc18ec10 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -56,16 +56,15 @@ jobs: pytest --cov-report term-missing:skip-covered scripts/tests/test_bids_validate.py - name: Upload coverage reports to Codecov - uses: codecov/codecov-action@v3 - env: - CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} + uses: codecov/codecov-action@v4 with: + token: ${{ secrets.CODECOV_TOKEN }} flags: unittests name: scilpy-unittests-${{ github.run_id }} verbose: true directory: .test_reports/ fail_ci_if_error: true - root_dir: $GITHUB_WORKSPACE/scilpy/ + plugin: pycoverage - name: Upload test reports and coverage to artifacts uses: actions/upload-artifact@v4.3.1 From 964ab4a172e17eb0344e3d894b638cf287e0222b Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Thu, 15 Feb 2024 13:21:53 -0500 Subject: [PATCH 72/97] Add test files with empty tests --- scilpy/dwi/tests/__init__.py | 0 scilpy/dwi/tests/test_operations.py | 9 +++++++++ scilpy/dwi/tests/test_utils.py | 9 +++++++++ 3 files changed, 18 insertions(+) create mode 100644 scilpy/dwi/tests/__init__.py create mode 100644 scilpy/dwi/tests/test_operations.py create mode 100644 scilpy/dwi/tests/test_utils.py diff --git a/scilpy/dwi/tests/__init__.py b/scilpy/dwi/tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/scilpy/dwi/tests/test_operations.py b/scilpy/dwi/tests/test_operations.py new file mode 100644 index 000000000..91d7ccb02 --- /dev/null +++ b/scilpy/dwi/tests/test_operations.py @@ -0,0 +1,9 @@ +# -*- coding: utf-8 -*- + + +def test_apply_bias_field(): + pass + + +def test_compute_dwi_attenuation(): + pass diff --git a/scilpy/dwi/tests/test_utils.py b/scilpy/dwi/tests/test_utils.py new file mode 100644 index 000000000..f6b4d5d51 --- /dev/null +++ b/scilpy/dwi/tests/test_utils.py @@ -0,0 +1,9 @@ +# -*- coding: utf-8 -*- + + +def test_extract_dwi_shell(): + pass + + +def test_extract_b0(): + pass From 09323a1de806ef97602eff71240f2b1cd32ffbbe Mon Sep 17 00:00:00 2001 From: AlexVCaron Date: Thu, 15 Feb 2024 13:35:50 -0500 Subject: [PATCH 73/97] use source pkgs instead of simple source --- .coveragerc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.coveragerc b/.coveragerc index e748ff743..9f2015df3 100644 --- a/.coveragerc +++ b/.coveragerc @@ -3,9 +3,9 @@ branch = True concurrency = multiprocessing data_file = .test_reports/.coverage relative_files = True -source = - scilpy/ - scripts/ +source_pkgs = + scilpy + scripts omit = scripts/tests/*.py scilpy/tests/**/*.py From 4ac0ac309e1436cb553a90a17c988710bd46b376 Mon Sep 17 00:00:00 2001 From: AlexVCaron Date: Thu, 15 Feb 2024 13:36:49 -0500 Subject: [PATCH 74/97] check reports before codecov run --- .github/workflows/test.yml | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 5fc18ec10..7820777ff 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -55,6 +55,13 @@ jobs: export C_INCLUDE_PATH=$pythonLocation/include/python${{ steps.python-selector.outputs.python-version }}:$C_INCLUDE_PATH pytest --cov-report term-missing:skip-covered scripts/tests/test_bids_validate.py + - name: Upload test reports and coverage to artifacts before codecov + uses: actions/upload-artifact@v4.3.1 + with: + name: test-reports-before + path: | + .test_reports/* + - name: Upload coverage reports to Codecov uses: codecov/codecov-action@v4 with: From 46ed33cc361023fd4a0f967c886af97a50aca326 Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Thu, 15 Feb 2024 13:48:23 -0500 Subject: [PATCH 75/97] Finish cleaning the main of dwi_ scripts, for eventual unit tests --- scilpy/dwi/operations.py | 87 ++++++++++++++++++++++ scilpy/dwi/tests/test_operations.py | 4 + scilpy/io/varian_fdf.py | 4 +- scripts/scil_dwi_detect_volume_outliers.py | 75 ++----------------- 4 files changed, 100 insertions(+), 70 deletions(-) diff --git a/scilpy/dwi/operations.py b/scilpy/dwi/operations.py index 7d6546075..29249e062 100644 --- a/scilpy/dwi/operations.py +++ b/scilpy/dwi/operations.py @@ -1,6 +1,12 @@ import logging +import math +import pprint + import numpy as np +from scilpy.gradients.bvec_bval_tools import identify_shells, \ + round_bvals_to_shell, DEFAULT_B0_THRESHOLD + def apply_bias_field(dwi_data, bias_field_data, mask_data): """ @@ -131,3 +137,84 @@ def compute_dwi_attenuation(dwi_weights: np.ndarray, b0: np.ndarray): dwi_attenuation[np.logical_not(np.isfinite(dwi_attenuation))] = 0. return dwi_attenuation + + +def detect_volume_outliers(data, bvecs, bvals, std_scale, verbose, + b0_thr=DEFAULT_B0_THRESHOLD): + """ + Parameters + ---------- + data: np.ndarray + The dwi data. + bvecs: np.ndarray + The bvecs + bvals: np.array + The b-values vector. + std_scale: float + How many deviation from the mean are required to be considered an + outlier. + verbose: bool + If True, print even more stuff. + b0_thr: float + Value below which b-values are considered as b0. + """ + results_dict = {} + shells_to_extract = identify_shells(bvals, b0_thr, sort=True)[0] + bvals = round_bvals_to_shell(bvals, shells_to_extract) + for bval in shells_to_extract[shells_to_extract > b0_thr]: + shell_idx = np.where(bvals == bval)[0] + shell = bvecs[shell_idx] + results_dict[bval] = np.ones((len(shell), 3)) * -1 + for i, vec in enumerate(shell): + if np.linalg.norm(vec) < 0.001: + continue + + dot_product = np.clip(np.tensordot(shell, vec, axes=1), -1, 1) + angle = np.arccos(dot_product) * 180 / math.pi + angle[np.isnan(angle)] = 0 + idx = np.argpartition(angle, 4).tolist() + idx.remove(i) + + avg_angle = np.average(angle[idx[:3]]) + corr = np.corrcoef([data[..., shell_idx[i]].ravel(), + data[..., shell_idx[idx[0]]].ravel(), + data[..., shell_idx[idx[1]]].ravel(), + data[..., shell_idx[idx[2]]].ravel()]) + results_dict[bval][i] = [shell_idx[i], avg_angle, + np.average(corr[0, 1:])] + + for key in results_dict.keys(): + avg_angle = np.round(np.average(results_dict[key][:, 1]), 4) + std_angle = np.round(np.std(results_dict[key][:, 1]), 4) + + avg_corr = np.round(np.average(results_dict[key][:, 2]), 4) + std_corr = np.round(np.std(results_dict[key][:, 2]), 4) + + outliers_angle = np.argwhere( + results_dict[key][:, 1] < avg_angle - (std_scale * std_angle)) + outliers_corr = np.argwhere( + results_dict[key][:, 2] < avg_corr - (std_scale * std_corr)) + + print('Results for shell {} with {} directions:' + .format(key, len(results_dict[key]))) + print('AVG and STD of angles: {} +/- {}' + .format(avg_angle, std_angle)) + print('AVG and STD of correlations: {} +/- {}' + .format(avg_corr, std_corr)) + + if len(outliers_angle) or len(outliers_corr): + print('Possible outliers ({} STD below or above average):' + .format(std_scale)) + print('Outliers based on angle [position (4D), value]') + for i in outliers_angle: + print(results_dict[key][i, :][0][0:2]) + print('Outliers based on correlation [position (4D), value]') + for i in outliers_corr: + print(results_dict[key][i, :][0][0::2]) + else: + print('No outliers detected.') + + if verbose: + print('Shell with b-value {}'.format(key)) + pprint.pprint(results_dict[key]) + print() diff --git a/scilpy/dwi/tests/test_operations.py b/scilpy/dwi/tests/test_operations.py index 91d7ccb02..7ea12e015 100644 --- a/scilpy/dwi/tests/test_operations.py +++ b/scilpy/dwi/tests/test_operations.py @@ -7,3 +7,7 @@ def test_apply_bias_field(): def test_compute_dwi_attenuation(): pass + + +def test_detect_volume_outliers(): + pass diff --git a/scilpy/io/varian_fdf.py b/scilpy/io/varian_fdf.py index 6deb6f546..7736d42a7 100644 --- a/scilpy/io/varian_fdf.py +++ b/scilpy/io/varian_fdf.py @@ -424,7 +424,9 @@ def correct_procpar_intensity(dwi_data, dwi_path, b0_path): P1 = 10**(Ldb/20) """ - + # Not really an io function: does not load and save! But it is the only + # method we have concerning varian_fdf. Kept here, as decided by guru + # Arnaud. dwi_gain = get_gain(dwi_path) b0_gain = get_gain(b0_path) diff --git a/scripts/scil_dwi_detect_volume_outliers.py b/scripts/scil_dwi_detect_volume_outliers.py index f32bbefce..c5af9833c 100755 --- a/scripts/scil_dwi_detect_volume_outliers.py +++ b/scripts/scil_dwi_detect_volume_outliers.py @@ -15,20 +15,16 @@ """ import argparse -import pprint from dipy.io.gradients import read_bvals_bvecs import nibabel as nib -import numpy as np +from scilpy.dwi.operations import detect_volume_outliers from scilpy.io.utils import (assert_inputs_exist, add_force_b0_arg, add_verbose_arg) from scilpy.gradients.bvec_bval_tools import (check_b0_threshold, - identify_shells, - normalize_bvecs, - round_bvals_to_shell) -import math + normalize_bvecs) def _build_arg_parser(): @@ -49,7 +45,7 @@ def _build_arg_parser(): 'diffusion weighting. [%(default)s]') p.add_argument('--std_scale', type=float, default=2.0, help='How many deviation from the mean are required to be ' - 'considered an outliers. [%(default)s]') + 'considered an outlier. [%(default)s]') add_force_b0_arg(p) add_verbose_arg(p) @@ -64,73 +60,14 @@ def main(): assert_inputs_exist(parser, [args.in_dwi, args.in_bval, args.in_bvec]) bvals, bvecs = read_bvals_bvecs(args.in_bval, args.in_bvec) - dwi = nib.load(args.in_dwi) - data = dwi.get_fdata() + data = nib.load(args.in_dwi).get_fdata() b0_thr = check_b0_threshold(args.force_b0_threshold, bvals.min(), args.b0_thr) bvecs = normalize_bvecs(bvecs) - results_dict = {} - shells_to_extract = identify_shells(bvals, b0_thr, sort=True)[0] - bvals = round_bvals_to_shell(bvals, shells_to_extract) - for bval in shells_to_extract[shells_to_extract > args.b0_thr]: - shell_idx = np.where(bvals == bval)[0] - shell = bvecs[shell_idx] - results_dict[bval] = np.ones((len(shell), 3)) * -1 - for i, vec in enumerate(shell): - if np.linalg.norm(vec) < 0.001: - continue - - dot_product = np.clip(np.tensordot(shell, vec, axes=1), -1, 1) - angle = np.arccos(dot_product) * 180 / math.pi - angle[np.isnan(angle)] = 0 - idx = np.argpartition(angle, 4).tolist() - idx.remove(i) - - avg_angle = np.average(angle[idx[:3]]) - corr = np.corrcoef([data[..., shell_idx[i]].ravel(), - data[..., shell_idx[idx[0]]].ravel(), - data[..., shell_idx[idx[1]]].ravel(), - data[..., shell_idx[idx[2]]].ravel()]) - results_dict[bval][i] = [shell_idx[i], avg_angle, - np.average(corr[0, 1:])] - - for key in results_dict.keys(): - avg_angle = np.round(np.average(results_dict[key][:, 1]), 4) - std_angle = np.round(np.std(results_dict[key][:, 1]), 4) - - avg_corr = np.round(np.average(results_dict[key][:, 2]), 4) - std_corr = np.round(np.std(results_dict[key][:, 2]), 4) - - outliers_angle = np.argwhere( - results_dict[key][:, 1] < avg_angle-(args.std_scale*std_angle)) - outliers_corr = np.argwhere( - results_dict[key][:, 2] < avg_corr-(args.std_scale*std_corr)) - - print('Results for shell {} with {} directions:'.format( - key, len(results_dict[key]))) - print('AVG and STD of angles: {} +/- {}'.format( - avg_angle, std_angle)) - print('AVG and STD of correlations: {} +/- {}'.format( - avg_corr, std_corr)) - - if len(outliers_angle) or len(outliers_corr): - print('Possible outliers ({} STD below or above average):'.format( - args.std_scale)) - print('Outliers based on angle [position (4D), value]') - for i in outliers_angle: - print(results_dict[key][i, :][0][0:2]) - print('Outliers based on correlation [position (4D), value]') - for i in outliers_corr: - print(results_dict[key][i, :][0][0::2]) - else: - print('No outliers detected.') - - if args.verbose: - print('Shell with b-value {}'.format(key)) - pprint.pprint(results_dict[key]) - print() + detect_volume_outliers(data, bvecs, bvals, args.std_scale, + args.verbose, b0_thr) if __name__ == "__main__": From a7bb40f485dc83c92cb5ad3ec44f7a35b680e697 Mon Sep 17 00:00:00 2001 From: AlexVCaron Date: Thu, 15 Feb 2024 14:04:48 -0500 Subject: [PATCH 76/97] use master as ref for PRs --- .github/workflows/test.yml | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 7820777ff..4410cd296 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -53,14 +53,7 @@ jobs: - name: Run tests run: | export C_INCLUDE_PATH=$pythonLocation/include/python${{ steps.python-selector.outputs.python-version }}:$C_INCLUDE_PATH - pytest --cov-report term-missing:skip-covered scripts/tests/test_bids_validate.py - - - name: Upload test reports and coverage to artifacts before codecov - uses: actions/upload-artifact@v4.3.1 - with: - name: test-reports-before - path: | - .test_reports/* + pytest --cov-report term-missing:skip-covered - name: Upload coverage reports to Codecov uses: codecov/codecov-action@v4 @@ -72,6 +65,7 @@ jobs: directory: .test_reports/ fail_ci_if_error: true plugin: pycoverage + commit_parent: ${{ github.event.pull_request.base.sha || github.event.before }} - name: Upload test reports and coverage to artifacts uses: actions/upload-artifact@v4.3.1 From 6ae0976636b493d3ed9297fb11955e4e41440a17 Mon Sep 17 00:00:00 2001 From: AlexVCaron Date: Thu, 15 Feb 2024 14:49:58 -0500 Subject: [PATCH 77/97] remove force base, as it breaks PR coverage counters --- .github/workflows/test.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 4410cd296..ff4f3bf28 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -65,7 +65,6 @@ jobs: directory: .test_reports/ fail_ci_if_error: true plugin: pycoverage - commit_parent: ${{ github.event.pull_request.base.sha || github.event.before }} - name: Upload test reports and coverage to artifacts uses: actions/upload-artifact@v4.3.1 From 773ed7f098a78d86a231daf7b09d1889d53a2cf2 Mon Sep 17 00:00:00 2001 From: AlexVCaron Date: Fri, 16 Feb 2024 16:14:25 -0500 Subject: [PATCH 78/97] update coveragerc --- .coveragerc | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/.coveragerc b/.coveragerc index 9f2015df3..21b501e1c 100644 --- a/.coveragerc +++ b/.coveragerc @@ -1,11 +1,11 @@ [run] branch = True concurrency = multiprocessing -data_file = .test_reports/.coverage -relative_files = True -source_pkgs = +data_file = .coverage +source_pkgs = scilpy scripts +relative_files = True omit = scripts/tests/*.py scilpy/tests/**/*.py @@ -16,6 +16,11 @@ omit = [report] skip_empty = True +skip_covered = True [html] title = Scilpy Coverage Report +directory = .test_reports/coverage.html + +[xml] +output = .test_reports/coverage.xml From d56551ee7f62033f2f9eae7209201480fe270716 Mon Sep 17 00:00:00 2001 From: AlexVCaron Date: Fri, 16 Feb 2024 16:14:56 -0500 Subject: [PATCH 79/97] simplify pytest.ini --- pytest.ini | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/pytest.ini b/pytest.ini index 5a7cbe175..b4b9224d1 100644 --- a/pytest.ini +++ b/pytest.ini @@ -33,8 +33,7 @@ junit_logging = out-err addopts = --html=.test_reports/pytest.html - --cov-report=html:.test_reports/coverage.html --junit-xml=.test_reports/junit.xml - --cov-report=xml:.test_reports/coverage.xml - --cov=scilpy/ - --cov=scripts/ + --cov + --cov-report html + --cov-report xml From 543d1c264b08b5567bf3d82aa5ddd65435706e03 Mon Sep 17 00:00:00 2001 From: AlexVCaron Date: Fri, 16 Feb 2024 16:15:37 -0500 Subject: [PATCH 80/97] tests should now run --- .github/workflows/test.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index ff4f3bf28..ded929612 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -62,7 +62,6 @@ jobs: flags: unittests name: scilpy-unittests-${{ github.run_id }} verbose: true - directory: .test_reports/ fail_ci_if_error: true plugin: pycoverage From d0410dd9b7021d198a5861b3863aeb6cce92f18a Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 10:19:08 -0500 Subject: [PATCH 81/97] remove unused StatefulTractogram --- scripts/scil_tractogram_project_map_to_streamlines.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index 18d0a64ed..0f322b429 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -13,7 +13,7 @@ import argparse import logging -from dipy.io.streamline import save_tractogram, StatefulTractogram +from dipy.io.streamline import save_tractogram from scilpy.io.image import load_img from scilpy.io.streamlines import load_tractogram_with_reference From ced021733b06e935cd4de60d46556603ee9dc0c5 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 10:19:59 -0500 Subject: [PATCH 82/97] removed extra spaces --- scripts/scil_tractogram_project_map_to_streamlines.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index 0f322b429..5db45c336 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -51,8 +51,8 @@ def _build_arg_parser(): p.add_argument('--endpoints_only', action='store_true', help='If set, will only project the map onto the \n' 'endpoints of the streamlines (all other values along \n' - ' streamlines will be NaN). If not set, will project \n' - ' the map onto all points of the streamlines.') + 'streamlines will be NaN). If not set, will project \n' + 'the map onto all points of the streamlines.') p.add_argument('--overwrite_data', action='store_true', default=False, help='If set, will overwrite data_per_point in the ' 'output tractogram, otherwise previous data will be ' From 0c8e07b057d04f9ff950549e30011e17d3dc32a2 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 10:20:25 -0500 Subject: [PATCH 83/97] removed default false --- scripts/scil_tractogram_project_map_to_streamlines.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index 5db45c336..c41e55da7 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -53,7 +53,7 @@ def _build_arg_parser(): 'endpoints of the streamlines (all other values along \n' 'streamlines will be NaN). If not set, will project \n' 'the map onto all points of the streamlines.') - p.add_argument('--overwrite_data', action='store_true', default=False, + p.add_argument('--overwrite_data', action='store_true', help='If set, will overwrite data_per_point in the ' 'output tractogram, otherwise previous data will be ' 'preserved in the output tractogram.') From e401f55886b59b95dbbd1ea4c931e4d691f6d1f1 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 10:32:43 -0500 Subject: [PATCH 84/97] changed to nib.load --- ...l_tractogram_project_map_to_streamlines.py | 24 ++++++++++--------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index c41e55da7..9a5eeecfa 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -13,9 +13,9 @@ import argparse import logging +import nibabel as nib from dipy.io.streamline import save_tractogram -from scilpy.io.image import load_img from scilpy.io.streamlines import load_tractogram_with_reference from scilpy.io.utils import (add_overwrite_arg, add_reference_arg, @@ -53,10 +53,15 @@ def _build_arg_parser(): 'endpoints of the streamlines (all other values along \n' 'streamlines will be NaN). If not set, will project \n' 'the map onto all points of the streamlines.') - p.add_argument('--overwrite_data', action='store_true', - help='If set, will overwrite data_per_point in the ' - 'output tractogram, otherwise previous data will be ' - 'preserved in the output tractogram.') + + p.add_argument('--keep_all_dpp', action='store_true', + help='If set, previous data_per_point will be preserved ' + 'in the output tractogram. Else, only --out_dpp_name ' + 'keys will be saved.') + p.add_argument('--overwrite_dpp', action='store_true', default=False, + help='If set, if --keep_all_dpp is set and some --out_dpp_name ' + 'keys already existed in your data_per_point, allow overwriting ' + 'old data_per_point.') add_reference_arg(p) add_overwrite_arg(p) @@ -102,7 +107,7 @@ def main(): data_per_point = {} for fmap, dpp_name in zip(args.in_maps, args.out_dpp_name): logging.debug("Loading the map...") - map_img, map_dtype = load_img(fmap) + map_img = nib.load(fmap) map_data = map_img.get_fdata(caching='unchanged', dtype=float) map_res = map_img.header.get_zooms()[:3] @@ -126,11 +131,8 @@ def main(): out_sft = sft.from_sft(sft.streamlines, sft, data_per_point=data_per_point) else: - old_data_per_point = sft.data_per_point - for dpp_name in data_per_point: - old_data_per_point[dpp_name] = data_per_point[dpp_name] - out_sft = sft.from_sft(sft.streamlines, sft, - data_per_point=old_data_per_point) + sft.data_per_point.update(data_per_point) + out_sft = sft save_tractogram(out_sft, args.out_tractogram) From 0f0abcce00ce25aa5989b2c865c8c819cffb2f0d Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 10:57:42 -0500 Subject: [PATCH 85/97] added dual keep_dpp and ovewrite args --- scripts/scil_tractogram_project_map_to_streamlines.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index 9a5eeecfa..9acea7f9f 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -58,7 +58,7 @@ def _build_arg_parser(): help='If set, previous data_per_point will be preserved ' 'in the output tractogram. Else, only --out_dpp_name ' 'keys will be saved.') - p.add_argument('--overwrite_dpp', action='store_true', default=False, + p.add_argument('--overwrite_dpp', action='store_true', help='If set, if --keep_all_dpp is set and some --out_dpp_name ' 'keys already existed in your data_per_point, allow overwriting ' 'old data_per_point.') @@ -127,12 +127,12 @@ def main(): data_per_point[dpp_name] = streamline_data - if args.overwrite_data: - out_sft = sft.from_sft(sft.streamlines, sft, - data_per_point=data_per_point) - else: + if args.keep_all_dpp: sft.data_per_point.update(data_per_point) out_sft = sft + else: + out_sft = sft.from_sft(sft.streamlines, sft, + data_per_point=data_per_point) save_tractogram(out_sft, args.out_tractogram) From 3fed0f1535f4f20729e89a9f8e6dea39f78b5f40 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 11:01:22 -0500 Subject: [PATCH 86/97] changed map to map_volume --- scilpy/tractograms/dps_and_dpp_management.py | 22 +++++++++---------- ...l_tractogram_project_map_to_streamlines.py | 4 ++-- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/scilpy/tractograms/dps_and_dpp_management.py b/scilpy/tractograms/dps_and_dpp_management.py index d919e4788..92cb3fb64 100644 --- a/scilpy/tractograms/dps_and_dpp_management.py +++ b/scilpy/tractograms/dps_and_dpp_management.py @@ -2,7 +2,7 @@ import numpy as np -def project_map_to_streamlines(sft, map, endpoints_only=False): +def project_map_to_streamlines(sft, map_volume, endpoints_only=False): """ Projects a map onto the points of streamlines. @@ -10,36 +10,36 @@ def project_map_to_streamlines(sft, map, endpoints_only=False): ---------- sft: StatefulTractogram Input tractogram. - map: DataVolume + map_volume: DataVolume Input map. Optional: --------- endpoints_only: bool - If True, will only project the map onto the endpoints of the + If True, will only project the map_volume onto the endpoints of the streamlines (all values along streamlines set to zero). If False, - will project the map onto all points of the streamlines. + will project the map_volume onto all points of the streamlines. Returns ------- streamline_data: - map projected to each point of the streamlines. + map_volume projected to each point of the streamlines. """ - if len(map.data.shape) == 4: - dimension = map.data.shape[3] + if len(map_volume.data.shape) == 4: + dimension = map_volume.data.shape[3] else: dimension = 1 streamline_data = [] if endpoints_only: for s in sft.streamlines: - p1_data = map.get_value_at_coordinate( + p1_data = map_volume.get_value_at_coordinate( s[0][0], s[0][1], s[0][2], space=sft.space, origin=sft.origin) - p2_data = map.get_value_at_coordinate( + p2_data = map_volume.get_value_at_coordinate( s[-1][0], s[-1][1], s[-1][2], space=sft.space, origin=sft.origin) - thisstreamline_data = [] + if dimension == 1: thisstreamline_data = np.ones((len(s), 1)) * np.nan else: @@ -57,7 +57,7 @@ def project_map_to_streamlines(sft, map, endpoints_only=False): for s in sft.streamlines: thisstreamline_data = [] for p in s: - thisstreamline_data.append(map.get_value_at_coordinate( + thisstreamline_data.append(map_volume.get_value_at_coordinate( p[0], p[1], p[2], space=sft.space, origin=sft.origin)) streamline_data.append( diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index 9acea7f9f..74468ccb0 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -116,11 +116,11 @@ def main(): else: interp = "nearest" - map = DataVolume(map_data, map_res, interp) + map_volume = DataVolume(map_data, map_res, interp) logging.debug("Projecting map onto streamlines") streamline_data = project_map_to_streamlines( - sft, map, + sft, map_volume, endpoints_only=args.endpoints_only) logging.debug("Saving the tractogram...") From 691c620c7e6fbed10dc37f5effa81bbba8de54b8 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 11:06:25 -0500 Subject: [PATCH 87/97] updated logging in project_map --- scripts/scil_tractogram_project_map_to_streamlines.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index 74468ccb0..7e8ad326c 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -77,7 +77,10 @@ def main(): assert_outputs_exist(parser, args, [args.out_tractogram]) - logging.debug("Loading the tractogram...") + if args.verbose: + logging.getLogger().setLevel(logging.INFO) + + logging.info("Loading the tractogram...") sft = load_tractogram_with_reference(parser, args, args.in_tractogram) sft.to_voxmm() sft.to_corner() @@ -106,7 +109,7 @@ def main(): data_per_point = {} for fmap, dpp_name in zip(args.in_maps, args.out_dpp_name): - logging.debug("Loading the map...") + logging.info("Loading the map...") map_img = nib.load(fmap) map_data = map_img.get_fdata(caching='unchanged', dtype=float) map_res = map_img.header.get_zooms()[:3] @@ -118,12 +121,12 @@ def main(): map_volume = DataVolume(map_data, map_res, interp) - logging.debug("Projecting map onto streamlines") + logging.info("Projecting map onto streamlines") streamline_data = project_map_to_streamlines( sft, map_volume, endpoints_only=args.endpoints_only) - logging.debug("Saving the tractogram...") + logging.info("Saving the tractogram...") data_per_point[dpp_name] = streamline_data From 172e821f05f59689673a25b4d47290003540e2fb Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 11:09:20 -0500 Subject: [PATCH 88/97] remove if statement --- scilpy/tractograms/dps_and_dpp_management.py | 6 +----- scripts/scil_tractogram_project_map_to_streamlines.py | 5 +++++ 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/scilpy/tractograms/dps_and_dpp_management.py b/scilpy/tractograms/dps_and_dpp_management.py index 92cb3fb64..4fcdd8e21 100644 --- a/scilpy/tractograms/dps_and_dpp_management.py +++ b/scilpy/tractograms/dps_and_dpp_management.py @@ -40,11 +40,7 @@ def project_map_to_streamlines(sft, map_volume, endpoints_only=False): s[-1][0], s[-1][1], s[-1][2], space=sft.space, origin=sft.origin) - if dimension == 1: - thisstreamline_data = np.ones((len(s), 1)) * np.nan - else: - thisstreamline_data = np.ones( - (len(s), p1_data.shape[0])) * np.nan + thisstreamline_data = np.ones((len(s), dimension)) * np.nan thisstreamline_data[0] = p1_data thisstreamline_data[-1] = p2_data diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index 7e8ad326c..02fa59051 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -137,6 +137,11 @@ def main(): out_sft = sft.from_sft(sft.streamlines, sft, data_per_point=data_per_point) + print("New data_per_point keys are: ") + for key in args.out_dpp_name: + print(" - {} with shape per point {}" + .format(key, out_sft.data_per_point[key][0].shape[1:])) + save_tractogram(out_sft, args.out_tractogram) From d4c4f572fe03e783d3f7cfe03f1158c43224f767 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 11:23:18 -0500 Subject: [PATCH 89/97] fixed comments in dpp_math --- scripts/scil_tractogram_dpp_math.py | 58 ++++++++++++++--------------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/scripts/scil_tractogram_dpp_math.py b/scripts/scil_tractogram_dpp_math.py index 12034bf8f..4cb7852fd 100755 --- a/scripts/scil_tractogram_dpp_math.py +++ b/scripts/scil_tractogram_dpp_math.py @@ -2,24 +2,24 @@ # -*- coding: utf-8 -*- """ -Performs an operation on data per point from input streamlines. - -Two modes of operation are supported: data per streamline (dps) and data per -point (dpp). - -In dps mode, the operation is performed across all dimensions of the data -resulting in a single value per streamline. - -In dpp mode the operation is performed on each point separately, resulting in -a single value per point. - -If endpoints_only and dpp mode is set operation will only +Performs an operation on data per point (dpp) from input streamlines. + +Although the input data always comes from the dpp, the output can be either +a dpp or a data_per_streamline (dps), depending on the chosen options. +Two modes of operation are supported: dpp and dps. + - In dps mode, the operation is performed on dpp across the dimension of + the streamlines resulting in a single value (or array in the 4D case) + per streamline, stored as dps. + - In dpp mode, the operation is performed on each point separately, resulting in + a new dpp. + +If endpoints_only and dpp mode is set the operation will only be calculated at the streamline endpoints the rest of the values along the streamline will be NaN If endpoints_only and dps mode is set operation will be calculated across the data at the endpoints and stored as a -single value per streamline. +single value (or array in the 4D case) per streamline. Endpoint only operation: correlation: correlation calculated between arrays extracted from @@ -53,31 +53,31 @@ def _build_arg_parser(): p.add_argument('operation', metavar='OPERATION', choices=['mean', 'sum', 'min', 'max', 'correlation'], - help='The type of operation to be performed on the ' - 'streamlines. Must\nbe one of the following: ' + help='The type of operation to be performed on the \n' + 'streamlines. Must\nbe one of the following: \n' '%(choices)s.') p.add_argument('dpp_or_dps', metavar='DPP_OR_DPS', choices=['dpp', 'dps'], - help='Set to dps if the operation is to be performed ' - 'across all dimensions resulting in a single value per ' - 'streamline. Set to dpp if the operation is to be ' - 'performed on each point separately resulting in a single ' + help='Set to dps if the operation is to be performed \n' + 'across all dimensions resulting in a single value per \n' + 'streamline. Set to dpp if the operation is to be \n' + 'performed on each point separately resulting in a single \n' 'value per point.') p.add_argument('in_tractogram', metavar='INPUT_FILE', - help='Input tractogram containing streamlines and ' + help='Input tractogram containing streamlines and \n' 'metadata.') p.add_argument('--in_dpp_name', nargs='+', required=True, - help='Name or list of names of the data_per_point for ' - 'operation to be performed on. If more than one dpp ' - 'is selected, the same operation will be applied ' + help='Name or list of names of the data_per_point for \n' + 'operation to be performed on. If more than one dpp \n' + 'is selected, the same operation will be applied \n' 'separately to each one.') p.add_argument('--out_name', nargs='+', required=True, - help='Name of the resulting data_per_point or ' - 'data_per_streamline to be saved in the output ' - 'tractogram. If more than one --in_dpp_name was used, ' + help='Name of the resulting data_per_point or \n' + 'data_per_streamline to be saved in the output \n' + 'tractogram. If more than one --in_dpp_name was used, \n' 'enter the same number of --out_name values.') p.add_argument('out_tractogram', metavar='OUTPUT_FILE', - help='The file where the remaining streamlines ' + help='The file where the remaining streamlines \n' 'are saved.') p.add_argument('--endpoints_only', action='store_true', default=False, @@ -85,8 +85,8 @@ def _build_arg_parser(): 'If not set, will perform operation on all streamline \n' 'points.') p.add_argument('--overwrite_data', action='store_true', default=False, - help='If set, will overwrite the data_per_point or ' - 'data_per_streamline in the output tractogram, otherwise ' + help='If set, will overwrite the data_per_point or \n' + 'data_per_streamline in the output tractogram, otherwise \n' 'previous data will be preserved in the output tractogram.') add_reference_arg(p) From 0a4966a17086299304e324b1d0c6c5defcc8d990 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 11:34:55 -0500 Subject: [PATCH 90/97] updated to mode required argument --- scripts/scil_tractogram_dpp_math.py | 9 ++++----- scripts/scil_tractogram_project_map_to_streamlines.py | 2 +- scripts/tests/test_tractogram_dpp_math.py | 4 ++-- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/scripts/scil_tractogram_dpp_math.py b/scripts/scil_tractogram_dpp_math.py index 4cb7852fd..fe01a0380 100755 --- a/scripts/scil_tractogram_dpp_math.py +++ b/scripts/scil_tractogram_dpp_math.py @@ -56,8 +56,7 @@ def _build_arg_parser(): help='The type of operation to be performed on the \n' 'streamlines. Must\nbe one of the following: \n' '%(choices)s.') - p.add_argument('dpp_or_dps', metavar='DPP_OR_DPS', - choices=['dpp', 'dps'], + p.add_argument('--mode', required=True, choices=['dpp', 'dps'], help='Set to dps if the operation is to be performed \n' 'across all dimensions resulting in a single value per \n' 'streamline. Set to dpp if the operation is to be \n' @@ -138,7 +137,7 @@ def main(): 'point. Exiting.') return - if args.operation == 'correlation' and args.dpp_or_dps == 'dpp': + if args.operation == 'correlation' and args.mode == 'dpp': logging.info('Correlation operation requires dps mode. Exiting.') return @@ -161,7 +160,7 @@ def main(): args.operation, sft, in_dpp_name) data_per_streamline[out_name] = new_dps - elif args.dpp_or_dps == 'dpp': + elif args.mode == 'dpp': # Results in new data per point logging.info( 'Performing {} on data from each streamine point ' @@ -170,7 +169,7 @@ def main(): new_dpp = perform_streamline_operation_per_point( args.operation, sft, in_dpp_name, args.endpoints_only) data_per_point[out_name] = new_dpp - elif args.dpp_or_dps == 'dps': + elif args.mode == 'dps': # Results in new data per streamline logging.info( 'Performing {} across each streamline and saving resulting ' diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index 02fa59051..e7578a80b 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -99,7 +99,7 @@ def main(): parser.error('The output names (out_dpp_names) must be unique.') # Check to see if the output names already exist in the input tractogram - if not args.overwrite_data: + if not args.overwrite_dpp: for out_dpp_name in args.out_dpp_name: if out_dpp_name in sft.data_per_point: logging.info('out_name {} already exists in input tractogram. ' diff --git a/scripts/tests/test_tractogram_dpp_math.py b/scripts/tests/test_tractogram_dpp_math.py index ef89d6eca..69387b53d 100644 --- a/scripts/tests/test_tractogram_dpp_math.py +++ b/scripts/tests/test_tractogram_dpp_math.py @@ -34,9 +34,9 @@ def test_execution_tractogram_point_math_mean_3D_defaults(script_runner): ret = script_runner.run('scil_tractogram_dpp_math.py', 'mean', - 'dps', t1_on_bundle, 't1_mean_on_streamlines.trk', + '--mode', 'dps', '--in_dpp_name', 't1', '--out_name', 't1_mean') @@ -59,9 +59,9 @@ def test_execution_tractogram_point_math_mean_4D_correlation(script_runner): ret = script_runner.run('scil_tractogram_dpp_math.py', 'correlation', - 'dps', fodf_on_bundle, 'fodf_correlation_on_streamlines.trk', + '--mode', 'dpp', '--in_dpp_name', 'fodf', '--out_name', 'fodf_correlation') From efa733a4af7e48d458416ba31cb27ec34f2b787a Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 11:53:17 -0500 Subject: [PATCH 91/97] added keep_all_dpp_dps option --- scripts/scil_tractogram_dpp_math.py | 46 +++++++++++++---------------- 1 file changed, 20 insertions(+), 26 deletions(-) diff --git a/scripts/scil_tractogram_dpp_math.py b/scripts/scil_tractogram_dpp_math.py index fe01a0380..5784e8e04 100755 --- a/scripts/scil_tractogram_dpp_math.py +++ b/scripts/scil_tractogram_dpp_math.py @@ -56,15 +56,15 @@ def _build_arg_parser(): help='The type of operation to be performed on the \n' 'streamlines. Must\nbe one of the following: \n' '%(choices)s.') + p.add_argument('in_tractogram', metavar='INPUT_FILE', + help='Input tractogram containing streamlines and \n' + 'metadata.') p.add_argument('--mode', required=True, choices=['dpp', 'dps'], help='Set to dps if the operation is to be performed \n' 'across all dimensions resulting in a single value per \n' 'streamline. Set to dpp if the operation is to be \n' 'performed on each point separately resulting in a single \n' 'value per point.') - p.add_argument('in_tractogram', metavar='INPUT_FILE', - help='Input tractogram containing streamlines and \n' - 'metadata.') p.add_argument('--in_dpp_name', nargs='+', required=True, help='Name or list of names of the data_per_point for \n' 'operation to be performed on. If more than one dpp \n' @@ -78,15 +78,19 @@ def _build_arg_parser(): p.add_argument('out_tractogram', metavar='OUTPUT_FILE', help='The file where the remaining streamlines \n' 'are saved.') - p.add_argument('--endpoints_only', action='store_true', default=False, help='If set, will only perform operation on endpoints \n' 'If not set, will perform operation on all streamline \n' 'points.') - p.add_argument('--overwrite_data', action='store_true', default=False, - help='If set, will overwrite the data_per_point or \n' - 'data_per_streamline in the output tractogram, otherwise \n' - 'previous data will be preserved in the output tractogram.') + p.add_argument('--keep_all_dpp_dps', action='store_true', + help='If set, previous data_per_point will be preserved \n' + 'in the output tractogram. Else, only --out_dpp_name \n' + 'keys will be saved.') + p.add_argument('--overwrite_dpp_dps', action='store_true', + help='If set, if --keep_all_dpp_dps is set and some \n' + '--out_dpp_name keys already existed in your \n' + ' data_per_point or data_per_streamline, allow \n' + ' overwriting old data_per_point.') add_reference_arg(p) add_verbose_arg(p) @@ -141,10 +145,10 @@ def main(): logging.info('Correlation operation requires dps mode. Exiting.') return - if not args.overwrite_data: + if not args.overwrite_dpp_dps: if in_dpp_name in args.out_name: logging.info('out_name {} already exists in input tractogram. ' - 'Set overwrite_data or choose a different ' + 'Set overwrite_dpp_dps or choose a different ' 'out_name. Exiting.'.format(in_dpp_name)) return @@ -178,25 +182,15 @@ def main(): args.operation, sft, in_dpp_name, args.endpoints_only) data_per_streamline[out_name] = new_data_per_streamline - if args.overwrite_data: + if args.keep_all_dpp_dps: + sft.data_per_streamline.update(data_per_streamline) + sft.data_per_point.update(data_per_point) + + new_sft = sft + else: new_sft = sft.from_sft(sft.streamlines, sft, data_per_point=data_per_point, data_per_streamline=data_per_streamline) - else: - old_data_per_streamline = sft.data_per_streamline - old_data_per_point = sft.data_per_point - - if data_per_point is not None: - for key, value in data_per_point.items(): - old_data_per_point[key] = value - - if data_per_streamline is not None: - for key, value in data_per_streamline.items(): - old_data_per_streamline[key] = value - - new_sft = sft.from_sft(sft.streamlines, sft, - data_per_point=old_data_per_point, - data_per_streamline=old_data_per_streamline) # Print DPP names logging.info('New data per point names:') From f4b4e394fb8a317f8b4668e03e27077f16b5d6b1 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 12:04:58 -0500 Subject: [PATCH 92/97] moved streamline ops to dps_and_dpp --- scilpy/tractograms/dps_and_dpp_management.py | 141 ++++++++++++++++++ scilpy/tractograms/streamline_operations.py | 143 +------------------ scripts/scil_tractogram_dpp_math.py | 27 ++-- 3 files changed, 160 insertions(+), 151 deletions(-) diff --git a/scilpy/tractograms/dps_and_dpp_management.py b/scilpy/tractograms/dps_and_dpp_management.py index 4fcdd8e21..1f85ef2dc 100644 --- a/scilpy/tractograms/dps_and_dpp_management.py +++ b/scilpy/tractograms/dps_and_dpp_management.py @@ -61,3 +61,144 @@ def project_map_to_streamlines(sft, map_volume, endpoints_only=False): (len(thisstreamline_data), dimension))) return streamline_data + + +def perform_streamline_operation_per_point(op_name, sft, dpp_name='metric', + endpoints_only=False): + """Peforms an operation per point for all streamlines. + + Parameters + ---------- + op_name: str + A callable that takes a list of streamline data per point (4D) and + returns a list of streamline data per point. + sft: StatefulTractogram + The streamlines used in the operation. + dpp_name: str + The name of the data per point to be used in the operation. + endpoints_only: bool + If True, will only perform operation on endpoints + + Returns + ------- + new_sft: StatefulTractogram + sft with data per streamline resulting from the operation. + """ + + # Performing operation + call_op = OPERATIONS[op_name] + if endpoints_only: + new_data_per_point = [] + for s in sft.data_per_point[dpp_name]: + this_data_per_point = np.nan * np.ones((len(s), 1)) + this_data_per_point[0] = call_op(s[0]) + this_data_per_point[-1] = call_op(s[-1]) + new_data_per_point.append( + np.reshape(this_data_per_point, (len(this_data_per_point), 1))) + else: + new_data_per_point = [] + for s in sft.data_per_point[dpp_name]: + this_data_per_point = [] + for p in s: + this_data_per_point.append(call_op(p)) + new_data_per_point.append( + np.reshape(this_data_per_point, (len(this_data_per_point), 1))) + + # Extracting streamlines + return new_data_per_point + + +def perform_operation_per_streamline(op_name, sft, dpp_name='metric', + endpoints_only=False): + """Performs an operation across all data points for each streamline. + + Parameters + ---------- + op_name: str + A callable that takes a list of streamline data per streamline and + returns a list of data per streamline. + sft: StatefulTractogram + The streamlines used in the operation. + dpp_name: str + The name of the data per point to be used in the operation. + endpoints_only: bool + If True, will only perform operation on endpoints + + Returns + ------- + new_sft: StatefulTractogram + sft with data per streamline resulting from the operation. + """ + # Performing operation + call_op = OPERATIONS[op_name] + if endpoints_only: + new_data_per_streamline = [] + for s in sft.data_per_point[dpp_name]: + start = s[0] + end = s[-1] + concat = np.concatenate((start[:], end[:])) + new_data_per_streamline.append(call_op(concat)) + else: + new_data_per_streamline = [] + for s in sft.data_per_point[dpp_name]: + s_np = np.asarray(s) + new_data_per_streamline.append(call_op(s_np)) + + return new_data_per_streamline + + +def perform_pairwise_streamline_operation_on_endpoints(op_name, sft, + dpp_name='metric'): + """Peforms an operation across endpoints for each streamline. + + Parameters + ---------- + op_name: str + A callable that takes a list of streamline data per streamline and + returns a list of data per streamline. + sft: StatefulTractogram + The streamlines used in the operation. + dpp_name: str + The name of the data per point to be used in the operation. + + Returns + ------- + new_sft: StatefulTractogram + sft with data per streamline resulting from the operation. + """ + # Performing operation + call_op = OPERATIONS[op_name] + new_data_per_streamline = [] + for s in sft.data_per_point[dpp_name]: + new_data_per_streamline.append(call_op(s[0], s[-1])[0, 1]) + + return new_data_per_streamline + + +def stream_mean(array): + return np.squeeze(np.mean(array, axis=0)) + + +def stream_sum(array): + return np.squeeze(np.sum(array, axis=0)) + + +def stream_min(array): + return np.squeeze(np.min(array, axis=0)) + + +def stream_max(array): + return np.squeeze(np.max(array, axis=0)) + + +def stream_correlation(array1, array2): + return np.corrcoef(array1, array2) + + +OPERATIONS = { + 'mean': stream_mean, + 'sum': stream_sum, + 'min': stream_min, + 'max': stream_max, + 'correlation': stream_correlation, +} diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index 3a98316e0..f51d56256 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -387,145 +387,4 @@ def smooth_line_spline(streamline, smoothing_parameter, nb_ctrl_points): smoothed_streamline[0] = streamline[0] smoothed_streamline[-1] = streamline[-1] - return smoothed_streamline - - -def perform_streamline_operation_per_point(op_name, sft, dpp_name='metric', - endpoints_only=False): - """Peforms an operation per point for all streamlines. - - Parameters - ---------- - op_name: str - A callable that takes a list of streamline data per point (4D) and - returns a list of streamline data per point. - sft: StatefulTractogram - The streamlines used in the operation. - dpp_name: str - The name of the data per point to be used in the operation. - endpoints_only: bool - If True, will only perform operation on endpoints - - Returns - ------- - new_sft: StatefulTractogram - sft with data per streamline resulting from the operation. - """ - - # Performing operation - call_op = OPERATIONS[op_name] - if endpoints_only: - new_data_per_point = [] - for s in sft.data_per_point[dpp_name]: - this_data_per_point = np.nan * np.ones((len(s), 1)) - this_data_per_point[0] = call_op(s[0]) - this_data_per_point[-1] = call_op(s[-1]) - new_data_per_point.append( - np.reshape(this_data_per_point, (len(this_data_per_point), 1))) - else: - new_data_per_point = [] - for s in sft.data_per_point[dpp_name]: - this_data_per_point = [] - for p in s: - this_data_per_point.append(call_op(p)) - new_data_per_point.append( - np.reshape(this_data_per_point, (len(this_data_per_point), 1))) - - # Extracting streamlines - return new_data_per_point - - -def perform_operation_per_streamline(op_name, sft, dpp_name='metric', - endpoints_only=False): - """Performs an operation across all data points for each streamline. - - Parameters - ---------- - op_name: str - A callable that takes a list of streamline data per streamline and - returns a list of data per streamline. - sft: StatefulTractogram - The streamlines used in the operation. - dpp_name: str - The name of the data per point to be used in the operation. - endpoints_only: bool - If True, will only perform operation on endpoints - - Returns - ------- - new_sft: StatefulTractogram - sft with data per streamline resulting from the operation. - """ - # Performing operation - call_op = OPERATIONS[op_name] - if endpoints_only: - new_data_per_streamline = [] - for s in sft.data_per_point[dpp_name]: - start = s[0] - end = s[-1] - concat = np.concatenate((start[:], end[:])) - new_data_per_streamline.append(call_op(concat)) - else: - new_data_per_streamline = [] - for s in sft.data_per_point[dpp_name]: - s_np = np.asarray(s) - new_data_per_streamline.append(call_op(s_np)) - - return new_data_per_streamline - - -def perform_pairwise_streamline_operation_on_endpoints(op_name, sft, - dpp_name='metric'): - """Peforms an operation across endpoints for each streamline. - - Parameters - ---------- - op_name: str - A callable that takes a list of streamline data per streamline and - returns a list of data per streamline. - sft: StatefulTractogram - The streamlines used in the operation. - dpp_name: str - The name of the data per point to be used in the operation. - - Returns - ------- - new_sft: StatefulTractogram - sft with data per streamline resulting from the operation. - """ - # Performing operation - call_op = OPERATIONS[op_name] - new_data_per_streamline = [] - for s in sft.data_per_point[dpp_name]: - new_data_per_streamline.append(call_op(s[0], s[-1])[0, 1]) - - return new_data_per_streamline - - -def stream_mean(array): - return np.squeeze(np.mean(array, axis=0)) - - -def stream_sum(array): - return np.squeeze(np.sum(array, axis=0)) - - -def stream_min(array): - return np.squeeze(np.min(array, axis=0)) - - -def stream_max(array): - return np.squeeze(np.max(array, axis=0)) - - -def stream_correlation(array1, array2): - return np.corrcoef(array1, array2) - - -OPERATIONS = { - 'mean': stream_mean, - 'sum': stream_sum, - 'min': stream_min, - 'max': stream_max, - 'correlation': stream_correlation, -} + return smoothed_streamline \ No newline at end of file diff --git a/scripts/scil_tractogram_dpp_math.py b/scripts/scil_tractogram_dpp_math.py index 5784e8e04..6d4bb8baf 100755 --- a/scripts/scil_tractogram_dpp_math.py +++ b/scripts/scil_tractogram_dpp_math.py @@ -39,12 +39,11 @@ add_verbose_arg, assert_inputs_exist, assert_outputs_exist) -from scilpy.tractograms.streamline_operations import ( +from scilpy.tractograms.dps_and_dpp_management import ( perform_pairwise_streamline_operation_on_endpoints, perform_streamline_operation_per_point, perform_operation_per_streamline) - def _build_arg_parser(): p = argparse.ArgumentParser( formatter_class=argparse.RawTextHelpFormatter, @@ -134,9 +133,15 @@ def main(): .format(in_dpp_name)) return + if args.dpp_or_dps == 'dpp' and (len(data_shape) == 1 or + data_shape[1] == 1): + logging.warning('dpp from key {} is a single number per point. ' + 'Performing a dpp-mode operation on this data ' + 'will not do anything. Continuing.') + # Check if first data_per_point is multivalued data_shape = sft.data_per_point[in_dpp_name][0].shape - if args.operation == 'correlation' and len(data_shape) == 1: + if args.operation == 'correlation' and len(data_shape) == 2 and data_shape[1] > 1: logging.info('Correlation operation requires multivalued data per ' 'point. Exiting.') return @@ -193,14 +198,18 @@ def main(): data_per_streamline=data_per_streamline) # Print DPP names - logging.info('New data per point names:') - for key in new_sft.data_per_point.keys(): - logging.info(key) + if data_per_point not in [None, {}]: + print("New data_per_point keys are: ") + for key in args.out_name: + print(" - {} with shape per point {}" + .format(key, new_sft.data_per_point[key][0].shape[1:])) # Print DPS names - logging.info('New data per streamline names:') - for key in new_sft.data_per_streamline.keys(): - logging.info(key) + if data_per_streamline not in [None, {}]: + print("New data_per_streamline keys are: ") + for key in args.out_name: + print(" - {} with shape per streamline {}" + .format(key, new_sft.data_per_streamline[key].shape[1:])) # Save the new streamlines (and metadata) logging.info('Saving {} streamlines to {}.'.format(len(new_sft), From e0e70a3e8c94043a784e330ace05f1a741119aca Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 12:15:36 -0500 Subject: [PATCH 93/97] blank tests to dps_dpp module --- .../tests/test_dps_and_dpp_management.py | 15 +++++++++++++++ .../tests/test_streamline_operations.py | 15 --------------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/scilpy/tractograms/tests/test_dps_and_dpp_management.py b/scilpy/tractograms/tests/test_dps_and_dpp_management.py index fb93d736a..1fc2a7e88 100644 --- a/scilpy/tractograms/tests/test_dps_and_dpp_management.py +++ b/scilpy/tractograms/tests/test_dps_and_dpp_management.py @@ -1,3 +1,18 @@ def test_project_map_to_streamlines(): # toDo pass + + +def test_perform_streamline_operation_per_point(): + # toDo + pass + + +def test_perform_operation_per_streamline(): + # toDo + pass + + +def test_perform_streamline_operation_on_endpoints(): + # toDo + pass diff --git a/scilpy/tractograms/tests/test_streamline_operations.py b/scilpy/tractograms/tests/test_streamline_operations.py index f8e0d4e8a..faeb4e57c 100644 --- a/scilpy/tractograms/tests/test_streamline_operations.py +++ b/scilpy/tractograms/tests/test_streamline_operations.py @@ -367,18 +367,3 @@ def test_smooth_line_spline(): dist_2 = np.linalg.norm(noisy_streamline - smoothed_streamline) assert dist_1 < dist_2 - - -def test_perform_streamline_operation_per_point(): - # toDo - pass - - -def test_perform_operation_per_streamline(): - # toDo - pass - - -def test_perform_streamline_operation_on_endpoints(): - # toDo - pass From 6d752e5640dfc0b20708b0a95a2756f1f0c762ee Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 12:23:54 -0500 Subject: [PATCH 94/97] pep8 --- scilpy/tractograms/streamline_operations.py | 2 +- scripts/scil_tractogram_dpp_math.py | 16 +++++++++------- ...scil_tractogram_project_map_to_streamlines.py | 10 +++++----- 3 files changed, 15 insertions(+), 13 deletions(-) diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index f51d56256..7df3a396c 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -387,4 +387,4 @@ def smooth_line_spline(streamline, smoothing_parameter, nb_ctrl_points): smoothed_streamline[0] = streamline[0] smoothed_streamline[-1] = streamline[-1] - return smoothed_streamline \ No newline at end of file + return smoothed_streamline diff --git a/scripts/scil_tractogram_dpp_math.py b/scripts/scil_tractogram_dpp_math.py index 6d4bb8baf..f391976ef 100755 --- a/scripts/scil_tractogram_dpp_math.py +++ b/scripts/scil_tractogram_dpp_math.py @@ -10,8 +10,8 @@ - In dps mode, the operation is performed on dpp across the dimension of the streamlines resulting in a single value (or array in the 4D case) per streamline, stored as dps. - - In dpp mode, the operation is performed on each point separately, resulting in - a new dpp. + - In dpp mode, the operation is performed on each point separately, + resulting in a new dpp. If endpoints_only and dpp mode is set the operation will only be calculated at the streamline endpoints the rest of the @@ -44,6 +44,7 @@ perform_streamline_operation_per_point, perform_operation_per_streamline) + def _build_arg_parser(): p = argparse.ArgumentParser( formatter_class=argparse.RawTextHelpFormatter, @@ -62,8 +63,8 @@ def _build_arg_parser(): help='Set to dps if the operation is to be performed \n' 'across all dimensions resulting in a single value per \n' 'streamline. Set to dpp if the operation is to be \n' - 'performed on each point separately resulting in a single \n' - 'value per point.') + 'performed on each point separately resulting in a \n' + 'single value per point.') p.add_argument('--in_dpp_name', nargs='+', required=True, help='Name or list of names of the data_per_point for \n' 'operation to be performed on. If more than one dpp \n' @@ -141,7 +142,8 @@ def main(): # Check if first data_per_point is multivalued data_shape = sft.data_per_point[in_dpp_name][0].shape - if args.operation == 'correlation' and len(data_shape) == 2 and data_shape[1] > 1: + if args.operation == 'correlation' and len(data_shape) == 2 \ + and data_shape[1] > 1: logging.info('Correlation operation requires multivalued data per ' 'point. Exiting.') return @@ -202,14 +204,14 @@ def main(): print("New data_per_point keys are: ") for key in args.out_name: print(" - {} with shape per point {}" - .format(key, new_sft.data_per_point[key][0].shape[1:])) + .format(key, new_sft.data_per_point[key][0].shape[1:])) # Print DPS names if data_per_streamline not in [None, {}]: print("New data_per_streamline keys are: ") for key in args.out_name: print(" - {} with shape per streamline {}" - .format(key, new_sft.data_per_streamline[key].shape[1:])) + .format(key, new_sft.data_per_streamline[key].shape[1:])) # Save the new streamlines (and metadata) logging.info('Saving {} streamlines to {}.'.format(len(new_sft), diff --git a/scripts/scil_tractogram_project_map_to_streamlines.py b/scripts/scil_tractogram_project_map_to_streamlines.py index e7578a80b..027515781 100755 --- a/scripts/scil_tractogram_project_map_to_streamlines.py +++ b/scripts/scil_tractogram_project_map_to_streamlines.py @@ -55,13 +55,13 @@ def _build_arg_parser(): 'the map onto all points of the streamlines.') p.add_argument('--keep_all_dpp', action='store_true', - help='If set, previous data_per_point will be preserved ' - 'in the output tractogram. Else, only --out_dpp_name ' + help='If set, previous data_per_point will be preserved \n' + 'in the output tractogram. Else, only --out_dpp_name \n' 'keys will be saved.') p.add_argument('--overwrite_dpp', action='store_true', - help='If set, if --keep_all_dpp is set and some --out_dpp_name ' - 'keys already existed in your data_per_point, allow overwriting ' - 'old data_per_point.') + help='If set, if --keep_all_dpp is set and some \n' + '--out_dpp_name keys already existed in your \n' + 'data_per_point, allow overwriting old data_per_point.') add_reference_arg(p) add_overwrite_arg(p) From a21aa16f221ee12576bd9c7c02367531cd2eaae1 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 12:47:02 -0500 Subject: [PATCH 95/97] 4D data stored as shape (#,) fixed checks --- scripts/scil_tractogram_dpp_math.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/scripts/scil_tractogram_dpp_math.py b/scripts/scil_tractogram_dpp_math.py index f391976ef..bdb54faf2 100755 --- a/scripts/scil_tractogram_dpp_math.py +++ b/scripts/scil_tractogram_dpp_math.py @@ -40,9 +40,9 @@ assert_inputs_exist, assert_outputs_exist) from scilpy.tractograms.dps_and_dpp_management import ( - perform_pairwise_streamline_operation_on_endpoints, - perform_streamline_operation_per_point, - perform_operation_per_streamline) + perform_pairwise_streamline_operation_on_endpoints, + perform_streamline_operation_per_point, + perform_operation_per_streamline) def _build_arg_parser(): @@ -134,16 +134,15 @@ def main(): .format(in_dpp_name)) return - if args.dpp_or_dps == 'dpp' and (len(data_shape) == 1 or - data_shape[1] == 1): + # warning if dpp mode and data in single number per point + data_shape = sft.data_per_point[in_dpp_name][0].shape + if args.mode == 'dpp' and data_shape[0] == 1: logging.warning('dpp from key {} is a single number per point. ' 'Performing a dpp-mode operation on this data ' 'will not do anything. Continuing.') # Check if first data_per_point is multivalued - data_shape = sft.data_per_point[in_dpp_name][0].shape - if args.operation == 'correlation' and len(data_shape) == 2 \ - and data_shape[1] > 1: + if args.operation == 'correlation' and data_shape > 1: logging.info('Correlation operation requires multivalued data per ' 'point. Exiting.') return From 8db77a91ecb8fe5d7913fecb2562fcc1341aede1 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 12:47:48 -0500 Subject: [PATCH 96/97] fixed dimension on tuple --- scripts/scil_tractogram_dpp_math.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/scil_tractogram_dpp_math.py b/scripts/scil_tractogram_dpp_math.py index bdb54faf2..80fbd091f 100755 --- a/scripts/scil_tractogram_dpp_math.py +++ b/scripts/scil_tractogram_dpp_math.py @@ -142,7 +142,7 @@ def main(): 'will not do anything. Continuing.') # Check if first data_per_point is multivalued - if args.operation == 'correlation' and data_shape > 1: + if args.operation == 'correlation' and data_shape[0] > 1: logging.info('Correlation operation requires multivalued data per ' 'point. Exiting.') return From 47f287bd224dbebf11e33ac063425d3772d0f796 Mon Sep 17 00:00:00 2001 From: grahamlittlephd Date: Mon, 19 Feb 2024 12:57:57 -0500 Subject: [PATCH 97/97] updated dimension on data_shape --- scripts/scil_tractogram_dpp_math.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/scil_tractogram_dpp_math.py b/scripts/scil_tractogram_dpp_math.py index 80fbd091f..bedaace47 100755 --- a/scripts/scil_tractogram_dpp_math.py +++ b/scripts/scil_tractogram_dpp_math.py @@ -142,7 +142,7 @@ def main(): 'will not do anything. Continuing.') # Check if first data_per_point is multivalued - if args.operation == 'correlation' and data_shape[0] > 1: + if args.operation == 'correlation' and data_shape[0] == 1: logging.info('Correlation operation requires multivalued data per ' 'point. Exiting.') return