kymotrack: add sample_from_channel

changelog.md

@@ -13,6 +13,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/kymotracking.rst

@@ -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.
+For every point on the track, the corresponding 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 dead time (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
 ---------------------------
 

lumicks/pylake/kymotracker/kymotrack.py

@@ -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,40 @@ def _split(self, node):
 
         return before, after
 
+    def sample_from_channel(
+        self, channel_slice, include_dead_time=True, *, reduce=np.mean
+    ) -> Slice:
+        """Downsample channel data corresponding to the time points of this track
+
+        For each time point in the track, downsample the channel data over the time limits of the
+        corresponding 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.
+        reduce : callable
+            The :mod:`numpy` function which is going to reduce multiple samples into one. The
+            default is :func:`np.mean <numpy.mean>`, but :func:`np.sum <numpy.sum>` could also be
+            appropriate for some cases, e.g. photon counts.
+        """
+        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], reduce=reduce
+            )
+        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/tests/test_image_sampling.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,53 @@ def test_origin_warning_sample_from_image():
         ),
     ):
         tracks[0].sample_from_image(0)
+
+
+@pytest.mark.parametrize(
+    "time_idx, ref_dead_included, ref_dead_excluded, ref_reduce_max",
+    [
+        ([0, 1, 2], [14.5, 24.5, 34.5], [12.5, 22.5, 32.5], [15.0, 25.0, 35.0]),
+        ([0, 2], [14.5, 34.5], [12.5, 32.5], [15.0, 35.0]),
+        ([], [], [], []),
+    ],
+)
+def test_sample_from_channel(time_idx, ref_dead_included, ref_dead_excluded, ref_reduce_max):
+    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)
+
+    sampled = kymotrack.sample_from_channel(data, include_dead_time=False, reduce=np.max)
+    np.testing.assert_allclose(sampled.data, ref_reduce_max)
+
+
+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)