diff --git a/docs/_static/img/slice_imgs_gpu.webm b/docs/_static/img/slice_imgs_gpu.webm new file mode 100644 index 000000000..fa21da866 Binary files /dev/null and b/docs/_static/img/slice_imgs_gpu.webm differ diff --git a/docs/tables/slice_imgs.csv b/docs/tables/slice_imgs.csv new file mode 100644 index 000000000..6097f771e --- /dev/null +++ b/docs/tables/slice_imgs.csv @@ -0,0 +1,9 @@ +FRAMES N,TIME (S) +500,2.92 +1000,2.134 +2000,3.881 +4000,7.358 +8000,15.066 +16000,35.156 +32000,57.9812 +IMG SIZE: 714 x 528, diff --git a/setup.py b/setup.py index 293145fd3..6172ab9d1 100644 --- a/setup.py +++ b/setup.py @@ -24,7 +24,7 @@ # Setup configuration setuptools.setup( name="Simba-UW-tf-dev", - version="2.1.2", + version="2.1.3", author="Simon Nilsson, Jia Jie Choong, Sophia Hwang", author_email="sronilsson@gmail.com", description="Toolkit for computer classification and analysis of behaviors in experimental animals", diff --git a/simba/data_processors/cuda/angle_3pt.py b/simba/data_processors/cuda/angle_3pt.py deleted file mode 100644 index a00d4ee6e..000000000 --- a/simba/data_processors/cuda/angle_3pt.py +++ /dev/null @@ -1,64 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -import math - -import numpy as np -from numba import cuda - -from simba.utils.read_write import read_df - -THREADS_PER_BLOCK = 256 -@cuda.jit -def _get_3pt_angle_kernel(x_dev, y_dev, z_dev, results): - i = cuda.grid(1) - - if i >= x_dev.shape[0]: - return - if i < x_dev.shape[0]: - x_x, x_y = x_dev[i][0], x_dev[i][1] - y_x, y_y = y_dev[i][0], y_dev[i][1] - z_x, z_y = z_dev[i][0], z_dev[i][1] - D = math.degrees(math.atan2(z_y - y_y, z_x - y_x) - math.atan2(x_y - y_y, x_x - y_x)) - if D < 0: - D += 360 - results[i] = D - -def get_3pt_angle(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> np.ndarray: - """ - Computes the angle formed by three points in 2D space for each corresponding row in the input arrays using - GPU. The points x, y, and z represent the coordinates of three points in space, and the angle is calculated - at point `y` between the line segments `xy` and `yz`. - - .. image:: _static/img/get_3pt_angle_cuda.png - :width: 500 - :align: center - - :param x: A numpy array of shape (n, 2) representing the first point (e.g., nose) coordinates. - :param y: A numpy array of shape (n, 2) representing the second point (e.g., center) coordinates, where the angle is computed. - :param z: A numpy array of shape (n, 2) representing the second point (e.g., center) coordinates, where the angle is computed. - :return: A numpy array of shape (n, 1) containing the calculated angles (in degrees) for each row. - - :example: - >>> video_path = r"/mnt/c/troubleshooting/mitra/project_folder/videos/501_MA142_Gi_CNO_0514.mp4" - >>> data_path = r"/mnt/c/troubleshooting/mitra/project_folder/csv/outlier_corrected_movement_location/501_MA142_Gi_CNO_0514 - test.csv" - >>> df = read_df(file_path=data_path, file_type='csv') - >>> y = df[['Center_x', 'Center_y']].values - >>> x = df[['Nose_x', 'Nose_y']].values - >>> z = df[['Tail_base_x', 'Tail_base_y']].values - >>> angle_x = get_3pt_angle(x=x, y=y, z=z) - """ - - - x = np.ascontiguousarray(x).astype(np.float32) - y = np.ascontiguousarray(y).astype(np.float32) - n, m = x.shape - x_dev = cuda.to_device(x) - y_dev = cuda.to_device(y) - z_dev = cuda.to_device(z) - results = cuda.device_array((n, m), dtype=np.int32) - bpg = (n + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK - _get_3pt_angle_kernel[bpg, THREADS_PER_BLOCK](x_dev, y_dev, z_dev, results) - results = results.copy_to_host() - cuda.current_context().memory_manager.deallocations.clear() - return results diff --git a/simba/data_processors/cuda/circular_statistics.py b/simba/data_processors/cuda/circular_statistics.py deleted file mode 100644 index fccd5fe8d..000000000 --- a/simba/data_processors/cuda/circular_statistics.py +++ /dev/null @@ -1,680 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -import math -from typing import Optional, Tuple - -try: - from typing import Literal -except: - from typing_extensions import Literal - -import numpy as np -from numba import cuda, int32 - -try: - import cupy as cp -except: - import numpy as cp - -from simba.utils.checks import check_float, check_int, check_valid_array -from simba.utils.enums import Formats - -THREADS_PER_BLOCK = 1024 - -@cuda.jit() -def _cuda_direction_from_two_bps(x, y, results): - i = cuda.grid(1) - if i > x.shape[0]: - return - else: - a = math.atan2(x[i][0] - y[i][0], y[i][1] - x[i][1]) * (180 / math.pi) - a = int32(a + 360 if a < 0 else a) - results[i] = a - - -def direction_from_two_bps(x: np.ndarray, y: np.ndarray) -> np.ndarray: - """ - Compute the directionality in degrees from two body-parts. E.g., ``nape`` and ``nose``, - or ``swim_bladder`` and ``tail`` with GPU acceleration. - - .. csv-table:: - :header: EXPECTED RUNTIMES - :file: ../../../docs/tables/direction_two_bps.csv - :widths: 10, 90 - :align: center - :header-rows: 1 - - .. seealso:: - For CPU function see :func:`~simba.mixins.circular_statistics.CircularStatisticsMixin.direction_two_bps`. - - :parameter np.ndarray x: Size len(frames) x 2 representing x and y coordinates for first body-part. - :parameter np.ndarray y: Size len(frames) x 2 representing x and y coordinates for second body-part. - :return: Frame-wise directionality in degrees. - :rtype: np.ndarray. - - """ - x = np.ascontiguousarray(x).astype(np.int32) - y = np.ascontiguousarray(y).astype(np.int32) - x_dev = cuda.to_device(x) - y_dev = cuda.to_device(y) - results = cuda.device_array((x.shape[0]), dtype=np.int32) - bpg = (x.shape[0] + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK - _cuda_direction_from_two_bps[bpg, THREADS_PER_BLOCK](x_dev, y_dev, results) - results = results.copy_to_host() - return results - - -def sliding_circular_hotspots(x: np.ndarray, - time_window: float, - sample_rate: float, - bins: np.ndarray, - batch_size: Optional[int] = int(3.5e+7)) -> np.ndarray: - """ - Calculate the proportion of data points falling within specified circular bins over a sliding time window using GPU - - This function processes time series data representing angles (in degrees) and calculates the proportion of data - points within specified angular bins over a sliding window. The calculations are performed in batches to - accommodate large datasets efficiently. - - .. csv-table:: - :header: EXPECTED RUNTIMES - :file: ../../../docs/tables/sliding_circular_hotspots.csv - :widths: 10, 45, 45 - :align: center - :header-rows: 1 - - .. seealso:: - For CPU function see :func:`~simba.mixins.circular_statistics.CircularStatisticsMixin.sliding_circular_hotspots`. - - - :param np.ndarray x: The input time series data in degrees. Should be a 1D numpy array. - :param float time_window: The size of the sliding window in seconds. - :param float sample_rate: The sample rate of the time series data (i.e., hz, fps). - :param ndarray bins: 2D array of shape representing circular bins defining [start_degree, end_degree] inclusive. - :param Optional[int] batch_size: The size of each batch for processing the data. Default is 5e+7 (50m). - :return: A 2D numpy array where each row corresponds to a time point in `data`, and each column represents a circular bin. The values in the array represent the proportion of data points within each bin at each time point. The first column represents the first bin. - :rtype: np.ndarray - """ - - n = x.shape[0] - x = cp.asarray(x, dtype=cp.float16) - results = cp.full((x.shape[0], bins.shape[0]), dtype=cp.float16, fill_value=-1) - window_size = int(cp.ceil(time_window * sample_rate)) - for cnt, left in enumerate(range(0, n, batch_size)): - right = int(min(left + batch_size, n)) - if cnt > 0: - left = left - window_size + 1 - x_batch = x[left:right] - x_batch = cp.lib.stride_tricks.sliding_window_view(x_batch, window_size).astype(cp.float16) - batch_results = cp.full((x_batch.shape[0], bins.shape[0]), dtype=cp.float16, fill_value=-1) - for bin_cnt in range(bins.shape[0]): - if bins[bin_cnt][0] > bins[bin_cnt][1]: - mask = ((x_batch >= bins[bin_cnt][0]) & (x_batch <= 360)) | ((x_batch >= 0) & (x_batch <= bins[bin_cnt][1])) - else: - mask = (x_batch >= bins[bin_cnt][0]) & (x_batch <= bins[bin_cnt][1]) - count_per_row = cp.array(mask.sum(axis=1) / window_size).reshape(-1, ) - batch_results[:, bin_cnt] = count_per_row - results[left + window_size - 1:right, ] = batch_results - return results.get() - - - - -def sliding_circular_mean(x: np.ndarray, - time_window: float, - sample_rate: int, - batch_size: Optional[int] = 3e+7) -> np.ndarray: - - """ - Calculate the sliding circular mean over a time window for a series of angles. - - This function computes the circular mean of angles in the input array `x` over a specified sliding window. - The circular mean is a measure of the average direction for angles, which is especially useful for angular data - where traditional averaging would not be meaningful due to the circular nature of angles (e.g., 359° and 1° should average to 0°). - - The calculation is performed using a sliding window approach, where the circular mean is computed for each window - of angles. The function leverages GPU acceleration via CuPy for efficiency when processing large datasets. - - The circular mean :math:`\\mu` for a set of angles is calculated using the following formula: - - .. math:: - - \\mu = \\text{atan2}\\left(\\frac{1}{N} \\sum_{i=1}^{N} \\sin(\\theta_i), \\frac{1}{N} \\sum_{i=1}^{N} \\cos(\\theta_i)\\right) - - - :math:`\\theta_i` are the angles in radians within the sliding window - - :math:`N` is the number of samples in the window - - - .. csv-table:: - :header: EXPECTED RUNTIMES - :file: ../../../docs/tables/sliding_circular_mean.csv - :widths: 10, 45, 45 - :align: center - :header-rows: 1 - - .. seealso:: - For CPU function see :func:`~simba.mixins.circular_statistics.CircularStatisticsMixin.sliding_circular_mean`. - - :param np.ndarray x: Input array containing angle values in degrees. The array should be 1-dimensional. - :param float time_window: Time duration for the sliding window, in seconds. This determines the number of samples in each window based on the `sample_rate`. - :param int sample_rate: The number of samples per second (i.e., FPS). This is used to calculate the window size in terms of array indices. - :param Optional[int] batch_size: The maximum number of elements to process in each batch. This is used to handle large arrays by processing them in chunks to avoid memory overflow. Defaults to 3e+7 (30 million elements). - :return np.ndarray: A 1D numpy array of the same length as `x`, containing the circular mean for each sliding window. Values before the window is fully populated will be set to -1. - - :example: - >>> x = np.random.randint(0, 361, (i, )).astype(np.int32) - >>> results = sliding_circular_mean(x, 1, 10) - """ - - - window_size = np.ceil(time_window * sample_rate).astype(np.int64) - n = x.shape[0] - results = cp.full(x.shape[0], -1, dtype=np.int32) - for cnt, left in enumerate(range(0, int(n), int(batch_size))): - right = np.int32(min(left + batch_size, n)) - if cnt > 0: - left = left - window_size+1 - x_batch = cp.asarray(x[left:right]) - x_batch = cp.lib.stride_tricks.sliding_window_view(x_batch, window_size) - x_batch = np.deg2rad(x_batch) - cos, sin = cp.cos(x_batch).astype(np.float32), cp.sin(x_batch).astype(np.float32) - r = cp.rad2deg(cp.arctan2(cp.mean(sin, axis=1), cp.mean(cos, axis=1))) - r = cp.where(r < 0, r + 360, r) - results[left + window_size - 1:right] = r - return results.get() - - - -def sliding_circular_range(x: np.ndarray, - time_window: float, - sample_rate: float, - batch_size: Optional[int] = int(5e+7)) -> np.ndarray: - """ - Computes the sliding circular range of a time series data array using GPU. - - This function calculates the circular range of a time series data array using a sliding window approach. - The input data is assumed to be in degrees, and the function handles the circular nature of the data - by considering the circular distance between angles. - - .. math:: - - R = \\min \\left( \\text{max}(\\Delta \\theta) - \\text{min}(\\Delta \\theta), \\, 360 - \\text{max}(\\Delta \\theta) + \\text{min}(\\Delta \\theta) \\right) - - where: - - - :math:`\\Delta \\theta` is the difference between angles within the window, - - :math:`360` accounts for the circular nature of the data (i.e., wrap-around at 360 degrees). - - .. csv-table:: - :header: EXPECTED RUNTIMES - :file: ../../../docs/tables/sliding_circular_range.csv - :widths: 10, 45, 45 - :align: center - :header-rows: 1 - - .. seealso:: - For CPU function see :func:`~simba.mixins.circular_statistics.CircularStatisticsMixin.sliding_circular_range`. - - :param np.ndarray x: The input time series data in degrees. Should be a 1D numpy array. - :param float time_window: The size of the sliding window in seconds. - :param float sample_rate: The sample rate of the time series data (i.e., hz, fps). - :param Optional[int] batch_size: The size of each batch for processing the data. Default is 5e+7 (50m). - :return: A numpy array containing the sliding circular range values. - :rtype: np.ndarray - - :example: - >>> x = np.random.randint(0, 361, (19, )).astype(np.int32) - >>> p = sliding_circular_range(x, 1, 10) - """ - - n = x.shape[0] - x = cp.asarray(x, dtype=cp.float16) - results = cp.zeros_like(x, dtype=cp.int16) - x = cp.deg2rad(x).astype(cp.float16) - window_size = int(cp.ceil(time_window * sample_rate)) - for cnt, left in enumerate(range(0, n, batch_size)): - right = int(min(left + batch_size, n)) - if cnt > 0: - left = left - window_size + 1 - x_batch = x[left:right] - x_batch = cp.lib.stride_tricks.sliding_window_view(x_batch, window_size).astype(cp.float16) - x_batch = cp.sort(x_batch) - results[left + window_size - 1:right] = cp.abs(cp.rint(cp.rad2deg(cp.amin(cp.vstack([x_batch[:, -1] - x_batch[:, 0], 2 * cp.pi - cp.max(cp.diff(x_batch), axis=1)]).T, axis=1)))) - return results.get() - - - - -def sliding_circular_std(x: np.ndarray, - time_window: float, - sample_rate: float, - batch_size: Optional[int] = int(5e+7)) -> np.ndarray: - - """ - Calculate the sliding circular standard deviation of a time series data on GPU. - - This function computes the circular standard deviation over a sliding window for a given time series array. - The time series data is assumed to be in degrees, and the function converts it to radians for computation. - The sliding window approach is used to handle large datasets efficiently, processing the data in batches. - - The circular standard deviation (σ) is computed using the formula: - - .. math:: - - \sigma = \sqrt{-2 \cdot \log \left|\text{mean}\left(\exp(i \cdot x_{\text{batch}})\right)\right|} - - where :math:`x_{\text{batch}}` is the data within the current sliding window, and :math:`\text{mean}` and - :math:`\log` are computed in the circular (complex plane) domain. - - .. csv-table:: - :header: EXPECTED RUNTIMES - :file: ../../../docs/tables/sliding_circular_std.csv - :widths: 10, 45, 45 - :align: center - :header-rows: 1 - - .. seealso:: - For CPU function see :func:`~simba.mixins.circular_statistics.CircularStatisticsMixin.sliding_circular_std`. - - :param np.ndarray x: The input time series data in degrees. Should be a 1D numpy array. - :param float time_window: The size of the sliding window in seconds. - :param float sample_rate: The sample rate of the time series data (i.e., hz, fps). - :param Optional[int] batch_size: The size of each batch for processing the data. Default is 5e+7 (50m). - - :return: A numpy array containing the sliding circular standard deviation values. - :rtype: np.ndarray - """ - - - n = x.shape[0] - x = cp.asarray(x, dtype=cp.float16) - results = cp.zeros_like(x, dtype=cp.float16) - x = np.deg2rad(x).astype(cp.float16) - window_size = int(np.ceil(time_window * sample_rate)) - for cnt, left in enumerate(range(0, n, batch_size)): - right = int(min(left + batch_size, n)) - if cnt > 0: - left = left - window_size + 1 - x_batch = x[left:right] - x_batch = cp.lib.stride_tricks.sliding_window_view(x_batch, window_size).astype(cp.float16) - m = cp.log(cp.abs(cp.mean(cp.exp(1j * x_batch), axis=1))) - stdev = cp.rad2deg(cp.sqrt(-2 * m)) - results[left + window_size - 1:right] = stdev - - return results.get() - - -def sliding_rayleigh_z(x: np.ndarray, - time_window: float, - sample_rate: float, - batch_size: Optional[int] = int(5e+7)) -> Tuple[np.ndarray, np.ndarray]: - - """ - Computes the Rayleigh Z-statistic over a sliding window for a given time series of angles - - This function calculates the Rayleigh Z-statistic, which tests the null hypothesis that the population of angles - is uniformly distributed around the circle. The calculation is performed over a sliding window across the input - time series, and results are computed in batches for memory efficiency. - - Data is processed using GPU acceleration via CuPy, which allows for faster computation compared to a CPU-based approach. - - .. note:: - Adapted from ``pingouin.circular.circ_rayleigh`` and ``pycircstat.tests.rayleigh``. - - - **Rayleigh Z-statistic:** - - The Rayleigh Z-statistic is given by: - - .. math:: - - R = \frac{1}{n} \sqrt{\left(\sum_{i=1}^{n} \cos(\theta_i)\right)^2 + \left(\sum_{i=1}^{n} \sin(\theta_i)\right)^2} - - where: - - :math:`\theta_i` are the angles in the window. - - :math:`n` is the number of angles in the window. - - - .. csv-table:: - :header: EXPECTED RUNTIMES - :file: ../../../docs/tables/sliding_rayleigh_z.csv - :widths: 10, 45, 45 - :align: center - :header-rows: 1 - - .. seealso:: - For CPU function see :func:`~simba.mixins.circular_statistics.CircularStatisticsMixin.sliding_rayleigh_z`. - - - :param np.ndarray x: Input array of angles in degrees. Should be a 1D numpy array. - :param float time_window: The size of the sliding window in time units (e.g., seconds). - :param float sample_rate: The sampling rate of the input time series in samples per time unit (e.g., Hz, fps). - :param Optional[int] batch_size: The number of samples to process in each batch. Default is 5e7 (50m). Reducing this value may save memory at the cost of longer computation time. - :return: - A tuple containing two numpy arrays: - - **z_results**: Rayleigh Z-statistics for each position in the input array where the window was fully applied. - - **p_results**: Corresponding p-values for the Rayleigh Z-statistics. - :rtype: Tuple[np.ndarray, np.ndarray] - """ - - n = x.shape[0] - x = cp.asarray(x, dtype=cp.float16) - z_results = cp.zeros_like(x, dtype=cp.float16) - p_results = cp.zeros_like(x, dtype=cp.float16) - x = np.deg2rad(x).astype(cp.float16) - window_size = int(np.ceil(time_window * sample_rate)) - for cnt, left in enumerate(range(0, n, batch_size)): - right = int(min(left + batch_size, n)) - if cnt > 0: - left = left - window_size + 1 - x_batch = x[left:right] - x_batch = cp.lib.stride_tricks.sliding_window_view(x_batch, window_size).astype(cp.float16) - cos_sums = cp.nansum(cp.cos(x_batch), axis=1) ** 2 - sin_sums = cp.nansum(cp.sin(x_batch), axis=1) ** 2 - R = cp.sqrt(cos_sums + sin_sums) / window_size - Z = window_size * (R**2) - P = cp.exp(np.sqrt(1 + 4 * window_size + 4 * (window_size ** 2 - R ** 2)) - (1 + 2 * window_size)) - z_results[left + window_size - 1:right] = Z - p_results[left + window_size - 1:right] = P - - return z_results.get(), p_results.get() - - -def sliding_resultant_vector_length(x: np.ndarray, - time_window: float, - sample_rate: int, - batch_size: Optional[int] = 3e+7) -> np.ndarray: - - """ - Calculate the sliding resultant vector length over a time window for a series of angles. - - This function computes the resultant vector length (R) for each window of angles in the input array `x`. - The resultant vector length is a measure of the concentration of angles, and it ranges from 0 to 1, where 1 - indicates all angles point in the same direction, and 0 indicates uniform distribution of angles. - - For a given sliding window of angles, the resultant vector length :math:`R` is calculated using the following formula: - - .. math:: - - R = \\frac{1}{N} \\sqrt{\\left(\\sum_{i=1}^{N} \\cos(\\theta_i)\\right)^2 + \\left(\\sum_{i=1}^{N} \\sin(\\theta_i)\\right)^2} - - where: - - - :math:`\\theta_i` are the angles in radians within the sliding window - - :math:`N` is the number of samples in the window - - The computation is performed in a sliding window manner over the entire array, utilizing GPU acceleration - with CuPy for efficiency, especially on large datasets. - - - .. csv-table:: - :header: EXPECTED RUNTIMES - :file: ../../../docs/tables/sliding_resultant_vector_length.csv - :widths: 10, 10, 80 - :align: center - :header-rows: 1 - - .. seealso:: - For CPU function see :func:`~simba.mixins.circular_statistics.CircularStatisticsMixin.sliding_resultant_vector_length`. - - :param np.ndarray x: Input array containing angle values in degrees. The array should be 1-dimensional. - :param float time_window: Time duration for the sliding window, in seconds. This determines the number of samples in each window based on the `sample_rate`. - :param int sample_rate: The number of samples per second (i.e., FPS). This is used to calculate the window size in terms of array indices. - :param Optional[int] batch_size: The maximum number of elements to process in each batch. This is used to handle large arrays by processing them in chunks to avoid memory overflow. Defaults to 3e+7 (30 million elements). - :return np.ndarray: A 1D numpy array of the same length as `x`, containing the resultant vector length for each sliding window. Values before the window is fully populated will be set to -1. - - :example: - >>> x = np.random.randint(0, 361, (5000, )).astype(np.int32) - >>> results = sliding_resultant_vector_length(x, 1, 10) - """ - - window_size = np.ceil(time_window * sample_rate).astype(np.int64) - n = x.shape[0] - results = cp.full(x.shape[0], -1, dtype=np.float32) - for cnt, left in enumerate(range(0, int(n), int(batch_size))): - right = np.int32(min(left + batch_size, n)) - if cnt > 0: - left = left - window_size+1 - x_batch = cp.asarray(x[left:right]) - x_batch = cp.lib.stride_tricks.sliding_window_view(x_batch, window_size) - x_batch = np.deg2rad(x_batch) - cos, sin = cp.cos(x_batch).astype(np.float32), cp.sin(x_batch).astype(np.float32) - cos_sum, sin_sum = cp.sum(cos, axis=1), cp.sum(sin, axis=1) - r = np.sqrt(cos_sum ** 2 + sin_sum ** 2) / window_size - results[left+window_size-1:right] = r - return results.get() - - -def direction_from_three_bps(x: np.ndarray, - y: np.ndarray, - z: np.ndarray, - batch_size: Optional[int] = int(1.5e+7)) -> np.ndarray: - - """ - Calculate the direction angle based on the coordinates of three body points using GPU acceleration. - - This function computes the mean direction angle (in degrees) for a batch of coordinates - provided in the form of NumPy arrays. The calculation is based on the arctangent of the - difference in x and y coordinates between pairs of points. The result is a value in - the range [0, 360) degrees. - - .. seealso:: - :func:`simba.mixins.circular_statistics.CircularStatisticsMixin.direction_three_bps` - - :param np.ndarray x: A 2D array of shape (N, 2) containing the x-coordinates of the first body part (nose) - :param np.ndarray y: A 2D array of shape (N, 2) containing the coordinates of the second body part (left ear). - :param np.ndarray z: A 2D array of shape (N, 2) containing the coordinates of the second body part (right ear). - :param Optional[int] batch_size: The size of the batch to be processed in each iteration. Default is 15 million. - :return: An array of shape (N,) containing the computed direction angles in degrees. - :rtype np.ndarray: - """ - - check_valid_array(data=x, source=direction_from_three_bps.__name__, accepted_ndims=(2,), accepted_dtypes=Formats.NUMERIC_DTYPES.value) - check_valid_array(data=y, source=direction_from_three_bps.__name__, accepted_shapes=(x.shape,), accepted_dtypes=Formats.NUMERIC_DTYPES.value) - check_valid_array(data=z, source=direction_from_three_bps.__name__, accepted_shapes=(x.shape,), accepted_dtypes=Formats.NUMERIC_DTYPES.value) - check_int(value=batch_size, name=direction_from_three_bps.__name__, min_value=1) - results = cp.full((x.shape[0]), fill_value=-1, dtype=np.int16) - - for l in range(0, x.shape[0], batch_size): - r = l + batch_size - x_batch = cp.array(x[l:r]) - y_batch = cp.array(y[l:r]) - z_batch = cp.array(z[l:r]) - left_ear_to_nose = cp.arctan2(x_batch[:, 0] - y_batch[:, 0], y_batch[:, 1] - x_batch[:,1]) - right_ear_nose = cp.arctan2(x_batch[:, 0] - z_batch[:, 0], z_batch[:, 1] - x_batch[:, 1]) - mean_angle_rad = cp.arctan2(cp.sin(left_ear_to_nose) + cp.sin(right_ear_nose), cp.cos(left_ear_to_nose) + cp.cos(right_ear_nose)) - results[l:r] = (cp.degrees(mean_angle_rad) + 360) % 360 - - return results.get() - - -@cuda.jit() -def _instantaneous_angular_velocity(x, stride, results): - r = cuda.grid(1) - l = np.int32(r - (stride[0])) - if (r > results.shape[0]) or (l < 0): - results[r] = -1 - else: - d = math.pi - (abs(math.pi - abs(x[l] - x[r]))) - results[r] = d * (180 / math.pi) - - -def instantaneous_angular_velocity(x: np.ndarray, stride: Optional[int] = 1) -> np.ndarray: - """ - Calculate the instantaneous angular velocity between angles in a given array. - - This function uses CUDA to perform parallel computations on the GPU. - - The angular velocity is computed using the difference in angles between - the current and previous values (with a specified stride) in the array. - The result is returned in degrees per unit time. - - .. csv-table:: - :header: EXPECTED RUNTIMES - :file: ../../../docs/tables/instantaneous_angular_velocity.csv - :widths: 10, 90 - :align: center - :header-rows: 1 - - .. math:: - \omega = \frac{{\Delta \theta}}{{\Delta t}} = \frac{{180}}{{\pi}} \times \left( \pi - \left| \pi - \left| \theta_r - \theta_l \right| \right| \right) - - where: - - \( \theta_r \) is the current angle. - - \( \theta_l \) is the angle at the specified stride before the current angle. - - \( \Delta t \) is the time difference between the two angles. - - - .. seealso:: - :func:`simba.mixins.circular_statistics.CircularStatisticsMixin.instantaneous_angular_velocity` - - :param np.ndarray x: Array of angles in degrees, for which the instantaneous angular velocity will be calculated. - :param Optional[int] stride: The stride or lag (in frames) to use when calculating the difference in angles. Defaults to 1. - :return: Array of instantaneous angular velocities corresponding to the input angles. Velocities are in degrees per unit time. - :rtype: np.ndarray - """ - - x = np.deg2rad(x).astype(np.int16) - stride = np.array([stride]).astype(np.int64) - bpg = (x.shape[0] + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK - x_dev = cuda.to_device(x) - stride_dev = cuda.to_device(stride) - results = cuda.device_array(x.shape[0], dtype=np.float32) - _instantaneous_angular_velocity[bpg, THREADS_PER_BLOCK](x_dev, stride_dev, results) - return results.copy_to_host() - - -@cuda.jit(device=True) -def _rad2deg(x): - return x * (180/math.pi) - - -@cuda.jit() -def _sliding_bearing(x, stride, results): - r = cuda.grid(1) - l = np.int32(r - (stride[0])) - if (r > results.shape[0]-1) or (l < 0): - results[r] = -1 - else: - x1, y1 = x[l, 0], x[l, 1] - x2, y2 = x[r, 0], x[r, 1] - bearing = _rad2deg(math.atan2(x2 - x1, y2 - y1)) - results[r] = (bearing + 360) % 360 - - -def sliding_bearing(x: np.ndarray, - stride: Optional[float] = 1, - sample_rate: Optional[float] = 1) -> np.ndarray: - """ - Compute the bearing between consecutive points in a 2D coordinate array using a sliding window approach using GPU acceleration. - - This function calculates the angle (bearing) in degrees between each point and a point a certain number of - steps ahead (defined by `stride`) in the 2D coordinate array `x`. The bearing is calculated using the - arctangent of the difference in coordinates, converted from radians to degrees. - - .. csv-table:: - :header: EXPECTED RUNTIMES - :file: ../../../docs/tables/sliding_bearing.csv - :widths: 10, 90 - :align: center - :header-rows: 1 - - .. seealso:: - :func:`simba.mixins.circular_statistics.CircularStatisticsMixin.sliding_bearing` - - :param np.ndarray x: A 2D array of shape `(n, 2)` where each row represents a point with `x` and `y` coordinates. The array must be numeric. - :param Optional[float] stride: The time (multiplied by `sample_rate`) to look ahead when computing the bearing in seconds. Defaults to 1. - :param Optional[float] sample_rate: A multiplier applied to the `stride` value to determine the actual step size for calculating the bearing. E.g., frames per second. Defaults to 1. If the resulting stride is less than 1, it is automatically set to 1. - :return:A 1D array of shape `(n,)` containing the calculated bearings in degrees. Values outside the valid range (i.e., where the stride exceeds array bounds) are set to -1. - :rtype: np.ndarray - """ - - check_valid_array(data=x, source=f'{sliding_bearing.__name__} x', accepted_ndims=(2,), accepted_axis_1_shape=(2,), accepted_dtypes=Formats.NUMERIC_DTYPES.value) - check_float(name=f'{sliding_bearing.__name__} stride', value=stride, min_value=10e-6, max_value=x.shape[0]-1) - check_float(name=f'{sliding_bearing.__name__} sample_rate', value=sample_rate, min_value=10e-6, max_value=x.shape[0]-1) - stride = int(stride * sample_rate) - if stride < 1: - stride = 1 - stride = np.array([stride]).astype(np.int64) - bpg = (x.shape[0] + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK - x_dev = cuda.to_device(x) - stride_dev = cuda.to_device(stride) - results = cuda.device_array(x.shape[0], dtype=np.float32) - _sliding_bearing[bpg, THREADS_PER_BLOCK](x_dev, stride_dev, results) - return results.copy_to_host() - - -@cuda.jit(device=True) -def _rad2deg(x): - return x * (180 / math.pi) - - -@cuda.jit() -def _sliding_angular_diff(data, strides, results): - x, y = cuda.grid(2) - if (x > data.shape[0] - 1) or (y > strides.shape[0] - 1): - return - else: - stride = int(strides[y]) - if x - stride < 0: - return - a_2 = data[x] - a_1 = data[x - stride] - distance = math.pi - abs(math.pi - abs(a_1 - a_2)) - distance = abs(int(_rad2deg(distance)) + 1) - results[x][y] = distance - - -def sliding_angular_diff(x: np.ndarray, - time_windows: np.ndarray, - fps: float) -> np.ndarray: - """ - Calculate the sliding angular differences for a given time window using GPU acceleration. - - - This function computes the angular differences between each angle in `x` - and the corresponding angle located at a distance determined by the time window - and frame rate (fps). The results are returned as a 2D array where each row corresponds - to a position in `x`, and each column corresponds to a different time window. - - .. csv-table:: - :header: EXPECTED RUNTIMES - :file: ../../../docs/tables/sliding_angular_diff.csv - :widths: 10, 90 - :align: center - :header-rows: 1 - - - .. seealso:: - :func:`simba.mixins.circular_statistics.CircularStatisticsMixin.sliding_angular_diff` - - .. math:: - \text{difference} = \pi - |\pi - |a_1 - a_2|| - - Where: - - \( a_1 \) is the angle at position `x`. - - \( a_2 \) is the angle at position `x - \text{stride}`. - - :param np.ndarray x: 1D array of angles in degrees. - :param np.ndarray time_windows: 1D array of time windows in seconds to determine the stride (distance in frames) between angles. - :param float fps: Frame rate (frames per second) used to convert time windows to strides. - :return: 2D array of angular differences. Each row corresponds to an angle in `x`, and each column corresponds to a time window. - :rtype: np.ndarray - """ - - x = np.deg2rad(x) - strides = np.zeros(time_windows.shape[0]) - for i in range(time_windows.shape[0]): - strides[i] = np.ceil(time_windows[i] * fps).astype(np.int32) - x_dev = cuda.to_device(x) - stride_dev = cuda.to_device(strides) - results = cuda.device_array((x.shape[0], time_windows.shape[0])) - grid_x = (x.shape[0] + THREADS_PER_BLOCK - 1) // THREADS_PER_BLOCK - grid_y = (strides.shape[0] + THREADS_PER_BLOCK - 1) - blocks_per_grid = (grid_x, grid_y) - _sliding_angular_diff[blocks_per_grid, THREADS_PER_BLOCK](x_dev, stride_dev, results) - results = results.copy_to_host().astype(np.int32) - return results - - diff --git a/simba/data_processors/cuda/convex_hull.py b/simba/data_processors/cuda/convex_hull.py deleted file mode 100644 index ba56f8f89..000000000 --- a/simba/data_processors/cuda/convex_hull.py +++ /dev/null @@ -1,111 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -import numpy as np -from numba import cuda, njit - -THREADS_PER_BLOCK = 128 - -@cuda.jit(device=True) -def _cross_test(x, y, x1, y1, x2, y2): - """Cross product test for determining whether left of line.""" - cross = (x - x1) * (y2 - y1) - (y - y1) * (x2 - x1) - return cross < 0 - -@cuda.jit -def _convex_hull_kernel(pts: np.ndarray, results: np.ndarray) -> np.ndarray: - """ - CUDA kernel for the Jarvis March algorithm. - - .. note:: - `Modified from Jacob Hultman `_ - - :param pts: M x N x 2 array where M is the number of frames, N is the number of body-parts, and 2 representing the x and y coordinates of the body-parts. - :param results: M x N array where M is the number of frames, and N is the indexes of the body-parts belonging to the hull. If -1, the body-part does not belong to the hull. - """ - row = cuda.grid(1) - if row >= pts.shape[0]: - return - - point_on_hull = 0 - min_x = pts[row, 0, 0] - for j in range(pts.shape[1]): - x = pts[row, j, 0] - if x < min_x: - min_x = x - point_on_hull = j - startpoint = point_on_hull - count = 0 - while True: - results[row, count] = point_on_hull - count += 1 - endpoint = 0 - for j in range(pts.shape[1]): - if endpoint == point_on_hull: - endpoint = j - elif _cross_test( - pts[row, j, 0], - pts[row, j, 1], - pts[row, point_on_hull, 0], - pts[row, point_on_hull, 1], - pts[row, endpoint, 0], - pts[row, endpoint, 1], - ): - endpoint = j - point_on_hull = endpoint - if endpoint == startpoint: - break - for j in range(count, pts.shape[1], 1): - results[row, j] = -1 - - -#@jit(nopython=True) -@njit("(int32[:, :, :], int32[:, :])") -def _slice_hull_idx(points: np.ndarray, point_idx: np.ndarray): - results = np.zeros_like(points) - for i in range(point_idx.shape[0]): - results[i] = points[i][point_idx[i]] - return results - - -def get_convex_hull(pts: np.ndarray) -> np.ndarray: - """ - Compute the convex hull for each set of 2D points in parallel using CUDA and the Jarvis March algorithm. - - This function processes a batch of 2D point sets (frames) and computes the convex hull for each set. The convex hull of a set of points is the smallest convex polygon that contains all the points. - - The function uses a variant of the Gift Wrapping algorithm (Jarvis March) to compute the convex hull. It finds the leftmost point, then iteratively determines the next point on the hull by checking the orientation of the remaining points. The results are stored in the `results` array, where each row corresponds to a frame and contains the indices of the points forming the convex hull. Points not on the hull are marked with `-1`. - - - .. image:: _static/img/get_convex_hull_cuda.png - :width: 300 - :align: center - - .. note:: - `Modified from Jacob Hultman `_ - - :param pts: A 3D numpy array of shape (M, N, 2) where: - - M is the number of frames. - - N is the number of points (body-parts) in each frame. - - The last dimension (2) represents the x and y coordinates of each point. - - :return: An upated 3D numpy array of shape (M, N, 2) consisting of the points in the hull. - - - :example: - >>> video_path = r"/mnt/c/troubleshooting/mitra/project_folder/videos/501_MA142_Gi_CNO_0514.mp4" - >>> data_path = r"/mnt/c/troubleshooting/mitra/project_folder/csv/outlier_corrected_movement_location/501_MA142_Gi_CNO_0514 - test.csv" - >>> df = read_df(file_path=data_path, file_type='csv') - >>> frame_data = df.values.reshape(len(df), -1, 2) - >>> x = get_convex_hull(frame_data) - """ - - pts = np.ascontiguousarray(pts).astype(np.int32) - n, m, _ = pts.shape - bpg = (n + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK - pts_dev = cuda.to_device(pts) - results = cuda.device_array((n, m), dtype=np.int32) - _convex_hull_kernel[bpg, THREADS_PER_BLOCK](pts_dev, results) - hull = results.copy_to_host().astype(np.int32) - hull = _slice_hull_idx(pts, hull) - return hull diff --git a/simba/data_processors/cuda/convex_hull_area.py b/simba/data_processors/cuda/convex_hull_area.py deleted file mode 100644 index 44fd775f2..000000000 --- a/simba/data_processors/cuda/convex_hull_area.py +++ /dev/null @@ -1,42 +0,0 @@ -from typing import Optional - -import cupy as cp -import numpy as np - -from simba.utils.checks import check_float, check_valid_array -from simba.utils.enums import Formats - - -def poly_area(data: np.ndarray, - pixels_per_mm: Optional[float] = 1.0, - batch_size: Optional[int] = int(0.5e+7)) -> np.ndarray: - - """ - Compute the area of a polygon using GPU acceleration. - - This function calculates the area of polygons defined by sets of points in a 3D array. - Each 2D slice along the first dimension represents a polygon, with each row corresponding - to a point in the polygon and each column representing the x and y coordinates. - - The computation is done in batches to handle large datasets efficiently. - - :param data: A 3D numpy array of shape (N, M, 2), where N is the number of polygons, M is the number of points per polygon, and 2 represents the x and y coordinates. - :param pixels_per_mm: Optional scaling factor to convert the area from pixels squared to square millimeters. Default is 1.0. - :param batch_size: Optional batch size for processing the data in chunks to fit in memory. Default is 0.5e+7. - :return: A 1D numpy array of shape (N,) containing the computed area of each polygon in square millimeters. - """ - - check_valid_array(data=data, source=f'{poly_area} data', accepted_ndims=(3,), accepted_dtypes=Formats.NUMERIC_DTYPES.value) - check_float(name=f'{poly_area} pixels_per_mm', min_value=10e-16, value=pixels_per_mm) - results = cp.full((data.shape[0]), fill_value=cp.nan, dtype=cp.int32) - for l in range(0, data.shape[0], batch_size): - r = l + batch_size - x = cp.asarray(data[l:r, :, 0]) - y = cp.asarray(data[l:r, :, 1]) - x_r = cp.roll(x, shift=1, axis=1) - y_r = cp.roll(y, shift=1, axis=1) - dot_xy_roll_y = cp.sum(x * y_r, axis=1) - dot_y_roll_x = cp.sum(y * x_r, axis=1) - results[l:r] = (0.5 * cp.abs(dot_xy_roll_y - dot_y_roll_x)) / pixels_per_mm - - return results.get() diff --git a/simba/data_processors/cuda/count_values_in_range.py b/simba/data_processors/cuda/count_values_in_range.py deleted file mode 100644 index c17c7e2f5..000000000 --- a/simba/data_processors/cuda/count_values_in_range.py +++ /dev/null @@ -1,53 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - - -import numpy as np -from numba import cuda - -THREADS_PER_BLOCK = 256 - -@cuda.jit -def _get_values_in_range_kernel(values, ranges, results): - i = cuda.grid(1) - if i >= values.shape[0]: - return - v = values[i] - for j in range(ranges.shape[0]): - l, u = ranges[j][0], ranges[j][1] - cnt = 0 - for k in v: - if k <= u and k >= l: - cnt += 1 - results[i, j] = cnt - - -def count_values_in_ranges(x: np.ndarray, r: np.ndarray) -> np.ndarray: - """ - Counts the number of values in each feature within specified ranges for each row in a 2D array using CUDA. - - .. image:: _static/img/get_euclidean_distance_cuda.png - :width: 500 - :align: center - - :param np.ndarray x: 2d array with feature values. - :param np.ndarray r: 2d array with lower and upper boundaries. - :return np.ndarray: 2d array of size len(x) x len(r) with the counts of values in each feature range (inclusive). - - :example: - >>> x = np.random.randint(1, 11, (10, 10)).astype(np.int8) - >>> r = np.array([[1, 6], [6, 11]]) - >>> r_x = count_values_in_ranges(x=x, r=r) - """ - - x = np.ascontiguousarray(x).astype(np.float32) - r = np.ascontiguousarray(r).astype(np.float32) - n, m = x.shape[0], r.shape[0] - values_dev = cuda.to_device(x) - ranges_dev = cuda.to_device(r) - results = cuda.device_array((n, m), dtype=np.int32) - bpg = (n + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK - _get_values_in_range_kernel[bpg, THREADS_PER_BLOCK](values_dev, ranges_dev, results) - results = results.copy_to_host() - cuda.current_context().memory_manager.deallocations.clear() - return results \ No newline at end of file diff --git a/simba/data_processors/cuda/create_average_cupy.py b/simba/data_processors/cuda/create_average_cupy.py deleted file mode 100644 index 409473653..000000000 --- a/simba/data_processors/cuda/create_average_cupy.py +++ /dev/null @@ -1,108 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - - -import os -from typing import Optional, Union - -import cupy as cp -import cv2 -import numpy as np - -from simba.utils.checks import (check_file_exist_and_readable, - check_if_dir_exists, - check_if_string_value_is_valid_video_timestamp, - check_int, check_nvidea_gpu_available, - check_that_hhmmss_start_is_before_end) -from simba.utils.data import find_frame_numbers_from_time_stamp -from simba.utils.errors import FFMPEGCodecGPUError, InvalidInputError -from simba.utils.printing import stdout_success -from simba.utils.read_write import ( - check_if_hhmmss_timestamp_is_valid_part_of_video, get_fn_ext, - get_video_meta_data, read_img_batch_from_video_gpu) - - -def create_average_frm(video_path: Union[str, os.PathLike], - start_frm: Optional[int] = None, - end_frm: Optional[int] = None, - start_time: Optional[str] = None, - end_time: Optional[str] = None, - save_path: Optional[Union[str, os.PathLike]] = None, - batch_size: Optional[int] = 3000, - verbose: Optional[bool] = False) -> Union[None, np.ndarray]: - - """ - Computes the average frame using GPU acceleration from a specified range of frames or time interval in a video file. - This average frame typically used for background substraction. - - The function reads frames from the video, calculates their average, and optionally saves the result - to a specified file. If `save_path` is provided, the average frame is saved as an image file; - otherwise, the average frame is returned as a NumPy array. - - - :param Union[str, os.PathLike] video_path: The path to the video file from which to extract frames. - :param Optional[int] start_frm: The starting frame number (inclusive). Either `start_frm`/`end_frm` or `start_time`/`end_time` must be provided, but not both. - :param Optional[int] end_frm: The ending frame number (exclusive). - :param Optional[str] start_time: The start time in the format 'HH:MM:SS' from which to begin extracting frames. - :param Optional[str] end_time: The end time in the format 'HH:MM:SS' up to which frames should be extracted. - :param Optional[Union[str, os.PathLike]] save_path: The path where the average frame image will be saved. If `None`, the average frame is returned as a NumPy array. - :param Optional[int] batch_size: The number of frames to process in each batch. Default is 3000. Increase if your RAM allows it. - :param Optional[bool] verbose: If `True`, prints progress and informational messages during execution. - :return: Returns `None` if the result is saved to `save_path`. Otherwise, returns the average frame as a NumPy array. - - :example: - >>> create_average_frm(video_path=r"C:\troubleshooting\RAT_NOR\project_folder\videos\2022-06-20_NOB_DOT_4_downsampled.mp4", verbose=True, start_frm=0, end_frm=9000) - """ - - def average_3d_stack(image_stack: np.ndarray) -> np.ndarray: - num_frames, height, width, _ = image_stack.shape - image_stack = cp.array(image_stack).astype(cp.float32) - img = cp.clip(cp.sum(image_stack, axis=0) / num_frames, 0, 255).astype(cp.uint8) - return img.get() - - if not check_nvidea_gpu_available(): - raise FFMPEGCodecGPUError(msg="No GPU found (as evaluated by nvidea-smi returning None)", source=create_average_frm.__name__) - - - if ((start_frm is not None) or (end_frm is not None)) and ((start_time is not None) or (end_time is not None)): - raise InvalidInputError(msg=f'Pass start_frm and end_frm OR start_time and end_time', source=create_average_frm.__name__) - elif type(start_frm) != type(end_frm): - raise InvalidInputError(msg=f'Pass start frame and end frame', source=create_average_frm.__name__) - elif type(start_time) != type(end_time): - raise InvalidInputError(msg=f'Pass start time and end time', source=create_average_frm.__name__) - if save_path is not None: - check_if_dir_exists(in_dir=os.path.dirname(save_path), source=create_average_frm.__name_) - check_file_exist_and_readable(file_path=video_path) - video_meta_data = get_video_meta_data(video_path=video_path) - video_name = get_fn_ext(filepath=video_path)[1] - if verbose: - print(f'Getting average frame from {video_name}...') - if (start_frm is not None) and (end_frm is not None): - check_int(name='start_frm', value=start_frm, min_value=0, max_value=video_meta_data['frame_count']) - check_int(name='end_frm', value=end_frm, min_value=0, max_value=video_meta_data['frame_count']) - if start_frm > end_frm: - raise InvalidInputError(msg=f'Start frame ({start_frm}) has to be before end frame ({end_frm}).', source=create_average_frm.__name__) - frame_ids = list(range(start_frm, end_frm)) - elif (start_time is not None) and (end_time is not None): - check_if_string_value_is_valid_video_timestamp(value=start_time, name=create_average_frm.__name__) - check_if_string_value_is_valid_video_timestamp(value=end_time, name=create_average_frm.__name__) - check_that_hhmmss_start_is_before_end(start_time=start_time, end_time=end_time, name=create_average_frm.__name__) - check_if_hhmmss_timestamp_is_valid_part_of_video(timestamp=start_time, video_path=video_path) - frame_ids = find_frame_numbers_from_time_stamp(start_time=start_time, end_time=end_time, fps=video_meta_data['fps']) - else: - frame_ids = list(range(0, video_meta_data['frame_count'])) - frame_ids = [frame_ids[i:i+batch_size] for i in range(0,len(frame_ids),batch_size)] - avg_imgs = [] - for batch_cnt in range(len(frame_ids)): - start_idx, end_idx = frame_ids[batch_cnt][0], frame_ids[batch_cnt][-1] - if start_idx == end_idx: - continue - imgs = read_img_batch_from_video_gpu(video_path=video_path, start_frm=start_idx, end_frm=end_idx, verbose=verbose) - avg_imgs.append(average_3d_stack(image_stack=imgs)) - avg_img = average_3d_stack(image_stack=np.stack(avg_imgs, axis=0)) - if save_path is not None: - cv2.imwrite(save_path, avg_img) - if verbose: - stdout_success(msg=f'Saved average frame at {save_path}', source=create_average_frm.__name__) - else: - return avg_img \ No newline at end of file diff --git a/simba/data_processors/cuda/create_average_frame_cuda.py b/simba/data_processors/cuda/create_average_frame_cuda.py deleted file mode 100644 index 6820841bf..000000000 --- a/simba/data_processors/cuda/create_average_frame_cuda.py +++ /dev/null @@ -1,152 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -import os -from typing import Optional, Union - -try: - from typing import Literal -except: - from typing_extensions import Literal - -from copy import deepcopy - -import cupy as cp -import cv2 -import numpy as np -from numba import cuda - -from simba.utils.checks import (check_file_exist_and_readable, - check_if_dir_exists, - check_if_string_value_is_valid_video_timestamp, - check_if_valid_img, check_instance, check_int, - check_nvidea_gpu_available, - check_that_hhmmss_start_is_before_end) -from simba.utils.data import find_frame_numbers_from_time_stamp -from simba.utils.errors import FFMPEGCodecGPUError, InvalidInputError -from simba.utils.printing import stdout_success -from simba.utils.read_write import ( - check_if_hhmmss_timestamp_is_valid_part_of_video, get_fn_ext, - get_video_meta_data, read_img_batch_from_video_gpu) - - -def average_3d_stack_cupy(image_stack: np.ndarray) -> np.ndarray: - num_frames, height, width, _ = image_stack.shape - image_stack = cp.array(image_stack).astype(cp.float32) - img = cp.clip(cp.sum(image_stack, axis=0) / num_frames, 0, 255).astype(cp.uint8) - return img.get() - -@cuda.jit() -def _average_3d_stack_cuda(data, results): - y, x, i = cuda.grid(3) - if i < 0 or x < 0 or y < 0: - return - if i > data.shape[0] - 1 or x > data.shape[1] - 1 or y > data.shape[2] - 1: - return - else: - sum_value = 0.0 - for n in range(data.shape[0]): - sum_value += data[n, y, x, i] - results[y, x, i] = sum_value / data.shape[0] - -def _average_3d_stack_cuda(image_stack: np.ndarray) -> np.ndarray: - check_instance(source=_average_3d_stack_cuda.__name__, instance=image_stack, accepted_types=(np.ndarray,)) - check_if_valid_img(data=image_stack[0], source=_average_3d_stack_cuda.__name__) - if image_stack.ndim != 4: - return image_stack - x = np.ascontiguousarray(image_stack) - x_dev = cuda.to_device(x) - results = cuda.device_array((x.shape[1], x.shape[2], x.shape[3]), dtype=np.float32) - grid_x = (x.shape[1] + 16 - 1) // 16 - grid_y = (x.shape[2] + 16 - 1) // 16 - grid_z = 3 - threads_per_block = (16, 16, 1) - blocks_per_grid = (grid_y, grid_x, grid_z) - _average_3d_stack_cuda[blocks_per_grid, threads_per_block](x_dev, results) - results = results.copy_to_host() - return results - - - -def create_average_frm_cuda(video_path: Union[str, os.PathLike], - start_frm: Optional[int] = None, - end_frm: Optional[int] = None, - start_time: Optional[str] = None, - end_time: Optional[str] = None, - save_path: Optional[Union[str, os.PathLike]] = None, - batch_size: Optional[int] = 6000, - verbose: Optional[bool] = False) -> Union[None, np.ndarray]: - """ - Computes the average frame using GPU acceleration from a specified range of frames or time interval in a video file. - This average frame typically used for background substraction. - - - The function reads frames from the video, calculates their average, and optionally saves the result - to a specified file. If `save_path` is provided, the average frame is saved as an image file; - otherwise, the average frame is returned as a NumPy array. - - - :param Union[str, os.PathLike] video_path: The path to the video file from which to extract frames. - :param Optional[int] start_frm: The starting frame number (inclusive). Either `start_frm`/`end_frm` or `start_time`/`end_time` must be provided, but not both. - :param Optional[int] end_frm: The ending frame number (exclusive). - :param Optional[str] start_time: The start time in the format 'HH:MM:SS' from which to begin extracting frames. - :param Optional[str] end_time: The end time in the format 'HH:MM:SS' up to which frames should be extracted. - :param Optional[Union[str, os.PathLike]] save_path: The path where the average frame image will be saved. If `None`, the average frame is returned as a NumPy array. - :param Optional[int] batch_size: The number of frames to process in each batch. Default is 3000. Increase if your RAM allows it. - :param Optional[bool] verbose: If `True`, prints progress and informational messages during execution. - :return: Returns `None` if the result is saved to `save_path`. Otherwise, returns the average frame as a NumPy array. - - :example: - >>> create_average_frm(video_path=r"C:\troubleshooting\RAT_NOR\project_folder\videos\2022-06-20_NOB_DOT_4_downsampled.mp4", verbose=True, start_frm=0, end_frm=9000) - """ - - if not check_nvidea_gpu_available(): - raise FFMPEGCodecGPUError(msg="No GPU found (as evaluated by nvidea-smi returning None)", - source=create_average_frm_cuda.__name__) - - if ((start_frm is not None) or (end_frm is not None)) and ((start_time is not None) or (end_time is not None)): - raise InvalidInputError(msg=f'Pass start_frm and end_frm OR start_time and end_time', - source=create_average_frm_cuda.__name__) - elif type(start_frm) != type(end_frm): - raise InvalidInputError(msg=f'Pass start frame and end frame', source=create_average_frm_cuda.__name__) - elif type(start_time) != type(end_time): - raise InvalidInputError(msg=f'Pass start time and end time', source=create_average_frm_cuda.__name__) - if save_path is not None: - check_if_dir_exists(in_dir=os.path.dirname(save_path), source=create_average_frm_cuda.__name_) - check_file_exist_and_readable(file_path=video_path) - video_meta_data = get_video_meta_data(video_path=video_path) - video_name = get_fn_ext(filepath=video_path)[1] - if verbose: - print(f'Getting average frame from {video_name}...') - if (start_frm is not None) and (end_frm is not None): - check_int(name='start_frm', value=start_frm, min_value=0, max_value=video_meta_data['frame_count']) - check_int(name='end_frm', value=end_frm, min_value=0, max_value=video_meta_data['frame_count']) - if start_frm > end_frm: - raise InvalidInputError(msg=f'Start frame ({start_frm}) has to be before end frame ({end_frm}).', - source=create_average_frm_cuda.__name__) - frame_ids = list(range(start_frm, end_frm)) - elif (start_time is not None) and (end_time is not None): - check_if_string_value_is_valid_video_timestamp(value=start_time, name=create_average_frm_cuda.__name__) - check_if_string_value_is_valid_video_timestamp(value=end_time, name=create_average_frm_cuda.__name__) - check_that_hhmmss_start_is_before_end(start_time=start_time, end_time=end_time, - name=create_average_frm_cuda.__name__) - check_if_hhmmss_timestamp_is_valid_part_of_video(timestamp=start_time, video_path=video_path) - frame_ids = find_frame_numbers_from_time_stamp(start_time=start_time, end_time=end_time, - fps=video_meta_data['fps']) - else: - frame_ids = list(range(0, video_meta_data['frame_count'])) - frame_ids = [frame_ids[i:i + batch_size] for i in range(0, len(frame_ids), batch_size)] - avg_imgs = [] - for batch_cnt in range(len(frame_ids)): - start_idx, end_idx = frame_ids[batch_cnt][0], frame_ids[batch_cnt][-1] - if start_idx == end_idx: - continue - imgs = read_img_batch_from_video_gpu(video_path=video_path, start_frm=start_idx, end_frm=end_idx, verbose=verbose) - avg_imgs.append(_average_3d_stack_cuda(image_stack=np.stack(list(imgs.values()), axis=0))) - avg_img = _average_3d_stack_cuda(image_stack=np.stack(avg_imgs, axis=0)) - if save_path is not None: - cv2.imwrite(save_path, avg_img) - if verbose: - stdout_success(msg=f'Saved average frame at {save_path}', source=create_average_frm.__name__) - else: - return avg_img \ No newline at end of file diff --git a/simba/data_processors/cuda/create_shap_log.py b/simba/data_processors/cuda/create_shap_log.py deleted file mode 100644 index f58177a53..000000000 --- a/simba/data_processors/cuda/create_shap_log.py +++ /dev/null @@ -1,154 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -import os -from typing import List, Optional, Tuple, Union - -import numpy as np -import pandas as pd -import shap -from sklearn.ensemble import RandomForestClassifier - -from simba.mixins.train_model_mixin import TrainModelMixin -from simba.utils.checks import (check_if_dir_exists, check_instance, check_int, - check_nvidea_gpu_available, check_str, - check_valid_array, check_valid_dataframe, - check_valid_lst) -from simba.utils.enums import Formats -from simba.utils.errors import FFMPEGCodecGPUError -from simba.utils.printing import SimbaTimer, stdout_success -from simba.utils.read_write import read_df, write_df -from simba.utils.warnings import NotEnoughDataWarning - - -def create_shap_log(rf_clf: Union[str, os.PathLike, RandomForestClassifier], - x: Union[pd.DataFrame, np.ndarray], - y: Union[pd.DataFrame, pd.Series, np.ndarray], - cnt_present: int, - cnt_absent: int, - x_names: Optional[List[str]] = None, - clf_name: Optional[str] = None, - save_dir: Optional[Union[str, os.PathLike]] = None, - verbose: Optional[bool] = True) -> Union[None, Tuple[pd.DataFrame, pd.DataFrame, int]]: - """ - Computes SHAP (SHapley Additive exPlanations) values using a GPU for a RandomForestClassifier, - based on specified counts of positive and negative samples, and optionally saves the results. - - .. image:: _static/img/create_shap_log_cuda.png - :width: 500 - :align: center - - :param Union[str, os.PathLike, RandomForestClassifier] rf_clf: Trained RandomForestClassifier model or path to the saved model. Can be a string, os.PathLike object, or an instance of RandomForestClassifier. - :param Union[pd.DataFrame, np.ndarray] x: Input features used for SHAP value computation. Can be a pandas DataFrame or numpy ndarray. - :param Union[pd.DataFrame, pd.Series, np.ndarray] y: Target labels corresponding to the input features. Can be a pandas DataFrame, pandas Series, or numpy ndarray with 0 and 1 values. - :param int cnt_present: Number of positive samples (label=1) to include in the SHAP value computation. - :param int cnt_absent: Number of negative samples (label=0) to include in the SHAP value computation. - :param Optional[List[str]] x_names: Optional list of feature names corresponding to the columns in `x`. If `x` is a DataFrame, this is extracted automatically. - :param Optional[str] clf_name: Optional name for the classifier, used in naming output files. If not provided, it is extracted from the `y` labels if possible. - :param Optional[Union[str, os.PathLike]] save_dir: Optional directory path where the SHAP values and corresponding raw features are saved as CSV files. - :param Optional[bool] verbose: Optional boolean flag indicating whether to print progress messages. Defaults to True. - :return Union[None, Tuple[pd.DataFrame, pd.DataFrame, int]]: If `save_dir` is None, returns a tuple containing: - - V: DataFrame with SHAP values, expected value, sum of SHAP values, prediction probability, and target labels. - - R: DataFrame containing the raw feature values for the selected samples. - - expected_value: The expected value from the SHAP explainer. - If `save_dir` is provided, the function returns None and saves the output to CSV files in the specified directory. - - :example: - >>> x = np.random.random((1000, 501)).astype(np.float32) - >>> y = np.random.randint(0, 2, size=(len(x), 1)).astype(np.int32) - >>> clf_names = [str(x) for x in range(501)] - >>> results = create_shap_log(rf_clf=MODEL_PATH, x=x, y=y, cnt_present=int(i/2), cnt_absent=int(i/2), clf_name='TEST', x_names=clf_names, verbose=False) - """ - - timer = SimbaTimer(start=True) - if verbose: - print('Computing SHAP values (GPU)...') - if not check_nvidea_gpu_available(): - raise FFMPEGCodecGPUError(msg="No GPU found (as evaluated by nvidea-smi returning None)", - source=create_shap_log.__name__) - check_instance(source=f'{create_shap_log.__name__} rf_clf', instance=rf_clf, - accepted_types=(str, RandomForestClassifier)) - if isinstance(rf_clf, (str, os.PathLike)): - rf_clf = TrainModelMixin().read_pickle(file_path=rf_clf) - check_instance(source=f'{create_shap_log.__name__} x', instance=x, accepted_types=(pd.DataFrame, np.ndarray)) - if isinstance(x, np.ndarray): - check_valid_lst(data=x_names, source=f'{create_shap_log.__name__} x_names', valid_dtypes=(str,), - exact_len=x.shape[1]) - check_valid_array(data=x, source=f'{create_shap_log.__name__} x', accepted_ndims=[2, ], - accepted_dtypes=Formats.NUMERIC_DTYPES.value) - else: - check_valid_dataframe(df=x, source=f'{create_shap_log.__name__} x', - valid_dtypes=Formats.NUMERIC_DTYPES.value) - x_names = list(x.columns) - x = x.values - check_instance(source=f'{create_shap_log.__name__} y', instance=y, - accepted_types=(pd.DataFrame, np.ndarray, pd.Series)) - if isinstance(y, np.ndarray): - check_str(name=f'{create_shap_log.__name__} clf_name', value=clf_name) - y = y.flatten() - elif isinstance(y, pd.Series): - clf_name = y.name - y = y.values.flatten() - else: - check_valid_dataframe(df=y, source=f'{create_shap_log.__name__} y', - valid_dtypes=Formats.NUMERIC_DTYPES.value, max_axis_1=1) - clf_name = list(y.columns)[0] - y = y.values.flatten() - save_shap_path, save_raw_path = None, None - if save_dir is not None: - check_if_dir_exists(in_dir=save_dir) - save_shap_path = os.path.join(save_dir, f"SHAP_values_{clf_name}.csv") - save_raw_path = os.path.join(save_dir, f"RAW_SHAP_feature_values_{clf_name}.csv") - check_valid_array(data=y, source=f'{create_shap_log.__name__} y', accepted_values=[0, 1]) - check_int(name=f'{create_shap_log.__name__} cnt_present', value=cnt_present, min_value=1) - check_int(name=f'{create_shap_log.__name__} cnt_absent', value=cnt_absent, min_value=1) - target_cnt = np.sum(y) - absent_cnt = y.shape[0] - target_cnt - - if cnt_present > target_cnt: - NotEnoughDataWarning( - msg=f"Data contains {target_cnt} behavior-present annotations. This is less the number of frames you specified to calculate shap values for ({cnt_present}). SimBA will calculate shap scores for the {target_cnt} behavior-present frames available", - source=create_shap_log.__name__) - cnt_present = target_cnt - if absent_cnt < cnt_absent: - NotEnoughDataWarning( - msg=f"Data contains {absent_cnt} behavior-absent annotations. This is less the number of frames you specified to calculate shap values for ({cnt_absent}). SimBA will calculate shap scores for the {absent_cnt} behavior-absent frames available", - source=create_shap_log.__name__) - cnt_absent = absent_cnt - - target_idx = np.argwhere(y == 1).flatten() - absent_idx = np.argwhere(y == 0).flatten() - target_idx = np.sort(np.random.choice(target_idx, cnt_present)) - absent_idx = np.sort(np.random.choice(absent_idx, cnt_absent)) - target_x = x[target_idx] - absent_x = x[absent_idx] - X = np.vstack([target_x, absent_x]).astype(np.float32) - Y = np.hstack([np.ones(target_x.shape[0]), np.zeros(absent_x.shape[0])]).astype(np.int32) - explainer = shap.explainers.GPUTree(model=rf_clf, data=None, model_output='raw', - feature_names='tree_path_dependent') - shap_values = explainer.shap_values(X, check_additivity=True) - V = pd.DataFrame(shap_values[1], columns=x_names).astype(np.float32) - sum = V.sum(axis=1) - expected_value = explainer.expected_value[1] - p = TrainModelMixin().clf_predict_proba(clf=rf_clf, x_df=X) - - V['EXPECTED_VALUE'] = expected_value.round(4) - V['SUM'] = sum + V['EXPECTED_VALUE'] - V['PREDICTION_PROBABILITY'] = p.round(4) - V['SUM'] = V['SUM'].round(4) - V[clf_name] = Y - x_idx = np.hstack([target_idx, absent_idx]) - R = pd.DataFrame(x[x_idx, :], columns=x_names) - timer.stop_timer() - if save_dir is None: - if verbose: - stdout_success(msg=f'Shap values compute complete (GPU) for {len(V)} observations.', elapsed_time=timer.elapsed_time_str) - return (V, R, expected_value) - else: - write_df(df=V, file_type='csv', save_path=save_shap_path) - write_df(df=R, file_type='csv', save_path=save_raw_path) - if verbose: - stdout_success(msg=f'Shap values compute complete (GPU) for {len(V)} observations, and saved in {save_dir}', elapsed_time=timer.elapsed_time_str) - - - diff --git a/simba/data_processors/cuda/direction_two_bps.py b/simba/data_processors/cuda/direction_two_bps.py deleted file mode 100644 index 4dee71a60..000000000 --- a/simba/data_processors/cuda/direction_two_bps.py +++ /dev/null @@ -1,45 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -import math - -import numpy as np -from numba import cuda, int32 - -THREADS_PER_BLOCK = 1024 - -@cuda.jit() -def _cuda_direction_from_two_bps(x, y, results): - i = cuda.grid(1) - if i > x.shape[0]: - return - else: - a = math.atan2(x[i][0] - y[i][0], y[i][1] - x[i][1]) * (180 / math.pi) - a = int32(a + 360 if a < 0 else a) - results[i] = a - - -def direction_from_two_bps(x: np.ndarray, y: np.ndarray) -> np.ndarray: - """ - Compute the directionality in degrees from two body-parts. E.g., ``nape`` and ``nose``, - or ``swim_bladder`` and ``tail`` with GPU acceleration. - - .. image:: _static/img/direction_from_two_bps_cuda.png - :width: 1200 - :align: center - - - :parameter np.ndarray x: Size len(frames) x 2 representing x and y coordinates for first body-part. - :parameter np.ndarray y: Size len(frames) x 2 representing x and y coordinates for second body-part. - :return np.ndarray: Frame-wise directionality in degrees. - - """ - x = np.ascontiguousarray(x).astype(np.int32) - y = np.ascontiguousarray(y).astype(np.int32) - x_dev = cuda.to_device(x) - y_dev = cuda.to_device(y) - results = cuda.device_array((x.shape[0]), dtype=np.int32) - bpg = (x.shape[0] + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK - _cuda_direction_from_two_bps[bpg, THREADS_PER_BLOCK](x_dev, y_dev, results) - results = results.copy_to_host() - return results diff --git a/simba/data_processors/cuda/euclidan_distance_cuda.py b/simba/data_processors/cuda/euclidan_distance_cuda.py deleted file mode 100644 index dae3f2880..000000000 --- a/simba/data_processors/cuda/euclidan_distance_cuda.py +++ /dev/null @@ -1,52 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -import math - -import numpy as np -from numba import cuda - -from simba.mixins.feature_extraction_mixin import FeatureExtractionMixin -from simba.utils.read_write import read_df - -THREADS_PER_BLOCK = 128 - -@cuda.jit -def _euclidean_distance_kernel(x_dev, y_dev, results): - i = cuda.grid(1) - if i < x_dev.shape[0]: - p = (math.sqrt((x_dev[i][0] - y_dev[i][0]) ** 2 + (x_dev[i][1] - y_dev[i][1]) ** 2)) - results[i] = p - -def get_euclidean_distance_cuda(x: np.ndarray, y: np.ndarray) -> np.ndarray: - """ - Computes the Euclidean distance between two sets of points using CUDA for GPU acceleration. - - .. image:: _static/img/get_euclidean_distance_cuda.png - :width: 500 - :align: center - - :param np.ndarray x: A 2D array of shape (n, m) representing n points in m-dimensional space. Each row corresponds to a point. - :param np.ndarray y: A 2D array of shape (n, m) representing n points in m-dimensional space. Each row corresponds to a point. - :return np.ndarray: A 1D array of shape (n,) where each element represents the Euclidean distance between the corresponding points in `x` and `y`. - - :example: - >>> video_path = r"/mnt/c/troubleshooting/mitra/project_folder/videos/501_MA142_Gi_CNO_0514.mp4" - >>> data_path = r"/mnt/c/troubleshooting/mitra/project_folder/csv/outlier_corrected_movement_location/501_MA142_Gi_CNO_0514 - test.csv" - >>> df = read_df(file_path=data_path, file_type='csv')[['Center_x', 'Center_y']] - >>> shifted_df = FeatureExtractionMixin.create_shifted_df(df=df, periods=1) - >>> x = shifted_df[['Center_x', 'Center_y']].values - >>> y = shifted_df[['Center_x_shifted', 'Center_y_shifted']].values - >>> get_euclidean_distance_cuda(x=x, y=y) - """ - - x = np.ascontiguousarray(x).astype(np.int32) - y = np.ascontiguousarray(y).astype(np.int32) - n, m = x.shape - x_dev = cuda.to_device(x) - y_dev = cuda.to_device(y) - results = cuda.device_array((n, m), dtype=np.int32) - bpg = (n + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK - _euclidean_distance_kernel[bpg, THREADS_PER_BLOCK](x_dev, y_dev, results) - results = results.copy_to_host().astype(np.int32) - return results diff --git a/simba/data_processors/cuda/euclidean_distances_cuda.py b/simba/data_processors/cuda/euclidean_distances_cuda.py deleted file mode 100644 index b4bc909d4..000000000 --- a/simba/data_processors/cuda/euclidean_distances_cuda.py +++ /dev/null @@ -1,38 +0,0 @@ -from typing import Optional - -import cupy as cp -import numpy as np - -from simba.utils.checks import check_int, check_valid_array -from simba.utils.enums import Formats - - -def get_euclidean_distance_cupy(x: np.ndarray, - y: np.ndarray, - batch_size: Optional[int] = int(3.5e10+7)) -> np.ndarray: - """ - Computes the Euclidean distance between corresponding pairs of points in two 2D arrays - using CuPy for GPU acceleration. The computation is performed in batches to handle large - datasets efficiently. - - :param np.ndarray x: A 2D NumPy array with shape (n, 2), where each row represents a point in a 2D space. - :param np.ndarray y: A 2D NumPy array with shape (n, 2), where each row represents a point in a 2D space. The shape of `y` must match the shape of `x`. - :param Optional[int] batch_size: The number of points to process in a single batch. This parameter controls memory usage and can be adjusted based on available GPU memory. The default value is large (`3.5e10 + 7`) to maximize GPU utilization, but it can be lowered if memory issues arise. - :return: A 1D NumPy array of shape (n,) containing the Euclidean distances between corresponding points in `x` and `y`. - :return: A 1D NumPy array of shape (n,) containing the Euclidean distances between corresponding points in `x` and `y`. - :rtype: np.ndarray - - :example: - >>> x = np.array([[1, 2], [3, 4], [5, 6]]) - >>> y = np.array([[7, 8], [9, 10], [11, 12]]) - >>> distances = get_euclidean_distance_cupy(x, y) - """ - check_valid_array(data=x, source=check_valid_array.__name__, accepted_ndims=[2,], accepted_dtypes=Formats.NUMERIC_DTYPES.value) - check_valid_array(data=y, source=check_valid_array.__name__, accepted_ndims=[2, ], accepted_dtypes=Formats.NUMERIC_DTYPES.value, accepted_shapes=(x.shape,)) - check_int(name='batch_size', value=batch_size, min_value=1) - results = cp.full((x.shape[0]), fill_value=cp.nan, dtype=cp.float32) - for l in range(0, x.shape[0], batch_size): - r = l + batch_size - batch_x, batch_y = cp.array(x[l:r]), cp.array(y[l:r]) - results[l:r] = (cp.sqrt((batch_x[:, 0] - batch_y[:, 0]) ** 2 + (batch_x[:, 1] - batch_y[:, 1]) ** 2)) - return results.get() \ No newline at end of file diff --git a/simba/data_processors/cuda/image.py b/simba/data_processors/cuda/image.py index af1916896..39cbc1e9d 100644 --- a/simba/data_processors/cuda/image.py +++ b/simba/data_processors/cuda/image.py @@ -4,7 +4,7 @@ import os from typing import Optional, Union - +import math try: from typing import Literal except: @@ -30,10 +30,9 @@ from simba.utils.data import find_frame_numbers_from_time_stamp from simba.utils.enums import Formats from simba.utils.errors import FFMPEGCodecGPUError, InvalidInputError -from simba.utils.printing import stdout_success -from simba.utils.read_write import ( - check_if_hhmmss_timestamp_is_valid_part_of_video, get_fn_ext, - get_video_meta_data, read_img_batch_from_video_gpu) +from simba.utils.printing import stdout_success, SimbaTimer +from simba.utils.read_write import (check_if_hhmmss_timestamp_is_valid_part_of_video, get_fn_ext, get_video_meta_data, read_img_batch_from_video_gpu) +from simba.mixins.image_mixin import ImageMixin PHOTOMETRIC = 'photometric' DIGITAL = 'digital' @@ -670,4 +669,170 @@ def segment_img_stack_horizontal(imgs: np.ndarray, else: imgs = imgs[:, int(px_crop/2):int((imgs.shape[0] - px_crop) / 2), :] - return imgs.get() \ No newline at end of file + return imgs.get() + + + +@cuda.jit(device=True) +def _cuda_is_inside_polygon(x, y, polygon_vertices): + """ + Checks if the pixel location is inside the polygon. + + :param int x: Pixel x location. + :param int y: Pixel y location. + :param np.ndarray polygon_vertices: 2-dimentional array representing the x and y coordinates of the polygon vertices. + :return: Boolean representing if the x and y are located in the polygon. + """ + + n = len(polygon_vertices) + p2x, p2y, xints, inside = 0.0, 0.0, 0.0, False + p1x, p1y = polygon_vertices[0] + for j in range(n + 1): + p2x, p2y = polygon_vertices[j % n] + if ( + (y > min(p1y, p2y)) + and (y <= max(p1y, p2y)) + and (x <= max(p1x, p2x)) + ): + if p1y != p2y: + xints = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x + if p1x == p2x or x <= xints: + inside = not inside + p1x, p1y = p2x, p2y + return inside + + + +@cuda.jit(device=True) +def _cuda_is_inside_circle(x, y, circle_x, circle_y, circle_r): + """ + Device func to check if the pixel location is inside a circle. + + :param int x: Pixel x location. + :param int y: Pixel y location. + :param int circle_x: Center of circle x coordinate. + :param int circle_y: Center of circle y coordinate. + :param int y: Circle radius. + :return: Boolean representing if the x and y are located in the circle. + """ + + p = (math.sqrt((x - circle_x) ** 2 + (y - circle_y) ** 2)) + if p <= circle_r: + return True + else: + return False + +@cuda.jit() +def _cuda_create_rectangle_masks(shapes, imgs): + x, y, n = cuda.grid(3) + if n < 0 or n > (imgs.shape[0] -1): + return + if y < 0 or y > (imgs.shape[1] -1): + return + if x < 0 or x > (imgs.shape[2] -1): + return + else: + polygon = shapes[n] + inside = _cuda_is_inside_polygon(x, y, polygon) + if not inside: + imgs[n, y, x] = 0 + +@cuda.jit() +def _cuda_create_circle_masks(shapes, imgs): + x, y, n = cuda.grid(3) + if n < 0 or n > (imgs.shape[0] -1): + return + if y < 0 or y > (imgs.shape[1] -1): + return + if x < 0 or x > (imgs.shape[2] -1): + return + else: + circle_x, circle_y, circle_r = shapes[n][0], shapes[n][1], shapes[n][2] + inside = _cuda_is_inside_circle(x, y, circle_x, circle_y, circle_r) + if not inside: + imgs[n, y, x] = 0 + +def slice_imgs(video_path: Union[str, os.PathLike], + shapes: np.ndarray, + batch_size: Optional[int] = 1000, + verbose: Optional[bool] = True): + + """ + Slice frames from a video based on given shape coordinates (rectangles or circles) and return the cropped regions using GPU acceleration. + + + .. video:: _static/img/slice_imgs_gpu.webm + :width: 800 + :autoplay: + :loop: + + .. csv-table:: + :header: EXPECTED RUNTIMES + :file: ../../../docs/tables/slice_imgs.csv + :widths: 10, 90 + :align: center + :class: simba-table + :header-rows: 1 + + :param Union[str, os.PathLike] video_path: Path to the video file. + :param np.ndarray shapes: A NumPy array of shape `(n, m, 2)` where `n` is the number of frames. Each frame contains 4 (x, y) or one points representing the bounding shapes (e.g., rectangles edges) or centroid and radius (if circles) to slice from each frame. + :param Optional[int] batch_size: Optional; default is 500. The number of frames to process in each batch for memory efficiency. Larger batches are faster but use more memory. + :param Optional[bool] verbose: If True, prints progress during the slicing process. + :return: A NumPy array of sliced images with shape `(n, h, w, 3)` if the video is in color, or `(n, h, w)` if the video is grayscale. Here, `n` is the number of frames, and `h`, `w` are the height and width of the frames, respectively. + :rtype: np.ndarray + + :example I rectangles: + >>> data_path = r"/mnt/c/troubleshooting/mitra/project_folder/csv/outlier_corrected_movement_location/FRR_gq_Saline_0624.csv" # PATH TO A DATA FILE + >>> video_path = r'/mnt/c/troubleshooting/mitra/project_folder/videos/FRR_gq_Saline_0624.mp4' # PATH TO AN ASSOCIATED VIDEO FILE + >>> nose_arr = read_df(file_path=data_path, file_type='csv', usecols=['Nose_x', 'Nose_y', 'Tail_base_x', 'Tail_base_y', 'Left_side_x', 'Left_side_y', 'Right_side_x', 'Right_side_y']).values.reshape(-1, 4, 2)[0:1000] ## READ THE BODY-PART THAT DEFINES THE HULL AND CONVERT TO ARRAY + >>> polygons = GeometryMixin().multiframe_bodyparts_to_polygon(data=nose_arr, parallel_offset=60) ## CONVERT THE BODY-PART TO POLYGONS WITH A LITTLE BUFFER + >>> polygons = GeometryMixin().multiframe_minimum_rotated_rectangle(shapes=polygons) # CONVERT THE POLYGONS TO RECTANGLES (I.E., WITH 4 UNIQUE POINTS). + >>> polygon_lst = [] # GET THE POINTS OF THE RECTANGLES + >>> for i in polygons: polygon_lst.append(np.array(i.exterior.coords)) + >>> polygons = np.stack(polygon_lst, axis=0) + >>> sliced_imgs = slice_imgs(video_path=video_path, shapes=polygons) #SLICE THE RECTANGLES IN THE VIDEO. + + + :example II circles: + >>> data_path = r"/mnt/c/troubleshooting/mitra/project_folder/csv/outlier_corrected_movement_location/FRR_gq_Saline_0624.csv" # PATH TO A DATA FILE + >>> video_path = r'/mnt/c/troubleshooting/mitra/project_folder/videos/FRR_gq_Saline_0624.mp4' # PATH TO AN ASSOCIATED VIDEO FILE + >>> nose_arr = read_df(file_path=data_path, file_type='csv', usecols=['Nose_x', 'Nose_y']).values[0:6000] ## READ THE BODY-PART THAT DEFINES THE CENTER OF CIRCLE + >>> nose_arr = np.hstack((nose_arr, np.full((nose_arr.shape[0], 1), fill_value=50))).astype(np.int32) ## APPEND THE RADIUS OF THE CIRCLE TO THE DATA + >>> sliced_imgs = slice_imgs(video_path=video_path, shapes=nose_arr) #SLICE THE CIRCLE IN THE VIDEO. + """ + THREADS_PER_BLOCK = (32, 32, 1) + + video_meta_data = get_video_meta_data(video_path=video_path) + video_meta_data['frame_count'] = shapes.shape[0] + n, w, h = video_meta_data['frame_count'], video_meta_data['width'], video_meta_data['height'] + is_color = ImageMixin.is_video_color(video=video_path) + timer = SimbaTimer(start=True) + if is_color: + results = np.zeros((n, h, w, 3), dtype=np.uint8) + else: + results = np.zeros((n, h, w), dtype=np.uint8) + for start_img_idx in range(0, n, batch_size): + end_img_idx = start_img_idx + batch_size + if end_img_idx > video_meta_data['frame_count']: + end_img_idx = video_meta_data['frame_count'] + if verbose: + print(f'Processing images {start_img_idx} to {end_img_idx} (of {n})...') + batch_n = end_img_idx - start_img_idx + batch_imgs = read_img_batch_from_video_gpu(video_path=video_path, start_frm=start_img_idx, end_frm=end_img_idx) + batch_imgs = np.stack(list(batch_imgs.values()), axis=0) + batch_shapes = shapes[start_img_idx:end_img_idx].astype(np.int32) + x_dev = cuda.to_device(batch_shapes) + batch_img_dev = cuda.to_device(batch_imgs) + grid_x = math.ceil(w / THREADS_PER_BLOCK[0]) + grid_y = math.ceil(h / THREADS_PER_BLOCK[1]) + grid_z = math.ceil(batch_n / THREADS_PER_BLOCK[2]) + bpg = (grid_x, grid_y, grid_z) + if batch_shapes.shape[1] == 3: + _cuda_create_circle_masks[bpg, THREADS_PER_BLOCK](x_dev, batch_img_dev) + else: + _cuda_create_rectangle_masks[bpg, THREADS_PER_BLOCK](x_dev, batch_img_dev) + results[start_img_idx: end_img_idx] = batch_img_dev.copy_to_host() + timer.stop_timer() + if verbose: + stdout_success(msg='Shapes sliced in video.', elapsed_time=timer.elapsed_time_str) + return results diff --git a/simba/data_processors/cuda/img_stack_brightness.py b/simba/data_processors/cuda/img_stack_brightness.py deleted file mode 100644 index 8922d8863..000000000 --- a/simba/data_processors/cuda/img_stack_brightness.py +++ /dev/null @@ -1,95 +0,0 @@ -import time -from copy import deepcopy -from typing import Optional, Union - -import numpy as np - -try: - from typing import Literal -except: - from typing_extensions import Literal - -from numba import cuda - -from simba.utils.checks import check_if_valid_img, check_instance -from simba.utils.read_write import read_img_batch_from_video_gpu - -PHOTOMETRIC = 'photometric' -DIGITAL = 'digital' - -@cuda.jit() -def _photometric(data, results): - y, x, i = cuda.grid(3) - if i < 0 or x < 0 or y < 0: - return - if i > results.shape[0] - 1 or x > results.shape[1] - 1 or y > results.shape[2] - 1: - return - else: - r, g, b = data[i][x][y][0], data[i][x][y][1], data[i][x][y][2] - results[i][x][y] = (0.2126 * r) + (0.7152 * g) + (0.0722 * b) - -@cuda.jit() -def _digital(data, results): - y, x, i = cuda.grid(3) - if i < 0 or x < 0 or y < 0: - return - if i > results.shape[0] - 1 or x > results.shape[1] - 1 or y > results.shape[2] - 1: - return - else: - r, g, b = data[i][x][y][0], data[i][x][y][1], data[i][x][y][2] - results[i][x][y] = (0.299 * r) + (0.587 * g) + (0.114 * b) - -def img_stack_brightness(x: np.ndarray, - method: Optional[Literal['photometric', 'digital']] = 'digital', - ignore_black: Optional[bool] = True) -> np.ndarray: - """ - Calculate the average brightness of a stack of images using a specified method. - - - - **Photometric Method**: The brightness is calculated using the formula: - - .. math:: - \text{brightness} = 0.2126 \cdot R + 0.7152 \cdot G + 0.0722 \cdot B - - - **Digital Method**: The brightness is calculated using the formula: - - .. math:: - \text{brightness} = 0.299 \cdot R + 0.587 \cdot G + 0.114 \cdot B - - - :param np.ndarray x: A 4D array of images with dimensions (N, H, W, C), where N is the number of images, H and W are the height and width, and C is the number of channels (RGB). - :param Optional[Literal['photometric', 'digital']] method: The method to use for calculating brightness. It can be 'photometric' for the standard luminance calculation or 'digital' for an alternative set of coefficients. Default is 'digital'. - :param Optional[bool] ignore_black: If True, black pixels (i.e., pixels with brightness value 0) will be ignored in the calculation of the average brightness. Default is True. - :return np.ndarray: A 1D array of average brightness values for each image in the stack. If `ignore_black` is True, black pixels are ignored in the averaging process. - - - :example: - >>> imgs = read_img_batch_from_video_gpu(video_path=r"/mnt/c/troubleshooting/RAT_NOR/project_folder/videos/2022-06-20_NOB_DOT_4_downsampled.mp4", start_frm=0, end_frm=5000) - >>> imgs = np.stack(list(imgs.values()), axis=0) - >>> x = img_stack_brightness(x=imgs) - """ - - check_instance(source=img_stack_brightness.__name__, instance=x, accepted_types=(np.ndarray,)) - check_if_valid_img(data=x[0], source=img_stack_brightness.__name__) - x = np.ascontiguousarray(x).astype(np.uint8) - if x.ndim == 4: - grid_x = (x.shape[1] + 16 - 1) // 16 - grid_y = (x.shape[2] + 16 - 1) // 16 - grid_z = x.shape[0] - threads_per_block = (16, 16, 1) - blocks_per_grid = (grid_y, grid_x, grid_z) - x_dev = cuda.to_device(x) - results = cuda.device_array((x.shape[0], x.shape[1], x.shape[2]), dtype=np.uint8) - if method == PHOTOMETRIC: - _photometric[blocks_per_grid, threads_per_block](x_dev, results) - else: - _digital[blocks_per_grid, threads_per_block](x_dev, results) - results = results.copy_to_host() - if ignore_black: - masked_array = np.ma.masked_equal(results, 0) - results = np.mean(masked_array, axis=(1, 2)).filled(0) - else: - results = deepcopy(x) - results = np.mean(results, axis=(1, 2)) - - return results \ No newline at end of file diff --git a/simba/data_processors/cuda/img_stack_mse.py b/simba/data_processors/cuda/img_stack_mse.py deleted file mode 100644 index 78b18b56d..000000000 --- a/simba/data_processors/cuda/img_stack_mse.py +++ /dev/null @@ -1,106 +0,0 @@ -from typing import Optional - -import numpy as np -from numba import cuda - -from simba.utils.checks import (check_if_valid_img, check_instance, - check_valid_array) - - -@cuda.jit() -def _grey_mse(data, ref_img, stride, batch_cnt, mse_arr): - y, x, i = cuda.grid(3) - stride = stride[0] - batch_cnt = batch_cnt[0] - if batch_cnt == 0: - if (i - stride) < 0 or x < 0 or y < 0: - return - else: - if i < 0 or x < 0 or y < 0: - return - if i > mse_arr.shape[0] - 1 or x > mse_arr.shape[1] - 1 or y > mse_arr.shape[2] - 1: - return - else: - img_val = data[i][x][y] - if i == 0: - prev_val = ref_img[x][y] - else: - img_val = data[i][x][y] - prev_val = data[i - stride][x][y] - mse_arr[i][x][y] = (img_val - prev_val) ** 2 - - -@cuda.jit() -def _rgb_mse(data, ref_img, stride, batch_cnt, mse_arr): - y, x, i = cuda.grid(3) - stride = stride[0] - batch_cnt = batch_cnt[0] - if batch_cnt == 0: - if (i - stride) < 0 or x < 0 or y < 0: - return - else: - if i < 0 or x < 0 or y < 0: - return - if i > mse_arr.shape[0] - 1 or x > mse_arr.shape[1] - 1 or y > mse_arr.shape[2] - 1: - return - else: - img_val = data[i][x][y] - if i != 0: - prev_val = data[i - stride][x][y] - else: - prev_val = ref_img[x][y] - r_diff = (img_val[0] - prev_val[0]) ** 2 - g_diff = (img_val[1] - prev_val[1]) ** 2 - b_diff = (img_val[2] - prev_val[2]) ** 2 - mse_arr[i][x][y] = r_diff + g_diff + b_diff - -def stack_sliding_mse(x: np.ndarray, - stride: Optional[int] = 1, - batch_size: Optional[int] = 1000) -> np.ndarray: - """ - Computes the Mean Squared Error (MSE) between each image in a stack and a reference image, - where the reference image is determined by a sliding window approach with a specified stride. - The function is optimized for large image stacks by processing them in batches. - - :param np.ndarray x: Input array of images, where the first dimension corresponds to the stack of images. The array should be either 3D (height, width, channels) or 4D (batch, height, width, channels). - :param Optional[int] stride: The stride or step size for the sliding window that determines the reference image. Defaults to 1, meaning the previous image in the stack is used as the reference. - :param Optional[int] batch_size: The number of images to process in a single batch. Larger batch sizes may improve performance but require more GPU memory. Defaults to 1000. - :return: A 1D NumPy array containing the MSE for each image in the stack compared to its corresponding reference image. The length of the array is equal to the number of images in the input stack. - :rtype: np.ndarray - - """ - - check_instance(source=stack_sliding_mse.__name__, instance=x, accepted_types=(np.ndarray,)) - check_if_valid_img(data=x[0], source=stack_sliding_mse.__name__) - check_valid_array(data=x, source=stack_sliding_mse.__name__, accepted_ndims=[3, 4]) - stride = np.array([stride], dtype=np.int32) - stride_dev = cuda.to_device(stride) - out = np.full((x.shape[0]), fill_value=0.0, dtype=np.float32) - for batch_cnt, l in enumerate(range(0, x.shape[0], batch_size)): - r = l + batch_size - batch_x = x[l:r] - if batch_cnt != 0: - if x.ndim == 3: - ref_img = x[l-stride].astype(np.uint8).reshape(x.shape[1], x.shape[2]) - else: - ref_img = x[l-stride].astype(np.uint8).reshape(x.shape[1], x.shape[2], 3) - else: - ref_img = np.full_like(x[l], dtype=np.uint8, fill_value=0) - ref_img = ref_img.astype(np.uint8) - grid_x = (batch_x.shape[1] + 16 - 1) // 16 - grid_y = (batch_x.shape[2] + 16 - 1) // 16 - grid_z = batch_x.shape[0] - threads_per_block = (16, 16, 1) - blocks_per_grid = (grid_y, grid_x, grid_z) - ref_img_dev = cuda.to_device(ref_img) - x_dev = cuda.to_device(batch_x) - results = cuda.device_array((batch_x.shape[0], batch_x.shape[1], batch_x.shape[2]), dtype=np.uint8) - batch_cnt_dev = np.array([batch_cnt], dtype=np.int32) - if x.ndim == 3: - _grey_mse[blocks_per_grid, threads_per_block](x_dev, ref_img_dev, stride_dev, batch_cnt_dev, results) - else: - _rgb_mse[blocks_per_grid, threads_per_block](x_dev, ref_img_dev, stride_dev, batch_cnt_dev, results) - results = results.copy_to_host() - results = np.mean(results, axis=(1, 2)) - out[l:r] = results - return out \ No newline at end of file diff --git a/simba/data_processors/cuda/imgs_to_grayscale_cupy.py b/simba/data_processors/cuda/imgs_to_grayscale_cupy.py deleted file mode 100644 index 5815174b6..000000000 --- a/simba/data_processors/cuda/imgs_to_grayscale_cupy.py +++ /dev/null @@ -1,44 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -from typing import Optional - -import cupy as cp -import numpy as np - -from simba.utils.checks import check_if_valid_img, check_instance -from simba.utils.read_write import read_img_batch_from_video_gpu - - -def img_stack_to_grayscale_cupy(imgs: np.ndarray, - batch_size: Optional[int] = 250) -> np.ndarray: - """ - Converts a stack of color images to grayscale using GPU acceleration with CuPy. - - :param np.ndarray imgs: A 4D NumPy array representing a stack of images with shape (num_images, height, width, channels). The images are expected to have 3 channels (RGB). - :param Optional[int] batch_size: The number of images to process in each batch. Defaults to 250. Adjust this parameter to fit your GPU's memory capacity. - :return np.ndarray: m A 3D NumPy array of shape (num_images, height, width) containing the grayscale images. If the input array is not 4D, the function returns the input as is. - - :example: - >>> imgs = read_img_batch_from_video_gpu(video_path=r"/mnt/c/troubleshooting/RAT_NOR/project_folder/videos/2022-06-20_NOB_IOT_1_cropped.mp4", verbose=False, start_frm=0, end_frm=i) - >>> imgs = np.stack(list(imgs.values()), axis=0).astype(np.uint8) - >>> gray_imgs = img_stack_to_grayscale_cupy(imgs=imgs) - """ - - - check_instance(source=img_stack_to_grayscale_cupy.__name__, instance=imgs, accepted_types=(np.ndarray,)) - check_if_valid_img(data=imgs[0], source=img_stack_to_grayscale_cupy.__name__) - if imgs.ndim != 4: - return imgs - results = cp.zeros((imgs.shape[0], imgs.shape[1], imgs.shape[2]), dtype=np.uint8) - n = int(np.ceil((imgs.shape[0] / batch_size))) - imgs = np.array_split(imgs, n) - start = 0 - for i in range(len(imgs)): - img_batch = cp.array(imgs[i]) - batch_cnt = img_batch.shape[0] - end = start + batch_cnt - vals = (0.07 * img_batch[:, :, :, 2] + 0.72 * img_batch[:, :, :, 1] + 0.21 * img_batch[:, :, :, 0]) - results[start:end] = vals.astype(cp.uint8) - start = end - return results.get() \ No newline at end of file diff --git a/simba/data_processors/cuda/imgs_to_greyscale_cuda.py b/simba/data_processors/cuda/imgs_to_greyscale_cuda.py deleted file mode 100644 index fea95a7bd..000000000 --- a/simba/data_processors/cuda/imgs_to_greyscale_cuda.py +++ /dev/null @@ -1,53 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -import numpy as np -from numba import cuda - -from simba.utils.checks import check_if_valid_img, check_instance -from simba.utils.read_write import read_img_batch_from_video_gpu - - -@cuda.jit() -def _img_stack_to_grayscale(data, results): - y, x, i = cuda.grid(3) - if i < 0 or x < 0 or y < 0: - return - if i > results.shape[0] - 1 or x > results.shape[1] - 1 or y > results.shape[2] - 1: - return - else: - b = 0.07 * data[i][x][y][2] - g = 0.72 * data[i][x][y][1] - r = 0.21 * data[i][x][y][0] - val = b + g + r - results[i][x][y] = val - -def img_stack_to_grayscale_cuda(x: np.ndarray) -> np.ndarray: - """ - Convert image stack to grayscale using CUDA. - - :param np.ndarray x: 4d array of color images in numpy format. - :return np.ndarray: 3D array of greyscaled images. - - :example: - >>> imgs = read_img_batch_from_video_gpu(video_path=r"/mnt/c/troubleshooting/mitra/project_folder/videos/temp_2/592_MA147_Gq_Saline_0516_downsampled.mp4", verbose=False, start_frm=0, end_frm=i) - >>> imgs = np.stack(list(imgs.values()), axis=0).astype(np.uint8) - >>> grey_images = img_stack_to_grayscale_cuda(x=imgs) - """ - check_instance(source=img_stack_to_grayscale_cuda.__name__, instance=imgs, accepted_types=(np.ndarray,)) - check_if_valid_img(data=x[0], source=img_stack_to_grayscale_cuda.__name__) - if x.ndim != 4: - return x - x = np.ascontiguousarray(x).astype(np.uint8) - x_dev = cuda.to_device(x) - results = cuda.device_array((x.shape[0], x.shape[1], x.shape[2]), dtype=np.uint8) - grid_x = (x.shape[1] + 16 - 1) // 16 - grid_y = (x.shape[2] + 16 - 1) // 16 - grid_z = x.shape[0] - threads_per_block = (16, 16, 1) - blocks_per_grid = (grid_y, grid_x, grid_z) - _img_stack_to_grayscale[blocks_per_grid, threads_per_block](x_dev, results) - results = results.copy_to_host() - return results - - diff --git a/simba/data_processors/cuda/is_inside_circle.py b/simba/data_processors/cuda/is_inside_circle.py deleted file mode 100644 index 6d18445c1..000000000 --- a/simba/data_processors/cuda/is_inside_circle.py +++ /dev/null @@ -1,39 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -import math - -import numpy as np -from numba import cuda - -THREADS_PER_BLOCK = 1024 -@cuda.jit -def _cuda_is_inside_circle(x, y, r, results): - i = cuda.grid(1) - if i > results.shape[0]: - return - else: - p = (math.sqrt((x[i][0] - y[0][0]) ** 2 + (x[i][1] - y[0][1]) ** 2)) - if p <= r[0]: - results[i] = 1 -def is_inside_circle(x: np.ndarray, y: np.ndarray, r: float) -> np.ndarray: - """ - Determines whether points in array `x` are inside the rectangle defined by the top left and bottom right vertices in array `y`. - - :param np.ndarray x: 2d numeric np.ndarray size (N, 2). - :param np.ndarray y: 2d numeric np.ndarray size (2, 2) (top left[x, y], bottom right[x, y]) - :return np.ndarray: 2d numeric boolean (N, 1) with 1s representing the point being inside the rectangle and 0 if the point is outside the rectangle. - """ - - x = np.ascontiguousarray(x).astype(np.int32) - y = np.ascontiguousarray(y).astype(np.int32) - x_dev = cuda.to_device(x) - y_dev = cuda.to_device(y) - r = np.array([r]).astype(np.float32) - r_dev = cuda.to_device(r) - results = cuda.device_array((x.shape[0]), dtype=np.int8) - bpg = (x.shape[0] + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK - del x, y - _cuda_is_inside_circle[bpg, THREADS_PER_BLOCK](x_dev, y_dev, r_dev, results) - results = results.copy_to_host() - return results diff --git a/simba/data_processors/cuda/is_inside_polygon.py b/simba/data_processors/cuda/is_inside_polygon.py deleted file mode 100644 index b106a05fa..000000000 --- a/simba/data_processors/cuda/is_inside_polygon.py +++ /dev/null @@ -1,65 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -import numpy as np -from numba import cuda - -THREADS_PER_BLOCK = 1024 - -@cuda.jit -def _cuda_is_inside_polygon(x, p, r): - i = cuda.grid(1) - if i > r.shape[0]: - return - else: - x, y, n = x[i][0], x[i][1], len(p) - p2x, p2y, xints, inside = 0.0, 0.0, 0.0, False - p1x, p1y = p[0] - for j in range(n + 1): - p2x, p2y = p[j % n] - if ( - (y > min(p1y, p2y)) - and (y <= max(p1y, p2y)) - and (x <= max(p1x, p2x)) - ): - if p1y != p2y: - xints = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x - if p1x == p2x or x <= xints: - inside = not inside - p1x, p1y = p2x, p2y - if inside: - r[i] = 1 - - -def is_inside_polygon(x: np.ndarray, y: np.ndarray) -> np.ndarray: - """ - Determines whether points in array `x` are inside the polygon defined by the vertices in array `y`. - - This function uses GPU acceleration to perform the point-in-polygon test. The points in `x` are tested against - the polygon defined by the vertices in `y`. The result is an array where each element indicates whether - the corresponding point is inside the polygon. - - .. image:: _static/img/is_inside_polygon_cuda.webp - :width: 500 - :align: center - - :param np.ndarray x: An array of shape (N, 2) where each row represents a point in 2D space. The points are checked against the polygon. - :param np.ndarray y: An array of shape (M, 2) where each row represents a vertex of the polygon in 2D space. - :return np.ndarray: An array of shape (N,) where each element is 1 if the corresponding point in `x` is inside the polygon defined by `y`, and 0 otherwise. - - :example: - >>> x = np.random.randint(0, 200, (i, 2)).astype(np.int8) - >>> y = np.random.randint(0, 200, (4, 2)).astype(np.int8) - >>> results = is_inside_polygon(x=x, y=y) - >>> print(results) - >>> [1 0 1 0 1 1 0 0 1 0] - """ - x = np.ascontiguousarray(x).astype(np.int32) - y = np.ascontiguousarray(y).astype(np.int32) - x_dev = cuda.to_device(x) - y_dev = cuda.to_device(y) - results = cuda.device_array((x.shape[0]), dtype=np.int8) - bpg = (x.shape[0] + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK - _cuda_is_inside_polygon[bpg, THREADS_PER_BLOCK](x_dev, y_dev, results) - results = results.copy_to_host() - return results \ No newline at end of file diff --git a/simba/data_processors/cuda/is_inside_rectangle.py b/simba/data_processors/cuda/is_inside_rectangle.py deleted file mode 100644 index 3920380a8..000000000 --- a/simba/data_processors/cuda/is_inside_rectangle.py +++ /dev/null @@ -1,45 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -import numpy as np -from numba import cuda - -THREADS_PER_BLOCK = 1024 -@cuda.jit -def _cuda_is_inside_rectangle(x, y, r): - i = cuda.grid(1) - if i > r.shape[0]: - return - else: - if (x[i][0] >= y[0][0]) and (x[i][0] <= y[1][0]): - if (x[i][1] >= y[0][1]) and (x[i][1] <= y[1][1]): - r[i] = 1 - -def is_inside_rectangle(x: np.ndarray, y: np.ndarray) -> np.ndarray: - """ - Determines whether points in array `x` are inside the rectangle defined by the top left and bottom right vertices in array `y`. - - .. csv-table:: - :header: EXPECTED RUNTIMES - :file: ../../../docs/tables/is_inside_rectangle.csv - :widths: 10, 45, 45 - :align: center - :class: simba-table - :header-rows: 1 - - :param np.ndarray x: 2d numeric np.ndarray size (N, 2). - :param np.ndarray y: 2d numeric np.ndarray size (2, 2) (top left[x, y], bottom right[x, y]) - :return np.ndarray: 2d numeric boolean (N, 1) with 1s representing the point being inside the rectangle and 0 if the point is outside the rectangle. - - """ - - - x = np.ascontiguousarray(x).astype(np.int32) - y = np.ascontiguousarray(y).astype(np.int32) - x_dev = cuda.to_device(x) - y_dev = cuda.to_device(y) - results = cuda.device_array((x.shape[0]), dtype=np.int8) - bpg = (x.shape[0] + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK - _cuda_is_inside_rectangle[bpg, THREADS_PER_BLOCK](x_dev, y_dev, results) - results = results.copy_to_host() - return results \ No newline at end of file diff --git a/simba/data_processors/cuda/sliding_circular_hotspots.py b/simba/data_processors/cuda/sliding_circular_hotspots.py deleted file mode 100644 index fafebc3ae..000000000 --- a/simba/data_processors/cuda/sliding_circular_hotspots.py +++ /dev/null @@ -1,55 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -from typing import Optional - -try: - from typing import Literal -except: - from typing_extensions import Literal - -import cupy as cp -import numpy as np - - -def sliding_circular_hotspots(x: np.ndarray, - time_window: float, - sample_rate: float, - bins: np.ndarray, - batch_size: Optional[int] = int(3.5e+7)) -> np.ndarray: - """ - Calculate the proportion of data points falling within specified circular bins over a sliding time window using GPU - - This function processes time series data representing angles (in degrees) and calculates the proportion of data - points within specified angular bins over a sliding window. The calculations are performed in batches to - accommodate large datasets efficiently. - - :param np.ndarray x: The input time series data in degrees. Should be a 1D numpy array. - :param float time_window: The size of the sliding window in seconds. - :param float sample_rate: The sample rate of the time series data (i.e., hz, fps). - :param ndarray bins: 2D array of shape representing circular bins defining [start_degree, end_degree] inclusive. - :param Optional[int] batch_size: The size of each batch for processing the data. Default is 5e+7 (50m). - :return: A 2D numpy array where each row corresponds to a time point in `data`, and each column represents a circular bin. The values in the array represent the proportion of data points within each bin at each time point. The first column represents the first bin. - :rtype: np.ndarray - """ - - n = x.shape[0] - x = cp.asarray(x, dtype=cp.float16) - results = cp.full((x.shape[0], bins.shape[0]), dtype=cp.float16, fill_value=-1) - window_size = int(cp.ceil(time_window * sample_rate)) - for cnt, left in enumerate(range(0, n, batch_size)): - right = int(min(left + batch_size, n)) - if cnt > 0: - left = left - window_size + 1 - x_batch = x[left:right] - x_batch = cp.lib.stride_tricks.sliding_window_view(x_batch, window_size).astype(cp.float16) - batch_results = cp.full((x_batch.shape[0], bins.shape[0]), dtype=cp.float16, fill_value=-1) - for bin_cnt in range(bins.shape[0]): - if bins[bin_cnt][0] > bins[bin_cnt][1]: - mask = ((x_batch >= bins[bin_cnt][0]) & (x_batch <= 360)) | ((x_batch >= 0) & (x_batch <= bins[bin_cnt][1])) - else: - mask = (x_batch >= bins[bin_cnt][0]) & (x_batch <= bins[bin_cnt][1]) - count_per_row = cp.array(mask.sum(axis=1) / window_size).reshape(-1, ) - batch_results[:, bin_cnt] = count_per_row - results[left + window_size - 1:right, ] = batch_results - return results.get() \ No newline at end of file diff --git a/simba/data_processors/cuda/sliding_circular_mean.py b/simba/data_processors/cuda/sliding_circular_mean.py deleted file mode 100644 index 0d670ae8a..000000000 --- a/simba/data_processors/cuda/sliding_circular_mean.py +++ /dev/null @@ -1,61 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -from typing import Optional - -import cupy -import numpy as np - - -def sliding_circular_mean(x: np.ndarray, - time_window: float, - sample_rate: int, - batch_size: Optional[int] = 3e+7) -> np.ndarray: - - """ - Calculate the sliding circular mean over a time window for a series of angles. - - This function computes the circular mean of angles in the input array `x` over a specified sliding window. - The circular mean is a measure of the average direction for angles, which is especially useful for angular data - where traditional averaging would not be meaningful due to the circular nature of angles (e.g., 359° and 1° should average to 0°). - - The calculation is performed using a sliding window approach, where the circular mean is computed for each window - of angles. The function leverages GPU acceleration via CuPy for efficiency when processing large datasets. - - The circular mean :math:`\\mu` for a set of angles is calculated using the following formula: - - .. math:: - - \\mu = \\text{atan2}\\left(\\frac{1}{N} \\sum_{i=1}^{N} \\sin(\\theta_i), \\frac{1}{N} \\sum_{i=1}^{N} \\cos(\\theta_i)\\right) - - - :math:`\\theta_i` are the angles in radians within the sliding window - - :math:`N` is the number of samples in the window - - - :param np.ndarray x: Input array containing angle values in degrees. The array should be 1-dimensional. - :param float time_window: Time duration for the sliding window, in seconds. This determines the number of samples in each window based on the `sample_rate`. - :param int sample_rate: The number of samples per second (i.e., FPS). This is used to calculate the window size in terms of array indices. - :param Optional[int] batch_size: The maximum number of elements to process in each batch. This is used to handle large arrays by processing them in chunks to avoid memory overflow. Defaults to 3e+7 (30 million elements). - :return np.ndarray: A 1D numpy array of the same length as `x`, containing the circular mean for each sliding window. Values before the window is fully populated will be set to -1. - - :example: - >>> x = np.random.randint(0, 361, (i, )).astype(np.int32) - >>> results = sliding_circular_mean(x, 1, 10) - """ - - - window_size = np.ceil(time_window * sample_rate).astype(np.int64) - n = x.shape[0] - results = cupy.full(x.shape[0], -1, dtype=np.int32) - for cnt, left in enumerate(range(0, int(n), int(batch_size))): - right = np.int32(min(left + batch_size, n)) - if cnt > 0: - left = left - window_size+1 - x_batch = cupy.asarray(x[left:right]) - x_batch = cupy.lib.stride_tricks.sliding_window_view(x_batch, window_size) - x_batch = np.deg2rad(x_batch) - cos, sin = cupy.cos(x_batch).astype(np.float32), cupy.sin(x_batch).astype(np.float32) - r = cupy.rad2deg(cupy.arctan2(cupy.mean(sin, axis=1), cupy.mean(cos, axis=1))) - r = cupy.where(r < 0, r + 360, r) - results[left + window_size - 1:right] = r - return results.get() diff --git a/simba/data_processors/cuda/sliding_circular_range.py b/simba/data_processors/cuda/sliding_circular_range.py deleted file mode 100644 index ecdd4aba7..000000000 --- a/simba/data_processors/cuda/sliding_circular_range.py +++ /dev/null @@ -1,60 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -from typing import Optional - -try: - from typing import Literal -except: - from typing_extensions import Literal - -import cupy as cp -import numpy as np - - -def sliding_circular_range(x: np.ndarray, - time_window: float, - sample_rate: float, - batch_size: Optional[int] = int(5e+7)) -> np.ndarray: - """ - Computes the sliding circular range of a time series data array using GPU. - - This function calculates the circular range of a time series data array using a sliding window approach. - The input data is assumed to be in degrees, and the function handles the circular nature of the data - by considering the circular distance between angles. - - .. math:: - - R = \\min \\left( \\text{max}(\\Delta \\theta) - \\text{min}(\\Delta \\theta), \\, 360 - \\text{max}(\\Delta \\theta) + \\text{min}(\\Delta \\theta) \\right) - - where: - - - :math:`\\Delta \\theta` is the difference between angles within the window, - - :math:`360` accounts for the circular nature of the data (i.e., wrap-around at 360 degrees). - - :param np.ndarray x: The input time series data in degrees. Should be a 1D numpy array. - :param float time_window: The size of the sliding window in seconds. - :param float sample_rate: The sample rate of the time series data (i.e., hz, fps). - :param Optional[int] batch_size: The size of each batch for processing the data. Default is 5e+7 (50m). - :return: A numpy array containing the sliding circular range values. - :rtype: np.ndarray - - :example: - >>> x = np.random.randint(0, 361, (19, )).astype(np.int32) - >>> p = sliding_circular_range(x, 1, 10) - """ - - n = x.shape[0] - x = cp.asarray(x, dtype=cp.float16) - results = cp.zeros_like(x, dtype=cp.int16) - x = cp.deg2rad(x).astype(cp.float16) - window_size = int(cp.ceil(time_window * sample_rate)) - for cnt, left in enumerate(range(0, n, batch_size)): - right = int(min(left + batch_size, n)) - if cnt > 0: - left = left - window_size + 1 - x_batch = x[left:right] - x_batch = cp.lib.stride_tricks.sliding_window_view(x_batch, window_size).astype(cp.float16) - x_batch = cp.sort(x_batch) - results[left + window_size - 1:right] = cp.abs(cp.rint(cp.rad2deg(cp.amin(cp.vstack([x_batch[:, -1] - x_batch[:, 0], 2 * cp.pi - cp.max(cp.diff(x_batch), axis=1)]).T, axis=1)))) - return results.get() \ No newline at end of file diff --git a/simba/data_processors/cuda/sliding_circular_std.py b/simba/data_processors/cuda/sliding_circular_std.py deleted file mode 100644 index e32cb1fe3..000000000 --- a/simba/data_processors/cuda/sliding_circular_std.py +++ /dev/null @@ -1,60 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -from typing import Optional - -try: - from typing import Literal -except: - from typing_extensions import Literal - -import cupy as cp -import numpy as np - - -def sliding_circular_std(x: np.ndarray, - time_window: float, - sample_rate: float, - batch_size: Optional[int] = int(5e+7)) -> np.ndarray: - """ - Calculate the sliding circular standard deviation of a time series data on GPU. - - This function computes the circular standard deviation over a sliding window for a given time series array. - The time series data is assumed to be in degrees, and the function converts it to radians for computation. - The sliding window approach is used to handle large datasets efficiently, processing the data in batches. - - The circular standard deviation (σ) is computed using the formula: - - .. math:: - - \sigma = \sqrt{-2 \cdot \log \left|\text{mean}\left(\exp(i \cdot x_{\text{batch}})\right)\right|} - - where :math:`x_{\text{batch}}` is the data within the current sliding window, and :math:`\text{mean}` and - :math:`\log` are computed in the circular (complex plane) domain. - - :param np.ndarray x: The input time series data in degrees. Should be a 1D numpy array. - :param float time_window: The size of the sliding window in seconds. - :param float sample_rate: The sample rate of the time series data (i.e., hz, fps). - :param Optional[int] batch_size: The size of each batch for processing the data. Default is 5e+7 (50m). - - :return: A numpy array containing the sliding circular standard deviation values. - :rtype: np.ndarray - """ - - - n = x.shape[0] - x = cp.asarray(x, dtype=cp.float16) - results = cp.zeros_like(x, dtype=cp.float16) - x = np.deg2rad(x).astype(cp.float16) - window_size = int(np.ceil(time_window * sample_rate)) - for cnt, left in enumerate(range(0, n, batch_size)): - right = int(min(left + batch_size, n)) - if cnt > 0: - left = left - window_size + 1 - x_batch = x[left:right] - x_batch = cp.lib.stride_tricks.sliding_window_view(x_batch, window_size).astype(cp.float16) - m = cp.log(cp.abs(cp.mean(cp.exp(1j * x_batch), axis=1))) - stdev = cp.rad2deg(cp.sqrt(-2 * m)) - results[left + window_size - 1:right] = stdev - - return results.get() diff --git a/simba/data_processors/cuda/sliding_mean.py b/simba/data_processors/cuda/sliding_mean.py deleted file mode 100644 index ff00fcd5e..000000000 --- a/simba/data_processors/cuda/sliding_mean.py +++ /dev/null @@ -1,54 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -import numpy as np -from numba import cuda - -THREADS_PER_BLOCK = 1024 - -@cuda.jit(device=True) -def _cuda_sum(x: np.ndarray): - s = 0 - for i in range(x.shape[0]): - s += x[i] - return s - -@cuda.jit -def _cuda_sliding_mean(x: np.ndarray, d: np.ndarray, results: np.ndarray): - r = cuda.grid(1) - l = np.int32(r - (d[0] - 1)) - if (r >= results.shape[0]) or (l < 0): - results[r] = -1 - else: - x_i = x[l:r+1] - s = _cuda_sum(x_i) - results[r] = s / x_i.shape[0] - -def sliding_mean(x: np.ndarray, time_window: float, sample_rate: int) -> np.ndarray: - """ - Computes the mean of values within a sliding window over a 1D numpy array `x` using CUDA for acceleration. - - .. image:: _static/img/sliding_mean_cuda.png - :width: 500 - :align: center - - :param np.ndarray x: The input 1D numpy array of floats. The array over which the sliding window sum is computed. - :param float time_window:The size of the sliding window in seconds. This window slides over the array `x` to compute the sum. - :param int sample_rate: The number of samples per second in the array `x`. This is used to convert the time-based window size into the number of samples. - :return np.ndarray: A numpy array containing the sum of values within each position of the sliding window. - - :example: - >>> x = np.random.randint(1, 11, (100, )).astype(np.float32) - >>> time_window = 1 - >>> sample_rate = 10 - >>> r_x = sliding_mean(x=x, time_window=time_window, sample_rate=10) - """ - x = np.ascontiguousarray(x).astype(np.int32) - window_size = np.array([np.ceil(time_window * sample_rate)]) - x_dev = cuda.to_device(x) - delta_dev = cuda.to_device(window_size) - results = cuda.device_array(x.shape, dtype=np.float32) - bpg = (x.shape[0] + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK - _cuda_sliding_mean[bpg, THREADS_PER_BLOCK](x_dev, delta_dev, results) - results = results.copy_to_host() - return results \ No newline at end of file diff --git a/simba/data_processors/cuda/sliding_min.py b/simba/data_processors/cuda/sliding_min.py deleted file mode 100644 index a3db6769d..000000000 --- a/simba/data_processors/cuda/sliding_min.py +++ /dev/null @@ -1,52 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -import numpy as np -from numba import cuda - -THREADS_PER_BLOCK = 1024 - -@cuda.jit -def _cuda_sliding_min(x: np.ndarray, d: np.ndarray, results: np.ndarray): - def _cuda_min(a, b): - return a if a < b else b - r = cuda.grid(1) - l = np.int32(r - (d[0]-1)) - if (r > results.shape[0]) or (l < 0): - results[r] = -1 - else: - x_i = x[l:r-1] - local_min = x_i[0] - for k in range(x_i.shape[0]): - local_min = _cuda_min(local_min, x_i[k]) - results[r] = local_min - -def sliding_min(x: np.ndarray, time_window: float, sample_rate: int) -> np.ndarray: - """ - Computes the minimum value within a sliding window over a 1D numpy array `x` using CUDA for acceleration. - - .. image:: _static/img/sliding_min_cuda.png - :width: 500 - :align: center - - :param np.ndarray x: Input 1D numpy array of floats. The array over which the sliding window minimum is computed. - :param float time_window: The size of the sliding window in seconds. - :param intsample_rate: The sampling rate of the data, which determines the number of samples per second. - :return: A numpy array containing the minimum value for each position of the sliding window. - - :example: - >>> x = np.arange(0, 10000000) - >>> time_window = 1 - >>> sample_rate = 10 - >>> sliding_min(x=x, time_window=time_window, sample_rate=sample_rate) - """ - - x = np.ascontiguousarray(x).astype(np.float32) - window_size = np.array([np.ceil(time_window * sample_rate)]) - x_dev = cuda.to_device(x) - delta_dev = cuda.to_device(window_size) - results = cuda.device_array(x.shape, dtype=np.float32) - bpg = (x.shape[0] + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK - _cuda_sliding_min[bpg, THREADS_PER_BLOCK](x_dev, delta_dev, results) - results = results.copy_to_host() - return results diff --git a/simba/data_processors/cuda/sliding_rayleigh_z.py b/simba/data_processors/cuda/sliding_rayleigh_z.py deleted file mode 100644 index ff0a65fc0..000000000 --- a/simba/data_processors/cuda/sliding_rayleigh_z.py +++ /dev/null @@ -1,76 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -from typing import Optional, Tuple - -try: - from typing import Literal -except: - from typing_extensions import Literal - -import cupy as cp -import numpy as np - - -def sliding_rayleigh_z(x: np.ndarray, - time_window: float, - sample_rate: float, - batch_size: Optional[int] = int(5e+7)) -> Tuple[np.ndarray, np.ndarray]: - - """ - Computes the Rayleigh Z-statistic over a sliding window for a given time series of angles - - This function calculates the Rayleigh Z-statistic, which tests the null hypothesis that the population of angles - is uniformly distributed around the circle. The calculation is performed over a sliding window across the input - time series, and results are computed in batches for memory efficiency. - - Data is processed using GPU acceleration via CuPy, which allows for faster computation compared to a CPU-based approach. - - .. note:: - Adapted from ``pingouin.circular.circ_rayleigh`` and ``pycircstat.tests.rayleigh``. - - - **Rayleigh Z-statistic:** - - The Rayleigh Z-statistic is given by: - - .. math:: - - R = \frac{1}{n} \sqrt{\left(\sum_{i=1}^{n} \cos(\theta_i)\right)^2 + \left(\sum_{i=1}^{n} \sin(\theta_i)\right)^2} - - where: - - :math:`\theta_i` are the angles in the window. - - :math:`n` is the number of angles in the window. - - :param np.ndarray x: Input array of angles in degrees. Should be a 1D numpy array. - :param float time_window: The size of the sliding window in time units (e.g., seconds). - :param float sample_rate: The sampling rate of the input time series in samples per time unit (e.g., Hz, fps). - :param Optional[int] batch_size: The number of samples to process in each batch. Default is 5e7 (50m). Reducing this value may save memory at the cost of longer computation time. - :return: - A tuple containing two numpy arrays: - - **z_results**: Rayleigh Z-statistics for each position in the input array where the window was fully applied. - - **p_results**: Corresponding p-values for the Rayleigh Z-statistics. - :rtype: Tuple[np.ndarray, np.ndarray] - """ - - n = x.shape[0] - x = cp.asarray(x, dtype=cp.float16) - z_results = cp.zeros_like(x, dtype=cp.float16) - p_results = cp.zeros_like(x, dtype=cp.float16) - x = np.deg2rad(x).astype(cp.float16) - window_size = int(np.ceil(time_window * sample_rate)) - for cnt, left in enumerate(range(0, n, batch_size)): - right = int(min(left + batch_size, n)) - if cnt > 0: - left = left - window_size + 1 - x_batch = x[left:right] - x_batch = cp.lib.stride_tricks.sliding_window_view(x_batch, window_size).astype(cp.float16) - cos_sums = cp.nansum(cp.cos(x_batch), axis=1) ** 2 - sin_sums = cp.nansum(cp.sin(x_batch), axis=1) ** 2 - R = cp.sqrt(cos_sums + sin_sums) / window_size - Z = window_size * (R**2) - P = cp.exp(np.sqrt(1 + 4 * window_size + 4 * (window_size ** 2 - R ** 2)) - (1 + 2 * window_size)) - z_results[left + window_size - 1:right] = Z - p_results[left + window_size - 1:right] = P - - return z_results.get(), p_results.get() \ No newline at end of file diff --git a/simba/data_processors/cuda/sliding_resultant_vector_length.py b/simba/data_processors/cuda/sliding_resultant_vector_length.py deleted file mode 100644 index 4d85e21cf..000000000 --- a/simba/data_processors/cuda/sliding_resultant_vector_length.py +++ /dev/null @@ -1,63 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -from typing import Optional - -import cupy -import numpy as np - - -def sliding_resultant_vector_length(x: np.ndarray, - time_window: float, - sample_rate: int, - batch_size: Optional[int] = 3e+7) -> np.ndarray: - - """ - Calculate the sliding resultant vector length over a time window for a series of angles. - - This function computes the resultant vector length (R) for each window of angles in the input array `x`. - The resultant vector length is a measure of the concentration of angles, and it ranges from 0 to 1, where 1 - indicates all angles point in the same direction, and 0 indicates uniform distribution of angles. - - For a given sliding window of angles, the resultant vector length :math:`R` is calculated using the following formula: - - .. math:: - - R = \\frac{1}{N} \\sqrt{\\left(\\sum_{i=1}^{N} \\cos(\\theta_i)\\right)^2 + \\left(\\sum_{i=1}^{N} \\sin(\\theta_i)\\right)^2} - - where: - - - :math:`\\theta_i` are the angles in radians within the sliding window - - :math:`N` is the number of samples in the window - - The computation is performed in a sliding window manner over the entire array, utilizing GPU acceleration - with CuPy for efficiency, especially on large datasets. - - - :param np.ndarray x: Input array containing angle values in degrees. The array should be 1-dimensional. - :param float time_window: Time duration for the sliding window, in seconds. This determines the number of samples in each window based on the `sample_rate`. - :param int sample_rate: The number of samples per second (i.e., FPS). This is used to calculate the window size in terms of array indices. - :param Optional[int] batch_size: The maximum number of elements to process in each batch. This is used to handle large arrays by processing them in chunks to avoid memory overflow. Defaults to 3e+7 (30 million elements). - :return np.ndarray: A 1D numpy array of the same length as `x`, containing the resultant vector length for each sliding window. Values before the window is fully populated will be set to -1. - - - :example: - >>> x = np.random.randint(0, 361, (5000, )).astype(np.int32) - >>> results = sliding_resultant_vector_length(x, 1, 10) - """ - - window_size = np.ceil(time_window * sample_rate).astype(np.int64) - n = x.shape[0] - results = cupy.full(x.shape[0], -1, dtype=np.float32) - for cnt, left in enumerate(range(0, int(n), int(batch_size))): - right = np.int32(min(left + batch_size, n)) - if cnt > 0: - left = left - window_size+1 - x_batch = cupy.asarray(x[left:right]) - x_batch = cupy.lib.stride_tricks.sliding_window_view(x_batch, window_size) - x_batch = np.deg2rad(x_batch) - cos, sin = cupy.cos(x_batch).astype(np.float32), cupy.sin(x_batch).astype(np.float32) - cos_sum, sin_sum = cupy.sum(cos, axis=1), cupy.sum(sin, axis=1) - r = np.sqrt(cos_sum ** 2 + sin_sum ** 2) / window_size - results[left+window_size-1:right] = r - return results.get() \ No newline at end of file diff --git a/simba/data_processors/cuda/sliding_spearmans_rank.py b/simba/data_processors/cuda/sliding_spearmans_rank.py deleted file mode 100644 index 85dae6af6..000000000 --- a/simba/data_processors/cuda/sliding_spearmans_rank.py +++ /dev/null @@ -1,67 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -from typing import Optional - -import cupy as cp -import numpy as np - - -def sliding_spearmans_rank(x: np.ndarray, - y: np.ndarray, - time_window: float, - sample_rate: int, - batch_size: Optional[int] = int(1.6e+7)) -> np.ndarray: - """ - Computes the Spearman's rank correlation coefficient between two 1D arrays `x` and `y` - over sliding windows of size `time_window * sample_rate`. The computation is performed - in batches to optimize memory usage, leveraging GPU acceleration with CuPy. - - .. math:: - \rho = 1 - \frac{6 \sum d_i^2}{n_w(n_w^2 - 1)} - - .. math:: - The function uses CuPy to perform GPU-accelerated calculations. Ensure that your environment - supports GPU computation with CuPy installed. - - Where: - - \( \rho \) is the Spearman's rank correlation coefficient. - - \( d_i \) is the difference between the ranks of corresponding elements in the sliding window. - - \( n_w \) is the size of the sliding window. - - :param np.ndarray x: The first 1D array containing the values for Feature 1. - :param np.ndarray y: The second 1D array containing the values for Feature 2. - :param float time_window: The size of the sliding window in seconds. - :param int sample_rate: The sampling rate (samples per second) of the data. - :param Optional[int] batch_size: The size of each batch to process at a time for memory efficiency. Defaults to 1.6e7. - :return: A 1D numpy array containing the Spearman's rank correlation coefficient for each sliding window. - :rtype: np.ndarray - - :example: - >>> x = np.array([9, 10, 13, 22, 15, 18, 15, 19, 32, 11]) - >>> y = np.array([11, 12, 15, 19, 21, 26, 19, 20, 22, 19]) - >>> sliding_spearmans_rank(x, y, time_window=0.5, sample_rate=2) - """ - - - window_size = int(np.ceil(time_window * sample_rate)) - n = x.shape[0] - results = cp.full(n, -1, dtype=cp.float32) - - for cnt, left in enumerate(range(0, n, batch_size)): - right = int(min(left + batch_size, n)) - if cnt > 0: - left = left - window_size + 1 - x_batch = cp.asarray(x[left:right]) - y_batch = cp.asarray(y[left:right]) - x_batch = cp.lib.stride_tricks.sliding_window_view(x_batch, window_size) - y_batch = cp.lib.stride_tricks.sliding_window_view(y_batch, window_size) - rank_x = cp.argsort(cp.argsort(x_batch, axis=1), axis=1) - rank_y = cp.argsort(cp.argsort(y_batch, axis=1), axis=1) - d_squared = cp.sum((rank_x - rank_y) ** 2, axis=1) - n_w = window_size - s = 1 - (6 * d_squared) / (n_w * (n_w ** 2 - 1)) - - results[left + window_size - 1:right] = s - - return cp.asnumpy(results) \ No newline at end of file diff --git a/simba/data_processors/cuda/sliding_std.py b/simba/data_processors/cuda/sliding_std.py deleted file mode 100644 index d33ccb451..000000000 --- a/simba/data_processors/cuda/sliding_std.py +++ /dev/null @@ -1,58 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -import numpy as np -from numba import cuda - -THREADS_PER_BLOCK = 1024 - -@cuda.jit(device=True) -def _cuda_sum(x: np.ndarray): - s = 0 - for i in range(x.shape[0]): - s += x[i] - return s - -@cuda.jit(device=True) -def _cuda_std(x: np.ndarray, x_hat: float): - std = 0 - for i in range(x.shape[0]): - std += (x[0] - x_hat) ** 2 - return std - -@cuda.jit(device=False) -def _cuda_sliding_std(x: np.ndarray, d: np.ndarray, results: np.ndarray): - r = cuda.grid(1) - l = np.int32(r - (d[0] - 1)) - if (r >= results.shape[0]) or (l < 0): - results[r] = -1 - else: - x_i = x[l:r + 1] - s = _cuda_sum(x_i) - m = s / x_i.shape[0] - std = _cuda_std(x_i, m) - results[r] = std - -def sliding_std(x: np.ndarray, time_window: float, sample_rate: int) -> np.ndarray: - """ - - :param np.ndarray x: The input 1D numpy array of floats. The array over which the sliding window sum is computed. - :param float time_window:The size of the sliding window in seconds. This window slides over the array `x` to compute the sum. - :param int sample_rate: The number of samples per second in the array `x`. This is used to convert the time-based window size into the number of samples. - :return np.ndarray: A numpy array containing the sum of values within each position of the sliding window. - - :example: - >>> x = np.random.randint(1, 11, (100, )).astype(np.float32) - >>> time_window = 1 - >>> sample_rate = 10 - >>> r_x = sliding_sum(x=x, time_window=time_window, sample_rate=10) - """ - x = np.ascontiguousarray(x).astype(np.int32) - window_size = np.array([np.ceil(time_window * sample_rate)]) - x_dev = cuda.to_device(x) - delta_dev = cuda.to_device(window_size) - results = cuda.device_array(x.shape, dtype=np.float32) - bpg = (x.shape[0] + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK - _cuda_sliding_std[bpg, THREADS_PER_BLOCK](x_dev, delta_dev, results) - results = results.copy_to_host() - return results \ No newline at end of file diff --git a/simba/data_processors/cuda/sliding_sum.py b/simba/data_processors/cuda/sliding_sum.py deleted file mode 100644 index 89f15045e..000000000 --- a/simba/data_processors/cuda/sliding_sum.py +++ /dev/null @@ -1,47 +0,0 @@ -__author__ = "Simon Nilsson" -__email__ = "sronilsson@gmail.com" - -import numpy as np -from numba import cuda - -THREADS_PER_BLOCK = 1024 - -@cuda.jit -def _cuda_sliding_sum(x: np.ndarray, d: np.ndarray, results: np.ndarray): - r = cuda.grid(1) - l = np.int32(r - (d[0]-1)) - if (r > results.shape[0]) or (l < 0): - results[r] = -1 - else: - x_i = x[l:r] - local_sum = 0 - for k in range(x_i.shape[0]): - local_sum += x_i[k] - results[r-1] = local_sum - -def sliding_sum(x: np.ndarray, time_window: float, sample_rate: int) -> np.ndarray: - """ - Computes the sum of values within a sliding window over a 1D numpy array `x` using CUDA for acceleration. - - :param np.ndarray x: The input 1D numpy array of floats. The array over which the sliding window sum is computed. - :param float time_window:The size of the sliding window in seconds. This window slides over the array `x` to compute the sum. - :param int sample_rate: The number of samples per second in the array `x`. This is used to convert the time-based window size into the number of samples. - :return np.ndarray: A numpy array containing the sum of values within each position of the sliding window. - - :example: - >>> x = np.random.randint(1, 11, (100, )).astype(np.float32) - >>> time_window = 1 - >>> sample_rate = 10 - >>> r_x = sliding_sum(x=x, time_window=time_window, sample_rate=10) - """ - x = np.ascontiguousarray(x).astype(np.float32) - window_size = np.array([np.ceil(time_window * sample_rate)]) - x_dev = cuda.to_device(x) - delta_dev = cuda.to_device(window_size) - results = cuda.device_array(x.shape, dtype=np.float32) - bpg = (x.shape[0] + (THREADS_PER_BLOCK - 1)) // THREADS_PER_BLOCK - _cuda_sliding_sum[bpg, THREADS_PER_BLOCK](x_dev, delta_dev, results) - results = results.copy_to_host() - return results - - diff --git a/simba/mixins/geometry_mixin.py b/simba/mixins/geometry_mixin.py index 6056d9355..11994b233 100644 --- a/simba/mixins/geometry_mixin.py +++ b/simba/mixins/geometry_mixin.py @@ -3548,12 +3548,34 @@ def filter_low_p_bps_for_shapes(x: np.ndarray, p: np.ndarray, threshold: float): results[i][j] = new_val return results - - - - - - + @staticmethod + def get_shape_lengths_widths(shapes: Union[List[Polygon], Polygon]) -> Dict[str, Any]: + """ + Calculate the lengths and widths of the minimum bounding rectangles of polygons. + + :param Union[List[Polygon], Polygon] shapes: A single Polygon or a list of Polygons for which the MBR dimensions are calculated. Each polygon is assumed to be a valid shapely Polygon object. If a single Polygon is provided, it is internally converted to a list + :return: A dictionary containing: + - 'lengths': A list of the lengths of the MBR for each polygon. + - 'widths': A list of the widths of the MBR for each polygon. + - 'max_length': The maximum length found among all polygons. + - 'min_length': The minimum length found among all polygons. + - 'max_width': The maximum width found among all polygons. + - 'min_width': The minimum width found among all polygons. + :rtype: Dict[str, Any] + """ + widths, lengths, max_length, max_width, min_length, min_width = [], [], -np.inf, -np.inf, np.inf, np.inf + if isinstance(shapes, Polygon): + shapes = [shapes] + for shape in shapes: + shape_cords = list(zip(*shape.exterior.coords.xy)) + mbr_lengths = [LineString((shape_cords[i], shape_cords[i + 1])).length for i in range(len(shape_cords) - 1)] + width, length = min(mbr_lengths), max(mbr_lengths) + min_length, max_length = min(min_length, length), max(max_length, length) + min_width, max_width = min(min_width, width), max(max_width, width) + lengths.append(length); + widths.append(width) + return {'lengths': lengths, 'widths': widths, 'max_length': max_length, 'min_length': min_length, + 'min_width': min_width, 'max_width': max_width} # data = np.array([[[364, 308], [383, 323], [403, 335], [423, 351]], # [[356, 307], [376, 319], [396, 331], [419, 347]], diff --git a/simba/mixins/image_mixin.py b/simba/mixins/image_mixin.py index 3cbd9755b..295f3a42e 100644 --- a/simba/mixins/image_mixin.py +++ b/simba/mixins/image_mixin.py @@ -1689,6 +1689,25 @@ def get_blob_locations(video_path: Union[str, os.PathLike], stdout_success(f'Video {video_meta["video_name"]} blob detection complete', elapsed_time=timer.elapsed_time_str) return results.astype(np.int32) + @staticmethod + def is_video_color(video: Union[str, os.PathLike, cv2.VideoCapture]): + """ + Determines whether a video is in color or greyscale. + + :param Union[str, os.PathLike, cv2.VideoCapture] video: The video source, either a cv2.VideoCapture object or a path to a file on disk. + :return: Returns `True` if the video is in color (has more than one channel), and `False` if the video is greyscale (single channel). + :rtype: bool + """ + + frm = ImageMixin.find_first_non_uniform_clr_frm(video_path=video) + if frm.ndim > 2: + if frm.shape[2] != 1: + return True + else: + return False + else: + return False + #x = ImageMixin.get_blob_locations(video_path=r"C:\troubleshooting\RAT_NOR\project_folder\videos\2022-06-20_NOB_DOT_4_downsampled_bg_subtracted.mp4", gpu=True) # imgs = ImageMixin().read_all_img_in_dir(dir='/Users/simon/Desktop/envs/simba/troubleshooting/RAT_NOR/project_folder/videos/examples') diff --git a/simba/roi_tools/ROI_define.py b/simba/roi_tools/ROI_define.py index 18468f0f8..886021b96 100644 --- a/simba/roi_tools/ROI_define.py +++ b/simba/roi_tools/ROI_define.py @@ -11,9 +11,8 @@ from simba.roi_tools.ROI_image import ROI_image_class from simba.roi_tools.ROI_move_shape import move_edge, update_all_tags from simba.roi_tools.ROI_multiply import create_emty_df -from simba.roi_tools.ROI_size_calculations import (circle_size_calc, - polygon_size_calc, - rectangle_size_calc) +from simba.roi_tools.ROI_size_calculations import (circle_size_calc, polygon_size_calc, rectangle_size_calc) +from simba.mixins.plotting_mixin import PlottingMixin from simba.ui.tkinter_functions import SimbaButton, hxtScrollbar from simba.utils.checks import check_file_exist_and_readable from simba.utils.enums import Formats, TagNames @@ -83,8 +82,8 @@ def __init__(self, config_path: str, video_path: str): self.img_no = 1 self.duplicate_jump_size = 20 self.click_sens = 10 - self.text_size = 5 - self.text_thickness = 3 + self.text_size, _, _ = PlottingMixin().get_optimal_font_scales(text='TEN DIGITS', accepted_px_width=int(self.video_info['Resolution_width']/10), accepted_px_height=int(self.video_info['Resolution_height']/10), text_thickness=2, font=cv2.FONT_HERSHEY_SIMPLEX) + self.text_thickness = 2 self.line_type = -1 self.named_shape_colors = get_color_dict() self.window_menus() @@ -318,7 +317,7 @@ def interact_menus(self): self.zoom_pct.insert(0, 10) self.pan = Button(self.interact_frame, text="Pan", fg=self.non_select_color, state=DISABLED, command=lambda: self.set_interact_state("pan")) - self.shape_info_btn = SimbaButton(parent=self.interact_frame, txt="SHOW SHAPE INFO", img='info', txt_clr=self.non_select_color, enabled=False, cmd=self.show_shape_information) + self.shape_info_btn = SimbaButton(parent=self.interact_frame, txt="SHOW SHAPE INFO", img='info', txt_clr=self.non_select_color, enabled=True, cmd=self.show_shape_information) self.interact_frame.grid(row=6, sticky=W) @@ -377,27 +376,21 @@ def show_shape_information(self): self.rectangle_size_dict = {} self.rectangle_size_dict["Rectangles"] = {} for rectangle in self.image_data.out_rectangles: - self.rectangle_size_dict["Rectangles"][rectangle["Name"]] = ( - rectangle_size_calc(rectangle, self.curr_px_mm) - ) + self.rectangle_size_dict["Rectangles"][rectangle["Name"]] = (rectangle_size_calc(rectangle, self.curr_px_mm)) self.image_data.rectangle_size_dict = self.rectangle_size_dict if len(self.image_data.out_circles) > 0: self.circle_size_dict = {} self.circle_size_dict["Circles"] = {} for circle in self.image_data.out_circles: - self.circle_size_dict["Circles"][circle["Name"]] = circle_size_calc( - circle, self.curr_px_mm - ) + self.circle_size_dict["Circles"][circle["Name"]] = circle_size_calc(circle, self.curr_px_mm) self.image_data.circle_size_dict = self.circle_size_dict if len(self.image_data.out_polygon) > 0: self.polygon_size_dict = {} self.polygon_size_dict["Polygons"] = {} for polygon in self.image_data.out_polygon: - self.polygon_size_dict["Polygons"][polygon["Name"]] = ( - polygon_size_calc(polygon, self.curr_px_mm) - ) + self.polygon_size_dict["Polygons"][polygon["Name"]] = (polygon_size_calc(polygon, self.curr_px_mm)) self.image_data.polygon_size_dict = self.polygon_size_dict self.image_data.insert_all_ROIs_into_image(show_size_info=True) diff --git a/simba/roi_tools/ROI_image.py b/simba/roi_tools/ROI_image.py index 3a17c0a38..9cbe3ea2f 100644 --- a/simba/roi_tools/ROI_image.py +++ b/simba/roi_tools/ROI_image.py @@ -1,7 +1,6 @@ import itertools import os import re -import threading from copy import deepcopy import cv2 @@ -14,6 +13,7 @@ from simba.utils.enums import ConfigKey, Keys, Paths from simba.utils.errors import InvalidInputError from simba.utils.read_write import get_fn_ext, read_config_file +from simba.mixins.plotting_mixin import PlottingMixin class ROI_image_class: @@ -601,78 +601,20 @@ def insert_all_ROIs_into_image( if show_size_info is True: area_cm = self.rectangle_size_dict["Rectangles"][e["Name"]]["area_cm"] width_cm = self.rectangle_size_dict["Rectangles"][e["Name"]]["width_cm"] - height_cm = self.rectangle_size_dict["Rectangles"][e["Name"]][ - "height_cm" - ] - cv2.putText( - self.working_frame, - str(height_cm), - ( - int(e["Tags"]["Left tag"][0] + e["Thickness"]), - e["Tags"]["Left tag"][1], - ), - cv2.FONT_HERSHEY_SIMPLEX, - self.text_size / 10, - self.colors[e["Color name"]], - self.text_thickness, - lineType=self.line_type, - ) - cv2.putText( - self.working_frame, - str(width_cm), - ( - e["Tags"]["Bottom tag"][0], - int(e["Tags"]["Bottom tag"][1] - e["Thickness"]), - ), - cv2.FONT_HERSHEY_SIMPLEX, - self.text_size / 10, - self.colors[e["Color name"]], - self.text_thickness, - lineType=self.line_type, - ) - cv2.putText( - self.working_frame, - str(area_cm), - (e["Tags"]["Center tag"][0], e["Tags"]["Center tag"][1]), - cv2.FONT_HERSHEY_SIMPLEX, - self.text_size / 10, - self.colors[e["Color name"]], - self.text_thickness, - lineType=self.line_type, - ) + height_cm = self.rectangle_size_dict["Rectangles"][e["Name"]]["height_cm"] + self.working_frame = PlottingMixin().put_text(img=self.working_frame, text=f'AREA: {str(area_cm)}, HEIGHT: {str(height_cm)}, WIDTH: {str(width_cm)}', pos=(int(e["Tags"]["Left tag"][0] + e["Thickness"]), e["Tags"]["Left tag"][1]), font_size=self.text_size, font_thickness=self.text_thickness, font=cv2.FONT_HERSHEY_SIMPLEX, text_color=self.colors[e["Color name"]]) + print(f'{e["Name"]} area: {area_cm}') for c in self.out_circles: self.no_shapes += 1 - cv2.circle( - self.working_frame, - (c["centerX"], c["centerY"]), - c["radius"], - c["Color BGR"], - int(c["Thickness"]), - lineType=self.line_type, - ) + cv2.circle( self.working_frame, (c["centerX"], c["centerY"]), c["radius"], c["Color BGR"], int(c["Thickness"]), lineType=self.line_type) if ROI_ear_tags is True: for t in c["Tags"]: - cv2.circle( - self.working_frame, - c["Tags"][t], - c["Ear_tag_size"], - self.colors[c["Color name"]], - -1, - ) + cv2.circle( self.working_frame, c["Tags"][t], c["Ear_tag_size"], self.colors[c["Color name"]], -1) if show_size_info is True: area_cm = self.circle_size_dict["Circles"][c["Name"]]["area_cm"] - cv2.putText( - self.working_frame, - str(area_cm), - (c["Tags"]["Center tag"][0], c["Tags"]["Center tag"][1]), - cv2.FONT_HERSHEY_SIMPLEX, - self.text_size / 10, - self.colors[c["Color name"]], - self.text_thickness, - lineType=self.line_type, - ) - + self.working_frame = PlottingMixin().put_text(img=self.working_frame, text=f'AREA: {str(area_cm)}', pos=(c["Tags"]["Center tag"][0], c["Tags"]["Center tag"][1]), font_size=self.text_size, font_thickness=self.text_thickness, font=cv2.FONT_HERSHEY_SIMPLEX, text_color=self.colors[c["Color name"]]) + print(f'{c["Name"]} area: {area_cm}') for pg in self.out_polygon: self.no_shapes += 1 pts = np.array(pg["vertices"]).reshape((-1, 1, 2)) @@ -715,17 +657,8 @@ def insert_all_ROIs_into_image( ) if show_size_info is True: area_cm = self.polygon_size_dict["Polygons"][pg["Name"]]["area_cm"] - cv2.putText( - self.working_frame, - str(area_cm), - (pg["Center_X"], pg["Center_Y"]), - cv2.FONT_HERSHEY_SIMPLEX, - self.text_size / 10, - self.colors[pg["Color name"]], - self.text_thickness, - lineType=self.line_type, - ) - + self.working_frame = PlottingMixin().put_text(img=self.working_frame, text=f'AREA: {str(area_cm)}', pos=(pg["Center_X"], pg["Center_Y"]), font_size=self.text_size, font_thickness=self.text_thickness, font=cv2.FONT_HERSHEY_SIMPLEX, text_color=self.colors[pg["Color name"]]) + print(f'{pg["Name"]} area: {area_cm}') cv2.namedWindow("Define shape", cv2.WINDOW_NORMAL) cv2.imshow("Define shape", self.working_frame) cv2.waitKey(100) diff --git a/simba/roi_tools/ROI_size_calculations.py b/simba/roi_tools/ROI_size_calculations.py index 4993f5ace..e4f98447b 100644 --- a/simba/roi_tools/ROI_size_calculations.py +++ b/simba/roi_tools/ROI_size_calculations.py @@ -1,6 +1,7 @@ import math - import numpy as np +from shapely.geometry import Polygon +from scipy.spatial import ConvexHull def rectangle_size_calc(rectangle_dict: dict, px_mm: float) -> dict: @@ -18,9 +19,7 @@ def rectangle_size_calc(rectangle_dict: dict, px_mm: float) -> dict: rectangle_dict["height_cm"] = round((rectangle_dict["height"] / px_mm) / 10, 2) rectangle_dict["width_cm"] = round((rectangle_dict["width"] / px_mm) / 10, 2) - rectangle_dict["area_cm"] = round( - rectangle_dict["width_cm"] * rectangle_dict["height_cm"], 2 - ) + rectangle_dict["area_cm"] = round(rectangle_dict["width_cm"] * rectangle_dict["height_cm"], 2) return rectangle_dict @@ -53,18 +52,14 @@ def polygon_size_calc(polygon_dict, px_mm) -> dict: >>> polygon_size_calc(polygon_dict={'vertices': np.array([[0, 2], [200, 98], [100, 876], [10, 702]])}, px_mm=5) >>> {'vertices': [[ 0, 2], [200, 98], [100, 876], [ 10, 702]], 'area_cm': 45.29} """ - - y_vals = polygon_dict["vertices"][:, 0] - x_vals = polygon_dict["vertices"][:, 1] - poly_area_px = 0.5 * np.abs( - np.dot(x_vals, np.roll(y_vals, 1)) - np.dot(y_vals, np.roll(x_vals, 1)) - ) - polygon_dict["area_cm"] = round((poly_area_px / px_mm) / 500, 2) + polygon = polygon_dict["vertices"] + area = round((ConvexHull(polygon).area / px_mm) / 10, 2) + polygon_dict["area_cm"] = area return polygon_dict -polygon_size_calc( - polygon_dict={"vertices": np.array([[0, 2], [200, 98], [100, 876], [10, 702]])}, - px_mm=5, -) +# polygon_size_calc( +# polygon_dict={"vertices": np.array([[0, 2], [200, 98], [100, 876], [10, 702]])}, +# px_mm=5, +# ) diff --git a/simba/utils/read_write.py b/simba/utils/read_write.py index 317912b11..045e6cd33 100644 --- a/simba/utils/read_write.py +++ b/simba/utils/read_write.py @@ -39,7 +39,8 @@ check_if_string_value_is_valid_video_timestamp, check_instance, check_int, check_nvidea_gpu_available, check_str, - check_valid_dataframe, check_valid_lst) + check_valid_dataframe, check_valid_lst, + check_valid_array, check_valid_boolean) from simba.utils.enums import ConfigKey, Dtypes, Formats, Keys, Options from simba.utils.errors import (DataHeaderError, DuplicationError, FFMPEGCodecGPUError, FileExistError, @@ -2336,3 +2337,69 @@ def read_boris_file(file_path: Union[str, os.PathLike], return results else: write_pickle(data=results, save_path=save_path) + + +def img_stack_to_video(x: np.ndarray, + save_path: Union[str, os.PathLike], + fps: float, + gpu: Optional[bool] = False, + bitrate: Optional[int] = 5000) -> None: + + """ + Converts a NumPy image stack to a video file, with optional GPU acceleration and configurable bitrate. + + :param np.ndarray x: A NumPy array representing the image stack. The array should have shape (N, H, W) for greyscale or (N, H, W, 3) for RGB images, where N is the number of frames, H is the height, and W is the width. + :param Union[str, os.PathLike] save_path: Path to the output video file where the video will be saved. + :param float fps: Frames per second for the output video. Should be a positive floating-point number. + :param Optional[bool] gpu: Whether to use GPU acceleration for encoding. If True, the video encoding will use NVIDIA's NVENC encoder. Defaults to False. + :param Optional[int] bitrate: Bitrate for the video encoding in kilobits per second (kbps). Should be an integer between 1000 and 35000. Defaults to 5000. + :return: None + """ + + check_if_dir_exists(in_dir=os.path.dirname(save_path), source=img_stack_to_video.__name__) + check_valid_array(data=x, source=img_stack_to_video.__name__, accepted_ndims=(3, 4)) + check_float(name=f'{img_stack_to_video.__name__} fps', value=fps, min_value=10e-6) + check_valid_boolean(value=gpu, source=img_stack_to_video.__name__) + check_int(name=f'{img_stack_to_video.__name__} bitrate', value=bitrate, min_value=1000, max_value=35000) + if gpu and not check_nvidea_gpu_available(): + raise FFMPEGCodecGPUError('No GPU found but GPU flag is True') + is_color = (x.ndim == 4 and x.shape[3] == 3) + N, H, W = x.shape[:3] + pix_fmt = 'gray' if not is_color else 'rgb24' + timer = SimbaTimer(start=True) + vcodec = 'mpeg4' + if gpu: + vcodec = 'h264_nvenc' + + cmd = [ + 'ffmpeg', + '-loglevel', 'error', + '-stats', + '-hide_banner', + '-y', + '-f', 'rawvideo', + '-vcodec', 'rawvideo', + '-s', f'{W}x{H}', + '-pix_fmt', pix_fmt, + '-r', str(fps), + '-i', '-', + '-an', + '-vcodec', f'{vcodec}', + '-b:v', f'{bitrate}k', + save_path + ] + + process = subprocess.Popen(cmd, stdin=subprocess.PIPE) + for frame in x: + if frame.dtype != np.uint8: + frame = (255 * np.clip(frame, 0, 1)).astype(np.uint8) + process.stdin.write(frame.tobytes()) + + process.stdin.close() + process.wait() + timer.stop_timer() + stdout_success(msg=f'Video complete. Saved at {save_path}', elapsed_time=timer.elapsed_time_str) + + + +