Compute locomotion features

[lochhh]

Description

What is this PR

• Bug fix
• Addition of a new feature
• Other

Why is this PR needed?
This PR resolves #3.

What does this PR do?

• This PR introduces the analysis.kinematics module with basic functions to compute xarray.DataArrays:
  • displacement
  • velocity (1st order derivative)
  • acceleration (2nd order derivative)
• These functions are used to add kinematic properties via PosesAccessor, allowing the corresponding kinematic properties to be computed and stored in the dataset (if not already). This means that calling ds.poses.<kinematic property> (e.g. ds.poses.velocity) will return the velocity DataArray and store this in ds. Since velocity is now stored in ds, one can retrieve this using ds.velocity.
• This PR also adds checks for expected xarray.Dataset.data_vars and dims in PosesAccessor.validate()

References

#3

How has this PR been tested?

New tests have been written and run locally and through CI.

Does this PR require an update to the documentation?

We may want to add this to the "getting started" guide or add to "examples".

Checklist:

• The code has been tested locally
• Tests have been added to cover all new functionality
• The documentation has been updated to reflect any changes
• The code has been formatted with pre-commit

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 99.17%. Comparing base (33e09ed) to head (a642280).
Report is 3 commits behind head on main.

❗ Current head a642280 differs from pull request most recent head e5cb36c. Consider uploading reports for the commit e5cb36c to get more accurate results

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #106      +/-   ##
==========================================
+ Coverage   99.13%   99.17%   +0.04%     
==========================================
  Files           8        9       +1     
  Lines         461      487      +26     
==========================================
+ Hits          457      483      +26     
  Misses          4        4              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[niksirbi]

Hey @lochhh thanks for this!

I have some high-level questions/comments regarding this implementation.

User interface

With the current implementation, the use-case would be something like:

from movement import sample_data
from movement.analysis import kinematics

ds = sample_data.fetch_sample_data(
    "SLEAP_three-mice_Aeon_proofread.predictions.slp"
)

ds["displacement_norm"] = kinematics.displacement_vector(ds["pose_tracks"])["magnitude"]
ds["velocity"] = kinematics.velocity(ds["pose_tracks"])
ds["acceleration"] = kinematics.acceleration(ds["pose_tracks"])
ds["speed"] = kinematics.velocity_vector(ds["pose_tracks"])["magnitude"]

I've played around a bit with using the new module in a notebook, in a separate branch).

This is fine for now, but are we planning to provide a more convenient API for the users later? Perhaps by adding convenience methods/properties to that accessor object? Have I understood correctly that we are leaving that for a later PR?

Regarding the interface as is, there is a bit of repetition because each function appears "twice", once in xy form and once in magnitude/angle form. What do you think about merging the two into something like below?

def velocity(data: xr.DataArray, format: Literal['xy', 'vector']):

There is a bit of an inconsistency I guess, because the xy functions return a DataArray, while the vector methods return a xr.Dataset. Perhaps this will become less consequential if we provide sth like kinematics.compute_all() which will return an xr.Dataset containing all relevant variables as DataArrays.

If we don't end up merging the two formats as one function, we may have to rename them as velocity_xy and velocity_vector, just to be very explicit and to avoid confusion.

A different approach would be to have say one velocity() function returning 3 DataArrays in the same Dataset: xy, magnitude, angle.

I'm thinking out loud here, none of these opinions are held strongly, so I'll open the discussion to @sfmig and @b-peri as well. Feel free to chime in on this one (no rush).

Testing

Currently, you test by generating some toy pose data, where the kinematic variables are known or very easy to derive, and assert that our functions give the expected results. To be on the extra safe side, shouldn't we also test with some more "real-world" data - i.e. data that contains NaN values perhaps and doesn't have acceleration=0? I wouldn't run the tests on all our sample datasets (that seems excessive and inefficient), but maybe a few of them?

[sfmig]

Looking good! ✨

I have two comments but with a caveat: I am not super familiar with the data structure so feel free to take / leave whatever makes most / least sense with the agreed structure. This could also very well be a separate PRs.

• The first comment is about nomenclature. What do you guys think of using commands for the names of the methods, so that the distinction between method and attribute is more clear? Just like we do with fetch_sample_data. In @niksirbi's example, I would suggest something like:

  from movement import sample_data
  from movement.analysis import kinematics
  
  ds = sample_data.fetch_sample_data(
      "SLEAP_three-mice_Aeon_proofread.predictions.slp"
  )
  ds["velocity"] = kinematics.compute_velocity(ds["pose_tracks"])
  ds["acceleration"] = kinematics.compute_acceleration(ds["pose_tracks"])

  I am assuming that compute_velocity and compute_acceleration return a cartesian vector (x,y, and maybe z).

• The second comment is about making slightly more generic methods. I would suggest that instead of computing the norm of the velocity vector like this:

  ds["speed"] = kinematics.velocity_vector(ds["pose_tracks"])["magnitude"]

  we instead have methods that operate on vectors (for example, computing the norm for any input vector). So maybe with the nomenclature from above:

  ds["speed"] = kinematics.compute_velocity(ds["pose_tracks"]).compute_norm()
  ds["heading"] = kinematics.compute_velocity(ds["pose_tracks"]).compute_theta()

  For computing "heading"/theta, another option could be to have methods that transform vectors between coordinate systems: cartesian (x,y) and polar (rho, theta) in 2D, and between cartesian (x,y,z), spherical and cylindrical in 3D. I would personally prefer this option. So something like:

  ds["velocity_polar"] = kinematics.compute_velocity(ds["pose_tracks"]).cart2polar()

  where ds["velocity_polar"] would be a vector whose first component is rho (aka magnitude) and its second component is theta (aka angle with respect to the positive x-axis). IIRC Matlab works like this and I'd say it is intuitive. For reference there was a discussion about including methods like this in numpy, and we may be able to use cmath module.

Happy to discuss in a behaviour meeting or do some pair programming on this!

[niksirbi]

ds["velocity_polar"] = kinematics.compute_velocity(ds["pose_tracks"]).cart2polar()

Thanks Sofia! I would also find this quite intuitive, and I like the idea of making more generic methods for operating on vectors, as I foresee much reuse for those.

We would have to see how to implement these in the context of xarray objects. One idea off the top of my head would be to write an accessor object for data arrays containing cartesian space coords (x,y or x,y,z) and to implement the transformations as methods of that accessor class. But there are also other ways of achieving the same goal.

One thing that this discussion highlights for me is the importance of defining clear and consistent terminology for the various coordinate systems. We should discuss this in our meeting + hackathon this week.

[lochhh]

Thanks @niksirbi and @sfmig for the great suggestions and feedback! I will incorporate the following small changes for now:

• use the more efficient dt = time1 - time
• rename methods to include verb, e.g. compute_velocity
• refactor _vector methods to more generic compute_norm and compute_theta

The idea is to eventually have users call ds.velocity, ds.displacement, etc. and create these on the fly if they do not yet exist. This can be done for instance by overriding xr.Dataset.__getattr__ - but since xarray recommends against subclassing, it makes more sense to use the existing PosesAccessor and do something like:

@property
    def displacement(self) -> xr.DataArray:
        """Return the displacement between consecutive x, y
        locations of each keypoint of each individual.
        """
        pose_tracks = self._obj[self.var_names[0]]
        self._obj["displacement"] = kinematics.compute_displacement(pose_tracks)
        return self._obj["displacement"]

that allows a user to call ds.poses.displacement and have displacement stored as a data variable in ds. We can definitely add these to the current PR.

The vectors require further discussion this week. I like the idea of having more generic methods, e.g. ds.poses.displacement.cart2polar() - whilst noting that the current _vector norm and theta are different from polar coordinates, as the former lack a fixed origin - but it may not be as straightforward as compute_theta(ds.poses.displacement), compute_norm(ds.poses.displacement), or compute_vector_magnitude_direction(ds.poses.displacement) which returns both, as @niksirbi has also pointed out. Also I'm not a fan of the long name compute_vector_magnitude_direction. We should certainly add real-world test cases and can do so during our Friday hackathon.

[niksirbi]

The idea is to eventually have users call ds.velocity, ds.displacement, etc. and create these on the fly if they do not yet exist. This can be done for instance by overriding xr.Dataset.__getattr__ - but since xarray recommends against subclassing, it makes more sense to use the existing PosesAccessor and do something like:

That all sounds very reasonable to me, but one note here: since we will be able access many properties through that object, I wonder whether PosesAccessor is an appropriate name for it. Perhaps we should name it MovementAccessor or MoveAccessor, and it's identifier should simply be move, i.e ds.move.displacement, ds.move.velocity, etc. Just one more terminology-related issue to decide on.

This goes alongside my other idea of renaming pose_tracks to simply position, to harmonise with the new kinematic properties #100.

[sfmig]

Conclusions from chat w/ @niksirbi and @lochhh

Agreed terminology

• default (cartesian) coordinate system:
  • origin top left of the image,
  • x-axis is positive across columns (pointing right),
  • y-axis resulting from the right hand rule; i.e., positive across rows (pointing down).
    This convention matches DLC, napari and most image processing tools.
• position vector for keypoint i at time t: vector from the origin of the coordinate system to keypoint i $$\vec{r_i}(t) = (x_i(t), y_i(t))$$
• velocity vector for keypoint i at time t: instant velocity of keypoint i at time t $$\vec{v_i}(t) = (\dot{x}_i(t), \dot{y}_i(t))$$
• acceleration vector for keypoint i at time t: instant acceleration of keypoint i at time t $$\vec{a_i}(t) = (\ddot{x}_i(t), \ddot{y}_i(t))$$
• displacement vector for keypoint i at time t: vector representing the displacement of a keypoint i with respect to its position at the previous frame t-1 $$\vec{d_i}(t) = (x_i(t) - x_i(t-1), y_i(t) - y_i(t-1))$$

In this PR:

• Rename the accessor to use the keyword move.
• Implement just the cartesian version of the above vectors.

Note that:

• We will assume the data variables position, velocity and acceleration are by default expressed in the above cartesian coordinate system .
• But the same variables can be expressed in polar coordinates; these will be called position_polar, velocity_polar, accelaration_polar. Instead of a space axis, these polar variables will have a space_polar axis, with two coordinates rho and phi.
• We plan to implement vector utils, that for now will have functions that convert cartesian to polar coordinates cart2polar and polar2cart - this is now logged in Implement polar versions of kinematic vectors #115

[lochhh]

• This PR also adds checks for expected xarray.Dataset.data_vars and dims in PosesAccessor.validate()

I'm leaning towards adding this to the PR for #118, and then rebase the current PR on top. What are your thoughts?

[niksirbi]

• This PR also adds checks for expected xarray.Dataset.data_vars and dims in PosesAccessor.validate()

I'm leaning towards adding this to the PR for #118, and then rebase the current PR on top. What are your thoughts?

I think that makes sense, it better fits the scope of #118

[niksirbi]

Well done @lochhh!

I think this is very close to being merged. Apart from trivial comments, I have two more substantial points to raise:

• Now that the accessor has some useful functionality, it should be added to the API docs (in api_index.rst). The class docstring should probably be modiefied to reflect that "pose_tracks" and "confidence" are not the only possible data variables. We should mention that if one of the properties is called, a data variable with the same name will also be created.
• the way these properties are currently implemented, they will be computed the first time a user asks for them, and subsequent calls will skip the computation and not overwrite the property. I am wondering about the use-case where a user computes say velocity, then does some additional filtering of pose tracks (perhaps dropping more outliers) and then asks for velocity again. I'm not sure how to deal with that case. But we can discuss tomorrow altogether, since it's also relevant for PR Implement functions to drop and interpolate over low-confidence datapoints #116

[sonarcloud]

Quality Gate passed

Issues
0 New issues

Measures
0 Security Hotspots
No data about Coverage
0.0% Duplication on New Code

See analysis details on SonarCloud

[lochhh]

Well done @lochhh!

I think this is very close to being merged. Apart from trivial comments, I have two more substantial points to raise:

• Now that the accessor has some useful functionality, it should be added to the API docs (in api_index.rst). The class docstring should probably be modiefied to reflect that "pose_tracks" and "confidence" are not the only possible data variables. We should mention that if one of the properties is called, a data variable with the same name will also be created.
• the way these properties are currently implemented, they will be computed the first time a user asks for them, and subsequent calls will skip the computation and not overwrite the property. I am wondering about the use-case where a user computes say velocity, then does some additional filtering of pose tracks (perhaps dropping more outliers) and then asks for velocity again. I'm not sure how to deal with that case. But we can discuss tomorrow altogether, since it's also relevant for PR Implement functions to drop and interpolate over low-confidence datapoints #116

Thanks @niksirbi for reviewing this! I've incorporated the changes you suggested.

[niksirbi]

Update API docs render well on my machine!

[niksirbi]

LGTM now 🚀