diff --git a/examples/advanced/compute_kinematics.ipynb b/examples/advanced/compute_kinematics.ipynb new file mode 100644 index 000000000..b00b4450d --- /dev/null +++ b/examples/advanced/compute_kinematics.ipynb @@ -0,0 +1,661 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "eff2e446f20126d1", + "metadata": {}, + "outputs": [], + "source": [ + "\"\"\"Compute and visualise kinematics\n", + "====================================\n", + "\n", + "Compute displacement, velocity and acceleration, and\n", + "visualise the results.\n", + "\"\"\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "95fd0e6e143c9313", + "metadata": {}, + "outputs": [], + "source": [ + "# Imports\n", + "# -------\n", + "\n", + "# For interactive plots: install ipympl with `pip install ipympl` and uncomment\n", + "# the following line in your notebook\n", + "%matplotlib widget\n", + "import numpy as np\n", + "from matplotlib import pyplot as plt\n", + "\n", + "import movement.kinematics as kin\n", + "from movement import sample_data\n", + "from movement.plots import plot_centroid_trajectory\n", + "from movement.utils.vector import compute_norm" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cfa097444531721c", + "metadata": {}, + "outputs": [], + "source": [ + "# Load sample dataset\n", + "# ------------------------\n", + "# First, we load an example dataset. In this case, we select the\n", + "# ``SLEAP_three-mice_Aeon_proofread`` sample data.\n", + "ds = sample_data.fetch_dataset(\n", + " \"SLEAP_three-mice_Aeon_proofread.analysis.h5\",\n", + ")\n", + "\n", + "print(ds)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "54669c1656427214", + "metadata": {}, + "outputs": [], + "source": [ + "# We can see in the printed description of the dataset ``ds`` that\n", + "# the data was acquired at 50 fps, and the time axis is expressed in seconds.\n", + "# It includes data for three individuals(``AEON3B_NTP``, ``AEON3B_TP1``,\n", + "# and ``AEON3B_TP2``), and only one keypoint called ``centroid`` was tracked\n", + "# in ``x`` and ``y`` dimensions.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a7bff37ac3163f17", + "metadata": {}, + "outputs": [], + "source": [ + "# The loaded dataset ``ds`` contains two data arrays:\n", + "# ``position`` and ``confidence``.\n", + "# To compute displacement, velocity and acceleration, we will need the\n", + "# ``position`` one:\n", + "position = ds.position" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27af8f02b26ab6a9", + "metadata": {}, + "outputs": [], + "source": [ + "# Visualise the data\n", + "# ---------------------------\n", + "# First, let's visualise the trajectories of the mice in the XY plane,\n", + "# colouring them by individual.\n", + "# For this we can use :func:`movement.plots.plot_centroid_trajectory`\n", + "# which is a wrapper around :func:`matplotlib.pyplot.scatter`.\n", + "# The fig and ax objects returned can be used to further customise the plot.\n", + "\n", + "# Create a single figure and axes\n", + "fig, ax = plt.subplots(1, 1)\n", + "# Invert y-axis so (0,0) is in the top-left,\n", + "# matching typical image coordinate systems\n", + "ax.invert_yaxis()\n", + "# Plot trajectories for each mouse on the same axes\n", + "for mouse_name, col in zip(\n", + " position.individuals.values,\n", + " [\"r\", \"g\", \"b\"], # colours\n", + " strict=False,\n", + "):\n", + " plot_centroid_trajectory(\n", + " position,\n", + " individual=mouse_name,\n", + " ax=ax, # Use the same axes for all plots\n", + " c=col,\n", + " marker=\"o\",\n", + " s=10,\n", + " alpha=0.2,\n", + " label=mouse_name,\n", + " )\n", + " ax.legend().set_alpha(1)\n", + "ax.set_xlabel(\"x (pixels)\")\n", + "ax.set_ylabel(\"y (pixels)\")\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d72e4c9a7643b6d0", + "metadata": {}, + "outputs": [], + "source": [ + "# We can see that the trajectories of the three mice are close to a circular\n", + "# arc. Notice that the x and y axes are set to equal scales, and that the\n", + "# origin of the coordinate system is at the top left of the image. This\n", + "# follows the convention for SLEAP and most image processing tools.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "899b58068e43f1da", + "metadata": {}, + "outputs": [], + "source": [ + "# We can also plot the trajectories of the mice in the XY plane independently,\n", + "# colouring the data points based on their timestamps. This is the default\n", + "# behaviour of\n", + "# :func:`plot_centroid_trajectory()`\n", + "# when the ``c`` argument is not provided:\n", + "fig, axes = plt.subplots(2, 2, sharey=True)\n", + "for mouse_name, ax in zip(\n", + " position.individuals.values, axes.flat, strict=False\n", + "):\n", + " ax.invert_yaxis()\n", + " fig, ax = plot_centroid_trajectory(\n", + " position,\n", + " individual=mouse_name,\n", + " ax=ax,\n", + " s=2,\n", + " )\n", + " ax.set_aspect(\"equal\")\n", + " ax.set_xlim(150, 1250)\n", + " ax.set_ylim(500, 1100)\n", + " ax.set_title(f\"Trajectory {mouse_name}\")\n", + " ax.set_xlabel(\"x (pixels)\")\n", + " ax.set_ylabel(\"y (pixels)\")\n", + " ax.collections[0].colorbar.set_label(\"Time (frames)\")\n", + "# Hide the unused subplot (4th one)\n", + "axes[1, 1].set_visible(False)\n", + "fig.tight_layout()\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ca71f2c8ce1e8c96", + "metadata": {}, + "outputs": [], + "source": [ + "# These plots show that for this snippet of the data,\n", + "# two of the mice (``AEON3B_NTP`` and ``AEON3B_TP1``)\n", + "# moved around the circle in clockwise direction, and\n", + "# the third mouse (``AEON3B_TP2``) followed an anti-clockwise direction.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa9590f1ddcff50d", + "metadata": {}, + "outputs": [], + "source": [ + "# We can also inspect the components of the position vector against time\n", + "# using ``xarray``'s built-in plotting methods. We use\n", + "# :meth:`xarray.DataArray.squeeze` to\n", + "# remove the dimension of length 1 from the data (the ``keypoints`` dimension).\n", + "position.squeeze().plot.line(x=\"time\", row=\"individuals\", aspect=2, size=2.5)\n", + "plt.gcf().show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7442cd10792dc504", + "metadata": {}, + "outputs": [], + "source": [ + "# If we use ``xarray``'s plotting function, the axes units are automatically\n", + "# taken from the data array. In our case, ``time`` is expressed in seconds,\n", + "# and the ``x`` and ``y`` coordinates of the ``position`` are in pixels.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b6389d2f455ec3ca", + "metadata": {}, + "outputs": [], + "source": [ + "# Compute displacement vectors\n", + "# ----------------------------\n", + "# The :mod:`movement.kinematics` module\n", + "# provides functions to compute various kinematic variables,\n", + "# such as displacement, velocity, and acceleration. Below we showcase\n", + "# how these functions can be used.\n", + "#\n", + "# We can compute the forward displacement vectors as follows:\n", + "forward_displacement = kin.compute_forward_displacement(position)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f9aaa4b7d358a35", + "metadata": {}, + "outputs": [], + "source": [ + "# The :func:`movement.kinematics.compute_forward_displacement`\n", + "# function will return a data array equivalent to the ``position`` one,\n", + "# but holding displacement data along the ``space`` axis.\n", + "#\n", + "# The ``forward_displacement`` data array holds, for a given individual and\n", + "# keypoint at timestep ``t``, the vector that goes from its current position\n", + "# at time ``t`` to its next position at time ``t+1``.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fecc78a0e6b09228", + "metadata": {}, + "outputs": [], + "source": [ + "# And what happens in the last timestep, when there is no next timepoint?\n", + "# We define the forward displacement vector then to be the\n", + "# zero vector. This way the shape of the ``forward_displacement`` data array\n", + "# is the same as the ``position`` array:\n", + "print(f\"Shape of position: {position.shape}\")\n", + "print(f\"Shape of displacement: {forward_displacement.shape}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "abcab41d78507719", + "metadata": {}, + "outputs": [], + "source": [ + "# We can visualise the forward displacement vectors with a quiver plot. In\n", + "# this case we focus on the mouse ``AEON3B_TP2``:\n", + "mouse_name = \"AEON3B_TP2\"\n", + "\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot()\n", + "\n", + "# plot position data\n", + "sc = ax.scatter(\n", + " position.sel(individuals=mouse_name, space=\"x\"),\n", + " position.sel(individuals=mouse_name, space=\"y\"),\n", + " s=15,\n", + " c=position.time,\n", + " cmap=\"viridis\",\n", + ")\n", + "\n", + "# plot forward displacement vectors: at t, vector from t to t+1\n", + "ax.quiver(\n", + " position.sel(individuals=mouse_name, space=\"x\"),\n", + " position.sel(individuals=mouse_name, space=\"y\"),\n", + " forward_displacement.sel(individuals=mouse_name, space=\"x\"),\n", + " forward_displacement.sel(individuals=mouse_name, space=\"y\"),\n", + " angles=\"xy\",\n", + " scale=1,\n", + " scale_units=\"xy\",\n", + " headwidth=7,\n", + " headlength=9,\n", + " headaxislength=9,\n", + ")\n", + "\n", + "ax.set_xlim(480, 600)\n", + "ax.set_ylim(980, 1080)\n", + "ax.set_xlabel(\"x (pixels)\")\n", + "ax.set_ylabel(\"y (pixels)\")\n", + "ax.set_title(f\"Zoomed in forward trajectory of {mouse_name}\")\n", + "ax.invert_yaxis()\n", + "sc.set_clim(8.8, 9.2)\n", + "fig.colorbar(sc, ax=ax, label=\"time (s)\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2b404bd7ae7bbc7", + "metadata": {}, + "outputs": [], + "source": [ + "# We can visually verify that indeed the forward displacement vector\n", + "# connects the previous and current positions as expected.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "df006963be556e", + "metadata": {}, + "outputs": [], + "source": [ + "# Similarly, with :func:`movement.kinematics.compute_backward_displacement`\n", + "# we can compute the backward displacement vectors, which connect the current\n", + "# position to the previous one:\n", + "backward_displacement = kin.compute_backward_displacement(position)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15ad238cc567f502", + "metadata": {}, + "outputs": [], + "source": [ + "# In this case, the backward displacement vector at the first timestep\n", + "# is defined as the zero vector, since there is no previous position.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28c3fbd36a0905e6", + "metadata": {}, + "outputs": [], + "source": [ + "# Adapting the code snippet from above, we can visually check that the\n", + "# backward displacement vector is indeed the reverse of the forward\n", + "# displacement vector.\n", + "\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot()\n", + "\n", + "sc = ax.scatter(\n", + " position.sel(individuals=mouse_name, space=\"x\"),\n", + " position.sel(individuals=mouse_name, space=\"y\"),\n", + " s=15,\n", + " c=position.time,\n", + " cmap=\"viridis\",\n", + ")\n", + "\n", + "ax.quiver(\n", + " position.sel(individuals=mouse_name, space=\"x\"),\n", + " position.sel(individuals=mouse_name, space=\"y\"),\n", + " backward_displacement.sel(individuals=mouse_name, space=\"x\"),\n", + " backward_displacement.sel(individuals=mouse_name, space=\"y\"),\n", + " angles=\"xy\",\n", + " scale=1,\n", + " scale_units=\"xy\",\n", + " headwidth=7,\n", + " headlength=9,\n", + " headaxislength=9,\n", + ")\n", + "\n", + "ax.set_xlim(480, 600)\n", + "ax.set_ylim(980, 1080)\n", + "ax.set_xlabel(\"x (pixels)\")\n", + "ax.set_ylabel(\"y (pixels)\")\n", + "ax.set_title(f\"Zoomed in backward trajectory of {mouse_name}\")\n", + "ax.invert_yaxis()\n", + "sc.set_clim(8.8, 9.2)\n", + "fig.colorbar(sc, ax=ax, label=\"time (s)\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "61ed02cba59293db", + "metadata": {}, + "outputs": [], + "source": [ + "# Compute path length\n", + "# --------------------\n", + "# We can compute the distance travelled by the\n", + "# mouse as the sum of the lengths of all\n", + "# displacement vectors along its trajectory.\n", + "# Both backward and forward displacement vectors\n", + "# should give the same result:\n", + "\n", + "# length of each forward displacement vector\n", + "forward_displacement_lengths = compute_norm(\n", + " forward_displacement.sel(individuals=mouse_name)\n", + ")\n", + "\n", + "# length of each backward displacement vector\n", + "backward_displacement_lengths = compute_norm(\n", + " backward_displacement.sel(individuals=mouse_name)\n", + ")\n", + "\n", + "# check their lengths are the same\n", + "np.testing.assert_almost_equal(\n", + " forward_displacement_lengths.values[:-1], # exclude last timestep\n", + " backward_displacement_lengths.values[1:], # exclude first timestep\n", + ")\n", + "\n", + "# sum the lengths of all displacement vectors (in pixels)\n", + "total_displacement_fwd = forward_displacement_lengths.sum(dim=\"time\").values[0]\n", + "total_displacement_bwd = backward_displacement_lengths.sum(dim=\"time\").values[\n", + " 0\n", + "]\n", + "\n", + "print(\n", + " f\"The mouse {mouse_name}'s path length is {total_displacement_fwd:.2f} \"\n", + " \"pixels long (using forward displacement)\"\n", + ")\n", + "print(\n", + " f\"The mouse {mouse_name}'s path length is {total_displacement_bwd:.2f} \"\n", + " \"pixels long (using backward displacement)\"\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6065379c3a4334a", + "metadata": {}, + "outputs": [], + "source": [ + "# We provide a convenience function\n", + "# :func:`movement.kinematics.compute_path_length`\n", + "# to compute the path length for all individuals and keypoints in a position\n", + "# data array. We can verify that using this function gives the same result\n", + "# as before for the ``AEON3B_TP2`` mouse:\n", + "\n", + "path_lengths = kin.compute_path_length(ds.position)\n", + "\n", + "for mouse_name in path_lengths.individuals.values:\n", + " print(\n", + " f\"Path length for {mouse_name}: \"\n", + " f\"{path_lengths.sel(individuals=mouse_name).values[0]:.2f} pixels\"\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "245f3310f4556cd5", + "metadata": {}, + "outputs": [], + "source": [ + "# Compute velocity\n", + "# ----------------\n", + "# We can also compute the velocity vectors for all individuals in our data\n", + "# array:\n", + "velocity = kin.compute_velocity(position)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "96c6aefd9862261", + "metadata": {}, + "outputs": [], + "source": [ + "# The :func:`movement.kinematics.compute_velocity`\n", + "# function will return a data array equivalent to\n", + "# the ``position`` one, but holding velocity data along the ``space`` axis,\n", + "# rather than position data. Notice how ``xarray`` nicely deals with the\n", + "# different individuals and spatial dimensions for us! ✨\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fb41a97a7d54d858", + "metadata": {}, + "outputs": [], + "source": [ + "# We can plot the components of the velocity vector against time\n", + "# using ``xarray``'s built-in plotting methods. We use\n", + "# :meth:`xarray.DataArray.squeeze` to\n", + "# remove the dimension of length 1 from the data (the ``keypoints`` dimension).\n", + "\n", + "velocity.squeeze().plot.line(x=\"time\", row=\"individuals\", aspect=2, size=2.5)\n", + "plt.gcf().show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "baedbc6a0ab6b7e", + "metadata": {}, + "outputs": [], + "source": [ + "# The components of the velocity vector seem noisier than the components of\n", + "# the position vector.\n", + "# This is expected, since we are estimating the velocity using differences in\n", + "# position (which is somewhat noisy), over small stepsizes.\n", + "# More specifically, we use :func:`numpy.gradient` internally, which\n", + "# uses second order central differences.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "408e47b1e721a427", + "metadata": {}, + "outputs": [], + "source": [ + "# We can also visualise the speed, as the magnitude (norm)\n", + "# of the velocity vector:\n", + "fig, axes = plt.subplots(3, 1, sharex=True, sharey=True)\n", + "for mouse_name, ax in zip(velocity.individuals.values, axes, strict=False):\n", + " # compute the magnitude of the velocity vector for one mouse\n", + " speed_one_mouse = compute_norm(velocity.sel(individuals=mouse_name))\n", + " # plot speed against time\n", + " ax.plot(speed_one_mouse)\n", + " ax.set_title(mouse_name)\n", + " ax.set_xlabel(\"time (s)\")\n", + " ax.set_ylabel(\"speed (px/s)\")\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "92173d472b4a7988", + "metadata": {}, + "outputs": [], + "source": [ + "# To visualise the direction of the velocity vector at each timestep, we can\n", + "# again use a quiver plot:\n", + "mouse_name = \"AEON3B_TP2\"\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot()\n", + "# plot trajectory (position data)\n", + "sc = ax.scatter(\n", + " position.sel(individuals=mouse_name, space=\"x\"),\n", + " position.sel(individuals=mouse_name, space=\"y\"),\n", + " s=15,\n", + " c=position.time,\n", + " cmap=\"viridis\",\n", + ")\n", + "# plot velocity vectors\n", + "ax.quiver(\n", + " position.sel(individuals=mouse_name, space=\"x\"),\n", + " position.sel(individuals=mouse_name, space=\"y\"),\n", + " velocity.sel(individuals=mouse_name, space=\"x\"),\n", + " velocity.sel(individuals=mouse_name, space=\"y\"),\n", + " angles=\"xy\",\n", + " scale=2,\n", + " scale_units=\"xy\",\n", + " color=\"r\",\n", + ")\n", + "ax.axis(\"equal\")\n", + "ax.set_xlabel(\"x (pixels)\")\n", + "ax.set_ylabel(\"y (pixels)\")\n", + "ax.set_title(f\"Velocity quiver plot for {mouse_name}\")\n", + "ax.invert_yaxis()\n", + "fig.colorbar(sc, ax=ax, label=\"time (s)\")\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41a64eea0e6ad0c0", + "metadata": {}, + "outputs": [], + "source": [ + "# Here we scaled the length of vectors to half of their actual value\n", + "# (``scale=2``) for easier visualisation.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eab3e8cae5abbd07", + "metadata": {}, + "outputs": [], + "source": [ + "# Compute acceleration\n", + "# ---------------------\n", + "# Let's now compute the acceleration for all individuals in our data\n", + "# array:\n", + "accel = kin.compute_acceleration(position)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f401fc6ac0118fe2", + "metadata": {}, + "outputs": [], + "source": [ + "# and plot of the components of the acceleration vector ``ax``, ``ay`` per\n", + "# individual:\n", + "fig, axes = plt.subplots(3, 1, sharex=True, sharey=True)\n", + "for mouse_name, ax in zip(accel.individuals.values, axes, strict=False):\n", + " # plot x-component of acceleration vector\n", + " ax.plot(\n", + " accel.sel(individuals=mouse_name, space=[\"x\"]).squeeze(),\n", + " label=\"ax\",\n", + " )\n", + " # plot y-component of acceleration vector\n", + " ax.plot(\n", + " accel.sel(individuals=mouse_name, space=[\"y\"]).squeeze(),\n", + " label=\"ay\",\n", + " )\n", + " ax.set_title(mouse_name)\n", + " ax.set_xlabel(\"time (s)\")\n", + " ax.set_ylabel(\"acceleration (px/s**2)\")\n", + " ax.legend(loc=\"center right\", bbox_to_anchor=(1.07, 1.07))\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2369424047919858", + "metadata": {}, + "outputs": [], + "source": [ + "# We can also compute and visualise the magnitude (norm) of the\n", + "# acceleration vector for each individual:\n", + "fig, axes = plt.subplots(3, 1, sharex=True, sharey=True)\n", + "for mouse_name, ax in zip(accel.individuals.values, axes, strict=False):\n", + " # compute magnitude of the acceleration vector for one mouse\n", + " accel_one_mouse = compute_norm(accel.sel(individuals=mouse_name))\n", + "\n", + " # plot acceleration against time\n", + " ax.plot(accel_one_mouse)\n", + " ax.set_title(mouse_name)\n", + " ax.set_xlabel(\"time (s)\")\n", + " ax.set_ylabel(\"accel (px/s**2)\")\n", + "fig.tight_layout()" + ] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/compute_kinematics.py b/examples/compute_kinematics.py index 1c7e08e0c..c71cff68f 100644 --- a/examples/compute_kinematics.py +++ b/examples/compute_kinematics.py @@ -419,7 +419,7 @@ ) ax.set_title(mouse_name) ax.set_xlabel("time (s)") - ax.set_ylabel("speed (px/s**2)") + ax.set_ylabel("accel (px/s**2)") ax.legend(loc="center right", bbox_to_anchor=(1.07, 1.07)) fig.tight_layout() diff --git a/movement/validators/files.py b/movement/validators/files.py index e75baf8ac..5d3c91f12 100644 --- a/movement/validators/files.py +++ b/movement/validators/files.py @@ -113,9 +113,7 @@ def _file_is_readable(value: Path) -> None: def _file_is_writable(value: Path) -> None: - """Ensure the file does not exist and parent directory is writable.""" - if value.exists(): - raise logger.error(FileExistsError(f"File {value} already exists.")) + """Ensure parent directory is writable.""" if not os.access(value.parent, os.W_OK): raise logger.error( PermissionError(