improve resulting displacement vectors (#657)

* cart2pol now always sets phi=0 when rho=0

* Improve displace vectors computation

* deprecate movement.kinematics.compute_displacement
* add movement.kinematics.compute_forward_displacement
* add movement.kinematics.compute_backward_displacement
* adapt examples and docstrings accordingly
* adapt tests to the implementation changes

* update compute_kinematics example

shows the existance of compute_path_length

* Add backward displacement to the example. Small edits to wording and flow

* Fix double ticks

* Deprecation warning after short summary following numpydoc style

https://numpydoc.readthedocs.io/en/latest/format.html#deprecation-warning

* Factor out deprecation test for kinematics as a general test

* Clarify handling of unsigned zeros in cart2pol

* Fix rendering bits

* Fix broken SLEAP link

---------

Co-authored-by: sfmig <[email protected]>

Author: carlocastoldi

examples/compute_kinematics.py

@@ -1,4 +1,4 @@
-"""Compute and visualise kinematics.
+"""Compute and visualise kinematics
 ====================================
 
 Compute displacement, velocity and acceleration, and
@@ -12,6 +12,7 @@
 # For interactive plots: install ipympl with `pip install ipympl` and uncomment
 # the following line in your notebook
 # %matplotlib widget
+import numpy as np
 from matplotlib import pyplot as plt
 
 import movement.kinematics as kin
@@ -50,9 +51,8 @@
 # ---------------------------
 # First, let's visualise the trajectories of the mice in the XY plane,
 # colouring them by individual.
-# We use :func:`movement.plots.plot_centroid_trajectory` which is a wrapper
-# around :func:`matplotlib.pyplot.scatter` that simplifies plotting the
-# trajectories of individuals in the dataset.
+# For this we can use :func:`movement.plots.plot_centroid_trajectory`
+# which is a wrapper around :func:`matplotlib.pyplot.scatter`.
 # The fig and ax objects returned can be used to further customise the plot.
 
 # Create a single figure and axes
@@ -77,6 +77,8 @@
         label=mouse_name,
     )
     ax.legend().set_alpha(1)
+ax.set_xlabel("x (pixels)")
+ax.set_ylabel("y (pixels)")
 fig.show()
 
 # %%
@@ -86,21 +88,31 @@
 # follows the convention for SLEAP and most image processing tools.
 
 # %%
-# By default :func:`plot_centroid_trajectory()<movement.plots.\
-# plot_centroid_trajectory>` colours data points based on their timestamps:
-fig, axes = plt.subplots(3, 1, sharey=True)
-for mouse_name, ax in zip(position.individuals.values, axes, strict=False):
+# We can also plot the trajectories of the mice in the XY plane independently,
+# colouring the data points based on their timestamps. This is the default
+# behaviour of
+# :func:`plot_centroid_trajectory()<movement.plots.plot_centroid_trajectory>`
+# when the ``c`` argument is not provided:
+fig, axes = plt.subplots(2, 2, sharey=True)
+for mouse_name, ax in zip(
+    position.individuals.values, axes.flat, strict=False
+):
     ax.invert_yaxis()
     fig, ax = plot_centroid_trajectory(
         position,
         individual=mouse_name,
         ax=ax,
         s=2,
     )
+    ax.set_aspect("equal")
+    ax.set_xlim(150, 1250)
+    ax.set_ylim(500, 1100)
     ax.set_title(f"Trajectory {mouse_name}")
     ax.set_xlabel("x (pixels)")
     ax.set_ylabel("y (pixels)")
     ax.collections[0].colorbar.set_label("Time (frames)")
+# Hide the unused subplot (4th one)
+axes[1, 1].set_visible(False)
 fig.tight_layout()
 fig.show()
 
@@ -111,7 +123,7 @@
 # the third mouse (``AEON3B_TP2``) followed an anti-clockwise direction.
 
 # %%
-# We can also easily plot the components of the position vector against time
+# We can also inspect the components of the position vector against time
 # using ``xarray``'s built-in plotting methods. We use
 # :meth:`xarray.DataArray.squeeze` to
 # remove the dimension of length 1 from the data (the ``keypoints`` dimension).
@@ -124,37 +136,36 @@
 # and the ``x`` and ``y`` coordinates of the ``position`` are in pixels.
 
 # %%
-# Compute displacement
-# ---------------------
+# Compute displacement vectors
+# ----------------------------
 # The :mod:`movement.kinematics` module
-# provides functions to compute various kinematic quantities,
-# such as displacement, velocity, and acceleration.
-# We can start off by computing the distance travelled by the mice along
-# their trajectories:
-
-displacement = kin.compute_displacement(position)
+# provides functions to compute various kinematic variables,
+# such as displacement, velocity, and acceleration. Below we showcase
+# how these functions can be used.
+#
+# We can compute the forward displacement vectors as follows:
+forward_displacement = kin.compute_forward_displacement(position)
 
 # %%
-# The :func:`movement.kinematics.compute_displacement`
+# The :func:`movement.kinematics.compute_forward_displacement`
 # function will return a data array equivalent to the ``position`` one,
-# but holding displacement data along the ``space`` axis, rather than
-# position data.
+# but holding displacement data along the ``space`` axis.
 #
-# The ``displacement`` data array holds, for a given individual and keypoint
-# at timestep ``t``, the vector that goes from its previous position at time
-# ``t-1`` to its current position at time ``t``.
+# The ``forward_displacement`` data array holds, for a given individual and
+# keypoint at timestep ``t``, the vector that goes from its current position
+# at time ``t`` to its next position at time ``t+1``.
 
 # %%
-# And what happens at ``t=0``, since there is no previous timestep?
-# We define the displacement vector at time ``t=0`` to be the zero vector.
-# This way the shape of the ``displacement`` data array is the
-# same as the  ``position`` array:
+# And what happens in the last timestep, when there is no next timepoint?
+# We define the forward displacement vector then to be the
+# zero vector. This way the shape of the ``forward_displacement`` data array
+# is the same as the ``position`` array:
 print(f"Shape of position: {position.shape}")
-print(f"Shape of displacement: {displacement.shape}")
+print(f"Shape of displacement: {forward_displacement.shape}")
 
 # %%
-# We can visualise these displacement vectors with a quiver plot. In this case
-# we focus on the mouse ``AEON3B_TP2``:
+# We can visualise the forward displacement vectors with a quiver plot. In
+# this case we focus on the mouse ``AEON3B_TP2``:
 mouse_name = "AEON3B_TP2"
 
 fig = plt.figure()
@@ -169,12 +180,12 @@
     cmap="viridis",
 )
 
-# plot displacement vectors: at t, vector from t-1 to t
+# plot forward displacement vectors: at t, vector from t to t+1
 ax.quiver(
     position.sel(individuals=mouse_name, space="x"),
     position.sel(individuals=mouse_name, space="y"),
-    displacement.sel(individuals=mouse_name, space="x"),
-    displacement.sel(individuals=mouse_name, space="y"),
+    forward_displacement.sel(individuals=mouse_name, space="x"),
+    forward_displacement.sel(individuals=mouse_name, space="y"),
     angles="xy",
     scale=1,
     scale_units="xy",
@@ -183,39 +194,37 @@
     headaxislength=9,
 )
 
-ax.axis("equal")
-ax.set_xlim(450, 575)
-ax.set_ylim(950, 1075)
+ax.set_xlim(480, 600)
+ax.set_ylim(980, 1080)
 ax.set_xlabel("x (pixels)")
 ax.set_ylabel("y (pixels)")
 ax.set_title(f"Zoomed in trajectory of {mouse_name}")
 ax.invert_yaxis()
 fig.colorbar(sc, ax=ax, label="time (s)")
 
+
 # %%
-# Notice that this figure is not very useful as a visual check:
-# we can see that there are vectors defined for each point in
-# the trajectory, but we have no easy way to verify they are indeed
-# the displacement vectors from ``t-1`` to ``t``.
+# We can visually verify that indeed the forward displacement vector
+# connects the previous and current positions as expected.
 
 # %%
-# If instead we plot
-# the opposite of the displacement vector, we will see that at every time
-# ``t``, the vectors point to the position at ``t-1``.
-# Remember that the displacement vector is defined as the vector at
-# time ``t``, that goes from the previous position ``t-1`` to the
-# current position at ``t``. Therefore, the opposite vector will point
-# from the position point at ``t``, to the position point at ``t-1``.
+# Similarly, with :func:`movement.kinematics.compute_backward_displacement`
+# we can compute the backward displacement vectors, which connect the current
+# position to the previous one:
+backward_displacement = kin.compute_backward_displacement(position)
 
 # %%
-# We can easily do this by flipping the sign of the displacement vector in
-# the plot above:
-mouse_name = "AEON3B_TP2"
+# In this case, the backward displacement vector at the first timestep
+# is defined as the zero vector, since there is no previous position.
+
+# %%
+# Adapting the code snippet from above, we can visually check that the
+# backward displacement vector is indeed the reverse of the forward
+# displacement vector.
 
 fig = plt.figure()
 ax = fig.add_subplot()
 
-# plot position data
 sc = ax.scatter(
     position.sel(individuals=mouse_name, space="x"),
     position.sel(individuals=mouse_name, space="y"),
@@ -224,22 +233,21 @@
     cmap="viridis",
 )
 
-# plot displacement vectors: at t, vector from t-1 to t
 ax.quiver(
     position.sel(individuals=mouse_name, space="x"),
     position.sel(individuals=mouse_name, space="y"),
-    -displacement.sel(individuals=mouse_name, space="x"),  # flipped sign
-    -displacement.sel(individuals=mouse_name, space="y"),  # flipped sign
+    backward_displacement.sel(individuals=mouse_name, space="x"),
+    backward_displacement.sel(individuals=mouse_name, space="y"),
     angles="xy",
     scale=1,
     scale_units="xy",
     headwidth=7,
     headlength=9,
     headaxislength=9,
 )
-ax.axis("equal")
-ax.set_xlim(450, 575)
-ax.set_ylim(950, 1075)
+
+ax.set_xlim(480, 600)
+ax.set_ylim(980, 1080)
 ax.set_xlabel("x (pixels)")
 ax.set_ylabel("y (pixels)")
 ax.set_title(f"Zoomed in trajectory of {mouse_name}")
@@ -248,30 +256,65 @@
 
 
 # %%
-# Now we can visually verify that indeed the displacement vector
-# connects the previous and current positions as expected.
+# Compute path length
+# --------------------
+# We can compute the distance travelled by the
+# mouse as the sum of the lengths of all
+# displacement vectors along its trajectory.
+# Both backward and forward displacement vectors
+# should give the same result:
+
+# length of each forward displacement vector
+forward_displacement_lengths = compute_norm(
+    forward_displacement.sel(individuals=mouse_name)
+)
 
-# %%
-# With the displacement data we can compute the distance travelled by the
-# mouse along its trajectory.
+# length of each backward displacement vector
+backward_displacement_lengths = compute_norm(
+    backward_displacement.sel(individuals=mouse_name)
+)
 
-# length of each displacement vector
-displacement_vectors_lengths = compute_norm(
-    displacement.sel(individuals=mouse_name)
+# check their lengths are the same
+np.testing.assert_almost_equal(
+    forward_displacement_lengths.values[:-1],  # exclude last timestep
+    backward_displacement_lengths.values[1:],  # exclude first timestep
 )
 
 # sum the lengths of all displacement vectors (in pixels)
-total_displacement = displacement_vectors_lengths.sum(dim="time").values[0]
+total_displacement_fwd = forward_displacement_lengths.sum(dim="time").values[0]
+total_displacement_bwd = backward_displacement_lengths.sum(dim="time").values[
+    0
+]
 
 print(
-    f"The mouse {mouse_name}'s trajectory is {total_displacement:.2f} "
-    "pixels long"
+    f"The mouse {mouse_name}'s path length is {total_displacement_fwd:.2f} "
+    "pixels long (using forward displacement)"
 )
+print(
+    f"The mouse {mouse_name}'s path length is {total_displacement_bwd:.2f} "
+    "pixels long (using backward displacement)"
+)
+
+
+# %%
+# We provide a convenience function
+# :func:`movement.kinematics.compute_path_length`
+# to compute the path length for all individuals and keypoints in a position
+# data array. We can verify that using this function gives the same result
+# as before for the ``AEON3B_TP2`` mouse:
+
+path_lengths = kin.compute_path_length(ds.position)
+
+for mouse_name in path_lengths.individuals.values:
+    print(
+        f"Path length for {mouse_name}: "
+        f"{path_lengths.sel(individuals=mouse_name).values[0]:.2f} pixels"
+    )
 
 # %%
 # Compute velocity
 # ----------------
-# We can easily compute the velocity vectors for all individuals in our data
+# We can also compute the velocity vectors for all individuals in our data
 # array:
 velocity = kin.compute_velocity(position)
 
@@ -294,9 +337,9 @@
 # %%
 # The components of the velocity vector seem noisier than the components of
 # the position vector.
-# This is expected, since we are deriving the velocity using differences in
+# This is expected, since we are estimating the velocity using differences in
 # position (which is somewhat noisy), over small stepsizes.
-# More specifically, we use numpy's gradient implementation, which
+# More specifically, we use :func:`numpy.gradient` internally, which
 # uses second order central differences.
 
 # %%
@@ -315,7 +358,7 @@
 
 # %%
 # To visualise the direction of the velocity vector at each timestep, we can
-# use a quiver plot:
+# again use a quiver plot:
 mouse_name = "AEON3B_TP2"
 fig = plt.figure()
 ax = fig.add_subplot()

movement/io/load_poses.py

@@ -287,7 +287,7 @@ def from_sleap_file(
     ----------
     .. [1] https://docs.sleap.ai/latest/learnings/export-analysis/
     .. [2] https://github.com/talmolab/sleap/blob/v1.3.3/sleap/info/write_tracking_h5.py#L59
-    .. [3] https://docs.sleap.ai/latest/how-to-guides/tracking-and-proofreading/
+    .. [3] https://docs.sleap.ai/latest/guides/tracking-and-proofreading/
 
     Examples
     --------

movement/kinematics/__init__.py

@@ -4,6 +4,8 @@
 from movement.kinematics.kinematics import (
     compute_acceleration,
     compute_displacement,
+    compute_forward_displacement,
+    compute_backward_displacement,
     compute_path_length,
     compute_speed,
     compute_time_derivative,
@@ -18,6 +20,8 @@
 
 __all__ = [
     "compute_displacement",
+    "compute_forward_displacement",
+    "compute_backward_displacement",
     "compute_velocity",
     "compute_acceleration",
     "compute_speed",