code: def _plot_mri_contours(
      *,
      mri_fname,
      surfaces,
      src,
      orientation="coronal",
      slices=None,
      show=True,
      show_indices=False,
      show_orientation=False,
      width=512,
      slices_as_subplots=True,
  ):
      """Plot BEM contours on anatomical MRI slices.
  
      Parameters
      ----------
      slices_as_subplots : bool
          Whether to add all slices as subplots to a single figure, or to
          create a new figure for each slice. If ``False``, return NumPy
          arrays instead of Matplotlib figures.
  
      Returns
      -------
      matplotlib.figure.Figure | list of array
          The plotted slices.
      """
      import matplotlib.pyplot as plt
      from matplotlib import patheffects
  
      from ..source_space._source_space import _ensure_src
  
      # For ease of plotting, we will do everything in voxel coordinates.
      _validate_type(show_orientation, (bool, str), "show_orientation")
      if isinstance(show_orientation, str):
          _check_option(
              "show_orientation",
              show_orientation,
              ("always",),
              extra="when str",
          )
      _check_option(
          "orientation", orientation, ("coronal", "axial", "sagittal")
      )
  
      # Load the T1 data
      _, _, _, _, _, nim = _read_mri_info(
          mri_fname, units="mm", return_img=True
      )
  
      data, rasvox_mri_t = _reorient_image(nim)
      mri_rasvox_t = np.linalg.inv(rasvox_mri_t)
      axis, x, y = _mri_orientation(orientation)
  
      n_slices = data.shape[axis]
  
      # if no slices were specified, pick some equally-spaced ones automatically
      if slices is None:
          slices = np.round(
              np.linspace(start=0, stop=n_slices - 1, num=14)
          ).astype(int)
  
          # omit first and last one (not much brain visible there anyway…)
          slices = slices[1:-1]
  
      slices = np.atleast_1d(slices).copy()
      slices[slices < 0] += n_slices  # allow negative indexing
      if (
          not np.array_equal(np.sort(slices), slices)
          or slices.ndim != 1
          or slices.size < 1
          or slices[0] < 0
          or slices[-1] >= n_slices
          or slices.dtype.kind not in "iu"
      ):
          raise ValueError(
              "slices must be a sorted 1D array of int with unique "
              "elements, at least one element, and no elements "
              "greater than %d, got %s" % (n_slices - 1, slices)
          )
  
      # create of list of surfaces
      surfs = list()
      for file_name, color in surfaces:
          surf = dict()
          surf["rr"], surf["tris"] = read_surface(file_name)
          # move surface to voxel coordinate system
          surf["rr"] = apply_trans(mri_rasvox_t, surf["rr"])
          surfs.append((surf, color))
  
      sources = list()
      if src is not None:
          _ensure_src(src, extra=" or None")
          # Eventually we can relax this by allowing ``trans`` if need be
          if src[0]["coord_frame"] != FIFF.FIFFV_COORD_MRI:
              raise ValueError(
                  "Source space must be in MRI coordinates, got "
                  f'{_frame_to_str[src[0]["coord_frame"]]}'
              )
          for src_ in src:
              points = src_["rr"][src_["inuse"].astype(bool)]
              sources.append(apply_trans(mri_rasvox_t, points * 1e3))
          sources = np.concatenate(sources, axis=0)
  
      # get the figure dimensions right
      if slices_as_subplots:
          n_col = 4
          fig, axs, _, _ = _prepare_trellis(len(slices), n_col)
          fig.set_facecolor("k")
          dpi = fig.get_dpi()
          n_axes = len(axs)
      else:
          n_col = n_axes = 1
          dpi = 96
          # 2x standard MRI resolution is probably good enough for the
          # traces
          w = width / dpi
          figsize = (w, w / data.shape[x] * data.shape[y])
  
      bounds = np.concatenate(
          [[-np.inf], slices[:-1] + np.diff(slices) / 2.0, [np.inf]]
      )  # float
      slicer = [slice(None)] * 3
      ori_labels = dict(R="LR", A="PA", S="IS")
      xlabels, ylabels = ori_labels["RAS"[x]], ori_labels["RAS"[y]]
      path_effects = [
          patheffects.withStroke(
              linewidth=4, foreground="k", alpha=0.75
          )
      ]
      figs = []
      for ai, (sl, lower, upper) in enumerate(
          zip(slices, bounds[:-1], bounds[1:])
      ):
          if slices_as_subplots:
              ax = axs[ai]
          else:
              # No need for constrained layout here because we make our axes fill the
              # entire figure
              fig = _figure_agg(figsize=figsize, dpi=dpi, facecolor="k")
              ax = fig.add_axes(
                  [0, 0, 1, 1], frame_on=False, facecolor="k"
              )
  
          # adjust the orientations for good view
          slicer[axis] = sl
          dat = data[tuple(slicer)].T
  
          # First plot the anatomical data
          ax.imshow(dat, cmap=plt.cm.gray, origin="lower")
          ax.set_autoscale_on(False)
          ax.axis("off")
          ax.set_aspect("equal")  # XXX eventually could deal with zooms
  
          # and then plot the contours on top
          for surf, color in surfs:
              with warnings.catch_warnings(
                  record=True
              ):  # ignore contour warn
                  warnings.simplefilter("ignore")
                  ax.tricontour(
                      surf["rr"][:, x],
                      surf["rr"][:, y],
                      surf["tris"],
                      surf["rr"][:, axis],
                      levels=[sl],
                      colors=color,
                      linewidths=1.0,
                      zorder=1,
                  )
  
          if len(sources):
              in_slice = (sources[:, axis] >= lower) & (
                  sources[:, axis] < upper
              )
              ax.scatter(
                  sources[in_slice, x],
                  sources[in_slice, y],
                  marker=".",
                  color="#FF00FF",
                  s=1,
                  zorder=2,
              )
          if show_indices:
              ax.text(
                  dat.shape[1] // 8 + 0.5,
                  0.5,
                  str(sl),
                  color="w",
                  fontsize="x-small",
                  va="bottom",
                  ha="left",
              )
          # label the axes
          kwargs = dict(
              color="#66CCEE",
              fontsize="medium",
              path_effects=path_effects,
              family="monospace",
              clip_on=False,
              zorder=5,
              weight="bold",
          )
          always = show_orientation == "always"
          if show_orientation:
              if ai % n_col == 0 or always:  # left
                  ax.text(
                      0,
                      dat.shape[0] / 2.0,
                      xlabels[0],
                      va="center",
                      ha="left",
                      **kwargs,
                  )
              if (
                  ai % n_col == n_col - 1 or ai == n_axes - 1 or always
              ):  # right
                  ax.text(
                      dat.shape[1] - 1,
                      dat.shape[0] / 2.0,
                      xlabels[1],
                      va="center",
                      ha="right",
                      **kwargs,
                  )
              if ai >= n_axes - n_col or always:  # bottom
                  ax.text(
                      dat.shape[1] / 2.0,
                      0,
                      ylabels[0],
                      ha="center",
                      va="bottom",
                      **kwargs,
                  )
              if ai < n_col or n_col == 1 or always:  # top
                  ax.text(
                      dat.shape[1] / 2.0,
                      dat.shape[0] - 1,
                      ylabels[1],
                      ha="center",
                      va="top",
                      **kwargs,
                  )
  
          if not slices_as_subplots:
              # convert to NumPy array
              with io.BytesIO() as buff:
                  fig.savefig(
                      buff,
                      format="raw",
                      bbox_inches="tight",
                      pad_inches=0,
                      dpi=dpi,
                  )
                  w_, h_ = fig.canvas.get_width_height()
                  plt.close(fig)
                  buff.seek(0)
                  fig_array = np.frombuffer(
                      buff.getvalue(), dtype=np.uint8
                  )
  
              fig = fig_array.reshape((int(h_), int(w_), -1))
              figs.append(fig)
  
      if slices_as_subplots:
          plt_show(show, fig=fig)
          return fig
      else:
          return figs
  
imports: import copy
  import io
  import os
  import os.path
  import warnings
  from collections import defaultdict
  from glob import glob
  from itertools import cycle
  from pathlib import Path
  import numpy
  from scipy.signal import filtfilt, freqz, group_delay, lfilter, sosfilt, sosfiltfilt
  from .._fiff.constants import FIFF
  from .._fiff.pick import _DATA_CH_TYPES_SPLIT, _picks_by_type, pick_channels, pick_info, pick_types
  from .._fiff.proj import make_projector
  from .._freesurfer import _check_mri, _mri_orientation, _read_mri_info, _reorient_image
  from ..defaults import DEFAULTS
  from ..filter import estimate_ringing_samples
  from ..fixes import _safe_svd
  from ..rank import compute_rank
  from ..surface import read_surface
  from ..transforms import _frame_to_str, apply_trans
  from ..utils import _check_option, _mask_to_onsets_offsets, _on_missing, _pl, fill_doc, get_subjects_dir, logger, verbose, warn
  from .utils import _figure_agg, _get_color_list, _prepare_trellis, _validate_type, plt_show
label_desc: Use image node_id 4 for calling a global function with example usage: _plot_mri_contours() and returns: fig, figs