Repository: mne-tools/mne-python

Issue:
add function to plot heatmap of eye-tracking gaze position
### Describe the new feature or enhancement

I think plotting eyegaze data as a heatmap can be more intuitive than a timeseries plot. Having a function in MNE for this could be nice. Maybe we could add an `eyetracking` module to `mne.viz` for this?

### Describe your proposed implementation

Ultimately, It would be great if this function could be called from either a `Raw` or `Epochs` object, (and for `Raw`, there could be  `start` and `duration` parameters), but as a start we could just implement the function to take an `Epochs` object.

I have some scratch code for this, but it could be improved before merging into mne. For example maybe we could implement some gaussian smoothing to make the heatmap traces look nicer.

Here's a really simple example from a short saccadic pursuit tasK:

![heatmap](https://github.com/mne-tools/mne-python/assets/52462026/fed82496-6594-414c-b30c-5c0dad5b0ee6)


### Describe possible alternatives

Open to ideas!

### Additional context

Adding this function is one of the tasks I listed in my GSoC proposal.

Discussion:
larsoner: Looks like a good start, let me know when I should really look @scott-huberty !
scott-huberty: > Looks like a good start, let me know when I should really look @scott-huberty !

Alright! It's not ready yet and I'm a bit stuck on a few points, can I meet with you @britta-wstnr or one of the other devs tomorrow or next week?
scott-huberty: > Looks like a good start, let me know when I should really look @scott-huberty !

Alright! It's not ready yet and I'm a bit stuck on a few points, can I meet with you @britta-wstnr or one of the other devs tomorrow or next week?
larsoner: I can't meet tomorrow morning my time, which would be my only reasonable overlap with @britta-wstnr . So feel free to meet with her in my morning if you can and if not ping me on Discord and we can meet afternoon EDT
britta-wstnr: @scott-huberty we can meet during dev meeting today?
scott-huberty: > @scott-huberty we can meet during dev meeting today?

Sure!
scott-huberty: Sorry for the delay on this! I think this is finally ready for review. @larsoner / @drammock I assume I'll mainly be corresponding with you next week for the code sprint. Could we carve out some time next week to look at this?

Also CC @britta-wstnr as my GSOC mentor and @mscheltienne .

### Tutorials 

* [Here](https://output.circle-artifacts.com/output/job/ce8c15b8-f4c6-4afe-b8e9-697a30e7502d/artifacts/0/html/auto_tutorials/visualization/30_eyetracking_heatmap.html) is a new visualization tutorial demonstrating the new functionality.

* I've also added a section that uses the funciton to the [current eye-tracking tutorial](https://output.circle-artifacts.com/output/job/ce8c15b8-f4c6-4afe-b8e9-697a30e7502d/artifacts/0/html/auto_tutorials/preprocessing/90_eyetracking_data.html).

### Tests
Per the [failing](https://dev.azure.com/mne-tools/mne-python/_build/results?buildId=27011&view=ms.vss-test-web.build-test-results-tab&runId=117982&resultId=100060&paneView=debug) tests: I had a quick look and there it seems related to a (new?) deprecation warning related to Numba?


larsoner: For the test can you add `ignore:Code using Numba extension API.*:` to mne/conftest.py ? Agreed it's unrelated to this PR but should get things green. I can't look at the code today but agree we should next week unless everyone else is happy today, in which case feel free to merge without me looking
scott-huberty: > For the test can you add `ignore:Code using Numba extension API.*:` to mne/conftest.py ? 

Ok will do!



> I can't look at the code today but agree we should next week unless everyone else is happy today, in which case feel free to merge without me looking

No rush on my side
britta-wstnr: Back at work after being sick - will have a thorough look soon! Tutorial page looks really cool on first glimpse, @scott-huberty !
scott-huberty: @larsoner FYI I added a new entry into the `doc_dict` for `cmap` docstring (it's more general than the doct_dict entry for `cmap_topomap`, which specifically uses `"Reds"` as a default). I think this could be re-used elsewhere, for example in [plot_csd](https://github.com/mne-tools/mne-python/blob/8051f6ddf68a797ec0c9e1801aea3c73808bba03/mne/viz/misc.py#L1441-L1444). Just let me know if you see a problem with this.


@britta-wstnr / @mscheltienne, can I ask you to also review when you have the time?


scott-huberty: @drammock since you've reviewed the tutorial, a question for you on tutorial organization; The heatmap tutorial I've added currently sits in the [Tutorials](https://mne.tools/stable/auto_tutorials/index.html) section. It could probably instead live in the (and maybe would be a better fit for?) [How-to examples](https://mne.tools/stable/auto_examples/index.html) section (and maybe it would be a better fit there?). Do you have any opinion on where this tutorial should live?
drammock: > The heatmap tutorial I've added currently sits in the [Tutorials](https://mne.tools/stable/auto_tutorials/index.html) section. It could probably instead live in the [How-to examples](https://mne.tools/stable/auto_examples/index.html) section (and maybe it would be a better fit there?). Do you have any opinion on where this tutorial should live?

Good question.  It's how-to like in that it only really demonstrates how to do one thing. It's tutorial like in that the code and explanations are interleaved and somewhat verbose and there are demarcated sections. If this were in examples, it would be better IMO to omit details about the experimental paradigm, and probably just show 2 figures each with only 1 subplot: one full heatmap and one heatmap overlayed on an image. That would also let you get rid of the custom image plotting function, which would shorten the recipe. Basically for the recipes I advise to distill as much as possible down to the core commands necessary to complete the task; users can extrapolate to their use case from there.

If you're willing to do the distillations and move to the examples folder I think that would make sense. If you're sick of it LMK and I can probably find time to do it this week.
scott-huberty: > If you're willing to do the distillations and move to the examples folder I think that would make sense. If you're sick of it LMK and I can probably find time to do it this week.

No problem I can do that this week. I agree It will probably be a net-positive to simplify that tutorial so users can quickly get what they need.
scott-huberty: Hey @drammock , another question. In `mne.viz.eyetracking.heatmap` I import `from scipy.ndimage import gaussian_filter`  at the top of the file. should I instead be using `lazy_loader` to do:

```
import lazy
sp = lazy.load('scipy')
...
sp.ndimage.gaussian_filter
```

I'm still a little confused on how to use `lazy_loader`..
drammock: > should I instead be using `lazy_loader`

Nope. It's possible to do it that way, but not necessary. For us, lazy loader only ever gets used in the init files 
drammock: @larsoner I'm happy with this one but you have requested changes here.

Patch:
diff --git a/doc/api/visualization.rst b/doc/api/visualization.rst
index 03a2d22a5..359f13731 100644
--- a/doc/api/visualization.rst
+++ b/doc/api/visualization.rst
@@ -88,6 +88,21 @@ Visualization
    get_browser_backend
    use_browser_backend
 
+Eyetracking
+-----------
+
+.. currentmodule:: mne.viz.eyetracking
+
+:py:mod:`mne.viz.eyetracking`:
+
+.. automodule:: mne.viz.eyetracking
+   :no-members:
+   :no-inherited-members:
+.. autosummary::
+   :toctree: generated/
+
+   plot_gaze
+
 UI Events
 ---------
 
diff --git a/doc/changes/devel.rst b/doc/changes/devel.rst
index 912af2555..6474a846f 100644
--- a/doc/changes/devel.rst
+++ b/doc/changes/devel.rst
@@ -38,7 +38,7 @@ Enhancements
 - Add support for writing forward solutions to HDF5 and convenience function :meth:`mne.Forward.save` (:gh:`12036` by `Eric Larson`_)
 - Refactored internals of :func:`mne.read_annotations` (:gh:`11964` by `Paul Roujansky`_)
 - Add support for drawing MEG sensors in :ref:`mne coreg` (:gh:`12098` by `Eric Larson`_)
-- By default MNE-Python creates matplotlib figures with ``layout='constrained'`` rather than the default ``layout='tight'`` (:gh:`12050` by `Mathieu Scheltienne`_ and `Eric Larson`_)
+- By default MNE-Python creates matplotlib figures with ``layout='constrained'`` rather than the default ``layout='tight'`` (:gh:`12050`, :gh:`12103` by `Mathieu Scheltienne`_ and `Eric Larson`_)
 - Enhance :func:`~mne.viz.plot_evoked_field` with a GUI that has controls for time, colormap, and contour lines (:gh:`11942` by `Marijn van Vliet`_)
 - Add :class:`mne.viz.ui_events.UIEvent` linking for interactive colorbars, allowing users to link figures and change the colormap and limits interactively. This supports :func:`~mne.viz.plot_evoked_topomap`, :func:`~mne.viz.plot_ica_components`, :func:`~mne.viz.plot_tfr_topomap`, :func:`~mne.viz.plot_projs_topomap`, :meth:`~mne.Evoked.plot_image`, and :meth:`~mne.Epochs.plot_image` (:gh:`12057` by `Santeri Ruuskanen`_)
 
diff --git a/doc/documentation/datasets.rst b/doc/documentation/datasets.rst
index 3946ef64b..348da772e 100644
--- a/doc/documentation/datasets.rst
+++ b/doc/documentation/datasets.rst
@@ -481,16 +481,34 @@ EYELINK
 =======
 :func:`mne.datasets.eyelink.data_path`
 
-A small example dataset from a pupillary light reflex experiment. Both EEG (EGI) and
-eye-tracking (SR Research EyeLink; ASCII format) data were recorded and stored in
-separate files. 1 participant fixated on the screen while short light flashes appeared.
-Event onsets were recorded by a photodiode attached to the screen and were
-sent to both the EEG and eye-tracking systems.
+Two small example datasets of eye-tracking data from SR Research EyeLink.
+
+EEG-Eyetracking
+^^^^^^^^^^^^^^^
+:func:`mne.datasets.eyelink.data_path`. Data exists at ``/eeg-et/``.
+
+Contains both EEG (EGI) and eye-tracking (ASCII format) data recorded from a
+pupillary light reflex experiment, stored in separate files. 1 participant fixated
+on the screen while short light flashes appeared. Event onsets were recorded by a
+photodiode attached to the screen and were sent to both the EEG and eye-tracking
+systems. 
 
 .. topic:: Examples
 
     * :ref:`tut-eyetrack`
 
+Freeviewing
+^^^^^^^^^^^
+:func:`mne.datasets.eyelink.data_path`. Data exists at ``/freeviewing/``.
+
+Contains eye-tracking data (ASCII format) from 1 participant who was free-viewing a
+video of a natural scene. In some videos, the natural scene was pixelated such that
+the people in the scene were unrecognizable.
+
+.. topic:: Examples
+
+    * :ref:`tut-eyetrack-heatmap`
+
 References
 ==========
 
diff --git a/examples/visualization/eyetracking_plot_heatmap.py b/examples/visualization/eyetracking_plot_heatmap.py
new file mode 100644
index 000000000..00c9fee66
--- /dev/null
+++ b/examples/visualization/eyetracking_plot_heatmap.py
@@ -0,0 +1,89 @@
+# -*- coding: utf-8 -*-
+"""
+.. _tut-eyetrack-heatmap:
+
+=============================================
+Plotting eye-tracking heatmaps in MNE-Python
+=============================================
+
+This tutorial covers plotting eye-tracking position data as a heatmap.
+
+.. seealso::
+
+    :ref:`tut-importing-eyetracking-data`
+    :ref:`tut-eyetrack`
+
+"""
+
+# %%
+# Data loading
+# ------------
+#
+# As usual we start by importing the modules we need and loading some
+# :ref:`example data <eyelink-dataset>`: eye-tracking data recorded from SR research's
+# ``'.asc'`` file format.
+
+
+import matplotlib.pyplot as plt
+
+import mne
+from mne.viz.eyetracking import plot_gaze
+
+task_fpath = mne.datasets.eyelink.data_path() / "freeviewing"
+et_fpath = task_fpath / "sub-01_task-freeview_eyetrack.asc"
+stim_fpath = task_fpath / "stim" / "naturalistic.png"
+
+raw = mne.io.read_raw_eyelink(et_fpath)
+
+# %%
+# Process and epoch the data
+# --------------------------
+#
+# First we will interpolate missing data during blinks and epoch the data.
+
+mne.preprocessing.eyetracking.interpolate_blinks(raw, interpolate_gaze=True)
+raw.annotations.rename({"dvns": "natural"})  # more intuitive
+event_ids = {"natural": 1}
+events, event_dict = mne.events_from_annotations(raw, event_id=event_ids)
+
+epochs = mne.Epochs(
+    raw, events=events, event_id=event_dict, tmin=0, tmax=20, baseline=None
+)
+
+
+# %%
+# Plot a heatmap of the eye-tracking data
+# ---------------------------------------
+#
+# To make a heatmap of the eye-tracking data, we can use the function
+# :func:`~mne.viz.eyetracking.plot_gaze`. We will need to define the dimensions of our
+# canvas; for this file, the eye position data are reported in pixels, so we'll use the
+# screen resolution of the participant screen (1920x1080) as the width and height. We
+# can also use the sigma parameter to smooth the plot.
+
+px_width, px_height = 1920, 1080
+cmap = plt.get_cmap("viridis")
+plot_gaze(epochs["natural"], width=px_width, height=px_height, cmap=cmap, sigma=50)
+
+# %%
+# Overlaying plots with images
+# ----------------------------
+#
+# We can use matplotlib to plot gaze heatmaps on top of stimuli images. We'll
+# customize a :class:`~matplotlib.colors.Colormap` to make some values of the heatmap
+# completely transparent. We'll then use the ``vlim`` parameter to force the heatmap to
+# start at a value greater than the darkest value in our previous heatmap, which will
+# make the darkest colors of the heatmap transparent.
+
+cmap.set_under("k", alpha=0)  # make the lowest values transparent
+ax = plt.subplot()
+ax.imshow(plt.imread(stim_fpath))
+plot_gaze(
+    epochs["natural"],
+    width=px_width,
+    height=px_height,
+    vlim=(0.0003, None),
+    sigma=50,
+    cmap=cmap,
+    axes=ax,
+)
diff --git a/mne/datasets/config.py b/mne/datasets/config.py
index 76ea0934e..2cf4d50a0 100644
--- a/mne/datasets/config.py
+++ b/mne/datasets/config.py
@@ -345,9 +345,9 @@ MNE_DATASETS["fake"] = dict(
 
 # eyelink dataset
 MNE_DATASETS["eyelink"] = dict(
-    archive_name="eeg-eyetrack_data.zip",
-    hash="md5:c4fc788fe01737e08e9086c90cab642d",
-    url=("https://osf.io/63fjm/download?version=1"),
-    folder_name="eyelink-example-data",
+    archive_name="MNE-eyelink-data.zip",
+    hash="md5:68a6323ef17d655f1a659c3290ee1c3f",
+    url=("https://osf.io/xsu4g/download?version=1"),
+    folder_name="MNE-eyelink-data",
     config_key="MNE_DATASETS_EYELINK_PATH",
 )
diff --git a/mne/utils/docs.py b/mne/utils/docs.py
index 9e72685e7..40308a907 100644
--- a/mne/utils/docs.py
+++ b/mne/utils/docs.py
@@ -743,6 +743,14 @@ clim : str | dict
     ``pos_lims``, as the surface plot must show the magnitude.
 """
 
+docdict[
+    "cmap"
+] = """
+cmap : matplotlib colormap | str | None
+        The :class:`~matplotlib.colors.Colormap` to use. Defaults to ``None``, which
+        will use the matplotlib default colormap.
+"""
+
 docdict[
     "cmap_topomap"
 ] = """
diff --git a/mne/viz/_mpl_figure.py b/mne/viz/_mpl_figure.py
index 2149aa73d..b0a059c97 100644
--- a/mne/viz/_mpl_figure.py
+++ b/mne/viz/_mpl_figure.py
@@ -2331,38 +2331,16 @@ class MNEBrowseFigure(BrowserBase, MNEFigure):
 class MNELineFigure(MNEFigure):
     """Interactive figure for non-scrolling line plots."""
 
-    def __init__(self, inst, n_axes, figsize, *, layout=None, **kwargs):
-        super().__init__(figsize=figsize, inst=inst, layout=layout, **kwargs)
-
-        # AXES: default margins (inches)
-        l_margin = 0.8
-        r_margin = 0.2
-        b_margin = 0.65
-        t_margin = 0.35
-        # AXES: default margins (figure-relative coordinates)
-        left = self._inch_to_rel(l_margin)
-        right = 1 - self._inch_to_rel(r_margin)
-        bottom = self._inch_to_rel(b_margin, horiz=False)
-        top = 1 - self._inch_to_rel(t_margin, horiz=False)
-        # AXES: make subplots
-        axes = [self.add_subplot(n_axes, 1, 1)]
-        for ix in range(1, n_axes):
-            axes.append(self.add_subplot(n_axes, 1, ix + 1, sharex=axes[0]))
-        self.subplotpars.update(
-            left=left, bottom=bottom, top=top, right=right, hspace=0.4
-        )
-        # save useful things
-        self.mne.ax_list = axes
-
-    def _resize(self, event):
-        """Handle resize event."""
-        old_width, old_height = self.mne.fig_size_px
-        new_width, new_height = self._get_size_px()
-        new_margins = _calc_new_margins(
-            self, old_width, old_height, new_width, new_height
+    def __init__(self, inst, n_axes, figsize, *, layout="constrained", **kwargs):
+        super().__init__(
+            figsize=figsize,
+            inst=inst,
+            layout=layout,
+            sharex=True,
+            **kwargs,
         )
-        self.subplots_adjust(**new_margins)
-        self.mne.fig_size_px = (new_width, new_height)
+        for ix in range(n_axes):
+            self.add_subplot(n_axes, 1, ix + 1)
 
 
 def _close_all():
@@ -2426,11 +2404,10 @@ def _line_figure(inst, axes=None, picks=None, **kwargs):
             FigureClass=MNELineFigure,
             figsize=figsize,
             n_axes=n_axes,
-            layout=None,
             **kwargs,
         )
         fig.mne.fig_size_px = fig._get_size_px()  # can't do in __init__
-        axes = fig.mne.ax_list
+        axes = fig.axes
     return fig, axes
 
 
diff --git a/mne/viz/epochs.py b/mne/viz/epochs.py
index 5918d6d6a..613f3a2b6 100644
--- a/mne/viz/epochs.py
+++ b/mne/viz/epochs.py
@@ -302,9 +302,11 @@ def plot_epochs_image(
 
     # check for compatible `fig` / `axes`; instantiate figs if needed; add
     # fig(s) and axes into group_by
+    needs_colorbar = colorbar and (axes is not None or fig is not None)
     group_by = _validate_fig_and_axes(
-        fig, axes, group_by, evoked, colorbar, clear=clear
+        fig, axes, group_by, evoked, colorbar=needs_colorbar, clear=clear
     )
+    del fig, axes, needs_colorbar, clear
 
     # prepare images in advance to get consistent vmin/vmax.
     # At the same time, create a subsetted epochs object for each group
@@ -649,20 +651,26 @@ def _plot_epochs_image(
         ax["evoked"].xaxis.set_major_locator(loc)
         ax["evoked"].yaxis.set_major_locator(AutoLocator())
 
+    fig = ax_im.get_figure()
+
     # draw the colorbar
     if colorbar:
         from matplotlib.pyplot import colorbar as cbar
 
-        this_colorbar = cbar(im, cax=ax["colorbar"])
-        this_colorbar.ax.set_ylabel(unit, rotation=270, labelpad=12)
+        if "colorbar" in ax:  # axes supplied by user
+            this_colorbar = cbar(im, cax=ax["colorbar"])
+            this_colorbar.ax.set_ylabel(unit, rotation=270, labelpad=12)
+        else:  # we created them
+            this_colorbar = fig.colorbar(im, ax=ax_im)
+            this_colorbar.ax.set_title(unit)
         if cmap[1]:
             ax_im.CB = DraggableColorbar(
                 this_colorbar, im, kind="epochs_image", ch_type=unit
             )
 
     # finish
-    plt_show(show)
-    return ax_im.get_figure()
+    plt_show(show, fig=fig)
+    return fig
 
 
 def plot_drop_log(
diff --git a/mne/viz/evoked.py b/mne/viz/evoked.py
index 88340295d..6abcbcc0d 100644
--- a/mne/viz/evoked.py
+++ b/mne/viz/evoked.py
@@ -2027,7 +2027,7 @@ def plot_evoked_joint(
     if topomap_args.get("colorbar", True):
         from matplotlib import ticker
 
-        cbar = fig.colorbar(map_ax[0].images[0], ax=map_ax, cax=cbar_ax)
+        cbar = fig.colorbar(map_ax[0].images[0], ax=map_ax, cax=cbar_ax, shrink=0.8)
         cbar.ax.grid(False)  # auto-removal deprecated as of 2021/10/05
         if isinstance(contours, (list, np.ndarray)):
             cbar.set_ticks(contours)
diff --git a/mne/viz/eyetracking/__init__.py b/mne/viz/eyetracking/__init__.py
new file mode 100644
index 000000000..7de13fd89
--- /dev/null
+++ b/mne/viz/eyetracking/__init__.py
@@ -0,0 +1,5 @@
+"""Eye-tracking visualization routines."""
+#
+# License: BSD-3-Clause
+
+from .heatmap import plot_gaze
diff --git a/mne/viz/eyetracking/heatmap.py b/mne/viz/eyetracking/heatmap.py
new file mode 100644
index 000000000..d3ff4756d
--- /dev/null
+++ b/mne/viz/eyetracking/heatmap.py
@@ -0,0 +1,155 @@
+# Authors: Scott Huberty <seh33@uw.edu>
+#
+# License: BSD-3-Clause
+
+import numpy as np
+from scipy.ndimage import gaussian_filter
+
+from ...utils import _ensure_int, _validate_type, fill_doc, logger
+from ..utils import plt_show
+
+
+@fill_doc
+def plot_gaze(
+    epochs,
+    width,
+    height,
+    *,
+    sigma=25,
+    cmap=None,
+    alpha=1.0,
+    vlim=(None, None),
+    axes=None,
+    show=True,
+):
+    """Plot a heatmap of eyetracking gaze data.
+
+    Parameters
+    ----------
+    epochs : instance of Epochs
+        The :class:`~mne.Epochs` object containing eyegaze channels.
+    width : int
+        The width dimension of the plot canvas. For example, if the eyegaze data units
+        are pixels, and the participant screen resolution was 1920x1080, then the width
+        should be 1920.
+    height : int
+        The height dimension of the plot canvas. For example, if the eyegaze data units
+        are pixels, and the participant screen resolution was 1920x1080, then the height
+        should be 1080.
+    sigma : float | None
+        The amount of Gaussian smoothing applied to the heatmap data (standard
+        deviation in pixels). If ``None``, no smoothing is applied. Default is 25.
+    %(cmap)s
+    alpha : float
+        The opacity of the heatmap (default is 1).
+    %(vlim_plot_topomap)s
+    %(axes_plot_topomap)s
+    %(show)s
+
+    Returns
+    -------
+    fig : instance of Figure
+        The resulting figure object for the heatmap plot.
+
+    Notes
+    -----
+    .. versionadded:: 1.6
+    """
+    from mne import BaseEpochs
+    from mne._fiff.pick import _picks_to_idx
+
+    _validate_type(epochs, BaseEpochs, "epochs")
+    _validate_type(alpha, "numeric", "alpha")
+    _validate_type(sigma, ("numeric", None), "sigma")
+    width = _ensure_int(width, "width")
+    height = _ensure_int(height, "height")
+
+    pos_picks = _picks_to_idx(epochs.info, "eyegaze")
+    gaze_data = epochs.get_data(picks=pos_picks)
+    gaze_ch_loc = np.array([epochs.info["chs"][idx]["loc"] for idx in pos_picks])
+    x_data = gaze_data[:, np.where(gaze_ch_loc[:, 4] == -1)[0], :]
+    y_data = gaze_data[:, np.where(gaze_ch_loc[:, 4] == 1)[0], :]
+
+    if x_data.shape[1] > 1:  # binocular recording. Average across eyes
+        logger.info("Detected binocular recording. Averaging positions across eyes.")
+        x_data = np.nanmean(x_data, axis=1)  # shape (n_epochs, n_samples)
+        y_data = np.nanmean(y_data, axis=1)
+    canvas = np.vstack((x_data.flatten(), y_data.flatten()))  # shape (2, n_samples)
+
+    # Create 2D histogram
+    # Bin into image-like format
+    hist, _, _ = np.histogram2d(
+        canvas[1, :],
+        canvas[0, :],
+        bins=(height, width),
+        range=[[0, height], [0, width]],
+    )
+    # Convert density from samples to seconds
+    hist /= epochs.info["sfreq"]
+    # Smooth the heatmap
+    if sigma:
+        hist = gaussian_filter(hist, sigma=sigma)
+
+    return _plot_heatmap_array(
+        hist,
+        width=width,
+        height=height,
+        cmap=cmap,
+        alpha=alpha,
+        vmin=vlim[0],
+        vmax=vlim[1],
+        axes=axes,
+        show=show,
+    )
+
+
+def _plot_heatmap_array(
+    data,
+    width,
+    height,
+    cmap=None,
+    alpha=None,
+    vmin=None,
+    vmax=None,
+    axes=None,
+    show=True,
+):
+    """Plot a heatmap of eyetracking gaze data from a numpy array."""
+    import matplotlib.pyplot as plt
+
+    # Prepare axes
+    if axes is not None:
+        from matplotlib.axes import Axes
+
+        _validate_type(axes, Axes, "axes")
+        ax = axes
+        fig = ax.get_figure()
+    else:
+        fig, ax = plt.subplots(constrained_layout=True)
+
+    ax.set_title("Gaze heatmap")
+    ax.set_xlabel("X position")
+    ax.set_ylabel("Y position")
+
+    # Prepare the heatmap
+    alphas = 1 if alpha is None else alpha
+    vmin = np.nanmin(data) if vmin is None else vmin
+    vmax = np.nanmax(data) if vmax is None else vmax
+    extent = [0, width, height, 0]  # origin is the top left of the screen
+
+    # Plot heatmap
+    im = ax.imshow(
+        data,
+        aspect="equal",
+        cmap=cmap,
+        alpha=alphas,
+        extent=extent,
+        origin="upper",
+        vmin=vmin,
+        vmax=vmax,
+    )
+
+    # Prepare the colorbar
+    fig.colorbar(im, ax=ax, shrink=0.6, label="Dwell time (seconds)")
+    plt_show(show)
+    return fig
diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py
index a778c72dc..8629525ac 100644
--- a/mne/viz/topomap.py
+++ b/mne/viz/topomap.py
@@ -304,15 +304,13 @@ def _add_colorbar(
     im,
     cmap,
     *,
-    side="right",
     title=None,
     format=None,
-    size="5%",
     kind=None,
     ch_type=None,
 ):
     """Add a colorbar to an axis."""
-    cbar = ax.figure.colorbar(im, format=format)
+    cbar = ax.figure.colorbar(im, format=format, shrink=0.6)
     if cmap is not None and cmap[1]:
         ax.CB = DraggableColorbar(cbar, im, kind, ch_type)
     cax = cbar.ax
@@ -1712,7 +1710,6 @@ def plot_ica_components(
                     im,
                     cmap,
                     title="AU",
-                    side="right",
                     format=cbar_fmt,
                     kind="ica_comp_topomap",
                     ch_type=ch_type,
@@ -2564,7 +2561,7 @@ def _plot_topomap_multi_cbar(
     )
 
     if colorbar:
-        cbar, cax = _add_colorbar(ax, im, cmap, title=None, size="10%", format=cbar_fmt)
+        cbar, cax = _add_colorbar(ax, im, cmap, title=None, format=cbar_fmt)
         cbar.set_ticks(_vlim)
         if unit is not None:
             cbar.ax.set_ylabel(unit, fontsize=8)
@@ -3744,7 +3741,7 @@ def plot_bridged_electrodes(
     if title is not None:
         im.axes.set_title(title)
     if colorbar:
-        cax = fig.colorbar(im)
+        cax = fig.colorbar(im, shrink=0.6)
         cax.set_label(r"Electrical Distance ($\mu$$V^2$)")
     return fig
 
diff --git a/tutorials/epochs/20_visualize_epochs.py b/tutorials/epochs/20_visualize_epochs.py
index ef16211c8..4a1606afe 100644
--- a/tutorials/epochs/20_visualize_epochs.py
+++ b/tutorials/epochs/20_visualize_epochs.py
@@ -12,7 +12,6 @@ plot the sensor locations and projectors stored in `~mne.Epochs` objects.
 We'll start by importing the modules we need, loading the continuous (raw)
 sample data, and cropping it to save memory:
 """
-
 # %%
 
 import mne
diff --git a/tutorials/evoked/20_visualize_evoked.py b/tutorials/evoked/20_visualize_evoked.py
index 180f995d8..e07b1e2e6 100644
--- a/tutorials/evoked/20_visualize_evoked.py
+++ b/tutorials/evoked/20_visualize_evoked.py
@@ -10,6 +10,7 @@ This tutorial shows the different visualization methods for
 
 As usual we'll start by importing the modules we need:
 """
+
 # %%
 
 import numpy as np
diff --git a/tutorials/preprocessing/90_eyetracking_data.py b/tutorials/preprocessing/90_eyetracking_data.py
index 2e71f3d00..a788988ba 100644
--- a/tutorials/preprocessing/90_eyetracking_data.py
+++ b/tutorials/preprocessing/90_eyetracking_data.py
@@ -31,9 +31,10 @@ the ERP and pupil response to the light flashes (i.e. the pupillary light reflex
 import mne
 from mne.datasets.eyelink import data_path
 from mne.preprocessing.eyetracking import read_eyelink_calibration
+from mne.viz.eyetracking import plot_gaze
 
-et_fpath = data_path() / "sub-01_task-plr_eyetrack.asc"
-eeg_fpath = data_path() / "sub-01_task-plr_eeg.mff"
+et_fpath = data_path() / "eeg-et" / "sub-01_task-plr_eyetrack.asc"
+eeg_fpath = data_path() / "eeg-et" / "sub-01_task-plr_eeg.mff"
 
 raw_et = mne.io.read_raw_eyelink(et_fpath, create_annotations=["blinks"])
 raw_eeg = mne.io.read_raw_egi(eeg_fpath, preload=True, verbose="warning")
@@ -123,7 +124,9 @@ raw_et.plot(scalings=dict(eyegaze=1e3))
 # window 50 ms before and 200 ms after the blink, so that the noisy data surrounding
 # the blink is also interpolated.
 
-mne.preprocessing.eyetracking.interpolate_blinks(raw_et, buffer=(0.05, 0.2))
+mne.preprocessing.eyetracking.interpolate_blinks(
+    raw_et, buffer=(0.05, 0.2), interpolate_gaze=True
+)
 
 # %%
 # .. important:: By default, :func:`~mne.preprocessing.eyetracking.interpolate_blinks`,
@@ -176,6 +179,7 @@ mne.preprocessing.realign_raw(
 )
 # Add EEG channels to the eye-tracking raw object
 raw_et.add_channels([raw_eeg], force_update_info=True)
+del raw_eeg  # free up some memory
 
 # Define a few channel groups of interest and plot the data
 frontal = ["E19", "E11", "E4", "E12", "E5"]
@@ -194,11 +198,25 @@ raw_et.plot(events=et_events, event_id=event_dict, event_color="g", order=picks_
 # Now let's extract epochs around our flash events. We should see a clear pupil
 # constriction response to the flashes.
 
-epochs = mne.Epochs(raw_et, events=et_events, event_id=event_dict, tmin=-0.3, tmax=3)
+# Skip baseline correction for now. We will apply baseline correction later.
+epochs = mne.Epochs(
+    raw_et, events=et_events, event_id=event_dict, tmin=-0.3, tmax=3, baseline=None
+)
+del raw_et  # free up some memory
 epochs[:8].plot(events=et_events, event_id=event_dict, order=picks_idx)
 
+# %%
+# For this experiment, the participant was instructed to fixate on a crosshair in the
+# center of the screen. Let's plot the gaze position data to confirm that the
+# participant primarily kept their gaze fixated at the center of the screen.
+
+plot_gaze(epochs, width=1920, height=1080)
+
+# %%
+# .. seealso:: :ref:`tut-eyetrack-heatmap`
+
 # %%
 # Finally, let's plot the evoked responses to the light flashes to get a sense of the
 # average pupillary light response, and the associated ERP in the EEG data.
 
-epochs.average().plot(picks=occipital + pupil)
+epochs.apply_baseline().average().plot(picks=occipital + pupil)
diff --git a/tutorials/raw/40_visualize_raw.py b/tutorials/raw/40_visualize_raw.py
index 8ce101a2e..28fbfa8fa 100644
--- a/tutorials/raw/40_visualize_raw.py
+++ b/tutorials/raw/40_visualize_raw.py
@@ -13,7 +13,6 @@ As usual we'll start by importing the modules we need, loading some
 :ref:`example data <sample-dataset>`, and cropping the `~mne.io.Raw`
 object to just 60 seconds before loading it into RAM to save memory:
 """
-
 # %%
 
 import os


Test patch:
diff --git a/mne/viz/eyetracking/tests/__init__.py b/mne/viz/eyetracking/tests/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/mne/viz/eyetracking/tests/test_heatmap.py b/mne/viz/eyetracking/tests/test_heatmap.py
new file mode 100644
index 000000000..99103c552
--- /dev/null
+++ b/mne/viz/eyetracking/tests/test_heatmap.py
@@ -0,0 +1,35 @@
+# Authors: Scott Huberty <seh33@uw.edu>
+#
+# License: Simplified BSD
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pytest
+
+import mne
+
+
+@pytest.mark.parametrize("axes", [None, True])
+def test_plot_heatmap(axes):
+    """Test plot_gaze."""
+    # Create a toy epochs instance
+    info = info = mne.create_info(
+        ch_names=["xpos", "ypos"], sfreq=100, ch_types="eyegaze"
+    )
+    # simulate a steady fixation at the center of the screen
+    width, height = (1920, 1080)
+    shape = (1, 100)  # x or y, time
+    data = np.vstack([np.full(shape, width / 2), np.full(shape, height / 2)])
+    epochs = mne.EpochsArray(data[None, ...], info)
+    epochs.info["chs"][0]["loc"][4] = -1
+    epochs.info["chs"][1]["loc"][4] = 1
+
+    if axes:
+        axes = plt.subplot()
+    fig = mne.viz.eyetracking.plot_gaze(
+        epochs, width=width, height=height, axes=axes, cmap="Greys", sigma=None
+    )
+    img = fig.axes[0].images[0].get_array()
+    # We simulated a 2D histogram where only the central pixel (960, 540) was active
+    assert img.T[width // 2, height // 2] == 1  # central pixel is active
+    assert np.sum(img) == 1  # only the central pixel should be active