Merge pull request #1087 from EmmaRenauld/next_connectivity2

Cleaning the connectivity pairwise agreement script

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

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)