We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
AxisError when computing `reporting.get_cluster_table` 
This issue concerns what seems to me like a bug in the `reporting.get_cluster_table` method. I have been encoutering a `AxisError` raised by numpy when trying to compute the table with certain `stat_threshold` values.

With a high threshold value (e.g `3.`) the computation works fine, but when using smaller threshold values (1e-2), the error is thrown.

The error happens because the `ijk` value passed to a numpy `apply_along_axis` [here](https://github.com/nilearn/nilearn/blob/ac1a934e3b2b4061f894b518c960412df9ea4f11/nilearn/reporting/_get_clusters_table.py#L83) is an empty list. I don't know why yet. You can reproduce the error with the snippet below.

I guess it makes less sense to extract peaks with a lower threshold, but it still worth investigating why the error occurs.

### Expected behavior

`get_clusters_table` return the clusters table whatever the `stat_threshold` value

### Actual behavior

`get_clusters_table` raises an error when `stat_threshold` values are small.

### Steps and code to reproduce bug

```python
import tempfile

from nilearn import reporting
from nilearn.datasets import fetch_neurovault_ids

DATA_PATH = tempfile.TemporaryDirectory()

collection_id = 307
neurovault_element = fetch_neurovault_ids(
    collection_ids=[collection_id],
    mode="download_new",
    data_dir=str(DATA_PATH),
)
sample_brain_map = neurovault_element.images[0]

# Working fine with a stat threshold of 3.
neurovault_peaks = reporting.get_clusters_table(
    sample_brain_map,
    stat_threshold=3.,
    min_distance=12.,
)

# Raising error with a stat threshold of 1e-2
neurovault_peaks = reporting.get_clusters_table(
    sample_brain_map,
    stat_threshold=1e-2,
    min_distance=12.,
)
```


--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `nilearn/reporting/_get_clusters_table.py`
Content:
```
1 """
2 This module implements plotting functions useful to report analysis results.
3 
4 Author: Martin Perez-Guevara, Elvis Dohmatob, 2017
5 """
6 
7 import warnings
8 from string import ascii_lowercase
9 
10 import numpy as np
11 import pandas as pd
12 import nibabel as nib
13 from scipy import ndimage
14 
15 from nilearn.image import get_data
16 from nilearn.image.resampling import coord_transform
17 from nilearn._utils import check_niimg_3d
18 
19 
20 def _local_max(data, affine, min_distance):
21     """Find all local maxima of the array, separated by at least min_distance.
22     Adapted from https://stackoverflow.com/a/22631583/2589328
23 
24     Parameters
25     ----------
26     data : array_like
27         3D array of with masked values for cluster.
28 
29     affine : np.ndarray
30         Square matrix specifying the position of the image array data
31         in a reference space.
32 
33     min_distance : int
34         Minimum distance between local maxima in ``data``, in terms of mm.
35 
36     Returns
37     -------
38     ijk : `numpy.ndarray`
39         (n_foci, 3) array of local maxima indices for cluster.
40 
41     vals : `numpy.ndarray`
42         (n_foci,) array of values from data at ijk.
43 
44     """
45     ijk, vals = _identify_subpeaks(data)
46     xyz, ijk, vals = _sort_subpeaks(ijk, vals, affine)
47     ijk, vals = _pare_subpeaks(xyz, ijk, vals, min_distance)
48     return ijk, vals
49 
50 
51 def _identify_subpeaks(data):
52     """Identify cluster peak and subpeaks based on minimum distance.
53 
54     Parameters
55     ----------
56     data : `numpy.ndarray`
57         3D array of with masked values for cluster.
58 
59     Returns
60     -------
61     ijk : `numpy.ndarray`
62         (n_foci, 3) array of local maxima indices for cluster.
63     vals : `numpy.ndarray`
64         (n_foci,) array of values from data at ijk.
65 
66     Notes
67     -----
68     When a cluster's local maxima correspond to contiguous voxels with the
69     same values (as in a binary cluster), this function determines the center
70     of mass for those voxels.
71     """
72     # Initial identification of subpeaks with minimal minimum distance
73     data_max = ndimage.filters.maximum_filter(data, 3)
74     maxima = data == data_max
75     data_min = ndimage.filters.minimum_filter(data, 3)
76     diff = (data_max - data_min) > 0
77     maxima[diff == 0] = 0
78 
79     labeled, n_subpeaks = ndimage.label(maxima)
80     labels_index = range(1, n_subpeaks + 1)
81     ijk = np.array(ndimage.center_of_mass(data, labeled, labels_index))
82     ijk = np.round(ijk).astype(int)
83     vals = np.apply_along_axis(
84         arr=ijk, axis=1, func1d=_get_val, input_arr=data
85     )
86     # Determine if all subpeaks are within the cluster
87     # They may not be if the cluster is binary and has a shape where the COM is
88     # outside the cluster, like a donut.
89     cluster_idx = np.vstack(np.where(labeled)).T.tolist()
90     subpeaks_outside_cluster = [
91         i
92         for i, peak_idx in enumerate(ijk.tolist())
93         if peak_idx not in cluster_idx
94     ]
95     vals[subpeaks_outside_cluster] = np.nan
96     if subpeaks_outside_cluster:
97         warnings.warn(
98             "Attention: At least one of the (sub)peaks falls outside of the "
99             "cluster body."
100         )
101     return ijk, vals
102 
103 
104 def _sort_subpeaks(ijk, vals, affine):
105     # Sort subpeaks in cluster in descending order of stat value
106     order = (-vals).argsort()
107     vals = vals[order]
108     ijk = ijk[order, :]
109     xyz = nib.affines.apply_affine(affine, ijk)  # Convert to xyz in mm
110     return xyz, ijk, vals
111 
112 
113 def _pare_subpeaks(xyz, ijk, vals, min_distance):
114     # Reduce list of subpeaks based on distance
115     keep_idx = np.ones(xyz.shape[0]).astype(bool)
116     for i in range(xyz.shape[0]):
117         for j in range(i + 1, xyz.shape[0]):
118             if keep_idx[i] == 1:
119                 dist = np.linalg.norm(xyz[i, :] - xyz[j, :])
120                 keep_idx[j] = dist > min_distance
121     ijk = ijk[keep_idx, :]
122     vals = vals[keep_idx]
123     return ijk, vals
124 
125 
126 def _get_val(row, input_arr):
127     """Small function for extracting values from array based on index.
128     """
129     i, j, k = row
130     return input_arr[i, j, k]
131 
132 
133 def get_clusters_table(stat_img, stat_threshold, cluster_threshold=None,
134                        two_sided=False, min_distance=8.):
135     """Creates pandas dataframe with img cluster statistics.
136 
137     Parameters
138     ----------
139     stat_img : Niimg-like object,
140        Statistical image (presumably in z- or p-scale).
141 
142     stat_threshold : `float`
143         Cluster forming threshold in same scale as `stat_img` (either a
144         p-value or z-scale value).
145 
146     cluster_threshold : `int` or `None`, optional
147         Cluster size threshold, in voxels.
148 
149     two_sided : `bool`, optional
150         Whether to employ two-sided thresholding or to evaluate positive values
151         only. Default=False.
152 
153     min_distance : `float`, optional
154         Minimum distance between subpeaks in mm. Default=8mm.
155 
156     Returns
157     -------
158     df : `pandas.DataFrame`
159         Table with peaks and subpeaks from thresholded `stat_img`. For binary
160         clusters (clusters with >1 voxel containing only one value), the table
161         reports the center of mass of the cluster,
162         rather than any peaks/subpeaks.
163 
164     """
165     cols = ['Cluster ID', 'X', 'Y', 'Z', 'Peak Stat', 'Cluster Size (mm3)']
166 
167     # check that stat_img is niimg-like object and 3D
168     stat_img = check_niimg_3d(stat_img)
169     # If cluster threshold is used, there is chance that stat_map will be
170     # modified, therefore copy is needed
171     if cluster_threshold is None:
172         stat_map = get_data(stat_img)
173     else:
174         stat_map = get_data(stat_img).copy()
175 
176     # Define array for 6-connectivity, aka NN1 or "faces"
177     conn_mat = np.zeros((3, 3, 3), int)
178     conn_mat[1, 1, :] = 1
179     conn_mat[1, :, 1] = 1
180     conn_mat[:, 1, 1] = 1
181     voxel_size = np.prod(stat_img.header.get_zooms())
182 
183     signs = [1, -1] if two_sided else [1]
184     no_clusters_found = True
185     rows = []
186     for sign in signs:
187         # Flip map if necessary
188         temp_stat_map = stat_map * sign
189 
190         # Binarize using CDT
191         binarized = temp_stat_map > stat_threshold
192         binarized = binarized.astype(int)
193 
194         # If the stat threshold is too high simply return an empty dataframe
195         if np.sum(binarized) == 0:
196             warnings.warn(
197                 'Attention: No clusters with stat {0} than {1}'.format(
198                     'higher' if sign == 1 else 'lower',
199                     stat_threshold * sign,
200                 )
201             )
202             continue
203 
204         # Extract connected components above cluster size threshold
205         label_map = ndimage.measurements.label(binarized, conn_mat)[0]
206         clust_ids = sorted(list(np.unique(label_map)[1:]))
207         for c_val in clust_ids:
208             if cluster_threshold is not None and np.sum(
209                     label_map == c_val) < cluster_threshold:
210                 temp_stat_map[label_map == c_val] = 0
211                 binarized[label_map == c_val] = 0
212 
213         # If the cluster threshold is too high simply return an empty dataframe
214         # this checks for stats higher than threshold after small clusters
215         # were removed from temp_stat_map
216         if np.sum(temp_stat_map > stat_threshold) == 0:
217             warnings.warn(
218                 'Attention: No clusters with more than {0} voxels'.format(
219                     cluster_threshold,
220                 )
221             )
222             continue
223 
224         # Now re-label and create table
225         label_map = ndimage.measurements.label(binarized, conn_mat)[0]
226         clust_ids = sorted(list(np.unique(label_map)[1:]))
227         peak_vals = np.array(
228             [np.max(temp_stat_map * (label_map == c)) for c in clust_ids])
229         # Sort by descending max value
230         clust_ids = [clust_ids[c] for c in (-peak_vals).argsort()]
231 
232         for c_id, c_val in enumerate(clust_ids):
233             cluster_mask = label_map == c_val
234             masked_data = temp_stat_map * cluster_mask
235 
236             cluster_size_mm = int(np.sum(cluster_mask) * voxel_size)
237 
238             # Get peaks, subpeaks and associated statistics
239             subpeak_ijk, subpeak_vals = _local_max(
240                 masked_data,
241                 stat_img.affine,
242                 min_distance=min_distance,
243             )
244             subpeak_vals *= sign  # flip signs if necessary
245             subpeak_xyz = np.asarray(
246                 coord_transform(
247                     subpeak_ijk[:, 0],
248                     subpeak_ijk[:, 1],
249                     subpeak_ijk[:, 2],
250                     stat_img.affine,
251                 )
252             ).tolist()
253             subpeak_xyz = np.array(subpeak_xyz).T
254 
255             # Only report peak and, at most, top 3 subpeaks.
256             n_subpeaks = np.min((len(subpeak_vals), 4))
257             for subpeak in range(n_subpeaks):
258                 if subpeak == 0:
259                     row = [
260                         c_id + 1,
261                         subpeak_xyz[subpeak, 0],
262                         subpeak_xyz[subpeak, 1],
263                         subpeak_xyz[subpeak, 2],
264                         subpeak_vals[subpeak],
265                         cluster_size_mm,
266                     ]
267                 else:
268                     # Subpeak naming convention is cluster num+letter:
269                     # 1a, 1b, etc
270                     sp_id = '{0}{1}'.format(
271                         c_id + 1,
272                         ascii_lowercase[subpeak - 1],
273                     )
274                     row = [
275                         sp_id,
276                         subpeak_xyz[subpeak, 0],
277                         subpeak_xyz[subpeak, 1],
278                         subpeak_xyz[subpeak, 2],
279                         subpeak_vals[subpeak],
280                         '',
281                     ]
282                 rows += [row]
283 
284         # If we reach this point, there are clusters in this sign
285         no_clusters_found = False
286 
287     if no_clusters_found:
288         df = pd.DataFrame(columns=cols)
289     else:
290         df = pd.DataFrame(columns=cols, data=rows)
291 
292     return df
293
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/nilearn/reporting/_get_clusters_table.py b/nilearn/reporting/_get_clusters_table.py
--- a/nilearn/reporting/_get_clusters_table.py
+++ b/nilearn/reporting/_get_clusters_table.py
@@ -12,7 +12,7 @@
 import nibabel as nib
 from scipy import ndimage
 
-from nilearn.image import get_data
+from nilearn._utils.niimg import _safe_get_data
 from nilearn.image.resampling import coord_transform
 from nilearn._utils import check_niimg_3d
 
@@ -168,10 +168,8 @@
     stat_img = check_niimg_3d(stat_img)
     # If cluster threshold is used, there is chance that stat_map will be
     # modified, therefore copy is needed
-    if cluster_threshold is None:
-        stat_map = get_data(stat_img)
-    else:
-        stat_map = get_data(stat_img).copy()
+    stat_map = _safe_get_data(stat_img, ensure_finite=True,
+                              copy_data=(cluster_threshold is not None))
 
     # Define array for 6-connectivity, aka NN1 or "faces"
     conn_mat = np.zeros((3, 3, 3), int)

{"golden_diff": "diff --git a/nilearn/reporting/_get_clusters_table.py b/nilearn/reporting/_get_clusters_table.py\n--- a/nilearn/reporting/_get_clusters_table.py\n+++ b/nilearn/reporting/_get_clusters_table.py\n@@ -12,7 +12,7 @@\n import nibabel as nib\n from scipy import ndimage\n \n-from nilearn.image import get_data\n+from nilearn._utils.niimg import _safe_get_data\n from nilearn.image.resampling import coord_transform\n from nilearn._utils import check_niimg_3d\n \n@@ -168,10 +168,8 @@\n     stat_img = check_niimg_3d(stat_img)\n     # If cluster threshold is used, there is chance that stat_map will be\n     # modified, therefore copy is needed\n-    if cluster_threshold is None:\n-        stat_map = get_data(stat_img)\n-    else:\n-        stat_map = get_data(stat_img).copy()\n+    stat_map = _safe_get_data(stat_img, ensure_finite=True,\n+                              copy_data=(cluster_threshold is not None))\n \n     # Define array for 6-connectivity, aka NN1 or \"faces\"\n     conn_mat = np.zeros((3, 3, 3), int)\n", "issue": "AxisError when computing `reporting.get_cluster_table` \nThis issue concerns what seems to me like a bug in the `reporting.get_cluster_table` method. I have been encoutering a `AxisError` raised by numpy when trying to compute the table with certain `stat_threshold` values.\r\n\r\nWith a high threshold value (e.g `3.`) the computation works fine, but when using smaller threshold values (1e-2), the error is thrown.\r\n\r\nThe error happens because the `ijk` value passed to a numpy `apply_along_axis` [here](https://github.com/nilearn/nilearn/blob/ac1a934e3b2b4061f894b518c960412df9ea4f11/nilearn/reporting/_get_clusters_table.py#L83) is an empty list. I don't know why yet. You can reproduce the error with the snippet below.\r\n\r\nI guess it makes less sense to extract peaks with a lower threshold, but it still worth investigating why the error occurs.\r\n\r\n### Expected behavior\r\n\r\n`get_clusters_table` return the clusters table whatever the `stat_threshold` value\r\n\r\n### Actual behavior\r\n\r\n`get_clusters_table` raises an error when `stat_threshold` values are small.\r\n\r\n### Steps and code to reproduce bug\r\n\r\n```python\r\nimport tempfile\r\n\r\nfrom nilearn import reporting\r\nfrom nilearn.datasets import fetch_neurovault_ids\r\n\r\nDATA_PATH = tempfile.TemporaryDirectory()\r\n\r\ncollection_id = 307\r\nneurovault_element = fetch_neurovault_ids(\r\n    collection_ids=[collection_id],\r\n    mode=\"download_new\",\r\n    data_dir=str(DATA_PATH),\r\n)\r\nsample_brain_map = neurovault_element.images[0]\r\n\r\n# Working fine with a stat threshold of 3.\r\nneurovault_peaks = reporting.get_clusters_table(\r\n    sample_brain_map,\r\n    stat_threshold=3.,\r\n    min_distance=12.,\r\n)\r\n\r\n# Raising error with a stat threshold of 1e-2\r\nneurovault_peaks = reporting.get_clusters_table(\r\n    sample_brain_map,\r\n    stat_threshold=1e-2,\r\n    min_distance=12.,\r\n)\r\n```\r\n\n", "before_files": [{"content": "\"\"\"\nThis module implements plotting functions useful to report analysis results.\n\nAuthor: Martin Perez-Guevara, Elvis Dohmatob, 2017\n\"\"\"\n\nimport warnings\nfrom string import ascii_lowercase\n\nimport numpy as np\nimport pandas as pd\nimport nibabel as nib\nfrom scipy import ndimage\n\nfrom nilearn.image import get_data\nfrom nilearn.image.resampling import coord_transform\nfrom nilearn._utils import check_niimg_3d\n\n\ndef _local_max(data, affine, min_distance):\n    \"\"\"Find all local maxima of the array, separated by at least min_distance.\n    Adapted from https://stackoverflow.com/a/22631583/2589328\n\n    Parameters\n    ----------\n    data : array_like\n        3D array of with masked values for cluster.\n\n    affine : np.ndarray\n        Square matrix specifying the position of the image array data\n        in a reference space.\n\n    min_distance : int\n        Minimum distance between local maxima in ``data``, in terms of mm.\n\n    Returns\n    -------\n    ijk : `numpy.ndarray`\n        (n_foci, 3) array of local maxima indices for cluster.\n\n    vals : `numpy.ndarray`\n        (n_foci,) array of values from data at ijk.\n\n    \"\"\"\n    ijk, vals = _identify_subpeaks(data)\n    xyz, ijk, vals = _sort_subpeaks(ijk, vals, affine)\n    ijk, vals = _pare_subpeaks(xyz, ijk, vals, min_distance)\n    return ijk, vals\n\n\ndef _identify_subpeaks(data):\n    \"\"\"Identify cluster peak and subpeaks based on minimum distance.\n\n    Parameters\n    ----------\n    data : `numpy.ndarray`\n        3D array of with masked values for cluster.\n\n    Returns\n    -------\n    ijk : `numpy.ndarray`\n        (n_foci, 3) array of local maxima indices for cluster.\n    vals : `numpy.ndarray`\n        (n_foci,) array of values from data at ijk.\n\n    Notes\n    -----\n    When a cluster's local maxima correspond to contiguous voxels with the\n    same values (as in a binary cluster), this function determines the center\n    of mass for those voxels.\n    \"\"\"\n    # Initial identification of subpeaks with minimal minimum distance\n    data_max = ndimage.filters.maximum_filter(data, 3)\n    maxima = data == data_max\n    data_min = ndimage.filters.minimum_filter(data, 3)\n    diff = (data_max - data_min) > 0\n    maxima[diff == 0] = 0\n\n    labeled, n_subpeaks = ndimage.label(maxima)\n    labels_index = range(1, n_subpeaks + 1)\n    ijk = np.array(ndimage.center_of_mass(data, labeled, labels_index))\n    ijk = np.round(ijk).astype(int)\n    vals = np.apply_along_axis(\n        arr=ijk, axis=1, func1d=_get_val, input_arr=data\n    )\n    # Determine if all subpeaks are within the cluster\n    # They may not be if the cluster is binary and has a shape where the COM is\n    # outside the cluster, like a donut.\n    cluster_idx = np.vstack(np.where(labeled)).T.tolist()\n    subpeaks_outside_cluster = [\n        i\n        for i, peak_idx in enumerate(ijk.tolist())\n        if peak_idx not in cluster_idx\n    ]\n    vals[subpeaks_outside_cluster] = np.nan\n    if subpeaks_outside_cluster:\n        warnings.warn(\n            \"Attention: At least one of the (sub)peaks falls outside of the \"\n            \"cluster body.\"\n        )\n    return ijk, vals\n\n\ndef _sort_subpeaks(ijk, vals, affine):\n    # Sort subpeaks in cluster in descending order of stat value\n    order = (-vals).argsort()\n    vals = vals[order]\n    ijk = ijk[order, :]\n    xyz = nib.affines.apply_affine(affine, ijk)  # Convert to xyz in mm\n    return xyz, ijk, vals\n\n\ndef _pare_subpeaks(xyz, ijk, vals, min_distance):\n    # Reduce list of subpeaks based on distance\n    keep_idx = np.ones(xyz.shape[0]).astype(bool)\n    for i in range(xyz.shape[0]):\n        for j in range(i + 1, xyz.shape[0]):\n            if keep_idx[i] == 1:\n                dist = np.linalg.norm(xyz[i, :] - xyz[j, :])\n                keep_idx[j] = dist > min_distance\n    ijk = ijk[keep_idx, :]\n    vals = vals[keep_idx]\n    return ijk, vals\n\n\ndef _get_val(row, input_arr):\n    \"\"\"Small function for extracting values from array based on index.\n    \"\"\"\n    i, j, k = row\n    return input_arr[i, j, k]\n\n\ndef get_clusters_table(stat_img, stat_threshold, cluster_threshold=None,\n                       two_sided=False, min_distance=8.):\n    \"\"\"Creates pandas dataframe with img cluster statistics.\n\n    Parameters\n    ----------\n    stat_img : Niimg-like object,\n       Statistical image (presumably in z- or p-scale).\n\n    stat_threshold : `float`\n        Cluster forming threshold in same scale as `stat_img` (either a\n        p-value or z-scale value).\n\n    cluster_threshold : `int` or `None`, optional\n        Cluster size threshold, in voxels.\n\n    two_sided : `bool`, optional\n        Whether to employ two-sided thresholding or to evaluate positive values\n        only. Default=False.\n\n    min_distance : `float`, optional\n        Minimum distance between subpeaks in mm. Default=8mm.\n\n    Returns\n    -------\n    df : `pandas.DataFrame`\n        Table with peaks and subpeaks from thresholded `stat_img`. For binary\n        clusters (clusters with >1 voxel containing only one value), the table\n        reports the center of mass of the cluster,\n        rather than any peaks/subpeaks.\n\n    \"\"\"\n    cols = ['Cluster ID', 'X', 'Y', 'Z', 'Peak Stat', 'Cluster Size (mm3)']\n\n    # check that stat_img is niimg-like object and 3D\n    stat_img = check_niimg_3d(stat_img)\n    # If cluster threshold is used, there is chance that stat_map will be\n    # modified, therefore copy is needed\n    if cluster_threshold is None:\n        stat_map = get_data(stat_img)\n    else:\n        stat_map = get_data(stat_img).copy()\n\n    # Define array for 6-connectivity, aka NN1 or \"faces\"\n    conn_mat = np.zeros((3, 3, 3), int)\n    conn_mat[1, 1, :] = 1\n    conn_mat[1, :, 1] = 1\n    conn_mat[:, 1, 1] = 1\n    voxel_size = np.prod(stat_img.header.get_zooms())\n\n    signs = [1, -1] if two_sided else [1]\n    no_clusters_found = True\n    rows = []\n    for sign in signs:\n        # Flip map if necessary\n        temp_stat_map = stat_map * sign\n\n        # Binarize using CDT\n        binarized = temp_stat_map > stat_threshold\n        binarized = binarized.astype(int)\n\n        # If the stat threshold is too high simply return an empty dataframe\n        if np.sum(binarized) == 0:\n            warnings.warn(\n                'Attention: No clusters with stat {0} than {1}'.format(\n                    'higher' if sign == 1 else 'lower',\n                    stat_threshold * sign,\n                )\n            )\n            continue\n\n        # Extract connected components above cluster size threshold\n        label_map = ndimage.measurements.label(binarized, conn_mat)[0]\n        clust_ids = sorted(list(np.unique(label_map)[1:]))\n        for c_val in clust_ids:\n            if cluster_threshold is not None and np.sum(\n                    label_map == c_val) < cluster_threshold:\n                temp_stat_map[label_map == c_val] = 0\n                binarized[label_map == c_val] = 0\n\n        # If the cluster threshold is too high simply return an empty dataframe\n        # this checks for stats higher than threshold after small clusters\n        # were removed from temp_stat_map\n        if np.sum(temp_stat_map > stat_threshold) == 0:\n            warnings.warn(\n                'Attention: No clusters with more than {0} voxels'.format(\n                    cluster_threshold,\n                )\n            )\n            continue\n\n        # Now re-label and create table\n        label_map = ndimage.measurements.label(binarized, conn_mat)[0]\n        clust_ids = sorted(list(np.unique(label_map)[1:]))\n        peak_vals = np.array(\n            [np.max(temp_stat_map * (label_map == c)) for c in clust_ids])\n        # Sort by descending max value\n        clust_ids = [clust_ids[c] for c in (-peak_vals).argsort()]\n\n        for c_id, c_val in enumerate(clust_ids):\n            cluster_mask = label_map == c_val\n            masked_data = temp_stat_map * cluster_mask\n\n            cluster_size_mm = int(np.sum(cluster_mask) * voxel_size)\n\n            # Get peaks, subpeaks and associated statistics\n            subpeak_ijk, subpeak_vals = _local_max(\n                masked_data,\n                stat_img.affine,\n                min_distance=min_distance,\n            )\n            subpeak_vals *= sign  # flip signs if necessary\n            subpeak_xyz = np.asarray(\n                coord_transform(\n                    subpeak_ijk[:, 0],\n                    subpeak_ijk[:, 1],\n                    subpeak_ijk[:, 2],\n                    stat_img.affine,\n                )\n            ).tolist()\n            subpeak_xyz = np.array(subpeak_xyz).T\n\n            # Only report peak and, at most, top 3 subpeaks.\n            n_subpeaks = np.min((len(subpeak_vals), 4))\n            for subpeak in range(n_subpeaks):\n                if subpeak == 0:\n                    row = [\n                        c_id + 1,\n                        subpeak_xyz[subpeak, 0],\n                        subpeak_xyz[subpeak, 1],\n                        subpeak_xyz[subpeak, 2],\n                        subpeak_vals[subpeak],\n                        cluster_size_mm,\n                    ]\n                else:\n                    # Subpeak naming convention is cluster num+letter:\n                    # 1a, 1b, etc\n                    sp_id = '{0}{1}'.format(\n                        c_id + 1,\n                        ascii_lowercase[subpeak - 1],\n                    )\n                    row = [\n                        sp_id,\n                        subpeak_xyz[subpeak, 0],\n                        subpeak_xyz[subpeak, 1],\n                        subpeak_xyz[subpeak, 2],\n                        subpeak_vals[subpeak],\n                        '',\n                    ]\n                rows += [row]\n\n        # If we reach this point, there are clusters in this sign\n        no_clusters_found = False\n\n    if no_clusters_found:\n        df = pd.DataFrame(columns=cols)\n    else:\n        df = pd.DataFrame(columns=cols, data=rows)\n\n    return df\n", "path": "nilearn/reporting/_get_clusters_table.py"}], "after_files": [{"content": "\"\"\"\nThis module implements plotting functions useful to report analysis results.\n\nAuthor: Martin Perez-Guevara, Elvis Dohmatob, 2017\n\"\"\"\n\nimport warnings\nfrom string import ascii_lowercase\n\nimport numpy as np\nimport pandas as pd\nimport nibabel as nib\nfrom scipy import ndimage\n\nfrom nilearn._utils.niimg import _safe_get_data\nfrom nilearn.image.resampling import coord_transform\nfrom nilearn._utils import check_niimg_3d\n\n\ndef _local_max(data, affine, min_distance):\n    \"\"\"Find all local maxima of the array, separated by at least min_distance.\n    Adapted from https://stackoverflow.com/a/22631583/2589328\n\n    Parameters\n    ----------\n    data : array_like\n        3D array of with masked values for cluster.\n\n    affine : np.ndarray\n        Square matrix specifying the position of the image array data\n        in a reference space.\n\n    min_distance : int\n        Minimum distance between local maxima in ``data``, in terms of mm.\n\n    Returns\n    -------\n    ijk : `numpy.ndarray`\n        (n_foci, 3) array of local maxima indices for cluster.\n\n    vals : `numpy.ndarray`\n        (n_foci,) array of values from data at ijk.\n\n    \"\"\"\n    ijk, vals = _identify_subpeaks(data)\n    xyz, ijk, vals = _sort_subpeaks(ijk, vals, affine)\n    ijk, vals = _pare_subpeaks(xyz, ijk, vals, min_distance)\n    return ijk, vals\n\n\ndef _identify_subpeaks(data):\n    \"\"\"Identify cluster peak and subpeaks based on minimum distance.\n\n    Parameters\n    ----------\n    data : `numpy.ndarray`\n        3D array of with masked values for cluster.\n\n    Returns\n    -------\n    ijk : `numpy.ndarray`\n        (n_foci, 3) array of local maxima indices for cluster.\n    vals : `numpy.ndarray`\n        (n_foci,) array of values from data at ijk.\n\n    Notes\n    -----\n    When a cluster's local maxima correspond to contiguous voxels with the\n    same values (as in a binary cluster), this function determines the center\n    of mass for those voxels.\n    \"\"\"\n    # Initial identification of subpeaks with minimal minimum distance\n    data_max = ndimage.filters.maximum_filter(data, 3)\n    maxima = data == data_max\n    data_min = ndimage.filters.minimum_filter(data, 3)\n    diff = (data_max - data_min) > 0\n    maxima[diff == 0] = 0\n\n    labeled, n_subpeaks = ndimage.label(maxima)\n    labels_index = range(1, n_subpeaks + 1)\n    ijk = np.array(ndimage.center_of_mass(data, labeled, labels_index))\n    ijk = np.round(ijk).astype(int)\n    vals = np.apply_along_axis(\n        arr=ijk, axis=1, func1d=_get_val, input_arr=data\n    )\n    # Determine if all subpeaks are within the cluster\n    # They may not be if the cluster is binary and has a shape where the COM is\n    # outside the cluster, like a donut.\n    cluster_idx = np.vstack(np.where(labeled)).T.tolist()\n    subpeaks_outside_cluster = [\n        i\n        for i, peak_idx in enumerate(ijk.tolist())\n        if peak_idx not in cluster_idx\n    ]\n    vals[subpeaks_outside_cluster] = np.nan\n    if subpeaks_outside_cluster:\n        warnings.warn(\n            \"Attention: At least one of the (sub)peaks falls outside of the \"\n            \"cluster body.\"\n        )\n    return ijk, vals\n\n\ndef _sort_subpeaks(ijk, vals, affine):\n    # Sort subpeaks in cluster in descending order of stat value\n    order = (-vals).argsort()\n    vals = vals[order]\n    ijk = ijk[order, :]\n    xyz = nib.affines.apply_affine(affine, ijk)  # Convert to xyz in mm\n    return xyz, ijk, vals\n\n\ndef _pare_subpeaks(xyz, ijk, vals, min_distance):\n    # Reduce list of subpeaks based on distance\n    keep_idx = np.ones(xyz.shape[0]).astype(bool)\n    for i in range(xyz.shape[0]):\n        for j in range(i + 1, xyz.shape[0]):\n            if keep_idx[i] == 1:\n                dist = np.linalg.norm(xyz[i, :] - xyz[j, :])\n                keep_idx[j] = dist > min_distance\n    ijk = ijk[keep_idx, :]\n    vals = vals[keep_idx]\n    return ijk, vals\n\n\ndef _get_val(row, input_arr):\n    \"\"\"Small function for extracting values from array based on index.\n    \"\"\"\n    i, j, k = row\n    return input_arr[i, j, k]\n\n\ndef get_clusters_table(stat_img, stat_threshold, cluster_threshold=None,\n                       two_sided=False, min_distance=8.):\n    \"\"\"Creates pandas dataframe with img cluster statistics.\n\n    Parameters\n    ----------\n    stat_img : Niimg-like object,\n       Statistical image (presumably in z- or p-scale).\n\n    stat_threshold : `float`\n        Cluster forming threshold in same scale as `stat_img` (either a\n        p-value or z-scale value).\n\n    cluster_threshold : `int` or `None`, optional\n        Cluster size threshold, in voxels.\n\n    two_sided : `bool`, optional\n        Whether to employ two-sided thresholding or to evaluate positive values\n        only. Default=False.\n\n    min_distance : `float`, optional\n        Minimum distance between subpeaks in mm. Default=8mm.\n\n    Returns\n    -------\n    df : `pandas.DataFrame`\n        Table with peaks and subpeaks from thresholded `stat_img`. For binary\n        clusters (clusters with >1 voxel containing only one value), the table\n        reports the center of mass of the cluster,\n        rather than any peaks/subpeaks.\n\n    \"\"\"\n    cols = ['Cluster ID', 'X', 'Y', 'Z', 'Peak Stat', 'Cluster Size (mm3)']\n\n    # check that stat_img is niimg-like object and 3D\n    stat_img = check_niimg_3d(stat_img)\n    # If cluster threshold is used, there is chance that stat_map will be\n    # modified, therefore copy is needed\n    stat_map = _safe_get_data(stat_img, ensure_finite=True,\n                              copy_data=(cluster_threshold is not None))\n\n    # Define array for 6-connectivity, aka NN1 or \"faces\"\n    conn_mat = np.zeros((3, 3, 3), int)\n    conn_mat[1, 1, :] = 1\n    conn_mat[1, :, 1] = 1\n    conn_mat[:, 1, 1] = 1\n    voxel_size = np.prod(stat_img.header.get_zooms())\n\n    signs = [1, -1] if two_sided else [1]\n    no_clusters_found = True\n    rows = []\n    for sign in signs:\n        # Flip map if necessary\n        temp_stat_map = stat_map * sign\n\n        # Binarize using CDT\n        binarized = temp_stat_map > stat_threshold\n        binarized = binarized.astype(int)\n\n        # If the stat threshold is too high simply return an empty dataframe\n        if np.sum(binarized) == 0:\n            warnings.warn(\n                'Attention: No clusters with stat {0} than {1}'.format(\n                    'higher' if sign == 1 else 'lower',\n                    stat_threshold * sign,\n                )\n            )\n            continue\n\n        # Extract connected components above cluster size threshold\n        label_map = ndimage.measurements.label(binarized, conn_mat)[0]\n        clust_ids = sorted(list(np.unique(label_map)[1:]))\n        for c_val in clust_ids:\n            if cluster_threshold is not None and np.sum(\n                    label_map == c_val) < cluster_threshold:\n                temp_stat_map[label_map == c_val] = 0\n                binarized[label_map == c_val] = 0\n\n        # If the cluster threshold is too high simply return an empty dataframe\n        # this checks for stats higher than threshold after small clusters\n        # were removed from temp_stat_map\n        if np.sum(temp_stat_map > stat_threshold) == 0:\n            warnings.warn(\n                'Attention: No clusters with more than {0} voxels'.format(\n                    cluster_threshold,\n                )\n            )\n            continue\n\n        # Now re-label and create table\n        label_map = ndimage.measurements.label(binarized, conn_mat)[0]\n        clust_ids = sorted(list(np.unique(label_map)[1:]))\n        peak_vals = np.array(\n            [np.max(temp_stat_map * (label_map == c)) for c in clust_ids])\n        # Sort by descending max value\n        clust_ids = [clust_ids[c] for c in (-peak_vals).argsort()]\n\n        for c_id, c_val in enumerate(clust_ids):\n            cluster_mask = label_map == c_val\n            masked_data = temp_stat_map * cluster_mask\n\n            cluster_size_mm = int(np.sum(cluster_mask) * voxel_size)\n\n            # Get peaks, subpeaks and associated statistics\n            subpeak_ijk, subpeak_vals = _local_max(\n                masked_data,\n                stat_img.affine,\n                min_distance=min_distance,\n            )\n            subpeak_vals *= sign  # flip signs if necessary\n            subpeak_xyz = np.asarray(\n                coord_transform(\n                    subpeak_ijk[:, 0],\n                    subpeak_ijk[:, 1],\n                    subpeak_ijk[:, 2],\n                    stat_img.affine,\n                )\n            ).tolist()\n            subpeak_xyz = np.array(subpeak_xyz).T\n\n            # Only report peak and, at most, top 3 subpeaks.\n            n_subpeaks = np.min((len(subpeak_vals), 4))\n            for subpeak in range(n_subpeaks):\n                if subpeak == 0:\n                    row = [\n                        c_id + 1,\n                        subpeak_xyz[subpeak, 0],\n                        subpeak_xyz[subpeak, 1],\n                        subpeak_xyz[subpeak, 2],\n                        subpeak_vals[subpeak],\n                        cluster_size_mm,\n                    ]\n                else:\n                    # Subpeak naming convention is cluster num+letter:\n                    # 1a, 1b, etc\n                    sp_id = '{0}{1}'.format(\n                        c_id + 1,\n                        ascii_lowercase[subpeak - 1],\n                    )\n                    row = [\n                        sp_id,\n                        subpeak_xyz[subpeak, 0],\n                        subpeak_xyz[subpeak, 1],\n                        subpeak_xyz[subpeak, 2],\n                        subpeak_vals[subpeak],\n                        '',\n                    ]\n                rows += [row]\n\n        # If we reach this point, there are clusters in this sign\n        no_clusters_found = False\n\n    if no_clusters_found:\n        df = pd.DataFrame(columns=cols)\n    else:\n        df = pd.DataFrame(columns=cols, data=rows)\n\n    return df\n", "path": "nilearn/reporting/_get_clusters_table.py"}]}