magtrack.core ============= .. py:module:: magtrack.core Functions --------- .. autoapisummary:: magtrack.core.binmean magtrack.core.pearson magtrack.core.gaussian magtrack.core.gaussian_2d magtrack.core.crop_stack_to_rois magtrack.core.parabolic_vertex magtrack.core._qi_sample_axis_profiles magtrack.core._qi_quadratic_offsets magtrack.core.qi magtrack.core.center_of_mass magtrack.core.auto_conv magtrack.core.auto_conv_sub_pixel magtrack.core.auto_conv_multiline magtrack.core.auto_conv_multiline_sub_pixel magtrack.core.radial_profile magtrack.core.fft_profile magtrack.core.lookup_z magtrack.core.stack_to_xyzp_advanced magtrack.core.stack_to_xyzp Module Contents --------------- .. py:function:: binmean(x, weights, n_bins: int) Compute mean values per bin for 2D arrays, similar to ``numpy.bincount``. Note: CPU or GPU: The code is agnostic of CPU and GPU usage. If the first parameter is on the GPU the computation/result will be on the GPU. Otherwise, the computation/result will be on the CPU. The input ``x`` is clipped *in place* so that values above ``n_bins`` fall back within the valid bin range; entries clipped to ``n_bins`` are ignored when returning the binned means. :param x: Input array to bin. :type x: 2D int array, shape (n_values, n_datasets) :param weights: Weights associated with ``x``; should be floating point to allow averaging. :type weights: 2D float array, shape (n_values, n_datasets) :param n_bins: The number of bins to be used. Values will be binned as integers between 0 and n_bins. :type n_bins: int :returns: **bin_means** -- Binned average values of weights. :rtype: 2D float array, shape (n_bins, n_datasets) .. py:function:: pearson(x, y) Calculate the Pearson correlation coefficient between each row of x and y. Note: CPU or GPU: The code is agnostic of CPU and GPU usage. If the first parameter is on the GPU the computation/result will be on the GPU. Otherwise, the computation/result will be on the CPU. :param x: 2D array whose columns are correlated with the columns of ``y``. :type x: array, shape (n, m) :param y: 2D array whose columns are correlated with the columns of ``x``. :type y: array, shape (n, k) :returns: **r** -- Pearson correlation coefficients between each column of ``y`` and each column of ``x``. :rtype: array, shape (k, m) .. py:function:: gaussian(x, mu, sigma) Calculate a 1D Gaussian function. :param x: x coordinates where to evaluate the gaussian :type x: array_like :param mu: Mean (center) :type mu: float :param sigma: Standard deviation :type sigma: float :returns: 1D array containing the gaussian evaluated at x coordinates :rtype: array_like .. py:function:: gaussian_2d(x, y, mu_x, mu_y, sigma) Calculate a 2D Gaussian image. Calculates a 2D Gaussian image for each center (mu_x, mu_y) provided along the grid (x, y) all sharing the same sigma in x and y (sigma). Note: CPU or GPU: The code is agnostic of CPU and GPU usage. If the first parameter is on the GPU the computation/result will be on the GPU. Otherwise, the computation/result will be on the CPU. :param x: x coordinates where to evaluate the gaussian :type x: 1D array :param y: y coordinates where to evaluate the gaussian :type y: 1D array :param mu_x: Mean (center) in x direction for each image (one center per image) :type mu_x: 1D array :param mu_y: Mean (center) in y direction for each image (one center per image) :type mu_y: 1D array :param sigma: Standard deviation in x and y direction :type sigma: float :returns: 3D array of shape ``(len(x), len(y), n_images)`` containing the gaussian evaluated at ``(x, y)`` coordinates for each image center :rtype: array .. py:function:: crop_stack_to_rois(stack, rois) Takes a 3D image-stack and crops it to the region of interests (ROIs). Given a 3D image-stack and a list of ROIs, this function will crop around each ROI and return a 4D array. Note the ROIs must be squares. Note: CPU or GPU: The code is agnostic of CPU and GPU usage. If the first parameter is on the GPU the computation/result will be on the GPU. Otherwise, the computation/result will be on the CPU. However, it is recommended to use the CPU and then transfer the result to the GPU and perform downstream analysis on the GPU. :param stack: Note the images must be square. :type stack: 3D ndarray of any type, shape (stack_width, stack_height, n_images) :param rois: Each row is an ROI. The columns are [top, bottom, left, right]. :type rois: 2D int ndarray, shape (n_roi, 4) :returns: **cropped_stack** -- Same type as input stack :rtype: 4D ndarray, shape (cropped_width, cropped_width, n_images, n_roi) .. py:function:: parabolic_vertex(data, vertex_est, n_local: int, weighted=True) Refine local min/max using parabolic interpolation. Given an estimated location of a local minimum or maximum, this function fits the surrounding datapoints to a parabola and interpolates the vertex. Note: CPU or GPU: The code is agnostic of CPU and GPU usage. If the first parameter is on the GPU the computation/result will be on the GPU. Otherwise, the computation/result will be on the CPU. :param data: Sequence of datasets arranged row-wise for fitting. :type data: array of float, shape (n_datasets, n_datapoints) :param vertex_est: Initial vertex estimates corresponding to each dataset. :type vertex_est: array of float, shape (n_datasets,) :param n_local: The number of local datapoints to be fit. Must be an odd integer >= 3. :type n_local: int :param weighted: Whether to apply a simple weighting procedure to emphasize the more central points in the fit. Default is True. :type weighted: bool, optional :returns: **vertex** -- Refined vertex locations :rtype: array of float, shape (n_datasets,) .. py:function:: _qi_sample_axis_profiles(stack, x, y, axis) Sample 1D profiles along ``axis`` using quadratic interpolation support. :param stack: Image stack where the first two axes correspond to ``y`` and ``x`` and the third axis indexes frames. :type stack: 3D array_like :param x: Approximate center coordinates for each frame. They must have the same length as the number of frames in ``stack``. :type x: 1D array_like :param y: Approximate center coordinates for each frame. They must have the same length as the number of frames in ``stack``. :type y: 1D array_like :param axis: Axis along which to collect the three-point profile. ``0`` samples the column profile (varying ``y``); ``1`` samples the row profile (varying ``x``). :type axis: int :returns: Samples of shape ``(n_frames, 3)`` corresponding to offsets of ``-1``, ``0`` and ``+1`` pixels along the chosen axis. :rtype: array_like .. py:function:: _qi_quadratic_offsets(samples) Compute quadratic-interpolation offsets from three-point samples. .. py:function:: qi(stack, x_old, y_old) Refine centers using quadratic interpolation along x and y axes. This routine samples the intensity profiles along the horizontal and vertical axes through the supplied center estimates and performs quadratic interpolation to recover sub-pixel offsets. Note: CPU or GPU: The code is agnostic of CPU and GPU usage. If the first parameter is on the GPU the computation/result will be on the GPU. Otherwise, the computation/result will be on the CPU. Intermediate values remain on the originating device and the function respects the caller's backend. :param stack: Image stack containing square frames to refine. :type stack: 3D float array, shape (n_pixels, n_pixels, n_images) :param x_old: Initial estimates of the x coordinates. :type x_old: 1D float array, shape (n_images) :param y_old: Initial estimates of the y coordinates. :type y_old: 1D float array, shape (n_images) :returns: Refined ``(x, y)`` coordinates with sub-pixel precision. :rtype: tuple of array_like .. py:function:: center_of_mass(stack, background='none') Calculate x and y by center-of-mass For each 2D image of a 3D image-stack compute the center-of-mass along the x- and y-axes. To avoid bias from the images' background, a pre-processing step can be taken to remove the background with the ``background`` keyword argument. The default, ``background='none'`` leaves the data unchanged. When ``background='mean'`` the per-frame mean is subtracted before calculating the center-of-mass. Likewise, ``background='median'`` subtracts the per-frame median. This function is faster than the version from ``scipy`` or ``cupyx.scipy``. Note: CPU or GPU: The code is agnostic of CPU and GPU usage. If the first parameter is on the GPU the computation/result will be on the GPU. Otherwise, the computation/result will be on the CPU. :param stack: The image-stack. The images must be square. :type stack: 3D float array, shape (n_pixels, n_pixels, n_images) :param background: Background pre-processing. ``'none'`` (default) uses the raw data. ``'mean'`` subtracts the per-image mean. ``'median'`` subtracts the per-image median. :type background: str, optional :returns: * **x** (*1D float array, shape (n_images,)*) -- The x coordinates of the center * **y** (*1D float array, shape (n_images,)*) -- The y coordinates of the center .. py:function:: auto_conv(stack, x_old, y_old, return_conv=False) Recalculate the center of a symmetric object using auto-convolution. For each 2D image of a 3D image stack, use the previous center to select the central row and column. Convolve these against reversed versions of themselves (auto-convolution). Then take the maximum as the new center. Optionally, by setting ``return_conv`` to ``True`` the convolution results can be returned directly, which is useful for sub-pixel fitting. Note: CPU or GPU: The code is agnostic of CPU and GPU usage. If the first parameter is on the GPU the computation/result will be on the GPU. Otherwise, the computation/result will be on the CPU. :param stack: The image stack. The images must be square. :type stack: 3D float array, shape (n_pixels, n_pixels, n_images) :param x_old: Estimated x coordinates near the true centers. :type x_old: 1D float array, shape (n_images,) :param y_old: Estimated y coordinates near the true centers. :type y_old: 1D float array, shape (n_images,) :param return_conv: Whether to return the convolutions instead of the updated centers. The default is ``False``. :type return_conv: bool, optional :returns: Return values differ depending on ``return_conv``: If ``return_conv`` is ``False`` x : 1D float array, shape (n_images,) The x coordinates of the center. y : 1D float array, shape (n_images,) The y coordinates of the center. If ``return_conv`` is ``True`` col_max : 1D int array, shape (n_images,) Indices of the maxima of the column convolutions. row_max : 1D int array, shape (n_images,) Indices of the maxima of the row convolutions. col_con : 2D float array, shape (n_pixels, n_images) Column convolutions (unchanged orientation). row_con : 2D float array, shape (n_images, n_pixels) Row convolutions; note the axes are ordered ``(n_images, n_pixels)``. :rtype: tuple of ndarray .. py:function:: auto_conv_sub_pixel(stack, x_old, y_old, n_local=5) Re-calculate center of symmetric object by auto-convolution sub-pixel fit For each 2D image of a 3D image-stack: use the previous center to select the central row and column. Convolve these against themselves. Use several points around the maximum of the convolution to fit a parabola and use the vertex of the parabola as the center to find the sub-pixel coordinates. Note: CPU or GPU: The code is agnostic of CPU and GPU usage. If the parameters are on the GPU the computation/result will be on the GPU. Otherwise, the computation/result will be on the CPU. :param stack: The image-stack. Note, the images must be square. :type stack: 3D float array, shape (n_pixels, n_pixels, n_images) :param x_old: Estimated x coordinates near the true centers. :type x_old: 1D float array, shape (n_images) :param y_old: Estimated y coordinates near the true centers. :type y_old: 1D float array, shape (n_images) :param n_local: The number of local points around the vertex to be used in parabolic fitting. Must be an odd int >=3. :type n_local: int :returns: * **x** (*1D float array, shape (n_images,)*) -- The x coordinates of the center. * **y** (*1D float array, shape (n_images,)*) -- The y coordinates of the center. .. py:function:: auto_conv_multiline(stack, x_old, y_old, line_ratio=0.05, return_conv=False) Re-calculate center of symmetric object by multi-line auto-convolution For each 2D image of a 3D image-stack: use the previous center to select multiple rows and columns (determined by ``line_ratio``). Average the resulting signals, convolve them against themselves (auto-convolution) Then take the maximum as the new center. Optionally, by setting return_conv to True the convolution results can be returned directly. This is useful for sub-pixel fitting. Note: CPU or GPU: The code is agnostic of CPU and GPU usage. If the first parameter is on the GPU the computation/result will be on the GPU. Otherwise, the computation/result will be on the CPU. :param stack: The image-stack. The images must be square. :type stack: 3D float array, shape (n_pixels, n_pixels, n_images) :param x_old: Estimated x coordinates near the true centers. :type x_old: 1D float array, shape (n_images) :param y_old: Estimated y coordinates near the true centers. :type y_old: 1D float array, shape (n_images) :param line_ratio: Fraction of the frame width that determines how many neighbouring lines are averaged before convolution. :type line_ratio: float, optional :param return_conv: Whether to return the convolution or return the new center. The default is False. :type return_conv: bool, optional :returns: * *tuple* -- see information below * *If return_conv is False* -- x : 1D float array, shape (n_images,) The x coordinates of the center y : 1D float array, shape (n_images,) The y coordinates of the center * *If return_conv is True* -- col_max : 1D int array, shape (n_images,) The index of the maximum of the column convolution row_max : 1D int array, shape (n_images,) The index of the maximum of the row convolution col_con : 2D float array, shape (n_pixels, n_images) The column convolution row_con : 2D float array, shape (n_images, n_pixels) The row convolution .. py:function:: auto_conv_multiline_sub_pixel(stack, x_old, y_old, line_ratio=0.1, n_local=5) Re-calculate center of symmetric object by multi-line auto-convolution with sub-pixel fit For each 2D image of a 3D image-stack: use the previous center to select multiple rows and columns (determined by ``line_ratio``). Average the resulting signals, convolve them against themselves (auto-convolution) and use several points around the maximum to fit a parabola. The vertex of the parabola is used to determine the sub-pixel coordinates of the center. Note: CPU or GPU: The code is agnostic of CPU and GPU usage. If the first parameter is on the GPU the computation/result will be on the GPU. Otherwise, the computation/result will be on the CPU. :param stack: The image-stack. The images must be square. :type stack: 3D float array, shape (n_pixels, n_pixels, n_images) :param x_old: Estimated x coordinates near the true centers. :type x_old: 1D float array, shape (n_images) :param y_old: Estimated y coordinates near the true centers. :type y_old: 1D float array, shape (n_images) :param line_ratio: The ratio relative to the total image width of lines to be used in the convolutions. :type line_ratio: float, optional :param n_local: The number of local points around the vertex to be used in parabolic fitting. Must be an odd int >=3. :type n_local: int, optional :returns: * **x** (*1D float array, shape (n_images,)*) -- The x coordinates of the center. * **y** (*1D float array, shape (n_images,)*) -- The y coordinates of the center. .. py:function:: radial_profile(stack, x, y, oversample=1) Calculate the average radial profile about a center For each 2D image of a 3D image-stack: calculate the average radial profile about the corresponding center (x and y). The profile is calculated by binning. For each pixel in an image the Euclidean distance from the center is calculated. The distance is then used to bin each pixel. When ``oversample`` equals 1 the bin widths are 1 pixel wide; higher values split each native bin into finer ``1 / oversample`` pixel slices. The bins are then normalized by the number of pixels in each bin to find the average intensity in each bin. The number of bins (n_bins) is ((stack.shape[0] // 2) * oversample). Note: CPU or GPU: The code is agnostic of CPU and GPU usage. If the first parameter is on the GPU the computation/result will be on the GPU. Otherwise, the computation/result will be on the CPU. :param stack: The image-stack. Note, the images must be square. :type stack: 3D float array, shape (n_pixels, n_pixels, n_images) :param x: x-coordinates of the center. :type x: 1D float array, shape (n_images) :param y: y-coordinates of the center. :type y: 1D float array, shape (n_images) :param oversample: Oversampling factor applied to the radial distances before binning. Increasing the factor multiplies the number of radial bins and thus the resolution of the profile by the same amount. Must be an integer greater than or equal to 1. :type oversample: int, optional :returns: **profiles** -- The average radial profile of each image about the center :rtype: 2D float array, shape (n_bins, n_images) .. py:function:: fft_profile(stack, x, y, oversample=4, rmin=0.0, rmax=0.5, gaus_factor=6.0) Compute FFT-based radial intensity profiles for a stack of images. The images are first weighted in-place by a 2D Gaussian centered at the requested locations. A real 2D FFT is then evaluated for each weighted image, and the magnitude spectrum is azimuthally averaged into oversampled radial bins. :param stack: Image stack to profile. The images must be square with an even width. This array is modified in-place by Gaussian weighting prior to the FFT step. :type stack: array_like, shape (n_pixels, n_pixels, n_images) :param x: X-coordinates of the Gaussian centers in pixel units. :type x: array_like, shape (n_images,) :param y: Y-coordinates of the Gaussian centers in pixel units. :type y: array_like, shape (n_images,) :param oversample: Radial oversampling factor (>=1) applied when binning FFT magnitudes. :type oversample: int, default=4 :param rmin: Minimum normalized radial frequency (0–0.5 Nyquist) to keep in the returned profile. :type rmin: float, default=0.0 :param rmax: Maximum normalized radial frequency (0–0.5 Nyquist) considered when building the radial profile. :type rmax: float, default=0.5 :param gaus_factor: Divisor controlling the Gaussian width relative to the image size. :type gaus_factor: float, default=6.0 :returns: **profile** -- Oversampled radial magnitude profiles for each image, sliced to the bins corresponding to the radial range ``[rmin, rmax]``. :rtype: array_like, shape (n_selected_bins, n_images) .. py:function:: lookup_z(profiles, zlut, n_local=5) Calculate the corresponding sub-planar z-coordinate of each profile by LUT For each image's profile in ``profiles``: find the best matching profile in the Z-LUT (lookup table). The lookup table stores a leading row of z-coordinates and radial profiles below; the first radial bin corresponds to the central pixel and is ignored during correlation (hence ``zlut[2:, :]``). Fits the local points around the best matching profile to find sub-planar fit in between columns of the LUT. Note: CPU or GPU: The code is agnostic of CPU and GPU usage. If the first parameter is on the GPU the computation/result will be on the GPU. Otherwise, the computation/result will be on the CPU. :param profiles: The average radial profile of each image about the center :type profiles: 2D float array, shape (n_bins, n_images) :param zlut: The reference radial profiles and corresponding z-coordinates. The first row (``zlut[0, :]``) holds the z-axis values. The remaining rows contain the reference radial profiles produced by :func:`radial_profile`; their first bin (``zlut[1, :]``) is skipped to avoid the central pixel when matching. :type zlut: 2D float array, shape (1+n_bins, n_ref) :param n_local: The number of local points around the vertex to be used in parabolic fitting. Must be an odd int >=3. Default is 5. :type n_local: int, optional :returns: **z** -- z-coordinates :rtype: 1D float array, shape (n_images) .. py:function:: stack_to_xyzp_advanced(stack, zlut=None, **kwargs) Calculate image-stack XYZ and profiles (Z is nan if Z-LUT is None) :param stack: The image-stack. Note, the images must be square. It is expected it is in the regular CPU memory. It will be transferred to the GPU. :type stack: 3D float array, shape (n_pixels, n_pixels, n_images) :param zlut: The reference radial profiles and corresponding z-coordinates. The first row (zlut[0, :]) are the z-coordinates. The rest of the rows are the corresponding profiles as generated by radial_profile. It is expected it is already in the GPU memory. The defualt is None. :type zlut: 2D float array, shape (1+n_bins, n_ref), optional :param \*\*kwargs: Additional keyword arguments controlling individual processing stages. The following keys are recognised: * ``"center_of_mass"`` (dict, default ``{}``): forwarded to :func:`center_of_mass`. * ``"auto_conv"`` (dict, default ``{}``): forwarded to :func:`auto_conv`. * ``"n auto_conv_multiline_sub_pixel"`` (int, default ``5``): number of :func:`auto_conv_multiline_sub_pixel` refinement iterations. * ``"auto_conv_multiline_sub_pixel"`` (dict, default ``{}``): forwarded to :func:`auto_conv_multiline_sub_pixel`. * ``"use fft_profile"`` (bool, default ``True``): when ``True`` compute profiles with :func:`fft_profile`; when ``False`` compute profiles with :func:`radial_profile`. * ``"fft_profile"`` (dict, default ``{}``): forwarded to :func:`fft_profile` when ``use fft_profile`` is ``True``. * ``"radial_profile"`` (dict, default ``{}``): forwarded to :func:`radial_profile` when ``use fft_profile`` is ``False``. * ``"lookup_z"`` (dict, default ``{}``): forwarded to :func:`lookup_z` when ``zlut`` is provided. :type \*\*kwargs: dict, optional :returns: * **x** (*1D float array, shape (n_images)*) -- x-coordinates * **y** (*1D float array, shape (n_images)*) -- y-coordinates * **z** (*1D float array, shape (n_images)*) -- z-coordinates * **profiles** (*2D float array, shape (n_bins, n_images)*) -- The average radial profile of each image about the center .. py:function:: stack_to_xyzp(stack, zlut=None) Estimate XYZ coordinates and radial profiles from an image stack. This convenience wrapper orchestrates the CPU/GPU-agnostic pipeline used throughout MagTrack: the x and y are first estimated with :func:`center_of_mass`, refined with :func:`auto_conv`, further refined by five iterations of :func:`auto_conv_multiline_sub_pixel`, and then converted into radial profiles via :func:`radial_profile`. When a Z-look-up table is provided, :func:`lookup_z` translates those profiles into axial coordinates; otherwise, NaNs are returned for the z. Note: CPU or GPU: The code is agnostic of CPU and GPU usage. If the first parameter is on the GPU the computation/result will be on the GPU. Otherwise, the computation/result will be on the CPU. :param stack: 3-D image stack containing square images. The array can reside on the CPU (NumPy) or GPU (CuPy). :type stack: array-like, shape (n_pixels, n_pixels, n_images) :param zlut: Radial-profile look-up table whose first row stores the reference z-positions and remaining rows contain the corresponding template profiles. If omitted, the axial coordinate output is filled with NaNs. :type zlut: array-like, shape (1 + n_bins, n_ref), optional :returns: * **x** (*1D float array, shape (n_images)*) -- x-coordinates * **y** (*1D float array, shape (n_images)*) -- y-coordinates * **z** (*1D float array, shape (n_images)*) -- z-coordinates or a NaN array when ``zlut`` is None * **profiles** (*2D float array, shape (n_bins, n_images)*) -- The average radial profile of each image about the center