Skip to content
Snippets Groups Projects
Commit 78883e38 authored by Olivier Bertrand's avatar Olivier Bertrand
Browse files

Adapt doc into ipython notebooks

parent 4d5d49aa
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 944 additions and 339 deletions
doc/source/camera_3dtrajectory/cam_calibration.rst deleted 100644 → 0
+ 0
− 27
View file @ 4d5d49aa
Camera Calibration
==================
A camera is an optical intrument for aquiring images. Cameras are composed of a sensors (converting photon to electrical signal) and a lens (foccussing light rays on the sensors). \
Many experimental paradigm require the recording of animal motion with one or more camera. \
The experimenter is more interested by the position of the animal in his or her arena (world \
coordinate system) than by the position of the animal within the image (camera-coordinate system). Therefore, we need to transform the position of the animal in the image to its position in the world. To do so, we need to calibrate our camera. This calibration require two steps:
* remove distortion from the lens (intrinsic calibration)
* determine the position and orientation (i.e. pose) of the camera in the environment (extrinsic calibration)
Intrinsic calibration
---------------------
`With Matlab <https://de.mathworks.com/help/vision/ug/single-camera-calibrator-app.html>`_
`With Opencv <https://docs.opencv.org/3.3.1/dc/dbb/tutorial_py_calibration.html>`_
Extrinsic calibration
---------------------
.. automodule:: bodorientpy.calib
Alternatively you can use Matlab:
`With Matlab <https://de.mathworks.com/help/vision/ref/estimateworldcamerapose.html>`
doc/source/camera_3dtrajectory/cam_triangulation.rst deleted 100644 → 0
+ 0
− 11
View file @ 4d5d49aa
Triangulation
=============
In computer vision triangulation refers to the process of determining a point in 3D space given its projections onto two, or more, images. In order to solve this problem it is necessary to know the parameters of the camera projection function from 3D to 2D for the cameras involved, in the simplest case represented by the camera matrices. Triangulation is sometimes also referred to as reconstruction.
The triangulation problem is in theory trivial. Since each point in an image corresponds to a line in 3D space, all points on the line in 3D are projected to the point in the image. If a pair of corresponding points in two, or more images, can be found it must be the case that they are the projection of a common 3D point x. The set of lines generated by the image points must intersect at x (3D point).
Richard Hartley and Andrew Zisserman (2003). Multiple View Geometry in computer vision. Cambridge University Press
.. automodule:: bodorientpy.triangulate
doc/source/camera_3dtrajectory/examples/extrinsic_calibration.py deleted 100644 → 0
+ 0
− 50
View file @ 4d5d49aa
"""
Script to calibrate your camera system, knowing the intrinsic \
calibration of the camera-lenses
"""
import pandas as pd
import pkg_resources
from bodorientpy.io import opencv as opencv_io
from bodorientpy.calib import calibrates_extrinsic
# The filename used for the manhttan
manhattan3d_filename = pkg_resources.resource_filename(
'bodorientpy', 'resources/calib/manhattan.hdf')
manhattan3d_key = 'manhattan'
# where to save the calibration
filename_fullcalib = 'full_calibration.xml'
# a filelist of intrinsic calibration
filenames_intrinsics = dict()
filenames_intrinsics[0] = pkg_resources.resource_filename(
'bodorientpy', 'resources/calib/calib_20180117_cam_0.xml')
filenames_intrinsics[1] = pkg_resources.resource_filename(
'bodorientpy', 'resources/calib/calib_20180117_cam_2.xml')
# a filelist for manhattan
filenames_manhattans = dict()
filenames_manhattans[0] = pkg_resources.resource_filename(
'bodorientpy', 'resources/calib/manhattan_cam_0_xypts.csv')
filenames_manhattans[1] = pkg_resources.resource_filename(
'bodorientpy', 'resources/calib/manhattan_cam_2_xypts.csv')
# number of cameras
ncams = 3
# threshold to exclude point in top-left corner of manhattan
corner_th = 60 # in pixel
# Load the manhttanx 3d
manhattan_3d = pd.read_hdf(manhattan3d_filename, manhattan3d_key)
# Load the intrinsics camlibraitons
cameras_intrinsics = opencv_io.cameras_intrinsic_calibration(
filenames_intrinsics)
# Load the manhattans
cameras_extrinsics = calibrates_extrinsic(filenames_manhattans,
manhattan_3d,
cameras_intrinsics,
corner_th)
# Save the full calibration
opencv_io.save_cameras_calibration(cameras_intrinsics,
cameras_extrinsics,
filename_fullcalib,
overwrite=True)
doc/source/camera_3dtrajectory/examples/manhattan.py deleted 100644 → 0
+ 0
− 52
View file @ 4d5d49aa
"""
Building your own manhattan
"""
import pandas as pd
import bodorientpy.calib.plot as cplot
manhattan_filename = 'manhattan.hdf'
manhattan_key = 'manhattan'
num_points = 27 # Number of points
manhattan_3d = pd.DataFrame(index=range(num_points),
columns=['x', 'y', 'z'])
# Then you manually fill the manhattan measurement
manhattan_3d.loc[0, :] = [-10, -10, 4.66]
manhattan_3d.loc[1, :] = [-10, -5, 4.39]
manhattan_3d.loc[2, :] = [-10, -0, 4.63]
manhattan_3d.loc[3, :] = [-10, +5, 4.38]
manhattan_3d.loc[4, :] = [-10, +10, 4.85]
manhattan_3d.loc[5, :] = [-8.09, -0.25, 10.13]
manhattan_3d.loc[6, :] = [-5, -10, 4.52]
manhattan_3d.loc[7, :] = [-5, -5, 4.57]
manhattan_3d.loc[8, :] = [-5, -0, 4.36]
manhattan_3d.loc[9, :] = [-5, +5, 4.45]
manhattan_3d.loc[10, :] = [-5, +10, 4.43]
manhattan_3d.loc[11, :] = [0, -10, 4.70]
manhattan_3d.loc[12, :] = [0, -5, 4.57]
manhattan_3d.loc[13, :] = [0, +5, 4.16]
manhattan_3d.loc[14, :] = [0, +10, 4.43]
manhattan_3d.loc[15, :] = [+5, -10, 4.70]
manhattan_3d.loc[16, :] = [+5, -5, 4.56]
manhattan_3d.loc[17, :] = [+5, -0, 4.27]
manhattan_3d.loc[18, :] = [+5, +5, 4.49]
manhattan_3d.loc[19, :] = [+5, +10, 4.50]
manhattan_3d.loc[20, :] = [+8.62, -8.38, 10.19]
manhattan_3d.loc[21, :] = [+8.55, +8.72, 10.09]
manhattan_3d.loc[22, :] = [+10, -10, 4.51]
manhattan_3d.loc[23, :] = [+10, -5, 4.33]
manhattan_3d.loc[24, :] = [+10, -0, 4.69]
manhattan_3d.loc[25, :] = [+10, +5, 4.69]
manhattan_3d.loc[26, :] = [+10, +10, 4.71]
manhattan_3d *= 10 # Because millimeters are nicer to use than centimeters
# And save it
manhattan_3d.to_hdf(manhattan_filename, manhattan_key)
cplot.manhattan(manhattan_3d, unit='mm')
doc/source/camera_3dtrajectory/examples/manhattan2d.py deleted 100644 → 0
+ 0
− 44
View file @ 4d5d49aa
import pkg_resources
import matplotlib.pyplot as plt
import pandas as pd
manhattan_imfile = list()
manhattan_imfile.append(pkg_resources.resource_filename(
'bodorientpy', 'resources/calib/20180117_cam_0.tif'))
manhattan_imfile.append(pkg_resources.resource_filename(
'bodorientpy', 'resources/calib/20180117_cam_2.tif'))
# a filelist for manhattan
filenames_manhattans = list()
filenames_manhattans.append(pkg_resources.resource_filename(
'bodorientpy', 'resources/calib/manhattan_cam_0_xypts.csv'))
filenames_manhattans.append(pkg_resources.resource_filename(
'bodorientpy', 'resources/calib/manhattan_cam_2_xypts.csv'))
f, axarr = plt.subplots(1, 2, figsize=(15, 8))
for i, ax in enumerate(axarr):
image = plt.imread(manhattan_imfile[i])
manhattan_2d = pd.read_csv(filenames_manhattans[i],
names=['x', 'y'],
header=0)
# Show the manhattan image
ax.imshow(image, cmap='gray')
toplotx = manhattan_2d.x
# Because inverted y axis in image
toploty = image.shape[0]-manhattan_2d.y
markersize = 10
# Plot marker
ax.plot(toplotx,
toploty, 'o',
markersize=markersize,
markerfacecolor="None",
markeredgecolor='red',
markeredgewidth=2)
# Plot marker label
for mi, xy in enumerate(zip(toplotx, toploty)):
ax.text(xy[0]+2*markersize,
xy[1]+2*markersize,
'{}'.format(mi), color='r',
horizontalalignment='left')
doc/source/camera_3dtrajectory/examples/manhattan_0001.hdf deleted 100644 → 0
+ 0
− 0
View file @ 4d5d49aa
File deleted
doc/source/camera_3dtrajectory/examples/npairwise_triangulation.py deleted 100644 → 0
+ 0
− 29
View file @ 4d5d49aa
import pkg_resources
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # noqa F401
from bodorientpy.io import opencv
from bodorientpy.io import marker_cam_fromfiles, marker_cam2pts_cam
from bodorientpy.triangulate import triangulate_ncam_pairwise
# Load full calibration
full_calibfile = pkg_resources.resource_filename(
'bodorientpy', 'resources/calib/full_calibration_20180117.xml')
cameras_calib = opencv.load_cameras_calibration(full_calibfile)
# Load markers on camera
filenames = [pkg_resources.resource_filename(
'bodorientpy', 'resources/calib/bee_005_cam_0_bodyxypts.csv'),
pkg_resources.resource_filename(
'bodorientpy', 'resources/calib/bee_005_cam_2_bodyxypts.csv')]
marker_cam = marker_cam_fromfiles(filenames)
nmarkers = len(marker_cam.columns.levels[0])
# Plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
for mark_i in range(3):
pts_cam = marker_cam2pts_cam(marker_cam, mark_i)
pts3d, _, _ = triangulate_ncam_pairwise(cameras_calib, pts_cam)
ax.plot(pts3d[:, 0], pts3d[:, 1], pts3d[:, 2])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
doc/source/camera_3dtrajectory/index.rst deleted 100644 → 0
+ 0
− 9
View file @ 4d5d49aa
===================
Video to trajectory
===================
.. toctree::
:maxdepth: 2
cam_calibration.rst
cam_triangulation.rst
doc/source/classification/examples/saccade_extraction.py deleted 100644 → 0
+ 0
− 38
View file @ 4d5d49aa
from bodorientpy.trajectory.random import saccadic_traj
from bodorientpy.classification.sacintersaccade import extract
import matplotlib.pyplot as plt
import numpy as np
# fix the seed for constant example
np.random.seed(0)
# generate fake saccadic data
trajectory, saccade_df = saccadic_traj()
#
thresholds = [2, 3]
sacintersac, angvel = extract(trajectory, thresholds)
f, axarr = plt.subplots(3, 1)
ax = axarr[0]
trajectory.loc[:, ['x', 'y', 'z']].plot(ax=ax)
ax = axarr[1]
trajectory.loc[:, ['angle_0', 'angle_1', 'angle_2']].plot(ax=ax)
ax = axarr[2]
angvel.plot(ax=ax)
ax.axhline(min(thresholds), color='r', linestyle=':')
ax.axhline(max(thresholds), color='r', linestyle='-.')
ylim = ax.get_ylim()
ax.fill_between(sacintersac.index, min(ylim), max(ylim),
where=sacintersac.saccade >= 0,
facecolor='green',
alpha=0.2)
ax.fill_between(sacintersac.index, min(ylim), max(ylim),
where=sacintersac.intersaccade >= 0,
facecolor='black',
alpha=0.2)
f.show()
print(sacintersac)
import time
time.sleep(6)
doc/source/classification/examples/saccade_loliplot.py deleted 100644 → 0
+ 0
− 26
View file @ 4d5d49aa
import numpy as np
from bodorientpy.trajectory.random import saccadic_traj
from bodorientpy.trajectory.plot import get_color_frame_dataframe
from bodorientpy.trajectory.plot import lollipops
import matplotlib.pyplot as plt
# fix the seed for constant example
np.random.seed(0)
# generate fake saccadic data
trajectory, saccade_df = saccadic_traj()
# for colors
colors, sm = get_color_frame_dataframe(frame_range=[
trajectory.index.min(),
trajectory.index.max()], cmap=plt.get_cmap('jet'))
f = plt.figure(figsize=(15, 10))
ax = f.add_subplot(111, projection='3d')
lollipops(ax, trajectory,
colors, step_lollipop=20,
offset_lollipop=1, lollipop_marker='o')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
cb = f.colorbar(sm)
cb.set_label('frame number')
f.show()
doc/source/classification/index.rst deleted 100644 → 0
+ 0
− 5
View file @ 4d5d49aa
==============
Classification
==============
.. include:: sac_intersaccade.rst
doc/source/classification/sac_intersaccade.rst deleted 100644 → 0
+ 0
− 34
View file @ 4d5d49aa
Saccade and intersaccade
========================
The trajectories of flying insects are often composed of saccade, \
a rapide rotation, and intersaccade, period of no rapid turn. During \
a saccade the angular velocity can reach several hundred or thousands of \
degrees per second.
An example of a saccadic flight
-------------------------------
.. plot:: classification/examples/saccade_loliplot.py
Extraction
----------
To extract the saccade from a trajectory, we need to use the angular \
velocity, and two thresholds. When the angular velocity is above the \
highest threshold, we have a saccade. When the angular velocity is \
below the lowest threshold, we have an intersaccade. Part of the \
trajectory between the threshold are unclassified. Indeed it could part\
of the saccade or intersaccade, depending on the noise presents in the \
trajectory.
.. literalinclude:: examples/saccade_extraction.py
:lines: 2,13,14
.. plot:: classification/examples/saccade_extraction.py
Prototypical mouvement
======================
See Elke Braun paper (under construction)
doc/source/index.rst
+ 0
− 2
View file @ 78883e38
...@@ -24,8 +24,6 @@ Content ...@@ -24,8 +24,6 @@ Content
gettingstarted gettingstarted
overview/index overview/index
classification/index
camera_3dtrajectory/index
orientation/index orientation/index
tutorials/index tutorials/index
references/index references/index
......
doc/source/tutorials/02-recording-animal-trajectory.ipynb
+ 283
− 0
View file @ 78883e38
Source diff could not be displayed: it is too large. Options to address this: view the blob.
doc/source/tutorials/04-error-propagation.ipynb 0 → 100644
+ 248
− 0
View file @ 78883e38
%% Cell type:markdown id: tags:
# Numerical error propagation
## Sources of observational errors
Observational error (or measurement error) is the difference between \
a measured value of quantity and its true value. An error is not \
necessarly a "mistake" (e.g. in Statistics). Variability is an \
inherent part of things being measured and of the measurement process.
Measurement errors can be divided into two components:
* random error
* systematic error (bias).
Types of bias (non exhaustiv):
* *Selection bias* (Berksonian bias) involves individuals being more \
likely to be selected for study than others.
* *The bias of an estimator* is the difference between an estimator's \
expectations and the true value of the parameter being estimated.
* *Omitted-variable bias* is the bias that appears in estimates of \
parameters in a regression analysis when the assumed specification \
omits an independent variable that should be in the model.
* *Detection bias* occurs when a phenomenon is more likely to be \
observed for a particular set of study subjects.
* *Funding bias* may lead to selection of outcomes, test samples, \
or test procedures that favor a study's financial sponsor.
* *Reporting bias* involves a skew in the availability of data, such \
that observations of a certain kind are more likely to be reported.
* *Analytical bias* arise due to the way that the results are evaluated.
* *Exclusion bias* arise due to the systematic exclusion of certain \
individuals from the study.
* *Observer bias* arises when the researcher unconsciously influences \
the experiment due to cognitive bias where judgement may alter how an \
experiment is carried out / how results are recorded.
## Propagation of uncertainty
In statistics, propagation of uncertainty (or propagation of error) \
is the effect of variables' uncertainties (or errors) on the \
uncertainty of a function based on them. When the variables are \
the values of experimental measurements they have uncertainties \
due to measurement limitations (e.g., instrument precision) which \
propagate to the combination of variables in the function.
## Non-linear function of one variable
Let's be $f(x)$ a non linear function of variables $x$, and \
$\Delta x$ the measurement error of $x$.
We see that the measurement error of $x$ (gray) propagate on $y$ \
(red), via $f(x)$ (black line). The propagation depends on the \
measured value $x_{o}$, the slope of $f(x)$, the curvature of \
$f(x)$, ... , i.e. the derivatives of $f(x)$ relativ to $x$.
The propagated error $\Delta y$ can therefore be express as a \
Taylor-expansion of $f(x)$ around $x_0$ the measured value, i.e. \
calculated the value $f(x_{o}+\Delta x)$ for $\Delta x \approx 0$
$$f(x_{o}+\Delta x) = \sum_{n=0}^{\inf} \
\frac{f^{(n)}(x_{o})}{n!}(\Delta x)^{n}$$
Using only the first-order of the equation:
$$f(x_{o}+\Delta x) \approx f(x_{o})+\frac{df}{dx}(x_{o})\Delta x$$
The propagated error is defined has \
$\Delta y=\left|f(x_{o}+\Delta x)-f(x_{o})\right|$. From the \
equation above we find that \
$\Delta y=\left|\frac{df}{dx}(x_{o})\right|\Delta x $
Instead of deriving the equation, one can approximate the derivitive \
of the function numerically. For example, to propagate the error on \
the function $f(x)=(x+1)(x-2)(x+3)$, one need to \
define the function, the error $\delta_x$, and the point of \
interest $x_0$.
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
from navipy.errorprop import propagate_error
%matplotlib inline
```
%% Cell type:code id: tags:
``` python
def func_f(x):
return (x+1)*(x-2)*(x+3)
delta_x = 0.5
```
%% Cell type:markdown id: tags:
The the error can be calculated as follow:
%% Cell type:code id: tags:
``` python
x0 = -1
delta_y0 = propagate_error(func_f, x0, delta_x)
# at another point
x1 = 3
delta_y1 = propagate_error(func_f, x1, delta_x)
```
%% Cell type:markdown id: tags:
In the following figure, the of error on x has been propagated at two \
different locations:
%% Cell type:code id: tags:
``` python
def plot_error_f(ax, func_f, x, delta_x, delta_y, x0):
y0 = func_f(x0)
xdelta = np.linspace(x0-delta_x, x0+delta_x, len(x))
ydelta = np.linspace(y0-delta_y, y0+delta_y, len(x))
y = func_f(x)
ax.fill_between(xdelta, min(y), max(y), facecolor='k',
alpha=0.3, edgecolor='None')
ax.fill_between(x, min(ydelta), max(ydelta),
facecolor='r', alpha=0.5, edgecolor='None')
ax.plot(x, y, 'k', linewidth=2)
ax.set_xlim([min(x), max(x)])
ax.set_ylim([min(y), max(y)])
ax.set_xlabel('x')
ax.set_ylabel('y=f(x)')
# Assin the x values to plot
x = np.linspace(-5, 4.5, 100)
# Create figure and plot
fig, axarr = plt.subplots(1, 2, figsize=(15, 5))
plot_error_f(axarr[0], func_f, x, delta_x, delta_y0, x0=x0)
plot_error_f(axarr[1], func_f, x, delta_x, delta_y1, x0=x1)
fig.show()
```
%% Output
/home/bolirev/.virtualenv/toolbox-navigation/lib/python3.6/site-packages/matplotlib-2.1.0-py3.6-linux-x86_64.egg/matplotlib/figure.py:418: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
"matplotlib is currently using a non-GUI backend, "
%% Cell type:markdown id: tags:
## Non-linear function of multiple variables
We can extend the framework developped above to function $f$ of \
multiple variables, $x_0,x_1,...x_i,...x_n$. When f is a set of \
non-linear combination of the variables x, an interval propagation \
could be performed in order to compute intervals which contain all \
consistent values for the variables. In a probabilistic approach, \
the function f must usually be linearized by approximation to a \
first-order Taylor series expansion, though in some cases, exact \
formulas can be derived that do not depend on the expansion as is \
the case for the exact variance of products. The Taylor expansion \
would be:
$$f_k \approx f^0_k+ \sum_i^n \frac{\partial f_k}{\partial {x_i}} x_i$$
where \partial f_k/\partial x_i denotes the partial derivative of fk \
with respect to the i-th variable, evaluated at the mean value of \
all components of vector x. Or in matrix notation,
$$\mathrm{f} \approx \mathrm{f}^0 + \mathrm{J} \mathrm{x}$$
where J is the Jacobian matrix. Since f_0 is a constant it does not \
contribute to the error on f. Therefore, the propagation of error \
follows the linear case, above, but replacing the linear \
coefficients, Aik and Ajk by the partial derivatives, \
$\frac{\partial f_k}{\partial x_i}$ and \
$\frac{\partial f_k}{\partial x_j}$. In matrix notation,
$$\mathrm{\Sigma}^\mathrm{f} = \mathrm{J}
\mathrm{\Sigma}^\mathrm{x} \mathrm{J}^\top$$
That is, the Jacobian of the function is used to transform the rows \
and columns of the variance-covariance matrix of the argument.
doc/source/tutorials/05-classification.ipynb 0 → 100644
+ 172
− 0
View file @ 78883e38
%% Cell type:markdown id: tags:
# Classification
## Saccade and intersaccade
The trajectories of flying insects are often composed of saccade, \
a rapide rotation, and intersaccade, period of no rapid turn. During \
a saccade the angular velocity can reach several hundred or thousands of \
degrees per second.
An example of a saccadic flight
%% Cell type:code id: tags:
``` python
import numpy as np
from navipy.trajectories.random import saccadic_traj
from navipy.tools.plots import get_color_frame_dataframe
import matplotlib.pyplot as plt
%matplotlib inline
```
%% Cell type:code id: tags:
``` python
# fix the seed for constant example
np.random.seed(0)
# generate fake saccadic data
trajectory, saccade_df = saccadic_traj()
# for colors
colors, sm = get_color_frame_dataframe(frame_range=[
trajectory.index.min(),
trajectory.index.max()], cmap=plt.get_cmap('jet'))
# Plot he trajectory in lollipops
f = plt.figure(figsize=(15, 10))
ax = f.add_subplot(111, projection='3d')
trajectory.lollipops(ax,
colors, step_lollipop=20,
offset_lollipop=1, lollipop_marker='o')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
cb = f.colorbar(sm)
cb.set_label('frame number')
```
%% Output
%% Cell type:markdown id: tags:
To extract the saccade from a trajectory, we need to use the angular
velocity, and two thresholds. When the angular velocity is above the
highest threshold, we have a saccade. When the angular velocity is
below the lowest threshold, we have an intersaccade. Part of the
trajectory between the threshold are unclassified. Indeed it could part
of the saccade or intersaccade, depending on the noise presents in the
trajectory.
%% Cell type:code id: tags:
``` python
from navipy.trajectories.tools import extract_sacintersac
thresholds = [0.02, 0.03]
sacintersac, angvel = extract_sacintersac(trajectory, thresholds)
```
%% Cell type:markdown id: tags:
And then plot for example the result of the saccade intersaccade
extraction as follow:
%% Cell type:code id: tags:
``` python
f, axarr = plt.subplots(3, 1)
ax = axarr[0]
trajectory.loc[:, 'location'].plot(ax=ax)
ax = axarr[1]
trajectory.loc[:, trajectory.rotation_mode].plot(ax=ax)
ax = axarr[2]
angvel.plot(ax=ax)
ax.axhline(min(thresholds), color='r', linestyle=':')
ax.axhline(max(thresholds), color='r', linestyle='-.')
ylim = ax.get_ylim()
ax.fill_between(sacintersac.index, min(ylim), max(ylim),
where=sacintersac.saccade >= 0,
facecolor='green',
alpha=0.2)
ax.fill_between(sacintersac.index, min(ylim), max(ylim),
where=sacintersac.intersaccade >= 0,
facecolor='black',
alpha=0.2)
f.show()
```
%% Output
/home/bolirev/.virtualenv/toolbox-navigation/lib/python3.6/site-packages/matplotlib-2.1.0-py3.6-linux-x86_64.egg/matplotlib/figure.py:418: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
"matplotlib is currently using a non-GUI backend, "
doc/source/tutorials/index.rst
+ 1
− 0
View file @ 78883e38
...@@ -8,6 +8,7 @@ Tutorials ...@@ -8,6 +8,7 @@ Tutorials
02-recording-animal-trajectory.ipynb 02-recording-animal-trajectory.ipynb
03-rendering-along-trajectory.ipynb 03-rendering-along-trajectory.ipynb
04-error-propagation.ipynb 04-error-propagation.ipynb
05-classification.ipynb
Average place-code vector homing Average place-code vector homing
-------------------------------- --------------------------------
......
navipy/maths/euler.py
+ 6
− 6
View file @ 78883e38
...@@ -132,13 +132,13 @@ def angle_rate_matrix(ai, aj, ak, axes='sxyz'): ...@@ -132,13 +132,13 @@ def angle_rate_matrix(ai, aj, ak, axes='sxyz'):
veck = np.zeros(3) veck = np.zeros(3)
veck[k] = 1 veck[k] = 1
if frame == 0: # static if frame == 0: # static
roti = matrix(ai, 0, 0, axes=axes) roti = matrix(ai, 0, 0, axes=axes)[:3, :3]
rotj = matrix(0, aj, 0, axes=axes) rotj = matrix(0, aj, 0, axes=axes)[:3, :3]
mat = np.hstack([veci, roti.dot(vecj), roti.dot(rotj).dot(veck)]) mat = np.vstack([veci, roti.dot(vecj), roti.dot(rotj).dot(veck)])
else: else:
rotk = matrix(0, 0, ak, axes=axes).transpose() rotk = matrix(0, 0, ak, axes=axes)[:3, :3].transpose()
rotj = matrix(0, aj, 0, axes=axes).transpose() rotj = matrix(0, aj, 0, axes=axes)[:3, :3].transpose()
mat = np.hstack([rotk.dot(rotj.dot(veci)), rotk.dot(vecj), veck]) mat = np.vstack([rotk.dot(rotj.dot(veci)), rotk.dot(vecj), veck])
return mat return mat
......
navipy/tools/__init__.py
+ 71
− 0
View file @ 78883e38
"""
Some tools
"""
import pandas as pd
import numpy as np
def extract_block(data, thresholds):
"""
data is a pandas series with integer consecutive indexes
return a dataframe, index by block containing block start and end
a block is defined as:
- data are greater than threshold_1
- extend to the last value of data below threshold_2 \
before the block
- extend to the first value of data below threshold_2 \
after the block
threshold_1 should be higher than threshold_2
"""
threshold_2 = min(thresholds)
threshold_1 = max(thresholds)
# create a dataframe containing:
# in column th_1: 1 for data above th_1, nan other wise
# in column th_2: 1 for data below th_2, nan other wise
treshold_df = pd.DataFrame(index=data.index,
columns=['th_1', 'th_2'])
treshold_df.th_1[data > threshold_1] = 1
treshold_df.th_2[data <= threshold_2] = 1
treshold_df['frame'] = treshold_df.index
# Calculate unextended block
subdf = treshold_df[['th_1', 'frame']].dropna(how='any')
start_block = subdf.frame[(subdf.frame.diff() > 1)
| (subdf.frame.diff().isnull())]
end_block = subdf.frame[(subdf.frame[::-1].diff() < -1)
| (subdf.frame[::-1].diff().isnull())]
block_df = pd.DataFrame({'start_th1': start_block.reset_index().frame,
'end_th1': end_block.reset_index().frame})
# extend block based on th_2
block_df['end_th2'] = np.nan
block_df['start_th2'] = np.nan
for i, row in block_df.iterrows():
point_below_th2 = treshold_df.loc[:row.start_th1,
'th_2'].dropna()
# check if a point is below th2
if len(point_below_th2) > 0:
start_th2 = point_below_th2.index[-1]
else:
start_th2 = row.start_th1
# Check if the point is before th1
if start_th2 <= row.start_th1:
block_df.loc[i, 'start_th2'] = start_th2
point_below_th2 = treshold_df.loc[row.start_th1:, 'th_2'].dropna()
# check if a point is below th2
if len(point_below_th2) > 0:
end_th2 = point_below_th2.index[0]
else:
end_th2 = row.end_th1
# Check if the point is after th1
if end_th2 >= row.end_th1:
block_df.loc[i, 'end_th2'] = end_th2
return block_df
def extract_block_nonans(data):
return extract_block(data.isnull() == 0, thresholds=[0.4, 0.5])
navipy/trajectories/__init__.py
+ 163
− 6
View file @ 78883e38
...@@ -6,8 +6,10 @@ import numpy as np ...@@ -6,8 +6,10 @@ import numpy as np
import navipy.maths.constants as mconst import navipy.maths.constants as mconst
import navipy.maths.quaternion as htq import navipy.maths.quaternion as htq
import navipy.maths.euler as hte import navipy.maths.euler as hte
import navipy.maths.homogeneous_transformations as htf
from navipy.errorprop import propagate_error from navipy.errorprop import propagate_error
from transformations import markers2translate, markers2euler from .transformations import markers2translate, markers2euler
from navipy.tools.plots import get_color_dataframe
def _markers2position(x, kwargs): def _markers2position(x, kwargs):
...@@ -403,14 +405,14 @@ class Trajectory(pd.DataFrame): ...@@ -403,14 +405,14 @@ class Trajectory(pd.DataFrame):
if error is not None: if error is not None:
euclidian_error = error.loc[index_i] euclidian_error = error.loc[index_i]
if not np.isnan(euclidian_error): if not np.isnan(euclidian_error):
covar = euclidian_error*np.eye(9) covar = euclidian_error * np.eye(9)
err_pos = propagate_error(_markers2position, x, err_pos = propagate_error(_markers2position, x,
covar, args=kwargs, covar, args=kwargs,
epsilon=euclidian_error/10) epsilon=euclidian_error / 10)
err_angle = propagate_error(_markers2angles, x, err_angle = propagate_error(_markers2angles, x,
covar, args=kwargs, covar, args=kwargs,
epsilon=euclidian_error/10) epsilon=euclidian_error / 10)
self.trajectory_error.loc[index_i, 'location'] = \ self.trajectory_error.loc[index_i, 'location'] = \
np.diagonal(err_pos) np.diagonal(err_pos)
self.trajectory_error.loc[index_i, self.rotation_mode] \ self.trajectory_error.loc[index_i, self.rotation_mode] \
...@@ -419,5 +421,160 @@ class Trajectory(pd.DataFrame): ...@@ -419,5 +421,160 @@ class Trajectory(pd.DataFrame):
self.loc[index_i, 'location'] = position self.loc[index_i, 'location'] = position
self.loc[index_i, self.rotation_mode] = orientation self.loc[index_i, self.rotation_mode] = orientation
def lollipops(self): def world2body(self, markers):
raise NameError('Not implemented') """ Transform markers in world coordinate to body coordinate
"""
transformed_markers = pd.DataFrame(index=self.index,
columns=markers.index,
dtype=float)
for index_i, row in self.iterrows():
angles = row.loc[self.rotation_mode].values
translate = row.loc['location'].values
trans_mat = htf.compose_matrix(angles=angles,
translate=translate,
axes=self.rotation_mode)
for marki in markers.index.levels[0]:
markpos = [markers.loc[marki, 'x'],
markers.loc[marki, 'y'],
markers.loc[marki, 'z'],
1]
markpos = trans_mat.dot(markpos)[:3]
transformed_markers.loc[index_i,
(marki, ['x', 'y', 'z'])] = markpos
return transformed_markers
def velocity(self):
""" Calculate the velocity on a trajectory
"""
velocity = pd.DataFrame(index=self.index,
columns=['dx', 'dy', 'dz',
'dalpha_0',
'dalpha_1',
'dalpha_2'],
dtype=float)
diffrow = self.diff()
velocity.loc[:, ['dx', 'dy', 'dz']] = \
diffrow.loc[:, 'location'].values
for index_i, row in self.iterrows():
if self.rotation_mode == 'quaternion':
raise NameError('Not implemented')
else:
rot = hte.angular_velocity(
ai=row.loc[(self.rotation_mode, 'alpha_0')],
aj=row.loc[(self.rotation_mode, 'alpha_1')],
ak=row.loc[(self.rotation_mode, 'alpha_2')],
dai=diffrow.loc[index_i, (self.rotation_mode,
'alpha_0')],
daj=diffrow.loc[index_i, (self.rotation_mode,
'alpha_1')],
dak=diffrow.loc[index_i, (self.rotation_mode,
'alpha_2')],
axes=self.rotation_mode)
velocity.loc[index_i, ['dalpha_0',
'dalpha_1',
'dalpha_2']] = rot.squeeze()
return velocity
def lollipops(self, ax,
colors=None, step_lollipop=1,
offset_lollipop=0, lollipop_marker='o',
linewidth=1, lollipop_tail_width=1,
lollipop_head_size=1
):
""" lollipops plot
create a lollipop plot for the trajectory with its associated \
direction. Handels missing frames by leaving
gaps in the lollipop plot. However indices of the colors, trajectory
**Note** Gap in the index of the trajectory dataframe, will \
lead to gap in the plotted trajectory and the color bar will \
therefore miss some values (otherwise the color-coded frames \
will not be linear)
:param ax: is a matplotlib axes with 3d projection
:param trajectory: is a pandas dataframe with columns \
['x','y','z', 'euler_0','euler_1','euler_2']
:param euler_axes: the axes rotation convention \
(see homogenous.transformations for details)
:param colors: is a pandas dataframe with columns \
['r','g','b','a'] and indexed as trajectory. \
(default: time is color coded)
:param step_lollipop: number of frames between two lollipops
:param offset_lollipop: the first lollipop to be plotted
:param lollipop_marker: the head of the lollipop
:param linewidth:
:param lollipop_tail_width:
:param lollipop_head_size:
"""
direction = self.facing_direction()
if colors is None:
timeseries = pd.Series(data=self.index,
index=self.index)
colors, sm = get_color_dataframe(timeseries)
# Create a continuous index from colors
frames = np.arange(colors.index.min(), colors.index.max() + 1)
# Calculate agent tail direction
# use function in trajectory to get any point bodyref to worldref
tailmarker = pd.Series(data=0,
index=pd.MultiIndex.from_product(
[[0],
['x', 'y', 'z']]))
tailmarker.loc[0, 'x'] = 1
tail = self.world2body(tailmarker)
tail = tail.loc[:, 0]
# Plot the agent trajectory
# - loop through consecutive point
# - Two consecutive point are associated with the color
# of the first point
x = self.loc[:, ('location', 'x')]
y = self.loc[:, ('location', 'y')]
z = self.loc[:, ('location', 'z')]
for frame_i, frame_j in zip(frames[:-1], frames[1:]):
# Frames may be missing in trajectory,
# and therefore can not be plotted
if (frame_i in self.index) and \
(frame_j in self.index) and \
(frame_i in self.index):
color = [colors.r[frame_i], colors.g[frame_i],
colors.b[frame_i], colors.a[frame_i]]
# Create the line to plot
line_x = [x[frame_i], x[frame_j]]
line_y = [y[frame_i], y[frame_j]]
line_z = [z[frame_i], z[frame_j]]
# Actual plot command
ax.plot(xs=line_x, ys=line_y, zs=line_z,
color=color, linewidth=linewidth)
# Plot the lollipop
# - loop through the frames with a step of step_lollipop
# - The lollipop is colored with the color of this frame
# - Each lollipop is composed of a marker,
# a point on the agent trajectory
# and a line representing the body (anti facing direction)
for frame_i in frames[offset_lollipop::step_lollipop]:
# Frames may be missing in trajectory,
# and therefore can not be plotted
if (frame_i in self.index) and \
(frame_i in direction.index) and \
(frame_i in colors.index):
color = [colors.r[frame_i], colors.g[frame_i],
colors.b[frame_i], colors.a[frame_i]]
# Create the line to plot
line_x = [self.x[frame_i],
tail.x[frame_i]]
line_y = [self.y[frame_i],
tail.y[frame_i]]
line_z = [self.z[frame_i],
tail.z[frame_i]]
# Actual plot command
ax.plot(xs=line_x, ys=line_y, zs=line_z,
color=color, linewidth=lollipop_tail_width)
ax.plot(xs=[line_x[0]],
ys=[line_y[0]],
zs=[line_z[0]],
color=color,
marker=lollipop_marker,
markersize=lollipop_head_size)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment