fix: off-grid laser lines, dropped emission sample, and QE without an emission filter

src/microsim/schema/modality/_simple_psf.py

@@ -348,7 +348,7 @@ def bin_spectrum(
         bins = np.linspace(w_min, w_max, num_bins + 1)
 
     # Use groupby_bins to bin the data within the filtered region
-    binned = masked.groupby_bins(Axis.W, bins=bins)
+    binned = masked.groupby_bins(Axis.W, bins=bins, include_lowest=True)
 
     # Create a new DataArray with the summed intensities and centroid wavelengths
     binned_spectrum = binned.sum(Axis.W)

src/microsim/schema/optical_config/config.py

@@ -71,7 +71,10 @@ def absorption_rate(self, fluorophore: Fluorophore) -> xr.DataArray:
         and singleton dimensions F and C.
         """
         illum_flux_density = self.illumination_flux_density  # photons/cm^2/s
-        cross_section = fluorophore.absorption_cross_section  # cm^2/fluorophore
+        # cm^2/fluorophore, at each wavelength of the illumination
+        cross_section = fluorophore.absorption_cross_section.interp(
+            {Axis.W: illum_flux_density.coords[Axis.W]}, kwargs={"fill_value": 0}
+        )
 
         # absorption rate in photons/s/fluorophore
         abs_rate = illum_flux_density * cross_section
@@ -120,18 +123,40 @@ def total_emission_rate(
             # - pulsed excitation: assumes CW (steady state), not pulses << lifetime
             s = self.saturation_parameter(fluorophore)
             tot_absorption_rate = tot_absorption_rate / (1 + s)
-        em_rate = fluorophore.emission_spectrum.as_xarray()
-        # norm area to 1
-        em_rate = em_rate / em_rate.sum()
-        # multiply by quantum yield and total absorption rate
-        em_rate = em_rate * (fluorophore.quantum_yield or 1) * tot_absorption_rate
+        em_rate = fluorophore.emission_per_absorbed_photon * tot_absorption_rate
         # add singleton coordinates for fluorophore and self
         em_rate = em_rate.expand_dims({Axis.F: [fluorophore], Axis.C: [self]})
         em_rate.name = "emission_rate"
         em_rate.attrs["long_name"] = "Emission rate"
         em_rate.attrs["units"] = "photons/s"
         return em_rate
 
+    def detected_fraction(
+        self, fluorophore: Fluorophore, detector_qe: float | Spectrum | None = None
+    ) -> xr.DataArray:
+        """Return the fraction of absorbed photons that are detected in this channel.
+
+        This is the number of photons detected per photon absorbed by `fluorophore`,
+        as a function of emission wavelength: the quantum yield, times the emission
+        spectrum (normalized to an area of 1), times the transmission of the emission
+        filters and the detector QE.  It does *not* include the fraction of emitted
+        light collected by the objective (see `ObjectiveLens.collection_efficiency`).
+        """
+        if detector_qe is None and self.detector is not None:
+            detector_qe = self.detector.qe
+        detected = fluorophore.emission_per_absorbed_photon
+        wavelength = detected.coords[Axis.W]
+        if self.emission is not None:
+            detected = detected * self.emission.spectrum.at(wavelength)
+        if isinstance(detector_qe, Spectrum):
+            detected = detected * detector_qe.at(wavelength)
+        elif detector_qe is not None:
+            detected = detected * detector_qe
+        detected.name = "detected_fraction"
+        detected.attrs["long_name"] = "Detected photons per absorbed photon"
+        detected.attrs["units"] = ""
+        return detected
+
     def filtered_emission_rate(
         self,
         fluorophore: Fluorophore,
@@ -156,17 +181,13 @@ def filtered_emission_rate(
             - the emission filter/beamsplitter spectra
             - camera QE, if passed
         """
-        if self.emission is None:
-            em_array: xr.DataArray | float = 1.0
-        else:
-            em_spectrum = self.emission.spectrum
-            if detector_qe is None and self.detector is not None:
-                detector_qe = self.detector.qe
-            if detector_qe is not None:
-                em_spectrum = em_spectrum * detector_qe
-            em_array = em_spectrum.as_xarray()
-
-        final = self.total_emission_rate(fluorophore, saturate=saturate) * em_array
+        absorption_rate = self.absorption_rate(fluorophore).sum()
+        if saturate:
+            absorption_rate /= 1 + self.saturation_parameter(fluorophore)
+        detected = self.detected_fraction(fluorophore, detector_qe)
+        final = absorption_rate * detected
+        # add singleton coordinates for fluorophore and self
+        final = final.expand_dims({Axis.F: [fluorophore], Axis.C: [self]})
         final.name = "filtered_emission_rate"
         final.attrs["long_name"] = "Filtered Emission rate"
         final.attrs["units"] = "photons/s"
@@ -200,7 +221,7 @@ def illumination(self) -> Spectrum | None:
                 for light in rest:
                     illum_spect = illum_spect + light.scaled_spectrum()
             if exc:
-                return illum_spect * exc.spectrum
+                return illum_spect.transmitted(exc.spectrum)
             return illum_spect
         return exc.spectrum if exc else None
 

src/microsim/schema/sample/fluorophore.py

@@ -37,6 +37,15 @@ def absorption_cross_section(self) -> xr.DataArray:
         out.name = "cross_section"
         return cast("xr.DataArray", out)
 
+    @property
+    def emission_per_absorbed_photon(self) -> xr.DataArray:
+        """Photons emitted per absorbed photon, at each wavelength.
+
+        The emission spectrum (normalized to an area of 1) times the quantum yield.
+        """
+        em = self.emission_spectrum.as_xarray()
+        return em / em.sum() * (self.quantum_yield or 1)
+
     @classmethod
     def from_fpbase(cls, name: str) -> "Fluorophore":
         from fpbase import get_fluorophore

src/microsim/schema/spectrum.py

@@ -13,6 +13,7 @@
 
 if TYPE_CHECKING:
     from fpbase.models import Spectrum as FPbaseSpectrum
+    from numpy.typing import ArrayLike
     from xarray.plot.accessor import DataArrayPlotAccessor
 
 
@@ -122,7 +123,42 @@ def __truediv__(self, other: "float | pint.Quantity | Spectrum") -> "Spectrum":
             return self._intensity_op(other, np.true_divide)
         return self.model_copy(update={"intensity": self.intensity / other})
 
+    def at(self, wavelength: "ArrayLike") -> np.ndarray:
+        """Return the intensity at `wavelength` (nm).
+
+        Values between samples are linearly interpolated; outside the range of this
+        spectrum, the intensity is 0.
+        """
+        order = np.argsort(self.wavelength)
+        return np.interp(  # type: ignore [no-any-return]
+            np.asarray(wavelength, dtype=float),
+            self.wavelength[order],
+            self.intensity[order],
+            left=0,
+            right=0,
+        )
+
+    def transmitted(self, transmission: "Spectrum") -> "Spectrum":
+        """Return this (light) spectrum after passing through `transmission`.
+
+        Unlike `self * transmission`, the result keeps this spectrum's wavelengths
+        (e.g. a laser line stays a laser line), and `transmission` is interpolated.
+        """
+        return self.model_copy(
+            update={"intensity": self.intensity * transmission.at(self.wavelength)}
+        )
+
     def _intensity_op(self, other: "Spectrum", op: np.ufunc) -> "Spectrum":
+        if op is np.multiply:
+            # e.g. combining filters, or a filter and detector QE: evaluate both
+            # at every wavelength either one is sampled at.
+            wavelength = np.union1d(self.wavelength, other.wavelength)
+            return self.model_copy(
+                update={
+                    "intensity": self.at(wavelength) * other.at(wavelength),
+                    "wavelength": wavelength,
+                }
+            )
         selfx, otherx = xr.align(self.as_xarray(), other.as_xarray(), join="outer")
         selfx = selfx.fillna(0)
         otherx = otherx.fillna(0)

tests/test_photometry.py

@@ -2,11 +2,14 @@
 
 import numpy as np
 import pytest
+import xarray as xr
 from pydantic import ValidationError
 
 from microsim import schema as ms
 from microsim.psf import make_confocal_psf
+from microsim.schema.dimensions import Axis
 from microsim.schema.lens import ObjectiveLens
+from microsim.schema.modality._simple_psf import bin_spectrum
 from microsim.schema.optical_config.lib import FITC
 from tests._util import skipif_no_internet
 
@@ -152,6 +155,47 @@ def test_confocal_saturation_in_psf_not_rates() -> None:
     assert wf.modality._saturation_parameter(em_spectrum) == 0
 
 
+def test_off_grid_laser_wavelength() -> None:
+    """Lasers between the sampled wavelengths of filters/spectra still excite."""
+    fluor = _green_fluor()
+    rates = [
+        ms.OpticalConfig(lights=[ms.LightSource.laser(w)], filters=FITC.filters)
+        .absorption_rate(fluor)
+        .sum()
+        .item()
+        for w in (488, 488.5)
+    ]
+    assert rates[0] > 0
+    assert rates[1] == pytest.approx(rates[0], rel=1e-3)
+
+
+def test_spectrum_product_on_different_grids() -> None:
+    """Samples present in only one spectrum are interpolated, not zeroed."""
+    wvl = np.arange(500, 600)
+    fine = ms.Spectrum(wavelength=wvl, intensity=np.ones(wvl.size))
+    for other in (wvl[::2], wvl + 0.5):  # coarser, or offset by half a sample
+        half = ms.Spectrum(wavelength=other, intensity=np.full(other.size, 0.5))
+        product = fine * half
+        assert isinstance(product.intensity, np.ndarray)
+        inside = np.arange(501, 598)
+        np.testing.assert_allclose(product.at(inside), 0.5)
+
+
+def test_detector_qe_without_emission_filter() -> None:
+    oc = ms.OpticalConfig(lights=[ms.LightSource.laser(488)])
+    full = oc.filtered_emission_rate(_green_fluor()).sum().item()
+    half = oc.filtered_emission_rate(_green_fluor(), detector_qe=0.5).sum().item()
+    assert half == pytest.approx(full / 2)
+
+
+def test_bin_spectrum_keeps_all_samples() -> None:
+    spectrum = xr.DataArray(
+        [1.0, 2, 3, 4], dims=[Axis.W], coords={Axis.W: [500, 501, 502, 503]}
+    )
+    for bins in (1, 2, 3):
+        assert float(bin_spectrum(spectrum, bins=bins).sum()) == 10
+
+
 # flat spectra across the visible, so the fluorophore is excited and detected in FITC
 _WAVES = np.arange(350, 750)
 FLAT = ms.Fluorophore(