diff --git a/scilpy/stats/matrix_stats.py b/scilpy/stats/matrix_stats.py index d77b5e1fc..fe15d5dbd 100644 --- a/scilpy/stats/matrix_stats.py +++ b/scilpy/stats/matrix_stats.py @@ -1,4 +1,5 @@ # -*- coding: utf-8 -*- +import itertools import logging import bct @@ -7,6 +8,8 @@ from scipy.stats import t as stats_t from statsmodels.stats.multitest import multipletests +from scilpy.tractanalysis.reproducibility_measures import compute_dice_voxel + def _ttest_stat_only(x, y, tail): t = np.mean(x) - np.mean(y) @@ -180,3 +183,59 @@ def omega_sigma(matrix): (path_length / path_length_rand) return float(omega), float(sigma) + + +def pairwise_agreement(matrices, ref_matrix=None, normalize=False): + """ + The similarity measures will be computed for each pair. Alternatively, you + can compare all matrices to a single reference, ref_matrix. + + Parameters + ---------- + matrices: list[np.ndarray] + Input matrices + ref_matrix: Optional[np.ndarray] + Optional reference matrix. + normalize: bool + If true, will normalize all matrices from zero to one. + + Returns + ------- + output_measures_dict: dict + A dict with list of values for each pair of matrices: + { + 'RMSE': root-mean-square error + 'correlation': correlation + 'w_dice_voxels': weighted dice, agreement of the values. + 'dice_voxels': agreement of the binarized matrices + } + """ + def _prepare_matrix(tmp_mat): + # Removing the min now simplifies computations + tmp_mat -= np.min(tmp_mat) + if normalize: + return tmp_mat / np.max(tmp_mat) + return tmp_mat + + matrices = [_prepare_matrix(m) for m in matrices] + if ref_matrix is not None: + ref_matrix = _prepare_matrix(ref_matrix) + + output_measures_dict = {'RMSE': [], 'correlation': [], + 'w_dice_voxels': [], 'dice_voxels': []} + + if ref_matrix is not None: + pairs = list(itertools.product(matrices, [ref_matrix])) + else: + pairs = list(itertools.combinations(matrices, r=2)) + + for i in pairs: + rmse = np.sqrt(np.mean((i[0] - i[1]) ** 2)) + output_measures_dict['RMSE'].append(rmse) + corrcoef = np.corrcoef(i[0].ravel(), i[1].ravel()) + output_measures_dict['correlation'].append(corrcoef[0][1]) + dice, w_dice = compute_dice_voxel(i[0], i[1]) + output_measures_dict['dice_voxels'].append(dice) + output_measures_dict['w_dice_voxels'].append(w_dice) + + return output_measures_dict diff --git a/scripts/scil_connectivity_pairwise_agreement.py b/scripts/scil_connectivity_pairwise_agreement.py index 41e62687d..79eeb2023 100755 --- a/scripts/scil_connectivity_pairwise_agreement.py +++ b/scripts/scil_connectivity_pairwise_agreement.py @@ -2,10 +2,17 @@ # -*- coding: utf-8 -*- """ -Evaluate pair-wise similarity measures of connectivity matrix. +Evaluate pair-wise similarity measures of connectivity matrices. The computed similarity measures are: -sum of square difference and pearson correlation coefficent +- RMSE: root-mean-square difference +- Pearson correlation coefficent +- w-dice: weighted dice, agreement of the values. +- dice: agreement of the binarized matrices + +If more than two matrices are given in input, the similarity measures will be +computed for each pair. Alternatively, you can compare all matrices to a +single reference, using --single_compare. Formerly: scil_evaluate_connectivity_pairwaise_agreement_measures.py """ @@ -23,6 +30,7 @@ assert_inputs_exist, assert_outputs_exist, load_matrix_in_any_format) +from scilpy.stats.matrix_stats import pairwise_agreement from scilpy.tractanalysis.reproducibility_measures import compute_dice_voxel @@ -53,48 +61,27 @@ def main(): args = parser.parse_args() logging.getLogger().setLevel(logging.getLevelName(args.verbose)) - assert_inputs_exist(parser, args.in_matrices) + assert_inputs_exist(parser, args.in_matrices, args.single_compare) assert_outputs_exist(parser, args, args.out_json) + # Check that --single_compare is not loaded twice + if args.single_compare and args.single_compare in args.in_matrices: + id = args.in_matrices.index(args.single_compare) + args.in_matrices.pop(id) + + # Load matrices all_matrices = [] for filename in args.in_matrices: - tmp_mat = load_matrix_in_any_format(filename) - tmp_mat = tmp_mat.astype(float) - tmp_mat -= np.min(tmp_mat) - if args.normalize: - all_matrices.append(tmp_mat / np.max(tmp_mat)) - else: - all_matrices.append(tmp_mat) - - if args.single_compare: - tmp_mat = load_matrix_in_any_format(args.single_compare) - tmp_mat = tmp_mat.astype(float) - tmp_mat -= np.min(tmp_mat) - if args.normalize: - all_matrices.append(tmp_mat / np.max(tmp_mat)) - else: - all_matrices.append(tmp_mat) - - output_measures_dict = {'RMSE': [], 'correlation': [], - 'w_dice_voxels': [], 'dice_voxels': []} + all_matrices.append(load_matrix_in_any_format(filename).astype(float)) + ref_matrix = None if args.single_compare: - if args.single_compare in args.in_matrices: - id = args.in_matrices.index(args.single_compare) - all_matrices.pop(id) - pairs = list(itertools.product(all_matrices[:-1], [all_matrices[-1]])) - else: - pairs = list(itertools.combinations(all_matrices, r=2)) - - for i in pairs: - rmse = np.sqrt(np.mean((i[0]-i[1])**2)) - output_measures_dict['RMSE'].append(rmse) - corrcoef = np.corrcoef(i[0].ravel(), i[1].ravel()) - output_measures_dict['correlation'].append(corrcoef[0][1]) - dice, w_dice = compute_dice_voxel(i[0], i[1]) - output_measures_dict['dice_voxels'].append(dice) - output_measures_dict['w_dice_voxels'].append(w_dice) + ref_matrix = load_matrix_in_any_format( + args.single_compare).astype(float) + # Compute and save + output_measures_dict = pairwise_agreement(all_matrices, ref_matrix, + normalize=args.normalize) with open(args.out_json, 'w') as outfile: json.dump(output_measures_dict, outfile, indent=args.indent, sort_keys=args.sort_keys)