Add a lifetime calculator module

[phasorpy]

Add a lifetime calculator module

[Mariochem92]

Thank you for the update! I appreciate your consideration of the input and code.

Regarding the use of masks, it sounds like a useful approach, especially for handling images. In this case I try to handle segments of the contour plot on the phasor cloud. I know it's customary to go with normal distribution approaches, but I had in mind a fully non-parametric contour plot in 2D. In case you suggest to mask the phasor plot that actually sounds interesting!

Separating the histogramming, analysis of histograms/distributions, and plotting into modules sounds like a great idea for better organization and modularity.

Please take your time to address these priorities, and I'll be available to eventually discuss when the time is right.

[Mariochem92]

Hello everyone,

I'd like to share a set of Python functions in this code snippet for the proposed phasorpy.lifetime module, with the hope that they can be of assistance in addressing some of the features described in the project request.

The code relies on a few external libraries:

• NumPy: Used for various mathematical operations, including array handling and calculations.

• Shapely: Required for working with polygons and geometric operations. Install it using pip install shapely.

• Matplotlib: Utilized for generating plots and visualizing data. Install it using pip install matplotlib.

Here's an overview of the key functions and their roles in the context of your project:

1. phasor_coordinates_FD(phi, mod):

   • Purpose: Converts phase and modulation data into phasor coordinates in the frequency domain.

2. ph_coordinates2phmd(g, s):

   • Purpose: Converts phasor coordinates back into phase and module.

3. polished_hist(g, s, hist_bins):

   • Purpose: Generates a polished 2D histogram from provided phasor data.

4. contour_plot(histogram, levels, limits):

   • Purpose: Creates a contour plot from a histogram with specified contour levels and limits.

5. analyze_contourplot(contour_segments):

   • Purpose: Analyzes contour segments of the phasor plot histogram and identifies both phasors and regions of interest.

6. simple_lt(g, s, flim_frequency):

   • Purpose: Calculates lifetime values from phasor coordinates and FLIM frequency, as τ_m and τ _φ

7. lifetime_computation(g, s, flim_frequency):

   • Purpose: Computes lifetime values and standard deviations from phasor coordinates and FLIM frequency with error propagation

8. refplot(names, lifetimes, flim_frequency):

   • Purpose: Generates reference phasor coordinates based on sample lifetimes and FLIM frequency.

9. Segment_IntensityFractions(RefA, RefB, target, labels):

   • Purpose: Calculates intensity fractions between two reference points and a target phasor.

I've left the FRET-related functionalities as they are, as I'm not experienced in that domain. If you have any specific guidance or if you'd like to explore an object-oriented approach, please let me know.

Warm regards,
Mario

---

import numpy as np
from shapely.geometry import Polygon
import matplotlib.pyplot as plt


def phasor_coordinates_FD(phi, mod):
# Convert phase and module information into phasor coordinates in the frequency domain.
# Parameters:
# - phi: Phase angle in radians.
# - mod: Module.
# Returns:
# - g: Real part of the phasor coordinates.
# - s: Imaginary part of the phasor coordinates.
#
# This function takes phase and moduleinformation (phi and mod) and computes the corresponding phasor
# coordinates (g, s) in the phasor plot. It is valuable for converting phase and modulation data
# into phasor representations.
  g = mod * np.cos(phi)
  s = mod * np.sin(phi)
  return g, s

def ph_coordinates2phmd(g, s):
# Convert phasor coordinates into phase and module
# Parameters:
# - g: Real part of the phasor coordinates.
# - s: Imaginary part of the phasor coordinates.
# Returns:
# - phi: Phase angle in radians.
# - mod: Module.
#
# This function takes phasor coordinates (g, s) and computes the corresponding phase angle (phi) in radians
# and module (mod). It is useful for transforming phasor data into more interpretable phase and
# module information for further analysis or visualization.
  phi = np.arctan2(s, g)
  mod = np.sqrt(g**2 + s**2)
  return phi, mod


def polished_hist(g,s,hist_bins):
# Generate a 2D histogram from provided phasor data.
# Parameters:
# - g: Real part of data points.
# - s: Imaginary part of data points.
# - hist_bins: Number of bins for histogram generation.
# Returns:
# - hist: The 2D histogram data.
# - filtered_array: A binary mask indicating valid data points.
# - limits: Extent limits defining the data range.
# - hist_color: Histogram values at specific data points.

# This function creates a polished 2D histogram using provided data points, with a specified number of bins.
# It returns the histogram data, a binary mask indicating valid data points, extent limits defining the data range,
# and histogram values at specific data points. This polished histogram is useful for visualizing data distribution
# and further analysis.
  hist, xedges, yedges = np.histogram2d(g, s, bins=hist_bins, range=[[-0.1, 1.1], [-0.1, 0.5]])
  xidx = np.clip(np.digitize(g, xedges), 0, hist.shape[0] - 1)
  yidx = np.clip(np.digitize(s, yedges), 0, hist.shape[1] - 1)
  hist_color = hist[xidx, yidx]
  hist=hist.T
  histmask = np.where(hist > 0, 1, 0)
  filtered_array = np.copy(histmask)
  limits =	[xedges[0],xedges[-1],yedges[0],yedges[-1]]
  return hist,filtered_array,limits, hist_color

def contour_plot(histogram, levels,limits):
# Create a contour plot from a histogram with specified contour levels and limits.
# Parameters:
# - histogram: The input histogram for contour plotting.
# - levels: Contour levels to be displayed/analyzed.
# - limits: Extent limits defining the data range.
# Returns:
# - contour_segments: List of contour segments extracted from the contour plot.

# This function generates a contour plot from the phasor plot histogram using specified contour levels
# and extent limits. It returns a list of contour segments extracted from the generated contour plot,
# allowing further analysis or visualization of the data.
  cs=plt.contour(histogram,levels,extent=limits,linewidths=1,zorder=1)
  contour_segments=cs.allsegs
  return contour_segments

def analyze_contourplot(contour_segments):
# Analyze contour segments of phasor plot histofram and identify regions of interest.
# Parameters:
# - contour_segments: List of contour segments.
# Returns:
# - perimeters: List of identified ROIs segments.
# - phasors: List of identified phasor segments within the ROIs.
#The result is two lists: `phasors` (regions within existing ROIs) and `perimeters` (regions considered to be ROIs).
#- It checks if each contour is part of an existing perimeter (ROI) or a new region of interest.
#- If a contour is within a perimeter and has a valid shape, it's marked as part of an ROI.
#- If the ROI already contains a phasor (distinct subregion), the function updates the phasor with the new contour or creates a new phasor if none exists.
#- If the contour is not part of any ROI, and it's a valid shape, it's added to the list of perimeters.
  perimeters = []
  phasors = []


  for level, contours in enumerate(contour_segments):
      for contour in contours:
          polygon = Polygon(contour)
          roi_found = False
          phasor_found = False

          for idx, perimeter in enumerate(perimeters):
              perimeter_polygon = Polygon(perimeter)
              if perimeter_polygon.contains(polygon) and len(contour)>2:
                  roi_found = True
                  if phasors:
                    for idxp,phasor in enumerate(phasors):
                      phasor_polygon = Polygon(phasor)
                      if phasor_polygon.contains(polygon):
                          phasor_polygon = Polygon(phasor)
                          phasors[idxp]=contour
                          phasor_found = True
                    if phasor_found == False:
                      phasors.append(contour)
                  else:
                    phasors.append(contour)
          if not roi_found:
              if len(contour) > 2:
                  perimeters.append(contour)                   
  return phasors, perimeters

def simple_lt(g, s, flim_frequency):
# Calculate lifetime values from phasor coordinates and FLIM frequency.
# Parameters:
# - g: Real part of the phasor coordinates.
# - s: Imaginary part of the phasor coordinates.
# - flim_frequency: FLIM (Fluorescence Lifetime Imaging Microscopy) frequency.
# Returns:
# - tau_m: Lifetime modulation.
# - tau_phi: Lifetime phase.
#
# This function computes lifetime values (tau_m and tau_phi) based on the provided phasor coordinates (g, s)
# and the FLIM frequency. It calculates lifetime modulation and phase components using the given formulas.
# In the case of input as a NumPy array, the function handles array-specific calculations and ensures
# that any potential NaN values are set to zero to prevent errors.
  omega = flim_frequency * np.pi * 2
  if isinstance(g, np.ndarray):
    tau_m = np.zeros(np.size(g))
    tau_phi = np.zeros(np.size(g))
  tau_m = np.sqrt((1 - (g**2 + s**2)) / (omega**2 * (g**2 + s**2))) * 10e2
  tau_phi = 1 / omega * s / g * 10e2
  if isinstance(g, np.ndarray):
    tau_m[np.isnan(tau_m) > 0] = 0
    tau_phi[np.isnan(tau_phi) > 0] = 0
  return tau_m, tau_phi

def lifetime_computation(g, s, flim_frequency):
# Compute lifetime values and standard deviations from phasor coordinates and FLIM frequency.
# Parameters:
# - g: Real part of the phasor coordinates.
# - s: Imaginary part of the phasor coordinates.
# - flim_frequency: FLIM (Fluorescence Lifetime Imaging Microscopy) frequency.
# Returns:
# - lt: List containing lifetime values, [lt_m, lt_phi].
# - lt_std: List containing standard deviations for lifetime values, [lt_m_std, lt_phi_std].
# - phasor_position: List with the average phasor coordinates, [g_av, s_av].
# - phasor_std: List with standard deviations for phasor coordinates, [g_std, s_std].
#
# This function calculates lifetime values and their standard deviations based on phasor coordinates (g, s)
# and the FLIM frequency. It computes average phasor coordinates and their standard deviations,
# as well as lifetime values and their associated standard deviations. The function also checks
# for cases where the phasor position may prevent the determination of a finite lifetime value.
  omega = flim_frequency * 2 * np.pi
  lt = [0, 0]
  lt_std = [0, 0]
  lt_m = []
  lt_phi = []

  g_av = np.average(g)
  s_av = np.average(s)
  phasor_position = [g_av, s_av]
  g_var = np.average((g - g_av)**2)
  s_var = np.average((s - s_av)**2)
  g_std = np.sqrt(g_var)
  s_std = np.sqrt(s_var)
  phasor_position = [g_av, s_av]
  phasor_std = [g_std, s_std]
  lt_m, lt_phi = simple_lt(g_av, s_av, flim_frequency)
  lt_phi_std = np.sqrt((g_std / g_av)**2 + (s_std / s_av)**2)
  lt_m_std = 2*g_std/np.sqrt(g_av**2)+ 2*s_std/np.sqrt(s_av**2)
  lt = [lt_m, lt_phi]
  lt_std = [lt_m_std, lt_phi_std]
  if not lt_m:
      print('Due to the phasor position it was not possible to determine a lifetime finite value')
  return lt, lt_std, phasor_position, phasor_std

def refplot(names, lifetimes, flim_frequency):
# Generates reference phasor coordinates based on sample lifetimes and FLIM frequency.
# Parameters:
# - names: List of sample names.
# - lifetimes: List of lifetime values corresponding to the samples.
# - flim_frequency: FLIM (Fluorescence Lifetime Imaging Microscopy) frequency.
# Returns:
# - g_ref: Real part of the reference phasor coordinates.
# - s_ref: Imaginary part of the reference phasor coordinates.
#
# This function computes reference phasor coordinates for each sample based on their respective
# lifetime values and the FLIM frequency. The phasor coordinates are computed as (g_ref, s_ref)
# and represent the real and imaginary parts of the phasor, respectively. These coordinates can
# be used for comparative analysis, work with lifetime mixtures and visualization of the samples in the context of FLIM data.

  g_ref = []
  s_ref = []
  omega = flim_frequency * np.pi * 10e-3
  for tau in lifetimes:
    M = 1 / np.sqrt(1 + (omega * tau)**2)
    phi = np.arctan(omega * tau)
    g_ref.append(M * np.cos(phi))
    s_ref.append(M * np.sin(phi))
  return g_ref, s_ref

def Segment_IntensityFractions(RefA, RefB, target, labels):
    # Calculate intensity fractions between two reference points and a target point.
    # Parameters:
    # - RefA: Coordinates of the first reference point.
    # - RefB: Coordinates of the second reference point.
    # - target: Coordinates of the target phasor.
    # - labels: A dictionary for labeling categories.
    # Returns a list of intensity fractions.

    Tot_distance = ((RefA[0] - RefB[0])**2 + (RefA[1] - RefB[1])**2)**0.5
    distanceA = ((RefA[0] - target[0])**2 + (RefA[1] - target[1])**2)**0.5
    distanceB = ((RefB[0] - target[0])**2 + (RefB[1] - target[1])**2)**0.5

    fractionA = distanceB / (distanceA + distanceB)
    fractionB = distanceA / (distanceA + distanceB)

    ratios = [fractionA, fractionB]
    return ratios

[cgohlke]

Thanks for the input and code. That functionality will certainly be integrated in this library. Histogramming, analysis of histograms/distributions, and plotting will be separated into modules.

Numpy and matplotlib are already core requirements. Rather than working directly with shapes, we might want to work via masks.

I'm currently focusing on I/O, #10, #11, and other formats. It will take a while before I can start working on this issue.