.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/compute_angular_derivatives.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code or to run this example in your browser via Binder. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_compute_angular_derivatives.py: Compute angular velocity and acceleration ============================================ Compute the head angular velocity and acceleration. .. GENERATED FROM PYTHON SOURCE LINES 8-10 Imports ------- .. GENERATED FROM PYTHON SOURCE LINES 10-30 .. code-block:: Python # For interactive plots: install ipympl with `pip install ipympl` and uncomment # the following line in your notebook # %matplotlib widget import xarray as xr from matplotlib import pyplot as plt from movement import sample_data from movement.filtering import ( filter_by_confidence, interpolate_over_time, rolling_filter, ) from movement.kinematics import ( compute_angular_acceleration, compute_angular_velocity, compute_forward_vector, ) .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 34-45 Before you start ---------------- This example builds on the :ref:`sphx_glr_examples_compute_head_direction.py` example, which explains how ``movement`` represents orientations as vectors and as angles. We recommend reading that one first. Here we focus on how fast an orientation changes over time, i.e. its *angular velocity* and *angular acceleration*. While this example uses the head direction of a freely moving mouse, the same functions can be used with any 2D orientation represented as vectors. .. GENERATED FROM PYTHON SOURCE LINES 47-52 Preparing the data ------------------ As in :ref:`sphx_glr_examples_compute_head_direction.py`, we use a single-mouse dataset consisting of a top-down video of a mouse exploring an elevated plus maze. .. GENERATED FROM PYTHON SOURCE LINES 52-57 .. code-block:: Python ds = sample_data.fetch_dataset("DLC_single-mouse_EPM.predictions.h5") print(f"fps: {ds.fps}") ds .. rst-class:: sphx-glr-script-out .. code-block:: none fps: 30.0 .. raw:: html 
<xarray.Dataset> Size: 4MB
    Dimensions:     (time: 18485, space: 2, keypoint: 8, individual: 1)
    Coordinates:
      * time        (time) float64 148kB 0.0 0.03333 0.06667 ... 616.1 616.1 616.1
      * space       (space) <U1 8B 'x' 'y'
      * keypoint    (keypoint) <U13 416B 'snout' 'left_ear' ... 'tail_end'
      * individual  (individual) <U4 16B 'id_0'
    Data variables:
        position    (time, space, keypoint, individual) float64 2MB 508.4 ... 390.5
        confidence  (time, keypoint, individual) float64 1MB 0.0002829 ... 0.9978
    Attributes: (6)


.. GENERATED FROM PYTHON SOURCE LINES 58-65 Next we :ref:`remove low-confidence predictions `, apply a :ref:`rolling median filter ` to suppress tracking outliers, and interpolate to fill short gaps in the data. Each operation here overwrites the ``position`` variable in the dataset. We select that variable and drop the "individual" dimension for convenience, as we are only working with one mouse here. .. GENERATED FROM PYTHON SOURCE LINES 65-85 .. code-block:: Python # preprocessing ds["position"] = filter_by_confidence( ds.position, ds.confidence, threshold=0.9, ) ds["position"] = rolling_filter( ds.position, window=5, min_periods=2, statistic="median", ) ds["position"] = interpolate_over_time(ds.position, max_gap=2) # remove the "individual" dimension position = ds.position.drop_vars("individual").squeeze("individual") position .. raw:: html 
<xarray.DataArray 'position' (time: 18485, space: 2, keypoint: 8)> Size: 2MB
    nan nan nan nan nan nan nan nan ... 479.3 489.0 443.2 436.7 456.8 381.6 390.4
    Coordinates:
      * time      (time) float64 148kB 0.0 0.03333 0.06667 0.1 ... 616.1 616.1 616.1
      * space     (space) <U1 8B 'x' 'y'
      * keypoint  (keypoint) <U13 416B 'snout' 'left_ear' ... 'tailbase' 'tail_end'
    Attributes: (1)


.. GENERATED FROM PYTHON SOURCE LINES 86-89 We compute the head vector from the left and right ear keypoints, as in the :ref:`head direction example `. .. GENERATED FROM PYTHON SOURCE LINES 89-98 .. code-block:: Python head_vector = compute_forward_vector( position, left_keypoint="left_ear", right_keypoint="right_ear", camera_view="top_down", ) head_vector .. raw:: html 
<xarray.DataArray 'forward_vector' (time: 18485, space: 2)> Size: 296kB
    nan nan nan nan nan nan nan ... 0.9424 0.4128 0.9108 0.4108 0.9117 0.4108 0.9117
    Coordinates:
      * time     (time) float64 148kB 0.0 0.03333 0.06667 0.1 ... 616.1 616.1 616.1
      * space    (space) <U1 8B 'x' 'y'
    Attributes: (1)


.. GENERATED FROM PYTHON SOURCE LINES 99-102 The ``head_vector`` array contains unit vectors (i.e. of length 1). It has the same number of time points as the ``position`` array it was derived from, and two spatial dimensions (x and y). .. GENERATED FROM PYTHON SOURCE LINES 104-108 Compute angular velocity ------------------------ We compute the angular head velocity (AHV) by passing the ``head_vector`` array to :func:`~movement.kinematics.compute_angular_velocity()`. .. GENERATED FROM PYTHON SOURCE LINES 108-112 .. code-block:: Python ahv = compute_angular_velocity(head_vector) ahv .. raw:: html 
<xarray.DataArray 'angular_velocity' (time: 18485)> Size: 148kB
    nan nan nan nan nan nan nan ... 0.4809 1.054 1.582 -0.4557 -1.234 0.03325 0.0
    Coordinates:
      * time     (time) float64 148kB 0.0 0.03333 0.06667 0.1 ... 616.1 616.1 616.1


.. GENERATED FROM PYTHON SOURCE LINES 113-118 We get a 1D time-varying array in radians per time unit (here seconds). To return angular velocity in degrees per unit time instead, pass ``in_degrees=True``. We can also assign metadata attributes to describe the variable and its units; these are useful for plotting. .. GENERATED FROM PYTHON SOURCE LINES 118-124 .. code-block:: Python ahv_deg = compute_angular_velocity(head_vector, in_degrees=True) ahv_deg.attrs["long_name"] = "Angular head velocity" ahv_deg.attrs["units"] = "deg/s" ahv_deg .. raw:: html 
<xarray.DataArray 'angular_velocity' (time: 18485)> Size: 148kB
    nan nan nan nan nan nan nan nan ... 27.55 60.38 90.63 -26.11 -70.68 1.905 0.0
    Coordinates:
      * time     (time) float64 148kB 0.0 0.03333 0.06667 0.1 ... 616.1 616.1 616.1
    Attributes: (2)


.. GENERATED FROM PYTHON SOURCE LINES 125-129 Let's plot two histograms of angular head velocities: one with a linear y-axis and one with a logarithmic y-axis. We expect the distribution to be symmetric about zero and heavy-tailed. The logarithmic scale makes the tails easier to visualise. .. GENERATED FROM PYTHON SOURCE LINES 129-140 .. code-block:: Python N_BINS = 100 fig, axes = plt.subplots(1, 2, figsize=(8, 3), sharex=True) ahv_deg.plot.hist(bins=N_BINS, ax=axes[0], color="grey") ahv_deg.plot.hist(bins=N_BINS, ax=axes[1], yscale="log", color="grey") axes[0].set_title("Linear y-axis") axes[1].set_title("Logarithmic y-axis") plt.tight_layout() plt.show() .. image-sg:: /examples/images/sphx_glr_compute_angular_derivatives_001.png :alt: Linear y-axis, Logarithmic y-axis :srcset: /examples/images/sphx_glr_compute_angular_derivatives_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 141-153 .. admonition:: Sign convention :class: note Positive angular velocities correspond to rotations from the positive x-axis towards the positive y-axis. In image coordinates (y pointing down) with a top-down camera, that is a *clockwise* rotation on screen. We see a few values beyond :math:`\pm` 1000 deg/s, and some even beyond :math:`\pm` 2000 deg/s. Mice can turn their heads very quickly, but some of these extreme values are likely due to noise. Next we'll look at ways to smooth the AHV estimate. .. GENERATED FROM PYTHON SOURCE LINES 155-190 Smooth angular velocity ----------------------- The angular velocity signal is noisy, because differentiation amplifies small frame-to-frame jitter. The preprocessing we have already applied to the ``position`` data helps with that, but may not be sufficient. Here we look at two further places to smooth the data: **A. Pre-smoothing.** Smooth ``head_vector`` with a rolling mean before differentiating. Since ``head_vector`` contains unit vectors, averaging them is the correct way to average directions: the averaged vector points in the mean direction of the window, also known as the *circular mean* of the head angles. The averaged vectors are generally shorter than 1, but this does not matter, since only their direction is used. Averaging the angles themselves would fail near :math:`\pm\pi`; for example, 179Ā° and -179Ā° would average to 0Ā° (pointing the opposite way) instead of 180Ā°. **B. Windowed estimate.** Passing ``window`` (an odd number of frames) to :func:`~movement.kinematics.compute_angular_velocity()` fits a straight line to the head angle within a centred window of that length and returns its slope. This **least-squares estimate** requires data to be uniformly sampled in time. The estimate is only computed where the full window contains valid (non-NaN) values. As a result, the first and last ``window // 2`` frames, and any frames whose window contains a missing value, are assigned NaN. .. note:: Both approaches widen existing gaps in the data. By default, :func:`~movement.filtering.rolling_filter` returns NaN wherever its window contains a NaN, and so does the least-squares estimate. Each gap therefore grows by ``window // 2`` frames on each side, per step. Pass ``min_periods`` to ``rolling_filter`` to relax this. Let's compare the two approaches, and their combination. Feel free to experiment with ``WINDOW_VEC``, ``WINDOW_LSQ`` and ``TIME_SLICE`` to see how they affect the results. .. GENERATED FROM PYTHON SOURCE LINES 190-228 .. code-block:: Python WINDOW_VEC = 5 # in frames, used for pre-smoothing the vector WINDOW_LSQ = 5 # in frames, used for the least-squares slope estimation TIME_SLICE = slice(80, 90) # in seconds, controls what's plotted head_vector_smooth = rolling_filter(head_vector, WINDOW_VEC, statistic="mean") ahv_presmoothed = compute_angular_velocity(head_vector_smooth, in_degrees=True) ahv_lsq = compute_angular_velocity( head_vector, window=WINDOW_LSQ, in_degrees=True ) ahv_presmoothed_lsq = compute_angular_velocity( head_vector_smooth, window=WINDOW_LSQ, in_degrees=True ) label_to_da = { "raw": ahv_deg, "A: smoothed vector": ahv_presmoothed, "B: least-squares": ahv_lsq, "A + B": ahv_presmoothed_lsq, } colors = ["grey", "tab:blue", "tab:red", "tab:purple"] fig, axes = plt.subplots(3, 1, figsize=(8, 6), sharex=True, sharey=True) for i, ax in enumerate(axes): title, da = list(label_to_da.items())[i + 1] color = colors[i + 1] ahv_deg.sel(time=TIME_SLICE).plot(ax=ax, color=colors[0], alpha=0.5) da.sel(time=TIME_SLICE).plot(ax=ax, color=color, alpha=0.8) ax.set_xlim(TIME_SLICE.start, TIME_SLICE.stop) ax.set_title(title, color=color) ax.set_xlabel("") ax.set_ylabel(ahv_deg.attrs["units"]) axes[-1].set_xlabel("Time (s)") fig.suptitle(f"{ahv_deg.attrs['long_name']} (raw vs smoothed)") fig.tight_layout() plt.show() .. image-sg:: /examples/images/sphx_glr_compute_angular_derivatives_002.png :alt: Angular head velocity (raw vs smoothed), A: smoothed vector, B: least-squares, A + B :srcset: /examples/images/sphx_glr_compute_angular_derivatives_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 229-232 All three approaches damp extreme values, while preserving the timing of the main peaks and troughs. To compare how much they damp, let's look at the distribution of angular head velocities before and after smoothing. .. GENERATED FROM PYTHON SOURCE LINES 232-250 .. code-block:: Python raw_range = (ahv_deg.min().item(), ahv_deg.max().item()) # same bins for all fig, axes = plt.subplots(1, 4, figsize=(8, 3), sharex=True, sharey=True) for i, ax in enumerate(axes): title, da = list(label_to_da.items())[i] da.plot.hist( ax=ax, bins=N_BINS, range=raw_range, yscale="log", color=colors[i] ) p99 = abs(da).quantile(0.99, skipna=True).item() ax.set_title( f"{title}\n99% within Ā±{p99:.0f} deg/s", color=colors[i], fontsize=10 ) ax.set_xlabel(ahv_deg.attrs["units"]) fig.suptitle(f"Distribution of {ahv_deg.attrs['long_name']}") fig.tight_layout() plt.show() .. image-sg:: /examples/images/sphx_glr_compute_angular_derivatives_003.png :alt: Distribution of Angular head velocity, raw 99% within Ā±627 deg/s, A: smoothed vector 99% within Ā±474 deg/s, B: least-squares 99% within Ā±607 deg/s, A + B 99% within Ā±542 deg/s :srcset: /examples/images/sphx_glr_compute_angular_derivatives_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 251-271 The two approaches reduce the overall spread of the distribution by a similar amount, but differ in the tails. Pre-smoothing (A) damps typical values more (see the range containing 99% of values in each title), but some extreme values beyond :math:`\pm` 2000 deg/s survive. The least-squares estimate (B) removes these, but damps typical values less. Combining both (A + B) removes the extreme values too, while damping typical values more than B alone. .. admonition:: Centred vs trailing window estimate :class: note ``movement`` uses *centred* windows, so smoothing does not shift the signal in time. Some applications, however, require a derivative estimate based on a window ending at each time point. We can obtain this trailing estimate by shifting the centred estimate forward by ``window // 2`` frames: .. code-block:: python ahv_lsq_trailing = ahv_lsq.shift(time=WINDOW_LSQ // 2) .. GENERATED FROM PYTHON SOURCE LINES 273-289 Compare to published values --------------------------- `Hennestad et al. (2021) `_ analysed published data of freely moving mice (`Laurens et al., 2019 `_) and found that about 90% of all head rotations are slower than 180 deg/s. As a sanity check, let's compute the same fraction for the A + B estimate, ignoring missing values. .. note:: The experimental setups differ, so we only expect a rough match. In Laurens et al., mice foraged in a circular arena (50 cm diameter), and head direction was tracked from two LEDs mounted on the headstage. Here, a mouse explores an elevated plus maze, and head direction is derived from ear keypoints estimated by DeepLabCut and prepocessed with ``movement``. .. GENERATED FROM PYTHON SOURCE LINES 289-294 .. code-block:: Python abs_ahv = abs(ahv_presmoothed_lsq) frac_below_180 = (abs_ahv < 180).sum() / abs_ahv.notnull().sum() print(f"AHV below 180 deg/s: {frac_below_180.item():.1%}") .. rst-class:: sphx-glr-script-out .. code-block:: none AHV below 180 deg/s: 91.7% .. GENERATED FROM PYTHON SOURCE LINES 295-300 This closely matches the published value, which suggests that the smoothed estimate is physiologically plausible. Let's visualise this on a histogram of the A + B estimate, this time with a linear y-axis. To focus on the bulk of the distribution, we only plot values within :math:`\pm` 500 deg/s. .. GENERATED FROM PYTHON SOURCE LINES 300-313 .. code-block:: Python fig, ax = plt.subplots(figsize=(5, 3)) ahv_presmoothed_lsq.plot.hist( ax=ax, bins=N_BINS, range=(-500, 500), color=colors[3] ) for x in (-180, 180): ax.axvline(x, color="black", linestyle="--") ax.set_title(f"A + B: {frac_below_180.item():.1%} within Ā±180 deg/s") ax.set_xlabel(ahv_deg.attrs["units"]) fig.tight_layout() plt.show() .. image-sg:: /examples/images/sphx_glr_compute_angular_derivatives_004.png :alt: A + B: 91.7% within Ā±180 deg/s :srcset: /examples/images/sphx_glr_compute_angular_derivatives_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 315-329 Compute angular acceleration ---------------------------- :func:`~movement.kinematics.compute_angular_acceleration()` takes the same arguments as :func:`~movement.kinematics.compute_angular_velocity()`. If ``window`` is supplied, it fits a quadratic to the head angle over the window and returns its second derivative. For this, ``window`` must be odd and greater than 2. Differentiating twice amplifies frame-to-frame jitter even more than differentiating once, so the raw angular acceleration is often dominated by noise. Smoothing is therefore even more important here than for angular velocity. Let's compare the raw estimate with the least-squares estimate (B) and its combination with pre-smoothing (A + B). Feel free to experiment with ``WINDOW_LSQ_ACC`` to see how it affects the results. .. GENERATED FROM PYTHON SOURCE LINES 329-364 .. code-block:: Python WINDOW_LSQ_ACC = 7 # in frames, used for the least-squares quadratic fit aha_deg = compute_angular_acceleration(head_vector, in_degrees=True) aha_deg.attrs["long_name"] = "Angular head acceleration" aha_deg.attrs["units"] = "deg/sĀ²" aha_lsq = compute_angular_acceleration( head_vector, window=WINDOW_LSQ_ACC, in_degrees=True ) aha_presmoothed_lsq = compute_angular_acceleration( head_vector_smooth, window=WINDOW_LSQ_ACC, in_degrees=True ) label_to_da_acc = { "B: least-squares": (aha_lsq, colors[2]), "A + B": (aha_presmoothed_lsq, colors[3]), } fig, axes = plt.subplots(2, 1, figsize=(8, 4), sharex=True, sharey=True) for ax, (title, (da, color)) in zip( axes, label_to_da_acc.items(), strict=True ): aha_deg.sel(time=TIME_SLICE).plot(ax=ax, color=colors[0], alpha=0.5) da.sel(time=TIME_SLICE).plot(ax=ax, color=color, alpha=0.8) ax.set_xlim(TIME_SLICE.start, TIME_SLICE.stop) ax.set_title(title, color=color) ax.set_xlabel("") ax.set_ylabel(aha_deg.attrs["units"]) axes[-1].set_xlabel("Time (s)") fig.suptitle(f"{aha_deg.attrs['long_name']} (raw vs smoothed)") fig.tight_layout() plt.show() .. image-sg:: /examples/images/sphx_glr_compute_angular_derivatives_005.png :alt: Angular head acceleration (raw vs smoothed), B: least-squares, A + B :srcset: /examples/images/sphx_glr_compute_angular_derivatives_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 365-372 As for angular velocity, both estimates preserve the timing of the main peaks while strongly reducing their amplitude. Combining pre-smoothing with the least-squares estimate (A + B) damps the amplitude considerably more than the least-squares estimate alone (B). Bear in mind that smoothing damps genuine fast accelerations as well as tracking noise, so the window sizes should be chosen with the expected timescale of head movements in mind. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.248 seconds) .. _sphx_glr_download_examples_compute_angular_derivatives.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/neuroinformatics-unit/movement/main?filepath=notebooks/examples/compute_angular_derivatives.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: compute_angular_derivatives.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: compute_angular_derivatives.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: compute_angular_derivatives.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_