diff --git a/movement/kinematics.py b/movement/kinematics.py index 2880a4888..77af99e22 100644 --- a/movement/kinematics.py +++ b/movement/kinematics.py @@ -1,6 +1,7 @@ """Compute kinematic variables like velocity and acceleration.""" import itertools +from collections.abc import Hashable from typing import Literal import numpy as np @@ -205,8 +206,8 @@ def compute_speed(data: xr.DataArray) -> xr.DataArray: def compute_forward_vector( data: xr.DataArray, - left_keypoint: str, - right_keypoint: str, + left_keypoint: Hashable, + right_keypoint: Hashable, camera_view: Literal["top_down", "bottom_up"] = "top_down", ) -> xr.DataArray: """Compute a 2D forward vector given two left-right symmetric keypoints. @@ -223,9 +224,9 @@ def compute_forward_vector( The input data representing position. This must contain the two symmetrical keypoints located on the left and right sides of the body, respectively. - left_keypoint : str + left_keypoint : Hashable Name of the left keypoint, e.g., "left_ear" - right_keypoint : str + right_keypoint : Hashable Name of the right keypoint, e.g., "right_ear" camera_view : Literal["top_down", "bottom_up"], optional The camera viewing angle, used to determine the upwards @@ -357,8 +358,8 @@ def compute_head_direction_vector( def compute_forward_vector_angle( data: xr.DataArray, - left_keypoint: str, - right_keypoint: str, + left_keypoint: Hashable, + right_keypoint: Hashable, reference_vector: xr.DataArray | ArrayLike = (1, 0), camera_view: Literal["top_down", "bottom_up"] = "top_down", in_radians: bool = False, @@ -371,7 +372,7 @@ def compute_forward_vector_angle( Forward vector angle is the :func:`signed angle\ ` between the reference vector and the animal's :func:`forward vector\ - `). + `. The returned angles are in degrees, spanning the range :math:`(-180, 180]`, unless ``in_radians`` is set to ``True``. @@ -381,10 +382,10 @@ def compute_forward_vector_angle( The input data representing position. This must contain the two symmetrical keypoints located on the left and right sides of the body, respectively. - left_keypoint : str + left_keypoint : Hashable Name of the left keypoint, e.g., "left_ear", used to compute the forward vector. - right_keypoint : str + right_keypoint : Hashable Name of the right keypoint, e.g., "right_ear", used to compute the forward vector. reference_vector : xr.DataArray | ArrayLike, optional diff --git a/movement/roi/base.py b/movement/roi/base.py index 2f93877d0..4bfe139d7 100644 --- a/movement/roi/base.py +++ b/movement/roi/base.py @@ -2,43 +2,59 @@ from __future__ import annotations -from collections.abc import Sequence +from collections.abc import Callable, Hashable, Sequence from typing import Literal, TypeAlias +import numpy as np import shapely +import xarray as xr from numpy.typing import ArrayLike from shapely.coords import CoordinateSequence from movement.utils.broadcasting import broadcastable_method from movement.utils.logging import log_error +from movement.utils.vector import compute_signed_angle_2d LineLike: TypeAlias = shapely.LinearRing | shapely.LineString -PointLike: TypeAlias = tuple[float, float] +PointLike: TypeAlias = list[float] | tuple[float, ...] PointLikeList: TypeAlias = Sequence[PointLike] RegionLike: TypeAlias = shapely.Polygon SupportedGeometry: TypeAlias = LineLike | RegionLike class BaseRegionOfInterest: - """Base class for representing regions of interest (RoIs). + """Base class for regions of interest (RoIs). Regions of interest can be either 1 or 2 dimensional, and are represented - by appropriate ``shapely.Geometry`` objects depending on which. Note that - there are a number of discussions concerning subclassing ``shapely`` - objects; - - - https://github.com/shapely/shapely/issues/1233. - - https://stackoverflow.com/questions/10788976/how-do-i-properly-inherit-from-a-superclass-that-has-a-new-method - - To avoid the complexities of subclassing ourselves, we simply elect to wrap - the appropriate ``shapely`` object in the ``_shapely_geometry`` attribute, - accessible via the property ``region``. This also has the benefit of - allowing us to 'forbid' certain operations (that ``shapely`` would - otherwise interpret in a set-theoretic sense, giving confusing answers to - users). - - This class is not designed to be instantiated directly. It can be - instantiated, however its primary purpose is to reduce code duplication. + by corresponding ``shapely.Geometry`` objects. + + To avoid the complexities of subclassing ``shapely`` objects (due to them + relying on ``__new__()`` methods, see + https://github.com/shapely/shapely/issues/1233), we simply wrap the + relevant ``shapely`` object in the ``_shapely_geometry`` attribute of the + class. This is accessible via the property ``region``. This also allows us + to forbid certain operations and make the manipulation of ``shapely`` + objects more user friendly. + + Although this class can be instantiated directly, it is not designed for + this. Its primary purpose is to reduce code duplication. + + Notes + ----- + A region of interest includes the points that make up its boundary and the + points contained in its interior. This convention means that points inside + the region will be treated as having zero distance to the region, and the + approach vector from these points to the region will be the null vector. + + This may be undesirable in certain situations, when we explicitly want to + find the distance of a point to the boundary of a region, for example. To + accommodate this, most methods of this class accept a keyword argument that + forces the method to perform all computations using only the boundary of + the region, rather than the region itself. For polygons, this will force + the method in question to only consider distances or closest points on the + segments that form the (interior and exterior) boundaries. For segments, + the boundary is considered to be just the two endpoints of the segment. + """ __default_name: str = "Un-named region" @@ -90,6 +106,103 @@ def region(self) -> SupportedGeometry: """``shapely.Geometry`` representation of the region.""" return self._shapely_geometry + @staticmethod + def _boundary_angle_computation( + position: xr.DataArray, + reference_vector: xr.DataArray | np.ndarray, + how_to_compute_vector_to_region: Callable[ + [xr.DataArray], xr.DataArray + ], + angle_rotates: Literal["vec to ref", "ref to vec"] = "vec to ref", + in_degrees: bool = False, + ) -> xr.DataArray: + """Perform a boundary angle computation. + + Intended for internal use when conducting boundary angle computations, + to reduce code duplication. All boundary angle computations involve two + parts: + + - From some given spatial position data, compute the "vector towards + the region". This is typically the approach vector, but might also be + the normal vector if we are dealing with a segment or the plane + supported by a segment. + - Compute the signed angle between the "vector towards the region" and + some given reference vector, which may be constant or varying in time + (such as an animal's heading or forward vector). + + As such, we generalise the process into this internal method, and + provide more explicit wrappers to the user to make working with the + methods easier. + + Parameters + ---------- + position : xarray.DataArray + Spatial position data, that is passed to + ``how_to_compute_vector_to_region`` and used to compute the + "vector to the region". + reference_vector : xarray.DataArray | np.ndarray + Constant or time-varying vector to take signed angle with the + "vector to the region". + how_to_compute_vector_to_region : Callable + How to compute the "vector to the region" from ``position``. + angle_rotates : Literal["vec to ref", "ref to vec"] + Convention for the returned signed angle. ``"vec to ref"`` returns + the signed angle between the "vector to the region" and the + ``reference_vector``. ``"ref to vec"`` returns the opposite. + Default is ``"vec to ref"`` + in_degrees : bool + If ``True``, angles are returned in degrees. Otherwise angles are + returned in radians. Default ``False``. + + """ + if angle_rotates == "vec to ref": + ref_is_left_operand = False + elif angle_rotates == "ref to vec": + ref_is_left_operand = True + else: + raise ValueError(f"Unknown angle convention: {angle_rotates}") + + vec_to_segment = how_to_compute_vector_to_region(position) + + angles = compute_signed_angle_2d( + vec_to_segment, + reference_vector, + v_as_left_operand=ref_is_left_operand, + ) + if in_degrees: + angles = np.rad2deg(angles) + return angles + + @staticmethod + def _reassign_space_dim( + da: xr.DataArray, + old_dimension: Hashable, + ) -> xr.DataArray: + """Rename a computed dimension 'space' and assign coordinates. + + Intended for internal use when chaining ``DataArray``-broadcastable + operations together. In some instances, the outputs drop the spatial + coordinates, or the "space" axis is returned under a different name. + This needs to be corrected before further computations can be + performed. + + Parameters + ---------- + da : xarray.DataArray + ``DataArray`` lacking a "space" dimension, that is to be assigned. + old_dimension : Hashable + The dimension that should be renamed to "space", and reassigned + coordinates. + + """ + return da.rename({old_dimension: "space"}).assign_coords( + { + "space": ["x", "y"] + if len(da[old_dimension]) == 2 + else ["x", "y", "z"] + } + ) + def __init__( self, points: PointLikeList, @@ -200,3 +313,258 @@ def contains_point( point ) return point_is_inside + + @broadcastable_method(only_broadcastable_along="space") + def compute_distance_to( + self, point: ArrayLike, boundary_only: bool = False + ) -> float: + """Compute the distance from the region to a point. + + Parameters + ---------- + point : ArrayLike + Coordinates of a point, from which to find the nearest point in the + region defined by ``self``. + boundary_only : bool, optional + If ``True``, compute the distance from ``point`` to the boundary of + the region, rather than the closest point belonging to the region. + Default ``False``. + + Returns + ------- + float + Euclidean distance from the ``point`` to the (closest point on the) + region. + + """ + region_to_consider = ( + self.region.boundary if boundary_only else self.region + ) + return shapely.distance(region_to_consider, shapely.Point(point)) + + @broadcastable_method( + only_broadcastable_along="space", new_dimension_name="nearest point" + ) + def compute_nearest_point_to( + self, /, position: ArrayLike, boundary_only: bool = False + ) -> np.ndarray: + """Compute (one of) the nearest point(s) in the region to ``position``. + + If there are multiple equidistant points, one of them is returned. + + Parameters + ---------- + position : ArrayLike + Coordinates of a point, from which to find the nearest point in the + region. + boundary_only : bool, optional + If ``True``, compute the nearest point to ``position`` that is on + the boundary of ``self``. Default ``False``. + + Returns + ------- + np.ndarray + Coordinates of the point on ``self`` that is closest to + ``position``. + + """ + region_to_consider = ( + self.region.boundary if boundary_only else self.region + ) + # shortest_line returns a line from 1st arg to 2nd arg, + # therefore the point on self is the 0th coordinate + return np.array( + shapely.shortest_line( + region_to_consider, shapely.Point(position) + ).coords[0] + ) + + @broadcastable_method( + only_broadcastable_along="space", new_dimension_name="vector to" + ) + def compute_approach_vector( + self, + point: ArrayLike, + boundary_only: bool = False, + unit: bool = False, + ) -> np.ndarray: + """Compute the approach vector from a ``point`` to the region. + + The approach vector is defined as the vector directed from the + ``point`` provided, to the closest point that belongs to the region. + + Parameters + ---------- + point : ArrayLike + Coordinates of a point to compute the vector to (or from) the + region. + boundary_only : bool + If ``True``, the approach vector to the boundary of the region is + computed. Default ``False``. + unit : bool + If ``True``, the approach vector is returned normalised, otherwise + it is not normalised. Default is ``False``. + + Returns + ------- + np.ndarray + Approach vector from the point to the region. + + See Also + -------- + compute_allocentric_angle_to_nearest_point : + Relies on this method to compute the approach vector. + compute_egocentric_angle_to_nearest_point : + Relies on this method to compute the approach vector. + + """ + region_to_consider = ( + self.region.boundary if boundary_only else self.region + ) + + # "point to region" by virtue of order of arguments to shapely call + directed_line = shapely.shortest_line( + shapely.Point(point), region_to_consider + ) + + approach_vector = np.array(directed_line.coords[1]) - np.array( + directed_line.coords[0] + ) + if unit: + norm = np.sqrt(np.sum(approach_vector**2)) + # Cannot normalise the 0 vector + if not np.isclose(norm, 0.0): + approach_vector /= norm + return approach_vector + + def compute_allocentric_angle_to_nearest_point( + self, + position: xr.DataArray, + angle_rotates: Literal[ + "approach to ref", "ref to approach" + ] = "approach to ref", + boundary_only: bool = False, + in_degrees: bool = False, + reference_vector: np.ndarray | xr.DataArray = None, + ) -> float: + """Compute the allocentric angle to the nearest point in the region. + + With the term "allocentric", we indicate that we are measuring angles + with respect to a reference frame that is fixed relative to the + experimental/camera setup. By default, we assume this is the positive + x-axis of the coordinate system in which ``position`` is. + + The allocentric angle is the :func:`signed angle\ + ` between the approach + vector and a given reference vector. + + Parameters + ---------- + position : xarray.DataArray + ``DataArray`` of spatial positions. + angle_rotates : Literal["approach to ref", "ref to approach"] + Direction of the signed angle returned. Default is + ``"approach to ref"``. + boundary_only : bool + If ``True``, the allocentric angle to the closest boundary point of + the region is computed. Default ``False``. + in_degrees : bool + If ``True``, angles are returned in degrees. Otherwise angles are + returned in radians. Default ``False``. + reference_vector : ArrayLike | xr.DataArray + The reference vector to be used. Dimensions must be compatible with + the argument of the same name that is passed to + :func:`compute_signed_angle_2d`. Default ``(1., 0.)``. + + See Also + -------- + compute_approach_vector : + The method used to compute the approach vector. + compute_egocentric_angle_to_nearest_point : + Related class method for computing the egocentric angle to the + region. + movement.utils.vector.compute_signed_angle_2d : + The underlying function used to compute the signed angle between + the approach vector and the reference vector. + + """ + if reference_vector is None: + reference_vector = np.array([[1.0, 0.0]]) + + return self._boundary_angle_computation( + position=position, + reference_vector=reference_vector, + how_to_compute_vector_to_region=lambda p: self._reassign_space_dim( + self.compute_approach_vector( + p, boundary_only=boundary_only, unit=False + ), + "vector to", + ), + angle_rotates=angle_rotates.replace("approach", "vec"), # type: ignore + in_degrees=in_degrees, + ) + + def compute_egocentric_angle_to_nearest_point( + self, + direction: xr.DataArray, + position: xr.DataArray, + angle_rotates: Literal[ + "approach to direction", "direction to approach" + ] = "approach to direction", + boundary_only: bool = False, + in_degrees: bool = False, + ) -> xr.DataArray: + """Compute the egocentric angle to the nearest point in the region. + + With the term "egocentric", we indicate that we are measuring angles + with respect to a reference frame that is varying in time relative to + the experimental/camera setup. + + The egocentric angle is the signed angle between the approach vector + and a ``direction`` vector (examples include the forward vector of + a given individual, or the velocity vector of a given point). + + Parameters + ---------- + direction : xarray.DataArray + An array of vectors representing a given direction, + e.g., the forward vector(s). + position : xarray.DataArray + `DataArray` of spatial positions, considered the origin of the + ``direction`` vector. + angle_rotates : {"approach to direction", "direction to approach"} + Direction of the signed angle returned. Default is + ``"approach to direction"``. + boundary_only : bool + Passed to ``compute_approach_vector`` (see Notes). Default + ``False``. + in_degrees : bool + If ``True``, angles are returned in degrees. Otherwise angles are + returned in radians. Default ``False``. + + See Also + -------- + compute_allocentric_angle_to_nearest_point : + Related class method for computing the egocentric angle to the + region. + compute_approach_vector : + The method used to compute the approach vector. + movement.utils.vector.compute_signed_angle_2d : + The underlying function used to compute the signed angle between + the approach vector and the reference vector. + + """ + return self._boundary_angle_computation( + position=position, + reference_vector=direction, + how_to_compute_vector_to_region=lambda p: self._reassign_space_dim( + self.compute_approach_vector( + p, boundary_only=boundary_only, unit=False + ), + "vector to", + ), + angle_rotates=angle_rotates.replace("approach", "vec").replace( # type: ignore + "direction", "ref" + ), + in_degrees=in_degrees, + ) diff --git a/movement/roi/line.py b/movement/roi/line.py index 5359c830e..64553fb84 100644 --- a/movement/roi/line.py +++ b/movement/roi/line.py @@ -1,6 +1,16 @@ """1-dimensional lines of interest.""" -from movement.roi.base import BaseRegionOfInterest, PointLikeList +from typing import Literal + +import numpy as np +import xarray as xr +from numpy.typing import ArrayLike + +from movement.roi.base import ( + BaseRegionOfInterest, + PointLikeList, +) +from movement.utils.broadcasting import broadcastable_method class LineOfInterest(BaseRegionOfInterest): @@ -56,3 +66,88 @@ def __init__( """ super().__init__(points, dimensions=1, closed=loop, name=name) + + @broadcastable_method( + only_broadcastable_along="space", new_dimension_name="normal" + ) + def normal(self, on_same_side_as: ArrayLike = (0.0, 0.0)) -> np.ndarray: + """Compute the unit normal to this line. + + The unit normal is a vector perpendicular to the input line + whose norm is equal to 1. The direction of the normal vector + is not fully defined: the line divides the 2D plane in two + halves, and the normal could be pointing to either of the half-planes. + For example, a horizontal line divides the 2D plane in a + bottom and a top half-plane, and we can choose whether + the normal points "upwards" or "downwards". We use a sample + point to define the half-plane the normal vector points to. + + If this is a multi-segment line, the method raises an error. + + Parameters + ---------- + on_same_side_as : ArrayLike + A sample point in the (x,y) plane the normal is in. If multiple + points are given, one normal vector is returned for each point + given. By default, the origin is used. + + Raises + ------ + ValueError : When the normal is requested for a multi-segment geometry. + + """ + # A multi-segment geometry always has at least 3 coordinates. + if len(self.coords) > 2: + raise ValueError( + "Normal is not defined for multi-segment geometries." + ) + + on_same_side_as = np.array(on_same_side_as) + + parallel_to_line = np.array(self.region.coords[1]) - np.array( + self.region.coords[0] + ) + normal = np.array([parallel_to_line[1], -parallel_to_line[0]]) + normal /= np.sqrt(np.sum(normal**2)) + + if np.dot((on_same_side_as - self.region.coords[0]), normal) < 0: + normal *= -1.0 + return normal + + def compute_angle_to_normal( + self, + direction: xr.DataArray, + position: xr.DataArray, + angle_rotates: Literal[ + "direction to normal", "normal to direction" + ] = "normal to direction", + in_degrees: bool = False, + ) -> xr.DataArray: + """Compute the angle between the normal to the segment and a direction. + + Parameters + ---------- + direction : xarray.DataArray + An array of vectors representing a given direction, + e.g., the forward vector(s). + position : xr.DataArray + Spatial positions, considered the origin of the ``direction``. + angle_rotates : Literal["direction to normal", "normal to direction"] + Sign convention of the angle returned. Default is + ``"normal to direction"``. + in_degrees : bool + If ``True``, angles are returned in degrees. Otherwise angles are + returned in radians. Default ``False``. + + """ + return self._boundary_angle_computation( + position=position, + reference_vector=direction, + how_to_compute_vector_to_region=lambda p: self._reassign_space_dim( + -1.0 * self.normal(p), "normal" + ), + angle_rotates=angle_rotates.replace("direction", "ref").replace( # type: ignore + "normal", "vec" + ), + in_degrees=in_degrees, + ) diff --git a/movement/validators/arrays.py b/movement/validators/arrays.py index 24417c502..20e09df13 100644 --- a/movement/validators/arrays.py +++ b/movement/validators/arrays.py @@ -1,5 +1,7 @@ """Validators for data arrays.""" +from collections.abc import Hashable + import numpy as np import xarray as xr @@ -8,7 +10,7 @@ def validate_dims_coords( data: xr.DataArray, - required_dim_coords: dict[str, list[str]], + required_dim_coords: dict[str, list[str] | list[Hashable]], exact_coords: bool = False, ) -> None: """Validate dimensions and coordinates in a data array. @@ -23,7 +25,7 @@ def validate_dims_coords( ---------- data : xarray.DataArray The input data array to validate. - required_dim_coords : dict of {str: list of str} + required_dim_coords : dict of {str: list of str | list of Hashable} A dictionary mapping required dimensions to a list of required coordinate values along each dimension. exact_coords : bool, optional diff --git a/tests/fixtures/roi.py b/tests/fixtures/roi.py new file mode 100644 index 000000000..0b72f7d38 --- /dev/null +++ b/tests/fixtures/roi.py @@ -0,0 +1,56 @@ +import numpy as np +import pytest +import xarray as xr + +from movement.roi import LineOfInterest +from movement.roi.polygon import PolygonOfInterest + + +@pytest.fixture +def segment_of_y_equals_x() -> LineOfInterest: + """Line segment from (0,0) to (1,1).""" + return LineOfInterest([(0, 0), (1, 1)]) + + +@pytest.fixture() +def unit_square_pts() -> np.ndarray: + """Points that define the 4 corners of a unit-length square. + + The points have the lower-left corner positioned at (0,0). + """ + return np.array( + [ + [0.0, 0.0], + [1.0, 0.0], + [1.0, 1.0], + [0.0, 1.0], + ], + dtype=float, + ) + + +@pytest.fixture() +def unit_square_hole(unit_square_pts: np.ndarray) -> np.ndarray: + """Hole in the shape of a 0.5 side-length square centred on (0.5, 0.5).""" + return 0.25 + (unit_square_pts.copy() * 0.5) + + +@pytest.fixture +def unit_square(unit_square_pts: xr.DataArray) -> PolygonOfInterest: + """Square of unit side-length centred on (0.5, 0.5).""" + return PolygonOfInterest(unit_square_pts, name="Unit square") + + +@pytest.fixture +def unit_square_with_hole( + unit_square_pts: xr.DataArray, unit_square_hole: xr.DataArray +) -> PolygonOfInterest: + """Square of unit side length with an internal hole. + + The "outer" square is centred on (0.5, 0.5) and has side length 1. + The "inner" square, or hole, is centred on (0.5, 0.5) and has side length + 0.5. + """ + return PolygonOfInterest( + unit_square_pts, holes=[unit_square_hole], name="Unit square with hole" + ) diff --git a/tests/test_unit/conftest.py b/tests/test_unit/conftest.py new file mode 100644 index 000000000..f1cabffb9 --- /dev/null +++ b/tests/test_unit/conftest.py @@ -0,0 +1,62 @@ +from collections.abc import Callable + +import numpy as np +import pytest +import xarray as xr + + +@pytest.fixture(scope="session") +def push_into_range() -> Callable[ + [xr.DataArray | np.ndarray, float, float], xr.DataArray | np.ndarray +]: + """Return a function for wrapping angles. + + This is a factory fixture that returns a method for wrapping angles + into a user-specified range. + """ + + def _push_into_range( + numeric_values: xr.DataArray | np.ndarray, + lower: float = -180.0, + upper: float = 180.0, + ) -> xr.DataArray | np.ndarray: + """Coerce values into the range (lower, upper]. + + Primarily used to wrap returned angles into a particular range, + such as (-pi, pi]. + + The interval width is the value ``upper - lower``. + Each element in ``values`` that starts less than or equal to the + ``lower`` bound has multiples of the interval width added to it, + until the result lies in the desirable interval. + + Each element in ``values`` that starts greater than the ``upper`` + bound has multiples of the interval width subtracted from it, + until the result lies in the desired interval. + """ + translated_values = ( + numeric_values.values.copy() + if isinstance(numeric_values, xr.DataArray) + else numeric_values.copy() + ) + + interval_width = upper - lower + if interval_width <= 0: + raise ValueError( + f"Upper bound ({upper}) must be strictly " + f"greater than lower bound ({lower})" + ) + + while np.any( + (translated_values <= lower) | (translated_values > upper) + ): + translated_values[translated_values <= lower] += interval_width + translated_values[translated_values > upper] -= interval_width + + if isinstance(numeric_values, xr.DataArray): + translated_values = numeric_values.copy( + deep=True, data=translated_values + ) + return translated_values + + return _push_into_range diff --git a/tests/test_unit/test_kinematics.py b/tests/test_unit/test_kinematics.py index 57fd2479a..18200cd3e 100644 --- a/tests/test_unit/test_kinematics.py +++ b/tests/test_unit/test_kinematics.py @@ -707,48 +707,6 @@ class TestForwardVectorAngle: y_axis = np.array([0.0, 1.0]) sqrt_2 = np.sqrt(2.0) - @staticmethod - def push_into_range( - numeric_values: xr.DataArray | np.ndarray, - lower: float = -180.0, - upper: float = 180.0, - ) -> xr.DataArray | np.ndarray: - """Translate values into the range (lower, upper]. - - The interval width is the value ``upper - lower``. - Each ``v`` in ``values`` that starts less than or equal to the - ``lower`` bound has multiples of the interval width added to it, - until the result lies in the desirable interval. - - Each ``v`` in ``values`` that starts greater than the ``upper`` - bound has multiples of the interval width subtracted from it, - until the result lies in the desired interval. - """ - translated_values = ( - numeric_values.values.copy() - if isinstance(numeric_values, xr.DataArray) - else numeric_values.copy() - ) - - interval_width = upper - lower - if interval_width <= 0: - raise ValueError( - f"Upper bound ({upper}) must be strictly " - f"greater than lower bound ({lower})" - ) - - while np.any( - (translated_values <= lower) | (translated_values > upper) - ): - translated_values[translated_values <= lower] += interval_width - translated_values[translated_values > upper] -= interval_width - - if isinstance(numeric_values, xr.DataArray): - translated_values = numeric_values.copy( - deep=True, data=translated_values - ) - return translated_values - @pytest.fixture def spinning_on_the_spot(self) -> xr.DataArray: """Simulate data for an individual's head spinning on the spot. @@ -794,6 +752,7 @@ def spinning_on_the_spot(self) -> xr.DataArray: ) def test_antisymmetry_properties( self, + push_into_range, spinning_on_the_spot: xr.DataArray, swap_left_right: bool, swap_camera_view: bool, @@ -843,16 +802,16 @@ def test_antisymmetry_properties( expected_orientations = without_orientations_swapped.copy(deep=True) if swap_left_right: - expected_orientations = self.push_into_range( + expected_orientations = push_into_range( expected_orientations + 180.0 ) if swap_camera_view: - expected_orientations = self.push_into_range( + expected_orientations = push_into_range( expected_orientations + 180.0 ) if swap_angle_rotation: expected_orientations *= -1.0 - expected_orientations = self.push_into_range(expected_orientations) + expected_orientations = push_into_range(expected_orientations) xr.testing.assert_allclose( with_orientations_swapped, expected_orientations diff --git a/tests/test_unit/test_roi/conftest.py b/tests/test_unit/test_roi/conftest.py deleted file mode 100644 index 8ba15cb94..000000000 --- a/tests/test_unit/test_roi/conftest.py +++ /dev/null @@ -1,21 +0,0 @@ -import numpy as np -import pytest - - -@pytest.fixture() -def unit_square_pts() -> np.ndarray: - return np.array( - [ - [0.0, 0.0], - [1.0, 0.0], - [1.0, 1.0], - [0.0, 1.0], - ], - dtype=float, - ) - - -@pytest.fixture() -def unit_square_hole(unit_square_pts: np.ndarray) -> np.ndarray: - """Hole in the shape of a 0.5 side-length square centred on (0.5, 0.5).""" - return 0.25 + (unit_square_pts.copy() * 0.5) diff --git a/tests/test_unit/test_roi/test_angles.py b/tests/test_unit/test_roi/test_angles.py new file mode 100644 index 000000000..3ee41e8a8 --- /dev/null +++ b/tests/test_unit/test_roi/test_angles.py @@ -0,0 +1,391 @@ +import re +from collections.abc import Iterable +from typing import Any + +import numpy as np +import pytest +import xarray as xr + +from movement.kinematics import compute_forward_vector +from movement.roi import LineOfInterest +from movement.roi.base import BaseRegionOfInterest + + +def sample_position_array() -> xr.DataArray: + """Return a simulated position array to test the egocentric angle. + + The data has time, space, and keypoints dimensions. + + The keypoints are left, right, midpt (midpoint), and wild. + The midpt is the mean of the left and right keypoints. + The wild keypoint is a point in the plane distinct from the midpt, and is + used to test that the function respects the origin of the forward vector + that the user provides (even if it is physically non-nonsensical). + + time 1: + left @ (1.25, 0.), right @ (1., -0.25), wild @ (-0.25, -0.25) + + Fwd vector is (1, -1). + time 2: + left @ (-0.25, 0.5), right @ (-0.25, 0.25), wild @ (-0.5, 0.375) + + Fwd vector is (-1, 0). + time 3: + left @ (0.375, 0.375), right @ (0.5, 0.375), wild @ (1.25, 1.25) + + Fwd vector is (0, -1). + time 4: + left @ (1., -0.25), right @ (1.25, 0.), wild @ (-0.25, -0.25) + + This is time 1 but with left and right swapped. + Fwd vector is (-1, 1). + time 5: + left @ (0.25, 0.5), right @ (0.375, 0.25), wild @ (0.5, 0.65) + + Fwd vector is (-2, -1). + acos(2/sqrt(5)) is the expected angle for the midpt. + acos(1/sqrt(5)) should be that for the wild point IF going to the + boundary. + """ + points = np.zeros(shape=(5, 2, 4)) + points[:, :, 0] = [ + [1.25, 0.0], + [-0.25, 0.5], + [0.375, 0.375], + [1.0, -0.25], + [0.25, 0.5], + ] + points[:, :, 1] = [ + [1.0, -0.25], + [-0.25, 0.25], + [0.5, 0.375], + [1.25, 0.0], + [0.375, 0.25], + ] + points[:, :, 3] = [ + [-0.25, 1.25], + [-0.5, 0.375], + [1.25, 1.25], + [-0.25, -0.25], + [0.5, 0.65], + ] + points[:, :, 2] = np.mean(points[:, :, 0:2], axis=2) + return xr.DataArray( + data=points, + dims=["time", "space", "keypoints"], + coords={ + "space": ["x", "y"], + "keypoints": ["left", "right", "midpt", "wild"], + }, + ) + + +@pytest.mark.parametrize( + ["region", "fn_args", "fn_kwargs", "expected_output", "egocentric"], + [ + pytest.param( + "unit_square_with_hole", + [ + compute_forward_vector( + sample_position_array(), "left", "right" + ), + sample_position_array().sel(keypoints="midpt", drop=True), + ], + {"angle_rotates": "elephant to region"}, + ValueError("Unknown angle convention: elephant to region"), + True, + id="[Egocentric] Unknown angle convention", + ), + pytest.param( + "unit_square_with_hole", + [sample_position_array()], + {"angle_rotates": "elephant to region"}, + ValueError("Unknown angle convention: elephant to region"), + False, + id="[Allocentric] Unknown angle convention", + ), + pytest.param( + "unit_square_with_hole", + [ + compute_forward_vector( + sample_position_array(), "left", "right" + ), + sample_position_array().sel(keypoints="midpt", drop=True), + ], + {"in_degrees": True}, + np.array( + [ + 0.0, + 180.0, + 0.0, + 180.0, + np.rad2deg(np.arccos(2.0 / np.sqrt(5.0))), + ] + ), + True, + id="[Egocentric] Default args", + ), + pytest.param( + "unit_square_with_hole", + [ + compute_forward_vector( + sample_position_array(), "left", "right" + ), + sample_position_array().sel(keypoints="wild", drop=True), + ], + {"in_degrees": True}, + np.array( + [ + 180.0, + 180.0, + 45.0, + -90.0, + np.rad2deg(np.pi / 2.0 + np.arcsin(1.0 / np.sqrt(5.0))), + ] + ), + True, + id="[Egocentric] Non-default position", + ), + pytest.param( + "unit_square", + [ + compute_forward_vector( + sample_position_array(), "left", "right" + ), + sample_position_array().sel(keypoints="midpt", drop=True), + ], + {"in_degrees": True}, + np.array( + [ + 0.0, + 180.0, + float("nan"), + 180.0, + float("nan"), + ] + ), + True, + id="[Egocentric] 0-approach vectors (nan returns)", + ), + pytest.param( + "unit_square", + [ + compute_forward_vector( + sample_position_array(), "left", "right" + ), + sample_position_array().sel(keypoints="midpt", drop=True), + ], + { + "boundary_only": True, + "in_degrees": True, + }, + np.array( + [ + 0.0, + 180.0, + 0.0, + 180.0, + np.rad2deg(np.arccos(2.0 / np.sqrt(5.0))), + ] + ), + True, + id="[Egocentric] Force boundary calculations", + ), + pytest.param( + "unit_square_with_hole", + [sample_position_array().sel(keypoints="midpt", drop=True)], + {}, + np.deg2rad( + np.array( + [ + -135.0, + 0.0, + 90.0, + -135.0, + 180.0, + ] + ) + ), + False, + id="[Allocentric] Default args", + ), + pytest.param( + "unit_square", + [sample_position_array().sel(keypoints="midpt", drop=True)], + {}, + np.deg2rad( + np.array( + [ + -135.0, + 0.0, + float("nan"), + -135.0, + float("nan"), + ] + ) + ), + False, + id="[Allocentric] 0-approach vectors", + ), + pytest.param( + "unit_square", + [sample_position_array().sel(keypoints="midpt", drop=True)], + {"boundary_only": True}, + np.deg2rad( + np.array( + [ + -135.0, + 0.0, + 90.0, + -135.0, + 180.0, + ] + ) + ), + False, + id="[Allocentric] Force boundary calculation", + ), + ], +) +def test_ego_and_allocentric_angle_to_region( + push_into_range, + region: BaseRegionOfInterest, + fn_args: Iterable[Any], + fn_kwargs: dict[str, Any], + expected_output: xr.DataArray | Exception, + egocentric: bool, + request, +) -> None: + """Test computation of the egocentric and allocentric angle. + + Note, we only test functionality explicitly introduced in this method. + Input arguments that are just handed to other functions (i.e., + ``in_degrees``), are not explicitly tested here. + + The ``angle_rotates`` argument is tested in all cases (signs should be + reversed when toggling the argument). + """ + if isinstance(region, str): + region = request.getfixturevalue(region) + if isinstance(expected_output, np.ndarray): + expected_output = xr.DataArray(data=expected_output, dims=["time"]) + + if egocentric: + which_method = "compute_egocentric_angle_to_nearest_point" + other_vector_name = "direction" + else: + which_method = "compute_allocentric_angle_to_nearest_point" + other_vector_name = "ref" + method = getattr(region, which_method) + + if isinstance(expected_output, Exception): + with pytest.raises( + type(expected_output), match=re.escape(str(expected_output)) + ): + method(*fn_args, **fn_kwargs) + else: + angles = method(*fn_args, **fn_kwargs) + xr.testing.assert_allclose(angles, expected_output) + + # Check reversal of the angle convention, + # which should just reverse the sign of the angles, + # except for 180-degrees -> 180-degrees. + if fn_kwargs.get("in_degrees", False): + lower = -180.0 + upper = 180.0 + else: + lower = -np.pi + upper = np.pi + if ( + fn_kwargs.get("angle_rotates", f"approach to {other_vector_name}") + == f"approach to {other_vector_name}" + ): + fn_kwargs["angle_rotates"] = f"{other_vector_name} to approach" + else: + fn_kwargs["angle_rotates"] = f"approach to {other_vector_name}" + + reverse_angles = push_into_range( + method(*fn_args, **fn_kwargs), lower=lower, upper=upper + ) + xr.testing.assert_allclose( + angles, push_into_range(-reverse_angles, lower=lower, upper=upper) + ) + + +@pytest.fixture +def points_around_segment() -> xr.DataArray: + """Sample points on either side of the segment y=x. + + Data has (time, space, keypoints) dimensions, shape (, 2, 2). + + Keypoints are "left" and "right". + + time 1: + left @ (0., 1.), right @ (0.05, 0.95). + Fwd vector is (-1, -1). + time 2: + left @ (1., 0.), right @ (0.95, 0.05). + Fwd vector is (-1, -1). + time 3: + left @ (1., 2.), right @ (1.05, 1.95). + Fwd vector is (-1, -1). + The egocentric angle will differ when using this point. + """ + points = np.zeros(shape=(3, 2, 2)) + points[:, :, 0] = [ + [0.0, 1.0], + [1.0, 0.0], + [1.0, 2.0], + ] + points[:, :, 1] = [ + [0.05, 0.95], + [0.95, 0.05], + [1.05, 1.95], + ] + return xr.DataArray( + data=points, + dims=["time", "space", "keypoints"], + coords={ + "space": ["x", "y"], + "keypoints": ["left", "right"], + }, + ) + + +def test_angle_to_normal( + segment_of_y_equals_x: LineOfInterest, + points_around_segment: xr.DataArray, +) -> None: + """Test the angle_to_support_plane method. + + This method checks two things: + + - The compute_angle_to_normal method returns the correct angle, and + - The method agrees with the egocentric angle computation, in the cases + that the two calculations should return the same value (IE when the + approach vector is the normal to the segment). And that the + returned angles are different otherwise. + """ + expected_output = xr.DataArray( + data=np.deg2rad([-90.0, -90.0, -90.0]), dims=["time"] + ) + should_be_same_as_egocentric = expected_output.copy( + data=[True, True, False], deep=True + ) + + fwd_vector = compute_forward_vector(points_around_segment, "left", "right") + positions = points_around_segment.mean(dim="keypoints") + angles_to_support = segment_of_y_equals_x.compute_angle_to_normal( + fwd_vector, positions + ) + xr.testing.assert_allclose(expected_output, angles_to_support) + + egocentric_angles = ( + segment_of_y_equals_x.compute_egocentric_angle_to_nearest_point( + fwd_vector, positions + ) + ) + values_are_close = egocentric_angles.copy( + data=np.isclose(egocentric_angles, angles_to_support), deep=True + ) + xr.testing.assert_equal(should_be_same_as_egocentric, values_are_close) diff --git a/tests/test_unit/test_roi/test_nearest_points.py b/tests/test_unit/test_roi/test_nearest_points.py new file mode 100644 index 000000000..35ce28495 --- /dev/null +++ b/tests/test_unit/test_roi/test_nearest_points.py @@ -0,0 +1,365 @@ +import re +from typing import Any + +import numpy as np +import pytest +import xarray as xr + +from movement.roi.base import BaseRegionOfInterest +from movement.roi.line import LineOfInterest + + +@pytest.fixture +def sample_target_points() -> dict[str, np.ndarray]: + """Sample 2D trajectory data.""" + return xr.DataArray( + np.array( + [ + [-0.5, 0.50], + [0.00, 0.50], + [0.40, 0.45], + [2.00, 1.00], + [0.40, 0.75], + [0.95, 0.90], + [0.80, 0.76], + ] + ), + dims=["time", "space"], + coords={"space": ["x", "y"]}, + ) + + +@pytest.fixture() +def unit_line_in_x() -> LineOfInterest: + return LineOfInterest([[0.0, 0.0], [1.0, 0.0]]) + + +@pytest.mark.parametrize( + ["region", "fn_kwargs", "expected_distances"], + [ + pytest.param( + "unit_square_with_hole", + {"boundary_only": True}, + np.array( + [ + 0.5, + 0.0, + 0.15, + 1.0, + 0.0, + 0.05, + np.sqrt((1.0 / 20.0) ** 2 + (1.0 / 100.0) ** 2), + ] + ), + id="Unit square w/ hole, boundary", + ), + pytest.param( + "unit_square_with_hole", + {}, + np.array([0.5, 0.0, 0.15, 1.0, 0.0, 0.0, 0.0]), + id="Unit square w/ hole, whole region", + ), + pytest.param( + "unit_line_in_x", + {}, + np.array( + [0.5 * np.sqrt(2.0), 0.5, 0.45, np.sqrt(2.0), 0.75, 0.9, 0.76] + ), + id="Unit line in x", + ), + pytest.param( + "unit_line_in_x", + {"boundary_only": True}, + np.array( + [ + 0.5 * np.sqrt(2.0), + 0.5, + np.sqrt(0.4**2 + 0.45**2), + np.sqrt(2.0), + np.sqrt(0.4**2 + 0.75**2), + np.sqrt(0.05**2 + 0.9**2), + np.sqrt(0.2**2 + 0.76**2), + ] + ), + id="Unit line in x, endpoints only", + ), + ], +) +def test_distance_point_to_region( + region: BaseRegionOfInterest, + sample_target_points: xr.DataArray, + fn_kwargs: dict[str, Any], + expected_distances: xr.DataArray, + request, +) -> None: + if isinstance(region, str): + region = request.getfixturevalue(region) + if isinstance(expected_distances, np.ndarray): + expected_distances = xr.DataArray( + data=expected_distances, dims=["time"] + ) + + computed_distances = region.compute_distance_to( + sample_target_points, **fn_kwargs + ) + + xr.testing.assert_allclose(computed_distances, expected_distances) + + +@pytest.mark.parametrize( + ["region", "other_fn_args", "expected_output"], + [ + pytest.param( + "unit_square", + {"boundary_only": True}, + np.array( + [ + [0.00, 0.50], + [0.00, 0.50], + [0.00, 0.45], + [1.00, 1.00], + [0.40, 1.00], + [1.00, 0.90], + [1.00, 0.76], + ] + ), + id="Unit square, boundary only", + ), + pytest.param( + "unit_square", + {}, + np.array( + [ + [0.00, 0.50], + [0.00, 0.50], + [0.40, 0.45], + [1.00, 1.00], + [0.40, 0.75], + [0.95, 0.90], + [0.80, 0.76], + ] + ), + id="Unit square, whole region", + ), + pytest.param( + "unit_square_with_hole", + {"boundary_only": True}, + np.array( + [ + [0.00, 0.50], + [0.00, 0.50], + [0.25, 0.45], + [1.00, 1.00], + [0.40, 0.75], + [1.00, 0.90], + [0.75, 0.75], + ] + ), + id="Unit square w/ hole, boundary only", + ), + pytest.param( + "unit_square_with_hole", + {}, + np.array( + [ + [0.00, 0.50], + [0.00, 0.50], + [0.25, 0.45], + [1.00, 1.00], + [0.40, 0.75], + [0.95, 0.90], + [0.80, 0.76], + ] + ), + id="Unit square w/ hole, whole region", + ), + pytest.param( + "unit_line_in_x", + {}, + np.array( + [ + [0.00, 0.00], + [0.00, 0.00], + [0.40, 0.00], + [1.00, 0.00], + [0.40, 0.00], + [0.95, 0.00], + [0.80, 0.00], + ] + ), + id="Line, whole region", + ), + pytest.param( + "unit_line_in_x", + {"boundary_only": True}, + np.array( + [ + [0.00, 0.00], + [0.00, 0.00], + [0.00, 0.00], + [1.00, 0.00], + [0.00, 0.00], + [1.00, 0.00], + [1.00, 0.00], + ] + ), + id="Line, boundary only", + ), + ], +) +def test_nearest_point_to( + region: BaseRegionOfInterest, + sample_target_points: xr.DataArray, + other_fn_args: dict[str, Any], + expected_output: xr.DataArray, + request, +) -> None: + if isinstance(region, str): + region = request.getfixturevalue(region) + if isinstance(expected_output, str): + expected_output = request.get(expected_output) + elif isinstance(expected_output, np.ndarray): + expected_output = xr.DataArray( + expected_output, + dims=["time", "nearest point"], + ) + + nearest_points = region.compute_nearest_point_to( + sample_target_points, **other_fn_args + ) + + xr.testing.assert_allclose(nearest_points, expected_output) + + +@pytest.mark.parametrize( + ["region", "position", "fn_kwargs", "possible_nearest_points"], + [ + pytest.param( + "unit_square", + [0.5, 0.5], + {"boundary_only": True}, + [ + np.array([0.0, 0.5]), + np.array([0.5, 0.0]), + np.array([1.0, 0.5]), + np.array([0.5, 1.0]), + ], + id="Centre of the unit square", + ), + pytest.param( + "unit_line_in_x", + [0.5, 0.0], + {"boundary_only": True}, + [ + np.array([0.0, 0.0]), + np.array([1.0, 0.0]), + ], + id="Boundary of a line", + ), + ], +) +def test_nearest_point_to_tie_breaks( + region: BaseRegionOfInterest, + position: np.ndarray, + fn_kwargs: dict[str, Any], + possible_nearest_points: list[np.ndarray], + request, +) -> None: + """Check behaviour when points are tied for nearest. + + This can only occur when we have a Polygonal region, or a multi-line 1D + region. In this case, there may be multiple points in the region of + interest that are tied for closest. ``shapely`` does not actually document + how it breaks ties here, but we can at least check that it identifies one + of the possible correct points. + """ + if isinstance(region, str): + region = request.getfixturevalue(region) + if not isinstance(position, np.ndarray | xr.DataArray): + position = np.array(position) + + nearest_point_found = region.compute_nearest_point_to( + position, **fn_kwargs + ) + + sq_dist_to_nearest_pt = np.sum((nearest_point_found - position) ** 2) + + n_matches = 0 + for possibility in possible_nearest_points: + # All possibilities should be approximately the same distance away + # from the position + assert np.isclose( + np.sum((possibility - position) ** 2), sq_dist_to_nearest_pt + ) + # We should match at least one possibility, + # track to see if we do. + if np.isclose(nearest_point_found, possibility).all(): + n_matches += 1 + assert n_matches == 1 + + +@pytest.mark.parametrize( + ["region", "point", "other_fn_args", "expected_output"], + [ + pytest.param( + "unit_square", + (-0.5, 0.0), + {"unit": True}, + np.array([1.0, 0.0]), + id="(-0.5, 0.0) -> unit square", + ), + pytest.param( + LineOfInterest([(0.0, 0.0), (1.0, 0.0)]), + (0.1, 0.5), + {"unit": True}, + np.array([0.0, -1.0]), + id="(0.1, 0.5) -> +ve x ray", + ), + pytest.param( + "unit_square", + (-0.5, 0.0), + {"unit": False}, + np.array([0.5, 0.0]), + id="Don't normalise output", + ), + pytest.param( + "unit_square", + (0.5, 0.5), + {"unit": True}, + np.array([0.0, 0.0]), + id="Interior point returns 0 vector", + ), + pytest.param( + "unit_square", + (0.25, 0.35), + {"boundary_only": True, "unit": True}, + np.array([-1.0, 0.0]), + id="Boundary, polygon", + ), + pytest.param( + LineOfInterest([(0.0, 0.0), (1.0, 0.0)]), + (0.1, 0.5), + {"boundary_only": True, "unit": True}, + np.array([-0.1, -0.5]) / np.sqrt(0.1**2 + 0.5**2), + id="Boundary, line", + ), + ], +) +def test_approach_vector( + region: BaseRegionOfInterest, + point: xr.DataArray, + other_fn_args: dict[str, Any], + expected_output: np.ndarray | Exception, + request, +) -> None: + if isinstance(region, str): + region = request.getfixturevalue(region) + + if isinstance(expected_output, Exception): + with pytest.raises( + type(expected_output), match=re.escape(str(expected_output)) + ): + region.compute_approach_vector(point, **other_fn_args) + else: + vector_to = region.compute_approach_vector(point, **other_fn_args) + assert np.allclose(vector_to, expected_output) diff --git a/tests/test_unit/test_roi/test_normal.py b/tests/test_unit/test_roi/test_normal.py new file mode 100644 index 000000000..ec789a75c --- /dev/null +++ b/tests/test_unit/test_roi/test_normal.py @@ -0,0 +1,63 @@ +import re + +import numpy as np +import pytest +from numpy.typing import ArrayLike + +from movement.roi import LineOfInterest + +SQRT_2 = np.sqrt(2.0) + + +@pytest.mark.parametrize( + ["segment", "point", "expected_normal"], + [ + pytest.param( + "segment_of_y_equals_x", + (0.0, 1.0), + (-1.0 / SQRT_2, 1.0 / SQRT_2), + id="Normal pointing to half-plane with point (0, 1)", + ), + pytest.param( + "segment_of_y_equals_x", + (1.0, 0.0), + (1.0 / SQRT_2, -1.0 / SQRT_2), + id="Normal pointing to half-plane with point (1, 0)", + ), + pytest.param( + LineOfInterest([(0.5, 0.5), (1.0, 1.0)]), + (1.0, 0.0), + (1.0 / SQRT_2, -1.0 / SQRT_2), + id="Segment does not start at origin", + ), + pytest.param( + "segment_of_y_equals_x", + (1.0, 2.0), + (-1.0 / SQRT_2, 1.0 / SQRT_2), + id="Necessary to extend segment to compute normal.", + ), + pytest.param( + LineOfInterest([(0.0, 0.0), (1.0, 0.0), (2.0, 0.0)]), + (0.5, 0.5), + ValueError("Normal is not defined for multi-segment geometries."), + id="Multi-segment lines do not have normals.", + ), + ], +) +def test_normal( + segment: LineOfInterest, + point: ArrayLike, + expected_normal: np.ndarray | Exception, + request, +) -> None: + if isinstance(segment, str): + segment = request.getfixturevalue(segment) + + if isinstance(expected_normal, Exception): + with pytest.raises( + type(expected_normal), match=re.escape(str(expected_normal)) + ): + segment.normal(point) + else: + computed_normal = segment.normal(point) + assert np.allclose(computed_normal, expected_normal)