diff --git a/scilpy/image/labels.py b/scilpy/image/labels.py index 90d73c3c1..4a8702102 100644 --- a/scilpy/image/labels.py +++ b/scilpy/image/labels.py @@ -24,7 +24,7 @@ def get_data_as_labels(in_img): curr_type = in_img.get_data_dtype() if np.issubdtype(curr_type, np.signedinteger) or \ - np.issubdtype(curr_type, np.unsignedinteger): + np.issubdtype(curr_type, np.unsignedinteger): return np.asanyarray(in_img.dataobj).astype(np.uint16) else: basename = os.path.basename(in_img.get_filename()) @@ -300,3 +300,48 @@ def dilate_labels(data, vox_size, distance, nbr_processes, data = data.reshape(img_shape) return data + + +def get_stats_in_label(map_data, label_data, label_lut): + """ + Get statistics about a map for each label in an atlas. + + Parameters + ---------- + map_data: np.ndarray + The map from which to get statistics. + label_data: np.ndarray + The loaded atlas. + label_lut: dict + The loaded label LUT (look-up table). + + Returns + ------- + out_dict: dict + A dict with one key per label name, and its values are the computed + statistics. + """ + (label_indices, label_names) = zip(*label_lut.items()) + + out_dict = {} + for label, name in zip(label_indices, label_names): + label = int(label) + if label != 0: + curr_data = (map_data[np.where(label_data == label)]) + nb_vx_roi = np.count_nonzero(label_data == label) + nb_seed_vx = np.count_nonzero(curr_data) + + if nb_seed_vx != 0: + mean_seed = np.sum(curr_data) / nb_seed_vx + max_seed = np.max(curr_data) + std_seed = np.sqrt(np.mean(abs(curr_data[curr_data != 0] - + mean_seed) ** 2)) + + out_dict[name] = {'ROI-idx': label, + 'ROI-name': str(name), + 'nb-vx-roi': int(nb_vx_roi), + 'nb-vx-seed': int(nb_seed_vx), + 'max': int(max_seed), + 'mean': float(mean_seed), + 'std': float(std_seed)} + return out_dict diff --git a/scilpy/image/tests/test_labels.py b/scilpy/image/tests/test_labels.py index bb7d1129a..e30685580 100644 --- a/scilpy/image/tests/test_labels.py +++ b/scilpy/image/tests/test_labels.py @@ -157,3 +157,8 @@ def test_split_labels(): assert len(out_labels) == 3 assert_equal(np.unique(out_labels[0]), [0, 6]) assert_equal(np.unique(out_labels[1]), [0]) + + +def test_stats_in_labels(): + # toDO. Will need to create a fake LUT. + pass diff --git a/scilpy/image/tests/test_volume_operations.py b/scilpy/image/tests/test_volume_operations.py index 02d52e3b1..2c7e7fb47 100644 --- a/scilpy/image/tests/test_volume_operations.py +++ b/scilpy/image/tests/test_volume_operations.py @@ -119,6 +119,16 @@ def test_compute_snr(): assert np.allclose(snr[0]['snr'], target_val, atol=0.00005) +def test_remove_outliers_ransac(): + # Could test, but uses mainly sklearn. Not testing again. + pass + + +def smooth_to_fwhm(): + # toDo + pass + + def test_resample_volume(): moving3d = np.pad(np.ones((4, 4, 4)), pad_width=1, mode='constant', constant_values=0) diff --git a/scilpy/image/volume_b0_synthesis.py b/scilpy/image/volume_b0_synthesis.py new file mode 100644 index 000000000..303ce8ef7 --- /dev/null +++ b/scilpy/image/volume_b0_synthesis.py @@ -0,0 +1,96 @@ +# -*- coding: utf-8 -*- +import logging +import os +import warnings + +# Disable tensorflow warnings +with warnings.catch_warnings(): + os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3' + warnings.simplefilter("ignore") + from dipy.nn.synb0 import Synb0 + +import numpy as np +from dipy.align.imaffine import AffineMap +from dipy.segment.tissue import TissueClassifierHMRF +from scipy.ndimage import gaussian_filter + +from scilpy.image.volume_operations import register_image + + +def compute_b0_synthesis(t1_data, t1_bet_data, b0_data, b0_bet_data, + template_data, t1_affine, b0_affine, template_affine, + verbose): + """ + Note. Tensorflow is required here, through dipy.Synb0. If not installed, + dipy will raise an error, like: + >> dipy.utils.tripwire.TripWireError: We need package tensorflow_addons for + these functions, but ``import tensorflow_addons`` raised an ImportError. + + Parameters + ---------- + t1_data: np.ndarary + The T1 (wholebrain, with the skull). + t1_bet_data: np.ndarray + The mask. + b0_data: np.ndarray + The b0 (wholebrain, with the skull). + b0_bet_data: np.ndarray + The mask. + template_data: np.ndarray + The template for typical usage. + t1_affine: np.ndarray + The T1's affine. + b0_affine: np.ndarray + The b0's affine. + template_affine: np.ndarray + The template's affine + verbose: bool + Whether to make dipy's Synb0 verbose. + + Returns + ------- + rev_b0: np.ndarray + The synthetized b0. + """ + logging.info('The usage of synthetic b0 is not fully tested.' + 'Be careful when using it.') + + # Crude estimation of the WM mean intensity and normalization + logging.info('Estimating WM mean intensity') + hmrf = TissueClassifierHMRF() + t1_bet_data = gaussian_filter(t1_bet_data, 2) + _, final_segmentation, _ = hmrf.classify(t1_bet_data, 3, 0.25, + tolerance=1e-4, max_iter=5) + avg_wm = np.mean(t1_data[final_segmentation == 3]) + + # Modifying t1 + t1_data /= avg_wm + t1_data *= 110 + + # SyNB0 works only in a standard space, so we need to register the images + logging.info('Registering images') + # Use the BET image for registration + t1_bet_to_b0, t1_bet_to_b0_transform = register_image( + b0_bet_data, b0_affine, t1_bet_data, t1_affine, fine=True) + affine_map = AffineMap(t1_bet_to_b0_transform, + b0_data.shape, b0_affine, + t1_data.shape, t1_affine) + t1_skull_to_b0 = affine_map.transform(t1_data.astype(np.float64)) + + # Then register to MNI (using the BET again) + _, t1_bet_to_b0_to_mni_transform = register_image( + template_data, template_affine, t1_bet_to_b0, b0_affine, fine=True) + affine_map = AffineMap(t1_bet_to_b0_to_mni_transform, + template_data.shape, template_affine, + b0_data.shape, b0_affine) + + # But for prediction, we want the skull + b0_skull_to_mni = affine_map.transform(b0_data.astype(np.float64)) + t1_skull_to_mni = affine_map.transform(t1_skull_to_b0.astype(np.float64)) + + logging.info('Running SyN-B0') + SyNb0 = Synb0(verbose) + rev_b0 = SyNb0.predict(b0_skull_to_mni, t1_skull_to_mni) + rev_b0 = affine_map.transform_inverse(rev_b0.astype(np.float64)) + + return rev_b0 diff --git a/scilpy/image/volume_operations.py b/scilpy/image/volume_operations.py index e8ede37ad..8ec3a2826 100644 --- a/scilpy/image/volume_operations.py +++ b/scilpy/image/volume_operations.py @@ -16,6 +16,7 @@ import numpy as np from numpy import ma from scipy.ndimage import binary_dilation, gaussian_filter +from sklearn import linear_model from scilpy.image.reslice import reslice # Don't use Dipy's reslice. Buggy. from scilpy.io.image import get_data_as_mask @@ -152,6 +153,10 @@ def apply_transform(transfo, reference, moving, elif moving_data.ndim == 4: if isinstance(moving_data[0, 0, 0], np.void): raise ValueError('Does not support TrackVis RGB') + else: + logging.warning('You are applying a transform to a 4D volume. If' + 'it is a DWI volume, make sure to rotate your ' + 'bvecs with scil_gradients_apply_transform.py') affine_map = AffineMap(np.linalg.inv(transfo), dim[0:3], grid2world, @@ -368,6 +373,47 @@ def compute_snr(dwi, bval, bvec, b0_thr, mask, return val +def remove_outliers_ransac(in_data, min_fit, fit_thr, max_iter): + """ + Remove outliers from image using the RANSAC algorithm. + + Parameters + ---------- + in_data: np.ndarray + The input. + min_fit: int + The minimum number of data values required to fit the model. + fit_thr: float + Threshold value for determining when a data point fits a model. + max_iter: int + The maximum number of iterations allowed in the algorithm. + + Returns + ------- + out_data: np.ndarray + Data without outliers. + """ + init_shape = in_data.shape + in_data_flat = in_data.flatten() + in_nzr_ind = np.nonzero(in_data_flat) + in_nzr_val = np.array(in_data_flat[in_nzr_ind]) + + X = in_nzr_ind[0][:, np.newaxis] + model_ransac = linear_model.RANSACRegressor( + base_estimator=linear_model.LinearRegression(), min_samples=min_fit, + residual_threshold=fit_thr, max_trials=max_iter) + model_ransac.fit(X, in_nzr_val) + + outlier_mask = np.logical_not(model_ransac.inlier_mask_) + outliers = X[outlier_mask] + logging.info('# outliers: {}'.format(len(outliers))) + + in_data_flat[outliers] = 0 + out_data = np.reshape(in_data_flat, init_shape) + + return out_data + + def smooth_to_fwhm(data, fwhm): """ Smooth a volume to given FWHM. diff --git a/scilpy/utils/metrics_tools.py b/scilpy/utils/metrics_tools.py index 66860549d..d4e5be865 100644 --- a/scilpy/utils/metrics_tools.py +++ b/scilpy/utils/metrics_tools.py @@ -377,30 +377,3 @@ def plot_metrics_stats(means, stds, title=None, xlabel=None, plt.close(fig) return fig - - -def get_roi_metrics_mean_std(density_map, metrics_files): - """ - Returns the mean and standard deviation of each metric, using the - provided density map. This can be a binary mask, - or contain weighted values between 0 and 1. - - Parameters - ------------ - density_map : ndarray - 3D numpy array containing a density map. - metrics_files : sequence - list of nibabel objects representing the metrics files. - - Returns - --------- - stats : list - list of tuples where the first element of the tuple is the mean - of a metric, and the second element is the standard deviation. - - """ - - return map(lambda metric_file: - weighted_mean_std(density_map, - metric_file.get_fdata(dtype=np.float64)), - metrics_files) diff --git a/scripts/scil_volume_apply_transform.py b/scripts/scil_volume_apply_transform.py index 6a6085dde..7ec3c058b 100755 --- a/scripts/scil_volume_apply_transform.py +++ b/scripts/scil_volume_apply_transform.py @@ -53,14 +53,11 @@ def main(): args = parser.parse_args() logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + # Verifications assert_inputs_exist(parser, [args.in_file, args.in_target_file, args.in_transfo]) assert_outputs_exist(parser, args, args.out_name) - transfo = load_matrix_in_any_format(args.in_transfo) - if args.inverse: - transfo = np.linalg.inv(transfo) - _, ref_extension = split_name_with_nii(args.in_target_file) _, in_extension = split_name_with_nii(args.in_file) if ref_extension not in ['.nii', '.nii.gz']: @@ -69,16 +66,14 @@ def main(): if in_extension not in ['.nii', '.nii.gz']: parser.error('{} is an unsupported format.'.format(args.in_file)) - # Load images and validate input type. + # Loading + transfo = load_matrix_in_any_format(args.in_transfo) + if args.inverse: + transfo = np.linalg.inv(transfo) moving = nib.load(args.in_file) - - if moving.get_fdata().ndim == 4: - logging.warning('You are applying a transform to a 4D dwi volume, ' - 'make sure to rotate your bvecs with ' - 'scil_gradients_apply_transform.py') - reference = nib.load(args.in_target_file) + # Processing, saving warped_img = apply_transform( transfo, reference, moving, keep_dtype=args.keep_dtype) diff --git a/scripts/scil_volume_b0_synthesis.py b/scripts/scil_volume_b0_synthesis.py index 79516b8a5..c12ba1a91 100644 --- a/scripts/scil_volume_b0_synthesis.py +++ b/scripts/scil_volume_b0_synthesis.py @@ -8,11 +8,11 @@ template, run SyNb0 and then transform the result back to the original space. SyNb0 is a deep learning model that predicts a synthetic a distortion-free -b0 image from a distorted b0 and T1w +b0 image from a distorted b0 and T1w. -This script must be used carefully, as it is not meant to be used in an -environment with the following dependencies already installed (not default -in Scilpy): +This script must be used carefully, as it is meant to be used in an +environment with the following dependencies already installed (not installed by +default in Scilpy): - tensorflow-addons - tensorrt - tensorflow @@ -20,33 +20,22 @@ import argparse import logging -import os -import warnings -# Disable tensorflow warnings -with warnings.catch_warnings(): - os.environ['TF_CPP_MIN_LOG_LEVEL'] = '3' - warnings.simplefilter("ignore") - from dipy.nn.synb0 import Synb0 - -from dipy.align.imaffine import AffineMap -from dipy.segment.tissue import TissueClassifierHMRF import nibabel as nib -import numpy as np -from scipy.ndimage import gaussian_filter +from scilpy.image.volume_b0_synthesis import compute_b0_synthesis +from scilpy.io.image import get_data_as_mask from scilpy.io.fetcher import get_synb0_template_path -from scilpy.io.utils import (add_overwrite_arg, - add_verbose_arg, - assert_inputs_exist, - assert_outputs_exist) -from scilpy.image.volume_operations import register_image +from scilpy.io.utils import (add_overwrite_arg, add_verbose_arg, + assert_inputs_exist, assert_outputs_exist, + assert_headers_compatible) EPILOG = """ [1] Schilling, Kurt G., et al. "Synthesized b0 for diffusion distortion correction (Synb0-DisCo)." Magnetic resonance imaging 64 (2019): 62-70. """ + def _build_arg_parser(): p = argparse.ArgumentParser( description=__doc__, @@ -72,73 +61,38 @@ def _build_arg_parser(): def main(): parser = _build_arg_parser() args = parser.parse_args() - assert_inputs_exist(parser, [args.in_b0, args.in_t1]) + assert_inputs_exist(parser, [args.in_b0, args.in_b0_mask, + args.in_t1, args.in_t1_mask]) assert_outputs_exist(parser, args, args.out_b0) + assert_headers_compatible(parser, [args.in_b0, args.in_b0_mask]) + assert_headers_compatible(parser, [args.in_t1, args.in_t1_mask]) logging.getLogger().setLevel(logging.getLevelName(args.verbose)) logging.info('The usage of synthetic b0 is not fully tested.' 'Be careful when using it.') + # Loading template_img = nib.load(get_synb0_template_path()) template_data = template_img.get_fdata() b0_img = nib.load(args.in_b0) - b0_skull_data = b0_img.get_fdata() - b0_mask_img = nib.load(args.in_b0_mask) - b0_mask_data = b0_mask_img.get_fdata() + b0_data = b0_img.get_fdata() + b0_mask_data = get_data_as_mask(nib.load(args.in_b0_mask)) t1_img = nib.load(args.in_t1) - t1_skull_data = t1_img.get_fdata() - t1_mask_img = nib.load(args.in_t1_mask) - t1_mask_data = t1_mask_img.get_fdata() - - b0_bet_data = np.zeros(b0_skull_data.shape) - b0_bet_data[b0_mask_data > 0] = b0_skull_data[b0_mask_data > 0] - t1_bet_data = np.zeros(t1_skull_data.shape) - t1_bet_data[t1_mask_data > 0] = t1_skull_data[t1_mask_data > 0] - - # Crude estimation of the WM mean intensity and normalization - logging.info('Estimating WM mean intensity') - hmrf = TissueClassifierHMRF() - t1_bet_data = gaussian_filter(t1_bet_data, 2) - _, final_segmentation, _ = hmrf.classify(t1_bet_data, 3, 0.25, - tolerance=1e-4, max_iter=5) - avg_wm = np.mean(t1_skull_data[final_segmentation == 3]) - t1_skull_data /= avg_wm - t1_skull_data *= 110 - - # SyNB0 works only in a standard space, so we need to register the images - logging.info('Registering images') - # Use the BET image for registration - t1_bet_to_b0, t1_bet_to_b0_transform = register_image(b0_bet_data, - b0_img.affine, - t1_bet_data, - t1_img.affine, - fine=True) - affine_map = AffineMap(t1_bet_to_b0_transform, - b0_skull_data.shape, b0_img.affine, - t1_skull_data.shape, t1_img.affine) - t1_skull_to_b0 = affine_map.transform(t1_skull_data.astype(np.float64)) - - # Then register to MNI (using the BET again) - _, t1_bet_to_b0_to_mni_transform = register_image(template_data, - template_img.affine, - t1_bet_to_b0, - b0_img.affine, - fine=True) - affine_map = AffineMap(t1_bet_to_b0_to_mni_transform, - template_data.shape, template_img.affine, - b0_skull_data.shape, b0_img.affine) - - # But for prediction, we want the skull - b0_skull_to_mni = affine_map.transform(b0_skull_data.astype(np.float64)) - t1_skull_to_mni = affine_map.transform(t1_skull_to_b0.astype(np.float64)) - - logging.info('Running SyN-B0') - SyNb0 = Synb0(args.verbose) - rev_b0 = SyNb0.predict(b0_skull_to_mni, t1_skull_to_mni) - rev_b0 = affine_map.transform_inverse(rev_b0.astype(np.float64)) + t1_data = t1_img.get_fdata() + t1_mask_data = get_data_as_mask(nib.load(args.in_t1_mask)) + + b0_bet_data = b0_data * b0_mask_data + t1_bet_data = t1_data * t1_mask_data + + # Processing + verbose = True if args.verbose in ['INFO', 'DEBUG'] else False + rev_b0 = compute_b0_synthesis(t1_data, t1_bet_data, b0_data, b0_bet_data, + template_data, t1_img.affine, b0_img.affine, + template_img.affine, verbose) + # Saving dtype = b0_img.get_data_dtype() nib.save(nib.Nifti1Image(rev_b0.astype(dtype), b0_img.affine), args.out_b0) diff --git a/scripts/scil_volume_count_non_zero_voxels.py b/scripts/scil_volume_count_non_zero_voxels.py index c8b984df9..e58cb7184 100755 --- a/scripts/scil_volume_count_non_zero_voxels.py +++ b/scripts/scil_volume_count_non_zero_voxels.py @@ -20,7 +20,8 @@ import numpy as np from scilpy.image.volume_operations import count_non_zero_voxels -from scilpy.io.utils import assert_inputs_exist, add_verbose_arg +from scilpy.io.utils import assert_inputs_exist, add_verbose_arg, \ + assert_outputs_exist, add_overwrite_arg def _build_arg_parser(): @@ -31,13 +32,14 @@ def _build_arg_parser(): p.add_argument( '--out', metavar='OUT_FILE', dest='out_filename', - help='name of the output file, which will be saved as a text file.') + help='Name of the output file, which will be saved as a text file.') p.add_argument( '--stats', action='store_true', dest='stats_format', - help='output the value using a stats format. Using this syntax will ' - 'append\na line to the output file, instead of creating a file ' - 'with only one line.\nThis is useful to create a file to be used ' - 'as the source of data for a graph.\nCan be combined with --id') + help='If set, output the value using a stats format. Using this ' + 'synthax will append\na line to the output file, instead of ' + 'creating a file with only one line.\nThis is useful to create ' + 'a file to be used as the source of data for a graph.\nCan be ' + 'combined with --id') p.add_argument( '--id', dest='value_id', help='Id of the current count. If used, the value of this argument ' @@ -45,6 +47,7 @@ def _build_arg_parser(): 'Mostly useful with --stats.') add_verbose_arg(p) + add_overwrite_arg(p) return p @@ -55,7 +58,12 @@ def main(): logging.getLogger().setLevel(logging.getLevelName(args.verbose)) assert_inputs_exist(parser, args.in_image) - # out_filename can exist or not + if not args.stats_format: + assert_outputs_exist(parser, args, [], args.out_filename) + elif args.overwrite: + parser.error("Using -f together with --stats is unclear. Please " + "either remove --stats to start anew, or remove -f to " + "append to existing file.") # Load image file im = nib.load(args.in_image).get_fdata(dtype=np.float32) diff --git a/scripts/scil_volume_crop.py b/scripts/scil_volume_crop.py index 71230fc08..b4982d813 100755 --- a/scripts/scil_volume_crop.py +++ b/scripts/scil_volume_crop.py @@ -50,7 +50,7 @@ def _build_arg_parser(): 'the bounding box to crop input file.') g1.add_argument('--output_bbox', help='Path of the pickle file where to write the ' - 'computed bounding box.') + 'computed bounding box. (.pickle extension)') return p diff --git a/scripts/scil_volume_remove_outliers_ransac.py b/scripts/scil_volume_remove_outliers_ransac.py index 7c3377b35..14da34c21 100755 --- a/scripts/scil_volume_remove_outliers_ransac.py +++ b/scripts/scil_volume_remove_outliers_ransac.py @@ -15,8 +15,8 @@ import nibabel as nib import numpy as np -from sklearn import linear_model +from scilpy.image.volume_operations import remove_outliers_ransac from scilpy.io.utils import (add_overwrite_arg, add_verbose_arg, assert_inputs_exist, @@ -33,14 +33,14 @@ def _build_arg_parser(): help='Corrected Nifti image.') p.add_argument('--min_fit', type=int, default=50, - help='The minimum number of data values required ' + - 'to fit the model. [%(default)s]') + help='The minimum number of data values required to fit ' + 'the model. [%(default)s]') p.add_argument('--max_iter', type=int, default=1000, - help='The maximum number of iterations allowed ' + - 'in the algorithm. [%(default)s]') + help='The maximum number of iterations allowed in the ' + 'algorithm. [%(default)s]') p.add_argument('--fit_thr', type=float, default=1e-2, - help='Threshold value for determining when a data ' + - 'point fits a model. [%(default)s]') + help='Threshold value for determining when a data point ' + 'fits a model. [%(default)s]') add_verbose_arg(p) add_overwrite_arg(p) @@ -53,6 +53,7 @@ def main(): args = parser.parse_args() logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + # Verifications assert_inputs_exist(parser, args.in_image) assert_outputs_exist(parser, args, args.out_image) @@ -66,6 +67,7 @@ def main(): parser.error('--fit_thr should be greater than 0. Current value: {}' .format(args.fit_thr)) + # Loading in_img = nib.load(args.in_image) in_data = in_img.get_fdata(dtype=np.float32) @@ -73,26 +75,11 @@ def main(): logging.warning('Be careful, your image doesn\'t seem to be an ad, ' 'md or rd.') - in_data_flat = in_data.flatten() - in_nzr_ind = np.nonzero(in_data_flat) - in_nzr_val = np.array(in_data_flat[in_nzr_ind]) + # Processing + out_data = remove_outliers_ransac(in_data, args.min_fit, args.fit_thr, + args.max_iter) - X = in_nzr_ind[0][:, np.newaxis] - model_ransac = linear_model.RANSACRegressor( - base_estimator=linear_model.LinearRegression(), - min_samples=args.min_fit, - residual_threshold=args.fit_thr, - max_trials=args.max_iter) - model_ransac.fit(X, in_nzr_val) - - outlier_mask = np.logical_not(model_ransac.inlier_mask_) - outliers = X[outlier_mask] - - logging.info('# outliers: {}'.format(len(outliers))) - - in_data_flat[outliers] = 0 - - out_data = np.reshape(in_data_flat, in_img.shape) + # Saving nib.save(nib.Nifti1Image(out_data, in_img.affine, in_img.header), args.out_image) diff --git a/scripts/scil_volume_resample.py b/scripts/scil_volume_resample.py index 42ed7d88a..e65e5b853 100755 --- a/scripts/scil_volume_resample.py +++ b/scripts/scil_volume_resample.py @@ -14,7 +14,8 @@ import nibabel as nib from scilpy.io.utils import (add_verbose_arg, add_overwrite_arg, - assert_inputs_exist, assert_outputs_exist) + assert_inputs_exist, assert_outputs_exist, + assert_headers_compatible) from scilpy.image.volume_operations import resample_volume @@ -66,6 +67,8 @@ def main(): # Checking args assert_inputs_exist(parser, args.in_image, args.ref) assert_outputs_exist(parser, args, args.out_image) + assert_headers_compatible(parser, args.in_image, args.ref) + if args.enforce_dimensions and not args.ref: parser.error("Cannot enforce dimensions without a reference image") diff --git a/scripts/scil_volume_stats_in_ROI.py b/scripts/scil_volume_stats_in_ROI.py index 9de905a4c..831373fdc 100755 --- a/scripts/scil_volume_stats_in_ROI.py +++ b/scripts/scil_volume_stats_in_ROI.py @@ -3,14 +3,12 @@ """ Compute the statistics (mean, std) of scalar maps, which can represent -diffusion metrics, in a ROI. +diffusion metrics, in a ROI. Prints the results. The mask can either be a binary mask, or a weighting mask. If the mask is a weighting mask it should either contain floats between 0 and 1 or should be -normalized with --normalize_weights. - -IMPORTANT: if the mask contains weights (and not 0 and 1 exclusively), the -standard deviation will also be weighted. +normalized with --normalize_weights. IMPORTANT: if the mask contains weights +(and not 0 and 1 exclusively), the standard deviation will also be weighted. """ import argparse @@ -28,7 +26,7 @@ assert_inputs_exist, assert_headers_compatible) from scilpy.utils.filenames import split_name_with_nii -from scilpy.utils.metrics_tools import get_roi_metrics_mean_std +from scilpy.utils.metrics_tools import weighted_mean_std def _build_arg_parser(): @@ -39,20 +37,20 @@ def _build_arg_parser(): help='Mask volume filename.\nCan be a binary mask or a ' 'weighted mask.') - p_metric = p.add_argument_group('Metrics input options') - g_metric = p_metric.add_mutually_exclusive_group(required=True) - g_metric.add_argument('--metrics_dir', - help='Metrics files directory. Name of the ' - 'directory containing the metrics files.') - g_metric.add_argument('--metrics', dest='metrics_file_list', nargs='+', - help='Metrics nifti filename. List of the names of ' - 'the metrics file, in nifti format.') + g = p.add_argument_group('Metrics input options') + gg = g.add_mutually_exclusive_group(required=True) + gg.add_argument('--metrics_dir', metavar='dir', + help='Name of the directory containing metrics files: ' + 'we will \nload all nifti files.') + gg.add_argument('--metrics', dest='metrics_file_list', nargs='+', + metavar='file', + help='Metrics nifti filename. List of the names of the ' + 'metrics file, \nin nifti format.') p.add_argument('--bin', action='store_true', - help='If set, will consider every value of the mask ' - 'higher than 0 to be part of the mask, and set to 1 ' - '(equivalent weighting for every voxel).') - + help='If set, will consider every value of the mask higher' + 'than 0 to be \npart of the mask (equivalent ' + 'weighting for every voxel).') p.add_argument('--normalize_weights', action='store_true', help='If set, the weights will be normalized to the [0,1] ' 'range.') @@ -69,81 +67,68 @@ def main(): args = parser.parse_args() logging.getLogger().setLevel(logging.getLevelName(args.verbose)) - if args.metrics_dir and os.path.exists(args.metrics_dir): - list_metrics_files = glob.glob(os.path.join(args.metrics_dir, - '*nii.gz')) - assert_inputs_exist(parser, [args.in_mask] + list_metrics_files) - assert_headers_compatible(parser, [args.in_mask] + list_metrics_files) - elif args.metrics_file_list: - assert_inputs_exist(parser, [args.in_mask] + args.metrics_file_list) - assert_headers_compatible(parser, - [args.in_mask] + args.metrics_file_list) - - # Load mask and validate content depending on flags - mask_img = nib.load(args.in_mask) + # Verifications - if len(mask_img.shape) > 3: - logging.error('Mask should be a 3D image.') - - # Can be a weighted image - mask_data = mask_img.get_fdata(dtype=np.float32) + # Get input list. Either all files in dir or given list. + if args.metrics_dir: + if not os.path.exists(args.metrics_dir): + parser.error("Metrics directory does not exist: {}" + .format(args.metrics_dir)) + assert_inputs_exist(parser, args.in_mask) + + tmp_file_list = glob.glob(os.path.join(args.metrics_dir, '*nii.gz')) + args.metrics_file_list = [os.path.join(args.metrics_dir, f) + for f in tmp_file_list] + else: + assert_inputs_exist(parser, [args.in_mask] + args.metrics_file_list) + assert_headers_compatible(parser, + [args.in_mask] + args.metrics_file_list) + # Loading + mask_data = nib.load(args.in_mask).get_fdata(dtype=np.float32) + if len(mask_data.shape) > 3: + parser.error('Mask should be a 3D image.') if np.min(mask_data) < 0: - logging.error('Mask should not contain negative values.') + parser.error('Mask should not contain negative values.') + roi_name = split_name_with_nii(os.path.basename(args.in_mask))[0] # Discussion about the way the normalization is done. # https://github.com/scilus/scilpy/pull/202#discussion_r411355609 + # Summary: + # 1) We don't want to normalize with data = (data-min) / (max-min) because + # it zeroes out the minimal values of the array. This is not a large error + # source, but not preferable. + # 2) data = data / max(data) or data = data / sum(data): in practice, when + # we use them in numpy using their weights argument, leads to the same + # result. if args.normalize_weights: mask_data /= np.max(mask_data) - - if np.min(mask_data) < 0.0 or np.max(mask_data) > 1.0: + elif args.bin: + mask_data[np.where(mask_data > 0.0)] = 1.0 + elif np.min(mask_data) < 0.0 or np.max(mask_data) > 1.0: parser.error('Mask data should only contain values between 0 and 1. ' 'Try --normalize_weights.') - if args.bin: - mask_data[np.where(mask_data > 0.0)] = 1.0 - - # Load all metrics files. - metrics_files = [] - if args.metrics_dir: - for f in sorted(os.listdir(args.metrics_dir)): - metric_img = nib.load(os.path.join(args.metrics_dir, f)) - if len(metric_img.shape)==3: - # Check if NaNs in metrics - if np.any(np.isnan(metric_img.get_fdata(dtype=np.float64))): - logging.warning('Metric \"{}\" contains some NaN.'.format(metric_img.get_filename()) + - ' Ignoring voxels with NaN.') - metrics_files.append(metric_img) - else: - parser.error('Metric {} is not compatible ({}D image).'.format(os.path.join(args.metrics_dir, f), - len(metric_img.shape))) - elif args.metrics_file_list: - assert_headers_compatible(parser, [args.in_mask] + - args.metrics_file_list) - for f in args.metrics_file_list: - metric_img = nib.load(f) - if len(metric_img.shape)==3: - # Check if NaNs in metrics - if np.any(np.isnan(metric_img.get_fdata(dtype=np.float64))): - logging.warning('Metric \"{}\" contains some NaN.'.format(metric_img.get_filename()) + - ' Ignoring voxels with NaN.') - metrics_files.append(metric_img) - else: - parser.error('Metric {} is not compatible ({}D image).'.format(f, - len(metric_img.shape))) - # Compute the mean values and standard deviations - stats = get_roi_metrics_mean_std(mask_data, metrics_files) - - roi_name = split_name_with_nii(os.path.basename(args.in_mask))[0] + # Load and process all metrics files. json_stats = {roi_name: {}} - for metric_file, (mean, std) in zip(metrics_files, stats): - metric_name = split_name_with_nii( - os.path.basename(metric_file.get_filename()))[0] - json_stats[roi_name][metric_name] = { - 'mean': mean.item(), - 'std': std.item() - } - + for f in args.metrics_file_list: + metric_img = nib.load(f) + metric_name = split_name_with_nii(os.path.basename(f))[0] + if len(metric_img.shape) == 3: + data = metric_img.get_fdata(dtype=np.float64) + if np.any(np.isnan(data)): + logging.warning("Metric '{}' contains some NaN. Ignoring " + "voxels with NaN." + .format(os.path.basename(f))) + mean, std = weighted_mean_std(mask_data, data) + json_stats[roi_name][metric_name] = {'mean': mean, + 'std': std} + else: + parser.error( + 'Metric {} is not a 3D image ({}D shape).' + .format(f, len(metric_img.shape))) + + # Print results print(json.dumps(json_stats, indent=args.indent, sort_keys=args.sort_keys)) diff --git a/scripts/scil_volume_stats_in_labels.py b/scripts/scil_volume_stats_in_labels.py index c4b926afb..07fc290b7 100755 --- a/scripts/scil_volume_stats_in_labels.py +++ b/scripts/scil_volume_stats_in_labels.py @@ -1,10 +1,11 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ -Computes the information from the seeding map for each cortical region -(corresponding to an atlas) associated with a specific bundle. -Here we want to estimate the seeding attribution to cortical area -affected by the bundle +Computes the information from the input map for each cortical region +(corresponding to an atlas). + +Hint: For instance, this script could be useful if you have a seed map from a +specific bundle, to know from which regions it originated. Formerly: scil_compute_seed_by_labels.py """ @@ -16,23 +17,21 @@ import nibabel as nib import numpy as np -from scilpy.image.labels import get_data_as_labels -from scilpy.io.utils import (add_overwrite_arg, - add_verbose_arg, - assert_inputs_exist) +from scilpy.image.labels import get_data_as_labels, get_stats_in_label +from scilpy.io.utils import (add_overwrite_arg, add_verbose_arg, + assert_inputs_exist, assert_headers_compatible) def _build_arg_parser(): - p = argparse.ArgumentParser( - description=__doc__, - formatter_class=argparse.RawTextHelpFormatter) + p = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawTextHelpFormatter) p.add_argument('in_labels', help='Path of the input label file.') p.add_argument('in_labels_lut', help='Path of the LUT file corresponding to labels,' 'used to name the regions of interest.') - p.add_argument('in_seed_maps', - help='Path of the input seed map file.') + p.add_argument('in_map', + help='Path of the input map file. Expecting a 3D file.') add_verbose_arg(p) add_overwrite_arg(p) @@ -45,42 +44,22 @@ def main(): args = parser.parse_args() logging.getLogger().setLevel(logging.getLevelName(args.verbose)) - required = args.in_labels, args.in_seed_maps, args.in_labels_lut - assert_inputs_exist(parser, required) - - # Load atlas image - label_img = nib.load(args.in_labels) - label_img_data = get_data_as_labels(label_img) + # Verifications + assert_inputs_exist(parser, + [args.in_labels, args.in_map, args.in_labels_lut]) + assert_headers_compatible(parser, [args.in_labels, args.in_map]) - # Load atlas lut + # Loading + label_data = get_data_as_labels(nib.load(args.in_labels)) with open(args.in_labels_lut) as f: label_dict = json.load(f) - (label_indices, label_names) = zip(*label_dict.items()) - - # Load seed image - seed_img = nib.load(args.in_seed_maps) - seed_img_data = seed_img.get_fdata(dtype=np.float32) - - for label, name in zip(label_indices, label_names): - label = int(label) - if label != 0: - curr_data = (seed_img_data[np.where(label_img_data == label)]) - nb_vx_roi = np.count_nonzero(label_img_data == label) - nb_seed_vx = np.count_nonzero(curr_data) - - if nb_seed_vx != 0: - mean_seed = np.sum(curr_data)/nb_seed_vx - max_seed = np.max(curr_data) - std_seed = np.sqrt(np.mean(abs(curr_data[curr_data != 0] - - mean_seed)**2)) - - print(json.dumps({'ROI-idx': label, - 'ROI-name': str(name), - 'nb-vx-roi': int(nb_vx_roi), - 'nb-vx-seed': int(nb_seed_vx), - 'max': int(max_seed), - 'mean': float(mean_seed), - 'std': float(std_seed)})) + map_data = nib.load(args.in_map).get_fdata(dtype=np.float32) + if len(map_data.shape) > 3: + parser.error('Mask should be a 3D image.') + + # Process + out_dict = get_stats_in_label(map_data, label_data, label_dict) + print(json.dumps(out_dict)) if __name__ == "__main__": diff --git a/scripts/tests/test_volume_b0_synthesis.py b/scripts/tests/test_volume_b0_synthesis.py index 003083740..a9c811dc7 100644 --- a/scripts/tests/test_volume_b0_synthesis.py +++ b/scripts/tests/test_volume_b0_synthesis.py @@ -7,12 +7,16 @@ import nibabel as nib import numpy as np import pytest +from _pytest.outcomes import Skipped from scilpy import SCILPY_HOME from scilpy.io.fetcher import fetch_data, get_testing_files_dict -tensorflow = pytest.importorskip("tensorflow") - +try: + tensorflow = pytest.importorskip("tensorflow") + have_tensorflow = True +except Skipped: + have_tensorflow = False # If they already exist, this only takes 5 seconds (check md5sum) fetch_data(get_testing_files_dict(), keys=['others.zip', 'processing.zip']) @@ -24,27 +28,25 @@ def test_help_option(script_runner): assert ret.success -@pytest.mark.skipif(tensorflow is None, reason="Tensorflow not installed") def test_synthesis(script_runner, monkeypatch): - monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) - in_t1 = os.path.join(SCILPY_HOME, 'others', - 't1.nii.gz') - in_b0 = os.path.join(SCILPY_HOME, 'processing', - 'b0_mean.nii.gz') - - t1_img = nib.load(in_t1) - b0_img = nib.load(in_b0) - t1_data = t1_img.get_fdata() - b0_data = b0_img.get_fdata() - t1_data[t1_data > 0] = 1 - b0_data[b0_data > 0] = 1 - nib.save(nib.Nifti1Image(t1_data.astype(np.uint8), t1_img.affine), - 't1_mask.nii.gz') - nib.save(nib.Nifti1Image(b0_data.astype(np.uint8), b0_img.affine), - 'b0_mask.nii.gz') - - ret = script_runner.run('scil_volume_b0_synthesis.py', - in_t1, 't1_mask.nii.gz', - in_b0, 'b0_mask.nii.gz', - 'b0_synthesized.nii.gz', '-v') - assert ret.success + if have_tensorflow: + monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) + in_t1 = os.path.join(SCILPY_HOME, 'others', 't1.nii.gz') + in_b0 = os.path.join(SCILPY_HOME, 'processing', 'b0_mean.nii.gz') + + t1_img = nib.load(in_t1) + b0_img = nib.load(in_b0) + t1_data = t1_img.get_fdata() + b0_data = b0_img.get_fdata() + t1_data[t1_data != 0] = 1 + b0_data[b0_data > 0] = 1 + nib.save(nib.Nifti1Image(t1_data.astype(np.uint8), t1_img.affine), + 't1_mask.nii.gz') + nib.save(nib.Nifti1Image(b0_data.astype(np.uint8), b0_img.affine), + 'b0_mask.nii.gz') + + ret = script_runner.run('scil_volume_b0_synthesis.py', + in_t1, 't1_mask.nii.gz', + in_b0, 'b0_mask.nii.gz', + 'b0_synthesized.nii.gz', '-v') + assert ret.success diff --git a/scripts/tests/test_volume_count_non_zero_voxels.py b/scripts/tests/test_volume_count_non_zero_voxels.py index 6b000e08d..cec9c23bc 100644 --- a/scripts/tests/test_volume_count_non_zero_voxels.py +++ b/scripts/tests/test_volume_count_non_zero_voxels.py @@ -17,9 +17,21 @@ def test_help_option(script_runner): assert ret.success -def test_execution_others(script_runner, monkeypatch): +def test_execution_simple_print(script_runner, monkeypatch): monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) - in_img = os.path.join(SCILPY_HOME, 'others', - 'rgb.nii.gz') + in_img = os.path.join(SCILPY_HOME, 'others', 'rgb.nii.gz') ret = script_runner.run('scil_volume_count_non_zero_voxels.py', in_img) assert ret.success + + +def test_execution_save_in_file(script_runner): + os.chdir(os.path.expanduser(tmp_dir.name)) + in_img = os.path.join(SCILPY_HOME, 'others', 'rgb.nii.gz') + ret = script_runner.run('scil_volume_count_non_zero_voxels.py', in_img, + '--out', 'printed.txt') + assert ret.success + + # Then re-use the same out file with --stats + ret = script_runner.run('scil_volume_count_non_zero_voxels.py', in_img, + '--out', 'printed.txt', '--stats') + assert ret.success diff --git a/scripts/tests/test_volume_crop.py b/scripts/tests/test_volume_crop.py index f4ad89a8d..145cd528b 100644 --- a/scripts/tests/test_volume_crop.py +++ b/scripts/tests/test_volume_crop.py @@ -20,5 +20,11 @@ def test_help_option(script_runner): def test_execution_processing(script_runner, monkeypatch): monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) in_dwi = os.path.join(SCILPY_HOME, 'processing', 'dwi.nii.gz') - ret = script_runner.run('scil_volume_crop.py', in_dwi, 'dwi_crop.nii.gz') + ret = script_runner.run('scil_volume_crop.py', in_dwi, 'dwi_crop.nii.gz', + '--output_bbox', 'bbox.pickle') + assert ret.success + + # Then try to load back the same box + ret = script_runner.run('scil_crop_volume.py', in_dwi, 'dwi_crop2.nii.gz', + '--input_bbox', 'bbox.pickle') assert ret.success diff --git a/scripts/tests/test_volume_stats_in_labels.py b/scripts/tests/test_volume_stats_in_labels.py index e0e1c4bd8..6190f3d46 100644 --- a/scripts/tests/test_volume_stats_in_labels.py +++ b/scripts/tests/test_volume_stats_in_labels.py @@ -1,7 +1,19 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- +import os + +from scilpy import SCILPY_HOME def test_help_option(script_runner): ret = script_runner.run('scil_volume_stats_in_labels.py', '--help') assert ret.success + + +def test_execution(script_runner): + in_map = os.path.join(SCILPY_HOME, 'plot', 'fa.nii.gz') + in_atlas = os.path.join(SCILPY_HOME, 'plot', 'atlas_brainnetome.nii.gz') + atlas_lut = os.path.join(SCILPY_HOME, 'plot', 'atlas_brainnetome.json') + ret = script_runner.run('scil_volume_stats_in_labels.py', + in_atlas, atlas_lut, in_map) + assert ret.success