diff --git a/.buildinfo b/.buildinfo index e9b701c2..fd9192ac 100644 --- a/.buildinfo +++ b/.buildinfo @@ -1,4 +1,4 @@ # Sphinx build info version 1 # This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done. -config: 1f7f213a712a67423f3eefd33baed626 +config: a6ac25651f02e78abf3a4ac84ba4923e tags: 645f666f9bcd5a90fca523b33c5a78b7 diff --git a/.doctrees/api/movement.io.load_poses.from_dlc_df.doctree b/.doctrees/api/movement.io.load_poses.from_dlc_df.doctree deleted file mode 100644 index 40a5cc8f..00000000 Binary files a/.doctrees/api/movement.io.load_poses.from_dlc_df.doctree and /dev/null differ diff --git a/.doctrees/api/movement.io.load_poses.from_dlc_file.doctree b/.doctrees/api/movement.io.load_poses.from_dlc_file.doctree index cb631d40..14ef626b 100644 Binary files a/.doctrees/api/movement.io.load_poses.from_dlc_file.doctree and b/.doctrees/api/movement.io.load_poses.from_dlc_file.doctree differ diff --git a/.doctrees/api/movement.io.load_poses.from_dlc_style_df.doctree b/.doctrees/api/movement.io.load_poses.from_dlc_style_df.doctree new file mode 100644 index 00000000..188c5cf0 Binary files /dev/null and b/.doctrees/api/movement.io.load_poses.from_dlc_style_df.doctree differ diff --git a/.doctrees/api/movement.io.load_poses.from_file.doctree b/.doctrees/api/movement.io.load_poses.from_file.doctree index 6ff6af21..eaa91318 100644 Binary files a/.doctrees/api/movement.io.load_poses.from_file.doctree and b/.doctrees/api/movement.io.load_poses.from_file.doctree differ diff --git a/.doctrees/api/movement.io.load_poses.from_lp_file.doctree b/.doctrees/api/movement.io.load_poses.from_lp_file.doctree index 7a1dd3d7..a9e22daa 100644 Binary files a/.doctrees/api/movement.io.load_poses.from_lp_file.doctree and b/.doctrees/api/movement.io.load_poses.from_lp_file.doctree differ diff --git a/.doctrees/api/movement.io.load_poses.from_numpy.doctree b/.doctrees/api/movement.io.load_poses.from_numpy.doctree new file mode 100644 index 00000000..9e4996b6 Binary files /dev/null and b/.doctrees/api/movement.io.load_poses.from_numpy.doctree differ diff --git a/.doctrees/api/movement.io.load_poses.from_sleap_file.doctree b/.doctrees/api/movement.io.load_poses.from_sleap_file.doctree index 49c36db4..1aaf85b0 100644 Binary files a/.doctrees/api/movement.io.load_poses.from_sleap_file.doctree and b/.doctrees/api/movement.io.load_poses.from_sleap_file.doctree differ diff --git a/.doctrees/api/movement.io.save_poses.to_dlc_df.doctree b/.doctrees/api/movement.io.save_poses.to_dlc_df.doctree deleted file mode 100644 index 705909bf..00000000 Binary files a/.doctrees/api/movement.io.save_poses.to_dlc_df.doctree and /dev/null differ diff --git a/.doctrees/api/movement.io.save_poses.to_dlc_file.doctree b/.doctrees/api/movement.io.save_poses.to_dlc_file.doctree index 48bfcb8d..1d555985 100644 Binary files a/.doctrees/api/movement.io.save_poses.to_dlc_file.doctree and b/.doctrees/api/movement.io.save_poses.to_dlc_file.doctree differ diff --git a/.doctrees/api/movement.io.save_poses.to_dlc_style_df.doctree b/.doctrees/api/movement.io.save_poses.to_dlc_style_df.doctree new file mode 100644 index 00000000..8ea93860 Binary files /dev/null and b/.doctrees/api/movement.io.save_poses.to_dlc_style_df.doctree differ diff --git a/.doctrees/api/movement.io.save_poses.to_lp_file.doctree b/.doctrees/api/movement.io.save_poses.to_lp_file.doctree index fe5d50dd..7559ba1a 100644 Binary files a/.doctrees/api/movement.io.save_poses.to_lp_file.doctree and b/.doctrees/api/movement.io.save_poses.to_lp_file.doctree differ diff --git a/.doctrees/api/movement.io.save_poses.to_sleap_analysis_file.doctree b/.doctrees/api/movement.io.save_poses.to_sleap_analysis_file.doctree index 02b6c944..7220ec30 100644 Binary files a/.doctrees/api/movement.io.save_poses.to_sleap_analysis_file.doctree and b/.doctrees/api/movement.io.save_poses.to_sleap_analysis_file.doctree differ diff --git a/.doctrees/api/movement.io.validators.ValidDeepLabCutCSV.doctree b/.doctrees/api/movement.io.validators.ValidDeepLabCutCSV.doctree new file mode 100644 index 00000000..9c77cdb5 Binary files /dev/null and b/.doctrees/api/movement.io.validators.ValidDeepLabCutCSV.doctree differ diff --git a/.doctrees/api/movement.io.validators.ValidFile.doctree b/.doctrees/api/movement.io.validators.ValidFile.doctree index 96fad33f..ba3c0773 100644 Binary files a/.doctrees/api/movement.io.validators.ValidFile.doctree and b/.doctrees/api/movement.io.validators.ValidFile.doctree differ diff --git a/.doctrees/api/movement.io.validators.ValidPosesCSV.doctree b/.doctrees/api/movement.io.validators.ValidPosesCSV.doctree deleted file mode 100644 index 0b9fff6f..00000000 Binary files a/.doctrees/api/movement.io.validators.ValidPosesCSV.doctree and /dev/null differ diff --git a/.doctrees/api/movement.io.validators.ValidPosesDataset.doctree b/.doctrees/api/movement.io.validators.ValidPosesDataset.doctree index 1647e8a0..f1d0d310 100644 Binary files a/.doctrees/api/movement.io.validators.ValidPosesDataset.doctree and b/.doctrees/api/movement.io.validators.ValidPosesDataset.doctree differ diff --git a/.doctrees/api_index.doctree b/.doctrees/api_index.doctree index 04475529..73ac3be2 100644 Binary files a/.doctrees/api_index.doctree and b/.doctrees/api_index.doctree differ diff --git a/.doctrees/community/contributing.doctree b/.doctrees/community/contributing.doctree index f54b778a..c1f616fe 100644 Binary files a/.doctrees/community/contributing.doctree and b/.doctrees/community/contributing.doctree differ diff --git a/.doctrees/environment.pickle b/.doctrees/environment.pickle index 2d4d6352..e79c3521 100644 Binary files a/.doctrees/environment.pickle and b/.doctrees/environment.pickle differ diff --git a/.doctrees/examples/compute_kinematics.doctree b/.doctrees/examples/compute_kinematics.doctree index bc8c0e25..7cc16de7 100644 Binary files a/.doctrees/examples/compute_kinematics.doctree and b/.doctrees/examples/compute_kinematics.doctree differ diff --git a/.doctrees/examples/compute_polar_coordinates.doctree b/.doctrees/examples/compute_polar_coordinates.doctree new file mode 100644 index 00000000..dad3ad19 Binary files /dev/null and b/.doctrees/examples/compute_polar_coordinates.doctree differ diff --git a/.doctrees/examples/filter_and_interpolate.doctree b/.doctrees/examples/filter_and_interpolate.doctree index 0c6d6e2e..399b412f 100644 Binary files a/.doctrees/examples/filter_and_interpolate.doctree and b/.doctrees/examples/filter_and_interpolate.doctree differ diff --git a/.doctrees/examples/index.doctree b/.doctrees/examples/index.doctree index adc01c07..078235b5 100644 Binary files a/.doctrees/examples/index.doctree and b/.doctrees/examples/index.doctree differ diff --git a/.doctrees/examples/load_and_explore_poses.doctree b/.doctrees/examples/load_and_explore_poses.doctree index 87de2e74..1ffaf1b1 100644 Binary files a/.doctrees/examples/load_and_explore_poses.doctree and b/.doctrees/examples/load_and_explore_poses.doctree differ diff --git a/.doctrees/examples/sg_execution_times.doctree b/.doctrees/examples/sg_execution_times.doctree index bd6444f5..520b2bab 100644 Binary files a/.doctrees/examples/sg_execution_times.doctree and b/.doctrees/examples/sg_execution_times.doctree differ diff --git a/.doctrees/examples/smooth.doctree b/.doctrees/examples/smooth.doctree new file mode 100644 index 00000000..3a3956b8 Binary files /dev/null and b/.doctrees/examples/smooth.doctree differ diff --git a/.doctrees/sg_execution_times.doctree b/.doctrees/sg_execution_times.doctree index 07b645bd..2e91e951 100644 Binary files a/.doctrees/sg_execution_times.doctree and b/.doctrees/sg_execution_times.doctree differ diff --git a/_downloads/07f4b1b69497c93530bbe772d7eb535a/smooth.py b/_downloads/07f4b1b69497c93530bbe772d7eb535a/smooth.py new file mode 100644 index 00000000..3bb44162 --- /dev/null +++ b/_downloads/07f4b1b69497c93530bbe772d7eb535a/smooth.py @@ -0,0 +1,314 @@ +"""Smooth pose tracks +===================== + +Smooth pose tracks using the median and Savitzky-Golay filters. +""" + +# %% +# Imports +# ------- + +from matplotlib import pyplot as plt +from scipy.signal import welch + +from movement import sample_data +from movement.filtering import ( + interpolate_over_time, + median_filter, + savgol_filter, +) + +# %% +# Load a sample dataset +# --------------------- +# Let's load a sample dataset and print it to inspect its contents. +# Note that if you are running this notebook interactively, you can simply +# type the variable name (here ``ds_wasp``) in a cell to get an interactive +# display of the dataset's contents. + +ds_wasp = sample_data.fetch_dataset("DLC_single-wasp.predictions.h5") +print(ds_wasp) + +# %% +# We see that the dataset contains a single individual (a wasp) with two +# keypoints tracked in 2D space. The video was recorded at 40 fps and lasts for +# ~27 seconds. + +# %% +# Define a plotting function +# -------------------------- +# Let's define a plotting function to help us visualise the effects smoothing +# both in the time and frequency domains. +# The function takes as inputs two datasets containing raw and smooth data +# respectively, and plots the position time series and power spectral density +# (PSD) for a given individual and keypoint. The function also allows you to +# specify the spatial coordinate (``x`` or ``y``) and a time range to focus on. + + +def plot_raw_and_smooth_timeseries_and_psd( + ds_raw, + ds_smooth, + individual="individual_0", + keypoint="stinger", + space="x", + time_range=None, +): + # If no time range is specified, plot the entire time series + if time_range is None: + time_range = slice(0, ds_raw.time[-1]) + + selection = { + "time": time_range, + "individuals": individual, + "keypoints": keypoint, + "space": space, + } + + fig, ax = plt.subplots(2, 1, figsize=(10, 6)) + + for ds, color, label in zip( + [ds_raw, ds_smooth], ["k", "r"], ["raw", "smooth"] + ): + # plot position time series + pos = ds.position.sel(**selection) + ax[0].plot( + pos.time, + pos, + color=color, + lw=2, + alpha=0.7, + label=f"{label} {space}", + ) + + # generate interpolated dataset to avoid NaNs in the PSD calculation + ds_interp = interpolate_over_time(ds, max_gap=None, print_report=False) + pos_interp = ds_interp.position.sel(**selection) + # compute and plot the PSD + freq, psd = welch(pos_interp, fs=ds.fps, nperseg=256) + ax[1].semilogy( + freq, + psd, + color=color, + lw=2, + alpha=0.7, + label=f"{label} {space}", + ) + + ax[0].set_ylabel(f"{space} position (px)") + ax[0].set_xlabel("Time (s)") + ax[0].set_title("Time Domain") + ax[0].legend() + + ax[1].set_ylabel("PSD (px$^2$/Hz)") + ax[1].set_xlabel("Frequency (Hz)") + ax[1].set_title("Frequency Domain") + ax[1].legend() + + plt.tight_layout() + fig.show() + + +# %% +# Smoothing with a median filter +# ------------------------------ +# Here we use the :py:func:`movement.filtering.median_filter` function to +# apply a rolling window median filter to the wasp dataset. +# The ``window_length`` parameter is defined in seconds (according to the +# ``time_unit`` dataset attribute). + +ds_wasp_medfilt = median_filter(ds_wasp, window_length=0.1) + +# %% +# We see from the printed report that the dataset has no missing values +# neither before nor after smoothing. Let's visualise the effects of the +# median filter in the time and frequency domains. + +plot_raw_and_smooth_timeseries_and_psd( + ds_wasp, ds_wasp_medfilt, keypoint="stinger" +) + +# %% +# We see that the median filter has removed the "spikes" present around the +# 14 second mark in the raw data. However, it has not dealt the big shift +# occurring during the final second. In the frequency domain, we can see that +# the filter has reduced the power in the high-frequency components, without +# affecting the low frequency components. +# +# This illustrates what the median filter is good at: removing brief "spikes" +# (e.g. a keypoint abruptly jumping to a different location for a frame or two) +# and high-frequency "jitter" (often present due to pose estimation +# working on a per-frame basis). + +# %% +# Choosing parameters for the median filter +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# You can control the behaviour of :py:func:`movement.filtering.median_filter` +# via two parameters: ``window_length`` and ``min_periods``. +# To better understand the effect of these parameters, let's use a +# dataset that contains missing values. + +ds_mouse = sample_data.fetch_dataset("SLEAP_single-mouse_EPM.analysis.h5") +print(ds_mouse) + +# %% +# The dataset contains a single mouse with six keypoints tracked in +# 2D space. The video was recorded at 30 fps and lasts for ~616 seconds. We can +# see that there are some missing values, indicated as "nan" in the +# printed dataset. Let's apply the median filter to this dataset, with +# the ``window_length`` set to 0.1 seconds. + +ds_mouse_medfilt = median_filter(ds_mouse, window_length=0.1) + +# %% +# The report informs us that the raw data contains NaN values, most of which +# occur at the ``snout`` and ``tail_end`` keypoints. After filtering, the +# number of NaNs has increased. This is because the default behaviour of the +# median filter is to propagate NaN values, i.e. if any value in the rolling +# window is NaN, the output will also be NaN. +# +# To modify this behaviour, you can set the value of the ``min_periods`` +# parameter to an integer value. This parameter determines the minimum number +# of non-NaN values required in the window for the output to be non-NaN. +# For example, setting ``min_periods=2`` means that two non-NaN values in the +# window are sufficient for the median to be calculated. Let's try this. + +ds_mouse_medfilt = median_filter(ds_mouse, window_length=0.1, min_periods=2) + +# %% +# We see that this time the number of NaN values has decreased +# across all keypoints. +# Let's visualise the effects of the median filter in the time and frequency +# domains. Here we focus on the first 80 seconds for the ``snout`` keypoint. +# You can adjust the ``keypoint`` and ``time_range`` arguments to explore other +# parts of the data. + +plot_raw_and_smooth_timeseries_and_psd( + ds_mouse, ds_mouse_medfilt, keypoint="snout", time_range=slice(0, 80) +) + +# %% +# The smoothing once again reduces the power of high-frequency components, but +# the resulting time series stays quite close to the raw data. +# +# What happens if we increase the ``window_length`` to 2 seconds? + +ds_mouse_medfilt = median_filter(ds_mouse, window_length=2, min_periods=2) + +# %% +# The number of NaN values has decreased even further. That's because the +# chance of finding at least 2 valid values within a 2 second window is +# quite high. Let's plot the results for the same keypoint and time range +# as before. + +plot_raw_and_smooth_timeseries_and_psd( + ds_mouse, ds_mouse_medfilt, keypoint="snout", time_range=slice(0, 80) +) +# %% +# We see that the filtered time series is much smoother and it has even +# "bridged" over some small gaps. That said, it often deviates from the raw +# data, in ways that may not be desirable, depending on the application. +# That means that our choice of ``window_length`` may be too large. +# In general, you should choose a ``window_length`` that is small enough to +# preserve the original data structure, but large enough to remove +# "spikes" and high-frequency noise. Always inspect the results to ensure +# that the filter is not removing important features. + +# %% +# Smoothing with a Savitzky-Golay filter +# -------------------------------------- +# Here we use the :py:func:`movement.filtering.savgol_filter` function, +# which is a wrapper around :py:func:`scipy.signal.savgol_filter`. +# The Savitzky-Golay filter is a polynomial smoothing filter that can be +# applied to time series data on a rolling window basis. A polynomial of +# degree ``polyorder`` is fitted to the data in each window of length +# ``window_length``, and the value of the polynomial at the center of the +# window is used as the output value. +# +# Let's try it on the mouse dataset. + +ds_mouse_savgol = savgol_filter(ds_mouse, window_length=0.2, polyorder=2) + + +# %% +# We see that the number of NaN values has increased after filtering. This is +# for the same reason as with the median filter (in its default mode), i.e. +# if there is at least one NaN value in the window, the output will be NaN. +# Unlike the median filter, the Savitzky-Golay filter does not provide a +# ``min_periods`` parameter to control this behaviour. Let's visualise the +# effects in the time and frequency domains. + +plot_raw_and_smooth_timeseries_and_psd( + ds_mouse, ds_mouse_savgol, keypoint="snout", time_range=slice(0, 80) +) +# %% +# Once again, the power of high-frequency components has been reduced, but more +# missing values have been introduced. + +# %% +# Now let's take a look at the wasp dataset. + +ds_wasp_savgol = savgol_filter(ds_wasp, window_length=0.2, polyorder=2) + +# %% +plot_raw_and_smooth_timeseries_and_psd( + ds_wasp, + ds_wasp_savgol, + keypoint="stinger", +) +# %% +# This example shows two important limitations of the Savitzky-Golay filter. +# First, the filter can introduce artefacts around sharp boundaries. For +# example, focus on what happens around the sudden drop in position +# during the final second. Second, the PSD appears to have large periodic +# drops at certain frequencies. Both of these effects vary with the +# choice of ``window_length`` and ``polyorder``. You can read more about these +# and other limitations of the Savitzky-Golay filter in +# `this paper `_. + + +# %% +# Combining multiple smoothing filters +# ------------------------------------ +# You can also combine multiple smoothing filters by applying them +# sequentially. For example, we can first apply the median filter with a small +# ``window_length`` to remove "spikes" and then apply the Savitzky-Golay filter +# with a larger ``window_length`` to further smooth the data. +# Between the two filters, we can interpolate over small gaps to avoid the +# excessive proliferation of NaN values. Let's try this on the mouse dataset. +# First, let's apply the median filter. + +ds_mouse_medfilt = median_filter(ds_mouse, window_length=0.1, min_periods=2) + +# %% +# Next, let's linearly interpolate over gaps smaller than 1 second. + +ds_mouse_medfilt_interp = interpolate_over_time(ds_mouse_medfilt, max_gap=1) + +# %% +# Finally, let's apply the Savitzky-Golay filter. + +ds_mouse_medfilt_interp_savgol = savgol_filter( + ds_mouse_medfilt_interp, window_length=0.4, polyorder=2 +) + +# %% +# A record of all applied operations is stored in the dataset's ``log`` +# attribute. Let's inspect it to summarise what we've done. + +for entry in ds_mouse_medfilt_interp_savgol.log: + print(entry) + +# %% +# Now let's visualise the difference between the raw data and the final +# smoothed result. + +plot_raw_and_smooth_timeseries_and_psd( + ds_mouse, + ds_mouse_medfilt_interp_savgol, + keypoint="snout", + time_range=slice(0, 80), +) + +# %% +# Feel free to play around with the parameters of the applied filters and to +# also look at other keypoints and time ranges. diff --git a/_downloads/085258163d1f42f74bf76da18a848b11/compute_kinematics.ipynb b/_downloads/085258163d1f42f74bf76da18a848b11/compute_kinematics.ipynb index 250e4a8c..600b53fa 100644 --- a/_downloads/085258163d1f42f74bf76da18a848b11/compute_kinematics.ipynb +++ b/_downloads/085258163d1f42f74bf76da18a848b11/compute_kinematics.ipynb @@ -4,7 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Compute and visualise kinematics.\n\nCompute displacement, velocity and acceleration data on an example dataset and\nvisualise the results.\n" + "# Compute and visualise kinematics.\n\nCompute displacement, velocity and acceleration, and\nvisualise the results.\n" ] }, { diff --git a/_downloads/122338c6db2328ed9eeb5e9961344704/filter_and_interpolate.ipynb b/_downloads/122338c6db2328ed9eeb5e9961344704/filter_and_interpolate.ipynb index 5faa78fc..bd2d7e2e 100644 --- a/_downloads/122338c6db2328ed9eeb5e9961344704/filter_and_interpolate.ipynb +++ b/_downloads/122338c6db2328ed9eeb5e9961344704/filter_and_interpolate.ipynb @@ -4,7 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Filtering and interpolation\n\nFilter out points with low confidence scores and interpolate over\nmissing values.\n" + "# Drop outliers and interpolate\n\nFilter out points with low confidence scores and interpolate over\nmissing values.\n" ] }, { diff --git a/_downloads/844e74b700fac354e52a59dbe7329a84/smooth.ipynb b/_downloads/844e74b700fac354e52a59dbe7329a84/smooth.ipynb new file mode 100644 index 00000000..c9d60c82 --- /dev/null +++ b/_downloads/844e74b700fac354e52a59dbe7329a84/smooth.ipynb @@ -0,0 +1,427 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Smooth pose tracks\n\nSmooth pose tracks using the median and Savitzky-Golay filters.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Imports\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt\nfrom scipy.signal import welch\n\nfrom movement import sample_data\nfrom movement.filtering import (\n interpolate_over_time,\n median_filter,\n savgol_filter,\n)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Load a sample dataset\nLet's load a sample dataset and print it to inspect its contents.\nNote that if you are running this notebook interactively, you can simply\ntype the variable name (here ``ds_wasp``) in a cell to get an interactive\ndisplay of the dataset's contents.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "ds_wasp = sample_data.fetch_dataset(\"DLC_single-wasp.predictions.h5\")\nprint(ds_wasp)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see that the dataset contains a single individual (a wasp) with two\nkeypoints tracked in 2D space. The video was recorded at 40 fps and lasts for\n~27 seconds.\n\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Define a plotting function\nLet's define a plotting function to help us visualise the effects smoothing\nboth in the time and frequency domains.\nThe function takes as inputs two datasets containing raw and smooth data\nrespectively, and plots the position time series and power spectral density\n(PSD) for a given individual and keypoint. The function also allows you to\nspecify the spatial coordinate (``x`` or ``y``) and a time range to focus on.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def plot_raw_and_smooth_timeseries_and_psd(\n ds_raw,\n ds_smooth,\n individual=\"individual_0\",\n keypoint=\"stinger\",\n space=\"x\",\n time_range=None,\n):\n # If no time range is specified, plot the entire time series\n if time_range is None:\n time_range = slice(0, ds_raw.time[-1])\n\n selection = {\n \"time\": time_range,\n \"individuals\": individual,\n \"keypoints\": keypoint,\n \"space\": space,\n }\n\n fig, ax = plt.subplots(2, 1, figsize=(10, 6))\n\n for ds, color, label in zip(\n [ds_raw, ds_smooth], [\"k\", \"r\"], [\"raw\", \"smooth\"]\n ):\n # plot position time series\n pos = ds.position.sel(**selection)\n ax[0].plot(\n pos.time,\n pos,\n color=color,\n lw=2,\n alpha=0.7,\n label=f\"{label} {space}\",\n )\n\n # generate interpolated dataset to avoid NaNs in the PSD calculation\n ds_interp = interpolate_over_time(ds, max_gap=None, print_report=False)\n pos_interp = ds_interp.position.sel(**selection)\n # compute and plot the PSD\n freq, psd = welch(pos_interp, fs=ds.fps, nperseg=256)\n ax[1].semilogy(\n freq,\n psd,\n color=color,\n lw=2,\n alpha=0.7,\n label=f\"{label} {space}\",\n )\n\n ax[0].set_ylabel(f\"{space} position (px)\")\n ax[0].set_xlabel(\"Time (s)\")\n ax[0].set_title(\"Time Domain\")\n ax[0].legend()\n\n ax[1].set_ylabel(\"PSD (px$^2$/Hz)\")\n ax[1].set_xlabel(\"Frequency (Hz)\")\n ax[1].set_title(\"Frequency Domain\")\n ax[1].legend()\n\n plt.tight_layout()\n fig.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Smoothing with a median filter\nHere we use the :py:func:`movement.filtering.median_filter` function to\napply a rolling window median filter to the wasp dataset.\nThe ``window_length`` parameter is defined in seconds (according to the\n``time_unit`` dataset attribute).\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "ds_wasp_medfilt = median_filter(ds_wasp, window_length=0.1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see from the printed report that the dataset has no missing values\nneither before nor after smoothing. Let's visualise the effects of the\nmedian filter in the time and frequency domains.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "plot_raw_and_smooth_timeseries_and_psd(\n ds_wasp, ds_wasp_medfilt, keypoint=\"stinger\"\n)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see that the median filter has removed the \"spikes\" present around the\n14 second mark in the raw data. However, it has not dealt the big shift\noccurring during the final second. In the frequency domain, we can see that\nthe filter has reduced the power in the high-frequency components, without\naffecting the low frequency components.\n\nThis illustrates what the median filter is good at: removing brief \"spikes\"\n(e.g. a keypoint abruptly jumping to a different location for a frame or two)\nand high-frequency \"jitter\" (often present due to pose estimation\nworking on a per-frame basis).\n\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Choosing parameters for the median filter\nYou can control the behaviour of :py:func:`movement.filtering.median_filter`\nvia two parameters: ``window_length`` and ``min_periods``.\nTo better understand the effect of these parameters, let's use a\ndataset that contains missing values.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "ds_mouse = sample_data.fetch_dataset(\"SLEAP_single-mouse_EPM.analysis.h5\")\nprint(ds_mouse)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The dataset contains a single mouse with six keypoints tracked in\n2D space. The video was recorded at 30 fps and lasts for ~616 seconds. We can\nsee that there are some missing values, indicated as \"nan\" in the\nprinted dataset. Let's apply the median filter to this dataset, with\nthe ``window_length`` set to 0.1 seconds.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "ds_mouse_medfilt = median_filter(ds_mouse, window_length=0.1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The report informs us that the raw data contains NaN values, most of which\noccur at the ``snout`` and ``tail_end`` keypoints. After filtering, the\nnumber of NaNs has increased. This is because the default behaviour of the\nmedian filter is to propagate NaN values, i.e. if any value in the rolling\nwindow is NaN, the output will also be NaN.\n\nTo modify this behaviour, you can set the value of the ``min_periods``\nparameter to an integer value. This parameter determines the minimum number\nof non-NaN values required in the window for the output to be non-NaN.\nFor example, setting ``min_periods=2`` means that two non-NaN values in the\nwindow are sufficient for the median to be calculated. Let's try this.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "ds_mouse_medfilt = median_filter(ds_mouse, window_length=0.1, min_periods=2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see that this time the number of NaN values has decreased\nacross all keypoints.\nLet's visualise the effects of the median filter in the time and frequency\ndomains. Here we focus on the first 80 seconds for the ``snout`` keypoint.\nYou can adjust the ``keypoint`` and ``time_range`` arguments to explore other\nparts of the data.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "plot_raw_and_smooth_timeseries_and_psd(\n ds_mouse, ds_mouse_medfilt, keypoint=\"snout\", time_range=slice(0, 80)\n)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The smoothing once again reduces the power of high-frequency components, but\nthe resulting time series stays quite close to the raw data.\n\nWhat happens if we increase the ``window_length`` to 2 seconds?\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "ds_mouse_medfilt = median_filter(ds_mouse, window_length=2, min_periods=2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The number of NaN values has decreased even further. That's because the\nchance of finding at least 2 valid values within a 2 second window is\nquite high. Let's plot the results for the same keypoint and time range\nas before.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "plot_raw_and_smooth_timeseries_and_psd(\n ds_mouse, ds_mouse_medfilt, keypoint=\"snout\", time_range=slice(0, 80)\n)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see that the filtered time series is much smoother and it has even\n\"bridged\" over some small gaps. That said, it often deviates from the raw\ndata, in ways that may not be desirable, depending on the application.\nThat means that our choice of ``window_length`` may be too large.\nIn general, you should choose a ``window_length`` that is small enough to\npreserve the original data structure, but large enough to remove\n\"spikes\" and high-frequency noise. Always inspect the results to ensure\nthat the filter is not removing important features.\n\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Smoothing with a Savitzky-Golay filter\nHere we use the :py:func:`movement.filtering.savgol_filter` function,\nwhich is a wrapper around :py:func:`scipy.signal.savgol_filter`.\nThe Savitzky-Golay filter is a polynomial smoothing filter that can be\napplied to time series data on a rolling window basis. A polynomial of\ndegree ``polyorder`` is fitted to the data in each window of length\n``window_length``, and the value of the polynomial at the center of the\nwindow is used as the output value.\n\nLet's try it on the mouse dataset.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "ds_mouse_savgol = savgol_filter(ds_mouse, window_length=0.2, polyorder=2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see that the number of NaN values has increased after filtering. This is\nfor the same reason as with the median filter (in its default mode), i.e.\nif there is at least one NaN value in the window, the output will be NaN.\nUnlike the median filter, the Savitzky-Golay filter does not provide a\n``min_periods`` parameter to control this behaviour. Let's visualise the\neffects in the time and frequency domains.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "plot_raw_and_smooth_timeseries_and_psd(\n ds_mouse, ds_mouse_savgol, keypoint=\"snout\", time_range=slice(0, 80)\n)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Once again, the power of high-frequency components has been reduced, but more\nmissing values have been introduced.\n\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now let's take a look at the wasp dataset.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "ds_wasp_savgol = savgol_filter(ds_wasp, window_length=0.2, polyorder=2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "plot_raw_and_smooth_timeseries_and_psd(\n ds_wasp,\n ds_wasp_savgol,\n keypoint=\"stinger\",\n)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This example shows two important limitations of the Savitzky-Golay filter.\nFirst, the filter can introduce artefacts around sharp boundaries. For\nexample, focus on what happens around the sudden drop in position\nduring the final second. Second, the PSD appears to have large periodic\ndrops at certain frequencies. Both of these effects vary with the\nchoice of ``window_length`` and ``polyorder``. You can read more about these\nand other limitations of the Savitzky-Golay filter in\n[this paper](https://pubs.acs.org/doi/10.1021/acsmeasuresciau.1c00054).\n\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Combining multiple smoothing filters\nYou can also combine multiple smoothing filters by applying them\nsequentially. For example, we can first apply the median filter with a small\n``window_length`` to remove \"spikes\" and then apply the Savitzky-Golay filter\nwith a larger ``window_length`` to further smooth the data.\nBetween the two filters, we can interpolate over small gaps to avoid the\nexcessive proliferation of NaN values. Let's try this on the mouse dataset.\nFirst, let's apply the median filter.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "ds_mouse_medfilt = median_filter(ds_mouse, window_length=0.1, min_periods=2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, let's linearly interpolate over gaps smaller than 1 second.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "ds_mouse_medfilt_interp = interpolate_over_time(ds_mouse_medfilt, max_gap=1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, let's apply the Savitzky-Golay filter.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "ds_mouse_medfilt_interp_savgol = savgol_filter(\n ds_mouse_medfilt_interp, window_length=0.4, polyorder=2\n)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A record of all applied operations is stored in the dataset's ``log``\nattribute. Let's inspect it to summarise what we've done.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "for entry in ds_mouse_medfilt_interp_savgol.log:\n print(entry)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now let's visualise the difference between the raw data and the final\nsmoothed result.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "plot_raw_and_smooth_timeseries_and_psd(\n ds_mouse,\n ds_mouse_medfilt_interp_savgol,\n keypoint=\"snout\",\n time_range=slice(0, 80),\n)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Feel free to play around with the parameters of the applied filters and to\nalso look at other keypoints and time ranges.\n\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.9" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file diff --git a/_downloads/8616797cd8df925412b19674f364ffea/filter_and_interpolate.py b/_downloads/8616797cd8df925412b19674f364ffea/filter_and_interpolate.py index b3461dc2..d3e52ca5 100644 --- a/_downloads/8616797cd8df925412b19674f364ffea/filter_and_interpolate.py +++ b/_downloads/8616797cd8df925412b19674f364ffea/filter_and_interpolate.py @@ -1,5 +1,5 @@ -"""Filtering and interpolation -============================ +"""Drop outliers and interpolate +================================ Filter out points with low confidence scores and interpolate over missing values. diff --git a/_downloads/8c764d9f054b54b3d9a974597939bbdb/compute_polar_coordinates.py b/_downloads/8c764d9f054b54b3d9a974597939bbdb/compute_polar_coordinates.py new file mode 100644 index 00000000..bad25b79 --- /dev/null +++ b/_downloads/8c764d9f054b54b3d9a974597939bbdb/compute_polar_coordinates.py @@ -0,0 +1,404 @@ +"""Express 2D vectors in polar coordinates +============================================ + +Compute a vector representing head direction and express it in polar +coordinates. +""" +# %% +# Imports +# ------- + +# 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 + +from movement import sample_data +from movement.io import load_poses +from movement.utils.vector import cart2pol, pol2cart + +# %% +# Load sample dataset +# ------------------------ +# In this tutorial, we will use a sample dataset with a single individual +# (a mouse) and six keypoints. + +ds_path = sample_data.fetch_dataset_paths( + "SLEAP_single-mouse_EPM.analysis.h5" +)["poses"] +ds = load_poses.from_sleap_file(ds_path, fps=None) # force time_unit = frames + +print(ds) +print("-----------------------------") +print(f"Individuals: {ds.individuals.values}") +print(f"Keypoints: {ds.keypoints.values}") + + +# %% +# The loaded dataset ``ds`` contains two data variables:``position`` and +# ``confidence``. Both are stored as data arrays. In this tutorial, we will +# use only ``position``: +position = ds.position + + +# %% +# Compute head vector +# --------------------- +# To demonstrate how polar coordinates can be useful in behavioural analyses, +# we will compute the head vector of the mouse. +# +# We define it as the vector from the midpoint between the ears to the snout. + +# compute the midpoint between the ears +midpoint_ears = position.sel(keypoints=["left_ear", "right_ear"]).mean( + dim="keypoints" +) + +# compute the head vector +head_vector = position.sel(keypoints="snout") - midpoint_ears + +# drop the keypoints dimension +# (otherwise the `head_vector` data array retains a `snout` keypoint from the +# operation above) +head_vector = head_vector.drop_vars("keypoints") + +# %% +# Visualise the head trajectory +# -------------------------------------- +# We can plot the data to check that our computation of the head vector is +# correct. +# +# We can start by plotting the trajectory of the midpoint between the ears. We +# will refer to this as the head trajectory. + +fig, ax = plt.subplots(1, 1) +mouse_name = ds.individuals.values[0] + +sc = ax.scatter( + midpoint_ears.sel(individuals=mouse_name, space="x"), + midpoint_ears.sel(individuals=mouse_name, space="y"), + s=15, + c=midpoint_ears.time, + cmap="viridis", + marker="o", +) + +ax.axis("equal") +ax.set_xlabel("x (pixels)") +ax.set_ylabel("y (pixels)") +ax.invert_yaxis() +ax.set_title(f"Head trajectory ({mouse_name})") +fig.colorbar(sc, ax=ax, label=f"time ({ds.attrs['time_unit']})") +fig.show() + +# %% +# We can see that the majority of the head trajectory data is within a +# cruciform shape. This is because the dataset is of a mouse moving on an +# `Elevated Plus Maze `_. +# We can actually verify this is the case by overlaying the head +# trajectory on the sample frame of the dataset. + +# read sample frame +frame_path = sample_data.fetch_dataset_paths( + "SLEAP_single-mouse_EPM.analysis.h5" +)["frame"] +im = plt.imread(frame_path) + + +# plot sample frame +fig, ax = plt.subplots(1, 1) +ax.imshow(im) + +# plot head trajectory with semi-transparent markers +sc = ax.scatter( + midpoint_ears.sel(individuals=mouse_name, space="x"), + midpoint_ears.sel(individuals=mouse_name, space="y"), + s=15, + c=midpoint_ears.time, + cmap="viridis", + marker="o", + alpha=0.05, # transparency +) + +ax.axis("equal") +ax.set_xlabel("x (pixels)") +ax.set_ylabel("y (pixels)") +# No need to invert the y-axis now, since the image is plotted +# using a pixel coordinate system with origin on the top left of the image +ax.set_title(f"Head trajectory ({mouse_name})") + +fig.show() + +# %% +# The overlaid plot suggests the mouse spends most of its time in the +# covered arms of the maze. + +# %% +# Visualise the head vector +# --------------------------- +# To visually check our computation of the head vector, it is easier to select +# a subset of the data. We can focus on the trajectory of the head when the +# mouse is within a small rectangular area and time window. + +# area of interest +xmin, ymin = 600, 665 # pixels +x_delta, y_delta = 125, 100 # pixels + +# time window +time_window = range(1650, 1671) # frames + + +# %% +# For that subset of the data, we now plot the head vector. + +fig, ax = plt.subplots(1, 1) +mouse_name = ds.individuals.values[0] + +# plot midpoint between the ears, and color based on time +sc = ax.scatter( + midpoint_ears.sel(individuals=mouse_name, space="x", time=time_window), + midpoint_ears.sel(individuals=mouse_name, space="y", time=time_window), + s=50, + c=midpoint_ears.time[time_window], + cmap="viridis", + marker="*", +) + +# plot snout, and color based on time +sc = ax.scatter( + position.sel( + individuals=mouse_name, space="x", time=time_window, keypoints="snout" + ), + position.sel( + individuals=mouse_name, space="y", time=time_window, keypoints="snout" + ), + s=50, + c=position.time[time_window], + cmap="viridis", + marker="o", +) + +# plot the computed head vector +ax.quiver( + midpoint_ears.sel(individuals=mouse_name, space="x", time=time_window), + midpoint_ears.sel(individuals=mouse_name, space="y", time=time_window), + head_vector.sel(individuals=mouse_name, space="x", time=time_window), + head_vector.sel(individuals=mouse_name, space="y", time=time_window), + angles="xy", + scale=1, + scale_units="xy", + headwidth=7, + headlength=9, + headaxislength=9, + color="gray", +) + +ax.axis("equal") +ax.set_xlim(xmin, xmin + x_delta) +ax.set_ylim(ymin, ymin + y_delta) +ax.set_xlabel("x (pixels)") +ax.set_ylabel("y (pixels)") +ax.set_title(f"Zoomed in head vector ({mouse_name})") +ax.invert_yaxis() +fig.colorbar( + sc, + ax=ax, + label=f"time ({ds.attrs['time_unit']})", + ticks=list(time_window)[0::2], +) + +ax.legend( + [ + "midpoint_ears", + "snout", + "head_vector", + ], + loc="best", +) + +fig.show() + +# %% +# From the plot we can confirm the head vector goes from the midpoint between +# the ears to the snout, as we defined it. + + +# %% +# Express the head vector in polar coordinates +# ------------------------------------------------------------- +# A convenient way to inspect the orientation of a vector in 2D is by +# expressing it in polar coordinates. We can do this with the vector function +# ``cart2pol``: +head_vector_polar = cart2pol(head_vector) + +print(head_vector_polar) + +# %% +# Notice how the resulting array has a ``space_pol`` dimension with two +# coordinates: ``rho`` and ``phi``. These are the polar coordinates of the +# head vector. +# +# The coordinate ``rho`` is the norm (i.e., magnitude, length) of the vector. +# In our case, the distance from the midpoint between the ears to the snout. +# The coordinate ``phi`` is the orientation of the head vector relative to the +# positive x-axis, and ranges from -``pi`` to ``pi`` +# (following the `atan2 `_ convention). +# +# In our coordinate system, ``phi`` will be +# positive if the shortest path from the positive x-axis to the vector is +# clockwise. Conversely, ``phi`` will be negative if the shortest path from +# the positive x-axis to the vector is anti-clockwise. + +# %% +# Histogram of ``rho`` values +# ---------------------------- +# Since ``rho`` is the distance between the ears' midpoint and the snout, +# we would expect ``rho`` to be approximately constant in this data. We can +# check this by plotting a histogram of its values across the whole clip. + +fig, ax = plt.subplots(1, 1) + +# plot histogram using xarray's built-in histogram function +rho_data = head_vector_polar.sel(individuals=mouse_name, space_pol="rho") +rho_data.plot.hist(bins=50, ax=ax, edgecolor="lightgray", linewidth=0.5) + +# add mean +ax.axvline( + x=rho_data.mean().values, + c="b", + linestyle="--", +) + + +# add median +ax.axvline( + x=rho_data.median().values, + c="r", + linestyle="-", +) + +# add legend +ax.legend( + [ + f"mean = {np.nanmean(rho_data):.2f} pixels", + f"median = {np.nanmedian(rho_data):.2f} pixels", + ], + loc="best", +) +ax.set_ylabel("count") +ax.set_xlabel("rho (pixels)") +fig.show() + +# %% +# We can see that there is some spread in the value of ``rho`` in this +# dataset. This may be due to noise in the detection of the head keypoints, +# or due to the mouse tipping its snout upwards during the recording. + +# %% +# Histogram of ``phi`` values +# ------------------------------------- +# We can also explore which ``phi`` values are most common in the dataset with +# a circular histogram. + +# sphinx_gallery_thumbnail_number = 5 + +# compute number of bins +bin_width_deg = 5 # width of the bins in degrees +n_bins = int(360 / bin_width_deg) + +# initialise figure with polar projection +fig = plt.figure() +ax = fig.add_subplot(projection="polar") + +# plot histogram using xarray's built-in histogram function +head_vector_polar.sel(individuals=mouse_name, space_pol="phi").plot.hist( + bins=np.linspace(-np.pi, np.pi, n_bins + 1), ax=ax +) + +# axes settings +ax.set_title("phi histogram") +ax.set_theta_direction(-1) # theta increases in clockwise direction +ax.set_theta_offset(0) # set zero at the right +ax.set_xlabel("") # remove default x-label from xarray's plot.hist() + +# set xticks to match the phi values in degrees +n_xtick_edges = 9 +ax.set_xticks(np.linspace(0, 2 * np.pi, n_xtick_edges)[:-1]) +xticks_in_deg = ( + list(range(0, 180 + 45, 45)) + list(range(0, -180, -45))[-1:0:-1] +) +ax.set_xticklabels([str(t) + "\N{DEGREE SIGN}" for t in xticks_in_deg]) + +fig.show() + +# %% +# The ``phi`` circular histogram shows that the head vector appears at a +# variety of orientations in this dataset. + +# %% +# Polar plot of the head vector within a time window +# --------------------------------------------------- +# We can also use a polar plot to represent the head vector in time. This way +# we can visualise the head vector in a coordinate system that translates with +# the mouse but is always parallel to the pixel coordinate system. +# Again, this will be easier to visualise if we focus on a +# small time window. + +# select phi values within a time window +phi = head_vector_polar.sel( + individuals=mouse_name, + space_pol="phi", + time=time_window, +).values + +# plot tip of the head vector within that window, and color based on time +fig = plt.figure() +ax = fig.add_subplot(projection="polar") +sc = ax.scatter( + phi, + np.ones_like(phi), # assign a constant value rho=1 for visualization + c=time_window, + cmap="viridis", + s=50, +) + +# axes settings +ax.set_theta_direction(-1) # theta increases in clockwise direction +ax.set_theta_offset(0) # set zero at the right +cax = fig.colorbar( + sc, + ax=ax, + label=f"time ({ds.attrs['time_unit']})", + ticks=list(time_window)[0::2], +) + +# set xticks to match the phi values in degrees +n_xtick_edges = 9 +ax.set_xticks(np.linspace(0, 2 * np.pi, n_xtick_edges)[:-1]) +xticks_in_deg = ( + list(range(0, 180 + 45, 45)) + list(range(0, -180, -45))[-1:0:-1] +) +ax.set_xticklabels([str(t) + "\N{DEGREE SIGN}" for t in xticks_in_deg]) + +fig.show() + +# %% +# In the polar plot above, the midpoint between the ears is at the centre of +# the plot. The tip of the head vector (the ``snout``) is represented with +# color markers at a constant ``rho`` value of 1. Markers are colored by frame. +# The polar plot shows how in this small time window of 20 frames, +# the head of the mouse turned anti-clockwise. + +# %% +# Convert polar coordinates to cartesian +# ------------------------------------------ +# ``movement`` also provides a ``pol2cart`` convenience function to transform +# a vector in polar coordinates back to cartesian. +head_vector_cart = pol2cart(head_vector_polar) + +print(head_vector_cart) + +# %% +# Note that the resulting `head_vector_cart` array has a ``space`` dimension +# with two coordinates: ``x`` and ``y``. diff --git a/_downloads/b14755a378a13dd501d6f42af59fe770/compute_kinematics.py b/_downloads/b14755a378a13dd501d6f42af59fe770/compute_kinematics.py index a2da2645..b2e05c7b 100644 --- a/_downloads/b14755a378a13dd501d6f42af59fe770/compute_kinematics.py +++ b/_downloads/b14755a378a13dd501d6f42af59fe770/compute_kinematics.py @@ -2,7 +2,7 @@ """Compute and visualise kinematics. ==================================== -Compute displacement, velocity and acceleration data on an example dataset and +Compute displacement, velocity and acceleration, and visualise the results. """ diff --git a/_downloads/bc82bea3a5dd7bdba60b65220891d9e5/examples_python.zip b/_downloads/bc82bea3a5dd7bdba60b65220891d9e5/examples_python.zip index 80f3a66c..eb7bb6b0 100644 Binary files a/_downloads/bc82bea3a5dd7bdba60b65220891d9e5/examples_python.zip and b/_downloads/bc82bea3a5dd7bdba60b65220891d9e5/examples_python.zip differ diff --git a/_downloads/f6a82ada996eeecdfd55e392bcd9ac44/compute_polar_coordinates.ipynb b/_downloads/f6a82ada996eeecdfd55e392bcd9ac44/compute_polar_coordinates.ipynb new file mode 100644 index 00000000..4dd930a8 --- /dev/null +++ b/_downloads/f6a82ada996eeecdfd55e392bcd9ac44/compute_polar_coordinates.ipynb @@ -0,0 +1,315 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Express 2D vectors in polar coordinates\n\nCompute a vector representing head direction and express it in polar\ncoordinates.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Imports\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# For interactive plots: install ipympl with `pip install ipympl` and uncomment\n# the following line in your notebook\n# %matplotlib widget\nimport numpy as np\nfrom matplotlib import pyplot as plt\n\nfrom movement import sample_data\nfrom movement.io import load_poses\nfrom movement.utils.vector import cart2pol, pol2cart" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Load sample dataset\nIn this tutorial, we will use a sample dataset with a single individual\n(a mouse) and six keypoints.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "ds_path = sample_data.fetch_dataset_paths(\n \"SLEAP_single-mouse_EPM.analysis.h5\"\n)[\"poses\"]\nds = load_poses.from_sleap_file(ds_path, fps=None) # force time_unit = frames\n\nprint(ds)\nprint(\"-----------------------------\")\nprint(f\"Individuals: {ds.individuals.values}\")\nprint(f\"Keypoints: {ds.keypoints.values}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The loaded dataset ``ds`` contains two data variables:``position`` and\n``confidence``. Both are stored as data arrays. In this tutorial, we will\nuse only ``position``:\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "position = ds.position" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Compute head vector\nTo demonstrate how polar coordinates can be useful in behavioural analyses,\nwe will compute the head vector of the mouse.\n\nWe define it as the vector from the midpoint between the ears to the snout.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# compute the midpoint between the ears\nmidpoint_ears = position.sel(keypoints=[\"left_ear\", \"right_ear\"]).mean(\n dim=\"keypoints\"\n)\n\n# compute the head vector\nhead_vector = position.sel(keypoints=\"snout\") - midpoint_ears\n\n# drop the keypoints dimension\n# (otherwise the `head_vector` data array retains a `snout` keypoint from the\n# operation above)\nhead_vector = head_vector.drop_vars(\"keypoints\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualise the head trajectory\nWe can plot the data to check that our computation of the head vector is\ncorrect.\n\nWe can start by plotting the trajectory of the midpoint between the ears. We\nwill refer to this as the head trajectory.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(1, 1)\nmouse_name = ds.individuals.values[0]\n\nsc = ax.scatter(\n midpoint_ears.sel(individuals=mouse_name, space=\"x\"),\n midpoint_ears.sel(individuals=mouse_name, space=\"y\"),\n s=15,\n c=midpoint_ears.time,\n cmap=\"viridis\",\n marker=\"o\",\n)\n\nax.axis(\"equal\")\nax.set_xlabel(\"x (pixels)\")\nax.set_ylabel(\"y (pixels)\")\nax.invert_yaxis()\nax.set_title(f\"Head trajectory ({mouse_name})\")\nfig.colorbar(sc, ax=ax, label=f\"time ({ds.attrs['time_unit']})\")\nfig.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see that the majority of the head trajectory data is within a\ncruciform shape. This is because the dataset is of a mouse moving on an\n[Elevated Plus Maze](https://en.wikipedia.org/wiki/Elevated_plus_maze).\nWe can actually verify this is the case by overlaying the head\ntrajectory on the sample frame of the dataset.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# read sample frame\nframe_path = sample_data.fetch_dataset_paths(\n \"SLEAP_single-mouse_EPM.analysis.h5\"\n)[\"frame\"]\nim = plt.imread(frame_path)\n\n\n# plot sample frame\nfig, ax = plt.subplots(1, 1)\nax.imshow(im)\n\n# plot head trajectory with semi-transparent markers\nsc = ax.scatter(\n midpoint_ears.sel(individuals=mouse_name, space=\"x\"),\n midpoint_ears.sel(individuals=mouse_name, space=\"y\"),\n s=15,\n c=midpoint_ears.time,\n cmap=\"viridis\",\n marker=\"o\",\n alpha=0.05, # transparency\n)\n\nax.axis(\"equal\")\nax.set_xlabel(\"x (pixels)\")\nax.set_ylabel(\"y (pixels)\")\n# No need to invert the y-axis now, since the image is plotted\n# using a pixel coordinate system with origin on the top left of the image\nax.set_title(f\"Head trajectory ({mouse_name})\")\n\nfig.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The overlaid plot suggests the mouse spends most of its time in the\ncovered arms of the maze.\n\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualise the head vector\nTo visually check our computation of the head vector, it is easier to select\na subset of the data. We can focus on the trajectory of the head when the\nmouse is within a small rectangular area and time window.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# area of interest\nxmin, ymin = 600, 665 # pixels\nx_delta, y_delta = 125, 100 # pixels\n\n# time window\ntime_window = range(1650, 1671) # frames" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For that subset of the data, we now plot the head vector.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(1, 1)\nmouse_name = ds.individuals.values[0]\n\n# plot midpoint between the ears, and color based on time\nsc = ax.scatter(\n midpoint_ears.sel(individuals=mouse_name, space=\"x\", time=time_window),\n midpoint_ears.sel(individuals=mouse_name, space=\"y\", time=time_window),\n s=50,\n c=midpoint_ears.time[time_window],\n cmap=\"viridis\",\n marker=\"*\",\n)\n\n# plot snout, and color based on time\nsc = ax.scatter(\n position.sel(\n individuals=mouse_name, space=\"x\", time=time_window, keypoints=\"snout\"\n ),\n position.sel(\n individuals=mouse_name, space=\"y\", time=time_window, keypoints=\"snout\"\n ),\n s=50,\n c=position.time[time_window],\n cmap=\"viridis\",\n marker=\"o\",\n)\n\n# plot the computed head vector\nax.quiver(\n midpoint_ears.sel(individuals=mouse_name, space=\"x\", time=time_window),\n midpoint_ears.sel(individuals=mouse_name, space=\"y\", time=time_window),\n head_vector.sel(individuals=mouse_name, space=\"x\", time=time_window),\n head_vector.sel(individuals=mouse_name, space=\"y\", time=time_window),\n angles=\"xy\",\n scale=1,\n scale_units=\"xy\",\n headwidth=7,\n headlength=9,\n headaxislength=9,\n color=\"gray\",\n)\n\nax.axis(\"equal\")\nax.set_xlim(xmin, xmin + x_delta)\nax.set_ylim(ymin, ymin + y_delta)\nax.set_xlabel(\"x (pixels)\")\nax.set_ylabel(\"y (pixels)\")\nax.set_title(f\"Zoomed in head vector ({mouse_name})\")\nax.invert_yaxis()\nfig.colorbar(\n sc,\n ax=ax,\n label=f\"time ({ds.attrs['time_unit']})\",\n ticks=list(time_window)[0::2],\n)\n\nax.legend(\n [\n \"midpoint_ears\",\n \"snout\",\n \"head_vector\",\n ],\n loc=\"best\",\n)\n\nfig.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "From the plot we can confirm the head vector goes from the midpoint between\nthe ears to the snout, as we defined it.\n\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Express the head vector in polar coordinates\nA convenient way to inspect the orientation of a vector in 2D is by\nexpressing it in polar coordinates. We can do this with the vector function\n``cart2pol``:\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "head_vector_polar = cart2pol(head_vector)\n\nprint(head_vector_polar)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice how the resulting array has a ``space_pol`` dimension with two\ncoordinates: ``rho`` and ``phi``. These are the polar coordinates of the\nhead vector.\n\nThe coordinate ``rho`` is the norm (i.e., magnitude, length) of the vector.\nIn our case, the distance from the midpoint between the ears to the snout.\nThe coordinate ``phi`` is the orientation of the head vector relative to the\npositive x-axis, and ranges from -``pi`` to ``pi``\n(following the [atan2](https://en.wikipedia.org/wiki/Atan2) convention).\n\nIn our coordinate system, ``phi`` will be\npositive if the shortest path from the positive x-axis to the vector is\nclockwise. Conversely, ``phi`` will be negative if the shortest path from\nthe positive x-axis to the vector is anti-clockwise.\n\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Histogram of ``rho`` values\nSince ``rho`` is the distance between the ears' midpoint and the snout,\nwe would expect ``rho`` to be approximately constant in this data. We can\ncheck this by plotting a histogram of its values across the whole clip.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(1, 1)\n\n# plot histogram using xarray's built-in histogram function\nrho_data = head_vector_polar.sel(individuals=mouse_name, space_pol=\"rho\")\nrho_data.plot.hist(bins=50, ax=ax, edgecolor=\"lightgray\", linewidth=0.5)\n\n# add mean\nax.axvline(\n x=rho_data.mean().values,\n c=\"b\",\n linestyle=\"--\",\n)\n\n\n# add median\nax.axvline(\n x=rho_data.median().values,\n c=\"r\",\n linestyle=\"-\",\n)\n\n# add legend\nax.legend(\n [\n f\"mean = {np.nanmean(rho_data):.2f} pixels\",\n f\"median = {np.nanmedian(rho_data):.2f} pixels\",\n ],\n loc=\"best\",\n)\nax.set_ylabel(\"count\")\nax.set_xlabel(\"rho (pixels)\")\nfig.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see that there is some spread in the value of ``rho`` in this\ndataset. This may be due to noise in the detection of the head keypoints,\nor due to the mouse tipping its snout upwards during the recording.\n\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Histogram of ``phi`` values\nWe can also explore which ``phi`` values are most common in the dataset with\na circular histogram.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# compute number of bins\nbin_width_deg = 5 # width of the bins in degrees\nn_bins = int(360 / bin_width_deg)\n\n# initialise figure with polar projection\nfig = plt.figure()\nax = fig.add_subplot(projection=\"polar\")\n\n# plot histogram using xarray's built-in histogram function\nhead_vector_polar.sel(individuals=mouse_name, space_pol=\"phi\").plot.hist(\n bins=np.linspace(-np.pi, np.pi, n_bins + 1), ax=ax\n)\n\n# axes settings\nax.set_title(\"phi histogram\")\nax.set_theta_direction(-1) # theta increases in clockwise direction\nax.set_theta_offset(0) # set zero at the right\nax.set_xlabel(\"\") # remove default x-label from xarray's plot.hist()\n\n# set xticks to match the phi values in degrees\nn_xtick_edges = 9\nax.set_xticks(np.linspace(0, 2 * np.pi, n_xtick_edges)[:-1])\nxticks_in_deg = (\n list(range(0, 180 + 45, 45)) + list(range(0, -180, -45))[-1:0:-1]\n)\nax.set_xticklabels([str(t) + \"\\N{DEGREE SIGN}\" for t in xticks_in_deg])\n\nfig.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The ``phi`` circular histogram shows that the head vector appears at a\nvariety of orientations in this dataset.\n\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Polar plot of the head vector within a time window\nWe can also use a polar plot to represent the head vector in time. This way\nwe can visualise the head vector in a coordinate system that translates with\nthe mouse but is always parallel to the pixel coordinate system.\nAgain, this will be easier to visualise if we focus on a\nsmall time window.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# select phi values within a time window\nphi = head_vector_polar.sel(\n individuals=mouse_name,\n space_pol=\"phi\",\n time=time_window,\n).values\n\n# plot tip of the head vector within that window, and color based on time\nfig = plt.figure()\nax = fig.add_subplot(projection=\"polar\")\nsc = ax.scatter(\n phi,\n np.ones_like(phi), # assign a constant value rho=1 for visualization\n c=time_window,\n cmap=\"viridis\",\n s=50,\n)\n\n# axes settings\nax.set_theta_direction(-1) # theta increases in clockwise direction\nax.set_theta_offset(0) # set zero at the right\ncax = fig.colorbar(\n sc,\n ax=ax,\n label=f\"time ({ds.attrs['time_unit']})\",\n ticks=list(time_window)[0::2],\n)\n\n# set xticks to match the phi values in degrees\nn_xtick_edges = 9\nax.set_xticks(np.linspace(0, 2 * np.pi, n_xtick_edges)[:-1])\nxticks_in_deg = (\n list(range(0, 180 + 45, 45)) + list(range(0, -180, -45))[-1:0:-1]\n)\nax.set_xticklabels([str(t) + \"\\N{DEGREE SIGN}\" for t in xticks_in_deg])\n\nfig.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the polar plot above, the midpoint between the ears is at the centre of\nthe plot. The tip of the head vector (the ``snout``) is represented with\ncolor markers at a constant ``rho`` value of 1. Markers are colored by frame.\nThe polar plot shows how in this small time window of 20 frames,\nthe head of the mouse turned anti-clockwise.\n\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Convert polar coordinates to cartesian\n``movement`` also provides a ``pol2cart`` convenience function to transform\na vector in polar coordinates back to cartesian.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "head_vector_cart = pol2cart(head_vector_polar)\n\nprint(head_vector_cart)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that the resulting `head_vector_cart` array has a ``space`` dimension\nwith two coordinates: ``x`` and ``y``.\n\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.9" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file diff --git a/_downloads/fb625db3c50d423b1b7881136ffdeec8/examples_jupyter.zip b/_downloads/fb625db3c50d423b1b7881136ffdeec8/examples_jupyter.zip index 0cd38f48..1c09c3d2 100644 Binary files a/_downloads/fb625db3c50d423b1b7881136ffdeec8/examples_jupyter.zip and b/_downloads/fb625db3c50d423b1b7881136ffdeec8/examples_jupyter.zip differ diff --git a/_images/sphx_glr_compute_polar_coordinates_001.png b/_images/sphx_glr_compute_polar_coordinates_001.png new file mode 100644 index 00000000..f94db162 Binary files /dev/null and b/_images/sphx_glr_compute_polar_coordinates_001.png differ diff --git a/_images/sphx_glr_compute_polar_coordinates_002.png b/_images/sphx_glr_compute_polar_coordinates_002.png new file mode 100644 index 00000000..4e69f393 Binary files /dev/null and b/_images/sphx_glr_compute_polar_coordinates_002.png differ diff --git a/_images/sphx_glr_compute_polar_coordinates_003.png b/_images/sphx_glr_compute_polar_coordinates_003.png new file mode 100644 index 00000000..428b5940 Binary files /dev/null and b/_images/sphx_glr_compute_polar_coordinates_003.png differ diff --git a/_images/sphx_glr_compute_polar_coordinates_004.png b/_images/sphx_glr_compute_polar_coordinates_004.png new file mode 100644 index 00000000..ea66e4ff Binary files /dev/null and b/_images/sphx_glr_compute_polar_coordinates_004.png differ diff --git a/_images/sphx_glr_compute_polar_coordinates_005.png b/_images/sphx_glr_compute_polar_coordinates_005.png new file mode 100644 index 00000000..2a1a5adf Binary files /dev/null and b/_images/sphx_glr_compute_polar_coordinates_005.png differ diff --git a/_images/sphx_glr_compute_polar_coordinates_006.png b/_images/sphx_glr_compute_polar_coordinates_006.png new file mode 100644 index 00000000..91fbb753 Binary files /dev/null and b/_images/sphx_glr_compute_polar_coordinates_006.png differ diff --git a/_images/sphx_glr_compute_polar_coordinates_thumb.png b/_images/sphx_glr_compute_polar_coordinates_thumb.png new file mode 100644 index 00000000..8ee72c2e Binary files /dev/null and b/_images/sphx_glr_compute_polar_coordinates_thumb.png differ diff --git a/_images/sphx_glr_smooth_001.png b/_images/sphx_glr_smooth_001.png new file mode 100644 index 00000000..0bffbe3f Binary files /dev/null and b/_images/sphx_glr_smooth_001.png differ diff --git a/_images/sphx_glr_smooth_002.png b/_images/sphx_glr_smooth_002.png new file mode 100644 index 00000000..263458a5 Binary files /dev/null and b/_images/sphx_glr_smooth_002.png differ diff --git a/_images/sphx_glr_smooth_003.png b/_images/sphx_glr_smooth_003.png new file mode 100644 index 00000000..7ba87d87 Binary files /dev/null and b/_images/sphx_glr_smooth_003.png differ diff --git a/_images/sphx_glr_smooth_004.png b/_images/sphx_glr_smooth_004.png new file mode 100644 index 00000000..a9cda89d Binary files /dev/null and b/_images/sphx_glr_smooth_004.png differ diff --git a/_images/sphx_glr_smooth_005.png b/_images/sphx_glr_smooth_005.png new file mode 100644 index 00000000..c42c1caf Binary files /dev/null and b/_images/sphx_glr_smooth_005.png differ diff --git a/_images/sphx_glr_smooth_006.png b/_images/sphx_glr_smooth_006.png new file mode 100644 index 00000000..a53cbc90 Binary files /dev/null and b/_images/sphx_glr_smooth_006.png differ diff --git a/_images/sphx_glr_smooth_thumb.png b/_images/sphx_glr_smooth_thumb.png new file mode 100644 index 00000000..718beeee Binary files /dev/null and b/_images/sphx_glr_smooth_thumb.png differ diff --git a/_modules/index.html b/_modules/index.html index c0f66036..0bf76026 100644 --- a/_modules/index.html +++ b/_modules/index.html @@ -13,19 +13,19 @@ - - - + + + - - - - + + + + @@ -36,11 +36,11 @@ - - - + + + - + @@ -59,27 +59,23 @@ - +
+ Back to top - + id="pst-primary-sidebar-checkbox"/> + - + id="pst-secondary-sidebar-checkbox"/> +
@@ -101,12 +97,17 @@ Ctrl+K
+ +
+ +
+ -