From f2f8f52054f6badbaea74d87c836f9a6f2aded4d Mon Sep 17 00:00:00 2001 From: Talley Lambert Date: Tue, 13 Oct 2026 16:37:42 +0200 Subject: [PATCH] feat: acquisition tracks (sequential vs. simultaneous channels) - New `Simulation.tracks`: "sequential" (default, each channel excited only by its own light sources, as before), "simultaneous", or a list of `Track`s (bare lists of channel indices or names are cast, e.g. [[0, 2], [1, 3]]). All light sources of a track excite every fluorophore, so fluorophores can bleed into every channel of their track; saturation uses the total. - A light source used by several channels of a track is counted once; those channels must agree on its power and excitation filters (ValueError otherwise, or if they share only some light sources). Lights and filters are identified by their spectra (so e.g. `laser(488)` and `laser(488.0)` match, in any order); channels without lights (a lamp) by their excitation filters. - Point scanners: channels in a track must share one pixel dwell time. - New `Confocal.shared_pinhole` (default True): channels in a track share one pinhole, with `pinhole_au` at the average of their emission band centers. - `ChannelFluorescence` gets an optional `track` (None: acquired on its own, e.g. with its own pinhole), and `create(..., absorption_rate=, track=)`. - A rejected `Simulation` field assignment no longer leaves the new value set. - `Track` is exported from `microsim.schema`. Co-Authored-By: Claude Opus 5.5 --- docs/stages.md | 45 +++- examples/multi_channel_confocal.py | 9 +- src/microsim/schema/__init__.py | 3 +- src/microsim/schema/modality/_simple_psf.py | 62 +++-- .../schema/optical_config/_fluorescence.py | 84 +++++- src/microsim/schema/optical_config/config.py | 37 +++ src/microsim/schema/simulation.py | 144 ++++++++++- tests/test_photometry.py | 240 +++++++++++++++++- 8 files changed, 575 insertions(+), 49 deletions(-) diff --git a/docs/stages.md b/docs/stages.md index e2d5e11..8f1345a 100644 --- a/docs/stages.md +++ b/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 diff --git a/examples/multi_channel_confocal.py b/examples/multi_channel_confocal.py index 49676f6..4c318b2 100644 --- a/examples/multi_channel_confocal.py +++ b/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", ) diff --git a/src/microsim/schema/__init__.py b/src/microsim/schema/__init__.py index 42544c0..d5c0aa0 100644 --- a/src/microsim/schema/__init__.py +++ b/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", ] diff --git a/src/microsim/schema/modality/_simple_psf.py b/src/microsim/schema/modality/_simple_psf.py index 1222f11..e08e6e0 100644 --- a/src/microsim/schema/modality/_simple_psf.py +++ b/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 diff --git a/src/microsim/schema/optical_config/_fluorescence.py b/src/microsim/schema/optical_config/_fluorescence.py index 77acf1d..b1bd44c 100644 --- a/src/microsim/schema/optical_config/_fluorescence.py +++ b/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") diff --git a/src/microsim/schema/optical_config/config.py b/src/microsim/schema/optical_config/config.py index 9adc6e4..ea7d97f 100644 --- a/src/microsim/schema/optical_config/config.py +++ b/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(), + ) diff --git a/src/microsim/schema/simulation.py b/src/microsim/schema/simulation.py index 54de93b..a5f8cb2 100644 --- a/src/microsim/schema/simulation.py +++ b/src/microsim/schema/simulation.py @@ -5,12 +5,18 @@ from contextlib import suppress from dataclasses import replace from pathlib import Path -from typing import TYPE_CHECKING, Annotated, Any +from typing import TYPE_CHECKING, Annotated, Any, Literal import numpy as np import pandas as pd import xarray as xr -from pydantic import AfterValidator, Field, model_validator +from pydantic import ( + AfterValidator, + Field, + StrictInt, + ValidationError, + model_validator, +) from pydantic_core import PydanticSerializationError from microsim._data_array import ArrayProtocol, from_cache, to_cache @@ -23,7 +29,11 @@ from .lens import ObjectiveLens from .modality import Modality, Widefield from .optical_config import OpticalConfig, Placement -from .optical_config._fluorescence import ChannelFluorescence +from .optical_config._fluorescence import ( + ChannelFluorescence, + track_absorption_rate, + track_light_sources, +) from .optical_config.lib import FITC from .sample import Fluorophore, FluorophoreDistribution, Sample from .settings import Settings @@ -42,6 +52,7 @@ class SimulationKwargs(TypedDict, total=False): detector: Detector | None modality: Modality settings: Settings + tracks: "Tracks" output_path: "OutPath" | None @@ -52,6 +63,33 @@ def _check_extensions(path: Path) -> Path: OutPath = Annotated[Path, AfterValidator(_check_extensions)] + + +class Track(SimBaseModel): + """Channels acquired at the same time, as on a point scanner. + + All light sources of all channels in a track are on, and each channel's detector + records at once. May also be given as a bare list of channels, e.g. `[0, 2]`. + """ + + channels: list[StrictInt | str] = Field( + min_length=1, + description="Channels in this track: by index in `Simulation.channels`, or " + "by channel name.", + ) + + @model_validator(mode="before") + @classmethod + def _from_list(cls, value: Any) -> Any: + if isinstance(value, list | tuple): + return {"channels": value} + return value + + +Tracks = Literal["sequential", "simultaneous"] | list[Track] +_TRACKS_HINT = ( + " (To change both `channels` and `tracks`, first set `tracks='sequential'`.)" +) xr.set_options(keep_attrs=True) # type: ignore [no-untyped-call] # LRU cache of `Simulation._channel_fluorescence` results (small C x F x W arrays) @@ -73,6 +111,16 @@ class Simulation(SimBaseModel): # Per-pixel integration time: camera exposure, or pixel dwell time for # point-scanning modalities (e.g. Confocal). May be renamed `integration_time`. exposure_ms: float = 100 + tracks: Tracks = Field( + default="sequential", + description="Which channels are acquired at the same time (a *track*, as on " + "a point scanner): all light sources of all channels in a track are on, and " + "each channel's detector records at once. `'sequential'` (the default): " + "each channel is its own track, excited only by its own light sources. " + "`'simultaneous'`: one track with all channels. Or a list of tracks, each a " + "`Track` or a list of channels given by index in `channels` or by name, " + "e.g. `[[0, 2], [1, 3]]`. Every channel must be in exactly one track.", + ) settings: Settings = Field(default_factory=Settings) output_path: OutPath | None = None @@ -106,8 +154,67 @@ def _resolve_spaces(self) -> "Self": self.output_space.reference = self.truth_space return self + @model_validator(mode="after") + def _validate_tracks(self) -> "Self": + point_scanning = not self.modality.detection_type.sums_over_pixel_area + for track in self._track_indices(): + channels = [self.channels[c] for c in track] + track_light_sources(channels) # checks shared light sources + dwells = {oc.exposure_ms or self.exposure_ms for oc in channels} + if point_scanning and len(dwells) > 1: + raise ValueError( + f"Channels {[oc.name for oc in channels]} are acquired together " + "(in the same track), so they must have the same pixel dwell time " + f"(`exposure_ms`); got {sorted(dwells)}." + ) + return self + + def _track_indices(self) -> list[list[int]]: + """Return `tracks` as lists of channel indices.""" + n = len(self.channels) + if self.tracks == "sequential": + return [[i] for i in range(n)] + if self.tracks == "simultaneous": + return [list(range(n))] + names = [oc.name for oc in self.channels] + tracks: list[list[int]] = [] + for track in self.tracks: + indices = [] + for ch in track.channels: + if isinstance(ch, str): + if names.count(ch) != 1: + found = "matches several" if ch in names else "matches no" + raise ValueError( + f"Track channel {ch!r} {found} channel names " + f"({names}). Use channel indices instead." + ) + ch = names.index(ch) + elif not 0 <= ch < n: + raise ValueError( + f"Track channel index {ch} is out of range.{_TRACKS_HINT}" + ) + indices.append(ch) + tracks.append(indices) + listed = sorted(i for track in tracks for i in track) + if listed != list(range(n)): + raise ValueError( + "Every channel must be in exactly one track; got tracks " + f"{tracks} for {n} channels.{_TRACKS_HINT}" + ) + return tracks + def __setattr__(self, name: str, value: Any) -> None: - super().__setattr__(name, value) + old = vars(self).get(name) + was_set = name in self.model_fields_set + try: + super().__setattr__(name, value) + except ValidationError: + # model validators run after the new value is set: restore the old one + if name in type(self).model_fields: + vars(self)[name] = old + if not was_set: + self.model_fields_set.discard(name) + raise if name in type(self).model_fields: # any field may affect the ground truth (sample, spaces, settings...) vars(self).pop("_ground_truth", None) @@ -204,7 +311,8 @@ def _channel_fluorescence(self) -> list[list[ChannelFluorescence]]: Each item describes one fluorophore in one channel, without spatial information: how many photons it absorbs (and at which wavelengths) at the reference illumination intensity, and what fraction of them are detected (and - at which wavelengths). Fluorophores are the unique fluorophores in + at which wavelengths). It is excited by the light sources of all channels in + the same track (see `tracks`). Fluorophores are the unique fluorophores in `sample.labels`. """ fluors = list({lbl.fluorophore: None for lbl in self.sample.labels}) @@ -220,15 +328,22 @@ def _channel_fluorescence(self) -> list[list[ChannelFluorescence]]: for row, oc in zip(cached, self.channels, strict=True) ] - fluorescence = [ - [ - ChannelFluorescence.create( - oc, fluor, self.objective_lens, self.detector - ) - for fluor in fluors - ] - for oc in self.channels - ] + fluorescence: list[list[ChannelFluorescence]] = [[] for _ in self.channels] + for t_idx, track in enumerate(self._track_indices()): + channels = [self.channels[c] for c in track] + for fluor in fluors: + # all light sources of the track excite the fluorophore + absorption = track_absorption_rate(channels, fluor) + for c_idx in track: + fl = ChannelFluorescence.create( + self.channels[c_idx], + fluor, + self.objective_lens, + self.detector, + absorption_rate=absorption, + track=t_idx, + ) + fluorescence[c_idx].append(fl) if key is not None: _FLUORESCENCE_CACHE[key] = fluorescence if len(_FLUORESCENCE_CACHE) > _FLUORESCENCE_CACHE_SIZE: @@ -279,6 +394,7 @@ def _fluorescence_cache_key(self, fluors: list[Fluorophore]) -> tuple | None: try: return ( tuple(oc.model_dump_json() for oc in self.channels), + str(self._track_indices()), tuple(f.model_dump_json() for f in fluors), self.detector.model_dump_json(include={"qe"}) if self.detector else "", self.objective_lens.collection_efficiency, diff --git a/tests/test_photometry.py b/tests/test_photometry.py index 1c6edb8..6887e38 100644 --- a/tests/test_photometry.py +++ b/tests/test_photometry.py @@ -29,7 +29,12 @@ excitation_lines, ) from microsim.schema.optical_config import lib -from microsim.schema.optical_config.lib import FITC +from microsim.schema.optical_config.lib import ( + FITC, + LASER_488, + LASER_561, + spectral_detector, +) from tests._util import skipif_no_internet @@ -227,8 +232,9 @@ def _confocal_sim(oc: ms.OpticalConfig, **kwargs: Any) -> ms.Simulation: def _confocal_psf(sim: ms.Simulation) -> tuple[Any, float]: modality = cast("ms.Confocal", sim.modality) - [channel] = sim._channel_fluorescence() - detection = modality._detection(channel, sim.objective_lens) + fluorescence = sim._channel_fluorescence() + [channel] = fluorescence + detection = modality._detection(fluorescence, 0, sim.objective_lens) _key, psf, brightness = modality._system_psf( channel[0], detection, @@ -393,8 +399,8 @@ def _pinhole(oc: ms.OpticalConfig, **kwargs: Any) -> PinholeDetection: objective_lens=objective, modality=ms.Confocal(pinhole_au=2, **kwargs), ) - [channel] = sim._channel_fluorescence() - det = cast("ms.Confocal", sim.modality)._detection(channel, objective) + fluorescence = sim._channel_fluorescence() + det = cast("ms.Confocal", sim.modality)._detection(fluorescence, 0, objective) return cast("PinholeDetection", det) # one pinhole for both fluorophores, 2 AU at the band center @@ -425,6 +431,44 @@ def _pinhole(oc: ms.OpticalConfig, **kwargs: Any) -> PinholeDetection: ms.Confocal(pinhole_au=closed) +def test_confocal_shared_pinhole() -> None: + """Channels in the same track share one pinhole, sized for all their bands.""" + objective = ObjectiveLens(numerical_aperture=1.4) + channels = spectral_detector(4, 500, 600, [488]) # 25 nm bins, 500-600 nm + + def _radii(tracks: Any = "simultaneous", **kwargs: Any) -> list[float]: + sim = ms.Simulation( + truth_space={"shape": (4, 4, 4), "scale": (1, 1, 1)}, + sample=[{"distribution": ms.MatsLines(), "fluorophore": FLAT}], + channels=channels, + objective_lens=objective, + modality=ms.Confocal(**kwargs), + tracks=tracks, + ) + modality = cast("ms.Confocal", sim.modality) + fluorescence = sim._channel_fluorescence() + return [ + cast( + "PinholeDetection", modality._detection(fluorescence, c, objective) + ).radius_um + for c in range(len(channels)) + ] + + per_channel = _radii(shared_pinhole=False) + assert len(set(per_channel)) == len(channels) # one pinhole per band + assert _radii("sequential") == per_channel # each channel is its own track + shared = _radii() # shared within the track by default + assert len(set(shared)) == 1 + centers = [oc.emission_band_center() for oc in channels] + expected = airy_radius_um(float(np.mean(centers)) / 1000, 1.4) # type: ignore + assert shared[0] == pytest.approx(expected, rel=1e-3) + assert shared[0] == pytest.approx(airy_radius_um(0.55, 1.4), rel=0.01) + + # one pinhole per track + two_tracks = _radii([[0, 1], [2, 3]]) + assert two_tracks[0] == two_tracks[1] < two_tracks[2] == two_tracks[3] + + def test_confocal_channel_without_signal() -> None: """A channel that detects nothing is dark (no pinhole size is needed).""" far_red = ms.OpticalConfig( @@ -512,3 +556,189 @@ def test_laser_channels(laser: int) -> None: illum = ch.illumination assert illum.wavelength[np.argmax(illum.intensity)] == laser assert ch.emission.spectrum.at([laser])[0] < 1e-6 + + +def _track_sim( + channels: list[ms.OpticalConfig], tracks: Any, **kwargs: Any +) -> ms.Simulation: + return ms.Simulation( + truth_space={"shape": (4, 4, 4), "scale": (1, 1, 1)}, + sample=[{"distribution": ms.MatsLines(), "fluorophore": FLAT}], + channels=channels, + tracks=tracks, + **kwargs, + ) + + +def test_tracks_by_index_or_name() -> None: + channels = [LASER_488, LASER_561] + assert _track_sim(channels, "sequential")._track_indices() == [[0], [1]] + assert _track_sim(channels, "simultaneous")._track_indices() == [[0, 1]] + assert _track_sim(channels, [["561"], [0]])._track_indices() == [[1], [0]] + sim = _track_sim(channels, [[0, "561"]]) + assert sim.tracks == [ms.Track(channels=[0, "561"])] # bare lists are cast + assert ms.Simulation.model_validate_json(sim.model_dump_json()).tracks == sim.tracks + + bad = [ + ([[0]], "exactly one track"), # channel 1 missing + ([[0, 0], [1]], "exactly one track"), # channel 0 twice + ([["640"], [0, 1]], "matches no"), + ([[2], [0, 1]], "out of range"), + ] + for tracks, match in bad: + with pytest.raises(ValidationError, match=match): + _track_sim(channels, tracks) + with pytest.raises(ValidationError, match="matches several"): + _track_sim([LASER_488, LASER_488], [["488"], [1]]) + + +def test_simultaneous_track_excites_with_all_lasers() -> None: + channels = [LASER_488, LASER_561] + seq = _track_sim(channels, "sequential")._channel_fluorescence() + sim = _track_sim(channels, "simultaneous")._channel_fluorescence() + both = sum(row[0].total_absorption_rate(saturate=False) for row in seq) + for c in range(2): + assert sim[c][0].total_absorption_rate(saturate=False) == pytest.approx(both) + # detection is still each channel's own + assert sim[c][0].detected_fraction.equals(seq[c][0].detected_fraction) + assert [row[0].track for row in seq] == [0, 1] + assert [row[0].track for row in sim] == [0, 0] + + +def test_track_counts_a_shared_laser_once() -> None: + # two detectors (emission windows) behind the same 488 nm laser + em2 = ms.Bandpass(bandcenter=600, bandwidth=50, placement="EM") + other_window = LASER_488.model_copy( + update={"filters": [*LASER_488.filters[:-1], em2]} + ) + channels = [LASER_488, other_window] + seq = _track_sim(channels, "sequential")._channel_fluorescence() + sim = _track_sim(channels, "simultaneous")._channel_fluorescence() + for c in range(2): + assert sim[c][0].total_absorption_rate(saturate=False) == pytest.approx( + seq[c][0].total_absorption_rate(saturate=False) + ) + # ... so the channels must agree on how the laser is set up + brighter = other_window.model_copy(update={"power": 5}) + with pytest.raises(ValidationError, match="same power and excitation filters"): + _track_sim([LASER_488, brighter], "simultaneous") + _track_sim([LASER_488, brighter], "sequential") # fine, if acquired separately + two_lasers = LASER_488.model_copy( + update={"lights": [*LASER_488.lights, *LASER_561.lights]} + ) + with pytest.raises(ValidationError, match="some, but not all"): + _track_sim([LASER_488, two_lasers], "simultaneous") + + +def test_track_shares_pixel_dwell() -> None: + slow = LASER_561.model_copy(update={"exposure_ms": 0.01}) + with pytest.raises(ValidationError, match="same pixel dwell time"): + _track_sim([LASER_488, slow], "simultaneous", modality=ms.Confocal()) + _track_sim([LASER_488, slow], "sequential", modality=ms.Confocal()) + # cameras can have different exposures + _track_sim([LASER_488, slow], "simultaneous", modality=ms.Widefield()) + + +def test_track_light_identity() -> None: + # the same laser, written differently, is still one laser + as_float = LASER_488.model_copy(update={"lights": [ms.LightSource.laser(488.0)]}) + em2 = ms.Bandpass(bandcenter=600, bandwidth=50, placement="EM") + as_float = as_float.model_copy(update={"filters": [*LASER_488.filters[:-1], em2]}) + seq = _track_sim([LASER_488, as_float], "sequential")._channel_fluorescence() + sim = _track_sim([LASER_488, as_float], "simultaneous")._channel_fluorescence() + assert sim[0][0].total_absorption_rate(saturate=False) == pytest.approx( + seq[0][0].total_absorption_rate(saturate=False) + ) + # ... and the order of a channel's lights doesn't matter + lasers = [*LASER_488.lights, *LASER_561.lights] + a = LASER_488.model_copy(update={"lights": lasers}) + b = as_float.model_copy(update={"lights": lasers[::-1]}) + _track_sim([a, b], "simultaneous") + + +def test_track_lamp_channels() -> None: + # channels without lights (a lamp through excitation filters) behave the same + em2 = ms.Bandpass(bandcenter=600, bandwidth=50, placement="EM") + other_window = FITC.model_copy(update={"filters": [*FITC.filters[:-1], em2]}) + assert other_window.excitation == FITC.excitation # same excitation side + seq = _track_sim([FITC, other_window], "sequential")._channel_fluorescence() + sim = _track_sim([FITC, other_window], "simultaneous")._channel_fluorescence() + assert sim[0][0].total_absorption_rate(saturate=False) == pytest.approx( + seq[0][0].total_absorption_rate(saturate=False) + ) + brighter = other_window.model_copy(update={"power": 5}) + with pytest.raises(ValidationError, match="same power and excitation filters"): + _track_sim([FITC, brighter], "simultaneous") + + +def test_simultaneous_confocal_adds_up_lasers() -> None: + """Without saturation, a track's image is the sum of what each laser excites.""" + sim = ms.Simulation( + truth_space={"shape": (8, 32, 32), "scale": (0.1, 0.05, 0.05)}, + sample=[{"distribution": ms.MatsLines(), "fluorophore": FLAT}], # no lifetime + channels=[LASER_488, LASER_561], + tracks="simultaneous", + modality=ms.Confocal(), + settings=ms.Settings(float_dtype="float64", cache=False), + ) + both = sim.optical_image_per_fluor() + by_laser = [ + ms.ChannelFluorescence.from_channel(oc, FLAT).absorption_rate + for oc in sim.channels + ] + + def excited_by(absorption: xr.DataArray) -> xr.DataArray: + # every channel's own detection, excited by `absorption` alone + fluorescence = [ + [ + ms.ChannelFluorescence.create( + oc, FLAT, sim.objective_lens, absorption_rate=absorption, track=0 + ) + ] + for oc in sim.channels + ] + args = (sim.objective_lens, sim.settings, ms.NumpyAPI()) + return sim.modality.render(sim.ground_truth(), fluorescence, *args) + + total = sum(excited_by(a) for a in by_laser) + np.testing.assert_allclose(both, total, rtol=1e-5) + assert not np.allclose(both, excited_by(by_laser[0])) + + +def test_channel_fluorescence_without_track_has_own_pinhole() -> None: + channels = spectral_detector(2, 500, 600, [488]) + objective = ObjectiveLens(numerical_aperture=1.4) + fluorescence = [ + [ms.ChannelFluorescence.create(oc, FLAT, objective)] for oc in channels + ] + modality = ms.Confocal() + radii = [ + cast("PinholeDetection", modality._detection(fluorescence, c, objective)) + for c in range(2) + ] + assert radii[0].radius_um < radii[1].radius_um + + +def test_rejected_assignment_is_undone() -> None: + sim = _track_sim([LASER_488, LASER_561], [[0], [1]]) + with pytest.raises(ValidationError, match="out of range"): + sim.channels = [LASER_488] + assert sim.channels == [LASER_488, LASER_561] + with pytest.raises(ValidationError, match="exactly one track"): + sim.tracks = [[0]] + assert sim.tracks == [ms.Track(channels=[0]), ms.Track(channels=[1])] + sim = ms.Simulation( + truth_space={"shape": (4, 4, 4), "scale": (1, 1, 1)}, + sample=[ms.MatsLines()], + channels=[LASER_488, LASER_561], + ) + with pytest.raises(ValidationError): + sim.tracks = [[0]] + assert "tracks" not in sim.model_fields_set + + +def test_track_validation() -> None: + channels = [LASER_488, LASER_561] + for bad in ([[0, 1], []], [[True], [0]]): + with pytest.raises(ValidationError): + _track_sim(channels, bad)