diff --git a/docs/_static/img/cumsum_animal_geometries_grid.webp b/docs/_static/img/cumsum_animal_geometries_grid.webp new file mode 100644 index 000000000..0d0f242cd Binary files /dev/null and b/docs/_static/img/cumsum_animal_geometries_grid.webp differ diff --git a/docs/_static/img/path_aspect_ratio.webp b/docs/_static/img/path_aspect_ratio.webp new file mode 100644 index 000000000..5464b354f Binary files /dev/null and b/docs/_static/img/path_aspect_ratio.webp differ diff --git a/docs/tables/sliding_linearity_index_cuda.csv b/docs/tables/sliding_linearity_index_cuda.csv new file mode 100644 index 000000000..70b7d99d9 --- /dev/null +++ b/docs/tables/sliding_linearity_index_cuda.csv @@ -0,0 +1,13 @@ +FRAMES (MILLIONS),GPU TIME (S),GPU TIME (STEV) +2,0.03031,0.002 +4,0.05208,0.004 +8,0.08882,0.0056 +16,0.17612,0.0108 +32,0.37699,0.01587 +64,0.70194,0.01848 +128,1.36778,0.06145 +256,3.09965,0.21796 +512,17.9707,10.5061 +NVIDIA GeForce RTX 4070,, +time window = 2.5s @ 30 FPS,, +3 ITERATIONS,, diff --git a/simba/data_processors/cuda/timeseries.py b/simba/data_processors/cuda/timeseries.py index 5d03ce6cc..434a7f74f 100644 --- a/simba/data_processors/cuda/timeseries.py +++ b/simba/data_processors/cuda/timeseries.py @@ -3,10 +3,10 @@ import numpy as np from numba import cuda -from simba.data_processors.cuda.utils import _euclid_dist -from simba.data_processors.cuda.utils import _cuda_mean, _cuda_std +from simba.data_processors.cuda.utils import _euclid_dist, _cuda_mean, _cuda_std, _cuda_available from simba.utils.checks import check_float, check_valid_array from simba.utils.enums import Formats +from simba.utils.errors import SimBAGPUError THREADS_PER_BLOCK = 1024 @@ -158,6 +158,9 @@ def sliding_spatial_density_cuda(x: np.ndarray, :align: center :header-rows: 1 + .. seealso:: + :func:`simba.mixins.timeseries_features_mixin.TimeseriesFeatureMixin.spatial_density`, :func:`simba.mixins.timeseries_features_mixin.TimeseriesFeatureMixin.sliding_spatial_density` + :param np.ndarray x: A 2D array of shape (N, 2), where N is the number of points and each point has two spatial coordinates (x, y). The array represents the trajectory path of points in a 2D space (e.g., x and y positions in space). :param float radius: The radius (in millimeters) within which to count neighboring points around each trajectory point. Defines the area of interest around each point. :param float pixels_per_mm: The scaling factor that converts the physical radius (in millimeters) to pixel units for spatial density calculations. @@ -189,3 +192,69 @@ def sliding_spatial_density_cuda(x: np.ndarray, results = cuda.device_array(shape=x.shape[0], dtype=np.float16) _sliding_spatial_density_kernel[bpg, THREADS_PER_BLOCK](x_dev, time_window_frames_dev, radius_dev, results) return results.copy_to_host() + + +@cuda.jit() +def _sliding_linearity_index_kernel(x, time_frms, results): + r = cuda.grid(1) + if r >= x.shape[0] or r < 0: + return + l = int(r - time_frms[0]) + if l < 0 or l >= r: + return + sample_x = x[l:r] + straight_line_distance = _euclid_dist(sample_x[0], sample_x[-1]) + path_dist = 0 + for i in range(1, sample_x.shape[0]): + path_dist += _euclid_dist(sample_x[i - 1], sample_x[i]) + if path_dist == 0: + results[r] = 0.0 + else: + results[r] = straight_line_distance / path_dist + + +def sliding_linearity_index_cuda(x: np.ndarray, + window_size: float, + sample_rate: float) -> np.ndarray: + """ + + Calculates the straightness (linearity) index of a path using CUDA acceleration. + + The output is a value between 0 and 1, where 1 indicates a perfectly straight path. + + .. csv-table:: + :header: EXPECTED RUNTIMES + :file: ../../../docs/tables/sliding_linearity_index_cuda.csv + :widths: 10, 45, 45 + :align: center + :header-rows: 1 + + .. seealso:: + :func:`simba.mixins.timeseries_features_mixin.TimeseriesFeatureMixin.sliding_linearity_index`, :func:`simba.mixins.timeseries_features_mixin.TimeseriesFeatureMixin.linearity_index` + + + :param np.ndarray x: An (N, M) array representing the path, where N is the number of points and M is the number of spatial dimensions (e.g., 2 for 2D or 3 for 3D). Each row represents the coordinates of a point along the path. + :param float x: The size of the sliding window in seconds. This defines the time window over which the linearity index is calculated. The window size should be specified in seconds. + :param float sample_rate: The sample rate in Hz (samples per second), which is used to convert the window size from seconds to frames. + :return: A 1D array of length N, where each element represents the linearity index of the path within a sliding window. The value is a ratio between the straight-line distance and the actual path length for each window. Values range from 0 to 1, with 1 indicating a perfectly straight path. + :rtype: np.ndarray + + :example: + >>> x = np.random.randint(0, 500, (100, 2)).astype(np.float32) + >>> q = sliding_linearity_index_cuda(x=x, window_size=2, sample_rate=30) + """ + + check_valid_array(data=x, source=f'{sliding_linearity_index_cuda.__name__} x', accepted_ndims=(2,), + accepted_axis_1_shape=[2, ], accepted_dtypes=Formats.NUMERIC_DTYPES.value) + check_float(name=f'{sliding_linearity_index_cuda.__name__} window_size', value=window_size) + check_float(name=f'{sliding_linearity_index_cuda.__name__} sample_rate', value=sample_rate) + x = np.ascontiguousarray(x) + time_window_frames = np.array([max(1.0, np.ceil(window_size * sample_rate))]) + if not _cuda_available()[0]: + SimBAGPUError(msg='No GPU found', source=sliding_linearity_index_cuda.__name__) + x_dev = cuda.to_device(x) + time_window_frames_dev = cuda.to_device(time_window_frames) + bpg = (x.shape[0] + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK + results = cuda.device_array(shape=x.shape[0], dtype=np.float16) + _sliding_linearity_index_kernel[bpg, THREADS_PER_BLOCK](x_dev, time_window_frames_dev, results) + return results.copy_to_host() diff --git a/simba/mixins/geometry_mixin.py b/simba/mixins/geometry_mixin.py index 3ce6d768d..28b87650e 100644 --- a/simba/mixins/geometry_mixin.py +++ b/simba/mixins/geometry_mixin.py @@ -3071,7 +3071,7 @@ def cumsum_coord_geometries(self, Compute the cumulative time a body-part has spent inside a grid of geometries using multiprocessing. :param np.ndarray data: Input data array where rows represent frames and columns represent body-part x and y coordinates. - :param Dict[Tuple[int, int], Polygon] geometries: Dictionary of polygons representing spatial regions. Created by ``GeometryMixin.bucket_img_into_squares``. + :param Dict[Tuple[int, int], Polygon] geometries: Dictionary of polygons representing spatial regions. E.g., created by :func:`simba.mixins.geometry_mixin.GeometryMixin.bucket_img_into_grid_square` or :func:`simba.mixins.geometry_mixin.GeometryMixin.bucket_img_into_grid_hexagon`. :param Optional[int] fps: Frames per second (fps) for time normalization. If None, cumulative sum of frame count is returned. :example: @@ -3154,7 +3154,7 @@ def cumsum_bool_geometries( E.g., compute the cumulative time of classified events within spatial locations at all time-points of the video. :param np.ndarray data: Array containing spatial data with shape (n, 2). E.g., 2D-array with body-part coordinates. - :param Dict[Tuple[int, int], Polygon] geometries: Dictionary of polygons representing spatial regions. Created by ``GeometryMixin.bucket_img_into_squares``. + :param Dict[Tuple[int, int], Polygon] geometries: Dictionary of polygons representing spatial regions. E.g., created by :func:`simba.mixins.geometry_mixin.GeometryMixin.bucket_img_into_grid_square` or :func:`simba.mixins.geometry_mixin.GeometryMixin.bucket_img_into_grid_hexagon`. :param np.ndarray bool_data: Boolean array with shape (data.shape[0],) or (data.shape[0], 1) indicating the presence or absence in each frame. :param Optional[float] fps: Frames per second. If provided, the result is normalized by the frame rate. :param Optional[float] core_cnt: Number of CPU cores to use for parallel processing. Default is -1, which means using all available cores. @@ -3246,54 +3246,45 @@ def _cumsum_animal_geometries_grid_helper( results[k[0], k[1]] = 1 return results - def cumsum_animal_geometries_grid( - self, - data: List[Polygon], - grid: Dict[Tuple[int, int], Polygon], - fps: Optional[int] = None, - core_cnt: Optional[int] = -1, - verbose: Optional[bool] = True, - ): + def cumsum_animal_geometries_grid(self, + data: List[Polygon], + grid: Dict[Tuple[int, int], Polygon], + fps: Optional[int] = None, + core_cnt: Optional[int] = -1, + verbose: Optional[bool] = True) -> np.ndarray: + """ + .. image:: _static/img/cumsum_animal_geometries_grid.webp + :width: 400 + :align: center + + :param List[Polygon] data: List of polygons where every index represent a frame and every value the animal convex hull + :param Dict[Tuple[int, int], Polygon] grid: Dictionary of polygons representing spatial regions. E.g., created by :func:`simba.mixins.geometry_mixin.GeometryMixin.bucket_img_into_grid_square` or :func:`simba.mixins.geometry_mixin.GeometryMixin.bucket_img_into_grid_hexagon`. + :param Optional[float] fps: Frames per second. If provided, the result is normalized by the frame rate. + :param Optional[float] core_cnt: Number of CPU cores to use for parallel processing. Default is -1, which means using all available cores. + :param Optional[bool] verbose: If True, then prints progress. + :returns np.ndarray: Array of size (frames x horizontal bins x verical bins) with values representing time in seconds (if fps passed) or frames (if fps not passed) + """ + timer = SimbaTimer(start=True) - check_valid_lst( - data=data, - source=GeometryMixin.cumsum_animal_geometries_grid.__name__, - valid_dtypes=(Polygon,), - ) - check_instance( - source=GeometryMixin.cumsum_animal_geometries_grid.__name__, - instance=grid, - accepted_types=(dict,), - ) - if fps is not None: - check_int(name="fps", value=fps, min_value=1) + check_valid_lst(data=data,source=GeometryMixin.cumsum_animal_geometries_grid.__name__,valid_dtypes=(Polygon,)) + check_instance( source=GeometryMixin.cumsum_animal_geometries_grid.__name__, instance=grid, accepted_types=(dict,)) + if fps is not None: check_int(name="fps", value=fps, min_value=1) check_int(name="core_cnt", value=core_cnt, min_value=-1) - if core_cnt == -1: - core_cnt = find_core_cnt()[0] + if core_cnt == -1: core_cnt = find_core_cnt()[0] w, h = 0, 0 for k in grid.keys(): w, h = max(w, k[0]), max(h, k[1]) frm_id = np.arange(0, len(data)).reshape(-1, 1) data = np.hstack((frm_id, np.array(data).reshape(-1, 1))) img_arr = np.zeros((data.shape[0], h + 1, w + 1)) - with multiprocessing.Pool( - core_cnt, maxtasksperchild=Defaults.LARGE_MAX_TASK_PER_CHILD.value - ) as pool: - constants = functools.partial( - self._cumsum_animal_geometries_grid_helper, - grid=grid, - size=(h, w), - verbose=verbose, - ) + with multiprocessing.Pool(core_cnt, maxtasksperchild=Defaults.LARGE_MAX_TASK_PER_CHILD.value) as pool: + constants = functools.partial(self._cumsum_animal_geometries_grid_helper, grid=grid, size=(h, w), verbose=verbose) for cnt, result in enumerate(pool.imap(constants, data, chunksize=1)): img_arr[cnt] = result timer.stop_timer() - stdout_success( - msg="Cumulative animal geometries in grid complete", - elapsed_time=timer.elapsed_time_str, - ) + stdout_success(msg="Cumulative animal geometries in grid complete", elapsed_time=timer.elapsed_time_str,) if fps is None: return np.cumsum(img_arr, axis=0) else: diff --git a/simba/mixins/timeseries_features_mixin.py b/simba/mixins/timeseries_features_mixin.py index 9717bd3dd..140c41f97 100644 --- a/simba/mixins/timeseries_features_mixin.py +++ b/simba/mixins/timeseries_features_mixin.py @@ -1942,7 +1942,7 @@ def linearity_index(x: np.ndarray) -> float: :align: center .. seealso:: - :func:`simba.mixins.timeseries_features_mixin.TimeseriesFeatureMixin.sliding_linearity_index` + :func:`simba.mixins.timeseries_features_mixin.TimeseriesFeatureMixin.sliding_linearity_index`, :func:`simba.data_processors.cuda.timeseries.sliding_linearity_index_cuda` .. math:: \text{linearity\_index} = \frac{\text{straight\_line\_distance}}{\text{path\_length}} @@ -1982,7 +1982,7 @@ def sliding_linearity_index(x: np.ndarray, The Linearity Index measures how straight a path is by comparing the straight-line distance between the start and end points of each window to the total distance traveled along the path. .. seealso:: - :func:`simba.mixins.timeseries_features_mixin.TimeseriesFeatureMixin.linearity_index` + :func:`simba.mixins.timeseries_features_mixin.TimeseriesFeatureMixin.linearity_index`, :func:`simba.data_processors.cuda.timeseries.sliding_linearity_index_cuda` :param np.ndarray x: An (N, M) array representing the path, where N is the number of points and M is the number of spatial dimensions (e.g., 2 for 2D or 3 for 3D). Each row represents the coordinates of a point along the path. :param float x: The size of the sliding window in seconds. This defines the time window over which the linearity index is calculated. The window size should be specified in seconds. @@ -2217,7 +2217,7 @@ def spatial_density(x: np.ndarray, :align: center .. seealso:: - :func:`simba.mixins.timeseries_features_mixin.TimeseriesFeatureMixin.sliding_spatial_density` + :func:`simba.mixins.timeseries_features_mixin.TimeseriesFeatureMixin.sliding_spatial_density`, :func:`simba.data_processors.cuda.timeseries.sliding_spatial_density_cuda` :param np.ndarray x: A 2D array of shape (N, 2), where N is the number of points and each point has two spatial coordinates. :param float radius: The radius within which to count neighboring points around each point. Defines the area of interest around each trajectory point. @@ -2261,7 +2261,7 @@ def sliding_spatial_density(x: np.ndarray, of the trajectory. .. seealso:: - :func:`simba.mixins.timeseries_features_mixin.TimeseriesFeatureMixin.spatial_density` + :func:`simba.mixins.timeseries_features_mixin.TimeseriesFeatureMixin.spatial_density`, :func:`simba.data_processors.cuda.timeseries.sliding_spatial_density_cuda` :param np.ndarray x: A 2D array of shape (N, 2), where N is the number of points and each point has two spatial coordinates (x, y). The array represents the trajectory path of points in a 2D space (e.g., x and y positions in space). :param float radius: The radius (in millimeters) within which to count neighboring points around each trajectory point. Defines the area of interest around each point. @@ -2290,3 +2290,71 @@ def sliding_spatial_density(x: np.ndarray, results[r - 1] = total_neighbors / n_points return results + + @staticmethod + def path_aspect_ratio(x: np.ndarray, px_per_mm: float) -> float: + """ + Calculates the aspect ratio of the bounding box that encloses a given path. + + .. image:: _static/img/path_aspect_ratio.webp + :width: 400 + :align: center + + :param np.ndarray x: A 2D array of shape (N, 2) representing the path, where N is the number of points and each point has two spatial coordinates (e.g., x and y for 2D space). The path should be in the form of an array of consecutive (x, y) points. + :param float px_per_mm: Convertion factor representing the number of pixels per millimeter + :return: The aspect ratio of the bounding box enclosing the path. If the width or height of the bounding box is zero (e.g., if all points are aligned vertically or horizontally), returns -1. + :rtype: float + + :example: + >>> x = np.random.randint(0, 500, (10, 2)) + >>> TimeseriesFeatureMixin.path_aspect_ratio(x=x) + """ + + check_valid_array(data=x, source=TimeseriesFeatureMixin.path_aspect_ratio.__name__, accepted_ndims=(2,), accepted_axis_1_shape=[2, ], accepted_dtypes=Formats.NUMERIC_DTYPES.value) + check_float(name=TimeseriesFeatureMixin.path_aspect_ratio.__name__, value=px_per_mm) + xmin, ymin = np.min(x[:, 0]), np.min(x[:, 1]) + xmax, ymax = np.max(x[:, 0]), np.max(x[:, 1]) + w, h = (xmax - xmin), (ymax - ymin) + if w == 0 or h == 0: + return -1 + else: + return (w / h) * px_per_mm + + @staticmethod + def sliding_path_aspect_ratio(x: np.ndarray, + window_size: float, + sample_rate: float, + px_per_mm: float) -> np.ndarray: + + """ + Computes the aspect ratio of the bounding box for a sliding window along a path. + + This function calculates the aspect ratio (width/height) of the smallest bounding box that encloses a sequence of points within a sliding window over a 2D path. The path is defined by consecutive (x, y) coordinates. The sliding window moves forward by one point at each step, and the aspect ratio is computed for each position of the window. + + :param np.ndarray x: A 2D array of shape (N, 2) representing the path, where N is the number of points, and each point has two spatial coordinates (x and y). + :param float window_size: The size of the sliding window in seconds. + :param float px_per_mm: Convertion factor representing the number of pixels per millimeter + :param float sample_rate: The sample rate of the path data in points per second. + + :return: An array of aspect ratios for each position of the sliding window. If the window contains a path segment that is aligned vertically or horizontally (leading to a zero width or height), the function returns -1.0 for that position. NaN values are used for the initial positions where the window cannot be fully applied. + :rtype: np.ndarray + + :example: + >>> x = np.random.randint(0, 500, (10, 2)) + >>> TimeseriesFeatureMixin.(x=x, window_size=1, sample_rate=2) + """ + + window_frm = np.ceil(window_size * sample_rate).astype(np.int32) + results = np.full(x.shape[0], dtype=np.float32, fill_value=np.nan) + for r in range(window_frm, x.shape[0] + 1): + l = r - window_frm + sample_x = x[l:r, :] + xmin, ymin = np.min(sample_x[:, 0]), np.min(sample_x[:, 1]) + xmax, ymax = np.max(sample_x[:, 0]), np.max(sample_x[:, 1]) + w, h = (xmax - xmin), (ymax - ymin) + if w == 0 or h == 0: + results[r - 1] = -1.0 + else: + results[r - 1] = (w / h) * px_per_mm + + return results \ No newline at end of file diff --git a/simba/utils/errors.py b/simba/utils/errors.py index fc84e5e8e..69e865767 100644 --- a/simba/utils/errors.py +++ b/simba/utils/errors.py @@ -415,6 +415,12 @@ def __init__(self, msg: str, source: str = "", show_window: bool = False): super().__init__(msg=msg, source=source, show_window=show_window) +class SimBAGPUError(SimbaError): + def __init__(self, msg: str, source: str = "", show_window: bool = False): + msg = f"SIMBA GPU ERROR: {msg}" + super().__init__(msg=msg, source=source, show_window=show_window) + + # test = NoSpecifiedOutputError(msg='test', source='test.method') # test = FFMPEGNotFoundError(msg='323')