kymotrack: add sample_from_channel

Author: JoepVanlier

changelog.md

@@ -11,6 +11,7 @@
 * Added improved printing of calibration items under `channel.calibration` providing a more convenient overview of the items associated with a `Slice`.
 * Added improved printing of calibrations performed with `Pylake`.
 * Added parameter `titles` to customize title of each subplot in [`Kymo.plot_with_channels()`](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.kymo.Kymo.html#lumicks.pylake.kymo.Kymo.plot_with_channels).
+* Added [`KymoTrack.sample_from_channel()`](https://lumicks-pylake.readthedocs.io/en/latest/_api/lumicks.pylake.kymotracker.kymotrack.KymoTrack.html#lumicks.pylake.kymotracker.kymotrack.KymoTrack.sample_from_channel) to downsample channel data to the time points of a kymotrack.
 
 ## v1.5.2 | 2024-07-24
 

docs/tutorial/figures/kymotracking/sample_from_channel.png

@@ -397,6 +397,32 @@ Here `num_pixels` is the number of pixels to sum on either side of the track.
 .. note::
     For tracks obtained from tracking or :func:`~lumicks.pylake.refine_tracks_centroid`, the photon counts found in the attribute :attr:`~lumicks.pylake.kymotracker.kymotrack.KymoTrack.photon_counts` are computed by :func:`~lumicks.pylake.kymotracker.kymotrack.KymoTrack.sample_from_image` using `num_pixels=np.ceil(track_width / pixelsize) // 2` where `track_width` is the track width used for tracking or refinement.
 
+Averaging channel data over tracks
+----------------------------------
+
+It is also possible to average channel data over the track using :meth:`~lumicks.pylake.kymotracker.Kymotrack.sample_from_channel()`.
+For example, let's find out what the force was during a particular track::
+
+    force_slice = file.force1x
+    track_force = longest_track.sample_from_channel(force_slice, include_dead_time=True)
+
+When you call this function with a :class:`~lumicks.pylake.Slice`, it returns another :class:`~lumicks.pylake.Slice` with the downsampled channel data.
+What happens is that for every point on the track, the correct kymograph scan line is looked up and the channel data is averaged over the entire duration of that scan line.
+The parameter `include_dead_time` specifies whether the time it takes the mirror to return to its initial position after each scan line should be included in this average.
+For tracks which have not been refined (see :ref:`localization_refinement`) the result from this function may skip some scan lines entirely.
+
+We can plot these slices just like any other.
+Plotting this slice, we can see that the protein detaches shortly after the force drops::
+
+    plt.figure()
+    force_slice.plot(label="force (whole file)")
+    track_force.plot(start=force_slice.start, marker=".", label="force (longest track)")
+    plt.legend(loc="upper left")
+
+.. image:: figures/kymotracking/sample_from_channel.png
+
+We plotted the track here putting the time zero at the start time of the entire force plot by passing its `start` time.
+
 Plotting binding histograms
 ---------------------------
 

docs/tutorial/kymotracking.rst

@@ -3,6 +3,7 @@
 import itertools
 from copy import copy
 
+from ..channel import Slice, empty_slice
 from ..__about__ import __version__
 from ..detail.utilities import replace_key_aliases
 from .detail.peakfinding import _sum_track_signal
@@ -592,7 +593,22 @@ def in_rect(self, rect, all_points=False):
         return criterion(np.logical_and(time_match, position_match))
 
     def interpolate(self):
-        """Interpolate KymoTrack to whole pixel values"""
+        """Linearly Interpolates :class:`KymoTrack` to include all time points between the start
+        and end time of the track.
+
+        By default, the kymotracker returns tracks that only include points where the sigal
+        is above a detection thresholds. Consequently, the returned track may miss particular
+        time points where the signal dropped below a certain level. This function interpolates
+        the :class:`KymoTrack` such that *all* intermediate time points are present.
+
+        .. note::
+
+            This function *only* linearly interpolates and does not attempt to estimate where the
+            peak signal intensity is located for those interpolated points. If this is desired
+            instead please use a refinement method such as
+            :func:`~lumicks.pylake.refine_tracks_centroid()` or
+            :func:`~lumicks.pylake.refine_tracks_gaussian()`.
+        """
         interpolated_time = np.arange(int(np.min(self.time_idx)), int(np.max(self.time_idx)) + 1, 1)
         interpolated_coord = np.interp(interpolated_time, self.time_idx, self.coordinate_idx)
         return self._with_coordinates(interpolated_time, interpolated_coord)
@@ -621,6 +637,34 @@ def _split(self, node):
 
         return before, after
 
+    def sample_from_channel(self, channel_slice, include_dead_time=True) -> Slice:
+        """Sample channel data using the time points corresponding to this track
+
+        For each time point in the track, sample the channel data corresponding to that kymograph
+        scan line.
+
+        .. note::
+
+            Tracks may skip kymograph frames when the signal drops below a certain level. If you
+            intend to sample every time point on the track, remember to interpolate the track
+            first using :meth:`KymoTrack.interpolate()`
+
+        Parameters
+        ----------
+        channel_slice : Slice
+            Sample from this channel data.
+        include_dead_time : bool
+            Include the time that the mirror returns to its start position after each kymograph
+            line.
+        """
+        ts_ranges = self._kymo.line_timestamp_ranges(include_dead_time=include_dead_time)
+        try:
+            return channel_slice.downsampled_over(
+                [ts_ranges[time_idx] for time_idx in self.time_idx]
+            )
+        except ValueError:
+            return empty_slice
+
     def sample_from_image(self, num_pixels, reduce=np.sum, *, correct_origin=None):
         """Sample from image using coordinates from this KymoTrack.
 

lumicks/pylake/kymotracker/kymotrack.py

@@ -4,6 +4,8 @@
 import pytest
 
 from lumicks.pylake.kymo import _kymo_from_array
+from lumicks.pylake.channel import Slice, Continuous
+from lumicks.pylake.detail.imaging_mixins import _FIRST_TIMESTAMP
 from lumicks.pylake.kymotracker.kymotrack import KymoTrack
 from lumicks.pylake.kymotracker.kymotracker import track_greedy
 from lumicks.pylake.tests.data.mock_confocal import generate_kymo
@@ -114,3 +116,50 @@ def test_origin_warning_sample_from_image():
         ),
     ):
         tracks[0].sample_from_image(0)
+
+
+@pytest.mark.parametrize(
+    "time_idx, ref_dead_included, ref_dead_excluded",
+    [
+        ([0, 1, 2], [14.5, 24.5, 34.5], [12.5, 22.5, 32.5]),
+        ([0, 2], [14.5, 34.5], [12.5, 32.5]),
+        ([], [], []),
+    ],
+)
+def test_sample_from_channel(time_idx, ref_dead_included, ref_dead_excluded):
+    img = np.zeros((5, 5))
+    kymo = _kymo_from_array(
+        img,
+        "r",
+        line_time_seconds=1.0,
+        exposure_time_seconds=0.6,
+        start=_FIRST_TIMESTAMP + int(1e9),
+    )
+
+    data = Slice(Continuous(np.arange(100), start=_FIRST_TIMESTAMP, dt=int(1e8)))  # 10 Hz
+    kymotrack = KymoTrack(time_idx, time_idx, kymo, "red", 0)
+
+    sampled = kymotrack.sample_from_channel(data)
+    np.testing.assert_allclose(sampled.data, ref_dead_included)
+
+    sampled = kymotrack.sample_from_channel(data, include_dead_time=False)
+    np.testing.assert_allclose(sampled.data, ref_dead_excluded)
+
+
+def test_sample_from_channel_out_of_bounds():
+    kymo = _kymo_from_array(np.zeros((5, 5)), "r", line_time_seconds=1.0)
+    data = Slice(Continuous(np.arange(100), start=0, dt=int(1e8)))
+    kymotrack = KymoTrack([0, 6], [0, 6], kymo, "red", 0)
+
+    with pytest.raises(IndexError):
+        kymotrack.sample_from_channel(data, include_dead_time=False)
+
+
+def test_sample_from_channel_no_overlap():
+    img = np.zeros((5, 5))
+    kymo = _kymo_from_array(img, "r", start=_FIRST_TIMESTAMP, line_time_seconds=int(1e8))
+    data = Slice(Continuous(np.arange(100), start=kymo.stop + 100, dt=int(1e8)))
+    kymotrack = KymoTrack([0, 1, 2], [0, 1, 2], kymo, "red", 0)
+
+    with pytest.raises(RuntimeError, match="No overlap"):
+        _ = kymotrack.sample_from_channel(data)