diff --git a/docs/index.rst b/docs/index.rst index 49caa701..4a948282 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -33,6 +33,8 @@ Features * Sequential raytrace modeling of an optical system * Stratified random sampling of input rays for faster convergence * Image simulation of a given scene using an optical system +* Fast linear forward model approximating a raytraced system, for imaging many + scenes without raytracing each one * Spherical, conical, and toroidal surface sag profiles * Circular, rectangular, and polygonal apertures * Support for mirrors and arbitrary multilayer coatings diff --git a/optika/distortion/_distortion.py b/optika/distortion/_distortion.py index 49a7fcaf..3116c4e9 100644 --- a/optika/distortion/_distortion.py +++ b/optika/distortion/_distortion.py @@ -24,6 +24,7 @@ @dataclasses.dataclass(eq=False, repr=False) class AbstractDistortionModel( optika.mixins.Printable, + optika.mixins.Shaped, ): """ An interface describing an arbitrary distortion model, @@ -189,6 +190,15 @@ class SimpleDistortionModel( """The reference wavelength and the sensor position that the field center maps to at that wavelength.""" + @property + def shape(self) -> dict[str, int]: + return na.broadcast_shapes( + optika.shape(self.plate_scale), + optika.shape(self.dispersion), + optika.shape(self.angle), + optika.shape(self.reference), + ) + @functools.cached_property def matrix(self) -> na.SpectralPositionalMatrixArray: cos = np.cos(self.angle) @@ -344,6 +354,14 @@ class PolynomialDistortionModel( where: bool | na.AbstractScalar = True """A boolean mask selecting which calibration points to use for fitting.""" + @property + def shape(self) -> dict[str, int]: + shape = na.broadcast_shapes( + optika.shape(self.coordinates_scene), + optika.shape(self.coordinates_sensor), + ) + return {ax: n for ax, n in shape.items() if ax not in self._axis_scene} + @property def _axis_scene(self) -> tuple[str, ...]: """The logical axes over which the calibration points are distributed.""" @@ -358,6 +376,7 @@ def fit(self) -> na.PolynomialFitFunctionArray: outputs=self.coordinates_sensor, center=scene.mean(self._axis_scene), degree=self.degree, + axis_polynomial=self._axis_scene, where_polynomial=self.where, ) @@ -374,6 +393,7 @@ def fit_inverse(self) -> na.PolynomialFitFunctionArray: outputs=scene.position, center=inputs.mean(self._axis_scene), degree=self.degree, + axis_polynomial=self._axis_scene, where_polynomial=self.where, ) diff --git a/optika/distortion/_distortion_test.py b/optika/distortion/_distortion_test.py index f10185ed..182656d2 100644 --- a/optika/distortion/_distortion_test.py +++ b/optika/distortion/_distortion_test.py @@ -22,6 +22,7 @@ def _scene() -> na.SpectralPositionalVectorArray: class AbstractTestAbstractDistortionModel( test_mixins.AbstractTestPrintable, + test_mixins.AbstractTestShaped, ): def test_distort(self, a: optika.distortion.AbstractDistortionModel): coordinates = _scene() @@ -151,3 +152,43 @@ def test_plot_residual( assert isinstance(ax, na.ScalarArray) assert a.axis_wavelength in na.shape(ax) plt.close(fig) + + +def test_polynomial_distortion_model_channel(): + """ + Calibration points that vary along an axis orthogonal to the scene axes + (e.g. the channel axis of a multi-channel instrument) must be fit with an + independent polynomial per channel, not one polynomial averaged over all + the channels. + """ + scene = _scene() + + scale = na.ScalarArray([10, 12, 8] * u.mm / u.deg, axes="channel") + angle = na.ScalarArray([0, 10, -15] * u.deg, axes="channel") + + cos, sin = np.cos(angle), np.sin(angle) + sensor = na.Cartesian2dVectorArray( + x=scale * (cos * scene.position.x - sin * scene.position.y), + y=scale * (sin * scene.position.x + cos * scene.position.y), + ) + + a = optika.distortion.PolynomialDistortionModel( + coordinates_scene=scene, + coordinates_sensor=sensor, + axis_wavelength="wavelength", + axis_field=("field_x", "field_y"), + degree=1, + ) + + distorted = a.distort(scene).position + assert "channel" in distorted.shape + assert np.all((distorted - sensor).length < 1e-9 * u.mm) + + undistorted = a.undistort( + na.SpectralPositionalVectorArray( + wavelength=scene.wavelength, + position=sensor, + ) + ).position + assert "channel" in undistorted.shape + assert np.all((undistorted - scene.position).length < 1e-9 * u.deg) diff --git a/optika/radiometry/_effective_area.py b/optika/radiometry/_effective_area.py index ab1f6117..31b618b2 100644 --- a/optika/radiometry/_effective_area.py +++ b/optika/radiometry/_effective_area.py @@ -12,6 +12,7 @@ @dataclasses.dataclass(eq=False, repr=False) class AbstractEffectiveAreaModel( optika.mixins.Printable, + optika.mixins.Shaped, ): """ An interface describing the effective area of an optical system as a @@ -90,6 +91,14 @@ class InterpolatedEffectiveAreaModel( axis_wavelength: str = dataclasses.MISSING """The logical axis corresponding to changing wavelength.""" + @property + def shape(self) -> dict[str, int]: + shape = na.broadcast_shapes( + optika.shape(self.wavelength), + optika.shape(self.area), + ) + return {ax: n for ax, n in shape.items() if ax != self.axis_wavelength} + def __call__( self, wavelength: na.AbstractScalar, diff --git a/optika/radiometry/_effective_area_test.py b/optika/radiometry/_effective_area_test.py index 8d8144d3..7e41ac66 100644 --- a/optika/radiometry/_effective_area_test.py +++ b/optika/radiometry/_effective_area_test.py @@ -16,6 +16,7 @@ def _area() -> na.AbstractScalar: class AbstractTestAbstractEffectiveAreaModel( test_mixins.AbstractTestPrintable, + test_mixins.AbstractTestShaped, ): def test__call__(self, a: optika.radiometry.AbstractEffectiveAreaModel): wavelength = na.linspace(200, 900, axis="wavelength", num=5) * u.AA diff --git a/optika/radiometry/_vignetting.py b/optika/radiometry/_vignetting.py index 33221af8..2098ce2a 100644 --- a/optika/radiometry/_vignetting.py +++ b/optika/radiometry/_vignetting.py @@ -20,6 +20,7 @@ @dataclasses.dataclass(eq=False, repr=False) class AbstractVignettingModel( optika.mixins.Printable, + optika.mixins.Shaped, ): """ An interface describing an arbitrary vignetting model, which maps scene @@ -161,6 +162,14 @@ class PolynomialVignettingModel( where: bool | na.AbstractScalar = True """A boolean mask selecting which calibration points to use for fitting.""" + @property + def shape(self) -> dict[str, int]: + shape = na.broadcast_shapes( + optika.shape(self.coordinates_scene), + optika.shape(self.illumination), + ) + return {ax: n for ax, n in shape.items() if ax not in self._axis_scene} + @property def _axis_scene(self) -> tuple[str, ...]: """The logical axes over which the calibration points are distributed.""" @@ -175,6 +184,7 @@ def fit(self) -> na.PolynomialFitFunctionArray: outputs=self.illumination, center=scene.mean(self._axis_scene), degree=self.degree, + axis_polynomial=self._axis_scene, where_polynomial=self.where, ) diff --git a/optika/radiometry/_vignetting_test.py b/optika/radiometry/_vignetting_test.py index cbef3679..43f7a9db 100644 --- a/optika/radiometry/_vignetting_test.py +++ b/optika/radiometry/_vignetting_test.py @@ -26,6 +26,7 @@ def _illumination() -> na.AbstractScalar: class AbstractTestAbstractVignettingModel( test_mixins.AbstractTestPrintable, + test_mixins.AbstractTestShaped, ): def test__call__(self, a: optika.radiometry.AbstractVignettingModel): scene = _scene() @@ -125,3 +126,28 @@ def test_plot_residual( assert isinstance(ax, na.ScalarArray) assert a.axis_wavelength in na.shape(ax) plt.close(fig) + + +def test_polynomial_vignetting_model_channel(): + """ + Calibration points that vary along an axis orthogonal to the scene axes + (e.g. the channel axis of a multi-channel instrument) must be fit with an + independent polynomial per channel, not one polynomial averaged over all + the channels. + """ + scene = _scene() + + coefficient = na.ScalarArray([0.1, 0.2, 0.05] / u.deg**2, axes="channel") + illumination = 1 - coefficient * scene.position.length**2 + + a = optika.radiometry.PolynomialVignettingModel( + coordinates_scene=scene, + illumination=illumination, + axis_wavelength="wavelength", + axis_field=("field_x", "field_y"), + degree=2, + ) + + result = a(scene) + assert "channel" in result.shape + assert np.all(np.abs(result - illumination) < 1e-9) diff --git a/optika/rays/_ray_vectors.py b/optika/rays/_ray_vectors.py index 49fb0bfc..b165fcc0 100644 --- a/optika/rays/_ray_vectors.py +++ b/optika/rays/_ray_vectors.py @@ -1,5 +1,6 @@ from __future__ import annotations from typing import TypeVar, Generic +from collections.abc import Sequence import abc import dataclasses import numpy as np @@ -98,6 +99,13 @@ def type_explicit(self) -> type[na.AbstractExplicitArray]: def type_matrix(self) -> type[na.AbstractMatrixArray]: raise NotImplementedError + def volume_cell(self, axis: None | str | Sequence[str]) -> na.AbstractScalar: + """ + A ray bundle is a scattered collection of rays rather than a + logically-rectangular grid, so the per-voxel volume is undefined. + """ + raise NotImplementedError + @property def explicit(self) -> RayVectorArray: return super().explicit diff --git a/optika/sensors/_sensors.py b/optika/sensors/_sensors.py index 15ba6501..9b7000e8 100644 --- a/optika/sensors/_sensors.py +++ b/optika/sensors/_sensors.py @@ -89,6 +89,27 @@ def aperture(self): half_width=self.width_pixel * self.num_pixel / 2, ) + def pixels( + self, + position: na.AbstractCartesian2dVectorArray, + ) -> na.AbstractCartesian2dVectorArray: + """ + Convert an in-plane position on the sensor plane into fractional pixel + coordinates. + + Pixel ``0`` is the lower edge of the light-sensitive area (the same grid + that :meth:`collect` bins onto), so an integer value lands on a pixel + boundary and a half-integer on a pixel center. + + Parameters + ---------- + position + The in-plane (``xy``) position on the sensor plane, in physical + units. + """ + lower = self.aperture.bound_lower.xy + return (position - lower) / self.width_pixel * u.pix << u.pix + def collect( self, rays: optika.rays.RayVectorArray, @@ -170,6 +191,43 @@ def collect( return image, direction + @staticmethod + def _collapse_wavelength( + inputs: na.SpectralPositionalVectorArray, + axis_wavelength: str, + ) -> na.SpectralPositionalVectorArray: + """Collapse a wavelength axis to its two band edges.""" + wavelength = inputs.wavelength + return inputs.replace( + wavelength=na.stack( + arrays=[ + wavelength.min(axis_wavelength), + wavelength.max(axis_wavelength), + ], + axis=axis_wavelength, + ) + ) + + def _integrate( + self, + image: na.FunctionArray[ + na.SpectralPositionalVectorArray, + na.AbstractScalar, + ], + axis_wavelength: str, + noise: bool, + ) -> na.FunctionArray[ + na.SpectralPositionalVectorArray, + na.AbstractScalar, + ]: + """Sum electrons over wavelength into one readout, adding read noise.""" + electrons = image.outputs.sum(axis_wavelength) + if noise: + # add zero-mean Gaussian read noise once per readout + electrons = na.random.normal(loc=electrons, scale=self.read_noise) + inputs = self._collapse_wavelength(image.inputs, axis_wavelength) + return dataclasses.replace(image, inputs=inputs, outputs=electrons) + def expose( self, image: na.FunctionArray[ @@ -180,6 +238,8 @@ def expose( axis_wavelength: None | str = None, timedelta: None | u.Quantity | na.AbstractScalar = None, noise: bool = True, + integrate: bool = True, + uncertainty: bool = False, ) -> na.FunctionArray[ na.SpectralPositionalVectorArray, na.AbstractScalar, @@ -195,7 +255,10 @@ def expose( The photon flux is multiplied by the exposure time and converted to electrons using :meth:`~optika.sensors.materials.AbstractSensorMaterial.signal`, which - applies the quantum efficiency, noise, and charge-diffusion models. + applies the quantum efficiency and the shot, Fano, and charge-diffusion + noise per wavelength. If `integrate` is :obj:`True`, the electrons are + then summed over wavelength into a single readout and the sensor's + :attr:`read_noise` is added once. Parameters ---------- @@ -216,8 +279,21 @@ def expose( If :obj:`None` (the default), the value in :attr:`timedelta_exposure` will be used. noise - Whether to add shot noise, intrinsic sensor noise, and read noise - to the result. + Whether to add shot, Fano, and charge-diffusion noise per wavelength + and (if `integrate`) read noise once per readout. + integrate + Whether to integrate the electrons over wavelength into a single + readout, applying :attr:`read_noise` once. + A real imaging sensor cannot resolve the individual wavelengths, so + this defaults to :obj:`True`; :obj:`False` keeps the wavelengths + separate for demonstration. + uncertainty + Whether to attach the standard deviation of the measurement noise to + the result as a + :class:`~named_arrays.NormalUncertainScalarArray`, computed with + :meth:`uncertainty`. + The width uses the *expected* electrons, so the noiseless signal + model is only re-evaluated when `noise` is also :obj:`True`. """ if axis_wavelength is None: shape_wavelength = na.shape(image.inputs.wavelength) @@ -232,24 +308,44 @@ def expose( timedelta = self.timedelta_exposure photons = image.outputs * timedelta + wavelength = image.inputs.wavelength.cell_centers(axis_wavelength) + + def signal(noise: bool) -> na.AbstractScalar: + return self.material.signal( + photons=photons, + wavelength=wavelength, + direction=direction, + width_pixel=self.width_pixel, + axis_xy=(self.axis_pixel.x, self.axis_pixel.y), + noise=noise, + ) - electrons = self.material.signal( - photons=photons, - wavelength=image.inputs.wavelength.cell_centers(axis_wavelength), - direction=direction, - width_pixel=self.width_pixel, - axis_xy=(self.axis_pixel.x, self.axis_pixel.y), - noise=noise, - ) + result = dataclasses.replace(image, outputs=signal(noise)) + + if uncertainty: + # the width uses the expected electrons, which are the same as the + # nominal result unless `noise` added a realization + expected = result if not noise else image.replace(outputs=signal(False)) + width = self.uncertainty( + expected, + direction=direction, + axis_wavelength=axis_wavelength, + integrate=integrate, + ) - if noise: - # add zero-mean Gaussian read noise to each pixel - electrons = na.random.normal( - loc=electrons, - scale=self.read_noise, + if integrate: + # a real sensor reads out all wavelengths at once + result = self._integrate(result, axis_wavelength, noise) + + if uncertainty: + result = result.replace( + outputs=na.NormalUncertainScalarArray( + nominal=result.outputs, + width=width.outputs, + ), ) - return dataclasses.replace(image, outputs=electrons) + return result def photons_absorbed( self, @@ -260,6 +356,7 @@ def photons_absorbed( direction: float | na.AbstractScalar = 1, axis_wavelength: None | str = None, timedelta: None | u.Quantity | na.AbstractScalar = None, + integrate: bool = True, ) -> na.FunctionArray[ na.SpectralPositionalVectorArray, na.AbstractScalar, @@ -293,6 +390,12 @@ def photons_absorbed( The exposure time of the measurement. If :obj:`None` (the default), the value in :attr:`timedelta_exposure` will be used. + integrate + Whether `image` is a single wavelength-integrated readout (as + produced by :meth:`expose` with ``integrate=True``). + If :obj:`True` (the default), the readout is spread uniformly across + the wavelength bins before the per-wavelength inverse, mirroring the + integration performed by :meth:`expose`. """ if axis_wavelength is None: shape_wavelength = na.shape(image.inputs.wavelength) @@ -306,8 +409,15 @@ def photons_absorbed( if timedelta is None: timedelta = self.timedelta_exposure + electrons = image.outputs + + if integrate: + # spread the integrated readout uniformly across the wavelength bins + num_wavelength = na.shape(image.inputs.wavelength)[axis_wavelength] - 1 + electrons = electrons / num_wavelength + photons = self.material.photons_absorbed( - electrons=image.outputs, + electrons=electrons, wavelength=image.inputs.wavelength.cell_centers(axis_wavelength), direction=direction, ) @@ -322,6 +432,7 @@ def uncertainty( ], direction: float | na.AbstractScalar = 1, axis_wavelength: None | str = None, + integrate: bool = True, ) -> na.FunctionArray[ na.SpectralPositionalVectorArray, na.AbstractScalar, @@ -330,10 +441,11 @@ def uncertainty( Compute the standard deviation of the noise in an image of electrons measured by the sensor. - This combines the material's analytic noise model + This uses the material's analytic per-wavelength noise model (:meth:`~optika.sensors.materials.AbstractSensorMaterial.uncertainty`), - which accounts for shot, Fano, and partial-charge-collection noise, with - the sensor's :attr:`read_noise` (added in quadrature), and so is the + which accounts for shot, Fano, and partial-charge-collection noise. If + `integrate` is :obj:`True`, the per-wavelength variances are summed in + quadrature and the sensor's :attr:`read_noise` is added once, giving the deterministic counterpart of the noise added by :meth:`expose`. Parameters @@ -350,6 +462,11 @@ def uncertainty( The logical axis of `image` corresponding to changing wavelength. If :obj:`None` (the default), ``image.inputs.wavelength`` must have only one logical axis. + integrate + Whether to integrate the noise over wavelength into a single + readout: the per-wavelength variances are summed in quadrature and + :attr:`read_noise` is added once. + Defaults to :obj:`True`, matching :meth:`expose`. """ if axis_wavelength is None: shape_wavelength = na.shape(image.inputs.wavelength) @@ -366,10 +483,16 @@ def uncertainty( direction=direction, ) - # add the read noise in quadrature - uncertainty = np.sqrt(np.square(uncertainty) + np.square(self.read_noise)) + inputs = image.inputs + + if integrate: + # sum the per-wavelength variances and add the read noise once + variance = np.square(uncertainty).sum(axis_wavelength) + variance = variance + np.square(self.read_noise) + uncertainty = np.sqrt(variance) + inputs = self._collapse_wavelength(inputs, axis_wavelength) - return dataclasses.replace(image, outputs=uncertainty) + return dataclasses.replace(image, inputs=inputs, outputs=uncertainty) def measure( self, @@ -380,6 +503,7 @@ def measure( where: bool | na.AbstractScalar = True, timedelta: None | u.Quantity | na.AbstractScalar = None, noise: bool = True, + integrate: bool = True, ) -> na.FunctionArray[ na.SpectralPositionalVectorArray, na.AbstractScalar, @@ -425,6 +549,7 @@ def measure( axis_wavelength=axis_wavelength, timedelta=timedelta, noise=noise, + integrate=integrate, ) diff --git a/optika/sensors/_sensors_test.py b/optika/sensors/_sensors_test.py index 2de337b4..4ab660c3 100644 --- a/optika/sensors/_sensors_test.py +++ b/optika/sensors/_sensors_test.py @@ -24,6 +24,15 @@ def test_read_noise(self, a: optika.sensors.AbstractImagingSensor): result = a.read_noise assert result >= 0 * u.electron + def test_pixels(self, a: optika.sensors.AbstractImagingSensor): + # the corners of the light-sensitive area map to pixel 0 and + # `num_pixel` + pixels_lower = a.pixels(a.aperture.bound_lower.xy) + pixels_upper = a.pixels(a.aperture.bound_upper.xy) + assert na.unit(pixels_lower).is_equivalent(u.pix) + assert np.allclose(pixels_lower, 0 * u.pix) + assert np.allclose(pixels_upper, a.num_pixel * u.pix) + @pytest.mark.parametrize( argnames="rays", argvalues=[ @@ -121,9 +130,9 @@ def test_photons_absorbed(self, a: optika.sensors.AbstractImagingSensor): outputs=rate, ) - # `photons_absorbed` is the deterministic inverse of `expose` - electrons = a.expose(image, noise=False) - result = a.photons_absorbed(electrons) + # per wavelength, `photons_absorbed` is the exact inverse of `expose` + electrons = a.expose(image, noise=False, integrate=False) + result = a.photons_absorbed(electrons, integrate=False) assert isinstance(result, na.FunctionArray) assert isinstance(result.inputs, na.SpectralPositionalVectorArray) @@ -133,6 +142,23 @@ def test_photons_absorbed(self, a: optika.sensors.AbstractImagingSensor): rate.to_value(u.photon / u.s), ) + # an integrated readout is spread back over the wavelength bins + integrated = na.FunctionArray( + inputs=na.SpectralPositionalVectorArray( + wavelength=wavelength, + position=position, + ), + outputs=na.random.uniform( + low=0, + high=1000, + shape_random={a.axis_pixel.x: 5, a.axis_pixel.y: 5}, + ) + * u.electron, + ) + result = a.photons_absorbed(integrated, integrate=True) + assert result.outputs.unit.is_equivalent(u.photon / u.s) + assert np.all(np.isfinite(result.outputs.to_value(u.photon / u.s))) + def test_uncertainty(self, a: optika.sensors.AbstractImagingSensor): # electrons measured in a few pixels, as a function of the wavelength # bin edges @@ -157,14 +183,22 @@ def test_uncertainty(self, a: optika.sensors.AbstractImagingSensor): outputs=electrons, ) + # integrated over wavelength (the default), with read noise once result = a.uncertainty(image) assert isinstance(result, na.FunctionArray) assert isinstance(result.inputs, na.SpectralPositionalVectorArray) assert result.outputs.unit.is_equivalent(u.electron) + assert "wavelength" not in na.shape(result.outputs) # the total noise is at least the read noise (added in quadrature) assert np.all(result.outputs >= a.read_noise) + # per wavelength (no integration, no read noise) + result = a.uncertainty(image, integrate=False) + assert result.outputs.unit.is_equivalent(u.electron) + assert "wavelength" in na.shape(result.outputs) + assert np.all(result.outputs >= 0 * u.electron) + @pytest.mark.parametrize( argnames="a", diff --git a/optika/systems/__init__.py b/optika/systems/__init__.py index ce20b15b..08c58064 100644 --- a/optika/systems/__init__.py +++ b/optika/systems/__init__.py @@ -3,10 +3,13 @@ """ from ._systems import AbstractSystem +from ._linear import AbstractLinearSystem, LinearSystem from ._sequential import AbstractSequentialSystem, SequentialSystem __all__ = [ "AbstractSystem", + "AbstractLinearSystem", + "LinearSystem", "AbstractSequentialSystem", "SequentialSystem", ] diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py new file mode 100644 index 00000000..4380c606 --- /dev/null +++ b/optika/systems/_linear.py @@ -0,0 +1,857 @@ +from typing import Any +import abc +import dataclasses +import numpy as np +import astropy.units as u +import astropy.constants +import named_arrays as na +import optika +from . import AbstractSystem + +__all__ = [ + "AbstractLinearSystem", + "LinearSystem", +] + + +def _radiance_to_unit( + radiance: na.AbstractScalar, + wavelength: na.AbstractScalar, + unit: None | u.UnitBase, +) -> na.AbstractScalar: + """ + Express a photon spectral radiance in `unit`, scaling by the energy per + photon (:math:`hc/\\lambda`) when `unit` is an energy rather than a photon + unit. A `unit` compatible with neither raises a + :class:`~astropy.units.UnitConversionError`. + """ + if unit is None: + return radiance + + if unit.is_equivalent(na.unit_normalized(radiance)): + return radiance.to(unit) + + # the natural (photon) radiance is scaled by the energy per photon to + # express it in energy units. + energy_per_photon = ( + astropy.constants.h * astropy.constants.c / wavelength / u.photon + ) + return (radiance * energy_per_photon).to(unit) + + +@dataclasses.dataclass(eq=False, repr=False) +class AbstractLinearSystem( + AbstractSystem, +): + """ + An interface for a linear forward model of an optical system. + + A linear system approximates an exact + :class:`~optika.systems.SequentialSystem` by characterizing it with a set + of precomputed models instead of raytracing each scene: a + :attr:`distortion` model mapping object-plane coordinates onto the detector, + an :attr:`area_effective` model, and optional :attr:`vignetting` and + :attr:`field_stop` models, together with a :attr:`sensor`. + The forward model, :meth:`image`, conservatively regrids the scene through + the distortion model, weighting each cell by the effective area, vignetting, + and field stop, and converts the result into detector electrons. + + Concrete subclasses supply these models as attributes. + """ + + @property + @abc.abstractmethod + def distortion(self) -> optika.distortion.AbstractDistortionModel: + """ + A distortion model which maps coordinates on the object plane to + positions on the detector plane. + """ + + @property + @abc.abstractmethod + def area_effective(self) -> optika.radiometry.AbstractEffectiveAreaModel: + """A model of the effective area of the system's entrance aperture.""" + + @property + @abc.abstractmethod + def vignetting(self) -> optika.radiometry.AbstractVignettingModel: + """ + A vignetting model which only transmits a fraction of the light as a + function of coordinates on the object plane. + """ + + @property + @abc.abstractmethod + def field_stop(self) -> optika.apertures.AbstractAperture: + """ + A model of the field stop which blocks light on the object plane. + """ + + @property + @abc.abstractmethod + def sensor(self) -> optika.sensors.AbstractImagingSensor: + """ + A model of the sensor which converts incident light intensity to an + electrical signal. + """ + + @property + @abc.abstractmethod + def direction(self): + """The cosine of the incidence angle on the sensor surface.""" + + @property + def coordinates_sensor(self) -> na.AbstractCartesian2dVectorArray: + """ + The vertices of the sensor pixel grid onto which the scene is + regridded, derived from :attr:`sensor`. + + These are the pixel edges in fractional pixel coordinates + (:meth:`~optika.sensors.AbstractImagingSensor.pixels`), matching the + pixel units of :attr:`distortion`. + """ + sensor = self.sensor + return na.Cartesian2dVectorLinearSpace( + start=sensor.pixels(sensor.aperture.bound_lower.xy), + stop=sensor.pixels(sensor.aperture.bound_upper.xy), + axis=sensor.axis_pixel, + num=sensor.num_pixel + 1, + ) + + @property + def shape(self) -> dict[str, int]: + """ + The broadcasted shape of the component models, which is the shape of + any array-valued parametrization of this system. + """ + return na.broadcast_shapes( + optika.shape(self.distortion), + optika.shape(self.area_effective), + optika.shape(self.vignetting), + optika.shape(self.field_stop), + optika.shape(self.sensor), + optika.shape(self.direction), + ) + + @property + def transformation(self) -> None: + """ + A linear system has no geometric placement, so it has no + transformation into the global coordinate system. + """ + return None + + def weights( + self, + coordinates: na.SpectralPositionalVectorArray, + axis_wavelength: str, + axis_field: tuple[str, str], + ) -> tuple[na.AbstractScalar, dict[str, int], dict[str, int]]: + """ + Compute the weights which map the overlap of each pixel on the object + plane to each pixel on the sensor plane. + + Parameters + ---------- + coordinates + The vertices of each pixel on the object plane. + axis_wavelength + The logical axis corresponding to changing wavelength coordinate. + axis_field + The logical axes corresponding to changing field coordinate. + """ + + coordinates = coordinates.spectral_positional + + coordinates = coordinates.cell_centers(axis_wavelength) + + position_sensor = self.distortion.distort(coordinates).position + + # the conservative-regridding weights apply per input *cell*, but + # `coordinates` are the cell vertices, so evaluate the radiometric + # factors at the field cell centers (one fewer point along each axis). + coordinates_cell = coordinates.cell_centers(axis_field) + + vignetting = self.vignetting + if vignetting is not None: + weights_vignetting = vignetting(coordinates_cell) + else: + weights_vignetting = 1 + + field_stop = self.field_stop + if field_stop is not None: + weights_stop = field_stop( + position=na.Cartesian3dVectorArray( + x=coordinates_cell.position.x, + y=coordinates_cell.position.y, + ), + ) + else: + weights_stop = 1 + + weights_area = self.area_effective(coordinates_cell.wavelength) + + weights_input = ( + weights_vignetting * weights_stop * weights_area.to_value(self.weights_unit) + ) + + axis_pixel = self.sensor.axis_pixel + + result = na.regridding.weights( + coordinates_input=position_sensor, + coordinates_output=self.coordinates_sensor, + axis_input=axis_field, + axis_output=(axis_pixel.x, axis_pixel.y), + weights_input=weights_input, + method="conservative", + ) + + return result + + def weights_transposed( + self, + weights: tuple[na.AbstractScalar, dict[str, int], dict[str, int]], + coordinates: na.SpectralPositionalVectorArray, + axis_wavelength: str, + axis_field: tuple[str, str], + ) -> tuple[na.AbstractScalar, dict[str, int], dict[str, int]]: + """ + Compute the weights which map the overlap of each pixel on the object + plane to each pixel on the detector plane. + + Parameters + ---------- + coordinates + The vertices of each pixel on the object plane. + axis_wavelength + The logical axis corresponding to changing wavelength coordinate. + axis_field + The logical axes corresponding to changing field coordinate. + """ + + coordinates = coordinates.spectral_positional + + coordinates = coordinates.cell_centers(axis_wavelength) + + position_sensor = self.distortion.distort(coordinates).position + + # the conservative-regridding weights apply per input *cell*, but + # `coordinates` are the cell vertices, so evaluate the radiometric + # factors at the field cell centers (one fewer point along each axis). + coordinates_cell = coordinates.cell_centers(axis_field) + + vignetting = self.vignetting + if vignetting is not None: + weights_vignetting = vignetting(coordinates_cell) + else: + weights_vignetting = 1 + + weights_area = self.area_effective(coordinates_cell.wavelength) + + weights_input = weights_vignetting * weights_area.to_value(self.weights_unit) + + axis_pixel = self.sensor.axis_pixel + + result = na.regridding.transpose_weights_conservative( + weights=weights, + coordinates_input=position_sensor, + coordinates_output=self.coordinates_sensor, + axis_input=axis_field, + axis_output=(axis_pixel.x, axis_pixel.y), + weights_input=weights_input, + ) + + return result + + @property + def weights_unit(self): + """The units associated with :attr:`weights`.""" + return u.cm**2 + + def image_from_weights( + self, + weights: tuple[na.AbstractScalar, dict[str, int], dict[str, int]], + scene: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar], + *, + axis_wavelength: None | str = None, + axis_field: None | tuple[str, str] = None, + integrate: bool = True, + noise: bool = True, + uncertainty: bool = False, + **kwargs: Any, + ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: + """ + Apply a precomputed set of regridding weights to a scene, returning the + electrons measured by the sensor. + + This reuses the (expensive) conservative-regridding operator built by + :meth:`weights` for any scene: it integrates the spectral radiance over + each object-plane voxel, regrids onto the sensor plane, and applies the + sensor response + (:meth:`~optika.sensors.AbstractImagingSensor.expose`). + :meth:`image` is the special case that builds `weights` on the fly. + + Parameters + ---------- + weights + The conservative-regridding weights computed by :meth:`weights`. + scene + The spectral radiance of the observed scene, sampled on the + vertices of each pixel on the object plane. + The radiance may be given in either energy or photon units; the + sensor converts energy to photons if necessary. + axis_wavelength + The logical axis of `scene` corresponding to changing wavelength. + If :obj:`None` (the default), the single axis of + ``scene.inputs.wavelength`` that is not a field axis is used; a + :class:`ValueError` is raised if there is more than one. + axis_field + The logical axes of `scene` corresponding to changing position on + the object plane. + If :obj:`None` (the default), the axes of ``scene.inputs.position`` + are used. + integrate + Whether to integrate the electrons over wavelength into a single + sensor readout, applying the read noise once + (see :meth:`~optika.sensors.AbstractImagingSensor.expose`). + Defaults to :obj:`True`. + noise + Whether to include sensor noise in the result. + uncertainty + Whether to attach the standard deviation of the measurement noise + to the result, as a + :class:`~named_arrays.NormalUncertainScalarArray`, using the + sensor's :meth:`~optika.sensors.AbstractImagingSensor.uncertainty`. + kwargs + Additional keyword arguments passed to the sensor's + :meth:`~optika.sensors.AbstractImagingSensor.expose` method, such + as `timedelta`. + """ + + scene = scene.explicit + coordinates = scene.inputs + + if axis_field is None: + axis_field = tuple(na.shape(coordinates.position)) + + if axis_wavelength is None: + axis = set(na.shape(coordinates.wavelength)) - set(axis_field) + if len(axis) != 1: # pragma: nocover + raise ValueError( + f"unable to infer `axis_wavelength`: expected exactly one " + f"axis of `scene.inputs.wavelength` " + f"({na.shape(coordinates.wavelength)}) that is not a field " + f"axis ({axis_field}), got {axis}." + ) + (axis_wavelength,) = axis + + # volume of each voxel on the object plane: the spectral bin width + # times the solid angle (or area) subtended by each field pixel. + # `.spectral_positional` accepts a scene defined on a Doppler grid (as + # produced by the ctis `IdealInstrument`). + volume = coordinates.spectral_positional.volume_cell( + (axis_wavelength, *axis_field) + ) + + # integrate the spectral radiance over each voxel into a flux per unit + # collecting area. + rate = scene.outputs * volume + + rate_sensor = na.regridding.regrid_from_weights( + *weights, + values_input=rate, + ) + + # restore the unit stripped from `weights_input` inside `weights` + # (the effective collecting area), turning the rate per unit area into + # a rate per sensor pixel. + rate_sensor = rate_sensor * self.weights_unit + + # the flux incident on the sensor cannot be negative + rate_sensor = np.maximum(rate_sensor, 0) + + image = na.FunctionArray( + inputs=na.SpectralPositionalVectorArray( + wavelength=coordinates.wavelength, + position=self.coordinates_sensor, + ), + outputs=rate_sensor, + ) + + return self.sensor.expose( + image=image, + direction=self.direction, + axis_wavelength=axis_wavelength, + noise=noise, + integrate=integrate, + uncertainty=uncertainty, + **kwargs, + ) + + def backproject_from_weights( + self, + weights: tuple[na.AbstractScalar, dict[str, int], dict[str, int]], + image: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar], + coordinates: na.SpectralPositionalVectorArray, + *, + axis_wavelength: None | str = None, + axis_field: None | tuple[str, str] = None, + integrate: bool = True, + unit: None | u.UnitBase = None, + ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: + """ + Apply a precomputed set of transposed regridding weights to a + detector-plane image of electrons, returning the backprojected spectral + radiance. + + This is the transpose of :meth:`image_from_weights`: + :meth:`weights_transposed` builds the (expensive) transposed operator + once, and this method reuses it to invert the sensor response + (:meth:`~optika.sensors.AbstractImagingSensor.photons_absorbed`), + spread the result back onto the object plane, and divide out the voxel + volume. :meth:`backproject` is the special case that builds the weights + on the fly. + + Parameters + ---------- + weights + The transposed regridding weights computed by + :meth:`weights_transposed`. + image + The detector-plane image of electrons to project onto the object + plane, as produced by :meth:`image`. + coordinates + The vertices of each pixel on the object plane to project onto. + axis_wavelength + The logical axis of `coordinates` corresponding to changing + wavelength. + If :obj:`None` (the default), the single axis of + ``coordinates.wavelength`` that is not a field axis is used; a + :class:`ValueError` is raised if there is more than one. + axis_field + The logical axes of `coordinates` corresponding to changing position + on the object plane. + If :obj:`None` (the default), the axes of ``coordinates.position`` + are used. + integrate + Whether `image` is a single wavelength-integrated readout, to be + spread back over the wavelength bins before the transpose. + Defaults to :obj:`True`. + unit + The unit of the backprojected spectral radiance. + :meth:`image_from_weights` accepts a scene in either photon or + energy units, so the backprojection is expressed in whichever the + caller requests, converting between photon and energy units using + the energy per photon (:math:`hc/\\lambda`). + If :obj:`None` (the default), the radiance is left in the natural + (photon) units of the backprojection and is not converted. + """ + + coordinates = coordinates.explicit + + if axis_field is None: + axis_field = tuple(na.shape(coordinates.position)) + + if axis_wavelength is None: + axis = set(na.shape(coordinates.wavelength)) - set(axis_field) + if len(axis) != 1: # pragma: nocover + raise ValueError( + f"unable to infer `axis_wavelength`: expected exactly one " + f"axis of `coordinates.wavelength` " + f"({na.shape(coordinates.wavelength)}) that is not a field " + f"axis ({axis_field}), got {axis}." + ) + (axis_wavelength,) = axis + + # volume of each voxel on the object plane: the spectral bin width + # times the solid angle (or area) subtended by each field pixel. + # `.spectral_positional` accepts a Doppler grid (as used by the ctis + # `IdealInstrument`); the radiance is returned on the original grid. + volume = coordinates.spectral_positional.volume_cell( + (axis_wavelength, *axis_field) + ) + + # an integrated readout only carries the two band edges, but the sensor + # inverse spreads it across the wavelength bins it is reconstructed + # onto; give it the full target grid so it spreads over the right number + # of bins (and evaluates the quantum efficiency per bin) instead of + # treating the readout as a single bin. + if integrate: + image = image.replace( + inputs=image.inputs.replace( + wavelength=coordinates.spectral_positional.wavelength, + ), + ) + + # invert the detector response, mapping the measured electrons back into + # the photon rate per pixel produced by `image_from_weights`. + image = self.sensor.photons_absorbed( + image, + direction=self.direction, + axis_wavelength=axis_wavelength, + integrate=integrate, + ) + + radiance = na.regridding.regrid_from_weights( + *weights, + values_input=image.outputs, + ) + + # divide out the effective collecting area, the inverse of the + # multiplication performed by `image_from_weights`, converting the + # per-pixel rate back into a rate per unit collecting area. + radiance = radiance / self.weights_unit + + # recover the spectral radiance by undoing the integration over each + # object-plane voxel performed by `image`. + radiance = radiance / volume + + # express the (photon) radiance in the requested unit, converting to + # energy units with the energy per photon if necessary. + radiance = _radiance_to_unit( + radiance=radiance, + wavelength=coordinates.wavelength.cell_centers(axis_wavelength), + unit=unit, + ) + + return na.FunctionArray( + inputs=coordinates, + outputs=radiance, + ) + + def image( + self, + scene: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar], + *, + axis_wavelength: None | str = None, + axis_field: None | tuple[str, str] = None, + integrate: bool = True, + noise: bool = True, + uncertainty: bool = False, + **kwargs: Any, + ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: + """ + Linear forward model of the optical system. + + Maps the spectral radiance of a scene to the electrons measured by the + sensor. + + Parameters + ---------- + scene + The spectral radiance of the observed scene, sampled on the + vertices of each pixel on the object plane. + The radiance may be given in either energy units (such as + :math:`W / cm^2 / arcsec^2 / nm`) or photon units (such as + :math:`photon / s / cm^2 / arcsec^2 / nm`); the sensor converts + energy to photons if necessary. + axis_wavelength + The logical axis of `scene` corresponding to changing wavelength. + If :obj:`None` (the default), the single axis of + ``scene.inputs.wavelength`` that is not a field axis is used; a + :class:`ValueError` is raised if there is more than one. + axis_field + The logical axes of `scene` corresponding to changing position on + the object plane. + If :obj:`None` (the default), the axes of ``scene.inputs.position`` + are used. + integrate + Whether to integrate the electrons over wavelength into a single + sensor readout, applying the read noise once + (see :meth:`~optika.sensors.AbstractImagingSensor.expose`). + Defaults to :obj:`True`. + noise + Whether to include sensor noise in the result. + uncertainty + Whether to attach the standard deviation of the measurement noise + to the result, as a + :class:`~named_arrays.NormalUncertainScalarArray`, using the + sensor's :meth:`~optika.sensors.AbstractImagingSensor.uncertainty`. + kwargs + Additional keyword arguments passed to the sensor's + :meth:`~optika.sensors.AbstractImagingSensor.expose` method, such + as `timedelta`. + """ + + scene = scene.explicit + coordinates = scene.inputs + + if axis_field is None: + axis_field = tuple(na.shape(coordinates.position)) + + if axis_wavelength is None: + axis = set(na.shape(coordinates.wavelength)) - set(axis_field) + if len(axis) != 1: # pragma: nocover + raise ValueError( + f"unable to infer `axis_wavelength`: expected exactly one " + f"axis of `scene.inputs.wavelength` " + f"({na.shape(coordinates.wavelength)}) that is not a field " + f"axis ({axis_field}), got {axis}." + ) + (axis_wavelength,) = axis + + # `weights` folds in the effective area, vignetting, and field stop and + # is the expensive part of the forward model; `image_from_weights` + # reuses it to integrate the radiance, regrid, and expose the sensor. + weights = self.weights( + coordinates=coordinates, + axis_wavelength=axis_wavelength, + axis_field=axis_field, + ) + + return self.image_from_weights( + weights, + scene, + axis_wavelength=axis_wavelength, + axis_field=axis_field, + noise=noise, + uncertainty=uncertainty, + integrate=integrate, + **kwargs, + ) + + def backproject( + self, + image: ( + na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar] + | na.AbstractScalar + ), + coordinates: na.SpectralPositionalVectorArray, + *, + axis_wavelength: None | str = None, + axis_field: None | tuple[str, str] = None, + integrate: bool = True, + unit: None | u.UnitBase = None, + weights: None | tuple[na.AbstractScalar, dict[str, int], dict[str, int]] = None, + ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: + """ + Transpose of the linear forward model, :meth:`image`. + + Inverts the detector response (:meth:`~optika.sensors.AbstractImagingSensor.expose`) + and then applies the transpose of the optical regridding, projecting a + detector-plane image of electrons back onto the object plane by + spreading each pixel's value across every object-plane cell that could + have contributed to it. This is the transpose of :meth:`image`, not its + inverse (the geometric spreading is not undone, and the sensor noise is + not recovered). + + Parameters + ---------- + image + The detector-plane image of electrons to project onto the object + plane, as produced by :meth:`image`. + coordinates + The vertices of each pixel on the object plane to project onto. + axis_wavelength + The logical axis of `coordinates` corresponding to changing + wavelength. + If :obj:`None` (the default), the single axis of + ``coordinates.wavelength`` that is not a field axis is used; a + :class:`ValueError` is raised if there is more than one. + axis_field + The logical axes of `coordinates` corresponding to changing position + on the object plane. + If :obj:`None` (the default), the axes of ``coordinates.position`` + are used. + integrate + Whether `image` is a single wavelength-integrated readout, to be + spread back over the wavelength bins before the transpose. + Defaults to :obj:`True`. + unit + The unit of the backprojected spectral radiance. + :meth:`image` accepts a scene in either photon or energy units, so + the backprojection is expressed in whichever the caller requests, + converting between photon and energy units using the energy per + photon (:math:`hc/\\lambda`). + If :obj:`None` (the default), the radiance is left in the natural + (photon) units of the backprojection and is not converted. + weights + The forward regridding weights computed by :meth:`weights`. + If :obj:`None` (the default), they are computed from `coordinates`. + Supplying a precomputed value (for example one already built for + :meth:`image`) avoids rebuilding the expensive regridding operator. + """ + + coordinates = coordinates.explicit + + if axis_field is None: + axis_field = tuple(na.shape(coordinates.position)) + + if axis_wavelength is None: + axis = set(na.shape(coordinates.wavelength)) - set(axis_field) + if len(axis) != 1: # pragma: nocover + raise ValueError( + f"unable to infer `axis_wavelength`: expected exactly one " + f"axis of `coordinates.wavelength` " + f"({na.shape(coordinates.wavelength)}) that is not a field " + f"axis ({axis_field}), got {axis}." + ) + (axis_wavelength,) = axis + + # `weights_transposed` is the expensive transposed operator; + # `backproject_from_weights` reuses it to invert the sensor response, + # regrid onto the object plane, and recover the spectral radiance. + if weights is None: + weights = self.weights( + coordinates=coordinates, + axis_wavelength=axis_wavelength, + axis_field=axis_field, + ) + weights_transposed = self.weights_transposed( + weights=weights, + coordinates=coordinates, + axis_wavelength=axis_wavelength, + axis_field=axis_field, + ) + + return self.backproject_from_weights( + weights_transposed, + image, + coordinates=coordinates, + axis_wavelength=axis_wavelength, + axis_field=axis_field, + integrate=integrate, + unit=unit, + ) + + +@dataclasses.dataclass(eq=False, repr=False) +class LinearSystem( + AbstractLinearSystem, +): + """ + A linear forward model of an optical system, assembled from precomputed + distortion, effective area, and (optionally) vignetting and field stop + models. + + This is a fast approximation to + :class:`~optika.systems.SequentialSystem`. Once a sequential system has been + characterized (for example, its :meth:`~optika.systems.SequentialSystem.distortion`, + :meth:`~optika.systems.SequentialSystem.vignetting`, and + :meth:`~optika.systems.SequentialSystem.area_effective` models have been fit), + those models can be reused here to image many scenes without raytracing each + one. + + Examples + -------- + + Simulate an image of an airforce target using a simple spectrograph model. + + .. jupyter-execute:: + + import matplotlib.pyplot as plt + import astropy.units as u + import astropy.visualization + import named_arrays as na + import optika + + # The distortion, effective area, and sensor models that + # define the system. + distortion = optika.distortion.SimpleDistortionModel( + plate_scale=0.75 * u.arcsec / u.pix, + dispersion=3.75 * u.nm / u.pix, + angle=0 * u.deg, + reference=na.SpectralPositionalVectorArray( + wavelength=550 * u.nm, + position=na.Cartesian2dVectorArray(32, 32) * u.pix, + ), + ) + area_effective = optika.radiometry.InterpolatedEffectiveAreaModel( + wavelength=na.linspace(500, 600, axis="wavelength", num=10) * u.nm, + area=na.linspace(1, 2, axis="wavelength", num=10) * u.cm ** 2, + axis_wavelength="wavelength", + ) + sensor = optika.sensors.ImagingSensor( + width_pixel=15 * u.um, + axis_pixel=na.Cartesian2dVectorArray("detector_x", "detector_y"), + timedelta_exposure=1 * u.s, + num_pixel=na.Cartesian2dVectorArray(64, 64), + ) + + # Assemble the linear system from the component models. + system = optika.systems.LinearSystem( + area_effective=area_effective, + distortion=distortion, + sensor=sensor, + ) + + # Define the number of field points to sample. + num_field = 2 * system.sensor.num_pixel + + # Define the scene as an airforce target. The coordinates (inputs) + # are defined on cell vertices and the spectral radiance (outputs) + # on cell centers. + scene = na.FunctionArray( + inputs=na.SpectralPositionalVectorArray( + wavelength=na.linspace(549, 551, axis="wavelength", num=4) * u.nm, + position=na.Cartesian2dVectorLinearSpace( + start=-15 * u.arcsec, + stop=+15 * u.arcsec, + axis=na.Cartesian2dVectorArray("field_x", "field_y"), + num=num_field + 1, + ), + ), + outputs=optika.targets.airforce( + axis_x="field_x", + axis_y="field_y", + num_x=num_field.x, + num_y=num_field.y, + ) * 1e-16 * u.W / u.cm ** 2 / u.arcsec ** 2 / u.nm, + ) + + # Simulate an image of the scene using the linear forward model. + image = system.image(scene, noise=False) + + # Plot the original scene and the simulated image. + with astropy.visualization.quantity_support(): + fig, ax = plt.subplots( + ncols=2, + figsize=(8, 5), + constrained_layout=True, + ) + na.plt.pcolormesh( + scene.inputs.position, + C=scene.outputs.value, + ax=ax[0], + ) + na.plt.pcolormesh( + image.inputs.position, + C=image.outputs.value, + ax=ax[1], + ) + ax[0].set_title("scene") + ax[1].set_title("image") + """ + + area_effective: optika.radiometry.AbstractEffectiveAreaModel = dataclasses.MISSING + """A model of the effective area of the system's entrance aperture.""" + + distortion: optika.distortion.AbstractDistortionModel = dataclasses.MISSING + """ + A distortion model which maps coordinates on the object plane to + positions on the detector plane. + """ + + sensor: optika.sensors.AbstractImagingSensor = dataclasses.MISSING + """ + A model of the sensor which converts incident light intensity to an + electrical signal. + """ + + direction: float | na.AbstractScalar = 1 + """ + The cosine of the incidence angle of the light striking the sensor surface. + """ + + vignetting: None | optika.radiometry.AbstractVignettingModel = None + """ + A vignetting model which only transmits a fraction of the light as a + function of coordinates on the object plane. + If :obj:`None` (the default), the system will have no vignetting. + """ + + field_stop: None | optika.apertures.AbstractAperture = None + """ + A model of the field stop which blocks light on the object plane. + If :obj:`None` (the default), the sensor will be the field stop. + """ diff --git a/optika/systems/_linear_test.py b/optika/systems/_linear_test.py new file mode 100644 index 00000000..c35d8b9a --- /dev/null +++ b/optika/systems/_linear_test.py @@ -0,0 +1,290 @@ +import dataclasses + +import pytest +import numpy as np +import astropy.units as u +import named_arrays as na +import optika +from ._systems_test import AbstractTestAbstractSystem + + +def _distortion() -> optika.distortion.SimpleDistortionModel: + return optika.distortion.SimpleDistortionModel( + plate_scale=0.75 * u.arcsec / u.pix, + dispersion=3.75 * u.nm / u.pix, + angle=0 * u.deg, + reference=na.SpectralPositionalVectorArray( + wavelength=550 * u.nm, + position=na.Cartesian2dVectorArray(16, 16) * u.pix, + ), + ) + + +def _area_effective() -> optika.radiometry.InterpolatedEffectiveAreaModel: + return optika.radiometry.InterpolatedEffectiveAreaModel( + wavelength=na.linspace(400, 700, axis="wavelength", num=10) * u.nm, + area=na.linspace(1, 2, axis="wavelength", num=10) * u.cm**2, + axis_wavelength="wavelength", + ) + + +def _vignetting() -> optika.radiometry.PolynomialVignettingModel: + scene = na.SpectralPositionalVectorArray( + wavelength=na.linspace(500, 600, axis="wavelength", num=3) * u.nm, + position=na.Cartesian2dVectorLinearSpace( + start=-10 * u.arcsec, + stop=+10 * u.arcsec, + axis=na.Cartesian2dVectorArray("field_x", "field_y"), + num=5, + ), + ) + return optika.radiometry.PolynomialVignettingModel( + coordinates_scene=scene, + illumination=1 - 0.001 * (scene.position.length / u.arcsec) ** 2, + axis_wavelength="wavelength", + axis_field=("field_x", "field_y"), + degree=2, + ) + + +def _sensor() -> optika.sensors.ImagingSensor: + return optika.sensors.ImagingSensor( + width_pixel=15 * u.um, + axis_pixel=na.Cartesian2dVectorArray("detector_x", "detector_y"), + timedelta_exposure=1 * u.s, + num_pixel=na.Cartesian2dVectorArray(32, 32), + ) + + +def _scene(radiance: u.Quantity | na.AbstractScalar) -> na.FunctionArray: + return na.FunctionArray( + inputs=na.SpectralPositionalVectorArray( + wavelength=na.linspace(500, 600, axis="wavelength", num=4) * u.nm, + position=na.Cartesian2dVectorLinearSpace( + start=-10 * u.arcsec, + stop=+10 * u.arcsec, + axis=na.Cartesian2dVectorArray("field_x", "field_y"), + num=11, + ), + ), + outputs=na.random.uniform( + low=0 * radiance, + high=radiance, + shape_random=dict(field_x=10, field_y=10), + ), + ) + + +def _scene_doppler(radiance: u.Quantity | na.AbstractScalar) -> na.FunctionArray: + return na.FunctionArray( + inputs=na.DopplerPositionalVectorArray( + wavelength=na.linspace(500, 600, axis="wavelength", num=4) * u.nm, + wavelength_rest=550 * u.nm, + position=na.Cartesian2dVectorLinearSpace( + start=-10 * u.arcsec, + stop=+10 * u.arcsec, + axis=na.Cartesian2dVectorArray("field_x", "field_y"), + num=11, + ), + ), + outputs=na.random.uniform( + low=0 * radiance, + high=radiance, + shape_random=dict(field_x=10, field_y=10), + ), + ) + + +class AbstractTestAbstractLinearSystem( + AbstractTestAbstractSystem, +): + def test_distortion(self, a: optika.systems.AbstractLinearSystem): + assert isinstance(a.distortion, optika.distortion.AbstractDistortionModel) + + def test_area_effective(self, a: optika.systems.AbstractLinearSystem): + assert isinstance( + a.area_effective, + optika.radiometry.AbstractEffectiveAreaModel, + ) + + def test_sensor(self, a: optika.systems.AbstractLinearSystem): + assert isinstance(a.sensor, optika.sensors.AbstractImagingSensor) + + def test_coordinates_sensor(self, a: optika.systems.AbstractLinearSystem): + result = a.coordinates_sensor + assert isinstance(result, na.AbstractCartesian2dVectorArray) + assert na.unit_normalized(result).is_equivalent(u.pix) + + @pytest.mark.parametrize( + argnames="radiance", + argvalues=[ + 1e-18 * u.W / u.cm**2 / u.arcsec**2 / u.nm, + 1e3 * u.photon / u.s / u.cm**2 / u.arcsec**2 / u.nm, + ], + ) + @pytest.mark.parametrize("noise", [False, True]) + def test_image( + self, + a: optika.systems.AbstractLinearSystem, + noise: bool, + radiance: u.Quantity, + ): + scene = _scene(radiance) + result = a.image(scene, noise=noise) + assert isinstance(result, na.FunctionArray) + assert isinstance(result.inputs, na.SpectralPositionalVectorArray) + assert na.unit(result.outputs).is_equivalent(u.electron) + assert np.all(np.isfinite(result.outputs.value)) + + # the scene is binned onto the sensor pixel grid + axis_pixel = a.sensor.axis_pixel + assert axis_pixel.x in na.shape(result.outputs) + assert axis_pixel.y in na.shape(result.outputs) + + if not noise: + assert np.all(result.outputs >= 0 * u.electron) + assert result.outputs.sum() > 0 * u.electron + + def test_image_doppler(self, a: optika.systems.AbstractLinearSystem): + # a scene defined on a Doppler grid (as produced by the ctis + # `IdealInstrument`) is accepted and imaged onto the sensor. + scene = _scene_doppler(1e3 * u.photon / u.s / u.cm**2 / u.arcsec**2 / u.nm) + result = a.image(scene, noise=False) + assert isinstance(result, na.FunctionArray) + assert isinstance(result.inputs, na.SpectralPositionalVectorArray) + assert na.unit(result.outputs).is_equivalent(u.electron) + assert np.all(np.isfinite(result.outputs.value)) + assert result.outputs.sum() > 0 * u.electron + + @pytest.mark.parametrize( + argnames="radiance", + argvalues=[ + 1e3 * u.photon / u.s / u.cm**2 / u.arcsec**2 / u.nm, + ], + ) + def test_backproject( + self, + a: optika.systems.AbstractLinearSystem, + radiance: u.Quantity, + ): + scene = _scene(radiance) + image = a.image(scene, noise=False) + result = a.backproject(image, scene.inputs) + assert isinstance(result, na.FunctionArray) + assert isinstance(result.inputs, na.SpectralPositionalVectorArray) + + # the detector response is inverted, so the backprojection recovers a + # spectral radiance with the same units as the original scene. + assert na.unit(result.outputs).is_equivalent(na.unit(radiance)) + assert np.all(np.isfinite(result.outputs.value)) + assert result.outputs.sum() > 0 * na.unit(radiance) + + def test_backproject_doppler(self, a: optika.systems.AbstractLinearSystem): + # a Doppler object-plane grid is accepted, and the backprojected + # radiance is returned on that same Doppler grid. + scene = _scene_doppler(1e3 * u.photon / u.s / u.cm**2 / u.arcsec**2 / u.nm) + image = a.image(scene, noise=False) + result = a.backproject(image, scene.inputs) + assert isinstance(result, na.FunctionArray) + assert isinstance(result.inputs, na.DopplerPositionalVectorArray) + assert na.unit(result.outputs).is_equivalent( + u.photon / u.s / u.cm**2 / u.arcsec**2 / u.nm + ) + assert np.all(np.isfinite(result.outputs.value)) + + def test_backproject_unit(self, a: optika.systems.AbstractLinearSystem): + scene = _scene(1e3 * u.photon / u.s / u.cm**2 / u.arcsec**2 / u.nm) + image = a.image(scene, noise=False) + + # the backprojection can be expressed in photon or energy units + unit_photon = u.photon / u.s / u.cm**2 / u.arcsec**2 / u.nm + unit_energy = u.W / u.cm**2 / u.arcsec**2 / u.nm + + result_photon = a.backproject(image, scene.inputs, unit=unit_photon) + result_energy = a.backproject(image, scene.inputs, unit=unit_energy) + + assert na.unit_normalized(result_photon.outputs).is_equivalent(unit_photon) + assert na.unit_normalized(result_energy.outputs).is_equivalent(unit_energy) + assert np.all(np.isfinite(result_energy.outputs.value)) + + # a unit compatible with neither is rejected + with pytest.raises(u.UnitConversionError): + a.backproject(image, scene.inputs, unit=u.s) + + def test_from_weights_infer_axes(self, a: optika.systems.AbstractLinearSystem): + # `image_from_weights`/`backproject_from_weights` infer `axis_wavelength` + # and `axis_field` from the scene when they are not given explicitly. + scene = _scene(1e3 * u.photon / u.s / u.cm**2 / u.arcsec**2 / u.nm) + axis_wavelength = "wavelength" + axis_field = ("field_x", "field_y") + + weights = a.weights(scene.inputs, axis_wavelength, axis_field) + + image = a.image_from_weights(weights, scene, noise=False, integrate=False) + assert isinstance(image, na.FunctionArray) + assert na.unit(image.outputs).is_equivalent(u.electron) + assert np.all(np.isfinite(image.outputs.value)) + + weights_transposed = a.weights_transposed( + weights, scene.inputs, axis_wavelength, axis_field + ) + result = a.backproject_from_weights( + weights_transposed, image, scene.inputs, integrate=False + ) + assert isinstance(result, na.FunctionArray) + assert na.unit(result.outputs).is_equivalent( + u.photon / u.s / u.cm**2 / u.arcsec**2 / u.nm + ) + assert np.all(np.isfinite(result.outputs.value)) + + def test_image_area_unit_invariance( + self, + a: optika.systems.AbstractLinearSystem, + ): + # the result must not depend on the unit the effective area model + # happens to be expressed in, since `weights` strips the unit and + # `image_from_weights` restores it as `weights_unit` + area = a.area_effective + if not hasattr(area, "area"): # pragma: nocover + pytest.skip("requires an area model with an `area` attribute") + b = dataclasses.replace( + a, + area_effective=dataclasses.replace(area, area=area.area.to(u.mm**2)), + ) + scene = _scene(1e3 * u.photon / u.s / u.cm**2 / u.arcsec**2 / u.nm) + result_a = a.image(scene, noise=False) + result_b = b.image(scene, noise=False) + assert np.allclose(result_a.outputs, result_b.outputs) + + def test_image_uncertainty(self, a: optika.systems.AbstractLinearSystem): + scene = _scene(1e3 * u.photon / u.s / u.cm**2 / u.arcsec**2 / u.nm) + result = a.image(scene, noise=False, uncertainty=True) + + # the measurement noise is attached as a normal uncertain array + assert isinstance(result.outputs, na.NormalUncertainScalarArray) + assert na.unit(result.outputs.nominal).is_equivalent(u.electron) + assert na.unit(result.outputs.width).is_equivalent(u.electron) + assert np.all(result.outputs.width >= 0 * u.electron) + + +@pytest.mark.parametrize( + argnames="a", + argvalues=[ + optika.systems.LinearSystem( + area_effective=_area_effective(), + distortion=_distortion(), + sensor=_sensor(), + ), + optika.systems.LinearSystem( + area_effective=_area_effective(), + distortion=_distortion(), + sensor=_sensor(), + vignetting=_vignetting(), + field_stop=optika.apertures.RectangularAperture(half_width=15 * u.arcsec), + ), + ], +) +class TestLinearSystem( + AbstractTestAbstractLinearSystem, +): + pass diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index a81c8f1c..a1e617a6 100644 --- a/optika/systems/_sequential.py +++ b/optika/systems/_sequential.py @@ -1,4 +1,3 @@ -from __future__ import annotations from typing import Sequence, Callable, Any, ClassVar import abc import dataclasses @@ -14,6 +13,7 @@ import named_arrays as na import optika from . import AbstractSystem +from . import LinearSystem __all__ = [ "AbstractSequentialSystem", @@ -23,6 +23,8 @@ @dataclasses.dataclass(eq=False, repr=False) class AbstractSequentialSystem( + optika.mixins.DxfWritable, + optika.mixins.Plottable, AbstractSystem, ): """ @@ -1088,12 +1090,13 @@ def _rayfunction_from_vertices( def image( self, scene: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar], - pupil: None | na.AbstractCartesian2dVectorArray = None, + *, axis_wavelength: None | str = None, axis_field: None | tuple[str, str] = None, - axis_pupil: None | tuple[str, str] = None, integrate: bool = True, noise: bool = True, + pupil: None | na.AbstractCartesian2dVectorArray = None, + axis_pupil: None | tuple[str, str] = None, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: """ Forward model of the optical system. @@ -1105,10 +1108,6 @@ def image( The spectral radiance of the scene as a function of wavelength and field position. The inputs must be cell vertices. - pupil - The vertices of the pupil grid in either normalized or physical - coordinates. - If :obj:`None` (the default), the pupil grid will only have one cell. axis_wavelength The logical axis of `scene` corresponding to changing wavelength coordinate. If :obj:`None`, @@ -1119,12 +1118,6 @@ def image( If :obj:`None`, ``set(scene.inputs.position.shape) - set(self.shape) - {axis_wavelength}``, should have exactly two elements. - axis_pupil - The two logical axes of `pupil` corresponding to changing pupil coordinate. - If :obj:`None`, - ``set(pupil.shape) - set(self.shape) - {axis_wavelength,} - set(axis_field)``, - should have exactly two elements. - If `pupil` is :obj:`None`, this parameter is ignored. integrate Whether to integrate the wavelength axis. Real images usually have the wavelength axis integrated since they use @@ -1132,6 +1125,16 @@ def image( separate for introspective purposes. noise Whether to add noise to the result. + pupil + The vertices of the pupil grid in either normalized or physical + coordinates. + If :obj:`None` (the default), the pupil grid will only have one cell. + axis_pupil + The two logical axes of `pupil` corresponding to changing pupil coordinate. + If :obj:`None`, + ``set(pupil.shape) - set(self.shape) - {axis_wavelength,} - set(axis_field)``, + should have exactly two elements. + If `pupil` is :obj:`None`, this parameter is ignored. """ scene = scene.explicit @@ -1188,21 +1191,34 @@ def image( normalized_pupil=normalized_pupil, ) - if integrate: - wavelength = na.stack( - arrays=[ - wavelength.min(axis_wavelength), - wavelength.max(axis_wavelength), - ], - axis=axis_wavelength, - ) - return self.sensor.measure( rays=rayfunction.outputs, wavelength=wavelength, axis=(axis_wavelength,) + axis_field + axis_pupil, axis_wavelength=axis_wavelength, noise=noise, + integrate=integrate, + ) + + def backproject( + self, + image: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar], + coordinates: na.SpectralPositionalVectorArray, + *, + axis_wavelength: None | str = None, + axis_field: None | tuple[str, str] = None, + integrate: bool = True, + **kwargs: Any, + ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: + """ + Not implemented: a ray-traced :class:`SequentialSystem` is not a linear + operator, so it has no transpose. Build a + :class:`~optika.systems.LinearSystem` (for example with + :meth:`linearize`) to backproject. + """ + raise NotImplementedError( # pragma: nocover + "`SequentialSystem` does not support backprojection; " + "use a `LinearSystem` instead." ) def distortion( @@ -1268,9 +1284,10 @@ def distortion( # average only the unvignetted rays, falling back to all of the rays # for field points excluded from the fit so that the mean is never - # empty + # empty. The sensor-plane positions are expressed in pixel coordinates, + # so the fitted distortion maps the scene onto the pixel grid. coordinates_sensor = np.mean( - rays.outputs.position.xy, + self.sensor.pixels(rays.outputs.position.xy), axis=axis_pupil, where=unvignetted | ~where, ) @@ -1511,6 +1528,96 @@ def area_effective( axis_wavelength=axis_wavelength, ) + def linearize( + self, + wavelength: None | u.Quantity | na.AbstractScalar = None, + field: None | na.AbstractCartesian2dVectorArray = None, + pupil: None | na.AbstractCartesian2dVectorArray = None, + normalized_field: bool = True, + normalized_pupil: bool = True, + degree: int = 2, + ) -> LinearSystem: + """ + Construct a linear approximation of this system by fitting its + distortion, vignetting, and effective area models. + + The result is an :class:`~optika.systems.LinearSystem`, a fast forward + model which images scenes by conservative regridding instead of + raytracing each one. + + The resulting system's + :attr:`~optika.systems.LinearSystem.field_stop` is left as :obj:`None`; + the field stop is not modeled here, since field points outside it are + excluded when fitting the vignetting model rather than represented as a + falloff. + + Parameters + ---------- + wavelength + The wavelengths at which to sample the system. + If :obj:`None` (the default), ``self.grid_input.wavelength`` + will be used. + field + The field positions at which to sample the system, in either + normalized or physical units. + If :obj:`None` (the default), ``self.grid_input.field`` + will be used. + pupil + The **vertices** of the pupil grid, in either normalized or physical + units. The effective area fit uses these vertices to compute the + pupil cell areas, while the distortion and vignetting fits trace at + the corresponding cell centers. + If :obj:`None` (the default), the effective area fit uses its own + default pupil grid and the distortion and vignetting fits use + ``self.grid_input.pupil``. + normalized_field + A boolean flag indicating whether the `field` parameter is given + in normalized or physical units. + normalized_pupil + A boolean flag indicating whether the `pupil` parameter is given + in normalized or physical units. + degree + The degree of the polynomial distortion and vignetting models. + """ + + # `area_effective` interprets `pupil` as cell vertices (it needs them to + # compute the pupil cell areas), while `distortion` and `vignetting` + # trace at sample points, so give them the cell centers of the vertices. + if pupil is not None: + pupil_centers = pupil.cell_centers(axis=tuple(na.shape(pupil))) + else: + pupil_centers = None + + kwargs = dict( + wavelength=wavelength, + field=field, + normalized_field=normalized_field, + normalized_pupil=normalized_pupil, + ) + # the cosine of the refracted angle at which light strikes the sensor, + # computed the same way as + # :meth:`~optika.sensors.AbstractImagingSensor.collect` and averaged + # over the grid. + rays = self.rayfunction_default.outputs + direction = self.sensor.material.direction_refracted( + wavelength=rays.wavelength, + direction=rays.direction, + n=rays.n, + normal=self.sensor.sag.normal(rays.position), + ) + axis_grid = self.axis_wavelength_ + self.axis_field_ + self.axis_pupil_ + direction = direction.mean( + axis=tuple(ax for ax in axis_grid if ax in na.shape(direction)), + ) + + return LinearSystem( + area_effective=self.area_effective(pupil=pupil, **kwargs), + distortion=self.distortion(pupil=pupil_centers, degree=degree, **kwargs), + sensor=self.sensor, + direction=direction, + vignetting=self.vignetting(pupil=pupil_centers, degree=degree, **kwargs), + ) + def _rayfunction_and_axes( self, wavelength: None | u.Quantity | na.AbstractScalar = None, @@ -2012,7 +2119,7 @@ class SequentialSystem( ) mappable_image = na.plt.pcolormesh( image.inputs.position, - C=image.outputs.value.sum("wavelength"), + C=image.outputs.value, ax=ax[1], ) cbar_0 = fig.colorbar( diff --git a/optika/systems/_sequential_test.py b/optika/systems/_sequential_test.py index bd4fa4ce..306278df 100644 --- a/optika/systems/_sequential_test.py +++ b/optika/systems/_sequential_test.py @@ -10,6 +10,7 @@ class AbstractTestAbstractSequentialSystem( test_mixins.AbstractTestDxfWritable, + test_mixins.AbstractTestPlottable, AbstractTestAbstractSystem, ): def test_object(self, a: optika.systems.AbstractSequentialSystem): @@ -401,6 +402,52 @@ def test_area_effective( assert na.unit(result.area).is_equivalent(u.deg**2) assert np.all(result.area >= 0) + @pytest.mark.parametrize( + argnames="wavelength,field,pupil", + argvalues=[ + ( + None, + None, + None, + ), + ( + na.linspace(500, 600, axis="wavelength", num=3) * u.nm, + na.Cartesian2dVectorLinearSpace( + start=-1, + stop=1, + axis=na.Cartesian2dVectorArray("field_x", "field_y"), + num=5, + ), + na.Cartesian2dVectorLinearSpace( + start=-1, + stop=1, + axis=na.Cartesian2dVectorArray("pupil_x", "pupil_y"), + num=5, + ), + ), + ], + ) + def test_linearize( + self, + a: optika.systems.AbstractSequentialSystem, + wavelength: None | u.Quantity | na.AbstractScalar, + field: None | na.AbstractCartesian2dVectorArray, + pupil: None | na.AbstractCartesian2dVectorArray, + ): + if wavelength is None and not a.axis_wavelength_: + with pytest.raises(ValueError): + a.linearize(wavelength=wavelength, field=field, pupil=pupil) + return + result = a.linearize(wavelength=wavelength, field=field, pupil=pupil) + assert isinstance(result, optika.systems.LinearSystem) + assert isinstance(result.distortion, optika.distortion.AbstractDistortionModel) + assert isinstance(result.vignetting, optika.radiometry.AbstractVignettingModel) + assert isinstance( + result.area_effective, optika.radiometry.AbstractEffectiveAreaModel + ) + assert result.sensor is a.sensor + assert result.field_stop is None + def test_spot_diagram(self, a: optika.systems.AbstractSequentialSystem): fig, axs = a.spot_diagram() assert isinstance(fig, plt.Figure) diff --git a/optika/systems/_systems.py b/optika/systems/_systems.py index 5a779ea2..19da5dd3 100644 --- a/optika/systems/_systems.py +++ b/optika/systems/_systems.py @@ -2,6 +2,7 @@ from typing import Any import abc import dataclasses +import astropy.units as u import named_arrays as na import optika @@ -12,8 +13,6 @@ @dataclasses.dataclass(eq=False, repr=False) class AbstractSystem( - optika.mixins.DxfWritable, - optika.mixins.Plottable, optika.mixins.Printable, optika.mixins.Transformable, optika.mixins.Shaped, @@ -28,18 +27,81 @@ class AbstractSystem( def image( self, scene: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar], + *, + axis_wavelength: None | str = None, + axis_field: None | tuple[str, str] = None, + integrate: bool = True, + noise: bool = True, **kwargs: Any, - ) -> na.SpectralPositionalVectorArray: + ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: """ Forward model of the optical system. - Maps the given spectral radiance of a scene to detector counts. + Maps the given spectral radiance of a scene to the electrons measured + by the sensor. Parameters ---------- scene The spectral radiance of the scene as a function of wavelength and field position. + axis_wavelength + The logical axis of `scene` corresponding to changing wavelength. + axis_field + The logical axes of `scene` corresponding to changing field position. + integrate + Whether to integrate the electrons over wavelength into a single + sensor readout. + noise + Whether to include sensor noise in the result. kwargs Additional keyword arguments used by subclass implementations of this method. """ + + def backproject( + self, + image: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar], + coordinates: na.SpectralPositionalVectorArray, + *, + axis_wavelength: None | str = None, + axis_field: None | tuple[str, str] = None, + integrate: bool = True, + unit: None | u.UnitBase = None, + **kwargs: Any, + ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: + """ + Transpose of the forward model, :meth:`image`. + Maps an image of measured electrons back onto the object plane. + + Backprojection is an optional capability: only systems that are linear + operators support it, so the base implementation raises + :class:`NotImplementedError`. Subclasses that support it (such as + :class:`~optika.systems.LinearSystem`) override this method. + + Parameters + ---------- + image + The detector-plane image of electrons to project onto the object + plane. + coordinates + The vertices of each pixel on the object plane to project onto. + axis_wavelength + The logical axis of `coordinates` corresponding to changing + wavelength. + axis_field + The logical axes of `coordinates` corresponding to changing field + position. + integrate + Whether `image` is a single wavelength-integrated readout. + unit + The unit of the backprojected spectral radiance. + The forward model accepts a scene in either photon or energy units, + so the backprojection is expressed in whichever the caller requests. + If :obj:`None` (the default), the radiance is not converted. + kwargs + Additional keyword arguments used by subclass implementations + of this method. + """ + raise NotImplementedError( # pragma: nocover + f"{type(self).__name__} does not support backprojection." + ) diff --git a/optika/systems/_systems_test.py b/optika/systems/_systems_test.py index 11c63004..3f57b7ed 100644 --- a/optika/systems/_systems_test.py +++ b/optika/systems/_systems_test.py @@ -7,7 +7,6 @@ class AbstractTestAbstractSystem( - test_mixins.AbstractTestPlottable, test_mixins.AbstractTestPrintable, test_mixins.AbstractTestTransformable, test_mixins.AbstractTestShaped, diff --git a/pyproject.toml b/pyproject.toml index 66005d0e..63cc8109 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -15,7 +15,7 @@ classifiers = [ ] dependencies = [ "astropy!=6.1.5", - "named-arrays~=2.1", + "named-arrays~=2.4", "pymupdf", "joblib", "ezdxf",