phasorpy.filter

Filter signals and phasor coordinates.

The phasorpy.filter module provides functions to filter

• phasor coordinates:

  • phasor_filter_gaussian() (not implemented yet)

  • phasor_filter_median()

  • phasor_filter_pawflim()

  • phasor_threshold()

• signals:

  • signal_filter_svd() (spectral vector denoise)

  • signal_filter_median() (not implemented yet)

  • signal_filter_pca() (not implemented yet)

phasorpy.filter.phasor_filter_median(mean, real, imag, /, *, repeat=1, size=3, skip_axis=None, use_scipy=False, num_threads=None, **kwargs)

Return median-filtered phasor coordinates.

By default, apply a NaN-aware median filter independently to the real and imaginary components of phasor coordinates once with a kernel size of 3 multiplied by the number of dimensions of the input arrays. Return the intensity unchanged.

Parameters:

• mean (array_like) – Intensity of phasor coordinates.

• real (array_like) – Real component of phasor coordinates to be filtered.

• imag (array_like) – Imaginary component of phasor coordinates to be filtered.

• repeat (int, optional) – Number of times to apply median filter. The default is 1.

• size (int, optional) – Size of median filter kernel. The default is 3.

• skip_axis (int or sequence of int, optional) – Axes in mean to exclude from filter. By default, all axes except harmonics are included.

• use_scipy (bool, optional) – Use scipy.ndimage.median_filter(). This function has undefined behavior if the input arrays contain NaN values but is faster when filtering more than 2 dimensions. See issue #87.

• num_threads (int, optional) – Number of OpenMP threads to use for parallelization. Applies to filtering in two dimensions when not using scipy. By default, multi-threading is disabled. If zero, up to half of logical CPUs are used. OpenMP may not be available on all platforms.

• **kwargs – Optional arguments passed to scipy.ndimage.median_filter().

Returns:

• mean (ndarray) – Unchanged intensity of phasor coordinates.

• real (ndarray) – Filtered real component of phasor coordinates.

• imag (ndarray) – Filtered imaginary component of phasor coordinates.

Raises:

ValueError – If repeat is less than 0. If size is less than 1. The array shapes of mean, real, and imag do not match.

Examples

Apply three times a median filter with a kernel size of three:

>>> mean, real, imag = phasor_filter_median(
...     [[1, 2, 3], [4, 5, 6], [7, 8, 9]],
...     [[0.0, 0.0, 0.0], [0.5, 0.5, 0.5], [0.2, 0.2, 0.2]],
...     [[0.3, 0.3, 0.3], [0.6, math.nan, 0.6], [0.4, 0.4, 0.4]],
...     size=3,
...     repeat=3,
... )
>>> mean
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])
>>> real
array([[0, 0, 0],
       [0.2, 0.2, 0.2],
       [0.2, 0.2, 0.2]])
>>> imag
array([[0.3, 0.3, 0.3],
       [0.4, nan, 0.4],
       [0.4, 0.4, 0.4]])

phasorpy.filter.phasor_filter_pawflim(mean, real, imag, /, *, sigma=2.0, levels=1, harmonic=None, skip_axis=None)

Return pawFLIM wavelet-filtered phasor coordinates.

This function must only be used with calibrated, unprocessed phasor coordinates obtained from FLIM data. The coordinates must not be filtered, thresholded, or otherwise pre-processed.

The pawFLIM wavelet filter is described in [1].

Parameters:

• mean (array_like) – Intensity of phasor coordinates.

• real (array_like) – Real component of phasor coordinates to be filtered. Must have at least two harmonics in the first axis.

• imag (array_like) – Imaginary component of phasor coordinates to be filtered. Must have at least two harmonics in the first axis.

• sigma (float, optional) – Significance level to test difference between two phasors. Given in terms of the equivalent 1D standard deviations. sigma=2 corresponds to ~95% (or 5%) significance.

• levels (int, optional) – Number of levels for wavelet decomposition. Controls the maximum averaging area, which has a length of \(2^level\).

• harmonic (sequence of int or None, optional) – Harmonics included in first axis of real and imag. If None (default), the first axis of real and imag contains lower harmonics starting at and increasing by one. All harmonics must have a corresponding half or double harmonic.

• skip_axis (int or sequence of int, optional) – Axes in mean to exclude from filter. By default, all axes except harmonics are included.

Returns:

• mean (ndarray) – Unchanged intensity of phasor coordinates.

• real (ndarray) – Filtered real component of phasor coordinates.

• imag (ndarray) – Filtered imaginary component of phasor coordinates.

Raises:

ValueError – If level is less than 0. The array shapes of mean, real, and imag do not match. If real and imag have no harmonic axis. Number of harmonics in harmonic is less than 2 or does not match the first axis of real and imag. Not all harmonics in harmonic have a corresponding half or double harmonic.

References

[1]

Silberberg M, and Grecco H. pawFLIM: reducing bias and uncertainty to enable lower photon count in FLIM experiments. Methods Appl Fluoresc, 5(2): 024016 (2017)

Examples

Apply a pawFLIM wavelet filter with four significance levels (sigma) and three decomposition levels:

>>> mean, real, imag = phasor_filter_pawflim(
...     [[1, 1], [1, 1]],
...     [[[0.5, 0.8], [0.5, 0.8]], [[0.2, 0.4], [0.2, 0.4]]],
...     [[[0.5, 0.4], [0.5, 0.4]], [[0.4, 0.5], [0.4, 0.5]]],
...     sigma=4,
...     levels=3,
...     harmonic=[1, 2],
... )
>>> mean
array([[1, 1],
       [1, 1]])
>>> real
array([[[0.65, 0.65],
        [0.65, 0.65]],
       [[0.3, 0.3],
        [0.3, 0.3]]])
>>> imag
array([[[0.45, 0.45],
        [0.45, 0.45]],
       [[0.45, 0.45],
        [0.45, 0.45]]])

phasorpy.filter.phasor_threshold(mean, real, imag, /, mean_min=None, mean_max=None, *, real_min=None, real_max=None, imag_min=None, imag_max=None, phase_min=None, phase_max=None, modulation_min=None, modulation_max=None, open_interval=False, detect_harmonics=True, **kwargs)

Return phasor coordinates with values outside interval replaced by NaN.

Interval thresholds can be set for mean intensity, real and imaginary coordinates, and phase and modulation. Phasor coordinates smaller than minimum thresholds or larger than maximum thresholds are replaced with NaN. No threshold is applied by default. NaNs in mean or any real and imag harmonic are propagated to mean and all harmonics in real and imag.

Parameters:

• mean (array_like) – Intensity of phasor coordinates.

• real (array_like) – Real component of phasor coordinates.

• imag (array_like) – Imaginary component of phasor coordinates.

• mean_min (array_like, optional) – Lower threshold for mean intensity.

• mean_max (array_like, optional) – Upper threshold for mean intensity.

• real_min (array_like, optional) – Lower threshold for real coordinates.

• real_max (array_like, optional) – Upper threshold for real coordinates.

• imag_min (array_like, optional) – Lower threshold for imaginary coordinates.

• imag_max (array_like, optional) – Upper threshold for imaginary coordinates.

• phase_min (array_like, optional) – Lower threshold for phase angle.

• phase_max (array_like, optional) – Upper threshold for phase angle.

• modulation_min (array_like, optional) – Lower threshold for modulation.

• modulation_max (array_like, optional) – Upper threshold for modulation.

• open_interval (bool, optional) – If true, the interval is open, and the threshold values are not included in the interval. If false (default), the interval is closed, and the threshold values are included in the interval.

• detect_harmonics (bool, optional) – By default, detect presence of multiple harmonics from array shapes. If false, no harmonics are assumed to be present, and the function behaves like a numpy universal function.

• **kwargs – Optional arguments passed to numpy universal functions.

Returns:

• mean (ndarray) – Thresholded intensity of phasor coordinates.

• real (ndarray) – Thresholded real component of phasor coordinates.

• imag (ndarray) – Thresholded imaginary component of phasor coordinates.

Examples

Set phasor coordinates to NaN if mean intensity is smaller than 1.1:

>>> phasor_threshold([1, 2, 3], [0.1, 0.2, 0.3], [0.4, 0.5, 0.6], 1.1)
(array([nan, 2, 3]), array([nan, 0.2, 0.3]), array([nan, 0.5, 0.6]))

Set phasor coordinates to NaN if real component is smaller than 0.15 or larger than 0.25:

>>> phasor_threshold(
...     [1.0, 2.0, 3.0],
...     [0.1, 0.2, 0.3],
...     [0.4, 0.5, 0.6],
...     real_min=0.15,
...     real_max=0.25,
... )
(array([nan, 2, nan]), array([nan, 0.2, nan]), array([nan, 0.5, nan]))

Apply NaNs to other input arrays:

>>> phasor_threshold(
...     [numpy.nan, 2, 3], [0.1, 0.2, 0.3], [0.4, 0.5, numpy.nan]
... )
(array([nan, 2, nan]), array([nan, 0.2, nan]), array([nan, 0.5, nan]))

phasorpy.filter.signal_filter_svd(signal, /, spectral_vector=None, *, axis=-1, harmonic=None, sigma=0.05, vmin=None, dtype=None, num_threads=None)

Return spectral-vector-denoised signal.

The spectral vector denoising algorithm is based on a Gaussian weighted average calculation, with weights obtained in n-dimensional Chebyshev or Fourier space [2].

Parameters:

• signal (array_like) – Hyperspectral data to be denoised. A minimum of three samples are required along axis. The samples must be uniformly spaced.

• spectral_vector (array_like, optional) – Spectral vector. For example, phasor coordinates, PCA projected phasor coordinates, or Chebyshev coefficients. Must be of the same shape as signal with axis removed and an axis containing spectral space appended. If None (default), phasor coordinates are calculated at specified harmonic.

• axis (int, optional, default: -1) – Axis over which spectral_vector is computed if not provided. The default is the last axis (-1).

• harmonic (int, sequence of int, or 'all', optional) – Harmonics to include in calculating spectral_vector. If ‘all’, include all harmonics for signal samples along axis. Else, harmonics must be at least one and no larger than half the number of signal samples along axis. The default is the first harmonic (fundamental frequency). A minimum of harmonic * 2 + 1 samples are required along axis to calculate correct phasor coordinates at harmonic.

• sigma (float, default: 0.05) – Width of Gaussian filter in spectral vector space. Weighted averages are calculated using the spectra of signal items within a spectral vector Euclidean distance of 3 * sigma and intensity above vmin.

• vmin (float, optional) – Signal intensity along axis below which spectra are excluded from denoising.

• dtype (dtype_like, optional) – Data type of output arrays. Either float32 or float64. The default is float64 unless the signal is float32.

• num_threads (int, optional) – Number of OpenMP threads to use for parallelization. By default, multi-threading is disabled. If zero, up to half of logical CPUs are used. OpenMP may not be available on all platforms.

Returns:

Denoised signal of dtype. Spectra with integrated intensity below vmin are unchanged.

Return type:

ndarray

References

[2]

Harman RC, Lang RT, Kercher EM, Leven P, and Spring BQ. Denoising multiplexed microscopy images in n-dimensional spectral space. Biomed Opt Express, 13(8): 4298-4309 (2022)

Notes

This implementation is considered experimental. It is not validated against the reference implementation and may not be practical with real-world data. See discussion in issue #142.

Examples

Denoise a hyperspectral image with a Gaussian filter width of 0.1 in spectral vector space using first and second harmonic:

>>> signal = numpy.random.randint(0, 255, (8, 16, 16))
>>> signal_filter_svd(signal, axis=0, sigma=0.1, harmonic=[1, 2])
array([[[...]]])