We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
find_xyz_cut_coords crashes on some compact binary images
...trace back is:

``` python
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-44-802cef648458> in <module>()
----> 1 plotting.find_xyz_cut_coords(nib.load(roi_paths[2]))

/git/nilearn/nilearn/plotting/find_cuts.py in find_xyz_cut_coords(img, mask, activation_threshold)
     90         return .5 * np.array(data.shape)
     91     mask = largest_connected_component(mask)
---> 92     slice_x, slice_y, slice_z = ndimage.find_objects(mask)[0]
     93     my_map = my_map[slice_x, slice_y, slice_z]
     94     mask = mask[slice_x, slice_y, slice_z]
```

The image on which this crashes is attached.

[dmpfc-3.nii.zip](https://github.com/nilearn/nilearn/files/66476/dmpfc-3.nii.zip)


--- 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/plotting/find_cuts.py`
Content:
```
1 """
2 Tools to find activations and cut on maps
3 """
4 
5 # Author: Gael Varoquaux
6 # License: BSD
7 
8 import warnings
9 import numbers
10 import numpy as np
11 from scipy import ndimage
12 
13 # Local imports
14 from .._utils.ndimage import largest_connected_component
15 from ..image import new_img_like
16 from .._utils.extmath import fast_abs_percentile
17 from .._utils.numpy_conversions import as_ndarray
18 from ..image.resampling import get_mask_bounds, coord_transform
19 from ..image.image import _smooth_array
20 
21 ################################################################################
22 # Functions for automatic choice of cuts coordinates
23 ################################################################################
24 
25 
26 def find_xyz_cut_coords(img, mask=None, activation_threshold=None):
27     """ Find the center of the largest activation connected component.
28 
29         Parameters
30         -----------
31         img : 3D Nifti1Image
32             The brain map.
33         mask : 3D ndarray, boolean, optional
34             An optional brain mask.
35         activation_threshold : float, optional
36             The lower threshold to the positive activation. If None, the
37             activation threshold is computed using the 80% percentile of
38             the absolute value of the map.
39 
40         Returns
41         -------
42         x : float
43             the x world coordinate.
44         y : float
45             the y world coordinate.
46         z : float
47             the z world coordinate.
48     """
49     data = img.get_data()
50     # To speed up computations, we work with partial views of the array,
51     # and keep track of the offset
52     offset = np.zeros(3)
53 
54     # Deal with masked arrays:
55     if hasattr(data, 'mask'):
56         not_mask = np.logical_not(data.mask)
57         if mask is None:
58             mask = not_mask
59         else:
60             mask *= not_mask
61         data = np.asarray(data)
62 
63     # Get rid of potential memmapping
64     data = as_ndarray(data)
65     my_map = data.copy()
66     if mask is not None:
67         # check against empty mask
68         if mask.sum() == 0.:
69             warnings.warn(
70                 "Provided mask is empty. Returning center of mass instead.")
71             cut_coords = ndimage.center_of_mass(np.abs(my_map)) + offset
72             x_map, y_map, z_map = cut_coords
73             return np.asarray(coord_transform(x_map, y_map, z_map,
74                                               img.get_affine())).tolist()
75         slice_x, slice_y, slice_z = ndimage.find_objects(mask)[0]
76         my_map = my_map[slice_x, slice_y, slice_z]
77         mask = mask[slice_x, slice_y, slice_z]
78         my_map *= mask
79         offset += [slice_x.start, slice_y.start, slice_z.start]
80 
81     # Testing min and max is faster than np.all(my_map == 0)
82     if (my_map.max() == 0) and (my_map.min() == 0):
83         return .5 * np.array(data.shape)
84     if activation_threshold is None:
85         activation_threshold = fast_abs_percentile(my_map[my_map != 0].ravel(),
86                                                    80)
87     mask = np.abs(my_map) > activation_threshold - 1.e-15
88     # mask may be zero everywhere in rare cases
89     if mask.max() == 0:
90         return .5 * np.array(data.shape)
91     mask = largest_connected_component(mask)
92     slice_x, slice_y, slice_z = ndimage.find_objects(mask)[0]
93     my_map = my_map[slice_x, slice_y, slice_z]
94     mask = mask[slice_x, slice_y, slice_z]
95     my_map *= mask
96     offset += [slice_x.start, slice_y.start, slice_z.start]
97 
98     # For the second threshold, we use a mean, as it is much faster,
99     # althought it is less robust
100     second_threshold = np.abs(np.mean(my_map[mask]))
101     second_mask = (np.abs(my_map) > second_threshold)
102     if second_mask.sum() > 50:
103         my_map *= largest_connected_component(second_mask)
104     cut_coords = ndimage.center_of_mass(np.abs(my_map))
105     x_map, y_map, z_map = cut_coords + offset
106 
107     # Return as a list of scalars
108     return np.asarray(coord_transform(x_map, y_map, z_map,
109                                       img.get_affine())).tolist()
110 
111 
112 def _get_auto_mask_bounds(img):
113     """ Compute the bounds of the data with an automaticaly computed mask
114     """
115     data = img.get_data().copy()
116     affine = img.get_affine()
117     if hasattr(data, 'mask'):
118         # Masked array
119         mask = np.logical_not(data.mask)
120         data = np.asarray(data)
121     else:
122         # The mask will be anything that is fairly different
123         # from the values in the corners
124         edge_value = float(data[0, 0, 0] + data[0, -1, 0]
125                             + data[-1, 0, 0] + data[0, 0, -1]
126                             + data[-1, -1, 0] + data[-1, 0, -1]
127                             + data[0, -1, -1] + data[-1, -1, -1]
128                         )
129         edge_value /= 6
130         mask = np.abs(data - edge_value) > .005*data.ptp()
131     xmin, xmax, ymin, ymax, zmin, zmax = \
132             get_mask_bounds(new_img_like(img, mask, affine))
133     return (xmin, xmax), (ymin, ymax), (zmin, zmax)
134 
135 
136 def _transform_cut_coords(cut_coords, direction, affine):
137     """Transforms cut_coords back in image space
138 
139     Parameters
140     ----------
141     cut_coords: 1D array of length n_cuts
142         The coordinates to be transformed.
143 
144     direction: string, optional (default "z")
145         sectional direction; possible values are "x", "y", or "z"
146 
147     affine: 2D array of shape (4, 4)
148         The affine for the image.
149 
150     Returns
151     -------
152     cut_coords: 1D array of length n_cuts
153        The original cut_coords transformed image space.
154     """
155     # make kwargs
156     axis = 'xyz'.index(direction)
157     kwargs = {}
158     for name in 'xyz':
159         kwargs[name] = np.zeros(len(cut_coords))
160     kwargs[direction] = cut_coords
161     kwargs['affine'] = affine
162 
163     # We need atleast_1d to make sure that when n_cuts is 1 we do
164     # get an iterable
165     cut_coords = coord_transform(**kwargs)[axis]
166     return np.atleast_1d(cut_coords)
167 
168 
169 def find_cut_slices(img, direction='z', n_cuts=7, spacing='auto'):
170     """ Find 'good' cross-section slicing positions along a given axis.
171 
172     Parameters
173     ----------
174     img: 3D Nifti1Image
175         the brain map
176     direction: string, optional (default "z")
177         sectional direction; possible values are "x", "y", or "z"
178     n_cuts: int, optional (default 7)
179         number of cuts in the plot
180     spacing: 'auto' or int, optional (default 'auto')
181         minimum spacing between cuts (in voxels, not milimeters)
182         if 'auto', the spacing is .5 / n_cuts * img_length
183 
184     Returns
185     -------
186     cut_coords: 1D array of length n_cuts
187         the computed cut_coords
188 
189     Notes
190     -----
191     This code works by iteratively locating peak activations that are
192     separated by a distance of at least 'spacing'. If n_cuts is very
193     large and all the activated regions are covered, cuts with a spacing
194     less than 'spacing' will be returned.
195     """
196 
197     # misc
198     if not direction in 'xyz':
199         raise ValueError(
200             "'direction' must be one of 'x', 'y', or 'z'. Got '%s'" % (
201                 direction))
202     axis = 'xyz'.index(direction)
203     affine = img.get_affine()
204     orig_data = np.abs(img.get_data())
205     this_shape = orig_data.shape[axis]
206 
207     if not isinstance(n_cuts, numbers.Number):
208         raise ValueError("The number of cuts (n_cuts) must be an integer "
209                          "greater than or equal to 1. "
210                          "You provided a value of n_cuts=%s. " % n_cuts)
211 
212     # BF issue #575: Return all the slices along and axis if this axis
213     # is the display mode and there are at least as many requested
214     # n_slices as there are slices.
215     if n_cuts > this_shape:
216         warnings.warn('Too many cuts requested for the data: '
217                       'n_cuts=%i, data size=%i' % (n_cuts, this_shape))
218         return _transform_cut_coords(np.arange(this_shape), direction, affine)
219 
220     data = orig_data.copy()
221     if data.dtype.kind == 'i':
222         data = data.astype(np.float)
223 
224     data = _smooth_array(data, affine, fwhm='fast')
225 
226     # to control floating point error problems
227     # during given input value "n_cuts"
228     epsilon = np.finfo(np.float32).eps
229     difference = abs(round(n_cuts) - n_cuts)
230     if round(n_cuts) < 1. or difference > epsilon:
231         message = ("Image has %d slices in direction %s. "
232                    "Therefore, the number of cuts must be between 1 and %d. "
233                    "You provided n_cuts=%s " % (
234                        this_shape, direction, this_shape, n_cuts))
235         raise ValueError(message)
236     else:
237         n_cuts = int(round(n_cuts))
238 
239     if spacing == 'auto':
240         spacing = max(int(.5 / n_cuts * data.shape[axis]), 1)
241 
242     slices = [slice(None, None), slice(None, None), slice(None, None)]
243 
244     cut_coords = list()
245 
246     for _ in range(n_cuts):
247         # Find a peak
248         max_along_axis = np.unravel_index(np.abs(data).argmax(),
249                                           data.shape)[axis]
250 
251         # cancel out the surroundings of the peak
252         start = max(0, max_along_axis - spacing)
253         stop = max_along_axis + spacing
254         slices[axis] = slice(start, stop)
255         # We don't actually fully zero the neighborhood, to avoid ending
256         # up with fully zeros if n_cuts is too big: we can do multiple
257         # passes on the data
258         data[slices] *= 1.e-3
259 
260         cut_coords.append(max_along_axis)
261 
262     # We sometimes get duplicated cuts, so we add cuts at the beginning
263     # and the end
264     cut_coords = np.unique(cut_coords).tolist()
265     while len(cut_coords) < n_cuts:
266         # Candidates for new cuts:
267         slice_below = min(cut_coords) - 2
268         slice_above = max(cut_coords) + 2
269         candidates = [slice_above]
270         # One slice where there is the biggest gap in the existing
271         # cut_coords
272         if len(cut_coords) > 1:
273             middle_idx = np.argmax(np.diff(cut_coords))
274             slice_middle = int(.5 * (cut_coords[middle_idx]
275                                     + cut_coords[middle_idx + 1]))
276             if not slice_middle in cut_coords:
277                 candidates.append(slice_middle)
278         if slice_below >= 0:
279             # We need positive slice to avoid having negative
280             # indices, which would work, but not the way we think of them
281             candidates.append(slice_below)
282         best_weight = -10
283         for candidate in candidates:
284             if candidate >= this_shape:
285                 this_weight = 0
286             else:
287                 this_weight = np.sum(np.rollaxis(orig_data, axis)[candidate])
288             if this_weight > best_weight:
289                 best_candidate = candidate
290                 best_weight = this_weight
291 
292         cut_coords.append(best_candidate)
293         cut_coords = np.unique(cut_coords).tolist()
294 
295     cut_coords = np.array(cut_coords)
296     cut_coords.sort()
297 
298     return _transform_cut_coords(cut_coords, direction, affine)
299
```
--- 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/plotting/find_cuts.py b/nilearn/plotting/find_cuts.py
--- a/nilearn/plotting/find_cuts.py
+++ b/nilearn/plotting/find_cuts.py
@@ -15,6 +15,8 @@
 from ..image import new_img_like
 from .._utils.extmath import fast_abs_percentile
 from .._utils.numpy_conversions import as_ndarray
+from .._utils import check_niimg_3d
+from .._utils.niimg import _safe_get_data
 from ..image.resampling import get_mask_bounds, coord_transform
 from ..image.image import _smooth_array
 
@@ -46,7 +48,11 @@
         z : float
             the z world coordinate.
     """
-    data = img.get_data()
+    # if a pseudo-4D image or several images were passed (cf. #922),
+    # we reduce to a single 3D image to find the coordinates
+    img = check_niimg_3d(img)
+    data = _safe_get_data(img)
+
     # To speed up computations, we work with partial views of the array,
     # and keep track of the offset
     offset = np.zeros(3)

{"golden_diff": "diff --git a/nilearn/plotting/find_cuts.py b/nilearn/plotting/find_cuts.py\n--- a/nilearn/plotting/find_cuts.py\n+++ b/nilearn/plotting/find_cuts.py\n@@ -15,6 +15,8 @@\n from ..image import new_img_like\n from .._utils.extmath import fast_abs_percentile\n from .._utils.numpy_conversions import as_ndarray\n+from .._utils import check_niimg_3d\n+from .._utils.niimg import _safe_get_data\n from ..image.resampling import get_mask_bounds, coord_transform\n from ..image.image import _smooth_array\n \n@@ -46,7 +48,11 @@\n         z : float\n             the z world coordinate.\n     \"\"\"\n-    data = img.get_data()\n+    # if a pseudo-4D image or several images were passed (cf. #922),\n+    # we reduce to a single 3D image to find the coordinates\n+    img = check_niimg_3d(img)\n+    data = _safe_get_data(img)\n+\n     # To speed up computations, we work with partial views of the array,\n     # and keep track of the offset\n     offset = np.zeros(3)\n", "issue": "find_xyz_cut_coords crashes on some compact binary images\n...trace back is:\n\n``` python\n---------------------------------------------------------------------------\nValueError                                Traceback (most recent call last)\n<ipython-input-44-802cef648458> in <module>()\n----> 1 plotting.find_xyz_cut_coords(nib.load(roi_paths[2]))\n\n/git/nilearn/nilearn/plotting/find_cuts.py in find_xyz_cut_coords(img, mask, activation_threshold)\n     90         return .5 * np.array(data.shape)\n     91     mask = largest_connected_component(mask)\n---> 92     slice_x, slice_y, slice_z = ndimage.find_objects(mask)[0]\n     93     my_map = my_map[slice_x, slice_y, slice_z]\n     94     mask = mask[slice_x, slice_y, slice_z]\n```\n\nThe image on which this crashes is attached.\n\n[dmpfc-3.nii.zip](https://github.com/nilearn/nilearn/files/66476/dmpfc-3.nii.zip)\n\n", "before_files": [{"content": "\"\"\"\nTools to find activations and cut on maps\n\"\"\"\n\n# Author: Gael Varoquaux\n# License: BSD\n\nimport warnings\nimport numbers\nimport numpy as np\nfrom scipy import ndimage\n\n# Local imports\nfrom .._utils.ndimage import largest_connected_component\nfrom ..image import new_img_like\nfrom .._utils.extmath import fast_abs_percentile\nfrom .._utils.numpy_conversions import as_ndarray\nfrom ..image.resampling import get_mask_bounds, coord_transform\nfrom ..image.image import _smooth_array\n\n################################################################################\n# Functions for automatic choice of cuts coordinates\n################################################################################\n\n\ndef find_xyz_cut_coords(img, mask=None, activation_threshold=None):\n    \"\"\" Find the center of the largest activation connected component.\n\n        Parameters\n        -----------\n        img : 3D Nifti1Image\n            The brain map.\n        mask : 3D ndarray, boolean, optional\n            An optional brain mask.\n        activation_threshold : float, optional\n            The lower threshold to the positive activation. If None, the\n            activation threshold is computed using the 80% percentile of\n            the absolute value of the map.\n\n        Returns\n        -------\n        x : float\n            the x world coordinate.\n        y : float\n            the y world coordinate.\n        z : float\n            the z world coordinate.\n    \"\"\"\n    data = img.get_data()\n    # To speed up computations, we work with partial views of the array,\n    # and keep track of the offset\n    offset = np.zeros(3)\n\n    # Deal with masked arrays:\n    if hasattr(data, 'mask'):\n        not_mask = np.logical_not(data.mask)\n        if mask is None:\n            mask = not_mask\n        else:\n            mask *= not_mask\n        data = np.asarray(data)\n\n    # Get rid of potential memmapping\n    data = as_ndarray(data)\n    my_map = data.copy()\n    if mask is not None:\n        # check against empty mask\n        if mask.sum() == 0.:\n            warnings.warn(\n                \"Provided mask is empty. Returning center of mass instead.\")\n            cut_coords = ndimage.center_of_mass(np.abs(my_map)) + offset\n            x_map, y_map, z_map = cut_coords\n            return np.asarray(coord_transform(x_map, y_map, z_map,\n                                              img.get_affine())).tolist()\n        slice_x, slice_y, slice_z = ndimage.find_objects(mask)[0]\n        my_map = my_map[slice_x, slice_y, slice_z]\n        mask = mask[slice_x, slice_y, slice_z]\n        my_map *= mask\n        offset += [slice_x.start, slice_y.start, slice_z.start]\n\n    # Testing min and max is faster than np.all(my_map == 0)\n    if (my_map.max() == 0) and (my_map.min() == 0):\n        return .5 * np.array(data.shape)\n    if activation_threshold is None:\n        activation_threshold = fast_abs_percentile(my_map[my_map != 0].ravel(),\n                                                   80)\n    mask = np.abs(my_map) > activation_threshold - 1.e-15\n    # mask may be zero everywhere in rare cases\n    if mask.max() == 0:\n        return .5 * np.array(data.shape)\n    mask = largest_connected_component(mask)\n    slice_x, slice_y, slice_z = ndimage.find_objects(mask)[0]\n    my_map = my_map[slice_x, slice_y, slice_z]\n    mask = mask[slice_x, slice_y, slice_z]\n    my_map *= mask\n    offset += [slice_x.start, slice_y.start, slice_z.start]\n\n    # For the second threshold, we use a mean, as it is much faster,\n    # althought it is less robust\n    second_threshold = np.abs(np.mean(my_map[mask]))\n    second_mask = (np.abs(my_map) > second_threshold)\n    if second_mask.sum() > 50:\n        my_map *= largest_connected_component(second_mask)\n    cut_coords = ndimage.center_of_mass(np.abs(my_map))\n    x_map, y_map, z_map = cut_coords + offset\n\n    # Return as a list of scalars\n    return np.asarray(coord_transform(x_map, y_map, z_map,\n                                      img.get_affine())).tolist()\n\n\ndef _get_auto_mask_bounds(img):\n    \"\"\" Compute the bounds of the data with an automaticaly computed mask\n    \"\"\"\n    data = img.get_data().copy()\n    affine = img.get_affine()\n    if hasattr(data, 'mask'):\n        # Masked array\n        mask = np.logical_not(data.mask)\n        data = np.asarray(data)\n    else:\n        # The mask will be anything that is fairly different\n        # from the values in the corners\n        edge_value = float(data[0, 0, 0] + data[0, -1, 0]\n                            + data[-1, 0, 0] + data[0, 0, -1]\n                            + data[-1, -1, 0] + data[-1, 0, -1]\n                            + data[0, -1, -1] + data[-1, -1, -1]\n                        )\n        edge_value /= 6\n        mask = np.abs(data - edge_value) > .005*data.ptp()\n    xmin, xmax, ymin, ymax, zmin, zmax = \\\n            get_mask_bounds(new_img_like(img, mask, affine))\n    return (xmin, xmax), (ymin, ymax), (zmin, zmax)\n\n\ndef _transform_cut_coords(cut_coords, direction, affine):\n    \"\"\"Transforms cut_coords back in image space\n\n    Parameters\n    ----------\n    cut_coords: 1D array of length n_cuts\n        The coordinates to be transformed.\n\n    direction: string, optional (default \"z\")\n        sectional direction; possible values are \"x\", \"y\", or \"z\"\n\n    affine: 2D array of shape (4, 4)\n        The affine for the image.\n\n    Returns\n    -------\n    cut_coords: 1D array of length n_cuts\n       The original cut_coords transformed image space.\n    \"\"\"\n    # make kwargs\n    axis = 'xyz'.index(direction)\n    kwargs = {}\n    for name in 'xyz':\n        kwargs[name] = np.zeros(len(cut_coords))\n    kwargs[direction] = cut_coords\n    kwargs['affine'] = affine\n\n    # We need atleast_1d to make sure that when n_cuts is 1 we do\n    # get an iterable\n    cut_coords = coord_transform(**kwargs)[axis]\n    return np.atleast_1d(cut_coords)\n\n\ndef find_cut_slices(img, direction='z', n_cuts=7, spacing='auto'):\n    \"\"\" Find 'good' cross-section slicing positions along a given axis.\n\n    Parameters\n    ----------\n    img: 3D Nifti1Image\n        the brain map\n    direction: string, optional (default \"z\")\n        sectional direction; possible values are \"x\", \"y\", or \"z\"\n    n_cuts: int, optional (default 7)\n        number of cuts in the plot\n    spacing: 'auto' or int, optional (default 'auto')\n        minimum spacing between cuts (in voxels, not milimeters)\n        if 'auto', the spacing is .5 / n_cuts * img_length\n\n    Returns\n    -------\n    cut_coords: 1D array of length n_cuts\n        the computed cut_coords\n\n    Notes\n    -----\n    This code works by iteratively locating peak activations that are\n    separated by a distance of at least 'spacing'. If n_cuts is very\n    large and all the activated regions are covered, cuts with a spacing\n    less than 'spacing' will be returned.\n    \"\"\"\n\n    # misc\n    if not direction in 'xyz':\n        raise ValueError(\n            \"'direction' must be one of 'x', 'y', or 'z'. Got '%s'\" % (\n                direction))\n    axis = 'xyz'.index(direction)\n    affine = img.get_affine()\n    orig_data = np.abs(img.get_data())\n    this_shape = orig_data.shape[axis]\n\n    if not isinstance(n_cuts, numbers.Number):\n        raise ValueError(\"The number of cuts (n_cuts) must be an integer \"\n                         \"greater than or equal to 1. \"\n                         \"You provided a value of n_cuts=%s. \" % n_cuts)\n\n    # BF issue #575: Return all the slices along and axis if this axis\n    # is the display mode and there are at least as many requested\n    # n_slices as there are slices.\n    if n_cuts > this_shape:\n        warnings.warn('Too many cuts requested for the data: '\n                      'n_cuts=%i, data size=%i' % (n_cuts, this_shape))\n        return _transform_cut_coords(np.arange(this_shape), direction, affine)\n\n    data = orig_data.copy()\n    if data.dtype.kind == 'i':\n        data = data.astype(np.float)\n\n    data = _smooth_array(data, affine, fwhm='fast')\n\n    # to control floating point error problems\n    # during given input value \"n_cuts\"\n    epsilon = np.finfo(np.float32).eps\n    difference = abs(round(n_cuts) - n_cuts)\n    if round(n_cuts) < 1. or difference > epsilon:\n        message = (\"Image has %d slices in direction %s. \"\n                   \"Therefore, the number of cuts must be between 1 and %d. \"\n                   \"You provided n_cuts=%s \" % (\n                       this_shape, direction, this_shape, n_cuts))\n        raise ValueError(message)\n    else:\n        n_cuts = int(round(n_cuts))\n\n    if spacing == 'auto':\n        spacing = max(int(.5 / n_cuts * data.shape[axis]), 1)\n\n    slices = [slice(None, None), slice(None, None), slice(None, None)]\n\n    cut_coords = list()\n\n    for _ in range(n_cuts):\n        # Find a peak\n        max_along_axis = np.unravel_index(np.abs(data).argmax(),\n                                          data.shape)[axis]\n\n        # cancel out the surroundings of the peak\n        start = max(0, max_along_axis - spacing)\n        stop = max_along_axis + spacing\n        slices[axis] = slice(start, stop)\n        # We don't actually fully zero the neighborhood, to avoid ending\n        # up with fully zeros if n_cuts is too big: we can do multiple\n        # passes on the data\n        data[slices] *= 1.e-3\n\n        cut_coords.append(max_along_axis)\n\n    # We sometimes get duplicated cuts, so we add cuts at the beginning\n    # and the end\n    cut_coords = np.unique(cut_coords).tolist()\n    while len(cut_coords) < n_cuts:\n        # Candidates for new cuts:\n        slice_below = min(cut_coords) - 2\n        slice_above = max(cut_coords) + 2\n        candidates = [slice_above]\n        # One slice where there is the biggest gap in the existing\n        # cut_coords\n        if len(cut_coords) > 1:\n            middle_idx = np.argmax(np.diff(cut_coords))\n            slice_middle = int(.5 * (cut_coords[middle_idx]\n                                    + cut_coords[middle_idx + 1]))\n            if not slice_middle in cut_coords:\n                candidates.append(slice_middle)\n        if slice_below >= 0:\n            # We need positive slice to avoid having negative\n            # indices, which would work, but not the way we think of them\n            candidates.append(slice_below)\n        best_weight = -10\n        for candidate in candidates:\n            if candidate >= this_shape:\n                this_weight = 0\n            else:\n                this_weight = np.sum(np.rollaxis(orig_data, axis)[candidate])\n            if this_weight > best_weight:\n                best_candidate = candidate\n                best_weight = this_weight\n\n        cut_coords.append(best_candidate)\n        cut_coords = np.unique(cut_coords).tolist()\n\n    cut_coords = np.array(cut_coords)\n    cut_coords.sort()\n\n    return _transform_cut_coords(cut_coords, direction, affine)\n", "path": "nilearn/plotting/find_cuts.py"}], "after_files": [{"content": "\"\"\"\nTools to find activations and cut on maps\n\"\"\"\n\n# Author: Gael Varoquaux\n# License: BSD\n\nimport warnings\nimport numbers\nimport numpy as np\nfrom scipy import ndimage\n\n# Local imports\nfrom .._utils.ndimage import largest_connected_component\nfrom ..image import new_img_like\nfrom .._utils.extmath import fast_abs_percentile\nfrom .._utils.numpy_conversions import as_ndarray\nfrom .._utils import check_niimg_3d\nfrom .._utils.niimg import _safe_get_data\nfrom ..image.resampling import get_mask_bounds, coord_transform\nfrom ..image.image import _smooth_array\n\n################################################################################\n# Functions for automatic choice of cuts coordinates\n################################################################################\n\n\ndef find_xyz_cut_coords(img, mask=None, activation_threshold=None):\n    \"\"\" Find the center of the largest activation connected component.\n\n        Parameters\n        -----------\n        img : 3D Nifti1Image\n            The brain map.\n        mask : 3D ndarray, boolean, optional\n            An optional brain mask.\n        activation_threshold : float, optional\n            The lower threshold to the positive activation. If None, the\n            activation threshold is computed using the 80% percentile of\n            the absolute value of the map.\n\n        Returns\n        -------\n        x : float\n            the x world coordinate.\n        y : float\n            the y world coordinate.\n        z : float\n            the z world coordinate.\n    \"\"\"\n    # if a pseudo-4D image or several images were passed (cf. #922),\n    # we reduce to a single 3D image to find the coordinates\n    img = check_niimg_3d(img)\n    data = _safe_get_data(img)\n\n    # To speed up computations, we work with partial views of the array,\n    # and keep track of the offset\n    offset = np.zeros(3)\n\n    # Deal with masked arrays:\n    if hasattr(data, 'mask'):\n        not_mask = np.logical_not(data.mask)\n        if mask is None:\n            mask = not_mask\n        else:\n            mask *= not_mask\n        data = np.asarray(data)\n\n    # Get rid of potential memmapping\n    data = as_ndarray(data)\n    my_map = data.copy()\n    if mask is not None:\n        # check against empty mask\n        if mask.sum() == 0.:\n            warnings.warn(\n                \"Provided mask is empty. Returning center of mass instead.\")\n            cut_coords = ndimage.center_of_mass(np.abs(my_map)) + offset\n            x_map, y_map, z_map = cut_coords\n            return np.asarray(coord_transform(x_map, y_map, z_map,\n                                              img.get_affine())).tolist()\n        slice_x, slice_y, slice_z = ndimage.find_objects(mask)[0]\n        my_map = my_map[slice_x, slice_y, slice_z]\n        mask = mask[slice_x, slice_y, slice_z]\n        my_map *= mask\n        offset += [slice_x.start, slice_y.start, slice_z.start]\n\n    # Testing min and max is faster than np.all(my_map == 0)\n    if (my_map.max() == 0) and (my_map.min() == 0):\n        return .5 * np.array(data.shape)\n    if activation_threshold is None:\n        activation_threshold = fast_abs_percentile(my_map[my_map != 0].ravel(),\n                                                   80)\n    mask = np.abs(my_map) > activation_threshold - 1.e-15\n    # mask may be zero everywhere in rare cases\n    if mask.max() == 0:\n        return .5 * np.array(data.shape)\n    mask = largest_connected_component(mask)\n    slice_x, slice_y, slice_z = ndimage.find_objects(mask)[0]\n    my_map = my_map[slice_x, slice_y, slice_z]\n    mask = mask[slice_x, slice_y, slice_z]\n    my_map *= mask\n    offset += [slice_x.start, slice_y.start, slice_z.start]\n\n    # For the second threshold, we use a mean, as it is much faster,\n    # althought it is less robust\n    second_threshold = np.abs(np.mean(my_map[mask]))\n    second_mask = (np.abs(my_map) > second_threshold)\n    if second_mask.sum() > 50:\n        my_map *= largest_connected_component(second_mask)\n    cut_coords = ndimage.center_of_mass(np.abs(my_map))\n    x_map, y_map, z_map = cut_coords + offset\n\n    # Return as a list of scalars\n    return np.asarray(coord_transform(x_map, y_map, z_map,\n                                      img.get_affine())).tolist()\n\n\ndef _get_auto_mask_bounds(img):\n    \"\"\" Compute the bounds of the data with an automaticaly computed mask\n    \"\"\"\n    data = img.get_data().copy()\n    affine = img.get_affine()\n    if hasattr(data, 'mask'):\n        # Masked array\n        mask = np.logical_not(data.mask)\n        data = np.asarray(data)\n    else:\n        # The mask will be anything that is fairly different\n        # from the values in the corners\n        edge_value = float(data[0, 0, 0] + data[0, -1, 0]\n                            + data[-1, 0, 0] + data[0, 0, -1]\n                            + data[-1, -1, 0] + data[-1, 0, -1]\n                            + data[0, -1, -1] + data[-1, -1, -1]\n                        )\n        edge_value /= 6\n        mask = np.abs(data - edge_value) > .005*data.ptp()\n    xmin, xmax, ymin, ymax, zmin, zmax = \\\n            get_mask_bounds(new_img_like(img, mask, affine))\n    return (xmin, xmax), (ymin, ymax), (zmin, zmax)\n\n\ndef _transform_cut_coords(cut_coords, direction, affine):\n    \"\"\"Transforms cut_coords back in image space\n\n    Parameters\n    ----------\n    cut_coords: 1D array of length n_cuts\n        The coordinates to be transformed.\n\n    direction: string, optional (default \"z\")\n        sectional direction; possible values are \"x\", \"y\", or \"z\"\n\n    affine: 2D array of shape (4, 4)\n        The affine for the image.\n\n    Returns\n    -------\n    cut_coords: 1D array of length n_cuts\n       The original cut_coords transformed image space.\n    \"\"\"\n    # make kwargs\n    axis = 'xyz'.index(direction)\n    kwargs = {}\n    for name in 'xyz':\n        kwargs[name] = np.zeros(len(cut_coords))\n    kwargs[direction] = cut_coords\n    kwargs['affine'] = affine\n\n    # We need atleast_1d to make sure that when n_cuts is 1 we do\n    # get an iterable\n    cut_coords = coord_transform(**kwargs)[axis]\n    return np.atleast_1d(cut_coords)\n\n\ndef find_cut_slices(img, direction='z', n_cuts=7, spacing='auto'):\n    \"\"\" Find 'good' cross-section slicing positions along a given axis.\n\n    Parameters\n    ----------\n    img: 3D Nifti1Image\n        the brain map\n    direction: string, optional (default \"z\")\n        sectional direction; possible values are \"x\", \"y\", or \"z\"\n    n_cuts: int, optional (default 7)\n        number of cuts in the plot\n    spacing: 'auto' or int, optional (default 'auto')\n        minimum spacing between cuts (in voxels, not milimeters)\n        if 'auto', the spacing is .5 / n_cuts * img_length\n\n    Returns\n    -------\n    cut_coords: 1D array of length n_cuts\n        the computed cut_coords\n\n    Notes\n    -----\n    This code works by iteratively locating peak activations that are\n    separated by a distance of at least 'spacing'. If n_cuts is very\n    large and all the activated regions are covered, cuts with a spacing\n    less than 'spacing' will be returned.\n    \"\"\"\n\n    # misc\n    if not direction in 'xyz':\n        raise ValueError(\n            \"'direction' must be one of 'x', 'y', or 'z'. Got '%s'\" % (\n                direction))\n    axis = 'xyz'.index(direction)\n    affine = img.get_affine()\n    orig_data = np.abs(img.get_data())\n    this_shape = orig_data.shape[axis]\n\n    if not isinstance(n_cuts, numbers.Number):\n        raise ValueError(\"The number of cuts (n_cuts) must be an integer \"\n                         \"greater than or equal to 1. \"\n                         \"You provided a value of n_cuts=%s. \" % n_cuts)\n\n    # BF issue #575: Return all the slices along and axis if this axis\n    # is the display mode and there are at least as many requested\n    # n_slices as there are slices.\n    if n_cuts > this_shape:\n        warnings.warn('Too many cuts requested for the data: '\n                      'n_cuts=%i, data size=%i' % (n_cuts, this_shape))\n        return _transform_cut_coords(np.arange(this_shape), direction, affine)\n\n    data = orig_data.copy()\n    if data.dtype.kind == 'i':\n        data = data.astype(np.float)\n\n    data = _smooth_array(data, affine, fwhm='fast')\n\n    # to control floating point error problems\n    # during given input value \"n_cuts\"\n    epsilon = np.finfo(np.float32).eps\n    difference = abs(round(n_cuts) - n_cuts)\n    if round(n_cuts) < 1. or difference > epsilon:\n        message = (\"Image has %d slices in direction %s. \"\n                   \"Therefore, the number of cuts must be between 1 and %d. \"\n                   \"You provided n_cuts=%s \" % (\n                       this_shape, direction, this_shape, n_cuts))\n        raise ValueError(message)\n    else:\n        n_cuts = int(round(n_cuts))\n\n    if spacing == 'auto':\n        spacing = max(int(.5 / n_cuts * data.shape[axis]), 1)\n\n    slices = [slice(None, None), slice(None, None), slice(None, None)]\n\n    cut_coords = list()\n\n    for _ in range(n_cuts):\n        # Find a peak\n        max_along_axis = np.unravel_index(np.abs(data).argmax(),\n                                          data.shape)[axis]\n\n        # cancel out the surroundings of the peak\n        start = max(0, max_along_axis - spacing)\n        stop = max_along_axis + spacing\n        slices[axis] = slice(start, stop)\n        # We don't actually fully zero the neighborhood, to avoid ending\n        # up with fully zeros if n_cuts is too big: we can do multiple\n        # passes on the data\n        data[slices] *= 1.e-3\n\n        cut_coords.append(max_along_axis)\n\n    # We sometimes get duplicated cuts, so we add cuts at the beginning\n    # and the end\n    cut_coords = np.unique(cut_coords).tolist()\n    while len(cut_coords) < n_cuts:\n        # Candidates for new cuts:\n        slice_below = min(cut_coords) - 2\n        slice_above = max(cut_coords) + 2\n        candidates = [slice_above]\n        # One slice where there is the biggest gap in the existing\n        # cut_coords\n        if len(cut_coords) > 1:\n            middle_idx = np.argmax(np.diff(cut_coords))\n            slice_middle = int(.5 * (cut_coords[middle_idx]\n                                    + cut_coords[middle_idx + 1]))\n            if not slice_middle in cut_coords:\n                candidates.append(slice_middle)\n        if slice_below >= 0:\n            # We need positive slice to avoid having negative\n            # indices, which would work, but not the way we think of them\n            candidates.append(slice_below)\n        best_weight = -10\n        for candidate in candidates:\n            if candidate >= this_shape:\n                this_weight = 0\n            else:\n                this_weight = np.sum(np.rollaxis(orig_data, axis)[candidate])\n            if this_weight > best_weight:\n                best_candidate = candidate\n                best_weight = this_weight\n\n        cut_coords.append(best_candidate)\n        cut_coords = np.unique(cut_coords).tolist()\n\n    cut_coords = np.array(cut_coords)\n    cut_coords.sort()\n\n    return _transform_cut_coords(cut_coords, direction, affine)\n", "path": "nilearn/plotting/find_cuts.py"}]}