Merge pull request #4 from InstitutoDOr/stelzer_stats

speed up calculation

Author: tbortolini

Code/permutation_stats.py

@@ -35,37 +35,33 @@ def calculate_threshold_statistics_map(
     # Initialize an empty map to store threshold statistics
     threshold_statistics_map = np.zeros_like(parcellated_masks[0], dtype=float)
 
-    # Iterate over each parcellated mask
     for mask_idx, mask in enumerate(parcellated_masks):
+        import time
+        t  = time.time()
+        
         # Get voxel coordinates where mask is True using np.nonzero
         voxel_indices = np.nonzero(mask)
         print(f"Mask {mask_idx} with {len(voxel_indices[0])} voxels.")
 
-        # Get memory usage in bytes
-        import psutil
-        memory_usage_bytes = psutil.Process().memory_info().rss
-        print(f"Memory usage: {memory_usage_bytes/1024**2} MB")
-
-        # Iterate over voxel coordinates
-        for voxel_idx in zip(*voxel_indices):
-            voxel_values = []
+        voxel_values = np.zeros((len(accuracy_maps), len(voxel_indices[0])))
 
-            # Load accuracy maps from image directory
-            for accuracy_map_path in accuracy_maps:
-                accuracy_map_img = nib.load(accuracy_map_path)
-                accuracy_map_data = accuracy_map_img.get_fdata()
+        # Load accuracy maps and extract voxel values
+        for map_idx, accuracy_map in enumerate(accuracy_maps):
+            accuracy_map_data = nib.load(accuracy_map).get_fdata()
+            voxel_values[map_idx, :] = accuracy_map_data[voxel_indices]
 
-                # Get voxel value from accuracy map at current voxel index
-                voxel_value = accuracy_map_data[voxel_idx]
-                voxel_values.append(voxel_value)
+        # Calculate the threshold statistic for each voxel using np.percentile
+        percentile_value = 100 * (1 - p_value_threshold)
+        threshold_statistic = np.percentile(voxel_values, percentile_value, axis=0)
 
-            # Calculate the threshold statistic using np.percentile
-            threshold_statistic = np.percentile(
-                voxel_values, 100 * (1 - p_value_threshold)
-            )
+        # Assign the computed threshold statistic to the map at the voxel indices
+        threshold_statistics_map[voxel_indices] = threshold_statistic
 
-            # Store the threshold statistic in the map at the current voxel index
-            threshold_statistics_map[voxel_idx] = threshold_statistic
+        import psutil
+        # Get memory usage in bytes after processing each mask
+        memory_usage_bytes = psutil.Process().memory_info().rss
+        print(f"Memory usage: {memory_usage_bytes/1024**2} MB")
+        print(f"Time for parcellation: {time.time()-t} s")
 
     return threshold_statistics_map
 
@@ -361,6 +357,7 @@ def keep_significant_clusters(accuracy_map, accuracy_cluster_map, cluster_size):
     permutation_dir = Path("/Users/sebastian.hoefle/projects/idor/brain-analysis/fixtures/permutations/random_accuracy_maps")
     brain_mask = Path("/Users/sebastian.hoefle/projects/idor/brain-analysis/fixtures/permutations/brain_mask.nii.gz")
     # This represents the true original accuracy map. Here we use a map generated by generate_fake_clusters
+    # This map represents the group accuracy map, that is the mean of the single-subject accuracy maps
     accuracy_map = Path("/Users/sebastian.hoefle/projects/idor/brain-analysis/fixtures/permutations/clustered_accuracy_map.nii.gz")
     base_output_path = Path("/Users/sebastian.hoefle/projects/idor/brain-analysis/fixtures/stat_results")
     perm_output_path = Path(base_output_path, "permuted")