fix: confocal excitation at the excitation wavelength; redesign PSF pipeline, add acquisition tracks and laser channels

README.md

@@ -72,15 +72,28 @@ object.
 
 ```python
 from microsim import schema as ms
+from microsim.schema.detectors.lib import HYD_SP8
+from microsim.schema.optical_config.lib import LASER_488
 from microsim.util import ortho_plot
 
 # define the parameters of the simulation
 sim = ms.Simulation(
     truth_space=ms.ShapeScaleSpace(shape=(128, 512, 512), scale=(0.02, 0.01, 0.01)),
     output_space={"downscale": 8},
-    sample=ms.Sample(labels=[ms.MatsLines(density=0.5, length=30, azimuth=5, max_r=1)]),
-    modality=ms.Confocal(pinhole_au=0.2),
-    output_path="au02.tiff",
+    sample=ms.Sample(
+        labels=[
+            ms.FluorophoreDistribution(
+                distribution=ms.MatsLines(density=0.5, length=30, azimuth=5, max_r=1),
+                fluorophore="EGFP",
+            )
+        ]
+    ),
+    # 488 nm laser; power is the peak irradiance at the focus (W/cm^2)
+    channels=[LASER_488.model_copy(update={"power": 1e5})],
+    modality=ms.Confocal(pinhole_au=0.6),
+    detector=HYD_SP8,
+    exposure_ms=0.002,  # 2 µs pixel dwell time
+    output_path="au06.tiff",
 )
 
 # run it

docs/concept.md

@@ -84,11 +84,11 @@ can be used to make the simulation reproducible, and other general parameters.
 
 ### Point Spread Functions
 
-There is no top-level point spread function object in the simulation: most
-modalities will be responsible for generating their own point spread function
-based on the objective lens and other parameter that are passed into their
-`render()` method.  The `microsim.psf` module contains some common point spread
-functions, with `make_psf()` being the most common entrypoint.
+There is no top-level point spread function object in the simulation: each
+modality builds a *system PSF* for every fluorophore in every channel (the
+excitation spot times the detection PSF), based on the objective lens and on how bright each fluorophore is in each
+channel (and at which wavelengths).  The `microsim.psf` module contains
+`make_psf()` and `make_confocal_psf()`, for making PSFs directly.
 
 ### FPbase
 

docs/index.md

@@ -42,13 +42,26 @@ Construct and run a [`microsim.Simulation`][] object.
 
     ```python
     from microsim import schema as ms
+    from microsim.schema.detectors.lib import HYD_SP8
+    from microsim.schema.optical_config.lib import LASER_488
     from microsim.util import ortho_plot
 
     sim = ms.Simulation(
         truth_space=ms.ShapeScaleSpace(shape=(128, 512, 512), scale=(0.02, 0.01, 0.01)),
         output_space={"downscale": 8},
-        sample=ms.Sample(labels=[ms.MatsLines(density=0.5, length=30, azimuth=5, max_r=1)]),
-        modality=ms.Confocal(pinhole_au=0.2),
+        sample=ms.Sample(
+            labels=[
+                ms.FluorophoreDistribution(
+                    distribution=ms.MatsLines(density=0.5, length=30, azimuth=5, max_r=1),
+                    fluorophore="EGFP",
+                )
+            ]
+        ),
+        # 488 nm laser; power is the peak irradiance at the focus (W/cm^2)
+        channels=[LASER_488.model_copy(update={"power": 1e5})],
+        modality=ms.Confocal(pinhole_au=0.6),
+        detector=HYD_SP8,
+        exposure_ms=0.002,  # 2 µs pixel dwell time
     )
     result = sim.run()
     ortho_plot(result)
@@ -70,15 +83,38 @@ Construct and run a [`microsim.Simulation`][] object.
         sample=dict(
             labels=[
                 {
-                    "type": "matslines",
-                    "density": 0.5,
-                    "length": 30,
-                    "azimuth": 5,
-                    "max_r": 1,
+                    "distribution": {
+                        "type": "matslines",
+                        "density": 0.5,
+                        "length": 30,
+                        "azimuth": 5,
+                        "max_r": 1,
+                    },
+                    "fluorophore": "EGFP",
                 }
             ]
         ),
-        modality={"type": "confocal", "pinhole_au": 0.2},
+        channels=[
+            {
+                "name": "488",
+                "lights": [
+                    {"spectrum": {"wavelength": [487, 488, 489], "intensity": [0, 1, 0]}}
+                ],
+                "filters": [
+                    {"type": "longpass", "cuton": 495, "placement": "BS"},
+                    {
+                        "type": "bandpass",
+                        "bandcenter": 525,
+                        "bandwidth": 50,
+                        "placement": "EM",
+                    },
+                ],
+                "power": 1e5,
+            }
+        ],
+        modality={"type": "confocal", "pinhole_au": 0.6},
+        detector={"camera_type": "HyD", "qe": 0.45, "dark_current": 560},
+        exposure_ms=0.002,
     )
     result = sim.run()
     ortho_plot(result)
@@ -102,19 +138,37 @@ the result of the simulation to a file.
     "sample": {
         "labels": [
         {
-            "type": "matslines",
-            "density": 0.5,
-            "length": 30,
-            "azimuth": 5,
-            "max_r": 1
+            "distribution": {
+                "type": "matslines",
+                "density": 0.5,
+                "length": 30,
+                "azimuth": 5,
+                "max_r": 1
+            },
+            "fluorophore": "EGFP"
         }
         ]
     },
+    "channels": [
+        {
+            "name": "488",
+            "lights": [
+                {"spectrum": {"wavelength": [487, 488, 489], "intensity": [0, 1, 0]}}
+            ],
+            "filters": [
+                {"type": "longpass", "cuton": 495, "placement": "BS"},
+                {"type": "bandpass", "bandcenter": 525, "bandwidth": 50, "placement": "EM"}
+            ],
+            "power": 100000
+        }
+    ],
     "modality": {
         "type": "confocal",
-        "pinhole_au": 0.2
+        "pinhole_au": 0.6
     },
-    "output_path": "au02.tiff"
+    "detector": {"camera_type": "HyD", "qe": 0.45, "dark_current": 560},
+    "exposure_ms": 0.002,
+    "output_path": "au06.tiff"
 }
 ```
 

docs/stages.md

@@ -213,6 +213,47 @@ where:
 - $T_{\text{em}}(\lambda)$ is the transmission of the emission filters at each wavelength.
 - $\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
+light source `power`): 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:
@@ -345,6 +386,26 @@ which retains the full `(C, F, Z, Y, X)` dimensions.
     In microsim, the PSF largely be determined by the [`ObjectiveLens`][microsim.schema.ObjectiveLens] and the
     [`Modality`](api.md#modality) of the simulation
 
+!!! info "Widefield vs. confocal"
+
+    For every modality, the signal from a fluorophore at some position (relative
+    to the focus) is *how strongly it is excited there* times *the fraction of its
+    emission that reaches the detector*: the excitation spot times the detection
+    PSF, or the **system PSF**.
+
+    - **Widefield**: the whole field is evenly illuminated, so only the detection
+      PSF (the image of a point source on the camera) shapes the image.
+    - **Confocal**: the system PSF is the excitation spot times the detection PSF
+      through the pinhole.  The excitation spot is computed at the excitation
+      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.  As on most 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
 that each channel dimension will contain the emission of *all* fluorophores that
@@ -356,15 +417,16 @@ bin in the detector.
 !!! example
 
     This example sets up a simulation with a confocal microscope with a 1.4 NA
-    objective lens and a 1.2 AU pinhole. The simulation is set up to use the "491"
-    optical configuration from the [Example Yokogawa
-    Setup](https://www.fpbase.org/microscope/4yL4ggAo/?c=491) (spinning disc confocal)
-    microscope on FPbase.
+    objective lens, a 1.2 AU pinhole, and a 488 nm laser (`LASER_488`: the laser,
+    a dichroic, and a 525/50 emission filter).  For point-scanning modalities, the
+    laser `power` is the peak irradiance at the focus.
 
     ```python
+    from microsim.schema.optical_config.lib import LASER_488
+
     sim = Simulation(
         # ...,
-        channels=["4yL4ggAo::491"],
+        channels=[LASER_488.model_copy(update={"power": 1e5})],  # W/cm² at the focus
         objective_lens={
             "numerical_aperture": 1.4,
             "immersion_medium_ri": 1.515,

docs/tutorial.ipynb

@@ -972,7 +972,8 @@
     "how confocal imaging reduces the axial extent of the PSF and increases the\n",
     "contrast of the image.\n",
     "\n",
-    "For point-scanning modalities like `Confocal`, the illumination `power` is the\n",
+    "Point-scanning confocals illuminate with lasers and detect with point detectors\n",
+    "(PMTs or hybrid detectors) rather than cameras.  The illumination `power` is the\n",
     "*peak* irradiance at the focus of the scanned spot (typically kW/cm², much higher\n",
     "than widefield), and `exposure_ms` is the integration (dwell) time *per pixel*\n",
     "(typically a few µs):"
@@ -1013,9 +1014,15 @@
     }
    ],
    "source": [
+    "from microsim.schema.detectors.lib import HYD_SP8\n",
+    "from microsim.schema.optical_config.lib import LASER_488, LASER_561\n",
+    "\n",
     "sim.modality = ms.Confocal(pinhole_au=1)\n",
-    "sim.channels[0].power = 5_000  # W/cm² at the focus\n",
-    "sim.channels[1].power = 5_000\n",
+    "sim.channels = [\n",
+    "    LASER_488.model_copy(update={\"power\": 5_000}),  # W/cm² at the focus\n",
+    "    LASER_561.model_copy(update={\"power\": 5_000}),\n",
+    "]\n",
+    "sim.detector = HYD_SP8\n",
     "\n",
     "ortho_plot(sim.digital_image(exposure_ms=0.01), z=9)  # 10 µs pixel dwell"
    ]

examples/basic_confocal.py

@@ -1,4 +1,6 @@
 from microsim import schema as ms
+from microsim.schema.detectors.lib import HYD_SP8
+from microsim.schema.optical_config.lib import LASER_488
 from microsim.util import ortho_plot
 
 sim = ms.Simulation(
@@ -11,10 +13,12 @@
             concentration=5,
         )
     ],
-    exposure_ms=200,
+    # 488 nm laser; power is the peak irradiance at the focus (~50 µW at NA 1.4)
+    channels=[LASER_488.model_copy(update={"power": 1e5})],  # W/cm^2
     modality=ms.Confocal(pinhole_au=0.5),
-    settings=ms.Settings(random_seed=100, max_psf_radius_aus=8),
-    detector=ms.CameraCCD(qe=0.82, read_noise=2, bit_depth=12),
+    detector=HYD_SP8,
+    exposure_ms=0.002,  # 2 µs pixel dwell time
+    settings=ms.Settings(random_seed=100),
 )
 
 result = sim.run()

examples/fpbase_scope.py

@@ -11,9 +11,9 @@
             )
         ]
     ),
-    modality=ms.Confocal(pinhole_au=0.5),
+    modality=ms.Widefield(),
     channels=[ms.OpticalConfig.from_fpbase("i6WL2W", "Widefield Dual Green")],
-    settings=ms.Settings(random_seed=100, max_psf_radius_aus=8),
+    settings=ms.Settings(random_seed=100),
     detector=ms.CameraCCD(qe=0.82, read_noise=2, bit_depth=12),
     # output_path="au1.tif",
 )

examples/multi_channel_confocal.py

@@ -1,11 +1,14 @@
 import time
 
 from microsim import schema as ms
-from microsim.schema.optical_config.lib import EYFP, FITC
+from microsim.schema.detectors.lib import HYD_SP8
+from microsim.schema.optical_config.lib import LASER_488, LASER_515
 from microsim.util import ndview
 
-# short stubby lines labeled in EGFP and long thin lines labeled in EYFP
-# bleedthrough of the EGFP signal in to the EYFP is more visible than vice versa
+# short stubby lines labeled in EGFP and long thin lines labeled in EYFP, imaged
+# with 488 and 515 nm lasers at the same time (one track).  With both lasers on,
+# each fluorophore bleeds into the other's channel.  Try `tracks="sequential"`:
+# 515 nm barely excites EGFP, so its bleedthrough into the EYFP channel vanishes
 
 sim = ms.Simulation(
     truth_space=ms.ShapeScaleSpace(shape=(64, 256, 256), scale=(0.04, 0.02, 0.02)),
@@ -23,9 +26,15 @@
         ]
     ),
     modality=ms.Confocal(pinhole_au=1),
-    settings=ms.Settings(random_seed=100, max_psf_radius_aus=4),
-    detector=ms.CameraCCD(qe=0.82, read_noise=4, bit_depth=12),
-    channels=[FITC, EYFP],
+    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
+    # power is the peak irradiance at the focus (W/cm^2)
+    channels=[
+        LASER_488.model_copy(update={"power": 1e5}),
+        LASER_515.model_copy(update={"power": 1e5}),
+    ],
     # output_path="au1.tif",
 )