Some tweaks to dF/F code

Author: StephenThornquist

siffpy/core/flim/flimparams.py

@@ -1022,6 +1022,36 @@ def sample(
 
         return np.histogram(samples, bins = x_range)[0]
 
+    def fraction_to_empirical(self, fractions : np.ndarray) -> np.ndarray:
+        """
+        Converts the `fractions` argument into the equivalent empirical lifetime
+        using the exponential parameter fits.
+
+        ## Arguments
+
+        - `fractions : np.ndarray`
+
+            An array of fractions of the total photons in each bin. Should be
+            dimensions (..., self.n_exps - 1), with the presumption that the final
+            bin adds up to 1 (i.e. fractions.sum(axis=-1) <= 1)
+
+        ## Returns
+
+        - `empirical : np.ndarray`
+
+            An array of the empirical lifetimes corresponding to the fractions
+            in the input array.
+        """
+        if self.n_exp == 2:
+            fractions = np.expand_dims(fractions, axis=-1)
+
+        return (
+            fractions @ np.array([exp.tau for exp in self.exps[:-1]]) 
+            + ((1-fractions.sum(axis=-1)) * self.exps[-1].tau)
+        )
+
+
+    
     def __eq__(self, other)->bool:
         """
         Two FLIMParams objects are equal if they have the same parameters

siffpy/core/siffreader.py

@@ -192,10 +192,12 @@ def open(
             self._time_axis_epoch = None
 
         #self.events = io.find_events(self.im_params, self.get_frames_metadata())
-        flim_params = io.load_flim_params(filename)
-        if any(flim_params):
-            self._flim_params = flim_params
-        #print("Finished opening and reading file")
+        try:
+            flim_params = io.load_flim_params(filename)
+            if any(flim_params):
+                self._flim_params = flim_params
+        except Exception as e:
+            print(f"Failed to load FLIM parameters with error: {e}")
 
     def close(self) -> None:
         """ Closes opened file """
@@ -518,16 +520,16 @@ def sec_to_frames(self, seconds : float, base : str = 'volume')->int:
 
         # Arguments
 
-        * seconds : float
-            Time in seconds to convert
+        * `seconds` : *float*
+            * Time in seconds to convert
 
-        * base : str
-            Base to convert to. Can be 'volume' or 'frame'. Defaults to 'volume'.
+        * `base` : *str*
+            * Base to convert to. Can be 'volume' or 'frame'. Defaults to 'volume'.
 
         # Returns
         
-        * frames : int
-            Number of frames that time represents.
+        * `frames` : *int*
+            * Number of frames that time represents.
 
         # Example
 
@@ -544,7 +546,7 @@ def sec_to_frames(self, seconds : float, base : str = 'volume')->int:
         >> 10
         ```
         """
-
+        base = base.lower()
         if base == 'volume':
             return int(seconds/self.dt_volume)
         if base == 'frame':
@@ -558,6 +560,20 @@ def get_frames_metadata(self, frames : Optional[List[int]] = None) -> List[io.Fr
         return [io.FrameMetaData(meta_dict)
             for meta_dict in io.frame_metadata_to_dict(self.siffio.get_frame_metadata(frames=frames))
         ]
+    
+    def epoch_time_to_volume_index(self, epoch_time : int) -> int:
+        """
+        Returns the index of the volume containing the epoch time provided.
+        """
+        return np.searchsorted(self.t_axis(reference_time = 'epoch'), epoch_time)
+
+    def epoch_time_to_frame_index(self, epoch_time : int, correct_flyback : bool = False) -> int:
+        """
+        Returns the index of the frame containing the epoch time provided. Note:
+        by default frames are numbered by _frame triggers_,
+        meaning they include flyback frames in their indexing!
+        """
+        return np.searchsorted(self.get_time(reference_time = 'epoch'), epoch_time)
 
     def epoch_to_frame_time(self, epoch_time : int) -> float:
         """ Converts epoch time to frame time for this experiment (returned in seconds) """
@@ -1326,11 +1342,6 @@ def get_frames_flim(
 
         method = FlimMethod(method)
 
-        if method != FlimMethod.EMPIRICAL:
-            raise NotImplementedError(
-                "Only empirical lifetime method is implemented in `siffio` backend so far"
-            )
-
         registration_dict = self.registration_dict if (
                 registration_dict is None
                 and hasattr(self, 'registration_dict') 
@@ -1354,8 +1365,6 @@ def get_frames_flim(
             method = method.value,
             units = 'countbins',
         )
-        
-
 
     def _get_frames_flim_srm(
         self,
@@ -1382,10 +1391,6 @@ def _get_frames_flim_srm(
             frames = list(range(self.im_params.num_frames))
 
         method = FlimMethod(method)
-        if method != FlimMethod.EMPIRICAL:
-            raise NotImplementedError(
-                "Only empirical lifetime method is implemented in `SiffIO` backend so far"
-            )
 
         registration_dict = self.registration_dict if (
                 registration_dict is None
@@ -1471,6 +1476,10 @@ def sum_mask_flim(
         * `return_framewise` : bool
             If True, does not sum across timepoints, and returns a flat array
             corresponding to the frames.
+
+        * `flim_method` : Union[str, FlimMethod]
+            The method to use for the FLIM analysis. Default is 'empirical'.
+            Options are any of the `FlimMethod` enum values (`siffpy.siffmath.flim.FlimMethod`).
             
         # Returns
 
@@ -1501,10 +1510,6 @@ def sum_mask_flim(
             )
 
         flim_method = FlimMethod(flim_method)
-        if flim_method != FlimMethod.EMPIRICAL:
-            raise NotImplementedError(
-                "Only empirical lifetime method is implemented in `siffio` backend so far"
-            )
 
         timepoint_end = (
             self.im_params.num_timepoints
@@ -1545,7 +1550,7 @@ def sum_mask_flim(
                 summed_flim_data, 
                 intensity = summed_intensity_data,
                 FLIMParams = params,
-                method = 'empirical lifetime',
+                method = flim_method.value,
                 info_string = "ROI",
                 units = FlimUnits.COUNTBINS,
             )
@@ -1554,7 +1559,7 @@ def sum_mask_flim(
             summed_flim_data, 
             intensity = summed_intensity_data,
             FLIMParams = params,
-            method = 'empirical lifetime',
+            method = flim_method.value,
             info_string = "ROI",
             units = FlimUnits.COUNTBINS,
         ).reshape(
@@ -1701,6 +1706,25 @@ def sum_masks_flim(
             If True, returns a flat array of all frames, rather than summing across
             timepoints.
 
+        * `flim_method` : Union[str, FlimMethod]
+            The method to use for the FLIM analysis. Default is 'empirical'.
+            Options are any of the `FlimMethod` enum values (`siffpy.siffmath.flim.FlimMethod`).
+        
+        # Returns
+
+        * FlimTrace
+            TODO :DOCUMENT
+
+        # Example
+
+        TODO
+        ```python
+
+            from siffpy import SiffReader
+            reader = SiffReader('example.siff')
+
+        ```
+
         """
         if self.backend == 'siffreadermodule':
             return self._sum_masks_flim_srm(
@@ -1715,10 +1739,10 @@ def sum_masks_flim(
             )
 
         flim_method = FlimMethod(flim_method)
-        if flim_method != FlimMethod.EMPIRICAL:
-            raise NotImplementedError(
-                "Only empirical lifetime method is implemented in `siffio` backend so far"
-            )
+        # if flim_method != FlimMethod.EMPIRICAL:
+        #     raise NotImplementedError(
+        #         "Only empirical lifetime method is implemented in `siffio` backend so far"
+        #     )
 
         if isinstance(masks, list):
             masks = np.array(masks).squeeze()
@@ -1765,7 +1789,7 @@ def sum_masks_flim(
                 flim_summed,
                 intensity = intensity_summed,
                 FLIMParams = params,
-                method = 'empirical lifetime',
+                method = flim_method.value,
                 info_string = "Multi-ROIs",
                 units = FlimUnits.COUNTBINS,
             )
@@ -1776,7 +1800,7 @@ def sum_masks_flim(
             flim_summed,
             intensity = intensity_summed,
             FLIMParams = params,
-            method = 'empirical lifetime',
+            method = flim_method.value,
             info_string = "Multi-ROIs",
             units = FlimUnits.COUNTBINS,
         ).reshape(
@@ -1985,4 +2009,6 @@ def _rinfo_safe_convert(registration_info : Union[RegistrationInfo, Dict]) -> Di
     if isinstance(registration_info, RegistrationInfo):
         return registration_info.yx_shifts
     if isinstance(registration_info, dict) and 'yx_shifts' in registration_info:
-        return registration_info['yx_shifts']
+        return registration_info['yx_shifts']
+    if isinstance(registration_info, dict):
+        return registration_info

siffpy/siffmath/flim/phasor.py

@@ -1,12 +1,15 @@
 """
 Generic functions for dealing with phasors.
 """
-from typing import Any, Union
+from typing import Any, Union, Tuple
 import numpy as np
 
 from siffpy.core.flim import FLIMParams
 
-def tau_to_phasor(tau : Union[np.ndarray,float], rep_period: float) -> Union[np.ndarray,complex]:
+def tau_to_phasor(
+        tau : Union[np.ndarray,float],
+        rep_period: float,
+    ) -> Union[np.ndarray,complex]:
     """
     Returns the corresponding phasor value for a single exponential
     with the given lifetime and laser repetition period. Tau
@@ -47,7 +50,10 @@ def tau_to_phasor(tau : Union[np.ndarray,float], rep_period: float) -> Union[np.
     s = (2*np.pi*tau/rep_period)/(1 + (2*np.pi*tau/rep_period)**2)
     return g + 1j*s
 
-def phasor_to_tau(phasor : Union[np.ndarray,complex], rep_period : float) -> Union[np.ndarray,float]:
+def phasor_to_tau(
+        phasor : Union[np.ndarray,complex],
+        rep_period : float
+    ) -> Union[np.ndarray,float]:
     """
     Converts a phasor to a tau value, assuming it corresponds to a
     single exponential with the given laser repetition period.
@@ -70,6 +76,43 @@ def phasor_to_tau(phasor : Union[np.ndarray,complex], rep_period : float) -> Uni
     s = phasor.imag
     return rep_period*s/(g*2*np.pi)
 
+def phasor_to_lifetimes(
+        phasor : Union[np.ndarray, complex],
+        rep_period : float,
+    ) -> Tuple[Union[np.ndarray,float], Union[np.ndarray,float]]:
+    """
+    Returns the modulation and phase lifetime values for a given
+    phasor using the relations:
+    - tau_m = sqrt[(1/M)^2 - 1]/(2*pi*f)
+    - tau_phi = tan(phi)/(2*pi*f)
+
+    where M = np.abs(phasor) and phi = np.angle(phasor) and f
+    is the frequency of laser pulses, with the phasor
+    value defined as
+
+    phasor = M * exp( 1j * phi )
+
+    # Arguments
+
+    - `phasor : Union[np.ndarray, complex]`
+        The phasors to transform (as an array or a single value)
+
+    - `rep_period : float`
+        The repetition rate of the laser pulses (in Hz)
+
+    # Returns
+
+    - `lifetimes : (tau_m : float, tau_phi : float)`
+        A tuple of the modulation and phase lifetimes in seconds.
+
+    """
+    m = np.abs(phasor)
+    phi = np.angle(phasor)
+    return (
+        np.sqrt(m**(-2) - 1)*rep_period/(2*np.pi),
+        np.tan(phi)*rep_period/(2*np.pi)
+    )
+
 def phasor_to_fraction(
         phasor : Union[np.ndarray[Any, np.complex128], complex],
         params : FLIMParams,
@@ -101,7 +144,7 @@ def phasor_to_fraction(
 
     `np.ndarray[Any, np.float64]`
         The fraction in each exponential state corresponding to the given phasor.
-        Shape is size(input_phasor) x (params.n_exps - 1), with the last state in
+        Shape is (params.n_exps - 1) x size(input_phasor), with the last state in
         the `FLIMParams` object's `exps` attribute being implicit (1 - sum(fractions)).
     """
 
@@ -116,7 +159,7 @@ def phasor_to_fraction(
     return np.array([
         ((phasor - startpt)/(endpt - startpt)).real
         for endpt, startpt in zip(bounds[:-1], bounds[1:])
-    ]).T.squeeze()
+    ]).T.squeeze().T
 
 def universal_circle() -> np.ndarray[Any, np.complex128]:
     """
@@ -195,6 +238,7 @@ def correct_phasor(
         phasor : np.ndarray[Any, np.complex128],
         params : 'FLIMParams',
         hist_length : int,
+        rotate_by_offset : bool = True,
         subtract_noise : bool = False,
     ) -> np.ndarray[Any, np.complex128]:
     """
@@ -220,6 +264,11 @@ def correct_phasor(
         The length of the histogram used to compute the phasor -- i.e., the
         number of histogram bins between each laser pulse.
 
+    `rotate_by_offset : bool = True`
+        Whether or not to rotate the phasor by the offset in the `FLIMParams`
+        object. If `False`, the phasor will not be corrected for the instrument
+        response function.
+
     `subtract_noise : bool = False`
         Whether or not to additionally project the phasor onto the line
         connecting the exponentials of the `FLIMParams` object.
@@ -265,9 +314,11 @@ def correct_phasor(
         >>> [0.660594+0.36399072j]
     ```
     """
+    rotated_phasor = phasor
     with params.as_units('countbins'):
-        offset = params.tau_offset
-        rotated_phasor = phasor*np.exp(-1j*2*np.pi*offset/hist_length)
+        if rotate_by_offset:
+            offset = params.tau_offset
+            rotated_phasor *= np.exp(-1j*2*np.pi*offset/hist_length)
 
         if subtract_noise:
             if len(params.exps) > 2: