Add photons_absorbed, the inverse of the sensor signal/expose

optika/sensors/_sensors.py

@@ -235,6 +235,69 @@ def expose(
 
         return dataclasses.replace(image, outputs=electrons)
 
+    def photons_absorbed(
+        self,
+        image: na.FunctionArray[
+            na.SpectralPositionalVectorArray,
+            na.AbstractScalar,
+        ],
+        direction: float | na.AbstractScalar = 1,
+        axis_wavelength: None | str = None,
+        timedelta: None | u.Quantity | na.AbstractScalar = None,
+    ) -> na.FunctionArray[
+        na.SpectralPositionalVectorArray,
+        na.AbstractScalar,
+    ]:
+        """
+        Invert :meth:`expose`, mapping the electrons measured in each pixel back
+        into a photon flux absorbed by the light-sensitive region.
+
+        The absorbance is *not* restored, since :meth:`expose` runs the
+        detector with an absorbance of one (the absorbance is usually accounted
+        for elsewhere, such as in the effective area of an optical system), so
+        this only divides out the quantum yield, the charge collection
+        efficiency, and the exposure time. It is the deterministic inverse of
+        :meth:`expose`; the sensor noise is not undone.
+
+        Parameters
+        ----------
+        image
+            The electrons measured in each pixel, as a function of wavelength
+            and pixel position.
+            The wavelength inputs (``image.inputs.wavelength``) must be the
+            bin *edges*, not the centers.
+        direction
+            The cosine of the refracted angle inside the light-sensitive region,
+            matching the value passed to :meth:`expose`.
+        axis_wavelength
+            The logical axis of `image` corresponding to changing wavelength.
+            If :obj:`None` (the default), ``image.inputs.wavelength`` must have
+            only one logical axis.
+        timedelta
+            The exposure time of the measurement.
+            If :obj:`None` (the default), the value in :attr:`timedelta_exposure`
+            will be used.
+        """
+        if axis_wavelength is None:
+            shape_wavelength = na.shape(image.inputs.wavelength)
+            if len(shape_wavelength) != 1:  # pragma: nocover
+                raise ValueError(
+                    f"if `axis_wavelength` is `None`, `image.inputs.wavelength` "
+                    f"must have exactly one logical axis, got {shape_wavelength}."
+                )
+            (axis_wavelength,) = shape_wavelength
+
+        if timedelta is None:
+            timedelta = self.timedelta_exposure
+
+        photons = self.material.photons_absorbed(
+            electrons=image.outputs,
+            wavelength=image.inputs.wavelength.cell_centers(axis_wavelength),
+            direction=direction,
+        )
+
+        return dataclasses.replace(image, outputs=photons / timedelta)
+
     def measure(
         self,
         rays: optika.rays.RayVectorArray,

optika/sensors/_sensors_test.py

@@ -1,4 +1,5 @@
 import pytest
+import dataclasses
 import numpy as np
 import astropy.units as u
 import named_arrays as na
@@ -88,6 +89,46 @@ def test_measure(
         assert a.axis_pixel.x in result_lines.outputs.shape
         assert a.axis_pixel.y in result_lines.outputs.shape
 
+    def test_photons_absorbed(self, a: optika.sensors.AbstractImagingSensor):
+        # use a nonzero exposure time so the default `timedelta` is invertible
+        a = dataclasses.replace(a, timedelta_exposure=10 * u.s)
+
+        # a photon rate incident on a few pixels, as a function of the
+        # wavelength bin edges
+        wavelength = na.linspace(500, 600, axis="wavelength", num=4) * u.nm
+        position = na.Cartesian2dVectorArray(
+            x=na.arange(0, 5, axis=a.axis_pixel.x) * u.pix,
+            y=na.arange(0, 5, axis=a.axis_pixel.y) * u.pix,
+        )
+        rate = (
+            na.random.uniform(
+                low=0,
+                high=100,
+                shape_random={"wavelength": 3, a.axis_pixel.x: 5, a.axis_pixel.y: 5},
+            )
+            * u.photon
+            / u.s
+        )
+        image = na.FunctionArray(
+            inputs=na.SpectralPositionalVectorArray(
+                wavelength=wavelength,
+                position=position,
+            ),
+            outputs=rate,
+        )
+
+        # `photons_absorbed` is the deterministic inverse of `expose`
+        electrons = a.expose(image, noise=False)
+        result = a.photons_absorbed(electrons)
+
+        assert isinstance(result, na.FunctionArray)
+        assert isinstance(result.inputs, na.SpectralPositionalVectorArray)
+        assert result.outputs.unit.is_equivalent(u.photon / u.s)
+        assert np.allclose(
+            result.outputs.to_value(u.photon / u.s),
+            rate.to_value(u.photon / u.s),
+        )
+
 
 @pytest.mark.parametrize(
     argnames="a",

optika/sensors/materials/_materials.py

@@ -1513,6 +1513,33 @@ def photons_incident(
             The vector perpendicular to the surface of the sensor.
         """
 
+    @abc.abstractmethod
+    def photons_absorbed(
+        self,
+        electrons: u.Quantity | na.AbstractScalar,
+        wavelength: u.Quantity | na.AbstractScalar,
+        direction: float | na.AbstractScalar = 1,
+    ) -> na.AbstractScalar:
+        """
+        Given the number of electrons measured by the sensor, compute the
+        expected number of photons *absorbed* by the light-sensitive region.
+
+        This is the inverse of :meth:`signal`: it divides out only the quantum
+        yield and the charge collection efficiency, not the absorbance. The
+        absorbance is deliberately excluded because it is usually accounted for
+        elsewhere (for example in the effective area of the optical system).
+
+        Parameters
+        ----------
+        electrons
+            The number of electrons measured by each pixel.
+        wavelength
+            The vacuum wavelength of the absorbed photons.
+        direction
+            The cosine of the refracted angle inside the light-sensitive region,
+            as produced by :meth:`direction_refracted`.
+        """
+
 
 @dataclasses.dataclass(eq=False, repr=False)
 class IdealSensorMaterial(
@@ -1575,6 +1602,14 @@ def photons_incident(
     ) -> na.AbstractScalar:
         return electrons * u.photon / u.electron
 
+    def photons_absorbed(
+        self,
+        electrons: u.Quantity | na.AbstractScalar,
+        wavelength: u.Quantity | na.AbstractScalar,
+        direction: float | na.AbstractScalar = 1,
+    ) -> na.AbstractScalar:
+        return electrons * u.photon / u.electron
+
 
 @dataclasses.dataclass(eq=False, repr=False)
 class AbstractSiliconSensorMaterial(
@@ -2106,6 +2141,37 @@ def photons_incident(
 
         return electrons / qe
 
+    def photons_absorbed(
+        self,
+        electrons: u.Quantity | na.AbstractScalar,
+        wavelength: u.Quantity | na.AbstractScalar,
+        direction: float | na.AbstractScalar = 1,
+    ) -> na.AbstractScalar:
+        # `direction` is the cosine of the refracted angle *inside* the
+        # substrate (as passed to `signal`), so compute the substrate
+        # absorption directly from it and divide out only the quantum yield and
+        # charge collection efficiency, matching `signal` at ``absorbance=1``.
+        n_substrate = self._chemical.n(wavelength)
+
+        absorption = absorption_effective(
+            wavelength=wavelength,
+            n_substrate=n_substrate,
+            direction_substrate=direction,
+        )
+
+        iqy = quantum_yield_ideal(
+            wavelength=wavelength,
+            temperature=self.temperature,
+        )
+
+        cce = charge_collection_efficiency(
+            absorption=absorption,
+            thickness_implant=self.thickness_implant,
+            cce_backsurface=self.cce_backsurface,
+        )
+
+        return electrons / (iqy * cce)
+
     def efficiency(
         self,
         rays: optika.rays.AbstractRayVectorArray,

optika/sensors/materials/_materials_test.py

@@ -417,6 +417,54 @@ def test_photons_incident(
         assert isinstance(na.as_named_array(result), na.AbstractScalar)
         assert result.unit.is_equivalent(u.photon)
 
+    @pytest.mark.parametrize(
+        argnames="photons",
+        argvalues=[
+            100 * u.photon,
+        ],
+    )
+    @pytest.mark.parametrize(
+        argnames="wavelength",
+        argvalues=[
+            100 * u.AA,
+        ],
+    )
+    @pytest.mark.parametrize(
+        argnames="direction",
+        argvalues=[
+            1,
+            0.5,
+        ],
+    )
+    def test_photons_absorbed(
+        self,
+        a: optika.sensors.materials.AbstractSensorMaterial,
+        photons: u.Quantity | na.AbstractScalar,
+        wavelength: u.Quantity | na.AbstractScalar,
+        direction: float | na.AbstractScalar,
+    ):
+        # `photons_absorbed` inverts the noiseless `signal` (which uses unit
+        # absorbance), recovering the number of absorbed photons.
+        electrons = a.signal(
+            photons=photons,
+            wavelength=wavelength,
+            direction=direction,
+            noise=False,
+        )
+        result = a.photons_absorbed(
+            electrons=electrons,
+            wavelength=wavelength,
+            direction=direction,
+        )
+        assert isinstance(na.as_named_array(result), na.AbstractScalar)
+        assert result.unit.is_equivalent(u.photon)
+        assert np.allclose(
+            na.as_named_array(result / photons).ndarray.to_value(
+                u.dimensionless_unscaled
+            ),
+            1,
+        )
+
     @pytest.mark.parametrize(
         argnames="wavelength",
         argvalues=[