We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
get_clusters_table does not allow file input and does not check that image is 3D
<!--Provide a brief description of the bug.-->
get_clusters_table does not let me pass in a filename, and does not ensure the image is 3D, which can raise an error of there is an "extra dimension" in the data.

<!--Please fill in the following information, to the best of your ability.-->
Nilearn version: 0.7.0

### Expected behavior

`get_clusters_table` should take filenames as input and ensure/check inputs are 3D

### Actual behavior

raises error if passed filename or 4D input (but effectively 3D) is passed.

### Steps and code to reproduce bug

```python
>>> from nilearn.reporting import get_clusters_table
>>> import numpy as np
>>> import nibabel as nib
>>> shape = (9, 10, 11)
>>> data = np.zeros(shape)
>>> data[2:4, 5:7, 6:8] = 5.
>>> stat_img = nib.Nifti1Image(data[...,np.newaxis], np.eye(4))
>>> get_clusters_table(stat_img, 0.0)

```

Error:
```
  File "/home/jdkent/projects/nilearn/nilearn/reporting/_get_clusters_table.py", line 186, in get_clusters_table
    label_map = ndimage.measurements.label(binarized, conn_mat)[0]
  File "/home/jdkent/projects/nilearn/.nilearn/lib/python3.8/site-packages/scipy/ndimage/measurements.py", line 183, in label
    raise RuntimeError('structure and input must have equal rank')
RuntimeError: structure and input must have equal rank
```


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

{"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@@ -14,6 +14,7 @@\n \n from nilearn.image import get_data\n from nilearn.image.resampling import coord_transform\n+from nilearn._utils import check_niimg_3d\n \n \n def _local_max(data, affine, min_distance):\n@@ -159,6 +160,8 @@\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", "issue": "get_clusters_table does not allow file input and does not check that image is 3D\n<!--Provide a brief description of the bug.-->\r\nget_clusters_table does not let me pass in a filename, and does not ensure the image is 3D, which can raise an error of there is an \"extra dimension\" in the data.\r\n\r\n<!--Please fill in the following information, to the best of your ability.-->\r\nNilearn version: 0.7.0\r\n\r\n### Expected behavior\r\n\r\n`get_clusters_table` should take filenames as input and ensure/check inputs are 3D\r\n\r\n### Actual behavior\r\n\r\nraises error if passed filename or 4D input (but effectively 3D) is passed.\r\n\r\n### Steps and code to reproduce bug\r\n\r\n```python\r\n>>> from nilearn.reporting import get_clusters_table\r\n>>> import numpy as np\r\n>>> import nibabel as nib\r\n>>> shape = (9, 10, 11)\r\n>>> data = np.zeros(shape)\r\n>>> data[2:4, 5:7, 6:8] = 5.\r\n>>> stat_img = nib.Nifti1Image(data[...,np.newaxis], np.eye(4))\r\n>>> get_clusters_table(stat_img, 0.0)\r\n\r\n```\r\n\r\nError:\r\n```\r\n  File \"/home/jdkent/projects/nilearn/nilearn/reporting/_get_clusters_table.py\", line 186, in get_clusters_table\r\n    label_map = ndimage.measurements.label(binarized, conn_mat)[0]\r\n  File \"/home/jdkent/projects/nilearn/.nilearn/lib/python3.8/site-packages/scipy/ndimage/measurements.py\", line 183, in label\r\n    raise RuntimeError('structure and input must have equal rank')\r\nRuntimeError: structure and input must have equal rank\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\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                       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    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    # 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    conn_mat = np.zeros((3, 3, 3), int)  # 6-connectivity, aka NN1 or \"faces\"\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    # Binarize using CDT\n    binarized = 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('Attention: No clusters with stat higher than %f' %\n                      stat_threshold)\n        return pd.DataFrame(columns=cols)\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            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 stat_map\n    if np.sum(stat_map > stat_threshold) == 0:\n        warnings.warn('Attention: No clusters with more than %d voxels' %\n                      cluster_threshold)\n        return pd.DataFrame(columns=cols)\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(stat_map * (label_map == c)) for c in clust_ids])\n    clust_ids = [clust_ids[c] for c in\n                 (-peak_vals).argsort()]  # Sort by descending max value\n\n    rows = []\n    for c_id, c_val in enumerate(clust_ids):\n        cluster_mask = label_map == c_val\n        masked_data = 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(masked_data, stat_img.affine,\n                                               min_distance=min_distance)\n        subpeak_xyz = np.asarray(coord_transform(subpeak_ijk[:, 0],\n                                                 subpeak_ijk[:, 1],\n                                                 subpeak_ijk[:, 2],\n                                                 stat_img.affine)).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 = [c_id + 1, subpeak_xyz[subpeak, 0],\n                       subpeak_xyz[subpeak, 1], subpeak_xyz[subpeak, 2],\n                       subpeak_vals[subpeak], cluster_size_mm]\n            else:\n                # Subpeak naming convention is cluster num+letter: 1a, 1b, etc\n                sp_id = '{0}{1}'.format(c_id + 1, ascii_lowercase[subpeak - 1])\n                row = [sp_id, subpeak_xyz[subpeak, 0], subpeak_xyz[subpeak, 1],\n                       subpeak_xyz[subpeak, 2], subpeak_vals[subpeak], '']\n            rows += [row]\n    df = pd.DataFrame(columns=cols, data=rows)\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.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                       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    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    conn_mat = np.zeros((3, 3, 3), int)  # 6-connectivity, aka NN1 or \"faces\"\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    # Binarize using CDT\n    binarized = 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('Attention: No clusters with stat higher than %f' %\n                      stat_threshold)\n        return pd.DataFrame(columns=cols)\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            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 stat_map\n    if np.sum(stat_map > stat_threshold) == 0:\n        warnings.warn('Attention: No clusters with more than %d voxels' %\n                      cluster_threshold)\n        return pd.DataFrame(columns=cols)\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(stat_map * (label_map == c)) for c in clust_ids])\n    clust_ids = [clust_ids[c] for c in\n                 (-peak_vals).argsort()]  # Sort by descending max value\n\n    rows = []\n    for c_id, c_val in enumerate(clust_ids):\n        cluster_mask = label_map == c_val\n        masked_data = 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(masked_data, stat_img.affine,\n                                               min_distance=min_distance)\n        subpeak_xyz = np.asarray(coord_transform(subpeak_ijk[:, 0],\n                                                 subpeak_ijk[:, 1],\n                                                 subpeak_ijk[:, 2],\n                                                 stat_img.affine)).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 = [c_id + 1, subpeak_xyz[subpeak, 0],\n                       subpeak_xyz[subpeak, 1], subpeak_xyz[subpeak, 2],\n                       subpeak_vals[subpeak], cluster_size_mm]\n            else:\n                # Subpeak naming convention is cluster num+letter: 1a, 1b, etc\n                sp_id = '{0}{1}'.format(c_id + 1, ascii_lowercase[subpeak - 1])\n                row = [sp_id, subpeak_xyz[subpeak, 0], subpeak_xyz[subpeak, 1],\n                       subpeak_xyz[subpeak, 2], subpeak_vals[subpeak], '']\n            rows += [row]\n    df = pd.DataFrame(columns=cols, data=rows)\n    return df\n", "path": "nilearn/reporting/_get_clusters_table.py"}]}