feat: acquisition tracks (sequential vs. simultaneous channels)

docs/stages.md

@@ -214,11 +214,47 @@ where:
 - $\text{QE}(\lambda)$ is the quantum efficiency of the detector at each wavelength.
 
 These rates are for a fluorophore at the *reference* illumination intensity (the
-channel's `power`, i.e. `OpticalConfig.power`, in W/cm²): anywhere in a widefield
+channel's `power`, i.e. `OpticalConfig.power`, in W/cm²; for channels acquired
+together, the sum over the track, see below): anywhere in a widefield
 field of view, or at the center of the focused spot in a confocal.  By default they include excitation saturation
 at that intensity; use `filtered_emission_rates(saturate=False)` for the rates
 without it.
 
+!!! info "Sequential vs. simultaneous acquisition (tracks)"
+
+    By default, channels are acquired one after another, so each channel is
+    excited only by its own light sources.  `Simulation.tracks` groups channels
+    that are acquired at the same time (a *track*): all lasers of all channels in
+    a track are on together, and each channel's detector records at once, so
+    every fluorophore can bleed into every channel of its track.  A laser listed
+    in several channels of a track is one laser, counted once, so those channels
+    must agree on its `power` and excitation filters.  On a point scanner, the
+    channels of a track also share one pixel dwell time (`exposure_ms`).
+
+    Each channel still uses its own filters (including its dichroic) to detect
+    light.  On a real scanner, the channels of a track share one (multiband)
+    dichroic, so give them matching dichroics if that matters for your
+    simulation.
+
+    ```python
+    from microsim.schema.optical_config.lib import (
+        LASER_405,
+        LASER_488,
+        LASER_561,
+        LASER_640,
+    )
+
+    sim = Simulation(
+        # ...,
+        channels=[LASER_405, LASER_488, LASER_561, LASER_640],
+        tracks=[[0, 2], [1, 3]],  # or by name: [["405", "561"], ["488", "640"]]
+        # tracks="simultaneous",  # all channels at once
+    )
+    ```
+
+    Each list is shorthand for a [`Track`][microsim.schema.Track], e.g.
+    `Track(channels=[0, 2])`.
+
 The output of this stage (given by
 [`Simulation.filtered_emission_rates`][microsim.Simulation.filtered_emission_rates])
 is a 3D array with dimensions `(C, F, W)`, but lacking spatial information:
@@ -365,8 +401,11 @@ which retains the full `(C, F, Z, Y, X)` dimensions.
       wavelengths (with several lasers or a broadband source, each wavelength is
       weighted by how strongly the fluorophore absorbs it).  Excitation
       saturation is applied locally, so it is strongest at the center of the
-      spot.  `pinhole_au` is in Airy units at the center of the channel's
-      emission filter band (or at `airy_unit_wavelength_nm`, if set).
+      spot.  As on many point scanners, channels acquired together (in the same
+      track, see below) are detected through one pinhole: `pinhole_au` is in Airy
+      units at the average of their emission filter band centers (or at
+      `airy_unit_wavelength_nm`, if set).  With `shared_pinhole=False`, each
+      channel instead has its own pinhole, matched to its emission band.
 
 The output of this stage is a 4D array with dimensions `(C, Z, Y, X)` where the
 wavelength and fluorophore dimensions have been collapsed into the channel. Note

examples/multi_channel_confocal.py

@@ -6,9 +6,11 @@
 from microsim.util import ndview
 
 # short stubby lines labeled in EGFP and long thin lines labeled in EYFP, imaged
-# with 488 and 515 nm lasers.  EYFP is also excited at 488 nm and emits in the
-# 525/50 band, so it bleeds through into the 488 channel (about half as bright as
-# EGFP there); much less EGFP bleeds through into the 515 channel.
+# with 488 and 515 nm lasers at the same time (one track), so each fluorophore
+# bleeds into the other's channel (in the 488 channel, EYFP even outshines EGFP).
+# Try `tracks="sequential"`: 515 nm barely excites EGFP, so its bleedthrough into
+# the EYFP channel mostly vanishes (EYFP still bleeds into the 488 channel, since
+# 488 nm excites it too).
 
 sim = ms.Simulation(
     truth_space=ms.ShapeScaleSpace(shape=(64, 256, 256), scale=(0.04, 0.02, 0.02)),
@@ -29,6 +31,7 @@
     settings=ms.Settings(random_seed=100),
     detector=HYD_SP8,
     exposure_ms=0.002,  # 2 µs pixel dwell time
+    tracks="simultaneous",  # both lasers on, both detectors recording
     channels=[LASER_488, LASER_515],
     # output_path="au1.tif",
 )

src/microsim/schema/__init__.py

@@ -16,7 +16,7 @@
 )
 from .sample import Fluorophore, FluorophoreDistribution, MatsLines, Sample
 from .settings import Settings
-from .simulation import Simulation
+from .simulation import Simulation, Track
 from .space import DownscaledSpace, ExtentScaleSpace, ShapeExtentSpace, ShapeScaleSpace
 from .spectrum import Spectrum
 
@@ -56,6 +56,7 @@
     "Simulation",
     "Spectrum",
     "SpectrumFilter",
+    "Track",
     "Widefield",
 ]
 

src/microsim/schema/modality/_simple_psf.py

@@ -43,9 +43,12 @@ def _illumination(self) -> Illumination:
         return WidefieldIllumination()
 
     def _detection(
-        self, channel: Sequence[ChannelFluorescence], objective: ObjectiveLens
+        self,
+        fluorescence: Sequence[Sequence[ChannelFluorescence]],
+        c_idx: int,
+        objective: ObjectiveLens,
     ) -> Detection:
-        """Detection for one channel (`channel` has one item per fluorophore)."""
+        """Detection for channel `c_idx` (`fluorescence` is `[channel][fluor]`)."""
         return CameraDetection()
 
     def render(
@@ -75,7 +78,7 @@ def render(
                     if fl.detected_rate(saturate=False) <= 1e-12:
                         continue  # this fluorophore is not seen in this channel
                     if detection is None:  # (only needed if something is seen)
-                        detection = self._detection(channel, objective_lens)
+                        detection = self._detection(fluorescence, c_idx, objective_lens)
                     logger.info(f"Fluor: {fluor}")
                     with logging_indented():
                         key, psf, bright = self._system_psf(
@@ -216,7 +219,16 @@ class Confocal(_PSFModality):
         default=None,
         gt=0,
         description="Wavelength (nm) at which `pinhole_au` is defined.  By default, "
-        "1 AU is matched to the center of each channel's emission filter band.",
+        "1 AU is matched to the average emission filter band center of the channels "
+        "sharing the pinhole (or, if not `shared_pinhole`, to each channel's band "
+        "center).",
+    )
+    shared_pinhole: bool = Field(
+        default=True,
+        description="If True (the default), channels acquired together (in the same "
+        "track, see `Simulation.tracks`) are detected through the same pinhole, as "
+        "on many point scanners; the pinhole is set anew for each track.  If False, "
+        "each channel has its own pinhole, matched to its emission band.",
     )
 
     detection_type = PinholeDetection
@@ -225,20 +237,33 @@ def _illumination(self) -> Illumination:
         return FocusedIllumination()
 
     def _detection(
-        self, channel: Sequence[ChannelFluorescence], objective: ObjectiveLens
+        self,
+        fluorescence: Sequence[Sequence[ChannelFluorescence]],
+        c_idx: int,
+        objective: ObjectiveLens,
     ) -> Detection:
-        if (wvl_nm := self.airy_unit_wavelength_nm) is None:
-            wvl_nm = channel[0].channel.emission_band_center()
+        # the channels sharing this channel's pinhole: those in the same track
+        track = fluorescence[c_idx][0].track
+        if self.shared_pinhole and track is not None:
+            group = [row for row in fluorescence if row[0].track == track]
+        else:
+            group = [fluorescence[c_idx]]
+        wvl_nm = self.airy_unit_wavelength_nm
+        if wvl_nm is None:
+            centers = [row[0].channel.emission_band_center() for row in group]
+            if None not in centers:
+                wvl_nm = float(np.mean(centers))  # type: ignore [arg-type]
         if wvl_nm is None:
             # no band-limited emission filter (e.g. a longpass): use the center of
-            # the light detected in this channel instead
-            wvl_nm = _detected_light_center(channel)
+            # the light detected in these channel(s) instead
+            wvl_nm = _detected_light_center([fl for row in group for fl in row])
+            names = ", ".join(repr(row[0].channel.name) for row in group)
             # (warnings are shown once per distinct message)
             warnings.warn(
-                f"Channel {channel[0].channel.name!r} has no band-limited emission "
-                f"filter to match the pinhole size to; using 1 AU at {wvl_nm:.0f} nm "
-                "(the center of the detected emission).  Set "
-                "`Confocal.airy_unit_wavelength_nm` to choose.",
+                f"Channel(s) {names} have no band-limited emission filter to match "
+                f"the pinhole size to; using 1 AU at {wvl_nm:.0f} nm (the center of "
+                "the detected emission).  Set `Confocal.airy_unit_wavelength_nm` to "
+                "choose.",
                 stacklevel=2,
             )
         radius_um = self.pinhole_au * airy_radius_um(
@@ -305,18 +330,19 @@ def render(
         )
 
 
-def _detected_light_center(channel: Sequence[ChannelFluorescence]) -> float:
-    """Mean wavelength of all light detected in a channel (all fluorophores)."""
+def _detected_light_center(fluorescence: Sequence[ChannelFluorescence]) -> float:
+    """Mean wavelength of all light detected (all channels and fluorophores given)."""
     total: float = 0
     weighted: float = 0
-    for fl in channel:
+    for fl in fluorescence:
         detected = fl.detected_spectrum(saturate=False)
         total += float(detected.sum())
         weighted += float((detected * detected.coords[Axis.W]).sum())
     if total <= 0:
+        names = sorted({fl.channel.name for fl in fluorescence})
         raise ValueError(
-            "Cannot determine the pinhole size: no light is detected in channel "
-            f"{channel[0].channel.name!r}.  Set `Confocal.airy_unit_wavelength_nm`."
+            "Cannot determine the pinhole size: no light is detected in channel(s) "
+            f"{names}.  Set `Confocal.airy_unit_wavelength_nm`."
         )
     return weighted / total
 

src/microsim/schema/optical_config/_fluorescence.py

@@ -3,10 +3,12 @@
 from dataclasses import dataclass, replace
 from typing import TYPE_CHECKING
 
+import xarray as xr
+
 from microsim.schema.dimensions import Axis
 
 if TYPE_CHECKING:
-    import xarray as xr
+    from collections.abc import Sequence
 
     from microsim.schema.detectors import Detector
     from microsim.schema.lens import ObjectiveLens
@@ -21,7 +23,8 @@ class ChannelFluorescence:
     """How bright one fluorophore is in one channel, and at which wavelengths.
 
     This has no spatial information: it describes a single fluorophore at the
-    *reference* illumination intensity (`OpticalConfig.power`), which is the intensity
+    *reference* illumination intensity (`OpticalConfig.power`, summed over the
+    channels of a track that are acquired together), which is the intensity
     everywhere in a widefield field of view, or at the center of a focused spot.  The
     rates are given with or without saturation at that intensity: modalities that
     don't illuminate evenly (e.g. confocal) apply saturation locally instead.
@@ -35,17 +38,24 @@ class ChannelFluorescence:
     fluorophore : Fluorophore
         The fluorophore.
     absorption_rate : xr.DataArray
-        Photons absorbed per second by one fluorophore, at each excitation wavelength.
+        Photons absorbed per second by one fluorophore, at each excitation wavelength,
+        from all light sources that are on while this channel is acquired (i.e. those
+        of every channel in its track).
     detected_fraction : xr.DataArray
         Photons detected per photon absorbed, at each emission wavelength: quantum
         yield x emission spectrum x emission filters x detector QE x the fraction of
         emitted light collected by the objective.
+    track : int | None
+        Index of the track this channel belongs to: the group of channels acquired
+        at the same time (see `Simulation.tracks`).  None if unknown: the channel is
+        then treated as acquired on its own (e.g. with its own confocal pinhole).
     """
 
     channel: OpticalConfig
     fluorophore: Fluorophore
     absorption_rate: xr.DataArray  # (W,)
     detected_fraction: xr.DataArray  # (W,)
+    track: int | None = None
 
     @classmethod
     def create(
@@ -54,14 +64,26 @@ def create(
         fluorophore: Fluorophore,
         objective: ObjectiveLens,
         detector: Detector | None = None,
+        *,
+        absorption_rate: xr.DataArray | None = None,
+        track: int | None = None,
     ) -> ChannelFluorescence:
-        """Describe `fluorophore` in `channel`."""
+        """Describe `fluorophore` in `channel`.
+
+        `absorption_rate` is the fluorophore's absorption from all light sources of
+        the channel's track (see `track_absorption_rate`).  By default, only from
+        `channel` itself.
+        """
         qe = detector.qe if detector else None
         fl = cls.from_channel(channel, fluorophore, detector_qe=qe)
+        if absorption_rate is None:
+            absorption_rate = fl.absorption_rate
         # fraction of isotropic emission entering the objective. PSFs are normalized
         # to collected light, so this is not double-counted (see `detection_psf`)
         detected = fl.detected_fraction * objective.collection_efficiency
-        return replace(fl, detected_fraction=detected)
+        return replace(
+            fl, absorption_rate=absorption_rate, detected_fraction=detected, track=track
+        )
 
     @classmethod
     def from_channel(
@@ -129,3 +151,55 @@ def detected_spectrum(self, *, saturate: bool = True) -> xr.DataArray:
         out.name = "filtered_emission_rate"
         out.attrs.update(long_name="Filtered Emission rate", units="photons/s")
         return out
+
+
+def track_light_sources(track: Sequence[OpticalConfig]) -> list[OpticalConfig]:
+    """Return the channels whose light sources excite `track`, one per light source.
+
+    A light source used by several channels of a track (e.g. two detection windows
+    behind one laser) is on once, so those channels must agree on its power and on
+    the excitation-side filters (e.g. the dichroic).
+    """
+    sources: dict[tuple, OpticalConfig] = {}
+    light_owner: dict[tuple, OpticalConfig] = {}
+    for oc in track:
+        # channels without lights are illuminated through their filters (a lamp)
+        key = (
+            tuple(sorted(oc._light_keys()))
+            if oc.lights
+            else oc._excitation_filter_key()
+        )
+        first = sources.setdefault(key, oc)
+        if first is not oc:
+            if first._excitation_key() != oc._excitation_key():
+                raise ValueError(
+                    f"Channels {first.name!r} and {oc.name!r} are acquired together "
+                    "(in the same track) with the same light source, so they must "
+                    "use the same power and excitation filters (e.g. dichroic)."
+                )
+            continue
+        for light in oc._light_keys():
+            if (other := light_owner.setdefault(light, oc)) is not oc:
+                raise ValueError(
+                    f"Channels {other.name!r} and {oc.name!r} are acquired together "
+                    "(in the same track) and share some, but not all, light sources. "
+                    "Give each light source to one channel (or the same set of light "
+                    "sources to both)."
+                )
+    return list(sources.values())
+
+
+def track_absorption_rate(
+    track: Sequence[OpticalConfig], fluorophore: Fluorophore
+) -> xr.DataArray:
+    """Photons absorbed per second by `fluorophore`, from all light sources of `track`.
+
+    Returns one value per excitation wavelength (dimension `W`).
+    """
+    rates = [
+        oc.absorption_rate(fluorophore).squeeze([Axis.F, Axis.C], drop=True)
+        for oc in track_light_sources(track)
+    ]
+    if len(rates) == 1:
+        return rates[0]
+    return xr.concat(rates, dim="source", join="outer", fill_value=0).sum("source")

src/microsim/schema/optical_config/config.py

@@ -221,6 +221,35 @@ def _crossing(i: int, j: int) -> float:
             return None  # runs off the end of the data: open-ended
         return float(lo + hi) / 2
 
+    def _light_keys(self) -> tuple[tuple[bytes, bytes], ...]:
+        """Identifies each light source, by its spectrum."""
+        return tuple(_spectrum_key(light.spectrum) for light in self.lights)
+
+    def _excitation_filter_key(self) -> tuple:
+        """Identifies the excitation-path filters (in any order)."""
+        ex_side = {Placement.EX_PATH, Placement.BS, Placement.BS_INV, Placement.ALL}
+        return tuple(
+            sorted(
+                (f.placement.value, _spectrum_key(f.spectrum))
+                for f in self.filters
+                if f.placement in ex_side
+            )
+        )
+
+    def _excitation_key(self) -> tuple:
+        """Identifies the excitation side: lights, power and excitation-path filters.
+
+        Channels with equal keys illuminate the sample identically (e.g. two
+        detectors behind the same laser).
+        """
+        lights = sorted(
+            zip(
+                self._light_keys(), (light.power for light in self.lights), strict=True
+            ),
+            key=lambda item: item[0],
+        )
+        return (tuple(lights), self.power, self._excitation_filter_key())
+
     @property
     def excitation(self) -> Filter | None:
         """Combine all excitation filters into a single spectrum."""
@@ -489,3 +518,11 @@ def plot_emission(
         ax.legend()
         ax.set_xlim(400, 800)
         ax.set_xlabel("transmission [%]")
+
+
+def _spectrum_key(spectrum: Spectrum) -> tuple[bytes, bytes]:
+    """Identifies a spectrum by its values (so that e.g. `488` and `488.0` match)."""
+    return (
+        np.asarray(spectrum.wavelength, dtype=float).tobytes(),
+        np.asarray(spectrum.intensity, dtype=float).tobytes(),
+    )