We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
nistats: dependence between function calls – strange bug may affect reproducibility
Hi,

I've encountered very strange behavior of two functions – `map_threshold` from `nistats.thresholding` module and `get_clusters_table` from `nistats.reporting` module. It seems like the threshold value returned by `map_threshold` depends on the last output value(s) of `get_clusters_table` function or some internal function variables. Impact seems to be severe (the higher the `cluster_threshold` value passed to `get_cluster_table`, the higher the threshold returned by subsequent call of `map_threshold`). 

Here is a simple demonstration: https://github.com/kbonna/decidenet/blob/master/activation_analysis/nistats_map_threshold_bug.ipynb

I've tried to find a source of the problem studying code of both functions, but I everything seems to be ok at first glance. 

Am I missing something simple?

--- 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     ijk, vals = _identify_subpeaks(data)
44     xyz, ijk, vals = _sort_subpeaks(ijk, vals, affine)
45     ijk, vals = _pare_subpeaks(xyz, ijk, vals, min_distance)
46     return ijk, vals
47 
48 
49 def _identify_subpeaks(data):
50     # Initial identification of subpeaks with minimal minimum distance
51     data_max = ndimage.filters.maximum_filter(data, 3)
52     maxima = (data == data_max)
53     data_min = ndimage.filters.minimum_filter(data, 3)
54     diff = ((data_max - data_min) > 0)
55     maxima[diff == 0] = 0
56 
57     labeled, n_subpeaks = ndimage.label(maxima)
58     labels_index = range(1, n_subpeaks + 1)
59     ijk = np.array(ndimage.center_of_mass(data, labeled, labels_index))
60     ijk = np.round(ijk).astype(int)
61     vals = np.apply_along_axis(arr=ijk, axis=1, func1d=_get_val,
62                                input_arr=data)
63     return ijk, vals
64 
65 
66 def _sort_subpeaks(ijk, vals, affine):
67     # Sort subpeaks in cluster in descending order of stat value
68     order = (-vals).argsort()
69     vals = vals[order]
70     ijk = ijk[order, :]
71     xyz = nib.affines.apply_affine(affine, ijk)  # Convert to xyz in mm
72     return xyz, ijk, vals
73 
74 
75 def _pare_subpeaks(xyz, ijk, vals, min_distance):
76     # Reduce list of subpeaks based on distance
77     keep_idx = np.ones(xyz.shape[0]).astype(bool)
78     for i in range(xyz.shape[0]):
79         for j in range(i + 1, xyz.shape[0]):
80             if keep_idx[i] == 1:
81                 dist = np.linalg.norm(xyz[i, :] - xyz[j, :])
82                 keep_idx[j] = dist > min_distance
83     ijk = ijk[keep_idx, :]
84     vals = vals[keep_idx]
85     return ijk, vals
86 
87 
88 def _get_val(row, input_arr):
89     """Small function for extracting values from array based on index.
90     """
91     i, j, k = row
92     return input_arr[i, j, k]
93 
94 
95 def get_clusters_table(stat_img, stat_threshold, cluster_threshold=None,
96                        min_distance=8.):
97     """Creates pandas dataframe with img cluster statistics.
98 
99     Parameters
100     ----------
101     stat_img : Niimg-like object,
102        Statistical image (presumably in z- or p-scale).
103 
104     stat_threshold: `float`
105         Cluster forming threshold in same scale as `stat_img` (either a
106         p-value or z-scale value).
107 
108     cluster_threshold : `int` or `None`, optional
109         Cluster size threshold, in voxels.
110 
111     min_distance: `float`, optional
112         Minimum distance between subpeaks in mm. Default is 8 mm.
113 
114     Returns
115     -------
116     df : `pandas.DataFrame`
117         Table with peaks and subpeaks from thresholded `stat_img`. For binary
118         clusters (clusters with >1 voxel containing only one value), the table
119         reports the center of mass of the cluster,
120         rather than any peaks/subpeaks.
121     """
122     cols = ['Cluster ID', 'X', 'Y', 'Z', 'Peak Stat', 'Cluster Size (mm3)']
123     stat_map = get_data(stat_img)
124     conn_mat = np.zeros((3, 3, 3), int)  # 6-connectivity, aka NN1 or "faces"
125     conn_mat[1, 1, :] = 1
126     conn_mat[1, :, 1] = 1
127     conn_mat[:, 1, 1] = 1
128     voxel_size = np.prod(stat_img.header.get_zooms())
129 
130     # Binarize using CDT
131     binarized = stat_map > stat_threshold
132     binarized = binarized.astype(int)
133 
134     # If the stat threshold is too high simply return an empty dataframe
135     if np.sum(binarized) == 0:
136         warnings.warn('Attention: No clusters with stat higher than %f' %
137                       stat_threshold)
138         return pd.DataFrame(columns=cols)
139 
140     # Extract connected components above cluster size threshold
141     label_map = ndimage.measurements.label(binarized, conn_mat)[0]
142     clust_ids = sorted(list(np.unique(label_map)[1:]))
143     for c_val in clust_ids:
144         if cluster_threshold is not None and np.sum(
145             label_map == c_val) < cluster_threshold:
146             stat_map[label_map == c_val] = 0
147             binarized[label_map == c_val] = 0
148 
149     # If the cluster threshold is too high simply return an empty dataframe
150     # this checks for stats higher than threshold after small clusters
151     # were removed from stat_map
152     if np.sum(stat_map > stat_threshold) == 0:
153         warnings.warn('Attention: No clusters with more than %d voxels' %
154                       cluster_threshold)
155         return pd.DataFrame(columns=cols)
156 
157     # Now re-label and create table
158     label_map = ndimage.measurements.label(binarized, conn_mat)[0]
159     clust_ids = sorted(list(np.unique(label_map)[1:]))
160     peak_vals = np.array(
161         [np.max(stat_map * (label_map == c)) for c in clust_ids])
162     clust_ids = [clust_ids[c] for c in
163                  (-peak_vals).argsort()]  # Sort by descending max value
164 
165     rows = []
166     for c_id, c_val in enumerate(clust_ids):
167         cluster_mask = label_map == c_val
168         masked_data = stat_map * cluster_mask
169 
170         cluster_size_mm = int(np.sum(cluster_mask) * voxel_size)
171 
172         # Get peaks, subpeaks and associated statistics
173         subpeak_ijk, subpeak_vals = _local_max(masked_data, stat_img.affine,
174                                                min_distance=min_distance)
175         subpeak_xyz = np.asarray(coord_transform(subpeak_ijk[:, 0],
176                                                  subpeak_ijk[:, 1],
177                                                  subpeak_ijk[:, 2],
178                                                  stat_img.affine)).tolist()
179         subpeak_xyz = np.array(subpeak_xyz).T
180 
181         # Only report peak and, at most, top 3 subpeaks.
182         n_subpeaks = np.min((len(subpeak_vals), 4))
183         for subpeak in range(n_subpeaks):
184             if subpeak == 0:
185                 row = [c_id + 1, subpeak_xyz[subpeak, 0],
186                        subpeak_xyz[subpeak, 1], subpeak_xyz[subpeak, 2],
187                        subpeak_vals[subpeak], cluster_size_mm]
188             else:
189                 # Subpeak naming convention is cluster num+letter: 1a, 1b, etc
190                 sp_id = '{0}{1}'.format(c_id + 1, ascii_lowercase[subpeak - 1])
191                 row = [sp_id, subpeak_xyz[subpeak, 0], subpeak_xyz[subpeak, 1],
192                        subpeak_xyz[subpeak, 2], subpeak_vals[subpeak], '']
193             rows += [row]
194     df = pd.DataFrame(columns=cols, data=rows)
195     return df
196
```
--- 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
@@ -120,7 +120,14 @@
         rather than any peaks/subpeaks.
     """
     cols = ['Cluster ID', 'X', 'Y', 'Z', 'Peak Stat', 'Cluster Size (mm3)']
-    stat_map = get_data(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()
+
     conn_mat = np.zeros((3, 3, 3), int)  # 6-connectivity, aka NN1 or "faces"
     conn_mat[1, 1, :] = 1
     conn_mat[1, :, 1] = 1

{"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@@ -120,7 +120,14 @@\n         rather than any peaks/subpeaks.\n     \"\"\"\n     cols = ['Cluster ID', 'X', 'Y', 'Z', 'Peak Stat', 'Cluster Size (mm3)']\n-    stat_map = get_data(stat_img)\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", "issue": "nistats: dependence between function calls \u2013 strange bug may affect reproducibility\nHi,\r\n\r\nI've encountered very strange behavior of two functions \u2013 `map_threshold` from `nistats.thresholding` module and `get_clusters_table` from `nistats.reporting` module. It seems like the threshold value returned by `map_threshold` depends on the last output value(s) of `get_clusters_table` function or some internal function variables. Impact seems to be severe (the higher the `cluster_threshold` value passed to `get_cluster_table`, the higher the threshold returned by subsequent call of `map_threshold`). \r\n\r\nHere is a simple demonstration: https://github.com/kbonna/decidenet/blob/master/activation_analysis/nistats_map_threshold_bug.ipynb\r\n\r\nI've tried to find a source of the problem studying code of both functions, but I everything seems to be ok at first glance. \r\n\r\nAm I missing something simple?\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    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    # 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(arr=ijk, axis=1, func1d=_get_val,\n                               input_arr=data)\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 is 8 mm.\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    cols = ['Cluster ID', 'X', 'Y', 'Z', 'Peak Stat', 'Cluster Size (mm3)']\n    stat_map = get_data(stat_img)\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\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    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    # 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(arr=ijk, axis=1, func1d=_get_val,\n                               input_arr=data)\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 is 8 mm.\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    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"}]}